Reproducing Peyrache et al. (2015) from dandiset 000056¶
Paper. Peyrache A, Lacroix MM, Petersen PC, Buzsaki G (2015). Internally organized mechanisms of the head direction sense. Nature Neuroscience 18, 569-575. doi:10.1038/nn.3968. This attempt used the full text of the NIH author manuscript in PubMed Central (PMC4376557): Results, figure legends and Methods. The supplementary files (Supplementary Table 1 and Figures) could not be downloaded.
Dandiset. DANDI:000056, version
0.250624.0430. The paper (2015) cites no DANDI version. This is the most recent of the two
published versions; the earlier one, 0.210812.1518, holds the same 40 assets (same asset ids
and sizes, shown below) and differs only in metadata (contributors, affiliations and a
relatedResource entry for the paper), none of which the analysis uses.
What the dandiset holds. 40 sessions of 7 mice: spike times of sorted units, LED positions, sleep states and LFP. The units carry no electrode, shank or brain region, so the antero-dorsal thalamic nucleus (ADn) and the post-subiculum (PoS) cannot be separated. The tested claims pool the HD cells of each session; the claims that compare the two structures are blocked.
No analysis code is linked by the paper or the dandiset. No data from outside the dandiset are used. Random numbers are used only for the Poisson surrogates of fig2f-sws-tenfold-faster and fig4b-noise-correlations, each section with its own fixed seed.
| Claim | Tested |
|---|---|
| hd-cell-fraction: 353 of 1077 cells are HD cells | tested (pooled over regions) |
| fig2a-rates-preserved: rates correlated across states | tested |
| fig2b-pair-correlations-preserved: pair correlations similar across states | tested |
| fig3a-packet-60-degrees: correlation crosses zero at about 60 degrees | tested |
| fig2f-sws-tenfold-faster: internal HD signal about tenfold faster in SWS | tested |
| fig2f-rem-similar-to-wake: REM speed similar to waking | tested |
| fig4b-noise-correlations: excess at 0 degrees, dip at 60 degrees | tested |
| decoding-tracks-hd: decoded HD tracks the head in waking | tested |
| ADn against PoS: peak rate, HD information, ensemble sizes, correlation strength, coherence, noise correlations, OLE, decoded tuning, GLM prediction, synchrony, monosynaptic connections, shared preferred direction of connected pairs, LFP (13 claims) | blocked, no region labels |
In [1]:
import os
import re
os.environ.setdefault("OMP_NUM_THREADS", "1")
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from scipy.ndimage import gaussian_filter1d
import hd
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started_all = started = time.time()
sessions = hd.load_all() # LOAD
print(f"loaded {len(sessions)} files from dandiset {hd.DANDISET_ID} version {hd.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
loaded 40 files from dandiset 000056 version 0.250624.0430 in 237 s
Data check¶
Each check of the metadata against the values is stated under its output.
In [2]:
ids_new, ids_old = hd.asset_ids(hd.DANDISET_VERSION), hd.asset_ids(hd.EARLIER_VERSION)
print(f"Assets: {len(ids_new)} in {hd.DANDISET_VERSION}, {len(ids_old)} in {hd.EARLIER_VERSION}; "
f"same paths: {set(ids_new) == set(ids_old)}; same asset ids and sizes: {ids_new == ids_old}")
info = pd.DataFrame([{
"session": s["session"], "mouse": s["mouse"], "units": len(s["units"]),
"units_0_spikes": int((s["n_spikes"] == 0).sum()), "spikes_M": round(s["n_spikes"].sum() / 1e6, 2),
"hours": round(s["duration_s"] / 3600, 2), "states": "sws" in s["states"], "position": s["has_position"],
"raw_ephys": s["has_raw"], "groups": len(s["electrode_groups"]),
"sites": "+".join(str(v) for v in s["electrode_groups"].values()),
"annotations": ",".join(s["annotations"])} for s in sessions])
print("\nSessions per mouse:")
print(info.groupby("mouse").agg(sessions=("session", "size"), units=("units", "sum")).to_string())
print(f"\nMice {info.mouse.nunique()} (paper 7); sessions {len(info)} (paper 42); unit rows {info.units.sum()} "
f"(paper: 720 thalamic + 357 PoS = 1077 recorded cells); rows with no spikes {info.units_0_spikes.sum()}")
print("Species:", sorted({s["species"] for s in sessions}), "; electrode locations:",
sorted({l for s in sessions for l in s["locations"]}))
print("Sessions without a sleep-state table:", info[~info.states].session.tolist())
print("Sessions without processed position:", info[~info.position].session.tolist())
print("\n" + info.to_string(index=False))
Assets: 40 in 0.250624.0430, 40 in 0.210812.1518; same paths: True; same asset ids and sizes: True
Sessions per mouse:
sessions units
mouse
Mouse12 5 615
Mouse17 9 382
Mouse20 5 111
Mouse24 4 126
Mouse25 10 301
Mouse28 6 359
Mouse32 1 65
Mice 7 (paper 7); sessions 40 (paper 42); unit rows 1959 (paper: 720 thalamic + 357 PoS = 1077 recorded cells); rows with no spikes 103
Species: ['Mus musculus'] ; electrode locations: ['unknown']
Sessions without a sleep-state table: ['Mouse32-140820']
Sessions without processed position: ['Mouse17-130201', 'Mouse32-140820']
session mouse units units_0_spikes spikes_M hours states position raw_ephys groups sites annotations
Mouse12-120806 Mouse12 75 3 18.70 6.22 True True True 13 8+8+8+8+8+8+8+8+4+4+4+4+5 dow,hdD,rip
Mouse12-120807 Mouse12 77 2 23.84 6.17 True True False 13 8+8+8+8+8+8+8+8+4+4+4+4+8 dow,hdD,rip
Mouse12-120808 Mouse12 81 2 27.45 6.95 True True False 13 8+8+8+8+8+8+8+8+4+4+4+4+8 dow,hdD,rip
Mouse12-120809 Mouse12 292 1 24.79 6.77 True True True 14 8+8+8+8+8+8+8+8+4+4+4+4+4+5 dow,hdD,rip
Mouse12-120810 Mouse12 90 1 19.31 5.52 True True True 14 8+8+8+8+8+8+8+8+4+4+4+4+4+5 dow,hdD,rip
Mouse17-130125 Mouse17 16 4 9.07 5.99 True True False 9 8+8+8+8+8+8+8+8+6 dow,rip
Mouse17-130128 Mouse17 39 2 10.35 5.23 True True False 9 8+8+8+8+8+8+8+8+6 dow,rip
Mouse17-130129 Mouse17 52 3 17.71 6.45 True True False 9 8+8+8+8+8+8+8+8+6 dow,hdD,rip
Mouse17-130130 Mouse17 46 2 19.95 6.31 True True False 9 8+8+8+8+8+8+8+8+6 dow,hdD,rip
Mouse17-130131 Mouse17 42 1 14.71 6.19 True True False 9 8+8+8+8+8+8+8+8+6 dow,rip
Mouse17-130201 Mouse17 46 2 18.60 5.46 True False False 9 8+8+8+8+8+8+8+8+6 dow,rip
Mouse17-130202 Mouse17 48 2 22.21 6.25 True True False 9 8+8+8+8+8+8+8+8+6 dow,rip
Mouse17-130203 Mouse17 50 0 9.64 6.93 True True False 9 8+8+8+8+8+8+8+8+6 dow,rip
Mouse17-130204 Mouse17 43 0 29.12 7.71 True True False 9 8+8+8+8+8+8+8+8+6 dow,rip
Mouse20-130514 Mouse20 17 5 2.39 4.34 True True False 9 8+8+8+8+8+8+7+6+14 dow,rip
Mouse20-130515 Mouse20 18 5 2.56 6.12 True True False 9 8+8+8+8+8+8+7+6+14 dow,rip
Mouse20-130516 Mouse20 20 4 4.17 6.85 True True False 10 8+8+8+8+8+8+7+6+8+6 dow,rip
Mouse20-130517 Mouse20 28 3 7.58 7.30 True True False 10 8+8+8+8+8+8+7+6+8+6 dow,rip
Mouse20-130520 Mouse20 28 1 7.93 5.37 True True False 10 8+8+8+8+8+8+7+6+8+6 dow,rip
Mouse24-131213 Mouse24 22 1 4.35 2.99 True True False 8 8+8+8+8+8+8+8+8
Mouse24-131216 Mouse24 34 1 12.37 6.01 True True False 8 8+8+8+8+8+8+8+8 dow
Mouse24-131217 Mouse24 33 1 8.97 5.08 True True False 8 8+8+8+8+8+8+8+8 dow
Mouse24-131218 Mouse24 37 2 11.84 5.72 True True False 8 8+8+8+8+8+8+8+8 dow
Mouse25-140123 Mouse25 34 1 8.42 5.85 True True False 8 8+8+8+8+8+8+8+8
Mouse25-140124 Mouse25 30 2 4.32 5.60 True True False 8 8+8+8+8+8+8+8+8
Mouse25-140128 Mouse25 27 3 5.16 5.77 True True False 8 8+8+8+8+8+8+8+8
Mouse25-140129 Mouse25 28 3 9.13 5.95 True True False 8 8+8+8+8+8+8+8+8
Mouse25-140130 Mouse25 33 3 9.33 6.52 True True False 8 8+8+8+8+8+8+8+8 dow,hdD
Mouse25-140131 Mouse25 34 2 6.59 5.18 True True False 8 8+8+8+8+8+8+8+8 dow
Mouse25-140203 Mouse25 26 2 3.50 5.41 True True False 8 8+8+8+8+8+8+8+8 dow
Mouse25-140204 Mouse25 34 3 8.64 6.54 True True False 8 8+8+8+8+8+8+8+8 dow
Mouse25-140205 Mouse25 32 3 5.19 4.72 True True False 8 8+8+8+8+8+8+8+8 dow
Mouse25-140206 Mouse25 23 3 6.50 5.29 True True False 8 8+8+8+8+8+8+8+8 dow
Mouse28-140310 Mouse28 47 2 8.73 5.94 True True False 11 10+10+10+10+10+10+4+8+8+8+8 dow
Mouse28-140311 Mouse28 54 3 14.55 7.32 True True False 11 10+10+10+10+10+10+4+8+8+8+8 dow,hdD
Mouse28-140312 Mouse28 72 0 23.09 8.48 True True False 11 10+10+10+10+10+10+4+8+8+8+8 dow
Mouse28-140313 Mouse28 73 2 17.85 7.85 True True False 11 10+10+10+10+10+10+4+8+8+8+8 dow
Mouse28-140317 Mouse28 54 0 19.41 6.70 True True False 11 10+10+10+10+10+10+4+8+8+8+8 dow
Mouse28-140318 Mouse28 59 0 17.58 6.20 True True False 11 10+10+10+10+10+10+4+8+8+8+8 dow,hdD
Mouse32-140820 Mouse32 65 23 5.14 2.71 False False False 9 8+8+8+8+8+8+8+8+5 dow,hdD,rip
Both published versions hold the same 40 assets (same asset ids and sizes), so the choice of
version does not change the data. The files hold 7 mice, as the paper has, and 40 sessions
against the paper's 42. The units table of every file holds only spike times: no electrode,
shank, brain region, waveform or quality measure. Every electrode's location is "unknown".
So no unit can be assigned to the ADn or the PoS from the files. There are many more unit rows than
the 1077 cells of the paper (printed above), and some rows have no spikes at all. Mouse32 has no
sleep-state table, and Mouse32 and Mouse17-130201 have no processed position, so they are
left out of every analysis that needs those (both have only raw LED positions in
acquisition).
In [3]:
# What the unit rows are: their ids, their rates, and what the electrode groups say.
print("Unit ids are 0 to N-1 in table order in every file:",
all(np.array_equal(s["unit_ids"], np.arange(len(s["units"]))) for s in sessions))
print("Units table columns: spike times only. Electrode group attributes (digits replaced by N):",
sorted({(re.sub(r"\d+", "N", a.get("description", "")), a.get("location", "")) for s in sessions for a in s["group_attrs"].values()}))
rate_bins = [-1, 0, 0.05, 0.5, 5, 20, 1e9]
rate_labels = ["0 spikes", "<0.05 Hz", "0.05-0.5", "0.5-5", "5-20", ">20 Hz"]
ur = pd.DataFrame([{"session": s["session"], "mouse": s["mouse"], "rate": n / s["duration_s"]} for s in sessions for n in s["n_spikes"]])
ur["rate_bin"] = pd.cut(ur.rate, rate_bins, labels=rate_labels)
tab_units = ur.pivot_table(index="mouse", columns="rate_bin", values="rate", aggfunc="size", observed=False)
tab_units["all"] = tab_units.sum(axis=1)
print("\nUnit rows by mean rate over the recording, per mouse:")
print(tab_units.to_string())
print("\nMouse12-120809:", ur[ur.session == "Mouse12-120809"].rate_bin.value_counts().reindex(rate_labels).to_dict())
print(f"Unit rows with a mean rate of at least 0.05 Hz: {int((ur.rate >= 0.05).sum())}; at least 0.5 Hz: {int((ur.rate >= 0.5).sum())} (paper: 1077 cells)")
Unit ids are 0 to N-1 in table order in every file: True
Units table columns: spike times only. Electrode group attributes (digits replaced by N): [('shankN electrodes', 'unknown')]
Unit rows by mean rate over the recording, per mouse:
rate_bin 0 spikes <0.05 Hz 0.05-0.5 0.5-5 5-20 >20 Hz all
mouse
Mouse12 9 30 176 201 135 64 615
Mouse17 16 0 46 151 96 73 382
Mouse20 18 0 4 43 31 15 111
Mouse24 5 0 6 64 22 29 126
Mouse25 25 0 29 132 68 47 301
Mouse28 7 0 37 168 84 63 359
Mouse32 23 0 2 27 7 6 65
Mouse12-120809: {'0 spikes': 1, '<0.05 Hz': 30, '0.05-0.5': 136, '0.5-5': 80, '5-20': 32, '>20 Hz': 13}
Unit rows with a mean rate of at least 0.05 Hz: 1826; at least 0.5 Hz: 1526 (paper: 1077 cells)
Unit ids run from 0 to N-1 in table order in every file, so they do not encode a shank and cluster number as Klusters output would. The units table has no column besides spike times, and every electrode group is described only as "shankN electrodes" with location "unknown", so neither the ids, the table order nor the groups map units to probes or regions. The unit rows cannot be brought down to the paper's 1077 cells from the files: 1826 rows have a mean rate of at least 0.05 Hz and 1526 of at least 0.5 Hz. Mouse12-120809 holds 292 rows, 167 of them below 0.5 Hz, against 75 to 90 in the other Mouse12 sessions.
In [4]:
# Which mice were recorded in the PoS as well? The files do not say. The Methods say 4 mice had
# the thalamus and CA1 wires for sleep scoring, and 3 had the thalamus and the PoS (21 sessions).
groups = info.groupby("mouse").agg(sites=("sites", lambda x: sorted(set(x))), ripples=("annotations", lambda x: sum("rip" in a for a in x)),
sessions=("session", "size"))
print(groups.to_string())
THALAMUS_ONLY = ["Mouse12", "Mouse17", "Mouse20", "Mouse32"]
DUAL = ["Mouse24", "Mouse25", "Mouse28"]
print("\nSessions in mice without ripple annotations or small electrode groups:",
int(info[info.mouse.isin(DUAL)].shape[0]), "(paper: 21 sessions in the 3 dual-site mice)")
sites ripples sessions mouse Mouse12 [8+8+8+8+8+8+8+8+4+4+4+4+4+5, 8+8+8+8+8+8+8+8+4+4+4+4+5, 8+8+8+8+8+8+8+8+4+4+4+4+8] 5 5 Mouse17 [8+8+8+8+8+8+8+8+6] 9 9 Mouse20 [8+8+8+8+8+8+7+6+14, 8+8+8+8+8+8+7+6+8+6] 5 5 Mouse24 [8+8+8+8+8+8+8+8] 0 4 Mouse25 [8+8+8+8+8+8+8+8] 0 10 Mouse28 [10+10+10+10+10+10+4+8+8+8+8] 0 6 Mouse32 [8+8+8+8+8+8+8+8+5] 1 1 Sessions in mice without ripple annotations or small electrode groups: 20 (paper: 21 sessions in the 3 dual-site mice)
The electrode groups and annotation series fit the Methods if Mouse12, Mouse17, Mouse20 and
Mouse32 are the four mice with the thalamic probe and CA1 wires. Each has small extra
groups of 4 to 14 sites besides 8-site shanks, and ripple annotations (rip), which need a CA1 signal.
Mouse24, Mouse25 and Mouse28 would then be the three dual-site mice, with 20 sessions here
against 21 in the paper. Mouse28 has a 6-shank probe with 10 sites per shank and 5 more groups of
4 to 8 sites, as two probes would give. This assignment is consistent with the data but not
shown by them, and even in these mice the files do not say which units are on which probe.
It is used only for variants. The annotation series dow ("DOWN ON 5", "UP ON ALL1"),
hdD ("HD DOWN ON", "HD UP ON") and rip (ripple start, peak and stop) have no description.
In [5]:
# Spike times: each unit's times should be sorted.
rows = []
for s in sessions:
f_unsorted = not s["spikes_sorted"]
rows.append({"session": s["session"], "unsorted": f_unsorted,
"unsorted_units": len(s["unsorted_units"]),
"units_over_2^20_spikes": int((s["n_spikes"] > 2 ** 20).sum()),
"max_mean_rate_hz": round(s["n_spikes"].max() / s["duration_s"], 1)})
order = pd.DataFrame(rows)
print("Files with a unit whose spike times are not sorted:", int(order.unsorted.sum()), "of", len(order))
print(order.to_string(index=False))
# Each unsorted unit: its mean rate, where its times step back (index of the last spike before
# the step, and that index + 1 divided by 65536), and the span of each sorted run.
rows = []
for s in sessions:
for q in s["unsorted_units"]:
rows.append({"session": s["session"], "unit": q["unit"], "spikes": q["n"],
"mean_rate_hz": round(q["n"] / s["duration_s"], 1), "steps_back": len(q["breaks"]),
"break_index": q["breaks"][0], "(index+1)/65536": round((q["breaks"][0] + 1) / 65536, 4),
"runs_first_last_s": "; ".join(f"{a:.0f}-{b:.0f}" for _, a, b in q["runs"]),
"recording_s": round(s["duration_s"]), "times_repeated": q["repeated"]})
uns = pd.DataFrame(rows)
print("\nUnsorted units:")
print(uns.to_string(index=False))
print(f"\n{len(uns)} unsorted units in {uns.session.nunique()} files; with exactly two runs: {int((uns.steps_back == 1).sum())}; "
f"break at an index of the form k x 65536 - 1: {int((((uns.break_index + 1) % 65536) == 0).sum())}; "
f"mean rate {uns.mean_rate_hz.min()} to {uns.mean_rate_hz.max()} Hz (median {uns.mean_rate_hz.median()})")
Files with a unit whose spike times are not sorted: 36 of 40
session unsorted unsorted_units units_over_2^20_spikes max_mean_rate_hz
Mouse12-120806 True 1 6 81.6
Mouse12-120807 True 1 4 135.7
Mouse12-120808 True 1 7 133.4
Mouse12-120809 True 1 4 114.7
Mouse12-120810 True 1 5 93.6
Mouse17-130125 True 1 3 167.0
Mouse17-130128 True 1 5 123.9
Mouse17-130129 True 1 6 137.4
Mouse17-130130 True 1 6 148.2
Mouse17-130131 True 1 5 81.6
Mouse17-130201 True 1 6 160.6
Mouse17-130202 True 1 8 144.0
Mouse17-130203 True 1 1 78.4
Mouse17-130204 True 1 8 207.9
Mouse20-130514 True 1 0 36.8
Mouse20-130515 True 1 0 30.1
Mouse20-130516 False 0 1 66.7
Mouse20-130517 True 1 2 110.8
Mouse20-130520 False 0 1 86.1
Mouse24-131213 True 1 1 99.8
Mouse24-131216 True 1 4 117.3
Mouse24-131217 True 1 3 84.2
Mouse24-131218 True 1 5 99.4
Mouse25-140123 True 1 3 83.9
Mouse25-140124 True 1 0 29.4
Mouse25-140128 True 1 0 49.0
Mouse25-140129 True 1 3 82.2
Mouse25-140130 True 1 2 96.2
Mouse25-140131 False 0 0 55.5
Mouse25-140203 False 0 0 30.5
Mouse25-140204 True 1 3 114.7
Mouse25-140205 True 1 2 105.5
Mouse25-140206 True 1 3 90.2
Mouse28-140310 True 1 2 55.9
Mouse28-140311 True 1 6 81.5
Mouse28-140312 True 1 7 109.6
Mouse28-140313 True 1 5 54.0
Mouse28-140317 True 1 5 109.1
Mouse28-140318 True 1 3 98.1
Mouse32-140820 True 1 1 112.6
Unsorted units:
session unit spikes mean_rate_hz steps_back break_index (index+1)/65536 runs_first_last_s recording_s times_repeated
Mouse12-120806 1 1255273 56.1 1 1245183 19.0 64-22374; 1-22374 22374 0
Mouse12-120807 1 765700 34.5 1 655359 10.0 1421-22223; 0-22223 22223 0
Mouse12-120808 1 2208479 88.2 1 2162687 33.0 619-25036; 0-25036 25036 0
Mouse12-120809 1 1315460 53.9 1 1179647 18.0 556-24389; 1-24387 24389 0
Mouse12-120810 1 851079 42.8 1 786431 12.0 1919-19865; 0-19865 19865 0
Mouse17-130125 3 3603827 167.0 1 2359295 36.0 143-21580; 0-21580 21580 22281
Mouse17-130128 3 921983 49.0 1 917503 14.0 342-18815; 36-18810 18815 0
Mouse17-130129 3 1825528 78.6 1 1769471 27.0 240-23238; 0-23237 23238 0
Mouse17-130130 3 2325304 102.4 1 2293759 35.0 189-23802; 3-23802 22701 0
Mouse17-130131 2 216599 9.7 1 196607 3.0 6122-22284; 0-22277 22284 0
Mouse17-130201 3 2301372 117.1 1 2293759 35.0 756-19652; 0-19637 19652 0
Mouse17-130202 2 1738353 77.2 1 1376255 21.0 887-22516; 0-22516 22516 0
Mouse17-130203 1 280850 11.3 1 262143 4.0 360-24959; 0-24959 24960 0
Mouse17-130204 1 78901 2.8 1 65535 1.0 13704-27744; 1-27737 27744 0
Mouse20-130514 2 483681 31.0 1 458751 7.0 2730-15626; 0-15624 15626 0
Mouse20-130515 2 663625 30.1 1 524287 8.0 1505-22017; 0-22017 22017 0
Mouse20-130517 1 2909404 110.8 1 196607 3.0 2833-26269; 0-26269 26269 2508
Mouse24-131213 1 530417 49.3 1 262143 4.0 703-10754; 0-10754 10754 0
Mouse24-131216 1 1445283 66.7 1 1376255 21.0 2179-21653; 0-21653 21653 0
Mouse24-131217 1 1034464 56.6 1 983039 15.0 1388-18271; 1-18270 18271 0
Mouse24-131218 1 1252354 60.8 1 851967 13.0 242-20591; 0-20591 20591 0
Mouse25-140123 1 1413734 67.1 1 393215 6.0 2989-21062; 0-21062 21062 3562
Mouse25-140124 1 405648 20.1 1 131071 2.0 5627-20165; 0-20165 20165 821
Mouse25-140128 2 89170 4.3 1 65535 1.0 4311-20778; 1-20778 20778 0
Mouse25-140129 2 876857 40.9 1 851967 13.0 1497-21435; 2-21433 21435 0
Mouse25-140130 2 2258573 96.2 1 851967 13.0 1466-23478; 0-23478 23478 11600
Mouse25-140204 1 2701678 114.7 1 1572863 24.0 162-23554; 0-23554 23554 16342
Mouse25-140205 1 1791112 105.5 1 1114111 17.0 24-16976; 0-16976 16976 11950
Mouse25-140206 1 1717877 90.2 1 524287 8.0 1900-19037; 0-19037 19037 6613
Mouse28-140310 1 1194576 55.9 1 1179647 18.0 287-21371; 2-21369 21371 0
Mouse28-140311 2 2062469 78.2 1 2031615 31.0 67-26365; 3-26363 26366 0
Mouse28-140312 1 3343931 109.6 1 524287 8.0 2420-30519; 0-30520 30520 7500
Mouse28-140313 1 1224314 43.3 1 1048575 16.0 812-28276; 1-28276 28276 0
Mouse28-140317 1 1929693 80.0 1 1900543 29.0 305-24115; 1-24115 24115 0
Mouse28-140318 1 769457 34.5 1 524287 8.0 1156-22334; 0-22334 22334 0
Mouse32-140820 23 1099621 112.6 1 720895 11.0 104-9763; 0-9763 9763 1775
36 unsorted units in 36 files; with exactly two runs: 36; break at an index of the form k x 65536 - 1: 36; mean rate 2.8 to 167.0 Hz (median 58.7)
In 36 of 40 files one unit's spike times are not sorted. Each of these 36 units holds exactly two sorted runs, and in all 36 the first run ends at an index of the form k x 65536 - 1. The second run starts within 36 s of the start of the recording, the first between 24 and 13704 s, and the two runs of a unit end within 15 s of each other. In 26 of the 36 units no time of the second run appears in the first; in the other 10, 821 to 22281 times do. This breaks the NWB rule that spike times are sorted and looks like an error in the conversion. The helper sorts every unit's times when it loads them. These units have mean rates of 2.8 to 167.0 Hz (median 58.7). None of the 34 in sessions with an open-field epoch is an HD cell (shown with the HD cells below).
In [6]:
pos = pd.DataFrame([{
"session": s["session"], "rate_hz": s["pos_rate"], "start_s": s["pos_start"], "unit": s["pos_unit"],
"conversion": s["pos_conversion"], "raw_unit": s["raw_led"]["unit"], "raw_conversion": s["raw_led"]["conversion"],
"raw_max": s["raw_led"]["max_first_10000"], "led_missing": round(s["led_missing"], 3),
"pos_minus_lfp_s": round(s["n_pos"] / s["pos_rate"] - s["duration_s"], 1),
"wake_start_min": s["wake_epoch"][0] / 60 if s["wake_epoch"] else np.nan,
"wake_min": (s["wake_epoch"][1] - s["wake_epoch"][0]) / 60 if s["wake_epoch"] else np.nan,
"extent_x_cm": round(s["wake_extent_cm"][0], 1) if s["wake_epoch"] else np.nan,
"extent_y_cm": round(s["wake_extent_cm"][1], 1) if s["wake_epoch"] else np.nan,
"candidate_runs_min": [(a, b) for a, b in s["expl_runs_min"] if b - a >= 15]} for s in sessions if s["has_position"]])
print(pos.to_string(index=False))
print(f"\nPosition rate {pos.rate_hz.unique()} Hz (Methods: resampled at 39 Hz).")
print(f"Open-field epoch, minutes: median {pos.wake_min.median():.1f}, range {pos.wake_min.min():.0f} to {pos.wake_min.max():.0f} "
f"(Methods: 30 to 45 min of foraging); within 27 to 49.5 min: {int(pos.wake_min.between(27, 49.5).sum())} of {pos.wake_min.notna().sum()}")
print(f"1st-99th percentile extent of the head in it (cm): x {pos.extent_x_cm.min():.0f} to {pos.extent_x_cm.max():.0f}, "
f"y {pos.extent_y_cm.min():.0f} to {pos.extent_y_cm.max():.0f} (arena 53 x 46 cm)")
session rate_hz start_s unit conversion raw_unit raw_conversion raw_max led_missing pos_minus_lfp_s wake_start_min wake_min extent_x_cm extent_y_cm candidate_runs_min Mouse12-120806 39.0625 0.0 meters 0.01 meters NaN 321.00 0.081 0.1 114.0 36.0 58.8 51.3 [(112, 152), (262, 290), (325, 372)] Mouse12-120807 39.0625 0.0 meters 0.01 meters NaN 382.00 0.134 0.2 95.0 33.0 60.6 51.8 [(93, 130), (253, 285), (328, 370)] Mouse12-120808 39.0625 0.0 meters 0.01 meters NaN 316.00 0.116 0.1 111.0 31.0 57.7 47.6 [(109, 144), (292, 315), (370, 417)] Mouse12-120809 39.0625 0.0 meters 0.01 meters NaN -1.00 0.294 0.1 123.0 36.0 59.5 50.0 [(121, 161), (309, 338), (356, 406)] Mouse12-120810 39.0625 0.0 meters 0.01 meters NaN 370.76 0.232 0.1 98.0 31.0 58.7 51.3 [(96, 131), (248, 276)] Mouse17-130125 39.0625 0.0 meters 0.01 meters NaN 283.00 0.010 0.1 132.0 31.0 52.0 59.5 [(130, 165)] Mouse17-130128 39.0625 0.0 meters 0.01 meters NaN 329.00 0.051 0.1 116.0 32.0 51.5 59.6 [(114, 150)] Mouse17-130129 39.0625 0.0 meters 0.01 meters NaN 304.00 0.007 0.0 204.0 57.0 51.6 59.6 [(202, 263)] Mouse17-130130 39.0625 0.0 meters 0.01 meters NaN 293.00 0.021 0.1 140.0 35.0 51.0 61.9 [(138, 177), (341, 378)] Mouse17-130131 39.0625 0.0 meters 0.01 meters NaN 316.00 0.008 0.0 140.0 41.0 58.4 49.0 [(138, 183)] Mouse17-130202 39.0625 0.0 meters 0.01 meters NaN 312.00 0.015 0.0 128.0 30.0 60.5 51.1 [(126, 160), (351, 375)] Mouse17-130203 39.0625 0.0 meters 0.01 meters NaN 310.00 0.006 0.1 159.0 41.0 53.8 62.7 [(157, 202), (382, 411)] Mouse17-130204 39.0625 0.0 meters 0.01 meters NaN 290.70 0.012 0.1 166.0 31.0 53.8 61.8 [(164, 199), (383, 427)] Mouse20-130514 39.0625 0.0 meters 0.01 meters NaN 282.68 0.023 0.1 102.0 31.0 53.1 62.5 [(100, 135)] Mouse20-130515 39.0625 0.0 meters 0.01 meters NaN 283.00 0.188 -665.3 142.0 30.0 59.3 51.5 [(140, 174)] Mouse20-130516 39.0625 0.0 meters 0.01 meters NaN 314.00 0.063 -4905.2 120.0 32.0 61.1 52.4 [(118, 154), (311, 329)] Mouse20-130517 39.0625 0.0 meters 0.01 meters NaN 281.00 0.097 -6538.0 133.0 31.0 52.9 61.5 [(131, 166), (309, 328)] Mouse20-130520 39.0625 0.0 meters 0.01 meters NaN 277.00 0.311 0.0 120.0 31.0 58.0 50.0 [(118, 153)] Mouse24-131213 39.0625 0.0 meters 0.01 meters NaN 496.00 0.011 0.0 2.0 23.0 60.8 54.6 [(0, 27)] Mouse24-131216 39.0625 0.0 meters 0.01 meters NaN 377.00 0.229 0.0 108.0 27.0 61.8 53.2 [(106, 137), (311, 360)] Mouse24-131217 39.0625 0.0 meters 0.01 meters NaN 358.00 0.142 0.0 105.0 30.0 51.7 61.8 [(103, 137), (286, 304)] Mouse24-131218 39.0625 0.0 meters 0.01 meters NaN 371.90 0.059 0.0 115.0 28.0 60.3 51.8 [(113, 145), (302, 343)] Mouse25-140123 39.0625 0.0 meters 0.01 meters NaN 376.60 0.004 0.1 102.0 36.0 58.4 55.8 [(100, 140), (300, 350)] Mouse25-140124 39.0625 0.0 meters 0.01 meters NaN 355.66 0.001 0.0 112.0 29.0 60.9 57.5 [(110, 143)] Mouse25-140128 39.0625 0.0 meters 0.01 meters NaN 365.00 0.004 0.0 131.0 31.0 77.4 59.2 [(129, 164)] Mouse25-140129 39.0625 0.0 meters 0.01 meters NaN 350.94 0.002 0.0 88.0 23.0 79.4 69.3 [(86, 113)] Mouse25-140130 39.0625 0.0 meters 0.01 meters NaN 373.00 0.001 0.1 142.0 32.0 57.2 61.2 [(140, 176), (355, 391)] Mouse25-140131 39.0625 0.0 meters 0.01 meters NaN 387.80 0.000 0.0 168.0 32.0 62.0 57.7 [(166, 202)] Mouse25-140203 39.0625 0.0 meters 0.01 meters NaN 355.00 0.001 0.0 162.0 32.0 64.9 57.8 [(160, 196)] Mouse25-140204 39.0625 0.0 meters 0.01 meters NaN 379.00 0.002 0.0 132.0 31.0 66.9 62.0 [(130, 165)] Mouse25-140205 39.0625 0.0 meters 0.01 meters NaN 320.00 0.003 0.0 138.0 31.0 66.6 61.8 [(136, 171)] Mouse25-140206 39.0625 0.0 meters 0.01 meters NaN 356.73 0.005 0.0 148.0 31.0 63.2 60.9 [(146, 181)] Mouse28-140310 39.0625 0.0 meters 0.01 meters NaN 474.00 0.001 0.1 103.0 24.0 57.3 58.7 [(101, 129), (322, 356)] Mouse28-140311 39.0625 0.0 meters 0.01 meters NaN 392.80 0.007 0.1 104.0 29.0 54.3 59.9 [(102, 135), (399, 422)] Mouse28-140312 39.0625 0.0 meters 0.01 meters NaN 451.25 0.000 0.1 152.0 32.0 55.5 61.1 [(150, 186)] Mouse28-140313 39.0625 0.0 meters 0.01 meters NaN 480.00 0.012 0.1 92.0 39.0 53.8 58.8 [(90, 133)] Mouse28-140317 39.0625 0.0 meters 0.01 meters NaN 485.60 0.029 0.0 114.0 30.0 56.3 61.3 [(112, 146)] Mouse28-140318 39.0625 0.0 meters 0.01 meters NaN 463.40 0.134 0.0 124.0 32.0 56.0 59.0 [(122, 158)] Position rate [39.0625] Hz (Methods: resampled at 39 Hz). Open-field epoch, minutes: median 31.0, range 23 to 57 (Methods: 30 to 45 min of foraging); within 27 to 49.5 min: 34 of 38 1st-99th percentile extent of the head in it (cm): x 51 to 79, y 48 to 69 (arena 53 x 46 cm)
Processed positions are stored with conversion 0.01 to meters, so the values are in cm, and
the extent of the head during the detected open-field epoch (about 50 to 80 cm per axis between the 1st and 99th
percentiles) is consistent with a 53 x 46 cm arena, given that the LEDs sit above the head and
the camera views at an angle. The raw LED series in acquisition are stated in meters with a
conversion of NaN, but hold values in the hundreds, which are camera pixels. The sampling rate
(39.06 Hz) matches the Methods. The files do not mark the open-field epoch: the "Awake" state
also covers waking in the home cage. It is found here from position (see the helper). In
Mouse12, later runs with a much larger extent (over 100 cm) are another environment, and in
several mice there is a second run in the same arena late in the session; the Methods describe
a single exploration, so the first run is used. Detected epochs last 23 to 57 min, around the
Methods' 30 to 45 min (they include 2 min trimmed at each end of the detection window). In three
Mouse20 sessions the position series stops 11 to 109 min before the LFP, after the exploration.
In [7]:
st = []
for s in sessions:
if "sws" not in s["states"]:
continue
t0, t1, lab = s["state_table"]
wake = s.get("wake_epoch")
row = {"session": s["session"], "sws_h": hd.epochs_total(s["states"]["sws"]) / 3600,
"rem_h": hd.epochs_total(s["states"]["rem"]) / 3600, "awake_h": hd.epochs_total(s["states"]["awake_label"]) / 3600,
"table_start_s": float(t0[0]), "table_end_h": float(t1[-1]) / 3600, "recording_h": s["duration_s"] / 3600,
"gaps_s": float(np.sum(t0[1:] - t1[:-1]))}
if wake:
inside = [(max(a, wake[0]), min(b, wake[1])) for a, b in zip(t0, t1)]
frac = {L: sum(max(0, b - a) for (a, b), l in zip(inside, lab) if l == L) / (wake[1] - wake[0]) for L in ("Awake", "Non-REM", "REM")}
row.update({f"wake_epoch_{k}": round(v, 3) for k, v in frac.items()})
row["outside_wake_h"] = (s["duration_s"] - (wake[1] - wake[0])) / 3600
st.append(row)
st = pd.DataFrame(st)
print(st.round(2).to_string(index=False))
print(f"\nSleep (SWS + REM) per session: mean {(st.sws_h + st.rem_h).mean():.2f} h, s.d. {(st.sws_h + st.rem_h).std():.2f}, "
f"range {(st.sws_h + st.rem_h).min():.2f} to {(st.sws_h + st.rem_h).max():.2f}")
print(f"Recording outside the open-field epoch: mean {st.outside_wake_h.mean():.2f} h, s.d. {st.outside_wake_h.std():.2f} "
"(Results: sleep sessions of 5 h, s.d. 1 h)")
print(f"Fraction of the open-field epoch labelled Awake: median {st['wake_epoch_Awake'].median():.3f}, min {st['wake_epoch_Awake'].min():.3f}")
session sws_h rem_h awake_h table_start_s table_end_h recording_h gaps_s wake_epoch_Awake wake_epoch_Non-REM wake_epoch_REM outside_wake_h Mouse12-120806 2.29 0.36 3.55 1.0 6.21 6.22 65.0 0.98 0.00 0.02 5.62 Mouse12-120807 2.99 0.36 2.80 1.0 6.17 6.17 81.0 1.00 0.00 0.00 5.62 Mouse12-120808 3.11 0.33 3.49 1.0 6.95 6.95 59.0 1.00 0.00 0.00 6.44 Mouse12-120809 3.37 0.36 3.01 1.0 6.77 6.77 96.0 1.00 0.00 0.00 6.17 Mouse12-120810 3.18 0.27 2.05 1.0 5.52 5.52 68.0 1.00 0.00 0.00 5.00 Mouse17-130125 3.21 0.52 2.24 1.0 5.99 5.99 99.0 0.99 0.01 0.00 5.48 Mouse17-130128 2.62 0.32 2.26 1.0 5.23 5.23 65.0 1.00 0.00 0.00 4.69 Mouse17-130129 3.58 0.38 2.47 1.0 6.45 6.45 91.0 1.00 0.00 0.00 5.50 Mouse17-130130 3.45 0.38 2.45 1.0 6.31 6.31 88.0 0.98 0.01 0.00 5.72 Mouse17-130131 3.34 0.39 2.43 1.0 6.19 6.19 91.0 1.00 0.00 0.00 5.51 Mouse17-130201 2.94 0.31 2.19 1.0 5.46 5.46 86.0 NaN NaN NaN NaN Mouse17-130202 3.06 0.38 2.77 1.0 6.25 6.25 145.0 1.00 0.00 0.00 5.75 Mouse17-130203 3.60 0.25 3.06 1.0 6.93 6.93 106.0 1.00 0.00 0.00 6.25 Mouse17-130204 3.50 0.44 3.74 1.0 7.71 7.71 88.0 0.99 0.00 0.00 7.19 Mouse20-130514 2.20 0.27 1.84 1.0 4.34 4.34 77.0 1.00 0.00 0.00 3.82 Mouse20-130515 3.39 0.57 2.13 1.0 6.12 6.12 87.0 1.00 0.00 0.00 5.62 Mouse20-130516 3.54 0.54 2.74 1.0 6.85 6.85 69.0 0.99 0.01 0.00 6.31 Mouse20-130517 4.02 0.56 2.69 1.0 7.30 7.30 70.0 0.99 0.00 0.01 6.78 Mouse20-130520 2.74 0.34 2.28 1.0 5.37 5.37 61.0 1.00 0.00 0.00 4.86 Mouse24-131213 1.49 0.20 1.29 1.0 2.99 2.99 36.0 1.00 0.00 0.00 2.60 Mouse24-131216 2.87 0.31 2.81 1.0 6.01 6.01 93.0 0.99 0.01 0.00 5.56 Mouse24-131217 2.93 0.36 1.76 1.0 5.07 5.08 85.0 1.00 0.00 0.00 4.58 Mouse24-131218 2.70 0.33 2.66 61.0 5.72 5.72 61.0 1.00 0.00 0.00 5.25 Mouse25-140123 3.10 0.32 2.41 1.0 5.85 5.85 69.0 0.99 0.00 0.01 5.25 Mouse25-140124 3.45 0.30 1.82 1.0 5.60 5.60 87.0 1.00 0.00 0.00 5.12 Mouse25-140128 3.61 0.33 1.81 1.0 5.77 5.77 104.0 0.99 0.00 0.00 5.26 Mouse25-140129 3.64 0.38 1.90 1.0 5.95 5.95 110.0 1.00 0.00 0.00 5.57 Mouse25-140130 3.31 0.33 2.86 1.0 6.52 6.52 90.0 1.00 0.00 0.00 5.99 Mouse25-140131 3.20 0.32 1.64 1.0 5.18 5.18 90.0 0.98 0.02 0.00 4.65 Mouse25-140203 3.12 0.28 1.99 1.0 5.41 5.41 78.0 1.00 0.00 0.00 4.88 Mouse25-140204 4.30 0.46 1.75 1.0 6.54 6.54 117.0 1.00 0.00 0.00 6.03 Mouse25-140205 2.71 0.20 1.79 1.0 4.72 4.72 55.0 0.97 0.03 0.00 4.20 Mouse25-140206 2.79 0.27 2.21 1.0 5.29 5.29 64.0 1.00 0.00 0.00 4.77 Mouse28-140310 2.66 0.27 2.98 1.0 5.94 5.94 99.0 1.00 0.00 0.00 5.54 Mouse28-140311 2.54 0.16 4.58 1.0 7.32 7.32 164.0 1.00 0.00 0.00 6.84 Mouse28-140312 3.60 0.34 4.48 1.0 8.48 8.48 193.0 1.00 0.00 0.00 7.94 Mouse28-140313 2.86 0.22 4.73 1.0 7.85 7.85 142.0 1.00 0.00 0.00 7.20 Mouse28-140317 2.46 0.05 4.13 1.0 6.70 6.70 178.0 1.00 0.00 0.00 6.20 Mouse28-140318 2.52 0.14 3.49 1.0 6.20 6.20 190.0 1.00 0.00 0.00 5.67 Sleep (SWS + REM) per session: mean 3.41 h, s.d. 0.60, range 1.69 to 4.76 Recording outside the open-field epoch: mean 5.56 h, s.d. 0.98 (Results: sleep sessions of 5 h, s.d. 1 h) Fraction of the open-field epoch labelled Awake: median 1.000, min 0.973
The state tables run from about 1 s to the end of the recording, in whole seconds with 1 s between rows. SWS and REM total 1.7 to 4.8 hours per session (mean 3.4). The time outside the open-field epoch averages 5.6 hours (s.d. 1.0), close to the Results' 5 hours (s.d. 1) for the sleep sessions. Nearly all of the detected open-field epoch is labelled "Awake", which supports the detection.
In [8]:
# Clock checks. (1) Position against spikes: the head-direction tuning of tuned units, with head
# direction shifted against the spikes. (2) The state table against position: head speed in
# "Awake" and "Non-REM" seconds, with the labels shifted in time.
shifts = np.round(np.arange(-3, 3.01, 0.1), 2)
clock = []
for s in sessions:
if not s.get("wake_epoch"):
continue
tuned = sorted([u for u in s["units"] if "spikes" in u], key=lambda u: -u["peak"] * u["r"])[:6]
w0, w1 = s["wake_epoch"]
R = np.array([[hd.von_mises_fit(hd.tuning_curve(hd.spike_angles(u["spikes"][(u["spikes"] >= w0 - 5) & (u["spikes"] < w1 + 5)], s, sh),
s["occupancy"]))[2] for sh in shifts] for u in tuned]).mean(0)
clock.append({"session": s["session"], "units": len(tuned), "best_shift_s": shifts[np.argmax(R)],
"r_at_best": R.max(), "r_at_0": R[shifts == 0][0], "r_at_3s": R[-1]})
clock = pd.DataFrame(clock)
print("Mean resultant length of tuning (6 most tuned units) against the shift of head direction:")
print(clock.round(2).to_string(index=False))
Mean resultant length of tuning (6 most tuned units) against the shift of head direction:
session units best_shift_s r_at_best r_at_0 r_at_3s
Mouse12-120806 6 -0.1 0.78 0.77 0.27
Mouse12-120807 6 -0.1 0.81 0.81 0.25
Mouse12-120808 6 -2.6 0.71 0.35 0.14
Mouse12-120809 6 0.0 0.91 0.91 0.25
Mouse12-120810 6 -0.1 0.92 0.90 0.28
Mouse17-130125 3 0.0 0.61 0.61 0.15
Mouse17-130128 6 -0.1 0.73 0.73 0.20
Mouse17-130129 6 -0.2 0.75 0.73 0.22
Mouse17-130130 6 -0.2 0.75 0.72 0.24
Mouse17-130131 6 -0.2 0.85 0.83 0.26
Mouse17-130202 6 0.0 0.60 0.60 0.19
Mouse17-130203 6 -0.1 0.68 0.67 0.15
Mouse17-130204 6 0.0 0.76 0.76 0.17
Mouse20-130514 2 -0.1 0.45 0.44 0.20
Mouse20-130515 5 0.0 0.48 0.48 0.19
Mouse20-130516 3 0.0 0.82 0.82 0.38
Mouse20-130517 6 -0.1 0.74 0.74 0.28
Mouse20-130520 1 -0.1 0.89 0.88 0.26
Mouse24-131213 5 -0.1 0.56 0.56 0.20
Mouse24-131216 6 -0.2 0.72 0.69 0.18
Mouse24-131217 6 0.0 0.81 0.81 0.24
Mouse24-131218 6 -0.1 0.84 0.83 0.20
Mouse25-140123 6 0.0 0.74 0.74 0.36
Mouse25-140124 6 -0.1 0.66 0.66 0.36
Mouse25-140128 6 -0.1 0.67 0.66 0.22
Mouse25-140129 6 -0.1 0.67 0.67 0.26
Mouse25-140130 6 0.0 0.91 0.91 0.33
Mouse25-140131 6 0.0 0.88 0.88 0.23
Mouse25-140203 6 -0.1 0.86 0.85 0.27
Mouse25-140204 6 0.0 0.79 0.79 0.28
Mouse25-140205 6 -0.1 0.84 0.84 0.35
Mouse25-140206 6 -0.1 0.86 0.85 0.32
Mouse28-140310 6 0.0 0.82 0.82 0.51
Mouse28-140311 6 0.0 0.75 0.75 0.33
Mouse28-140312 6 0.0 0.76 0.76 0.35
Mouse28-140313 6 -0.1 0.86 0.86 0.41
Mouse28-140317 6 -0.1 0.83 0.83 0.44
Mouse28-140318 6 -0.1 0.83 0.83 0.53
In [9]:
s808 = next(s for s in sessions if s["session"] == "Mouse12-120808")
tuned = sorted([u for u in s808["units"] if "spikes" in u], key=lambda u: -u["peak"] * u["r"])[:6]
w0, w1 = s808["wake_epoch"]
shifts808 = np.round(np.arange(-5, 1.01, 0.1), 2)
q_times, q_shifts = [], []
for q in range(4):
a, b = w0 + q * (w1 - w0) / 4, w0 + (q + 1) * (w1 - w0) / 4
R = np.array([[hd.von_mises_fit(hd.tuning_curve(hd.spike_angles(u["spikes"][(u["spikes"] >= a) & (u["spikes"] < b)], s808, sh),
s808["occupancy"]))[2] for sh in shifts808] for u in tuned]).mean(0)
q_times.append((a + b) / 2); q_shifts.append(shifts808[np.argmax(R)])
print(f"Mouse12-120808, quarter {q + 1} of the open-field epoch ({a:.0f} to {b:.0f} s): best shift {shifts808[np.argmax(R)]:+.1f} s, "
f"r {R.max():.2f} (at 0 s {R[shifts808 == 0][0]:.2f})")
Mouse12-120808, quarter 1 of the open-field epoch (6660 to 7125 s): best shift -1.5 s, r 0.71 (at 0 s 0.49) Mouse12-120808, quarter 2 of the open-field epoch (7125 to 7590 s): best shift -2.5 s, r 0.84 (at 0 s 0.31) Mouse12-120808, quarter 3 of the open-field epoch (7590 to 8055 s): best shift -2.8 s, r 0.86 (at 0 s 0.35)
Mouse12-120808, quarter 4 of the open-field epoch (8055 to 8520 s): best shift -3.4 s, r 0.80 (at 0 s 0.25)
In [10]:
lag_s = [-600, -120, -30, -10, 0, 10, 30, 120, 600]
rows = []
for s in sessions:
if "sws" not in s["states"] or not s["has_position"]:
continue
t0, t1, lab = s["state_table"]
sp = s["speed_1s"]
for L in lag_s:
row = {"session": s["session"], "shift_s": L}
for name, label in (("awake", "Awake"), ("sws", "Non-REM")):
sec = np.concatenate([np.arange(int(a), int(b)) for a, b, l in zip(t0, t1, lab) if l == label]) + L
sec = sec[(sec >= 0) & (sec < len(sp))]
row[name] = np.nanmedian(sp[sec])
rows.append(row)
speed = pd.DataFrame(rows)
speed["contrast"] = speed.awake - speed.sws
table = speed.pivot(index="session", columns="shift_s", values="contrast")
print("Median head speed in Awake minus Non-REM seconds (cm/s), state labels shifted by the given seconds:")
print(table.round(2).to_string())
print("\nSessions where the contrast is largest at shift 0:", int((table.idxmax(axis=1) == 0).sum()), "of", len(table),
"; shift of the largest contrast in the others:", table.idxmax(axis=1)[table.idxmax(axis=1) != 0].to_dict())
print("Median over sessions:", table.median().round(2).to_dict())
Median head speed in Awake minus Non-REM seconds (cm/s), state labels shifted by the given seconds:
shift_s -600 -120 -30 -10 0 10 30 120 600
session
Mouse12-120806 3.20 4.24 4.33 4.36 4.37 4.35 4.30 4.06 2.70
Mouse12-120807 3.85 4.86 4.96 5.00 5.01 5.01 4.98 4.81 3.26
Mouse12-120808 3.77 5.26 5.46 5.51 5.51 5.51 5.50 5.35 3.53
Mouse12-120809 5.15 6.02 6.21 6.25 6.27 6.25 6.19 6.07 4.70
Mouse12-120810 4.35 4.94 5.09 5.15 5.17 5.19 5.18 5.14 4.57
Mouse17-130125 2.74 3.04 3.20 3.24 3.25 3.21 3.17 2.98 2.89
Mouse17-130128 2.24 2.74 2.86 2.88 2.88 2.84 2.77 2.70 2.09
Mouse17-130129 1.05 3.02 3.26 3.35 3.36 3.30 3.29 3.14 1.33
Mouse17-130130 2.86 3.60 3.86 4.00 3.99 3.93 3.86 3.64 2.81
Mouse17-130131 2.36 3.14 3.45 3.56 3.60 3.57 3.51 3.25 2.61
Mouse17-130202 1.76 2.80 3.04 3.11 3.20 3.18 3.10 2.97 2.03
Mouse17-130203 2.22 2.84 2.94 3.00 3.01 2.98 2.98 2.83 2.14
Mouse17-130204 2.52 3.41 3.55 3.60 3.63 3.60 3.57 3.45 2.48
Mouse20-130514 1.46 2.02 2.18 2.20 2.15 2.11 2.00 1.79 1.05
Mouse20-130515 1.37 2.65 2.80 2.90 2.93 2.86 2.81 2.67 1.14
Mouse20-130516 2.47 2.87 2.93 2.98 3.01 2.98 2.96 2.87 2.22
Mouse20-130517 1.96 2.83 3.07 3.17 3.22 3.20 3.15 3.02 2.12
Mouse20-130520 1.37 2.31 2.46 2.53 2.57 2.56 2.55 2.39 1.26
Mouse24-131213 3.01 3.45 3.54 3.58 3.59 3.56 3.55 3.42 2.07
Mouse24-131216 2.14 3.58 3.82 3.89 3.93 3.92 3.94 3.79 2.34
Mouse24-131217 1.60 3.42 3.74 3.82 3.83 3.82 3.81 3.63 1.78
Mouse24-131218 2.30 3.39 3.65 3.70 3.70 3.66 3.67 3.54 2.81
Mouse25-140123 1.59 3.54 3.71 3.83 3.84 3.84 3.76 3.45 2.09
Mouse25-140124 0.30 1.93 2.47 2.56 2.56 2.47 2.37 2.01 0.40
Mouse25-140128 0.50 2.59 3.47 3.71 3.71 3.55 3.26 2.09 0.60
Mouse25-140129 0.32 2.88 3.69 4.07 4.07 3.85 3.61 2.47 0.26
Mouse25-140130 2.51 5.10 5.54 5.62 5.62 5.57 5.50 5.10 3.07
Mouse25-140131 0.12 4.21 5.07 5.24 5.26 5.18 4.99 4.32 0.22
Mouse25-140203 0.33 2.93 3.68 3.76 3.74 3.62 3.29 2.27 0.18
Mouse25-140204 0.69 4.09 4.69 4.85 4.83 4.71 4.57 3.93 0.54
Mouse25-140205 1.19 4.64 4.97 5.03 5.02 4.98 4.85 4.47 0.92
Mouse25-140206 3.45 5.02 5.25 5.39 5.41 5.35 5.25 5.04 3.08
Mouse28-140310 1.55 2.05 2.23 2.28 2.28 2.25 2.22 2.06 1.44
Mouse28-140311 1.19 1.60 2.03 2.23 2.24 2.11 1.94 1.51 1.17
Mouse28-140312 0.86 1.46 1.58 1.75 1.76 1.69 1.55 1.39 1.01
Mouse28-140313 0.41 0.96 1.27 1.36 1.38 1.29 1.15 0.76 0.36
Mouse28-140317 0.67 0.82 0.98 1.12 1.15 1.05 0.95 0.71 0.60
Mouse28-140318 0.31 0.57 0.83 0.97 0.99 0.86 0.72 0.46 0.28
Sessions where the contrast is largest at shift 0: 27 of 38 ; shift of the largest contrast in the others: {'Mouse12-120810': 10, 'Mouse17-130128': -10, 'Mouse17-130130': -10, 'Mouse20-130514': -10, 'Mouse24-131216': 30, 'Mouse24-131218': -10, 'Mouse25-140130': -10, 'Mouse25-140203': -10, 'Mouse25-140204': -10, 'Mouse25-140205': -10, 'Mouse28-140310': -10}
Median over sessions: {-600: 1.68, -120: 3.03, -30: 3.5, -10: 3.59, 0: 3.61, 10: 3.57, 30: 3.4, 120: 3.08, 600: 2.05}
In every session but one, head-direction tuning is strongest with the head direction shifted by -0.2 to 0 s against the spikes. This fits the Methods' note that the first video frame is misaligned by 0 to 60 ms, and shows that spikes and position share a clock. Mouse12-120808 is the exception: its best shift is -2.6 s overall and drifts from about -1.5 to -3.4 s across the epoch, so in this session position and spikes do not share a clock. It is kept as deposited, and with this offset no unit in it passes the HD criteria; with position realigned by the best shift of each quarter (a variant under hd-cell-fraction), 13 do. The state labels share the clock of position: the contrast in head speed between Awake and Non-REM seconds is largest at zero shift in 27 of 38 sessions, at -10 s in 9 others, at +10 s in one and at +30 s in one, and falls as the labels are shifted further (table above). Whether the state table and the spikes share a clock is consistent with this, as spikes and position do, but is not tested on its own.
In [11]:
# Recordings deposited twice: compare unit spike counts, recording lengths and position lengths.
sig = pd.DataFrame([{"session": s["session"], "n_units": len(s["units"]), "duration_s": round(s["duration_s"], 1),
"n_pos": s.get("n_pos", -1), "counts": tuple(sorted(s["n_spikes"]))} for s in sessions])
dups = [(a.session, b.session) for i, a in sig.iterrows() for j, b in sig.iterrows()
if i < j and (a.counts == b.counts or a.duration_s == b.duration_s or (a.n_pos == b.n_pos and a.n_pos > 0))]
print("Pairs of files sharing spike counts, recording length or position length:", dups)
print("Units are not matched across sessions; every analysis is within a session.")
Pairs of files sharing spike counts, recording length or position length: [] Units are not matched across sessions; every analysis is within a session.
No two files share their unit spike counts, recording length or number of position frames, so no recording appears to be deposited twice. The analyses do not treat a unit as the same across sessions.
HD cells¶
HD cells are classified in each session with an open-field epoch, as in the helper (Methods' criteria: concentration > 1, peak > 1 Hz, Rayleigh p < 0.001).
In [12]:
cells = pd.DataFrame([{
"session": s["session"], "mouse": s["mouse"], "unit": u["unit"], "n_total": u["n"], "n_wake": u.get("n_wake", 0),
"kappa": u.get("kappa", np.nan), "peak": u.get("peak", np.nan), "rayleigh_p": u.get("rayleigh_p", np.nan),
"pref": u.get("pref", np.nan), "is_hd": u.get("is_hd", False),
"rate_wake": u.get("rate_wake", np.nan), "rate_rem": u.get("rate_rem", np.nan), "rate_sws": u.get("rate_sws", np.nan),
"mean_rate": u["n"] / s["duration_s"], "unsorted": u["unsorted"]}
for s in sessions if s.get("wake_epoch") for u in s["units"]])
hdc = cells[cells.is_hd]
print(f"Units in sessions with an open-field epoch: {len(cells)}; with spikes in it: {(cells.n_wake > 0).sum()}; HD cells: {len(hdc)}")
print(f"Unsorted units in these sessions: {int(cells.unsorted.sum())}; HD cells among them: {int(hdc.unsorted.sum())}")
print(f"HD cells: peak rate {hdc.peak.min():.1f} to {hdc.peak.max():.1f} Hz, concentration {hdc.kappa.min():.2f} to {hdc.kappa.max():.2f}")
Units in sessions with an open-field epoch: 1848; with spikes in it: 1731; HD cells: 378 Unsorted units in these sessions: 34; HD cells among them: 0 HD cells: peak rate 1.0 to 181.7 Hz, concentration 1.01 to 11.60
hd-cell-fraction¶
Of the recorded cells, 242 of 720 in the anterior thalamus and 111 of 357 in the PoS were classified as HD cells, 353 of 1077 in all.
In [13]:
readings = {"(a) units with spikes in the open-field epoch": cells.n_wake > 0, "(b) every unit row": cells.n_wake >= 0}
for name, m in readings.items():
print(f"{name}: HD {int(cells[m].is_hd.sum())} of {int(m.sum())} = {cells[m].is_hd.mean():.3f} "
f"(expected 353 +- 35.3 and 0.328 +- 0.033)")
by_mouse = cells[cells.n_wake > 0].groupby("mouse").agg(sessions=("session", "nunique"), units=("unit", "size"), hd=("is_hd", "sum"))
by_mouse["fraction"] = (by_mouse.hd / by_mouse.units).round(3)
print("\nBy mouse (reading a):")
print(by_mouse.to_string())
per_session = cells[cells.n_wake > 0].groupby("session").agg(units=("unit", "size"), hd=("is_hd", "sum"))
print("\nHD cells per session: median", per_session.hd.median(), ", range", per_session.hd.min(), "to", per_session.hd.max())
(a) units with spikes in the open-field epoch: HD 378 of 1731 = 0.218 (expected 353 +- 35.3 and 0.328 +- 0.033)
(b) every unit row: HD 378 of 1848 = 0.205 (expected 353 +- 35.3 and 0.328 +- 0.033)
By mouse (reading a):
sessions units hd fraction
mouse
Mouse12 5 595 64 0.108
Mouse17 8 315 55 0.175
Mouse20 5 90 14 0.156
Mouse24 4 118 32 0.271
Mouse25 10 267 103 0.386
Mouse28 6 346 110 0.318
HD cells per session: median 10.0 , range 0 to 26
In [14]:
print("Variants (reading a):")
m = cells.n_wake > 0
th = m & cells.mouse.isin(THALAMUS_ONLY)
print(f" mice that appear thalamus only: HD {int(cells[th].is_hd.sum())} of {int(th.sum())} = {cells[th].is_hd.mean():.3f} "
"(paper, thalamus: 242 of 720 = 0.336)")
du = m & cells.mouse.isin(DUAL)
print(f" mice that appear dual-site: HD {int(cells[du].is_hd.sum())} of {int(du.sum())} = {cells[du].is_hd.mean():.3f}")
no808 = m & (cells.session != "Mouse12-120808")
print(f" without Mouse12-120808 (clock offset): HD {int(cells[no808].is_hd.sum())} of {int(no808.sum())} = {cells[no808].is_hd.mean():.3f}")
for name, cond in {
"concentration > 0.8": (cells.kappa > 0.8) & (cells.peak > 1) & (cells.rayleigh_p < 0.001),
"concentration > 1.2": (cells.kappa > 1.2) & (cells.peak > 1) & (cells.rayleigh_p < 0.001),
"concentration > 1, peak > 5 Hz": (cells.kappa > 1) & (cells.peak > 5) & (cells.rayleigh_p < 0.001),
"unsorted units left out": cells.is_hd & ~cells.unsorted}.items():
print(f" {name}: HD {int(cond[m].sum())} of {int(m.sum())} = {cond[m].mean():.3f}")
fig, ax = plt.subplots(1, 2, figsize=(9, 3))
ax[0].hist(np.log10(cells[m].kappa.clip(1e-3)), 60, color="0.6")
ax[0].axvline(0, color="k", lw=1); ax[0].set_xlabel("log10 concentration"); ax[0].set_ylabel("units")
ax[1].hist(hdc.pref, np.arange(0, 361, 15), color="0.4"); ax[1].set_xlabel("preferred direction of HD cells (deg)")
plt.tight_layout(); plt.show()
print(f"Concentration plotted from {cells[m].kappa.min():.3g} to {cells[m].kappa.max():.3g} (log scale, values below 1e-3 at 1e-3)")
Variants (reading a): mice that appear thalamus only: HD 133 of 1000 = 0.133 (paper, thalamus: 242 of 720 = 0.336) mice that appear dual-site: HD 245 of 731 = 0.335 without Mouse12-120808 (clock offset): HD 378 of 1653 = 0.229 concentration > 0.8: HD 424 of 1731 = 0.245 concentration > 1.2: HD 357 of 1731 = 0.206 concentration > 1, peak > 5 Hz: HD 285 of 1731 = 0.165 unsorted units left out: HD 378 of 1731 = 0.218
Concentration plotted from 0.00477 to 91.8 (log scale, values below 1e-3 at 1e-3)
In [15]:
# Variant: Mouse12-120808 with position realigned to the spikes, by the best shift of each
# quarter of the open-field epoch (data check above), interpolated linearly between the
# quarters' centres; and with one constant shift (the median of the quarters). The shifts were
# chosen on the 6 most tuned units, so this favours those units.
url808 = dict(hd.asset_urls())[s808["path"]]
n_other = int(cells[m & (cells.session != "Mouse12-120808")].is_hd.sum())
for name, shift in (("per-quarter shift", (np.array(q_times), np.array(q_shifts))),
(f"constant shift {np.median(q_shifts):+.1f} s", (np.array([0.0]), np.array([np.median(q_shifts)])))):
s808r = hd.load_session(s808["path"], url808, hd_shift=shift)
hd808 = [u for u in s808r["units"] if u.get("is_hd")]
n_wake808 = sum(u.get("n_wake", 0) > 0 for u in s808r["units"])
print(f"Mouse12-120808, {name}: {len(hd808)} HD cells of {n_wake808} units with spikes in the epoch; with the other sessions "
f"{n_other + len(hd808)} HD cells (expected 353 +- 35.3, upper bound 388.3), fraction (a) {(n_other + len(hd808)) / int(m.sum()):.3f}")
print("HD cells in the other Mouse12 sessions:", per_session.loc[[x for x in per_session.index if x.startswith("Mouse12")], "hd"].to_dict())
# HD cells per session by mouse group, for the blocked claim ensemble-sizes (ADn 8.4 +- 5.1, PoS 5 +- 2.8).
ps = per_session.join(cells.groupby("session").mouse.first())
for name, mice in (("thalamus only", THALAMUS_ONLY), ("dual-site", DUAL)):
d = ps[ps.mouse.isin(mice)]
d1 = d.drop(index="Mouse12-120808", errors="ignore")
print(f"HD cells per session, mice that appear {name}: mean {d.hd.mean():.1f}, s.d. {d.hd.std():.1f} over {len(d)} sessions; "
f"without Mouse12-120808 {d1.hd.mean():.1f}, s.d. {d1.hd.std():.1f} over {len(d1)}")
Mouse12-120808, per-quarter shift: 13 HD cells of 78 units with spikes in the epoch; with the other sessions 391 HD cells (expected 353 +- 35.3, upper bound 388.3), fraction (a) 0.226
Mouse12-120808, constant shift -2.6 s: 12 HD cells of 78 units with spikes in the epoch; with the other sessions 390 HD cells (expected 353 +- 35.3, upper bound 388.3), fraction (a) 0.225
HD cells in the other Mouse12 sessions: {'Mouse12-120806': 21, 'Mouse12-120807': 13, 'Mouse12-120808': 0, 'Mouse12-120809': 20, 'Mouse12-120810': 10}
HD cells per session, mice that appear thalamus only: mean 7.4, s.d. 6.2 over 18 sessions; without Mouse12-120808 7.8, s.d. 6.1 over 17
HD cells per session, mice that appear dual-site: mean 12.2, s.d. 5.9 over 20 sessions; without Mouse12-120808 12.2, s.d. 5.9 over 20
Result: partially reproduced. By the Methods' criteria 378 units are HD cells, within the tolerance of the paper's 353. They are 0.218 of the 1731 units with spikes in the open-field epoch (0.205 of all 1848 rows of the 38 sessions with position), outside 0.328 plus or minus 0.033, because the files hold many more units than the paper's 1077 recorded cells; even the 1526 rows with a mean rate of 0.5 Hz or more exceed it. The fraction is 0.108 to 0.175 in the mice that appear to be thalamus only (Mouse12, Mouse17, Mouse20) and 0.271 to 0.386 in the others. The count passes only with Mouse12-120808 as deposited, where position and spikes do not share a clock and no unit passes the criteria. With position realigned to the spikes that session has 13 HD cells (12 with one constant shift of -2.6 s), and the count becomes 391 (390), above the upper bound of 388.3, so the match of the count is fragile. As the unit rows cannot be brought down to the paper's 1077 cells, the agreement of the count may also be coincidental. The regional split (242 of 720 thalamic, 111 of 357 PoS) cannot be tested.
fig2a-rates-preserved¶
The firing rates of HD neurons are strongly correlated between waking and SWS and between waking and REM sleep.
In [16]:
SLEEP = [s for s in sessions if s.get("wake_epoch") and "sws" in s["states"]]
sleep_names = {s["session"] for s in SLEEP}
r2a = hdc[hdc.session.isin(sleep_names)]
print(f"HD cells in sessions with sleep states: {len(r2a)} (expected 277 +- 27.7) in {r2a.session.nunique()} sessions")
print("Rates (Hz): " + ", ".join(f"{k} {r2a[f'rate_{k}'].min():.2f} to {r2a[f'rate_{k}'].max():.2f}" for k in hd.STATES))
for k in ("sws", "rem"):
r, p = stats.pearsonr(r2a.rate_wake, r2a[f"rate_{k}"])
pos_ = (r2a.rate_wake > 0) & (r2a[f"rate_{k}"] > 0)
rl, pl = stats.pearsonr(np.log10(r2a.rate_wake[pos_]), np.log10(r2a[f"rate_{k}"][pos_]))
rs, ps = stats.spearmanr(r2a.rate_wake, r2a[f"rate_{k}"])
print(f"wake vs {k.upper()}: Pearson r = {r:.3f} (p = {p:.2g}); log rates r = {rl:.3f} (n = {pos_.sum()}, p = {pl:.2g}); Spearman {rs:.3f}")
print("\nBy mouse (Pearson r of rates):")
print(r2a.groupby("mouse").apply(lambda d: pd.Series({"cells": len(d), "r_sws": stats.pearsonr(d.rate_wake, d.rate_sws)[0],
"r_rem": stats.pearsonr(d.rate_wake, d.rate_rem)[0]}), include_groups=False).round(3).to_string())
print("\nMedian rate ratio SWS/wake", round((r2a.rate_sws / r2a.rate_wake).median(), 3), ", REM/wake", round((r2a.rate_rem / r2a.rate_wake).median(), 3))
fig, ax = plt.subplots(1, 2, figsize=(8, 3.5))
for a, k in zip(ax, ("sws", "rem")):
a.loglog(r2a.rate_wake, r2a[f"rate_{k}"], ".", ms=4, color="0.3")
a.plot([0.05, 100], [0.05, 100], "k:", lw=1); a.set_xlabel("wake rate (Hz)"); a.set_ylabel(f"{k.upper()} rate (Hz)")
plt.tight_layout(); plt.show()
print("Points with a zero rate are not drawn on the log axes:", int(((r2a.rate_sws == 0) | (r2a.rate_rem == 0)).sum()))
HD cells in sessions with sleep states: 378 (expected 277 +- 27.7) in 37 sessions
Rates (Hz): wake 0.09 to 29.49, rem 0.12 to 29.40, sws 0.05 to 20.34
wake vs SWS: Pearson r = 0.838 (p = 3.7e-101); log rates r = 0.806 (n = 378, p = 9e-88); Spearman 0.840
wake vs REM: Pearson r = 0.891 (p = 4.5e-131); log rates r = 0.883 (n = 378, p = 2.5e-125); Spearman 0.901
By mouse (Pearson r of rates):
cells r_sws r_rem
mouse
Mouse12 64.0 0.873 0.915
Mouse17 55.0 0.855 0.906
Mouse20 14.0 0.889 0.942
Mouse24 32.0 0.718 0.928
Mouse25 103.0 0.832 0.866
Mouse28 110.0 0.836 0.907
Median rate ratio SWS/wake 0.756 , REM/wake 0.892
Points with a zero rate are not drawn on the log axes: 0
Result: reproduced. Across 378 HD cells, waking rates correlate with SWS rates (r = 0.838) and with REM rates (r = 0.891), both p < 1e-100, and r is 0.72 to 0.94 in every mouse. The paper has 277 cells in this comparison; there are 378 here, as all HD cells of the sessions with sleep states are used.
fig2b-pair-correlations-preserved¶
The pairwise correlations between HD neurons, within and across structures, are similar across brain states.
In [17]:
def hd_list(s):
return [u for u in s["units"] if u.get("is_hd")]
pair_rows, session_counts = [], {}
for s in SLEEP:
us = hd_list(s)
if len(us) < 2:
continue
corr = {}
for k in hd.STATES:
b = hd.PAIR_BIN_S[k]
mask = hd.bin_mask(s["states"][k], s["duration_s"], b)
counts = np.array([hd.binned(u["spikes"], mask, b) for u in us])
pairs = [(i, j) for i in range(len(us)) for j in range(i + 1, len(us))]
corr[k] = hd.pearson_pairs(counts, pairs)
for n, (i, j) in enumerate(pairs):
rates = [[us[i][f"rate_{k}"], us[j][f"rate_{k}"]] for k in hd.STATES]
pair_rows.append({"session": s["session"], "mouse": s["mouse"], "i": us[i]["unit"], "j": us[j]["unit"],
"offset": abs(hd.circ_diff(us[i]["pref"], us[j]["pref"])),
"all_states": bool(np.all(np.array(rates) >= hd.MIN_RATE_HZ)),
"any_state": bool(np.any(np.all(np.array(rates) >= hd.MIN_RATE_HZ, axis=1))),
**{f"c_{k}": corr[k][n] for k in hd.STATES}})
pairs_df = pd.DataFrame(pair_rows)
for name, col in (("(a) both cells >= 0.5 Hz in every state", "all_states"), ("(b) both cells >= 0.5 Hz in at least one state", "any_state")):
d = pairs_df[pairs_df[col]].dropna(subset=["c_wake", "c_sws", "c_rem"])
out = [f"{name}: {len(d)} pairs in {d.session.nunique()} sessions"]
for k in ("sws", "rem"):
r, p = stats.pearsonr(d.c_wake, d[f"c_{k}"])
out.append(f"wake vs {k.upper()} r = {r:.3f} (p = {p:.2g})")
print("; ".join(out))
pa = pairs_df[pairs_df.all_states].dropna(subset=["c_wake", "c_sws", "c_rem"])
print("Correlation ranges: " + ", ".join(f"{k} {pa[f'c_{k}'].min():.3f} to {pa[f'c_{k}'].max():.3f}" for k in hd.STATES))
print("Pairs dropped under (a) for low rates:", int((~pairs_df.all_states).sum()), "of", len(pairs_df))
print("\nBy mouse, reading (a):")
print(pa.groupby("mouse").apply(lambda d: pd.Series({"pairs": len(d), "r_sws": stats.pearsonr(d.c_wake, d.c_sws)[0] if len(d) > 2 else np.nan,
"r_rem": stats.pearsonr(d.c_wake, d.c_rem)[0] if len(d) > 2 else np.nan}), include_groups=False).round(3).to_string())
th_pairs = pa[pa.mouse.isin(THALAMUS_ONLY)]
print(f"\nVariant, mice that appear thalamus only: {len(th_pairs)} pairs; r SWS {stats.pearsonr(th_pairs.c_wake, th_pairs.c_sws)[0]:.3f}, "
f"r REM {stats.pearsonr(th_pairs.c_wake, th_pairs.c_rem)[0]:.3f}")
fig, ax = plt.subplots(1, 2, figsize=(8, 3.5))
for a, k in zip(ax, ("sws", "rem")):
a.plot(pa.c_wake, pa[f"c_{k}"], ".", ms=3, color="0.3"); a.axhline(0, color="0.8", lw=0.5); a.axvline(0, color="0.8", lw=0.5)
a.set_xlabel("wake correlation (300 ms)"); a.set_ylabel(f"{k.upper()} correlation ({int(hd.PAIR_BIN_S[k] * 1000)} ms)")
plt.tight_layout(); plt.show()
(a) both cells >= 0.5 Hz in every state: 1573 pairs in 33 sessions; wake vs SWS r = 0.938 (p = 0); wake vs REM r = 0.963 (p = 0)
(b) both cells >= 0.5 Hz in at least one state: 2246 pairs in 35 sessions; wake vs SWS r = 0.925 (p = 0); wake vs REM r = 0.952 (p = 0)
Correlation ranges: wake -0.377 to 0.977, rem -0.281 to 0.935, sws -0.188 to 0.840
Pairs dropped under (a) for low rates: 888 of 2461
By mouse, reading (a):
pairs r_sws r_rem
mouse
Mouse12 344.0 0.955 0.976
Mouse17 161.0 0.950 0.976
Mouse20 16.0 0.913 0.970
Mouse24 110.0 0.923 0.958
Mouse25 335.0 0.941 0.960
Mouse28 607.0 0.927 0.950
Variant, mice that appear thalamus only: 521 pairs; r SWS 0.952, r REM 0.976
Result: reproduced. Over 1573 pairs of HD cells (reading a), the waking pair correlations correlate with those in SWS (r = 0.938) and REM (r = 0.963); reading (b) gives 0.925 and 0.952 over 2246 pairs, and every mouse gives at least 0.91.
fig3a-packet-60-degrees¶
The zero-lag correlation between HD cell pairs falls with the offset of their preferred directions, positive below 60 degrees, negative above it and near zero at about 60 degrees, and this holds in waking, REM and SWS.
In [18]:
def crossing(d, k, width=10):
edges = np.arange(0, 181, width)
centers = (edges[:-1] + edges[1:]) / 2
means = np.array([d[(d.offset >= a) & (d.offset < b)][f"c_{k}"].mean() for a, b in zip(edges[:-1], edges[1:])])
for n in range(len(means) - 1):
if means[n] > 0 and means[n + 1] <= 0:
return centers[n] + (centers[n + 1] - centers[n]) * means[n] / (means[n] - means[n + 1]), means
return np.nan, means
for label, d in (("reading (a), 10 degree bins", pa), ("reading (b), 10 degree bins", pairs_df[pairs_df.any_state].dropna(subset=["c_wake", "c_sws", "c_rem"])),
("mice that appear thalamus only, (a)", th_pairs), ("mice that appear dual-site, (a)", pa[pa.mouse.isin(DUAL)])):
out = []
for k in hd.STATES:
x, _ = crossing(d, k)
r, p = stats.pearsonr(d.offset, d[f"c_{k}"])
out.append(f"{k} crossing {x:.1f} deg, r(offset) = {r:.3f}, p = {p:.2g}")
print(f"{label} ({len(d)} pairs): " + "; ".join(out))
print("\nVariant, 15 degree bins, reading (a): " + "; ".join(f"{k} {crossing(pa, k, 15)[0]:.1f} deg" for k in hd.STATES))
print("Mean correlation per 10 degree bin, reading (a):")
tab = pd.DataFrame({k: crossing(pa, k)[1] for k in hd.STATES}, index=[f"{a}-{a + 10}" for a in range(0, 180, 10)])
tab["pairs"] = [int(((pa.offset >= a) & (pa.offset < a + 10)).sum()) for a in range(0, 180, 10)]
print(tab.round(4).to_string())
fig, ax = plt.subplots(figsize=(5, 3.5))
for k, c in zip(hd.STATES, ("k", "tab:orange", "tab:green")):
ax.plot(np.arange(5, 180, 10), tab[k], "-o", ms=3, color=c, label=k)
ax.axhline(0, color="0.7", lw=0.5); ax.axvline(60, color="0.7", ls="--", lw=0.5)
ax.set_xlabel("offset of preferred directions (deg)"); ax.set_ylabel("mean correlation"); ax.legend(); plt.tight_layout(); plt.show()
reading (a), 10 degree bins (1573 pairs): wake crossing 54.2 deg, r(offset) = -0.751, p = 5.5e-286; rem crossing 56.9 deg, r(offset) = -0.741, p = 9.2e-274; sws crossing 55.9 deg, r(offset) = -0.671, p = 1.7e-206
reading (b), 10 degree bins (2246 pairs): wake crossing 54.1 deg, r(offset) = -0.728, p = 0; rem crossing 57.1 deg, r(offset) = -0.715, p = 0; sws crossing 56.5 deg, r(offset) = -0.632, p = 3.6e-251
mice that appear thalamus only, (a) (521 pairs): wake crossing 53.6 deg, r(offset) = -0.798, p = 5.5e-116; rem crossing 55.5 deg, r(offset) = -0.772, p = 2.3e-104; sws crossing 55.4 deg, r(offset) = -0.717, p = 1.5e-83
mice that appear dual-site, (a) (1052 pairs): wake crossing 54.6 deg, r(offset) = -0.727, p = 1e-173; rem crossing 57.4 deg, r(offset) = -0.725, p = 6.1e-172; sws crossing 56.1 deg, r(offset) = -0.649, p = 6.6e-127
Variant, 15 degree bins, reading (a): wake 55.4 deg; rem 57.9 deg; sws 57.7 deg
Mean correlation per 10 degree bin, reading (a):
wake rem sws pairs
0-10 0.5703 0.4960 0.3345 82
10-20 0.4831 0.4325 0.2870 80
20-30 0.3545 0.3348 0.2089 96
30-40 0.2302 0.2034 0.1263 89
40-50 0.0823 0.0808 0.0423 108
50-60 -0.0070 0.0129 0.0026 78
60-70 -0.0851 -0.0548 -0.0267 86
70-80 -0.1081 -0.0783 -0.0326 78
80-90 -0.1326 -0.0989 -0.0468 73
90-100 -0.1403 -0.1074 -0.0430 97
100-110 -0.1567 -0.1285 -0.0550 91
110-120 -0.1526 -0.1166 -0.0498 90
120-130 -0.1591 -0.1186 -0.0498 83
130-140 -0.1424 -0.1201 -0.0466 94
140-150 -0.1586 -0.1334 -0.0583 92
150-160 -0.1435 -0.1205 -0.0479 83
160-170 -0.1596 -0.1246 -0.0548 78
170-180 -0.1603 -0.1336 -0.0581 95
Result: reproduced. The mean zero-lag correlation crosses zero at 54.2 degrees in waking, 56.9 in REM and 55.9 in SWS, all within 54 to 66 degrees, and falls with the offset in each state (r = -0.75, -0.74 and -0.67, p < 1e-200). The waking crossing is at the edge of the tolerance; with 15 degree bins the crossings are 55.4 to 57.9 degrees, and under reading (b) the waking crossing is 54.1 degrees, just outside it.
fig2f-sws-tenfold-faster¶
The angular velocity of the internal HD signal, estimated from the temporal profile of pairwise cross-correlograms, is about tenfold faster during SWS than during waking.
In [19]:
started = time.time()
max_lag = int(round(hd.CCG_MAX_LAG_S / hd.CCG_BIN_S))
lags = np.arange(-max_lag, max_lag + 1) * hd.CCG_BIN_S
fits = []
for s in SLEEP:
us = hd_list(s)
keep = [(i, j) for i in range(len(us)) for j in range(i + 1, len(us))
if any(us[i][f"rate_{k}"] >= hd.MIN_RATE_HZ and us[j][f"rate_{k}"] >= hd.MIN_RATE_HZ for k in hd.STATES)]
if not keep:
continue
profiles = {}
for (i, j) in keep:
rho = hd.angular_correlation(us[i]["curve"], us[j]["curve"])
profiles[(i, j)] = hd.expected_profile(rho, hd.SPEED_GRID[-1] * hd.CCG_MAX_LAG_S, du=2.0)
for k in hd.STATES:
mask = hd.bin_mask(s["states"][k], s["duration_s"], hd.CCG_BIN_S)
counts = np.array([hd.binned(u["spikes"], mask, hd.CCG_BIN_S) for u in us])
ccgs = hd.ccg_all(counts, keep, max_lag)
for (i, j), ccg in zip(keep, ccgs):
v, edge, r2 = hd.fit_speed(ccg, lags, *profiles[(i, j)])
fits.append({"session": s["session"], "mouse": s["mouse"], "i": us[i]["unit"], "j": us[j]["unit"], "state": k,
"speed": v, "edge": edge, "r2": r2, "offset": abs(hd.circ_diff(us[i]["pref"], us[j]["pref"])),
"all_states": all(us[q][f"rate_{x}"] >= hd.MIN_RATE_HZ for q in (i, j) for x in hd.STATES)})
fits = pd.DataFrame(fits)
speeds_b = fits.pivot_table(index=["session", "mouse", "i", "j", "offset", "all_states"], columns="state", values="speed").reset_index()
r2 = fits.pivot_table(index=["session", "i", "j"], columns="state", values="r2").add_prefix("r2_").reset_index()
edges = fits.groupby(["session", "i", "j"]).edge.any().reset_index()
speeds_b = speeds_b.merge(edges, on=["session", "i", "j"]).merge(r2, on=["session", "i", "j"])
speeds = speeds_b[speeds_b.all_states]
print(f"Pairs fitted: {len(speeds)} in {speeds.session.nunique()} sessions ({time.time() - started:.0f} s); "
f"dropped with a fit at the edge of the speed range in any state: {int(speeds.edge.sum())}")
good = speeds[~speeds.edge].copy()
good["sws_wake"] = good.sws / good.wake
good["rem_wake"] = good.rem / good.wake
for k in hd.STATES:
print(f"{k}: median speed {good[k].median():.1f} deg/s (quartiles {good[k].quantile(0.25):.1f} to {good[k].quantile(0.75):.1f}; range {good[k].min():.1f} to {good[k].max():.1f})")
n_faster = int((good.sws_wake > 1).sum())
print(f"\nSWS/wake ratio: median {good.sws_wake.median():.2f} over {len(good)} pairs (expected 10 +- 1); quartiles "
f"{good.sws_wake.quantile(0.25):.2f} to {good.sws_wake.quantile(0.75):.2f}; SWS faster in {n_faster} of {len(good)} pairs "
f"(sign test p = {stats.binomtest(n_faster, len(good)).pvalue:.2g})")
print("\nBy mouse:")
print(good.groupby("mouse").agg(pairs=("sws_wake", "size"), sws_wake=("sws_wake", "median"), rem_wake=("rem_wake", "median"),
wake=("wake", "median"), sws=("sws", "median"), rem=("rem", "median")).round(2).to_string())
print("\nBy offset of preferred directions:")
good["offset_bin"] = pd.cut(good.offset, [0, 30, 60, 90, 180])
print(good.groupby("offset_bin", observed=True).agg(pairs=("sws_wake", "size"), sws_wake=("sws_wake", "median"), rem_wake=("rem_wake", "median")).round(2).to_string())
th = good[good.mouse.isin(THALAMUS_ONLY)]
print(f"\nVariant, mice that appear thalamus only: {len(th)} pairs, median SWS/wake {th.sws_wake.median():.2f}, REM/wake {th.rem_wake.median():.2f}")
good_b = speeds_b[~speeds_b.edge]
print(f"Variant, reading (b), both cells >= 0.5 Hz in at least one state: {len(speeds_b)} pairs, {int(speeds_b.edge.sum())} dropped at the edge; "
f"median SWS/wake {(good_b.sws / good_b.wake).median():.2f} and REM/wake {(good_b.rem / good_b.wake).median():.2f} over {len(good_b)} pairs")
print("Fit quality, R^2 of scale x expected + offset to the CCG, reading (a) pairs kept: " + ", ".join(
f"{k} median {good[f'r2_{k}'].median():.2f} (quartiles {good[f'r2_{k}'].quantile(0.25):.2f} to {good[f'r2_{k}'].quantile(0.75):.2f})" for k in hd.STATES))
for thr in (0.25, 0.5):
g = good[(good[[f"r2_{k}" for k in hd.STATES]] >= thr).all(axis=1)]
print(f"Variant, pairs with R^2 >= {thr} in every state: {len(g)} pairs, median SWS/wake {g.sws_wake.median():.2f}, REM/wake {g.rem_wake.median():.2f}")
print(f"Variant, all pairs including edge fits: median SWS/wake {(speeds.sws / speeds.wake).median():.2f}, REM/wake {(speeds.rem / speeds.wake).median():.2f}")
Pairs fitted: 1573 in 33 sessions (59 s); dropped with a fit at the edge of the speed range in any state: 64
wake: median speed 36.9 deg/s (quartiles 29.2 to 46.8; range 2.4 to 2139.2)
rem: median speed 33.5 deg/s (quartiles 20.7 to 49.0; range 2.1 to 1600.0)
sws: median speed 324.4 deg/s (quartiles 228.6 to 418.0; range 2.1 to 4025.7)
SWS/wake ratio: median 8.33 over 1509 pairs (expected 10 +- 1); quartiles 5.91 to 11.31; SWS faster in 1420 of 1509 pairs (sign test p = 4.1e-309)
By mouse:
pairs sws_wake rem_wake wake sws rem
mouse
Mouse12 338 8.88 0.84 41.48 404.35 35.65
Mouse17 154 8.51 0.85 34.61 298.09 30.10
Mouse20 16 13.52 1.07 34.80 466.18 31.38
Mouse24 104 8.41 0.69 41.11 356.72 31.09
Mouse25 330 8.01 0.79 38.02 333.51 30.93
Mouse28 567 7.89 0.99 34.36 291.33 35.28
By offset of preferred directions:
pairs sws_wake rem_wake
offset_bin
(0, 30] 256 7.69 0.88
(30, 60] 265 8.61 0.79
(60, 90] 230 4.92 0.68
(90, 180] 758 9.61 0.93
Variant, mice that appear thalamus only: 508 pairs, median SWS/wake 8.73, REM/wake 0.84
Variant, reading (b), both cells >= 0.5 Hz in at least one state: 2246 pairs, 117 dropped at the edge; median SWS/wake 8.02 and REM/wake 0.85 over 2129 pairs
Fit quality, R^2 of scale x expected + offset to the CCG, reading (a) pairs kept: wake median 0.84 (quartiles 0.62 to 0.94), rem median 0.55 (quartiles 0.28 to 0.78), sws median 0.78 (quartiles 0.51 to 0.90)
Variant, pairs with R^2 >= 0.25 in every state: 1063 pairs, median SWS/wake 8.46, REM/wake 0.87
Variant, pairs with R^2 >= 0.5 in every state: 748 pairs, median SWS/wake 8.56, REM/wake 0.88
Variant, all pairs including edge fits: median SWS/wake 8.33, REM/wake 0.86
In [20]:
# Reference for the waking estimate: the actual angular speed of the head in the open-field epoch.
actual = []
for s in SLEEP:
h = s["hd"].copy()
ok = ~np.isnan(h)
un = np.rad2deg(np.unwrap(np.deg2rad(np.interp(np.arange(len(h)), np.flatnonzero(ok), h[ok]))))
v = np.abs(np.diff(gaussian_filter1d(un, 2))) * s["pos_rate"]
actual.append({"session": s["session"], "actual_mean_deg_s": np.mean(v[ok[1:] & ok[:-1]])})
actual = pd.DataFrame(actual)
comp = good.groupby("session").wake.median().reset_index().merge(actual, on="session")
print("Actual mean angular head speed in waking (deg/s, head direction smoothed with a 2-frame s.d.): median over sessions "
f"{actual.actual_mean_deg_s.median():.1f} (range {actual.actual_mean_deg_s.min():.1f} to {actual.actual_mean_deg_s.max():.1f})")
print(f"Median CCG estimate of waking speed per session against the actual mean speed: Pearson r = "
f"{stats.pearsonr(comp.wake, comp.actual_mean_deg_s)[0]:.2f} over {len(comp)} sessions; median ratio estimate/actual {(comp.wake / comp.actual_mean_deg_s).median():.2f}")
fig, ax = plt.subplots(1, 2, figsize=(9, 3.2))
bins = np.geomspace(1, 5000, 50)
for k, c in zip(hd.STATES, ("k", "tab:orange", "tab:green")):
ax[0].hist(good[k], bins, histtype="step", color=c, label=k)
ax[0].set_xscale("log"); ax[0].set_xlabel("fitted mean angular speed (deg/s)"); ax[0].legend()
ax[1].hist(np.log10(good.sws_wake), 40, histtype="step", color="tab:green", label="SWS/wake")
ax[1].hist(np.log10(good.rem_wake), 40, histtype="step", color="tab:orange", label="REM/wake")
ax[1].axvline(1, color="0.6", ls="--"); ax[1].axvline(0, color="0.6", ls=":"); ax[1].set_xlabel("log10 ratio"); ax[1].legend()
plt.tight_layout(); plt.show()
print(f"Fitted speeds plotted from {good[list(hd.STATES)].min().min():.1f} to {good[list(hd.STATES)].max().max():.1f} deg/s, all within the axes")
Actual mean angular head speed in waking (deg/s, head direction smoothed with a 2-frame s.d.): median over sessions 62.3 (range 39.1 to 72.3) Median CCG estimate of waking speed per session against the actual mean speed: Pearson r = 0.20 over 33 sessions; median ratio estimate/actual 0.60
Fitted speeds plotted from 2.1 to 4025.7 deg/s, all within the axes
In [21]:
# Calibration on surrogates with a known speed. In each session, independent Poisson spike trains
# are drawn from each HD cell's waking tuning curve read at the head direction of the open-field
# epoch (10 ms bins, missing frames interpolated), and a second set from the same trajectory
# played ten times faster (wrapping around the epoch), with the same number of bins. Each
# reading (a) pair is fitted as above. The true mean speed is the mean absolute change of the
# interpolated head direction per bin.
rng_2f = np.random.default_rng(2)
cal = []
for s in SLEEP:
us = hd_list(s)
keep = [(i, j) for i in range(len(us)) for j in range(i + 1, len(us))
if all(us[i][f"rate_{k}"] >= hd.MIN_RATE_HZ and us[j][f"rate_{k}"] >= hd.MIN_RATE_HZ for k in hd.STATES)]
if not keep:
continue
w0, w1 = s["wake_epoch"]
h = s["hd"]
ok = ~np.isnan(h)
un = np.rad2deg(np.unwrap(np.deg2rad(np.interp(np.arange(len(h)), np.flatnonzero(ok), h[ok]))))
theta = np.interp(np.arange(w0, w1, hd.CCG_BIN_S) + hd.CCG_BIN_S / 2, s["hd_t"], un)
true_wake = np.mean(np.abs(np.diff(theta))) / hd.CCG_BIN_S
fast = theta[(np.arange(len(theta)) * 10) % len(theta)]
prof = {(i, j): hd.expected_profile(hd.angular_correlation(us[i]["curve"], us[j]["curve"]), hd.SPEED_GRID[-1] * hd.CCG_MAX_LAG_S, du=2.0)
for (i, j) in keep}
est = {}
for name, th_ in (("wake", theta), ("fast", fast)):
lam = np.array([u["curve"][np.floor(th_).astype(int) % 360] for u in us]) * hd.CCG_BIN_S
ccgs = hd.ccg_all(rng_2f.poisson(lam).astype(float), keep, max_lag)
est[name] = [hd.fit_speed(ccg, lags, *prof[p])[:2] for p, ccg in zip(keep, ccgs)]
for n, (i, j) in enumerate(keep):
cal.append({"session": s["session"], "offset": abs(hd.circ_diff(us[i]["pref"], us[j]["pref"])), "true_wake": true_wake,
"wake": est["wake"][n][0], "fast": est["fast"][n][0], "edge": est["wake"][n][1] or est["fast"][n][1]})
cal = pd.DataFrame(cal)
cg = cal[~cal.edge].copy()
cg["ratio"] = cg.fast / cg.wake
cg["est_over_true"] = cg.wake / cg.true_wake
print(f"Surrogate pairs: {len(cal)} in {cal.session.nunique()} sessions; dropped with a fit at the edge: {int(cal.edge.sum())}")
print(f"True mean speed of the surrogate trajectory (waking, per session): median {cal.groupby('session').true_wake.first().median():.1f} deg/s; "
f"surrogate waking estimate median {cg.wake.median():.1f} deg/s (data {good.wake.median():.1f}); estimate/true per pair: median {cg.est_over_true.median():.2f} (quartiles {cg.est_over_true.quantile(0.25):.2f} to {cg.est_over_true.quantile(0.75):.2f})")
print(f"Estimated ratio for a true ratio of 10: median {cg.ratio.median():.2f} over {len(cg)} pairs (quartiles {cg.ratio.quantile(0.25):.2f} to {cg.ratio.quantile(0.75):.2f})")
cg["offset_bin"] = pd.cut(cg.offset, [0, 30, 60, 90, 180])
print("By offset of preferred directions (surrogates):")
print(cg.groupby("offset_bin", observed=True).agg(pairs=("ratio", "size"), ratio=("ratio", "median"), est_over_true_wake=("est_over_true", "median")).round(2).to_string())
print(f"Data median SWS/wake ratio {good.sws_wake.median():.2f} against the surrogate median {cg.ratio.median():.2f} for a true 10: "
f"scaled by 10/{cg.ratio.median():.2f}, {good.sws_wake.median() * 10 / cg.ratio.median():.2f}")
Surrogate pairs: 1573 in 33 sessions; dropped with a fit at the edge: 11
True mean speed of the surrogate trajectory (waking, per session): median 99.9 deg/s; surrogate waking estimate median 37.0 deg/s (data 36.9); estimate/true per pair: median 0.37 (quartiles 0.31 to 0.44)
Estimated ratio for a true ratio of 10: median 4.99 over 1562 pairs (quartiles 2.76 to 7.02)
By offset of preferred directions (surrogates):
pairs ratio est_over_true_wake
offset_bin
(0, 30] 258 4.22 0.39
(30, 60] 274 5.75 0.37
(60, 90] 235 4.51 0.34
(90, 180] 795 5.27 0.36
Data median SWS/wake ratio 8.33 against the surrogate median 4.99 for a true 10: scaled by 10/4.99, 16.69
In [22]:
# Variant: angular speed of the Bayesian decoded HD (Methods: 10 ms bins, smoothed 300 ms s.d. in
# waking and REM and 50 ms in SWS), in sessions with at least 6 HD cells; mean absolute change
# of the decoded angle over 100 ms.
dec_rows = []
for s in SLEEP:
us = [u for u in hd_list(s)]
if len(us) < hd.MIN_HD_CELLS_DECODE:
continue
curves = np.array([u["curve"].reshape(-1, hd.DECODE_BIN_DEG).mean(1) for u in us])
row = {"session": s["session"], "mouse": s["mouse"], "hd_cells": len(us)}
for k in hd.STATES:
sd = 0.3 if k in ("wake", "rem") else 0.05
mask = hd.bin_mask(s["states"][k], s["duration_s"], hd.DECODE_BIN_S)
counts = np.array([gaussian_filter1d(hd.binned(u["spikes"], mask, hd.DECODE_BIN_S), sd / hd.DECODE_BIN_S) for u in us])
ang = hd.decode(counts, curves, hd.DECODE_BIN_S) * hd.DECODE_BIN_DEG
step = int(0.1 / hd.DECODE_BIN_S)
row[k] = np.mean(np.abs(hd.circ_diff(ang[step:], ang[:-step]))) / 0.1
dec_rows.append(row)
dec = pd.DataFrame(dec_rows)
dec["sws_wake"], dec["rem_wake"] = dec.sws / dec.wake, dec.rem / dec.wake
print(dec.round(2).to_string(index=False))
print(f"\nDecoder variant over {len(dec)} sessions: median SWS/wake {dec.sws_wake.median():.2f}, REM/wake {dec.rem_wake.median():.2f}")
session mouse hd_cells wake rem sws sws_wake rem_wake Mouse12-120806 Mouse12 21 37.20 37.26 499.31 13.42 1.00 Mouse12-120807 Mouse12 13 35.69 38.43 503.54 14.11 1.08 Mouse12-120809 Mouse12 20 40.42 30.67 405.09 10.02 0.76 Mouse12-120810 Mouse12 10 43.89 34.38 417.94 9.52 0.78 Mouse17-130128 Mouse17 6 61.36 77.18 515.18 8.40 1.26 Mouse17-130129 Mouse17 11 38.05 37.12 430.33 11.31 0.98 Mouse17-130130 Mouse17 12 36.19 28.79 420.86 11.63 0.80 Mouse17-130131 Mouse17 9 54.98 62.24 504.08 9.17 1.13 Mouse17-130203 Mouse17 6 58.42 61.05 446.93 7.65 1.04 Mouse20-130517 Mouse20 6 52.02 59.70 513.84 9.88 1.15 Mouse24-131216 Mouse24 9 54.04 56.48 475.40 8.80 1.05 Mouse24-131217 Mouse24 10 52.52 59.79 520.98 9.92 1.14 Mouse24-131218 Mouse24 11 52.02 43.37 497.62 9.57 0.83 Mouse25-140123 Mouse25 12 37.57 51.31 467.25 12.44 1.37 Mouse25-140124 Mouse25 7 40.99 36.76 361.77 8.83 0.90 Mouse25-140129 Mouse25 6 60.25 55.33 332.37 5.52 0.92 Mouse25-140130 Mouse25 13 34.58 43.19 391.98 11.33 1.25 Mouse25-140131 Mouse25 17 38.55 33.79 445.58 11.56 0.88 Mouse25-140203 Mouse25 10 35.58 33.71 450.07 12.65 0.95 Mouse25-140204 Mouse25 15 35.35 35.99 477.60 13.51 1.02 Mouse25-140205 Mouse25 11 38.76 41.54 463.68 11.96 1.07 Mouse25-140206 Mouse25 7 36.91 45.34 492.91 13.35 1.23 Mouse28-140310 Mouse28 15 27.35 38.98 433.62 15.86 1.43 Mouse28-140311 Mouse28 13 33.44 43.57 404.60 12.10 1.30 Mouse28-140312 Mouse28 26 29.93 36.60 356.07 11.90 1.22 Mouse28-140313 Mouse28 24 30.28 29.91 372.81 12.31 0.99 Mouse28-140317 Mouse28 16 29.90 31.45 344.03 11.51 1.05 Mouse28-140318 Mouse28 16 26.39 35.33 345.90 13.11 1.34 Decoder variant over 28 sessions: median SWS/wake 11.53, REM/wake 1.05
Result: partially reproduced. With the paper's cross-correlogram method, the internal HD signal is faster in SWS than in waking in 1420 of 1509 pairs, with a median ratio of 8.33, below the tolerance of 9 to 11 (quartiles 5.91 to 11.31; mice 7.89 to 13.52). Reading (b) gives 8.02 over 2129 pairs, and keeping only pairs whose fits reach R^2 of 0.25 or 0.5 in every state gives 8.46 and 8.56 (median R^2 0.84 in waking, 0.55 in REM, 0.78 in SWS). The estimator as implemented here does not recover known speeds. On Poisson surrogates driven by the measured head direction it returns 0.37 of the true mean speed in waking (median 37.0 deg/s, as on the data, against a true 99.9), and for the same trajectory played ten times faster it gives a median ratio of 4.99, not 10. Scaled by that compression (10/4.99), the data's 8.33 becomes 16.69. The data are therefore consistent with a speed-up below, at or above tenfold: the estimator shows that SWS is faster, but not by how much, and how far its bias in SWS matches that of the surrogates is not known. The decoder variant gives a median of 11.53 over 28 sessions.
fig2f-rem-similar-to-wake¶
The angular velocity of the internal HD signal is similar during REM sleep and waking.
In [23]:
w = stats.wilcoxon(np.log(good.rem_wake))
print(f"REM/wake ratio: median {good.rem_wake.median():.2f} over {len(good)} pairs (criterion 0.67 to 1.5); quartiles "
f"{good.rem_wake.quantile(0.25):.2f} to {good.rem_wake.quantile(0.75):.2f}; REM faster in {int((good.rem_wake > 1).sum())} pairs")
print(f"Wilcoxon signed-rank test of log ratio against 0: p = {w.pvalue:.2g}; median log10 ratio {np.log10(good.rem_wake).median():.3f}")
print(f"SWS/wake against REM/wake per pair: median of (SWS/wake)/(REM/wake) {(good.sws_wake / good.rem_wake).median():.2f}")
print(f"Decoder variant: median REM/wake {dec.rem_wake.median():.2f} over {len(dec)} sessions")
REM/wake ratio: median 0.87 over 1509 pairs (criterion 0.67 to 1.5); quartiles 0.60 to 1.26; REM faster in 581 pairs Wilcoxon signed-rank test of log ratio against 0: p = 6.3e-18; median log10 ratio -0.060 SWS/wake against REM/wake per pair: median of (SWS/wake)/(REM/wake) 8.90 Decoder variant: median REM/wake 1.05 over 28 sessions
Result: reproduced. The median REM to waking ratio is 0.87 over 1509 pairs (0.85 over 2129 pairs under reading b), inside 0.67 to 1.5; the decoder variant gives 1.05. By the cross-correlogram method the median REM speed is slightly below the waking one (Wilcoxon p = 6.3e-18, median log10 ratio -0.060), a test the paper does not make, but REM is faster in 581 of 1509 pairs, Mouse20 gives a ratio of 1.07 and the decoder variant 1.05.
fig4b-noise-correlations¶
During waking, the observed correlations of HD cell pairs exceed those expected from independent rate-coded neurons for pairs with overlapping fields (0 plus or minus 15 degrees), and fall below them for pairs at about 60 degrees.
In [24]:
def lagged(x, y, v, lag_bins, pearson=False):
"""Mean cross-product (or Pearson correlation) of x(t + lag) with y(t), lags -lag_bins..lag_bins,
over the times where both bins are valid (v, 1.0 or 0.0); x and y are zero in invalid bins."""
n, out = len(x), []
for L in range(-lag_bins, lag_bins + 1):
a, b, va, vb = (x[L:], y[:n - L], v[L:], v[:n - L]) if L >= 0 else (x[:n + L], y[-L:], v[:n + L], v[-L:])
c, sab = np.dot(va, vb), np.dot(a, b)
if not pearson:
out.append(sab / c)
continue
ma, mb = np.dot(a, vb) / c, np.dot(va, b) / c
out.append((sab / c - ma * mb) / np.sqrt((np.dot(a * a, vb) / c - ma ** 2) * (np.dot(va, b * b) / c - mb ** 2)))
return np.array(out)
rng_4b = np.random.default_rng(4) # Poisson surrogates of this section
noise_rows = []
lag_bins = int(round(0.05 / hd.CCG_BIN_S))
for s in SLEEP:
us = hd_list(s)
if len(us) < 2:
continue
w0, w1 = s["wake_epoch"]
edges = np.arange(w0, w1 + 1e-9, hd.CCG_BIN_S)
frame = np.round((edges[:-1] + hd.CCG_BIN_S / 2 - s["hd_t"][0]) * s["pos_rate"]).astype(int).clip(0, len(s["hd"]) - 1)
theta = s["hd"][frame]
valid = ~np.isnan(theta)
v = valid.astype(float)
spk = np.array([np.histogram(u["spikes"], edges)[0] * v for u in us])
inten = np.zeros_like(spk)
inten[:, valid] = np.array([u["curve"][np.floor(theta[valid]).astype(int) % 360] * hd.CCG_BIN_S for u in us])
# Surrogates: independent Poisson cells driven by the intensity alone. Their tuning curves are
# estimated again from the surrogate spikes (10 ms bins), and their own intensity is read from
# those, as for the data, so the bias includes that of estimating the curve from the same spikes.
sur = rng_4b.poisson(inten) * v
ang10 = np.floor(theta[valid]).astype(int) % 360
occ10 = np.bincount(ang10, minlength=360) * hd.CCG_BIN_S
inten_sur = np.zeros_like(spk)
inten_sur[:, valid] = np.array([hd.tuning_curve(np.repeat(theta[valid], sur[q, valid].astype(int)), occ10)[ang10] * hd.CCG_BIN_S
for q in range(len(us))])
for i in range(len(us)):
for j in range(i + 1, len(us)):
ok_k = [all(us[q][f"rate_{k}"] >= hd.MIN_RATE_HZ for q in (i, j)) for k in hd.STATES]
if not any(ok_k):
continue
obs, exp = lagged(spk[i], spk[j], v, lag_bins).mean(), lagged(inten[i], inten[j], v, lag_bins).mean()
noise_sur = lagged(sur[i], sur[j], v, lag_bins).mean() - lagged(inten_sur[i], inten_sur[j], v, lag_bins).mean()
chance = spk[i].sum() * spk[j].sum() / v.sum() ** 2
noise_rows.append({"session": s["session"], "mouse": s["mouse"], "offset": abs(hd.circ_diff(us[i]["pref"], us[j]["pref"])),
"all_states": all(ok_k), "noise_hz2": (obs - exp) / hd.CCG_BIN_S ** 2, "noise_rel": (obs - exp) / chance,
"obs_rel": obs / chance, "exp_rel": exp / chance, "sur_hz2": noise_sur / hd.CCG_BIN_S ** 2,
"noise_pearson": lagged(spk[i], spk[j], v, lag_bins, True).mean() - lagged(inten[i], inten[j], v, lag_bins, True).mean()})
noise_b = pd.DataFrame(noise_rows)
noise = noise_b[noise_b.all_states].copy()
print(f"Pairs (reading a): {len(noise)} in {noise.session.nunique()} sessions. Noise correlation = mean over lags -50..+50 ms of the spike-train "
f"cross-product minus the intensity cross-product, over bins where the head direction is known at both times, in Hz^2 "
f"(range {noise.noise_hz2.min():.2f} to {noise.noise_hz2.max():.2f}); also shown relative to the product of the mean counts (1 = chance)")
def groups_of(d):
return (("0 +- 15 deg", d.offset <= 15), ("60 +- 15 deg", d.offset.between(45, 75)))
for name, m in groups_of(noise):
d = noise[m]
wt = stats.wilcoxon(d.noise_hz2)
print(f"{name}: {len(d)} pairs, median noise {d.noise_hz2.median():.3f} Hz^2 (relative {d.noise_rel.median():+.3f}; observed "
f"{d.obs_rel.median():.3f}, expected {d.exp_rel.median():.3f}), {int((d.noise_hz2 > 0).sum())} positive; Wilcoxon p = {wt.pvalue:.2g}")
print("\nVariant, reading (b), both cells >= 0.5 Hz in at least one state:")
for name, m in groups_of(noise_b):
d = noise_b[m]
print(f" {name}: {len(d)} pairs, median {d.noise_hz2.median():.3f} Hz^2, {int((d.noise_hz2 > 0).sum())} positive, Wilcoxon p = {stats.wilcoxon(d.noise_hz2).pvalue:.2g}")
print("\nBias of the measure: the same measure on independent Poisson surrogates drawn from each cell's intensity, with tuning curves "
"estimated from the surrogate spikes (reading a):")
for name, m in groups_of(noise):
d = noise[m]
print(f" {name}: surrogate median {d.sur_hz2.median():+.3f} Hz^2 ({int((d.sur_hz2 > 0).sum())} of {len(d)} positive, Wilcoxon p = "
f"{stats.wilcoxon(d.sur_hz2).pvalue:.2g}); data minus surrogate median {(d.noise_hz2 - d.sur_hz2).median():+.3f} Hz^2, "
f"{int((d.noise_hz2 > d.sur_hz2).sum())} of {len(d)} above, Wilcoxon p = {stats.wilcoxon(d.noise_hz2 - d.sur_hz2).pvalue:.2g}")
grp = np.where(noise.offset <= 15, "0", np.where(noise.offset.between(45, 75), "60", "other"))
print("Per session, median surrogate (bias) and data noise at 0 and 60 deg (Hz^2):")
print(noise.assign(group=grp)[grp != "other"].pivot_table(index="session", columns="group", values=["sur_hz2", "noise_hz2"], aggfunc="median").round(2).to_string())
print("\nVariant, the first criterion (both cross-correlations as Pearson correlations, not comparable in scale):")
for name, m in groups_of(noise):
d = noise[m]
print(f" {name}: median {d.noise_pearson.median():+.4f}, Wilcoxon p = {stats.wilcoxon(d.noise_pearson).pvalue:.2g}")
noise["offset_bin"] = pd.cut(noise.offset, np.arange(0, 181, 15), include_lowest=True)
print("\nBy offset:")
print(noise.groupby("offset_bin", observed=True).agg(pairs=("noise_hz2", "size"), median_hz2=("noise_hz2", "median"), median_sur_hz2=("sur_hz2", "median"),
median_rel=("noise_rel", "median"), positive=("noise_hz2", lambda x: (x > 0).mean())).round(3).to_string())
print("\nBy mouse (median noise, Hz^2):")
print(noise.assign(group=grp).pivot_table(index="mouse", columns="group", values="noise_hz2", aggfunc=["median", "size"]).round(3).to_string())
th_n = noise[noise.mouse.isin(THALAMUS_ONLY)]
print("\nVariant, mice that appear thalamus only: " + "; ".join(
f"{name} {len(th_n[m])} pairs, median {th_n[m].noise_hz2.median():.3f} Hz^2" for name, m in groups_of(th_n)))
fig, ax = plt.subplots(figsize=(5, 3.5))
ax.plot(noise.offset, noise.noise_rel, ".", ms=3, color="0.4"); ax.axhline(0, color="k", lw=0.5)
ax.set_ylim(-0.5, 1.0); ax.set_xlabel("offset of preferred directions (deg)"); ax.set_ylabel("noise correlation (relative to chance)")
plt.tight_layout(); plt.show()
print("Points outside the y axis:", int(((noise.noise_rel < -0.5) | (noise.noise_rel > 1.0)).sum()), "of", len(noise))
Pairs (reading a): 1573 in 33 sessions. Noise correlation = mean over lags -50..+50 ms of the spike-train cross-product minus the intensity cross-product, over bins where the head direction is known at both times, in Hz^2 (range -229.10 to 503.06); also shown relative to the product of the mean counts (1 = chance)
0 +- 15 deg: 118 pairs, median noise 7.838 Hz^2 (relative +0.813; observed 3.782, expected 2.951), 110 positive; Wilcoxon p = 5.3e-19
60 +- 15 deg: 256 pairs, median noise -4.079 Hz^2 (relative -0.282; observed 0.657, expected 0.996), 35 positive; Wilcoxon p = 4.2e-32
Variant, reading (b), both cells >= 0.5 Hz in at least one state:
0 +- 15 deg: 171 pairs, median 4.222 Hz^2, 152 positive, Wilcoxon p = 2e-24
60 +- 15 deg: 378 pairs, median -1.626 Hz^2, 62 positive, Wilcoxon p = 6.8e-42
Bias of the measure: the same measure on independent Poisson surrogates drawn from each cell's intensity, with tuning curves estimated from the surrogate spikes (reading a):
0 +- 15 deg: surrogate median +0.644 Hz^2 (100 of 118 positive, Wilcoxon p = 4.7e-17); data minus surrogate median +6.967 Hz^2, 110 of 118 above, Wilcoxon p = 5e-19
60 +- 15 deg: surrogate median -0.365 Hz^2 (47 of 256 positive, Wilcoxon p = 8.1e-30); data minus surrogate median -3.580 Hz^2, 33 of 256 above, Wilcoxon p = 1.6e-31
Per session, median surrogate (bias) and data noise at 0 and 60 deg (Hz^2):
noise_hz2 sur_hz2
group 0 60 0 60
session
Mouse12-120806 14.28 -2.79 0.42 -0.26
Mouse12-120807 17.57 -4.67 1.13 -0.29
Mouse12-120809 9.08 -14.44 0.99 -1.97
Mouse12-120810 9.14 -16.78 0.98 -1.09
Mouse17-130128 6.38 -4.43 0.76 -0.27
Mouse17-130129 34.99 -8.57 2.59 -0.87
Mouse17-130130 10.06 -15.24 0.90 -1.19
Mouse17-130131 27.09 -9.62 2.25 -0.93
Mouse17-130202 NaN -0.86 NaN 0.28
Mouse17-130203 NaN 14.26 NaN 1.42
Mouse17-130204 7.22 -6.22 2.46 -1.78
Mouse20-130517 NaN -23.54 NaN -2.69
Mouse24-131216 1.53 4.38 -0.11 0.06
Mouse24-131217 1.66 -0.22 0.21 0.05
Mouse24-131218 3.59 -0.90 0.42 -0.25
Mouse25-140123 32.01 -0.28 6.12 -0.10
Mouse25-140129 NaN -0.60 NaN -0.11
Mouse25-140130 2.65 -6.26 0.56 -1.62
Mouse25-140131 2.64 -15.52 0.38 -4.06
Mouse25-140203 64.63 -9.99 5.76 -1.07
Mouse25-140204 6.22 -3.44 0.18 -0.42
Mouse25-140205 4.32 -6.18 0.76 -0.76
Mouse25-140206 22.02 -19.20 2.32 -3.01
Mouse28-140310 19.17 -0.82 0.74 -0.20
Mouse28-140311 8.71 -3.29 0.57 -0.30
Mouse28-140312 50.79 -3.24 2.94 -0.22
Mouse28-140313 4.19 -3.59 0.36 -0.47
Mouse28-140317 9.17 -0.83 0.95 -0.14
Mouse28-140318 11.80 -3.84 0.63 -0.06
Variant, the first criterion (both cross-correlations as Pearson correlations, not comparable in scale):
0 +- 15 deg: median -0.7808, Wilcoxon p = 4.2e-21
60 +- 15 deg: median -0.0117, Wilcoxon p = 0.024
By offset:
pairs median_hz2 median_sur_hz2 median_rel positive
offset_bin
(-0.001, 15.0] 118 7.838 0.644 0.813 0.932
(15.0, 30.0] 140 4.768 0.516 0.355 0.829
(30.0, 45.0] 145 -1.222 -0.082 -0.107 0.331
(45.0, 60.0] 130 -4.153 -0.299 -0.291 0.154
(60.0, 75.0] 126 -3.986 -0.476 -0.259 0.119
(75.0, 90.0] 111 -2.240 -0.340 -0.156 0.108
(90.0, 105.0] 139 -0.800 -0.171 -0.079 0.173
(105.0, 120.0] 139 -0.446 -0.070 -0.035 0.180
(120.0, 135.0] 126 -0.280 -0.066 -0.027 0.222
(135.0, 150.0] 143 -0.200 -0.033 -0.017 0.322
(150.0, 165.0] 118 -0.112 -0.010 -0.007 0.390
(165.0, 180.0] 138 -0.056 -0.051 -0.005 0.406
By mouse (median noise, Hz^2):
median size
group 0 60 other 0 60 other
mouse
Mouse12 13.170 -5.529 -0.662 32.0 51.0 261.0
Mouse17 14.898 -6.385 -0.575 12.0 29.0 120.0
Mouse20 NaN -23.536 -0.396 NaN 1.0 15.0
Mouse24 1.658 -0.691 -0.101 9.0 13.0 88.0
Mouse25 5.363 -4.937 -0.181 22.0 66.0 247.0
Mouse28 7.909 -3.258 -0.215 43.0 96.0 468.0
Variant, mice that appear thalamus only: 0 +- 15 deg 44 pairs, median 13.302 Hz^2; 60 +- 15 deg 81 pairs, median -6.385 Hz^2
Points outside the y axis: 115 of 1573
Result: reproduced. With the Methods' formula, over bins where the head direction is known at both times, the noise correlation is positive for pairs at 0 plus or minus 15 degrees (110 of 118 pairs, median 7.84 Hz^2, p = 5.3e-19) and negative for pairs at 60 plus or minus 15 degrees (35 of 256 positive, median -4.08 Hz^2, p = 4.2e-32). Reading (b) gives the same signs (medians 4.22 and -1.63 Hz^2, p = 2e-24 and 6.8e-42). On independent Poisson surrogates whose tuning curves are estimated from their own spikes, the measure is biased in the same direction (median +0.64 Hz^2 at 0 degrees and -0.37 Hz^2 at 60), about a tenth of the observed values; the data exceed their surrogates in 110 of 118 pairs at 0 degrees (p = 5e-19) and fall below them in 223 of 256 pairs at 60 degrees (p = 1.6e-31). The criterion was revised after a first run to follow the Methods' formula (see the claims file); under the first version, with Pearson scaling, the 0 degree pairs come out negative.
decoding-tracks-hd¶
The HD signal decoded with a Bayesian decoder from the HD cells faithfully tracks the head direction of the waking mouse.
In [25]:
def decode_session(s, select_on_train):
"""Errors (deg) of two-fold cross-validated decoding in the open-field epoch, and the cells used
per fold, or None when a fold has fewer than MIN_HD_CELLS_DECODE cells. With select_on_train,
HD cells are chosen by the Methods' criteria on the training half alone, from the tuned units
whose spike times the helper keeps (concentration above 0.5 and peak above 1 Hz over the epoch)."""
w0, w1 = s["wake_epoch"]
half = (w0 + w1) / 2
pool = [u for u in s["units"] if "spikes" in u] if select_on_train else hd_list(s)
errs, used = [], []
for train, test in (((w0, half), (half, w1)), ((half, w1), (w0, half))):
tr = (s["hd_t"] >= train[0]) & (s["hd_t"] < train[1]) & ~np.isnan(s["hd"])
occ = np.bincount(np.floor(s["hd"][tr]).astype(int) % 360, minlength=360) / s["pos_rate"]
curves, cells_ = [], []
for u in pool:
ang = hd.spike_angles(u["spikes"][(u["spikes"] >= train[0]) & (u["spikes"] < train[1])], s)
curve = hd.tuning_curve(ang, occ)
if select_on_train:
if curve.sum() <= 0:
continue
_, kappa, r = hd.von_mises_fit(curve)
if not (kappa > hd.KAPPA_MIN and curve.max() > hd.PEAK_MIN_HZ and hd.rayleigh_p(r, len(ang)) < hd.RAYLEIGH_P):
continue
curves.append(curve.reshape(-1, hd.DECODE_BIN_DEG).mean(1))
cells_.append(u)
if len(cells_) < hd.MIN_HD_CELLS_DECODE:
return None
edges = np.arange(test[0], test[1] + 1e-9, hd.DECODE_BIN_S)
counts = np.array([gaussian_filter1d(np.histogram(u["spikes"], edges)[0].astype(float), hd.DECODE_SMOOTH_S / hd.DECODE_BIN_S) for u in cells_])
decoded = (hd.decode(counts, np.array(curves), hd.DECODE_BIN_S) + 0.5) * hd.DECODE_BIN_DEG
frame = np.round((edges[:-1] + hd.DECODE_BIN_S / 2 - s["hd_t"][0]) * s["pos_rate"]).astype(int).clip(0, len(s["hd"]) - 1)
actual_hd = s["hd"][frame]
ok = ~np.isnan(actual_hd)
errs.append(np.abs(hd.circ_diff(decoded[ok], actual_hd[ok])))
used.append(len(cells_))
return np.concatenate(errs), used
dec_err, dec_sel = [], []
for s in [s for s in sessions if s.get("wake_epoch")]:
if len(hd_list(s)) >= hd.MIN_HD_CELLS_DECODE:
errs, used = decode_session(s, False)
dec_err.append({"session": s["session"], "mouse": s["mouse"], "hd_cells": len(hd_list(s)), "bins": len(errs),
"median_error_deg": np.median(errs), "errors": errs})
res = decode_session(s, True)
if res is not None:
dec_sel.append({"session": s["session"], "cells_per_fold": res[1], "median_error_deg": np.median(res[0]), "errors": res[0]})
dec_err = pd.DataFrame(dec_err)
dec_sel = pd.DataFrame(dec_sel)
pooled = np.median(np.concatenate(dec_err.errors.values))
meets = int((dec_err.median_error_deg <= 30).sum())
print(dec_err.drop(columns="errors").round(1).to_string(index=False))
print(f"\nSessions with at least {hd.MIN_HD_CELLS_DECODE} HD cells: {len(dec_err)}; pooled median absolute error {pooled:.1f} deg; "
f"sessions with median <= 30 deg: {meets} of {len(dec_err)} ({meets / len(dec_err):.0%})")
print(f"Chance: the median absolute error of a uniformly random angle is 90 deg. Median of session medians {dec_err.median_error_deg.median():.1f} deg.")
print("By mouse:", dec_err.groupby("mouse").median_error_deg.median().round(1).to_dict())
th_err = dec_err[dec_err.mouse.isin(THALAMUS_ONLY)]
print(f"Variant, mice that appear thalamus only: {len(th_err)} sessions, median of session medians {th_err.median_error_deg.median():.1f} deg")
sel_pooled = np.median(np.concatenate(dec_sel.errors.values))
print(f"\nVariant, HD cells selected on the training half only: {len(dec_sel)} sessions with at least {hd.MIN_HD_CELLS_DECODE} cells in both folds; "
f"pooled median error {sel_pooled:.1f} deg; session medians {dec_sel.median_error_deg.min():.1f} to {dec_sel.median_error_deg.max():.1f}; "
f"at or below 30 deg: {int((dec_sel.median_error_deg <= 30).sum())} of {len(dec_sel)}")
both = dec_err.merge(dec_sel, on="session", suffixes=("", "_sel"))
print(f"In the {len(both)} sessions in both: median of session medians {both.median_error_deg.median():.1f} deg (whole-epoch selection) "
f"against {both.median_error_deg_sel.median():.1f} deg (training-half selection)")
session mouse hd_cells bins median_error_deg
Mouse12-120806 Mouse12 21 209329 14.4
Mouse12-120807 Mouse12 13 188760 11.4
Mouse12-120809 Mouse12 20 208106 10.5
Mouse12-120810 Mouse12 10 180277 12.7
Mouse17-130128 Mouse17 6 188876 16.6
Mouse17-130129 Mouse17 11 334243 11.3
Mouse17-130130 Mouse17 12 198749 12.7
Mouse17-130131 Mouse17 9 243509 13.2
Mouse17-130203 Mouse17 6 244741 15.0
Mouse20-130517 Mouse20 6 185670 12.6
Mouse24-131216 Mouse24 9 161334 13.4
Mouse24-131217 Mouse24 10 179701 10.6
Mouse24-131218 Mouse24 11 167655 9.9
Mouse25-140123 Mouse25 12 215756 10.4
Mouse25-140124 Mouse25 7 172713 27.2
Mouse25-140129 Mouse25 6 135496 18.4
Mouse25-140130 Mouse25 13 191973 8.1
Mouse25-140131 Mouse25 17 191863 9.1
Mouse25-140203 Mouse25 10 191774 9.4
Mouse25-140204 Mouse25 15 183353 9.2
Mouse25-140205 Mouse25 11 181914 10.1
Mouse25-140206 Mouse25 7 182639 10.1
Mouse28-140310 Mouse28 15 143230 11.1
Mouse28-140311 Mouse28 13 171915 10.4
Mouse28-140312 Mouse28 26 191383 10.7
Mouse28-140313 Mouse28 24 233774 11.0
Mouse28-140317 Mouse28 16 168791 13.2
Mouse28-140318 Mouse28 16 175425 12.9
Sessions with at least 6 HD cells: 28; pooled median absolute error 11.7 deg; sessions with median <= 30 deg: 28 of 28 (100%)
Chance: the median absolute error of a uniformly random angle is 90 deg. Median of session medians 11.2 deg.
By mouse: {'Mouse12': 12.0, 'Mouse17': 13.2, 'Mouse20': 12.6, 'Mouse24': 10.6, 'Mouse25': 10.1, 'Mouse28': 11.0}
Variant, mice that appear thalamus only: 10 sessions, median of session medians 12.7 deg
Variant, HD cells selected on the training half only: 28 sessions with at least 6 cells in both folds; pooled median error 11.7 deg; session medians 8.1 to 24.7; at or below 30 deg: 28 of 28
In the 28 sessions in both: median of session medians 11.2 deg (whole-epoch selection) against 11.2 deg (training-half selection)
Result: reproduced. In the 28 sessions with at least 6 HD cells, the cross-validated decoder has a pooled median absolute error of 11.7 degrees (chance 90), and every session is at or below 30 degrees (8.1 to 27.2). Selecting the HD cells on the training half alone keeps the same 28 sessions and gives a pooled median of 11.7 degrees (session medians 8.1 to 24.7).
Summary¶
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.
| Claim | Outcome |
|---|---|
| hd-cell-fraction | partially reproduced (378 HD cells; fraction 0.218 against 0.328; 391 cells with Mouse12-120808 realigned) |
| fig2a-rates-preserved | reproduced (r = 0.838 SWS, 0.891 REM) |
| fig2b-pair-correlations-preserved | reproduced (r = 0.938 SWS, 0.963 REM) |
| fig3a-packet-60-degrees | reproduced (54.2, 56.9, 55.9 degrees) |
| fig2f-sws-tenfold-faster | partially reproduced (median ratio 8.33; the estimator gives 4.99 for a known tenfold, so the fold is not determined) |
| fig2f-rem-similar-to-wake | reproduced (median ratio 0.87) |
| fig4b-noise-correlations | reproduced |
| decoding-tracks-hd | reproduced (11.7 degrees) |
| 13 claims comparing ADn and PoS | blocked, metadata insufficient |
How the analysis differs from the paper's. HD cells of both structures are pooled, as the units carry no region. Waking is the open-field epoch found from position, as the files do not mark it. Choices the Methods leave open (the von Mises fit, the Rayleigh test, the scale between expected and observed cross-correlograms, the lag range, the decoder's angle bins, the measure of "similar" and "faithfully") are fixed in the claims file. Pairs on the same shank cannot be removed. Two sessions lack processed position and one lacks sleep states. Poisson surrogates driven by the measured head direction are used to check the bias of the angular velocity estimator and of the noise correlation measure; the paper reports no such check.
What was found about the dandiset. The units table holds only spike times, with no electrode, shank, region or waveform, and every electrode's location is "unknown", so the paper's comparisons of ADn and PoS cannot be made. The files hold 1959 unit rows (103 without spikes) in 40 sessions against the paper's 1077 cells in 42, with no flag for which units the paper used; unit ids are 0 to N-1 and carry no shank or cluster number, and the electrode groups give no location. The open-field epoch is not marked. In 36 files one unit's spike times are two sorted runs joined at an index of the form k x 65536 - 1. In Mouse12-120808 position and spikes drift apart by 1.5 to 3.4 s during the exploration; realigned, the session has 13 HD cells, and as deposited none. The raw LED positions are stated in meters with a NaN conversion but hold pixels. The annotation series have no description, and the dandiset's link to the paper names it "nternally organized..." with the title in place of the DOI.
In [26]:
print(f"Total run time {(time.time() - started_all) / 60:.1f} min")
Total run time 6.5 min