Perievent / Spike-triggered average

The perievent module allows for aligning timeseries and timestamps data around events, as well as computing event-triggered averages (e.g. spike-triggered averages).

Show code cell content

Hide code cell content

import pynapple as nap
import numpy as np
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)

Peri-Event Time Histograms (PETH)

Spike data

We will start with the common use-case of aligning the spiking activity of a unit to a set of events.

Single unit

Let’s simulate some uniform stimuli and a unit that has a gaussian firing field after the stimulus:

Show code cell source

Hide code cell source

stimuli = nap.Ts(t=np.arange(0, 1000, 1), time_units="s")

def generate_spiking_unit(offset):
    baseline = np.random.uniform(0, 1000, 500)
    burst = np.concatenate([
        np.random.normal(st + offset, 0.05, 3) for st in stimuli.times()
    ])
    return nap.Ts(t=np.sort(np.concatenate([baseline, burst])), time_units="s")

ts = generate_spiking_unit(offset=0.1)

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.set_xlabel("time (s)")
ax.legend(loc="upper left", bbox_to_anchor=(1, 1));

The compute_perievent function can align timestamps to a particular set of events:

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

Index    rate    events
-------  ------  --------
0        6.0     1.0
1        10.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      6.0     993.0
993      6.0     994.0
994      10.0    995.0
995      8.0     996.0
996      8.0     997.0
997      6.0     998.0
998      6.0     999.0

The returned object is a TsGroup containing the aligned timestamps per event. The event times are stored in the events metadata column.

If you want to get aligned counts, you can now easily call:

bin_size = 0.01
peth_counts = peth.count(bin_size)
peth_counts

Time (s)    0    1    2    3    4    ...
----------  ---  ---  ---  ---  ---  -----
-0.095      0    0    0    0    0    ...
-0.085      0    0    0    0    0    ...
-0.075      0    0    0    0    0    ...
-0.065      0    0    0    0    0    ...
-0.055      0    0    0    0    0    ...
-0.045      0    0    0    0    0    ...
-0.035      0    0    0    0    0    ...
...                                  ...
0.335       0    0    0    0    0    ...
0.345       0    0    0    0    0    ...
0.355       0    0    0    0    0    ...
0.365       0    0    0    0    0    ...
0.375       0    0    0    0    0    ...
0.385       0    0    0    0    0    ...
0.395       0    0    0    0    0    ...
Metadata
events      1.0  2.0  3.0  4.0  5.0  ...
dtype: int64, shape: (50, 999)

Visualizing the aligned timestamps can now be done by using the to_tsd function. This function flattens the TsGroup into one Tsd containing all the timestamps and storing the event ids as data. We can also take the mean of the counts and divide by the bin size to show an estimate of the aligned firing rate:

Show code cell source

Hide code cell source

def plot_peth(unit_peth, unit_peth_counts, ax_mean, ax_spikes, color=None):
    mean = np.mean(unit_peth_counts / bin_size, axis=1)
    ax_mean.plot(mean, color=color)
    ax_mean.set_ylabel("spikes/s")
    ax_mean.axvline(0.0, color="gray", linestyle="--")
    ax_spikes.plot(unit_peth.to_tsd(), "|", markersize=5, color=color)
    ax_mean.set_xlabel("time from event (s)")
    ax_spikes.set_ylabel("event")
    ax_spikes.axvline(0.0, color="gray", linestyle="--")
    

fig, (ax_spikes, ax_mean) = plt.subplots(
    2, 1, sharex=True, height_ratios=[1, 0.3], figsize=(6, 6), gridspec_kw={"hspace": 0.3}
)
plot_peth(peth, peth_counts, ax_mean, ax_spikes)
ax_spikes.set_title("peri = nap.compute_perievent(spikes, events)")
ax_mean.set_title("peri.count(bin_size) / bin_size")

Text(0.5, 1.0, 'peri.count(bin_size) / bin_size')

Multiple units

The same function can be applied to a group of units:

Show code cell source

Hide code cell source

tsgroup = nap.TsGroup(
    {1: ts, 2: generate_spiking_unit(offset=0.2), 3: generate_spiking_unit(0.3)}
)

