Decoding#

Hide code cell content
import pynapple as nap
import numpy as np
import pandas as pd
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)

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

Hint

Input to the bayesian decoding functions always include the tuning curves computed from nap.compute_1d_tuning_curves or nap.compute_2d_tuning_curves.

1-dimensional decoding#

Hide code cell content
from scipy.ndimage import gaussian_filter1d

# Fake Tuning curves
N = 6 # Number of neurons
bins = np.linspace(0, 2*np.pi, 61)
x = np.linspace(-np.pi, np.pi, len(bins)-1)
tmp = np.roll(np.exp(-(1.5*x)**2), (len(bins)-1)//2)
tc = np.array([np.roll(tmp, i*(len(bins)-1)//N) for i in range(N)]).T

tc_1d = pd.DataFrame(index=bins[0:-1], data=tc)

# Feature
T = 10000
dt = 0.01
timestep = np.arange(0, T)*dt
feature = nap.Tsd(
    t=timestep,
    d=gaussian_filter1d(np.cumsum(np.random.randn(T)*0.5), 20)%(2*np.pi)
    )
index = np.digitize(feature, bins)-1

# Spiking activity

count = np.random.poisson(tc[index])>0
tsgroup = nap.TsGroup({i:nap.Ts(timestep[count[:,i]]) for i in range(N)})
epoch = nap.IntervalSet(0, 10)

To decode, we need to compute tuning curves in 1D.

tcurves_1d = nap.compute_1d_tuning_curves(
    tsgroup, feature, nb_bins=61, minmax=(0, 2 * np.pi)
)

We can display the tuning curves of each neurons

Hide code cell source
plt.figure()
plt.plot(tcurves_1d)
plt.xlabel("Feature position")
plt.ylabel("Rate (Hz)")
plt.show()
../_images/2239e767922cd5c42b5dfb43e1671f43d75ea890aeadbe21a71c0f7d14715613.png

nap.decode_1d performs bayesian decoding:

decoded, proba_feature = nap.decode_1d(
    tuning_curves=tcurves_1d , # 2D tuning curves
    group=tsgroup, # Spiking activity
    ep=epoch, # Small epoch
    bin_size=0.06,  # How to bin the spike trains
    feature=feature, # Features
)

decoded is Tsd object containing the decoded feature value. proba_feature is a TsdFrame containing the probabilities of being in a particular feature bin over time.

Hide code cell source
plt.figure(figsize=(12, 6))
plt.subplot(211)
plt.plot(feature.restrict(epoch), label="True")
plt.plot(decoded, label="Decoded")
plt.legend()
plt.xlim(epoch[0,0], epoch[0,1])
plt.subplot(212)
plt.imshow(proba_feature.values.T, aspect="auto", origin="lower", cmap="viridis")
plt.xticks([0, len(decoded)], epoch.values[0])
plt.xlabel("Time (s)")
plt.show()
../_images/c9a58367b9ce4c82726738e804e6fe76571e681afc97b596cd49f4231d77095e.png

2-dimensional decoding#

Hide code cell content
dt = 0.1
epoch = nap.IntervalSet(start=0, end=1000, time_units="s")
features = np.vstack((np.cos(np.arange(0, 1000, dt)), np.sin(np.arange(0, 1000, dt)))).T
features = nap.TsdFrame(t=np.arange(0, 1000, dt),
    d=features,
    time_units="s",
    time_support=epoch,
    columns=["a", "b"],
)

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)

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

tcurves2d, binsxy = nap.compute_2d_tuning_curves(
    group=ts_group, # Spiking activity of 12 neurons
    features=features, # 2-dimensional features
    nb_bins=10,
    ep=epoch,
    minmax=(-1.0, 1.0, -1.0, 1.0), # Minmax of the features
)
/home/runner/.local/lib/python3.10/site-packages/pynapple/process/tuning_curves.py:269: RuntimeWarning: invalid value encountered in divide
  count = count / occupancy

We can display the tuning curves of each neuron

Hide code cell source
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()
../_images/e7475d13bf95bec464045f2e30be0ee4221eabe5a26104793f3fbe1a2891d4ff.png

nap.decode_2d performs bayesian decoding:

decoded, proba_feature = nap.decode_2d(
    tuning_curves=tcurves2d, # 2D tuning curves
    group=ts_group, # Spiking activity
    ep=epoch, # Epoch
    bin_size=0.1,  # How to bin the spike trains
    xy=binsxy, # Features position
    features=features, # Features
)
Hide code cell source
plt.figure(figsize=(15, 5))
plt.subplot(131)
plt.plot(features["a"].get(0,20), label="True")
plt.plot(decoded["a"].get(0,20), label="Decoded")
plt.legend()
plt.xlabel("Time (s)")
plt.ylabel("Feature a")
plt.subplot(132)
plt.plot(features["b"].get(0,20), label="True")
plt.plot(decoded["b"].get(0,20), label="Decoded")
plt.legend()
plt.xlabel("Time (s)")
plt.title("Feature b")
plt.subplot(133)
plt.plot(
    features["a"].get(0,20),
    features["b"].get(0,20),
    label="True",
)
plt.plot(
    decoded["a"].get(0,20),
    decoded["b"].get(0,20),
    label="Decoded",
)
plt.xlabel("Feature a")
plt.title("Feature b")
plt.legend()
plt.tight_layout()
plt.show()
../_images/92ca446e6d077c80e2ec54e93ce1cb92ef0c88bf9429dedb009729b38b4d8227.png