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]
0%| | 1.00/425 [00:00<01:06, 6.36MB/s]
0%| | 2.00/425 [00:00<01:01, 6.85MB/s]
1%| | 3.00/425 [00:00<00:57, 7.28MB/s]
1%| | 4.00/425 [00:00<00:53, 7.85MB/s]
1%| | 5.00/425 [00:00<00:54, 7.72MB/s]
1%|▏ | 6.00/425 [00:00<00:52, 7.95MB/s]
2%|▏ | 7.00/425 [00:00<00:54, 7.73MB/s]
2%|▏ | 8.00/425 [00:01<00:54, 7.61MB/s]
2%|▏ | 9.00/425 [00:01<00:57, 7.21MB/s]
2%|▏ | 10.0/425 [00:01<00:57, 7.16MB/s]
3%|▎ | 11.0/425 [00:01<00:59, 7.00MB/s]
3%|▎ | 12.0/425 [00:01<01:03, 6.48MB/s]
3%|▎ | 13.0/425 [00:01<01:08, 5.97MB/s]
3%|▎ | 14.0/425 [00:02<01:06, 6.14MB/s]
4%|▎ | 15.0/425 [00:02<01:06, 6.15MB/s]
4%|▍ | 16.0/425 [00:02<01:02, 6.57MB/s]
4%|▍ | 17.0/425 [00:02<00:57, 7.11MB/s]
4%|▍ | 18.0/425 [00:02<00:56, 7.20MB/s]
4%|▍ | 19.0/425 [00:02<00:52, 7.66MB/s]
5%|▍ | 20.0/425 [00:02<00:52, 7.71MB/s]
5%|▍ | 21.0/425 [00:02<00:57, 7.08MB/s]
5%|▌ | 22.0/425 [00:03<00:57, 7.01MB/s]
5%|▌ | 23.0/425 [00:03<00:55, 7.31MB/s]
6%|▌ | 24.0/425 [00:03<00:56, 7.05MB/s]
6%|▌ | 25.0/425 [00:03<00:53, 7.44MB/s]
6%|▌ | 26.0/425 [00:03<00:49, 8.05MB/s]
6%|▋ | 27.0/425 [00:03<00:54, 7.36MB/s]
7%|▋ | 28.0/425 [00:03<00:50, 7.82MB/s]
7%|▋ | 29.0/425 [00:04<00:48, 8.18MB/s]
7%|▋ | 30.0/425 [00:04<00:56, 7.00MB/s]
7%|▋ | 31.0/425 [00:04<00:54, 7.28MB/s]
8%|▊ | 32.0/425 [00:04<00:55, 7.08MB/s]
8%|▊ | 33.0/425 [00:04<00:55, 7.04MB/s]
8%|▊ | 34.0/425 [00:04<00:55, 7.06MB/s]
8%|▊ | 35.0/425 [00:04<00:54, 7.19MB/s]
8%|▊ | 36.0/425 [00:05<00:55, 7.07MB/s]
9%|▊ | 37.0/425 [00:05<00:54, 7.17MB/s]
9%|▉ | 38.0/425 [00:05<00:54, 7.16MB/s]
9%|▉ | 39.0/425 [00:05<00:53, 7.19MB/s]
9%|▉ | 40.0/425 [00:05<00:52, 7.31MB/s]
10%|▉ | 41.0/425 [00:05<00:49, 7.73MB/s]
10%|▉ | 42.0/425 [00:05<00:50, 7.51MB/s]
10%|█ | 43.0/425 [00:05<00:48, 7.95MB/s]
10%|█ | 44.0/425 [00:06<00:53, 7.07MB/s]
11%|█ | 45.0/425 [00:06<00:55, 6.82MB/s]
11%|█ | 46.0/425 [00:06<00:54, 6.90MB/s]
11%|█ | 47.0/425 [00:06<00:50, 7.56MB/s]
11%|█▏ | 48.0/425 [00:06<00:47, 7.99MB/s]
12%|█▏ | 49.0/425 [00:06<00:50, 7.46MB/s]
12%|█▏ | 50.0/425 [00:06<00:50, 7.36MB/s]
12%|█▏ | 52.0/425 [00:07<00:46, 8.05MB/s]
12%|█▏ | 53.0/425 [00:07<00:48, 7.67MB/s]
13%|█▎ | 54.0/425 [00:07<00:51, 7.27MB/s]
13%|█▎ | 55.0/425 [00:07<00:55, 6.70MB/s]
13%|█▎ | 56.0/425 [00:07<00:54, 6.78MB/s]
13%|█▎ | 57.0/425 [00:07<00:50, 7.27MB/s]
14%|█▎ | 58.0/425 [00:08<00:55, 6.60MB/s]
14%|█▍ | 59.0/425 [00:08<00:53, 6.84MB/s]
14%|█▍ | 60.0/425 [00:08<00:49, 7.32MB/s]
14%|█▍ | 61.0/425 [00:08<00:48, 7.53MB/s]
15%|█▍ | 62.0/425 [00:08<00:47, 7.61MB/s]
15%|█▍ | 63.0/425 [00:08<00:47, 7.68MB/s]
15%|█▌ | 64.0/425 [00:08<00:46, 7.79MB/s]
15%|█▌ | 65.0/425 [00:08<00:45, 7.90MB/s]
16%|█▌ | 66.0/425 [00:09<00:45, 7.97MB/s]
16%|█▌ | 67.0/425 [00:09<00:45, 7.94MB/s]
16%|█▌ | 68.0/425 [00:09<00:43, 8.14MB/s]
16%|█▌ | 69.0/425 [00:09<00:42, 8.35MB/s]
16%|█▋ | 70.0/425 [00:09<00:40, 8.70MB/s]
17%|█▋ | 72.0/425 [00:09<00:38, 9.14MB/s]
17%|█▋ | 73.0/425 [00:09<00:41, 8.52MB/s]
17%|█▋ | 74.0/425 [00:09<00:40, 8.66MB/s]
18%|█▊ | 75.0/425 [00:10<00:43, 8.01MB/s]
18%|█▊ | 76.0/425 [00:10<00:42, 8.16MB/s]
18%|█▊ | 77.0/425 [00:10<00:41, 8.36MB/s]
18%|█▊ | 78.0/425 [00:10<00:40, 8.59MB/s]
19%|█▊ | 79.0/425 [00:10<00:46, 7.43MB/s]
19%|█▉ | 80.0/425 [00:10<00:47, 7.21MB/s]
19%|█▉ | 81.0/425 [00:10<00:44, 7.67MB/s]
19%|█▉ | 82.0/425 [00:11<00:46, 7.43MB/s]
20%|█▉ | 83.0/425 [00:11<00:48, 7.07MB/s]
20%|█▉ | 84.0/425 [00:11<00:48, 7.10MB/s]
20%|██ | 85.0/425 [00:11<00:46, 7.29MB/s]
20%|██ | 86.0/425 [00:11<00:44, 7.63MB/s]
20%|██ | 87.0/425 [00:11<00:46, 7.26MB/s]
21%|██ | 88.0/425 [00:11<00:44, 7.59MB/s]
21%|██ | 89.0/425 [00:12<00:45, 7.31MB/s]
21%|██ | 90.0/425 [00:12<00:50, 6.61MB/s]
21%|██▏ | 91.0/425 [00:12<00:48, 6.87MB/s]
22%|██▏ | 92.0/425 [00:12<00:45, 7.28MB/s]
22%|██▏ | 93.0/425 [00:12<00:43, 7.55MB/s]
22%|██▏ | 94.0/425 [00:12<00:42, 7.71MB/s]
22%|██▏ | 95.0/425 [00:12<00:43, 7.61MB/s]
23%|██▎ | 96.0/425 [00:12<00:44, 7.43MB/s]
23%|██▎ | 97.0/425 [00:13<00:44, 7.35MB/s]
23%|██▎ | 98.0/425 [00:13<00:45, 7.21MB/s]
23%|██▎ | 99.0/425 [00:13<00:43, 7.44MB/s]
24%|██▎ | 100/425 [00:13<00:42, 7.69MB/s]
24%|██▍ | 101/425 [00:13<00:42, 7.60MB/s]
24%|██▍ | 102/425 [00:13<00:42, 7.65MB/s]
24%|██▍ | 103/425 [00:13<00:42, 7.64MB/s]
24%|██▍ | 104/425 [00:14<00:40, 7.91MB/s]
25%|██▍ | 105/425 [00:14<00:40, 7.81MB/s]
25%|██▍ | 106/425 [00:14<00:43, 7.29MB/s]
25%|██▌ | 107/425 [00:14<00:43, 7.33MB/s]
25%|██▌ | 108/425 [00:14<00:43, 7.37MB/s]
26%|██▌ | 109/425 [00:14<00:45, 6.93MB/s]
26%|██▌ | 110/425 [00:14<00:44, 7.01MB/s]
26%|██▌ | 111/425 [00:15<00:44, 7.04MB/s]
26%|██▋ | 112/425 [00:15<00:46, 6.80MB/s]
27%|██▋ | 113/425 [00:15<00:45, 6.80MB/s]
27%|██▋ | 114/425 [00:15<00:43, 7.21MB/s]
27%|██▋ | 115/425 [00:15<00:40, 7.62MB/s]
27%|██▋ | 116/425 [00:15<00:38, 8.11MB/s]
28%|██▊ | 117/425 [00:15<00:39, 7.81MB/s]
28%|██▊ | 118/425 [00:15<00:41, 7.33MB/s]
28%|██▊ | 119/425 [00:16<00:39, 7.68MB/s]
28%|██▊ | 120/425 [00:16<00:40, 7.51MB/s]
28%|██▊ | 121/425 [00:16<00:41, 7.37MB/s]
29%|██▊ | 122/425 [00:16<00:43, 6.92MB/s]
29%|██▉ | 123/425 [00:16<00:42, 7.18MB/s]
29%|██▉ | 124/425 [00:16<00:40, 7.34MB/s]
29%|██▉ | 125/425 [00:16<00:41, 7.20MB/s]
30%|██▉ | 126/425 [00:17<00:44, 6.77MB/s]
30%|██▉ | 127/425 [00:17<00:42, 6.95MB/s]
30%|███ | 128/425 [00:17<00:42, 6.99MB/s]
30%|███ | 129/425 [00:17<00:41, 7.17MB/s]
31%|███ | 130/425 [00:17<00:39, 7.46MB/s]
31%|███ | 131/425 [00:17<00:38, 7.61MB/s]
31%|███ | 132/425 [00:17<00:42, 6.90MB/s]
31%|███▏ | 133/425 [00:18<00:47, 6.16MB/s]
32%|███▏ | 134/425 [00:18<00:48, 6.06MB/s]
32%|███▏ | 135/425 [00:18<00:45, 6.35MB/s]
32%|███▏ | 136/425 [00:18<00:41, 6.89MB/s]
32%|███▏ | 137/425 [00:18<00:41, 7.02MB/s]
32%|███▏ | 138/425 [00:18<00:42, 6.71MB/s]
33%|███▎ | 139/425 [00:18<00:41, 6.92MB/s]
33%|███▎ | 140/425 [00:19<00:38, 7.37MB/s]
33%|███▎ | 141/425 [00:19<00:40, 7.04MB/s]
33%|███▎ | 142/425 [00:19<00:38, 7.44MB/s]
34%|███▎ | 143/425 [00:19<00:36, 7.68MB/s]
34%|███▍ | 144/425 [00:19<00:36, 7.59MB/s]
34%|███▍ | 145/425 [00:19<00:37, 7.38MB/s]
34%|███▍ | 146/425 [00:19<00:36, 7.69MB/s]
35%|███▍ | 147/425 [00:20<00:35, 7.87MB/s]
35%|███▍ | 148/425 [00:20<00:35, 7.74MB/s]
35%|███▌ | 149/425 [00:20<00:35, 7.81MB/s]
35%|███▌ | 150/425 [00:20<00:34, 7.93MB/s]
36%|███▌ | 151/425 [00:20<00:34, 8.03MB/s]
36%|███▌ | 152/425 [00:20<00:34, 7.83MB/s]
36%|███▌ | 153/425 [00:20<00:34, 8.00MB/s]
36%|███▌ | 154/425 [00:20<00:32, 8.22MB/s]
36%|███▋ | 155/425 [00:20<00:32, 8.32MB/s]
37%|███▋ | 156/425 [00:21<00:32, 8.20MB/s]
37%|███▋ | 157/425 [00:21<00:33, 8.07MB/s]
37%|███▋ | 158/425 [00:21<00:32, 8.17MB/s]
37%|███▋ | 159/425 [00:21<00:31, 8.35MB/s]
38%|███▊ | 160/425 [00:21<00:31, 8.48MB/s]
38%|███▊ | 161/425 [00:21<00:29, 8.80MB/s]
38%|███▊ | 162/425 [00:21<00:31, 8.43MB/s]
38%|███▊ | 163/425 [00:21<00:30, 8.46MB/s]
39%|███▊ | 164/425 [00:22<00:31, 8.31MB/s]
39%|███▉ | 165/425 [00:22<00:37, 6.96MB/s]
39%|███▉ | 166/425 [00:22<00:44, 5.89MB/s]
39%|███▉ | 167/425 [00:22<00:41, 6.22MB/s]
40%|███▉ | 168/425 [00:22<00:39, 6.51MB/s]
40%|███▉ | 169/425 [00:22<00:38, 6.61MB/s]
40%|████ | 170/425 [00:23<00:35, 7.10MB/s]
40%|████ | 171/425 [00:23<00:37, 6.76MB/s]
40%|████ | 172/425 [00:23<00:36, 6.97MB/s]
41%|████ | 173/425 [00:23<00:33, 7.55MB/s]
41%|████ | 174/425 [00:23<00:31, 7.95MB/s]
41%|████ | 175/425 [00:23<00:30, 8.16MB/s]
41%|████▏ | 176/425 [00:23<00:29, 8.34MB/s]
42%|████▏ | 177/425 [00:23<00:29, 8.49MB/s]
42%|████▏ | 178/425 [00:24<00:29, 8.25MB/s]
42%|████▏ | 179/425 [00:24<00:30, 8.03MB/s]
42%|████▏ | 180/425 [00:24<00:33, 7.34MB/s]
43%|████▎ | 181/425 [00:24<00:34, 7.17MB/s]
43%|████▎ | 182/425 [00:24<00:34, 7.04MB/s]
43%|████▎ | 183/425 [00:24<00:35, 6.78MB/s]
43%|████▎ | 184/425 [00:24<00:36, 6.59MB/s]
44%|████▎ | 185/425 [00:25<00:36, 6.53MB/s]
44%|████▍ | 186/425 [00:25<00:34, 6.95MB/s]
44%|████▍ | 187/425 [00:25<00:33, 7.18MB/s]
44%|████▍ | 188/425 [00:25<00:30, 7.73MB/s]
44%|████▍ | 189/425 [00:25<00:28, 8.24MB/s]
45%|████▍ | 190/425 [00:25<00:28, 8.11MB/s]
45%|████▍ | 191/425 [00:25<00:29, 8.02MB/s]
45%|████▌ | 192/425 [00:25<00:27, 8.47MB/s]
45%|████▌ | 193/425 [00:26<00:27, 8.58MB/s]
46%|████▌ | 194/425 [00:26<00:27, 8.43MB/s]
46%|████▌ | 195/425 [00:26<00:28, 8.14MB/s]
46%|████▌ | 196/425 [00:26<00:27, 8.19MB/s]
46%|████▋ | 197/425 [00:26<00:26, 8.54MB/s]
47%|████▋ | 198/425 [00:26<00:26, 8.70MB/s]
47%|████▋ | 199/425 [00:26<00:25, 8.88MB/s]
47%|████▋ | 200/425 [00:26<00:25, 8.92MB/s]
47%|████▋ | 201/425 [00:26<00:25, 8.92MB/s]
48%|████▊ | 202/425 [00:27<00:26, 8.47MB/s]
48%|████▊ | 203/425 [00:27<00:26, 8.23MB/s]
48%|████▊ | 204/425 [00:27<00:26, 8.20MB/s]
48%|████▊ | 205/425 [00:27<00:26, 8.26MB/s]
48%|████▊ | 206/425 [00:27<00:27, 8.02MB/s]
49%|████▊ | 207/425 [00:27<00:26, 8.25MB/s]
49%|████▉ | 208/425 [00:27<00:27, 7.91MB/s]
49%|████▉ | 209/425 [00:27<00:28, 7.53MB/s]
49%|████▉ | 210/425 [00:28<00:34, 6.25MB/s]
50%|████▉ | 211/425 [00:28<00:33, 6.32MB/s]
50%|████▉ | 211/425 [00:28<00:28, 7.43MB/s]
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[2], line 6
4 block_size = 1024 * 1024
5 with open(path, "wb") as f:
----> 6 for data in tqdm.tqdm(
7 r.iter_content(block_size),
8 unit="MB",
9 unit_scale=True,
10 total=math.ceil(int(r.headers.get("content-length", 0)) // block_size),
11 ):
12 f.write(data)
File ~/.local/lib/python3.12/site-packages/tqdm/std.py:1181, in tqdm.__iter__(self)
1178 time = self._time
1180 try:
-> 1181 for obj in iterable:
1182 yield obj
1183 # Update and possibly print the progressbar.
1184 # Note: does not call self.update(1) for speed optimisation.
File /usr/lib/python3/dist-packages/requests/models.py:816, in Response.iter_content.<locals>.generate()
814 if hasattr(self.raw, "stream"):
815 try:
--> 816 yield from self.raw.stream(chunk_size, decode_content=True)
817 except ProtocolError as e:
818 raise ChunkedEncodingError(e)
File /usr/lib/python3/dist-packages/urllib3/response.py:936, in HTTPResponse.stream(self, amt, decode_content)
934 else:
935 while not is_fp_closed(self._fp) or len(self._decoded_buffer) > 0:
--> 936 data = self.read(amt=amt, decode_content=decode_content)
938 if data:
939 yield data
File /usr/lib/python3/dist-packages/urllib3/response.py:879, in HTTPResponse.read(self, amt, decode_content, cache_content)
876 if len(self._decoded_buffer) >= amt:
877 return self._decoded_buffer.get(amt)
--> 879 data = self._raw_read(amt)
881 flush_decoder = amt is None or (amt != 0 and not data)
883 if not data and len(self._decoded_buffer) == 0:
File /usr/lib/python3/dist-packages/urllib3/response.py:814, in HTTPResponse._raw_read(self, amt)
811 fp_closed = getattr(self._fp, "closed", False)
813 with self._error_catcher():
--> 814 data = self._fp_read(amt) if not fp_closed else b""
815 if amt is not None and amt != 0 and not data:
816 # Platform-specific: Buggy versions of Python.
817 # Close the connection when no data is returned
(...)
822 # not properly close the connection in all cases. There is
823 # no harm in redundantly calling close.
824 self._fp.close()
File /usr/lib/python3/dist-packages/urllib3/response.py:799, in HTTPResponse._fp_read(self, amt)
796 return buffer.getvalue()
797 else:
798 # StringIO doesn't like amt=None
--> 799 return self._fp.read(amt) if amt is not None else self._fp.read()
File /usr/lib/python3.12/http/client.py:479, in HTTPResponse.read(self, amt)
476 if self.length is not None and amt > self.length:
477 # clip the read to the "end of response"
478 amt = self.length
--> 479 s = self.fp.read(amt)
480 if not s and amt:
481 # Ideally, we would raise IncompleteRead if the content-length
482 # wasn't satisfied, but it might break compatibility.
483 self._close_conn()
File /usr/lib/python3.12/socket.py:707, in SocketIO.readinto(self, b)
705 while True:
706 try:
--> 707 return self._sock.recv_into(b)
708 except timeout:
709 self._timeout_occurred = True
File /usr/lib/python3.12/ssl.py:1252, in SSLSocket.recv_into(self, buffer, nbytes, flags)
1248 if flags != 0:
1249 raise ValueError(
1250 "non-zero flags not allowed in calls to recv_into() on %s" %
1251 self.__class__)
-> 1252 return self.read(nbytes, buffer)
1253 else:
1254 return super().recv_into(buffer, nbytes, flags)
File /usr/lib/python3.12/ssl.py:1104, in SSLSocket.read(self, len, buffer)
1102 try:
1103 if buffer is not None:
-> 1104 return self._sslobj.read(len, buffer)
1105 else:
1106 return self._sslobj.read(len)
KeyboardInterrupt:
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)
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)")
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)
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()
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()
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()
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()
Authors
Kipp Freud (https://kippfreud.com/)
Guillaume Viejo