Neural analysis with Pynapple

This chapter is an introduction to Pynapple.
Pynapple is a python package for analysing neural data.
The package is actively maintained by the Flatiron Institute’s neuroRSE team.

Pynapple can do a lot more than what is covered in this chapter, take a look at the Pynapple documentation for more.
The only thing you’ll need to follow along is a working Python environment with pynapple installed. LINK TO PREVIOUS CHAPTERS.

Open your favourite way of running Python, and follow along!

Origins

2000: TSToolbox, MATLAB
McNaughton lab: David Redish & Francesco Battaglia
2016: TSToolbox2, MATLAB
Adrien Peyrache & Luke Sjulson
2018: neuroseries, Python
Francesco Battaglia
2021: pynapple, Python
Guillaume Viejo

presently with the Flatiron Institute’s neuroRSE team:

Guillaume Viejo
Sarah Jo Venditto
Edoardo Balzani
William Broderick
Wolf De Wulf, a 2025 neuroRSE intern

Goal

Pynapple was designed to lie in between pre- and postprocessing.
It contains functions facilitating alignment, wrangling, and performing basic neuroscientific analysis.

Recording

Processing
segmenting (calman)
spike-sorting (spikeinterface)

Analysis
model fitting (nemos)
interpretation (you)

1. Pynapple fundamentals

Time series (`Ts`)

We store a time series, for example a spike train, in a Ts object.

import pynapple as nap
ts = nap.Ts(t=[1,2,3,4,5])
ts

Time (s)
1.0
2.0
3.0
4.0
5.0
shape: 5

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(constrained_layout=True, figsize=(5,3))
for timepoint in ts.times():
    ax.axvline(timepoint)
ax.yaxis.set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_xlabel("time (s)")
plt.show()

Time series with data (`Tsd`)

We store a time series with corresponding data values, for example the speed of a mouse when it is running around, in a Tsd object.

import pynapple as nap
import numpy as np
tsd = nap.Tsd(t=[1,2,3,4,5], d=[1,1,1,2,1])
tsd

Time (s)
----------  --
1            1
2            1
3            1
4            2
5            1
dtype: int64, shape: (5,)

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(constrained_layout=True, figsize=(5,3))
ax.plot(tsd)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_xlabel("time (s)")
ax.set_ylabel("data (a.u.)")
plt.show()

Multiple time series

If we have multiple time series, they can be collected in a TsGroup.
This could be a collection of spike times for multiple spiking units, for example.

import pynapple as nap
tsgroup = nap.TsGroup({
    1: nap.Ts([1,2,3,4,5]), 
    2: nap.Ts([1.5, 2.2, 2.9, 4.2])
})
tsgroup

  Index    rate
-------  ------
      1    1.25
      2    1

Show plotting code

import matplotlib.pyplot as plt
import numpy as np
cmap = plt.get_cmap("tab10")
fig, ax = plt.subplots(constrained_layout=True, figsize=(5,3))
for i in tsgroup:
    ax.vlines(tsgroup[i].times(), i, i+1, label=i, color=cmap(i-1))
ax.yaxis.set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_xlabel("time (s)")
plt.ylim(1,3)
plt.legend(frameon=False)
plt.show()

Time series with multiple data columns (`TsdFrame`)

If we have a time series with multiple data associated with the same time points, they can be stored in a TsdFrame.
This could be a collection of activity traces from multiple cells, e.g. recorded using calcium-imaging.

import pynapple as nap
tsdframe = nap.TsdFrame(
    t=[1,2,3,4,5], 
    d=[[1,2],
       [1,2],
       [1,3],
       [2,1],
       [1,2]],
    columns=["A", "B"]
)
tsdframe

Time (s)      A    B
----------  ---  ---
1             1    2
2             1    2
3             1    3
4             2    1
5             1    2
dtype: int64, shape: (5, 2)

Show plotting code

import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots(constrained_layout=True, figsize=(5,3))
for i, col in enumerate(tsdframe.columns):
    ax.plot(tsdframe[:, i], label=col)
ax.set_ylabel("data (a.u.)")
ax.set_xlabel("time (s)")
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.legend(frameon=False)
plt.show()

Intervals (`IntervalSet`)

Often we want to store events that occur across a period rather than at single time points, for example a behavioural stimulus that lasts multiple seconds, or a period of sleep. To represent this, we store time intervals as start and end pairs in an IntervalSet.

import pynapple as nap
epochs = nap.IntervalSet(
    start=[1, 3, 7],
    end=[2, 5, 9]
)
epochs

  index    start    end
      0        1      2
      1        3      5
      2        7      9
shape: (3, 2), time unit: sec.

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(constrained_layout=True, figsize=(5,3))
for start, stop in epochs.values:
    ax.axvspan(start, stop)
