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). It can also compute average firing rate for different epochs (for example firing rate for different epochs of stimulus presentation).
Important
If you are using calcium imaging data with the activity of the cell as a continuous transient, the function to call ends with _continuous for continuous time series (e.g. compute_1d_tuning_curves_continuous).
Show code cell content
import pynapple as nap
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from pprint import pprint
custom_params = {"axes.spines.right": False, "axes.spines.top": False}
sns.set_theme(style="ticks", palette="colorblind", font_scale=1.5, rc=custom_params)
Show code cell content
group = {
0: nap.Ts(t=np.sort(np.random.uniform(0, 100, 10))),
1: nap.Ts(t=np.sort(np.random.uniform(0, 100, 20))),
2: nap.Ts(t=np.sort(np.random.uniform(0, 100, 30))),
}
tsgroup = nap.TsGroup(group)
from epochs#
The epochs should be stored in a dictionnary :
dict_ep = {
"stim0": nap.IntervalSet(start=0, end=20),
"stim1":nap.IntervalSet(start=30, end=70)
}
nap.compute_discrete_tuning_curves takes a TsGroup for spiking activity and a dictionary of epochs.
The output is a pandas DataFrame where each column is a unit in the TsGroup and each row is one IntervalSet type.
The value is the mean firing rate of the neuron during this set of intervals.
mean_fr = nap.compute_discrete_tuning_curves(tsgroup, dict_ep)
pprint(mean_fr)
0 1 2
stim0 0.200 0.150 0.25
stim1 0.075 0.275 0.30
from timestamps activity#
1-dimension tuning curves#
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)
generative_tc = np.array([np.roll(tmp, i*(len(bins)-1)//N) for i in range(N)]).T
# Feature
T = 50000
dt = 0.002
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(generative_tc[index])>0
tsgroup = nap.TsGroup(
{i:nap.Ts(timestep[count[:,i]]) for i in range(N)},
time_support = nap.IntervalSet(0, 100)
)
Mandatory arguments are TsGroup, Tsd (or TsdFrame with 1 column only)
and nb_bins for number of bins of the tuning curves.
If an IntervalSet is passed with ep, everything is restricted to ep
otherwise the time support of the feature is used.
The min and max of the tuning curve is by default the min and max of the feature. This can be tweaked with the argument minmax.
The output is a pandas DataFrame. Each column is a unit from TsGroup argument. The index of the DataFrame carries the center of the bin in feature space.
tuning_curve = nap.compute_1d_tuning_curves(
group=tsgroup,
feature=feature,
nb_bins=120,
minmax=(0, 2*np.pi)
)
print(tuning_curve)
0 1 2 3 4 5
0.026180 323.667188 45.759844 0.0 0.0 0.000000 34.598906
0.078540 322.154094 46.980805 0.0 0.0 0.000000 31.320537
0.130900 296.168993 65.948561 0.0 0.0 0.000000 15.587842
0.183260 314.776271 65.376610 0.0 0.0 0.000000 12.106780
0.235619 277.431709 100.212975 0.0 0.0 0.000000 7.384114
... ... ... ... ... ... ...
6.047566 280.411014 4.504595 0.0 0.0 2.252297 100.227230
6.099926 305.561667 24.573141 0.0 0.0 0.000000 72.651026
6.152286 304.835071 14.084789 0.0 0.0 0.000000 65.393662
6.204645 311.230714 29.337321 0.0 0.0 0.000000 39.541607
6.257005 305.699183 22.277673 0.0 0.0 0.000000 51.981238
[120 rows x 6 columns]
plt.figure()
plt.plot(tuning_curve)
plt.xlabel("Feature space")
plt.ylabel("Firing rate (Hz)")
plt.show()
Internally, the function is calling the method value_from which maps a timestamps
to its closest value in time from a Tsd object.
It is then possible to validate the tuning curves by displaying the
timestamps as well as their associated values.
Show code cell source
plt.figure()
plt.subplot(121)
plt.plot(tsgroup[3].value_from(feature), 'o')
plt.plot(feature, label="feature")
plt.ylabel("Feature")
plt.xlim(0, 2)
plt.xlabel("Time (s)")
plt.subplot(122)
plt.plot(tuning_curve[3].values, tuning_curve[3].index.values, label="Tuning curve (unit=3)")
plt.xlabel("Firing rate (Hz)")
plt.legend()
plt.show()
2-dimension tuning curves#
Show code cell content
dt = 0.01
T = 10
epoch = nap.IntervalSet(start=0, end=T, time_units="s")
features = np.vstack((np.cos(np.arange(0, T, dt)), np.sin(np.arange(0, T, dt)))).T
features = nap.TsdFrame(
t=np.arange(0, T, dt),
d=features,
time_units="s",
time_support=epoch,
columns=["a", "b"],
)
tsgroup = nap.TsGroup({
0: nap.Ts(t=np.sort(np.random.uniform(0, T, 10))),
1: nap.Ts(t=np.sort(np.random.uniform(0, T, 15))),
2: nap.Ts(t=np.sort(np.random.uniform(0, T, 20))),
}, time_support=epoch)
The group argument must be a TsGroup object.
The features argument must be a 2-columns TsdFrame object.
nb_bins can be an int or a tuple of 2 ints.
tcurves2d, binsxy = nap.compute_2d_tuning_curves(
group=tsgroup,
features=features,
nb_bins=(5,5),
minmax=(-1, 1, -1, 1)
)
pprint(tcurves2d)
{0: array([[7.14285714, 0. , 1.25 , 2.24719101, 0. ],
[0. , nan, nan, nan, 1.13636364],
[0. , nan, nan, nan, 2.4691358 ],
[0. , nan, nan, nan, 1.13636364],
[0. , 0. , 1.63934426, 0. , 0. ]]),
1: array([[3.57142857, 0. , 1.25 , 1.12359551, 3.50877193],
[0. , nan, nan, nan, 0. ],
[2.5 , nan, nan, nan, 2.4691358 ],
[0. , nan, nan, nan, 1.13636364],
[0. , 2.22222222, 1.63934426, 3.40909091, 1.75438596]]),
2: array([[3.57142857, 3.7037037 , 0. , 4.49438202, 1.75438596],
[0. , nan, nan, nan, 1.13636364],
[0. , nan, nan, nan, 2.4691358 ],
[2.27272727, nan, nan, nan, 3.40909091],
[3.57142857, 0. , 0. , 2.27272727, 1.75438596]])}
/home/runner/.local/lib/python3.12/site-packages/pynapple/process/tuning_curves.py:269: RuntimeWarning: invalid value encountered in divide
count = count / occupancy
tcurves2d is a dictionnary with each key a unit in TsGroup. binsxy is a numpy array representing the centers of the bins and is useful for plotting tuning curves. Bins that have never been visited by the feature have been assigned a NaN value.
Checking the accuracy of the tuning curves can be bone by displaying the spikes aligned to the features with the function value_from which assign to each spikes the corresponding features value for unit 0.
ts_to_features = tsgroup[0].value_from(features)
print(ts_to_features)
Time (s) a b
---------- -------- --------
1.55063 0.02079 0.99978
2.58047 -0.84641 0.53253
3.83257 -0.77226 -0.63531
3.95176 -0.69065 -0.72319
6.2142 0.99732 -0.07312
7.51446 0.33724 0.94142
7.92761 -0.07595 0.99711
8.24842 -0.38575 0.9226
8.89764 -0.86544 0.50102
9.61445 -0.9829 -0.18416
dtype: float64, shape: (10, 2)
tsgroup[0] which is a Ts object has been transformed to a TsdFrame object with each timestamps (spike times) being associated with a features value.
Show code cell source
plt.figure()
plt.subplot(121)
plt.plot(features["b"], features["a"], label="features")
plt.plot(ts_to_features["b"], ts_to_features["a"], "o", color="red", markersize=4, label="spikes")
plt.xlabel("feature b")
plt.ylabel("feature a")
[plt.axvline(b, linewidth=0.5, color='grey') for b in np.linspace(-1, 1, 6)]
[plt.axhline(b, linewidth=0.5, color='grey') for b in np.linspace(-1, 1, 6)]
plt.subplot(122)
extents = (
np.min(features["a"]),
np.max(features["a"]),
np.min(features["b"]),
np.max(features["b"]),
)
plt.imshow(tcurves2d[0],
origin="lower", extent=extents, cmap="viridis",
aspect='auto'
)
plt.title("Tuning curve unit 0")
plt.xlabel("feature b")
plt.ylabel("feature a")
plt.grid(False)
plt.colorbar()
plt.tight_layout()
plt.show()
from continuous activity#
Tuning curves compute with the following functions are usually made with data from calcium imaging activities.
1-dimension tuning curves#
Show code cell content
from scipy.ndimage import gaussian_filter1d
# Fake Tuning curves
N = 3 # 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)
generative_tc = np.array([np.roll(tmp, i*(len(bins)-1)//N) for i in range(N)]).T
# Feature
T = 50000
dt = 0.002
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
tmp = generative_tc[index]
tmp = tmp + np.random.randn(*tmp.shape)*1
# Calcium activity
tsdframe = nap.TsdFrame(
t=timestep,
d=tmp
)
Arguments are TsdFrame (for example continuous calcium data), Tsd or TsdFrame for the 1-d feature and nb_bins for the number of bins.
tuning_curves = nap.compute_1d_tuning_curves_continuous(
tsdframe=tsdframe,
feature=feature,
nb_bins=120,
minmax=(0, 2*np.pi)
)
print(tuning_curves)
0 1 2
0.026180 0.964677 0.038873 -0.067743
0.078540 0.948832 -0.004711 0.036663
0.130900 0.917935 0.097936 -0.013865
0.183260 0.993311 0.021283 0.028429
0.235619 0.843370 -0.042239 -0.059548
... ... ... ...
6.047566 0.814062 -0.068988 0.096578
6.099926 0.953426 -0.019140 0.012197
6.152286 0.958816 0.036384 -0.056039
6.204645 1.051622 -0.040729 -0.034881
6.257005 0.938252 -0.047950 0.048416
[120 rows x 3 columns]
plt.figure()
plt.plot(tuning_curves)
plt.xlabel("Feature space")
plt.ylabel("Mean activity")
plt.show()
2-dimension tuning curves#
Show code cell content
dt = 0.01
T = 10
epoch = nap.IntervalSet(start=0, end=T, time_units="s")
features = np.vstack((np.cos(np.arange(0, T, dt)), np.sin(np.arange(0, T, dt)))).T
features = nap.TsdFrame(
t=np.arange(0, T, dt),
d=features,
time_units="s",
time_support=epoch,
columns=["a", "b"],
)
# Calcium activity
tsdframe = nap.TsdFrame(
t=timestep,
d=np.random.randn(len(timestep), 2)
)
Arguments are TsdFrame (for example continuous calcium data), Tsd or TsdFrame for the 1-d feature and nb_bins for the number of bins.
tuning_curves, xy = nap.compute_2d_tuning_curves_continuous(
tsdframe=tsdframe,
features=features,
nb_bins=5,
)
print(tuning_curves)
{0: array([[ 0.0595357 , -0.02103981, 0.07091826, -0.03199865, 0.05891291],
[ 0.12010567, nan, nan, nan, -0.04951179],
[-0.05571585, nan, nan, nan, 0.02446444],
[ 0.01925542, nan, nan, nan, 0.00354843],
[ 0.06636559, -0.05387809, -0.07661865, 0.00594493, -0.02468248]]), 1: array([[ 0.10535301, -0.04457306, -0.00025789, -0.04025552, 0.0257059 ],
[ 0.06211647, nan, nan, nan, 0.03745972],
[ 0.05411783, nan, nan, nan, -0.04341233],
[-0.05377767, nan, nan, nan, 0.01578407],
[-0.00549912, -0.15478555, 0.01191716, 0.07132315, -0.01376528]])}
/home/runner/.local/lib/python3.12/site-packages/numpy/_core/fromnumeric.py:3860: RuntimeWarning: Mean of empty slice.
return _methods._mean(a, axis=axis, dtype=dtype,
/home/runner/.local/lib/python3.12/site-packages/numpy/_core/_methods.py:137: RuntimeWarning: invalid value encountered in divide
ret = um.true_divide(
plt.figure()
plt.subplot(121)
plt.plot(features["b"], features["a"], label="features")
plt.xlabel("feature b")
plt.ylabel("feature a")
[plt.axvline(b, linewidth=0.5, color='grey') for b in np.linspace(-1, 1, 6)]
[plt.axhline(b, linewidth=0.5, color='grey') for b in np.linspace(-1, 1, 6)]
plt.subplot(122)
extents = (
np.min(features["a"]),
np.max(features["a"]),
np.min(features["b"]),
np.max(features["b"]),
)
plt.imshow(tuning_curves[0],
origin="lower", extent=extents, cmap="viridis",
aspect='auto'
)
plt.title("Tuning curve unit 0")
plt.xlabel("feature b")
plt.ylabel("feature a")
plt.grid(False)
plt.colorbar()
plt.tight_layout()
plt.show()