fig, ax = plt.subplots(1, 1, constrained_layout=True)
unit_spacing = 0.15
for i, unit in enumerate(tsgroup):
    ax.vlines(
        tsgroup[unit].restrict(segment).times(),
        i * unit_spacing + 0.02,
        i * unit_spacing + 0.12,
        label=f"unit {unit}",
        color=plt.cm.tab10(i)
    )
ax.vlines(
    stimuli.restrict(segment).times(),
    0.0,
    0.44,
    color="gray",
    linestyle="--",
    label="stimulus",
)
ax.yaxis.set_visible(False)
ax.spines["left"].set_visible(False)
ax.set_xlabel("time (s)")
ax.legend(loc="upper left", bbox_to_anchor=(1, 1));

In this case, it returns a dict of TsGroup, containing the same object as before, but now per unit.

peth = nap.compute_perievent(data=tsgroup, events=stimuli, window=window)
peth

{1: Index    rate    events
 -------  ------  --------
 0        6.0     1.0
 1        10.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      6.0     993.0
 993      6.0     994.0
 994      10.0    995.0
 995      8.0     996.0
 996      8.0     997.0
 997      6.0     998.0
 998      6.0     999.0,
 2: Index    rate    events
 -------  ------  --------
 0        6.0     1.0
 1        8.0     2.0
 2        6.0     3.0
 3        6.0     4.0
 4        8.0     5.0
 5        6.0     6.0
 6        6.0     7.0
 ...      ...     ...
 992      6.0     993.0
 993      6.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,
 3: 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      4.0     993.0
 993      6.0     994.0
 994      8.0     995.0
 995      6.0     996.0
 996      8.0     997.0
 997      6.0     998.0
 998      6.0     999.0}

We can again visualize easily:

Show code cell source

Hide code cell source

fig, axs = plt.subplots(
    2,
    len(tsgroup),
    sharey="row",
    sharex=True,
    height_ratios=[0.3, 1.0],
    figsize=(15, 8),
)
fig.suptitle("nap.compute_perievent(tsgroup, stimuli, (-0.1, 0.4)", y=1.02)
for i, (unit, unit_axs) in enumerate(zip(tsgroup, axs.T)):
    plot_peth(peth[unit], peth[unit].count(bin_size), *unit_axs, color=plt.cm.tab10(i))
    unit_axs[0].set_title(f"unit {unit}")

Continuous data

If you have continuous data, e.g. calcium imaging traces, you can use the exact same function! compute_perievent is designed in such a way that it recognizes the input format and decides the correct computation. Hence, when given a Tsd, TsdFrame or even a TsdTensor, it will know what to do.

Note

compute_perievent only works with regularly sampled continuous data. If your data is irregularly sampled, try padding it with nan first.

Single unit

Let’s again start by simulating some data, but this time continuous traces:

Show code cell source

Hide code cell source

def generate_continuous_unit(burst_offset):
    t = np.arange(0, 1000, 0.02)
    d = np.random.uniform(0, 1, len(t))
    for st in stimuli.times():
        d += 4 * np.exp(-((t - (st + burst_offset)) ** 2) / (2 * 0.05**2))
    return nap.Tsd(t=t, d=np.clip(d, 0, 5), time_units="s")

tsd = generate_continuous_unit(burst_offset=0.1)

segment = nap.IntervalSet(100, 102.9)
fig, ax = plt.subplots(1, 1, constrained_layout=True)
ax.plot(tsd.restrict(segment), label="activity")
ax.vlines(
    stimuli.restrict(segment).times(),
    0.0,
    tsd.max(),
    color="gray",
    linestyle="--",
    label="stimulus",
)
ax.yaxis.set_visible(False)
ax.spines["left"].set_visible(False)
ax.set_xlabel("time (s)")
ax.legend(loc="upper left", bbox_to_anchor=(1, 1));

We can pass continuous units (as a Tsd) to the function in the exact same way:

peth = nap.compute_perievent(data=tsd, events=stimuli, window=window)
peth

