Reproducing Tingley and Buzsaki (2019) from dandiset 000218¶
Paper: Tingley D, Buzsaki G (2019). Routing of hippocampal ripples to subcortical structures via the lateral septum. Neuron 105, 138-149. doi:10.1016/j.neuron.2019.10.012. The full text was used: the accepted author manuscript in PubMed Central (PMC6952543). Its supplementary figures were not available, so claims on Figures S1 to S5 rest on the main text.
Dandiset: DANDI:000218 version 0.250624.0439, 98 NWB files (1.5 TB) from 6 rats, streamed from the archive. The paper cites no version (it says the data would go to CRCNS). There are two published versions, 0.220131.1608 and 0.250624.0439; they hold the same 98 files with the same blob ids (checked below) and differ only in contributor metadata, so the most recent one is pinned.
The files hold the LFP at 1250 Hz on every channel, sorted units with a region label in most files,
sleep states, epochs and trials. They hold no detected events. LS high-frequency oscillations
(HFOs) and CA1 ripples are detected here from the LFP, following the Methods and the buzcode
detector they name (see ls_hfo.py). The paper's own code repository holds no analysis of this
paper apart from the network model of Figure 8.
| Claim | Tested here |
|---|---|
| dataset-composition: 6 rats, 1989 LS neurons, 37 CA1+LS and 25 CA3+LS recordings | tested |
| fig1c-hfo-each-animal: HFOs with increased LS firing in each of 6 animals | tested |
| fig1c-event-count: 29,955 HFOs | tested |
| figs1-participation: 16% of LS neurons per HFO, 0.23 spikes per event | tested |
| fig1d-phase-locking: 92% of LS neurons phase-locked, trough preference | tested |
| fig1f-shorter-than-ripples: LS HFOs 27 ms against CA1 ripples 31 ms | tested |
| fig1g-slower-than-ripples: 147 against 159 Hz | tested |
| fig1h-no-amplitude-frequency-correlation | tested |
| fig2-confined-to-ls: HFOs only at LS sites | tested |
| fig2ef-intermediate-ls: power and locking peak in intermediate LS | blocked (metadata insufficient) |
| fig3-travelling-waves: propagation speeds | blocked (metadata insufficient) |
| fig4b-pre-post-unchanged: excess correlation PRE = POST | tested |
| fig4b-hfo-above-behav: HFO excess correlation above BEHAV and jitter | tested |
| fig4-sessions-few-change: 7 of 53 sessions differ | blocked after selection (too few recordings) |
| fig4c-temporal-bias: PRE and POST biases similar | blocked after selection (too few recordings) |
| fig4d-compression: above chance up to 20-fold | tested |
| figs2-post-rate-increase: rates rise in post-behavior NREM | tested |
| fig5b-nrem-and-rates: HFOs mostly in NREM, rarer than ripples | tested |
| fig5d-lag: CA1 leads LS by 6.4 ms | tested |
| fig5e-state-coupling: coupling 0.3 in NREM against 0.15 in WAKE | tested |
| fig5f-ripple-triggered-firing: higher in NREM | tested |
| fig5f-no-phase-coherence: CA1-LS coherence at chance | tested |
| fig6-optogenetic-hfos | blocked (data missing) |
| fig7e-ca1-spikes: 38% of LS neurons | tested |
| fig7e-state: 26% | tested |
| fig7e-duration: about 20% | tested |
| fig7e-ca1-depth: 23% | blocked (data missing) |
| fig7e-replay: 7.9, 6.6 and 7.1% | blocked (methods underspecified) |
| fig7e-seqnmf: 8.1% | blocked (methods underspecified) |
| fig7e-content-independent | blocked (depends on blocked claims) |
| fig8-model | blocked (other: a simulation, no recorded data) |
The criteria are in claims/000218/10.1016-j.neuron.2019.10.012.yaml; the outcomes and the
reasons are in attempts/000218-tingley-2019-01.yaml.
In [1]:
import collections
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
from scipy import stats
import ls_hfo
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
In [2]:
# Stream every file of the pinned version once and keep what the analyses need (about 17 min
# with 6 workers; most of it is loading each file's chunk index over the network).
started = time.time()
asset_list = sorted(ls_hfo.assets(), key=lambda a: a["path"])
results = ls_hfo.run_all([a["path"] for a in asset_list], workers=6)
print(f"dandiset {ls_hfo.DANDISET} version {ls_hfo.VERSION}: {len(asset_list)} files, "
f"{sum(a['size'] for a in asset_list) / 1e12:.2f} TB")
print(f"streamed and processed in {(time.time() - started) / 60:.1f} min")
errors = [r for r in results if "error" in r]
print("files that failed:", len(errors))
for r in errors:
print(" ", r["path"], r["error"])
skipped = [r for r in results if "skipped" in r]
print("files without LS electrodes (no LFP analysis):", len(skipped))
for r in skipped:
print(" ", r["path"], dict(collections.Counter(r["meta"]["el_location"])))
R = [r for r in results if "error" not in r and "skipped" not in r]
# The two published versions hold the same files: compare paths and blob ids.
older = {a["path"]: a["blob"] for a in ls_hfo.assets("0.220131.1608")}
pinned = {a["path"]: a["blob"] for a in asset_list}
print("0.220131.1608 and the pinned version hold the same paths and blobs:", older == pinned)
dandiset 000218 version 0.250624.0439: 98 files, 1.51 TB
streamed and processed in 23.9 min
files that failed: 0
files without LS electrodes (no LFP analysis): 4
sub-DT7/sub-DT7_ses-20170327-864um-360um-170327-115803_ecephys.nwb {np.str_('unknown'): 128}
sub-DT7/sub-DT7_ses-20170406-864um-1152um-merge_behavior+ecephys.nwb {np.str_('unknown'): 128}
sub-DT7/sub-DT7_ses-20170407-936um-1152um-merge_behavior+ecephys.nwb {np.str_('unknown'): 128}
sub-DT7/sub-DT7_ses-20170408-936um-1152um-merge_ecephys.nwb {np.str_('unknown'): 128}
0.220131.1608 and the pinned version hold the same paths and blobs: True
Data check¶
Each check of the description against the values is stated under its output.
In [3]:
# Inventory: what each file holds.
def region_counts(values):
return dict(sorted(collections.Counter(values).items()))
inv = []
for r in results:
m = r["meta"]
loc = m["units"]["location"]
inv.append({
"file": r["path"].split("/")[1].replace("_ecephys.nwb", "").replace("_behavior+ecephys.nwb", "")[:48],
"subject": m["subject"], "hours": round(m["duration"] / 3600, 2),
"electrodes": region_counts(m["el_location"]),
"units": len(m["units"]["spikes"]),
"unit_regions": "no location column" if loc is None else region_counts(loc),
"states": m["states"] is not None, "epochs": m["epochs"] is not None,
"trials": 0 if m["trials"] is None else len(m["trials"]["start"]),
})
inv = pd.DataFrame(inv)
print(inv.to_string())
print()
print("files per subject:", inv.subject.value_counts().sort_index().to_dict())
unit_regions = collections.Counter()
for r in results:
loc = r["meta"]["units"]["location"]
if loc is None:
unit_regions["(no location)"] += len(r["meta"]["units"]["spikes"])
else:
unit_regions.update(loc)
print("units by region label:", dict(unit_regions))
print("files whose units have no location column:",
[r["path"].split("/")[1][:40] for r in results if r["meta"]["units"]["location"] is None and len(r["meta"]["units"]["spikes"])])
print("LFP rates (Hz):", sorted({r["meta"]["lfp_rate"] for r in results}),
"| conversion:", sorted({r["meta"]["lfp_conversion"] for r in results}), sorted({r["meta"]["lfp_unit"] for r in results}))
print("electrode x, y, z all NaN in every file:", all(np.isnan(r["meta"]["el_xyz"]).all() for r in results))
print("LFP starting_time values (s):", sorted({r["meta"]["lfp_starting_time"] for r in results}))
print("units-table columns:", sorted({tuple(r["meta"]["units"]["columns"]) for r in results}))
print("first unit of the first file with a group_id column, by animal:")
for subj in sorted({r["meta"]["subject"] for r in results}):
r = next(r for r in results if r["meta"]["subject"] == subj and "group_id" in r["meta"]["units"]["first_row"])
print(" ", subj, r["meta"]["units"]["first_row"])
unit_groups_ok = all(set(r["meta"]["units"]["group"]) <= set(r["meta"]["el_group"]) for r in results if r["meta"]["units"]["group"] is not None)
print("every unit group_id is an electrode group of its file:", unit_groups_ok)
print("electrode group descriptions:", sorted({g["description"].split()[-1] for r in results for g in r["meta"]["electrode_groups"].values()}),
"| e.g.", sorted({g["description"] for r in results for g in r["meta"]["electrode_groups"].values()})[:2])
print("devices (name: description), over all files:")
for name, desc in sorted({(k, v) for r in results for k, v in r["meta"]["devices"].items()}):
print(" ", name, ":", desc)
file subject hours electrodes units unit_regions states epochs trials
0 sub-DT1_ses-DT1-rLS-20150723-1584um DT1 3.03 {'ls': 30} 74 {'ls': 74} False False 0
1 sub-DT1_ses-DT1-rLS-20150723-1728um DT1 2.19 {'ls': 30} 53 {'ls': 53} False False 0
2 sub-DT1_ses-DT1-rLS-20150724-1800um DT1 2.61 {'ls': 30} 51 {'ls': 51} False False 0
3 sub-DT1_ses-DT1-rLS-20150724-1872um DT1 3.51 {'ls': 30} 60 {'ls': 60} False False 0
4 sub-DT1_ses-DT1-rLS-20150725-1944um DT1 4.68 {'ls': 30} 43 {'ls': 43} False False 0
5 sub-DT1_ses-DT1-rLS-20150726-2016um DT1 3.51 {'ls': 30} 40 {'ls': 40} False False 0
6 sub-DT1_ses-DT1-rLS-20150727-2088um DT1 3.07 {'ls': 30} 44 {'ls': 44} False False 0
7 sub-DT1_ses-DT1-rLS-20150727-2232um DT1 3.49 {'ls': 30} 54 {'ls': 54} False False 0
8 sub-DT1_ses-DT1-rLS-20150728-2304um DT1 5.47 {'ls': 30} 43 {'ls': 43} False False 0
9 sub-DT1_ses-DT1-rLS-20150729-2448um DT1 2.92 {'ls': 30} 43 {'ls': 43} False False 0
10 sub-DT1_ses-DT1-rLS-20150729-2592um DT1 4.24 {'ls': 30} 50 {'ls': 50} False False 0
11 sub-DT1_ses-DT1-rLS-20150730-2736um DT1 3.12 {'ls': 30} 55 {'ls': 55} False False 0
12 sub-DT2_ses-DT2-rPPC-rCCG-3036um-362um-20160215- DT2 3.65 {'hpc': 64, 'ls': 64} 37 {'ls': 37} True True 77
13 sub-DT2_ses-DT2-rPPC-rCCG-3108um-434um-20160215- DT2 2.65 {'hpc': 64, 'ls': 64} 0 no location column False False 0
14 sub-DT2_ses-DT2-rPPC-rCCG-3324um-1000um-20160217 DT2 7.25 {'hpc': 64, 'ls': 64} 56 {'ls': 28, 'ppc': 28} True True 69
15 sub-DT2_ses-DT2-rPPC-rCCG-3324um-1144um-20160219 DT2 6.05 {'hpc': 64, 'ls': 64} 88 {'hpc': 67, 'ls': 21} True True 84
16 sub-DT2_ses-DT2-rPPC-rCCG-3396um-1180um-20160220 DT2 3.78 {'hpc': 64, 'ls': 64} 104 {'hpc': 82, 'ls': 22} True True 71
17 sub-DT2_ses-DT2-rPPC-rCCG-3396um-1180um-20160221 DT2 5.87 {'hpc': 64, 'ls': 64} 133 {'hpc': 112, 'ls': 21} True True 121
18 sub-DT2_ses-DT2-rPPC-rCCG-3468um-1216um-20160222 DT2 4.26 {'hpc': 64, 'ls': 64} 205 {'hpc': 154, 'ls': 51} True True 100
19 sub-DT2_ses-DT2-rPPC-rCCG-3468um-1216um-20160223 DT2 2.49 {'hpc': 64, 'ls': 64} 193 {'hpc': 148, 'ls': 45} True True 115
20 sub-DT2_ses-DT2-rPPC-rCCG-3468um-1216um-20160225 DT2 4.76 {'hpc': 64, 'ls': 64} 200 {'hpc': 159, 'ls': 41} True True 113
21 sub-DT2_ses-DT2-rPPC-rCCG-3468um-1252um-20160226 DT2 6.57 {'hpc': 64, 'ls': 64} 142 {'hpc': 99, 'ls': 43} True True 134
22 sub-DT2_ses-DT2-rPPC-rCCG-3540um-1288um-20160227 DT2 4.67 {'hpc': 64, 'ls': 64} 179 {'hpc': 126, 'ls': 53} True True 207
23 sub-DT2_ses-DT2-rPPC-rCCG-3540um-1288um-20160228 DT2 3.15 {'hpc': 64, 'ls': 64} 113 {'hpc': 79, 'ls': 34} True True 106
24 sub-DT2_ses-DT2-rPPC-rCCG-3612um-1360um-20160304 DT2 0.46 {'hpc': 64, 'ls': 64} 25 {'hpc': 9, 'ls': 16} True True 30
25 sub-DT5_ses-20161007-720um-720um-merge DT5 4.88 {'hpc': 64, 'ls': 32} 13 no location column False False 0
26 sub-DT5_ses-20161008-792um-864um-merge DT5 5.20 {'hpc': 64, 'ls': 32} 31 {'hpc': 6, 'ls': 25} False False 0
27 sub-DT5_ses-20161011-972um-1224um-merge DT5 9.18 {'hpc': 64, 'ls': 32} 52 {'hpc': 27, 'ls': 25} True True 186
28 sub-DT5_ses-20161014-1296um-1476um-merge DT5 5.56 {'hpc': 64, 'ls': 32} 82 {'hpc': 48, 'ls': 34} True True 151
29 sub-DT5_ses-20161015-1368um-1476um-merge DT5 4.84 {'hpc': 64, 'ls': 32} 67 {'hpc': 27, 'ls': 40} True True 115
30 sub-DT5_ses-20161016-1368um-1476um-merge DT5 0.89 {'hpc': 64, 'ls': 32} 32 {'hpc': 13, 'ls': 19} True True 33
31 sub-DT5_ses-20161017-1440um-2016um-merge DT5 6.33 {'hpc': 64, 'ls': 32} 83 {'hpc': 49, 'ls': 34} False False 0
32 sub-DT5_ses-20161018-1552um-2304um-merge DT5 8.53 {'ca3': 64, 'ls': 32} 100 {'ca3': 56, 'ls': 44} True True 213
33 sub-DT5_ses-20161019-1624um-2340um-161019-095311 DT5 6.52 {'ca3': 64, 'ls': 32} 97 {'ca3': 49, 'ls': 48} True True 238
34 sub-DT5_ses-20161020-1768um-2376um-161020-102642 DT5 7.33 {'hpc': 64, 'ls': 32} 102 {'hpc': 62, 'ls': 40} True True 245
35 sub-DT5_ses-20161021-1840um-2412um-merge DT5 7.95 {'ca3': 64, 'ls': 32} 93 {'ca3': 40, 'ls': 53} True True 54
36 sub-DT5_ses-20161023-1872um-2448um-161023-123504 DT5 1.25 {'hpc': 64, 'ls': 32} 33 {'hpc': 8, 'ls': 25} True True 138
37 sub-DT5_ses-20161025-1944um-2232um-merge DT5 7.37 {'hpc': 64, 'ls': 32} 68 {'hpc': 40, 'ls': 28} True True 117
38 sub-DT5_ses-20161026-1944um-2232um-merge DT5 7.68 {'hpc': 64, 'ls': 32} 67 {'hpc': 36, 'ls': 31} True True 57
39 sub-DT5_ses-20161027-1944um-2232um-161027-194125 DT5 0.63 {'hpc': 64, 'ls': 32} 39 {'hpc': 24, 'ls': 15} True True 190
40 sub-DT5_ses-20161029-2016um-2232um-merge DT5 5.25 {'hpc': 64, 'ls': 32} 68 {'hpc': 47, 'ls': 21} False False 0
41 sub-DT5_ses-20161030-2052um-2232um-161030-111142 DT5 6.92 {'hpc': 64, 'ls': 32} 56 {'hpc': 44, 'ls': 12} True True 225
42 sub-DT5_ses-20161031-2088um-2232um-161031-104001 DT5 8.71 {'ca3': 64, 'ls': 32} 70 {'ca3': 46, 'ls': 24} True True 185
43 sub-DT5_ses-20161101-2088um-2232um-merge DT5 6.41 {'ca3': 64, 'ls': 32} 69 {'ca3': 45, 'ls': 24} True True 76
44 sub-DT5_ses-20161104-2124um-2232um-merge DT5 6.68 {'ca3': 64, 'ls': 32} 61 {'ca3': 43, 'ls': 18} True True 150
45 sub-DT7_ses-20170314-72um-0um-170314-181434 DT7 1.55 {'hpc': 64, 'ls': 64} 13 no location column False False 0
46 sub-DT7_ses-20170315-216um-0um-merge DT7 7.22 {'hpc': 64, 'ls': 64} 0 no location column False False 0
47 sub-DT7_ses-20170315-288um-0um-170315-191240 DT7 0.60 {'acg': 23, 'cc': 18, 'hpc': 64, 'ls': 23} 61 {'acg': 25, 'cc': 5, 'hpc': 22, 'ls': 9} True True 123
48 sub-DT7_ses-20170316-288um-0um-merge DT7 6.76 {'hpc': 64, 'ls': 64} 45 no location column False False 0
49 sub-DT7_ses-20170317-288um-0um-merge DT7 6.87 {'hpc': 64, 'ls': 64} 43 no location column True True 76
50 sub-DT7_ses-20170318-288um-0um-170318-115923 DT7 7.97 {'acg': 38, 'hpc': 63, 'ls': 24, 'ms': 3} 40 {'hpc': 40} True True 252
51 sub-DT7_ses-20170319-288um-0um-170319-123708 DT7 7.00 {'hpc': 64, 'ls': 64} 27 no location column False False 0
52 sub-DT7_ses-20170320-432um-0um-merge DT7 8.32 {'acg': 34, 'hpc': 64, 'ls': 23, 'ms': 7} 70 {'hpc': 70} True True 326
53 sub-DT7_ses-20170321-216um-0um-170321-120259 DT7 10.01 {'acg': 34, 'hpc': 63, 'ls': 24, 'ms': 7} 100 {'hpc': 100} True True 249
54 sub-DT7_ses-20170324-576um-144um-170324-123932 DT7 2.73 {'acg': 31, 'hpc': 63, 'ls': 24, 'ms': 10} 0 no location column True True 121
55 sub-DT7_ses-20170325-576um-360um-170325-125511 DT7 8.28 {'acg': 29, 'hpc': 64, 'ls': 21, 'ms': 14} 81 {'acg': 6, 'hpc': 53, 'ls': 12, 'ms': 10} True True 239
56 sub-DT7_ses-20170327-864um-360um-170327-115803 DT7 9.35 {'unknown': 128} 0 no location column False False 0
57 sub-DT7_ses-20170401-864um-936um-merge DT7 4.05 {'cc': 20, 'hpc': 63, 'ls': 24, 'ms': 21} 47 {'cc': 1, 'hpc': 22, 'ls': 9, 'ms': 15} True True 81
58 sub-DT7_ses-20170405-864um-1152um-170405-102118 DT7 10.63 {'hpc': 64, 'ls': 64} 0 no location column False True 51
59 sub-DT7_ses-20170406-864um-1152um-merge DT7 10.92 {'unknown': 128} 0 no location column True True 27
60 sub-DT7_ses-20170407-936um-1152um-merge DT7 8.06 {'unknown': 128} 36 no location column False True 54
61 sub-DT7_ses-20170408-936um-1152um-merge DT7 4.18 {'unknown': 128} 0 no location column False False 0
62 sub-DT7_ses-20170409-1152um-1152um-170409-115803 DT7 9.13 {'hpc': 64, 'ls': 64} 0 no location column True True 31
63 sub-DT7_ses-20170410-1152um-1152um-170410-102443 DT7 8.23 {'hpc': 64, 'ls': 64} 61 no location column True True 82
64 sub-DT7_ses-20170412-1440um-1170um-170412-103153 DT7 7.55 {'hpc': 63, 'ls': 17, 'sFi': 48} 223 {'hpc': 139, 'ls': 30, 'sFi': 54} True True 76
65 sub-DT8_ses-20170125-0um-0um-merge DT8 6.50 {'hpc': 64, 'ls': 32} 28 {'ls': 28} False True 17
66 sub-DT8_ses-20170128-0um-72um-merge DT8 5.49 {'hpc': 64, 'ls': 32} 27 {'ls': 27} False True 71
67 sub-DT8_ses-20170201-72um-432um-merge DT8 3.90 {'hpc': 64, 'ls': 32} 45 {'ls': 13, 'ppc': 32} False True 137
68 sub-DT8_ses-20170202-72um-576um-170202-092033 DT8 5.48 {'hpc': 64, 'ls': 32} 11 {'ls': 11} False True 116
69 sub-DT8_ses-20170206-72um-1872um-merge DT8 9.07 {'hpc': 64, 'ls': 32} 15 {'ls': 15} False True 146
70 sub-DT8_ses-20170207-144um-1944um-merge DT8 7.82 {'hpc': 64, 'ls': 32} 19 {'ls': 19} False True 174
71 sub-DT8_ses-20170211-144um-1944um-merge DT8 5.50 {'hpc': 64, 'ls': 32} 20 {'ls': 20} False True 74
72 sub-DT8_ses-20170214-216um-1944um-170214-102804 DT8 5.56 {'hpc': 64, 'ls': 32} 8 {'ls': 8} False True 242
73 sub-DT8_ses-20170215-216um-1944um-170215-140739 DT8 6.24 {'hpc': 64, 'ls': 32} 21 {'ls': 21} False True 72
74 sub-DT8_ses-20170216-216um-1944um-merge DT8 6.05 {'hpc': 64, 'ls': 32} 22 {'ls': 22} False True 264
75 sub-DT8_ses-20170217-216um-1944um-merge DT8 5.80 {'hpc': 64, 'ls': 32} 17 {'ls': 17} False True 226
76 sub-DT8_ses-20170218-216um-1944um-170218-114724 DT8 5.05 {'hpc': 64, 'ls': 32} 16 {'ls': 16} False True 254
77 sub-DT8_ses-20170221-324um-1944um-170221-101929 DT8 5.57 {'hpc': 64, 'ls': 32} 18 {'ls': 18} False True 199
78 sub-DT8_ses-20170224-324um-2088um-170224-101710 DT8 8.74 {'hpc': 64, 'ls': 32} 20 {'ls': 20} False True 210
79 sub-DT8_ses-20170305-468um-2088um-170305-135233 DT8 4.59 {'hpc': 64, 'ls': 32} 15 {'ls': 15} False True 209
80 sub-DT8_ses-20170306-468um-2088um-170306-113628 DT8 8.03 {'hpc': 64, 'ls': 32} 11 {'ls': 11} False True 234
81 sub-DT8_ses-20170308-612um-2088um-merge DT8 8.30 {'hpc': 64, 'ls': 32} 15 {'ls': 15} False True 92
82 sub-DT9_ses-20170428-0um-0um-merge DT9 7.47 {'hpc': 64, 'ls': 32} 10 no location column False False 0
83 sub-DT9_ses-20170429-0um-0um-merge DT9 7.16 {'hpc': 64, 'ls': 32} 148 {'hpc': 134, 'ls': 14} True True 160
84 sub-DT9_ses-20170430-0um-0um-merge DT9 1.40 {'hpc': 64, 'ls': 32} 62 {'hpc': 46, 'ls': 16} True True 159
85 sub-DT9_ses-20170501-216um-0um-170501-112210 DT9 9.02 {'hpc': 64, 'ls': 32} 140 {'hpc': 90, 'ls': 50} True True 175
86 sub-DT9_ses-20170502-288um-0um-merge DT9 9.76 {'hpc': 64, 'ls': 32} 141 {'hpc': 89, 'ls': 52} True True 60
87 sub-DT9_ses-20170504-396um-0um-merge DT9 9.72 {'hpc': 64, 'ls': 32} 153 {'hpc': 101, 'ls': 52} True True 56
88 sub-DT9_ses-20170505-396um-0um-merge DT9 7.78 {'hpc': 64, 'ls': 32} 127 {'hpc': 81, 'ls': 46} True True 221
89 sub-DT9_ses-20170506-468um-36um-170506-111835 DT9 7.35 {'hpc': 64, 'ls': 32} 121 {'hpc': 80, 'ls': 41} True True 293
90 sub-DT9_ses-20170507-468um-36um-merge DT9 9.18 {'hpc': 64, 'ls': 32} 131 {'hpc': 88, 'ls': 43} True True 62
91 sub-DT9_ses-20170508-468um-36um-170508-102251 DT9 10.13 {'hpc': 64, 'ls': 32} 132 {'hpc': 83, 'ls': 49} True True 76
92 sub-DT9_ses-20170509-468um-36um-170509-103451 DT9 7.34 {'hpc': 64, 'ls': 32} 113 {'hpc': 74, 'ls': 39} True True 263
93 sub-DT9_ses-20170518-828um-864um-merge DT9 9.11 {'ca3': 64, 'ls': 32} 95 {'ca3': 62, 'ls': 33} True True 315
94 sub-DT9_ses-20170519-864um-900um-merge DT9 6.72 {'ca3': 64, 'ls': 32} 107 {'ca3': 76, 'ls': 31} True True 85
95 sub-DT9_ses-20170521-900um-936um-170521-125602 DT9 4.59 {'ca3': 64, 'ls': 32} 120 {'ca3': 92, 'ls': 28} True True 93
96 sub-DT9_ses-20170525-900um-936um-170525-150124 DT9 1.99 {'ca3': 64, 'ls': 32} 73 {'ca3': 57, 'ls': 16} True True 57
97 sub-DT9_ses-20170526-1116um-1000um-merge DT9 6.69 {'hpc': 64, 'ls': 32} 20 no location column False True 92
files per subject: {'DT1': 12, 'DT2': 13, 'DT5': 20, 'DT7': 20, 'DT8': 17, 'DT9': 16}
units by region label: {np.str_('ls'): 2448, '(no location)': 268, np.str_('ppc'): 60, np.str_('hpc'): 2778, np.str_('ca3'): 566, np.str_('acg'): 31, np.str_('cc'): 6, np.str_('ms'): 25, np.str_('sFi'): 54}
files whose units have no location column: ['sub-DT5_ses-20161007-720um-720um-merge_e', 'sub-DT7_ses-20170314-72um-0um-170314-181', 'sub-DT7_ses-20170316-288um-0um-merge_ece', 'sub-DT7_ses-20170317-288um-0um-merge_beh', 'sub-DT7_ses-20170319-288um-0um-170319-12', 'sub-DT7_ses-20170407-936um-1152um-merge_', 'sub-DT7_ses-20170410-1152um-1152um-17041', 'sub-DT9_ses-20170428-0um-0um-merge_eceph', 'sub-DT9_ses-20170526-1116um-1000um-merge']
LFP rates (Hz): [1250.0] | conversion: [1.9499999999999999e-07] ['volts']
electrode x, y, z all NaN in every file: True
LFP starting_time values (s): [0.0]
units-table columns: [(), ('clu_id', 'group_id', 'id', 'location', 'spike_times', 'spike_times_index'), ('id', 'spike_times', 'spike_times_index')]
first unit of the first file with a group_id column, by animal:
DT1 {'clu_id': '3', 'group_id': 'Group1', 'id': '1', 'location': 'ls'}
DT2 {'clu_id': '71', 'group_id': 'Group1', 'id': '1', 'location': 'ls'}
DT5 {'clu_id': '20', 'group_id': 'Group2', 'id': '1', 'location': 'ls'}
DT7 {'clu_id': '18', 'group_id': 'Group1', 'id': '1', 'location': 'hpc'}
DT8 {'clu_id': '2', 'group_id': 'Group1', 'id': '1', 'location': 'ls'}
DT9 {'clu_id': '36', 'group_id': 'Group1', 'id': '1', 'location': 'ls'}
every unit group_id is an electrode group of its file: True
electrode group descriptions: ['electrodes.'] | e.g. ['Group1 electrodes.', 'Group10 electrodes.']
devices (name: description), over all files:
Cambridge prob (1 x 64) : Silicon probe from Cambridge Neurotech. Electrophysiological data were acquired using an Intan RHD2000 system (Intan Technologies LLC) digitized with20 kHz rate.
Device_ecephys : Ecephys probe. Automatically generated.
Neuronexus probe (4 x 1) : according to author thse are reference sites a few millimeters dorsal from the restrecorded from an older neuronexus probe
Neuronexus probe (4 x 8) : A 4 (shanks) x 8 (electrodes) silicon probe from Neuronexus. Electrophysiological data were acquired using an Intan RHD2000 system (Intan Technologies LLC) digitized with 20 kHz rate.
Neuronexus probe (5 x 12) : A 5 (shanks) x 12 (electrodes) silicon probe from Neuronexus. Electrophysiological data were acquired using an Intan RHD2000 system (Intan Technologies LLC) digitized with 20 kHz rate.
Neuronexus probe (6 x 10) : A 6 (shanks) x 10 (electrodes) silicon probe from Neuronexus. Electrophysiological data were acquired using an Intan RHD2000 system (Intan Technologies LLC) digitized with 20 kHz rate.
Neuronexus probe (8 x 8) : A 8 (shanks) x 8 (electrodes) silicon probe from Neuronexus. Electrophysiological data were acquired using an Intan RHD2000 system (Intan Technologies LLC) digitized with 20 kHz rate.
The archive's summary gives 6 subjects, which matches the 6 rats (DT1, DT2, DT5, DT7, DT8, DT9) and the paper's 5 rats of Tingley and Buzsaki (2018) plus one tetrode rat. DT1 holds LS channels only, without sleep states, epochs or trials, which fits the tetrode rat used only in the phase-locking analysis; the files do not say which rat is which.
Unit region labels: ls, hpc, ca3, and a few cortical, septal and fimbria labels. The hpc
label does not separate CA1 from other hippocampal sites. Nine files have units without a location
column; their 268 units cannot be assigned to a region and are left out of every unit analysis.
The LFP is stored at 1250 Hz as int16 with a conversion of 1.95e-7 V (0.195 uV) per step. The electrode x, y and z are NaN in every file, so the probe geometry
and its orientation in the brain are not in the files.
Where the units table has a location column it holds clu_id, group_id and location besides
the spike times; the nine files without it hold only the ids and spike times. No file has unit
waveforms, a peak channel or a quality metric. group_id names an electrode group, and each group is one shank of
a probe (the electrode groups are described only as "GroupN electrodes."), so it is used as the
shank. The device names give the probe layouts (shanks by sites per shank), without site pitch or
the orientation of the probe in the brain. One device, "Neuronexus probe (4 x 1)", is described as
"reference sites a few millimeters dorsal" to the rest. The LFP starts at 0 s in every file, the
same origin as the spikes, states and trials.
In [4]:
# Duplicates, spike-time format, states, epochs and trials.
lengths = collections.Counter(r["meta"]["lfp_shape"][0] for r in results)
fps = collections.Counter(r["fingerprint"] for r in R)
print("files sharing an LFP length:", sum(v for v in lengths.values() if v > 1),
"| files sharing the checksum of the first 10 s of their LS channel:", sum(v for v in fps.values() if v > 1))
spike_counts = collections.Counter(tuple(len(s) for s in r["meta"]["units"]["spikes"]) for r in results if len(r["meta"]["units"]["spikes"]))
print("files sharing the spike counts of all their units:", sum(v for v in spike_counts.values() if v > 1))
bad = []
for r in results:
m = r["meta"]
for i, s in enumerate(m["units"]["spikes"]):
if len(s) and (np.any(np.diff(s) < 0) or s.min() < 0 or s.max() > m["duration"] + 1):
bad.append({"file": r["path"].split("/")[1][:40], "unit": i,
"region": None if m["units"]["location"] is None else m["units"]["location"][i],
"spikes": len(s), "out_of_order": int((np.diff(s) < 0).sum()),
"first": round(float(s.min()), 3), "last_minus_duration": round(float(s.max() - m["duration"]), 3)})
print(f"units with spike times out of order or outside the recording: {len(bad)}")
print(pd.DataFrame(bad).to_string())
oo_files = sorted({r["path"] for r in results for s in r["meta"]["units"]["spikes"] if len(s) and np.any(np.diff(s) < 0)})
print(f"units with spike times out of order: {sum(b['out_of_order'] > 0 for b in bad)}, in {len(oo_files)} files:")
for f_ in oo_files:
print(" ", f_)
labels = collections.Counter()
cover, epoch_ok, trial_dur = [], [], []
for r in results:
m = r["meta"]
if m["states"] is not None:
labels.update(m["states"]["label"])
cover.append((m["states"]["stop"] - m["states"]["start"]).sum() / m["duration"])
if m["epochs"] is not None and m["trials"] is not None:
ep, tr = m["epochs"], m["trials"]
epoch_ok.append(abs(ep["stop"][0] - tr["start"].min()) < 1 and abs(ep["stop"][1] - tr["stop"].max()) < 1)
trial_dur.append(tr["stop"] - tr["start"])
print("sleep-state labels (intervals):", dict(labels))
print(f"fraction of each recording covered by states: {np.min(cover):.2f} to {np.max(cover):.2f} over {len(cover)} files")
print("epoch tags:", sorted({tuple(r['meta']['epochs']['tags']) for r in results if r['meta']['epochs'] is not None}))
print(f"files where the BEHAV epoch runs from the first trial start to the last trial stop: {sum(epoch_ok)} of {len(epoch_ok)}")
td = np.concatenate(trial_dur)
print(f"trial durations: median {np.median(td):.1f} s, range {td.min():.2f} to {td.max():.1f} s, {len(td)} trials")
# Position data (for fig7e-replay): which files hold behavior series, and do they cover the trials?
rows = []
for r in results:
m = r["meta"]
pos = {k: v for k, v in m["behavior"].items() if k.endswith("position")}
if not pos and m["trials"] is None:
continue
k, v = next(iter(pos.items())) if pos else (None, None)
tr = m["trials"]
cov = np.nan
if v is not None and tr is not None:
cov = np.mean((tr["start"] >= v["t0"] - 1) & (tr["stop"] <= v["t1"] + 1))
rows.append({"subject": m["subject"], "position series": k, "samples": None if v is None else v["shape"][0],
"trials": 0 if tr is None else len(tr["start"]), "trials inside the position time": cov})
bp = pd.DataFrame(rows)
print("files with trials or a position series, by animal:")
print(bp.groupby("subject").agg(files=("trials", "size"), with_position=("position series", lambda x: int(x.notna().sum())),
trials=("trials", "sum"), mean_coverage=("trials inside the position time", "mean")).round(3).to_string())
print("position series names:", dict(collections.Counter(bp["position series"].dropna())))
files sharing an LFP length: 0 | files sharing the checksum of the first 10 s of their LS channel: 0 files sharing the spike counts of all their units: 0
units with spike times out of order or outside the recording: 25
file unit region spikes out_of_order first last_minus_duration
0 sub-DT2_ses-DT2-rPPC-rCCG-3324um-1144um- 3 ls 779907 1 0.000 -0.038
1 sub-DT2_ses-DT2-rPPC-rCCG-3324um-1144um- 16 ls 589905 1 0.000 -0.090
2 sub-DT2_ses-DT2-rPPC-rCCG-3324um-1144um- 29 hpc 392003 4 0.000 -0.210
3 sub-DT2_ses-DT2-rPPC-rCCG-3324um-1144um- 36 hpc 607301 7 0.000 -0.032
4 sub-DT2_ses-DT2-rPPC-rCCG-3324um-1144um- 58 hpc 782899 9 0.000 -0.001
5 sub-DT2_ses-DT2-rPPC-rCCG-3396um-1180um- 8 ls 582047 1 0.000 -0.054
6 sub-DT2_ses-DT2-rPPC-rCCG-3396um-1180um- 43 hpc 1318374 26 0.000 -0.002
7 sub-DT5_ses-20161007-720um-720um-merge_e 2 NaN 431966 2 0.000 -0.052
8 sub-DT5_ses-20161007-720um-720um-merge_e 11 NaN 512420 305 0.000 -0.005
9 sub-DT5_ses-20161017-1440um-2016um-merge 20 ls 1483149 16 0.000 -0.029
10 sub-DT5_ses-20161017-1440um-2016um-merge 33 ls 1503446 2782 0.000 -0.037
11 sub-DT5_ses-20161017-1440um-2016um-merge 41 hpc 1011538 6 0.000 -0.001
12 sub-DT5_ses-20161017-1440um-2016um-merge 49 hpc 696444 20 0.000 -0.000
13 sub-DT5_ses-20161017-1440um-2016um-merge 65 hpc 562440 16 0.000 -0.006
14 sub-DT5_ses-20161017-1440um-2016um-merge 82 hpc 1148051 2 0.000 -0.009
15 sub-DT5_ses-20161018-1552um-2304um-merge 30 ls 69323 1 0.724 -0.088
16 sub-DT5_ses-20161026-1944um-2232um-merge 31 hpc 550179 5 0.000 -0.004
17 sub-DT5_ses-20161029-2016um-2232um-merge 2 ls 1110215 11 0.000 -0.006
18 sub-DT5_ses-20161029-2016um-2232um-merge 14 ls 862447 1282 0.000 -0.044
19 sub-DT5_ses-20161029-2016um-2232um-merge 20 ls 940337 231 0.000 -0.038
20 sub-DT5_ses-20161029-2016um-2232um-merge 34 hpc 481326 11 0.000 -0.010
21 sub-DT5_ses-20161029-2016um-2232um-merge 60 hpc 840710 52 0.000 -0.029
22 sub-DT7_ses-20170314-72um-0um-170314-181 12 NaN 201038 1 0.000 -0.012
23 sub-DT7_ses-20170321-216um-0um-170321-12 3 hpc 184332 1 0.373 -0.042
24 sub-DT8_ses-20170201-72um-432um-merge_be 23 ppc 131838 9 0.144 -6.695
units with spike times out of order: 25, in 10 files:
sub-DT2/sub-DT2_ses-DT2-rPPC-rCCG-3324um-1144um-20160219-160219-130435_behavior+ecephys.nwb
sub-DT2/sub-DT2_ses-DT2-rPPC-rCCG-3396um-1180um-20160220-160220-141839_behavior+ecephys.nwb
sub-DT5/sub-DT5_ses-20161007-720um-720um-merge_ecephys.nwb
sub-DT5/sub-DT5_ses-20161017-1440um-2016um-merge_ecephys.nwb
sub-DT5/sub-DT5_ses-20161018-1552um-2304um-merge_behavior+ecephys.nwb
sub-DT5/sub-DT5_ses-20161026-1944um-2232um-merge_behavior+ecephys.nwb
sub-DT5/sub-DT5_ses-20161029-2016um-2232um-merge_ecephys.nwb
sub-DT7/sub-DT7_ses-20170314-72um-0um-170314-181434_ecephys.nwb
sub-DT7/sub-DT7_ses-20170321-216um-0um-170321-120259_behavior+ecephys.nwb
sub-DT8/sub-DT8_ses-20170201-72um-432um-merge_behavior+ecephys.nwb
sleep-state labels (intervals): {np.str_('Awake'): 1052, np.str_('Non-REM'): 4835, np.str_('MA'): 2043, np.str_('REM'): 940}
fraction of each recording covered by states: 0.82 to 1.00 over 54 files
epoch tags: [('before trials', 'during trials', 'after trials')]
files where the BEHAV epoch runs from the first trial start to the last trial stop: 74 of 74
trial durations: median 3.8 s, range 0.23 to 11.5 s, 10292 trials
files with trials or a position series, by animal:
files with_position trials mean_coverage
subject
DT2 12 12 1227 1.000
DT5 16 16 2373 1.000
DT7 14 14 1788 0.929
DT8 17 17 2737 1.000
DT9 15 15 2167 1.000
position series names: {'wheel_alternation_task/position': 21, 'central_alternation_task/position': 26, 'linear_task/position': 9, 'generic_position_tracking_task/position': 6, 'both_alternation_task/position': 11, 'linear-jump_task/position': 1}
No recording is deposited twice: no two files share an LFP length, the checksum of the first 10 s of their LS channel, or the spike counts of all their units.
25 units in 10 files have spike times out of order (1 to 2782 inversions each), a break of the NWB rule that
spike times are sorted. Spike times are sorted here before use. One ppc unit in DT8 ends 6.7 s
before the recording; all spikes fall inside the recordings.
Sleep states cover 82 to 100% of each recording that has them, with the labels Non-REM, Awake,
MA (microarousal) and REM. Epochs are named before trials, during trials and
after trials, and in every file the middle epoch runs from the first trial's start to the last
trial's stop, so the epochs are derived from the trials table. They are used as PRE, BEHAV and POST.
In [5]:
# Channels chosen, events at each detection step, and unit firing rates.
rows = []
for r in R:
m = r["meta"]
row = {"file": r["path"].split("/")[1][:40], "subject": m["subject"], "ls_ch": r["channels"]["ls"],
"hpc_ch": r["channels"]["hpc"], "emg_shanks": len(r["channels"]["emg"]),
"bad_channels": int((~r["channel_stats"]["good"]).sum()),
"emg>0.9 (%)": None if r["emg_summary"] is None else round(100 * r["emg_summary"]["frac_above"], 2)}
for k, v in r["ls_events"]["counts"].items():
row[f"ls_{k}"] = v
if "hpc_events" in r:
row["hpc_final"] = r["hpc_events"]["counts"]["emg"]
rows.append(row)
det = pd.DataFrame(rows)
print(det.groupby("subject")[[c for c in det.columns if c.startswith(("ls_t", "ls_m", "ls_p", "ls_d", "ls_e", "hpc_final"))]].sum().to_string())
print(f"LS HFOs over all {len(det)} recordings with LS electrodes: {det.ls_emg.sum()}; "
f"removed by the EMG check: {(det.ls_duration - det.ls_emg).sum()}")
print(f"recordings with a hippocampal (`hpc`) channel: {det.hpc_ch.notna().sum()}, CA1 ripples {int(det.hpc_final.sum())}")
print(f"channels dropped as bad (broadband SD outside 0.2 to 5 times the region median): {det.bad_channels.sum()} in {(det.bad_channels > 0).sum()} files")
print("EMG above 0.9, % of 0.1 s windows: median", det["emg>0.9 (%)"].median(), "max", det["emg>0.9 (%)"].max())
rates = []
for r in results:
m = r["meta"]
loc = m["units"]["location"]
for i, s in enumerate(m["units"]["spikes"]):
rates.append({"region": "(no location)" if loc is None else loc[i], "rate": len(s) / m["duration"]})
rates = pd.DataFrame(rates)
print("unit firing rates (Hz) by region:")
print(rates.groupby("region").rate.describe(percentiles=[0.5, 0.9]).round(2).to_string())
print("units above 50 Hz:", rates.groupby("region").rate.apply(lambda x: int((x > 50).sum())).to_dict())
ls_threshold ls_merged ls_peak ls_duration ls_emg hpc_final
subject
DT1 69170 44415 12371 11089 5292 0.0
DT2 15436 12111 4003 3501 1183 26531.0
DT5 31648 25279 8020 6991 4032 29858.0
DT7 476660 297919 48236 41570 38008 95114.0
DT8 25211 18640 5680 5063 2757 1854.0
DT9 126641 84954 20293 18656 15914 76804.0
LS HFOs over all 94 recordings with LS electrodes: 67186; removed by the EMG check: 19684
recordings with a hippocampal (`hpc`) channel: 72, CA1 ripples 230161
channels dropped as bad (broadband SD outside 0.2 to 5 times the region median): 5 in 3 files
EMG above 0.9, % of 0.1 s windows: median 0.55 max 17.13
unit firing rates (Hz) by region:
count mean std min 50% 90% max
region
(no location) 268.0 1.42 6.02 0.00 0.12 1.30 61.88
acg 31.0 2.11 2.70 0.00 0.79 6.10 10.00
ca3 566.0 1.19 3.88 0.00 0.18 1.73 37.80
cc 6.0 0.73 0.43 0.04 0.76 1.14 1.17
hpc 2778.0 1.78 5.60 0.00 0.21 3.29 96.88
ls 2448.0 1.64 4.00 0.00 0.63 3.56 65.99
ms 25.0 1.25 3.88 0.05 0.24 1.70 19.66
ppc 60.0 1.03 3.06 0.01 0.18 1.43 20.41
sFi 54.0 0.91 1.67 0.01 0.22 2.47 8.54
units above 50 Hz: {'(no location)': 1, np.str_('acg'): 0, np.str_('ca3'): 0, np.str_('cc'): 0, np.str_('hpc'): 5, np.str_('ls'): 5, np.str_('ms'): 0, np.str_('ppc'): 0, np.str_('sFi'): 0}
Detection counts are summed by animal for each step of the detector: threshold crossings, after merging, after the 5 SD peak test, after the duration limits, and after the EMG check. The EMG check removes 19,684 LS events (23% of those that pass the other steps). LS HFO counts differ widely between animals; DT7 contributes more than half of them (38,008 of 67,186). Only 5 channels in 3 files were dropped as bad.
Unit firing rates are low (median 0.18 to 0.63 Hz for the hpc, ca3 and ls units). 11 units
fire above 50 Hz on average; the files carry no isolation quality or single-unit label, so all
units are kept.
In [6]:
# Do the spikes, LFP, sleep states and trials share one clock?
# (1) Spikes against LFP: locking of spikes to the 120-200 Hz phase with the spikes shifted in time.
for key, name in (("clock_hpc", "hippocampal units, inside CA1 ripples"), ("clock_ls", "LS units, inside high-power epochs")):
curves = [r[key] for r in R if r.get(key) is not None and np.isfinite(r[key]["r"]).all()]
shifts = np.array(curves[0]["shifts_ms"])
M = np.array([c["r"] for c in curves])
best = shifts[M.argmax(axis=1)]
print(f"{name}: {len(curves)} recordings; median resultant by shift (ms):",
dict(zip(shifts.tolist(), np.round(np.median(M, axis=0), 3).tolist())))
print(f" recordings whose locking peaks within +-2 ms of zero shift: {np.sum(np.abs(best) <= 2)} of {len(best)};"
f" peaks elsewhere: {sorted(collections.Counter(best[np.abs(best) > 2].tolist()).items())}")
names = [r["path"].split("/")[1][:34] for r in R if r.get(key) is not None and np.isfinite(r[key]["r"]).all()]
print(" off zero:", [(n, int(b)) for n, b in zip(names, best) if abs(b) > 2])
# (2) Sleep states against LFP: CA1 ripples should fall in NREM far more than in WAKE.
def rate_in(times, states, label, shift=0.0):
sel = states["label"] == label
s0, s1 = states["start"][sel] + shift, states["stop"][sel] + shift
order = np.argsort(s0)
n = ls_hfo.counts_in(np.sort(times), s0[order], s1[order]).sum()
return n / (s1 - s0).sum()
rows = []
for r in R:
if "hpc_events" in r and r["meta"]["states"] is not None:
st, t = r["meta"]["states"], r["hpc_events"]["peak"]
row = {"file": r["path"].split("/")[1][:30]}
for sh in (-600, 0, 600):
row[sh] = rate_in(t, st, "Non-REM", sh) / max(rate_in(t, st, "Awake", sh), 1e-9)
rows.append(row)
cl = pd.DataFrame(rows)
print(f"\nCA1 ripple rate in NREM / WAKE over {len(cl)} recordings, median: true labels {cl[0].median():.2f}, "
f"labels shifted -10 min {cl[-600].median():.2f}, +10 min {cl[600].median():.2f}; "
f"true labels give the largest ratio in {(cl[0] >= cl[[-600, 600]].max(axis=1)).sum()} of {len(cl)}")
# (3) Trials against states: the animal is awake on the maze, so trial time should not be scored NREM.
rows = []
for r in R:
m = r["meta"]
if m["states"] is not None and m["trials"] is not None:
st = m["states"]
nrem = st["label"] == "Non-REM"
row = {}
for sh in (-3600, 0, 3600):
t = np.concatenate([np.arange(a, b, 0.1) for a, b in zip(m["trials"]["start"] + sh, m["trials"]["stop"] + sh)])
t = t[(t > 0) & (t < m["duration"])]
row[sh] = np.mean(ls_hfo.state_of(t, st) == "Non-REM") if len(t) else np.nan
rows.append(row)
tc = pd.DataFrame(rows)
print(f"fraction of trial time scored NREM over {len(tc)} recordings, mean: true {tc[0].mean():.3f}, "
f"trials shifted -1 h {tc[-3600].mean():.3f}, +1 h {tc[3600].mean():.3f}")
hippocampal units, inside CA1 ripples: 39 recordings; median resultant by shift (ms): {-20: 0.018, -10: 0.083, -5: 0.116, -2: 0.14, 0: 0.145, 2: 0.141, 5: 0.111, 10: 0.069, 20: 0.015}
recordings whose locking peaks within +-2 ms of zero shift: 32 of 39; peaks elsewhere: [(-20, 1), (-10, 2), (-5, 1), (5, 1), (10, 2)]
off zero: [('sub-DT5_ses-20161016-1368um-1476um', 10), ('sub-DT5_ses-20161020-1768um-2376um', 10), ('sub-DT5_ses-20161025-1944um-2232um', 5), ('sub-DT5_ses-20161026-1944um-2232um', -10), ('sub-DT5_ses-20161027-1944um-2232um', -20), ('sub-DT5_ses-20161029-2016um-2232um', -5), ('sub-DT5_ses-20161030-2052um-2232um', -10)]
LS units, inside high-power epochs: 78 recordings; median resultant by shift (ms): {-20: 0.018, -10: 0.036, -5: 0.087, -2: 0.123, 0: 0.156, 2: 0.136, 5: 0.096, 10: 0.037, 20: 0.019}
recordings whose locking peaks within +-2 ms of zero shift: 75 of 78; peaks elsewhere: [(5, 1), (10, 2)]
off zero: [('sub-DT2_ses-DT2-rPPC-rCCG-3036um-3', 10), ('sub-DT2_ses-DT2-rPPC-rCCG-3612um-1', 5), ('sub-DT9_ses-20170525-900um-936um-1', 10)]
CA1 ripple rate in NREM / WAKE over 43 recordings, median: true labels 8.98, labels shifted -10 min 1.90, +10 min 1.58; true labels give the largest ratio in 36 of 43
fraction of trial time scored NREM over 53 recordings, mean: true 0.017, trials shifted -1 h 0.306, +1 h 0.385
Three checks that the streams share one clock, each with a contrast that an offset would destroy.
- Spikes and LFP: the phase of the 120-200 Hz band at each spike is most consistent at zero shift and falls off within 5 to 10 ms, about one cycle. This holds for the hippocampal units inside CA1 ripples (32 of 39 recordings peak within +-2 ms) and for the LS units in high-power epochs (75 of 78). The 7 hippocampal recordings that peak elsewhere are all from DT5 (peaks at -20 to +10 ms), and the 3 LS recordings are two DT2 recordings and one DT9 recording. This check supports a common clock in most recordings and does not establish it in those 10.
- Sleep states and LFP: CA1 ripples are about 9 times as frequent in NREM as in WAKE with the states as stored, and about 2 times with the state table shifted by 10 min.
- Trials and states: 1.7% of trial time is scored NREM, against 31 to 39% with the trials shifted by 1 h.
The data are consistent with one clock for spikes, LFP, sleep states and trials.
dataset-composition¶
LS was recorded in 6 rats, with 1989 LS neurons; 37 recordings had CA1 and LS and 25 had CA3 and LS recorded simultaneously; LS HFO rates come from 89 recordings and CA1 ripple rates from 53.
In [7]:
def has_units(r, region):
loc = r["meta"]["units"]["location"]
return loc is not None and region in set(loc)
def has_el(r, region):
return region in set(r["meta"]["el_location"])
animals = sorted({r["meta"]["subject"] for r in R})
n_ls_units = sum(len(ls_hfo.sorted_spikes(r, "ls")) for r in R)
n_ls_tested = sum(sum(p["n"] > 50 for p in r["phase_locking"]) for r in R)
comp = pd.DataFrame({
"reading": ["(a) units", "(b) electrodes", "(c) electrodes and sleep states"],
"CA1+LS (37)": [sum(has_units(r, "ls") and has_units(r, "hpc") for r in results),
sum(has_el(r, "ls") and has_el(r, "hpc") for r in results), None],
"CA3+LS (25)": [sum(has_units(r, "ls") and has_units(r, "ca3") for r in results),
sum(has_el(r, "ls") and has_el(r, "ca3") for r in results), None],
"LS rate recordings (89)": [sum(has_units(r, "ls") for r in results), sum(has_el(r, "ls") for r in results),
sum(has_el(r, "ls") and r["meta"]["states"] is not None for r in results)],
"CA1 rate recordings (53)": [sum(has_units(r, "hpc") for r in results), sum(has_el(r, "hpc") for r in results),
sum(has_el(r, "hpc") and r["meta"]["states"] is not None for r in results)],
})
print(f"animals with LS recordings: {len(animals)} {animals}")
print(f"LS neurons: (a) all units labelled ls {n_ls_units}; (b) with more than 50 spikes in high-power epochs {n_ls_tested} (paper 1989)")
print(comp.to_string(index=False))
tol = {"CA1+LS (37)": (37, 3.7), "CA3+LS (25)": (25, 2.5), "LS rate recordings (89)": (89, 8.9), "CA1 rate recordings (53)": (53, 5.3)}
within = {k: [v is not None and abs(v - c) <= t for v in comp[k]] for k, (c, t) in tol.items()}
print("within tolerance by reading:", within)
print("LS neurons within tolerance:", {"(a)": abs(n_ls_units - 1989) <= 198.9, "(b)": abs(n_ls_tested - 1989) <= 198.9})
animals with LS recordings: 6 ['DT1', 'DT2', 'DT5', 'DT7', 'DT8', 'DT9']
LS neurons: (a) all units labelled ls 2448; (b) with more than 50 spikes in high-power epochs 2046 (paper 1989)
reading CA1+LS (37) CA3+LS (25) LS rate recordings (89) CA1 rate recordings (53)
(a) units 37.0 10.0 78 40
(b) electrodes 72.0 10.0 94 72
(c) electrodes and sleep states NaN NaN 53 43
within tolerance by reading: {'CA1+LS (37)': [True, False, False], 'CA3+LS (25)': [False, False, False], 'LS rate recordings (89)': [False, True, False], 'CA1 rate recordings (53)': [False, False, False]}
LS neurons within tolerance: {'(a)': False, '(b)': True}
Result: partially reproduced. The animal count matches (6). Three of the five other counts are
within tolerance under one reading: LS neurons with more than 50 spikes in high-power epochs
(2046 against 1989), files with LS and hpc units (37 against 37) and files with LS electrodes
(94 against 89). Files with LS and CA3 units number 10 against 25, and no reading gives 53 CA1
recordings (40, 72 or 43). All units labelled ls number 2448.
fig1c-hfo-each-animal¶
Transient 120-180 Hz HFOs occur in LS in each of the 6 animals, with a robust increase in the firing probability of LS neurons at the events.
Spike probability is the number of spikes in +-10 ms of the event time (the 20 bins of 1 ms from -10 to +10 ms), per event and neuron, against as many uniformly random times in the same recording. A variant compares the same window with each neuron's firing 50 to 100 ms before and after the same events.
In [8]:
# fig1c-hfo-each-animal: frequency and HFO-triggered LS firing by animal
rng = np.random.default_rng(102)
rows, curves, flank_rows = [], {}, []
for subj in animals:
freqs, ev_p, rnd_p, fl_p, n_events, peth_sum, peth_rnd = [], [], [], [], 0, 0, 0
for r in [r for r in R if r["meta"]["subject"] == subj]:
e = r["ls_events"]
freqs.append(e["frequency"])
n_events += len(e["peak"])
units = ls_hfo.sorted_spikes(r, "ls")
if len(e["peak"]) == 0 or not units:
continue
rnd_t = rng.uniform(0.2, r["meta"]["duration"] - 0.2, len(e["peak"]))
for s in units.values():
pe, pr = ls_hfo.peth(s, e["peak"]), ls_hfo.peth(s, rnd_t)
ev_p.append(pe[90:110].sum()) # +-10 ms around the event time (1 ms bins)
rnd_p.append(pr[90:110].sum())
fl_p.append((pe[:50].sum() + pe[150:].sum()) / 100 * 20) # 50-100 ms on each side, scaled to 20 bins
peth_sum = peth_sum + pe
peth_rnd = peth_rnd + pr
d = np.array(ev_p) - np.array(rnd_p)
dfl = np.array(ev_p) - np.array(fl_p)
flank_rows.append({"animal": subj, "LS neurons": len(dfl), "P(spike) +-10 ms at HFO": round(np.mean(ev_p), 4),
"flanks 50-100 ms": round(np.mean(fl_p), 4), "difference": round(dfl.mean(), 4),
"SEM": round(dfl.std(ddof=1) / np.sqrt(len(dfl)), 4),
"difference / SEM": round(dfl.mean() / (dfl.std(ddof=1) / np.sqrt(len(dfl))), 1)})
f = np.concatenate(freqs)
curves[subj] = (peth_sum / len(ev_p), peth_rnd / len(ev_p))
rows.append({"animal": subj, "HFOs": n_events, "median freq (Hz)": round(np.nanmedian(f), 1),
"LS neurons": len(ev_p), "P(spike) +-10 ms at HFO": round(np.mean(ev_p), 4),
"at random times": round(np.mean(rnd_p), 4), "difference": round(d.mean(), 4),
"SEM": round(d.std(ddof=1) / np.sqrt(len(d)), 4), "difference / SEM": round(d.mean() / (d.std(ddof=1) / np.sqrt(len(d))), 1)})
f1c = pd.DataFrame(rows)
f1c["freq in 120-180"] = f1c["median freq (Hz)"].between(120, 180)
f1c["above 3 SEM"] = f1c["difference / SEM"] > 3
print(f1c.to_string(index=False))
print("animals meeting both:", int((f1c["freq in 120-180"] & f1c["above 3 SEM"]).sum()), "of", len(f1c))
fig, axes = plt.subplots(1, len(animals), figsize=(15, 2.6), sharey=False)
lag = np.arange(-99.5, 100, 1)
for ax, subj in zip(axes, animals):
ax.plot(lag, curves[subj][0], "k", lw=1, label="HFO")
ax.plot(lag, curves[subj][1], "r", lw=1, label="random")
ax.set_title(subj); ax.set_xlabel("ms from HFO")
axes[0].set_ylabel("spikes per event per 1 ms bin"); axes[0].legend(fontsize=7)
plt.tight_layout(); plt.show()
print("plotted range (spikes per event per bin):", {s: (round(float(min(c[0].min(), c[1].min())), 4), round(float(max(c[0].max(), c[1].max())), 4)) for s, c in curves.items()})
# Variant (does not change the outcome): the same window against each neuron's own firing 50 to 100 ms
# before and after the same events, a baseline taken in the same moments as the events.
f1c_flank = pd.DataFrame(flank_rows)
print("\nvariant, baseline from the flanks of the PETH (50-100 ms either side of the event):")
print(f1c_flank.to_string(index=False))
print("animals above 3 SEM against the flanks:", int((f1c_flank["difference / SEM"] > 3).sum()), "of", len(f1c_flank))
animal HFOs median freq (Hz) LS neurons P(spike) +-10 ms at HFO at random times difference SEM difference / SEM freq in 120-180 above 3 SEM DT1 5292 155.5 610 0.0536 0.0350 0.0187 0.0031 6.0 True True DT2 1183 156.9 396 0.0563 0.0250 0.0314 0.0248 1.3 True False DT5 4032 155.6 560 0.0935 0.0592 0.0343 0.0252 1.4 True False DT7 38008 152.2 51 0.0045 0.0097 -0.0052 0.0015 -3.4 True False DT8 2757 156.2 296 0.0107 0.0172 -0.0065 0.0014 -4.6 True False DT9 15914 157.8 510 0.0151 0.0169 -0.0018 0.0013 -1.4 True False animals meeting both: 1 of 6
plotted range (spikes per event per bin): {'DT1': (0.0013, 0.0048), 'DT2': (0.0006, 0.0051), 'DT5': (0.0017, 0.0085), 'DT7': (0.0, 0.0008), 'DT8': (0.0003, 0.0013), 'DT9': (0.0004, 0.0013)}
variant, baseline from the flanks of the PETH (50-100 ms either side of the event):
animal LS neurons P(spike) +-10 ms at HFO flanks 50-100 ms difference SEM difference / SEM
DT1 610 0.0536 0.0386 0.0151 0.0029 5.2
DT2 396 0.0563 0.0321 0.0242 0.0155 1.6
DT5 560 0.0935 0.0565 0.0370 0.0218 1.7
DT7 51 0.0045 0.0056 -0.0012 0.0004 -2.6
DT8 296 0.0107 0.0133 -0.0027 0.0011 -2.4
DT9 510 0.0151 0.0145 0.0006 0.0013 0.5
animals above 3 SEM against the flanks: 1 of 6
Variants that do not change the outcome: only the events in NREM, and only the events with a CA1 ripple within 25 ms (animals without states or hippocampal channels drop out).
In [9]:
# Variant (does not change the outcome): NREM events only, and events with a CA1 ripple within 25 ms.
rng = np.random.default_rng(1021)
def triggered(sel_fn):
out = {}
for subj in animals:
d = []
for r in [r for r in R if r["meta"]["subject"] == subj]:
e = r["ls_events"]
t = e["peak"][sel_fn(r, e)]
units = ls_hfo.sorted_spikes(r, "ls")
if len(t) < 5 or not units:
continue
rnd_t = rng.uniform(0.2, r["meta"]["duration"] - 0.2, len(t))
d += [ls_hfo.peth(s, t)[90:110].sum() - ls_hfo.peth(s, rnd_t)[90:110].sum() for s in units.values()]
d = np.array(d)
out[subj] = (len(d), round(d.mean(), 4) if len(d) else None, round(d.mean() / (d.std(ddof=1) / np.sqrt(len(d))), 1) if len(d) > 2 else None)
return out
def coupled(r, e):
if "hpc_events" not in r:
return np.zeros(len(e["peak"]), bool)
rp = np.sort(r["hpc_events"]["peak"])
k = np.clip(np.searchsorted(rp, e["peak"]), 1, max(len(rp) - 1, 1))
near = np.minimum(np.abs(e["peak"] - rp[k - 1]), np.abs(e["peak"] - rp[np.minimum(k, len(rp) - 1)])) if len(rp) else np.full(len(e["peak"]), np.inf)
return near <= 0.025
print("NREM events only (neurons, mean difference, difference/SEM):", triggered(lambda r, e: e["state"] == "Non-REM"))
print("events with a CA1 ripple within 25 ms:", triggered(coupled))
NREM events only (neurons, mean difference, difference/SEM): {'DT1': (0, None, None), 'DT2': (158, np.float64(0.0384), np.float64(3.3)), 'DT5': (378, np.float64(-0.0026), np.float64(-0.4)), 'DT7': (51, np.float64(-0.0026), np.float64(-0.7)), 'DT8': (0, None, None), 'DT9': (478, np.float64(0.0092), np.float64(3.9))}
events with a CA1 ripple within 25 ms: {'DT1': (0, None, None), 'DT2': (396, np.float64(0.0286), np.float64(1.4)), 'DT5': (309, np.float64(0.0803), np.float64(2.0)), 'DT7': (51, np.float64(-0.0047), np.float64(-1.7)), 'DT8': (296, np.float64(-0.0087), np.float64(-5.7)), 'DT9': (352, np.float64(0.0117), np.float64(4.2))}
Result: not reproduced. Every animal's detected events have a median frequency of 152 to 158 Hz, inside 120 to 180 Hz, but LS firing within +-10 ms of the events exceeds firing at random times by more than 3 SEM in DT1 only (6.0 SEM). In DT2 and DT5 it is higher by 1.3 and 1.4 SEM, and in DT7, DT8 and DT9 it is lower (-3.4, -4.6 and -1.4 SEM). With NREM events only, DT2 (3.3 SEM) and DT9 (3.9 SEM) pass; with ripple-coupled events, DT9 passes (4.2 SEM). On the re-referenced channel (see the check after fig5f) no animal passes (largest 1.5 SEM, DT2; DT8 -4.4 SEM). Random times are a baseline that mixes states. Against each neuron's firing 50 to 100 ms before and after the same events, firing at the events is higher in DT1, DT2, DT5 and DT9 (DT9 0.0151 against 0.0145) and lower in DT7 and DT8, and exceeds 3 SEM only in DT1 (5.2 SEM; DT2 1.6, DT5 1.7, DT9 0.5). In DT2 and DT5 the differences are large (0.024 and 0.037) but vary widely across neurons. The test does not show raised LS firing at the events outside DT1; it does not show that the events lack it.
fig1c-event-count¶
Across 6 animals, 29,955 LS HFO events were analysed.
In [10]:
# fig1c-event-count
n_all = sum(len(r["ls_events"]["peak"]) for r in R)
n_units = sum(len(r["ls_events"]["peak"]) for r in R if has_units(r, "ls"))
print(f"LS HFOs: (a) all {len(R)} recordings with LS electrodes {n_all}; (b) the {sum(has_units(r, 'ls') for r in R)} recordings with LS units {n_units} (paper 29,955, tolerance 2,995.5)")
print("by animal (a):", {s: sum(len(r["ls_events"]["peak"]) for r in R if r["meta"]["subject"] == s) for s in animals})
by_state = collections.Counter(np.concatenate([r["ls_events"]["state"] for r in R]))
print("by sleep state at the event (all recordings; '' = no states or unscored):", dict(by_state))
LS HFOs: (a) all 94 recordings with LS electrodes 67186; (b) the 78 recordings with LS units 31984 (paper 29,955, tolerance 2,995.5)
by animal (a): {'DT1': 5292, 'DT2': 1183, 'DT5': 4032, 'DT7': 38008, 'DT8': 2757, 'DT9': 15914}
by sleep state at the event (all recordings; '' = no states or unscored): {'': 17366, 'Awake': 27676, 'Non-REM': 20019, 'MA': 1217, 'REM': 908}
Result: reproduced. In the recordings with LS units, 31,984 events are detected, within 10% of 29,955. Over all recordings with LS electrodes there are 67,186, 38,008 of them from DT7, and more fall in WAKE (27,676) than in NREM (20,019). The match of reading (b) does not mean the same events were found, and the count depends on the channel: on the re-referenced channel the recordings with LS units give 66,458 events.
figs1-participation¶
On average 16% of LS neurons take part in a given HFO, firing an average of 0.23 spikes per event.
In [11]:
# figs1-participation
frac_ev, spk_ev, per_rec = [], [], []
for r in R:
units = ls_hfo.sorted_spikes(r, "ls")
e = r["ls_events"]
if not units or len(e["peak"]) == 0:
continue
C = np.array([ls_hfo.counts_in(s, e["start"], e["stop"]) for s in units.values()]) # neurons x events
frac_ev.append((C > 0).mean(axis=0))
spk_ev.append(C.ravel())
per_rec.append({"file": r["path"].split("/")[1][:30], "neurons": len(units), "events": C.shape[1],
"fraction": (C > 0).mean(), "spikes": C.mean()})
per_rec = pd.DataFrame(per_rec)
fa, sa = np.concatenate(frac_ev).mean(), np.concatenate(spk_ev).mean()
print(f"(a) pooled over {sum(len(x) for x in frac_ev)} events: fraction of LS neurons firing {fa:.3f}, spikes per neuron per event {sa:.3f}")
print(f"(b) mean of {len(per_rec)} recordings: fraction {per_rec.fraction.mean():.3f}, spikes {per_rec.spikes.mean():.3f}")
print("paper: 0.16 (tolerance 0.016) and 0.23 (0.023)")
print(f"ratio to the paper's values: fraction {fa / 0.16:.2f} (pooled), {per_rec.fraction.mean() / 0.16:.2f} (recordings);"
f" spikes {sa / 0.23:.2f} (pooled), {per_rec.spikes.mean() / 0.23:.2f} (recordings)")
print(f"per-recording fraction: median {per_rec.fraction.median():.3f}, range {per_rec.fraction.min():.3f} to {per_rec.fraction.max():.3f}")
(a) pooled over 31984 events: fraction of LS neurons firing 0.027, spikes per neuron per event 0.051 (b) mean of 76 recordings: fraction 0.035, spikes 0.068 paper: 0.16 (tolerance 0.016) and 0.23 (0.023) ratio to the paper's values: fraction 0.17 (pooled), 0.22 (recordings); spikes 0.22 (pooled), 0.30 (recordings) per-recording fraction: median 0.017, range 0.000 to 0.182
Result: not reproduced. Pooled over events, 2.7% of LS neurons fire in an HFO, with 0.051 spikes per neuron per event; the mean over recordings gives 3.5% and 0.068. The fraction is 0.17 (pooled) to 0.22 (recordings) of the paper's value, and the spike count 0.22 to 0.30 of it. The highest recording reaches 18.2%. On the re-referenced channel the pooled values are 4.5% and 0.060.
fig1d-phase-locking¶
The vast majority of LS neurons (1840 of 1989, 92%) are significantly phase-locked to the 120-200 Hz band (Rayleigh P < 0.05), with a strong preference for the trough of each cycle.
The phase is taken on each recording's HFO detection channel; the paper does not say which channel it used.
In [12]:
# fig1d-phase-locking
pl = pd.DataFrame([{"animal": r["meta"]["subject"], **{k: p[k] for k in ("n", "p", "r", "phase")},
"r_null": np.nanmean(p["r_null"])} for r in R for p in r["phase_locking"]])
tested = pl[pl.n > 50]
sig = tested[tested.p < 0.05]
print(f"LS neurons: {len(pl)}; with more than 50 spikes in high-power epochs: {len(tested)}; significant (Rayleigh p < 0.05): {len(sig)}")
print(f"fraction significant: (a) of tested {len(sig) / len(tested):.3f}; (b) of all LS neurons {len(sig) / len(pl):.3f} (paper 0.925, tolerance 0.0925)")
print(f"resultant length of tested neurons: median {tested.r.median():.3f}; with the phase circularly shifted: median {tested.r_null.median():.3f}")
mean_dir = np.angle(np.exp(1j * sig.phase).mean())
rt = len(sig) * np.abs(np.exp(1j * sig.phase).mean())
p_pref = np.exp(np.sqrt(1 + 4 * len(sig) + 4 * (len(sig) ** 2 - rt ** 2)) - (1 + 2 * len(sig)))
print(f"preferred phases of significant neurons: circular mean {np.degrees(mean_dir) % 360:.1f} deg (trough = 180),"
f" resultant {np.abs(np.exp(1j * sig.phase).mean()):.3f}, Rayleigh p = {p_pref:.2g}")
print(f"within 45 deg of the trough: {abs(abs(np.degrees(mean_dir)) - 180) <= 45}")
print("by animal:")
print(tested.groupby("animal").apply(lambda d: pd.Series({"tested": len(d), "significant": int((d.p < 0.05).sum()),
"fraction": round((d.p < 0.05).mean(), 3), "median r": round(d.r.median(), 3),
"mean phase (deg)": round(np.degrees(np.angle(np.exp(1j * d.phase[d.p < 0.05]).mean())) % 360, 1)})).to_string())
fig, axes = plt.subplots(1, 2, figsize=(9, 3))
axes[0].hist(tested.r, bins=50, range=(0, 0.6), color="k", histtype="step", label="data")
axes[0].hist(tested.r_null, bins=50, range=(0, 0.6), color="r", histtype="step", label="shifted phase")
axes[0].set_xlabel("resultant length"); axes[0].legend()
axes[1].hist(np.degrees(sig.phase) % 360, bins=36, range=(0, 360), color="k")
axes[1].set_xlabel("preferred phase (deg; 180 = trough)")
plt.tight_layout(); plt.show()
print(f"resultant lengths above 0.6 (outside the axis): data {(tested.r > 0.6).sum()}, shifted {(tested.r_null > 0.6).sum()}")
LS neurons: 2448; with more than 50 spikes in high-power epochs: 2046; significant (Rayleigh p < 0.05): 1151
fraction significant: (a) of tested 0.563; (b) of all LS neurons 0.470 (paper 0.925, tolerance 0.0925)
resultant length of tested neurons: median 0.124; with the phase circularly shifted: median 0.048
preferred phases of significant neurons: circular mean 189.3 deg (trough = 180), resultant 0.755, Rayleigh p = 0
within 45 deg of the trough: True
by animal:
tested significant fraction median r mean phase (deg)
animal
DT1 559.0 391.0 0.699 0.151 190.5
DT2 304.0 102.0 0.336 0.085 201.1
DT5 463.0 326.0 0.704 0.160 187.2
DT7 35.0 2.0 0.057 0.053 296.4
DT8 241.0 136.0 0.564 0.119 185.8
DT9 444.0 194.0 0.437 0.098 187.1
resultant lengths above 0.6 (outside the axis): data 82, shifted 0
In [13]:
# Check: could the trough preference come from each unit's own spikes in the LFP? The phase is taken
# again on an LS channel on another shank than the unit's: the detection channel for units on other
# shanks, and for units on the detection shank the LS channel with the highest 120-200 Hz SD on
# another shank. Recordings with LS channels on one shank only drop out.
def circ_summary(ph):
c = np.exp(1j * np.asarray(ph)).mean()
n = len(ph)
R_ = n * np.abs(c)
p_ = np.exp(np.sqrt(1 + 4 * n + 4 * (n ** 2 - R_ ** 2)) - (1 + 2 * n))
return np.degrees(np.angle(c)) % 360, np.abs(c), p_
rows = []
for r in R:
grp = r["meta"]["units"]["group"]
if grp is None or "phase_locking_other" not in r:
continue
for k, i in enumerate(r["ls_units"]):
same = grp[i] == r["channel_groups"]["ls"]
d_, o_ = r["phase_locking"][k], r["phase_locking_other"][k]
off = o_ if same else d_
rows.append({"animal": r["meta"]["subject"], "on detection shank": same,
"n_det": d_["n"], "p_det": d_["p"], "r_det": d_["r"], "ph_det": d_["phase"],
"n_off": off["n"], "p_off": off["p"], "r_off": off["r"], "ph_off": off["phase"]})
ps = pd.DataFrame(rows)
print(f"LS units in recordings with LS channels on two or more shanks: {len(ps)} of {len(pl)}")
for label, d in (("units on the detection shank, detection channel", ps[ps["on detection shank"]].rename(columns=lambda c: c.replace("_det", ""))),
("units on other shanks, detection channel", ps[~ps["on detection shank"]].rename(columns=lambda c: c.replace("_det", ""))),
("all units, phase from another shank", ps.rename(columns=lambda c: c.replace("_off", "")))):
t_ = d[d.n > 50]
s_ = t_[t_.p < 0.05]
mdeg, res_, pp_ = circ_summary(s_.ph) if len(s_) else (np.nan, np.nan, np.nan)
print(f"{label}: tested {len(t_)}, significant {len(s_)} ({len(s_) / max(len(t_), 1):.3f}); median resultant {t_.r.median():.3f};"
f" above 0.6: {(t_.r > 0.6).sum()}; preferred phase {mdeg:.1f} deg (resultant {res_:.3f}, Rayleigh p = {pp_:.2g})")
t_ = ps[ps.n_off > 50]
print("phase from another shank, by animal:")
print(t_.groupby("animal").apply(lambda d: pd.Series({"tested": len(d), "fraction": round((d.p_off < 0.05).mean(), 3),
"mean phase (deg)": round(circ_summary(d.ph_off[d.p_off < 0.05])[0], 1) if (d.p_off < 0.05).any() else np.nan})).to_string())
LS units in recordings with LS channels on two or more shanks: 2388 of 2448
units on the detection shank, detection channel: tested 461, significant 352 (0.764); median resultant 0.191; above 0.6: 44; preferred phase 186.7 deg (resultant 0.800, Rayleigh p = 7.7e-123)
units on other shanks, detection channel: tested 1550, significant 797 (0.514); median resultant 0.110; above 0.6: 38; preferred phase 190.5 deg (resultant 0.737, Rayleigh p = 6.2e-225)
all units, phase from another shank: tested 2003, significant 1048 (0.523); median resultant 0.110; above 0.6: 44; preferred phase 189.7 deg (resultant 0.743, Rayleigh p = 8.4e-302)
phase from another shank, by animal:
tested fraction mean phase (deg)
animal
DT1 558.0 0.674 193.6
DT2 299.0 0.274 206.0
DT5 459.0 0.704 186.8
DT8 241.0 0.494 181.0
DT9 446.0 0.332 185.5
Result: partially reproduced. 1151 of 2046 tested LS neurons (56.3%) are significantly phase-locked, 47.0% of all 2448; the paper reports 92.5%. The trough preference holds: the preferred phases of significant neurons centre on 189 degrees (resultant 0.755; the Rayleigh p underflows to 0). Fractions range from 5.7% (DT7) to 70.4% (DT5) by animal. Median resultant lengths are 0.124 for the data and 0.048 with the phase circularly shifted.
The trough preference does not come from each unit's own spikes in the LFP. With the phase taken on an LS channel on another shank than the unit's, 1048 of 2003 tested neurons (52.3%) are locked, with preferred phases centred on 189.7 degrees (resultant 0.743), and 44 have a resultant above 0.6. Locking is stronger for the units on the detection shank (76.4% of 461, median resultant 0.191) than for those on other shanks (51.4% of 1550, 0.110), so part of the locking measured on the unit's own shank may come from its spikes or from its distance to the channel.
fig1f-shorter-than-ripples and fig1g-slower-than-ripples¶
LS HFOs are shorter than CA1 ripples from the same animals (27 +- 9 against 31 +- 14 ms; rank-sum P < 0.01).
LS HFOs are slower in frequency than CA1 ripples (147 +- 5.7 against 159 +- 8 Hz; rank-sum P < 0.01).
In [14]:
# fig1f-shorter-than-ripples and fig1g-slower-than-ripples: animals with both LS and hippocampal channels
both = [r for r in R if "hpc_events" in r]
ls_d = np.concatenate([r["ls_events"]["duration"] for r in both]) * 1000
ca_d = np.concatenate([r["hpc_events"]["duration"] for r in both]) * 1000
ls_f = np.concatenate([r["ls_events"]["frequency"] for r in both])
ca_f = np.concatenate([r["hpc_events"]["frequency"] for r in both])
ls_f, ca_f = ls_f[np.isfinite(ls_f)], ca_f[np.isfinite(ca_f)]
print(f"{len(both)} recordings with both channels, animals {sorted({r['meta']['subject'] for r in both})}")
for name, a, b, unit, ref in (("duration", ls_d, ca_d, "ms", "27 +- 9 against 31 +- 14"), ("frequency", ls_f, ca_f, "Hz", "147 +- 5.7 against 159 +- 8")):
p = stats.mannwhitneyu(a, b).pvalue
print(f"{name}: LS {a.mean():.1f} +- {a.std():.1f} {unit} (median {np.median(a):.1f}, n = {len(a)}); "
f"CA1 {b.mean():.1f} +- {b.std():.1f} (median {np.median(b):.1f}, n = {len(b)}); rank-sum p = {p:.2g} [paper {ref}]")
print("variant, all recordings with LS electrodes:",
f"LS duration {np.concatenate([r['ls_events']['duration'] for r in R]).mean() * 1000:.1f} ms,",
f"frequency {np.nanmean(np.concatenate([r['ls_events']['frequency'] for r in R])):.1f} Hz")
rows = []
for subj in sorted({r["meta"]["subject"] for r in both}):
rr = [r for r in both if r["meta"]["subject"] == subj]
a = np.concatenate([r["ls_events"]["duration"] for r in rr]) * 1000
b = np.concatenate([r["hpc_events"]["duration"] for r in rr]) * 1000
fa = np.concatenate([r["ls_events"]["frequency"] for r in rr])
fb = np.concatenate([r["hpc_events"]["frequency"] for r in rr])
rows.append({"animal": subj, "LS dur": round(a.mean(), 1), "CA1 dur": round(b.mean(), 1),
"p dur": stats.mannwhitneyu(a, b).pvalue, "LS freq": round(np.nanmean(fa), 1),
"CA1 freq": round(np.nanmean(fb), 1), "p freq": stats.mannwhitneyu(fa[np.isfinite(fa)], fb[np.isfinite(fb)]).pvalue})
print(pd.DataFrame(rows).to_string(index=False))
# NREM-only variant
st_ls = np.concatenate([r["ls_events"]["state"] for r in both]); st_ca = np.concatenate([r["hpc_events"]["state"] for r in both])
a, b = ls_d[st_ls == "Non-REM"], ca_d[st_ca == "Non-REM"]
fa = np.concatenate([r["ls_events"]["frequency"] for r in both])[st_ls == "Non-REM"]; fb = np.concatenate([r["hpc_events"]["frequency"] for r in both])[st_ca == "Non-REM"]
print(f"variant, NREM events only: duration LS {a.mean():.1f} +- {a.std():.1f}, CA1 {b.mean():.1f} +- {b.std():.1f} ms (p = {stats.mannwhitneyu(a, b).pvalue:.2g});"
f" frequency LS {np.nanmean(fa):.1f} +- {np.nanstd(fa):.1f}, CA1 {np.nanmean(fb):.1f} +- {np.nanstd(fb):.1f} Hz (p = {stats.mannwhitneyu(fa[np.isfinite(fa)], fb[np.isfinite(fb)]).pvalue:.2g})")
fig, axes = plt.subplots(1, 2, figsize=(9, 3))
axes[0].hist(ca_d, bins=60, range=(0, 150), density=True, histtype="step", color="k", label="CA1")
axes[0].hist(ls_d, bins=60, range=(0, 150), density=True, histtype="step", color="m", label="LS")
axes[0].set_xlabel("duration (ms)"); axes[0].legend()
axes[1].hist(ca_f, bins=60, range=(100, 220), density=True, histtype="step", color="k")
axes[1].hist(ls_f, bins=60, range=(100, 220), density=True, histtype="step", color="m")
axes[1].set_xlabel("peak frequency (Hz)")
plt.tight_layout(); plt.show()
print(f"outside the axes: durations above 150 ms LS {(ls_d > 150).sum()}, CA1 {(ca_d > 150).sum()};"
f" frequencies outside 100-220 Hz LS {((ls_f < 100) | (ls_f > 220)).sum()}, CA1 {((ca_f < 100) | (ca_f > 220)).sum()}")
72 recordings with both channels, animals ['DT2', 'DT5', 'DT7', 'DT8', 'DT9'] duration: LS 33.9 +- 25.0 ms (median 23.2, n = 56871); CA1 38.3 +- 23.7 (median 32.0, n = 230161); rank-sum p = 0 [paper 27 +- 9 against 31 +- 14] frequency: LS 153.6 +- 8.3 Hz (median 154.1, n = 56871); CA1 150.8 +- 9.2 (median 149.8, n = 230161); rank-sum p = 0 [paper 147 +- 5.7 against 159 +- 8] variant, all recordings with LS electrodes: LS duration 34.4 ms, frequency 153.9 Hz
animal LS dur CA1 dur p dur LS freq CA1 freq p freq DT2 25.6 38.5 7.598951e-201 156.8 154.8 1.745084e-22 DT5 22.5 35.1 3.447988e-201 156.1 147.8 0.000000e+00 DT7 35.2 40.2 0.000000e+00 151.8 150.3 0.000000e+00 DT8 27.5 28.0 8.320209e-01 155.6 156.2 4.893006e-02 DT9 33.8 37.4 1.204100e-304 157.8 150.9 0.000000e+00
variant, NREM events only: duration LS 38.0 +- 25.7, CA1 38.4 +- 22.3 ms (p = 8.3e-103); frequency LS 150.3 +- 8.0, CA1 151.1 +- 9.0 Hz (p = 5.4e-12)
outside the axes: durations above 150 ms LS 290, CA1 1011; frequencies outside 100-220 Hz LS 0, CA1 0
Result: fig1f-shorter-than-ripples reproduced. LS HFOs last 33.9 +- 25.0 ms and CA1 ripples 38.3 +- 23.7 ms, both inside the tolerances of the paper's SDs, and LS HFOs are shorter (medians 23.2 and 32.0 ms; the rank-sum p underflows to 0, as it does for any difference over 56,871 and 230,161 events). The means are 6 to 7 ms above the paper's. In four of five animals LS HFOs are significantly shorter; in DT8 they are not (27.5 against 28.0 ms, p = 0.83). With NREM events only, the means are 38.0 and 38.4 ms (p = 8.3e-103), and on the re-referenced channel 38.2 and 38.3 ms, so the 4 ms gap in the means is not seen in NREM events or on the re-referenced channel.
Result: fig1g-slower-than-ripples not reproduced. LS HFOs have a mean frequency of 153.6 +- 8.3 Hz and CA1 ripples 150.8 +- 9.2 Hz: the LS events are faster, the reverse of the claim, and both means are outside the tolerances (147 +- 5.7 and 159 +- 8). LS events are faster in four of five animals. With NREM events only, LS are slower by 0.8 Hz (150.3 against 151.1 Hz). On the re-referenced channel LS events are faster (156.3 against 150.8 Hz).
fig1h-no-amplitude-frequency-correlation¶
LS HFOs show no significant correlation between event amplitude and peak frequency.
In [15]:
# fig1h-no-amplitude-frequency-correlation
amp = np.concatenate([r["ls_events"]["amplitude"] for r in R])
frq = np.concatenate([r["ls_events"]["frequency"] for r in R])
ok = np.isfinite(frq)
pr = stats.pearsonr(amp[ok], frq[ok])
sr = stats.spearmanr(amp[ok], frq[ok])
print(f"all {ok.sum()} LS HFOs: Pearson r = {pr.statistic:.3f}, p = {pr.pvalue:.2g}; Spearman rho = {sr.statistic:.3f}, p = {sr.pvalue:.2g}")
rows = []
for subj in animals:
a = np.concatenate([r["ls_events"]["amplitude"] for r in R if r["meta"]["subject"] == subj])
f = np.concatenate([r["ls_events"]["frequency"] for r in R if r["meta"]["subject"] == subj])
k = np.isfinite(f)
t = stats.pearsonr(a[k], f[k])
rows.append({"animal": subj, "events": int(k.sum()), "r": round(t.statistic, 3), "p": t.pvalue})
print(pd.DataFrame(rows).to_string(index=False))
st = np.concatenate([r["ls_events"]["state"] for r in R])
k = ok & (st == "Non-REM")
t = stats.pearsonr(amp[k], frq[k])
print(f"variant, NREM events only ({k.sum()}): r = {t.statistic:.3f}, p = {t.pvalue:.2g}")
fig, ax = plt.subplots(figsize=(4, 3))
ax.plot(frq[ok], amp[ok], ".", ms=1, alpha=0.3, color="m")
ax.set_xlim(100, 220); ax.set_ylim(4, 30); ax.set_xlabel("peak frequency (Hz)"); ax.set_ylabel("amplitude (SD)")
plt.tight_layout(); plt.show()
print(f"points outside the axes: {((frq[ok] < 100) | (frq[ok] > 220) | (amp[ok] > 30)).sum()} of {ok.sum()}")
all 67186 LS HFOs: Pearson r = 0.056, p = 2.4e-47; Spearman rho = 0.123, p = 3.4e-223 animal events r p DT1 5292 0.038 5.274922e-03 DT2 1183 0.004 8.905490e-01 DT5 4032 0.034 2.880132e-02 DT7 38008 0.080 4.102187e-55 DT8 2757 0.066 4.836785e-04 DT9 15914 -0.009 2.776492e-01 variant, NREM events only (20019): r = 0.028, p = 7.3e-05
points outside the axes: 3241 of 67186
Result: not reproduced. Over 67,186 events the Pearson correlation is small but significant (r = 0.056, p = 2.4e-47; Spearman rho = 0.123). By animal, r is -0.009 to 0.080, significant in four of six; with NREM events only, r = 0.028 (p = 7.3e-05). The outcome follows the criterion's test; the correlation is negligible in size (r = 0.056 explains 0.3% of the variance) and is significant because of the event count. The paper does not name its test.
fig2-confined-to-ls¶
HFO events are seen only at recording sites within LS, not at sites in neighbouring structures.
DT7 has probes that span cingulate cortex (acg), corpus callosum (cc), medial septum (ms) and
fimbria (sFi) with LS. For each channel, the 120-200 Hz power within +-20 ms of the LS HFO times in
up to 30 min of all channels, divided by the power at random times.
In [16]:
# fig2-confined-to-ls
rows = []
for r in R:
if "localization" not in r:
continue
L, loc = r["localization"], r["meta"]["el_location"]
good = r["channel_stats"]["good"]
row = {"file": r["path"].split("/")[1][:34], "HFOs in block": L["n_events"],
"block (min)": round((L["window"][1] - L["window"][0]) / 60, 1)}
for g in sorted(set(loc)):
row[g] = round(float(np.nanmedian(L["ratio"][(loc == g) & good])), 2)
rows.append(row)
loc_df = pd.DataFrame(rows)
print("median over channels of HFO-triggered / random 120-200 Hz power, by region:")
print(loc_df.to_string(index=False))
others = [c for c in loc_df.columns if c not in ("file", "HFOs in block", "block (min)", "ls", "hpc", "ca3", "unknown")]
used = loc_df[loc_df["HFOs in block"] >= 20]
win = used.apply(lambda row: all(row["ls"] > row[c] for c in others if pd.notna(row[c])), axis=1)
print(f"recordings with at least 20 HFOs in the block: {len(used)}; LS above every other non-hippocampal region in {int(win.sum())}")
median over channels of HFO-triggered / random 120-200 Hz power, by region:
file HFOs in block block (min) acg cc hpc ls ms sFi
sub-DT7_ses-20170315-288um-0um-170 0 7.0 NaN NaN NaN NaN NaN NaN
sub-DT7_ses-20170318-288um-0um-170 332 30.0 5.31 NaN 2.37 5.38 5.05 NaN
sub-DT7_ses-20170320-432um-0um-mer 36 30.0 1.14 NaN 1.84 5.18 6.05 NaN
sub-DT7_ses-20170321-216um-0um-170 50 30.0 5.75 NaN 1.86 8.25 8.05 NaN
sub-DT7_ses-20170324-576um-144um-1 1 30.0 2.97 NaN 0.52 49.91 134.50 NaN
sub-DT7_ses-20170325-576um-360um-1 6 30.0 1.88 NaN 2.48 2.64 4.57 NaN
sub-DT7_ses-20170401-864um-936um-m 151 30.0 NaN 50.11 1.72 42.02 38.50 NaN
sub-DT7_ses-20170412-1440um-1170um 35 30.0 NaN NaN 2.31 64.80 NaN 72.09
recordings with at least 20 HFOs in the block: 5; LS above every other non-hippocampal region in 2
Result: not reproduced. Five recordings, all from DT7 (the only animal whose files label electrodes outside LS and the hippocampus), have at least 20 HFOs in the block. LS channels have the highest ratio in two: 2017-03-18 (5.38 against 5.31 for cingulate cortex) and 2017-03-21 (8.25 against 8.05 for medial septum). In the other three, medial septum (6.05 against 5.18), corpus callosum (50.1 against 42.0) or fimbria (72.1 against 64.8) channels are higher. With each channel re-referenced to the median of its shank and the re-referenced LS events as triggers, four recordings have at least 20 events and LS is highest in two; medial septum channels are higher in the other two (4.00 against 3.09 and 4.21 against 3.65). The detected events reach channels outside LS at similar or larger relative power. The test rests on one animal.
fig2ef-intermediate-ls (blocked)¶
HFO power and the proportion of phase-locked neurons are highest in intermediate LS and taper towards its dorsal and ventral borders.
The claim is blocked: it needs each recording's depth relative to the borders of LS, which the paper took from histology. The session names carry drive depths in um (one in DT1, two in the other animals, without saying which drive holds the LS probe). For the record, the HFO band SD of the detection channel, the LS HFO rate and the fraction of tested LS neurons that are phase-locked are shown against each depth within each animal.
In [17]:
# fig2ef-intermediate-ls (blocked): drive depth in the session name against HFO power and locking
import re
rows = []
for r in R:
name = r["path"].split("/")[1]
depths = [int(x) for x in re.findall(r"(\d+)um", name)]
t_ = [p_ for p_ in r["phase_locking"] if p_["n"] > 50]
rows.append({"animal": r["meta"]["subject"], "depth 1": depths[0], "depth 2": depths[1] if len(depths) > 1 else np.nan,
"band SD (uV)": r["channel_stats"]["band_sd"][r["channels"]["ls"]] * r["meta"]["lfp_conversion"] * 1e6,
"HFO rate (/min)": 60 * len(r["ls_events"]["peak"]) / r["meta"]["duration"],
"tested": len(t_), "locked fraction": np.mean([p_["p"] < 0.05 for p_ in t_]) if len(t_) >= 5 else np.nan})
dp = pd.DataFrame(rows)
out = []
for subj, d in dp.groupby("animal"):
for dcol in ("depth 1", "depth 2"):
if d[dcol].nunique() < 3:
continue
row = {"animal": subj, "depth": dcol, "recordings": len(d), "range (um)": f"{d[dcol].min():.0f}-{d[dcol].max():.0f}"}
for v in ("band SD (uV)", "HFO rate (/min)", "locked fraction"):
k = d[[dcol, v]].dropna()
row[f"rho {v}"] = round(stats.spearmanr(k[dcol], k[v]).statistic, 2) if len(k) >= 4 else np.nan
# where the largest value falls in the depth range (0 = shallowest, 1 = deepest)
row[f"argmax {v}"] = round((k.loc[k[v].idxmax(), dcol] - k[dcol].min()) / (k[dcol].max() - k[dcol].min()), 2) if len(k) >= 4 and k[dcol].max() > k[dcol].min() else np.nan
out.append(row)
print(pd.DataFrame(out).to_string(index=False))
print("recordings by animal with at least 5 tested neurons:", dp.dropna(subset=["locked fraction"]).groupby("animal").size().to_dict())
animal depth recordings range (um) rho band SD (uV) argmax band SD (uV) rho HFO rate (/min) argmax HFO rate (/min) rho locked fraction argmax locked fraction
DT1 depth 1 12 1584-2736 0.18 0.25 -0.30 0.00 -0.40 0.19
DT2 depth 1 13 3036-3612 -0.01 1.00 0.08 0.75 -0.32 0.57
DT2 depth 2 13 362-1360 0.03 1.00 0.03 0.86 -0.32 0.84
DT5 depth 1 20 720-2124 -0.14 0.05 -0.39 0.64 -0.33 0.81
DT5 depth 2 20 720-2448 -0.17 0.08 0.06 0.94 0.36 1.00
DT7 depth 1 16 72-1440 -0.58 0.79 0.11 0.79 NaN NaN
DT7 depth 2 16 0-1170 -0.61 0.98 0.26 0.98 NaN NaN
DT8 depth 1 17 0-612 -0.15 0.00 -0.07 0.12 -0.79 0.00
DT8 depth 2 17 0-2088 -0.16 0.00 -0.10 0.21 -0.84 0.03
DT9 depth 1 16 0-1116 0.37 1.00 0.16 0.00 -0.73 0.44
DT9 depth 2 16 0-1000 0.41 1.00 0.29 0.00 -0.76 0.00
recordings by animal with at least 5 tested neurons: {'DT1': 12, 'DT2': 11, 'DT5': 19, 'DT7': 2, 'DT8': 17, 'DT9': 14}
Result: blocked (metadata insufficient). Shown for the record. Within animals, the rank correlations with depth are -0.61 to 0.41 for the HFO band SD and -0.84 to 0.36 for the locked fraction. The band SD is largest at the shallowest or deepest recording (relative position at most 0.08 or at least 0.98) in 9 of 11 animal and depth series; the locked fraction peaks inside the range (relative position 0.19 to 0.84) in 5 of 9. Without the LS borders, and without knowing which drive holds the LS probe, these do not test the claim.
LS pairs: fig4b-pre-post-unchanged, fig4b-hfo-above-behav, fig4d-compression¶
The co-activation of LS neuron pairs during HFOs (excess correlation) is not modulated by waking experience on a familiar maze (PRE against POST, rank-sum P = 0.36).
Excess correlations during HFOs (PRE and POST) are higher than during maze traversal (rank-sum P < 0.001), and those of PRE and POST HFOs exceed chance from spikes jittered by +-200 ms (P < 0.01).
Cross-correlograms of LS pairs during HFOs correlate with those during traversal above chance up to about 20-fold temporal compression (P < 0.01), similarly for PRE and POST (P > 0.01 for all compression factors).
The pair selection of the criteria (pairs on different shanks with more than 5 spike pairs within +-200 ms in each epoch) is applied first; it decides which of these claims have enough data.
In [18]:
# Pairs of LS neurons on different shanks: CCGs of spikes inside LS HFOs (PRE, POST) and inside trials (BEHAV).
rng = np.random.default_rng(109)
FACTORS = [1, 2, 4, 6, 8, 10, 15, 20, 30, 40, 50]
NB = 400 # bins of 1 ms on each side
W200 = slice(NB - 200, NB + 201)
def warp_corr(c_hfo, c_beh, factor):
"""Correlation of the BEHAV CCG over +-120 ms with the HFO CCG over +-120/factor ms,
interpolated to the same 241 points."""
lags = np.arange(-NB, NB + 1)
x = np.linspace(-120 / factor, 120 / factor, 241)
a = np.interp(x, lags, c_hfo)
b = c_beh[NB - 120:NB + 121]
if a.std() == 0 or b.std() == 0:
return np.nan
return np.corrcoef(a, b)[0, 1]
pair_rows, comp_rows, kept_sessions = [], [], 0
t_start = time.time()
for r in R:
m = r["meta"]
units = ls_hfo.sorted_spikes(r, "ls")
if m["epochs"] is None or m["trials"] is None or len(units) < 2 or m["units"]["group"] is None:
continue
e = r["ls_events"]
order = np.argsort(e["start"])
es, ee, ep = e["start"][order], e["stop"][order], e["peak"][order]
pre = ep < m["epochs"]["stop"][0]
post = ep >= m["epochs"]["start"][2]
to = np.argsort(m["trials"]["start"])
ts, tt = m["trials"]["start"][to], m["trials"]["stop"][to]
sp = {"PRE": {}, "BEHAV": {}, "POST": {}}
for i, s in units.items():
sp["PRE"][i] = ls_hfo.in_intervals(s, es[pre], ee[pre])
sp["POST"][i] = ls_hfo.in_intervals(s, es[post], ee[post])
sp["BEHAV"][i] = ls_hfo.in_intervals(s, ts, tt)
grp = m["units"]["group"]
ids = list(units)
n_pairs = 0
for a_ in range(len(ids)):
for b_ in range(a_ + 1, len(ids)):
i, j = ids[a_], ids[b_]
if grp[i] == grp[j]:
continue
C = {k: ls_hfo.ccg(sp[k][i], sp[k][j], half=NB / 1000) for k in sp}
if any(C[k][W200].sum() <= 5 for k in C):
continue
row = {"file": r["path"].split("/")[1][:30], "subject": m["subject"]}
for k in C:
cn = C[k] / np.sqrt(len(sp[k][i]) * len(sp[k][j]))
row[k] = ls_hfo.excess_correlation(cn)
ji = np.sort(sp[k][i] + rng.uniform(-0.2, 0.2, len(sp[k][i])))
jj = np.sort(sp[k][j] + rng.uniform(-0.2, 0.2, len(sp[k][j])))
row[k + "_jitter"] = ls_hfo.excess_correlation(ls_hfo.ccg(ji, jj, half=NB / 1000) / np.sqrt(len(ji) * len(jj)))
row["bias_" + k] = cn[NB + 1:NB + 201].sum() - cn[NB - 200:NB].sum()
C[k] = cn
for k in ("PRE", "POST"):
shift = rng.integers(1, 2 * NB)
for fct in FACTORS:
comp_rows.append({"file": row["file"], "epoch": k, "factor": fct,
"corr": warp_corr(C[k], C["BEHAV"], fct),
"chance": warp_corr(np.roll(C[k], shift), C["BEHAV"], fct)})
pair_rows.append(row)
n_pairs += 1
kept_sessions += n_pairs > 0
pairs = pd.DataFrame(pair_rows)
comp = pd.DataFrame(comp_rows)
print(f"pairs kept: {len(pairs)} in {kept_sessions} recordings (paper 7,172 pairs in 53 sessions); {time.time() - t_start:.0f} s")
print("pairs by animal:", pairs.subject.value_counts().sort_index().to_dict())
print(pairs[["PRE", "BEHAV", "POST", "PRE_jitter", "BEHAV_jitter", "POST_jitter"]].describe().round(4).to_string())
pairs kept: 106 in 10 recordings (paper 7,172 pairs in 53 sessions); 6 s
pairs by animal: {'DT2': 4, 'DT5': 38, 'DT9': 64}
PRE BEHAV POST PRE_jitter BEHAV_jitter POST_jitter
count 106.0000 106.0000 106.0000 106.0000 106.0000 106.0000
mean 0.0111 0.0037 0.0109 0.0004 0.0001 0.0003
std 0.0272 0.0174 0.0314 0.0012 0.0007 0.0008
min -0.0001 -0.0056 0.0008 -0.0011 -0.0009 -0.0007
25% 0.0034 -0.0002 0.0033 -0.0003 -0.0002 -0.0002
50% 0.0062 0.0003 0.0059 0.0001 0.0000 0.0001
75% 0.0107 0.0015 0.0096 0.0008 0.0003 0.0006
max 0.2778 0.1625 0.3225 0.0059 0.0054 0.0049
Only 106 pairs from 10 recordings in 3 animals pass the selection, against 7,172 pairs in 53 sessions in the paper. Few LS spikes fall inside the detected HFOs (see figs1-participation), so most pairs have 5 or fewer spike pairs within +-200 ms in PRE or POST. The shortfall comes from this attempt's detection, not from the dandiset, which holds the spikes, epochs and trials of these recordings. The pair-level tests below use these 106 pairs and have far less power than the paper's. The per-session claims are blocked: 7 recordings have at least 5 pairs (fig4-sessions-few-change) and 4 have at least 10 (fig4c-temporal-bias); their outputs are printed for the record.
In [19]:
# fig4b-pre-post-unchanged and fig4b-hfo-above-behav
t = stats.mannwhitneyu(pairs.PRE, pairs.POST)
print(f"PRE against POST (rank-sum): medians {pairs.PRE.median():.4f} and {pairs.POST.median():.4f}, p = {t.pvalue:.3g};"
f" median of POST - PRE over pairs {np.median(pairs.POST - pairs.PRE):.4f}; Wilcoxon signed-rank p = {stats.wilcoxon(pairs.PRE, pairs.POST).pvalue:.3g} (variant)")
for k in ("PRE", "POST"):
t = stats.mannwhitneyu(pairs[k], pairs.BEHAV)
print(f"{k} against BEHAV: medians {pairs[k].median():.4f} and {pairs.BEHAV.median():.4f}, rank-sum p = {t.pvalue:.3g}")
for k in ("PRE", "BEHAV", "POST"):
t = stats.mannwhitneyu(pairs[k], pairs[k + "_jitter"])
print(f"{k} against jittered: medians {pairs[k].median():.4f} and {pairs[k + '_jitter'].median():.4f}, rank-sum p = {t.pvalue:.3g}")
print("by animal (median excess correlation):")
print(pairs.groupby("subject")[["PRE", "BEHAV", "POST"]].median().round(4).to_string())
print("by animal, PRE against POST rank-sum p:", pairs.groupby("subject").apply(lambda d: round(stats.mannwhitneyu(d.PRE, d.POST).pvalue, 3)).to_dict())
fig, ax = plt.subplots(figsize=(5, 3))
for k, c in (("PRE", "g"), ("BEHAV", "k"), ("POST", "b")):
ax.hist(pairs[k], bins=60, range=(-0.05, 0.25), histtype="step", color=c, label=k)
ax.set_xlabel("excess correlation"); ax.legend(); plt.tight_layout(); plt.show()
print("values outside the axis:", {k: int(((pairs[k] < -0.05) | (pairs[k] > 0.25)).sum()) for k in ("PRE", "BEHAV", "POST")})
PRE against POST (rank-sum): medians 0.0062 and 0.0059, p = 0.812; median of POST - PRE over pairs -0.0001; Wilcoxon signed-rank p = 0.617 (variant)
PRE against BEHAV: medians 0.0062 and 0.0003, rank-sum p = 9.12e-24
POST against BEHAV: medians 0.0059 and 0.0003, rank-sum p = 3.78e-25
PRE against jittered: medians 0.0062 and 0.0001, rank-sum p = 3.61e-31
BEHAV against jittered: medians 0.0003 and 0.0000, rank-sum p = 0.00223
POST against jittered: medians 0.0059 and 0.0001, rank-sum p = 7.78e-35
by animal (median excess correlation):
PRE BEHAV POST
subject
DT2 0.0307 0.0029 0.0209
DT5 0.0084 0.0014 0.0096
DT9 0.0039 0.0002 0.0037
by animal, PRE against POST rank-sum p: {'DT2': 0.486, 'DT5': 0.728, 'DT9': 0.877}
values outside the axis: {'PRE': 1, 'BEHAV': 0, 'POST': 1}
Result: fig4b-pre-post-unchanged reproduced. PRE and POST excess correlations do not differ (medians 0.0062 and 0.0059, rank-sum p = 0.81), on 106 pairs against the paper's 7,172.
Result: fig4b-hfo-above-behav reproduced. HFO excess correlations are higher than BEHAV (median 0.0003; PRE p = 9.1e-24, POST p = 3.8e-25), and above their jittered values (p = 3.6e-31 and 7.8e-35); BEHAV is also slightly above its jitter (p = 0.0022), as in the paper. Same 106 pairs.
In [20]:
# fig4-sessions-few-change
rows = []
for f, d in pairs.groupby("file"):
if len(d) >= 5:
p = stats.mannwhitneyu(d.PRE, d.POST).pvalue
rows.append({"file": f, "pairs": len(d), "p": p, "direction": "PRE < POST" if d.PRE.median() < d.POST.median() else "PRE > POST"})
sess = pd.DataFrame(rows)
sig = sess[sess.p < 0.01]
print(f"recordings with at least 5 pairs: {len(sess)}; PRE differs from POST (p < 0.01): {len(sig)} ({len(sig) / len(sess):.3f}; paper 7 of 53 = 0.132)")
print("directions among them:", sig.direction.value_counts().to_dict())
print(sess.sort_values("p").head(12).to_string(index=False))
recordings with at least 5 pairs: 7; PRE differs from POST (p < 0.01): 0 (0.000; paper 7 of 53 = 0.132)
directions among them: {}
file pairs p direction
sub-DT9_ses-20170518-828um-864 8 0.049883 PRE > POST
sub-DT9_ses-20170507-468um-36u 5 0.150794 PRE > POST
sub-DT9_ses-20170506-468um-36u 33 0.243200 PRE < POST
sub-DT5_ses-20161019-1624um-23 18 0.558337 PRE < POST
sub-DT5_ses-20161026-1944um-22 15 0.708923 PRE > POST
sub-DT9_ses-20170508-468um-36u 18 0.763745 PRE > POST
sub-DT5_ses-20161021-1840um-24 5 0.841270 PRE > POST
In [21]:
# fig4c-temporal-bias
rows = []
for f, d in pairs.groupby("file"):
if len(d) >= 10:
rows.append({"file": f, "pairs": len(d), "pre_behav": stats.pearsonr(d.bias_PRE, d.bias_BEHAV).statistic,
"behav_post": stats.pearsonr(d.bias_BEHAV, d.bias_POST).statistic})
tb = pd.DataFrame(rows)
t = stats.ttest_ind(tb.pre_behav, tb.behav_post)
print(f"{len(tb)} recordings with at least 10 pairs; correlation of biases PRE-BEHAV mean {tb.pre_behav.mean():.3f}, BEHAV-POST mean {tb.behav_post.mean():.3f}")
print(f"two-sample t-test p = {t.pvalue:.3g}; paired t-test p = {stats.ttest_rel(tb.pre_behav, tb.behav_post).pvalue:.3g} (variant)")
print(f"one-sample t-tests against 0: PRE-BEHAV p = {stats.ttest_1samp(tb.pre_behav, 0).pvalue:.3g}, BEHAV-POST p = {stats.ttest_1samp(tb.behav_post, 0).pvalue:.3g}")
fig, ax = plt.subplots(figsize=(3.5, 3.5))
ax.plot(tb.pre_behav, tb.behav_post, "k.")
ax.plot([-1, 1], [-1, 1], ":", color="grey"); ax.set_xlim(-1, 1); ax.set_ylim(-1, 1)
ax.set_xlabel("r(PRE, BEHAV)"); ax.set_ylabel("r(BEHAV, POST)"); plt.tight_layout(); plt.show()
4 recordings with at least 10 pairs; correlation of biases PRE-BEHAV mean 0.042, BEHAV-POST mean 0.084 two-sample t-test p = 0.805; paired t-test p = 0.805 (variant) one-sample t-tests against 0: PRE-BEHAV p = 0.696, BEHAV-POST p = 0.556
In [22]:
# fig4d-compression
rows = []
for (k, fct), d in comp.groupby(["epoch", "factor"]):
d = d.dropna()
rows.append({"epoch": k, "factor": fct, "pairs": len(d), "mean corr": round(d["corr"].mean(), 4),
"mean chance": round(d.chance.mean(), 4), "p vs chance": stats.wilcoxon(d["corr"], d.chance).pvalue})
cr = pd.DataFrame(rows)
pp = []
for fct in FACTORS:
a = comp[(comp.epoch == "PRE") & (comp.factor == fct)]["corr"].dropna()
b = comp[(comp.epoch == "POST") & (comp.factor == fct)]["corr"].dropna()
pp.append({"factor": fct, "PRE vs POST p": stats.mannwhitneyu(a, b).pvalue})
cr = cr.merge(pd.DataFrame(pp), on="factor")
print(cr.to_string(index=False))
above = {k: [f for f in FACTORS if cr[(cr.epoch == k) & (cr.factor == f)]["p vs chance"].iloc[0] < 0.01
and cr[(cr.epoch == k) & (cr.factor == f)]["mean corr"].iloc[0] > cr[(cr.epoch == k) & (cr.factor == f)]["mean chance"].iloc[0]] for k in ("PRE", "POST")}
print("factors where the correlation exceeds chance (p < 0.01):", above)
print("factors where PRE and POST differ (p <= 0.01):", [f for f, p in zip(FACTORS, cr.drop_duplicates("factor")["PRE vs POST p"]) if p <= 0.01])
epoch factor pairs mean corr mean chance p vs chance PRE vs POST p POST 1 77 0.0964 0.0006 0.000375 0.949116 POST 2 56 0.1210 -0.0068 0.000032 0.963387 POST 4 36 0.0811 -0.0042 0.002079 0.602634 POST 6 27 0.0647 -0.0032 0.028018 0.260503 POST 8 25 0.0604 -0.0075 0.133640 0.088583 POST 10 23 0.0651 -0.0040 0.044876 0.038438 POST 15 21 0.0438 -0.0233 0.137283 0.456961 POST 20 16 0.0730 -0.0313 0.007629 0.225755 POST 30 14 0.0780 -0.0100 0.118896 0.973508 POST 40 11 0.0441 -0.0240 0.101562 0.655348 POST 50 11 0.0575 -0.0145 0.053711 0.901660 PRE 1 59 0.1092 -0.0158 0.000412 0.949116 PRE 2 41 0.1098 -0.0073 0.001388 0.963387 PRE 4 24 0.0981 0.0198 0.177980 0.602634 PRE 6 19 0.0875 0.0390 0.768005 0.260503 PRE 8 18 0.0932 0.0239 0.246208 0.088583 PRE 10 15 0.0774 0.0399 0.678772 0.038438 PRE 15 12 0.0979 0.0307 0.176270 0.456961 PRE 20 10 0.0949 0.0414 0.193359 0.225755 PRE 30 7 0.0513 0.0193 0.578125 0.973508 PRE 40 5 0.0342 0.0294 1.000000 0.655348 PRE 50 5 0.0220 -0.0085 0.812500 0.901660
factors where the correlation exceeds chance (p < 0.01): {'PRE': [1, 2], 'POST': [1, 2, 4, 20]}
factors where PRE and POST differ (p <= 0.01): []
Result: fig4d-compression not reproduced. PRE correlations exceed chance at factors 1 and 2 only, POST at 1, 2, 4 and 20. PRE and POST do not differ at any factor (smallest p = 0.038). The number of pairs with a defined correlation falls from 59 (PRE) and 77 (POST) at factor 1 to 5 and 11 at factor 50, as the HFO correlogram is empty in the narrow windows; at factor 4, PRE rests on 24 pairs (p = 0.18). Chance uses one circular shift per pair for all factors.
figs2-post-rate-increase¶
Firing rates increase during post-behavior NREM sleep in CA1, CA3 and LS.
In [23]:
# figs2-post-rate-increase
rows = []
for r in R + skipped:
m = r["meta"]
if m["states"] is None or m["epochs"] is None or m["units"]["location"] is None:
continue
st = m["states"]
nrem = st["label"] == "Non-REM"
def nrem_within(a, b):
s0, s1 = np.clip(st["start"][nrem], a, b), np.clip(st["stop"][nrem], a, b)
keep = s1 > s0
return s0[keep], s1[keep]
pre_iv = nrem_within(m["epochs"]["start"][0], m["epochs"]["stop"][0])
post_iv = nrem_within(m["epochs"]["start"][2], m["epochs"]["stop"][2])
dpre, dpost = (pre_iv[1] - pre_iv[0]).sum(), (post_iv[1] - post_iv[0]).sum()
if dpre < 60 or dpost < 60:
continue
o1, o2 = np.argsort(pre_iv[0]), np.argsort(post_iv[0])
for i, s in enumerate(m["units"]["spikes"]):
s = np.sort(s)
rows.append({"region": m["units"]["location"][i], "subject": m["subject"],
"pre": ls_hfo.counts_in(s, pre_iv[0][o1], pre_iv[1][o1]).sum() / dpre,
"post": ls_hfo.counts_in(s, post_iv[0][o2], post_iv[1][o2]).sum() / dpost})
rates_pp = pd.DataFrame(rows)
out = []
for g in ("hpc", "ca3", "ls"):
d = rates_pp[(rates_pp.region == g) & ((rates_pp.pre > 0) | (rates_pp.post > 0))]
w = stats.wilcoxon(d.post, d.pre)
out.append({"region": g, "units": len(d), "recordings": None, "median PRE (Hz)": round(d.pre.median(), 3),
"median POST (Hz)": round(d.post.median(), 3),
"median POST/PRE": round(np.median((d.post + 1e-3) / (d.pre + 1e-3)), 3), "Wilcoxon p": w.pvalue})
print(pd.DataFrame(out).drop(columns="recordings").to_string(index=False))
print("by animal, median POST/PRE:")
print(rates_pp[(rates_pp.pre > 0) | (rates_pp.post > 0)].assign(ratio=lambda d: (d.post + 1e-3) / (d.pre + 1e-3))
.groupby(["region", "subject"]).ratio.agg(["count", "median"]).round(3).to_string())
region units median PRE (Hz) median POST (Hz) median POST/PRE Wilcoxon p
hpc 2552 0.228 0.227 0.998 5.212872e-01
ca3 466 0.233 0.232 1.030 2.224666e-01
ls 1326 0.618 0.642 1.060 8.355372e-07
by animal, median POST/PRE:
count median
region subject
acg DT7 6 1.050
ca3 DT5 236 1.075
DT9 230 0.973
cc DT7 1 1.509
hpc DT2 1024 0.994
DT5 284 1.027
DT7 424 0.881
DT9 820 1.027
ls DT2 395 1.027
DT5 402 1.101
DT7 51 1.504
DT9 478 1.028
ms DT7 25 0.755
ppc DT2 28 0.778
sFi DT7 54 0.973
Result: partially reproduced. LS units fire more in POST NREM than in PRE NREM (median ratio
1.06, Wilcoxon p = 8.4e-07, 1326 units). The hpc units (ratio 1.00, p = 0.52) and the ca3 units
(ratio 1.03, p = 0.22) do not.
fig5b-nrem-and-rates¶
LS HFOs occur primarily during NREM sleep, and their rates are lower than those of hippocampal SPW-Rs.
In [24]:
# fig5b-nrem-and-rates
rows = []
for r in R:
st = r["meta"]["states"]
if st is None:
continue
row = {"file": r["path"].split("/")[1][:30], "subject": r["meta"]["subject"],
"LS NREM": rate_in(r["ls_events"]["peak"], st, "Non-REM"), "LS WAKE": rate_in(r["ls_events"]["peak"], st, "Awake")}
if "hpc_events" in r:
row["CA1 NREM"] = rate_in(r["hpc_events"]["peak"], st, "Non-REM")
row["CA1 WAKE"] = rate_in(r["hpc_events"]["peak"], st, "Awake")
rows.append(row)
rt = pd.DataFrame(rows)
no_state = rt[rt[["LS NREM", "LS WAKE"]].isna().any(axis=1)]
print(f"recordings without NREM or without WAKE time, left out: {len(no_state)} {no_state.file.tolist()}")
rt = rt.dropna(subset=["LS NREM", "LS WAKE"])
w1 = stats.wilcoxon(rt["LS NREM"], rt["LS WAKE"])
print(f"{len(rt)} recordings with states. LS HFO rate (events/s): NREM median {rt['LS NREM'].median():.3f}, WAKE median {rt['LS WAKE'].median():.3f};"
f" NREM > WAKE in {(rt['LS NREM'] > rt['LS WAKE']).sum()} recordings; NREM < WAKE in {(rt['LS NREM'] < rt['LS WAKE']).sum()}; median NREM - WAKE {np.median(rt['LS NREM'] - rt['LS WAKE']):.4f}; Wilcoxon p = {w1.pvalue:.3g}")
b = rt.dropna(subset=["CA1 NREM"])
w2 = stats.wilcoxon(b["LS NREM"], b["CA1 NREM"])
print(f"{len(b)} recordings with both: NREM rate LS median {b['LS NREM'].median():.3f}, CA1 median {b['CA1 NREM'].median():.3f};"
f" LS < CA1 in {(b['LS NREM'] < b['CA1 NREM']).sum()}; Wilcoxon p = {w2.pvalue:.3g}")
print(f"WAKE rates: LS median {b['LS WAKE'].median():.3f}, CA1 median {b['CA1 WAKE'].median():.3f}")
print("by animal (median rates):")
print(rt.groupby("subject")[["LS NREM", "LS WAKE", "CA1 NREM", "CA1 WAKE"]].median().round(3).to_string())
fig, ax = plt.subplots(figsize=(4, 3))
for _, x in rt.iterrows():
ax.plot([0, 1], [x["LS WAKE"], x["LS NREM"]], "m-", lw=0.5)
for _, x in b.iterrows():
ax.plot([2, 3], [x["CA1 WAKE"], x["CA1 NREM"]], "k-", lw=0.5)
ax.set_xticks([0, 1, 2, 3], ["LS wake", "LS NREM", "CA1 wake", "CA1 NREM"]); ax.set_ylabel("events/s"); ax.set_ylim(0, 2)
plt.tight_layout(); plt.show()
print("rates above 2/s (outside the axis):", int((rt[["LS NREM", "LS WAKE"]] > 2).sum().sum() + (b[["CA1 NREM", "CA1 WAKE"]] > 2).sum().sum()))
recordings without NREM or without WAKE time, left out: 4 ['sub-DT5_ses-20161016-1368um-14', 'sub-DT5_ses-20161023-1872um-24', 'sub-DT5_ses-20161027-1944um-22', 'sub-DT9_ses-20170430-0um-0um-m']
49 recordings with states. LS HFO rate (events/s): NREM median 0.006, WAKE median 0.010; NREM > WAKE in 18 recordings; NREM < WAKE in 29; median NREM - WAKE -0.0016; Wilcoxon p = 0.44
39 recordings with both: NREM rate LS median 0.006, CA1 median 0.461; LS < CA1 in 34; Wilcoxon p = 3.07e-06
WAKE rates: LS median 0.009, CA1 median 0.036
by animal (median rates):
LS NREM LS WAKE CA1 NREM CA1 WAKE
subject
DT2 0.000 0.005 0.407 0.050
DT5 0.006 0.003 0.285 0.037
DT7 0.015 0.059 0.624 0.062
DT9 0.009 0.026 0.461 0.017
rates above 2/s (outside the axis): 0
Result: partially reproduced. LS HFO rates in NREM are far below CA1 ripple rates (medians 0.006 and 0.461 per s, lower in 34 of 39 recordings, Wilcoxon p = 3.1e-06). The detected LS events are not more frequent in NREM than in WAKE: of 49 recordings with both states, NREM is higher in 18 and lower in 29 (Wilcoxon p = 0.44). On the re-referenced channel the same holds: NREM higher in 17 and lower in 31 (p = 0.73), and LS below CA1 in 35 of 39 (p = 6.7e-06).
fig5d-lag and fig5e-state-coupling¶
LS HFOs are reliably coupled to CA1 ripples, with a median peak-to-peak lag of 6.4 +- 4.2 ms, CA1 leading LS.
HPC-LS coupling is stronger in NREM sleep than in waking; the fraction of LS HFOs within +-25 ms of a CA1 ripple is 0.3 in NREM and 0.15 in WAKE (P < 6.6e-6).
In [25]:
# fig5d-lag and fig5e-state-coupling
def nearest(t, ref):
ref = np.sort(ref)
if len(ref) == 0:
return np.full(len(t), np.nan)
k = np.clip(np.searchsorted(ref, t), 1, len(ref) - 1) if len(ref) > 1 else np.zeros(len(t), int)
lo, hi = ref[np.maximum(k - 1, 0)], ref[k]
return np.where(np.abs(t - lo) <= np.abs(t - hi), t - lo, t - hi) # t minus the nearest ref
lag_rows, all_lags, ccg_sum = [], [], 0
for r in R:
if "hpc_events" not in r:
continue
d = nearest(r["ls_events"]["peak"], r["hpc_events"]["peak"])
d = d[np.abs(d) <= 0.05] * 1000
all_lags.append(d)
if len(d):
lag_rows.append({"file": r["path"].split("/")[1][:30], "subject": r["meta"]["subject"], "pairs": len(d), "median lag": np.median(d)})
ccg_sum = ccg_sum + ls_hfo.ccg(np.sort(r["hpc_events"]["peak"]), np.sort(r["ls_events"]["peak"]), half=0.2, bin_s=0.005)
lags = pd.DataFrame(lag_rows)
pooled = np.concatenate(all_lags)
print(f"LS HFO minus nearest CA1 ripple within 50 ms: (a) median of {len(lags)} recording medians {lags['median lag'].median():.2f} ms"
f" (SD {lags['median lag'].std():.2f}); (b) pooled median {np.median(pooled):.2f} ms over {len(pooled)} pairs [paper 6.4 +- 4.2, CA1 leading]")
print(f"fraction of pairs with LS after CA1: {np.mean(pooled > 0):.3f}; at the same sample (|lag| < 0.5 ms): {np.mean(np.abs(pooled) < 0.5):.3f}")
nz = np.concatenate([d[np.abs(d) >= 1] for d in all_lags])
print(f"variant without pairs within 1 ms: pooled median {np.median(nz):.2f} ms over {len(nz)} pairs; LS after CA1 in {np.mean(nz > 0):.3f}")
print("fraction of pairs at the same sample, by animal:", {s_: round(float(np.mean(np.abs(np.concatenate([d for d, r in zip(all_lags, [r for r in R if 'hpc_events' in r]) if r['meta']['subject'] == s_])) < 0.5)), 3) for s_ in sorted({r['meta']['subject'] for r in R if 'hpc_events' in r})})
lagc = np.arange(-200, 201, 5)
flank = ccg_sum[np.abs(lagc) >= 100]
centre = ccg_sum[np.abs(lagc) <= 25]
print(f"pooled CA1-LS CCG (5 ms bins): peak within +-25 ms {centre.max():.0f} at {lagc[np.abs(lagc) <= 25][centre.argmax()]} ms;"
f" flanks (100-200 ms) mean {flank.mean():.1f}, SD {flank.std():.1f}; peak - flank = {(centre.max() - flank.mean()) / flank.std():.1f} SD")
print("by animal, median of recording medians, all pairs (ms):", lags.groupby("subject")["median lag"].median().round(2).to_dict())
subj_hpc = [r["meta"]["subject"] for r in R if "hpc_events" in r]
print("by animal, pooled median without pairs within 1 ms (ms):",
{s_: round(float(np.median(np.concatenate([d[np.abs(d) >= 1] for d, sj in zip(all_lags, subj_hpc) if sj == s_]))), 2) for s_ in sorted(set(subj_hpc))})
fig, axes = plt.subplots(1, 2, figsize=(9, 3))
axes[0].bar(lagc, ccg_sum, width=5, color="k"); axes[0].set_xlabel("LS HFO time - CA1 ripple time (ms)"); axes[0].set_ylabel("count")
axes[1].hist(pooled, bins=50, range=(-50, 50), color="m"); axes[1].set_xlabel("lag to nearest ripple (ms)")
plt.tight_layout(); plt.show()
rows = []
for r in R:
st = r["meta"]["states"]
if st is None or "hpc_events" not in r:
continue
e = r["ls_events"]
d = np.abs(nearest(e["peak"], r["hpc_events"]["peak"]))
row = {"file": r["path"].split("/")[1][:30], "subject": r["meta"]["subject"]}
for lab, name in (("Non-REM", "NREM"), ("Awake", "WAKE")):
sel = e["state"] == lab
row[name + " n"] = int(sel.sum())
row[name] = np.mean(d[sel] <= 0.025) if sel.sum() else np.nan
row[name + " coupled"] = int((d[sel] <= 0.025).sum())
rows.append(row)
cp = pd.DataFrame(rows)
cpk = cp[(cp["NREM n"] >= 10) & (cp["WAKE n"] >= 10)]
t = stats.ttest_ind(cpk.NREM, cpk.WAKE)
print(f"\n{len(cpk)} of {len(cp)} recordings with at least 10 LS HFOs in each state")
print(f"(a) mean over recordings: NREM {cpk.NREM.mean():.3f}, WAKE {cpk.WAKE.mean():.3f}; two-sample t-test p = {t.pvalue:.3g}; paired t-test p = {stats.ttest_rel(cpk.NREM, cpk.WAKE).pvalue:.3g} (variant)")
print(f"(b) pooled: NREM {cpk['NREM coupled'].sum() / cpk['NREM n'].sum():.3f}, WAKE {cpk['WAKE coupled'].sum() / cpk['WAKE n'].sum():.3f} [paper 0.3 and 0.15]")
# chance: the same fraction with the CA1 ripple times shifted by 1 s (keeps the ripple rate)
ch = []
for r in R:
st = r["meta"]["states"]
if st is None or "hpc_events" not in r:
continue
e = r["ls_events"]
d = np.abs(nearest(e["peak"], r["hpc_events"]["peak"] + 1.0))
for lab in ("Non-REM", "Awake"):
sel = e["state"] == lab
ch.append({"state": lab, "frac": np.mean(d[sel] <= 0.025) if sel.sum() >= 10 else np.nan})
print("reference, ripple times shifted by 1 s (mean over recordings):", pd.DataFrame(ch).groupby("state").frac.mean().round(3).to_dict())
print(cpk.groupby("subject")[["NREM", "WAKE"]].mean().round(3).to_string())
LS HFO minus nearest CA1 ripple within 50 ms: (a) median of 70 recording medians 0.00 ms (SD 2.35); (b) pooled median 0.00 ms over 5903 pairs [paper 6.4 +- 4.2, CA1 leading]
fraction of pairs with LS after CA1: 0.335; at the same sample (|lag| < 0.5 ms): 0.348
variant without pairs within 1 ms: pooled median 3.20 ms over 3168 pairs; LS after CA1 in 0.554
fraction of pairs at the same sample, by animal: {'DT2': 0.619, 'DT5': 0.387, 'DT7': 0.352, 'DT8': 0.711, 'DT9': 0.073}
pooled CA1-LS CCG (5 ms bins): peak within +-25 ms 2966 at 0 ms; flanks (100-200 ms) mean 95.9, SD 17.3; peak - flank = 166.3 SD
by animal, median of recording medians, all pairs (ms): {'DT2': 0.0, 'DT5': 0.0, 'DT7': 0.0, 'DT8': 0.0, 'DT9': 3.8}
by animal, pooled median without pairs within 1 ms (ms): {'DT2': -3.2, 'DT5': -5.6, 'DT7': -6.4, 'DT8': 1.6, 'DT9': 8.0}
25 of 43 recordings with at least 10 LS HFOs in each state
(a) mean over recordings: NREM 0.199, WAKE 0.116; two-sample t-test p = 0.0885; paired t-test p = 0.0957 (variant)
(b) pooled: NREM 0.088, WAKE 0.041 [paper 0.3 and 0.15]
reference, ripple times shifted by 1 s (mean over recordings): {'Awake': 0.005, 'Non-REM': 0.019}
NREM WAKE
subject
DT2 0.100 0.108
DT5 0.082 0.327
DT7 0.101 0.049
DT9 0.382 0.093
Result: fig5d-lag not reproduced. LS HFOs are coupled to CA1 ripples: the pooled correlogram peaks at 0 ms, 166 SD above its flanks. The lag is 0.00 ms under both readings, outside 6.4 +- 4.2 ms, and LS follows CA1 in only 33.5% of pairs; 34.8% of the coupled pairs fall on the same sample (71% in DT8, 7% in DT9). Without the pairs within 1 ms the pooled median is 3.2 ms, LS after CA1 in 55%, but by animal it is -3.2 (DT2), -5.6 (DT5), -6.4 (DT7), 1.6 (DT8) and 8.0 ms (DT9), so LS leads in three animals. Over all pairs, DT9's median of recording medians is 3.8 ms, inside the tolerance (2.2 to 10.6 ms); the other four animals give 0.0 ms. On the re-referenced channel, where 4.1% of pairs fall on the same sample, the coupling remains (53.4 SD at 0 ms) and the lag is 0.00 ms (a) and 0.80 ms (b), LS after CA1 in 52%, still outside the tolerance. The coupling part of the claim holds and the lag part does not. The criterion requires the lag within tolerance under one reading and defines no partial outcome, so the claim is not reproduced.
Result: fig5e-state-coupling not reproduced. Over 25 recordings, the fraction of LS HFOs within 25 ms of a CA1 ripple is 0.199 in NREM and 0.116 in WAKE (two-sample t-test p = 0.089); pooled, 0.088 and 0.041. The direction matches the paper, but the difference is not significant and the values are below 0.3 and 0.15. Chance (ripple times shifted by 1 s) is 0.019 and 0.005. On the re-referenced channel, over 16 recordings, 0.118 and 0.085 (p = 0.43).
fig5f-ripple-triggered-firing and fig5f-no-phase-coherence¶
LS firing around CA1 ripples is higher during NREM than during waking.
CA1-LS phase coherence in the 120-180 Hz band around CA1 ripples does not differ from that of randomly jittered windows, so LS does not inherit the ripple wave by wave.
In [26]:
# fig5f-ripple-triggered-firing and fig5f-no-phase-coherence
rows = []
for r in R:
st = r["meta"]["states"]
if st is None or "hpc_events" not in r:
continue
e = r["hpc_events"]
tn, tw = e["peak"][e["state"] == "Non-REM"], e["peak"][e["state"] == "Awake"]
if len(tn) < 20 or len(tw) < 20:
continue
for i, s in ls_hfo.sorted_spikes(r, "ls").items():
rn = ls_hfo.counts_in(s, tn - 0.025, tn + 0.025).mean() / 0.05
rw = ls_hfo.counts_in(s, tw - 0.025, tw + 0.025).mean() / 0.05
base_n = rate_in(s, st, "Non-REM")
base_w = rate_in(s, st, "Awake")
rows.append({"subject": r["meta"]["subject"], "NREM": rn, "WAKE": rw, "NREM gain": rn / base_n if base_n > 0 else np.nan,
"WAKE gain": rw / base_w if base_w > 0 else np.nan})
rf = pd.DataFrame(rows)
t = stats.ttest_rel(rf.NREM, rf.WAKE)
print(f"{len(rf)} LS neurons (paper 1,313). Rate within +-25 ms of CA1 ripples: NREM mean {rf.NREM.mean():.2f} Hz, WAKE mean {rf.WAKE.mean():.2f} Hz;"
f" NREM higher for {(rf.NREM > rf.WAKE).sum()} neurons; paired t-test p = {t.pvalue:.3g}")
g = rf.dropna(subset=["NREM gain", "WAKE gain"])
t2 = stats.ttest_rel(g["NREM gain"], g["WAKE gain"])
print(f"variant, rate divided by the neuron's mean rate in that state ({len(g)} neurons): NREM {g['NREM gain'].mean():.2f}, WAKE {g['WAKE gain'].mean():.2f}, paired t-test p = {t2.pvalue:.3g}")
print(rf.groupby("subject")[["NREM", "WAKE"]].agg(["count", "mean"]).round(2).to_string())
coh = [r["coherence"] for r in R if "coherence" in r and r["coherence"]["n"] > 0]
plv = np.array([c["plv"] for c in coh]); jit = np.array([c["jitter"] for c in coh]) # recordings x lags, recordings x 20 x lags
lagsc = coh[0]["lags"] * 1000
mid = len(lagsc) // 2
avg, javg = plv.mean(axis=0), jit.mean(axis=0) # lags; 20 x lags
mu, sd = javg.mean(axis=0), javg.std(axis=0)
print(f"\n{len(coh)} recordings. CA1-LS 120-180 Hz phase-locking value at the ripple time: {avg[mid]:.4f};"
f" jittered {mu[mid]:.4f} +- {sd[mid]:.4f} (SD of 20 sets); excess {(avg[mid] - mu[mid]) / sd[mid]:.1f} SD; bound mean + 3 SD = {mu[mid] + 3 * sd[mid]:.4f}")
print(f"lags within +-50 ms above the bound: {int((avg[np.abs(lagsc) <= 50] > (mu + 3 * sd)[np.abs(lagsc) <= 50]).sum())} of {int((np.abs(lagsc) <= 50).sum())}")
n_rip = np.array([c["n"] for c in coh])
print(f"ripples per recording in the coherence: median {np.median(n_rip):.0f}, range {n_rip.min()} to {n_rip.max()};"
f" a PLV of n random phase differences would be about {np.median(np.sqrt(np.pi / (4 * n_rip))):.3f} (median over recordings)")
print(f"jittered PLV at the ripple time, per recording: median {np.median(jit[:, :, mid].mean(axis=1)):.3f}, range {jit[:, :, mid].mean(axis=1).min():.3f} to {jit[:, :, mid].mean(axis=1).max():.3f}")
per = plv[:, mid] - jit[:, :, mid].mean(axis=1)
print(f"per recording, PLV minus mean jittered PLV: median {np.median(per):.4f}, above 0 in {(per > 0).sum()} of {len(per)}")
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(lagsc, avg, "m", label="CA1 ripples")
ax.fill_between(lagsc, mu - 3 * sd, mu + 3 * sd, color="k", alpha=0.3, label="jittered +-3 SD")
ax.set_xlabel("ms from CA1 ripple"); ax.set_ylabel("phase-locking value"); ax.legend(); plt.tight_layout(); plt.show()
print(f"plotted range: {min(avg.min(), (mu - 3 * sd).min()):.4f} to {max(avg.max(), (mu + 3 * sd).max()):.4f}")
978 LS neurons (paper 1,313). Rate within +-25 ms of CA1 ripples: NREM mean 1.39 Hz, WAKE mean 1.33 Hz; NREM higher for 648 neurons; paired t-test p = 0.238
variant, rate divided by the neuron's mean rate in that state (977 neurons): NREM 1.34, WAKE 1.56, paired t-test p = 0.000156
NREM WAKE
count mean count mean
subject
DT2 331 1.15 331 1.28
DT5 210 2.42 210 2.28
DT7 51 0.73 51 0.49
DT9 386 1.11 386 0.95
71 recordings. CA1-LS 120-180 Hz phase-locking value at the ripple time: 0.4051; jittered 0.2018 +- 0.0066 (SD of 20 sets); excess 30.7 SD; bound mean + 3 SD = 0.2217
lags within +-50 ms above the bound: 90 of 125
ripples per recording in the coherence: median 1722, range 8 to 19577; a PLV of n random phase differences would be about 0.021 (median over recordings)
jittered PLV at the ripple time, per recording: median 0.134, range 0.030 to 0.602
per recording, PLV minus mean jittered PLV: median 0.1277, above 0 in 70 of 71
plotted range: 0.1749 to 0.4089
Result: fig5f-ripple-triggered-firing not reproduced. Over 978 LS neurons, the rate within 25 ms of CA1 ripples is 1.39 Hz in NREM and 1.33 Hz in WAKE (paired t-test p = 0.24). Divided by each neuron's rate in the same state, WAKE is higher (1.56 against 1.34, p = 0.00016).
Result: fig5f-no-phase-coherence not reproduced. The CA1-LS phase-locking value at the ripple time is 0.405, against 0.202 +- 0.007 for jittered windows (31 SD above), and 90 of 125 lags within +-50 ms exceed the bound. It is higher at ripples in 70 of 71 recordings. The jittered value itself is high: with a median of 1722 ripples per recording, random phase differences would give about 0.021, and the jittered value per recording is 0.030 to 0.602 (median 0.134). The two channels keep a phase relation outside ripples too, a signal common to both. On the re-referenced LS channel the value at ripples is 0.164 against 0.074 +- 0.004 jittered (23.9 SD above, higher in 60 of 71 recordings): re-referencing on the shank lowers both, and the coherence at ripples stays above chance.
Check: LS events on a re-referenced channel¶
Several outputs above point to a 120-200 Hz signal shared by the LS detection channel and the hippocampal channel: a third of the coupled LS events fall on the same sample as the CA1 ripple, and the CA1-LS phase-locking value at jittered times is far above that of random phases. A signal common to every channel, such as the reference, would do this. Here the detection channel is re-referenced to the LS channel on its own shank with the lowest 120-200 Hz SD, which removes what is common to the shank, and events are detected again with the same parameters. The measures of the HFO claims are repeated on these events. This check is not in the paper's Methods and does not decide an outcome; it shows how much of each result rests on the shared signal.
In [27]:
# Re-referenced LS channel: the measures of the HFO claims on events detected on it.
RB = [r for r in R if "ls_events_bip" in r]
print(f"recordings with a second good LS channel on the detection shank: {len(RB)} of {len(R)}")
print("LS HFOs, as detected / re-referenced, by animal:",
{s_: (sum(len(r["ls_events"]["peak"]) for r in RB if r["meta"]["subject"] == s_),
sum(len(r["ls_events_bip"]["peak"]) for r in RB if r["meta"]["subject"] == s_)) for s_ in animals})
print("re-referenced events in recordings with LS units:", sum(len(r["ls_events_bip"]["peak"]) for r in RB if has_units(r, "ls")))
# fig1c-hfo-each-animal: firing within +-10 ms of the events against random times
rng = np.random.default_rng(1191)
rows = []
for subj in animals:
d = []
for r in [r for r in RB if r["meta"]["subject"] == subj]:
t = r["ls_events_bip"]["peak"]
units = ls_hfo.sorted_spikes(r, "ls")
if len(t) == 0 or not units:
continue
rnd_t = rng.uniform(0.2, r["meta"]["duration"] - 0.2, len(t))
d += [ls_hfo.peth(s, t)[90:110].sum() - ls_hfo.peth(s, rnd_t)[90:110].sum() for s in units.values()]
d = np.array(d)
f_ = np.concatenate([r["ls_events_bip"]["frequency"] for r in RB if r["meta"]["subject"] == subj] or [np.array([np.nan])])
rows.append({"animal": subj, "neurons": len(d), "median freq (Hz)": round(np.nanmedian(f_), 1) if np.isfinite(f_).any() else np.nan,
"difference": round(d.mean(), 4) if len(d) else np.nan,
"difference / SEM": round(d.mean() / (d.std(ddof=1) / np.sqrt(len(d))), 1) if len(d) > 2 else np.nan})
print("fig1c-hfo-each-animal, re-referenced:")
print(pd.DataFrame(rows).to_string(index=False))
# figs1-participation
fr, sp_ = [], []
for r in RB:
units = ls_hfo.sorted_spikes(r, "ls")
e = r["ls_events_bip"]
if units and len(e["peak"]):
C = np.array([ls_hfo.counts_in(s, e["start"], e["stop"]) for s in units.values()])
fr.append((C > 0).mean(axis=0)); sp_.append(C.ravel())
print(f"figs1-participation, re-referenced, pooled over {sum(len(x) for x in fr)} events: fraction {np.concatenate(fr).mean():.3f}, spikes {np.concatenate(sp_).mean():.3f}")
# fig1f and fig1g: duration and frequency against CA1 ripples
bb = [r for r in RB if "hpc_events" in r]
a_d = np.concatenate([r["ls_events_bip"]["duration"] for r in bb]) * 1000
a_f = np.concatenate([r["ls_events_bip"]["frequency"] for r in bb]); a_f = a_f[np.isfinite(a_f)]
c_d = np.concatenate([r["hpc_events"]["duration"] for r in bb]) * 1000
c_f = np.concatenate([r["hpc_events"]["frequency"] for r in bb]); c_f = c_f[np.isfinite(c_f)]
print(f"fig1f/fig1g, re-referenced, {len(bb)} recordings: duration LS {a_d.mean():.1f} against CA1 {c_d.mean():.1f} ms (p = {stats.mannwhitneyu(a_d, c_d).pvalue:.2g});"
f" frequency LS {a_f.mean():.1f} against CA1 {c_f.mean():.1f} Hz (p = {stats.mannwhitneyu(a_f, c_f).pvalue:.2g})")
# fig5b-nrem-and-rates
rows = []
for r in RB:
st = r["meta"]["states"]
if st is None:
continue
row = {"NREM": rate_in(r["ls_events_bip"]["peak"], st, "Non-REM"), "WAKE": rate_in(r["ls_events_bip"]["peak"], st, "Awake"),
"CA1 NREM": rate_in(r["hpc_events"]["peak"], st, "Non-REM") if "hpc_events" in r else np.nan}
rows.append(row)
rb = pd.DataFrame(rows).dropna(subset=["NREM", "WAKE"])
w = stats.wilcoxon(rb.NREM, rb.WAKE)
b2 = rb.dropna(subset=["CA1 NREM"])
print(f"fig5b, re-referenced, {len(rb)} recordings: LS rate NREM median {rb.NREM.median():.4f}/s, WAKE {rb.WAKE.median():.4f}/s,"
f" NREM higher in {(rb.NREM > rb.WAKE).sum()}, lower in {(rb.NREM < rb.WAKE).sum()}, Wilcoxon p = {w.pvalue:.3g};"
f" below the CA1 ripple rate in {(b2.NREM < b2['CA1 NREM']).sum()} of {len(b2)} (p = {stats.wilcoxon(b2.NREM, b2['CA1 NREM']).pvalue:.3g})")
# fig5d-lag
lag_rows, all_l, ccg_b = [], [], 0
for r in bb:
d = nearest(r["ls_events_bip"]["peak"], r["hpc_events"]["peak"])
d = d[np.abs(d) <= 0.05] * 1000
all_l.append(d)
if len(d):
lag_rows.append({"subject": r["meta"]["subject"], "pairs": len(d), "median lag": np.median(d)})
ccg_b = ccg_b + ls_hfo.ccg(np.sort(r["hpc_events"]["peak"]), np.sort(r["ls_events_bip"]["peak"]), half=0.2, bin_s=0.005)
lb = pd.DataFrame(lag_rows); pooled_b = np.concatenate(all_l)
fl = ccg_b[np.abs(lagc) >= 100]; ce = ccg_b[np.abs(lagc) <= 25]
print(f"fig5d, re-referenced: (a) median of {len(lb)} recording medians {lb['median lag'].median():.2f} ms; (b) pooled median {np.median(pooled_b):.2f} ms over {len(pooled_b)} pairs;"
f" LS after CA1 {np.mean(pooled_b > 0):.3f}; same sample {np.mean(np.abs(pooled_b) < 0.5):.3f};"
f" CCG peak within +-25 ms at {lagc[np.abs(lagc) <= 25][ce.argmax()]} ms, {(ce.max() - fl.mean()) / fl.std():.1f} SD above the flanks")
print(" pairs by animal:", lb.groupby("subject").pairs.sum().to_dict(), "| median lag by animal (ms):", lb.groupby("subject")["median lag"].median().round(2).to_dict())
# fig5e-state-coupling
rows = []
for r in bb:
if r["meta"]["states"] is None:
continue
e = r["ls_events_bip"]
d = np.abs(nearest(e["peak"], r["hpc_events"]["peak"]))
row = {}
for lab, name in (("Non-REM", "NREM"), ("Awake", "WAKE")):
sel = e["state"] == lab
row[name + " n"], row[name], row[name + " c"] = int(sel.sum()), np.mean(d[sel] <= 0.025) if sel.sum() else np.nan, int((d[sel] <= 0.025).sum())
rows.append(row)
cb = pd.DataFrame(rows)
cbk = cb[(cb["NREM n"] >= 10) & (cb["WAKE n"] >= 10)]
if len(cbk) >= 2:
print(f"fig5e, re-referenced, {len(cbk)} recordings with 10 events in each state: (a) NREM {cbk.NREM.mean():.3f}, WAKE {cbk.WAKE.mean():.3f},"
f" t-test p = {stats.ttest_ind(cbk.NREM, cbk.WAKE).pvalue:.3g}; (b) pooled NREM {cbk['NREM c'].sum() / cbk['NREM n'].sum():.3f}, WAKE {cbk['WAKE c'].sum() / cbk['WAKE n'].sum():.3f}")
else:
print(f"fig5e, re-referenced: {len(cbk)} recordings with 10 events in each state")
# fig5f-no-phase-coherence
cohb = [r["coherence_bip"] for r in RB if "coherence_bip" in r and r["coherence_bip"]["n"] > 0]
cohr = [r["coherence"] for r in RB if "coherence_bip" in r and r["coherence_bip"]["n"] > 0]
for label, cc in (("as detected, same recordings", cohr), ("re-referenced", cohb)):
P = np.array([c["plv"] for c in cc]); J = np.array([c["jitter"] for c in cc])
av, ja = P.mean(axis=0), J.mean(axis=0)
mu_, sd_ = ja.mean(axis=0), ja.std(axis=0)
per_ = P[:, mid] - J[:, :, mid].mean(axis=1)
print(f"fig5f coherence, {label}, {len(cc)} recordings: PLV at the ripple time {av[mid]:.4f}, jittered {mu_[mid]:.4f} +- {sd_[mid]:.4f},"
f" {(av[mid] - mu_[mid]) / sd_[mid]:.1f} SD; within mean +- 3 SD: {abs(av[mid] - mu_[mid]) <= 3 * sd_[mid]};"
f" lags within +-50 ms above the bound {int((av[np.abs(lagsc) <= 50] > (mu_ + 3 * sd_)[np.abs(lagsc) <= 50]).sum())} of {int((np.abs(lagsc) <= 50).sum())};"
f" higher than jittered in {(per_ > 0).sum()} recordings, median excess {np.median(per_):.4f}")
recordings with a second good LS channel on the detection shank: 94 of 94
LS HFOs, as detected / re-referenced, by animal: {'DT1': (5292, 8166), 'DT2': (1183, 3741), 'DT5': (4032, 34492), 'DT7': (38008, 47824), 'DT8': (2757, 1243), 'DT9': (15914, 14777)}
re-referenced events in recordings with LS units: 66458
fig1c-hfo-each-animal, re-referenced: animal neurons median freq (Hz) difference difference / SEM DT1 610 159.2 0.0052 1.3 DT2 396 156.3 0.0150 1.5 DT5 560 157.2 0.0226 1.3 DT7 60 155.4 0.0009 0.7 DT8 296 156.6 -0.0084 -4.4 DT9 510 157.2 -0.0032 -2.1
figs1-participation, re-referenced, pooled over 66458 events: fraction 0.045, spikes 0.060
fig1f/fig1g, re-referenced, 72 recordings: duration LS 38.2 against CA1 38.3 ms (p = 0); frequency LS 156.3 against CA1 150.8 Hz (p = 0)
fig5b, re-referenced, 49 recordings: LS rate NREM median 0.0012/s, WAKE 0.0021/s, NREM higher in 17, lower in 31, Wilcoxon p = 0.727; below the CA1 ripple rate in 35 of 39 (p = 6.69e-06)
fig5d, re-referenced: (a) median of 70 recording medians 0.00 ms; (b) pooled median 0.80 ms over 4509 pairs; LS after CA1 0.522; same sample 0.041; CCG peak within +-25 ms at 0 ms, 53.4 SD above the flanks
pairs by animal: {'DT2': 94, 'DT5': 352, 'DT7': 1616, 'DT8': 285, 'DT9': 2162} | median lag by animal (ms): {'DT2': 0.8, 'DT5': 0.8, 'DT7': 0.0, 'DT8': 0.0, 'DT9': 1.4}
fig5e, re-referenced, 16 recordings with 10 events in each state: (a) NREM 0.118, WAKE 0.085, t-test p = 0.426; (b) pooled NREM 0.092, WAKE 0.009
fig5f coherence, as detected, same recordings, 71 recordings: PLV at the ripple time 0.4051, jittered 0.2018 +- 0.0066, 30.7 SD; within mean +- 3 SD: False; lags within +-50 ms above the bound 90 of 125; higher than jittered in 70 recordings, median excess 0.1277
fig5f coherence, re-referenced, 71 recordings: PLV at the ripple time 0.1640, jittered 0.0742 +- 0.0038, 23.9 SD; within mean +- 3 SD: False; lags within +-50 ms above the bound 41 of 125; higher than jittered in 60 recordings, median excess 0.0291
In [28]:
# fig2-confined-to-ls with each channel re-referenced to the median of its shank, triggered by the
# re-referenced LS events (single-channel shanks give no value)
rows = []
for r in R:
if "localization_bip" not in r:
continue
L, loc = r["localization_bip"], r["meta"]["el_location"]
good = r["channel_stats"]["good"]
row = {"file": r["path"].split("/")[1][:34], "HFOs in block": L["n_events"]}
for g in sorted(set(loc)):
row[g] = round(float(np.nanmedian(L["ratio"][(loc == g) & good])), 2) if np.isfinite(L["ratio"][(loc == g) & good]).any() else np.nan
rows.append(row)
lb2 = pd.DataFrame(rows)
print(lb2.to_string(index=False))
u2 = lb2[lb2["HFOs in block"] >= 20]
oth = [c for c in lb2.columns if c not in ("file", "HFOs in block", "ls", "hpc", "ca3", "unknown")]
print(f"recordings with at least 20 re-referenced HFOs in the block: {len(u2)}; LS above every other non-hippocampal region in "
f"{int(u2.apply(lambda row: all(row['ls'] > row[c] for c in oth if pd.notna(row[c])), axis=1).sum()) if len(u2) else 0}")
file HFOs in block acg cc hpc ls ms sFi sub-DT7_ses-20170315-288um-0um-170 196 1.29 1.48 0.93 1.76 NaN NaN sub-DT7_ses-20170318-288um-0um-170 481 1.66 NaN 1.74 3.09 4.00 NaN sub-DT7_ses-20170320-432um-0um-mer 183 1.76 NaN 1.16 3.65 4.21 NaN sub-DT7_ses-20170321-216um-0um-170 0 NaN NaN NaN NaN NaN NaN sub-DT7_ses-20170324-576um-144um-1 2 8.59 NaN 1.20 14.34 16.72 NaN sub-DT7_ses-20170325-576um-360um-1 3 0.62 NaN 0.83 1.39 5.70 NaN sub-DT7_ses-20170401-864um-936um-m 2 NaN 1.35 0.26 1.66 0.95 NaN sub-DT7_ses-20170412-1440um-1170um 530 NaN NaN 1.07 2.83 NaN 2.0 recordings with at least 20 re-referenced HFOs in the block: 4; LS above every other non-hippocampal region in 2
On the re-referenced channel the share of coupled LS events on the same sample as the CA1 ripple falls from 34.8% to 4.1%, so most of those events came from a signal common to the shank. The other measures do not move towards the paper's: no animal shows LS firing at the events more than 3 SEM above random times, the lag stays at 0 to 0.8 ms, NREM rates are not above WAKE, and coherence at ripples stays above chance. Event counts change a lot (DT5 4,032 to 34,492; DT8 2,757 to 1,243), so the event set depends on the referencing, which the paper does not describe. No outcome is changed by this check.
fig6-optogenetic-hfos (blocked)¶
Optogenetic activation of LS GABAergic neurons with a Gaussian light waveform induces local oscillations of about 140 Hz.
The two rats used for optogenetics are not among the 6 subjects of the dandiset (see the data check), and no file holds an optogenetic series. The cell below searches the archive for other dandisets by the same authors.
In [29]:
# fig6-optogenetic-hfos (blocked): is there another dandiset of this paper on the archive?
for q in ("Tingley", "lateral septum"):
res_ = requests.get("https://api.dandiarchive.org/api/dandisets/", params={"search": q, "page_size": 100}, timeout=60)
res_.raise_for_status()
print(f"archive search {q!r}: {res_.json()['count']} dandisets")
for x in res_.json()["results"]:
v = x.get("most_recent_published_version") or x["draft_version"]
print(" ", x["identifier"], v["name"])
print("optogenetic series in the files:", sum("Optogenetic" in v["type"] for r in results for v in r["meta"]["acquisition"].values()),
"| subjects:", sorted({r["meta"]["subject"] for r in results}))
archive search 'Tingley': 3 dandisets 000213 Transformation of a Spatial Map across the Hippocampal-Lateral Septal Circuit 000218 Routing of Hippocampal Ripples to Subcortical Structures via the Lateral Septum 000233 A metabolic function of the hippocampal sharp wave-ripple
archive search 'lateral septum': 3 dandisets 000213 Transformation of a Spatial Map across the Hippocampal-Lateral Septal Circuit 000218 Routing of Hippocampal Ripples to Subcortical Structures via the Lateral Septum 000233 A metabolic function of the hippocampal sharp wave-ripple optogenetic series in the files: 0 | subjects: ['DT1', 'DT2', 'DT5', 'DT7', 'DT8', 'DT9']
Result: blocked (data missing). The 6 subjects are the recording rats, no file holds an optogenetic series, and archive searches for "Tingley" and "lateral septum" find only 000213, 000218 and 000233, none of which is described as holding these experiments.
fig7e-ca1-spikes, fig7e-state and fig7e-duration¶
The number of CA1 spikes in a SPW-R significantly predicts the number of spikes of 38% of LS neurons.
Brain state (WAKE or NREM) in which a SPW-R occurs significantly predicts the spike count of 26% of LS neurons.
The duration of hippocampal SPW-Rs significantly predicts the spike count of about 20% of LS neurons.
The model is fixed in the claims file (linear regression, 5-fold cross-validated MSE, 100 shuffles). Reading (a) calls a neuron significant when its MSE is below the 5th percentile of the shuffles, reading (b) when a t-test of the shuffles against it gives p < 0.05. The chance rate tests one shuffle against the others.
In [30]:
# fig7e-ca1-spikes, fig7e-state, fig7e-duration: does a ripple feature predict an LS neuron's spike count?
N_SHUF = 100
def fig7_rows(rng, fixed_window=None):
rows = []
for r in R:
if "hpc_events" not in r:
continue
ca1 = ls_hfo.sorted_spikes(r, "hpc")
lsu = ls_hfo.sorted_spikes(r, "ls")
e = r["hpc_events"]
if not ca1 or not lsu or len(e["peak"]) < 100:
continue
win = (e["start"], e["stop"]) if fixed_window is None else (e["peak"] - fixed_window, e["peak"] + fixed_window)
Y = np.array([ls_hfo.counts_in(s, *win) for s in lsu.values()]).T.astype(float) # ripples x neurons
keep_n = Y.sum(axis=0) >= 20
if not keep_n.any():
continue
Y = Y[:, keep_n]
n_ca1 = np.sum([ls_hfo.counts_in(s, e["start"], e["stop"]) for s in ca1.values()], axis=0).astype(float)
feats = {"ca1_spikes": (n_ca1, np.ones(len(n_ca1), bool)), "duration": (e["duration"], np.ones(len(n_ca1), bool)),
"state": ((e["state"] == "Non-REM").astype(float), np.isin(e["state"], ["Non-REM", "Awake"]))}
for name, (x, sel) in feats.items():
if sel.sum() < 100 or np.unique(x[sel]).size < 2:
continue
xs, ys = x[sel], Y[sel]
folds = rng.permutation(np.arange(len(xs)) % 5)
X = np.column_stack([xs] + [rng.permutation(xs) for _ in range(N_SHUF)])
mse = ls_hfo.cv_mse(ys, X, folds) # neurons x (1 + shuffles)
actual, shuf = mse[:, 0], mse[:, 1:]
for k in range(len(actual)):
s = shuf[k]
# chance: shuffle 0 tested against the other shuffles
c, cs = s[0], s[1:]
tt = stats.ttest_1samp(s, actual[k])
tc = stats.ttest_1samp(cs, c)
rows.append({"feature": name, "subject": r["meta"]["subject"], "file": r["path"].split("/")[1][:30],
"z": (s.mean() - actual[k]) / s.std(ddof=1) if s.std() > 0 else np.nan,
"sig_a": actual[k] < np.percentile(s, 5), "sig_b": (tt.pvalue < 0.05) and actual[k] < s.mean(),
"chance_a": c < np.percentile(cs, 5), "chance_b": (tc.pvalue < 0.05) and c < cs.mean()})
return pd.DataFrame(rows)
def fig7_summary(f7):
summary = []
for name, d in f7.groupby("feature"):
row = {"feature": name, "LS neurons": len(d), "recordings": d.file.nunique()}
for rd in ("a", "b"):
k, n, ch = d[f"sig_{rd}"].sum(), len(d), d[f"chance_{rd}"].mean()
row[f"% ({rd})"] = round(100 * k / n, 1)
row[f"chance % ({rd})"] = round(100 * ch, 1)
row[f"binomial p vs chance ({rd})"] = stats.binomtest(int(k), n, max(ch, 1e-6), alternative="greater").pvalue
summary.append(row)
return pd.DataFrame(summary)
f7 = fig7_rows(np.random.default_rng(118))
f7s = fig7_summary(f7)
print(f7s.to_string(index=False))
print("paper: CA1 spikes 38%, state 26%, duration about 20%")
print("by animal, % significant under reading (a):")
print(f7.groupby(["feature", "subject"]).sig_a.agg(["count", "mean"]).assign(mean=lambda d: (100 * d["mean"]).round(1)).to_string())
fig, ax = plt.subplots(figsize=(5, 3))
for j, (name, d) in enumerate(f7.groupby("feature")):
ax.plot(j + np.random.default_rng(1181).uniform(-0.2, 0.2, len(d)), d.z.clip(-5, 20), "k.", ms=2)
ax.set_xticks(range(f7.feature.nunique()), sorted(f7.feature.unique()))
ax.set_ylabel("SDs of shuffled MSE above actual (clipped -5 to 20)"); plt.tight_layout(); plt.show()
print("z outside the axis (clipped):", int(((f7.z < -5) | (f7.z > 20)).sum()), "of", len(f7))
# Variant: spikes counted in a fixed +-25 ms window around the ripple time, so the count does not grow with duration.
f7v = fig7_rows(np.random.default_rng(1182), fixed_window=0.025)
print("variant, fixed +-25 ms window:")
print(fig7_summary(f7v).to_string(index=False))
feature LS neurons recordings % (a) chance % (a) binomial p vs chance (a) % (b) chance % (b) binomial p vs chance (b)
ca1_spikes 880 27 79.4 4.7 0.000000e+00 88.4 40.5 1.009166e-193
duration 880 27 86.0 5.9 0.000000e+00 93.2 41.1 5.416100e-237
state 813 23 34.3 5.4 1.980314e-141 67.4 41.7 1.891534e-49
paper: CA1 spikes 38%, state 26%, duration about 20%
by animal, % significant under reading (a):
count mean
feature subject
ca1_spikes DT2 249 88.0
DT5 216 74.5
DT7 44 79.5
DT9 371 76.5
duration DT2 249 89.2
DT5 216 88.0
DT7 44 88.6
DT9 371 82.5
state DT2 249 43.0
DT5 158 34.8
DT7 44 63.6
DT9 362 24.6
z outside the axis (clipped): 616 of 2573
variant, fixed +-25 ms window:
feature LS neurons recordings % (a) chance % (a) binomial p vs chance (a) % (b) chance % (b) binomial p vs chance (b)
ca1_spikes 932 27 24.7 5.6 8.888078e-82 57.4 40.5 1.450791e-25
duration 932 27 17.8 7.2 5.603893e-27 51.4 42.4 1.897220e-08
state 862 23 34.9 6.0 2.953161e-142 65.5 42.6 4.866878e-42
Result: fig7e-ca1-spikes partially reproduced. The number of CA1 spikes predicts 79.4% of 880 LS neurons under reading (a) (chance 4.7%), above chance and far above the paper's 38%. Spikes are counted from ripple start to stop, so both the LS and the CA1 counts grow with the ripple's duration; with a fixed +-25 ms window, 24.7% under (a) (chance 5.6%). Reading (b) is not usable: its chance rate is 40.5% (88.4% significant), because a one-sample t-test of 99 shuffled MSEs against a single value is significant whenever that value lies a little away from their mean.
Result: fig7e-state partially reproduced. State predicts 34.3% of 813 LS neurons under (a) (chance 5.4%), above chance and above 26 +- 2.6%. With a fixed window, 34.9% (chance 6.0%). Under (b), 67.4% against a chance rate of 41.7%.
Result: fig7e-duration partially reproduced. Duration predicts 86.0% of 880 LS neurons under (a) (chance 5.9%), far above about 20%. Counting spikes from ripple start to stop makes the count grow with duration, which inflates this; with a fixed +-25 ms window it is 17.8% (chance 7.2%), just below the tolerance (18 to 22%). The fixed window is a variant and does not decide the outcome. Under (b), 93.2% against a chance rate of 41.1%.
fig7e-ca1-depth (blocked)¶
The mean depth within the CA1 pyramidal layer of the neurons active in a SPW-R predicts the spike count of 23% of LS neurons.
The depth of a CA1 neuron comes from its waveform amplitude across sites, and the units table holds
no waveforms or channels (see the data check). Waveforms could only be re-extracted from wide-band
data. The cell below lists the files that hold a series sampled above 5 kHz in acquisition, with
their size and the hippocampal units they hold.
In [31]:
# fig7e-ca1-depth (blocked): which files hold wide-band data from which waveforms could be re-extracted?
sizes = {a["path"]: a["size"] for a in asset_list}
rows = []
for r in results:
m = r["meta"]
raw = {k: v for k, v in m["acquisition"].items() if v["rate"] is not None and v["rate"] > 5000}
loc = m["units"]["location"]
rows.append({"subject": m["subject"], "file": r["path"].split("/")[1][:44], "wide-band series": ", ".join(raw) or None,
"rate (Hz)": max([v["rate"] for v in raw.values()], default=None),
"channels": max([v["shape"][1] if len(v["shape"]) > 1 else 1 for v in raw.values()], default=None),
"wide-band GB": sum(v["bytes"] for v in raw.values()) / 1e9, "file GB": sizes[r["path"]] / 1e9,
"hpc units": 0 if loc is None else int((loc == "hpc").sum())})
wb = pd.DataFrame(rows)
print("series types in acquisition:", dict(collections.Counter(v["type"] for r in results for v in r["meta"]["acquisition"].values())))
have = wb[wb["wide-band series"].notna()]
print(f"files with a wide-band series: {len(have)} of {len(wb)}; wide-band data {have['wide-band GB'].sum() / 1e3:.2f} TB,"
f" files {have['file GB'].sum() / 1e3:.2f} TB")
print(wb.groupby("subject").agg(files=("file", "size"), with_wide_band=("wide-band series", lambda x: int(x.notna().sum())),
hpc_units=("hpc units", "sum")).to_string())
print(have.drop(columns="file GB").round(1).to_string(index=False))
hv = have[have["hpc units"] > 0]
print(f"wide-band files with hpc units: {len(hv)}, {hv['wide-band GB'].sum() / 1e3:.2f} TB; streaming them within the"
f" 30 min run limit needs {hv['wide-band GB'].sum() * 1e3 / 1800:.0f} MB/s")
print(f"hpc units in files without wide-band data: {wb.loc[wb['wide-band series'].isna(), 'hpc units'].sum()} of {wb['hpc units'].sum()}")
series types in acquisition: {'ElectricalSeries': 27}
files with a wide-band series: 27 of 98; wide-band data 1.43 TB, files 1.15 TB
files with_wide_band hpc_units
subject
DT1 12 12 0
DT2 13 13 1035
DT5 20 0 431
DT7 20 0 446
DT8 17 0 0
DT9 16 2 866
subject file wide-band series rate (Hz) channels wide-band GB hpc units
DT1 sub-DT1_ses-DT1-rLS-20150723-1584um_ecephys. ElectricalSeries_raw 20000.0 30.0 13.1 0
DT1 sub-DT1_ses-DT1-rLS-20150723-1728um_ecephys. ElectricalSeries_raw 20000.0 30.0 9.5 0
DT1 sub-DT1_ses-DT1-rLS-20150724-1800um_ecephys. ElectricalSeries_raw 20000.0 30.0 11.3 0
DT1 sub-DT1_ses-DT1-rLS-20150724-1872um_ecephys. ElectricalSeries_raw 20000.0 30.0 15.2 0
DT1 sub-DT1_ses-DT1-rLS-20150725-1944um_ecephys. ElectricalSeries_raw 20000.0 30.0 20.2 0
DT1 sub-DT1_ses-DT1-rLS-20150726-2016um_ecephys. ElectricalSeries_raw 20000.0 30.0 15.2 0
DT1 sub-DT1_ses-DT1-rLS-20150727-2088um_ecephys. ElectricalSeries_raw 20000.0 30.0 13.3 0
DT1 sub-DT1_ses-DT1-rLS-20150727-2232um_ecephys. ElectricalSeries_raw 20000.0 30.0 15.1 0
DT1 sub-DT1_ses-DT1-rLS-20150728-2304um_ecephys. ElectricalSeries_raw 20000.0 30.0 23.6 0
DT1 sub-DT1_ses-DT1-rLS-20150729-2448um_ecephys. ElectricalSeries_raw 20000.0 30.0 12.6 0
DT1 sub-DT1_ses-DT1-rLS-20150729-2592um_ecephys. ElectricalSeries_raw 20000.0 30.0 18.3 0
DT1 sub-DT1_ses-DT1-rLS-20150730-2736um_ecephys. ElectricalSeries_raw 20000.0 30.0 13.5 0
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3036um-362um-20160 ElectricalSeries_raw 20000.0 128.0 67.3 0
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3108um-434um-20160 ElectricalSeries_raw 20000.0 128.0 48.8 0
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3324um-1000um-2016 ElectricalSeries_raw 20000.0 128.0 133.7 0
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3324um-1144um-2016 ElectricalSeries_raw 20000.0 128.0 111.5 67
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3396um-1180um-2016 ElectricalSeries_raw 20000.0 128.0 69.7 82
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3396um-1180um-2016 ElectricalSeries_raw 20000.0 128.0 108.2 112
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3468um-1216um-2016 ElectricalSeries_raw 20000.0 128.0 78.6 154
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3468um-1216um-2016 ElectricalSeries_raw 20000.0 128.0 45.8 148
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3468um-1216um-2016 ElectricalSeries_raw 20000.0 128.0 87.8 159
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3468um-1252um-2016 ElectricalSeries_raw 20000.0 128.0 121.1 99
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3540um-1288um-2016 ElectricalSeries_raw 20000.0 128.0 86.0 126
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3540um-1288um-2016 ElectricalSeries_raw 20000.0 128.0 58.0 79
DT2 sub-DT2_ses-DT2-rPPC-rCCG-3612um-1360um-2016 ElectricalSeries_raw 20000.0 128.0 8.4 9
DT9 sub-DT9_ses-20170501-216um-0um-170501-112210 ElectricalSeries_raw 20000.0 96.0 124.7 90
DT9 sub-DT9_ses-20170509-468um-36um-170509-10345 ElectricalSeries_raw 20000.0 96.0 101.5 74
wide-band files with hpc units: 12, 1.00 TB; streaming them within the 30 min run limit needs 556 MB/s
hpc units in files without wide-band data: 1579 of 2778
Result: blocked (data missing). 27 of 98 files hold a 20 kHz wide-band series: all 12 DT1 files, all 13 DT2 files and 2 DT9 files, which take 1.15 TB. Twelve of them (10 DT2, 2 DT9) hold hippocampal units, 1.00 TB of wide-band data, which would have to be streamed at about 556 MB/s to fit the 30 min run limit; for these recordings the claim is blocked by run time. The other 1579 of 2778 hippocampal units, including all of DT5 and DT7, have no wide-band data, and for them the waveforms are missing.
Summary¶
| Claim | Outcome |
|---|---|
| dataset-composition | partially reproduced |
| fig1c-hfo-each-animal | not reproduced |
| fig1c-event-count | reproduced |
| figs1-participation | not reproduced |
| fig1d-phase-locking | partially reproduced |
| fig1f-shorter-than-ripples | reproduced |
| fig1g-slower-than-ripples | not reproduced |
| fig1h-no-amplitude-frequency-correlation | not reproduced |
| fig2-confined-to-ls | not reproduced |
| fig2ef-intermediate-ls | blocked: metadata insufficient |
| fig3-travelling-waves | blocked: metadata insufficient |
| fig4b-pre-post-unchanged | reproduced |
| fig4b-hfo-above-behav | reproduced |
| fig4-sessions-few-change | blocked: other (7 recordings with 5 or more pairs) |
| fig4c-temporal-bias | blocked: other (4 recordings with 10 or more pairs) |
| fig4d-compression | not reproduced |
| figs2-post-rate-increase | partially reproduced |
| fig5b-nrem-and-rates | partially reproduced |
| fig5d-lag | not reproduced |
| fig5e-state-coupling | not reproduced |
| fig5f-ripple-triggered-firing | not reproduced |
| fig5f-no-phase-coherence | not reproduced |
| fig6-optogenetic-hfos | blocked: data missing |
| fig7e-ca1-spikes | partially reproduced |
| fig7e-state | partially reproduced |
| fig7e-duration | partially reproduced |
| fig7e-ca1-depth | blocked: data missing |
| fig7e-replay | blocked: methods underspecified |
| fig7e-seqnmf | blocked: methods underspecified |
| fig7e-content-independent | blocked: methods underspecified |
| fig8-model | blocked: other |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's:
- The files hold no detected events, and the paper does not say which channels it used. LS HFOs and CA1 ripples are detected here on one channel each, chosen by a fixed rule, with the Methods' parameters and the buzcode detector's open details. The paper's noise channel is replaced by the EMG estimated from the LFP. LS firing at the detected events exceeds firing at random times, or before and after the events, by more than 3 SEM in one animal of six, participation is a sixth to a third of the paper's, and a third of the events coupled to ripples fall on the same sample as the ripple. Neither the detection channel nor the re-referenced one is shown to capture the paper's HFOs; most of the outcomes on HFOs follow from that.
- A check added after the first verification: events detected on the LS channel re-referenced to another LS channel on its shank. It removes most LS events that fell on the same sample as a CA1 ripple (34.8% to 4.1%) but changes no outcome. The phase locking of fig1d was also measured with the phase from another shank than each unit's; the trough preference remains.
hpcunits and channels are taken as CA1. CA3 recordings are left out of ripple analyses, as their CA1 reference channel is not in the files.- Where the paper does not state a test or model (Figure 1H, Figure 7E, the pair inclusion of Figure 4), the choice was fixed in the claims file before the analysis.
What was found about the dandiset:
- No detected HFOs or ripples, no detection channels, and no unit waveforms, peak channels or
quality metrics are deposited (the units table holds
clu_id,group_idandlocation). Electrode coordinates are NaN in every file, and the probes are described only by their layout. - The LS and hippocampal channels share a 120-200 Hz signal: their phase difference is consistent even at random times, and many LS events fall on the same sample as CA1 ripples. Re-referencing within the LS shank reduces this. Some files list "reference sites a few millimeters dorsal" as a device; which reference each channel used is not stated.
- Nine files have units without a region label (268 units); the
hpclabel does not separate CA1. - 25 units in 10 files have spike times out of order.
- The CA3 recordings have no CA1 reference channel; the optogenetic rats are not in the dandiset.
- Wide-band data, from which unit waveforms could be re-extracted, are in 27 of 98 files (DT1, DT2 and two DT9 files).
- Spikes, LFP, sleep states and trials are consistent with one clock.