---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.4
kernelspec:
  display_name: Python 3
  language: python
  name: python3
---


Spikes-phase coupling
=====================

In this tutorial we will learn how to isolate phase information using band-pass filtering and combine it
with spiking data, to find phase preferences of spiking units.

Specifically, we will examine LFP and spiking data from a period of REM sleep, after traversal of a linear track.


```{code-cell} ipython3
import math
import os

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
import scipy
import seaborn as sns
import tqdm
import pynapple as nap

custom_params = {"axes.spines.right": False, "axes.spines.top": False}
sns.set_theme(style="ticks", palette="colorblind", font_scale=1.5, rc=custom_params)

```

***
Downloading the data
------------------
Let's download the data and save it locally


```{code-cell} ipython3
path = "Achilles_10252013_EEG.nwb"
if path not in os.listdir("."):
    r = requests.get(f"https://osf.io/2dfvp/download", stream=True)
    block_size = 1024 * 1024
    with open(path, "wb") as f:
        for data in tqdm.tqdm(
            r.iter_content(block_size),
            unit="MB",
            unit_scale=True,
            total=math.ceil(int(r.headers.get("content-length", 0)) // block_size),
        ):
            f.write(data)
```

***
Loading the data
------------------
Let's load and print the full dataset.


```{code-cell} ipython3
data = nap.load_file(path)
FS = 1250  # We know from the methods of the paper
print(data)
```

***
Selecting slices
-----------------------------------
For later visualization, we define an interval of 3 seconds of data during REM sleep.


```{code-cell} ipython3
ep_ex_rem = nap.IntervalSet(
    data["rem"]["start"][0] + 97.0,
    data["rem"]["start"][0] + 100.0,
)
```

Here we restrict the lfp to the REM epochs.


```{code-cell} ipython3
tsd_rem = data["eeg"][:,0].restrict(data["rem"])

# We will also extract spike times from all units in our dataset
# which occur during REM sleep
spikes = data["units"].restrict(data["rem"])
```

***
Plotting the LFP Activity
-----------------------------------
We should first plot our REM Local Field Potential data.


```{code-cell} ipython3
fig, ax = plt.subplots(1, constrained_layout=True, figsize=(10, 3))
ax.plot(tsd_rem.restrict(ep_ex_rem))
ax.set_title("REM Local Field Potential")
ax.set_ylabel("LFP (a.u.)")
ax.set_xlabel("time (s)")
```

***
Getting the Wavelet Decomposition
-----------------------------------
As we would expect, it looks like we have a very strong theta oscillation within our data
- this is a common feature of REM sleep. Let's perform a wavelet decomposition,
as we did in the last tutorial, to see get a more informative breakdown of the
frequencies present in the data.

We must define the frequency set that we'd like to use for our decomposition.


```{code-cell} ipython3
freqs = np.geomspace(5, 200, 25)
```

We compute the wavelet transform on our LFP data (only during the example interval).


```{code-cell} ipython3
cwt_rem = nap.compute_wavelet_transform(tsd_rem.restrict(ep_ex_rem), fs=FS, freqs=freqs)
```

***
Now let's plot the calculated wavelet scalogram.


```{code-cell} ipython3
# Define wavelet decomposition plotting function
def plot_timefrequency(freqs, powers, ax=None):
    im = ax.imshow(np.abs(powers), aspect="auto")
    ax.invert_yaxis()
    ax.set_xlabel("Time (s)")
    ax.set_ylabel("Frequency (Hz)")
    ax.get_xaxis().set_visible(False)
    ax.set(yticks=np.arange(len(freqs))[::2], yticklabels=np.rint(freqs[::2]))
    ax.grid(False)
    return im

fig = plt.figure(constrained_layout=True, figsize=(10, 6))
fig.suptitle("Wavelet Decomposition")
gs = plt.GridSpec(2, 1, figure=fig, height_ratios=[1.0, 0.3])

ax0 = plt.subplot(gs[0, 0])
im = plot_timefrequency(freqs, np.transpose(cwt_rem[:, :].values), ax=ax0)
cbar = fig.colorbar(im, ax=ax0, orientation="vertical")

ax1 = plt.subplot(gs[1, 0])
ax1.plot(tsd_rem.restrict(ep_ex_rem))
ax1.set_ylabel("LFP (a.u.)")
ax1.set_xlabel("Time (s)")
ax1.margins(0)
```

