Skip to content

Note

Click here to download the full example code

Quick start

The examplar data to replicate the figure in the jupyter notebook can be found here.

The data contains a sample recordings taken simultaneously from the anterodorsal thalamus and the hippocampus and contains both a sleep and wake session. It contains both head-direction cells (i.e. cells that fire for a particular head direction in the horizontal plane) and place cells (i.e. cells that fire for a particular position in the environment).

Preprocessing of the data was made with Kilosort 2.0 and spike sorting was made with Klusters.

Instructions for installing pynapple can be found here.


This notebook is meant to provide an overview of pynapple by going through:

  • Input output (IO). In this case, pynapple will load a NWB file using the NWBFile object within a project Folder that represent a dataset.
  • Core functions that handle time series, interval sets and groups of time series. See this notebook for a detailled usage of the core functions.
  • Process functions. A small collection of high-level functions widely used in system neuroscience. This notebook details those functions.

Warning

This tutorial uses seaborn and matplotlib for displaying the figure.

You can install both with pip install matplotlib seaborn

import numpy as np
import pandas as pd
import pynapple as nap
import matplotlib.pyplot as plt
import seaborn as sns

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

IO

The first step is to give the path to the data folder.

DATA_DIRECTORY = "../../your/path/to/MyProject/"

We can load the session with the function load_folder. Pynapple will walks throught the folder and collects every subfolders. We can use the attribute view or the function expand to display a tree view of the dataset. The treeview shows all the compatible data format (i.e npz files or NWBs files) and their equivalent pynapple type.

data = nap.load_folder(DATA_DIRECTORY)
data.view

Out:

๐Ÿ“‚ MyProject
โ””โ”€โ”€ ๐Ÿ“‚ sub-A2929
    โ””โ”€โ”€ ๐Ÿ“‚ A2929-200711
        โ”œโ”€โ”€ ๐Ÿ“‚ derivatives
        โ”‚   โ”œโ”€โ”€ spikes.npz      |        TsGroup
        โ”‚   โ”œโ”€โ”€ sleep_ep.npz    |        IntervalSet
        โ”‚   โ”œโ”€โ”€ position.npz    |        TsdFrame
        โ”‚   โ””โ”€โ”€ wake_ep.npz     |        IntervalSet
        โ”œโ”€โ”€ ๐Ÿ“‚ pynapplenwb
        โ”‚   โ””โ”€โ”€ A2929-200711    |        NWB file
        โ”œโ”€โ”€ x_plus_1.npz    |        Tsd
        โ””โ”€โ”€ stimulus-fish.npz       |        IntervalSet

The object data is a Folder object that allows easy navigation and interaction with a dataset. In this case, we want to load the NWB file in the folder /pynapplenwb. Data are always lazy loaded. No time series is loaded until it's actually called. When calling the NWB file, the object nwb is an interface to the NWB file. All the data inside the NWB file that are compatible with one of the pynapple objects are shown with their corresponding keys.

nwb = data["sub-A2929"]["A2929-200711"]["pynapplenwb"]["A2929-200711"]
print(nwb)

Out:

A2929-200711
โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”ฏโ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”‘
โ”‚ Keys                  โ”‚ Type        โ”‚
โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”ฟโ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”ฅ
โ”‚ units                 โ”‚ TsGroup     โ”‚
โ”‚ position_time_support โ”‚ IntervalSet โ”‚
โ”‚ epochs                โ”‚ IntervalSet โ”‚
โ”‚ z                     โ”‚ Tsd         โ”‚
โ”‚ y                     โ”‚ Tsd         โ”‚
โ”‚ x                     โ”‚ Tsd         โ”‚
โ”‚ rz                    โ”‚ Tsd         โ”‚
โ”‚ ry                    โ”‚ Tsd         โ”‚
โ”‚ rx                    โ”‚ Tsd         โ”‚
โ”•โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”ทโ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”โ”™

We can individually call each object and they are actually loaded.

units is a TsGroup object. It allows to group together time series with different timestamps and couple metainformation to each neuron. In this case, the location of where the neuron was recorded has been added when loading the session for the first time. We load units as spikes

