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.

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

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)

  0%|          | 0.00/425 [00:00<?, ?MB/s]

  1%|          | 3.00/425 [00:00<00:15, 27.7MB/s]

  2%|▏         | 8.00/425 [00:00<00:10, 40.0MB/s]

  4%|▎         | 15.0/425 [00:00<00:08, 50.6MB/s]

  5%|▍         | 21.0/425 [00:00<00:07, 53.9MB/s]

  7%|▋         | 28.0/425 [00:00<00:07, 55.7MB/s]

  8%|▊         | 34.0/425 [00:00<00:08, 47.1MB/s]

  9%|▉         | 40.0/425 [00:00<00:07, 49.7MB/s]

 11%|█         | 46.0/425 [00:00<00:07, 49.6MB/s]

 12%|█▏        | 52.0/425 [00:01<00:07, 50.2MB/s]

 14%|█▎        | 58.0/425 [00:01<00:07, 50.0MB/s]

 15%|█▌        | 64.0/425 [00:01<00:07, 45.7MB/s]

 16%|█▌        | 69.0/425 [00:01<00:08, 44.5MB/s]

 17%|█▋        | 74.0/425 [00:01<00:09, 38.7MB/s]

 19%|█▊        | 79.0/425 [00:01<00:09, 36.6MB/s]

 20%|█▉        | 84.0/425 [00:01<00:08, 38.4MB/s]

 21%|██        | 88.0/425 [00:02<00:08, 37.5MB/s]

 22%|██▏       | 92.0/425 [00:02<00:08, 37.6MB/s]

 23%|██▎       | 97.0/425 [00:02<00:08, 39.6MB/s]

 24%|██▍       | 102/425 [00:02<00:10, 31.6MB/s] 

 25%|██▍       | 106/425 [00:02<00:11, 27.8MB/s]

 26%|██▌       | 110/425 [00:02<00:11, 27.6MB/s]

 27%|██▋       | 114/425 [00:02<00:10, 29.8MB/s]

 28%|██▊       | 119/425 [00:03<00:09, 32.6MB/s]

 29%|██▉       | 123/425 [00:03<00:08, 33.6MB/s]

 30%|███       | 128/425 [00:03<00:07, 37.2MB/s]

 31%|███       | 132/425 [00:03<00:07, 37.3MB/s]

 33%|███▎      | 139/425 [00:03<00:06, 43.8MB/s]

 34%|███▍      | 146/425 [00:03<00:05, 50.7MB/s]

 36%|███▌      | 152/425 [00:03<00:06, 44.8MB/s]

 37%|███▋      | 158/425 [00:03<00:05, 47.7MB/s]

 39%|███▉      | 166/425 [00:03<00:04, 54.9MB/s]

 40%|████      | 172/425 [00:04<00:04, 55.9MB/s]

 42%|████▏     | 180/425 [00:04<00:03, 62.0MB/s]

 44%|████▍     | 187/425 [00:04<00:03, 63.8MB/s]

 46%|████▌     | 196/425 [00:04<00:03, 70.6MB/s]

 48%|████▊     | 205/425 [00:04<00:03, 72.7MB/s]

 50%|█████     | 213/425 [00:04<00:02, 72.8MB/s]

 52%|█████▏    | 221/425 [00:04<00:02, 71.9MB/s]

 54%|█████▍    | 230/425 [00:04<00:02, 76.4MB/s]

 56%|█████▌    | 239/425 [00:04<00:02, 79.0MB/s]

 58%|█████▊    | 248/425 [00:05<00:02, 79.7MB/s]

 60%|██████    | 257/425 [00:05<00:02, 82.1MB/s]

 63%|██████▎   | 266/425 [00:05<00:01, 80.1MB/s]

 65%|██████▍   | 275/425 [00:05<00:02, 72.2MB/s]

 67%|██████▋   | 283/425 [00:05<00:02, 62.5MB/s]

 68%|██████▊   | 290/425 [00:05<00:02, 62.0MB/s]

 70%|██████▉   | 297/425 [00:05<00:02, 57.9MB/s]

 71%|███████▏  | 303/425 [00:05<00:02, 50.2MB/s]

 73%|███████▎  | 309/425 [00:06<00:02, 49.9MB/s]

 75%|███████▍  | 317/425 [00:06<00:01, 55.5MB/s]

 76%|███████▌  | 324/425 [00:06<00:01, 57.1MB/s]

 78%|███████▊  | 331/425 [00:06<00:01, 59.3MB/s]

 80%|███████▉  | 339/425 [00:06<00:01, 64.7MB/s]

 81%|████████▏ | 346/425 [00:06<00:01, 60.0MB/s]

 84%|████████▎ | 355/425 [00:06<00:01, 67.2MB/s]

 86%|████████▌ | 364/425 [00:06<00:00, 73.1MB/s]

 88%|████████▊ | 372/425 [00:07<00:00, 73.3MB/s]

 90%|████████▉ | 381/425 [00:07<00:00, 77.5MB/s]

 92%|█████████▏| 389/425 [00:07<00:00, 75.6MB/s]

 93%|█████████▎| 397/425 [00:07<00:00, 76.1MB/s]

 95%|█████████▌| 405/425 [00:07<00:00, 72.9MB/s]

 97%|█████████▋| 413/425 [00:07<00:00, 69.4MB/s]

 99%|█████████▉| 421/425 [00:07<00:00, 66.0MB/s]

