Tuning curves#

With Pynapple you can easily compute n-dimensional tuning curves (for example, firing rate as a function of 1D angular direction or firing rate as a function of 2D position). It is also possible to compute average firing rate for different epochs (for example, firing rate for different epochs of stimulus presentation).

Hide code cell content

 
import pynapple as nap
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import xarray as xr
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)
xr.set_options(display_expand_attrs=False)

<xarray.core.options.set_options at 0x7f5b4b985910>

From timestamps or continuous activity#

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)
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)
    )

Computing tuning curves is done using compute_tuning_curves.

When computing from general time-series, mandatory arguments are:

  • data: a TsGroup (or single Ts) or TsdFrame (or single Tsd) containing the neural activity of one or more units.

  • features: a Tsd or TsdFrame containing one or more features.

By default, 10 bins are used for all features, but you can specify the number of bins, or the bin edges explicitly, using the bins argument.

The min and max of the tuning curves are by default the minima and maxima of the features. This can be tweaked with the range argument.

If an IntervalSet is passed with epochs, everything is restricted to epochs, otherwise the time support of the features is used.

If you do not want the sampling rate of the features to be estimated from the timestamps, you can pass it explicitly using the fs argument.

You can further also pass a list of strings to label each dimension via feature_names (by default the columns of the features are used).

The output is an xarray.DataArray in which the first dimension represents the units and further dimensions represent the features. The occupancy and bin edges are stored as attributes.

If you explicitly want a pd.DataFrame as output (which is only possible when you have just the one feature), you can set return_pandas=True. Note that this will not return the occupancy and bin edges.

1D tuning curves from spikes#

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

<xarray.DataArray (unit: 6, feature: 120)> Size: 6kB
array([[315.50218341, 321.66301969, 314.61864407, 302.42825607,
        286.20689655, 290.14989293, 250.54704595, 269.8019802 ,
        234.66981132, 218.42105263, 171.71717172, 171.79487179,
        133.81995134, 132.8320802 , 117.04545455, 110.97852029,
         75.17899761,  79.12087912,  35.86956522,  48.31460674,
         20.05012531,  18.56148492,  12.01923077,  17.5257732 ,
          9.25925926,  12.01923077,   4.95049505,   3.42465753,
          2.36966825,   4.47427293,   0.        ,   2.42718447,
          1.2345679 ,   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.        ,   0.        ,   0.        ,   0.        ,
          0.        ,   0.        ,   0.        ,   0.        ,
          0.        ,   0.        ,   0.        ,   0.        ,
          0.        ,   0.        ,   1.24688279,   0.        ,
          0.        ,   4.23728814,   1.44092219,   4.12087912,
          4.06504065,   3.04878049,  13.43283582,   9.09090909,
         20.83333333,  17.41293532,  21.86588921,  28.76712329,
         54.34782609,  62.65060241,  66.13756614,  86.73469388,
        112.82660333,  98.25870647, 135.33834586, 143.37349398,
        194.1031941 , 181.93717277, 220.69825436, 234.29951691,
        259.30232558, 256.09756098, 289.28571429, 304.03800475,
        284.09090909, 309.82367758, 339.37823834, 316.3507109 ,
        308.31099196, 322.54196643, 292.28855721, 306.97674419,
        289.93435449, 294.11764706, 273.69668246, 251.10619469,
        228.87323944, 224.05660377, 175.37313433, 184.7826087 ,
        144.41747573, 132.31552163, 108.59728507, 111.75115207,
         62.94536817,  75.37154989,  41.1622276 ,  34.68899522]])
Coordinates:
  * unit     (unit) int64 48B 0 1 2 3 4 5
  * feature  (feature) float64 960B 0.02618 0.07854 0.1309 ... 6.152 6.205 6.257
Attributes: (3)

The xarray.DataArray can be treated like a numpy array.

It has a shape:

tuning_curves_1d.shape

(6, 120)

It can be sliced:

tuning_curves_1d[1, 2:8]

<xarray.DataArray (feature: 6)> Size: 48B
array([ 73.09322034,  73.95143488, 103.44827586, 104.92505353,
       133.47921225, 139.85148515])