spikes = nwb["units"]
print(spikes)

Out:

  Index      rate  location      group
-------  --------  ----------  -------
      0   7.30358  adn               0
      1   5.73269  adn               0
      2   8.11944  adn               0
      3   6.67856  adn               0
      4  10.7712   adn               0
      5  11.0045   adn               0
      6  16.5164   adn               0
      7   2.19674  ca1               1
      8   2.0159   ca1               1
      9   1.06837  ca1               1
     10   3.91847  ca1               1
     11   3.30844  ca1               1
     12   1.0942   ca1               1
     13   1.27921  ca1               1
     14   1.32338  ca1               1

In this case, the TsGroup holds 15 neurons and it is possible to access, similar to a dictionary, the spike times of a single neuron:

neuron_0 = spikes[0]
print(neuron_0)

Out:

Time (s)
0.00845
0.03265
0.1323
0.3034
0.329
0.661
0.671
0.70995
1.5596
1.83285
1.8503
1.8595
1.88455
1.9147
1.9598
2.03655
2.088
2.108
2.12555
2.4989
...
1185.266
1185.27545
1185.28575
1185.29245
1185.30375
1185.31475
1185.3214
1185.33505
1185.34715
1185.3705
1185.3779
1185.3826
1185.39065
1185.4003
1185.84315
1186.12755
1189.384
1194.13475
1196.2075
1196.67675
shape: 8764

neuron_0 is a Ts object containing the times of the spikes.

The other information about the session is contained in nwb["epochs"]. In this case, the start and end of the sleep and wake epochs. If the NWB time intervals contains tags of the epochs, pynapple will try to group them together and return a dictionary of IntervalSet instead of IntervalSet.

epochs = nwb["epochs"]
print(epochs)

Out:

{'sleep':             start    end
       0        0    600
shape: (1, 2), time unit: sec., 'wake':             start    end
       0      600   1200
shape: (1, 2), time unit: sec.}

Finally this dataset contains tracking of the animal in the environment. rx, ry, rz represent respectively the roll, the yaw and the pitch of the head of the animal. x and z represent the position of the animal in the horizontal plane while y represents the elevation. Here we load only the head-direction as a Tsd object.

head_direction = nwb["ry"]
print(head_direction)

Out:

Time (s)
----------  -------
670.6407    5.20715
670.649     5.18103
670.65735   5.15551
670.66565   5.13654
670.674     5.12085
670.68235   5.10845
670.69065   5.09732
670.699     5.08613
670.70735   5.07296
670.71565   5.06046
670.724     5.04562
670.73235   5.02894
670.74065   5.00941
670.749     4.9876
670.75735   4.96359
670.76565   4.94159
670.774     4.92152
670.78235   4.90635
670.79065   4.89627
670.799     4.89221
...
1199.8366   4.29808
1199.84495  4.2702
1199.85325  4.24382
1199.8616   4.20729
1199.86995  4.16671
1199.87825  4.12178
1199.8866   4.07674
1199.89495  4.02607
1199.90325  3.97859
1199.9116   3.92706
1199.91995  3.87269
1199.92825  3.82848
1199.9366   3.78868
1199.94495  3.74765
1199.95325  3.70803
1199.9616   3.66595
1199.96995  3.63462
1199.97825  3.61785
1199.9866   3.60945
1199.99495  3.60938
dtype: float32, shape: (63527,)

Core

The core functions of pynapple provides many ways to manipulate time series. In this example, spike times are restricted to the wake epoch. Notice how the frequencies change from the original object.

wake_ep = epochs["wake"]

spikes_wake = spikes.restrict(wake_ep)
print(spikes_wake)

Out:

  Index      rate  location      group
-------  --------  ----------  -------
      0   4.84667  adn               0
      1   8.06333  adn               0
      2   7.10667  adn               0
      3   7.66333  adn               0
      4   7.96833  adn               0
      5  11.285    adn               0
      6  22.0833   adn               0
      7   1.82     ca1               1
      8   2.83667  ca1               1
      9   0.705    ca1               1
     10   4.775    ca1               1
     11   4.93     ca1               1
     12   1.71     ca1               1
     13   0.96833  ca1               1
     14   0.26167  ca1               1