Time (s)             0         1          2          3          4  ...
----------  ----------  --------  ---------  ---------  ---------  -----
-0.1        nan         0.67001   0.0352572  0.427523   0.108556   ...
-0.08       nan         0.235106  0.951892   0.239927   0.165249   ...
-0.06       nan         0.642842  0.918455   0.211255   0.783425   ...
-0.04       nan         0.789414  0.978445   0.386464   0.977083   ...
-0.02       nan         0.360932  0.748351   0.476296   1.1472     ...
0.0           1.08306   1.50333   0.678878   1.19077    0.99051    ...
0.02          2.00771   1.79345   1.71931    1.75735    1.48297    ...
...                                                                ...
0.28          0.600366  0.955442  0.673684   0.798336   0.554792   ...
0.3           0.819278  0.963498  0.655775   0.497935   0.914406   ...
0.32          0.202198  0.32481   0.433978   0.730133   0.102648   ...
0.34          0.147349  0.338948  0.37255    0.0411339  0.0875405  ...
0.36          0.785896  0.127397  0.386794   0.846991   0.646502   ...
0.38          0.6338    0.169052  0.894664   0.171691   0.257385   ...
0.4           0.551991  0.73567   0.575797   0.260073   0.93598    ...
dtype: float64, shape: (26, 1000)

This time, it will return a TsdFrame with a column per event.

We can again visualize, this time using a heatmap (i.e. using imshow):

Show code cell source

Hide code cell source

def plot_peth_continuous(unit_peth, ax_mean, ax, color=None):
    mean = np.nanmean(unit_peth, axis=1)
    ax_mean.plot(mean, color=color)
    ax_mean.set_ylabel("dF/F [a.u.]")
    ax_mean.axvline(0.0, color="gray", linestyle="--")
    im = ax.imshow(
        unit_peth.values.T,
        extent=(unit_peth.times()[0], unit_peth.times()[-1], 0, unit_peth.shape[1]),
        interpolation="none",
        aspect="auto",
        cmap="Grays",
    )
    ax.axvline(0.0, color="gray", linestyle="--")
    ax.set_xlabel("time from event (s)")
    ax.set_ylabel("event")
    return im

fig, (ax_mean, ax) = plt.subplots(
    2, 1, sharex=True, height_ratios=[0.3, 1.0], figsize=(6, 7)
)
fig.suptitle("nap.compute_perievent(tsd, stimuli, (-0.1, 0.4)", y=1.02)
im = plot_peth_continuous(peth, ax_mean, ax)
fig.colorbar(im, cax=fig.add_axes([0.92, 0.14, 0.02, 0.3]), label="dF/F [a.u.]");

Multiple units

The same function can also handle multiple continuous units, passed as a TsdFrame:

Show code cell source

Hide code cell source

tsdframe = np.stack(
    [tsd, generate_continuous_unit(0.2), generate_continuous_unit(0.3)], axis=1
)
fig, ax = plt.subplots(1, 1, constrained_layout=True)
for i in range(tsdframe.shape[1]):
    ax.plot(
        tsdframe[:, i].restrict(segment), color=plt.cm.tab10(i), label=f"unit {i+1}"
    )
ax.vlines(
    stimuli.restrict(segment).times(),
    0.0,
    tsd.max(),
    color="gray",
    linestyle="--",
    label="stimulus",
)
ax.yaxis.set_visible(False)
ax.spines["left"].set_visible(False)
ax.set_xlabel("time (s)")
ax.legend(loc="upper left", bbox_to_anchor=(1, 1));

peth = nap.compute_perievent(data=tsdframe, events=stimuli, window=window)
peth

Time (s)
----------  -----------------------------
-0.1        [[nan ... nan] ...]
-0.08       [[nan ... nan] ...]
-0.06       [[nan ... nan] ...]
-0.04       [[nan ... nan] ...]
-0.02       [[nan ... nan] ...]
0.0         [[1.083058 ... 0.885986] ...]
0.02        [[2.007713 ... 0.378317] ...]
...
0.28        [[0.600366 ... 4.078696] ...]
0.3         [[0.819278 ... 4.204869] ...]
0.32        [[0.202198 ... 4.117117] ...]
0.34        [[0.147349 ... 2.961957] ...]
0.36        [[0.785896 ... 2.018256] ...]
0.38        [[0.6338   ... 1.732266] ...]
0.4         [[0.551991 ... 0.686566] ...]
dtype: float64, shape: (26, 1000, 3)

In this case, it returns a TsdTensor, containing an added dimension for the events. We can again visualize by reusing our plotting function:

Show code cell source

Hide code cell source

fig, axs = plt.subplots(
    2, len(tsdframe.columns), sharex=True, height_ratios=[0.3, 1.0], figsize=(15, 8)
)