Coordinates:
  * feature  (feature) float64 48B 0.1309 0.1833 0.2356 0.288 0.3403 0.3927
    unit     int64 8B 1
Attributes: (3)

It can also be indexed using the coordinates:

tuning_curves_1d.sel(unit=1)

<xarray.DataArray (feature: 120)> Size: 960B
array([ 45.85152838,  56.89277899,  73.09322034,  73.95143488,
       103.44827586, 104.92505353, 133.47921225, 139.85148515,
       188.67924528, 178.94736842, 220.95959596, 250.        ,
       263.99026764, 258.14536341, 298.86363636, 272.07637232,
       317.42243437, 287.91208791, 345.65217391, 323.59550562,
       314.53634085, 300.46403712, 317.30769231, 275.25773196,
       282.40740741, 289.66346154, 261.13861386, 250.        ,
       215.63981043, 229.3064877 , 205.65552699, 174.75728155,
       120.98765432, 142.15686275,  98.12108559,  92.13483146,
        70.        ,  89.93576017,  44.11764706,  55.26992288,
        30.4136253 ,  32.15077605,  20.83333333,  11.21076233,
        10.66098081,   6.00961538,   0.        ,   6.83371298,
         0.        ,   1.12359551,   2.27272727,   1.2195122 ,
         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.1627907 ,
         0.        ,   0.        ,   2.36966825,   0.        ,
         1.17370892,   1.17924528,   3.73134328,   2.41545894,
         2.42718447,   5.08905852,  10.18099548,   8.06451613,
        14.25178147,  12.7388535 ,  13.31719128,  20.33492823])
Coordinates:
  * feature  (feature) float64 960B 0.02618 0.07854 0.1309 ... 6.152 6.205 6.257
    unit     int64 8B 1
Attributes: (3)

xarray further has matplotlib support, allowing for easy visualization:

tuning_curves_1d.plot.line(x="feature", add_legend=False)
plt.ylabel("Firing rate (Hz)")
plt.show()
../_images/d95541656d201763fe9147c2876691cb4c3725a8bb08d10b93d1132f6e83825d.png

You can either customize the plot labels yourself using matplotlib, or you can set them in the tuning curve object:

tuning_curves_1d.name = "Firing rate"
tuning_curves_1d.attrs["unit"] = "Hz"
tuning_curves_1d.coords["feature"].attrs["unit"] = "rad"
tuning_curves_1d.plot.line(x="feature", add_legend=False)
plt.show()
../_images/6ffa2a0a2dcf5eb163926179e493e976f7ebdc259b2e925d696c5fe026faa8b5.png

Internally, the compute_tuning_curves calls the value_from method which maps timestamps to their closest values 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.

Hide 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_curves_1d[3].values, tuning_curves_1d.coords["feature"], label="Tuning curve (unit=3)")
plt.xlabel("Firing rate (Hz)")
plt.legend()
plt.show()
../_images/ad0e0f20cd2eba5efcddad13b37745496bf7ca7f84dd5a411e7d9cf46c5eb794.png

It is also possible to just get the spike counts per bins. This can be done by setting the argument return_counts=True. The output is also a xarray.DataArray with the same dimensions as the tuning curves.

spike_counts = nap.compute_tuning_curves(
    data=tsgroup,
    features=feature,
    bins=30, 
    range=(0, 2*np.pi),
    feature_names=["feature"],
    return_counts=True
    )

Hide code cell source

 
plt.figure()
plt.subplot(131)
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(132)
plt.plot(tuning_curves_1d[3].values, tuning_curves_1d.coords["feature"])
plt.xlabel("Firing rate (Hz)")
plt.subplot(133)
plt.barh(spike_counts.coords["feature"], width=spike_counts[3].values, height=np.mean(np.diff(spike_counts.coords["feature"])))
plt.xlabel("Spike count")
plt.tight_layout()
plt.show()
../_images/39f5d089b423741e724eb1f8bf8b09e8863c6dfd83bc4ea982df1aa5fe30313f.png

2D tuning curves from spikes#

Hide 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)

If you pass more than 1 feature, a multi-dimensional tuning curve is computed:

tuning_curves_2d = nap.compute_tuning_curves(
    data=tsgroup, 
    features=features, 
    bins=(5,5),
    range=[(-1, 1), (-1, 1)],
    feature_names=["a", "b"]
)
tuning_curves_2d

