Decoding#
Show 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#
Show 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
Show code cell source
plt.figure()
plt.plot(tcurves_1d)
plt.xlabel("Feature position")
plt.ylabel("Rate (Hz)")
plt.show()
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.
Show 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()
2-dimensional decoding#
Show 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
Show 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()
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
)
Show 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()