***
Filtering Theta
---------------

As expected, there is a strong 8Hz component during REM sleep. We can filter it using the function `nap.apply_bandpass_filter`.


```{code-cell} ipython3
theta_band = nap.apply_bandpass_filter(tsd_rem, cutoff=(6.0, 10.0), fs=FS)
```

We can plot the original signal and the filtered signal.


```{code-cell} ipython3
plt.figure(constrained_layout=True, figsize=(12, 3))
plt.plot(tsd_rem.restrict(ep_ex_rem), alpha=0.5)
plt.plot(theta_band.restrict(ep_ex_rem))
plt.xlabel("Time (s)")
plt.show()
```

***
Computing phase
---------------

From the filtered signal, it is easy to get the phase using the Hilbert transform. Here we use scipy Hilbert method.


```{code-cell} ipython3
from scipy import signal

theta_phase = nap.Tsd(t=theta_band.t, d=np.angle(signal.hilbert(theta_band)))
```

Let's plot the phase.


```{code-cell} ipython3
plt.figure(constrained_layout=True, figsize=(12, 3))
plt.subplot(211)
plt.plot(tsd_rem.restrict(ep_ex_rem), alpha=0.5)
plt.plot(theta_band.restrict(ep_ex_rem))
plt.subplot(212)
plt.plot(theta_phase.restrict(ep_ex_rem), color='r')
plt.ylabel("Phase (rad)")
plt.xlabel("Time (s)")
plt.show()
```

***
Finding Phase of Spikes
-----------------------
Now that we have the phase of our theta wavelet, and our spike times, we can find the phase firing preferences
of each of the units using the `compute_1d_tuning_curves` function.

We will start by throwing away cells which do not have a high enough firing rate during our interval.


```{code-cell} ipython3
spikes = spikes[spikes.rate > 5.0]
```

The feature is the theta phase during REM sleep.


```{code-cell} ipython3
phase_modulation = nap.compute_1d_tuning_curves(
    group=spikes, feature=theta_phase, nb_bins=61, minmax=(-np.pi, np.pi)
)
```

Let's plot the first 3 neurons.


```{code-cell} ipython3
plt.figure(constrained_layout=True, figsize = (12, 3))
for i in range(3):
    plt.subplot(1,3,i+1)
    plt.plot(phase_modulation.iloc[:,i])
    plt.xlabel("Phase (rad)")
    plt.ylabel("Firing rate (Hz)")
plt.show()
```

There is clearly a strong modulation for the third neuron.
Finally, we can use the function `value_from` to align each spikes to the corresponding phase position and overlay
it with the LFP.


```{code-cell} ipython3
spike_phase = spikes[spikes.index[3]].value_from(theta_phase)
```

Let's plot it.


```{code-cell} ipython3
plt.figure(constrained_layout=True, figsize=(12, 3))
plt.subplot(211)
plt.plot(tsd_rem.restrict(ep_ex_rem), alpha=0.5)
plt.plot(theta_band.restrict(ep_ex_rem))
plt.subplot(212)
plt.plot(theta_phase.restrict(ep_ex_rem), alpha=0.5)
plt.plot(spike_phase.restrict(ep_ex_rem), 'o')
plt.ylabel("Phase (rad)")
plt.xlabel("Time (s)")
plt.show()
```

:::{card}
Authors
^^^
Kipp Freud (https://kippfreud.com/)

Guillaume Viejo

:::