<xarray.DataArray (unit: 3, a: 5, b: 5)> Size: 600B
array([[[0.        , 1.2345679 , 3.75      , 3.37078652, 0.        ],
        [0.        ,        nan,        nan,        nan, 1.13636364],
        [0.        ,        nan,        nan,        nan, 0.        ],
        [2.27272727,        nan,        nan,        nan, 0.        ],
        [0.        , 0.        , 1.63934426, 0.        , 0.        ]],

       [[0.        , 0.        , 1.25      , 2.24719101, 5.26315789],
        [0.        ,        nan,        nan,        nan, 2.27272727],
        [0.        ,        nan,        nan,        nan, 1.2345679 ],
        [0.        ,        nan,        nan,        nan, 2.27272727],
        [3.57142857, 4.44444444, 0.        , 1.13636364, 0.        ]],

       [[0.        , 1.2345679 , 1.25      , 2.24719101, 1.75438596],
        [0.        ,        nan,        nan,        nan, 0.        ],
        [2.5       ,        nan,        nan,        nan, 1.2345679 ],
        [2.27272727,        nan,        nan,        nan, 5.68181818],
        [0.        , 0.        , 1.63934426, 5.68181818, 1.75438596]]])
Coordinates:
  * unit     (unit) int64 24B 0 1 2
  * a        (a) float64 40B -0.8 -0.4 1.11e-16 0.4 0.8
  * b        (b) float64 40B -0.8 -0.4 1.11e-16 0.4 0.8
Attributes: (3)

tuning_curve_2d is a again an xarray.DataArray but now with three dimensions: one for the units of TsGroup and 2 for the features, the coordinates contain the centers of the bins. Bins that have never been visited by the feature have been assigned a NaN value.

Two-dimensional tuning curves can also easily be visualized:

tuning_curves_2d.name="Firing rate"
tuning_curves_2d.attrs["unit"]="Hz"
tuning_curves_2d.plot(col="unit")
plt.show()
../_images/056ba0d7255f55eccb5c81dbcd9e0d26ca450f465b7ca4490224b7909592deda.png

Verifying the accuracy of the tuning curves can once more be done 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
----------  ---------  ----------
2.65129     -0.881582   0.472031
2.83797     -0.954865   0.297041
3.10551     -0.999501   0.0315874
3.16223     -0.999831  -0.0184063
5.33695      0.587213  -0.809433
6.16614      0.993601  -0.112944
8.27051     -0.404122   0.914705
9.11071     -0.950865   0.309605
9.46675     -0.998978  -0.0452066
9.67952     -0.967607  -0.25246
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.

Hide code cell source

 
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8,4), sharey=True)
ax1.plot(features["b"], features["a"], label="features")
ax1.plot(ts_to_features["b"], ts_to_features["a"], "o", color="red", markersize=4, label="spikes")
ax1.set_xlabel("b")
ax1.set_ylabel("a")
[ax1.axvline(b, linewidth=0.5, color='grey') for b in np.linspace(-1, 1, 6)]
[ax1.axhline(b, linewidth=0.5, color='grey') for b in np.linspace(-1, 1, 6)]
extents = (
    np.min(features["a"]),
    np.max(features["a"]),
    np.min(features["b"]),
    np.max(features["b"]),
)
tuning_curves_2d[0].plot(ax=ax2)
ax2.set_ylabel("")
plt.tight_layout()
plt.show()
../_images/02cf2035ed8713446c9531c7c56704f6e85b5a988b97f15d523281d22c4f08d5.png

1D tuning curves from continuous activity#

Hide 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
    )

We do not always have spikes. Sometimes we are analysing continuous firing rates or calcium intensities. In that case, we can simply pass a Tsd or TsdFrame as group:

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