ax.set_xlabel("time (s)")
ax.spines['left'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()

2. Interactions

A key goal of Pynapple is to facilitate integrating neural activity and behavioural data for processing and analysis. Therefore, many Pynapple objects interact with each other.

Time supports (`time_support`)

Every time series object has a time support, an IntervalSet that gives the epochs over which the time series is defined.

import pynapple as nap
ts = nap.Ts(t=[1,2,3,4,5])
ts.time_support

  index    start    end
      0        1      5
shape: (1, 2), time unit: sec.

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(constrained_layout=True, figsize=(5,3))
for timepoint in ts.times():
    ax.axvline(timepoint)
for start, stop in ts.time_support.values:
    ax.axvspan(start, stop, alpha=0.2)
ax.yaxis.set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_xlabel("time (s)")
plt.show()

Restricting (`restrict`)

Any time series object can be restricted to certain IntervalSet, thereby updating its time support. This could be used, for example, to keep only the spike times that fall within the presentation of a behavioural stimulus.

import pynapple as nap
ts = nap.Ts(t=[1,2,3,4,5])
restriction = nap.IntervalSet(start=[1.3, 3.9], end=[3.5, 4.2])
restricted = ts.restrict(restriction)
restricted

Time (s)
2.0
3.0
4.0
shape: 3

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(constrained_layout=True, figsize=(5,3))
for timepoint in restricted.times():
    ax.axvline(timepoint)
for start, stop in restricted.time_support.values:
    ax.axvspan(start, stop, alpha=0.2)
ax.yaxis.set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_xlabel("time (s)")
ax.set_xlim(1,5)
plt.show()

Taking values (`value_from`)

Given one time series, you can value_from to take the value corresponding to the closest time point in another time series with data. This can be handy, for example, if you want to sample the position of an animal, at particular spike times.

import pynapple as nap
import numpy as np
ts = nap.Ts(t=[1,2,3,4,5])
tsd = nap.Tsd(t=np.arange(0, 10, 0.1), d=np.random.randn(100))
values = ts.value_from(tsd)
values

Time (s)
----------  ---------
1            0.483223
2           -0.423219
3           -0.446724
4            1.455
5            0.883494
dtype: float64, shape: (5,)

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(constrained_layout=True, figsize=(5,3))
ax.plot(tsd, label="tsd")
ax.plot(values, label="values", linestyle="none", marker=".")
minimum = tsd.values.min()
cmap = plt.get_cmap("tab10")
for timepoint, value in zip(values.times(), values.values, strict=True):
    ax.vlines(timepoint, ymin=minimum, ymax=value, color=cmap(1))
ax.set_ylabel("data (a.u.)")
ax.set_xlabel("time (s)")
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
fig.legend(frameon=False, bbox_to_anchor=(1.0, 1.1))
plt.show()

3. Wrangling

Pynapple objects are easy to transform, reshape, slice, and more.

Numpy

Most Numpy functions can be applied to Pynapple objects without hassle.

import pynapple as nap
import numpy as np
tsdframe = nap.TsdFrame(
    t=[1,2,3,4,5], 
    d=[[1,2],
       [1,2],
       [1,3],
       [2,1],
       [1,2]],
    columns=["A", "B"]
)
np.mean(tsdframe, axis=1)

Time (s)
----------  ---
1           1.5
2           1.5
3           2
4           1.5
5           1.5
dtype: float64, shape: (5,)

import pynapple as nap
import numpy as np
tsdframe = nap.TsdFrame(
    t=[1,2,3,4,5], 
    d=[[1,2],
       [1,2],
       [1,3],
       [2,1],
       [1,2]],
    columns=["A", "B"]
)
np.diff(tsdframe, axis=1)

Time (s)      0
----------  ---
1             1
2             1
3             2
4            -1
5             1
dtype: int64, shape: (5, 1)

Slicing

Pynapple objects can be sliced like Numpy arrays.

import pynapple as nap
import numpy as np
tsdframe = nap.TsdFrame(
    t=[1,2,3,4,5], 
    d=[[1,2],
       [1,2],
       [1,3],
       [2,1],
       [1,2]],
    columns=["A", "B"]
)
tsdframe[:, 1]

Time (s)
----------  --
1            2
2            2
3            3
4            1
5            2
dtype: int64, shape: (5,)

import pynapple as nap
import numpy as np
tsdframe = nap.TsdFrame(
    t=[1,2,3,4,5], 
    d=[[1,2],
       [1,2],
       [1,3],
       [2,1],
       [1,2]],
    columns=["A", "B"]
)
tsdframe[2:4]

Time (s)      A    B
----------  ---  ---
3             1    3
4             2    1
dtype: int64, shape: (2, 2)

4. Core functions

Pynapple has a core set of functions designed to help with neuroscientific analysis. The following are a subset of them, for an exhaustive list, see the Pynapple documentation.

Binning

The time points of any time series object can be counted.

import pynapple as nap
import numpy as np
time_step_s = 0.001
times = np.arange(0, 4, time_step_s)

# Create a firing pattern that increases and decreases as a sine wave
rate = 40 * (1 + np.sin(2*np.pi*times)) / 2

# Use the rate to create discrete spike times
spikes = nap.Ts(times[np.random.rand(len(times)) < rate * time_step_s])

# Bin the spike times
counts = spikes.count(bin_size=0.1)
counts

Time (s)
------------------  --
0.079                2
0.179                4
0.279                6
0.379                3
0.479                2
0.5790000000000001   2
0.679                0
...
3.379                4
3.4789999999999996   1
3.5789999999999997   0
3.679                0
3.779                0
3.879                0
3.9789999999999996   1
dtype: int64, shape: (40,)

Show plotting code

import matplotlib.pyplot as plt
fig, (ax_rate, ax_spikes, ax_count) = plt.subplots(
    3, 1, constrained_layout=True, figsize=(7.5,3), sharex=True
)
ax_rate.plot(times, rate)
ax_rate.spines['right'].set_visible(False)
ax_rate.spines['top'].set_visible(False)
ax_rate.set_ylabel("true rate [Hz]")
ax_spikes.vlines(spikes.times(), 0, 1)
ax_spikes.set_ylabel("spikes")
ax_spikes.set_yticks([])
ax_spikes.spines['left'].set_visible(False)
ax_spikes.spines['right'].set_visible(False)
ax_spikes.spines['top'].set_visible(False)
ax_count.plot(counts)
ax_count.set_ylabel("count")
ax_count.set_xlabel("time (s)")
ax_count.spines['right'].set_visible(False)
ax_count.spines['top'].set_visible(False)
plt.show()

Averaging

If you have continous data already, you would not count events in bins, but rather average values across bins. For example, if you wanted to downsample a signal without losing too much information.

import pynapple as nap
import numpy as np
time_step_s = 0.01
times = np.arange(0, 10, time_step_s)

# Create a sine wave signal
signal = np.sin(2 * np.pi * times) 

# Add random noise
signal += 0.2 * np.random.randn(len(times))

# Store as Tsd
tsd = nap.Tsd(t=times, d=signal)

# Average across 200ms bins
averaged = tsd.bin_average(bin_size=0.2)
averaged

Time (s)
-------------------  -----------
0.1                   0.572503
0.30000000000000004   0.88994
0.5                  -0.00437328
0.7                  -0.83483
0.9                  -0.716906
1.1                   0.506269
1.3                   0.901277
...
8.7                  -0.937301
8.9                  -0.593816
9.1                   0.465485
9.299999999999999     0.972677
9.5                   0.0892741
9.7                  -0.855876
9.9                  -0.573498
dtype: float64, shape: (50,)

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(constrained_layout=True, figsize=(5,3))
plt.plot(tsd, label="original")
plt.plot(averaged, label="averaged")
ax.set_ylabel("data (a.u.)")
ax.set_xlabel("time (s)")
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
fig.legend(frameon=False, bbox_to_anchor=(1.0, 1.1))
plt.show()

Thresholding

Time series with data can be thresholded, yielding new time series objects with updated time supports.

import pynapple as nap
import numpy as np
time_step_s = 0.01
times = np.arange(0, 10, time_step_s)

# Create a sine wave signal
signal = np.sin(2 * np.pi * times) 

# Add random noise
signal += 0.2 * np.random.randn(len(times))

# Store as Tsd
tsd = nap.Tsd(t=times, d=signal)

# Threshold
above = tsd.threshold(0.5, method="above")
below = tsd.threshold(-.5, method="below")
above.time_support

index    start                end
0        0.065                0.07500000000000001
1        0.08499999999999999  0.185
2        0.195                0.365
3        0.385                0.41500000000000004
4        1.0550000000000002   1.0750000000000002
5        1.0950000000000002   1.4249999999999998
6        2.075                2.375
...      ...                  ...
35       7.445                7.455
36       8.045000000000002    8.065000000000001
37       8.075                8.405000000000001
38       8.415                8.425
39       9.045000000000002    9.055
40       9.085                9.415
41       9.425                9.434999999999999
shape: (42, 2), time unit: sec.

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(constrained_layout=True, figsize=(5,3))
ax.plot(tsd, label="original")
ax.plot(above, label="above", marker=".", linestyle="none")
ax.plot(below, label="below", marker=".", linestyle="none")
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
fig.legend(frameon=False, bbox_to_anchor=(1.0, 1.2))
plt.show()

Smoothing

Time series with data can be smoothed. Smoothing makes a signal less noisy by averaging out small, fast changes while keeping the main pattern. The operation consists of convolving the signal with a Gaussian kernel.

import pynapple as nap
import numpy as np
time_step_s = 0.01
times = np.arange(0, 10, time_step_s)

# Create a sine wave signal
signal = np.sin(2 * np.pi * times) 

# Add random noise
signal += 0.2 * np.random.randn(len(times))

# Store as Tsd
tsd = nap.Tsd(t=times, d=signal)

# Gaussian smooth with a standard deviation of 100ms
smoothed = tsd.smooth(std=0.1)
smoothed

Time (s)
-----------------  ---------
0.0                 0.197228
0.01                0.223642
0.02                0.251785
0.03                0.281496
0.04                0.312576
0.05                0.34479
0.06                0.377874
...
9.93               -0.356862
9.94               -0.325914
9.950000000000001  -0.295973
9.96               -0.26724
9.97               -0.239888
9.98               -0.214056
9.99               -0.18985
dtype: float64, shape: (1000,)

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(constrained_layout=True, figsize=(5,3))
ax.plot(tsd, label="original")
ax.plot(smoothed, label="smoothed")
ax.set_ylabel("data (a.u.)")
ax.set_xlabel("time (s)")
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
fig.legend(frameon=False, bbox_to_anchor=(1.0, 1.1))
plt.show()

5. Advanced Functions

Beyond fundamental processing functions, Pynapple has loads of advanced functionality for neuroscientific analysis.

Autocorrelograms

Given one or a group of time series, you can compute autocorrelograms.
LINK TO PREVIOUS CHAPTERS.

import pynapple as nap
import numpy as np
time_step_s = 0.01
times = np.arange(0, 4, time_step_s)

# Create a firing pattern that increases and decreases as a sine wave
rate = 40 * (1 + np.sin(2*np.pi*times)) / 2

# Use the rate to create 2 spiking units, one with double the rate
units = nap.TsGroup({
    1: nap.Ts(times[np.random.rand(len(times)) < rate * time_step_s]),
    2: nap.Ts(times[np.random.rand(len(times)) < rate*2 * time_step_s]),
})

# Compute autocorrelograms
autocorrelogram = nap.compute_autocorrelogram(
    units, binsize=0.1, windowsize=1.0, norm=False
)
autocorrelogram

	1	2
-0.9	20.000000	42.603550
-0.8	17.037037	36.035503
-0.7	13.086420	30.414201
-0.6	10.123457	25.029586
-0.5	7.654321	23.136095
-0.4	10.370370	27.278107
-0.3	18.024691	36.272189
-0.2	21.728395	46.390533
-0.1	26.913580	55.443787
0.0	0.000000	0.000000
0.1	27.283951	55.621302
0.2	21.111111	46.153846
0.3	17.777778	36.035503
0.4	9.876543	27.159763
0.5	8.024691	22.662722
0.6	10.246914	24.911243
0.7	12.469136	29.585799
0.8	17.037037	35.621302
0.9	20.000000	42.721893

Show plotting code

import matplotlib.pyplot as plt
ax=autocorrelogram.plot()
plt.ylabel("rate (Hz)")
plt.xlabel("time (s)")
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.legend(frameon=False)
plt.show()

Crosscorrelograms

Given one or two groups of time series, you can also compute crosscorrelograms.

crosscorrelogram = nap.compute_crosscorrelogram(
    units, binsize=0.1, windowsize=1.0, norm=False
)
crosscorrelogram

	1
	2
-0.9	43.456790
-0.8	40.493827
-0.7	33.086420
-0.6	25.925926
-0.5	22.098765
-0.4	24.691358
-0.3	31.481481
-0.2	42.098765
-0.1	54.074074
0.0	60.493827
0.1	58.024691
0.2	48.518519
0.3	39.753086
0.4	28.765432
0.5	19.506173
0.6	19.259259
0.7	23.456790
0.8	33.086420
0.9	41.234568

Show plotting code

import matplotlib.pyplot as plt
ax=crosscorrelogram.plot()
plt.ylabel("rate (Hz)")
plt.xlabel("time (s)")
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.legend(frameon=False)
plt.show()

Signal processing

Pynapple wraps a lot of scipy functions for signal processing.

Let’s simulate a noisy multi-frequency signal (the combination of multiple sine waves oscillating at different frequencies):

import pynapple as nap
import numpy as np
sampling_rate_hz = 1000
times = np.linspace(0, 2, sampling_rate_hz * 2)

# Signals at different frequencies
signal_2Hz = np.cos(times*2*np.pi*2)
signal_10Hz = np.cos(times*2*np.pi*10)
signal_50Hz = np.cos(times*2*np.pi*50)

# Combine
signal = signal_2Hz+signal_10Hz+signal_50Hz

# Add noise
signal += np.random.normal(0, 0.5, len(times))

# Store in Tsd
tsd = nap.Tsd(t=times,d=signal)
tsd

Time (s)
---------------------  --------
0.0                    2.88351
0.0010005002501250625  2.50478
0.002001000500250125   3.3138
0.003001500750375187   2.98422
0.00400200100050025    2.77529
0.0050025012506253125  1.65997
0.006003001500750374   1.10145
...
1.9939969984992494     0.836853
1.9949974987493746     1.36542
1.9959979989994996     1.11617
1.9969984992496246     1.86567
1.9979989994997498     1.78296
1.9989994997498748     3.51441
2.0                    3.19231
dtype: float64, shape: (2000,)

Show plotting code

import matplotlib.pyplot as plt
fig, (ax2, ax10, ax50, ax) = plt.subplots(
    4,
    1,
    constrained_layout=True, 
    figsize=(5,6), 
    sharex=True
)
ax.plot(tsd)
ax2.plot(times,signal_2Hz, color="red")
ax2.spines['right'].set_visible(False)
ax2.spines['top'].set_visible(False)
ax2.set_ylabel("2Hz")
ax10.plot(times,signal_10Hz, color="red")
ax10.spines['right'].set_visible(False)
ax10.spines['top'].set_visible(False)
ax10.set_ylabel("10Hz")
ax50.plot(times,signal_50Hz, color="red")
ax50.spines['right'].set_visible(False)
ax50.spines['top'].set_visible(False)
ax50.set_ylabel("10Hz")
ax.set_ylabel("mixed signal (a.u.)")
ax.set_xlabel("time (s)")
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()

Bandpass filter

We can bandpass filter to extract certain frequency bands. For example, if we apply a bandpass filter with a 8Hz high-pass and 12Hz low-pass cutoff, we will keep only the part of the signal in the 8-12hz range.

filtered = nap.apply_bandpass_filter(
    data=tsd, 
    cutoff=(8, 12), 
    fs=sampling_rate_hz, 
    mode='butter'
)
filtered

Time (s)
---------------------  ----------
0.0                    0.26773
0.0010005002501250625  0.248751
0.002001000500250125   0.228765
0.003001500750375187   0.207815
0.00400200100050025    0.185952
0.0050025012506253125  0.163229
0.006003001500750374   0.139702
...
1.9939969984992494     0.00309984
1.9949974987493746     0.00279932
1.9959979989994996     0.00251099
1.9969984992496246     0.00223674
1.9979989994997498     0.00197809
1.9989994997498748     0.00173619
2.0                    0.00151186
dtype: float64, shape: (2000,)

Show plotting code

import matplotlib.pyplot as plt
fig, (ax_original, ax_filtered) = plt.subplots(
    2,
    1,
    constrained_layout=True, 
    figsize=(5,4), 
    sharex=True
)
ax_original.plot(tsd, color="red")
ax_original.set_ylabel("original (a.u.)")
ax_original.spines['right'].set_visible(False)
ax_original.spines['top'].set_visible(False)
ax_filtered.plot(filtered)
ax_filtered.set_xlabel("time (s)")
ax_filtered.set_ylabel("filtered (a.u.)")
ax_filtered.spines['right'].set_visible(False)
ax_filtered.spines['top'].set_visible(False)
plt.show()

Power-spectral density

We can compute the power-spectral density (PSD) to find frequencies of interest.
The PSD shows how strongly different frequencies contribute to a signal. You only ever compute it up to half the sampling rate (the Nyquist frequency), because higher frequencies cannot be reliably recovered from the recorded data.

psd = nap.compute_power_spectral_density(
    tsd, 
    fs=sampling_rate_hz
)
psd

	0
0.0	0.000005
0.5	0.000639
1.0	0.000033
1.5	0.000120
2.0	0.994257
...	...
497.5	0.000093
498.0	0.000155
498.5	0.000129
499.0	0.000934
499.5	0.000065

1000 rows × 1 columns

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(constrained_layout=True, figsize=(5,3))
ax.plot(psd)
ax.set_ylabel("power")
ax.set_xlabel("frequency (Hz)")
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()

Perievent

A typical neuroscientific analysis consists of aligning activity to a set of stimuli or events. For example, spiking activity around a stimulus presentation.

Let’s start by simulating a spiking unit that has clear stimulus-locked activity:

import pynapple as nap
import numpy as np
stimuli = nap.Ts(t=np.arange(0, 1000, 1), time_units="s")
baseline = np.random.uniform(0, 1000, 500)
burst = np.concatenate([
        np.random.normal(
            st + 0.1, 0.05, 3
        ) 
        for st in stimuli.times()
    ])
ts = nap.Ts(t=np.sort(np.concatenate([baseline, burst])))
ts

Time (s)
0.06696263198522495
0.11055034497211144
0.17988360058853112
0.5023373119573682
1.1244560582839225
1.1259295365660025
1.136160497284933
...
998.127294612097
998.1735347565713
999.0794552630998
999.1260268678623
999.162846014479
999.4782084278808
999.9329704556361
shape: 3500

Show plotting code

import matplotlib.pyplot as plt
segment = nap.IntervalSet(100, 103.9)
fig, ax = plt.subplots(1, 1, constrained_layout=True)
ax.vlines(ts.restrict(segment).times(), 0.04, 0.10, label="spikes")
ax.vlines(
    stimuli.restrict(segment).times(),
    0.0,
    0.14,
    color="gray",
    linestyle="--",
    label="stimulus",
)
ax.yaxis.set_visible(False)
ax.spines["left"].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_xlabel("time (s)")
plt.show()

Now, we can easily compute a peri-event time histogram (PETH) as follows:

peth = nap.compute_perievent(
    data=ts, 
    events=stimuli, 
    window=(-0.1, 0.4))
peth

Index    rate    events
-------  ------  --------
0        6.0     1.0
1        6.0     2.0
2        6.0     3.0
3        6.0     4.0
4        6.0     5.0
5        6.0     6.0
6        6.0     7.0
...      ...     ...
992      8.0     993.0
993      8.0     994.0
994      6.0     995.0
995      6.0     996.0
996      6.0     997.0
997      6.0     998.0
998      6.0     999.0

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(peth.to_tsd(), "|", markersize=5)
ax.set_ylabel("event")
ax.set_xlabel("time from event (s)")
ax.axvline(0.0, color="gray", linestyle="--")
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()

A PETH aligns all activity (spikes) with respect to the event times.

Tuning curves

A typical neuroscientific analysis is that of computing tuning curves: the mean activity curve with respect to a variable of interest. Well known examples are the response of neurons in the visual cortex with respect to the orientation of a bar stimulus, the response of neurons in the postsubiculum with respect to head direction, or the response of hippocampal neurons with respect to position on a track.

Tuning curves can be computed from both discrete (spike times) and continuous activity. Pynapple’s compute_tuning_curves function handles both cases, depending on whether you pass a Ts or Tsd as data. In what follows, we’ll simulate computing tuning curves from a circular variable (e.g. bar angle, or head direction).

Let’s start by simulating some spiking units that are clearly modulated by a circular feature:

import pandas as pd
import pynapple as nap
import numpy as np
from scipy.ndimage import gaussian_filter1d

n_neurons = 6

# Divide the circular feature space (0 to 2π) into bins
bins = np.linspace(0, 2*np.pi, 61)

# Create a Gaussian-shaped tuning curve centred at 0
# This represents a neuron that fires most strongly for one preferred value
x = np.linspace(-np.pi, np.pi, len(bins)-1)
tmp = np.roll(
    np.exp(-(1.5 * x) ** 2),  # Gaussian bump
    (len(bins)-1)//2          # Shift peak to the centre of the array
)

# Create tuning curves for all neurons by shifting the preferred value
# Each neuron responds to a different part of the circular space
tc = np.array([
    np.roll(tmp, i * (len(bins)-1)//n_neurons)
    for i in range(n_neurons)
]).T

# Create a time vector
total_timesteps = 10000
time_step_sec = 0.01
timestep = np.arange(total_timesteps) * time_step_sec

# Simulate a random trajectory through the circular feature space
# The cumulative sum creates a smooth random walk
# Gaussian filtering removes abrupt jumps
# Modulo 2π keeps the feature within the circular range [0, 2π)
feature = nap.Tsd(
    t=timestep,
    d=(
        gaussian_filter1d(
            np.cumsum(np.random.randn(total_timesteps) * 0.5),
            20
        ) % (2*np.pi)
    )
)

# Find which tuning curve bin corresponds to each feature value
index = np.digitize(feature, bins) - 1

# Look up each neuron's activity at the current feature value
# Then sample spikes randomly based on this activity level
count = np.random.poisson(tc[index]) > 0

tsgroup = nap.TsGroup({
    i + 1: nap.Ts(timestep[count[:, i]])
    for i in range(n_neurons)
})
epochs = nap.IntervalSet(0, 10)

Let’s visualise what we generated.
We now have a random walk in the feature space:

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(constrained_layout=True, figsize=(6, 4))
ax.plot(feature.restrict(epochs), linestyle="none", marker="o")
ax.set_yticks([0, 2*np.pi], ["0", "2π"])
ax.set_ylim(0, 2*np.pi)
ax.set_xlabel("time(s)")
ax.set_ylabel("feature (rad)")
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()

As well as the activity of 6 neurons, modulated by the feature:

Show plotting code

import matplotlib.pyplot as plt
fig, ax = plt.subplots(constrained_layout=True, figsize=(6, 4))
ax.plot(tsgroup.restrict(epochs).to_tsd(), linestyle="none", marker="|")
ax.tick_params(axis="y", length=0)
ax.set_xlabel("time(s)")
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()

We can now use Pynapple to easily compute tuning curves with respect to the circular feature.
The output is an xarray.DataArray:

tuning_curves = nap.compute_tuning_curves(
    data=tsgroup,
    features=feature,
    bins=120, 
    range=(0, 2*np.pi),
    feature_names=["feature"]
    )
tuning_curves

<xarray.DataArray (unit: 6, feature: 120)> Size: 6kB
array([[66.15384615, 52.63157895, 66.66666667, 60.        , 62.5       ,
        68.23529412, 60.21505376, 52.74725275, 48.86363636, 48.75      ,
        34.24657534, 32.32323232, 27.18446602, 25.        , 20.45454545,
        26.73267327, 18.18181818, 10.46511628,  5.19480519,  3.61445783,
         6.59340659,  2.98507463,  7.57575758,  6.09756098,  1.13636364,
         1.51515152,  5.12820513,  1.2987013 ,  2.7027027 ,  1.33333333,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         1.25      ,  1.26582278,  1.44927536,  0.        ,  3.17460317,
         2.98507463,  2.53164557,  1.61290323,  3.22580645,  6.94444444,
...
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.80645161,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  1.0989011 ,
         0.93457944,  0.        ,  1.03092784,  5.10204082,  3.92156863,
         7.21649485,  7.57575758, 14.92537313, 11.53846154, 23.95833333,
        19.31818182, 29.41176471, 26.74418605, 25.49019608, 32.14285714,
        46.25      , 41.7721519 , 44.92753623, 44.21052632, 53.96825397,
        49.25373134, 59.49367089, 70.96774194, 59.67741935, 63.88888889,
        58.33333333, 62.16216216, 55.69620253, 58.69565217, 54.94505495,
        55.88235294, 50.        , 51.89873418, 53.125     , 54.41176471,
        36.23188406, 44.92753623, 28.04878049, 32.91139241, 22.72727273,
        20.96774194, 11.11111111,  9.63855422, 10.34482759, 12.72727273]])
Coordinates:
  * unit     (unit) int64 48B 1 2 3 4 5 6
  * feature  (feature) float64 960B 0.02618 0.07854 0.1309 ... 6.152 6.205 6.257
Attributes:
    occupancy:  [ 65.  57.  60.  65.  64.  85.  93.  91.  88.  80.  73.  99. ...
    bin_edges:  [array([0.        , 0.05235988, 0.10471976, 0.15707963, 0.209...
    fs:         100.0
    rates:      [12.28122812 13.16131613 13.23132313 14.97149715 14.51145115 ...

xarray.DataArray

unit: 6
feature: 120

66.15 52.63 66.67 60.0 62.5 68.24 ... 20.97 11.11 9.639 10.34 12.73

array([[66.15384615, 52.63157895, 66.66666667, 60.        , 62.5       ,
        68.23529412, 60.21505376, 52.74725275, 48.86363636, 48.75      ,
        34.24657534, 32.32323232, 27.18446602, 25.        , 20.45454545,
        26.73267327, 18.18181818, 10.46511628,  5.19480519,  3.61445783,
         6.59340659,  2.98507463,  7.57575758,  6.09756098,  1.13636364,
         1.51515152,  5.12820513,  1.2987013 ,  2.7027027 ,  1.33333333,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         1.25      ,  1.26582278,  1.44927536,  0.        ,  3.17460317,
         2.98507463,  2.53164557,  1.61290323,  3.22580645,  6.94444444,
...
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.80645161,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  1.0989011 ,
         0.93457944,  0.        ,  1.03092784,  5.10204082,  3.92156863,
         7.21649485,  7.57575758, 14.92537313, 11.53846154, 23.95833333,
        19.31818182, 29.41176471, 26.74418605, 25.49019608, 32.14285714,
        46.25      , 41.7721519 , 44.92753623, 44.21052632, 53.96825397,
        49.25373134, 59.49367089, 70.96774194, 59.67741935, 63.88888889,
        58.33333333, 62.16216216, 55.69620253, 58.69565217, 54.94505495,
        55.88235294, 50.        , 51.89873418, 53.125     , 54.41176471,
        36.23188406, 44.92753623, 28.04878049, 32.91139241, 22.72727273,
        20.96774194, 11.11111111,  9.63855422, 10.34482759, 12.72727273]])

Coordinates: (2)

unit
(unit)
int64
1 2 3 4 5 6
```
array([1, 2, 3, 4, 5, 6])
```

feature

(feature)

float64

0.02618 0.07854 ... 6.205 6.257

array([0.02618 , 0.07854 , 0.1309  , 0.18326 , 0.235619, 0.287979, 0.340339,
       0.392699, 0.445059, 0.497419, 0.549779, 0.602139, 0.654498, 0.706858,
       0.759218, 0.811578, 0.863938, 0.916298, 0.968658, 1.021018, 1.073377,
       1.125737, 1.178097, 1.230457, 1.282817, 1.335177, 1.387537, 1.439897,
       1.492257, 1.544616, 1.596976, 1.649336, 1.701696, 1.754056, 1.806416,
       1.858776, 1.911136, 1.963495, 2.015855, 2.068215, 2.120575, 2.172935,
       2.225295, 2.277655, 2.330015, 2.382374, 2.434734, 2.487094, 2.539454,
       2.591814, 2.644174, 2.696534, 2.748894, 2.801253, 2.853613, 2.905973,
       2.958333, 3.010693, 3.063053, 3.115413, 3.167773, 3.220132, 3.272492,
       3.324852, 3.377212, 3.429572, 3.481932, 3.534292, 3.586652, 3.639011,
       3.691371, 3.743731, 3.796091, 3.848451, 3.900811, 3.953171, 4.005531,
       4.057891, 4.11025 , 4.16261 , 4.21497 , 4.26733 , 4.31969 , 4.37205 ,
       4.42441 , 4.47677 , 4.529129, 4.581489, 4.633849, 4.686209, 4.738569,
       4.790929, 4.843289, 4.895649, 4.948008, 5.000368, 5.052728, 5.105088,
       5.157448, 5.209808, 5.262168, 5.314528, 5.366887, 5.419247, 5.471607,
       5.523967, 5.576327, 5.628687, 5.681047, 5.733407, 5.785766, 5.838126,
       5.890486, 5.942846, 5.995206, 6.047566, 6.099926, 6.152286, 6.204645,
       6.257005])

Attributes: (4)
occupancy :
[ 65. 57. 60. 65. 64. 85. 93. 91. 88. 80. 73. 99. 103. 80. 88. 101. 110. 86. 77. 83. 91. 67. 66. 82. 88. 66. 78. 77. 74. 75. 66. 85. 85. 85. 84. 78. 71. 75. 79. 78. 97. 79. 73. 79. 77. 94. 78. 78. 124. 90. 103. 89. 116. 86. 84. 101. 94. 119. 104. 114. 105. 92. 92. 79. 79. 89. 124. 93. 89. 94. 92. 120. 85. 107. 91. 107. 82. 97. 98. 102. 97. 66. 67. 78. 96. 88. 68. 86. 102. 84. 80. 79. 69. 95. 63. 67. 79. 62. 62. 72. 72. 74. 79. 92. 91. 68. 64. 79. 64. 68. 69. 69. 82. 79. 66. 62. 72. 83. 58. 55.]
bin_edges :
[array([0. , 0.05235988, 0.10471976, 0.15707963, 0.20943951, 0.26179939, 0.31415927, 0.36651914, 0.41887902, 0.4712389 , 0.52359878, 0.57595865, 0.62831853, 0.68067841, 0.73303829, 0.78539816, 0.83775804, 0.89011792, 0.9424778 , 0.99483767, 1.04719755, 1.09955743, 1.15191731, 1.20427718, 1.25663706, 1.30899694, 1.36135682, 1.41371669, 1.46607657, 1.51843645, 1.57079633, 1.6231562 , 1.67551608, 1.72787596, 1.78023584, 1.83259571, 1.88495559, 1.93731547, 1.98967535, 2.04203522, 2.0943951 , 2.14675498, 2.19911486, 2.25147474, 2.30383461, 2.35619449, 2.40855437, 2.46091425, 2.51327412, 2.565634 , 2.61799388, 2.67035376, 2.72271363, 2.77507351, 2.82743339, 2.87979327, 2.93215314, 2.98451302, 3.0368729 , 3.08923278, 3.14159265, 3.19395253, 3.24631241, 3.29867229, 3.35103216, 3.40339204, 3.45575192, 3.5081118 , 3.56047167, 3.61283155, 3.66519143, 3.71755131, 3.76991118, 3.82227106, 3.87463094, 3.92699082, 3.97935069, 4.03171057, 4.08407045, 4.13643033, 4.1887902 , 4.24115008, 4.29350996, 4.34586984, 4.39822972, 4.45058959, 4.50294947, 4.55530935, 4.60766923, 4.6600291 , 4.71238898, 4.76474886, 4.81710874, 4.86946861, 4.92182849, 4.97418837, 5.02654825, 5.07890812, 5.131268 , 5.18362788, 5.23598776, 5.28834763, 5.34070751, 5.39306739, 5.44542727, 5.49778714, 5.55014702, 5.6025069 , 5.65486678, 5.70722665, 5.75958653, 5.81194641, 5.86430629, 5.91666616, 5.96902604, 6.02138592, 6.0737458 , 6.12610567, 6.17846555, 6.23082543, 6.28318531])]
fs :
100.0
rates :
[12.28122812 13.16131613 13.23132313 14.97149715 14.51145115 11.62116212]

It has 2 dimensions, one for the neurons, and one for the feature bins.

These can be easily visualised:

Show plotting code

fig, ax = plt.subplots(constrained_layout=True, figsize=(6, 4))
tuning_curves.plot.line(ax=ax, x="feature")
plt.ylabel("firing rate [Hz]")
plt.xticks([0, 2*np.pi], ["0", "2π"])
legend=ax.get_legend()
legend.set_frame_on(False)
legend.set_bbox_to_anchor((1.0, 0.7))
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()

We can see that the mean activity of the neurons is clearly modulated by the feature, and that each neuron has a clear preference, matching our generation pattern. For more, take a look at the Pynapple user guide on tuning curves, as well as the real examples.

Decoding

If we know the mean activity of many neurons with respect to a variable, we should be able to predict the current value of that variable, given the current activity of the neurons. This is the premise of Bayesian decoding. The tuning curves of a population of neurons is used as a ‘template’ of what activity matches what feature value. Combined with Poisson statistics, Bayesian decoding provides a simple but principled way of neural decoding.

Using Pynapple, we can easily do Bayesian decoding, provided we have:

a set of tuning curves computed using compute_tuning_curves
the activity of the neurons matching the tuning curves
a (temporal) bin size, to bin activity in
a smoothing window, if we want to smooth the activity (this often helps for Bayesian decoding, by reducing noise)

decoded, proba_feature = nap.decode_bayes(
    tuning_curves=tuning_curves,
    data=tsgroup,
    epochs=epochs,
    bin_size=0.02,
    sliding_window_size=4,
)
decoded

Time (s)
-------------------  ---------
0.01                 0.0785398
0.03                 0.0785398
0.05                 0.1309
0.06999999999999999  0.1309
0.09                 0.0261799
0.11                 6.25701
0.13                 5.94285
...
9.87                 1.33518
9.89                 1.54462
9.91                 1.12574
9.93                 1.12574
9.95                 1.12574
9.97                 0.968658
9.99                 0.968658
dtype: float64, shape: (500,)

The output consists of a Tsd, containing the predicted feature values at every time step, and a TsdFrame containing the estimated probabilities of each feature value at every time step:

Show plotting code

import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(figsize=(6, 5), nrows=2, ncols=1, sharex=True)
ax1.plot(
    feature,
    label="true",
    linestyle="none",
    marker="o",
    markersize=4,
    alpha=0.5
)
ax1.plot(
    decoded,
    label="decoded",
    c="orange",
    linestyle="none",
    marker="o",
    markersize=3,
)
ax1.spines['right'].set_visible(False)
ax1.spines['top'].set_visible(False)
fig.legend(frameon=False, bbox_to_anchor=(1.1, 0.8))
ax1.set_ylabel("feature")
ax1.set_yticks([0, 2*np.pi], ["0", "2π"])
proba_feature = (proba_feature.values - np.min(proba_feature.values, axis=1, keepdims=True)) / np.ptp(proba_feature.values, axis=1,keepdims=True)
im = ax2.imshow(proba_feature.T, aspect="auto", origin="lower", cmap="viridis", extent=(0, 10.0, 0, 2*np.pi))
cbar_ax = fig.add_axes([0.93, 0.1, 0.015, 0.36])
fig.colorbar(im, cax=cbar_ax, label="probability")
ax2.set_xlabel("time (s)")
ax2.set_ylabel("feature")
ax2.set_yticks([0, 2*np.pi], ["0", "2π"])
plt.show()

For more, take a look at the Pynapple user guide on decoding, as well as the real examples.

6. Integration with other packages

Pynaviz

Pynaviz provides interactive, high-performance visualisations designed to work seamlessly with Pynapple time series and video data.
It allows synchronised exploration of neural signals and behavioral recordings.
It is build on top of pygfx, a modern GPU-based rendering engine.

To install it:

pip install pynaviz[qt]

NeMoS

NeMoS (Neural ModelS) is a statistical modeling framework optimized for systems neuroscience and powered by JAX.
It streamlines the process of defining and selecting models, through a collection of easy-to-use methods for feature design.

The core of nemos includes GPU-accelerated, well-tested implementations of standard statistical models, currently focusing on the Generalized Linear Model (GLM).
Check out this page for many examples of neural modelling using nemos and pynapple.

To install it:

pip install nemos

7. Data

Pynapple provides many ways to load your data.

NWB

The main input format Pynapple can read is Neurodata Without Borders (NWB).
Pynapple will parse all data in the file that can be represented as a Pynapple object and present a dictionary-like interface for selecting them.

Show plotting code

import os
import requests
nwb_path = 'A2929-200711.nwb'

if nwb_path not in os.listdir("."):
  r = requests.get(f"https://osf.io/fqht6/download", stream=True)
  block_size = 1024*1024
  with open(nwb_path, 'wb') as f:
    for data in r.iter_content(block_size):
        f.write(data)

import pynapple as nap
nwb_path = 'A2929-200711.nwb'
data = nap.load_file(nwb_path)
data

A2929-200711
┍━━━━━━━━━━━━━━━━━━━━━━━┯━━━━━━━━━━━━━┑
│ Keys                  │ Type        │
┝━━━━━━━━━━━━━━━━━━━━━━━┿━━━━━━━━━━━━━┥
│ units                 │ TsGroup     │
│ position_time_support │ IntervalSet │
│ epochs                │ IntervalSet │
│ z                     │ Tsd         │
│ y                     │ Tsd         │
│ x                     │ Tsd         │
│ rz                    │ Tsd         │
│ ry                    │ Tsd         │
│ rx                    │ Tsd         │
┕━━━━━━━━━━━━━━━━━━━━━━━┷━━━━━━━━━━━━━┙

Raw data

Raw data can be loaded with pynapple through the NEO library (also used by SpikeInterface). Many formats are supported, see here.

import pynapple as nap
data = nap.EphysReader("path/to/raw")

Saving as `.npz`

Pynapple objects can be saved (and loaded) as Numpy .npz files.

import pynapple as nap
import numpy as np
tsd = nap.Tsd(t=np.arange(10), d=np.arange(10))
tsd.save("my_tsd.npz")
loaded = nap.load_file("my_tsd.npz")
loaded

Time (s)
----------  --
0            0
1            1
2            2
3            3
4            4
5            5
6            6
7            7
8            8
9            9
dtype: int64, shape: (10,)

Streaming from DANDI

Data in the NWB format can be directly streamed from the DANDI archive into Pynapple.
Here is useful function to keep handy that does this for you:

import pynapple as nap
from pynwb import NWBHDF5IO
from dandi.dandiapi import DandiAPIClient
import fsspec
from fsspec.implementations.cached import CachingFileSystem
import h5py

def stream_dandiset(dandiset_id, filepath):
    with DandiAPIClient() as client:
        asset = client.get_dandiset(dandiset_id, "draft").get_asset_by_path(filepath)
        s3_url = asset.get_content_url(follow_redirects=1, strip_query=True)

    # first, create a virtual filesystem based on the http protocol
    fs = fsspec.filesystem("http")

    # create a cache to save downloaded data to disk (optional)
    fs = CachingFileSystem(
        fs=fs,
        cache_storage="nwb-cache",  # Local folder for the cache
    )

    # next, open the file
    file = h5py.File(fs.open(s3_url, "rb"))
    io = NWBHDF5IO(file=file, load_namespaces=True)

    return nap.NWBFile(io.read())

Your turn

Now it is your turn!

Start by streaming a single session from this dataset on the DANDI archive, it contains recordings from head direction cells in the postsubiculum. Use the streaming function we provided in the data section! LINK TO PREVIOUS CHAPTERS.

Then, try to compute head direction tuning curves, and use those to decode head direction from the population activity.

If you are confident, you can write your own code using the snippets provided in this chapter.
If you want a bit more help, follow along in this notebook.

Origins

Goal

1. Pynapple fundamentals

Time series (Ts)

Time series with data (Tsd)

Multiple time series

Time series with multiple data columns (TsdFrame)

Intervals (IntervalSet)

2. Interactions

Time supports (time_support)

Restricting (restrict)

Taking values (value_from)

3. Wrangling

Numpy

Slicing

4. Core functions

Binning

Averaging

Thresholding

Smoothing

5. Advanced Functions

Autocorrelograms

Crosscorrelograms

Signal processing

Bandpass filter

Power-spectral density

Perievent

Tuning curves

Decoding

6. Integration with other packages

Pynaviz

NeMoS

7. Data

NWB

Raw data

Saving as .npz

Streaming from DANDI

Your turn

Time series (`Ts`)

Time series with data (`Tsd`)

Time series with multiple data columns (`TsdFrame`)

Intervals (`IntervalSet`)

Time supports (`time_support`)

Restricting (`restrict`)

Taking values (`value_from`)

Saving as `.npz`