for i, unit_axs in zip(range(tsdframe.shape[1]), axs.T):
    im = plot_peth_continuous(peth[:, :, i], *unit_axs, color=plt.cm.tab10(i))
    unit_axs[0].set_title(f"unit {i+1}")
fig.suptitle("nap.compute_perievent(tsdframe, stimuli, (-0.1, 0.4)", y=1.02)
fig.colorbar(im, cax=fig.add_axes([0.92, 0.14, 0.02, 0.3]), label="dF/F [a.u.]");

In the most complex case, you can even pass a TsdTensor to compute_perievent. It will return a TsdTensor with an added dimension for the events. Visualizing such a tensor becomes a bit complex, but feel free to try!

Event-Triggered Average (ETA/STA)

The compute_event_triggered_average computes the average of a continuous feature aligned to a set of events. The classic use-case is to recover the stimulus feature that drives a neuron’s spiking.

Note

In neuroscience, this is commonly known as a spike-triggered average (STA). For convenience, compute_spike_triggered_average is provided as an alias.

Single unit

Let’s simulate a slowly varying stimulus and a neuron that fires preferentially at its peaks:

Show code cell source

Hide code cell source

t = np.arange(0, 1000, 0.02)
feature = nap.Tsd(t=t, d=np.sin(2 * np.pi * t / 10), time_units="s")

def generate_spiking_unit(phase):
    rate = np.clip(np.sin(2 * np.pi * t / 10 + phase) * 10 + 5, 0, None)
    return nap.Ts(t=np.sort(t[np.random.rand(len(t)) < rate * 0.01]), time_units="s")

ts = generate_spiking_unit(phase=0.0)

segment = nap.IntervalSet(100, 124.9)
fig, ax = plt.subplots(1, 1, constrained_layout=True)
ax.vlines(ts.restrict(segment).times(), 1.1, 1.5, label="spikes")
ax.plot(feature.restrict(segment), color="black", label="feature")
ax.yaxis.set_visible(False)
ax.spines["left"].set_visible(False)
ax.set_xlabel("time (s)")
ax.legend(loc='upper left', bbox_to_anchor=(1, 1));

We can pass this to the function, passing a bin size and a window:

eta = nap.compute_event_triggered_average(feature, ts, binsize=0.02, window=(-5,5))
eta

Time (s)            0
----------  ---------
-5.0        -0.657077
-4.98       -0.656894
-4.96       -0.656608
-4.94       -0.656218
-4.92       -0.655725
-4.9        -0.655127
-4.88       -0.654427
...
4.88        -0.664481
4.9         -0.664927
4.92        -0.665268
4.94        -0.665503
4.96        -0.665634
4.98        -0.665659
5.0         -0.665579
dtype: float64, shape: (501, 1)

The result is a TsdFrame with one column. If the neuron is driven by the stimulus, the ETA should recover the stimulus waveform preceding each spike:

Show code cell source

Hide code cell source

plt.plot(eta)
plt.axvline(0.0, color="gray", linestyle="--")
plt.xlabel("time from spike (s)")
plt.ylabel("stimulus [a.u.]")
plt.title("Spike triggered average \n nap.compute_event_triggered_average(feature, ts, binsize=0.02, window=(-5,5))");
plt.tight_layout();

Multiple units

The same function can be used for a group of units. When passing a TsGroup, the function returns one column per unit. Here, we simulate units driven by the stimulus at different phases:

Show code cell source

Hide code cell source

# simulation
tsgroup = nap.TsGroup({
    1: ts,
    2: generate_spiking_unit(phase=np.pi / 2),
    3: generate_spiking_unit(phase=np.pi),
})

# visualization
fig, ax = plt.subplots(1, 1, constrained_layout=True)
unit_spacing = 0.4
y_positions = np.linspace(1.24, 2.2, len(tsgroup))
for i, unit in enumerate(tsgroup):
    ax.vlines(
        tsgroup[unit].restrict(segment).times(),
        y_positions[i] - unit_spacing / 2,
        y_positions[i] + unit_spacing / 2,
        label=f"unit {unit}",
        color=plt.cm.tab10(i)
    )
ax.plot(feature.restrict(segment), color="black", label="feature")
ax.yaxis.set_visible(False)
ax.spines["left"].set_visible(False)
ax.set_xlabel("time (s)")
ax.legend(loc='upper left', bbox_to_anchor=(1, 1));