<xarray.DataArray (unit: 3, feature: 120)> Size: 3kB
array([[ 9.86157522e-01,  9.12894950e-01,  1.02917785e+00,
         9.36568217e-01,  9.32748301e-01,  8.12913728e-01,
         7.31943442e-01,  7.66946244e-01,  5.93173785e-01,
         5.73914098e-01,  3.83407412e-01,  3.98264196e-01,
         2.57429770e-01,  3.18848808e-01,  1.95130197e-01,
         2.93831580e-01,  1.87797297e-01,  1.79708718e-01,
         8.06631337e-02,  7.60499854e-02,  7.34629691e-02,
         5.08527582e-02,  3.00228219e-02,  2.08296228e-03,
         2.94358824e-03, -1.57981577e-02, -4.10941373e-02,
         3.75697557e-02,  4.78705320e-03, -3.43962395e-02,
        -1.70741407e-02, -8.54123420e-02,  6.43226104e-03,
        -1.15110635e-02, -1.64099029e-02, -8.16751202e-02,
        -1.33625201e-02, -2.73356684e-02,  1.77479920e-02,
         2.38781891e-02,  7.25649735e-02, -1.54870553e-01,
         2.83207908e-02,  4.45810516e-02, -7.41878069e-02,
         1.57975545e-02,  1.36160499e-02,  2.13802277e-02,
         4.19703626e-02,  1.38995829e-02, -3.86904221e-02,
         8.18219099e-02, -1.46206590e-02,  7.38805708e-02,
        -7.36581160e-03, -6.35111668e-02,  9.78018752e-03,
        -1.42664646e-01,  1.90335738e-02, -1.63731934e-02,
...
         1.44632306e-02,  1.56546566e-01,  1.91536056e-01,
         2.34762874e-01,  2.58923361e-01,  2.22623688e-01,
         3.13057726e-01,  3.34368697e-01,  5.15979494e-01,
         5.46292055e-01,  5.82710930e-01,  6.51420347e-01,
         7.89006626e-01,  7.74163225e-01,  9.30402991e-01,
         8.92839622e-01,  9.43730354e-01,  9.95088898e-01,
         9.84832708e-01,  9.59882204e-01,  9.72181846e-01,
         9.96994805e-01,  9.95076017e-01,  9.72209342e-01,
         7.76059664e-01,  8.59659983e-01,  7.10026950e-01,
         7.21771634e-01,  6.28968652e-01,  5.72144832e-01,
         4.38091668e-01,  4.77474957e-01,  3.42009271e-01,
         3.91889544e-01,  1.70109510e-01,  1.46810333e-01,
         1.72948177e-01,  1.42247099e-01,  1.32359302e-02,
         1.06725243e-01,  3.11522981e-02,  1.26289004e-01,
         4.80906344e-02,  6.40675911e-02, -5.08200487e-02,
         2.28327496e-02,  4.74275749e-02,  1.13854281e-02,
         6.25525727e-02,  6.22069197e-02, -5.06064101e-02,
         4.92245100e-02,  3.45255593e-02,  3.86506712e-02,
        -7.47408003e-02, -9.28849176e-02, -2.71762927e-03,
         8.96027501e-02, -3.24250616e-02,  6.52579167e-02]])
Coordinates:
  * unit     (unit) int64 24B 0 1 2
  * feature  (feature) float64 960B 0.02618 0.07854 0.1309 ... 6.152 6.205 6.257
Attributes: (3)
tuning_curves_1d.name="ΔF/F"
tuning_curves_1d.attrs["unit"]="a.u."
tuning_curves_1d.plot.line(x="feature", add_legend=False)
plt.show()
../_images/8d48735d2df313d853f7e8f8b38f059c782e674bea2ecce7d5bf2655ee730f33.png

2D tuning curves from continuous activity#

This also works with more than one feature:

Hide 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
ft = features.values
alpha = np.arctan2(ft[:, 1], ft[:, 0])
bin_centers = np.linspace(-np.pi, np.pi, 6)
kappa = 4.0
units=[]
for i, mu in enumerate(bin_centers):
    units.append(np.exp(kappa * np.cos(alpha - mu))) # wrapped Gaussian
units = np.stack(units, axis=1)
tsdframe = nap.TsdFrame(t=features.times(), d=units)

tuning_curves_2d = nap.compute_tuning_curves(
    data=tsdframe,
    features=features,
    bins=5,
    feature_names=["a", "b"]
    )
tuning_curves_2d

