Skip to content

Note

Click here to download the full example code

Advanced processing

The pynapple package provides a small set of high-level functions that are widely used in systems neuroscience.

  • Discrete correlograms
  • Tuning curves
  • Decoding
  • PETH
  • Randomization

This notebook provides few examples with artificial data.

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)

Discrete correlograms

First let's generate some data. Here we have two neurons recorded together. We can group them in a TsGroup.

ts1 = nap.Ts(t=np.sort(np.random.uniform(0, 1000, 2000)), time_units="s")
ts2 = nap.Ts(t=np.sort(np.random.uniform(0, 1000, 1000)), time_units="s")
epoch = nap.IntervalSet(start=0, end=1000, time_units="s")
ts_group = nap.TsGroup({0: ts1, 1: ts2}, time_support=epoch)

print(ts_group)

Out:

  Index    rate
-------  ------
      0       2
      1       1

First we can compute their autocorrelograms meaning the number of spikes of a neuron observed in a time windows centered around its own spikes. For this we can use the function compute_autocorrelogram. We need to specifiy the binsize and windowsize to bin the spike train.

autocorrs = nap.compute_autocorrelogram(
    group=ts_group, binsize=100, windowsize=1000, time_units="ms", ep=epoch  # ms
)
print(autocorrs, "\n")

Out:

           0     1
-0.9  1.0075  0.99
-0.8  1.0200  0.97
-0.7  0.9650  1.03
-0.6  0.9250  1.05
-0.5  1.0875  0.93
-0.4  1.0650  1.01
-0.3  1.0075  0.88
-0.2  1.0175  1.05
-0.1  1.0300  0.99
 0.0  0.0000  0.00
 0.1  1.0300  0.99
 0.2  1.0175  1.05
 0.3  1.0075  0.88
 0.4  1.0650  1.01
 0.5  1.0875  0.93
 0.6  0.9250  1.05
 0.7  0.9650  1.03
 0.8  1.0200  0.97
 0.9  1.0075  0.99 

The variable autocorrs is a pandas DataFrame with the center of the bins for the index and each columns is a neuron.

Similarly, we can compute crosscorrelograms meaning how many spikes of neuron 1 do I observe whenever neuron 0 fires. Here the function is called compute_crosscorrelogram and takes a TsGroup as well.

crosscorrs = nap.compute_crosscorrelogram(
    group=ts_group, binsize=100, windowsize=1000, time_units="ms"  # ms
)

print(crosscorrs, "\n")

Out:

          0
          1
-0.9  0.950
-0.8  0.995
-0.7  1.190
-0.6  1.110
-0.5  0.900
-0.4  1.015
-0.3  0.995
-0.2  0.935
-0.1  0.985
 0.0  0.950
 0.1  0.975
 0.2  1.105
 0.3  1.130
 0.4  1.035
 0.5  0.955
 0.6  0.955
 0.7  1.015
 0.8  1.020
 0.9  1.015 

Column name (0, 1) is read as cross-correlogram of neuron 0 and 1 with neuron 0 being the reference time.


Peri-Event Time Histogram (PETH)

A second way to examine the relationship between spiking and an event (i.e. stimulus) is to compute a PETH. pynapple uses the function compute_perievent to center spike time around the timestamps of an event within a given window.

stim = nap.Tsd(
    t=np.sort(np.random.uniform(0, 1000, 50)), d=np.random.rand(50), time_units="s"
)

peth0 = nap.compute_perievent(ts1, stim, minmax=(-0.1, 0.2), time_unit="s")

print(peth0)

Out:

Index    rate     ref_times
-------  -------  -----------
0        nan      29.40605
1        nan      63.01218
2        nan      63.0748
3        nan      77.28495
4        3.33333  104.30122
5        nan      117.75764
6        nan      130.47599
7        3.33333  143.45276
8        nan      183.43031
9        nan      184.06554
10       3.33333  186.16187
11       nan      194.2719
12       nan      220.44822
13       3.33333  237.7793
14       nan      243.841
15       nan      307.09021
16       3.33333  320.85718
17       6.66667  342.77955
18       nan      349.432
19       3.33333  379.01537
...      ...      ...
30       nan      606.62292
31       nan      613.98349
32       nan      636.69981
33       nan      641.26858
34       3.33333  672.17483
35       nan      690.24915
36       3.33333  705.10134
37       nan      732.03505
38       nan      761.7472
39       3.33333  762.27815
40       3.33333  762.41848
41       nan      783.45928
42       3.33333  810.42527
43       nan      811.15964
44       3.33333  818.97
45       nan      819.86238
46       nan      853.38755
47       3.33333  912.12068
48       nan      919.21089
49       3.33333  985.27014

It is then easy to create a raster plot around the times of the stimulation event by calling the to_tsd function of pynapple to "flatten" the TsGroup peth0.

mkdocs_gallery_thumbnail_number = 2

plt.figure(figsize=(10, 6))
plt.subplot(211)
plt.plot(np.sum(peth0.count(0.01), 1), linewidth=3, color="red")
plt.xlim(-0.1, 0.2)
plt.ylabel("Count")
plt.axvline(0.0)
plt.subplot(212)
plt.plot(peth0.to_tsd(), "|", markersize=20, color="red", mew=4)
plt.xlabel("Time from stim (s)")
plt.ylabel("Stimulus")
plt.xlim(-0.1, 0.2)
plt.axvline(0.0)

tutorial pynapple process

Out:

<matplotlib.lines.Line2D object at 0x7f5f38c4a110>

The same function can be applied to a group of neurons. In this case, it returns a dict of TsGroup

pethall = nap.compute_perievent(ts_group, stim, minmax=(-0.1, 0.2), time_unit="s")

print(pethall[1])

Out:

Index    rate     ref_times
-------  -------  -----------
0        nan      29.40605
1        nan      63.01218
2        nan      63.0748
3        3.33333  77.28495
4        nan      104.30122
5        6.66667  117.75764
6        nan      130.47599
7        nan      143.45276
8        nan      183.43031
9        3.33333  184.06554
10       nan      186.16187
11       nan      194.2719
12       nan      220.44822
13       nan      237.7793
14       3.33333  243.841
15       3.33333  307.09021
16       nan      320.85718
17       nan      342.77955
18       nan      349.432
19       nan      379.01537
...      ...      ...
30       6.66667  606.62292
31       nan      613.98349
32       nan      636.69981
33       nan      641.26858
34       nan      672.17483
35       3.33333  690.24915
36       nan      705.10134
37       3.33333  732.03505
38       3.33333  761.7472
39       3.33333  762.27815
40       3.33333  762.41848
41       nan      783.45928
42       nan      810.42527
43       nan      811.15964
44       3.33333  818.97
45       nan      819.86238
46       3.33333  853.38755
47       nan      912.12068
48       6.66667  919.21089
49       nan      985.27014

Tuning curves

pynapple can compute 1 dimension tuning curves (for example firing rate as a function of angular direction) and 2 dimension tuning curves ( for example firing rate as a function of position). In both cases, a TsGroup object can be directly passed to the function.

First we will create the 2D features:

dt = 0.1
features = np.vstack((np.cos(np.arange(0, 1000, dt)), np.sin(np.arange(0, 1000, dt)))).T
# features += np.random.randn(features.shape[0], features.shape[1])*0.05
features = nap.TsdFrame(
    t=np.arange(0, 1000, dt),
    d=features,
    time_units="s",
    time_support=epoch,
    columns=["a", "b"],
)

print(features)

plt.figure(figsize=(15, 7))
plt.subplot(121)
plt.plot(features[0:100])
plt.title("Features")
plt.xlabel("Time(s)")
plt.subplot(122)
plt.title("Features")
plt.plot(features["a"][0:100], features["b"][0:100])
plt.xlabel("Feature a")
plt.ylabel("Feature b")