The same operation can be applied to all time series.

# In this example, we want all the epochs for which position in `x` is above a certain threhsold. For this we use the function `threshold`.
posx = nwb["x"]

threshold = 0.08

posxpositive = posx.threshold(threshold)

plt.figure()
plt.plot(posx)
plt.plot(posxpositive, ".")
plt.axhline(threshold)
plt.xlabel("Time (s)")
plt.ylabel("x")
plt.title("x > {}".format(threshold))
plt.tight_layout()
plt.show()

x > 0.08

The epochs above the threshold can be accessed through the time support of the Tsd object. The time support is an important concept in the pynapple package. It helps the user to define the epochs for which the time serie should be defined. By default, Ts, Tsd and TsGroup objects possess a time support (defined as an IntervalSet). It is recommended to pass the time support when instantiating one of those objects.

epochs_above_thr = posxpositive.time_support
print(epochs_above_thr)

Out:

            start      end
       0  682.661  745.566
       1  752.24   752.44
       2  752.582  752.674
       3  757.498  758.998
       4  789.863  790.272
       5  875.225  876.067
       6  878.158  878.642
shape: (7, 2), time unit: sec.

Tuning curves

Let's do a more advanced analysis. Neurons from ADn (group 0 in the spikes group object) are know to fire for a particular direction. Therefore, we can compute their tuning curves, i.e. their firing rates as a function of the head-direction of the animal in the horizontal plane (ry). To do this, we can use the function compute_1d_tuning_curves. In this case, the tuning curves are computed over 120 bins and between 0 and 2$\pi$.

tuning_curves = nap.compute_1d_tuning_curves(
    group=spikes, feature=head_direction, nb_bins=121, minmax=(0, 2 * np.pi)
)

print(tuning_curves)

Out:

                 0    1   ...        13        14
0.025964  45.520459  0.0  ...  0.000000  2.069112
0.077891  55.049762  0.0  ...  0.000000  0.000000
0.129818  76.369034  0.0  ...  1.558552  0.000000
0.181745  82.179721  0.0  ...  2.608880  0.000000
0.233672  73.851374  0.0  ...  2.637549  0.000000
...             ...  ...  ...       ...       ...
6.049513  15.001060  0.0  ...  0.000000  1.363733
6.101440  22.327159  0.0  ...  0.000000  0.000000
6.153367  47.062150  0.0  ...  2.353107  0.000000
6.205295  56.003958  0.0  ...  0.000000  0.000000
6.257222  38.712414  0.0  ...  0.000000  0.000000

[121 rows x 15 columns]

We can plot tuning curves in polar plots.

neuron_location = spikes.get_info("location")  # to know where the neuron was recorded
plt.figure(figsize=(12, 9))

for i, n in enumerate(tuning_curves.columns):
    plt.subplot(3, 5, i + 1, projection="polar")
    plt.plot(tuning_curves[n])
    plt.title(neuron_location[n] + "-" + str(n), fontsize=18)

plt.tight_layout()
plt.show()

adn-0, adn-1, adn-2, adn-3, adn-4, adn-5, adn-6, ca1-7, ca1-8, ca1-9, ca1-10, ca1-11, ca1-12, ca1-13, ca1-14

While ADN neurons show obvious modulation for head-direction, it is not obvious for all CA1 cells. Therefore we want to restrict the remaining of the analyses to only ADN neurons. We can split the spikes group with the function getby_category.

spikes_by_location = spikes.getby_category("location")

print(spikes_by_location["adn"])
print(spikes_by_location["ca1"])

spikes_adn = spikes_by_location["adn"]

Out:

  Index      rate  location      group
-------  --------  ----------  -------
      0   7.30358  adn               0
      1   5.73269  adn               0
      2   8.11944  adn               0
      3   6.67856  adn               0
      4  10.7712   adn               0
      5  11.0045   adn               0
      6  16.5164   adn               0
  Index     rate  location      group