<xarray.DataArray (unit: 6, a: 5, b: 5)> Size: 1kB
array([[[1.72951050e+01, 3.91144332e+01, 5.31794858e+01, 3.81958293e+01,
         1.71975317e+01],
        [5.66125302e+00,            nan,            nan,            nan,
         5.64344925e+00],
        [1.09770223e+00,            nan,            nan,            nan,
         1.11244979e+00],
        [2.18237399e-01,            nan,            nan,            nan,
         2.18436481e-01],
        [6.17222795e-02, 2.75125688e-02, 1.88193311e-02, 2.74574544e-02,
         6.14029659e-02]],

       [[3.51467243e+01, 1.46084750e+01, 3.77319431e+00, 7.44172335e-01,
         1.70652443e-01],
        [5.16817365e+01,            nan,            nan,            nan,
         5.61114135e-02],
        [4.40798235e+01,            nan,            nan,            nan,
         2.31135316e-02],
        [1.99841156e+01,            nan,            nan,            nan,
         1.93967943e-02],
        [6.42751592e+00, 1.75046459e+00, 2.80401048e-01, 7.44271795e-02,
...
         3.52864865e+01],
        [5.62880827e-02,            nan,            nan,            nan,
         5.17474811e+01],
        [2.30117522e-02,            nan,            nan,            nan,
         4.42029261e+01],
        [1.93982667e-02,            nan,            nan,            nan,
         1.99958150e+01],
        [2.88192489e-02, 7.47912707e-02, 3.64106970e-01, 1.74279428e+00,
         6.38482631e+00]],

       [[1.72951050e+01, 3.91144332e+01, 5.31794858e+01, 3.81958293e+01,
         1.71975317e+01],
        [5.66125302e+00,            nan,            nan,            nan,
         5.64344925e+00],
        [1.09770223e+00,            nan,            nan,            nan,
         1.11244979e+00],
        [2.18237399e-01,            nan,            nan,            nan,
         2.18436481e-01],
        [6.17222795e-02, 2.75125688e-02, 1.88193311e-02, 2.74574544e-02,
         6.14029659e-02]]])
Coordinates:
  * unit     (unit) int64 48B 0 1 2 3 4 5
  * a        (a) float64 40B -0.8 -0.4 6.341e-07 0.4 0.8
  * b        (b) float64 40B -0.8 -0.4 1.268e-06 0.4 0.8
Attributes: (3)
tuning_curves_2d.name="ΔF/F"
tuning_curves_2d.attrs["unit"]="a.u."
tuning_curves_2d.plot(col="unit", col_wrap=3)
plt.show()
../_images/ca5081615f5b471c393e1623af47b808054c65ecff856554508a23e9d8571333.png

From epochs#

When computing from epochs, you should store them in a dictionary:

epochs_dict =  {
    "stim0": nap.IntervalSet(start=0, end=20),
    "stim1":nap.IntervalSet(start=30, end=70)
}

You can then compute the tuning curves using nap.compute_response_per_epoch. You can pass either a TsGroup for spikes, or a TsdFrame for rates/calcium activity.

The output is an xarray.DataArray with labeled dimensions:

tuning_curves = nap.compute_response_per_epoch(tsgroup, epochs_dict)
tuning_curves

<xarray.DataArray (unit: 3, epochs: 2)> Size: 48B
array([[0.5 , 0.  ],
       [0.75, 0.  ],
       [1.  , 0.  ]])
Coordinates:
  * unit     (unit) int64 24B 0 1 2
  * epochs   (epochs) <U5 40B 'stim0' 'stim1'

Mutual information#

Given a set of tuning curves, you can use compute_mutual_information to compute the mutual information between the activity of the neurons and the features, no matter what dimension. See the Skaggs et al. (1992) paper for more information on what mutual information computes.

MI = nap.compute_mutual_information(tuning_curves_1d)
MI
bits/sec bits/spike
0 0.304127 1.665517
1 0.317359 1.855041
2 0.307325 1.646332 
MI = nap.compute_mutual_information(tuning_curves_2d)
MI
bits/sec bits/spike
0 17.622510 1.333944
1 14.124129 1.809407
2 14.131191 1.776403
3 17.507148 1.308321
4 16.106702 1.135644
5 17.622510 1.333944

Take a look at the tutorial on head direction cells for a realistic example.