Features, Features

Out:

Time (s)           a         b
----------  --------  --------
0.0          1         0
0.1          0.995     0.09983
0.2          0.98007   0.19867
0.3          0.95534   0.29552
0.4          0.92106   0.38942
0.5          0.87758   0.47943
0.6          0.82534   0.56464
0.7          0.76484   0.64422
0.8          0.69671   0.71736
0.9          0.62161   0.78333
1.0          0.5403    0.84147
1.1          0.4536    0.89121
1.2          0.36236   0.93204
1.3          0.2675    0.96356
1.4          0.16997   0.98545
1.5          0.07074   0.99749
1.6         -0.0292    0.99957
1.7         -0.12884   0.99166
1.8         -0.2272    0.97385
1.9         -0.32329   0.9463
...
998.0        0.51785  -0.85547
998.1        0.60066  -0.7995
998.2        0.67748  -0.73554
998.3        0.74753  -0.66423
998.4        0.81011  -0.58628
998.5        0.86459  -0.50248
998.6        0.91043  -0.41365
998.7        0.94718  -0.3207
998.8        0.97447  -0.22453
998.9        0.99201  -0.12613
999.0        0.99965  -0.02646
999.1        0.9973    0.07347
999.2        0.98498   0.17267
999.3        0.96282   0.27014
999.4        0.93104   0.36491
999.5        0.88996   0.45604
999.6        0.83999   0.54261
999.7        0.78162   0.62375
999.8        0.71544   0.69867
999.9        0.64212   0.7666
dtype: float64, shape: (10000, 2)

Text(732.5909090909089, 0.5, 'Feature b')

Here we call the function compute_2d_tuning_curves. To check the accuracy of the tuning curves, we will display the spikes aligned to the features with the function value_from which assign to each spikes the corresponding feature value for neuron 0.

tcurves2d, binsxy = nap.compute_2d_tuning_curves(
    group=ts_group, features=features, nb_bins=10
)

ts_to_features = ts_group[1].value_from(features)

plt.figure()
plt.plot(ts_to_features["a"], ts_to_features["b"], "o", color="red", markersize=4)
extents = (
    np.min(features["b"]),
    np.max(features["b"]),
    np.min(features["a"]),
    np.max(features["a"]),
)
plt.imshow(tcurves2d[1].T, origin="lower", extent=extents, cmap="viridis")
plt.title("Tuning curve unit 0 2d")
plt.xlabel("feature a")
plt.ylabel("feature b")
plt.grid(False)
plt.show()

Tuning curve unit 0 2d

Out:

/mnt/home/gviejo/pynapple/pynapple/process/tuning_curves.py:223: RuntimeWarning: invalid value encountered in divide
  count = count / occupancy

Decoding

Pynapple supports 1 dimensional and 2 dimensional bayesian decoding. The function returns the decoded feature as well as the probabilities for each timestamps.

First we generate some artificial "place fields" in 2 dimensions based on the features.

This part is just to generate units with a relationship to the features (i.e. "place fields")

times = features.as_units("us").index.values
ft = features.values
alpha = np.arctan2(ft[:, 1], ft[:, 0])
bins = np.repeat(np.linspace(-np.pi, np.pi, 13)[::, np.newaxis], 2, 1)
bins += np.array([-2 * np.pi / 24, 2 * np.pi / 24])
ts_group = {}
for i in range(12):
    ts = times[(alpha >= bins[i, 0]) & (alpha <= bins[i + 1, 1])]
    ts_group[i] = nap.Ts(ts, time_units="us")

ts_group = nap.TsGroup(ts_group, time_support=epoch)
print(ts_group)

Out:

  Index    rate
-------  ------
      0   1.248
      1   1.665
      2   1.666
      3   1.664
      4   1.665
      5   1.67
      6   1.671
      7   1.673
      8   1.665
      9   1.666
     10   1.667
     11   1.248

