Skip to content

Note

Click here to download the full example code

Peyrache et al (2015) Tutorial

This tutorial demonstrates how we use Pynapple to generate Figure 4a in the publication. The NWB file for the example is hosted on OSF. We show below how to stream it. The entire dataset can be downloaded here.

See the documentation of Pynapple for instructions on installing the package.

This tutorial was made by Dhruv Mehrotra and Guillaume Viejo.

Warning

This tutorial uses seaborn and matplotlib for displaying the figure

You can install all with pip install matplotlib seaborn tqdm

mkdocs_gallery_thumbnail_number = 2

Now, import the necessary libraries:

import numpy as np
import pandas as pd
import pynapple as nap
import scipy.ndimage
import matplotlib.pyplot as plt
import seaborn as sns
import requests, math, os
import tqdm

custom_params = {"axes.spines.right": False, "axes.spines.top": False}
sns.set_theme(style="ticks", palette="colorblind", font_scale=1.5, rc=custom_params)

Downloading the data

It's a small NWB file.

path = "Mouse32-140822.nwb"
if path not in os.listdir("."):
    r = requests.get(f"https://osf.io/jb2gd/download", stream=True)
    block_size = 1024*1024
    with open(path, 'wb') as f:
        for data in tqdm.tqdm(r.iter_content(block_size), unit='MB', unit_scale=True,
            total=math.ceil(int(r.headers.get('content-length', 0))//block_size)):
            f.write(data)

Parsing the data

The first step is to load the data and other relevant variables of interest

data = nap.load_file(path)  # Load the NWB file for this dataset

What does this look like ?

print(data)

Out:

Mouse32-140822
┍━━━━━━━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━┑
 Keys                   Type        ┝━━━━━━━━━━━━━━━━━━━━━━━┿━━━━━━━━━━━━━┥
 units                  TsGroup      sws                    IntervalSet  rem                    IntervalSet  position_time_support  IntervalSet  epochs                 IntervalSet  ry                     Tsd         ┕━━━━━━━━━━━━━━━━━━━━━━━┷━━━━━━━━━━━━━┙

Head-Direction Tuning Curves

To plot Head-Direction Tuning curves, we need the spike timings and the orientation of the animal. These quantities are stored in the variables 'units' and 'ry'.

spikes = data["units"]  # Get spike timings
epochs = data["epochs"]  # Get the behavioural epochs (in this case, sleep and wakefulness)
angle = data["ry"]  # Get the tracked orientation of the animal

What does this look like ?

print(spikes)

Out:

Index    rate     location    group
-------  -------  ----------  -------
0        2.96981  thalamus    1
1        2.42638  thalamus    1
2        5.93417  thalamus    1
3        5.04432  thalamus    1
...      ...      ...         ...
45       0.55768  thalamus    6
46       1.15056  thalamus    6
47       0.46084  thalamus    6
48       0.19287  thalamus    7

Here, rate is the mean firing rate of the unit. Location indicates the brain region the unit was recorded from, and group refers to the shank number on which the cell was located.

This dataset contains units recorded from the anterior thalamus. Head-direction (HD) cells are found in the anterodorsal nucleus of the thalamus (henceforth referred to as ADn). Units were also recorded from nearby thalamic nuclei in this animal. For the purposes of our tutorial, we are interested in the units recorded in ADn. We can restrict ourselves to analysis of these units rather easily, using Pynapple.

spikes_adn = spikes.getby_category("location")["adn"]  # Select only those units that are in ADn

What does this look like ?

print(spikes_adn)

Out:

Index    rate      location    group
-------  --------  ----------  -------
4        0.30207   adn         2
5        0.87042   adn         2
6        0.36154   adn         2
7        10.51737  adn         2
...      ...       ...         ...
31       0.58829   adn         4
32       1.12899   adn         4
33       5.26316   adn         4
34       1.57122   adn         4

Let's compute some head-direction tuning curves. To do this in Pynapple, all you need is a single line of code!

Plot firing rate of ADn units as a function of heading direction, i.e. a head-direction tuning curve

tuning_curves = nap.compute_1d_tuning_curves(
    group=spikes_adn, 
    feature=angle, 
    nb_bins=61, 
    ep = epochs['wake'],
    minmax=(0, 2 * np.pi)
    )

What does this look like ?

print(tuning_curves)

Out:

                4         5         6          7         8         9         10        11        12         13         14        15        16        17  ...        21         22        23        24        25         26        27        28        29        30        31        32         33        34
0.051502  0.255172  0.127586  0.170115   5.316080  0.000000  0.722987  5.231023  0.467815  0.680458   9.313772   8.675842  0.000000  0.000000  1.190802  ...  0.680458   6.081595  1.233331  0.765516  0.722987   2.636776  0.000000  0.000000  0.297700  0.850573  0.552872  0.467815   5.273551  1.786203
0.154505  0.300635  0.000000  0.187897   8.305042  0.037579  0.338214  3.194247  0.338214  0.751587   8.718415  14.731115  0.000000  0.000000  1.127381  ...  0.150317   4.321628  0.638849  0.413373  0.338214   1.691072  0.000000  0.000000  0.375794  0.751587  1.014643  0.375794   8.380200  1.954127
0.257508  0.189885  0.094943  0.094943  11.867814  0.158238  0.094943  3.671110  0.316475  0.569655   9.430956  22.406433  0.000000  0.000000  0.886130  ...  0.063295   3.386283  0.949425  0.443065  0.632950   1.107663  0.031648  0.000000  0.348123  1.076015  1.677318  0.158238  10.506971  3.892643
0.360511  0.498062  0.052428  0.078641  16.724387  0.052428  0.314565  4.534983  0.340779  0.655344  10.170944  26.921545  0.052428  0.026214  1.494185  ...  0.183496   2.411667  0.812627  1.599040  0.498062   1.992247  0.078641  0.000000  0.366993  0.629131  1.782537  0.235924  12.556397  5.321396
0.463514  0.362941  0.111674  0.139593  20.631824  0.055837  0.335023  4.411134  0.279186  0.725883  11.502451  34.674865  0.027919  0.027919  0.642127  ...  0.251267   1.619277  0.307104  0.195430  0.335023   0.530453  0.055837  0.000000  0.418779  0.642127  2.456834  0.307104  15.271458  8.431408
...            ...       ...       ...        ...       ...       ...       ...       ...       ...        ...        ...       ...       ...       ...  ...       ...        ...       ...       ...       ...        ...       ...       ...       ...       ...       ...       ...        ...       ...
5.819672  0.063460  0.158650  0.190380   1.491313  0.793252  0.507681  5.584491  0.444221  1.554773   2.284565   1.269203  0.000000  0.000000  1.935534  ...  0.158650  18.339976  7.297914  4.283558  1.015362  12.565105  0.095190  0.000000  0.095190  2.221104  0.095190  0.253841   0.285571  0.507681
5.922675  0.024772  0.123861  0.123861   1.337703  0.569762  0.445901  4.558100  0.371584  0.990891   3.146080   1.015664  0.024772  0.000000  1.040436  ...  0.099089  17.984677  5.648080  1.213842  1.040436  10.032774  0.049545  0.074317  0.099089  1.535882  0.049545  0.396357   0.767941  0.272495
6.025678  0.000000  0.112276  0.028069   1.936754  0.533309  0.224551  4.406818  0.392965  1.263101   5.333091   1.571858  0.000000  0.028069  1.178894  ...  0.140345  13.978313  4.013853  2.189374  0.926274   6.427778  0.000000  0.084207  0.112276  1.038549  0.084207  0.364896   0.954343  0.196482
6.128681  0.000000  0.138009  0.165611   2.152940  0.220814  0.386425  4.223075  0.441629  1.076470   4.637101   1.932126  0.000000  0.027602  1.518099  ...  0.248416  13.773295  3.864251  3.119003  0.634841   6.265607  0.055204  0.027602  0.027602  0.938461  0.055204  0.331222   1.021266  0.441629
6.231684  0.067699  0.101548  0.135397   2.504850  0.000000  0.304644  3.791124  0.372343  0.541589   5.720536   3.147987  0.033849  0.000000  0.947781  ...  0.169247  10.256345  2.030959  0.304644  0.473891   3.892672  0.000000  0.033849  0.101548  0.812384  0.304644  0.541589   2.234055  0.609288

[61 rows x 31 columns]

Each row indicates an angular bin (in radians), and each column corresponds to a single unit. Let's compute the preferred angle quickly as follows:

pref_ang = tuning_curves.idxmax()

For easier visualization, we will colour our plots according to the preferred angle of the cell. To do so, we will normalize the range of angles we have, over a colourmap.

norm = plt.Normalize()  # Normalizes data into the range [0,1]
color = plt.cm.hsv(norm([i / (2 * np.pi) for i in pref_ang.values]))  # Assigns a colour in the HSV colourmap for each value of preferred angle
color = pd.DataFrame(index=pref_ang.index, data = color, columns = ['r', 'g', 'b', 'a'])

To make the tuning curves look nice, we will smooth them before plotting, using this custom function:

from scipy.ndimage import gaussian_filter1d
def smoothAngularTuningCurves(tuning_curves, sigma=2):

    tmp = np.concatenate((tuning_curves.values, tuning_curves.values, tuning_curves.values))
    tmp = gaussian_filter1d(tmp, sigma=sigma, axis=0)

    return pd.DataFrame(index = tuning_curves.index,
        data = tmp[tuning_curves.shape[0]:tuning_curves.shape[0]*2], 
        columns = tuning_curves.columns
        )

Therefore, we have:

smoothcurves = smoothAngularTuningCurves(tuning_curves, sigma=3)

What does this look like? Let's plot the tuning curves!

plt.figure(figsize=(12, 9))
for i, n in enumerate(pref_ang.sort_values().index.values):
    plt.subplot(8, 4, i + 1, projection='polar')  # Plot the curves in 8 rows and 4 columns
    plt.plot(
        smoothcurves[n], color=color.loc[n]
    )  # Colour of the curves determined by preferred angle    
    plt.xlabel("Angle (rad)")  # Angle in radian, on the X-axis
    plt.ylabel("Firing Rate (Hz)")  # Firing rate in Hz, on the Y-axis
    plt.xticks([])
plt.show()

tutorial HD dataset

Awesome!


Decoding

Now that we have HD tuning curves, we can go one step further. Using only the population activity of ADn units, we can decode the direction the animal is looking in. We will then compare this to the real head direction of the animal, and discover that population activity in the ADn indeed codes for HD.

To decode the population activity, we will be using a Bayesian Decoder as implemented in Pynapple. Just a single line of code!

decoded, proba_feature = nap.decode_1d(
    tuning_curves=tuning_curves,
    group=spikes_adn,
    ep=epochs["wake"],
    bin_size=0.1,  # second
    feature=angle,
)

What does this look like ?

print(decoded)

Out:

Time (s)
----------  --------
8812.35     0.772523
8812.45     0.66952
8812.55     0.463514
8812.65     0.66952
...
10770.95    5.92267
10771.05    5.92267
10771.15    5.92267
10771.25    5.92267
dtype: float64, shape: (19590,)

The variable 'decoded' indicates the most probable angle in which the animal was looking. There is another variable, 'proba_feature' that denotes the probability of a given angular bin at a given time point. We can look at it below:

print(proba_feature.as_dataframe())

Out:

              0.051502  0.154505      0.257508      0.360511      0.463514      0.566517      0.669520  0.772523      0.875526      0.978529  1.081532  ...      5.201653      5.304656      5.407659      5.510663      5.613666      5.716669      5.819672      5.922675      6.025678      6.128681      6.231684
8812.35   2.199077e-06  0.000223  3.717901e-03  1.769861e-02  7.205101e-02  1.209866e-01  1.601997e-01  0.312027  1.813195e-01  4.539781e-02  0.063692  ...  6.562027e-18  2.549308e-19  3.093103e-17  0.000000e+00  4.690809e-19  5.917759e-19  5.337070e-18  6.116697e-18  2.510509e-15  1.483408e-14  2.212419e-11
8812.45   8.561129e-08  0.000013  1.858017e-03  9.312326e-03  9.216920e-02  2.697128e-01  3.057500e-01  0.186276  1.164765e-01  9.345216e-03  0.008323  ...  1.019887e-19  8.619404e-22  1.694103e-18  6.042584e-21  6.914020e-21  1.348545e-20  4.031744e-20  2.920013e-20  1.609456e-18  1.628816e-16  2.464448e-13
8812.55   4.168300e-04  0.022715  1.916560e-01  1.919596e-01  2.489159e-01  1.441986e-01  1.160732e-01  0.055778  2.033150e-02  2.540693e-03  0.004145  ...  7.409675e-11  7.261910e-11  1.416640e-09  8.292585e-11  1.325661e-11  2.676119e-11  9.353997e-12  9.577405e-12  2.723566e-10  2.698994e-09  4.014847e-07
8812.65   1.082000e-05  0.000156  5.663501e-03  1.973657e-02  1.097107e-01  1.954540e-01  2.926565e-01  0.210710  1.295982e-01  1.561975e-02  0.018847  ...  1.755181e-17  1.724868e-19  1.374135e-16  3.578488e-18  4.528014e-18  1.448351e-17  1.244645e-16  1.160040e-15  2.430940e-14  6.175847e-13  4.426551e-10
8812.75   4.128198e-05  0.001369  2.088684e-02  3.927845e-02  1.338599e-01  1.742158e-01  4.097762e-01  0.170674  2.982752e-02  6.882122e-03  0.012688  ...  3.160851e-14  2.523405e-14  0.000000e+00  0.000000e+00  8.439360e-15  0.000000e+00  2.246921e-14  7.808235e-14  1.152593e-12  5.077673e-12  5.290852e-09
...                ...       ...           ...           ...           ...           ...           ...       ...           ...           ...       ...  ...           ...           ...           ...           ...           ...           ...           ...           ...           ...           ...           ...
10770.85  6.695624e-05  0.000003  1.968111e-07  1.057744e-07  5.481207e-11  2.365351e-14  9.349059e-18  0.000000  1.616283e-12  3.996825e-11  0.000000  ...  7.933647e-03  3.851417e-02  5.724365e-02  1.567618e-01  1.082532e-01  2.254218e-01  1.287529e-01  1.980321e-01  3.692355e-02  3.310523e-02  4.842028e-03
10770.95  2.924858e-04  0.000005  1.392026e-07  4.800910e-08  1.638920e-10  6.896839e-13  2.673364e-16  0.000000  1.276731e-12  3.111488e-11  0.000000  ...  4.482674e-03  1.260161e-02  2.855182e-02  1.289702e-01  9.672991e-02  1.741568e-01  1.688514e-01  2.808207e-01  3.751799e-02  5.668237e-02  8.429181e-03
10771.05  1.093979e-03  0.000115  1.300380e-05  3.998319e-06  1.464971e-08  5.651362e-11  1.128912e-12  0.000000  5.683691e-10  1.866787e-09  0.000000  ...  7.879675e-03  3.988051e-02  4.264079e-02  1.568604e-01  8.787573e-02  1.375497e-01  1.069793e-01  2.313362e-01  9.806256e-02  6.114471e-02  2.121581e-02
10771.15  5.537065e-03  0.001235  1.198087e-04  2.325507e-05  4.497914e-07  6.727892e-09  1.151520e-10  0.000000  6.206958e-09  1.253643e-08  0.000000  ...  2.924244e-03  1.250249e-02  1.941843e-02  5.342305e-02  4.553747e-02  1.063931e-01  7.923923e-02  2.470057e-01  1.578827e-01  1.547459e-01  1.115693e-01
10771.25  5.969857e-04  0.000058  2.970409e-06  2.588196e-06  3.024315e-09  3.310300e-12  3.821833e-15  0.000000  8.289259e-11  8.503706e-10  0.000000  ...  9.996361e-03  2.610720e-02  3.732318e-02  7.769811e-02  1.118871e-01  1.793572e-01  1.789009e-01  2.164274e-01  7.182830e-02  6.686085e-02  1.501718e-02

[19590 rows x 61 columns]

Each row of this pandas DataFrame is a time bin, and each column is an angular bin. The sum of all values in a row add up to 1.

Now, let's plot the raster plot for a given period of time, and overlay the actual and decoded HD on the population activity.

ep = nap.IntervalSet(
    start=10717, end=10730
)  # Select an arbitrary interval for plotting

plt.figure()
plt.rc("font", size=12)
for i, n in enumerate(spikes_adn.keys()):
    plt.plot(
        spikes[n].restrict(ep).fillna(pref_ang[n]), "|", color=color.loc[n]
    )  # raster plot for each cell
plt.plot(
    decoded.restrict(ep), "--", color="grey", linewidth=2, label="decoded HD"
)  # decoded HD
plt.legend(loc="upper left")

tutorial HD dataset

Out:

<matplotlib.legend.Legend object at 0x7f294cd2e9d0>

From this plot, we can see that the decoder is able to estimate the head-direction based on the population activity in ADn. Amazing!

What does the probability distribution in this example event look like? Ideally, the bins with the highest probability will correspond to the bins having the most spikes. Let's plot the probability matrix to visualize this.

smoothed = scipy.ndimage.gaussian_filter(
    proba_feature, 1
)  # Smoothening the probability distribution

# Create a DataFrame with the smoothed distribution
p_feature = pd.DataFrame(
    index=proba_feature.index.values,
    columns=proba_feature.columns.values,
    data=smoothed,
)
p_feature = nap.TsdFrame(p_feature)  # Make it a Pynapple TsdFrame

plt.figure()
plt.plot(
    angle.restrict(ep), "w", linewidth=2, label="actual HD", zorder=1
)  # Actual HD, in white
plt.plot(
    decoded.restrict(ep), "--", color="grey", linewidth=2, label="decoded HD", zorder=1
)  # Decoded HD, in grey

# Plot the smoothed probability distribution
plt.imshow(
    np.transpose(p_feature.restrict(ep).values),
    aspect="auto",
    interpolation="bilinear",
    extent=[ep["start"][0], ep["end"][0], 0, 2 * np.pi],
    origin="lower",
    cmap="viridis",
)

plt.xlabel("Time (s)")  # X-axis is time in seconds
plt.ylabel("Angle (rad)")  # Y-axis is the angle in radian
plt.colorbar(label="probability")

tutorial HD dataset

Out:

<matplotlib.colorbar.Colorbar object at 0x7f294c9961d0>

From this probability distribution, we observe that the decoded HD very closely matches the actual HD. Therefore, the population activity in ADn is a reliable estimate of the heading direction of the animal.

I hope this tutorial was helpful. If you have any questions, comments or suggestions, please feel free to reach out to the Pynapple Team!

Total running time of the script: ( 0 minutes 2.799 seconds)

Download Python source code: tutorial_HD_dataset.py

Download Jupyter notebook: tutorial_HD_dataset.ipynb

Gallery generated by mkdocs-gallery