426MB [00:07, 54.9MB/s]

Loading the data#

Let’s load and print the full dataset.

data = nap.load_file(path)
FS = 1250  # We know from the methods of the paper
print(data)

/home/runner/.local/lib/python3.10/site-packages/hdmf/spec/namespace.py:535: UserWarning: Ignoring cached namespace 'hdmf-common' version 1.7.0 because version 1.8.0 is already loaded.
  warn("Ignoring cached namespace '%s' version %s because version %s is already loaded."
/home/runner/.local/lib/python3.10/site-packages/hdmf/spec/namespace.py:535: UserWarning: Ignoring cached namespace 'core' version 2.6.0-alpha because version 2.7.0 is already loaded.
  warn("Ignoring cached namespace '%s' version %s because version %s is already loaded."
/home/runner/.local/lib/python3.10/site-packages/hdmf/spec/namespace.py:535: UserWarning: Ignoring cached namespace 'hdmf-experimental' version 0.4.0 because version 0.5.0 is already loaded.
  warn("Ignoring cached namespace '%s' version %s because version %s is already loaded."

Achilles_10252013_EEG
┍━━━━━━━━━━━━━┯━━━━━━━━━━━━━┑
│ Keys        │ Type        │
┝━━━━━━━━━━━━━┿━━━━━━━━━━━━━┥
│ units       │ TsGroup     │
│ rem         │ IntervalSet │
│ nrem        │ IntervalSet │
│ forward_ep  │ IntervalSet │
│ eeg         │ TsdFrame    │
│ theta_phase │ Tsd         │
│ position    │ Tsd         │
┕━━━━━━━━━━━━━┷━━━━━━━━━━━━━┙

Selecting slices#

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

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.

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.

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

Text(0.5, 0, 'time (s)')
../_images/07343401e036925fb433e9196e016dae559dc960a9bb469a39fd5c8901b13e40.png

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.

freqs = np.geomspace(5, 200, 25)

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

cwt_rem = nap.compute_wavelet_transform(tsd_rem.restrict(ep_ex_rem), fs=FS, freqs=freqs)

Now let’s plot the calculated wavelet scalogram.

# 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)
../_images/0cf06111c67ecf64c3a1ea6cc994a5930f01277c9a9666721e71b132d9873bb8.png

Filtering Theta#

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

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.

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()
../_images/83320f8d90e82d300483d90179e710a41d317679ed1ed637ac2b631f8a28076b.png

Computing phase#

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

from scipy import signal

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

Let’s plot the phase.

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()
../_images/fa13ee9b9d461999032f7e981be1bd33dc352fe7ad89422b1a905f0261191dd1.png

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.

spikes = spikes[spikes.rate > 5.0]

The feature is the theta phase during REM sleep.

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.

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()
../_images/c23a3be49dbf1eac1b5c0893ec060d1b5c5f254e0963092b06458a4a11591358.png

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.

spike_phase = spikes[spikes.index[3]].value_from(theta_phase)

Let’s plot it.

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()
../_images/48e16de1936b3df18b01dbe6418217dbe5d0d2a9aa7da2ef57c5455ab73a0005.png

Authors

Kipp Freud (https://kippfreud.com/)

Guillaume Viejo