To decode we need to compute tuning curves in 2D.

import warnings

warnings.filterwarnings("ignore")

tcurves2d, binsxy = nap.compute_2d_tuning_curves(
    group=ts_group,
    features=features,
    nb_bins=10,
    ep=epoch,
    minmax=(-1.0, 1.0, -1.0, 1.0),
)

Then we plot the "place fields".

plt.figure(figsize=(20, 9))
for i in ts_group.keys():
    plt.subplot(2, 6, i + 1)
    plt.imshow(
        tcurves2d[i], extent=(binsxy[1][0], binsxy[1][-1], binsxy[0][0], binsxy[0][-1])
    )
    plt.xticks()
plt.show()

tutorial pynapple process

Then we call the actual decoding function in 2d.

decoded, proba_feature = nap.decode_2d(
    tuning_curves=tcurves2d,
    group=ts_group,
    ep=epoch,
    bin_size=0.1,  # second
    xy=binsxy,
    features=features,
)


plt.figure(figsize=(15, 5))
plt.subplot(131)
plt.plot(features["a"].as_units("s").loc[0:20], label="True")
plt.plot(decoded["a"].as_units("s").loc[0:20], label="Decoded")
plt.legend()
plt.xlabel("Time (s)")
plt.ylabel("Feature a")
plt.subplot(132)
plt.plot(features["b"].as_units("s").loc[0:20], label="True")
plt.plot(decoded["b"].as_units("s").loc[0:20], label="Decoded")
plt.legend()
plt.xlabel("Time (s)")
plt.ylabel("Feature b")
plt.subplot(133)
plt.plot(
    features["a"].as_units("s").loc[0:20],
    features["b"].as_units("s").loc[0:20],
    label="True",
)
plt.plot(
    decoded["a"].as_units("s").loc[0:20],
    decoded["b"].as_units("s").loc[0:20],
    label="Decoded",
)
plt.xlabel("Feature a")
plt.ylabel("Feature b")
plt.legend()
plt.tight_layout()
plt.show()

tutorial pynapple process


Randomization

Pynapple provides some ready-to-use randomization methods to compute null distributions for statistical testing. Different methods preserve or destroy different features of the data, here's a brief overview.

shift_timestamps shifts all the timestamps in a Ts object by the same random amount, wrapping the end of the time support to its beginning. This randomization preserves the temporal structure in the data but destroys the temporal relationships with other quantities (e.g. behavioural data). When applied on a TsGroup object, each series in the group is shifted independently.

ts = nap.Ts(t=np.sort(np.random.uniform(0, 100, 10)), time_units="ms")
rand_ts = nap.shift_timestamps(ts, min_shift=1, max_shift=20)

shuffle_ts_intervals computes the intervals between consecutive timestamps, permutes them, and generates a new set of timestamps with the permuted intervals. This procedure preserve the distribution of intervals, but not their sequence.

ts = nap.Ts(t=np.sort(np.random.uniform(0, 100, 10)), time_units="s")
rand_ts = nap.shuffle_ts_intervals(ts)

jitter_timestamps shifts each timestamp in the data of an independent random amount. When applied with a small max_jitter, this procedure destroys the fine temporal structure of the data, while preserving structure on longer timescales.

ts = nap.Ts(t=np.sort(np.random.uniform(0, 100, 10)), time_units="s")
rand_ts = nap.jitter_timestamps(ts, max_jitter=1)

resample_timestamps uniformly re-draws the same number of timestamps in ts, in the same time support. This procedures preserve the total number of timestamps, but destroys any other feature of the original data.

ts = nap.Ts(t=np.sort(np.random.uniform(0, 100, 10)), time_units="s")
rand_ts = nap.resample_timestamps(ts)

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

Download Python source code: tutorial_pynapple_process.py

Download Jupyter notebook: tutorial_pynapple_process.ipynb

Gallery generated by mkdocs-gallery