-------  -------  ----------  -------
      7  2.19674  ca1               1
      8  2.0159   ca1               1
      9  1.06837  ca1               1
     10  3.91847  ca1               1
     11  3.30844  ca1               1
     12  1.0942   ca1               1
     13  1.27921  ca1               1
     14  1.32338  ca1               1

Correlograms

A classical question with head-direction cells is how pairs stay coordinated across brain states i.e. wake vs sleep (see Peyrache, A., Lacroix, M. M., Petersen, P. C., & Buzsรกki, G. (2015). Internally organized mechanisms of the head direction sense. Nature neuroscience, 18(4), 569-575.)

In this example, this coordination across brain states will be evaluated with cross-correlograms of pairs of neurons. We can call the function compute_crosscorrelogram during both sleep and wake epochs.

cc_wake = nap.compute_crosscorrelogram(
    group=spikes_adn,
    binsize=20,  # ms
    windowsize=4000,  # ms
    ep=epochs["wake"],
    norm=True,
    time_units="ms",
)

cc_sleep = nap.compute_crosscorrelogram(
    group=spikes_adn,
    binsize=5,  # ms
    windowsize=400,  # ms
    ep=epochs["sleep"],
    norm=True,
    time_units="ms",
)

From the previous figure, we can see that neurons 0 and 1 fires for opposite directions during wake. Therefore we expect their cross-correlograms to show a trough around 0 time lag, meaning those two neurons do not fire spikes together. A similar trough during sleep for the same pair thus indicates a persistence of their coordination even if the animal is not moving its head. mkdocs_gallery_thumbnail_number = 3

xtwake = cc_wake.index.values
xtsleep = cc_sleep.index.values

plt.figure(figsize=(15, 5))
plt.subplot(131, projection="polar")
plt.plot(tuning_curves[[0, 1]])  # The tuning curves of the pair [0,1]
plt.subplot(132)
plt.fill_between(
    xtwake, np.zeros_like(xtwake), cc_wake[(0, 1)].values, color="darkgray"
)
plt.title("wake")
plt.xlabel("Time (ms)")
plt.ylabel("CC")
plt.subplot(133)
plt.fill_between(
    xtsleep, np.zeros_like(xtsleep), cc_sleep[(0, 1)].values, color="lightgrey"
)
plt.title("sleep")
plt.xlabel("Time (ms)")
plt.ylabel("CC")
plt.tight_layout()
plt.show()

wake, sleep


Decoding

This last analysis shows how to use the pynapple's decoding function.

The previous result indicates a persistent coordination of head-direction cells during sleep. Therefore it is possible to decode a virtual head-direction signal even if the animal is not moving its head. This example uses the function decode_1d which implements bayesian decoding (see : Zhang, K., Ginzburg, I., McNaughton, B. L., & Sejnowski, T. J. (1998). Interpreting neuronal population activity by reconstruction: unified framework with application to hippocampal place cells. Journal of neurophysiology, 79(2), 1017-1044.)

First we can validate the decoding function with the real position of the head of the animal during wake.

tuning_curves_adn = nap.compute_1d_tuning_curves(
    spikes_adn, head_direction, nb_bins=61, minmax=(0, 2 * np.pi)
)

decoded, proba_angle = nap.decode_1d(
    tuning_curves=tuning_curves_adn,
    group=spikes_adn,
    ep=epochs["wake"],
    bin_size=0.3,  # second
    feature=head_direction,
)
print(decoded)

Out:

Time (s)
----------  -------
600.15      2.11156
600.45      2.11156
600.75      2.31757
601.05      2.00856
601.35      2.11156
601.65      2.11156
601.95      2.21457
602.25      2.21457
602.55      2.21457
602.85      2.21457
603.15      2.21457
603.45      2.42057
603.75      3.03859
604.05      3.03859
604.35      3.2446
604.65      3.4506
604.95      3.2446
605.25      3.2446
605.55      3.3476
605.85      2.93559
...
1194.15     4.89264
1194.45     4.78964
1194.75     4.48063
1195.05     4.48063
1195.35     4.48063
1195.65     4.99565
1195.95     4.99565
1196.25     4.89264
1196.55     4.89264
1196.85     4.78964
1197.15     4.48063
1197.45     4.27463
1197.75     4.27463
1198.05     4.58364
1198.35     4.99565
1198.65     4.58364
1198.95     4.27463
1199.25     4.58364
1199.55     4.58364
1199.85     3.86261
dtype: float64, shape: (2000,)

