Reproducing Duszkiewicz et al. (2024) from dandiset 000939¶
Paper. Duszkiewicz, Orhan, Skromne Carrasco, Brown, Owczarek, Vite, Wood and Peyrache
(2024), "Local origin of excitatory-inhibitory tuning equivalence in a cortical network",
Nature Neuroscience 27, 782-792, doi:10.1038/s41593-024-01588-5. The claims file
(claims/000939/10.1038-s41593-024-01588-5.yaml) was written from the version of record
(PMC11001581) with its Methods, Extended Data legends and Supplementary Information, and from
the paper's MATLAB code on figshare (doi:10.6084/m9.figshare.24921252, version 3).
Dataset. DANDI:000939, version 0.261016.2310 (doi:10.48324/dandi.000939/0.261016.2310), 31 NWB files, one PoSub session per mouse, streamed from the archive. The paper cites no version of the dandiset (its data statement names the figshare item), so the most recent published version is pinned. Its 31 files hold the same sessions as version 0.250207.0025, with the raw voltage added in 2026; the data check below shows that the unit, spike and head-direction sample counts of every file are the same in both versions.
Outside data. Claims marked data: elsewhere in the claims file use the figshare item
above, version 3: Dataset_5 (ADN recordings), Dataset_6 (optogenetic sessions with light over
ADN), Dataset_4 (Neuropixel session) and Dataset_3 (cue rotation), each streamed and read in
memory. Three uses serve dandiset claims without deciding them: the code's session list
(SessionLists_ForAnalysis.zip) and Dataset_4 identify the session the dandiset lacks
(cell-classes-counts), and the three PoSub ArchT mice of Dataset_6 that the dandiset lacks are
added in a variant of the optogenetic claims.
| Claim | Tested here |
|---|---|
| cell-classes-counts | tested |
| ed2b-excitatory-hd-information | tested |
| fig2a-fs-hd-information | tested |
| fig2b-decoding-from-fs | tested |
| fig2c-geometry-independent | tested |
| ed2hi-within-epoch-stability | tested |
| fig2f-cue-rotation | tested, figshare Dataset_3 |
| fig3c-adn-broader | tested, figshare Dataset_5 |
| fig3-ten-components-power | tested |
| fig3e-fourier-signatures-differ | tested, figshare Dataset_5 |
| fig3f-posub-hd-closer-to-fs | tested, figshare Dataset_5 |
| ed3k-robust-to-controls | tested, figshare Dataset_5 |
| fig3g-first-three-components | tested |
| fig3h-dominant-component | tested |
| ed3g-fs-spectra-heterogeneous | tested (ADN part reported) |
| ed3l-fs-phases-independent | tested |
| ed4a-autocorrelogram-symmetries | tested |
| ed4cd-isomap-posub-central | tested, figshare Dataset_5 |
| fig4-random-connectivity-simulation | blocked: simulation code and parameters not available |
| fig4f-synaptic-weights-dispersed | tested |
| fig4h-model-fit | blocked: fitting code not available, input numbers and initial values not given |
| ed7-adn-disinhibition | tested, figshare Dataset_6 |
| fig5b-posub-rates-increase | tested |
| ed7gj-tuning-preserved | tested |
| fig5e-hd-multiplicative-fs-additive | tested |
| ed8cd-multiplicative-detectable | tested |
| ed8e-adn-mixed-gain | tested, figshare Dataset_6 |
| fig6c-coupling-preserved-wake-rem | tested |
| ed9d-coupling-follows-tuning | tested |
| fig7b-wake-isomap-tuning | tested |
| ed10b-ring-topology | tested |
| fig7d-rem-isomap-tuning | tested |
| fig1ac-posub-border | tested, figshare Dataset_4 |
| ed6-trn-innervates-adn | blocked: anatomical images not deposited |
The notebook streams about 2.7 GB of processed data from the archive and 4 GB of tarballs from figshare; it takes about 37 minutes, fifteen of them in the GLM fits.
Data check¶
import time
import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
import posub
T0 = time.time()
print("dandiset", posub.DANDISET, "version", posub.VERSION)
S = posub.load_all()
print(f"streamed {len(S)} files in {time.time() - T0:.0f} s")
rows = []
for s in S:
sq = s["epochs"]["wake_square"][0]
rem = (s["sleep"].get("rem", np.zeros((0, 2)))[:, 1] - s["sleep"].get("rem", np.zeros((0, 2)))[:, 0]).sum()
rows.append(dict(subject=s["subject"], session=s["session_id"], orientation=s["orientation"],
units=len(s["spikes"]), epochs=",".join(s["epochs"]),
square_s=round(sq[1] - sq[0]), rem_s=round(float(rem)),
pulses=len(s.get("pulses", [])), virus=s["virus"][:40], genotype=s["genotype"][12:30]))
meta = pd.DataFrame(rows)
pd.set_option("display.width", 250, "display.max_columns", 20, "display.max_rows", 100)
print(meta.to_string())
print("units by orientation:", meta.groupby("orientation").units.agg(["count", "sum"]).to_dict())
# units and ranges of the variables used
hd_all = np.concatenate([s["hd"] for s in S]); pos_all = np.concatenate([s["pos"] for s in S])
print("head direction unit:", {s["hd_unit"] for s in S}, "range", np.nanmin(hd_all).round(3), np.nanmax(hd_all).round(3),
"NaN fraction", np.isnan(hd_all).mean().round(4))
print("position unit:", {s["pos_unit"] for s in S}, "range x", np.nanmin(pos_all[:, 0]).round(1), np.nanmax(pos_all[:, 0]).round(1),
"y", np.nanmin(pos_all[:, 1]).round(1), np.nanmax(pos_all[:, 1]).round(1))
dts = np.array([np.median(np.diff(s["hd_t"])) for s in S])
print("HD sampling interval (s), min and max over files:", dts.min().round(5), dts.max().round(5))
print("spike time resolution attribute:", {s["spike_res"] for s in S})
ttp = np.concatenate([s["ttp"] for s in S])
print("trough-to-peak (ms): n", len(ttp), "NaN", np.isnan(ttp).sum(), "range", np.nanmin(ttp), np.nanmax(ttp))
# do spikes and tracking share a clock: spikes fall inside the recording, tracking inside its epochs
for s in S:
last = max(sp.max() for sp in s["spikes"] if len(sp)); end = max(b for v in s["epochs"].values() for a, b in v)
assert last <= end + 1, s["subject"]
print("all spike times fall within the session's epochs")
# HD sampled only during exploration epochs?
for s in S:
fin = np.isfinite(s["hd"]); t = s["hd_t"][fin]
wake = [iv for k, v in s["epochs"].items() if k.startswith("wake") for iv in v]
s["_hd_out"] = (~posub.in_intervals(t, wake)).mean()
print("fraction of tracked HD samples outside wake epochs, max over files:", max(s["_hd_out"] for s in S).round(4))
# duplicate check: identical spike trains across files
sig = {}
for s in S:
for i, sp in enumerate(s["spikes"]):
key = (len(sp), round(float(sp[0]), 5) if len(sp) else None, round(float(sp[-1]), 5) if len(sp) else None)
sig.setdefault(key, []).append((s["subject"], i))
dups = {k: v for k, v in sig.items() if len(v) > 1}
print("units with identical spike count, first and last spike in two places:", len(dups))
hd_sig = {}
for s in S:
hd_sig.setdefault((len(s["hd"]), round(float(np.nansum(s["hd"])), 3)), []).append(s["subject"])
print("files with identical HD series:", [v for v in hd_sig.values() if len(v) > 1])
# the previous published version with the same sessions (0.250207.0025, processed data only):
# units and spike counts per file, against the pinned version
import h5py, remfile
from concurrent.futures import ThreadPoolExecutor
def _units(a):
f = h5py.File(remfile.File(posub.s3_url(a["asset_id"])), "r")
return a["path"].split("/")[0], len(f["units/id"]), f["units/spike_times"].shape[0], f["processing/behavior/CompassDirection/head-direction/data"].shape[0]
with ThreadPoolExecutor(8) as ex:
old = {k: (n, ns, nh) for k, n, ns, nh in ex.map(_units, posub.assets("0.250207.0025"))}
new = {s["path"].split("/")[0]: (len(s["spikes"]), sum(len(x) for x in s["spikes"]), len(s["hd"])) for s in S}
print("subjects in 0.250207.0025 and the pinned version:", len(old), len(new), "; same subjects:", set(old) == set(new))
print("files whose unit count, spike count or HD sample count differ between the versions:", [k for k in new if old.get(k) != new[k]])
# REM intervals (intervals/sleep_states): the same in 0.250207.0025 (intervals/rem there), and
# where they fall. The 2026 version says it fixed the indexing of sleep scoring epochs.
def _rem_old(a):
f = h5py.File(remfile.File(posub.s3_url(a["asset_id"])), "r")
if "rem" not in f["intervals"]:
return a["path"].split("/")[0], np.zeros((0, 2))
return a["path"].split("/")[0], np.c_[f["intervals/rem/start_time"][:], f["intervals/rem/stop_time"][:]]
with ThreadPoolExecutor(8) as ex:
rem_old = dict(ex.map(_rem_old, posub.assets("0.250207.0025")))
rem_rows = []
for s in S:
k = s["path"].split("/")[0]
r = s["sleep"].get("rem", np.zeros((0, 2)))
end = max(b for v in s["epochs"].values() for a, b in v)
smax = max(float(v.max()) for v in s["sleep"].values() if len(v))
tt = np.arange(0, max(end, smax), 0.1)
in_r = posub.in_intervals(tt, r)
in_home = in_r & posub.in_intervals(tt, s["epochs"]["home_cage"])
in_expl = in_r & posub.in_intervals(tt, [iv for t_, v in s["epochs"].items() if t_ != "home_cage" for iv in v])
o = rem_old[k]
rem_rows.append(dict(subject=s["subject"], session_end_s=round(end), sleep_states_end_s=round(smax), rem_s=round(float((r[:, 1] - r[:, 0]).sum())),
rem_in_home_cage_s=round(in_home.sum() * 0.1), rem_in_exploration_s=round(in_expl.sum() * 0.1),
rem_after_end_s=round((in_r & (tt >= end)).sum() * 0.1),
same_in_0_250207=bool(o.shape == r.shape and np.allclose(o, r))))
rem_tab = pd.DataFrame(rem_rows)
print(rem_tab.to_string())
print("files whose sleep states run past the session's end:", list(rem_tab[rem_tab.sleep_states_end_s > rem_tab.session_end_s + 5].subject),
"; with REM inside an exploration epoch:", list(rem_tab[rem_tab.rem_in_exploration_s > 0].subject))
# the exploration epochs that hold REM, with the seconds of REM in each
for s in S:
r = s["sleep"].get("rem", np.zeros((0, 2)))
for tag, v in s["epochs"].items():
ov = sum(max(0.0, min(b, y) - max(a, x)) for a, b in r for x, y in v)
if tag != "home_cage" and ov > 0:
print(f"{s['subject']}: {ov:.1f} s of REM inside the epoch '{tag}' {[tuple(np.round(e, 1)) for e in v]}")
# REM is used only inside the home-cage epochs (the sleep recordings, Methods); the rest is dropped
for s in S:
r = s["sleep"].get("rem", np.zeros((0, 2)))
clip = [(max(a, x), min(b, y)) for a, b in r for x, y in s["epochs"]["home_cage"] if min(b, y) > max(a, x)]
s["sleep"]["rem"] = np.array(clip, float).reshape(-1, 2)
print("REM kept (s) per session, inside the home-cage epochs:", {s["subject"]: round(float((s["sleep"]["rem"][:, 1] - s["sleep"]["rem"][:, 0]).sum())) for s in S})
dandiset 000939 version 0.261016.2310
streamed 31 files in 146 s
subject session orientation units epochs square_s rem_s pulses virus genotype
0 A1808 190807 vertical 101 home_cage,wake_square,wake_opto 912 39 240 Bilateral AAV9-emCBA-Flex-GFP virus inje VGAT-IRES-Cre mou
1 A1813 190408 vertical 70 home_cage,wake_square,wake_opto 1042 389 240 Bilateral AAV9-CAG-flex-ArchT-GFP virus VGAT-IRES-Cre mou
2 A1815 190523 vertical 65 home_cage,wake_square,wake_opto 907 285 240 Bilateral AAV9-CAG-flex-ArchT-GFP virus VGAT-IRES-Cre mou
3 A1821 190910 vertical 74 home_cage,wake_square,wake_opto 2707 324 240 Bilateral AAV9-emCBA-Flex-GFP virus inje VGAT-IRES-Cre mou
4 A1824 190827 vertical 70 home_cage,wake_square,wake_opto 1208 88 240 Bilateral AAV9-emCBA-Flex-GFP virus inje VGAT-IRES-Cre mou
5 A3701 191119 parallel 102 home_cage,wake_square,wake_triangle 2541 308 0 nan VGAT-IRES-Cre mou
6 A3702 191126 parallel 90 home_cage,wake_square,wake_triangle 2552 124 0 nan VGAT-IRES-Cre mou
7 A3703 191215 parallel 103 home_cage,wake_square,wake_triangle 2417 301 0 nan VGAT-IRES-Cre mou
8 A3705 200306 parallel 117 home_cage,wake_square,wake_triangle 2639 620 0 nan VGAT-IRES-Cre mou
9 A3706 200313 parallel 126 home_cage,wake_square,wake_triangle 2407 218 0 nan VGAT-IRES-Cre mou
10 A3707 200317 parallel 185 home_cage,wake_square,wake_triangle 1808 38 0 nan VGAT-IRES-Cre mou
11 A3708 200317b parallel 78 home_cage,wake_square,wake_triangle 1904 501 0 nan VGAT-IRES-Cre mou
12 A3709 200601 parallel 64 home_cage,wake_square 1886 978 0 nan VGAT-IRES-Cre mou
13 A3710 200609 parallel 107 home_cage,wake_square,wake_triangle 2415 288 0 nan VGAT-IRES-Cre mou
14 A3711 200810b parallel 142 home_cage,wake_square,wake_triangle 2409 591 0 nan VGAT-IRES-Cre mou
15 A3712 200903a parallel 126 home_cage,wake_square,wake_triangle 1807 248 0 nan VGAT-IRES-Cre mou
16 A3713 200909a parallel 88 home_cage,wake_square,wake_triangle 1822 200 0 nan VGAT-IRES-Cre mou
17 A3716 201015b parallel 147 home_cage,wake_square,wake_triangle 1826 284 0 nan VGAT-IRES-Cre mou
18 A3717 201021 parallel 98 home_cage,wake_square,wake_triangle 1808 1040 0 nan VGAT-IRES-Cre mou
19 A3723 201115 parallel 68 home_cage,wake_square,wake_triangle 2407 480 0 nan VGAT-IRES-Cre mou
20 A3728 210309b parallel 42 home_cage,wake_square,wake_triangle 1809 0 0 nan VGAT-IRES-Cre mou
21 A3730 210323b parallel 77 home_cage,wake_square,wake_triangle 2408 228 0 nan VGAT-IRES-Cre mou
22 A5505 200831 vertical 55 home_cage,wake_square,wake_triangle 1809 355 0 nan VGAT-IRES-Cre mou
23 A5506 200914a vertical 49 home_cage,wake_square,wake_triangle 1809 340 0 nan VGAT-IRES-Cre mou
24 A5507 200920 vertical 96 home_cage,wake_square,wake_triangle 1824 380 0 nan VGAT-IRES-Cre mou
25 A5508 200930 vertical 52 home_cage,wake_square,wake_triangle 1816 122 0 nan VGAT-IRES-Cre mou
26 A5801 200924 vertical 46 home_cage,wake_square 2425 139 0 nan VGAT-IRES-Cre mou
27 A6211 210407 vertical 61 home_cage,wake_square,wake_opto 1212 214 240 Unilateral AAV2/8-CAG-Flex-ArchT-eGFP vi Pvalb-IRES-Cre mo
28 A6215 210610 vertical 56 home_cage,wake_square,wake_opto 973 248 240 Unilateral AAV2/8-CAG-Flex-ArchT-eGFP vi SST-IRES-Cre mous
29 A6216 210623 vertical 75 home_cage,wake_square,wake_opto 1208 429 240 Unilateral AAV2/8-CAG-Flex-ArchT-eGFP vi SST-IRES-Cre mous
30 A9701 210805 vertical 61 home_cage,wake_square,wake_opto 1222 313 640 Unilateral AAV2/9-Ef1a-DIO-hChR2(H143R)- Pvalb-IRES-Cre mo
units by orientation: {'count': {'parallel': 17, 'vertical': 14}, 'sum': {'parallel': 1760, 'vertical': 931}}
head direction unit: {'radians'} range 0.0 6.283 NaN fraction 0.0095
position unit: {'centimeters'} range x 0.0 92.0 y 0.0 96.4
HD sampling interval (s), min and max over files: 0.01 0.01
spike time resolution attribute: {np.float64(5e-05)}
trough-to-peak (ms): n 2691 NaN 177 range 0.1250000000000001 1.17
all spike times fall within the session's epochs
fraction of tracked HD samples outside wake epochs, max over files: 0.0 units with identical spike count, first and last spike in two places: 0 files with identical HD series: []
subjects in 0.250207.0025 and the pinned version: 31 31 ; same subjects: True files whose unit count, spike count or HD sample count differ between the versions: []
subject session_end_s sleep_states_end_s rem_s rem_in_home_cage_s rem_in_exploration_s rem_after_end_s same_in_0_250207
0 A1808 10199 10197 39 39 0 0 True
1 A1813 11667 11665 389 389 0 0 True
2 A1815 9623 9621 285 285 0 0 True
3 A1821 11552 11550 324 267 57 0 True
4 A1824 9105 9103 88 88 0 0 True
5 A3701 14118 14116 308 308 0 0 True
6 A3702 12640 12638 124 124 0 0 True
7 A3703 13327 13325 301 301 0 0 True
8 A3705 16199 16197 620 620 0 0 True
9 A3706 10908 10906 218 218 0 0 True
10 A3707 11142 11140 38 38 0 0 True
11 A3708 13658 13658 501 501 0 0 True
12 A3709 13508 13506 978 978 0 0 True
13 A3710 12602 12600 288 288 0 0 True
14 A3711 12177 24351 591 293 0 298 True
15 A3712 11290 11288 248 236 12 0 True
16 A3713 11233 11231 200 200 0 0 True
17 A3716 11728 11726 284 284 0 0 True
18 A3717 12269 24536 1040 522 0 518 True
19 A3723 12241 24481 480 257 0 223 True
20 A3728 11271 11269 0 0 0 0 True
21 A3730 12455 12453 228 228 0 0 True
22 A5505 11190 11188 355 355 0 0 True
23 A5506 11193 11191 340 340 0 0 True
24 A5507 12133 12131 380 380 0 0 True
25 A5508 11228 11226 122 122 0 0 True
26 A5801 6114 6112 139 139 0 0 True
27 A6211 11164 11162 214 214 0 0 True
28 A6215 11258 11256 248 248 0 0 True
29 A6216 11159 11157 429 429 0 0 True
30 A9701 11387 11385 313 313 0 0 True
files whose sleep states run past the session's end: ['A3711', 'A3717', 'A3723'] ; with REM inside an exploration epoch: ['A1821', 'A3712']
A1821: 16.0 s of REM inside the epoch 'wake_square' [(np.float64(3042.6), np.float64(5749.6))]
A1821: 41.0 s of REM inside the epoch 'wake_opto' [(np.float64(8831.8), np.float64(11551.8))]
A3701: 0.3 s of REM inside the epoch 'wake_triangle' [(np.float64(11703.7), np.float64(14118.3))]
A3712: 12.0 s of REM inside the epoch 'wake_square' [(np.float64(3764.9), np.float64(5571.5))]
REM kept (s) per session, inside the home-cage epochs: {'A1808': 39, 'A1813': 389, 'A1815': 285, 'A1821': 267, 'A1824': 88, 'A3701': 308, 'A3702': 124, 'A3703': 301, 'A3705': 620, 'A3706': 218, 'A3707': 38, 'A3708': 501, 'A3709': 978, 'A3710': 288, 'A3711': 293, 'A3712': 236, 'A3713': 200, 'A3716': 284, 'A3717': 522, 'A3723': 257, 'A3728': 0, 'A3730': 228, 'A5505': 355, 'A5506': 340, 'A5507': 380, 'A5508': 122, 'A5801': 139, 'A6211': 214, 'A6215': 248, 'A6216': 429, 'A9701': 313}
Tuning curves of every unit¶
Computed once for the analyses below: square-arena curves (1 degree bins, 3 degree smoothing, speed above 2 cm/s), their time-reversed controls, the two halves of the epoch, the triangle-arena curves, the 2 degree smoothing variant (code), and curves at 10 cm/s (Extended Data Fig. 3k). The time-reversed control reverses the HD series in time over the used span (the tracked square epoch, or its half) before the speed mask is applied (claims file). Rates are spike counts over the whole square epoch (code, GetMeanRates: Rate(S, wake1Ep)).
U = [] # one row per unit
TC = {} # per session: dict of arrays (360, n units)
for si, s in enumerate(S):
t, ang = s["hd_t"], s["hd"]
a, b = s["epochs"]["wake_square"][0]
n = len(s["spikes"])
rate = np.array([((sp >= a) & (sp < b)).sum() / (b - a) for sp in s["spikes"]])
m = posub.good_samples(s, "wake_square")
d = {}
for sd in (3, 2):
d[f"sq{sd}"], occ = posub.tuning_curves(s["spikes"], t, ang, m, sd=sd)
d[f"sq{sd}_rev"], occ_r = posub.tuning_curves(s["spikes"], t, ang, m, sd=sd, reverse=True)
m1, m2 = posub.split_halves(t, m)
d[f"h1_{sd}"], _ = posub.tuning_curves(s["spikes"], t, ang, m1, sd=sd)
d[f"h2_{sd}"], _ = posub.tuning_curves(s["spikes"], t, ang, m2, sd=sd)
d[f"h1_{sd}_rev"], _ = posub.tuning_curves(s["spikes"], t, ang, m1, sd=sd, reverse=True)
d[f"h2_{sd}_rev"], _ = posub.tuning_curves(s["spikes"], t, ang, m2, sd=sd, reverse=True)
if sd == 3:
info = posub.hd_info(d["sq3"], occ)
info_rev = posub.hd_info(d["sq3_rev"], occ_r)
move_s = m.sum() * np.median(np.diff(t))
m10 = posub.good_samples(s, "wake_square", speed_th=10.0)
d["sq3_v10"], _ = posub.tuning_curves(s["spikes"], t, ang, m10)
rate_tri = np.full(n, np.nan)
if "wake_triangle" in s["epochs"]:
ta, tb = s["epochs"]["wake_triangle"][0]
rate_tri = np.array([((sp >= ta) & (sp < tb)).sum() / (tb - ta) for sp in s["spikes"]])
mt = posub.good_samples(s, "wake_triangle")
for sd in (3, 2):
d[f"tri{sd}"], _ = posub.tuning_curves(s["spikes"], t, ang, mt, sd=sd)
d[f"tri{sd}_rev"], _ = posub.tuning_curves(s["spikes"], t, ang, mt, sd=sd, reverse=True)
TC[si] = d
for i in range(n):
U.append(dict(sess=si, subject=s["subject"], unit=i, orientation=s["orientation"], rate=rate[i],
rate_tri=rate_tri[i], ttp=s["ttp"][i], info=info[i], info_rev=info_rev[i],
lab_ex=s["lab_ex"][i], lab_fs=s["lab_fs"][i], lab_hd=s["lab_hd"][i], move_s=move_s))
U = pd.DataFrame(U)
U["gd"] = (U.rate >= posub.RATE_MIN) & np.isfinite(U.ttp)
print(f"{len(U)} units; {U.gd.sum()} viable (rate at least 0.5 Hz in the square epoch and a trough-to-peak, i.e. a negative waveform)")
print("time moving above 2 cm/s in the square epoch (s), min and max over sessions:",
U.groupby("sess").move_s.first().min().round(0), U.groupby("sess").move_s.first().max().round(0))
print(f"tuning curves done at {time.time() - T0:.0f} s")
2691 units; 2420 viable (rate at least 0.5 Hz in the square epoch and a trough-to-peak, i.e. a negative waveform) time moving above 2 cm/s in the square epoch (s), min and max over sessions: 578.0 2094.0 tuning curves done at 288 s
cell-classes-counts¶
PoSub units split into 1,835 putative excitatory cells and 427 putative FS cells by firing rate and waveform, and 1,602 excitatory cells (87%) were HD cells with HD information above 0.2 bits per spike; the vertical recordings gave 931 units from 14 mice and the parallel ones 1,999 units from 18 mice.
Classes as the claims file fixes them (Methods): viable units (rate at least 0.5 Hz in the square epoch, negative waveform, read as a finite trough-to-peak), FS below 0.4 ms and above 10 Hz, excitatory above 0.4 ms and below 10 Hz; the code's 0.35 ms cut (with its >= 10 Hz and <= 0.35 ms) is a variant. HD cells are excitatory cells above 0.2 bits per spike, and, as a second reading, above the 99th percentile of the excitatory cells' time-reversed HD information.
The parallel-probe count is compared with the code's list of the paper's 32 PoSub sessions (figshare item, SessionLists_ForAnalysis.zip, dataset_All.list), and the session the dandiset lacks is counted from figshare Dataset_4.
classify = posub.classify
U["ex"], U["fs"] = classify(U)
U["hd"] = U.ex & (U["info"] > posub.HD_INFO_TH)
p99 = np.percentile(U["info_rev"][U.ex], 99)
rows = []
for name, cut, code in [("Methods, 0.4 ms", 0.4, False), ("code, 0.35 ms", 0.35, True)]:
ex, fs = classify(U, cut, code)
for thname, th in [("0.2 bits/spike", 0.2), (f"99th pct of reversed ({p99:.3f})", p99)]:
hd = ex & (U["info"] > th)
rows.append(dict(classes=name, hd_threshold=thname, excitatory=int(ex.sum()), FS=int(fs.sum()),
HD=int(hd.sum()), HD_fraction=round(hd.sum() / ex.sum(), 3)))
print(pd.DataFrame(rows).to_string())
print("units by orientation:", U.groupby("orientation").size().to_dict(),
"sessions:", U.groupby("orientation").sess.nunique().to_dict())
print("dandiset's own labels: excitatory", int(U.lab_ex.sum()), "FS", int(U.lab_fs.sum()), "HD", int(U.lab_hd.sum()))
print("agreement with the dandiset labels (main reading): ex", (U.ex == U.lab_ex).mean().round(4),
"fs", (U.fs == U.lab_fs).mean().round(4), "hd", (U.hd == U.lab_hd).mean().round(4))
exc, fsc = classify(U, 0.35, True)
print("agreement with the dandiset labels (code reading): ex", (exc == U.lab_ex).mean().round(4),
"fs", (fsc == U.lab_fs).mean().round(4), "hd", ((exc & (U["info"] > 0.2)) == U.lab_hd).mean().round(4))
print("units neither excitatory nor FS among viable ones (main reading):", int((U.gd & ~U.ex & ~U.fs).sum()))
print(U.groupby("subject")[["ex", "fs", "hd"]].sum().T.to_string())
fig, ax = plt.subplots(1, 2, figsize=(9, 3.5))
g = U[U.gd]
ax[0].scatter(g.ttp, g.rate, s=3, c=np.where(g.fs, "tab:orange", np.where(g.ex, "tab:purple", "grey")))
ax[0].set_yscale("log"); ax[0].set_xlabel("trough-to-peak (ms)"); ax[0].set_ylabel("rate in square (Hz)")
ax[1].hist(np.log10(U["info"][U.gd].clip(1e-4)), 80); ax[1].set_xlabel("log10 HD information (bits/spike), viable units")
plt.tight_layout(); plt.show()
print("plotted ranges: ttp", g.ttp.min(), g.ttp.max(), "rate", g.rate.min().round(2), g.rate.max().round(1),
"; units with information below 1e-4 drawn at -4:", int((U["info"][U.gd] < 1e-4).sum()))
import io, requests
# Which of the paper's PoSub sessions is missing: the code's session list (figshare item,
# SessionLists_ForAnalysis.zip, dataset_All.list, the 32 PoSub sessions) against the dandiset
import zipfile
zl = zipfile.ZipFile(io.BytesIO(requests.get(f"https://ndownloader.figshare.com/files/{posub.FIGSHARE_LISTS_FILE}", timeout=120).content))
paper_list = [x.strip() for x in zl.read("ListFiles/dataset_All.list").decode().splitlines() if x.strip()]
missing = [x for x in paper_list if x.split("-")[0] not in set(U.subject)]
print(f"dataset_All.list: {len(paper_list)} sessions; in the list and not in the dandiset: {missing}; in the dandiset and not in the list:",
sorted(set(U.subject) - {x.split('-')[0] for x in paper_list}))
# the missing session is the Neuropixel recording of figshare Dataset_4 (also used for fig1ac)
NP = posub.load_adn(file_id=posub.FIGSHARE_NPX_FILE)[0]
print(f"figshare Dataset_4 session {NP['session']}: {len(NP['spikes'])} units; the files' classes: excitatory {int(NP['file_types']['ex'].sum())}, "
f"FS {int(NP['file_types']['fs'].sum())}, HD {int(NP['file_types']['hd'].sum())}")
par = int((U.orientation == "parallel").sum())
print(f"parallel units, dandiset {par} + Neuropixel session {len(NP['spikes'])} = {par + len(NP['spikes'])} (paper 1,999 from 18 mice); "
f"largest dandiset recording {U.groupby('subject').size().max()} units")
classes hd_threshold excitatory FS HD HD_fraction
0 Methods, 0.4 ms 0.2 bits/spike 1661 430 1512 0.910
1 Methods, 0.4 ms 99th pct of reversed (0.272) 1661 430 1481 0.892
2 code, 0.35 ms 0.2 bits/spike 1699 415 1532 0.902
3 code, 0.35 ms 99th pct of reversed (0.272) 1699 415 1501 0.883
units by orientation: {'parallel': 1760, 'vertical': 931} sessions: {'parallel': 17, 'vertical': 14}
dandiset's own labels: excitatory 1698 FS 414 HD 1533
agreement with the dandiset labels (main reading): ex 0.9848 fs 0.9933 hd 0.9892
agreement with the dandiset labels (code reading): ex 0.9989 fs 0.9989 hd 0.9967
units neither excitatory nor FS among viable ones (main reading): 329
subject A1808 A1813 A1815 A1821 A1824 A3701 A3702 A3703 A3705 A3706 A3707 A3708 A3709 A3710 A3711 A3712 A3713 A3716 A3717 A3723 A3728 A3730 A5505 A5506 A5507 A5508 A5801 A6211 A6215 A6216 A9701
ex 32 31 35 38 40 66 63 67 70 90 123 47 24 72 95 80 37 104 60 48 34 54 33 27 60 38 29 36 39 46 43
fs 11 17 16 16 17 12 16 20 25 21 32 12 8 13 21 23 13 17 18 5 5 9 10 11 10 6 9 11 7 14 5
hd 30 31 33 37 39 66 61 65 66 86 117 43 22 72 89 72 25 91 57 37 27 46 30 24 42 33 26 35 29 44 37
plotted ranges: ttp 0.135 1.17 rate 0.5 97.8 ; units with information below 1e-4 drawn at -4: 0
dataset_All.list: 32 sessions; in the list and not in the dandiset: ['A8603-210602']; in the dandiset and not in the list: []
figshare Dataset_4 session A8603-210602: 239 units; the files' classes: excitatory 137, FS 13, HD 69 parallel units, dandiset 1760 + Neuropixel session 239 = 1999 (paper 1,999 from 18 mice); largest dandiset recording 185 units
Result: reproduced. Under the Methods' classes there are 1,661 excitatory cells (expected 1,835 +- 184), 430 FS cells (427 +- 43), 1,512 HD cells (1,602 +- 160) and an HD fraction of 0.910 (0.87 or 0.85 +- 0.087), all within tolerance, as are the code's classes (1,699, 415, 1,532, 0.902); with the HD threshold read as the 99th percentile of the time-reversed control (0.272 bits per spike) there are 1,481 and 1,501 HD cells. The vertical recordings hold 931 units from 14 mice, as the paper gives. The parallel recordings hold 1,760 units from 17 mice, below 1,999 - 200. The session missing from the dandiset is A8603-210602, the only session of the code's list of 32 (dataset_All.list) that the dandiset lacks: the Neuropixel recording of figshare Dataset_4, with 239 units, and 1,760 + 239 = 1,999, the paper's number. So the shortfall is one session's units, the case the criterion foresees; that session is larger than any dandiset recording (185 units at most), which is why it exceeds the H5 sessions' range. The dandiset's own labels agree with the code's classes for 99.9% of units for the excitatory and FS labels and 99.7% for the HD label.
Statistics used below¶
signrank (helper) is the two-sided Wilcoxon signed-rank test with its normal-approximation z (MATLAB's
signrank reports this z); ranksum (helper) is the two-sided Mann-Whitney U test with its z. All tests
are two-tailed (Methods).
signrank, ranksum, show, classify = posub.signrank, posub.ranksum, posub.show, posub.classify
ed2b-excitatory-hd-information¶
Putative excitatory PoSub cells carry more HD information than their time-reversed controls, and HD information across all PoSub cells is bimodal.
Log10 HD information of each excitatory cell against its time-reversed control, Wilcoxon signed-rank (criterion; code Figure_BasicProperties).
lg = lambda v: np.log10(np.clip(v, 1e-6, None))
for name, (exm, _) in {"Methods classes": (U.ex, U.fs), "code classes": classify(U, 0.35, True)}.items():
show(signrank(lg(U["info"][exm]), lg(U["info_rev"][exm])), f"excitatory, {name}:")
print("paper: n = 1835, z = 37.1")
print("per session, fraction of excitatory cells above their own control:",
U[U.ex].assign(up=lambda d: d["info"] > d["info_rev"]).groupby("subject").up.mean().round(2).describe()[["min", "50%", "max"]].to_dict())
h, e = np.histogram(lg(U["info"][U.gd]), 60)
fig, ax = plt.subplots(figsize=(5, 3)); ax.stairs(h, e); ax.set_xlabel("log10 HD info (bits/spike), all viable units"); plt.show()
# bimodality shown, not judged: the two largest local maxima of the smoothed histogram
from scipy.ndimage import gaussian_filter1d as _g
hs = _g(h.astype(float), 2); pk = [i for i in range(1, len(hs) - 1) if hs[i] >= hs[i - 1] and hs[i] >= hs[i + 1]]
print("local maxima of the smoothed histogram at log10 info =", [round((e[i] + e[i + 1]) / 2, 2) for i in pk], "range of values", e[0].round(2), e[-1].round(2))
excitatory, Methods classes: n=1661, median_a=0.1017, median_b=-1.414, z=-35.3, p=6.989e-273
excitatory, code classes: n=1699, median_a=0.09936, median_b=-1.418, z=-35.7, p=4.496e-279
paper: n = 1835, z = 37.1
per session, fraction of excitatory cells above their own control: {'min': 0.97, '50%': 1.0, 'max': 1.0}
local maxima of the smoothed histogram at log10 info = [np.float64(-1.67), np.float64(0.19)] range of values -3.08 0.52
Result: reproduced. Excitatory cells carry more HD information than their time-reversed controls (median log10 information 0.10 against -1.41, |z| = 35.3, p = 7e-273, n = 1,661; the paper gives z = 37.1 for 1,835); 97 to 100% of the excitatory cells of each session are above their own control. The histogram of HD information over all viable units has two maxima, near log10 = -1.67 and 0.19 (shown, not judged).
fig2a-fs-hd-information¶
The HD information carried by PoSub-FS cell tuning curves is higher than that of their time-reversed controls.
for name, (_, fsm) in {"Methods classes": (U.ex, U.fs), "code classes": classify(U, 0.35, True)}.items():
show(signrank(lg(U["info"][fsm]), lg(U["info_rev"][fsm])), f"FS, {name}:")
print("paper: n = 427, z = 17.9")
print("median HD information of PoSub-HD cells (bits/spike):", U["info"][U.hd].median().round(3),
"; of FS cells:", U["info"][U.fs].median().round(3), "; FS controls:", U["info_rev"][U.fs].median().round(3))
print("FS cells above their own control:", int((U["info"][U.fs] > U["info_rev"][U.fs]).sum()), "of", int(U.fs.sum()))
print(U[U.fs].assign(up=lambda d: d["info"] > d["info_rev"]).groupby("subject").up.mean().round(2).to_dict())
FS, Methods classes: n=430, median_a=-1.699, median_b=-2.633, z=-17.86, p=2.445e-71
FS, code classes: n=415, median_a=-1.699, median_b=-2.646, z=-17.54, p=6.842e-69
paper: n = 427, z = 17.9
median HD information of PoSub-HD cells (bits/spike): 1.348 ; of FS cells: 0.02 ; FS controls: 0.002
FS cells above their own control: 419 of 430
{'A1808': 1.0, 'A1813': 1.0, 'A1815': 1.0, 'A1821': 0.88, 'A1824': 1.0, 'A3701': 1.0, 'A3702': 1.0, 'A3703': 1.0, 'A3705': 1.0, 'A3706': 1.0, 'A3707': 0.97, 'A3708': 1.0, 'A3709': 1.0, 'A3710': 0.92, 'A3711': 0.9, 'A3712': 1.0, 'A3713': 0.92, 'A3716': 0.94, 'A3717': 0.94, 'A3723': 1.0, 'A3728': 1.0, 'A3730': 1.0, 'A5505': 1.0, 'A5506': 0.91, 'A5507': 1.0, 'A5508': 1.0, 'A5801': 1.0, 'A6211': 1.0, 'A6215': 1.0, 'A6216': 0.93, 'A9701': 1.0}
Result: reproduced. FS cells carry more HD information than their time-reversed controls (median log10 information -1.70 against -2.63, |z| = 17.86, p = 2e-71, n = 430; paper z = 17.9 for 427); 419 of 430 FS cells are above their own control. The median HD information of HD cells is 1.35 bits per spike, that of FS cells 0.02.
ed2hi-within-epoch-stability¶
Tuning curves of PoSub-FS and PoSub-HD cells computed from the two halves of the square-arena epoch are correlated above their time-reversed controls.
Zero-lag Pearson correlation of the two halves' curves of each cell, against that of the time-reversed halves (criterion; code GetHdTuning_xval and Figure_Stability). Smoothing 3 degrees (Methods) and 2 degrees (code).
def per_unit(key):
"""Stack a per-session curve array over all units, in the order of U."""
return np.concatenate([TC[si][key] for si in range(len(S))], axis=1)
for sd in (3, 2):
r = posub.zero_lag_r(per_unit(f"h1_{sd}"), per_unit(f"h2_{sd}"))
rr = posub.zero_lag_r(per_unit(f"h1_{sd}_rev"), per_unit(f"h2_{sd}_rev"))
U[f"r_half{sd}"], U[f"r_half{sd}_rev"] = r, rr
for cls in ("fs", "hd"):
show(signrank(r[U[cls]], rr[U[cls]]), f"{cls.upper()}, smoothing {sd} deg:")
print("paper: FS n = 427, z = 17.8; HD n = 1602, z = 34.7")
print("cells with an undefined correlation (a flat half):", int(np.isnan(U.r_half3[U.fs | U.hd]).sum()))
print("median per session of FS half-half r:", U[U.fs].groupby("subject").r_half3.median().round(2).describe()[["min", "50%", "max"]].round(2).to_dict())
FS, smoothing 3 deg: n=430, median_a=0.827, median_b=-0.01141, z=-17.91, p=1.039e-71
HD, smoothing 3 deg: n=1512, median_a=0.9721, median_b=-0.01075, z=-33.68, p=1.116e-248
FS, smoothing 2 deg: n=430, median_a=0.8, median_b=-0.004351, z=-17.92, p=8.43e-72
HD, smoothing 2 deg: n=1512, median_a=0.9658, median_b=-0.01389, z=-33.68, p=1.116e-248
paper: FS n = 427, z = 17.8; HD n = 1602, z = 34.7
cells with an undefined correlation (a flat half): 0
median per session of FS half-half r: {'min': 0.57, '50%': 0.83, 'max': 0.94}
/home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:376: RuntimeWarning: invalid value encountered in divide return (a * b).sum(0) / np.sqrt((a * a).sum(0) * (b * b).sum(0))
Result: reproduced. Correlations between the two halves of the square epoch are above the time-reversed controls for FS cells (median r 0.83 against -0.011, |z| = 17.91, p = 1e-71, n = 430) and HD cells (0.97 against -0.011, |z| = 33.7, p = 1e-248, n = 1,512), and the same with 2 degree smoothing (|z| = 17.92 and 33.7).
fig2c-geometry-independent¶
The tuning of PoSub-FS cells is preserved between the square and triangular arenas, as it is for PoSub-HD cells.
Sessions with a triangle epoch. FS cells above 10 Hz and HD cells above 0.5 Hz in both arenas (criterion; code Figure_Stability). Each session's triangle curves are shifted back by the circular mean of its HD cells' rotations (lag of maximum square-triangle cross-correlation; Methods). The control shifts the time-reversed triangle curves by the same amount and correlates them with the time-reversed square curves.
tri_rows = []
for sd in (3, 2):
rs, rc = [], []
for si, s in enumerate(S):
if f"tri{sd}" not in TC[si]:
continue
u = U[U.sess == si]
hdm = (u.hd & (u.rate > 0.5) & (u.rate_tri > 0.5)).values
fsm = (u.fs & (u.rate > 10) & (u.rate_tri > 10)).values
sq, tr = TC[si][f"sq{sd}"], TC[si][f"tri{sd}"]
sqr, trr = TC[si][f"sq{sd}_rev"], TC[si][f"tri{sd}_rev"]
lag = posub.circ_xcorr_lag(sq[:, hdm], tr[:, hdm])
rot = int(np.round(np.rad2deg(stats.circmean(np.deg2rad(lag))))) % 360
tr_s, trr_s = np.roll(tr, -rot, axis=0), np.roll(trr, -rot, axis=0)
for cls, m in (("HD", hdm), ("FS", fsm)):
r = posub.zero_lag_r(sq[:, m], tr_s[:, m]); c = posub.zero_lag_r(sqr[:, m], trr_s[:, m])
tri_rows += [dict(sd=sd, cls=cls, subject=s["subject"], rot=rot, r=a, r_rev=b) for a, b in zip(r, c)]
T2 = pd.DataFrame(tri_rows)
print("sessions with a triangle:", T2[T2.sd == 3].subject.nunique(), "; rotation (deg) per session:",
T2[(T2.sd == 3)].groupby("subject").rot.first().to_dict())
for sd in (3, 2):
for cls in ("FS", "HD"):
d = T2[(T2.sd == sd) & (T2.cls == cls)]
show(signrank(d.r, d.r_rev), f"{cls}, smoothing {sd} deg:")
print("paper: FS n = 264, z = 14.0; HD n = 1013, z = 27.6")
print(T2[T2.sd == 3].groupby(["cls", "subject"]).r.median().unstack(0).round(2).T.to_string())
sessions with a triangle: 20 ; rotation (deg) per session: {'A3701': 35, 'A3702': 5, 'A3703': 334, 'A3705': 10, 'A3706': 168, 'A3707': 38, 'A3708': 16, 'A3710': 4, 'A3711': 24, 'A3712': 5, 'A3713': 287, 'A3716': 19, 'A3717': 7, 'A3723': 353, 'A3728': 27, 'A3730': 217, 'A5505': 8, 'A5506': 355, 'A5507': 31, 'A5508': 177}
FS, smoothing 3 deg: n=273, median_a=0.7961, median_b=-0.01283, z=-14.18, p=1.264e-45
HD, smoothing 3 deg: n=1118, median_a=0.9565, median_b=-0.009054, z=-28.94, p=3.563e-184
FS, smoothing 2 deg: n=273, median_a=0.7743, median_b=-0.007338, z=-14.19, p=1.016e-45
HD, smoothing 2 deg: n=1118, median_a=0.9519, median_b=-0.009902, z=-28.94, p=3.332e-184
paper: FS n = 264, z = 14.0; HD n = 1013, z = 27.6
subject A3701 A3702 A3703 A3705 A3706 A3707 A3708 A3710 A3711 A3712 A3713 A3716 A3717 A3723 A3728 A3730 A5505 A5506 A5507 A5508
cls
FS 0.90 0.83 0.87 0.82 0.71 0.57 0.39 0.86 0.84 0.81 0.76 0.83 0.71 0.61 0.79 0.89 0.88 0.82 0.75 0.79
HD 0.95 0.96 0.96 0.98 0.97 0.77 0.96 0.98 0.95 0.98 0.97 0.98 0.95 0.93 0.88 0.96 0.92 0.93 0.96 0.93
Result: reproduced. In the 20 sessions with a triangle, square-triangle correlations after the session's rotation correction are above their time-reversed controls for FS cells (median r 0.80 against -0.01, |z| = 14.2, p = 1e-45, n = 273) and HD cells (0.96 against -0.01, |z| = 28.9, p = 4e-184, n = 1,118); with the code's 2 degree smoothing the same holds (|z| = 14.2 and 28.9).
fig2b-decoding-from-fs¶
The animal's HD can be decoded accurately from PoSub-FS cell spiking with a Bayesian decoder built on their tuning curves, and decoding error falls as more FS cells are used.
As GetBayesianDecoding.m (criterion): the square epoch from first to last tracked sample is split into two halves of equal time; tuning curves (360 bins, 3 degree smoothing) from the first half at speed above 2 cm/s; spike counts in 100 ms bins smoothed with a Gaussian of 2 bins s.d.; decoding over the second half, all times; error the absolute angular difference from the real HD sample nearest each bin centre, median per session. Chance: the same with tuning curves built from the HD time-reversed over the epoch, its error also taken from the real HD (criterion).
from scipy.ndimage import gaussian_filter1d
dec = []
for si, s in enumerate(S):
u = U[U.sess == si]
t, ang = s["hd_t"], s["hd"]
a, b = posub.tracked_epoch(s, "wake_square")
in_ep = (t >= a) & (t <= b)
ang_rev = ang.copy(); ang_rev[in_ep] = ang[in_ep][::-1]
mid = (a + b) / 2
m1 = posub.good_samples(s, "wake_square") & (t < mid)
edges = np.arange(a, b, 0.1); cen = edges[:-1] + 0.05
sec = cen >= mid
j = np.clip(np.searchsorted(t, cen[sec]), 1, len(t) - 1); j = np.where(np.abs(t[j - 1] - cen[sec]) < np.abs(t[j] - cen[sec]), j - 1, j)
row = dict(subject=s["subject"], n_fs=int(u.fs.sum()), n_hd=int(u.hd.sum()))
for cls in ("fs", "hd"):
idx = np.where(u[cls].values)[0]
if len(idx) == 0:
continue
sp = [s["spikes"][i] for i in idx]
Q = gaussian_filter1d(posub.bin_counts(sp, edges), 2, axis=0)
for kind, aa in (("real", ang), ("rev", ang_rev)):
tc, _ = posub.tuning_curves(sp, t, aa, m1 & np.isfinite(aa))
th = np.deg2rad(posub.bayes_decode(tc, Q[sec], 0.1))
err = np.abs(posub.angdiff(th, ang[j])) # both decoders scored against the real HD
ok = np.isfinite(ang[j])
row[f"{cls}_{kind}"] = np.rad2deg(np.median(err[ok]))
dec.append(row)
dec = pd.DataFrame(dec)
print(dec.round(1).to_string())
print(f"mean over {len(dec)} sessions of median error (deg): FS {dec.fs_real.mean():.1f}, FS chance {dec.fs_rev.mean():.1f}, "
f"HD {dec.hd_real.mean():.1f}, HD chance {dec.hd_rev.mean():.1f}")
show(signrank(dec.fs_real, dec.fs_rev), "FS error vs chance, per session:")
r, p = stats.pearsonr(dec.n_fs, dec.fs_real)
print(f"Pearson r between number of FS cells and FS median error = {r:.3f}, p = {p:.3g}, n = {len(dec)}")
rs, ps = stats.spearmanr(dec.n_fs, dec.fs_real)
print(f"(Spearman, shown only: rho = {rs:.3f}, p = {ps:.3g})")
fig, ax = plt.subplots(figsize=(4, 3)); ax.scatter(dec.n_fs, dec.fs_real); ax.set_xlabel("FS cells"); ax.set_ylabel("median error (deg)"); plt.show()
print("plotted ranges: cells", dec.n_fs.min(), dec.n_fs.max(), "error", dec.fs_real.min().round(1), dec.fs_real.max().round(1))
subject n_fs n_hd fs_real fs_rev hd_real hd_rev 0 A1808 11 30 52.8 97.0 9.8 80.8 1 A1813 17 31 67.7 67.5 10.7 78.6 2 A1815 16 33 23.7 105.5 9.2 123.5 3 A1821 16 37 29.4 93.5 10.0 91.1 4 A1824 17 39 18.8 84.4 9.2 72.9 5 A3701 12 66 21.3 106.3 11.6 86.9 6 A3702 16 61 19.4 94.1 8.4 113.4 7 A3703 20 65 26.0 73.5 8.1 96.7 8 A3705 25 66 16.7 105.3 9.2 82.7 9 A3706 21 86 26.7 84.3 7.5 84.6 10 A3707 32 117 35.2 90.4 10.3 78.5 11 A3708 12 43 55.3 89.5 8.2 109.1 12 A3709 8 22 36.4 98.6 11.3 112.9 13 A3710 13 72 32.3 99.4 8.6 119.3 14 A3711 21 89 33.9 109.0 10.6 95.7 15 A3712 23 72 21.1 89.8 9.2 72.3 16 A3713 13 25 31.8 113.1 11.1 82.7 17 A3716 17 91 30.4 90.5 9.6 97.5 18 A3717 18 57 23.7 101.3 10.4 130.2 19 A3723 5 37 45.3 112.0 9.6 84.1 20 A3728 5 27 74.0 76.9 14.4 105.2 21 A3730 9 46 24.4 87.1 9.0 98.3 22 A5505 10 30 22.1 117.1 8.9 84.6 23 A5506 11 24 54.8 98.0 12.5 72.0 24 A5507 10 42 32.3 67.1 8.6 50.3 25 A5508 6 33 34.7 85.2 8.8 96.5 26 A5801 9 26 47.9 82.9 11.7 79.8 27 A6211 11 35 48.0 90.9 11.5 79.1 28 A6215 7 29 30.9 79.0 8.8 112.0 29 A6216 14 44 32.5 94.5 13.5 73.7 30 A9701 5 37 66.9 92.7 7.9 127.8 mean over 31 sessions of median error (deg): FS 36.0, FS chance 92.8, HD 9.9, HD chance 92.7 FS error vs chance, per session: n=31, median_a=32.32, median_b=92.75, z=-4.84, p=1.296e-06 Pearson r between number of FS cells and FS median error = -0.449, p = 0.0113, n = 31 (Spearman, shown only: rho = -0.488, p = 0.00534)
plotted ranges: cells 5 32 error 16.7 74.0
Result: reproduced. Over 31 sessions the mean of the FS decoder's median error is 36.0 degrees against a chance level of 92.8 (both scored against the real HD; signed-rank over sessions p = 1.3e-6), and the error falls with the number of FS cells (Pearson r = -0.449, p = 0.011, n = 31). The HD-cell decoder's mean median error is 9.9 degrees.
fig2f-cue-rotation (figshare, data elsewhere)¶
In a cue rotation task PoSub-FS cell tuning curves follow the rotated distal landmark together with PoSub-HD cells.
The six cue rotation sessions are figshare Dataset_3.tar.gz (claims file), streamed here. Cell classes are the files' own CellTypes (hd, fs). Methods (Analysis of HD cell realignment after cue rotation): HD cell spikes in 50 ms bins smoothed with 100 ms s.d.; Bayesian decoding with tuning curves (1 degree bins, 3 degree smoothing) from the baseline period with a stable cue (LED1 first interval, as GetBayesianDecoding_CueRotation.m); decoder error the angular difference between real and decoded HD. Rotations at the starts of the 16 cue epochs (CueEpochs.mat). Exclusion (criterion): (1) the circular mean decoder error over the 200 s after the rotation differs by less than 45 degrees from that over the 200 s before it (the change, as the code computes it, since the cue returns to its first wall every second rotation); (2) mean speed below 2 cm/s in the preceding epoch (positions from Tracking_data.csv, metres). Main reading excludes a rotation failing either; the second reading only one failing both. A sigmoid is fitted to the decoder error (re-centred on the preceding epoch's mean error) from 50 s before to 200 s after each rotation; realignment ends where the normalized sigmoid reaches 0.99. For each cell, a tuning curve (1 degree, 3 degree smoothing, no speed threshold) over the rest of each epoch; its realignment is the lag of maximum cross-correlation with the curve of the preceding epoch. Per FS cell, the mean absolute difference between its realignment and the circular mean realignment of the session's HD cells, over the kept rotations; control, the same with the FS curves from the time-reversed HD; Wilcoxon signed-rank.
from scipy.optimize import curve_fit
CU = posub.load_adn(file_id=posub.FIGSHARE_CUE_FILE)
sigm = lambda x, a, b, c, d: a + (b - a) / (1 + 10 ** ((c - x) * d))
cue_rows, rot_rows = [], []
for c in CU:
t, ang = c["hd_t"], c["hd"]
hdc = np.where(c["file_types"]["hd"])[0]; fsc = np.where(c["file_types"]["fs"])[0]
b0, b1 = c["led1"][0]
mb = posub.in_intervals(t, [(b0, b1)]) & np.isfinite(ang)
tc_b, _ = posub.tuning_curves([c["spikes"][i] for i in hdc], t, ang, mb)
wa, wb = c["eps"]["wakeEp"][0]
edges = np.arange(wa, wb, 0.05); cen = edges[:-1] + 0.025
Q = gaussian_filter1d(posub.bin_counts([c["spikes"][i] for i in hdc], edges), 2, axis=0)
dec = np.deg2rad(posub.bayes_decode(tc_b, Q, 0.05))
j = np.clip(np.searchsorted(t, cen), 1, len(t) - 1); j = np.where(np.abs(t[j - 1] - cen) < np.abs(t[j] - cen), j - 1, j)
real = np.where(np.abs(t[j] - cen) < 0.02, ang[j], np.nan)
err = posub.angdiff(real, dec)
tr = c["track"]; _, uq = np.unique(tr[:, 0], return_index=True); tr = tr[uq] # repeated timestamps dropped
tt = tr[:, 0]; xy = tr[:, 1:3] * 100
sp_ = np.r_[0, np.hypot(*np.diff(gaussian_filter1d(xy, 10, axis=0), axis=0).T) / np.diff(tt)]
ce = c["cueEp"]
starts = ce[:, 0]
cmean = lambda x: np.angle(np.nanmean(np.exp(1j * x[np.isfinite(x)])))
for k, r0 in enumerate(starts):
pre = (cen >= r0 - 200) & (cen < r0); post = (cen >= r0) & (cen < r0 + 200)
change = abs(np.rad2deg(posub.angdiff(cmean(err[post]), cmean(err[pre]))))
prev_speed = sp_[(tt >= r0 - 200) & (tt < r0)].mean()
w = (cen >= r0 - 50) & (cen < r0 + 200) & np.isfinite(err)
x = cen[w]; y = np.rad2deg(posub.angdiff(err[w], cmean(err[pre])))
try:
p, _ = curve_fit(sigm, x, y, p0=[np.quantile(y, 0.05), np.quantile(y, 0.95), r0 + 10, 0.1], maxfev=20000)
yv = sigm(x, *p); yn = (yv - yv.min()) / (yv.max() - yv.min())
if yv[0] > yv[-1]:
yn = 1 - yn
end = x[np.argmin(np.abs(yn - 0.99))]
except Exception:
end = np.nan
# the code's rule (GetSigmoids_CueRotation.m): windows of the last 50 s of each epoch
w_pre = (cen >= r0 - 50) & (cen < r0); w_post = (cen >= r0 + 150) & (cen < r0 + 200)
code_change = abs(np.rad2deg(posub.angdiff(cmean(err[w_post]), cmean(err[w_pre]))))
e_ = err[w_pre][np.isfinite(err[w_pre])]
code_var = np.rad2deg(1 - np.abs(np.mean(np.exp(1j * e_))))
mh = (t >= r0 - 50) & (t < r0) & np.isfinite(ang)
ahv = np.nanmean(np.abs(np.rad2deg(np.diff(np.unwrap(ang[mh]))) / np.diff(t[mh])))
rot_rows.append(dict(session=c["session"], k=k, change=change, prev_speed=prev_speed, end_s=end - r0 if np.isfinite(end) else np.nan,
code_change=code_change, code_var=code_var, ahv=ahv))
R = pd.DataFrame([r for r in rot_rows if r["session"] == c["session"]])
# tuning curves over the rest of each epoch (from the end of realignment), and the epoch before the first rotation
segs = [(starts[0] - 200, starts[0])] + [((r0 + (e if np.isfinite(e) else 0)), ce[k, 1]) for k, (r0, e) in enumerate(zip(starts, R.end_s))]
cells = np.r_[hdc, fsc]
H = []; Hr = []
for a_, b_ in segs:
m = posub.in_intervals(t, [(a_, b_)]) & np.isfinite(ang)
H.append(posub.tuning_curves([c["spikes"][i] for i in cells], t, ang, m)[0])
Hr.append(posub.tuning_curves([c["spikes"][i] for i in cells], t, ang, m, reverse=True)[0])
nh = len(hdc)
for k in range(len(starts)):
lag = posub.circ_xcorr_lag(H[k], H[k + 1]); lag_r = posub.circ_xcorr_lag(Hr[k], Hr[k + 1])
hd_mean = stats.circmean(np.deg2rad(lag[:nh][np.isfinite(lag[:nh])]))
for jj, i in enumerate(fsc):
d = abs(np.rad2deg(posub.angdiff(np.deg2rad(lag[nh + jj]), hd_mean)))
dr = abs(np.rad2deg(posub.angdiff(np.deg2rad(lag_r[nh + jj]), hd_mean)))
cue_rows.append(dict(session=c["session"], cell=i, k=k, diff=d, diff_rev=dr,
hd_rot=np.rad2deg(hd_mean) % 360, change=R.change[k], prev_speed=R.prev_speed[k],
code_change=R.code_change[k], code_var=R.code_var[k], ahv=R.ahv[k]))
ROT = pd.DataFrame(rot_rows); CR = pd.DataFrame(cue_rows)
print(ROT.round(1).to_string())
print("rotations failing: change < 45:", int((ROT.change < 45).sum()), "; speed < 2 cm/s:", int((ROT.prev_speed < 2).sum()), "; realignment end (s) median", ROT.end_s.median().round(1), "range", ROT.end_s.min().round(1), ROT.end_s.max().round(1))
for name, excl in (("main: exclude a rotation failing either condition", lambda d: (d.change < 45) | (d.prev_speed < 2)),
("second: exclude only a rotation failing both", lambda d: (d.change < 45) & (d.prev_speed < 2)),
("change condition alone", lambda d: d.change < 45),
("code variant: change < 45, variance > 10 or AHV < 10", lambda d: (d.code_change < 45) | (d.code_var > 10) | (d.ahv < 10))):
kept = CR[~excl(CR)]
per = kept.groupby(["session", "cell"])[["diff", "diff_rev"]].mean()
print(f"{name}: rotations kept {int((~excl(ROT)).sum())} of {len(ROT)}; FS cells {len(per)}")
show(signrank(per["diff"], per["diff_rev"]), " FS minus HD realignment |difference| (deg), real vs control:")
print("paper: n = 99 FS cells, z = 8.25")
print("HD rotation (deg) at kept rotations, main reading:", CR[~((CR.change < 45) | (CR.prev_speed < 2))].groupby(["session", "k"]).hd_rot.first().round(0).values)
/tmp/ipykernel_2279/3067355471.py:3: RuntimeWarning: overflow encountered in power sigm = lambda x, a, b, c, d: a + (b - a) / (1 + 10 ** ((c - x) * d))
/home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:383: RuntimeWarning: invalid value encountered in divide z2 = (h2 - h2.mean(0)) / h2.std(0) /home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:382: RuntimeWarning: invalid value encountered in divide z1 = (h1 - h1.mean(0)) / h1.std(0)
session k change prev_speed end_s code_change code_var ahv 0 A3702-191128 0 88.8 3.6 17.8 90.3 1.2 98.0 1 A3702-191128 1 88.7 4.6 6.0 90.4 1.1 65.8 2 A3702-191128 2 82.5 5.0 10.2 86.7 0.6 67.4 3 A3702-191128 3 79.9 3.7 15.1 84.4 1.6 42.5 4 A3702-191128 4 82.8 4.3 33.0 89.4 0.3 65.8 5 A3702-191128 5 88.0 4.2 4.0 92.3 1.6 93.7 6 A3702-191128 6 81.6 4.3 49.6 92.9 1.0 94.1 7 A3702-191128 7 69.5 5.3 33.7 89.1 1.7 82.8 8 A3702-191128 8 75.4 6.3 9.8 81.4 0.9 82.4 9 A3702-191128 9 81.5 4.5 12.8 84.0 0.9 70.9 10 A3702-191128 10 81.7 4.1 13.1 84.9 1.0 99.0 11 A3702-191128 11 87.3 5.2 12.1 79.2 1.2 77.4 12 A3702-191128 12 88.6 4.2 23.9 82.4 1.0 62.7 13 A3702-191128 13 84.8 5.7 51.3 100.2 1.7 84.7 14 A3702-191128 14 79.8 4.1 14.8 95.9 1.4 83.9 15 A3702-191128 15 77.4 4.1 7.4 77.2 1.1 67.8 16 A3706-200315 0 91.1 2.3 18.7 87.8 0.7 59.8 17 A3706-200315 1 86.5 3.2 5.2 93.3 0.6 65.5 18 A3706-200315 2 85.0 2.6 13.1 97.3 0.6 39.2 19 A3706-200315 3 64.9 2.9 65.5 89.1 1.0 45.0 20 A3706-200315 4 75.6 2.3 12.6 91.5 0.5 53.9 21 A3706-200315 5 86.6 2.4 44.9 93.3 0.8 46.9 22 A3706-200315 6 86.8 2.0 27.1 109.2 0.8 34.5 23 A3706-200315 7 91.5 2.2 20.8 114.0 0.6 42.5 24 A3706-200315 8 85.1 1.9 52.6 96.8 1.9 51.7 25 A3706-200315 9 71.2 2.1 51.7 87.7 1.3 76.7 26 A3706-200315 10 66.1 2.3 81.8 90.4 0.7 50.8 27 A3706-200315 11 77.3 1.9 18.0 80.9 2.7 51.7 28 A3706-200315 12 81.8 2.3 50.6 81.2 1.2 78.7 29 A3706-200315 13 32.5 2.3 158.6 41.2 0.9 42.0 30 A3706-200315 14 31.3 2.6 0.2 40.5 2.0 44.7 31 A3706-200315 15 28.1 1.8 5.3 35.4 1.1 32.7 32 A3707-200318 0 89.1 3.2 3.9 97.2 4.1 54.6 33 A3707-200318 1 90.7 4.9 21.0 107.8 1.2 82.6 34 A3707-200318 2 89.8 3.3 33.2 98.8 2.4 51.0 35 A3707-200318 3 82.1 3.8 8.5 81.2 1.2 80.8 36 A3707-200318 4 88.2 4.0 6.0 87.9 1.1 40.6 37 A3707-200318 5 91.5 4.2 2.7 95.2 1.1 77.4 38 A3707-200318 6 86.5 4.2 5.7 92.6 2.0 67.2 39 A3707-200318 7 87.6 3.5 4.5 90.5 1.9 51.6 40 A3707-200318 8 86.1 2.4 14.6 96.2 2.9 42.0 41 A3707-200318 9 87.0 3.7 3.0 96.3 2.4 67.2 42 A3707-200318 10 90.6 3.9 2.3 97.3 0.6 74.4 43 A3707-200318 11 88.3 4.5 14.0 91.6 0.8 57.8 44 A3707-200318 12 87.5 4.5 9.9 79.6 1.3 88.0 45 A3707-200318 13 87.4 5.4 7.8 83.9 2.3 93.6 46 A3707-200318 14 83.8 5.3 8.0 86.9 0.7 77.2 47 A3707-200318 15 79.7 5.7 4.2 82.5 0.8 94.5 48 A3716-201015 0 81.0 4.7 16.1 85.6 2.0 98.8 49 A3716-201015 1 79.5 4.6 6.6 87.5 0.7 74.1 50 A3716-201015 2 83.1 4.7 3.2 87.2 3.6 134.8 51 A3716-201015 3 84.9 3.5 21.3 88.9 2.0 58.8 52 A3716-201015 4 84.3 4.0 11.2 92.8 0.8 107.3 53 A3716-201015 5 87.8 5.5 12.9 93.4 3.6 115.7 54 A3716-201015 6 79.7 4.0 33.7 76.5 0.9 97.5 55 A3716-201015 7 76.2 6.3 24.4 72.5 4.3 128.0 56 A3716-201015 8 85.0 3.8 24.8 84.1 1.2 90.1 57 A3716-201015 9 81.9 2.4 9.9 85.7 3.9 127.6 58 A3716-201015 10 75.1 3.0 24.9 78.7 4.4 68.6 59 A3716-201015 11 84.7 5.2 7.5 85.1 1.9 114.7 60 A3716-201015 12 87.1 4.3 32.6 98.6 1.9 74.6 61 A3716-201015 13 89.2 2.4 5.0 97.1 7.7 32.9 62 A3716-201015 14 85.8 2.4 35.3 85.7 3.2 137.8 63 A3716-201015 15 77.1 3.8 12.7 82.4 2.7 100.0 64 A3717-201022 0 92.8 2.8 22.2 87.9 1.7 60.3 65 A3717-201022 1 86.2 5.8 11.4 95.4 1.6 53.6 66 A3717-201022 2 86.8 3.8 9.6 97.4 2.3 113.6 67 A3717-201022 3 89.5 3.8 15.4 87.8 3.4 104.0 68 A3717-201022 4 82.9 2.9 30.9 89.3 1.3 64.0 69 A3717-201022 5 84.5 2.8 13.9 93.7 4.9 22.7 70 A3717-201022 6 86.8 3.6 25.3 84.4 5.6 78.1 71 A3717-201022 7 84.4 4.3 12.7 82.0 1.5 35.5 72 A3717-201022 8 86.7 2.9 7.9 88.6 3.4 101.1 73 A3717-201022 9 84.1 2.5 13.7 94.1 1.2 63.2 74 A3717-201022 10 83.9 3.2 8.1 96.7 5.0 28.9 75 A3717-201022 11 82.4 3.7 26.8 92.7 1.4 82.1 76 A3717-201022 12 87.9 3.0 19.5 99.3 7.4 65.5 77 A3717-201022 13 97.6 3.9 24.1 110.8 3.6 73.1 78 A3717-201022 14 94.7 1.0 5.7 97.9 14.3 10.8 79 A3717-201022 15 82.9 2.7 14.7 78.8 1.9 83.1 80 A3730-210324 0 89.5 4.8 13.8 94.6 0.4 63.0 81 A3730-210324 1 91.1 5.1 2.4 91.8 1.8 73.8 82 A3730-210324 2 88.0 3.4 6.1 82.3 2.5 47.9 83 A3730-210324 3 86.5 3.2 17.1 84.2 0.9 79.1 84 A3730-210324 4 81.9 3.5 12.9 86.6 0.7 53.6 85 A3730-210324 5 78.7 4.0 16.6 88.9 0.9 57.3 86 A3730-210324 6 81.0 3.1 24.5 90.9 0.7 70.5 87 A3730-210324 7 84.6 3.5 25.3 92.5 1.5 50.5 88 A3730-210324 8 79.5 3.0 12.1 87.1 0.9 40.7 89 A3730-210324 9 78.7 3.0 20.8 87.3 1.6 86.5 90 A3730-210324 10 79.9 3.8 54.2 89.5 1.1 53.1 91 A3730-210324 11 79.6 2.8 2.9 85.5 2.9 45.2 92 A3730-210324 12 84.1 3.5 19.3 94.4 0.6 51.1 93 A3730-210324 13 79.8 2.7 1.9 78.4 2.2 27.0 94 A3730-210324 14 75.3 1.6 55.7 66.2 22.3 6.4 95 A3730-210324 15 87.0 2.3 0.6 79.4 3.7 37.7 rotations failing: change < 45: 3 ; speed < 2 cm/s: 5 ; realignment end (s) median 13.8 range 0.2 158.6 main: exclude a rotation failing either condition: rotations kept 89 of 96; FS cells 99 FS minus HD realignment |difference| (deg), real vs control: n=99, median_a=36.17, median_b=91.79, z=-8.569, p=1.047e-17 second: exclude only a rotation failing both: rotations kept 95 of 96; FS cells 99 FS minus HD realignment |difference| (deg), real vs control: n=99, median_a=36.55, median_b=90.05, z=-8.59, p=8.727e-18 change condition alone: rotations kept 93 of 96; FS cells 99 FS minus HD realignment |difference| (deg), real vs control: n=99, median_a=35.97, median_b=90.43, z=-8.583, p=9.273e-18 code variant: change < 45, variance > 10 or AHV < 10: rotations kept 91 of 96; FS cells 99 FS minus HD realignment |difference| (deg), real vs control: n=99, median_a=35.97, median_b=91.6, z=-8.586, p=8.996e-18 paper: n = 99 FS cells, z = 8.25 HD rotation (deg) at kept rotations, main reading: [ 91. 269. 85. 275. 91. 267. 92. 272. 87. 272. 88. 269. 92. 264. 92. 277. 93. 270. 88. 274. 93. 267. 93. 267. 273. 76. 91. 271. 94. 267. 84. 273. 90. 273. 88. 273. 89. 270. 88. 270. 89. 273. 83. 86. 274. 85. 272. 90. 268. 90. 274. 87. 277. 79. 272. 88. 274. 87. 277. 263. 89. 269. 94. 272. 85. 269. 91. 273. 86. 272. 87. 270. 97. 83. 270. 94. 271. 89. 271. 85. 273. 89. 273. 87. 276. 81. 271. 88. 89.]
Result: reproduced. In the six figshare cue rotation sessions, the realignment of the 99 FS cells is closer to that of the HD cells than in the time-reversed control (median absolute difference 36.2 against 91.8 degrees, |z| = 8.57, p = 1e-17; paper z = 8.25), with 89 of 96 rotations kept under the main exclusion reading, and the same under the other readings and the code's rule (|z| 8.58 to 8.59). The HD cells rotate by about 90 and 270 degrees at the kept rotations.
Fourier spectra¶
Each 360-bin tuning curve (3 degree smoothing) is projected on harmonics 1 to 180 (code, Figure_Fourier). Power is read two ways (claims file): the squared modulus of the complex coefficient (code) and the root mean square of the sine and cosine projections (Methods' words). The first ten powers are normalized to their sum.
H = per_unit("sq3"); H_rev = per_unit("sq3_rev")
FT, FT_rev = posub.fourier(H), posub.fourier(H_rev)
P = {rd: posub.power(FT, rd) for rd in ("sq", "rms")}
P_rev = {rd: posub.power(FT_rev, rd) for rd in ("sq", "rms")}
N10 = {rd: posub.norm10(P[rd]) for rd in P}
N10_rev = {rd: posub.norm10(P_rev[rd]) for rd in P}
fsm, hdm = U.fs.values, U.hd.values
# Parseval check (Extended Data Fig. 3e): 2 x sum of squared moduli = variance of the curve
pv = 2 * P["sq"].sum(0) / H.var(0)
print("Parseval, 2 x sum |c_k|^2 over the curve's variance: median", np.nanmedian(pv[fsm | hdm]).round(4),
"range", np.nanmin(pv[fsm | hdm]).round(4), np.nanmax(pv[fsm | hdm]).round(4))
Parseval, 2 x sum |c_k|^2 over the curve's variance: median 1.0 range 1.0 1.0
fig3-ten-components-power¶
Fourier power of tuning curves decays quickly with frequency, and the first ten components hold on average 98% of the total power.
Expected 0.98, tolerance 0.02. Readings: population (a) FS and HD pooled, (b) each alone; power (a) squared modulus, (b) root mean square. Reproduced when one population reading under each power reading is within tolerance.
rows = []
for rd in ("sq", "rms"):
frac = P[rd][:10].sum(0) / P[rd].sum(0)
for pop, m in (("FS + HD", fsm | hdm), ("FS", fsm), ("HD", hdm)):
rows.append(dict(power=rd, population=pop, n=int(m.sum()), mean_fraction=frac[m].mean().round(4),
within=abs(frac[m].mean() - 0.98) <= 0.02))
print(pd.DataFrame(rows).to_string())
mp = P["sq"][:, fsm | hdm].mean(1)
print("mean squared-modulus power, harmonics 1 to 12, relative to harmonic 1:", (mp[:12] / mp[0]).round(4))
power population n mean_fraction within 0 sq FS + HD 1942 0.9888 True 1 sq FS 430 0.9691 True 2 sq HD 1512 0.9944 True 3 rms FS + HD 1942 0.8716 False 4 rms FS 430 0.7783 False 5 rms HD 1512 0.8981 False mean squared-modulus power, harmonics 1 to 12, relative to harmonic 1: [1. 0.6845 0.4009 0.2227 0.1005 0.0467 0.0216 0.0102 0.0051 0.0031 0.0019 0.0017]
Result: partially reproduced. With power as the squared modulus (the code, and the only reading under which the Parseval check holds, ratio 1.0000), the first ten harmonics hold on average 0.989 of the power of FS and HD cells pooled, 0.969 for FS and 0.994 for HD cells, all within 0.98 +- 0.02, and the mean power falls to 0.10 of the first harmonic's by the fifth. With power as the root mean square of the projections (the Methods' words) the fractions are 0.872, 0.778 and 0.898, none within tolerance, and the criterion asks for one population within tolerance under each power reading.
fig3g-first-three-components¶
Compared with time-reversed controls, the PoSub-FS spectrum has higher power for the first three Fourier components only.
Repeated-measures ANOVA, condition (real, time-reversed) by component (1 to 10), within cells; then a Wilcoxon signed-rank test per component with p multiplied by 10 (criterion; code).
from statsmodels.stats.anova import AnovaRM
for rd in ("sq", "rms"):
a, b = N10[rd][:, fsm], N10_rev[rd][:, fsm]
n = a.shape[1]
long = pd.DataFrame(dict(cell=np.tile(np.repeat(np.arange(n), 10), 2), comp=np.tile(np.tile(np.arange(1, 11), n), 2),
cond=np.repeat(["real", "rev"], 10 * n), y=np.r_[a.T.ravel(), b.T.ravel()]))
res = AnovaRM(long, "y", "cell", within=["cond", "comp"]).fit().anova_table
print(f"power {rd}: interaction F({int(res.loc['cond:comp', 'Num DF'])},{int(res.loc['cond:comp', 'Den DF'])}) = "
f"{res.loc['cond:comp', 'F Value']:.1f}, p = {res.loc['cond:comp', 'Pr > F']:.3g}")
out = []
for k in range(10):
d = signrank(a[k], b[k])
out.append(dict(component=k + 1, real_median=round(d["median_a"], 4), rev_median=round(d["median_b"], 4),
z=round(d["z"], 2), z_over_10=round(d["z"] / 10, 3), p_bonf=min(1, d["p"] * 10),
higher=d["median_a"] > d["median_b"]))
out = pd.DataFrame(out); print(out.to_string())
sig_hi = out[(out.p_bonf < 0.05) & out.higher].component.tolist()
sig_hi_mean = [k + 1 for k in range(10) if out.p_bonf[k] < 0.05 and a[k].mean() > b[k].mean()]
print(f"power {rd}: components significantly higher than control (medians): {sig_hi}; (means): {sig_hi_mean}")
print("paper: F(9,3834) = 64.1, p = 4e-110; components 1 to 3 higher (p = 8e-10, 3e-16, 4e-5)")
power sq: interaction F(9,3861) = 122.9, p = 1.02e-203 component real_median rev_median z z_over_10 p_bonf higher 0 1 0.3813 0.2013 -11.20 -1.120 4.003255e-28 True 1 2 0.2364 0.1190 -9.36 -0.936 7.838817e-20 True 2 3 0.1138 0.0926 -2.73 -0.273 6.407072e-02 True 3 4 0.0613 0.0761 -1.49 -0.149 1.000000e+00 False 4 5 0.0218 0.0586 -10.09 -1.009 5.872187e-23 False 5 6 0.0105 0.0448 -13.02 -1.302 9.547112e-38 False 6 7 0.0060 0.0472 -15.97 -1.597 2.033490e-56 False 7 8 0.0041 0.0394 -16.02 -1.602 9.163307e-57 False 8 9 0.0026 0.0395 -16.96 -1.696 1.540240e-63 False 9 10 0.0020 0.0388 -17.04 -1.704 3.941121e-64 False power sq: components significantly higher than control (medians): [1, 2]; (means): [1, 2]
power rms: interaction F(9,3861) = 227.2, p = 0 component real_median rev_median z z_over_10 p_bonf higher 0 1 0.2637 0.1582 -13.72 -1.372 7.252799e-42 True 1 2 0.2132 0.1252 -12.62 -1.262 1.569467e-35 True 2 3 0.1498 0.1101 -8.10 -0.810 5.416439e-15 True 3 4 0.1095 0.0984 -3.62 -0.362 2.993661e-03 True 4 5 0.0652 0.0891 -6.49 -0.649 8.628602e-10 False 5 6 0.0448 0.0784 -11.06 -1.106 1.925853e-27 False 6 7 0.0340 0.0776 -15.14 -1.514 8.933457e-51 False 7 8 0.0281 0.0717 -15.46 -1.546 6.019120e-53 False 8 9 0.0230 0.0726 -16.64 -1.664 3.788500e-61 False 9 10 0.0204 0.0728 -16.92 -1.692 3.287416e-63 False power rms: components significantly higher than control (medians): [1, 2, 3, 4]; (means): [1, 2, 3, 4] paper: F(9,3834) = 64.1, p = 4e-110; components 1 to 3 higher (p = 8e-10, 3e-16, 4e-5)
Result: partially reproduced. The condition by component interaction is significant under both power readings (squared modulus F(9,3861) = 122.9, p = 1e-203; root mean square F = 227.2; paper F(9,3834) = 64.1). After Bonferroni correction, with squared-modulus power components 1 and 2 are higher than control (p = 4e-28 and 8e-20) and component 3 is higher in median (0.114 against 0.093) but not significantly (p = 0.064); with root-mean-square power components 1 to 4 are higher (component 4, p = 0.003). The criterion's set, components 1, 2 and 3 and no other, holds under neither power reading; the FS spectrum is above control at the lowest components and below it from component 5 on under both.
fig3h-dominant-component¶
Of the 427 PoSub-FS cells, 231, 141 and 55 have their highest power in the first, second and third Fourier component, proportions similar to the relative mean power of those components.
Each FS cell goes to the largest of its first three powers (code, highComp). Readings of "similar": (a) the proportions are in the same order as the mean relative powers of components 1 to 3 renormalized over those three; (b) each proportion is within 10% of its renormalized mean power (read as a relative 10%; an absolute 0.10 is shown too).
high = np.argmax(N10["sq"][:3], axis=0) + 1
U["high"] = high
cnt = pd.Series(high[fsm]).value_counts().sort_index()
print("FS cells by dominant component (1, 2, 3):", cnt.to_dict(), "of", int(fsm.sum()), "; expected 231 +-23, 141 +-14, 55 +-6")
print("same assignment under the rms reading:", bool((np.argmax(N10["rms"][:3], 0)[fsm] == high[fsm] - 1).all()))
prop = cnt.values / cnt.sum()
for rd in ("sq", "rms"):
mpw = N10[rd][:3, fsm].mean(1); mpw = mpw / mpw.sum()
a_ok = list(np.argsort(-prop)) == list(np.argsort(-mpw))
rel = np.abs(prop - mpw) / mpw
print(f"power {rd}: proportions {prop.round(3)}, renormalized mean powers {mpw.round(3)}; (a) same order: {a_ok}; "
f"(b) relative differences {rel.round(3)}, all within 10%: {bool((rel <= 0.10).all())}; absolute differences {np.abs(prop - mpw).round(3)}")
print("per subject, FS cells by dominant component:")
print(U[U.fs].groupby("subject").high.value_counts().unstack(fill_value=0).T.to_string())
FS cells by dominant component (1, 2, 3): {1: 237, 2: 146, 3: 47} of 430 ; expected 231 +-23, 141 +-14, 55 +-6
same assignment under the rms reading: True
power sq: proportions [0.551 0.34 0.109], renormalized mean powers [0.486 0.326 0.187]; (a) same order: True; (b) relative differences [0.133 0.04 0.416], all within 10%: False; absolute differences [0.065 0.013 0.078]
power rms: proportions [0.551 0.34 0.109], renormalized mean powers [0.425 0.332 0.243]; (a) same order: True; (b) relative differences [0.297 0.023 0.551], all within 10%: False; absolute differences [0.126 0.008 0.134]
per subject, FS cells by dominant component:
subject A1808 A1813 A1815 A1821 A1824 A3701 A3702 A3703 A3705 A3706 A3707 A3708 A3709 A3710 A3711 A3712 A3713 A3716 A3717 A3723 A3728 A3730 A5505 A5506 A5507 A5508 A5801 A6211 A6215 A6216 A9701
high
1 6 9 8 7 7 5 8 12 18 12 19 6 2 8 12 12 7 8 8 5 4 4 7 7 5 2 6 4 6 9 4
2 5 4 7 4 8 6 4 7 7 6 11 3 5 4 6 9 2 8 8 0 1 4 2 4 4 4 3 4 0 5 1
3 0 4 1 5 2 1 4 1 0 3 2 3 1 1 3 2 4 1 2 0 0 1 1 0 1 0 0 3 1 0 0
Result: not reproduced. FS cells peak in component 1, 2 and 3 in 237, 146 and 47 cases; the first two are within tolerance of 231 and 141, the third is below 55 - 6. The criterion's reproduced and partial outcomes both need all counts within tolerance, so the claim is not reproduced. The proportions (0.551, 0.340, 0.109) are in the same order as the renormalized mean powers (0.486, 0.326, 0.187), reading (a), but differ by more than 10% for components 1 and 3 (13% and 42%), reading (b), under both power readings. The dandiset has 430 FS cells where the paper has 427.
ed3g-fs-spectra-heterogeneous¶
Fourier signatures are far more variable among PoSub-FS cells than among HD cells of the same region, so individual FS cells are narrowly tuned in Fourier space.
Distance of each cell's ten-component spectrum from its own population's mean: Euclidean (code, ftDist) and KL divergence D_KL(cell || mean) (legend, Methods direction), both power readings. FS against PoSub-HD, Mann-Whitney U, p multiplied by 3. The ANOVA and the ADN-HD comparisons are in the ADN section below.
DIST = {}
for rd in ("sq", "rms"):
for meas in ("euclid", "kl"):
d = np.full(len(U), np.nan)
for m in (fsm, hdm):
mu = N10[rd][:, m].mean(1, keepdims=True)
d[m] = np.sqrt(((N10[rd][:, m] - mu) ** 2).sum(0)) if meas == "euclid" else posub.kl(N10[rd][:, m], mu)
DIST[(rd, meas)] = d
r = ranksum(d[fsm], d[hdm]); r["p_x3"] = min(1, r["p"] * 3); r["z_over_3"] = r["z"] / 3
show(r, f"power {rd}, {meas}: FS vs PoSub-HD")
print("paper: legend z = 5.71 (a third of the test's z)")
power sq, euclid: FS vs PoSub-HD n_a=430, n_b=1512, median_a=0.311, median_b=0.1468, z=17.35, p=2.115e-67, p_x3=6.344e-67, z_over_3=5.782 power sq, kl: FS vs PoSub-HD n_a=430, n_b=1512, median_a=0.2589, median_b=0.08065, z=17.9, p=1.157e-71, p_x3=3.471e-71, z_over_3=5.967 power rms, euclid: FS vs PoSub-HD n_a=430, n_b=1512, median_a=0.168, median_b=0.09365, z=17.39, p=9.137e-68, p_x3=2.741e-67, z_over_3=5.798 power rms, kl: FS vs PoSub-HD n_a=430, n_b=1512, median_a=0.1105, median_b=0.04184, z=17.57, p=3.919e-69, p_x3=1.176e-68, z_over_3=5.858 paper: legend z = 5.71 (a third of the test's z)
ed3l-fs-phases-independent¶
Phases of the Fourier components are correlated with each other in individual HD cells but random relative to each other in PoSub-FS cells.
Circular correlation (circ_corrcc) across cells between twice the first phase (mod 2 pi) and the second phase, and three times the first and the third (criterion; code).
ph = np.angle(FT[:3])
for name, m in (("PoSub-HD", hdm), ("PoSub-FS", fsm)):
for k in (2, 3):
rho, p = posub.circ_corrcc(np.mod(k * ph[0, m], 2 * np.pi), np.mod(ph[k - 1, m], 2 * np.pi))
print(f"{name}: {k} x phase 1 against phase {k}: rho = {rho:.3f}, p = {p:.3g}, n = {int(m.sum())}")
# the FS result with code classes
exc_, fsc_ = classify(U, 0.35, True)
for k in (2, 3):
rho, p = posub.circ_corrcc(np.mod(k * ph[0, fsc_.values], 2 * np.pi), np.mod(ph[k - 1, fsc_.values], 2 * np.pi))
print(f"PoSub-FS, code classes: {k} x phase 1 against phase {k}: rho = {rho:.3f}, p = {p:.3g}")
PoSub-HD: 2 x phase 1 against phase 2: rho = 0.322, p = 0, n = 1512 PoSub-HD: 3 x phase 1 against phase 3: rho = 0.087, p = 9.79e-05, n = 1512 PoSub-FS: 2 x phase 1 against phase 2: rho = 0.014, p = 0.768, n = 430 PoSub-FS: 3 x phase 1 against phase 3: rho = 0.005, p = 0.913, n = 430 PoSub-FS, code classes: 2 x phase 1 against phase 2: rho = 0.006, p = 0.902 PoSub-FS, code classes: 3 x phase 1 against phase 3: rho = -0.016, p = 0.738
Result: reproduced. Across PoSub-HD cells twice the first phase is correlated with the second (rho = 0.32, p < 1e-15) and three times the first with the third (rho = 0.087, p = 1e-4); across FS cells neither is (rho = 0.014, p = 0.77; rho = 0.005, p = 0.91), nor with the code's classes (p = 0.90 and 0.74).
ed4a-autocorrelogram-symmetries¶
Tuning curve auto-correlograms of FS cells with the highest power in the second component peak at 180 degrees, and those with the highest power in the third at 120 and 240 degrees, above controls.
Cross-validated auto-correlograms from the two halves (3 degree smoothing), value at 180 degrees for twofold cells and at 120 degrees for threefold cells, against the same value from the time-reversed halves; Wilcoxon signed-rank (legend), Mann-Whitney U (code) as a variant. 240 degrees is shown and not judged.
AC = posub.xval_autocorr(per_unit("h1_3"), per_unit("h2_3"))
AC_rev = posub.xval_autocorr(per_unit("h1_3_rev"), per_unit("h2_3_rev"))
for k, lags in ((2, (180,)), (3, (120, 240))):
m = fsm & (high == k)
for lag in lags:
show(signrank(AC[lag, m], AC_rev[lag, m]), f"FS, dominant component {k}, {lag} deg, signed-rank:")
show(ranksum(AC[lag, m], AC_rev[lag, m]), f"FS, dominant component {k}, {lag} deg, rank-sum:")
m1 = fsm & (high == 1)
print("for comparison, onefold FS cells at 180 deg: median", np.median(AC[180, m1]).round(3), "control", np.median(AC_rev[180, m1]).round(3))
FS, dominant component 2, 180 deg, signed-rank: n=146, median_a=0.1934, median_b=0.01829, z=-7.3, p=2.885e-13 FS, dominant component 2, 180 deg, rank-sum: n_a=146, n_b=146, median_a=0.1934, median_b=0.01829, z=7.327, p=2.364e-13 FS, dominant component 3, 120 deg, signed-rank: n=47, median_a=0.1713, median_b=-0.002327, z=-3.683, p=0.0002309 FS, dominant component 3, 120 deg, rank-sum: n_a=47, n_b=47, median_a=0.1713, median_b=-0.002327, z=4.034, p=5.566e-05 FS, dominant component 3, 240 deg, signed-rank: n=47, median_a=0.182, median_b=-0.01674, z=-4.667, p=3.06e-06 FS, dominant component 3, 240 deg, rank-sum: n_a=47, n_b=47, median_a=0.182, median_b=-0.01674, z=4.866, p=1.16e-06 for comparison, onefold FS cells at 180 deg: median -0.4 control -0.022
/home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:362: RuntimeWarning: invalid value encountered in divide z1 = (h1 - h1.mean(0)) / h1.std(0, ddof=1) /home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:363: RuntimeWarning: invalid value encountered in divide z2 = (h2 - h2.mean(0)) / h2.std(0, ddof=1)
Result: reproduced. The cross-validated auto-correlogram at 180 degrees of the 146 twofold FS cells is above control (median 0.19 against 0.02, signed-rank p = 3e-13; rank-sum p = 2e-13), and at 120 degrees for the 47 threefold cells (0.17 against -0.002, p = 2e-4; rank-sum p = 6e-5). At 240 degrees, shown only, 0.18 against -0.02.
ADN recordings (figshare, claims with data elsewhere)¶
The ADN-HD cells of Fig. 3 are not in the dandiset. The claims file names their source, the
paper's figshare item doi:10.6084/m9.figshare.24921252 version 3, Dataset_5.tar.gz. It is
streamed and read in memory (helper load_adn). Tuning curves are computed from the spike
times, head direction and the file's speed as for PoSub; ADN-HD cells are units with rate
above 0.5 Hz and HD information at least 0.2 bits per spike (code, GetCellTypes_ADN), among
units the file's BrainArea marks as ADN (code, Figure_Fourier isAdn).
A = posub.load_adn()
print("ADN sessions:", len(A), "; units:", sum(len(a["spikes"]) for a in A), "; marked ADN:", sum(int(a["adn"].sum()) for a in A))
ADN = []
for a in A:
d = posub.adn_tuning(a)
d10 = posub.adn_tuning(a, speed_th=10.0)
hd = a["adn"] & (d["rate"] > 0.5) & (d["info"] >= 0.2)
d.update(hd=hd, h10=d10["h"], file_hd=a["file_hd"], session=a["session"])
ADN.append(d)
adn_hd = np.concatenate([d["hd"] for d in ADN]); adn_file = np.concatenate([d["file_hd"] for d in ADN])
print("ADN-HD cells recomputed:", int(adn_hd.sum()), "; in the files' CellTypes:", int(adn_file.sum()),
"; agreement:", (adn_hd == adn_file).mean().round(3), "; paper: 97")
HA = np.concatenate([d["h"][:, d["hd"]] for d in ADN], 1)
HA_rev = np.concatenate([d["h_rev"][:, d["hd"]] for d in ADN], 1)
HA10 = np.concatenate([d["h10"][:, d["hd"]] for d in ADN], 1)
HA1 = np.concatenate([d["h1"][:, d["hd"]] for d in ADN], 1); HA2 = np.concatenate([d["h2"][:, d["hd"]] for d in ADN], 1)
adn_sess = np.concatenate([[d["session"]] * int(d["hd"].sum()) for d in ADN])
# contamination scores, ADN-HD
adn_cq = np.concatenate([[posub.contamination(sp) for sp, k in zip(a["spikes"], d["hd"]) if k] for a, d in zip(A, ADN)])
ADN sessions: 19 ; units: 220 ; marked ADN: 220
ADN-HD cells recomputed: 97 ; in the files' CellTypes: 97 ; agreement: 1.0 ; paper: 97
fig3c-adn-broader¶
ADN-HD cells have broader tuning curves than PoSub-HD cells.
Width as the code computes it (full width at half maximum of the peak-aligned, normalized curve, its two sides averaged); Mann-Whitney U, PoSub-HD against ADN-HD.
w_pos, off_pos = posub.peak_width(H[:, hdm])
w_adn, off_adn = posub.peak_width(HA)
show(ranksum(w_adn, w_pos), "width (deg), ADN-HD vs PoSub-HD:")
print("paper: n = 1602 and 97, z = 11.4")
print("ADN-HD median width per session:", pd.Series(w_adn).groupby(adn_sess).median().to_dict())
width (deg), ADN-HD vs PoSub-HD: n_a=97, n_b=1512, median_a=62, median_b=42, z=11.58, p=3.681e-31
paper: n = 1602 and 97, z = 11.4
ADN-HD median width per session: {'A1104-180328': 77.0, 'A1104-180329': 80.0, 'A1106-180418': 80.0, 'A1106-180419': 63.0, 'A1109-180615': 48.0, 'A1109-180619': 54.0, 'A1110-180612': 53.0, 'A1110-180613': 49.0, 'A1110-180617': 61.0, 'A1110-180619': 72.0, 'A1110-180621': 68.0, 'A1110-180626': 59.0, 'A1112-180709': 77.0, 'A1112-180711': 82.0, 'A1113-180813': 46.0, 'A1113-180814': 48.0, 'A1113-180816': 52.0, 'B0201-210922': 61.0, 'B0703-211208': 59.0}
/home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:328: RuntimeWarning: Mean of empty slice off[i] = x[wi - 1: n - wi].mean() /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/numpy/_core/_methods.py:142: RuntimeWarning: invalid value encountered in scalar divide ret = ret.dtype.type(ret / rcount)
Result: reproduced. ADN-HD tuning curves (97 cells, the figshare count) are wider than PoSub-HD curves (median width 62 against 42 degrees, z = 11.58, p = 4e-31; paper z = 11.4).
fig3e-fourier-signatures-differ¶
The average Fourier signatures of PoSub-FS, PoSub-HD and ADN-HD populations differ, with ADN-HD cells having more power at low and less at high components than PoSub-HD and PoSub-FS cells, which are alike.
Repeated-measures ANOVA of the ten normalized powers, component by cell type interaction, cell type as a numeric covariate (code: FS 1, PoSub-HD 2, ADN-HD 3) and, as a variant, as a categorical factor. Direction readings (a) ADN mean above both PoSub means at component 1 and below both summed over components 3 to 10; (b) ADN cumulative signature above both at every component 1 to 9. Both power readings.
FTA = posub.fourier(HA)
NA = {rd: posub.norm10(posub.power(FTA, rd)) for rd in ("sq", "rms")}
for rd in ("sq", "rms"):
Y = np.r_[N10[rd][:, fsm].T, N10[rd][:, hdm].T, NA[rd].T]
g = np.r_[np.ones(fsm.sum()), 2 * np.ones(hdm.sum()), 3 * np.ones(NA[rd].shape[1])]
F, d1, d2, p = posub.rm_interaction(Y, g[:, None])
Fc, c1, c2, pc = posub.rm_interaction(Y, np.c_[g == 2, g == 3].astype(float))
Y2 = Y[g < 3]; F2, e1, e2, p2 = posub.rm_interaction(Y2, g[g < 3][:, None])
mf, mh, ma = N10[rd][:, fsm].mean(1), N10[rd][:, hdm].mean(1), NA[rd].mean(1)
ra = ma[0] > max(mf[0], mh[0]) and ma[2:].sum() < min(mf[2:].sum(), mh[2:].sum())
rb = bool(np.all(np.cumsum(ma)[:9] > np.maximum(np.cumsum(mf)[:9], np.cumsum(mh)[:9])))
print(f"power {rd}: numeric covariate F({d1},{d2}) = {F:.2f}, p = {p:.3g}; categorical F({c1},{c2}) = {Fc:.2f}, p = {pc:.3g}; "
f"PoSub-FS vs PoSub-HD only F({e1},{e2}) = {F2:.2f}, p = {p2:.3g}")
print(f" mean signatures FS {mf.round(3)}\n PoSub-HD {mh.round(3)}\n ADN-HD {ma.round(3)}")
print(f" reading (a) {ra}; reading (b) {rb}")
print("paper: F(9,19116) = 2.76, p = 0.003")
power sq: numeric covariate F(9,18333) = 1.09, p = 0.365; categorical F(18,18324) = 24.70, p = 4.96e-82; PoSub-FS vs PoSub-HD only F(9,17460) = 11.37, p = 6.15e-18 mean signatures FS [0.398 0.267 0.153 0.093 0.038 0.02 0.013 0.008 0.005 0.005] PoSub-HD [0.358 0.268 0.172 0.101 0.052 0.026 0.012 0.006 0.003 0.002] ADN-HD [0.509 0.292 0.123 0.049 0.016 0.006 0.003 0.001 0.001 0.001] reading (a) True; reading (b) True power rms: numeric covariate F(9,18333) = 10.68, p = 1.08e-16; categorical F(18,18324) = 44.63, p = 2.36e-155; PoSub-FS vs PoSub-HD only F(9,17460) = 21.09, p = 7.2e-36 mean signatures FS [0.271 0.212 0.155 0.116 0.073 0.051 0.039 0.032 0.026 0.024] PoSub-HD [0.248 0.213 0.167 0.126 0.087 0.059 0.04 0.027 0.019 0.013] ADN-HD [0.348 0.259 0.162 0.095 0.051 0.03 0.02 0.013 0.012 0.01 ] reading (a) True; reading (b) True paper: F(9,19116) = 2.76, p = 0.003
Result: partially reproduced. The direction holds under both direction readings and both power readings: the ADN-HD signature has more power at component 1 than both PoSub signatures (0.509 against 0.398 and 0.358, squared modulus) and less over components 3 to 10, and its cumulative signature is above both at components 1 to 9, the ones tested (at component 10 all three equal 1). The interaction is not significant under the main reading, squared-modulus power with cell type as the code's numeric covariate (F(9,18333) = 1.09, p = 0.365; paper F(9,19116) = 2.76, p = 0.003). It is significant with root-mean-square power (F = 10.68, p = 1e-16) and with cell type as a categorical factor (F(18,18324) = 24.7, p = 5e-82, and 44.6, p = 2e-155). The criterion's test thus holds under some readings and not under the one that follows the code. A numeric code tests a linear trend from FS through PoSub-HD to ADN-HD, and FS and ADN-HD both lie above PoSub-HD at component 1, which such a trend does not capture. The PoSub-only ANOVA, a variant, gives F(9,17460) = 11.4, p = 6e-18.
fig3f-posub-hd-closer-to-fs¶
Individual PoSub-HD cell Fourier spectra are closer to the average PoSub-FS signature than ADN-HD cell spectra are.
KL divergence between each HD cell's ten-component spectrum and the mean PoSub-FS spectrum, both directions and both power readings; Mann-Whitney U, PoSub-HD against ADN-HD.
def fs_distance(rd, Npos, Nadn, Nfs, direction):
S_ = Nfs.mean(1, keepdims=True)
f = (lambda c: posub.kl(c, S_)) if direction == "cell||FS" else (lambda c: posub.kl(np.broadcast_to(S_, c.shape), c))
return f(Npos), f(Nadn)
F3F = {}
for rd in ("sq", "rms"):
for direction in ("cell||FS", "FS||cell"):
dp, da = fs_distance(rd, N10[rd][:, hdm], NA[rd], N10[rd][:, fsm], direction)
r = ranksum(dp, da); F3F[(rd, direction)] = r
show(r, f"power {rd}, D_KL({direction}): PoSub-HD vs ADN-HD")
print("paper: p < 1e-9, n = 1602 and 97")
power sq, D_KL(cell||FS): PoSub-HD vs ADN-HD n_a=1512, n_b=97, median_a=0.08189, median_b=0.1332, z=-2.652, p=0.00801 power sq, D_KL(FS||cell): PoSub-HD vs ADN-HD n_a=1512, n_b=97, median_a=0.1057, median_b=0.2327, z=-5.33, p=9.845e-08 power rms, D_KL(cell||FS): PoSub-HD vs ADN-HD n_a=1512, n_b=97, median_a=0.04999, median_b=0.09001, z=-6.011, p=1.843e-09 power rms, D_KL(FS||cell): PoSub-HD vs ADN-HD n_a=1512, n_b=97, median_a=0.06313, median_b=0.1249, z=-6.971, p=3.149e-12 paper: p < 1e-9, n = 1602 and 97
Result: partially reproduced. PoSub-HD spectra are closer to the mean FS spectrum than ADN-HD spectra under all four readings (median KL 0.082 against 0.133, 0.106 against 0.233, 0.050 against 0.090 and 0.063 against 0.125), with p = 0.008 (squared modulus, D_KL(cell||FS), the Methods' direction), 1e-7, 1.8e-9 and 3e-12. The p-value is below 1e-9 under one reading only, the criterion's partial outcome.
ed3k-robust-to-controls¶
The smaller distance of PoSub-HD than ADN-HD spectra from the PoSub-FS signature holds after excluding multipeaked HD cells, after excluding noisy clusters, and with a 10 cm/s speed threshold; multipeak scores differ between ADN-HD and PoSub-HD cells.
The Fig. 3f test under each control and each reading. Multipeak score below 0.05 (HD cells; code), contamination score at most 0.2 (all cells, as the code's onlyGoodClusters filters HD, FS and viable cells), curves at speed above 10 cm/s (all cells).
The statement was amended after a reader's answer (questions/000939.yaml): the paper gives no direction for the multipeak test, so the criterion accepts it at p < 0.05 in either direction; the three Fourier distance tests keep the paper's direction (PoSub-HD closer, z < 0 here).
U["cq"] = np.nan
for si, s in enumerate(S):
idx = U.index[(U.sess == si) & (U.hd | U.fs)]
for j in idx:
U.loc[j, "cq"] = posub.contamination(s["spikes"][U.unit[j]])
cq_pos = U.cq.values
print("contamination score above 0.2: PoSub-HD", int((cq_pos[hdm] > 0.2).sum()), "of", int(hdm.sum()),
"; FS", int((cq_pos[fsm] > 0.2).sum()), "of", int(fsm.sum()), "; ADN-HD", int((adn_cq > 0.2).sum()), "of", len(adn_cq))
_, off_pos = posub.peak_width(H[:, hdm]); _, off_adn = posub.peak_width(HA)
H10 = per_unit("sq3_v10")
rows = []
for rd in ("sq", "rms"):
N10_v = posub.norm10(posub.power(posub.fourier(H10), rd)); NA_v = posub.norm10(posub.power(posub.fourier(HA10), rd))
for direction in ("cell||FS", "FS||cell"):
# multipeak
dp, da = fs_distance(rd, N10[rd][:, hdm][:, off_pos < 0.05], NA[rd][:, off_adn < 0.05], N10[rd][:, fsm], direction)
rows.append(dict(control="multipeak < 0.05", power=rd, direction=direction, **ranksum(dp, da)))
good_f = fsm & (cq_pos <= 0.2)
dp, da = fs_distance(rd, N10[rd][:, hdm & (cq_pos <= 0.2)], NA[rd][:, adn_cq <= 0.2], N10[rd][:, good_f], direction)
rows.append(dict(control="contamination <= 0.2", power=rd, direction=direction, **ranksum(dp, da)))
dp, da = fs_distance(rd, N10_v[:, hdm], NA_v, N10_v[:, fsm], direction)
rows.append(dict(control="speed > 10 cm/s", power=rd, direction=direction, **ranksum(dp, da)))
K = pd.DataFrame(rows)
print(K.to_string())
show(ranksum(off_adn, off_pos), "multipeak score, ADN-HD vs PoSub-HD:")
print("cells with no width (no side below half maximum): PoSub-HD", int(np.isnan(off_pos).sum()), "ADN-HD", int(np.isnan(off_adn).sum()))
mp = posub.ranksum(off_adn, off_pos)
print("criterion: Fourier tests p < 0.05 and PoSub-HD closer under every reading:", bool(((K.p < 0.05) & (K.z < 0)).all()),
"; multipeak test p < 0.05, either direction:", mp["p"] < 0.05)
print("paper: multipeak n = 789 and 69, z = 10.4; contamination n = 1472 and 83, z = 5.91; speed n = 1602 and 97, z = 6.46; multipeak scores z = 4.25")
contamination score above 0.2: PoSub-HD 113 of 1512 ; FS 11 of 430 ; ADN-HD 16 of 97
control power direction n_a n_b median_a median_b z p
0 multipeak < 0.05 sq cell||FS 762 69 0.036240 0.096549 -7.660942 1.849506e-14
1 contamination <= 0.2 sq cell||FS 1399 81 0.078335 0.129418 -2.347373 1.891312e-02
2 speed > 10 cm/s sq cell||FS 1512 97 0.092475 0.153157 -4.208819 2.568371e-05
3 multipeak < 0.05 sq FS||cell 762 69 0.050364 0.189815 -9.653824 4.747713e-22
4 contamination <= 0.2 sq FS||cell 1399 81 0.101140 0.230193 -5.044104 4.559709e-07
5 speed > 10 cm/s sq FS||cell 1512 97 0.137695 0.303498 -7.247645 4.244340e-13
6 multipeak < 0.05 rms cell||FS 762 69 0.031576 0.084124 -9.440657 3.713776e-21
7 contamination <= 0.2 rms cell||FS 1399 81 0.048760 0.093335 -5.781593 7.405522e-09
8 speed > 10 cm/s rms cell||FS 1512 97 0.057986 0.115558 -8.236841 1.769835e-16
9 multipeak < 0.05 rms FS||cell 762 69 0.036531 0.114327 -9.562691 1.150244e-21
10 contamination <= 0.2 rms FS||cell 1399 81 0.061566 0.123873 -6.745304 1.528479e-11
11 speed > 10 cm/s rms FS||cell 1512 97 0.073344 0.167866 -9.244523 2.365413e-20
multipeak score, ADN-HD vs PoSub-HD: n_a=97, n_b=1511, median_a=0.02171, median_b=0.04912, z=-4.175, p=2.986e-05
cells with no width (no side below half maximum): PoSub-HD 1 ADN-HD 0
criterion: Fourier tests p < 0.05 and PoSub-HD closer under every reading: True ; multipeak test p < 0.05, either direction: True
paper: multipeak n = 789 and 69, z = 10.4; contamination n = 1472 and 83, z = 5.91; speed n = 1602 and 97, z = 6.46; multipeak scores z = 4.25
/home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:328: RuntimeWarning: Mean of empty slice off[i] = x[wi - 1: n - wi].mean() /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/numpy/_core/_methods.py:142: RuntimeWarning: invalid value encountered in scalar divide ret = ret.dtype.type(ret / rcount) /home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:351: RuntimeWarning: invalid value encountered in divide return q / q.sum(0)
Result: reproduced. PoSub-HD spectra stay closer to the mean FS spectrum than ADN-HD spectra with multipeaked cells excluded (762 and 69 cells; paper 789 and 69), with noisy clusters excluded (1,399 and 81; paper 1,472 and 83) and at 10 cm/s (1,512 and 97), under all four readings (p from 0.019 to 5e-22). The multipeak scores differ (z = -4.18, p = 3e-5; paper z = 4.25), with ADN-HD cells scoring lower than PoSub-HD cells (median 0.022 against 0.049), as the paper's counts of cells kept below 0.05 (69 of 97 ADN-HD, 789 of 1,602 PoSub-HD) imply. The claim was amended after a reader's answer to say only that the scores differ, and the criterion accepts the multipeak test in either direction; before the amendment the statement had ADN-HD scores higher, and the outcome was partially reproduced.
ed3g-fs-spectra-heterogeneous, ADN part (reported, not judged)¶
for rd in ("sq", "rms"):
ma = NA[rd].mean(1, keepdims=True)
da_e = np.sqrt(((NA[rd] - ma) ** 2).sum(0)); da_k = posub.kl(NA[rd], ma)
for meas, da in (("euclid", da_e), ("kl", da_k)):
d = DIST[(rd, meas)]
F1 = stats.f_oneway(d[fsm], d[hdm], da)
r1, r2 = ranksum(d[fsm], da), ranksum(d[hdm], da)
print(f"power {rd}, {meas}: one-way ANOVA F(2,{fsm.sum() + hdm.sum() + len(da) - 3}) = {F1.statistic:.1f}, p = {F1.pvalue:.3g}; "
f"FS vs ADN z/3 = {r1['z'] / 3:.2f}, p x 3 = {min(1, 3 * r1['p']):.3g}; PoSub-HD vs ADN z/3 = {r2['z'] / 3:.2f}, p x 3 = {min(1, 3 * r2['p']):.3g}")
print("paper: F(2,2123) = 169, p = 5e-69; z = 3.40 and 0.45")
power sq, euclid: one-way ANOVA F(2,2036) = 184.3, p = 2.78e-74; FS vs ADN z/3 = 3.42, p x 3 = 2.92e-24; PoSub-HD vs ADN z/3 = 0.48, p x 3 = 0.457 power sq, kl: one-way ANOVA F(2,2036) = 171.2, p = 1.85e-69; FS vs ADN z/3 = 3.92, p x 3 = 2.14e-31; PoSub-HD vs ADN z/3 = 0.75, p x 3 = 0.0711 power rms, euclid: one-way ANOVA F(2,2036) = 174.3, p = 1.38e-70; FS vs ADN z/3 = 3.27, p x 3 = 3.2e-22; PoSub-HD vs ADN z/3 = -0.02, p x 3 = 1 power rms, kl: one-way ANOVA F(2,2036) = 143.8, p = 3.98e-59; FS vs ADN z/3 = 3.52, p x 3 = 1.32e-25; PoSub-HD vs ADN z/3 = 0.06, p x 3 = 1 paper: F(2,2123) = 169, p = 5e-69; z = 3.40 and 0.45
Result: reproduced. FS spectra are farther from their population mean than PoSub-HD spectra under all four readings (z/3 from 5.78 to 5.97, p x 3 below 1e-66; paper z = 5.71). With the figshare ADN-HD cells the one-way ANOVA gives F(2,2036) = 184 (Euclidean, squared modulus; paper F(2,2123) = 169), FS against ADN z/3 = 3.42 and PoSub-HD against ADN z/3 = 0.48 (paper 3.40 and 0.45), reported and not judged.
ed4cd-isomap-posub-central¶
In a joint Isomap embedding of tuning curve auto-correlograms, PoSub-HD cells lie nearer the centre of the PoSub-FS distribution than ADN-HD cells.
Isomap (12 neighbours, 2 dimensions) on the cross-validated auto-correlograms of FS, PoSub-HD and ADN-HD cells together; Euclidean distance of each HD cell to the mean FS position; Mann-Whitney U (criterion; Methods). Isomap from scikit-learn, not the MATLAB toolbox the paper used.
from sklearn.manifold import Isomap
ACA = posub.xval_autocorr(HA1, HA2)
Xac = np.c_[AC[:, fsm], AC[:, hdm], ACA].T
grp = np.r_[np.zeros(fsm.sum()), np.ones(hdm.sum()), 2 * np.ones(ACA.shape[1])]
okr = np.isfinite(Xac).all(1)
print("cells dropped for an undefined auto-correlogram (flat half):", int((~okr).sum()), "by group", np.bincount(grp[~okr].astype(int), minlength=3))
emb = Isomap(n_neighbors=12, n_components=2).fit_transform(Xac[okr]); g = grp[okr]
c = emb[g == 0].mean(0); dist = np.linalg.norm(emb - c, axis=1)
show(ranksum(dist[g == 1], dist[g == 2]), "distance to FS centre, PoSub-HD vs ADN-HD:")
print("paper: n = 1602 and 97, z = 4.26")
fig, ax = plt.subplots(figsize=(4.5, 4))
for k, col in ((0, "tab:orange"), (1, "tab:purple"), (2, "magenta")):
ax.scatter(*emb[g == k].T, s=3, c=col, label=["FS", "PoSub-HD", "ADN-HD"][k])
ax.legend(); plt.show()
print("embedding range", emb.min(0).round(2), emb.max(0).round(2))
cells dropped for an undefined auto-correlogram (flat half): 0 by group [0 0 0]
distance to FS centre, PoSub-HD vs ADN-HD: n_a=1512, n_b=97, median_a=2.337, median_b=3.394, z=-2.934, p=0.00335 paper: n = 1602 and 97, z = 4.26
embedding range [-9.61 -8.65] [10.41 10. ]
Result: reproduced. In the joint Isomap embedding of auto-correlograms, PoSub-HD cells lie nearer the mean FS position than ADN-HD cells (median distance 2.34 against 3.39, z = -2.93, p = 0.0034, n = 1,512 and 97; paper z = 4.26).
fig4f-synaptic-weights-dispersed¶
Putative synaptic weights from PoSub-HD to PoSub-FS cells, estimated from spike train cross-correlograms, are heavy-tailed with an s.d. larger than their mean, in the range where simulations predict the radial symmetries of FS tuning.
Cross-correlograms over the whole recording in 0.2 ms bins, HD cell as reference; baseline from a partially hollowed Gaussian (cch_conv, 41 bins); the 99th percentile of the Poisson at the largest predicted count within 0 to 2 ms; at least two consecutive bins above it; peak at positive lag, latency above 0.5 ms, no peak at negative lag (criterion; GetSynapticPairs.m, FindSynapse_Exc, Figure_SynapticConnections.m). Weight strZ (peak minus prediction over the square root of the prediction); strR is a variant. CCG span plus or minus 50 ms (Methods) and plus or minus 25 ms (code).
t_syn = time.time()
syn_rows = []
for si, s in enumerate(S):
u = U[U.sess == si]
hd_idx = np.where(u.hd.values)[0]; fs_idx = np.where(u.fs.values)[0]
for i in hd_idx:
for j in fs_idx:
row = dict(subject=s["subject"], hd=i, fs=j)
for span, nb in (("50", 500), ("25", 250)):
c = posub.ccg_counts(s["spikes"][i], s["spikes"][j], nbins=nb)
ln, zn, lp, zp, rp = posub.find_synapse_exc(c)
row.update({f"lat_neg{span}": ln, f"z_neg{span}": zn, f"lat{span}": lp, f"z{span}": zp, f"r{span}": rp})
syn_rows.append(row)
SY = pd.DataFrame(syn_rows)
print(f"HD-FS pairs: {len(SY)}; CCGs in {time.time() - t_syn:.0f} s")
for span in ("50", "25"):
syn = SY[SY[f"z{span}"].notna() & SY[f"z_neg{span}"].isna() & (SY[f"lat{span}"] > 0.5)]
for w in ("z", "r"):
x = syn[f"{w}{span}"]
print(f"span +-{span} ms, weight str{w.upper()}: {len(syn)} connections ({len(syn) / len(SY):.3%} of pairs); mean {x.mean():.2f}, "
f"s.d. {x.std():.2f}, CV {x.std() / x.mean():.2f}, skewness {stats.skew(x):.2f}, median {x.median():.2f}, max {x.max():.2f}")
print(" connections per subject:", syn.groupby("subject").size().to_dict())
fig, ax = plt.subplots(figsize=(4, 3))
x = SY[SY.z50.notna() & SY.z_neg50.isna() & (SY.lat50 > 0.5)].z50
ax.hist(x, 40); ax.set_xlabel("strZ"); plt.show()
print("plotted range:", x.min().round(2), x.max().round(2))
HD-FS pairs: 24498; CCGs in 144 s
span +-50 ms, weight strZ: 1176 connections (4.800% of pairs); mean 8.58, s.d. 16.32, CV 1.90, skewness 9.34, median 4.60, max 327.37
span +-50 ms, weight strR: 1176 connections (4.800% of pairs); mean 1.73, s.d. 1.23, CV 0.71, skewness 5.82, median 1.37, max 18.31
connections per subject: {'A1808': 21, 'A1813': 33, 'A1815': 65, 'A1821': 47, 'A1824': 54, 'A3701': 21, 'A3702': 25, 'A3703': 84, 'A3705': 58, 'A3706': 81, 'A3707': 95, 'A3708': 8, 'A3709': 26, 'A3710': 9, 'A3711': 87, 'A3712': 59, 'A3713': 13, 'A3716': 60, 'A3717': 45, 'A3723': 11, 'A3728': 7, 'A3730': 20, 'A5505': 30, 'A5506': 17, 'A5507': 35, 'A5508': 16, 'A5801': 18, 'A6211': 29, 'A6215': 20, 'A6216': 47, 'A9701': 35}
span +-25 ms, weight strZ: 1176 connections (4.800% of pairs); mean 8.58, s.d. 16.32, CV 1.90, skewness 9.34, median 4.60, max 327.37
span +-25 ms, weight strR: 1176 connections (4.800% of pairs); mean 1.73, s.d. 1.23, CV 0.71, skewness 5.82, median 1.37, max 18.31
connections per subject: {'A1808': 21, 'A1813': 33, 'A1815': 65, 'A1821': 47, 'A1824': 54, 'A3701': 21, 'A3702': 25, 'A3703': 84, 'A3705': 58, 'A3706': 81, 'A3707': 95, 'A3708': 8, 'A3709': 26, 'A3710': 9, 'A3711': 87, 'A3712': 59, 'A3713': 13, 'A3716': 60, 'A3717': 45, 'A3723': 11, 'A3728': 7, 'A3730': 20, 'A5505': 30, 'A5506': 17, 'A5507': 35, 'A5508': 16, 'A5801': 18, 'A6211': 29, 'A6215': 20, 'A6216': 47, 'A9701': 35}
plotted range: 2.55 327.37
Result: reproduced. 1,176 putative HD-to-FS connections are found among 24,498 pairs (4.8%). Their strZ weights have mean 8.58 and s.d. 16.32 (CV 1.90) and skewness 9.34, heavy-tailed with an s.d. above the mean. The strR weights, a variant, have CV 0.71. The CCG spans of 50 and 25 ms give identical results, as they must: the test window (0 to 2 ms) and the baseline kernel (41 bins of 0.2 ms) never reach either edge, so the span variant does not test anything.
Optogenetic sessions¶
Sessions with TRN injections in VGAT-Cre mice and light over ADN: ArchT (virus field names
ArchT and the thalamic reticular nucleus) and control (GFP only). Light ON is each 1 s pulse,
OFF the 4 s before each onset (Methods; code offStart = 4). Cells are the square-epoch classes
above, kept when their rates in ON, OFF and the square epoch (first to last tracked sample)
all exceed 0.5 Hz (code, Figure_Optogenetics). Tuning curves in ON and OFF use no speed
threshold (Methods): 1 degree bins with 3 degree smoothing for widths, 6 degree bins without
smoothing for gain. The light pulses are decoded from the files' optogenetic series (helper
light_pulses); their widths and intervals are printed, and they are checked against the
figshare pulse table further below.
def opto_group(s):
v = s["virus"]
if "thalamic reticular" not in v:
return None
return "ArchT" if "ArchT" in v else ("control" if "GFP" in v else None)
def opto_rows(spikes, t, ang, on, base, hd, fs, **meta):
"""One row per unit of an optogenetic session: ON and OFF rates, light modulation, tuning
curves in ON and OFF (no speed threshold), widths, gain factors. on is (240, 2) pulse
times, base the square epoch (first to last tracked sample)."""
off = np.c_[on[:, 0] - 4, on[:, 0]]
m_on = posub.in_intervals(t, on) & np.isfinite(ang); m_off = posub.in_intervals(t, off) & np.isfinite(ang)
h_on, _ = posub.tuning_curves(spikes, t, ang, m_on); h_off, _ = posub.tuning_curves(spikes, t, ang, m_off)
g_on, _ = posub.tuning_curves(spikes, t, ang, m_on, nbins=60, sd=0); g_off, _ = posub.tuning_curves(spikes, t, ang, m_off, nbins=60, sd=0)
w_on, _ = posub.peak_width(h_on); w_off, _ = posub.peak_width(h_off)
rows = []
for i, sp in enumerate(spikes):
fr_on = np.array([((sp >= a) & (sp < b)).sum() / (b - a) for a, b in on])
fr_off = np.array([((sp >= a) & (sp < b)).sum() / (b - a) for a, b in off])
x, y = g_off[:, i], g_on[:, i]
mx = x.max()
ok = mx > 0 and x.std() > 0 and y.std() > 0
r, p = stats.pearsonr(y, x) if ok else (np.nan, np.nan)
slope, icpt = np.polyfit(x / mx, y / mx, 1) if ok else (np.nan, np.nan)
rows.append(dict(meta, unit=i, hd=bool(hd[i]), fs=bool(fs[i]),
rate_on=fr_on.mean(), rate_off=fr_off.mean(), rate_base=((sp >= base[0]) & (sp < base[1])).sum() / (base[1] - base[0]),
p_light=stats.wilcoxon(fr_on, fr_off).pvalue if np.any(fr_on != fr_off) else 1.0,
r_tc=r, p_tc=p, mult=slope, add=icpt, w_on=w_on[i], w_off=w_off[i], g_off=x, g_on=y))
return rows
OP = []
for si, s in enumerate(S):
grp = opto_group(s)
if grp is None or "pulses" not in s:
continue
u = U[U.sess == si]
OP += opto_rows(s["spikes"], s["hd_t"], s["hd"], s["pulses"][:240], posub.tracked_epoch(s, "wake_square"),
u.hd.values, u.fs.values, subject=s["subject"], group=grp, sess=si)
OP = pd.DataFrame(OP)
OP["gd"] = (OP.rate_on > 0.5) & (OP.rate_off > 0.5) & (OP.rate_base > 0.5)
OP["mod"] = OP.rate_on / OP.rate_off
print("opto sessions:", OP.groupby("group").subject.unique().to_dict())
print("cells kept, by group and class:")
print(OP[OP.gd].groupby("group")[["hd", "fs"]].sum().to_string())
print("paper: ArchT 83 HD and 47 FS cells (5 mice); control 89 HD and 38 FS (3 mice)")
# the light pulses, decoded from the optogenetic series (helper light_pulses): count, width and
# interval between onsets, per session
for si in OP.sess.unique():
pl = S[si]["pulses"]; w = pl[:, 1] - pl[:, 0]; gap = np.diff(pl[:, 0])
print(f"{S[si]['subject']}: {len(pl)} pulses; width (s) min {w.min():.4f}, median {np.median(w):.4f}, max {w.max():.4f}; "
f"onset interval (s) min {gap.min():.3f}, median {np.median(gap):.3f}, max {gap.max():.3f}; first onset {pl[0, 0]:.1f} s, "
f"wake_opto epoch {[tuple(np.round(e, 1)) for e in S[si]['epochs'].get('wake_opto', [])]}")
# the optogenetic series itself (finding format_error): its length against the raw ElectricalSeries,
# its stated rate and unit, the channels read from it, and the values on each
for s in S:
if "ogen_info" not in s:
continue
oi = s["ogen_info"]
print(f"{s['subject']}: optogenetic_stim length {oi['length']:,} at {oi['rate']:.0f} Hz, unit {oi['unit']}; ElectricalSeries length {s['ephys_len']:,} "
f"at {s['ephys_rate']:.0f} Hz ({s['ephys_len'] / s['ephys_rate']:.0f} s); ratio {oi['length'] / s['ephys_len']:.3f}; "
f"channels {oi['channels']}, max per channel {oi['channel_max']}, min {oi['channel_min']}; light channel values {oi['light_values']}")
# the light times share the spike clock: the ON/OFF ratio with the true pulse times against the
# same ratio with every window moved 2.5 s later (ON in the middle of the 4 s gap, OFF across the
# pulse and the gap), over all kept cells
for grp in ("ArchT", "control"):
d = OP[OP.gd & (OP.group == grp)]
sh = []
for si in d.sess.unique():
s = S[si]; on = s["pulses"][:240] + 2.5; off = np.c_[on[:, 0] - 4, on[:, 0]]
for i in d[d.sess == si].unit:
sp = s["spikes"][i]
ron = posub.in_intervals(sp, on).sum() / 240; roff = posub.in_intervals(sp, off).sum() / 960
sh.append(ron / roff if roff > 0 else np.nan)
print(f"{grp}: median ON/OFF rate ratio over {len(d)} kept cells {d['mod'].median():.3f}; with windows moved 2.5 s {np.nanmedian(sh):.3f}")
/home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:328: RuntimeWarning: Mean of empty slice off[i] = x[wi - 1: n - wi].mean() /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/numpy/_core/_methods.py:142: RuntimeWarning: invalid value encountered in scalar divide ret = ret.dtype.type(ret / rcount)
opto sessions: {'ArchT': <StringArray>
['A1813', 'A1815']
Length: 2, dtype: str, 'control': <StringArray>
['A1808', 'A1821', 'A1824']
Length: 3, dtype: str}
cells kept, by group and class:
hd fs
group
ArchT 62 33
control 92 43
paper: ArchT 83 HD and 47 FS cells (5 mice); control 89 HD and 38 FS (3 mice)
A1808: 240 pulses; width (s) min 1.0000, median 1.0001, max 1.0011; onset interval (s) min 5.000, median 5.000, max 305.014; first onset 7795.0 s, wake_opto epoch [(np.float64(7489.7), np.float64(10198.7))]
A1813: 240 pulses; width (s) min 1.0000, median 1.0001, max 1.0011; onset interval (s) min 5.000, median 5.000, max 305.020; first onset 9212.7 s, wake_opto epoch [(np.float64(8908.4), np.float64(11667.1))]
A1815: 240 pulses; width (s) min 1.0000, median 1.0001, max 1.0011; onset interval (s) min 5.000, median 5.000, max 305.015; first onset 7224.7 s, wake_opto epoch [(np.float64(6915.2), np.float64(9622.9))]
A1821: 240 pulses; width (s) min 0.9991, median 1.0001, max 1.0101; onset interval (s) min 4.999, median 5.000, max 305.020; first onset 9149.5 s, wake_opto epoch [(np.float64(8831.8), np.float64(11551.8))]
A1824: 240 pulses; width (s) min 1.0000, median 1.0000, max 1.0011; onset interval (s) min 5.000, median 5.000, max 305.015; first onset 6698.3 s, wake_opto epoch [(np.float64(6393.1), np.float64(9105.5))]
A1808: optogenetic_stim length 203,973,840 at 20000 Hz, unit watts; ElectricalSeries length 203,973,840 at 20000 Hz (10199 s); ratio 1.000; channels 1, max per channel [1.0], min [0.0]; light channel values [0, 1]
A1813: optogenetic_stim length 933,367,680 at 20000 Hz, unit watts; ElectricalSeries length 233,341,920 at 20000 Hz (11667 s); ratio 4.000; channels 4, max per channel [0.0, 0.0, 0.0, 16368.0], min [0.0, 0.0, 0.0, 0.0]; light channel values [0, 16368]
A1815: optogenetic_stim length 769,829,760 at 20000 Hz, unit watts; ElectricalSeries length 192,457,440 at 20000 Hz (9623 s); ratio 4.000; channels 4, max per channel [0.0, 0.0, 0.0, 16368.0], min [0.0, 0.0, 0.0, 0.0]; light channel values [0, 16368]
A1821: optogenetic_stim length 924,145,920 at 20000 Hz, unit watts; ElectricalSeries length 231,036,480 at 20000 Hz (11552 s); ratio 4.000; channels 4, max per channel [0.0, 0.0, 0.0, 16368.0], min [0.0, 0.0, 0.0, 0.0]; light channel values [0, 16368]
A1824: optogenetic_stim length 728,438,400 at 20000 Hz, unit watts; ElectricalSeries length 182,109,600 at 20000 Hz (9105 s); ratio 4.000; channels 4, max per channel [0.0, 0.0, 0.0, 16368.0], min [0.0, 0.0, 0.0, 0.0]; light channel values [0, 16368]
A6211: optogenetic_stim length 893,102,400 at 20000 Hz, unit watts; ElectricalSeries length 223,275,600 at 20000 Hz (11164 s); ratio 4.000; channels 4, max per channel [0.0, 0.0, 0.0, 16368.0], min [0.0, 0.0, 0.0, 0.0]; light channel values [0, 16368]
A6215: optogenetic_stim length 900,662,400 at 20000 Hz, unit watts; ElectricalSeries length 225,165,600 at 20000 Hz (11258 s); ratio 4.000; channels 4, max per channel [0.0, 0.0, 0.0, 16368.0], min [0.0, 0.0, 0.0, 0.0]; light channel values [0, 16368]
A6216: optogenetic_stim length 892,716,480 at 20000 Hz, unit watts; ElectricalSeries length 223,179,120 at 20000 Hz (11159 s); ratio 4.000; channels 4, max per channel [0.0, 0.0, 0.0, 16368.0], min [0.0, 0.0, 0.0, 0.0]; light channel values [0, 16368]
A9701: optogenetic_stim length 910,967,040 at 20000 Hz, unit watts; ElectricalSeries length 227,741,760 at 20000 Hz (11387 s); ratio 4.000; channels 4, max per channel [0.0, 0.0, 0.0, 16368.0], min [0.0, 0.0, 0.0, 0.0]; light channel values [0, 16368]
ArchT: median ON/OFF rate ratio over 123 kept cells 1.402; with windows moved 2.5 s 0.911
control: median ON/OFF rate ratio over 191 kept cells 1.023; with windows moved 2.5 s 0.997
fig5b-posub-rates-increase¶
Thalamic disinhibition raises the firing rates of PoSub-HD and PoSub-FS cells in ArchT mice; the rise of FS cells is larger than in control-virus mice.
for cls in ("fs", "hd"):
a = OP[OP.gd & OP[cls] & (OP.group == "ArchT")]; c = OP[OP.gd & OP[cls] & (OP.group == "control")]
show(ranksum(a["mod"], c["mod"]), f"{cls.upper()} ON/OFF ratio, ArchT vs control:")
show(signrank(a.rate_on, a.rate_off), f"{cls.upper()} ArchT, ON vs OFF rate per cell:")
print(" by ArchT subject, median ratio:", a.groupby("subject")["mod"].median().round(3).to_dict(),
"; control:", c.groupby("subject")["mod"].median().round(3).to_dict())
print("paper: FS n = 47 and 38, z = 4.99; HD n = 83 and 89, z = 1.69")
FS ON/OFF ratio, ArchT vs control: n_a=33, n_b=43, median_a=1.331, median_b=1.04, z=5.413, p=6.388e-08
FS ArchT, ON vs OFF rate per cell: n=33, median_a=31.68, median_b=26.21, z=-4.458, p=8.272e-06
by ArchT subject, median ratio: {'A1813': 1.168, 'A1815': 1.614} ; control: {'A1808': 0.988, 'A1821': 1.042, 'A1824': 1.051}
HD ON/OFF ratio, ArchT vs control: n_a=62, n_b=92, median_a=1.463, median_b=0.9984, z=7.774, p=7.746e-15
HD ArchT, ON vs OFF rate per cell: n=62, median_a=3.514, median_b=2.795, z=-5.949, p=2.7e-09
by ArchT subject, median ratio: {'A1813': 1.37, 'A1815': 1.548} ; control: {'A1808': 0.949, 'A1821': 1.017, 'A1824': 1.039}
paper: FS n = 47 and 38, z = 4.99; HD n = 83 and 89, z = 1.69
Result: reproduced. In the dandiset's two ArchT mice, FS and HD cells raise their rates with light (median ON/OFF ratio 1.33 and 1.46; signed-rank p = 8e-6 and 3e-9), above the three control mice (1.04 and 1.00) and in each ArchT mouse; the FS rise exceeds control (z = 5.41, p = 6e-8, n = 33 and 43; paper z = 4.99 for 47 and 38). The HD comparison with control is also significant here (z = 7.77, n = 62 and 92), where the paper's is not (z = 1.69 for 83 and 89); the criterion does not judge it. The dandiset holds 2 of the paper's 5 ArchT mice; the variant below adds the other three from figshare.
ed7gj-tuning-preserved¶
Thalamic disinhibition preserves the HD tuning of PoSub-HD and PoSub-FS cells.
Two claims of absence on the legends' tests: PoSub-HD width OFF against ON (ArchT), Wilcoxon signed-rank (rank-sum, the code's test, as a variant); PoSub-FS ON-OFF tuning curve correlation (6 degree bins, unsmoothed), ArchT against control, Mann-Whitney U.
a = OP[OP.gd & OP.hd & (OP.group == "ArchT")]
show(signrank(a.w_off, a.w_on), "PoSub-HD width (deg), OFF vs ON, signed-rank:")
show(ranksum(a.w_off, a.w_on), "PoSub-HD width, OFF vs ON, rank-sum (code):")
af = OP[OP.gd & OP.fs & (OP.group == "ArchT")]; cf = OP[OP.gd & OP.fs & (OP.group == "control")]
show(ranksum(af.r_tc, cf.r_tc), "PoSub-FS ON-OFF tuning correlation, ArchT vs control:")
print("paper: width z = 0.30 (n read as 83); FS correlation n = 47 and 38, z = 1.14")
PoSub-HD width (deg), OFF vs ON, signed-rank: n=62, median_a=44, median_b=46, z=-0.6312, p=0.5279 PoSub-HD width, OFF vs ON, rank-sum (code): n_a=62, n_b=62, median_a=44, median_b=46, z=0, p=1 PoSub-FS ON-OFF tuning correlation, ArchT vs control: n_a=33, n_b=43, median_a=0.7199, median_b=0.6908, z=0.4244, p=0.6751 paper: width z = 0.30 (n read as 83); FS correlation n = 47 and 38, z = 1.14
Result: reproduced. Light does not change PoSub-HD tuning widths (44 against 46 degrees, signed-rank p = 0.53, n = 62; rank-sum p = 1.0), and FS ON-OFF tuning correlations do not differ between ArchT and control (median 0.72 against 0.69, p = 0.68, n = 33 and 43). Both are claims of absence tested on 62 and 33 ArchT cells where the paper had 83 and 47; the variant below with the five ArchT mice gives the same (p = 0.24 and 0.95).
fig5e-hd-multiplicative-fs-additive¶
Thalamic gain modulates PoSub-HD cell tuning multiplicatively and PoSub-FS cell tuning additively, so that FS cells receive thalamic input that is uniform across directions.
ON and OFF curves (6 degree bins, unsmoothed, no speed threshold) divided by the OFF maximum; linear fit of ON on OFF, slope multiplicative and intercept additive factor. ArchT cells are excluded (a) when both light modulation (Wilcoxon signed-rank of per-pulse ON against OFF rates) and the ON-OFF tuning correlation are not significant, (b) when either is not, (c) when the correlation is not (code).
sig_l, sig_r = OP.p_light < 0.05, OP.p_tc < 0.05
keep = {"(a) exclude when both fail": sig_l | sig_r, "(b) exclude when either fails": sig_l & sig_r, "(c) correlation only (code)": sig_r}
G = {}
for name, k in keep.items():
ah = OP[OP.gd & OP.hd & (OP.group == "ArchT") & k]; afs = OP[OP.gd & OP.fs & (OP.group == "ArchT") & k]
ch = OP[OP.gd & OP.hd & (OP.group == "control")]; cfs = OP[OP.gd & OP.fs & (OP.group == "control")]
excl = (int((OP.gd & OP.hd & (OP.group == "ArchT")).sum() - len(ah)), int((OP.gd & OP.fs & (OP.group == "ArchT")).sum() - len(afs)))
print(f"{name}: ArchT cells excluded HD {excl[0]}, FS {excl[1]}")
res = {"HD vs FS, mult": ranksum(ah.mult, afs.mult), "HD vs FS, add": ranksum(ah["add"], afs["add"]),
"HD vs control, mult": ranksum(ah.mult, ch.mult), "HD vs control, add": ranksum(ah["add"], ch["add"]),
"FS vs control, mult": ranksum(afs.mult, cfs.mult), "FS vs control, add": ranksum(afs["add"], cfs["add"])}
for k2, v in res.items():
show(v, " " + k2 + ":")
G[name] = res
print("paper: HD vs FS additive z = 7.72, multiplicative z = 5.99; HD vs control additive z = 1.57, multiplicative z = 6.80; FS vs control additive z = 4.16, multiplicative z = 0.36")
(a) exclude when both fail: ArchT cells excluded HD 0, FS 0 HD vs FS, mult: n_a=62, n_b=33, median_a=1.295, median_b=0.874, z=5.604, p=2.139e-08 HD vs FS, add: n_a=62, n_b=33, median_a=0.008556, median_b=0.2803, z=-6.323, p=2.624e-10 HD vs control, mult: n_a=62, n_b=92, median_a=1.295, median_b=0.9825, z=7.593, p=3.171e-14 HD vs control, add: n_a=62, n_b=92, median_a=0.008556, median_b=0.003254, z=1.323, p=0.1866 FS vs control, mult: n_a=33, n_b=43, median_a=0.874, median_b=0.8931, z=0.8541, p=0.396 FS vs control, add: n_a=33, n_b=43, median_a=0.2803, median_b=0.09411, z=3.747, p=0.0001831 (b) exclude when either fails: ArchT cells excluded HD 31, FS 4 HD vs FS, mult: n_a=31, n_b=29, median_a=1.439, median_b=0.8958, z=4.416, p=1.042e-05 HD vs FS, add: n_a=31, n_b=29, median_a=0.01657, median_b=0.2814, z=-5.022, p=5.311e-07 HD vs control, mult: n_a=31, n_b=92, median_a=1.439, median_b=0.9825, z=5.161, p=2.494e-07 HD vs control, add: n_a=31, n_b=92, median_a=0.01657, median_b=0.003254, z=2.447, p=0.01454 FS vs control, mult: n_a=29, n_b=43, median_a=0.8958, median_b=0.8931, z=1.269, p=0.2066 FS vs control, add: n_a=29, n_b=43, median_a=0.2814, median_b=0.09411, z=3.921, p=9.034e-05 (c) correlation only (code): ArchT cells excluded HD 0, FS 1 HD vs FS, mult: n_a=62, n_b=32, median_a=1.295, median_b=0.8787, z=5.474, p=4.507e-08 HD vs FS, add: n_a=62, n_b=32, median_a=0.008556, median_b=0.2794, z=-6.208, p=5.511e-10 HD vs control, mult: n_a=62, n_b=92, median_a=1.295, median_b=0.9825, z=7.593, p=3.171e-14 HD vs control, add: n_a=62, n_b=92, median_a=0.008556, median_b=0.003254, z=1.323, p=0.1866 FS vs control, mult: n_a=32, n_b=43, median_a=0.8787, median_b=0.8931, z=1.103, p=0.2722 FS vs control, add: n_a=32, n_b=43, median_a=0.2794, median_b=0.09411, z=3.599, p=0.0003258 paper: HD vs FS additive z = 7.72, multiplicative z = 5.99; HD vs control additive z = 1.57, multiplicative z = 6.80; FS vs control additive z = 4.16, multiplicative z = 0.36
Result: partially reproduced. Under all three exclusion readings HD multiplicative factors exceed FS ones (p from 1e-5 to 2e-8) and FS additive factors exceed HD ones (p from 5e-7 to 3e-10), HD multiplicative and FS additive factors exceed control (p below 3e-7 and 4e-4), and FS multiplicative factors do not differ from control (p 0.21 to 0.40). HD additive factors do not differ from control under readings (a) and (c) (p = 0.19) but do under reading (b) (p = 0.015), which excludes 31 of 62 HD cells; the criterion asks for every reading. The paper excluded 1 FS cell and no other cell, as readings (a) and (c) do here (0 HD and 0 or 1 FS cells); whether the claims file should keep reading (b) is asked of a reader (questions/000939.yaml).
ed8cd-multiplicative-detectable¶
Simulations of noisy gain on real FS tuning curves show that multiplicative gain would be detectable at the FS cells' tuning correlations, and the FS data are best explained by additive gain.
As Figure_GainSimulation.m: OFF curves (6 degree bins, unsmoothed) of all kept cells; for ArchT cells the simulated ON curve is (OFF + a x max(OFF)) x m with, per cell, no gain (m = 1, a = 0), multiplicative gain (m drawn from N(1.2, 0.33), a = 0) or additive gain (m = 1, a drawn from N(0.2, 0.1)); control cells get ON = OFF. Gaussian noise is added to OFF and ON independently, and factors are computed as for the data. The noise is read two ways: titrated, as the Methods and the criterion say, so that the no-gain simulation's median ON-OFF correlation matches the observed one (the main reading, second block below); and fixed at 1.8 x 1.8 Hz, as the code's noiseMag sets it (first block, a variant). The simulated ArchT FS factors are compared with the control group's, read two ways: the simulated control cells (ON = OFF plus noise, as Figure_GainSimulation.m builds them) and the recorded control factors. The code's choice of FS cells (ArchT cells with a significant tuning correlation, reading (c)) is used. Seeds fixed; the titrated simulation is repeated over 50 seeds.
rng = np.random.default_rng(8)
fs_a = OP[OP.gd & OP.fs & (OP.group == "ArchT") & (OP.p_tc < 0.05)]
fs_c = OP[OP.gd & OP.fs & (OP.group == "control")]
def sim_factors(rows, mult, add, noise=1.8 * 1.8):
out = []
for _, r in rows.iterrows():
x = r.g_off
m = rng.normal(*mult) if mult else 1.0
a = rng.normal(*add) if add else 0.0
y = (x + a * x.max()) * m
xo = x + rng.normal(0, noise, len(x)); yo = y + rng.normal(0, noise, len(x))
mx = xo.max()
sl, ic = np.polyfit(xo / mx, yo / mx, 1)
out.append((sl, ic, np.corrcoef(xo, yo)[0, 1]))
return np.array(out)
ctrl = sim_factors(fs_c, None, None)
print(f"observed ON-OFF correlation of ArchT FS cells: median {fs_a.r_tc.median():.3f}; simulated control: median {np.median(ctrl[:, 2]):.3f}")
for name, mult, add in (("no gain", None, None), ("multiplicative", (1.2, 0.33), None), ("additive", None, (0.2, 0.1))):
sim = sim_factors(fs_a, mult, add)
rm_, ra_ = ranksum(sim[:, 0], ctrl[:, 0]), ranksum(sim[:, 1], ctrl[:, 1])
print(f"{name}: simulated correlation median {np.median(sim[:, 2]):.3f}; multiplicative vs control z = {rm_['z']:.2f}, p = {rm_['p']:.3g} "
f"(medians {rm_['median_a']:.3f}, {rm_['median_b']:.3f}); additive vs control z = {ra_['z']:.2f}, p = {ra_['p']:.3g} "
f"(medians {ra_['median_a']:.3f}, {ra_['median_b']:.3f})")
vm, va = ranksum(sim[:, 0], fs_c.mult), ranksum(sim[:, 1], fs_c["add"])
print(f" variant, against the recorded control FS factors: multiplicative p = {vm['p']:.3g}, additive p = {va['p']:.3g}")
r1, r2 = ranksum(fs_a.mult, fs_c.mult), ranksum(fs_a["add"], fs_c["add"])
print(f"data: FS multiplicative vs control p = {r1['p']:.3g}, additive vs control p = {r2['p']:.3g}")
# Methods reading: the noise is titrated until the simulated correlations resemble the observed
# ones. Noise s.d. (Hz) chosen on a grid so that the no-gain simulation's median ON-OFF
# correlation is closest to the observed median of the ArchT FS cells.
grid = np.arange(1.0, 12.01, 0.25)
med = []
for nz in grid:
rng = np.random.default_rng(80)
med.append(np.median(sim_factors(fs_a, None, None, noise=nz)[:, 2]))
nz = grid[np.argmin(np.abs(np.array(med) - fs_a.r_tc.median()))]
print(f"titrated noise s.d. {nz:.2f} Hz (no-gain median r {med[int(np.argmin(np.abs(np.array(med) - fs_a.r_tc.median())))]:.3f}, observed {fs_a.r_tc.median():.3f})")
rng = np.random.default_rng(81)
ctrl_t = sim_factors(fs_c, None, None, noise=nz)
for name, mult, add in (("no gain", None, None), ("multiplicative", (1.2, 0.33), None), ("additive", None, (0.2, 0.1))):
sim = sim_factors(fs_a, mult, add, noise=nz)
rm_, ra_ = ranksum(sim[:, 0], ctrl_t[:, 0]), ranksum(sim[:, 1], ctrl_t[:, 1])
vm, va = ranksum(sim[:, 0], fs_c.mult), ranksum(sim[:, 1], fs_c["add"])
print(f"titrated, {name}: median r {np.median(sim[:, 2]):.3f}; vs simulated control: multiplicative p = {rm_['p']:.3g}, additive p = {ra_['p']:.3g}; "
f"vs recorded control: multiplicative p = {vm['p']:.3g}, additive p = {va['p']:.3g}")
# how often the criterion holds over 50 seeds (titrated noise, simulated control)
hits = []
for seed in range(50):
rng = np.random.default_rng(1000 + seed)
c_ = sim_factors(fs_c, None, None, noise=nz); mu_ = sim_factors(fs_a, (1.2, 0.33), None, noise=nz); ad_ = sim_factors(fs_a, None, (0.2, 0.1), noise=nz)
hits.append((ranksum(mu_[:, 0], c_[:, 0])["p"] < 0.05, ranksum(ad_[:, 1], c_[:, 1])["p"] < 0.05, ranksum(ad_[:, 0], c_[:, 0])["p"] >= 0.05))
hits = np.array(hits)
print("over 50 seeds: multiplicative detected", hits[:, 0].mean(), "; additive detected", hits[:, 1].mean(),
"; additive simulation without a multiplicative difference", hits[:, 2].mean(), "; all three", hits.all(1).mean())
observed ON-OFF correlation of ArchT FS cells: median 0.730; simulated control: median 0.657 no gain: simulated correlation median 0.818; multiplicative vs control z = 2.12, p = 0.0344 (medians 0.792, 0.656); additive vs control z = -1.86, p = 0.0631 (medians 0.124, 0.181) variant, against the recorded control FS factors: multiplicative p = 0.0572, additive p = 0.437 multiplicative: simulated correlation median 0.816; multiplicative vs control z = 2.82, p = 0.00492 (medians 0.790, 0.656); additive vs control z = -0.82, p = 0.413 (medians 0.121, 0.181) variant, against the recorded control FS factors: multiplicative p = 0.801, additive p = 0.138 additive: simulated correlation median 0.818; multiplicative vs control z = 2.66, p = 0.00802 (medians 0.859, 0.656); additive vs control z = 2.53, p = 0.0116 (medians 0.307, 0.181) variant, against the recorded control FS factors: multiplicative p = 0.162, additive p = 1.73e-06 data: FS multiplicative vs control p = 0.272, additive vs control p = 0.000326 titrated noise s.d. 4.25 Hz (no-gain median r 0.731, observed 0.730) titrated, no gain: median r 0.693; vs simulated control: multiplicative p = 0.114, additive p = 0.0937; vs recorded control: multiplicative p = 0.00257, additive p = 0.0392 titrated, multiplicative: median r 0.701; vs simulated control: multiplicative p = 0.0222, additive p = 0.49; vs recorded control: multiplicative p = 0.0235, additive p = 0.011 titrated, additive: median r 0.710; vs simulated control: multiplicative p = 0.127, additive p = 7.55e-05; vs recorded control: multiplicative p = 0.000755, additive p = 1.29e-09
over 50 seeds: multiplicative detected 0.98 ; additive detected 1.0 ; additive simulation without a multiplicative difference 0.82 ; all three 0.8
Result: reproduced. The criterion asks for noise titrated to the observed correlations, so the code's fixed noise (3.24 Hz) is a variant and not a reading of the criterion. With the noise titrated (4.25 Hz; no-gain median r 0.731 against the observed 0.730) and the simulated control cells, as Figure_GainSimulation.m builds them, simulated multiplicative gain is detected (p = 0.022) and simulated additive gain shows the data's pattern (additive p = 8e-5, multiplicative p = 0.13), as in the data (p = 3e-4 and 0.27); over 50 seeds all three hold in 80% of runs. The titration used seed 80; in the seed-81 run the simulated median correlations are 0.69 to 0.71, slightly below the observed 0.73. Against the recorded control factors the titrated additive simulation also changes the multiplicative factor (p = 0.0008), but this comparison fails its own check: the titrated no-gain simulation already differs from the recorded control (multiplicative p = 0.0026, additive p = 0.039), so it cannot tell one gain type from another. Under the one comparison that passes the no-gain check the criterion holds. The outcome was partially reproduced until round 3 of verification, as it counted the fixed-noise readings, which the criterion does not call for; the verifier agreed with the change to reproduced. Variant with the code's fixed noise (simulated correlations about 0.82, above the observed 0.73): against the simulated control the additive simulation changes the multiplicative factor (p = 0.008), and against the recorded control the multiplicative simulation is not detected (p = 0.80).
ADN optogenetic sessions (figshare, claims with data elsewhere)¶
The ADN recordings with light over ADN are not in the dandiset. The claims file names their source, figshare doi:10.6084/m9.figshare.24921252 version 3, Dataset_6.tar.gz, streamed here. The same tarball holds the PoSub ArchT sessions of three mice (A1801, A1804, A1806) that the dandiset lacks; the PoSub claims are judged on the dandiset, and these mice are added in a variant after the ADN claims. ON is each of the first 240 pulses (Opto_TS.csv), OFF the 4 s before each onset; cells are kept with rates above 0.5 Hz in ON, OFF and the square epoch (code, Figure_Optogenetics). ADN-HD cells as GetCellTypes_ADN (rate above 0.5 Hz, HD information at least 0.2 bits per spike in the square epoch at speed above 2 cm/s); non-HD cells the others above 0.5 Hz.
O = posub.load_adn(file_id=posub.FIGSHARE_OPTO_FILE)
print("sessions in Dataset_6:", [(o["session"], "ADN" if o["adn"].any() else "PoSub", {1: "ArchT", 2: "control"}.get(o.get("group")), len(o["spikes"])) for o in O])
dandi_subjects = {s["subject"] for s in S}
print("PoSub opto sessions in Dataset_6 whose mouse is not in the dandiset:",
[o["session"] for o in O if not o["adn"].any() and o["session"].split("-")[0] not in dandi_subjects])
AO = []
for o in O:
if not o["adn"].any():
continue
o = dict(o); o["eps"] = dict(o["eps"]); o["eps"].setdefault("wake1Ep", o["eps"]["wakeEp"])
tun = posub.adn_tuning(o)
hd = (tun["rate"] > 0.5) & (tun["info"] >= 0.2)
nhd = (tun["rate"] > 0.5) & (tun["info"] < 0.2)
on = o["pulses"][:240, :2]; off = np.c_[on[:, 0] - 4, on[:, 0]]
t, ang = o["hd_t"], o["hd"]
m_on = posub.in_intervals(t, on) & np.isfinite(ang); m_off = posub.in_intervals(t, off) & np.isfinite(ang)
h_on, _ = posub.tuning_curves(o["spikes"], t, ang, m_on); h_off, _ = posub.tuning_curves(o["spikes"], t, ang, m_off)
g_on, _ = posub.tuning_curves(o["spikes"], t, ang, m_on, nbins=60, sd=0); g_off, _ = posub.tuning_curves(o["spikes"], t, ang, m_off, nbins=60, sd=0)
w_on, _ = posub.peak_width(h_on); w_off, _ = posub.peak_width(h_off)
ep = o["eps"]["wakeEp"]; tt = t[posub.in_intervals(t, ep) & np.isfinite(ang)]; ba, bb = tt[0], tt[-1]
for i, sp in enumerate(o["spikes"]):
fr_on = np.array([((sp >= a) & (sp < b)).sum() / (b - a) for a, b in on])
fr_off = np.array([((sp >= a) & (sp < b)).sum() / (b - a) for a, b in off])
x, y = g_off[:, i], g_on[:, i]; mx = x.max()
ok = mx > 0 and x.std() > 0 and y.std() > 0
r, p = stats.pearsonr(y, x) if ok else (np.nan, np.nan)
slope, icpt = np.polyfit(x / mx, y / mx, 1) if ok else (np.nan, np.nan)
AO.append(dict(session=o["session"], mouse=o["session"].split("-")[0], hd=bool(hd[i]), nhd=bool(nhd[i]), file_hd=bool(o["file_hd"][i]),
rate_on=fr_on.mean(), rate_off=fr_off.mean(), rate_base=((sp >= ba) & (sp < bb)).sum() / (bb - ba),
p_light=stats.wilcoxon(fr_on, fr_off).pvalue if np.any(fr_on != fr_off) else 1.0, r_tc=r, p_tc=p,
mult=slope, add=icpt, w_on=w_on[i], w_off=w_off[i]))
AO = pd.DataFrame(AO)
AO["gd"] = (AO.rate_on > 0.5) & (AO.rate_off > 0.5) & (AO.rate_base > 0.5)
AO["mod"] = AO.rate_on / AO.rate_off
print("ADN opto sessions:", AO.session.nunique(), "mice:", AO.mouse.nunique(), "; units", len(AO), "; kept", int(AO.gd.sum()),
"; ADN-HD kept", int((AO.gd & AO.hd).sum()), "; non-HD kept", int((AO.gd & AO.nhd).sum()),
"; agreement of recomputed HD with the files' CellTypes", (AO.hd == AO.file_hd).mean().round(3))
print("paper: 127 thalamic cells including 52 HD cells from 3 mice")
sessions in Dataset_6: [('A1106-180418', 'ADN', 'ArchT', 4), ('A1106-180419', 'ADN', 'ArchT', 8), ('A1110-180612', 'ADN', 'ArchT', 14), ('A1110-180613', 'ADN', 'ArchT', 18), ('A1110-180617', 'ADN', 'ArchT', 12), ('A1110-180619', 'ADN', 'ArchT', 9), ('A1110-180621', 'ADN', 'ArchT', 23), ('A1110-180626', 'ADN', 'ArchT', 25), ('A1113-180813', 'ADN', 'ArchT', 8), ('A1113-180816', 'ADN', 'ArchT', 10), ('A1801-181004', 'PoSub', 'ArchT', 22), ('A1801-181015', 'PoSub', 'ArchT', 16), ('A1804-181225', 'PoSub', 'ArchT', 14), ('A1804-190104', 'PoSub', 'ArchT', 16), ('A1806-190307', 'PoSub', 'ArchT', 13), ('A1806-190314', 'PoSub', 'ArchT', 12), ('A1808-190807', 'PoSub', 'control', 101), ('A1813-190408', 'PoSub', 'ArchT', 70), ('A1815-190523', 'PoSub', 'ArchT', 65), ('A1821-190910', 'PoSub', 'control', 74), ('A1824-190827', 'PoSub', 'control', 70)]
PoSub opto sessions in Dataset_6 whose mouse is not in the dandiset: ['A1801-181004', 'A1801-181015', 'A1804-181225', 'A1804-190104', 'A1806-190307', 'A1806-190314']
/home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:328: RuntimeWarning: Mean of empty slice off[i] = x[wi - 1: n - wi].mean() /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/numpy/_core/_methods.py:142: RuntimeWarning: invalid value encountered in scalar divide ret = ret.dtype.type(ret / rcount)
ADN opto sessions: 10 mice: 3 ; units 131 ; kept 127 ; ADN-HD kept 53 ; non-HD kept 74 ; agreement of recomputed HD with the files' CellTypes 0.985 paper: 127 thalamic cells including 52 HD cells from 3 mice
ed7-adn-disinhibition¶
Light over ADN in VGAT-Cre mice with ArchT in TRN raises the firing of ADN-HD cells but not of non-HD thalamic cells, without broadening ADN-HD tuning curves.
Modulation ON over OFF rate, ADN-HD against non-HD, Mann-Whitney U, with the ADN-HD median above 1; ADN-HD width (1 degree bins, 3 degree smoothing) OFF against ON, Wilcoxon signed-rank (a claim of absence), the code's rank-sum as a variant.
a = AO[AO.gd & AO.hd]; n_ = AO[AO.gd & AO.nhd]
show(ranksum(a["mod"], n_["mod"]), "ON/OFF ratio, ADN-HD vs non-HD:")
show(signrank(a.rate_on, a.rate_off), "ADN-HD, ON vs OFF rate:")
show(signrank(n_.rate_on, n_.rate_off), "non-HD, ON vs OFF rate:")
show(signrank(a.w_off, a.w_on), "ADN-HD width OFF vs ON, signed-rank:")
show(ranksum(a.w_off, a.w_on), "ADN-HD width OFF vs ON, rank-sum (code):")
print("by mouse, ADN-HD median ratio:", a.groupby("mouse")["mod"].median().round(3).to_dict(), "; non-HD:", n_.groupby("mouse")["mod"].median().round(3).to_dict())
print("paper: n = 52 and 75, z = 6.57; width z = 0.67, p = 0.5")
ON/OFF ratio, ADN-HD vs non-HD: n_a=53, n_b=74, median_a=1.538, median_b=1.058, z=6.845, p=7.791e-12
ADN-HD, ON vs OFF rate: n=53, median_a=8.274, median_b=4.671, z=-6.334, p=2.386e-10
non-HD, ON vs OFF rate: n=74, median_a=6.358, median_b=5.605, z=-1.517, p=0.1294
ADN-HD width OFF vs ON, signed-rank: n=53, median_a=62, median_b=58, z=-0.5119, p=0.6087
ADN-HD width OFF vs ON, rank-sum (code): n_a=53, n_b=53, median_a=62, median_b=58, z=0.2591, p=0.7978
by mouse, ADN-HD median ratio: {'A1106': 1.529, 'A1110': 1.552, 'A1113': 1.51} ; non-HD: {'A1106': 1.026, 'A1110': 1.061, 'A1113': 1.06}
paper: n = 52 and 75, z = 6.57; width z = 0.67, p = 0.5
Result: reproduced. In the figshare ADN sessions (10 sessions, 3 mice, 127 cells kept, 53 ADN-HD), light raises ADN-HD rates more than non-HD rates (median ON/OFF ratio 1.54 against 1.06, z = 6.85, p = 8e-12; paper z = 6.57), non-HD rates do not change significantly (p = 0.13), and ADN-HD widths do not change (62 against 58 degrees, signed-rank p = 0.61; rank-sum p = 0.80).
ed8e-adn-mixed-gain¶
ADN-HD cells show multiplicative factors like those of PoSub-HD cells and larger additive factors.
Factors as in fig5e-hd-multiplicative-fs-additive; ADN-HD (figshare) against the dandiset's ArchT PoSub-HD cells, Mann-Whitney U, under the three exclusion readings.
for name in keep:
kk = {"(a) exclude when both fail": lambda d: (d.p_light < 0.05) | (d.p_tc < 0.05),
"(b) exclude when either fails": lambda d: (d.p_light < 0.05) & (d.p_tc < 0.05),
"(c) correlation only (code)": lambda d: d.p_tc < 0.05}[name]
ad = AO[AO.gd & AO.hd & kk(AO)]
ph = OP[OP.gd & OP.hd & (OP.group == "ArchT") & kk(OP)]
rm_, ra_ = ranksum(ad.mult, ph.mult), ranksum(ad["add"], ph["add"])
print(f"{name}: ADN-HD n = {len(ad)}, PoSub-HD n = {len(ph)}")
show(rm_, " multiplicative, ADN-HD vs PoSub-HD:"); show(ra_, " additive, ADN-HD vs PoSub-HD:")
print("paper: multiplicative z = 0.55, additive z = 5.07")
(a) exclude when both fail: ADN-HD n = 53, PoSub-HD n = 62 multiplicative, ADN-HD vs PoSub-HD: n_a=53, n_b=62, median_a=1.263, median_b=1.295, z=-0.7182, p=0.4744 additive, ADN-HD vs PoSub-HD: n_a=53, n_b=62, median_a=0.03174, median_b=0.008556, z=4.489, p=7.262e-06 (b) exclude when either fails: ADN-HD n = 36, PoSub-HD n = 31 multiplicative, ADN-HD vs PoSub-HD: n_a=36, n_b=31, median_a=1.304, median_b=1.439, z=-0.503, p=0.6194 additive, ADN-HD vs PoSub-HD: n_a=36, n_b=31, median_a=0.04181, median_b=0.01657, z=3.068, p=0.002199 (c) correlation only (code): ADN-HD n = 53, PoSub-HD n = 62 multiplicative, ADN-HD vs PoSub-HD: n_a=53, n_b=62, median_a=1.263, median_b=1.295, z=-0.7182, p=0.4744 additive, ADN-HD vs PoSub-HD: n_a=53, n_b=62, median_a=0.03174, median_b=0.008556, z=4.489, p=7.262e-06 paper: multiplicative z = 0.55, additive z = 5.07
Result: reproduced. Figshare ADN-HD cells and the dandiset's ArchT PoSub-HD cells have similar multiplicative factors (median 1.26 against 1.30, p = 0.47) and ADN-HD additive factors are larger (0.032 against 0.009, p = 7e-6; paper z = 0.55 and 5.07); the same holds under reading (b) (p = 0.62 and 0.002).
PoSub ArchT mice from figshare Dataset_6 (variant for fig5b, ed7gj, fig5e and ed8cd)¶
The dandiset holds 2 of the paper's 5 ArchT mice of the PoSub optogenetic experiment. The
other three (A1801, A1804 and A1806, two sessions each) are in figshare Dataset_6, streamed
above for the ADN claims. The four PoSub optogenetic claims are judged on the dandiset
(the claims file gives them data: dataset); this variant adds the three figshare mice to the
ArchT group, to show how the paper's sample size bears on the results. Cell classes for the
figshare sessions are recomputed as for the dandiset (Methods), from the files' trough-to-peak
and Velocity series; the files' own CellTypes are printed beside them. The dandiset's decoded
light pulses are checked against the figshare pulse table (Opto_TS.csv) for the five mice in
both.
# PoSub ArchT sessions of Dataset_6 whose mice are not in the dandiset; classes recomputed as for
# the dandiset (Methods: trough-to-peak from the file's WaveformFeatures, rate and HD information
# in the square epoch at speed above 2 cm/s from the file's Velocity series)
OPF, cls_rows = [], []
for o in O:
if o["adn"].any() or o.get("group") != 1 or o["session"].split("-")[0] in dandi_subjects:
continue
tun = posub.adn_tuning(o)
ttp = o["mat"]["tr2pk"].ravel().astype(float)
gd = (tun["rate"] >= posub.RATE_MIN) & np.isfinite(ttp)
ex = gd & (ttp > posub.TTP_CUT) & (tun["rate"] < posub.FS_RATE)
fs = gd & (ttp < posub.TTP_CUT) & (tun["rate"] > posub.FS_RATE)
hd = ex & (tun["info"] > posub.HD_INFO_TH)
cls_rows.append(dict(session=o["session"], units=len(o["spikes"]), hd=int(hd.sum()), fs=int(fs.sum()),
file_hd=int(o["file_types"]["hd"].sum()), file_fs=int(o["file_types"]["fs"].sum()),
agree_hd=(hd == o["file_types"]["hd"]).mean().round(3), agree_fs=(fs == o["file_types"]["fs"]).mean().round(3)))
t, ang = o["hd_t"], o["hd"]
ep = o["eps"]["wakeEp"]; tt = t[posub.in_intervals(t, ep) & np.isfinite(ang)]
OPF += opto_rows(o["spikes"], t, ang, o["pulses"][:240, :2], (tt[0], tt[-1]), hd, fs,
subject=o["session"].split("-")[0], group="ArchT", sess=o["session"])
print(pd.DataFrame(cls_rows).to_string())
OPF = pd.DataFrame(OPF)
OPF["gd"] = (OPF.rate_on > 0.5) & (OPF.rate_off > 0.5) & (OPF.rate_base > 0.5)
OPF["mod"] = OPF.rate_on / OPF.rate_off
OP5 = pd.concat([OP, OPF], ignore_index=True)
print("ArchT cells kept with the figshare mice added:", OP5[OP5.gd & (OP5.group == "ArchT")][["hd", "fs"]].sum().to_dict(),
"from", OP5[OP5.group == "ArchT"].subject.nunique(), "mice (paper 83 HD and 47 FS from 5 mice)")
# the dandiset's pulse decoding against the figshare pulse table, for the ArchT mice in both
for o in O:
m = [si for si, s in enumerate(S) if s["subject"] == o["session"].split("-")[0] and "pulses" in s]
if o["adn"].any() or not m:
continue
a_, b_ = S[m[0]]["pulses"][:240], o["pulses"][:240, :2]
print(f"{o['session']}: onset difference dandiset - figshare (s), median {np.median(a_[:, 0] - b_[:, 0]):.4f}, max |.| {np.abs(a_[:, 0] - b_[:, 0]).max():.4f}; "
f"widths {np.median(a_[:, 1] - a_[:, 0]):.4f} and {np.median(b_[:, 1] - b_[:, 0]):.4f}")
ah, ac = OP5[OP5.gd & OP5.hd & (OP5.group == "ArchT")], OP5[OP5.gd & OP5.hd & (OP5.group == "control")]
af, cf = OP5[OP5.gd & OP5.fs & (OP5.group == "ArchT")], OP5[OP5.gd & OP5.fs & (OP5.group == "control")]
print("fig5b, five ArchT mice:")
show(ranksum(af["mod"], cf["mod"]), " FS ON/OFF ratio, ArchT vs control:")
show(ranksum(ah["mod"], ac["mod"]), " HD ON/OFF ratio, ArchT vs control:")
print(" ArchT medians by mouse, FS:", af.groupby("subject")["mod"].median().round(3).to_dict(), "; HD:", ah.groupby("subject")["mod"].median().round(3).to_dict())
print("ed7gj, five ArchT mice:")
show(signrank(ah.w_off, ah.w_on), " PoSub-HD width, OFF vs ON, signed-rank:")
show(ranksum(ah.w_off, ah.w_on), " PoSub-HD width, OFF vs ON, rank-sum (code):")
show(ranksum(af.r_tc, cf.r_tc), " PoSub-FS ON-OFF tuning correlation, ArchT vs control:")
print("fig5e, five ArchT mice:")
sig5 = {"(a)": (OP5.p_light < 0.05) | (OP5.p_tc < 0.05), "(b)": (OP5.p_light < 0.05) & (OP5.p_tc < 0.05), "(c)": OP5.p_tc < 0.05}
for name, k in sig5.items():
a5h, a5f = OP5[OP5.gd & OP5.hd & (OP5.group == "ArchT") & k], OP5[OP5.gd & OP5.fs & (OP5.group == "ArchT") & k]
res = {"HD vs FS, mult": ranksum(a5h.mult, a5f.mult), "HD vs FS, add": ranksum(a5h["add"], a5f["add"]),
"HD vs control, mult": ranksum(a5h.mult, ac.mult), "HD vs control, add": ranksum(a5h["add"], ac["add"]),
"FS vs control, mult": ranksum(a5f.mult, cf.mult), "FS vs control, add": ranksum(a5f["add"], cf["add"])}
print(f" {name}: excluded HD {len(ah) - len(a5h)}, FS {len(af) - len(a5f)}; " +
"; ".join(f"{k2} z = {v['z']:.2f}, p = {v['p']:.2g}" for k2, v in res.items()))
print("ed8cd, five ArchT mice (FS cells with a significant tuning correlation, as the code):")
fs5 = af[af.p_tc < 0.05]
med5 = []
for nz in grid:
rng = np.random.default_rng(80)
med5.append(np.median(sim_factors(fs5, None, None, noise=nz)[:, 2]))
nz5 = grid[np.argmin(np.abs(np.array(med5) - fs5.r_tc.median()))]
rng = np.random.default_rng(82)
ctrl5 = sim_factors(cf, None, None, noise=nz5)
for name, mult, add in (("multiplicative", (1.2, 0.33), None), ("additive", None, (0.2, 0.1))):
sim = sim_factors(fs5, mult, add, noise=nz5)
rm_, ra_ = ranksum(sim[:, 0], ctrl5[:, 0]), ranksum(sim[:, 1], ctrl5[:, 1])
vm, va = ranksum(sim[:, 0], cf.mult), ranksum(sim[:, 1], cf["add"])
print(f" titrated noise {nz5:.2f} Hz, {name}: vs simulated control multiplicative p = {rm_['p']:.3g}, additive p = {ra_['p']:.3g}; "
f"vs recorded control multiplicative p = {vm['p']:.3g}, additive p = {va['p']:.3g}")
r1, r2 = ranksum(fs5.mult, cf.mult), ranksum(fs5["add"], cf["add"])
print(f" data: FS multiplicative vs control p = {r1['p']:.3g}, additive vs control p = {r2['p']:.3g}; n = {len(fs5)} and {len(cf)}")
/home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:328: RuntimeWarning: Mean of empty slice off[i] = x[wi - 1: n - wi].mean() /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/numpy/_core/_methods.py:142: RuntimeWarning: invalid value encountered in scalar divide ret = ret.dtype.type(ret / rcount)
session units hd fs file_hd file_fs agree_hd agree_fs
0 A1801-181004 22 13 5 9 5 0.818 1.000
1 A1801-181015 16 6 4 1 3 0.688 0.938
2 A1804-181225 14 3 3 1 4 0.857 0.929
3 A1804-190104 16 10 0 8 0 0.875 1.000
4 A1806-190307 13 6 3 6 2 0.846 0.923
5 A1806-190314 12 4 2 4 2 1.000 1.000
ArchT cells kept with the figshare mice added: {'hd': 102, 'fs': 50} from 5 mice (paper 83 HD and 47 FS from 5 mice)
A1808-190807: onset difference dandiset - figshare (s), median -0.0001, max |.| 0.0001; widths 1.0001 and 1.0000
A1813-190408: onset difference dandiset - figshare (s), median -0.0001, max |.| 0.0001; widths 1.0001 and 1.0000
A1815-190523: onset difference dandiset - figshare (s), median -0.0001, max |.| 0.0001; widths 1.0001 and 1.0000
A1821-190910: onset difference dandiset - figshare (s), median -0.0001, max |.| 0.0001; widths 1.0001 and 1.0000
A1824-190827: onset difference dandiset - figshare (s), median -0.0000, max |.| 0.0001; widths 1.0000 and 1.0000
fig5b, five ArchT mice:
FS ON/OFF ratio, ArchT vs control: n_a=50, n_b=43, median_a=1.288, median_b=1.04, z=5.602, p=2.167e-08
HD ON/OFF ratio, ArchT vs control: n_a=102, n_b=92, median_a=1.315, median_b=0.9984, z=7.091, p=1.345e-12
ArchT medians by mouse, FS: {'A1801': 1.46, 'A1804': 0.937, 'A1806': 1.194, 'A1813': 1.168, 'A1815': 1.614} ; HD: {'A1801': 1.298, 'A1804': 0.93, 'A1806': 1.144, 'A1813': 1.37, 'A1815': 1.548}
ed7gj, five ArchT mice:
PoSub-HD width, OFF vs ON, signed-rank: n=102, median_a=44, median_b=44, z=-1.186, p=0.2356
PoSub-HD width, OFF vs ON, rank-sum (code): n_a=102, n_b=102, median_a=44, median_b=44, z=1.171, p=0.2418
PoSub-FS ON-OFF tuning correlation, ArchT vs control: n_a=50, n_b=43, median_a=0.7042, median_b=0.6908, z=0.06164, p=0.9539
fig5e, five ArchT mice:
(a): excluded HD 2, FS 0; HD vs FS, mult z = 5.82, p = 5.9e-09; HD vs FS, add z = -8.04, p = 9.3e-16; HD vs control, mult z = 6.28, p = 3.5e-10; HD vs control, add z = 1.58, p = 0.11; FS vs control, mult z = 0.74, p = 0.46; FS vs control, add z = 4.09, p = 4.4e-05
(b): excluded HD 52, FS 7; HD vs FS, mult z = 4.11, p = 3.9e-05; HD vs FS, add z = -6.36, p = 2.1e-10; HD vs control, mult z = 3.90, p = 9.6e-05; HD vs control, add z = 2.40, p = 0.016; FS vs control, mult z = 1.29, p = 0.2; FS vs control, add z = 4.16, p = 3.3e-05
(c): excluded HD 2, FS 2; HD vs FS, mult z = 5.57, p = 2.6e-08; HD vs FS, add z = -7.85, p = 4.3e-15; HD vs control, mult z = 6.28, p = 3.5e-10; HD vs control, add z = 1.58, p = 0.11; FS vs control, mult z = 1.10, p = 0.27; FS vs control, add z = 3.88, p = 0.00011
ed8cd, five ArchT mice (FS cells with a significant tuning correlation, as the code):
titrated noise 4.00 Hz, multiplicative: vs simulated control multiplicative p = 0.0188, additive p = 0.14; vs recorded control multiplicative p = 0.126, additive p = 0.000157
titrated noise 4.00 Hz, additive: vs simulated control multiplicative p = 0.384, additive p = 2.64e-07; vs recorded control multiplicative p = 7.41e-05, additive p = 8.34e-12
data: FS multiplicative vs control p = 0.274, additive vs control p = 0.000106; n = 48 and 43
Result: variant, not judged. With the three figshare ArchT mice added (6 sessions; 102 HD and 50 FS ArchT cells kept from 5 mice, paper 83 and 47; the recomputed classes agree with the files' CellTypes for 69 to 100% of units per session for HD), the four claims come out as on the dandiset alone. fig5b: FS ArchT against control z = 5.60, p = 2e-8, ArchT medians 1.29 (FS) and 1.32 (HD) above the control medians, though one mouse (A1804) has medians below 1 (0.94 and 0.93). ed7gj: HD widths p = 0.24 (signed-rank), FS correlation p = 0.95. fig5e: every condition holds under readings (a) and (c) (HD additive against control p = 0.11), and under (b), which excludes 52 HD cells, HD additive against control gives p = 0.016. ed8cd: titrated noise (4.00 Hz) with the simulated control detects multiplicative gain (p = 0.019) and the additive simulation shows the data's pattern (additive p = 3e-7, multiplicative p = 0.38); against the recorded control the additive simulation changes the multiplicative factor (p = 7e-5). The dandiset's decoded light pulses match the figshare pulse table within 0.1 ms in all five PoSub opto mice in both.
fig6c-coupling-preserved-wake-rem¶
GLM cross-coupling of PoSub-HD pairs, PoSub-FS pairs and PoSub-HD to PoSub-FS pairs is preserved between wakefulness and REM sleep, less for FS pairs than for HD pairs.
Sessions with at least 100 s of REM (Methods). As GetGLM_Forward.m and GetGLM_Reverse.m: spike counts of every unit in 100 ms bins over the whole recording, smoothed with a Gaussian of 1 bin s.d.; the population rate (sum over all units) smoothed alike and z-scored; cell B's train z-scored; Poisson GLM of cell A's smoothed train on the population rate and cell B, within WAKE (square epoch) or REM; the coefficient on B at lag 0, averaged with that of the reversed pair. Cells as Figure_GLMRemWake.m selects them (ixHd = isHd & isGd, ixFs = isFs & isGd: the HD and FS classes, with no rate condition in WAKE or REM). REM is the scored REM inside the home-cage epochs (data check). Only lag 0 is fitted (the claim uses lag 0).
Pairs with fit warnings are removed (criterion). Figure_GLMRemWake.m removes a pair when
glmfit raised any warning (lastwarn) in any of its fits: WAKE or REM, either order, any of
the 201 lags. The fits here mirror glmfit's iteratively reweighted least squares (helper
glmfit_poisson_batch) and flag its two warnings for a Poisson fit: weights ill-conditioned
(a weight below eps^(2/3) of the largest at any iteration) and the iteration limit; a
singular system is flagged too. MATLAB is not run, so that this mirror flags the pairs glmfit
would is an assumption. Shown as variants: the round-1 rule (the final fitted means spanning
more than eps^(-4/3), or no convergence), no removal, and the round-1 selection of cells above
1 Hz in both WAKE and REM, which came from a sanity-check block of GetGLM_Forward.m that runs
after its return statement. The flagged pairs are described (coefficients and REM rates).
from scipy.ndimage import gaussian_filter1d
from concurrent.futures import ThreadPoolExecutor
glm_rows = []
t_glm = time.time()
EPS43 = (4 / 3) * np.log(np.finfo(float).eps)
for si, s in enumerate(S):
remiv = s["sleep"].get("rem", np.zeros((0, 2)))
if (remiv[:, 1] - remiv[:, 0]).sum() < 100:
continue
u = U[U.sess == si]
tmax = max(b for v in s["epochs"].values() for a, b in v)
edges = np.arange(0, tmax + 0.1, 0.1); cen = edges[:-1] + 0.05
Q = gaussian_filter1d(posub.bin_counts(s["spikes"], edges), 1, axis=0)
mua = Q.sum(1); mua = (mua - mua.mean()) / mua.std()
Qz = (Q - Q.mean(0)) / Q.std(0)
wa, wb = s["epochs"]["wake_square"][0]
states = {"WAKE": (cen >= wa) & (cen < wb), "REM": posub.in_intervals(cen, remiv)}
rate_w = np.array([((sp >= wa) & (sp < wb)).sum() / (wb - wa) for sp in s["spikes"]])
rem_dur = (remiv[:, 1] - remiv[:, 0]).sum()
rate_r = np.array([posub.in_intervals(sp, remiv).sum() / rem_dur for sp in s["spikes"]])
# the code's selection (Figure_GLMRemWake.m: ixHd = isHd & isGd, ixFs = isFs & isGd), with no
# rate condition in WAKE or REM; the 1 Hz condition is a variant below
cells = np.where((u.hd | u.fs).values)[0]
res = {}
for st, msk in states.items():
Qm, mm, Zm = Q[msk], mua[msk], Qz[msk]
def fit(a):
others = np.array([c for c in cells if c != a])
b, bad, lim = posub.glmfit_poisson_batch(Qm[:, a], mm, Zm[:, others])
with np.errstate(all="ignore"):
eta = b[0] + np.outer(mm, b[1]) + Zm[:, others] * b[2]
lr = eta.min(0) - eta.max(0) # log of the smallest over the largest fitted mean
return b[2], bad, lim, ~(lr > EPS43)
with ThreadPoolExecutor(4) as ex:
out = list(ex.map(fit, cells))
B = np.full((len(cells), len(cells)), np.nan); W = np.zeros_like(B, bool); L = W.copy(); F = W.copy()
for ia, (b2, bad, lim, fin) in enumerate(out):
col = [k for k in range(len(cells)) if k != ia]
B[ia, col] = b2; W[ia, col] = bad; L[ia, col] = lim; F[ia, col] = fin
res[st] = (B, W, L, F)
cls = np.where(u.hd.values[cells], "HD", "FS")
tcr = posub.zero_lag_r(np.repeat(TC[si]["sq3"][:, cells], len(cells), 1), np.tile(TC[si]["sq3"][:, cells], (1, len(cells)))).reshape(len(cells), len(cells))
for i in range(len(cells)):
for j in range(i + 1, len(cells)):
pair = "-".join(sorted([cls[i], cls[j]], key=lambda c: c != "HD"))
f = lambda st, k: bool(res[st][k][i, j] or res[st][k][j, i])
glm_rows.append(dict(subject=s["subject"], sess=si, ci=int(cells[i]), cj=int(cells[j]), pair=pair,
wake=(res["WAKE"][0][i, j] + res["WAKE"][0][j, i]) / 2, rem=(res["REM"][0][i, j] + res["REM"][0][j, i]) / 2,
tc_r=tcr[i, j], bad_scale=f("WAKE", 1) or f("REM", 1), iter_lim=f("WAKE", 2) or f("REM", 2),
final_ratio=f("WAKE", 3) or f("REM", 3),
rate_rem_min=min(rate_r[cells[i]], rate_r[cells[j]]), rate_wake_min=min(rate_w[cells[i]], rate_w[cells[j]])))
GL = pd.DataFrame(glm_rows)
GL["warn"] = GL.bad_scale | GL.iter_lim | ~np.isfinite(GL.wake) | ~np.isfinite(GL.rem) # glmfit's warnings, as mirrored
GL["old_rule"] = GL.iter_lim | GL.final_ratio | ~np.isfinite(GL.wake) | ~np.isfinite(GL.rem) # round-1 rule
GL["hz1"] = (GL.rate_rem_min > 1) & (GL.rate_wake_min > 1)
print(f"GLMs fitted in {time.time() - t_glm:.0f} s over {GL.subject.nunique()} sessions")
print("paper: HD-HD 34,270, FS-FS 2,492, HD-FS 8,992 (Extended Data Fig. 9d gives 18,992)")
print("pairs by type (code's selection):", GL.groupby("pair").size().to_dict())
print("pairs with a glmfit warning as mirrored (bad scaling at any iteration, or iteration limit, in WAKE or REM, either order):",
GL.groupby("pair").warn.sum().to_dict(), "; of which bad scaling", GL.groupby("pair").bad_scale.sum().to_dict(),
", iteration limit", GL.groupby("pair").iter_lim.sum().to_dict())
print("pairs flagged by the round-1 rule (final fitted means spanning more than eps^(-4/3), or no convergence):", GL.groupby("pair").old_rule.sum().to_dict())
print("both rules agree on", round(float((GL.warn == GL.old_rule).mean()), 4), "of pairs")
# what the flagged pairs are
for p_ in ("HD-HD", "FS-FS", "HD-FS"):
d = GL[GL.pair == p_]
for lab, m_ in (("flagged", d.warn), ("kept", ~d.warn)):
x = d[m_]
if not len(x):
print(f"{p_} {lab}: none"); continue
aw, ar = x.wake.abs(), x.rem.abs()
print(f"{p_} {lab}: n = {len(x)}; |coupling| WAKE median {np.nanmedian(aw):.3g}, max {np.nanmax(aw):.3g}; REM median {np.nanmedian(ar):.3g}, max {np.nanmax(ar):.3g}; "
f"lower REM rate of the pair median {x.rate_rem_min.median():.2f} Hz, below 1 Hz {(x.rate_rem_min < 1).mean():.2f}, below 0.1 Hz {(x.rate_rem_min < 0.1).mean():.2f}")
def corr_table(G, label):
for p_ in ("HD-HD", "FS-FS", "HD-FS"):
d = G[(G.pair == p_) & np.isfinite(G.wake) & np.isfinite(G.rem)]
r, p = stats.pearsonr(d.wake, d.rem)
print(f" {label}, {p_}: Pearson r = {r:.3f}, p = {p:.3g}, n = {len(d)} unordered pairs ({2 * len(d)} ordered); Spearman {stats.spearmanr(d.wake, d.rem)[0]:.3f}")
print("WAKE against REM coupling:")
corr_table(GL[~GL.warn], "main: code's selection, glmfit warnings removed")
corr_table(GL[~GL.old_rule], "variant: code's selection, round-1 rule")
corr_table(GL, "variant: code's selection, nothing removed")
corr_table(GL[GL.hz1 & ~GL.warn], "variant: cells above 1 Hz in WAKE and REM, glmfit warnings removed")
corr_table(GL[GL.hz1 & ~GL.old_rule], "variant: cells above 1 Hz in WAKE and REM, round-1 rule (round-1 main)")
GLc = GL[~GL.warn & np.isfinite(GL.wake) & np.isfinite(GL.rem)]
print("per session r (main), median and range:", {p_: GLc[GLc.pair == p_].groupby("subject").apply(lambda x: stats.pearsonr(x.wake, x.rem)[0] if len(x) > 2 else np.nan).describe()[["min", "50%", "max"]].round(3).to_dict() for p_ in ("HD-HD", "FS-FS", "HD-FS")})
GLMs fitted in 274 s over 27 sessions
paper: HD-HD 34,270, FS-FS 2,492, HD-FS 8,992 (Extended Data Fig. 9d gives 18,992)
pairs by type (code's selection): {'FS-FS': 2683, 'HD-FS': 19626, 'HD-HD': 36406}
pairs with a glmfit warning as mirrored (bad scaling at any iteration, or iteration limit, in WAKE or REM, either order): {'FS-FS': 0, 'HD-FS': 24, 'HD-HD': 1795} ; of which bad scaling {'FS-FS': 0, 'HD-FS': 24, 'HD-HD': 1794} , iteration limit {'FS-FS': 0, 'HD-FS': 22, 'HD-HD': 555}
pairs flagged by the round-1 rule (final fitted means spanning more than eps^(-4/3), or no convergence): {'FS-FS': 0, 'HD-FS': 24, 'HD-HD': 1794}
both rules agree on 1.0 of pairs
HD-HD flagged: n = 1795; |coupling| WAKE median 2.77, max 3.46e+73; REM median 2.05, max 1.03e+73; lower REM rate of the pair median 1.51 Hz, below 1 Hz 0.32, below 0.1 Hz 0.11
HD-HD kept: n = 34611; |coupling| WAKE median 0.281, max 4.17; REM median 0.24, max 10.7; lower REM rate of the pair median 2.10 Hz, below 1 Hz 0.15, below 0.1 Hz 0.01
FS-FS flagged: none
FS-FS kept: n = 2683; |coupling| WAKE median 0.0483, max 0.758; REM median 0.0528, max 1.18; lower REM rate of the pair median 22.32 Hz, below 1 Hz 0.02, below 0.1 Hz 0.00
HD-FS flagged: n = 24; |coupling| WAKE median 0.214, max 0.676; REM median 5.9, max 5.9; lower REM rate of the pair median 0.00 Hz, below 1 Hz 0.96, below 0.1 Hz 0.96
HD-FS kept: n = 19602; |coupling| WAKE median 0.115, max 3.16; REM median 0.11, max 21.9; lower REM rate of the pair median 2.97 Hz, below 1 Hz 0.09, below 0.1 Hz 0.01
WAKE against REM coupling:
main: code's selection, glmfit warnings removed, HD-HD: Pearson r = 0.641, p = 0, n = 34611 unordered pairs (69222 ordered); Spearman 0.762
main: code's selection, glmfit warnings removed, FS-FS: Pearson r = 0.364, p = 1.08e-84, n = 2683 unordered pairs (5366 ordered); Spearman 0.417
main: code's selection, glmfit warnings removed, HD-FS: Pearson r = 0.359, p = 0, n = 19602 unordered pairs (39204 ordered); Spearman 0.514
variant: code's selection, round-1 rule, HD-HD: Pearson r = 0.641, p = 0, n = 34612 unordered pairs (69224 ordered); Spearman 0.762
variant: code's selection, round-1 rule, FS-FS: Pearson r = 0.364, p = 1.08e-84, n = 2683 unordered pairs (5366 ordered); Spearman 0.417
variant: code's selection, round-1 rule, HD-FS: Pearson r = 0.359, p = 0, n = 19602 unordered pairs (39204 ordered); Spearman 0.514
variant: code's selection, nothing removed, HD-HD: Pearson r = -0.000, p = 0.995, n = 36228 unordered pairs (72456 ordered); Spearman 0.759
variant: code's selection, nothing removed, FS-FS: Pearson r = 0.364, p = 1.08e-84, n = 2683 unordered pairs (5366 ordered); Spearman 0.417
variant: code's selection, nothing removed, HD-FS: Pearson r = 0.356, p = 0, n = 19603 unordered pairs (39206 ordered); Spearman 0.514
variant: cells above 1 Hz in WAKE and REM, glmfit warnings removed, HD-HD: Pearson r = 0.698, p = 0, n = 26604 unordered pairs (53208 ordered); Spearman 0.786
variant: cells above 1 Hz in WAKE and REM, glmfit warnings removed, FS-FS: Pearson r = 0.380, p = 2.5e-91, n = 2639 unordered pairs (5278 ordered); Spearman 0.421
variant: cells above 1 Hz in WAKE and REM, glmfit warnings removed, HD-FS: Pearson r = 0.425, p = 0, n = 17006 unordered pairs (34012 ordered); Spearman 0.527
variant: cells above 1 Hz in WAKE and REM, round-1 rule (round-1 main), HD-HD: Pearson r = 0.698, p = 0, n = 26605 unordered pairs (53210 ordered); Spearman 0.786
variant: cells above 1 Hz in WAKE and REM, round-1 rule (round-1 main), FS-FS: Pearson r = 0.380, p = 2.5e-91, n = 2639 unordered pairs (5278 ordered); Spearman 0.421
variant: cells above 1 Hz in WAKE and REM, round-1 rule (round-1 main), HD-FS: Pearson r = 0.425, p = 0, n = 17006 unordered pairs (34012 ordered); Spearman 0.527
per session r (main), median and range: {'HD-HD': {'min': 0.408, '50%': 0.717, 'max': 0.837}, 'FS-FS': {'min': -0.257, '50%': 0.453, 'max': 0.732}, 'HD-FS': {'min': 0.145, '50%': 0.52, 'max': 0.824}}
# Warnings at other lags. GetGLM_Forward.m fits every pair at 201 lags (the target cell's smoothed
# counts circshifted by -100 to 100 bins of 0.1 s over the whole recording, then restricted to the
# state), and Figure_GLMRemWake.m removes a pair with a warning at any lag. Only lag 0 is fitted
# above. Here every third session of the 27 is also fitted at lags of -10, -1, 1 and 10 s (four of
# the 200 other lags, chosen for time), and the pairs with a warning at any of these lags are
# removed as well, to see how many more pairs the code's rule would remove and how the
# correlations of those sessions change.
t_lag = time.time()
lag_sess = sorted(GL.sess.unique())[::3]
LAGS = (-100, -10, 10, 100)
lag_warn = {lg: set() for lg in LAGS}
for si in lag_sess:
s = S[si]
remiv = s["sleep"]["rem"]
tmax = max(b for v in s["epochs"].values() for a, b in v)
edges = np.arange(0, tmax + 0.1, 0.1); cen = edges[:-1] + 0.05
Q = gaussian_filter1d(posub.bin_counts(s["spikes"], edges), 1, axis=0)
mua = Q.sum(1); mua = (mua - mua.mean()) / mua.std()
Qz = (Q - Q.mean(0)) / Q.std(0)
wa, wb = s["epochs"]["wake_square"][0]
states = {"WAKE": (cen >= wa) & (cen < wb), "REM": posub.in_intervals(cen, remiv)}
u = U[U.sess == si]
cells = np.where((u.hd | u.fs).values)[0]
for lg in LAGS:
Ql = np.roll(Q, lg, axis=0) # MATLAB's circshift(dQ, lag)
for st, msk in states.items():
Qm, mm, Zm = Ql[msk], mua[msk], Qz[msk]
def fit_lag(a):
others = np.array([c for c in cells if c != a])
b, bad, lim = posub.glmfit_poisson_batch(Qm[:, a], mm, Zm[:, others])
w = bad | lim | ~np.isfinite(b).all(0)
return {(si, int(min(a, o)), int(max(a, o))) for o in others[w]}
with ThreadPoolExecutor(4) as ex:
for w in ex.map(fit_lag, cells):
lag_warn[lg] |= w
print(f"fitted at {len(LAGS)} more lags in {len(lag_sess)} sessions ({[S[i]['subject'] for i in lag_sess]}) in {time.time() - t_lag:.0f} s")
GLs = GL[GL.sess.isin(lag_sess)].copy()
key = list(zip(GLs.sess, GLs.ci, GLs.cj))
for lg in LAGS:
GLs[f"warn_{lg}"] = [k in lag_warn[lg] for k in key]
GLs["warn_any"] = GLs.warn | GLs[[f"warn_{lg}" for lg in LAGS]].any(axis=1)
for p_ in ("HD-HD", "FS-FS", "HD-FS"):
d = GLs[GLs.pair == p_]
print(f"{p_}: {len(d)} pairs in these sessions; warning at lag 0 {int(d.warn.sum())}; "
+ "; ".join(f"at {lg / 10:+.0f} s {int(d[f'warn_{lg}'].sum())}" for lg in LAGS)
+ f"; at any of the five {int(d.warn_any.sum())}, of which not at lag 0 {int((d.warn_any & ~d.warn).sum())}")
print("WAKE against REM coupling in these sessions:")
corr_table(GLs[~GLs.warn], "warnings at lag 0 removed (as the main result)")
corr_table(GLs[~GLs.warn_any], "warnings at any of the five lags removed")
fitted at 4 more lags in 9 sessions (['A1813', 'A3701', 'A3705', 'A3709', 'A3712', 'A3717', 'A5505', 'A5508', 'A6215']) in 283 s HD-HD: 10507 pairs in these sessions; warning at lag 0 417; at -10 s 68; at -1 s 116; at +1 s 108; at +10 s 73; at any of the five 477, of which not at lag 0 60 FS-FS: 1017 pairs in these sessions; warning at lag 0 0; at -10 s 0; at -1 s 0; at +1 s 0; at +10 s 0; at any of the five 0, of which not at lag 0 0 HD-FS: 6528 pairs in these sessions; warning at lag 0 6; at -10 s 6; at -1 s 6; at +1 s 6; at +10 s 5; at any of the five 6, of which not at lag 0 0 WAKE against REM coupling in these sessions: warnings at lag 0 removed (as the main result), HD-HD: Pearson r = 0.638, p = 0, n = 10090 unordered pairs (20180 ordered); Spearman 0.747 warnings at lag 0 removed (as the main result), FS-FS: Pearson r = 0.282, p = 4.3e-20, n = 1017 unordered pairs (2034 ordered); Spearman 0.383 warnings at lag 0 removed (as the main result), HD-FS: Pearson r = 0.479, p = 0, n = 6522 unordered pairs (13044 ordered); Spearman 0.512 warnings at any of the five lags removed, HD-HD: Pearson r = 0.647, p = 0, n = 10030 unordered pairs (20060 ordered); Spearman 0.748 warnings at any of the five lags removed, FS-FS: Pearson r = 0.282, p = 4.3e-20, n = 1017 unordered pairs (2034 ordered); Spearman 0.383 warnings at any of the five lags removed, HD-FS: Pearson r = 0.479, p = 0, n = 6522 unordered pairs (13044 ordered); Spearman 0.512
Result: reproduced (resting on an assumption). With the code's cell selection (HD and FS cells, no rate condition), 27 sessions give 34,611 HD-HD, 2,683 FS-FS and 19,602 HD-FS pairs after the removal of fits with glmfit's warnings (paper 34,270, 2,492 and 18,992 in Extended Data Fig. 9d). WAKE and REM coupling coefficients are correlated for HD-HD pairs (r = 0.641), FS-FS pairs (r = 0.364, p = 1e-84) and HD-FS pairs (r = 0.359), all positive with p < 0.05, and FS-FS is below HD-HD; the Spearman correlations are 0.762, 0.417 and 0.514. The glmfit mirror flags 1,795 HD-HD pairs (4.9%), 24 HD-FS and no FS-FS pairs; the round-1 rule flags the same pairs (agreement 1.0). The flagged HD-HD fits are divergent (median |coupling| 2.77 in WAKE, up to 3e73, against 0.28 in kept pairs), and with them kept the pooled HD-HD Pearson r is 0.000 while the Spearman r is 0.759, so the outcome rests on the assumption that the mirror flags the pairs glmfit would. The round-1 selection of cells above 1 Hz in WAKE and REM gives r = 0.698, 0.380 and 0.425. Warnings at other lags (the code removes a pair with a warning at any of 201 lags) were checked in 9 of the 27 sessions at -10, -1, 1 and 10 s: they flag 60 more HD-HD pairs of 10,507 (0.6%; 417 at lag 0) and no more FS-FS or HD-FS pairs, and removing them moves the HD-HD r of those sessions from 0.638 to 0.647, FS-FS and HD-FS unchanged (0.282 and 0.479).
ed9d-coupling-follows-tuning¶
The polarity of the GLM coupling within a pair depends on the similarity of the two cells' HD tuning, during both WAKE and REM.
Pearson correlation across pairs between the zero-lag correlation of the two cells' square tuning curves and the lag-0 coupling, for HD-FS and FS-FS pairs, in WAKE and in REM (criterion; code Figure_GLMRemWake).
for p_ in ("HD-FS", "FS-FS"):
d = GLc[GLc.pair == p_]
for st in ("wake", "rem"):
r, p = stats.pearsonr(d.tc_r, d[st])
print(f"{p_}, {st.upper()}: r = {r:.3f}, p = {p:.3g}, n = {len(d)}")
d = GLc[GLc.pair == "HD-HD"]
print("HD-HD, shown only:", [f"{st}: r = {stats.pearsonr(d.tc_r, d[st])[0]:.3f}" for st in ("wake", "rem")])
HD-FS, WAKE: r = 0.758, p = 0, n = 19602 HD-FS, REM: r = 0.281, p = 0, n = 19602 FS-FS, WAKE: r = 0.451, p = 2.55e-134, n = 2683 FS-FS, REM: r = 0.123, p = 1.69e-10, n = 2683 HD-HD, shown only: ['wake: r = 0.669', 'rem: r = 0.503']
Result: reproduced. Coupling follows tuning similarity for HD-FS pairs (WAKE r = 0.758, REM r = 0.281, n = 19,602) and FS-FS pairs (WAKE r = 0.451, REM r = 0.123, n = 2,683), all p < 1e-9.
fig7b-wake-isomap-tuning¶
HD tuning curves of PoSub-HD and PoSub-FS cells computed from the virtual HD of the WAKE ring manifold are equivalent to those computed from the real HD.
Methods and GetIsomap_WAKE.m: HD cell spikes over the whole square epoch in 200 ms bins, smoothed with 400 ms s.d., square root; Isomap with 12 neighbours and 3 dimensions; virtual HD from the first two dimensions, its direction by the smaller circular variance of its difference from real HD, its origin by the mean offset of maximum cross-correlation between real and virtual curves over HD cells. Curves in 1 degree bins, 3 degree smoothing, from all times of the epoch (code); real HD taken at each bin centre. Control: the virtual HD reversed in time. Isomap is fitted on at most 4,000 vectors per session (helper; assumption, for time).
ISO = {}
iso_rows = []
for si, s in enumerate(S):
u = U[U.sess == si]
a, b = s["epochs"]["wake_square"][0]
edges = np.arange(a, b, 0.2); cen = edges[:-1] + 0.1
t = s["hd_t"]
j = np.clip(np.searchsorted(t, cen), 1, len(t) - 1); j = np.where(np.abs(t[j - 1] - cen) < np.abs(t[j] - cen), j - 1, j)
real = np.where(np.abs(t[j] - cen) < 0.05, s["hd"][j], np.nan)
hd_idx = np.where(u.hd.values)[0]
X = posub.pop_vectors([s["spikes"][i] for i in hd_idx], edges)
emb = posub.isomap_embed(X, 3, seed=700 + si)
v = posub.virtual_hd(emb, real)
allm = np.ones(len(cen), bool)
h_real = posub.binned_tuning(s["spikes"], edges, real, allm)
h_vir = posub.binned_tuning(s["spikes"], edges, v, allm)
lag = posub.circ_xcorr_lag(h_real[:, hd_idx], h_vir[:, hd_idx])
off = int(np.round(np.rad2deg(stats.circmean(np.deg2rad(lag))))) % 360
v = v - np.deg2rad(off) # the virtual curve rolled back by off bins to match the real one
h_vir = posub.binned_tuning(s["spikes"], edges, v, allm)
h_rev = posub.binned_tuning(s["spikes"], edges, v[::-1].copy(), allm)
r = posub.zero_lag_r(h_real, h_vir); rc = posub.zero_lag_r(h_real, h_rev)
rad = np.linalg.norm(emb[:, :2] - emb[:, :2].mean(0), axis=1)
ISO[si] = dict(emb=emb, v=v, real=real, edges=edges, rad=rad)
for i in np.where((u.hd | u.fs).values)[0]:
iso_rows.append(dict(subject=s["subject"], cls="HD" if u.hd.values[i] else "FS", r=r[i], r_rev=rc[i]))
IR = pd.DataFrame(iso_rows)
for cls in ("HD", "FS"):
d = IR[IR.cls == cls]
show(signrank(d.r, d.r_rev), f"{cls}: real vs virtual curve correlation, against reversed control:")
print("paper: HD n = 1602, z = 34.7; FS n = 427, z = 17.9")
print("per session median r (HD, FS):", IR.groupby(["subject", "cls"]).r.median().unstack().round(2).T.to_string())
print(f"Isomap done at {time.time() - T0:.0f} s")
HD: real vs virtual curve correlation, against reversed control: n=1512, median_a=0.9438, median_b=-0.001836, z=-33.68, p=1.127e-248 FS: real vs virtual curve correlation, against reversed control: n=430, median_a=0.842, median_b=0.01859, z=-17.95, p=4.437e-72 paper: HD n = 1602, z = 34.7; FS n = 427, z = 17.9 per session median r (HD, FS): subject A1808 A1813 A1815 A1821 A1824 A3701 A3702 A3703 A3705 A3706 A3707 A3708 A3709 A3710 A3711 A3712 A3713 A3716 A3717 A3723 A3728 A3730 A5505 A5506 A5507 A5508 A5801 A6211 A6215 A6216 A9701 cls FS 0.84 0.62 0.77 0.80 0.83 0.89 0.85 0.87 0.88 0.91 0.83 0.86 0.91 0.86 0.66 0.83 0.81 0.93 0.83 0.84 0.75 0.88 0.95 0.69 0.90 0.83 0.88 0.79 0.84 0.68 0.73 HD 0.95 0.91 0.90 0.94 0.94 0.95 0.89 0.95 0.97 0.96 0.96 0.97 0.98 0.96 0.77 0.89 0.93 0.97 0.97 0.88 0.91 0.89 0.95 0.84 0.97 0.87 0.97 0.92 0.96 0.91 0.92 Isomap done at 1368 s
Result: reproduced. Tuning curves from the virtual HD of the WAKE manifold match those from the real HD for HD cells (median r 0.94 against -0.002 for the reversed control, |z| = 33.7, p = 1e-248, n = 1,512) and FS cells (0.84 against 0.02, |z| = 17.95, p = 4e-72, n = 430).
ed10b-ring-topology¶
Population vectors of PoSub-HD cells lie on a ring during WAKE, farther from the centre of the manifold than shuffled vectors, and the ring has the same radius during REM.
Distance of each point to the centre of the first two Isomap dimensions, normalized to the mean distance of the real projection. Shuffle: each cell's binned train shifted circularly by a random number of bins, then embedded the same way. Real against shuffled mean distance, Wilcoxon signed-rank across sessions. REM against WAKE: the joint embedding below (WAKE vectors subsampled to the number of REM vectors), mean distance per session in each state, Wilcoxon signed-rank (a claim of absence, judged at p < 0.05).
rng = np.random.default_rng(10)
ring = []
for si, s in enumerate(S):
u = U[U.sess == si]; hd_idx = np.where(u.hd.values)[0]
e = ISO[si]["edges"]
C = posub.bin_counts([s["spikes"][i] for i in hd_idx], e)
C = np.stack([np.roll(C[:, k], rng.integers(len(C))) for k in range(C.shape[1])], 1)
Xs = np.sqrt(gaussian_filter1d(C, 2, axis=0))
es = posub.isomap_embed(Xs, 3, seed=900 + si)
rr = ISO[si]["rad"]; rs = np.linalg.norm(es[:, :2] - es[:, :2].mean(0), axis=1)
ring.append(dict(subject=s["subject"], real=1.0, shuffled=rs.mean() / rr.mean(), real_cv=rr.std() / rr.mean(), shuffled_cv=rs.std() / rs.mean()))
ring = pd.DataFrame(ring)
print(ring.round(3).to_string())
show(signrank(ring.real, ring.shuffled), "mean normalized distance, real (1 by construction) vs shuffled:")
print("paper: n = 32, z = 4.79")
subject real shuffled real_cv shuffled_cv 0 A1808 1.0 0.607 0.226 0.568 1 A1813 1.0 0.623 0.333 0.469 2 A1815 1.0 0.618 0.241 0.490 3 A1821 1.0 0.384 0.158 0.418 4 A1824 1.0 0.490 0.231 0.493 5 A3701 1.0 0.350 0.188 0.479 6 A3702 1.0 0.380 0.185 0.567 7 A3703 1.0 0.421 0.201 0.465 8 A3705 1.0 0.403 0.230 0.459 9 A3706 1.0 0.349 0.154 0.523 10 A3707 1.0 0.652 0.314 0.540 11 A3708 1.0 0.403 0.219 0.523 12 A3709 1.0 0.537 0.300 0.462 13 A3710 1.0 0.330 0.130 0.513 14 A3711 1.0 0.410 0.318 0.501 15 A3712 1.0 0.408 0.209 0.463 16 A3713 1.0 0.490 0.200 0.559 17 A3716 1.0 0.408 0.171 0.492 18 A3717 1.0 0.393 0.192 0.534 19 A3723 1.0 0.442 0.237 0.542 20 A3728 1.0 0.449 0.215 0.580 21 A3730 1.0 0.418 0.204 0.484 22 A5505 1.0 0.511 0.205 0.481 23 A5506 1.0 0.504 0.280 0.616 24 A5507 1.0 0.386 0.191 0.524 25 A5508 1.0 0.460 0.244 0.487 26 A5801 1.0 0.383 0.187 0.611 27 A6211 1.0 0.450 0.177 0.581 28 A6215 1.0 0.466 0.181 0.498 29 A6216 1.0 0.465 0.224 0.509 30 A9701 1.0 0.461 0.167 0.580 mean normalized distance, real (1 by construction) vs shuffled: n=31, median_a=1, median_b=0.4416, z=-4.86, p=1.174e-06 paper: n = 32, z = 4.79
Result: reproduced. The real population vectors lie farther from the manifold centre than shuffled ones in every session (shuffled mean distance 0.33 to 0.65 of the real, signed-rank p = 1.2e-6, n = 31; paper z = 4.79), and in the joint WAKE-REM embedding the mean REM distance does not differ from WAKE (median ratio 0.98, p = 0.46, n = 27; paper z = 0.17).
fig7d-rem-isomap-tuning¶
During REM sleep the HD tuning of PoSub-HD and PoSub-FS cells can be estimated from the virtual HD of the ring manifold alone and matches their WAKE tuning.
Sessions with at least 100 s of REM (Methods). HD-cell population vectors (200 ms, 400 ms s.d., square root) in REM and in the square epoch; WAKE vectors subsampled at random to the number of REM vectors; Isomap (12 neighbours, 3 dimensions) on both together; virtual HD set in direction and origin from the WAKE vectors' real HD. Curves in 6 degree bins with 6 degree smoothing. Cells with REM rate above 0.5 Hz (HD) or 10 Hz (FS). WAKE curve read (i) as the WAKE Isomap curve and (ii) as the real-HD WAKE curve; control, the REM curve from the REM virtual HD reversed in time.
rng = np.random.default_rng(7)
rem_rows, rad_rows = [], []
for si, s in enumerate(S):
remiv = s["sleep"].get("rem", np.zeros((0, 2)))
if (remiv[:, 1] - remiv[:, 0]).sum() < 100:
continue
u = U[U.sess == si]; hd_idx = np.where(u.hd.values)[0]
tmax = max(b for v in s["epochs"].values() for a, b in v)
edges = np.arange(0, tmax, 0.2); cen = edges[:-1] + 0.1
in_rem = posub.in_intervals(cen, remiv)
a, b = s["epochs"]["wake_square"][0]
in_wk = (cen >= a) & (cen < b)
wk_idx = np.sort(rng.choice(np.where(in_wk)[0], in_rem.sum(), replace=False))
rem_idx = np.where(in_rem)[0]
X = posub.pop_vectors([s["spikes"][i] for i in hd_idx], edges)
sel = np.r_[wk_idx, rem_idx]
emb = posub.isomap_embed(X[sel], 3, seed=1100 + si)
t = s["hd_t"]
j = np.clip(np.searchsorted(t, cen[sel]), 1, len(t) - 1); j = np.where(np.abs(t[j - 1] - cen[sel]) < np.abs(t[j] - cen[sel]), j - 1, j)
real = np.where(np.abs(t[j] - cen[sel]) < 0.05, s["hd"][j], np.nan)
nw = len(wk_idx); real[nw:] = np.nan
v = posub.virtual_hd(emb, real) # direction from the WAKE vectors (REM real HD is NaN)
ang_full = np.full(len(cen), np.nan); ang_full[sel] = v
real_full = np.full(len(cen), np.nan); real_full[sel[:nw]] = real[:nw]
mw = np.zeros(len(cen), bool); mw[wk_idx] = True
mr = np.zeros(len(cen), bool); mr[rem_idx] = True
kw = dict(nbins=60, sd=1)
hw_real = posub.binned_tuning(s["spikes"], edges, real_full, mw, **kw)
hw_vir = posub.binned_tuning(s["spikes"], edges, ang_full, mw, **kw)
lag = posub.circ_xcorr_lag(hw_real[:, hd_idx], hw_vir[:, hd_idx])
off = stats.circmean(np.deg2rad(lag * 6))
ang_full = ang_full - off
hw_vir = posub.binned_tuning(s["spikes"], edges, ang_full, mw, **kw)
hr_vir = posub.binned_tuning(s["spikes"], edges, ang_full, mr, **kw)
rev = ang_full.copy(); rev[rem_idx] = ang_full[rem_idx][::-1]
hr_rev = posub.binned_tuning(s["spikes"], edges, rev, mr, **kw)
rem_dur = (remiv[:, 1] - remiv[:, 0]).sum()
rate_rem = np.array([posub.in_intervals(sp, remiv).sum() / rem_dur for sp in s["spikes"]])
ri = posub.zero_lag_r(hw_vir, hr_vir); rii = posub.zero_lag_r(hw_real, hr_vir)
ci = posub.zero_lag_r(hw_vir, hr_rev); cii = posub.zero_lag_r(hw_real, hr_rev)
for i in range(len(u)):
cls = "HD" if u.hd.values[i] and rate_rem[i] > 0.5 else ("FS" if u.fs.values[i] and rate_rem[i] > 10 else None)
if cls:
rem_rows.append(dict(subject=s["subject"], cls=cls, r_iso=ri[i], r_real=rii[i], c_iso=ci[i], c_real=cii[i]))
rad = np.linalg.norm(emb[:, :2] - emb[:, :2].mean(0), axis=1); rad = rad / rad[:nw].mean()
wv = posub.zero_lag_r(hw_real[:, hd_idx], hw_vir[:, hd_idx])
rad_rows.append(dict(subject=s["subject"], rem_s=round(rem_dur), n_vectors=int(in_rem.sum()), wake=rad[:nw].mean(), rem=rad[nw:].mean(),
wake_real_vs_virtual_r=np.nanmedian(wv)))
RR = pd.DataFrame(rem_rows); RD = pd.DataFrame(rad_rows)
print(RD.round(3).to_string())
print("sessions with at least 100 s of REM:", len(RD), "(paper 26 of 32)")
for cls in ("HD", "FS"):
d = RR[RR.cls == cls]
show(signrank(d.r_iso, d.c_iso), f"{cls}, reading (i) WAKE Isomap curve vs REM, against control:")
show(signrank(d.r_real, d.c_real), f"{cls}, reading (ii) WAKE real-HD curve vs REM, against control:")
print("paper: HD n = 1148, z = 29.2; FS n = 317, z = 10.4")
print("ed10b, REM vs WAKE radius in the joint embedding:")
show(signrank(RD.rem, RD.wake), " mean normalized distance, REM vs WAKE:")
print("paper: n = 26, z = 0.17")
/home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:376: RuntimeWarning: invalid value encountered in divide return (a * b).sum(0) / np.sqrt((a * a).sum(0) * (b * b).sum(0))
/home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:376: RuntimeWarning: invalid value encountered in divide return (a * b).sum(0) / np.sqrt((a * a).sum(0) * (b * b).sum(0))
/home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:376: RuntimeWarning: invalid value encountered in divide return (a * b).sum(0) / np.sqrt((a * a).sum(0) * (b * b).sum(0))
/home/runner/work/encore/encore/entries/000939-duszkiewicz-2024/posub.py:376: RuntimeWarning: invalid value encountered in divide return (a * b).sum(0) / np.sqrt((a * a).sum(0) * (b * b).sum(0))
subject rem_s n_vectors wake rem wake_real_vs_virtual_r 0 A1813 389 1945 1.0 1.053 0.770 1 A1815 285 1425 1.0 1.275 0.628 2 A1821 267 1335 1.0 0.796 0.924 3 A3701 308 1538 1.0 0.919 0.906 4 A3702 124 620 1.0 1.032 0.928 5 A3703 301 1505 1.0 1.032 0.903 6 A3705 620 3100 1.0 0.931 0.976 7 A3706 218 1090 1.0 0.946 0.971 8 A3708 501 2505 1.0 0.849 0.975 9 A3709 978 4890 1.0 0.985 0.973 10 A3710 288 1440 1.0 0.905 0.951 11 A3711 293 1465 1.0 1.182 0.836 12 A3712 236 1180 1.0 0.771 0.960 13 A3713 200 1000 1.0 0.762 0.902 14 A3716 284 1420 1.0 1.133 0.918 15 A3717 522 2610 1.0 0.771 0.890 16 A3723 257 1285 1.0 0.841 0.806 17 A3730 228 1140 1.0 0.939 0.834 18 A5505 355 1775 1.0 1.002 0.875 19 A5506 340 1700 1.0 1.039 0.906 20 A5507 380 1900 1.0 1.038 0.932 21 A5508 122 610 1.0 1.107 0.902 22 A5801 139 695 1.0 1.252 0.919 23 A6211 214 1070 1.0 1.228 0.869 24 A6215 248 1240 1.0 1.045 0.904 25 A6216 429 2145 1.0 0.937 0.892 26 A9701 313 1565 1.0 0.754 0.952 sessions with at least 100 s of REM: 27 (paper 26 of 32) HD, reading (i) WAKE Isomap curve vs REM, against control: n=1259, median_a=0.8562, median_b=-0.09993, z=-30.4, p=6.218e-203 HD, reading (ii) WAKE real-HD curve vs REM, against control: n=1259, median_a=0.8199, median_b=-0.1102, z=-30.31, p=9.392e-202 FS, reading (i) WAKE Isomap curve vs REM, against control: n=336, median_a=0.4265, median_b=-0.003542, z=-11.36, p=6.942e-30 FS, reading (ii) WAKE real-HD curve vs REM, against control: n=336, median_a=0.361, median_b=-0.02412, z=-11.45, p=2.443e-30 paper: HD n = 1148, z = 29.2; FS n = 317, z = 10.4 ed10b, REM vs WAKE radius in the joint embedding: mean normalized distance, REM vs WAKE: n=27, median_a=0.9849, median_b=1, z=-0.7448, p=0.4564 paper: n = 26, z = 0.17
Result: reproduced. In 27 sessions with at least 100 s of REM inside the home-cage epochs, REM tuning curves from the virtual HD match WAKE curves above the time-reversed control for HD cells (median r 0.86, |z| = 30.4, n = 1,259) and FS cells (0.43, |z| = 11.4, n = 336), with the WAKE curve from the Isomap angle, and the same with the real-HD WAKE curve (0.82 and 0.36, |z| = 30.3 and 11.5); all p < 1e-29.
fig1ac-posub-border (figshare, data elsewhere)¶
Along a Neuropixel probe crossing pRSC and PoSub, HD tuning of excitatory cells rises in a step at the functional border between the two regions.
The Neuropixel session is figshare Dataset_4.tar.gz (claims file), streamed here. Excitatory cells as for PoSub (rate at least 0.5 Hz in the square epoch, trough-to-peak above 0.4 ms from the file's WaveformFeatures, rate below 10 Hz); HD information from the square-epoch curve at speed above 2 cm/s (file's Velocity). Depth from CellDepth.mat (micrometres), set to 0 at the shallowest cell as the code sets its depth. Running average as Figure_SingleRec_NPX.m: mean HD information in 50 um bins from 0 to -2.5 mm, empty bins interpolated, Gaussian smoothing of 2 bins. The judgement is qualitative (criterion).
# NP, figshare Dataset_4, was streamed in the cell-classes-counts section
tunN = posub.adn_tuning(NP)
ttpN = NP["mat"]["tr2pk"].ravel(); depN = NP["mat"]["cellDep"].ravel().astype(float)
depN = (depN - depN.max()) / 1000
exN = (tunN["rate"] >= 0.5) & np.isfinite(ttpN) & (ttpN > 0.4) & (tunN["rate"] < 10)
print("Neuropixel session", NP["session"], "; units", len(NP["spikes"]), "; excitatory", int(exN.sum()),
"; depth range (mm)", depN.min(), depN.max(), "; the files' own excitatory count", int(NP["file_types"]["ex"].sum()))
top = -0.05 * np.arange(50); br = np.c_[top, top - 0.05]
prof = np.array([tunN["info"][exN & (depN < a) & (depN >= b)].mean() if (exN & (depN < a) & (depN >= b)).any() else np.nan for a, b in br])
ok = np.isfinite(prof); tt = np.arange(len(prof))
prof_i = prof.copy(); prof_i[~ok] = np.interp(tt[~ok], tt[ok], prof[ok])
prof_s = gaussian_filter1d(prof_i, 2)
mid = br.mean(1)
print("running mean HD information (bits/spike) in every 50 um bin, by bin centre (mm):")
for k0 in range(0, len(mid), 10):
print(" " + " ".join(f"{d_:6.3f}: {v_:.3f}" for d_, v_ in zip(mid[k0:k0 + 10], prof_s[k0:k0 + 10])))
for lab, m_ in (("top 1.4 mm (centres above -1.4)", mid > -1.4), ("from 1.8 mm down (centres below -1.8)", mid < -1.8)):
print(f"{lab}: running mean {prof_s[m_].min():.3f} to {prof_s[m_].max():.3f}")
loc = np.where((prof_s[1:-1] > prof_s[:-2]) & (prof_s[1:-1] > prof_s[2:]))[0] + 1
print("local maxima of the running mean at (mm, bits/spike):", [(round(float(mid[k]), 3), round(float(prof_s[k]), 3)) for k in loc])
dpr = np.diff(prof_s)
k = int(np.argmax(dpr if prof_s[-1] > prof_s[0] else -dpr))
print(f"largest change between adjacent bins at {mid[k]:.3f} mm ({dpr[k]:+.3f}); profile range {prof_s.min():.3f} to {prof_s.max():.3f}; "
f"shallow third mean {prof_s[: len(prof_s) // 3].mean():.3f}, deep third mean {prof_s[-len(prof_s) // 3:].mean():.3f}")
half = (prof_s.max() + prof_s.min()) / 2
cross = np.where(np.diff(np.sign(prof_s - half)) != 0)[0]
print("crossings of the half-range level at depths (mm):", mid[cross].round(3))
fig, ax = plt.subplots(1, 2, figsize=(8, 4))
ax[0].scatter(tunN["info"][exN], depN[exN], s=6); ax[0].set_xlabel("HD info (bits/spike)"); ax[0].set_ylabel("depth (mm)")
ax[1].plot(prof_s, mid); ax[1].set_xlabel("running mean HD info"); plt.tight_layout(); plt.show()
print("plotted ranges: info", tunN["info"][exN].min().round(3), tunN["info"][exN].max().round(3))
# Which end is which. The files give the depth (CellDepth.mat, cellDep, um) and no probe
# geometry. Figure_SingleRec_NPX.m sets depth = depth - max(depth), so the largest value is the
# top of the profile (0 mm), as here. GetCellTypes_PoSub.m carries, commented out, the line
# "gd(cellDep > 0) = 0; % For the neuropixel rec", which drops the units above cellDep = 0 from
# the PoSub cells: the code's zero is the PoSub border, with PoSub at negative values.
rawN = NP["mat"]["cellDep"].ravel().astype(float)
print(f"cellDep (um): range {rawN.min():.0f} to {rawN.max():.0f}; units above 0: {(rawN > 0).sum()} of {len(rawN)}")
print(f"half-range crossing(s) in the file's units (um): {(mid[cross] * 1000 + rawN.max()).round(0)}")
for lab, m_ in (("cellDep above 0 (pRSC side by the code's line)", rawN > 0), ("cellDep at or below 0 (PoSub side)", rawN <= 0)):
x = tunN["info"][exN & m_]
print(f"{lab}: {len(x)} excitatory cells, mean HD information {x.mean():.3f}, median {np.median(x):.3f} bits/spike; "
f"above 0.2 bits/spike {np.mean(x > 0.2):.2f}")
Neuropixel session A8603-210602 ; units 239 ; excitatory 134 ; depth range (mm) -2.32 0.0 ; the files' own excitatory count 137 running mean HD information (bits/spike) in every 50 um bin, by bin centre (mm): -0.025: 0.008 -0.075: 0.008 -0.125: 0.008 -0.175: 0.008 -0.225: 0.008 -0.275: 0.009 -0.325: 0.014 -0.375: 0.026 -0.425: 0.050 -0.475: 0.088 -0.525: 0.134 -0.575: 0.174 -0.625: 0.191 -0.675: 0.182 -0.725: 0.155 -0.775: 0.127 -0.825: 0.110 -0.875: 0.106 -0.925: 0.107 -0.975: 0.105 -1.025: 0.098 -1.075: 0.092 -1.125: 0.097 -1.175: 0.117 -1.225: 0.152 -1.275: 0.192 -1.325: 0.226 -1.375: 0.240 -1.425: 0.237 -1.475: 0.238 -1.525: 0.270 -1.575: 0.349 -1.625: 0.463 -1.675: 0.590 -1.725: 0.715 -1.775: 0.828 -1.825: 0.909 -1.875: 0.934 -1.925: 0.896 -1.975: 0.827 -2.025: 0.777 -2.075: 0.788 -2.125: 0.859 -2.175: 0.959 -2.225: 1.052 -2.275: 1.119 -2.325: 1.159 -2.375: 1.180 -2.425: 1.189 -2.475: 1.192 top 1.4 mm (centres above -1.4): running mean 0.008 to 0.240 from 1.8 mm down (centres below -1.8): running mean 0.777 to 1.192 local maxima of the running mean at (mm, bits/spike): [(-0.625, 0.191), (-0.925, 0.107), (-1.375, 0.24), (-1.875, 0.934)] largest change between adjacent bins at -1.625 mm (+0.127); profile range 0.008 to 1.192; shallow third mean 0.074, deep third mean 0.940 crossings of the half-range level at depths (mm): [-1.675]
plotted ranges: info 0.001 2.248 cellDep (um): range -680 to 1640; units above 0: 145 of 239 half-range crossing(s) in the file's units (um): [-35.] cellDep above 0 (pRSC side by the code's line): 71 excitatory cells, mean HD information 0.145, median 0.080 bits/spike; above 0.2 bits/spike 0.14 cellDep at or below 0 (PoSub side): 63 excitatory cells, mean HD information 0.905, median 0.757 bits/spike; above 0.2 bits/spike 0.89
Result: reproduced (qualitative). Along the figshare Neuropixel probe the running mean HD information of 134 excitatory cells lies between 0.008 and 0.240 bits per spike over the top 1.4 mm (bin centres above -1.4 mm) and between 0.777 and 1.192 from 1.8 mm down, crossing its half-range level once, at -1.675 mm (top third mean 0.074, bottom third 0.940); the profile has a smaller local maximum at -0.625 mm (0.191), and two more at -0.925 mm (0.107) and -1.375 mm (0.240). The crossing lies at cellDep = -35 um in the file's units, next to the code's zero: the commented-out line of GetCellTypes_PoSub.m drops units above cellDep = 0 from the PoSub cells for this recording, so the low end is the pRSC side (71 excitatory cells, mean 0.145 bits per spike) and the high end the PoSub side (63 cells, mean 0.905), as the criterion requires. The files hold no probe geometry, so which end is shallower rests on the code's depth - max(depth), not on the data.
Claims not tested¶
- fig4-random-connectivity-simulation. A result of simulation. The paper's Python simulation code is not in the figshare item, and the Methods and Supplementary Methods do not give the input tuning widths, the number of inputs or the log-normal weight parameters of each point.
- fig4h-model-fit. The PyTorch model fit is not deposited, and the Supplementary Methods do not give the numbers of inputs or the initial values.
- ed6-trn-innervates-adn. Anatomical tracing images, which are neither in the dandiset nor
in the figshare item nor anywhere the paper names; recorded as blocked for
other, since the claims file marks the claim as data elsewhere and the protocol does not allowdata_missingfor such claims.
Summary¶
| Claim | Data | Outcome |
|---|---|---|
| cell-classes-counts | dandiset | reproduced (parallel units short by the one missing session, the Neuropixel recording) |
| ed2b-excitatory-hd-information | dandiset | reproduced |
| fig2a-fs-hd-information | dandiset | reproduced |
| fig2b-decoding-from-fs | dandiset | reproduced |
| fig2c-geometry-independent | dandiset | reproduced |
| ed2hi-within-epoch-stability | dandiset | reproduced |
| fig2f-cue-rotation | figshare Dataset_3 | reproduced |
| fig3c-adn-broader | figshare Dataset_5 | reproduced |
| fig3-ten-components-power | dandiset | partially reproduced (holds for squared-modulus power only) |
| fig3e-fourier-signatures-differ | figshare Dataset_5 | partially reproduced (direction holds; interaction significant under some readings only) |
| fig3f-posub-hd-closer-to-fs | figshare Dataset_5 | partially reproduced (p < 1e-9 under one of four readings) |
| ed3k-robust-to-controls | figshare Dataset_5 | reproduced (claim amended after a reader's answer) |
| fig3g-first-three-components | dandiset | partially reproduced (components 1-2 or 1-4 above control, not 1-3) |
| fig3h-dominant-component | dandiset | not reproduced (third count 47, outside 55 +- 6) |
| ed3g-fs-spectra-heterogeneous | dandiset | reproduced |
| ed3l-fs-phases-independent | dandiset | reproduced |
| ed4a-autocorrelogram-symmetries | dandiset | reproduced |
| ed4cd-isomap-posub-central | figshare Dataset_5 | reproduced |
| fig4-random-connectivity-simulation | none deposited | blocked |
| fig4f-synaptic-weights-dispersed | dandiset | reproduced |
| fig4h-model-fit | none deposited | blocked |
| ed7-adn-disinhibition | figshare Dataset_6 | reproduced |
| fig5b-posub-rates-increase | dandiset | reproduced |
| ed7gj-tuning-preserved | dandiset | reproduced |
| fig5e-hd-multiplicative-fs-additive | dandiset | partially reproduced (one control comparison fails under one exclusion reading) |
| ed8cd-multiplicative-detectable | dandiset | reproduced (titrated noise and the simulated control, the one control that passes the no-gain check; holds in 80% of 50 seeds) |
| ed8e-adn-mixed-gain | figshare Dataset_6 | reproduced |
| fig6c-coupling-preserved-wake-rem | dandiset | reproduced (assumed: the glmfit mirror flags glmfit's warnings) |
| ed9d-coupling-follows-tuning | dandiset | reproduced |
| fig7b-wake-isomap-tuning | dandiset | reproduced |
| ed10b-ring-topology | dandiset | reproduced |
| fig7d-rem-isomap-tuning | dandiset | reproduced |
| fig1ac-posub-border | figshare Dataset_4 | reproduced (qualitative) |
| ed6-trn-innervates-adn | none deposited | blocked |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. Python reimplementations of the MATLAB code (tuning curves, HD information, Bayesian decoding, Poisson GLMs by a mirror of glmfit's iteratively reweighted least squares, CCG synapse detection); Isomap from scikit-learn, fitted on at most 4,000 population vectors per session; speed computed from the position, as the files hold no speed series; only lag 0 of the GLM, so pairs are removed for warnings at lag 0 only (four more lags in 9 sessions remove 0.6% more HD-HD pairs); REM restricted to the home-cage epochs; the dandiset's 31 sessions in place of the paper's 32, and 2 of the 5 ArchT mice for the PoSub optogenetic claims (the other three added from figshare in a variant); cue rotation realignment from a Methods-based reimplementation of the code's sigmoid procedure.
What was found about the dataset. The dandiset holds 31 of the paper's 32 PoSub sessions: all 14 vertical ones (931 units, as the paper gives) and 17 of 18 parallel ones (1,760 units against 1,999). The missing one is the Neuropixel session A8603-210602 (239 units, figshare Dataset_4), which the paper's parallel count includes. In three files (A3711, A3717, A3723) the sleep states run to twice the session's length, so that about half of their REM lies after the recording ends, and in three some REM lies inside an exploration epoch (A1821, 16 s in the square and 41 s in the optogenetic epoch; A3712, 12 s in the square; A3701, 0.3 s in the triangle); the REM intervals are the same in version 0.250207.0025. Of the five ArchT mice of the PoSub optogenetic experiment it holds two (A1813, A1815); the other three (A1801, A1804, A1806, six sessions) are on figshare only, as are the ADN, cue rotation and Neuropixel recordings. The optogenetic stimulus series of eight files holds four interleaved channels in one dimension, so its length at the stated 20 kHz is four times the session's; in A1808 it holds one channel of 0 and 1. The dandiset's unit labels (is_excitatory, is_fast_spiking, is_head_direction) follow the code's 0.35 ms cut, not the Methods' 0.4 ms. The files hold no speed series and the probe orientation is given only in the device description.