eta = nap.compute_event_triggered_average(feature, tsgroup, binsize=0.02, window=10.0)
eta

Time (s)           1             2          3
----------  --------  ------------  ---------
-10.0       0.657175  -0.000487287  -0.649027
-9.98       0.656984   0.00773503   -0.649103
-9.96       0.65669    0.0159561    -0.649076
-9.94       0.656292   0.0241747    -0.648943
-9.92       0.65579    0.0323937    -0.648707
-9.9        0.655185   0.0406075    -0.648369
-9.88       0.654476   0.0488192    -0.647929
...
9.88        0.65782   -0.0517478    -0.647786
9.9         0.65824   -0.0435362    -0.648461
9.92        0.658557  -0.035326     -0.649035
9.94        0.65877   -0.0271103    -0.649506
9.96        0.658879  -0.0188903    -0.649874
9.98        0.658883  -0.0106674    -0.65014
10.0        0.658784  -0.00243851   -0.650303
dtype: float64, shape: (1001, 3)

Each unit recovers a phase-shifted version of the stimulus:

Show code cell source

Hide code cell source

fig, ax = plt.subplots()
for i, unit in enumerate(eta.columns):
    ax.plot(eta[:, i], label=f"unit {unit}", color=plt.cm.tab10(i))
ax.axvline(0.0, color="gray", linestyle="--")
ax.set_xlabel("time from spike (s)")
ax.set_ylabel("stimulus [a.u.]")
ax.set_title("Spike Triggered Average")
ax.legend(loc='upper left', bbox_to_anchor=(1, 1));

Multiple features

You can also pass multiple features at once. Here, we simulate two units each driven by a different frequency, and pass both stimulus features together:

Show code cell source

Hide code cell source

tsdframe = nap.TsdFrame(
    t=t,
    d=np.stack([
        np.sin(2 * np.pi * t / 10),
        np.sin(2 * np.pi * t / 5),
    ], axis=1),
    columns=["10s", "5s"],
)

def generate_spiking_unit_from_feature(feature):
    rate = np.clip(feature * 10 + 10, 0, None)
    return nap.Ts(t=np.sort(t[np.random.rand(len(t)) < rate * 0.02]), time_units="s")

tsgroup = nap.TsGroup({
    1: generate_spiking_unit_from_feature(tsdframe[:, 0].d),
    2: generate_spiking_unit_from_feature(tsdframe[:, 1].d),
})

eta = nap.compute_event_triggered_average(tsdframe, tsgroup, binsize=0.02, window=10.0)
eta

Time (s)
----------  ----------------------------
-10.0       [[ 0.501369, -0.008342] ...]
-9.98       [[ 0.501225, -0.008271] ...]
-9.96       [[ 0.501003, -0.008195] ...]
-9.94       [[ 0.500702, -0.008112] ...]
-9.92       [[ 0.500323, -0.008023] ...]
-9.9        [[ 0.499864, -0.007929] ...]
-9.88       [[ 0.499327, -0.00783 ] ...]
...
9.88        [[ 0.501959, -0.009343] ...]
9.9         [[ 0.502293, -0.009328] ...]
9.92        [[ 0.502547, -0.009308] ...]
9.94        [[ 0.502721, -0.009281] ...]
9.96        [[ 0.502818, -0.009246] ...]
9.98        [[ 0.502833, -0.009211] ...]
10.0        [[ 0.502768, -0.009169] ...]
dtype: float64, shape: (1001, 2, 2)

The result is a TsdTensor. Each unit recovers its own driving feature while the other averages to near zero:

Show code cell source

Hide code cell source

fig, axs = plt.subplots(1, len(tsgroup), sharey=True, figsize=(10, 4))
for i, (unit, ax) in enumerate(zip(tsgroup, axs)):
    for feat in range(len(tsdframe.columns)):
        ax.plot(eta.t, eta[:, i, feat], label=tsdframe.columns[feat])
    ax.axvline(0.0, color="gray", linestyle="--")
    ax.set_title(f"unit {unit}")
    ax.set_xlabel("time from spike (s)")
axs[0].set_ylabel("stimulus [a.u.]")
fig.suptitle("Spike Triggered Average", y=1.1)
axs[-1].legend(loc='upper left', bbox_to_anchor=(1, 1));