We can plot the decoded head-direction along with the true head-direction.

plt.figure(figsize=(20, 5))
plt.plot(head_direction.as_units("s"), label="True")
plt.plot(decoded.as_units("s"), label="Decoded")
plt.legend()
plt.xlabel("Time (s)")
plt.ylabel("Head-direction (rad)")
plt.show()

tutorial pynapple quick start


Raster

Finally we can decode activity during sleep and overlay spiking activity of ADN neurons as a raster plot (in this case only during the first 4 seconds). Pynapple return as well the probability of being in a particular state. We can display it next to the spike train.

First let's decode during sleep with a bin size of 40 ms.

decoded_sleep, proba_angle_Sleep = nap.decode_1d(
    tuning_curves=tuning_curves_adn,
    group=spikes_adn,
    ep=epochs["sleep"],
    bin_size=0.04,  # second
    feature=head_direction,
)

Here we are gonna chain the TsGroup function set_info and the function to_tsd to flatten the TsGroup and quickly assign to each spikes a corresponding value found in the metadata table. Any columns of the metadata table can be assigned to timestamps in a TsGroup.

Here the value assign to the spikes comes from the preferred firing direction of the neurons. The following line is a quick way to sort the neurons based on their preferred firing direction

order = np.argsort(np.argmax(tuning_curves_adn.values, 0))
print(order)

Out:

[0 4 2 6 1 3 5]

Assigning order as a metadata of TsGroup

spikes_adn.set_info(order=order)
print(spikes_adn)

Out:

  Index      rate  location      group    order
-------  --------  ----------  -------  -------
      0   7.30358  adn               0        0
      1   5.73269  adn               0        4
      2   8.11944  adn               0        2
      3   6.67856  adn               0        6
      4  10.7712   adn               0        1
      5  11.0045   adn               0        3
      6  16.5164   adn               0        5

"Flattening" the TsGroup to a Tsd based on order. It's then very easy to call plot on tsd_adn to display the raster

tsd_adn = spikes_adn.to_tsd("order")
print(tsd_adn)

Out:

Time (s)
----------  --
0.00845      0
0.03265      0
0.07745      6
0.1323       0
0.14045      5
0.1511       5
0.16955      5
0.17905      5
0.1935       5
0.19925      5
0.2039       5
0.21795      5
0.2286       5
0.23705      5
0.2524       5
0.2781       5
0.29455      5
0.3034       0
0.329        0
0.51505      2
...
1199.74275   4
1199.743     6
1199.74975   4
1199.7549    6
1199.761     4
1199.7736    6
1199.7756    4
1199.788     4
1199.79995   4
1199.81095   4
1199.8219    4
1199.83955   6
1199.8453    4
1199.856     4
1199.8897    4
1199.9065    5
1199.91745   5
1199.94065   5
1199.95035   5
1199.96795   5
dtype: float64, shape: (79349,)

Plotting everything

subep = nap.IntervalSet(start=0, end=10, time_units="s")

plt.figure(figsize=(19, 10))
plt.subplot(211)
plt.plot(tsd_adn.restrict(subep), "|", markersize=20)
plt.xlim(subep.start[0], subep.end[0])
plt.ylabel("Order")
plt.title("Decoding during sleep")
plt.subplot(212)
p = proba_angle_Sleep.restrict(subep)
plt.imshow(p.values.T, aspect="auto", origin="lower", cmap="viridis")
plt.title("Probability")
plt.xticks([0, p.shape[0] - 1], subep.values[0])
plt.yticks([0, p.shape[1]], ["0", "360"])
plt.legend()
plt.xlabel("Time (s)")
plt.ylabel("Head-direction (deg)")
plt.legend()
plt.show()

Decoding during sleep, Probability

Out:

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

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

Download Python source code: tutorial_pynapple_quick_start.py

Download Jupyter notebook: tutorial_pynapple_quick_start.ipynb

Gallery generated by mkdocs-gallery