Reproducing Girardeau, Inema & Buzsaki (2017) from dandiset 000061¶
Paper. Girardeau G, Inema I, Buzsaki G (2017). Reactivations of emotional memory in the hippocampus-amygdala system during sleep. Nature Neuroscience 20, 1634-1642. doi:10.1038/nn.4637.
Claims. claims/000061/10.1038-nn.4637.yaml, written by a reader from the version of record
(main text, figure legends, Online Methods, reporting summary) and the Supplementary Information,
with the paper's code (github.com/buzsakilab/papers, GGirardeau-BLAHpcInteractions-Package, at
commit cdcec89f), and checked by a second reader. Its criteria fix every method; the shared
methods M1 to M8 and the open details O1 to O7 are in its opening comment. This attempt
(000061-girardeau-2017-02) replaces 000061-girardeau-2017-01, whose claims were written from
the abstract alone. Three questions this attempt asked (questions/000061.yaml) were answered by
a reader and the claims file updated: laps are judged on O1 reading c alone when the dataset has
neither lap intervals nor position; the sleep firing claims are judged on the identified cells (O4
reading a) alone; and the all-cell tests of reactivation-peaks-at-ripples,
direction-all-cell-pairs and ripple-modulated-cells-larger-reactivation no longer depend on
claims that need putative cell types, so they are judged here.
Dandiset. DANDI:000061, version 0.250624.0426. The paper cites no DANDI version (it says the data will be on CRCNS). There are two published versions, 0.210812.1517 and 0.250624.0426, with the same 40 files, and the draft (modified 2026-10-05) holds the same 40 assets as 0.250624.0426; the data check prints this. The most recent published version is pinned, as the protocol asks when the paper cites none: it is the version the earlier attempt tested, no later published version exists, and the draft differs from it only in metadata the analysis does not use (keywords, anatomy terms, ethics approval, affiliations).
What the dandiset holds. 40 NWB files, one per recording day, from four rats (Rat08 to Rat11):
sorted units with a region label, sleep states, epochs, and event times of hippocampal ripples,
air puffs and rewards; some files also hold raw voltage, LFP and video frames. There is no
position tracking, and the putative cell types the paper uses (pyramidal cell or interneuron,
from monosynaptic connections and then k-means on the waveform) are not in it: five Rat08 files
have a cell_type column that marks only cells identified as excitatory or inhibitory. Every
claim that needs putative pyramidal cells is therefore blocked; where a claim's criterion judges
its all-cell test alone when the pyramidal test cannot be run, that test decides it.
Random numbers are drawn only where power and bounds are given for the claims of absence
(bla-int-sleep-firing, cen-no-reactivation, ripple-rate-unchanged); each of those sections
has its own generator with a fixed seed, and none of these numbers decides an outcome.
| Claim | Tested here |
|---|---|
unit-counts |
tested |
rats-learn-puff-location |
blocked (too large: no position; only raw video frames) |
each-rat-learns-puff-location |
blocked (too large, depends on rats-learn-puff-location) |
dz-speed-slower-post-run |
blocked (too large, depends on rats-learn-puff-location) |
puff-location-no-speed-bias |
blocked (too large, depends on rats-learn-puff-location) |
habituation-over-days |
tested (laps under O1 reading c) |
bla-pyr-sleep-firing-increase |
tested (identified cells, five Rat08 files) |
bla-int-sleep-firing |
tested (identified cells, five Rat08 files) |
hpc-bla-reactivation-post-sleep |
blocked (metadata insufficient: no putative cell types) |
reactivation-all-cell-pairs |
tested |
reactivation-decay |
blocked (metadata insufficient) |
no-rem-reactivation |
blocked (metadata insufficient); the all-cell variant is shown |
piriform-weaker-reactivation |
blocked (metadata insufficient) |
cen-no-reactivation |
tested |
intra-structure-reactivation |
blocked (metadata insufficient); the all-cell tests are shown |
bla-ripple-up-subgroup |
blocked (metadata insufficient); fractions over all BLA units shown |
bla-ripple-down-subgroup |
blocked (metadata insufficient) |
ripple-modulated-cells-larger-reactivation |
tested (all-cell test alone; the pyramidal test cannot be run) |
run-correlated-pair-fractions |
blocked (metadata insufficient) |
reactivation-involves-ripple-up-cells |
blocked (metadata insufficient) |
pair-group-ev |
blocked (metadata insufficient) |
contributing-up-cells-gain-increase |
blocked (metadata insufficient) |
ev-independent-of-firing-rate |
blocked (metadata insufficient) |
reactivation-peaks-at-ripples |
tested (all-cell test alone; the pyramidal test cannot be run) |
airpuff-trajectory-ripple-reactivation |
blocked (metadata insufficient) |
safe-trajectory-no-ripple-reactivation |
blocked (metadata insufficient) |
aversive-direction-stronger |
blocked (metadata insufficient) |
direction-all-cell-pairs |
tested (laps under O1 reading c) |
bla-rate-same-across-directions |
blocked (too large: O7 reading b needs position); the runnable readings are shown |
ripple-rate-unchanged |
tested |
contributing-pairs-run-sleep-coupling |
blocked (metadata insufficient) |
contributing-quartiles-run-sleep-coupling |
blocked (metadata insufficient) |
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import hpcbla
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
notebook_started = started = time.time()
sessions = sorted(hpcbla.load_all(), key=lambda s: s["path"])
print(f"loaded {len(sessions)} files from dandiset {hpcbla.DANDISET_ID} version "
f"{hpcbla.DANDISET_VERSION} in {time.time() - started:.0f} s")
loaded 40 files from dandiset 000061 version 0.250624.0426 in 244 s
Data check¶
Each check of the metadata against the values is stated under its output. The status of each
finding of the earlier attempt (000061-girardeau-2017-01) is summarised at the end of this
section.
from dandi.dandiapi import DandiAPIClient
with DandiAPIClient() as client:
ds = client.get_dandiset(hpcbla.DANDISET_ID)
print("Versions:", [(v.identifier, str(v.modified)[:19]) for v in ds.get_versions()])
meta = client.get_dandiset(hpcbla.DANDISET_ID, hpcbla.DANDISET_VERSION).get_raw_metadata()
draft_meta = client.get_dandiset(hpcbla.DANDISET_ID, "draft").get_raw_metadata()
pinned_assets = {a.path: a.identifier for a in client.get_dandiset(hpcbla.DANDISET_ID, hpcbla.DANDISET_VERSION).get_assets()}
draft_assets = {a.path: a.identifier for a in client.get_dandiset(hpcbla.DANDISET_ID, "draft").get_assets()}
summary = meta["assetsSummary"]
print("Archive summary: files", summary["numberOfFiles"], "; subjects", summary["numberOfSubjects"],
"; species", [s["name"] for s in summary["species"]], "; measured", summary["variableMeasured"])
print("Draft against pinned version: assets", len(draft_assets), "and", len(pinned_assets),
"; same asset identifiers:", sum(pinned_assets.get(k) == v for k, v in draft_assets.items()))
skip = {"@context", "citation", "id", "url", "version", "doi", "datePublished", "publishedBy", "manifestLocation",
"schemaVersion", "dateCreated", "dateModified", "assetsSummary"}
print("Metadata fields that differ in the draft:", sorted(k for k in set(meta) | set(draft_meta)
if k not in skip and meta.get(k) != draft_meta.get(k)))
print("Related resources (pinned):", [(r.get("relation"), r.get("identifier")) for r in meta.get("relatedResource", [])])
rows = []
for s in sessions:
u = s["units"]
pre, post, run = hpcbla.sleep_epochs(s)
region = pd.Series(u["region"]).value_counts() if u else pd.Series(dtype=int)
rows.append({
"session": s["session"], "rat": s["rat"], "units": 0 if u is None else len(u["region"]),
"hpc": int(sum(region.get(k, 0) for k in hpcbla.HPC_LABELS)),
"bla": int(sum(region.get(k, 0) for k in hpcbla.BLA_LABELS)),
"cen": int(sum(region.get(k, 0) for k in hpcbla.CEN_LABELS)),
"pir": int(sum(region.get(k, 0) for k in hpcbla.PIR_LABELS)),
"cell_type": bool(u and "cell_type" in u["columns"]),
"states": s["states"] is not None,
"ripples": 0 if s["ripples"] is None else len(s["ripples"]["peak"]),
"puffs": len(s["events"]["puf"][1]) if "puf" in s["events"] else 0,
"rewards": len(s["events"]["lrw"][1]) + len(s["events"]["rrw"][1]),
"pre_sleep_s": round(float(np.diff(pre).sum())), "post_sleep_s": round(float(np.diff(post).sum())),
"run_s": None if run is None else round(run[1] - run[0]),
"lfp": s["lfp"] is not None, "raw": s["raw"] is not None, "videos": len(s["videos"]),
"trials": None if s["trials"] is None else len(s["trials"][0])})
info = pd.DataFrame(rows)
print("\nSubjects:", sorted(info.rat.unique()), "; files per rat:", info.rat.value_counts().sort_index().to_dict())
print("Species, age, sex:", sorted({(s["species"], s["age"], s["sex"]) for s in sessions}))
print("\n", info.to_string(index=False))
print("\nFiles without a units table:", info.session[info.units == 0].tolist())
print("Files without sleep states:", info.session[~info.states].tolist())
print("Files without ripple events:", info.session[info.ripples == 0].tolist())
print("Rats with hippocampal units:", sorted(info.rat[info.hpc > 0].unique()))
print("Files with raw voltage:", int(info.raw.sum()), "; with LFP:", int(info.lfp.sum()), "; with video:", int((info.videos > 0).sum()))
Versions: [('draft', '2026-10-05 19:47:59'), ('0.210812.1517', '2021-08-12 15:16:41'), ('0.250624.0426', '2025-06-24 04:26:54')]
Archive summary: files 40 ; subjects 4 ; species ['Brown rat'] ; measured ['Units', 'LFP', 'ElectricalSeries', 'ProcessingModule']
Draft against pinned version: assets 40 and 40 ; same asset identifiers: 40
Metadata fields that differ in the draft: ['about', 'contributor', 'ethicsApproval', 'keywords']
Related resources (pinned): [('dcite:IsDescribedBy', 'https://doi.org/10.1038/nn.4637')]
Subjects: ['Rat08', 'Rat09', 'Rat10', 'Rat11'] ; files per rat: {'Rat08': 13, 'Rat09': 10, 'Rat10': 7, 'Rat11': 10}
Species, age, sex: [('Rattus norvegicus domestica - Long Evans', '3 months', 'Male')]
session rat units hpc bla cen pir cell_type states ripples puffs rewards pre_sleep_s post_sleep_s run_s lfp raw videos trials
Rat08-20130708 Rat08 166 9 15 64 9 False True 9473 3 9 9890 10664 3603 True True 2 215.0
Rat08-20130709 Rat08 201 34 78 25 18 True True 10713 2 22 9725 11239 3621 True True 3 185.0
Rat08-20130710 Rat08 190 31 87 17 13 True True 9251 9 24 10098 10224 3607 True True 3 113.0
Rat08-20130711 Rat08 201 24 92 18 23 True True 8193 5 36 8903 8425 3612 True True 3 289.0
Rat08-20130712 Rat08 199 59 40 0 11 True True 6118 14 69 10318 9222 3328 True True 3 266.0
Rat08-20130713 Rat08 252 60 94 0 98 True True 6835 20 44 9199 9766 3244 True True 3 183.0
Rat08-20130715 Rat08 166 20 39 0 96 False True 6945 34 82 9266 9334 2730 True False 3 NaN
Rat08-20130716 Rat08 187 19 34 0 95 False True 8934 30 81 10663 11500 2486 True False 3 NaN
Rat08-20130717 Rat08 197 42 17 0 79 False True 7857 31 99 9671 10985 1797 True False 3 NaN
Rat08-20130718 Rat08 219 50 0 0 118 False False 0 49 146 11154 11121 2425 True False 3 NaN
Rat08-20130719 Rat08 157 0 0 0 109 False False 0 29 83 7054 10068 2126 True False 3 NaN
Rat08-20130720 Rat08 182 0 0 0 138 False False 0 38 92 10801 9430 1819 True False 3 NaN
Rat08-20130722 Rat08 0 0 0 0 0 False False 0 40 108 8573 11377 1856 True False 3 NaN
Rat09-20140324 Rat09 136 0 35 70 0 False True 0 20 37 9887 10655 3065 True True 3 137.0
Rat09-20140325 Rat09 97 0 49 42 0 False False 0 37 62 6871 9671 3111 True True 3 167.0
Rat09-20140326 Rat09 140 0 108 32 0 False False 0 78 162 12038 8735 2708 True True 3 276.0
Rat09-20140328 Rat09 75 0 0 0 0 False False 0 54 129 0 9763 2027 True True 2 69.0
Rat09-20140329 Rat09 67 0 67 0 0 False False 0 50 137 8932 8349 1935 True True 3 219.0
Rat09-20140402 Rat09 0 0 0 0 0 False False 0 59 158 9800 10256 2226 True True 3 271.0
Rat09-20140403 Rat09 0 0 0 0 0 False False 0 40 98 8984 7443 1299 True True 3 190.0
Rat09-20140404 Rat09 0 0 0 0 0 False False 0 52 123 8982 10688 1347 True True 3 240.0
Rat09-20140405 Rat09 0 0 0 0 0 False False 0 44 116 5513 8597 908 True True 3 185.0
Rat09-20140407 Rat09 0 0 0 0 0 False False 0 0 28 8248 5062 621 True True 1 71.0
Rat10-20140624 Rat10 101 29 37 32 0 False True 7114 26 79 8141 9106 1890 True False 3 NaN
Rat10-20140626 Rat10 76 15 33 20 0 False True 9302 20 66 8126 7554 1364 True False 3 NaN
Rat10-20140701 Rat10 40 12 6 0 3 False True 9275 14 72 9766 12766 2618 True False 3 NaN
Rat10-20140702 Rat10 61 18 2 0 7 False True 4534 38 162 8551 9870 1333 True False 3 NaN
Rat10-20140704 Rat10 57 17 9 0 28 False True 5503 9 111 8103 9854 1388 True False 3 NaN
Rat10-20140705 Rat10 49 18 0 0 22 False True 5717 53 158 10876 9559 1492 True False 3 NaN
Rat10-20140708 Rat10 56 7 0 0 44 False True 5934 24 152 7925 9390 1353 True False 3 NaN
Rat11-20150313 Rat11 136 32 0 62 0 False True 3724 19 56 11423 9948 1461 True False 3 NaN
Rat11-20150315 Rat11 127 22 0 57 0 False True 4942 30 102 13153 9933 1067 True False 3 NaN
Rat11-20150316 Rat11 190 57 0 43 0 False True 3503 46 90 11336 12113 1246 True False 3 NaN
Rat11-20150325 Rat11 161 8 86 0 0 False True 4219 32 89 14305 10695 812 True False 3 NaN
Rat11-20150326 Rat11 202 27 74 0 0 False True 3823 32 85 11065 10850 903 True False 3 NaN
Rat11-20150327 Rat11 214 68 0 0 9 False True 4446 30 82 11965 11606 1042 True False 3 NaN
Rat11-20150328 Rat11 167 50 84 0 0 False True 4130 37 111 11595 11582 1023 True False 3 NaN
Rat11-20150331 Rat11 156 83 32 0 0 False False 5155 21 71 13324 7614 1025 True False 3 NaN
Rat11-20150401 Rat11 77 18 0 0 15 False False 4384 28 78 10826 10173 941 True False 3 NaN
Rat11-20150403 Rat11 0 0 0 0 0 False False 5694 0 62 0 8550 1064 True False 1 NaN
Files without a units table: ['Rat08-20130722', 'Rat09-20140402', 'Rat09-20140403', 'Rat09-20140404', 'Rat09-20140405', 'Rat09-20140407', 'Rat11-20150403']
Files without sleep states: ['Rat08-20130718', 'Rat08-20130719', 'Rat08-20130720', 'Rat08-20130722', 'Rat09-20140325', 'Rat09-20140326', 'Rat09-20140328', 'Rat09-20140329', 'Rat09-20140402', 'Rat09-20140403', 'Rat09-20140404', 'Rat09-20140405', 'Rat09-20140407', 'Rat11-20150331', 'Rat11-20150401', 'Rat11-20150403']
Files without ripple events: ['Rat08-20130718', 'Rat08-20130719', 'Rat08-20130720', 'Rat08-20130722', 'Rat09-20140324', 'Rat09-20140325', 'Rat09-20140326', 'Rat09-20140328', 'Rat09-20140329', 'Rat09-20140402', 'Rat09-20140403', 'Rat09-20140404', 'Rat09-20140405', 'Rat09-20140407']
Rats with hippocampal units: ['Rat08', 'Rat10', 'Rat11']
Files with raw voltage: 16 ; with LFP: 40 ; with video: 40
The archive's summary says four subjects, which matches the four rats in the files, all male Long Evans rats given as 3 months old. The draft holds the same assets as the pinned version (same identifiers) and differs in descriptive metadata only, so no finding of the earlier attempt can have been fixed in the data. Seven files have no units table, sixteen have no sleep states and fourteen have no ripple events; the Rat09 files have no hippocampal units and no ripples, which fits the paper's rat whose hippocampal probe failed (rat 2). The dandiset holds 40 recording days, where the paper reports 61 sessions.
labels = pd.Series(np.concatenate([s["units"]["label"] for s in sessions if s["units"]])).value_counts()
print("Unit location labels over all files (as stored):")
print(labels.to_string())
print("\nElectrode table locations:", sorted({x for s in sessions for x in s["electrode_locations"]}))
print("Files whose units table has no location column:",
[s["session"] for s in sessions if s["units"] and not s["units"]["has_location_column"]])
print("\nUnits table columns:", sorted({tuple(sorted(s["units"]["columns"])) for s in sessions if s["units"]}))
typed = [s for s in sessions if s["units"] and "cell_type" in s["units"]["columns"]]
print("Files with a cell_type column:", [s["session"] for s in typed])
print("Its description:", sorted({s["units"]["cell_type_description"] for s in typed}))
ct = pd.DataFrame([{"session": s["session"], "region": s["units"]["region"][i], "cell_type": s["units"]["cell_type"][i],
"rate_hz": s["units"]["n_spikes"][i] / s["units"]["last_spike"]}
for s in typed for i in range(len(s["units"]["region"]))])
print(ct.groupby(["region", "cell_type"]).size().unstack(fill_value=0).to_string())
print("\nWhole-recording rate (Hz) by cell_type, all regions: median",
ct.groupby("cell_type").rate_hz.median().round(2).to_dict())
ident = ct[ct.region.isin(hpcbla.BLA_LABELS)].groupby(["session", "cell_type"]).size().unstack(fill_value=0)
print("BLA units by cell_type and file:\n", ident.to_string())
print("Files with at least", hpcbla.MIN_CELLS, "identified excitatory BLA units (M4 minimum):", int((ident.get("excitatory", 0) >= hpcbla.MIN_CELLS).sum()))
hpc_exc = ct[ct.region.isin(hpcbla.HPC_LABELS) & (ct.cell_type == "excitatory")].groupby("session").size()
print("Hippocampal units typed excitatory per file:", hpc_exc.to_dict(), "; files with 6 or more:", int((hpc_exc >= 6).sum()))
print("\nUsed as hippocampus:", sorted(hpcbla.HPC_LABELS), "; BLA:", sorted(hpcbla.BLA_LABELS), "; central nuclei:",
sorted(hpcbla.CEN_LABELS), "; piriform:", sorted(hpcbla.PIR_LABELS), "; wider amygdala variant:",
sorted(hpcbla.BLA_COMPLEX_LABELS), "(case ignored)")
Unit location labels over all files (as stored):
bla 1101
hpc 779
pir 653
cen 482
PIR 282
blv 272
bmp 264
ladl 248
den 181
ven 98
(none) 75
BLV 51
HIP 50
no label 43
PLCo1 37
cpu 35
BMP 24
BLA 17
plc01 10
Electrode table locations: [np.str_('unknown')]
Files whose units table has no location column: ['Rat09-20140328']
Units table columns: [('cell_type', 'electrode_group', 'location', 'shank_id', 'spike_times'), ('electrode_group', 'location', 'shank_id', 'spike_times'), ('electrode_group', 'shank_id', 'spike_times')]
Files with a cell_type column: ['Rat08-20130709', 'Rat08-20130710', 'Rat08-20130711', 'Rat08-20130712', 'Rat08-20130713']
Its description: ['Type of cell this has been classified as.']
cell_type excitatory inhibitory unclassified
region
bla 152 36 203
cen 15 3 42
den 26 12 119
hpc 9 11 188
ladl 3 1 17
no label 14 2 27
pir 20 8 135
Whole-recording rate (Hz) by cell_type, all regions: median {np.str_('excitatory'): 0.27, np.str_('inhibitory'): 8.91, np.str_('unclassified'): 0.42}
BLA units by cell_type and file:
cell_type excitatory inhibitory unclassified
session
Rat08-20130709 34 10 34
Rat08-20130710 38 7 42
Rat08-20130711 47 12 33
Rat08-20130712 11 3 26
Rat08-20130713 22 4 68
Files with at least 6 identified excitatory BLA units (M4 minimum): 5
Hippocampal units typed excitatory per file: {'Rat08-20130710': 1, 'Rat08-20130711': 4, 'Rat08-20130712': 2, 'Rat08-20130713': 2} ; files with 6 or more: 0
Used as hippocampus: ['hip', 'hpc'] ; BLA: ['bla'] ; central nuclei: ['cen'] ; piriform: ['pir'] ; wider amygdala variant: ['bla', 'blv', 'bmp', 'ladl'] (case ignored)
The unit labels are short codes with no description in the files or the dandiset metadata. They
match the structure names of the paper's Supplementary Fig. 3 (BLA, BLV, LaDL, BMP; Pir, VEn,
DEn; and CeN for the central nuclei), so bla, cen, pir and hpc/HIP are read as the
paper's BLA, central nuclei, piriform cortex and hippocampus (M1); this is consistent with the
paper's names, not shown by the dandiset. Some labels appear in two cases (bla and BLA,
pir and PIR, plc01 and PLCo1); others (cpu, no label) are not defined anywhere. The
electrode table gives every channel the location unknown, and Rat09-20140328 has no location
column.
Cell types. Only five Rat08 files (20130709 to 20130713) have a cell_type column, described only
as "Type of cell this has been classified as", with the values excitatory, inhibitory and
unclassified. Most units are unclassified, and almost no hippocampal unit is typed (no file has six
typed excitatory hippocampal units), while the paper's k-means step types every remaining cell as
a putative pyramidal cell or interneuron. The column therefore holds at most the paper's
monosynaptic identification, and the putative types are missing: excitatory units are read as
monosynaptically identified pyramidal cells and inhibitory ones as identified interneurons (O4
reading a). This reading is consistent with the data (inhibitory units fire faster than
excitatory ones, the BLA ratio of excitatory to inhibitory units is close to the paper's 675 to
175, and the short-latency check below finds connections of the marked kind for nearly every
marked unit), not shown by it. The files hold no spike waveforms; the raw voltage of 16 files would be
needed to classify cells, which the claims file does not provide for. The earlier attempt did not
read this column. Every claim whose criterion uses putative pyramidal cells or interneurons is
blocked as metadata_insufficient, except bla-rate-same-across-directions, which position
blocks first (too large).
# What the cell_type column marks: for each BLA unit of the five typed files, the number of loaded
# units (hippocampus, amygdala, central nuclei) it excites or inhibits at 1 to 4 ms, at three
# thresholds. Identified cells should have such partners far more often than unclassified ones.
started = time.time()
rows = []
for s in typed:
for i in hpcbla.units_of(s, hpcbla.BLA_LABELS):
rows.append({"session": s["session"], "unit": i, "cell_type": s["units"]["cell_type"][i] or "unclassified",
"rate_hz": s["units"]["n_spikes"][i] / s["units"]["last_spike"], **hpcbla.mono_connections(s, i)})
mono = pd.DataFrame(rows)
print(f"computed in {time.time() - started:.0f} s; BLA units", len(mono), "; partners searched among",
sum(len(s["units"]["trains"]) for s in typed) - len(typed), "loaded units of the same file (all files)")
table = pd.DataFrame({f"{what} at p < {p:g}": (mono[f"{what}_{p:g}"] > 0).groupby(mono.cell_type).mean()
for p in hpcbla.MONO_P for what in ("excites", "inhibits")})
table.insert(0, "units", mono.groupby("cell_type").size())
print("Share of BLA units with at least one short-latency partner:")
print(table.round(3).to_string())
print("\nMedian whole-recording rate (Hz):", {str(k): v for k, v in mono.groupby("cell_type").rate_hz.median().round(2).items()})
computed in 93 s; BLA units 391 ; partners searched among 671 loaded units of the same file (all files)
Share of BLA units with at least one short-latency partner:
units excites at p < 0.001 inhibits at p < 0.001 excites at p < 1e-06 inhibits at p < 1e-06 excites at p < 1e-10 inhibits at p < 1e-10
cell_type
excitatory 152 1.000 0.099 0.993 0.020 0.941 0.007
inhibitory 36 0.833 1.000 0.806 0.944 0.722 0.944
unclassified 203 0.714 0.177 0.443 0.079 0.310 0.034
Median whole-recording rate (Hz): {'excitatory': 0.21, 'inhibitory': 10.62, 'unclassified': 0.35}
Among the 391 BLA units of the five typed files, every unit marked excitatory (152) excites at least one loaded unit at 1 to 4 ms at p < 0.001, and 99.3% do at p < 1e-6, against 71% and 44% of the 203 unclassified units. Of the 36 units marked inhibitory, 94.4% inhibit one at p < 1e-6, against 2.0% of the excitatory and 7.9% of the unclassified units, and they fire at a median of 10.6 Hz against 0.21 Hz. So the marked units are those with short-latency partners of the marked sign, which fits the paper's monosynaptic identification and supports reading the column as O4 reading a. The check is not the paper's procedure: its windows and thresholds are a standard choice, it searches partners only among the loaded units (hippocampus, amygdala labels, central nuclei), many unclassified units also have short-latency peaks at the looser thresholds, and 81% of the inhibitory units also show an excitatory peak at p < 1e-6, which this test cannot tell from shared input. It decides no claim.
rows = []
for s in sessions:
end_raw = max(b for _, _, b in s["epochs_raw"])
lfp_s = s["lfp"][0][0] / s["lfp"][1] if s["lfp"] else np.nan
rows.append({"session": s["session"], "epochs_last_stop_raw": end_raw,
"lfp_duration_s": round(lfp_s, 1), "ratio": round(end_raw / lfp_s, 3),
"last_spike_s": round(s["units"]["last_spike"], 1) if s["units"] else np.nan,
"lfp_rate_hz": s["lfp"][1] if s["lfp"] else np.nan,
"states_span_s": (round(min(a for _, a, _ in s["states"])), round(max(b for _, _, b in s["states"])))
if s["states"] else None,
"state_labels": ",".join(sorted({lab for lab, _, _ in s["states"]})) if s["states"] else None})
times = pd.DataFrame(rows)
print("Epoch table column description:", sorted({s["epoch_unit_stated"] for s in sessions}))
print(times.to_string(index=False))
print("\nRatio of the epoch table's last stop to the LFP duration: min", times.ratio.min(), "max", times.ratio.max())
print("Epoch tags:", sorted({t for s in sessions for t, _, _ in s["epochs"]}))
print("Files that tag sleep only as 'sleep':", {s["session"]: [t for t, _, _ in s["epochs"]] for s in sessions
if "presleep" not in [t for t, _, _ in s["epochs"]] and "postsleep" not in [t for t, _, _ in s["epochs"]]})
print("Files with more than one run epoch:", {s["session"]: [(round(a), round(b)) for t, a, b in s["epochs"] if t == "run"]
for s in sessions if sum(t == "run" for t, _, _ in s["epochs"]) > 1})
Epoch table column description: ['Start time of epoch, in seconds']
session epochs_last_stop_raw lfp_duration_s ratio last_spike_s lfp_rate_hz states_span_s state_labels
Rat08-20130708 24759200.0 24759.2 1000.0 24759.2 1250.0 (5, 24754) Awake,Non-REM,REM
Rat08-20130709 25510500.0 25510.5 1000.0 25510.5 1250.0 (5, 25505) Awake,Non-REM,REM
Rat08-20130710 25195700.0 25195.7 1000.0 25195.7 1250.0 (5, 25190) Awake,Non-REM,REM
Rat08-20130711 21867700.0 21867.7 1000.0 21867.7 1250.0 (5, 21862) Awake,Non-REM,REM
Rat08-20130712 23980425.0 23980.4 1000.0 23980.4 1250.0 (5, 23975) Awake,Non-REM,REM
Rat08-20130713 22680000.0 22680.0 1000.0 22680.0 1250.0 (5, 22675) Awake,Non-REM,REM
Rat08-20130715 21827500.0 21827.5 1000.0 21827.5 1250.0 (5, 21822) Awake,Non-REM,REM
Rat08-20130716 25112700.0 25112.7 1000.0 25112.7 1250.0 (5, 25107) Awake,Non-REM,REM
Rat08-20130717 23243600.0 23243.6 1000.0 23243.6 1250.0 (5, 23238) Awake,Non-REM,REM
Rat08-20130718 25651100.0 25651.1 1000.0 25651.1 1250.0 None NaN
Rat08-20130719 19752575.0 19752.6 1000.0 19752.6 1250.0 None NaN
Rat08-20130720 22498700.0 22498.7 1000.0 22498.7 1250.0 None NaN
Rat08-20130722 22295200.0 22295.2 1000.0 NaN 1250.0 None NaN
Rat09-20140324 23930200.0 23930.2 1000.0 23930.2 1250.0 (5, 23925) Awake,Non-REM,REM
Rat09-20140325 20033500.0 20033.5 1000.0 20033.5 1250.0 None NaN
Rat09-20140326 24306500.0 24306.5 1000.0 24306.5 1250.0 None NaN
Rat09-20140328 12660300.0 12660.3 1000.0 12660.3 1250.0 None NaN
Rat09-20140329 20108100.0 20108.1 1000.0 20108.1 1250.0 None NaN
Rat09-20140402 23074300.0 23074.3 1000.0 NaN 1250.0 None NaN
Rat09-20140403 18526600.0 18526.6 1000.0 NaN 1250.0 None NaN
Rat09-20140404 21850700.0 21850.7 1000.0 NaN 1250.0 None NaN
Rat09-20140405 15987200.0 15987.2 1000.0 NaN 1250.0 None NaN
Rat09-20140407 13930600.0 13930.6 1000.0 NaN 1250.0 None NaN
Rat10-20140624 20018200.0 20018.2 1000.0 20018.2 1250.0 (5, 20013) Awake,Non-REM,REM
Rat10-20140626 17839200.0 17839.2 1000.0 17839.2 1250.0 (5, 17834) Awake,Non-REM,REM
Rat10-20140701 26196300.0 26196.3 1000.0 26196.3 1250.0 (5, 26191) Awake,Non-REM,REM
Rat10-20140702 21061687.5 21061.7 1000.0 21061.7 1250.0 (5, 21056) Awake,Non-REM,REM
Rat10-20140704 20183000.0 20183.0 1000.0 20183.0 1250.0 (5, 20178) Awake,Non-REM,REM
Rat10-20140705 22905400.0 22905.4 1000.0 22905.4 1250.0 (5, 22900) Awake,Non-REM,REM
Rat10-20140708 19519200.0 19519.2 1000.0 19519.1 1250.0 (5, 19514) Awake,Non-REM,REM
Rat11-20150313 23213900.0 23213.9 1000.0 23213.9 1250.0 (5, 23208) Awake,Non-REM,REM
Rat11-20150315 24989700.0 24989.7 1000.0 24989.7 1250.0 (5, 24984) Awake,Non-REM,REM
Rat11-20150316 25062400.0 25062.4 1000.0 25062.4 1250.0 (5, 25057) Awake,Non-REM,REM
Rat11-20150325 26214300.0 26214.3 1000.0 26214.3 1250.0 (5, 26209) Awake,Non-REM,REM
Rat11-20150326 23205600.0 23205.6 1000.0 23205.6 1250.0 (5, 23200) Awake,Non-REM,REM
Rat11-20150327 25025187.5 25025.2 1000.0 25025.2 1250.0 (5, 25020) Awake,Non-REM,REM
Rat11-20150328 25004100.0 25004.1 1000.0 25004.1 1250.0 (5, 24999) Awake,Non-REM,REM
Rat11-20150331 22777000.0 22777.0 1000.0 22777.0 1250.0 None NaN
Rat11-20150401 22477800.0 22477.8 1000.0 22477.8 1250.0 None NaN
Rat11-20150403 9613500.0 9613.5 1000.0 NaN 1250.0 None NaN
Ratio of the epoch table's last stop to the LFP duration: min 1000.0 max 1000.0
Epoch tags: [np.str_('postrun'), np.str_('postsleep'), np.str_('prerun'), np.str_('presleep'), np.str_('run'), np.str_('sleep'), np.str_('water')]
Files that tag sleep only as 'sleep': {'Rat08-20130708': [np.str_('sleep'), np.str_('run'), np.str_('sleep'), np.str_('run')], 'Rat09-20140407': [np.str_('sleep'), np.str_('run'), np.str_('sleep')], 'Rat11-20150403': [np.str_('run'), np.str_('sleep')]}
Files with more than one run epoch: {'Rat08-20130708': [(9890, 13493), (24157, 24759)]}
The epoch table says its times are in seconds, but its last stop time is 1000 times the LFP
duration in every file: the epochs are stored in milliseconds. The helper divides them by 1000.
Sleep states, spikes, ripples, puffs and rewards are in seconds and span the session; the states
(Awake, Non-REM, REM) cover the whole recording, so wake includes the track epochs. Three files
tag their sleep epochs only as sleep; there the sleep before the first run is pre-sleep and the
sleep after it post-sleep. Rat08-20130708 has a second run of about 600 s after its post-sleep;
the first run is the training run used here.
rows = []
for s in sessions:
eps = {t: (a, b) for t, a, b in s["epochs"]}
for name, (n_frames, first, last, unit) in s["videos"].items():
tag = name.rsplit("-", 1)[-1]
a, b = eps.get(tag, (np.nan, np.nan))
rows.append({"session": s["session"], "video": tag, "frames": n_frames, "stated_unit": unit,
"first_ts": round(first, 1), "last_ts": round(last, 1),
"video_span_ts/1000": round((last - first) / 1000, 1), "epoch_start_s": round(a, 1),
"epoch_duration_s": round(b - a, 1)})
videos = pd.DataFrame(rows)
print(videos.head(12).to_string(index=False))
print("...", len(videos), "videos in all")
ok = np.isfinite(videos.epoch_start_s)
print("Videos whose first timestamp / 1000 is within 5 s of their epoch's start:",
int((np.abs(videos.first_ts[ok] / 1000 - videos.epoch_start_s[ok]) < 5).sum()), "of", int(ok.sum()))
print("Videos whose span / 1000 is within 5% of their epoch's duration:",
int((np.abs(videos["video_span_ts/1000"][ok] / videos.epoch_duration_s[ok] - 1) < 0.05).sum()), "of", int(ok.sum()))
rows = []
for s in sessions:
if s["trials"] is None:
continue
a, b, desc = s["trials"]
_, _, run = hpcbla.sleep_epochs(s)
rows.append({"session": s["session"], "trials": len(a), "description": desc,
"median_duration_s": round(float(np.median(b - a)), 1),
"fraction_in_run": round(float(np.mean((a >= run[0]) & (a < run[1]))), 2)})
print("\nTrials tables:")
print(pd.DataFrame(rows).to_string(index=False))
session video frames stated_unit first_ts last_ts video_span_ts/1000 epoch_start_s epoch_duration_s
Rat08-20130708 run 54431 seconds 10350766.7 12165066.7 1814.3 24156.8 602.4
Rat08-20130708 run 9101 seconds 12165066.7 12468366.7 303.3 24156.8 602.4
Rat08-20130709 prerun 4796 seconds 12376400.0 12536200.0 159.8 0.0 317.4
Rat08-20130709 run 54697 seconds 10246933.3 12070100.0 1823.2 10042.2 3620.7
Rat08-20130709 postrun 9191 seconds 12070100.0 12376400.0 306.3 24902.1 608.4
Rat08-20130710 prerun 4764 seconds 12523766.7 12682500.0 158.7 0.0 315.4
Rat08-20130710 run 54495 seconds 10228866.7 12045300.0 1816.4 10413.6 3607.2
Rat08-20130710 postrun 9937 seconds 12045300.0 12376466.7 331.2 24245.2 658.0
Rat08-20130711 prerun 4662 seconds 10856433.3 11011766.7 155.3 0.0 308.4
Rat08-20130711 run 54593 seconds 8725300.0 10545000.0 1819.7 9211.5 3612.3
Rat08-20130711 postrun 9345 seconds 10545000.0 10856433.3 311.4 21248.9 618.8
Rat08-20130712 prerun 4461 seconds 11928733.3 12077366.7 148.6 0.0 295.3
... 114 videos in all
Videos whose first timestamp / 1000 is within 5 s of their epoch's start: 1 of 114
Videos whose span / 1000 is within 5% of their epoch's duration: 93 of 114
Trials tables:
session trials description median_duration_s fraction_in_run
Rat08-20130708 215 experimental trials 3.2 0.75
Rat08-20130709 185 experimental trials 3.1 0.61
Rat08-20130710 113 experimental trials 3.0 0.53
Rat08-20130711 289 experimental trials 3.5 0.69
Rat08-20130712 266 experimental trials 3.6 0.64
Rat08-20130713 183 experimental trials 3.5 0.80
Rat09-20140324 137 experimental trials 3.2 0.85
Rat09-20140325 167 experimental trials 3.1 0.84
Rat09-20140326 276 experimental trials 4.7 0.78
Rat09-20140328 69 experimental trials 4.0 0.00
Rat09-20140329 219 experimental trials 3.9 0.66
Rat09-20140402 271 experimental trials 3.9 0.73
Rat09-20140403 190 experimental trials 3.6 0.59
Rat09-20140404 240 experimental trials 4.1 0.64
Rat09-20140405 185 experimental trials 3.9 0.54
Rat09-20140407 71 experimental trials 4.0 1.00
The video timestamps are stated in seconds but are in milliseconds, and they are not on the
session clock: the first timestamp matches the start of the epoch a video is named after in 1 of
114 videos, while the spans match the epochs' durations in 93 of 114. The videos are not used.
The trials table of the Rat08 and Rat09 files is described only as "experimental trials"; its
intervals last about 3 to 5 s and fall partly outside the run, so they are not laps (O1 reading
a has no lap intervals to use).
# Video storage: frame shape and compression of every video series, the size of the run videos,
# and the time taken to stream frames of one run video.
storage = pd.DataFrame([{"session": s["session"], "video": k.rsplit("-", 1)[-1], "frames": shape[0],
"frame_shape": shape[1:], "dtype": dt, "compression": comp,
"uncompressed_gb": np.prod(shape) / 1e9}
for s in sessions for k, (shape, dt, comp) in s["video_storage"].items()])
print("Video series:", len(storage), "; frame shapes:", sorted(set(storage.frame_shape)), "; dtypes:", sorted(set(storage.dtype)),
"; compression:", sorted(set(map(str, storage.compression))))
run_videos = storage[storage.video == "run"]
print("Run videos:", len(run_videos), "in", run_videos.session.nunique(), "files;", int(run_videos.frames.sum()), "frames,",
round(float(run_videos.uncompressed_gb.sum())), "GB uncompressed")
first = sessions[0]
name = next(k for k in first["video_storage"] if k.endswith("-run"))
seconds, frame_bytes = hpcbla.time_frame_reads(dict(hpcbla.asset_urls())[first["path"]], name, 20)
print(f"Streaming 20 frames of {name}: {seconds:.1f} s ({seconds / 20:.2f} s per frame of {frame_bytes} bytes)")
print(f"At that rate the run frames alone would take about {run_videos.frames.sum() * seconds / 20 / 3600:.0f} hours to stream")
raw = [(s["session"], s["raw"]) for s in sessions if s["raw"]]
print("Raw voltage: files", len(raw), "; shapes and chunks, first three:", raw[:3])
Video series: 114 ; frame shapes: [(360, 640, 3)] ; dtypes: ['uint8'] ; compression: ['gzip'] Run videos: 40 in 39 files; 1920573 frames, 1328 GB uncompressed
Streaming 20 frames of Video: Rat08-20130708-02-run: 5.4 s (0.27 s per frame of 691200 bytes)
At that rate the run frames alone would take about 143 hours to stream
Raw voltage: files 16 ; shapes and chunks, first three: [('Rat08-20130708', ((495184000, 160), (241790, 1))), ('Rat08-20130709', ((510210000, 160), (249126, 1))), ('Rat08-20130710', ((503914000, 160), (246052, 1)))]
The videos are stored as raw frames (360 by 640 pixels, three colours) inside the NWB files, compressed with gzip, with no position derived from them. Tracking the rat in the run videos would mean streaming the frames counted above and running a tracker, which does not fit the time limit, and their timestamps would first have to be aligned to the session clock. The behaviour claims that need position are blocked as too large. The raw voltage of 16 files (next check) is chunked one channel at a time over about four hours of samples, so reading spike waveforms from it would mean streaming whole channels.
# What each file holds: its groups (three levels, NWB type), against the variables the archive's
# asset metadata lists for it.
asset_vars = hpcbla.asset_variables()
kinds = sorted({k for s in sessions for _, k in s["groups"] if k})
print("NWB group types in any file:", kinds)
print("Groups of type Position, SpatialSeries, BehavioralTimeSeries or CompassDirection in any file:",
[(s["session"], g) for s in sessions for g, k in s["groups"] if k in ("Position", "SpatialSeries", "BehavioralTimeSeries", "CompassDirection")] or "none")
rows = []
for s in sessions:
groups = {g for g, _ in s["groups"]}
v = asset_vars[s["path"]]
rows.append({"session": s["session"],
"acquisition": ",".join(sorted(g.split("/")[1] for g in groups if g.startswith("acquisition/") and g.count("/") == 1)),
"processing": ",".join(sorted(g.split("/", 1)[1] for g in groups if g.startswith("processing/") and g.count("/") == 2)),
"raw_in_file": s["raw"] is not None, "raw_in_metadata": "ElectricalSeries" in v,
"behavior_module_in_file": "processing/behavior" in groups, "ProcessingModule_in_metadata": "ProcessingModule" in v,
"units_in_file": s["units"] is not None, "Units_in_metadata": "Units" in v})
tree = pd.DataFrame(rows)
print(tree.to_string(index=False))
print("\nFiles with raw voltage:", int(tree.raw_in_file.sum()), "; listed with ElectricalSeries in the metadata:", int(tree.raw_in_metadata.sum()),
"; agreeing:", int((tree.raw_in_file == tree.raw_in_metadata).sum()), "of", len(tree))
print("Files with a behavior module:", int(tree.behavior_module_in_file.sum()), "; listed with ProcessingModule:",
int(tree.ProcessingModule_in_metadata.sum()), "; agreeing:", int((tree.behavior_module_in_file == tree.ProcessingModule_in_metadata).sum()))
print("Files with units:", int(tree.units_in_file.sum()), "; listed with Units:", int(tree.Units_in_metadata.sum()),
"; agreeing:", int((tree.units_in_file == tree.Units_in_metadata).sum()))
NWB group types in any file: ['AnnotationSeries', 'ElectricalSeries', 'LFP', 'ProcessingModule', 'TimeIntervals']
Groups of type Position, SpatialSeries, BehavioralTimeSeries or CompassDirection in any file: none
session acquisition processing raw_in_file raw_in_metadata behavior_module_in_file ProcessingModule_in_metadata units_in_file Units_in_metadata
Rat08-20130708 ElectricalSeries behavior/states,ecephys/LFP True False True False True True
Rat08-20130709 ElectricalSeries behavior/states,ecephys/LFP True False True False True True
Rat08-20130710 ElectricalSeries behavior/states,ecephys/LFP True False True False True True
Rat08-20130711 ElectricalSeries behavior/states,ecephys/LFP True False True False True True
Rat08-20130712 ElectricalSeries behavior/states,ecephys/LFP True False True False True True
Rat08-20130713 ElectricalSeries behavior/states,ecephys/LFP True False True False True True
Rat08-20130715 behavior/states,ecephys/LFP False True True False True True
Rat08-20130716 behavior/states,ecephys/LFP False True True False True True
Rat08-20130717 behavior/states,ecephys/LFP False True True False True True
Rat08-20130718 ecephys/LFP False True False True True True
Rat08-20130719 ecephys/LFP False True False True True True
Rat08-20130720 ecephys/LFP False True False True True True
Rat08-20130722 ecephys/LFP False True False True False False
Rat09-20140324 ElectricalSeries behavior/states,ecephys/LFP True False True False True True
Rat09-20140325 ElectricalSeries ecephys/LFP True False False True True True
Rat09-20140326 ElectricalSeries ecephys/LFP True False False True True True
Rat09-20140328 ElectricalSeries ecephys/LFP True False False True True True
Rat09-20140329 ElectricalSeries ecephys/LFP True False False True True True
Rat09-20140402 ElectricalSeries ecephys/LFP True False False True False False
Rat09-20140403 ElectricalSeries ecephys/LFP True False False True False False
Rat09-20140404 ElectricalSeries ecephys/LFP True False False True False False
Rat09-20140405 ElectricalSeries ecephys/LFP True False False True False False
Rat09-20140407 ElectricalSeries ecephys/LFP True False False True False False
Rat10-20140624 behavior/states,ecephys/LFP False True True False True True
Rat10-20140626 behavior/states,ecephys/LFP False True True False True True
Rat10-20140701 behavior/states,ecephys/LFP False True True False True True
Rat10-20140702 behavior/states,ecephys/LFP False True True False True True
Rat10-20140704 behavior/states,ecephys/LFP False True True False True True
Rat10-20140705 behavior/states,ecephys/LFP False True True False True True
Rat10-20140708 behavior/states,ecephys/LFP False True True False True True
Rat11-20150313 behavior/states,ecephys/LFP False True True False True True
Rat11-20150315 behavior/states,ecephys/LFP False True True False True True
Rat11-20150316 behavior/states,ecephys/LFP False True True False True True
Rat11-20150325 behavior/states,ecephys/LFP False True True False True True
Rat11-20150326 behavior/states,ecephys/LFP False True True False True True
Rat11-20150327 behavior/states,ecephys/LFP False True True False True True
Rat11-20150328 behavior/states,ecephys/LFP False True True False True True
Rat11-20150331 ecephys/LFP False True False True True True
Rat11-20150401 ecephys/LFP False True False True True True
Rat11-20150403 ecephys/LFP False True False True False False
Files with raw voltage: 16 ; listed with ElectricalSeries in the metadata: 24 ; agreeing: 0 of 40
Files with a behavior module: 24 ; listed with ProcessingModule: 16 ; agreeing: 0
Files with units: 33 ; listed with Units: 33 ; agreeing: 40
No file holds position: the only NWB group types under acquisition (videos aside), processing,
stimulus and intervals are AnnotationSeries, ElectricalSeries, LFP, ProcessingModule and
TimeIntervals, and no file has a Position, SpatialSeries, BehavioralTimeSeries or
CompassDirection group. The processing modules are behavior (sleep states only, 24 files) and
ecephys (LFP). Raw voltage is acquisition/ElectricalSeries in 16 files (Rat08 20130708 to
20130713 and all ten Rat09 files). The archive's asset metadata (variableMeasured) disagrees with
the files for these two in all 40 files: it lists ElectricalSeries for the 24 files that have no
raw voltage and not for the 16 that have it, and ProcessingModule for the 16 files without a
behavior module and not for the 24 with one. Its Units entries agree with the files in all 40. The
files are read directly, so the asset metadata is wrong, not the loader; this is a new finding.
rows = []
for s in sessions:
u = s["units"]
if u and u["unsorted"]:
rows.append({"session": s["session"], "unsorted_units": u["unsorted"],
"falls_back_after_spike": [x for i in u["unsorted"] for x in u["unsorted_at"][i]],
"unit0_last_spike_s": round(u["unit0_last"], 2), "unit1_first_spike_s": round(u["unit1_first"], 2),
"units_left_out": u["mixed"]})
print(pd.DataFrame(rows).to_string(index=False))
print(hpcbla.MIXED_UNITS_NOTE)
session unsorted_units falls_back_after_spike unit0_last_spike_s unit1_first_spike_s units_left_out Rat08-20130710 [1] [65535] 12430.94 12431.15 [0, 1] Rat08-20130711 [1] [65535] 6613.77 6614.18 [0, 1] Rat08-20130716 [1] [327679] 182.71 183.31 [0, 1] Rat08-20130718 [1] [196607] 8651.72 8651.73 [0, 1] Rat09-20140325 [1] [131071] 5940.00 5940.50 [0, 1] Rat09-20140329 [1] [65535] 2129.97 2130.18 [0, 1] Rat10-20140701 [1] [458751] 3133.85 3133.86 [0, 1] Rat10-20140702 [1] [65535] 8973.54 8973.84 [0, 1] Rat11-20150326 [1] [65535] 9326.75 9326.84 [0, 1] Rat11-20150327 [1] [65535] 7018.67 7018.70 [0, 1] In some files the spike times of unit 1 fall back once, after 65536 times k spikes, and unit 0 stops part way through the session where unit 1 starts: the first part of unit 1 is the end of unit 0. Units 0 and 1 of those files are left out of the analyses.
In ten files the spike times of unit 1 are not sorted: they fall back once, always after a
multiple of 65536 spikes, and unit 0 stops part way through the session at the time where unit 1
starts. The first part of unit 1 is the end of unit 0, so the index that separates them is wrong.
Units 0 and 1 of these files are left out of the analyses (they are still counted in
unit-counts, which counts units as stored).
# Clock checks. (1) Ripples against the sleep states, with the states shifted by 600 s as the
# contrast. (2) Hippocampal multi-unit rate within 25 ms of ripple peaks against 1 s away.
rows = []
for s in sessions:
if s["states"] is None or s["ripples"] is None or s["units"] is None:
continue
peaks = s["ripples"]["peak"]
r = {"session": s["session"], "ripples_ordered": s["ripples"]["ordered"],
"channels": ",".join(s["ripples"]["channels"]),
"ripple_ms_median": round(1000 * float(np.median(s["ripples"]["stop"] - s["ripples"]["start"])), 1)}
for label in ("Non-REM", "REM", "Awake"):
iv = hpcbla.intervals_of(s, label)
dur = np.diff(iv).sum()
r[f"{label}_hz"] = round(np.count_nonzero(hpcbla.in_intervals(peaks, iv)) / dur, 3)
r[f"{label}_hz_shift"] = round(np.count_nonzero(hpcbla.in_intervals(peaks, iv + 600)) / dur, 3)
u = s["units"]
hpc = hpcbla.units_of(s, hpcbla.HPC_LABELS)
if hpc:
spikes = np.sort(np.concatenate([u["trains"][i] for i in hpc]))
def rate(t, w=0.025):
return (np.searchsorted(spikes, t + w) - np.searchsorted(spikes, t - w)).mean() / (2 * w)
r["hpc_hz_at_peak"] = round(rate(peaks), 1)
r["hpc_hz_peak+1s"] = round(rate(peaks + 1), 1)
rows.append(r)
clock = pd.DataFrame(rows)
print(clock.to_string(index=False))
print("\nRipples per s of Non-REM: min", clock["Non-REM_hz"].min(), "max", clock["Non-REM_hz"].max(), "; REM max", clock["REM_hz"].max())
print("Median over files: ripples per s in Non-REM", clock["Non-REM_hz"].median(),
"(states shifted:", clock["Non-REM_hz_shift"].median(), "); in REM", clock["REM_hz"].median(),
"(shifted:", clock["REM_hz_shift"].median(), ")")
print("Median ratio of hippocampal rate at ripple peaks to 1 s later:",
round(float((clock.hpc_hz_at_peak / clock["hpc_hz_peak+1s"]).median()), 2), "; files with hippocampal units:",
int(clock.hpc_hz_at_peak.notna().sum()))
session ripples_ordered channels ripple_ms_median Non-REM_hz Non-REM_hz_shift REM_hz REM_hz_shift Awake_hz Awake_hz_shift hpc_hz_at_peak hpc_hz_peak+1s Rat08-20130708 True 2 38.4 0.835 0.427 0.019 0.561 0.130 0.319 73.6 24.8 Rat08-20130709 True 11 40.8 0.881 0.462 0.019 0.554 0.129 0.366 291.1 75.9 Rat08-20130710 True 30 42.4 0.808 0.473 0.018 0.393 0.130 0.296 183.3 52.8 Rat08-20130711 True 20 48.0 0.893 0.397 0.049 0.324 0.136 0.362 171.1 50.2 Rat08-20130712 True 23 41.6 0.770 0.366 0.022 0.337 0.062 0.208 242.9 77.6 Rat08-20130713 True 23 43.2 0.626 0.368 0.024 0.353 0.125 0.242 243.2 73.7 Rat08-20130715 True 5 46.4 0.996 0.419 0.038 0.471 0.165 0.285 63.6 16.7 Rat08-20130716 True 3 45.6 0.928 0.452 0.037 0.530 0.119 0.293 95.8 39.0 Rat08-20130717 True 13 44.8 0.781 0.372 0.017 0.304 0.100 0.318 210.3 79.2 Rat10-20140624 True 130 44.0 0.583 0.381 0.031 0.283 0.171 0.327 221.7 61.4 Rat10-20140626 True 156 44.8 0.607 0.442 0.189 0.388 0.508 0.531 126.9 53.8 Rat10-20140701 True 151 47.2 0.615 0.374 0.038 0.421 0.185 0.317 124.3 38.2 Rat10-20140702 True 153 43.2 0.499 0.234 0.012 0.214 0.121 0.203 149.9 48.0 Rat10-20140704 True 137 45.6 0.456 0.301 0.014 0.323 0.121 0.224 88.7 17.9 Rat10-20140705 True 129 44.8 0.553 0.249 0.031 0.105 0.176 0.249 154.7 47.3 Rat10-20140708 True 157 48.8 0.630 0.379 0.048 0.366 0.199 0.259 74.1 21.6 Rat11-20150313 True 142 30.4 0.349 0.175 0.002 0.141 0.063 0.148 105.9 34.3 Rat11-20150315 True 136 45.6 0.465 0.224 0.012 0.178 0.113 0.183 162.3 64.5 Rat11-20150316 True 139 43.2 0.360 0.087 0.024 0.160 0.080 0.153 392.5 117.6 Rat11-20150325 True 152 34.4 0.462 0.210 0.018 0.029 0.070 0.151 89.3 20.9 Rat11-20150326 True 134 34.4 0.327 0.165 0.010 0.110 0.074 0.167 204.4 48.1 Rat11-20150327 True 153 33.6 0.426 0.214 0.009 0.178 0.071 0.158 378.5 112.0 Rat11-20150328 True 137 36.0 0.431 0.184 0.028 0.133 0.092 0.159 271.1 97.8 Ripples per s of Non-REM: min 0.327 max 0.996 ; REM max 0.189 Median over files: ripples per s in Non-REM 0.607 (states shifted: 0.368 ); in REM 0.022 (shifted: 0.323 ) Median ratio of hippocampal rate at ripple peaks to 1 s later: 3.3 ; files with hippocampal units: 23
In every file with states and ripples, ripples occur far more often per second of Non-REM than of REM, and with the states shifted by 600 s the two rates move toward each other. Hippocampal units fire several times faster within 25 ms of a ripple peak than 1 s later. The ripple annotations come in ordered start, peak and stop triplets on one channel per file. So the states, the ripples and the spikes share one clock, as far as these contrasts show.
# (3) Lap direction decoded from hippocampal rates in each lap (leave one lap out, nearest
# centroid), with the true laps and with laps shifted by half their own duration.
def decode_direction(s, half_shift=False):
u = s["units"]
ids = hpcbla.units_of(s, hpcbla.HPC_LABELS)
laps = hpcbla.run_laps(s)
if len(ids) < 5 or len(laps) < 6 or len({d for _, _, d in laps}) < 2:
return None
if half_shift:
laps = [(a + (b - a) / 2, b + (b - a) / 2, d) for a, b, d in laps]
x = np.array([[np.count_nonzero((u["trains"][i] >= a) & (u["trains"][i] < b)) / (b - a) for i in ids]
for a, b, _ in laps])
y = np.array([d == "LR" for _, _, d in laps])
x = np.log1p(x)
x = (x - x.mean(0)) / (x.std(0) + 1e-9)
hits = 0
for k in range(len(y)):
keep = np.arange(len(y)) != k
c1, c0 = x[keep & y].mean(0), x[keep & ~y].mean(0)
hits += (np.sum((x[k] - c1) ** 2) < np.sum((x[k] - c0) ** 2)) == y[k]
return hits / len(y), len(y)
rows = []
for s in sessions:
if s["units"] is None:
continue
true, shifted = decode_direction(s), decode_direction(s, True)
if true:
rows.append({"session": s["session"], "laps": true[1], "accuracy": round(true[0], 3),
"accuracy_half_lap_shift": round(shifted[0], 3)})
dec = pd.DataFrame(rows)
print(dec.to_string(index=False))
print("Median accuracy", dec.accuracy.median(), "; with laps shifted by half a lap", dec.accuracy_half_lap_shift.median(),
"; chance 0.5")
session laps accuracy accuracy_half_lap_shift Rat08-20130710 11 1.000 0.818 Rat08-20130711 8 1.000 0.500 Rat08-20130712 22 1.000 0.818 Rat08-20130713 18 1.000 0.722 Rat08-20130715 40 0.975 0.750 Rat08-20130716 48 1.000 0.833 Rat08-20130717 52 1.000 0.942 Rat08-20130718 87 1.000 0.897 Rat10-20140624 37 0.973 0.757 Rat10-20140626 34 1.000 0.853 Rat10-20140701 26 0.923 0.808 Rat10-20140702 69 1.000 0.841 Rat10-20140704 50 1.000 0.960 Rat10-20140705 92 0.978 0.880 Rat10-20140708 81 1.000 0.877 Rat11-20150313 38 0.974 0.842 Rat11-20150315 53 0.962 0.887 Rat11-20150316 65 0.954 0.969 Rat11-20150325 60 1.000 0.850 Rat11-20150326 61 1.000 0.967 Rat11-20150327 57 1.000 0.912 Rat11-20150328 63 1.000 0.921 Rat11-20150331 32 1.000 0.938 Rat11-20150401 53 0.925 0.868 Median accuracy 1.0 ; with laps shifted by half a lap 0.8605 ; chance 0.5
Hippocampal rates tell the lap direction apart in most laps built from the reward times, and less well when the laps are shifted by half their length. The reward events and the spikes share one clock, as far as this contrast shows; it does not test the timing to better than a fraction of a lap.
# (4) Air puffs: pooled rate in the 0.5 s after each puff against the 0.5 s before, in BLA
# (wider labels) and hippocampus, with the puff times shifted by 30 s as the contrast; and the lap
# direction that holds the puffs.
rows = []
for s in sessions:
u = s["units"]
if "puf" not in s["events"]:
continue
puffs = s["events"]["puf"][1]
direction, counts = hpcbla.puff_direction(s)
row = {"session": s["session"], "puffs": len(puffs), "direction": direction, **counts}
for name, labels_ in (("bla", hpcbla.BLA_COMPLEX_LABELS), ("hpc", hpcbla.HPC_LABELS)):
ids = hpcbla.units_of(s, labels_)
if len(ids) < 3 or len(puffs) < 5:
continue
spikes = np.sort(np.concatenate([u["trains"][i] for i in ids]))
def change(t):
after = (np.searchsorted(spikes, t + 0.5) - np.searchsorted(spikes, t)).sum()
before = (np.searchsorted(spikes, t) - np.searchsorted(spikes, t - 0.5)).sum()
return (after - before) / max(before, 1)
row[f"{name}_change"] = round(change(puffs), 3)
row[f"{name}_change_shift30"] = round(change(puffs + 30), 3)
rows.append(row)
puff = pd.DataFrame(rows)
print(puff.to_string(index=False))
print("\nMedian relative change after puffs, BLA:", puff.bla_change.median(), "(shifted:", puff.bla_change_shift30.median(),
"); hippocampus:", puff.hpc_change.median(), "(shifted:", puff.hpc_change_shift30.median(), ")")
print("Files with BLA change above 0.3:", int((puff.bla_change > 0.3).sum()), "; with shifted times:",
int((puff.bla_change_shift30 > 0.3).sum()))
has_dir = puff.direction.notna()
in_dir = int(sum(r[r.direction] for _, r in puff[has_dir].iterrows()))
print("Files with puffs:", len(puff), "; without a puff direction:", puff.session[puff.direction.isna()].tolist())
print("Puffs inside run laps of the air puff direction:", in_dir, "of", int(puff.puffs.sum()),
"; in laps of the other direction:", int(sum(r["RL" if r.direction == "LR" else "LR"] for _, r in puff[has_dir].iterrows())),
"; outside run laps:", int(puff.none.sum()))
session puffs direction LR RL none bla_change bla_change_shift30 hpc_change hpc_change_shift30 Rat08-20130708 3 LR 3 0 0 NaN NaN NaN NaN Rat08-20130709 2 NaN 0 0 2 NaN NaN NaN NaN Rat08-20130710 9 RL 1 8 0 0.182 -0.074 0.040 -0.219 Rat08-20130711 5 RL 0 4 1 0.572 0.221 0.485 -0.288 Rat08-20130712 14 LR 14 0 0 0.459 0.108 0.001 0.056 Rat08-20130713 20 LR 19 1 0 0.134 -0.037 0.089 -0.064 Rat08-20130715 34 RL 0 33 1 -0.300 -0.035 0.562 -0.169 Rat08-20130716 30 LR 28 2 0 0.773 -0.044 0.148 -0.049 Rat08-20130717 31 RL 0 31 0 0.180 0.045 0.468 0.033 Rat08-20130718 49 RL 2 47 0 0.069 0.050 -0.067 0.062 Rat08-20130719 29 LR 29 0 0 -0.246 0.033 NaN NaN Rat08-20130720 38 LR 36 2 0 0.069 0.168 NaN NaN Rat08-20130722 40 RL 0 40 0 NaN NaN NaN NaN Rat09-20140324 20 RL 1 19 0 0.141 0.084 NaN NaN Rat09-20140325 37 LR 36 0 1 0.001 -0.027 NaN NaN Rat09-20140326 78 RL 4 73 1 -0.032 0.029 NaN NaN Rat09-20140328 54 RL 1 53 0 NaN NaN NaN NaN Rat09-20140329 50 RL 1 47 2 -0.012 -0.045 NaN NaN Rat09-20140402 59 LR 56 3 0 NaN NaN NaN NaN Rat09-20140403 40 RL 3 36 1 NaN NaN NaN NaN Rat09-20140404 52 LR 49 2 1 NaN NaN NaN NaN Rat09-20140405 44 LR 42 2 0 NaN NaN NaN NaN Rat10-20140624 26 RL 2 24 0 0.108 -0.035 0.019 -0.079 Rat10-20140626 20 RL 0 20 0 -0.062 0.067 0.057 0.173 Rat10-20140701 14 RL 0 14 0 -0.005 0.360 -0.044 0.272 Rat10-20140702 38 LR 38 0 0 0.155 -0.006 0.026 -0.011 Rat10-20140704 9 LR 8 1 0 -0.667 -0.273 0.179 0.217 Rat10-20140705 53 RL 0 53 0 -0.375 -0.083 -0.037 -0.017 Rat10-20140708 24 RL 0 24 0 0.206 -0.054 -0.041 -0.198 Rat11-20150313 19 RL 0 19 0 -0.056 0.117 -0.100 0.021 Rat11-20150315 30 RL 0 26 4 0.264 0.127 -0.082 -0.010 Rat11-20150316 46 RL 0 44 2 0.753 0.108 0.142 0.074 Rat11-20150325 32 RL 0 32 0 0.008 0.055 0.047 -0.183 Rat11-20150326 32 LR 32 0 0 -0.040 0.029 0.505 -0.177 Rat11-20150327 30 RL 0 30 0 0.008 -0.039 0.140 -0.048 Rat11-20150328 37 RL 0 37 0 0.077 -0.001 -0.066 0.081 Rat11-20150331 21 RL 0 21 0 -0.105 -0.017 -0.023 0.084 Rat11-20150401 28 LR 28 0 0 -0.162 -0.011 0.238 -0.038 Median relative change after puffs, BLA: 0.038500000000000006 (shifted: 0.014 ); hippocampus: 0.0435 (shifted: -0.014 ) Files with BLA change above 0.3: 4 ; with shifted times: 1 Files with puffs: 38 ; without a puff direction: ['Rat08-20130709'] Puffs inside run laps of the air puff direction: 1153 of 1197 ; in laps of the other direction: 28 ; outside run laps: 16
In the 38 files with puffs, 1153 of 1197 puffs fall inside run laps of the direction that holds most of them, 28 in laps of the other direction and 16 outside run laps. This gives the day's air puff direction for the lap-based analyses (O1 reading c laps, M8 lap types). Rat08-20130709 has no puff inside a run lap (2 puffs, 5 laps) and so no direction; two files have no puffs. The pooled response to puffs is small in most files; this is consistent with the puffs sharing the spike clock but does not show it for every file.
# Recordings deposited twice: compare unit spike counts and first spikes, ripple times and puff
# times between every pair of files, and units within a file.
def signature(s):
u = s["units"]
return {"units": None if u is None else tuple(zip(u["n_spikes"].tolist(), np.round(u["first_spike"], 4).tolist())),
"ripples": None if s["ripples"] is None else tuple(np.round(s["ripples"]["peak"][:50], 4)),
"puffs": tuple(np.round(s["events"]["puf"][1], 4)) if "puf" in s["events"] else None,
"rewards": tuple(np.round(np.sort(np.concatenate([s["events"]["lrw"][1], s["events"]["rrw"][1]])), 4))}
sigs = {s["session"]: signature(s) for s in sessions}
names = list(sigs)
same = [(a, b, k) for i, a in enumerate(names) for b in names[i + 1:] for k in sigs[a]
if sigs[a][k] is not None and sigs[a][k] == sigs[b][k]]
print("Pairs of files sharing units, ripples, puffs or rewards:", same if same else "none")
within = {s["session"]: len(sigs[s["session"]]["units"]) - len(set(sigs[s["session"]]["units"]))
for s in sessions if s["units"]}
print("Units sharing spike count and first spike within a file:", {k: v for k, v in within.items() if v} or "none")
Pairs of files sharing units, ripples, puffs or rewards: none Units sharing spike count and first spike within a file: none
No two files share their units, ripples, puffs or rewards, and no two units of a file share a spike count and first spike, so no recording was deposited twice. Units are not matched across days; like the paper, each session is analysed on its own, and a cell recorded on several days may be counted more than once.
# Sessions entering each analysis: states, pre- and post-sleep, a run, and at least 6 cells in each
# structure (M4, M7); REM needs more than 180 s in both sleep epochs (M4).
rows = []
for s in sessions:
if s["units"] is None:
continue
pre, post, run = hpcbla.sleep_epochs(s)
has_sleep = s["states"] is not None and len(pre) > 0 and len(post) > 0 and run is not None
r = {"session": s["session"], "rat": s["rat"], "sleep": has_sleep, "ripples": s["ripples"] is not None,
"hpc": len(hpcbla.units_of(s, hpcbla.HPC_LABELS)), "bla": len(hpcbla.units_of(s, hpcbla.BLA_LABELS)),
"cen": len(hpcbla.units_of(s, hpcbla.CEN_LABELS)),
"bla_wide": len(hpcbla.units_of(s, hpcbla.BLA_COMPLEX_LABELS))}
if has_sleep:
rem = hpcbla.session_epoch_rows(s, "REM")
nrem = hpcbla.session_epoch_rows(s)
r.update({"pre_nrem_s": round(nrem["pre_s"]), "post_nrem_s": round(nrem["post_s"]),
"pre_rem_s": round(rem["pre_s"]), "post_rem_s": round(rem["post_s"]),
"rem_ok": rem["pre_s"] > hpcbla.MIN_REM_S and rem["post_s"] > hpcbla.MIN_REM_S})
rows.append(r)
use = pd.DataFrame(rows)
use["hpc_bla"] = use.sleep & (use.hpc >= hpcbla.MIN_CELLS) & (use.bla >= hpcbla.MIN_CELLS)
use["hpc_cen"] = use.sleep & (use.hpc >= hpcbla.MIN_CELLS) & (use.cen >= hpcbla.MIN_CELLS)
print(use.to_string(index=False))
for col in ("hpc_bla", "hpc_cen"):
d = use[use[col]]
print(f"{col}: {len(d)} sessions {d.rat.value_counts().sort_index().to_dict()}; with REM above 180 s in both: {int(d.rem_ok.sum())}")
print("Units entering hpc_bla (hippocampus, BLA):", int(use.hpc[use.hpc_bla].sum()), int(use.bla[use.hpc_bla].sum()))
session rat sleep ripples hpc bla cen bla_wide pre_nrem_s post_nrem_s pre_rem_s post_rem_s rem_ok hpc_bla hpc_cen
Rat08-20130708 Rat08 True True 9 15 64 60 4564.0 4502.0 844.0 784.0 True True True
Rat08-20130709 Rat08 True True 34 78 25 99 4638.0 5376.0 733.0 1095.0 True True True
Rat08-20130710 Rat08 True True 29 87 17 87 3956.0 4851.0 571.0 673.0 True True True
Rat08-20130711 Rat08 True True 22 92 18 92 3496.0 3201.0 719.0 354.0 True True True
Rat08-20130712 Rat08 True True 59 40 0 40 3281.0 3125.0 377.0 217.0 True True False
Rat08-20130713 Rat08 True True 60 94 0 94 4466.0 3852.0 1088.0 915.0 True True False
Rat08-20130715 Rat08 True True 20 39 0 39 1193.0 2935.0 166.0 599.0 False True False
Rat08-20130716 Rat08 True True 17 34 0 73 2639.0 4784.0 500.0 703.0 True True False
Rat08-20130717 Rat08 True True 42 17 0 76 4780.0 3487.0 958.0 421.0 True True False
Rat08-20130718 Rat08 False False 48 0 0 51 NaN NaN NaN NaN NaN False False
Rat08-20130719 Rat08 False False 0 0 0 48 NaN NaN NaN NaN NaN False False
Rat08-20130720 Rat08 False False 0 0 0 7 NaN NaN NaN NaN NaN False False
Rat09-20140324 Rat09 True False 0 35 70 66 5543.0 5317.0 657.0 730.0 True False False
Rat09-20140325 Rat09 False False 0 49 40 55 NaN NaN NaN NaN NaN False False
Rat09-20140326 Rat09 False False 0 108 32 108 NaN NaN NaN NaN NaN False False
Rat09-20140328 Rat09 False False 0 0 0 0 NaN NaN NaN NaN NaN False False
Rat09-20140329 Rat09 False False 0 65 0 65 NaN NaN NaN NaN NaN False False
Rat10-20140624 Rat10 True True 29 37 32 40 5058.0 4370.0 613.0 894.0 True True True
Rat10-20140626 Rat10 True True 15 33 20 33 2374.0 2870.0 368.0 480.0 True True True
Rat10-20140701 Rat10 True True 12 4 0 21 6092.0 4787.0 1007.0 624.0 True False False
Rat10-20140702 Rat10 True True 18 1 0 34 3018.0 2385.0 306.0 105.0 False False False
Rat10-20140704 Rat10 True True 17 9 0 9 3818.0 5514.0 313.0 847.0 True True False
Rat10-20140705 Rat10 True True 18 0 0 6 1973.0 2722.0 294.0 280.0 True False False
Rat10-20140708 Rat10 True True 7 0 0 3 1732.0 3382.0 313.0 712.0 True False False
Rat11-20150313 Rat11 True True 32 0 62 16 5000.0 3076.0 425.0 519.0 True False True
Rat11-20150315 Rat11 True True 22 0 57 39 2837.0 3325.0 221.0 443.0 True False True
Rat11-20150316 Rat11 True True 57 0 43 90 2337.0 3098.0 253.0 122.0 False False True
Rat11-20150325 Rat11 True True 8 86 0 153 3317.0 2853.0 456.0 564.0 True True False
Rat11-20150326 Rat11 True True 27 72 0 156 4604.0 4050.0 620.0 904.0 True True False
Rat11-20150327 Rat11 True True 68 0 0 102 2128.0 5508.0 176.0 764.0 False False False
Rat11-20150328 Rat11 True True 50 84 0 117 2222.0 3341.0 386.0 669.0 True True False
Rat11-20150331 Rat11 False True 83 32 0 68 NaN NaN NaN NaN NaN False False
Rat11-20150401 Rat11 False True 18 0 0 20 NaN NaN NaN NaN NaN False False
hpc_bla: 15 sessions {'Rat08': 9, 'Rat10': 3, 'Rat11': 3}; with REM above 180 s in both: 14
hpc_cen: 9 sessions {'Rat08': 4, 'Rat10': 2, 'Rat11': 3}; with REM above 180 s in both: 8
Units entering hpc_bla (hippocampus, BLA): 438 817
Under the claims file's rule (pre- and post-sleep NREM, a run, at least six cells of each
structure), the sessions and units printed above enter the hippocampus-BLA and
hippocampus-central nuclei analyses. Rat09 has no hippocampal units, so as in the paper only three
rats enter these analyses. Which of the paper's rats 1 to 4 each of Rat08 to Rat11 is, the
dandiset does not say; Rat09 is the rat without hippocampal units, as the paper's rat 2, and the
others are taken in order (Rat08 rat 1, Rat10 rat 3, Rat11 rat 4) for unit-counts, as an
assumption.
Status of the findings of 000061-girardeau-2017-01. The version tested is the same
(0.250624.0426), and the draft holds the same assets, so none of the data findings can have been
fixed; each check above was re-run and gives the same result.
| Earlier finding | Status | Shown by |
|---|---|---|
| Epochs table in milliseconds though it says seconds | still present | epoch check (ratio 1000 in every file) |
| Video timestamps in milliseconds, stated as seconds, not on the session clock | still present | video check (1 of 114 starts match) |
spike_times_index wrong at units 0 and 1 of ten files |
still present | unsorted-unit check (ten files) |
| Files without units (7), sleep states (16) or ripples (14); Rat09 without hippocampal units | still present | file table |
| No position tracking, only raw video frames | still present | video storage check |
Undescribed unit region labels; electrode table all unknown; one file without labels |
still present in the dandiset; the paper's Supplementary Fig. 3 now gives the structure names the codes match | label check |
trials table undescribed and not laps |
still present | trials check |
| Ripple labels carry an undescribed number; events stored as stimulus annotations | still present | clock check (one channel per file) |
| Methods unavailable (paper not open, no code linked) | fixed for this attempt: the full text was supplied and the paper links its code; the dandiset still links no code (its related resources name only the paper) | header, related resources above |
New in this attempt: the cell_type column of five files, and the missing putative cell types
and waveforms (label check above).
unit-counts¶
Over 61 sessions in 4 rats, 7,390 well-isolated units were recorded (rat 1, 2,444; rat 2, 1,294; rat 3, 1,138; rat 4, 2,514), of which 2,038 were in the BLA, 782 in the central nuclei, 1,560 in the piriform cortex and 1,210 in the hippocampus.
Units as stored in the units tables, counted per rat and per location label (M1, case ignored). Rats are matched to the paper's numbering as stated in the data check.
rat_number = {"Rat08": "rat 1", "Rat09": "rat 2", "Rat10": "rat 3", "Rat11": "rat 4"}
units_all = pd.DataFrame([{"session": s["session"], "rat": s["rat"], "region": r}
for s in sessions if s["units"] for r in s["units"]["region"]])
region_sets = {"BLA": hpcbla.BLA_LABELS, "central nuclei": hpcbla.CEN_LABELS, "piriform cortex": hpcbla.PIR_LABELS,
"hippocampus": hpcbla.HPC_LABELS}
expected_counts = [("sessions", 61, len(sessions)), ("units, all rats", 7390, len(units_all))]
for rat, paper_n in (("Rat08", 2444), ("Rat09", 1294), ("Rat10", 1138), ("Rat11", 2514)):
expected_counts.append((f"units, {rat_number[rat]} ({rat})", paper_n, int((units_all.rat == rat).sum())))
for name, paper_n in (("BLA", 2038), ("central nuclei", 782), ("piriform cortex", 1560), ("hippocampus", 1210)):
expected_counts.append((f"{name} units", paper_n, int(units_all.region.isin(region_sets[name]).sum())))
counts = pd.DataFrame(expected_counts, columns=["quantity", "paper", "dandiset"])
counts["tolerance"] = 0.1 * counts.paper
counts["within"] = (counts.dandiset - counts.paper).abs() <= counts.tolerance
print(counts.to_string(index=False))
print("\nFiles with a units table:", units_all.session.nunique(), "of", len(sessions))
print("Units by rat and region:")
print(units_all.assign(region=units_all.region.map(lambda r: next((k for k, v in region_sets.items() if r in v), "other: " + r)))
.groupby(["region", "rat"]).size().unstack(fill_value=0).to_string())
print("\nVariant, wider amygdala labels (bla, blv, bmp, ladl):", int(units_all.region.isin(hpcbla.BLA_COMPLEX_LABELS).sum()))
quantity paper dandiset tolerance within
sessions 61 40 6.1 False
units, all rats 7390 4702 739.0 False
units, rat 1 (Rat08) 2444 2317 244.4 True
units, rat 2 (Rat09) 1294 515 129.4 False
units, rat 3 (Rat10) 1138 440 113.8 False
units, rat 4 (Rat11) 2514 1430 251.4 False
BLA units 2038 1118 203.8 False
central nuclei units 782 482 78.2 False
piriform cortex units 1560 935 156.0 False
hippocampus units 1210 829 121.0 False
Files with a units table: 33 of 40
Units by rat and region:
rat Rat08 Rat09 Rat10 Rat11
region
BLA 496 259 87 276
central nuclei 124 144 52 162
hippocampus 348 0 116 365
other: (none) 0 75 0 0
other: blv 153 0 59 111
other: bmp 51 0 4 233
other: cpu 0 0 0 35
other: den 181 0 0 0
other: ladl 66 37 0 145
other: no label 43 0 0 0
other: plc01 0 0 10 0
other: plco1 37 0 0 0
other: ven 11 0 8 79
piriform cortex 807 0 104 24
Variant, wider amygdala labels (bla, blv, bmp, ladl): 1977
Result: not reproduced. The dandiset holds 40 sessions (33 with a units table) and 4702 units, against the paper's 61 sessions and 7,390 units, both outside the 10% tolerance. Only Rat08, taken as rat 1, is within tolerance (2317 units against 2,444); Rat09, Rat10 and Rat11 hold 515, 440 and 1430 units against 1,294, 1,138 and 2,514. By label the dandiset has 1118 BLA, 482 central nuclei, 935 piriform and 829 hippocampal units, against 2,038, 782, 1,560 and 1,210; with the wider amygdala labels it has 1977 amygdala units. The dandiset holds part of the paper's data.
rats-learn-puff-location, each-rat-learns-puff-location, dz-speed-slower-post-run, puff-location-no-speed-bias (blocked)¶
Rats learned the daily location of the air puff. The speed ratio between the previous and the current danger zone rises from pre-run to run and post-run, reversing between pre-run and post-run.
The four behaviour claims need the rat's speed in the 20 cm before the air puff location of the
day and of the previous day, which needs position on the track and the location of each puff.
The dandiset holds no position tracking and no puff location, only the puff times and the raw
video frames, and tracking the frames does not fit the time limit (data check: run frames and
streaming time). each-rat-learns-puff-location, dz-speed-slower-post-run and
puff-location-no-speed-bias depend on rats-learn-puff-location and are blocked with it.
Result: blocked (too large). No analysis is run.
habituation-over-days¶
The aversive valence of the air puff diminishes over training days; rats run more laps per unit time as they habituate.
Habituation index per session = laps / run epoch duration x 100 (Performance.m). Laps are O1 reading c (between consecutive rewards at opposite ends, both inside the run): readings a (lap intervals the dataset holds) and b (position) cannot be run, as the dandiset has neither. As answered (questions/000061.yaml), the criterion now judges the claim on reading c alone in that case, the back-and-forth runs between the two reward ends that the Methods count. The three lap counts of the criterion: all laps (back-and-forth), air puff laps, and left-to-right laps. A reading c lap includes the time on the starting platform after its reward, so lap counts are those of back-and-forth runs between rewards. Readings of "gradually diminished": (a) pooled Spearman correlation of day rank within the rat and index, positive with p < 0.05; (b) positive Spearman correlation in at least 3 of the 4 rats (shown as the sign alone and with p < 0.05).
rows = []
for s in sessions:
_, _, run = hpcbla.sleep_epochs(s)
laps_ = hpcbla.run_laps(s)
direction, _ = hpcbla.puff_direction(s)
duration = run[1] - run[0]
rows.append({"session": s["session"], "rat": s["rat"], "run_s": round(duration), "laps": len(laps_),
"back_and_forth": 100 * len(laps_) / duration,
"airpuff": 100 * sum(d == direction for _, _, d in laps_) / duration if direction else np.nan,
"left_to_right": 100 * sum(d == "LR" for _, _, d in laps_) / duration})
hab = pd.DataFrame(rows)
hab["day"] = hab.groupby("rat").session.rank().astype(int)
print(hab.round(3).to_string(index=False))
rows = []
for kind in ("back_and_forth", "airpuff", "left_to_right"):
d = hab.dropna(subset=[kind])
rho, p = stats.spearmanr(d.day, d[kind])
per_rat = {rat: stats.spearmanr(g.day, g[kind]) for rat, g in d.groupby("rat")}
rows.append({"laps": kind, "sessions": len(d), "pooled_rho": round(rho, 3), "pooled_p": p,
"reading_a": bool(rho > 0 and p < 0.05),
**{f"rho_{rat}": round(r.statistic, 2) for rat, r in per_rat.items()},
**{f"p_{rat}": round(r.pvalue, 4) for rat, r in per_rat.items()},
"rats_positive": sum(r.statistic > 0 for r in per_rat.values()),
"rats_positive_p05": sum(r.statistic > 0 and r.pvalue < 0.05 for r in per_rat.values())})
hab_tests = pd.DataFrame(rows)
hab_tests["reading_b"] = hab_tests.rats_positive >= 3
print(hab_tests.to_string(index=False))
session rat run_s laps back_and_forth airpuff left_to_right day
Rat08-20130708 Rat08 3603 3 0.083 0.056 0.056 1
Rat08-20130709 Rat08 3621 5 0.138 NaN 0.055 2
Rat08-20130710 Rat08 3607 11 0.305 0.166 0.139 3
Rat08-20130711 Rat08 3612 8 0.221 0.111 0.111 4
Rat08-20130712 Rat08 3328 22 0.661 0.331 0.331 5
Rat08-20130713 Rat08 3244 18 0.555 0.277 0.277 6
Rat08-20130715 Rat08 2730 40 1.465 0.733 0.733 7
Rat08-20130716 Rat08 2486 48 1.931 0.965 0.965 8
Rat08-20130717 Rat08 1797 52 2.894 1.447 1.447 9
Rat08-20130718 Rat08 2425 87 3.588 1.773 1.815 10
Rat08-20130719 Rat08 2126 55 2.587 1.270 1.270 11
Rat08-20130720 Rat08 1819 70 3.848 1.924 1.924 12
Rat08-20130722 Rat08 1856 80 4.309 2.155 2.155 13
Rat09-20140324 Rat09 3065 21 0.685 0.359 0.326 1
Rat09-20140325 Rat09 3111 43 1.382 0.707 0.707 2
Rat09-20140326 Rat09 2708 108 3.988 1.994 1.994 3
Rat09-20140328 Rat09 2027 81 3.996 1.973 2.022 4
Rat09-20140329 Rat09 1935 78 4.031 2.016 2.016 5
Rat09-20140402 Rat09 2226 101 4.538 2.246 2.246 6
Rat09-20140403 Rat09 1299 63 4.850 2.386 2.463 7
Rat09-20140404 Rat09 1347 84 6.238 3.119 3.119 8
Rat09-20140405 Rat09 908 54 5.946 2.973 2.973 9
Rat09-20140407 Rat09 621 27 4.348 NaN 2.254 10
Rat10-20140624 Rat10 1890 37 1.958 1.006 0.953 1
Rat10-20140626 Rat10 1364 34 2.493 1.246 1.246 2
Rat10-20140701 Rat10 2618 26 0.993 0.497 0.497 3
Rat10-20140702 Rat10 1333 69 5.175 2.550 2.550 4
Rat10-20140704 Rat10 1388 50 3.602 1.801 1.801 5
Rat10-20140705 Rat10 1492 92 6.165 3.082 3.082 6
Rat10-20140708 Rat10 1353 81 5.987 3.031 2.957 7
Rat11-20150313 Rat11 1461 38 2.600 1.300 1.300 1
Rat11-20150315 Rat11 1067 53 4.965 2.436 2.530 2
Rat11-20150316 Rat11 1246 65 5.218 2.569 2.649 3
Rat11-20150325 Rat11 812 60 7.386 3.693 3.693 4
Rat11-20150326 Rat11 903 61 6.752 3.431 3.431 5
Rat11-20150327 Rat11 1042 57 5.469 2.783 2.687 6
Rat11-20150328 Rat11 1023 63 6.159 3.128 3.031 7
Rat11-20150331 Rat11 1025 32 3.121 1.560 1.560 8
Rat11-20150401 Rat11 941 53 5.632 2.869 2.869 9
Rat11-20150403 Rat11 1064 61 5.734 NaN 2.914 10
laps sessions pooled_rho pooled_p reading_a rho_Rat08 rho_Rat09 rho_Rat10 rho_Rat11 p_Rat08 p_Rat09 p_Rat10 p_Rat11 rats_positive rats_positive_p05 reading_b
back_and_forth 40 0.383 0.014743 True 0.97 0.87 0.82 0.33 0.0 0.0012 0.0234 0.3466 4 3 True
airpuff 37 0.331 0.045660 True 0.97 0.97 0.82 0.32 0.0 0.0000 0.0234 0.4064 4 3 True
left_to_right 40 0.397 0.011277 True 0.97 0.90 0.82 0.33 0.0 0.0003 0.0234 0.3466 4 3 True
fig, ax = plt.subplots(figsize=(5, 3.5))
for rat, d in hab.groupby("rat"):
ax.plot(d.day, d.back_and_forth, "o-", label=rat)
ax.set_xlabel("training day within rat")
ax.set_ylabel("laps per 100 s of run")
ax.legend(fontsize=8)
plt.show()
print("Plotted range:", round(hab.back_and_forth.min(), 3), "to", round(hab.back_and_forth.max(), 3), "; no point outside the axes")
Plotted range: 0.083 to 7.386 ; no point outside the axes
Result: reproduced, under O1 reading c, the only lap reading the dandiset allows and, as answered, the one the criterion now judges the claim on. The habituation index rises with training day for each lap count: pooled Spearman rho 0.383 (p = 0.015, 40 sessions) for all laps, 0.331 (p = 0.046, 37 sessions) for air puff laps and 0.397 (p = 0.011, 40 sessions) for left-to-right laps. The correlation is positive in all four rats for each count, with p < 0.05 in Rat08, Rat09 and Rat10 and not in Rat11 (rho 0.32 to 0.33, p of 0.35 to 0.41). Both readings of "gradually diminished" hold for all three lap counts; the air puff count is the closest to its threshold (p = 0.046).
bla-pyr-sleep-firing-increase¶
BLA pyramidal cells fire faster in REM sleep and, to a lesser extent, in NREM sleep than in wakefulness.
Per BLA cell, rates over all REM, all NREM and all wake intervals of the session (StatesFiringRates.m), ratios (state - wake)/(state + wake), and a two-sided Wilcoxon signed rank test of each ratio against 0. Cells per O4 reading a: the dandiset's identified excitatory BLA units (five Rat08 files); reading b (all putative pyramidal cells) cannot be run, as the putative types are not in the dandiset. As answered (questions/000061.yaml), O4 now judges the claim on reading a alone: the Fig. 2c legend's cells are the monosynaptically identified ones. The paper's cells come from four rats; these come from one, a smaller sample.
rows = []
for s in typed:
for ctype in ("excitatory", "inhibitory"):
for r in hpcbla.state_rates(s, hpcbla.units_of(s, hpcbla.BLA_LABELS, ctype)):
rows.append({"session": s["session"], "type": ctype, **r})
state_fr = pd.DataFrame(rows)
for k in ("rem", "nrem"):
state_fr[f"{k}_ratio"] = (state_fr[k] - state_fr.wake) / (state_fr[k] + state_fr.wake)
def ratio_tests(d):
out = {}
for k in ("rem", "nrem"):
x = d[f"{k}_ratio"].dropna()
out.update({f"{k}_n": len(x), f"{k}_median": round(float(x.median()), 4),
f"{k}_positive": int((x > 0).sum()), f"{k}_p": stats.wilcoxon(x).pvalue if len(x) > 1 else np.nan})
return pd.Series(out)
print("Cells with an undefined ratio (no spike in either state):",
state_fr.groupby("type")[["rem_ratio", "nrem_ratio"]].apply(lambda d: d.isna().sum()).to_dict())
print("\nIdentified BLA pyramidal cells (excitatory):")
pyr = state_fr[state_fr.type == "excitatory"]
print(ratio_tests(pyr).to_string())
print("\nBy session:")
print(pyr.groupby("session").apply(ratio_tests).to_string())
print("\nMedian rates (Hz): REM", round(pyr.rem.median(), 3), "; NREM", round(pyr.nrem.median(), 3), "; wake", round(pyr.wake.median(), 3))
Cells with an undefined ratio (no spike in either state): {'rem_ratio': {'excitatory': 0, 'inhibitory': 0}, 'nrem_ratio': {'excitatory': 0, 'inhibitory': 0}}
Identified BLA pyramidal cells (excitatory):
rem_n 1.520000e+02
rem_median 3.915000e-01
rem_positive 1.200000e+02
rem_p 6.995842e-15
nrem_n 1.520000e+02
nrem_median 1.793000e-01
nrem_positive 1.080000e+02
nrem_p 7.191666e-12
By session:
rem_n rem_median rem_positive rem_p nrem_n nrem_median nrem_positive nrem_p
session
Rat08-20130709 34.0 0.3546 30.0 1.502573e-06 34.0 0.0929 23.0 0.002046
Rat08-20130710 38.0 0.4602 34.0 3.252353e-09 38.0 0.2660 28.0 0.000008
Rat08-20130711 47.0 0.4467 36.0 4.388972e-05 47.0 0.2327 37.0 0.000003
Rat08-20130712 11.0 0.3682 8.0 3.222656e-02 11.0 0.1403 8.0 0.067383
Rat08-20130713 22.0 0.1742 12.0 5.660677e-01 22.0 0.0265 12.0 0.502840
Median rates (Hz): REM 0.396 ; NREM 0.257 ; wake 0.15
Result: reproduced, on the 152 identified excitatory BLA units of five Rat08 sessions (O4 reading a; the paper had 675 cells from four rats). The REM/wake ratio has a median of 0.39 (positive in 120 of 152, p = 7.0e-15) and the NREM/wake ratio a median of 0.18 (positive in 108, p = 7.2e-12), so both are positive and the REM median exceeds the NREM one. By session the REM ratio is positive with p < 0.05 in four of five sessions and the NREM ratio in three; Rat08-20130713 (22 cells) reaches neither.
bla-int-sleep-firing¶
BLA interneurons do not change their firing rate between REM sleep and wake, and fire slightly less in NREM sleep than in wake.
As above, for the identified inhibitory BLA units (O4 reading a, which as answered decides the claim alone). The REM part is a claim of absence, judged on the two-sided Wilcoxon signed rank test at p < 0.05.
inter = state_fr[state_fr.type == "inhibitory"]
print("Identified BLA interneurons (inhibitory):")
print(ratio_tests(inter).to_string())
print("\nBy session:")
print(inter.groupby("session").apply(ratio_tests).to_string())
print("\nMedian rates (Hz): REM", round(inter.rem.median(), 3), "; NREM", round(inter.nrem.median(), 3), "; wake", round(inter.wake.median(), 3))
Identified BLA interneurons (inhibitory):
rem_n 36.000000
rem_median 0.028600
rem_positive 22.000000
rem_p 0.020054
nrem_n 36.000000
nrem_median -0.035000
nrem_positive 16.000000
nrem_p 0.264819
By session:
rem_n rem_median rem_positive rem_p nrem_n nrem_median nrem_positive nrem_p
session
Rat08-20130709 10.0 0.0196 6.0 0.431641 10.0 -0.0509 4.0 0.232422
Rat08-20130710 7.0 -0.0165 3.0 0.812500 7.0 -0.0767 2.0 0.578125
Rat08-20130711 12.0 0.0886 9.0 0.020996 12.0 0.0369 7.0 0.380371
Rat08-20130712 3.0 0.0235 2.0 0.750000 3.0 -0.0475 1.0 0.750000
Rat08-20130713 4.0 0.0279 2.0 0.625000 4.0 -0.0283 2.0 0.625000
Median rates (Hz): REM 12.003 ; NREM 10.994 ; wake 10.705
# Context (decides nothing): effect sizes against the paper's, and bootstrap intervals of the
# median ratios.
rng_int = np.random.default_rng(6102)
for k, paper_z in (("rem", 1.25), ("nrem", -6.00)):
x = inter[f"{k}_ratio"].dropna().to_numpy()
lo, hi = hpcbla.median_ci(x, rng_int)
p = stats.wilcoxon(x).pvalue
print(f"{k.upper()}/wake ratio, {len(x)} cells: median {np.median(x):.3f}, bootstrap 95% interval {lo:.3f} to {hi:.3f};"
f" effect size r here {hpcbla.effect_r(p, len(x), np.median(x)):.2f}, the paper {paper_z / np.sqrt(175):.2f} (175 cells);"
f" shift detected with 80% power {hpcbla.min_detectable_shift(x, len(x), rng_int):.3f}")
x = pyr.rem_ratio.dropna().to_numpy()
print(f"For scale, identified excitatory cells: REM/wake median {np.median(x):.3f}, interval {hpcbla.median_ci(x, rng_int)[0]:.3f} to {hpcbla.median_ci(x, rng_int)[1]:.3f}")
REM/wake ratio, 36 cells: median 0.029, bootstrap 95% interval -0.015 to 0.109; effect size r here 0.39, the paper 0.09 (175 cells); shift detected with 80% power 0.039 NREM/wake ratio, 36 cells: median -0.035, bootstrap 95% interval -0.077 to 0.032; effect size r here -0.19, the paper -0.45 (175 cells); shift detected with 80% power 0.068 For scale, identified excitatory cells: REM/wake median 0.392, interval 0.321 to 0.460
Result: not reproduced. On the 36 identified inhibitory BLA units of five Rat08 sessions, the REM/wake ratio has a median of 0.029, positive in 22 of 36, with p = 0.020: the paper's test finds a change between REM and wake, where the claim says there is none. The NREM/wake ratio has a negative median (-0.035, positive in 16 of 36) but p = 0.26, so the decrease is not significant. Neither part of the criterion holds. The cells are about a fifth of the paper's 175, from one rat, and by session only Rat08-20130711 (12 cells) has a REM ratio with p < 0.05. As context, the REM effect is small (median 0.029, against 0.39 for the excitatory cells), and the bootstrap 95% interval of the median (-0.015 to 0.109) includes zero while the signed rank test, which uses all the ranks, has p = 0.020. Its effect size r is 0.39 against the paper's 0.09; no test here compares the two samples, and the interval of the median includes zero. The NREM effect is weaker than the paper's (r -0.19 against -0.45), and with 36 cells only a shift of the median of about 0.068 would be detected with 80% power, larger than the observed -0.035.
hpc-bla-reactivation-post-sleep (blocked)¶
Hippocampus-BLA pyramidal cell pair correlations of the run are reinstated in post-sleep NREM, measured as explained variance (EV) exceeding the reverse explained variance (REV) across sessions.
The criterion uses putative pyramidal cells of the hippocampus and the BLA. The dandiset has no
putative types, and even its identified excitatory units leave no session with six hippocampal
cells (data check). The all-cell version of this test is the claim reactivation-all-cell-pairs,
below. reactivation-decay (pyramidal pairs, depends on this claim) and
piriform-weaker-reactivation (pyramidal pairs) are blocked for the same reason.
Result: blocked (metadata insufficient). No analysis is run here.
reactivation-all-cell-pairs¶
Hippocampus-BLA NREM reactivation holds when pyramidal cells and interneurons are both included.
EV and REV per session (M4) over all hippocampus-BLA pairs, from 50 ms bins of the run epoch and of the NREM of pre- and post-sleep; two-sided Wilcoxon signed rank test of EV against REV across sessions. The successive 20-min periods of post-NREM (24,000 bins each, with the whole pre-NREM and run) are reported, as the criterion says, and decide nothing.
started = time.time()
ev_rows = []
for s in sessions:
r = use[use.session == s["session"]]
if r.empty or not r.sleep.iloc[0]:
continue
hpc = hpcbla.units_of(s, hpcbla.HPC_LABELS)
nrem_rows, rem_rows = hpcbla.session_epoch_rows(s), hpcbla.session_epoch_rows(s, "REM")
row = {"session": s["session"], "rat": s["rat"]}
for name, labels_ in (("bla", hpcbla.BLA_LABELS), ("cen", hpcbla.CEN_LABELS), ("bla_wide", hpcbla.BLA_COMPLEX_LABELS)):
other = hpcbla.units_of(s, labels_)
if len(hpc) < hpcbla.MIN_CELLS or len(other) < hpcbla.MIN_CELLS:
continue
res = hpcbla.ev_session(s, hpc, other, nrem_rows, post_periods=(name == "bla"))
row.update({f"{name}_{k}": v for k, v in res.items()})
if bool(r.rem_ok.iloc[0]) and name != "bla_wide":
res = hpcbla.ev_session(s, hpc, other, rem_rows)
row.update({f"{name}_rem_{k}": v for k, v in res.items()})
ev_rows.append(row)
ev = pd.DataFrame(ev_rows)
print(f"computed in {time.time() - started:.0f} s")
def ev_test(df, prefix, label):
d = df.dropna(subset=[f"{prefix}_ev"])
diff = d[f"{prefix}_ev"] - d[f"{prefix}_rev"]
res = stats.wilcoxon(d[f"{prefix}_ev"], d[f"{prefix}_rev"])
print(f"{label}: {len(d)} sessions {d.rat.value_counts().sort_index().to_dict()}; EV median {d[f'{prefix}_ev'].median():.3f}%,"
f" REV median {d[f'{prefix}_rev'].median():.3f}%; EV - REV mean {diff.mean():.3f}% (s.e.m. {diff.sem():.3f}),"
f" median {diff.median():.3f}%, EV > REV in {int((diff > 0).sum())}; Wilcoxon two-sided p = {res.pvalue:.2g}")
for rat, g in d.groupby("rat"):
x = g[f"{prefix}_ev"] - g[f"{prefix}_rev"]
print(f" {rat}: {len(g)} sessions, EV - REV median {x.median():.3f}%, EV > REV in {int((x > 0).sum())}")
return d, diff
cols = ["session", "rat", "bla_ev", "bla_rev", "bla_pairs", "bla_ev_0", "bla_rev_0", "bla_ev_1", "bla_rev_1", "bla_ev_2", "bla_rev_2", "bla_bins_2"]
print(ev[cols].dropna(subset=["bla_ev"]).round(3).to_string(index=False))
d_all, diff_all = ev_test(ev, "bla", "Hippocampus-BLA, all cells, NREM")
computed in 16 s
session rat bla_ev bla_rev bla_pairs bla_ev_0 bla_rev_0 bla_ev_1 bla_rev_1 bla_ev_2 bla_rev_2 bla_bins_2
Rat08-20130708 Rat08 2.254 0.293 135.0 4.387 0.189 1.169 0.200 0.191 0.117 24000.0
Rat08-20130709 Rat08 6.115 0.124 2652.0 6.777 0.795 2.180 2.257 3.328 1.265 24000.0
Rat08-20130710 Rat08 5.403 0.548 2523.0 3.942 1.850 2.174 1.844 1.626 2.168 24000.0
Rat08-20130711 Rat08 5.606 0.000 2024.0 5.691 0.198 1.379 0.556 1.667 0.698 16020.0
Rat08-20130712 Rat08 1.076 0.004 2360.0 1.640 0.002 0.700 0.002 0.000 0.078 14500.0
Rat08-20130713 Rat08 1.434 0.224 5546.0 1.542 0.483 0.563 0.553 0.263 0.699 24000.0
Rat08-20130715 Rat08 8.444 0.084 780.0 6.640 0.162 6.305 0.139 0.401 0.507 10700.0
Rat08-20130716 Rat08 6.448 0.072 578.0 5.757 0.001 3.922 0.023 4.841 0.057 24000.0
Rat08-20130717 Rat08 13.226 0.097 714.0 9.220 0.199 7.124 0.381 6.334 0.485 21740.0
Rat10-20140624 Rat10 0.785 0.821 1073.0 0.008 2.434 0.526 1.424 1.113 1.179 24000.0
Rat10-20140626 Rat10 3.367 2.750 495.0 1.489 4.774 3.548 3.230 1.848 4.348 9400.0
Rat10-20140704 Rat10 1.765 1.788 153.0 0.000 0.696 1.178 1.556 0.101 0.805 24000.0
Rat11-20150325 Rat11 2.156 0.020 672.0 1.218 0.199 0.604 0.346 2.376 0.304 9060.0
Rat11-20150326 Rat11 0.904 0.536 1800.0 1.396 0.539 1.071 0.682 0.338 1.380 24000.0
Rat11-20150328 Rat11 1.768 1.150 4200.0 1.802 1.693 1.493 1.697 0.449 2.960 18820.0
Hippocampus-BLA, all cells, NREM: 15 sessions {'Rat08': 9, 'Rat10': 3, 'Rat11': 3}; EV median 2.254%, REV median 0.224%; EV - REV mean 3.483% (s.e.m. 0.985), median 1.960%, EV > REV in 13; Wilcoxon two-sided p = 0.00031
Rat08: 9 sessions, EV - REV median 5.605%, EV > REV in 9
Rat10: 3 sessions, EV - REV median -0.023%, EV > REV in 1
Rat11: 3 sessions, EV - REV median 0.618%, EV > REV in 3
print("Reported, not deciding: successive 20-min periods of post-NREM (all sessions of the test above)")
for k in range(3):
x = d_all[f"bla_ev_{k}"] - d_all[f"bla_rev_{k}"]
ok = x.dropna()
print(f" period {k + 1}: {len(ok)} sessions, EV - REV mean {ok.mean():.3f}%, EV > REV in {int((ok > 0).sum())},"
f" Wilcoxon p = {stats.wilcoxon(ok).pvalue:.2g}")
print("Sessions whose third period is shorter than 20 min:", int((d_all.bla_bins_2 < hpcbla.PERIOD_BINS).sum()),
"; post-NREM shorter than 40 min (no third period):", int(d_all.bla_ev_2.isna().sum()))
_ = ev_test(ev, "bla_wide", "Variant, wider amygdala labels (bla, blv, bmp, ladl)")
fig, ax = plt.subplots(figsize=(4, 4))
ax.scatter(d_all.bla_rev, d_all.bla_ev, s=16)
lim = [0, max(d_all.bla_ev.max(), d_all.bla_rev.max()) * 1.05]
ax.plot(lim, lim, "k--", lw=0.8)
ax.set_xlabel("REV (%)")
ax.set_ylabel("EV (%)")
plt.show()
print("Plotted ranges: EV", round(d_all.bla_ev.min(), 3), "to", round(d_all.bla_ev.max(), 3), "; REV",
round(d_all.bla_rev.min(), 3), "to", round(d_all.bla_rev.max(), 3), "; no point outside the axes")
Reported, not deciding: successive 20-min periods of post-NREM (all sessions of the test above)
period 1: 15 sessions, EV - REV mean 2.486%, EV > REV in 12, Wilcoxon p = 0.018
period 2: 15 sessions, EV - REV mean 1.270%, EV > REV in 11, Wilcoxon p = 0.035
period 3: 15 sessions, EV - REV mean 0.522%, EV > REV in 6, Wilcoxon p = 1
Sessions whose third period is shorter than 20 min: 7 ; post-NREM shorter than 40 min (no third period): 0
Variant, wider amygdala labels (bla, blv, bmp, ladl): 22 sessions {'Rat08': 9, 'Rat10': 6, 'Rat11': 7}; EV median 1.793%, REV median 0.145%; EV - REV mean 1.924% (s.e.m. 0.519), median 1.530%, EV > REV in 17; Wilcoxon two-sided p = 0.0005
Rat08: 9 sessions, EV - REV median 3.660%, EV > REV in 9
Rat10: 6 sessions, EV - REV median 0.080%, EV > REV in 3
Rat11: 7 sessions, EV - REV median 1.562%, EV > REV in 5
Plotted ranges: EV 0.785 to 13.226 ; REV 0.0 to 2.75 ; no point outside the axes
Result: reproduced. Over 15 sessions (9 Rat08, 3 Rat10, 3 Rat11; the paper had 27), EV exceeds REV in 13, with a mean EV - REV of 3.48% (s.e.m. 0.99) and a two-sided Wilcoxon p of 0.00031. By rat, EV exceeds REV in all 9 Rat08 and all 3 Rat11 sessions, and in 1 of 3 Rat10 sessions (median -0.023%). The wider amygdala labels give the same result (22 sessions, 17 with EV > REV, p = 0.0005). Reported only: EV - REV falls over the successive 20-min periods of post-NREM (means 2.49, 1.27 and 0.52%; p = 0.018, 0.035 and 1).
no-rem-reactivation (blocked)¶
Hippocampus-BLA and hippocampus-piriform cortex reactivations are not significant in REM sleep.
The criterion tests pyramidal pairs, for hippocampus-BLA and hippocampus-piriform, which needs the putative types. The criterion's variant, hippocampus-BLA with all cells in REM (sessions with more than 180 s of REM in each sleep epoch), is shown and decides nothing.
_ = ev_test(ev, "bla_rem", "Variant, hippocampus-BLA, all cells, REM")
Variant, hippocampus-BLA, all cells, REM: 14 sessions {'Rat08': 8, 'Rat10': 3, 'Rat11': 3}; EV median 0.070%, REV median 0.032%; EV - REV mean 0.339% (s.e.m. 0.315), median 0.040%, EV > REV in 8; Wilcoxon two-sided p = 0.43
Rat08: 8 sessions, EV - REV median 0.040%, EV > REV in 5
Rat10: 3 sessions, EV - REV median 0.048%, EV > REV in 2
Rat11: 3 sessions, EV - REV median -0.196%, EV > REV in 1
Result: blocked (metadata insufficient). The printed variant is a description only: with all cells, hippocampus-BLA EV does not differ from REV in REM (14 sessions, EV > REV in 8, p = 0.43).
cen-no-reactivation¶
There are no reactivations between the central amygdala nuclei and the hippocampus.
A claim of absence, judged on the paper's test: EV and REV (M4) over all hippocampus-central nuclei pairs, two-sided Wilcoxon signed rank of EV against REV in NREM and in REM (REM sessions with more than 180 s in each sleep epoch).
print(ev[["session", "rat", "cen_ev", "cen_rev", "cen_pairs", "cen_rem_ev", "cen_rem_rev"]].dropna(subset=["cen_ev"]).round(3).to_string(index=False))
d_cen, diff_cen = ev_test(ev, "cen", "Hippocampus-central nuclei, all cells, NREM")
d_cen_rem, diff_cen_rem = ev_test(ev, "cen_rem", "Hippocampus-central nuclei, all cells, REM")
session rat cen_ev cen_rev cen_pairs cen_rem_ev cen_rem_rev
Rat08-20130708 Rat08 0.088 0.011 567.0 0.289 0.000
Rat08-20130709 Rat08 2.232 1.372 850.0 1.418 1.210
Rat08-20130710 Rat08 6.853 1.024 493.0 0.271 0.168
Rat08-20130711 Rat08 0.564 2.978 396.0 0.805 2.477
Rat10-20140624 Rat10 4.096 0.044 928.0 0.070 0.003
Rat10-20140626 Rat10 1.067 4.929 300.0 0.676 0.136
Rat11-20150313 Rat11 1.279 0.099 1984.0 0.006 0.047
Rat11-20150315 Rat11 0.405 0.308 1254.0 0.045 0.000
Rat11-20150316 Rat11 0.958 0.355 2279.0 NaN NaN
Hippocampus-central nuclei, all cells, NREM: 9 sessions {'Rat08': 4, 'Rat10': 2, 'Rat11': 3}; EV median 1.067%, REV median 0.355%; EV - REV mean 0.714% (s.e.m. 0.979), median 0.603%, EV > REV in 7; Wilcoxon two-sided p = 0.3
Rat08: 4 sessions, EV - REV median 0.468%, EV > REV in 3
Rat10: 2 sessions, EV - REV median 0.095%, EV > REV in 1
Rat11: 3 sessions, EV - REV median 0.603%, EV > REV in 3
Hippocampus-central nuclei, all cells, REM: 8 sessions {'Rat08': 4, 'Rat10': 2, 'Rat11': 2}; EV median 0.280%, REV median 0.091%; EV - REV mean -0.058% (s.e.m. 0.239), median 0.085%, EV > REV in 6; Wilcoxon two-sided p = 0.25
Rat08: 4 sessions, EV - REV median 0.156%, EV > REV in 3
Rat10: 2 sessions, EV - REV median 0.304%, EV > REV in 2
Rat11: 2 sessions, EV - REV median 0.002%, EV > REV in 1
# Context for the claim of absence (decides nothing): the effect size against the paper's, a
# bootstrap interval of the median EV - REV, the shift these sessions would detect, and the power
# to detect a reactivation as large as the hippocampus-BLA one of reactivation-all-cell-pairs.
rng_cen = np.random.default_rng(6101)
for label, x, paper_z in (("NREM", diff_cen.to_numpy(), 0.36), ("REM", diff_cen_rem.to_numpy(), -0.40)):
n = len(x)
lo, hi = hpcbla.median_ci(x, rng_cen)
p = stats.wilcoxon(x).pvalue
bla_like = float((stats.wilcoxon(rng_cen.choice(diff_all.to_numpy(), size=(2000, n)), axis=1).pvalue < 0.05).mean())
print(f"{label}: {n} sessions; effect size r here {hpcbla.effect_r(p, n, np.median(x)):.2f}, the paper {paper_z / np.sqrt(19):.2f} (19 sessions);"
f" median EV - REV {np.median(x):.3f}%, bootstrap 95% interval {lo:.3f} to {hi:.3f}%;"
f" shift detected with 80% power {hpcbla.min_detectable_shift(x, n, rng_cen):.3f}%;"
f" power to detect the hippocampus-BLA EV - REV (median {diff_all.median():.3f}%, {len(diff_all)} sessions resampled): {bla_like:.2f}")
NREM: 9 sessions; effect size r here 0.34, the paper 0.08 (19 sessions); median EV - REV 0.603%, bootstrap 95% interval -2.414 to 4.052%; shift detected with 80% power 3.030%; power to detect the hippocampus-BLA EV - REV (median 1.960%, 15 sessions resampled): 0.92
REM: 8 sessions; effect size r here 0.41, the paper -0.09 (19 sessions); median EV - REV 0.085%, bootstrap 95% interval -0.041 to 0.289%; shift detected with 80% power 0.950%; power to detect the hippocampus-BLA EV - REV (median 1.960%, 15 sessions resampled): 0.93
Result: reproduced. Hippocampus-central nuclei EV does not differ from REV in NREM (9 sessions, median EV - REV 0.60%, EV > REV in 7, p = 0.30) or in REM (8 sessions, median 0.085%, EV > REV in 6, p = 0.25), so the paper's test finds no reactivation. The paper had 19 sessions in each test; with 9 and 8 the test has less power, and the NREM median is positive. As context: the effect sizes here (r 0.34 in NREM, 0.41 in REM) are larger than the paper's (0.08 and -0.09); with 9 sessions a shift of the median EV - REV of about 3.0% (NREM) or 0.95% (REM) would be detected with 80% power, and the bootstrap 95% interval of the NREM median (-2.4 to 4.1%) includes the hippocampus-BLA median of 1.96%. Resampling the hippocampus-BLA sessions at these sizes, a reactivation distributed as the BLA's would reach p < 0.05 in 92% (NREM) and 93% (REM) of samples, so such a reactivation would usually be seen; yet the interval of the NREM median does not exclude one of the BLA's median size. The absence holds on the paper's test and is weak evidence.
intra-structure-reactivation (blocked)¶
Within the hippocampus and within the BLA, pair correlations are reactivated in NREM but not in REM sleep.
The criterion needs all four NREM tests, two of them on putative pyramidal pairs, so it cannot be decided. The all-cell tests (distinct pairs within one structure, M4) are shown and decide nothing.
rows = []
for s in sessions:
r = use[use.session == s["session"]]
if r.empty or not r.sleep.iloc[0]:
continue
for name, labels_ in (("hippocampus", hpcbla.HPC_LABELS), ("BLA", hpcbla.BLA_LABELS)):
ids = hpcbla.units_of(s, labels_)
if len(ids) < hpcbla.MIN_CELLS:
continue
row = {"session": s["session"], "rat": s["rat"], "structure": name, "cells": len(ids)}
res = hpcbla.ev_session(s, ids, None, hpcbla.session_epoch_rows(s))
row.update({"nrem_ev": res["ev"], "nrem_rev": res["rev"]})
if bool(r.rem_ok.iloc[0]):
res = hpcbla.ev_session(s, ids, None, hpcbla.session_epoch_rows(s, "REM"))
row.update({"rem_ev": res["ev"], "rem_rev": res["rev"]})
rows.append(row)
intra = pd.DataFrame(rows)
for name, g in intra.groupby("structure"):
for state in ("nrem", "rem"):
d = g.dropna(subset=[f"{state}_ev"])
x = d[f"{state}_ev"] - d[f"{state}_rev"]
print(f"{name}, all pairs, {state.upper()}: {len(d)} sessions {d.rat.value_counts().sort_index().to_dict()},"
f" EV - REV median {x.median():.3f}%, EV > REV in {int((x > 0).sum())}, Wilcoxon p = {stats.wilcoxon(x).pvalue:.2g}")
BLA, all pairs, NREM: 16 sessions {'Rat08': 9, 'Rat09': 1, 'Rat10': 3, 'Rat11': 3}, EV - REV median 1.443%, EV > REV in 9, Wilcoxon p = 0.14
BLA, all pairs, REM: 15 sessions {'Rat08': 8, 'Rat09': 1, 'Rat10': 3, 'Rat11': 3}, EV - REV median 0.892%, EV > REV in 9, Wilcoxon p = 0.21
hippocampus, all pairs, NREM: 23 sessions {'Rat08': 9, 'Rat10': 7, 'Rat11': 7}, EV - REV median 6.701%, EV > REV in 19, Wilcoxon p = 0.0054
hippocampus, all pairs, REM: 19 sessions {'Rat08': 8, 'Rat10': 6, 'Rat11': 5}, EV - REV median 0.825%, EV > REV in 12, Wilcoxon p = 0.83
# Variant (description only): within the BLA, the identified excitatory units (the five typed
# Rat08 files), sessions with at least 6 of them.
rows = []
for s in typed:
r = use[use.session == s["session"]]
ids = hpcbla.units_of(s, hpcbla.BLA_LABELS, "excitatory")
row = {"session": s["session"], "excitatory_bla": len(ids), "sleep": bool(r.sleep.iloc[0])}
if len(ids) >= hpcbla.MIN_CELLS and row["sleep"]:
res = hpcbla.ev_session(s, ids, None, hpcbla.session_epoch_rows(s))
row.update({"nrem_ev": res["ev"], "nrem_rev": res["rev"], "pairs": res["pairs"]})
rows.append(row)
intra_pyr = pd.DataFrame(rows)
print(intra_pyr.round(3).to_string(index=False))
d = intra_pyr.dropna(subset=["nrem_ev"])
x = d.nrem_ev - d.nrem_rev
print(f"BLA identified excitatory pairs, NREM: {len(d)} sessions, EV - REV median {x.median():.3f}%, EV > REV in {int((x > 0).sum())},"
f" Wilcoxon p = {stats.wilcoxon(x).pvalue:.2g} (the smallest two-sided p with {len(d)} sessions is {2 / 2 ** len(d):.3g})")
print("The hippocampal pyramidal test cannot be formed: at most", int(hpc_exc.max()), "identified excitatory hippocampal units in a file")
session excitatory_bla sleep nrem_ev nrem_rev pairs Rat08-20130709 34 True 5.763 0.207 561 Rat08-20130710 38 True 0.722 0.158 703 Rat08-20130711 47 True 1.753 0.676 1081 Rat08-20130712 11 True 31.234 3.202 55 Rat08-20130713 22 True 0.886 7.405 231 BLA identified excitatory pairs, NREM: 5 sessions, EV - REV median 1.077%, EV > REV in 4, Wilcoxon p = 0.44 (the smallest two-sided p with 5 sessions is 0.0625) The hippocampal pyramidal test cannot be formed: at most 4 identified excitatory hippocampal units in a file
Result: blocked (metadata insufficient). The printed tests are a description only: with all pairs, the hippocampus reactivates in NREM (23 sessions, EV > REV in 19, p = 0.0054) and the BLA does not reach p < 0.05 (16 sessions, 9, p = 0.14); in REM neither does (p = 0.83 and 0.21). Since the criterion needs all four NREM tests and the BLA all-pairs test already fails, the claim could at best be partially reproduced even with putative types. A BLA-only variant on the identified excitatory units of the five typed files (all five have at least six) gives EV above REV in 4 of 5 sessions, p = 0.44; with five sessions the two-sided test cannot reach p < 0.05 (smallest p 0.0625), and the hippocampal pyramidal test cannot be formed (at most 4 identified excitatory hippocampal units in a file).
bla-ripple-up-subgroup, bla-ripple-down-subgroup (blocked)¶
A subset of BLA neurons is positively modulated during hippocampal SPW-Rs, 25.8% of interneurons (42 of 163) and 11.1% of pyramidal cells (137 of 1,233).
Both claims give fractions of putative pyramidal cells and interneurons, which the dandiset does
not type, and are blocked. The Poisson test of M6 (both readings of O2) is run on every BLA unit of
the sessions with states and ripples, and its fractions are shown over all BLA units and over the
identified units of the five typed files, as a description only. The ripple modulation is also
used in the description of ripple-modulated-cells-larger-reactivation, below.
rows = []
for s in sessions:
if s["units"] is None or s["states"] is None or s["ripples"] is None:
continue
for r in hpcbla.ripple_modulation(s, hpcbla.units_of(s, hpcbla.BLA_LABELS)):
rows.append({"session": s["session"], "rat": s["rat"], "cell_type": s["units"]["cell_type"][r["unit"]] or "(no column)", **r})
mod = pd.DataFrame(rows)
print("BLA units tested:", len(mod), "in", mod.session.nunique(), "sessions; untested (no spike in ripples or baseline):", int((~mod.tested).sum()))
summary = mod.groupby("cell_type")[["up_a", "down_a", "up_b", "down_b"]].agg(["sum", "mean"]).round(3)
summary.insert(0, "units", mod.groupby("cell_type").size())
print(summary.to_string())
print("All BLA units: up", int(mod.up_a.sum()), f"({mod.up_a.mean():.3f}) reading a,", int(mod.up_b.sum()), f"({mod.up_b.mean():.3f}) reading b;",
"down", int(mod.down_a.sum()), f"({mod.down_a.mean():.3f}) reading a,", int(mod.down_b.sum()), f"({mod.down_b.mean():.3f}) reading b")
print(mod.groupby("rat")[["up_a", "down_a"]].agg(["sum", "mean"]).round(3).to_string())
BLA units tested: 822 in 17 sessions; untested (no spike in ripples or baseline): 0
units up_a down_a up_b down_b
sum mean sum mean sum mean sum mean
cell_type
(no column) 431 65 0.151 28 0.065 69 0.160 30 0.070
excitatory 152 9 0.059 23 0.151 10 0.066 26 0.171
inhibitory 36 8 0.222 17 0.472 8 0.222 17 0.472
unclassified 203 23 0.113 17 0.084 23 0.113 18 0.089
All BLA units: up 105 (0.128) reading a, 110 (0.134) reading b; down 85 (0.103) reading a, 91 (0.111) reading b
up_a down_a
sum mean sum mean
rat
Rat08 58 0.117 62 0.125
Rat10 5 0.060 12 0.143
Rat11 42 0.174 11 0.045
Result: blocked (metadata insufficient). bla-ripple-down-subgroup depends on this claim and
is blocked with it. The printed fractions are a description only: over all 822 BLA units of 17
sessions, 105 (12.8%) are upmodulated and 85 (10.3%) downmodulated under O2 reading a (110 and 91
under reading b); among the identified units of five files, 8 of 36 inhibitory (22.2%) and 9 of
152 excitatory units (5.9%) are upmodulated.
ripple-modulated-cells-larger-reactivation¶
Reactivation (EV - REV) is larger for hippocampus-BLA pairs whose BLA cell is ripple-modulated than for pairs with a non-modulated BLA cell.
The modulated cells are taken per M6, over every BLA unit, with no cell type. The pyramidal test needs the putative types and cannot be run, so, as the criterion now says (answered, questions/000061.yaml), the all-cell test is judged alone: reproduced in part when it holds, not reproduced when it does not. Per session, EV - REV with the modulated BLA units (up or down, M6) and with the non-modulated ones, sessions with more than 5 hippocampal and more than 5 modulated BLA units; two-sided Wilcoxon rank sum test, and the signed rank test on the same sessions as the criterion's variant. Under both O2 readings.
rows = []
for s in sessions:
r = use[use.session == s["session"]]
if r.empty or not r.sleep.iloc[0] or s["ripples"] is None:
continue
hpc = hpcbla.units_of(s, hpcbla.HPC_LABELS)
m = mod[mod.session == s["session"]]
for reading in ("a", "b"):
is_mod = m[f"up_{reading}"] | m[f"down_{reading}"]
modulated, other = [int(u) for u in m.unit[is_mod]], [int(u) for u in m.unit[~is_mod]]
if len(hpc) <= 5 or len(modulated) <= 5:
continue
rows_ = hpcbla.session_epoch_rows(s)
a = hpcbla.ev_session(s, hpc, modulated, rows_)
b = hpcbla.ev_session(s, hpc, other, rows_)
rows.append({"session": s["session"], "rat": s["rat"], "reading": reading, "modulated": len(modulated), "other": len(other),
"mod_ev_rev": a["ev"] - a["rev"], "non_ev_rev": b["ev"] - b["rev"]})
modev = pd.DataFrame(rows)
print(modev.round(3).to_string(index=False))
for reading, d in modev.groupby("reading"):
d = d.dropna()
print(f"O2 reading {reading}: {len(d)} sessions; EV - REV mean modulated {d.mod_ev_rev.mean():.2f}%, non-modulated {d.non_ev_rev.mean():.2f}%;"
f" rank sum p = {stats.ranksums(d.mod_ev_rev, d.non_ev_rev).pvalue:.2g}; signed rank p = {stats.wilcoxon(d.mod_ev_rev, d.non_ev_rev).pvalue:.2g};"
f" modulated larger in {int((d.mod_ev_rev > d.non_ev_rev).sum())}")
print(" expected (all cells): modulated 6.07 +/- 1.3%, non-modulated 1.76 +/- 0.56%")
for rat, g in d.groupby("rat"):
print(f" {rat}: {len(g)} sessions, mean modulated {g.mod_ev_rev.mean():.2f}%, non-modulated {g.non_ev_rev.mean():.2f}%,"
f" modulated larger in {int((g.mod_ev_rev > g.non_ev_rev).sum())}")
n_eligible = int((use.sleep & use.ripples).sum())
print(f"Sessions entering: {modev.session.nunique()} of the {n_eligible} with sleep and ripples; the other"
f" {n_eligible - modev.session.nunique()} have 5 or fewer hippocampal or modulated BLA units")
session rat reading modulated other mod_ev_rev non_ev_rev Rat08-20130709 Rat08 a 30 48 8.321 1.795 Rat08-20130709 Rat08 b 32 46 8.190 1.850 Rat08-20130710 Rat08 a 18 69 3.845 4.764 Rat08-20130710 Rat08 b 18 69 3.845 4.764 Rat08-20130711 Rat08 a 26 66 9.191 2.752 Rat08-20130711 Rat08 b 29 63 9.137 2.719 Rat08-20130712 Rat08 a 10 30 1.374 0.896 Rat08-20130712 Rat08 b 10 30 1.374 0.896 Rat08-20130713 Rat08 a 13 81 -0.017 1.395 Rat08-20130713 Rat08 b 13 81 -0.017 1.395 Rat08-20130715 Rat08 a 6 33 19.363 2.363 Rat08-20130715 Rat08 b 6 33 19.363 2.363 Rat08-20130716 Rat08 a 10 24 3.628 7.416 Rat08-20130716 Rat08 b 10 24 3.628 7.416 Rat10-20140624 Rat10 a 9 28 -0.098 -0.016 Rat10-20140624 Rat10 b 9 28 -0.098 -0.016 Rat10-20140626 Rat10 a 6 27 -6.249 -3.374 Rat10-20140626 Rat10 b 6 27 -6.249 -3.374 Rat11-20150325 Rat11 a 14 72 6.219 2.214 Rat11-20150325 Rat11 b 17 69 5.955 1.773 Rat11-20150326 Rat11 a 21 51 1.526 -0.089 Rat11-20150326 Rat11 b 21 51 1.526 -0.089 Rat11-20150328 Rat11 a 18 66 3.017 0.048 Rat11-20150328 Rat11 b 20 64 2.637 0.057 O2 reading a: 12 sessions; EV - REV mean modulated 4.18%, non-modulated 1.68%; rank sum p = 0.27; signed rank p = 0.2; modulated larger in 7 expected (all cells): modulated 6.07 +/- 1.3%, non-modulated 1.76 +/- 0.56% Rat08: 7 sessions, mean modulated 6.53%, non-modulated 3.05%, modulated larger in 4 Rat10: 2 sessions, mean modulated -3.17%, non-modulated -1.70%, modulated larger in 0 Rat11: 3 sessions, mean modulated 3.59%, non-modulated 0.72%, modulated larger in 3 O2 reading b: 12 sessions; EV - REV mean modulated 4.11%, non-modulated 1.65%; rank sum p = 0.3; signed rank p = 0.23; modulated larger in 7 expected (all cells): modulated 6.07 +/- 1.3%, non-modulated 1.76 +/- 0.56% Rat08: 7 sessions, mean modulated 6.50%, non-modulated 3.06%, modulated larger in 4 Rat10: 2 sessions, mean modulated -3.17%, non-modulated -1.70%, modulated larger in 0 Rat11: 3 sessions, mean modulated 3.37%, non-modulated 0.58%, modulated larger in 3 Sessions entering: 12 of the 23 with sleep and ripples; the other 11 have 5 or fewer hippocampal or modulated BLA units
Result: not reproduced. With all cells, EV - REV is larger with the modulated BLA units in 7 of 12 sessions (means 4.18 and 1.68%), and the rank sum test does not reach p < 0.05 (p = 0.27 under O2 reading a, 0.3 under reading b). The pyramidal test cannot be run, so the all-cell test decides, and it does not hold. The non-modulated mean is within the expected 1.76 +/- 0.56%; the modulated mean is below the expected 6.07 +/- 1.3%. By rat, the modulated units give the larger EV - REV in 4 of 7 Rat08 sessions, 0 of 2 Rat10 sessions and 3 of 3 Rat11 sessions (reading a). 12 of the 23 sessions with sleep and ripples enter; the paper had 19.
run-correlated-pair-fractions, reactivation-involves-ripple-up-cells, pair-group-ev, contributing-up-cells-gain-increase, ev-independent-of-firing-rate, contributing-pairs-run-sleep-coupling, contributing-quartiles-run-sleep-coupling (blocked)¶
Among nine groups of hippocampus-BLA pyramidal pairs formed by run correlation and BLA ripple modulation, the pairs with a significant positive run correlation and an upmodulated BLA partner show the largest pre-NREM to post-NREM increase in correlation, and differ from every other group.
Each of these criteria pools hippocampus-BLA pyramidal pairs (putative types, M1), and all but
ev-independent-of-firing-rate also use the ripple modulation of bla-ripple-up-subgroup. The
dandiset has no putative types, so all seven are blocked. The earlier attempt tested
reactivation-involves-ripple-up-cells with joint-pattern weights over all BLA units, a measure the
claims file no longer uses.
Result: blocked (metadata insufficient). No analysis is run.
reactivation-peaks-at-ripples¶
Reactivation strength of the whole-run template around hippocampal ripple peaks is higher in post-NREM than in pre-NREM, for pyramidal pairs and for all cell pairs.
The first test uses pyramidal cells and cannot be run. As answered (questions/000061.yaml), the
claim no longer depends on hpc-bla-reactivation-post-sleep: sessions enter by the M7 minimum of 6
cells per structure for all cells, and the criterion judges the all-cell test alone, partially
reproduced when it has p < 0.05. R per M7 with the whole run as template, the ripple-peak value of
pre- and post-NREM per session, one-tailed Wilcoxon signed rank test (post greater than pre).
airpuff-trajectory-ripple-reactivation, safe-trajectory-no-ripple-reactivation and
aversive-direction-stronger test pyramidal cells and are blocked; direction-all-cell-pairs
repeats their tests with all cells, with the air puff and safe templates (laps per O1 reading c),
and is tested below from the same computation. A lap of reading c runs from the reward at one end
to the reward at the other, so it includes the time on the starting platform after the reward,
which the M8 laps (position inside a 150-cm window) leave out: the air puff and safe templates hold
some awake activity on the platforms. The whole-run template includes the platforms in the paper
too (M2).
started = time.time()
rows, curves = [], {}
for s in sessions:
r = use[use.session == s["session"]]
if r.empty or not r.hpc_bla.iloc[0] or s["ripples"] is None:
continue
_, _, run = hpcbla.sleep_epochs(s)
hpc, bla = hpcbla.units_of(s, hpcbla.HPC_LABELS), hpcbla.units_of(s, hpcbla.BLA_LABELS)
for template, iv in (("run", np.array([run])), ("airpuff", hpcbla.lap_intervals(s, "puff")), ("safe", hpcbla.lap_intervals(s, "safe"))):
if len(iv) == 0:
continue
out = hpcbla.reactivation_strength(s, hpc, bla, iv)
curves[(s["session"], template)] = (out["pre_curve"], out["post_curve"])
rows.append({"session": s["session"], "rat": s["rat"], "template": template, "hpc_total": len(hpc), "bla_total": len(bla),
**{k: v for k, v in out.items() if not k.endswith("curve")}})
rs = pd.DataFrame(rows)
rs["post_minus_pre"] = rs.post_peak - rs.pre_peak
print(f"computed in {time.time() - started:.0f} s")
print(rs.round(4).to_string(index=False))
def one_tailed(d, label):
res = stats.wilcoxon(d.post_peak, d.pre_peak, alternative="greater")
print(f"{label}: {len(d)} sessions {d.rat.value_counts().sort_index().to_dict()}; post minus pre median {d.post_minus_pre.median():.4f},"
f" positive in {int((d.post_minus_pre > 0).sum())}; one-tailed Wilcoxon p = {res.pvalue:.2g}")
by_rat = d.groupby("rat").post_minus_pre.agg(sessions="size", median="median", positive=lambda v: int((v > 0).sum()))
print(by_rat.round(4).to_string())
one_tailed(rs[rs.template == "run"], "Whole-run template, all cells")
run_r = rs[rs.template == "run"]
print("Units left out of R (silent in the template or in the NREM of a sleep epoch), per session:")
print(run_r.assign(hpc_left_out=run_r.hpc_total - run_r.hpc_used, bla_left_out=run_r.bla_total - run_r.bla_used)
[["session", "hpc_total", "hpc_used", "hpc_left_out", "bla_total", "bla_used", "bla_left_out"]].to_string(index=False))
centres = hpcbla.peri_centres()
mean_pre = np.mean([curves[k][0] for k in curves if k[1] == "run"], axis=0)
mean_post = np.mean([curves[k][1] for k in curves if k[1] == "run"], axis=0)
fig, ax = plt.subplots(figsize=(5.5, 3.5))
ax.plot(centres, mean_pre, label="pre-NREM")
ax.plot(centres, mean_post, label="post-NREM")
ax.axvspan(-0.25, 0.25, color="0.9")
ax.set_xlabel("time from ripple peak (s)")
ax.set_ylabel("mean z-scored R, all cells")
ax.legend(fontsize=8)
plt.show()
print("Plotted range:", round(float(min(mean_pre.min(), mean_post.min())), 4), "to", round(float(max(mean_pre.max(), mean_post.max())), 4))
computed in 19 s
session rat template hpc_total bla_total template_bins hpc_used bla_used pre_peak pre_ripples post_peak post_ripples post_minus_pre
Rat08-20130708 Rat08 run 9 15 72060 9 15 0.0009 3581 -0.0155 3993 -0.0163
Rat08-20130708 Rat08 airpuff 9 15 39915 9 15 -0.0009 3581 -0.0177 3993 -0.0168
Rat08-20130708 Rat08 safe 9 15 463 7 11 0.0021 3581 0.0022 3993 0.0001
Rat08-20130709 Rat08 run 34 78 72414 34 78 0.0304 3952 0.0432 4919 0.0129
Rat08-20130710 Rat08 run 29 87 72144 29 87 0.0191 3021 0.0293 4222 0.0102
Rat08-20130710 Rat08 airpuff 29 87 28766 29 87 0.0060 3021 0.0375 4222 0.0315
Rat08-20130710 Rat08 safe 29 87 39292 29 86 0.0156 3021 0.0116 4222 -0.0041
Rat08-20130711 Rat08 run 22 92 72246 22 92 0.0218 3152 0.0620 2987 0.0402
Rat08-20130711 Rat08 airpuff 22 92 55135 22 92 0.0094 3152 0.0421 2987 0.0327
Rat08-20130711 Rat08 safe 22 92 11194 22 91 0.0030 3152 0.0056 2987 0.0026
Rat08-20130712 Rat08 run 59 40 66554 59 40 0.0077 2753 0.0376 2241 0.0300
Rat08-20130712 Rat08 airpuff 59 40 24040 59 40 0.0077 2753 0.0310 2241 0.0233
Rat08-20130712 Rat08 safe 59 40 37948 58 40 0.0022 2753 0.0372 2241 0.0350
Rat08-20130713 Rat08 run 60 94 64870 59 94 -0.0123 2743 -0.0168 2494 -0.0046
Rat08-20130713 Rat08 airpuff 60 94 36380 59 94 -0.0117 2743 -0.0041 2494 0.0076
Rat08-20130713 Rat08 safe 60 94 18512 59 92 -0.0184 2743 -0.0131 2494 0.0053
Rat08-20130715 Rat08 run 20 39 54596 20 39 -0.0131 1116 0.0555 3008 0.0686
Rat08-20130715 Rat08 airpuff 20 39 28145 20 39 -0.0147 1116 0.0402 3008 0.0548
Rat08-20130715 Rat08 safe 20 39 16378 20 39 -0.0013 1116 0.0435 3008 0.0448
Rat08-20130716 Rat08 run 17 34 49724 17 34 0.0070 2243 0.0253 4649 0.0183
Rat08-20130716 Rat08 airpuff 17 34 25790 17 34 0.0093 2243 0.0281 4649 0.0188
Rat08-20130716 Rat08 safe 17 34 22282 17 34 0.0043 2243 0.0123 4649 0.0080
Rat08-20130717 Rat08 run 42 17 35940 42 17 0.0029 3451 0.1135 3006 0.1106
Rat08-20130717 Rat08 airpuff 42 17 19427 42 17 0.0139 3451 0.1139 3006 0.0999
Rat08-20130717 Rat08 safe 42 17 15939 42 17 -0.0096 3451 0.0804 3006 0.0900
Rat10-20140624 Rat10 run 29 37 37792 29 37 0.0274 3105 0.0252 2395 -0.0023
Rat10-20140624 Rat10 airpuff 29 37 16528 28 36 0.0342 3105 0.0164 2395 -0.0178
Rat10-20140624 Rat10 safe 29 37 20752 29 37 0.0199 3105 0.0337 2395 0.0138
Rat10-20140626 Rat10 run 15 33 27280 15 33 0.0120 1451 0.0344 1730 0.0224
Rat10-20140626 Rat10 airpuff 15 33 12509 15 31 0.0263 1451 0.0414 1730 0.0151
Rat10-20140626 Rat10 safe 15 33 14409 15 32 -0.0093 1451 0.0113 1730 0.0207
Rat10-20140704 Rat10 run 17 9 27760 17 9 0.0005 1666 0.0367 2661 0.0362
Rat10-20140704 Rat10 airpuff 17 9 11784 17 9 -0.0068 1666 0.0325 2661 0.0394
Rat10-20140704 Rat10 safe 17 9 10710 17 9 -0.0036 1666 0.0168 2661 0.0205
Rat11-20150325 Rat11 run 8 86 16246 8 84 0.0070 1558 0.0525 1291 0.0455
Rat11-20150325 Rat11 airpuff 8 86 7960 8 82 0.0091 1558 0.0414 1291 0.0323
Rat11-20150325 Rat11 safe 8 86 6947 8 82 0.0062 1558 0.0562 1291 0.0499
Rat11-20150326 Rat11 run 27 72 18068 25 72 0.0526 1526 0.0781 1306 0.0255
Rat11-20150326 Rat11 airpuff 27 72 6567 25 71 0.0688 1526 0.0685 1306 -0.0003
Rat11-20150326 Rat11 safe 27 72 10935 25 69 0.0417 1526 0.0763 1306 0.0347
Rat11-20150328 Rat11 run 50 84 20458 50 84 0.2179 819 0.1983 1576 -0.0197
Rat11-20150328 Rat11 airpuff 50 84 12479 50 84 0.1954 819 0.1784 1576 -0.0169
Rat11-20150328 Rat11 safe 50 84 7399 50 81 0.2183 819 0.1930 1576 -0.0253
Whole-run template, all cells: 15 sessions {'Rat08': 9, 'Rat10': 3, 'Rat11': 3}; post minus pre median 0.0224, positive in 11; one-tailed Wilcoxon p = 0.0042
sessions median positive
rat
Rat08 9 0.0183 7
Rat10 3 0.0224 2
Rat11 3 0.0255 2
Units left out of R (silent in the template or in the NREM of a sleep epoch), per session:
session hpc_total hpc_used hpc_left_out bla_total bla_used bla_left_out
Rat08-20130708 9 9 0 15 15 0
Rat08-20130709 34 34 0 78 78 0
Rat08-20130710 29 29 0 87 87 0
Rat08-20130711 22 22 0 92 92 0
Rat08-20130712 59 59 0 40 40 0
Rat08-20130713 60 59 1 94 94 0
Rat08-20130715 20 20 0 39 39 0
Rat08-20130716 17 17 0 34 34 0
Rat08-20130717 42 42 0 17 17 0
Rat10-20140624 29 29 0 37 37 0
Rat10-20140626 15 15 0 33 33 0
Rat10-20140704 17 17 0 9 9 0
Rat11-20150325 8 8 0 86 84 2
Rat11-20150326 27 25 2 72 72 0
Rat11-20150328 50 50 0 84 84 0
Plotted range: -0.0049 to 0.0755
Result: partially reproduced. The pyramidal test cannot be run, and the all-cell test holds: the ripple-peak strength of the whole-run template is higher in post-NREM than in pre-NREM in 11 of 15 sessions (median difference 0.0224, one-tailed p = 0.0042; the paper had p = 0.00189 over 27 sessions). By rat, post exceeds pre in 7 of 9 Rat08, 2 of 3 Rat10 and 2 of 3 Rat11 sessions. R leaves out units silent in the template or in a sleep epoch's NREM: 1 hippocampal unit in Rat08-20130713, 2 in Rat11-20150326 and 2 BLA units in Rat11-20150325, none elsewhere. As the criterion now says, this is partially reproduced. The pyramidal part of the statement is not tested.
print("direction-all-cell-pairs, its three tests with all cells:")
puff_r, safe_r = rs[rs.template == "airpuff"], rs[rs.template == "safe"]
one_tailed(puff_r, "Air puff template, post greater than pre")
one_tailed(safe_r, "Safe template, post greater than pre")
both = puff_r.set_index("session")[["post_minus_pre"]].join(safe_r.set_index("session")[["post_minus_pre"]], lsuffix="_airpuff", rsuffix="_safe").dropna()
res = stats.wilcoxon(both.post_minus_pre_airpuff, both.post_minus_pre_safe, alternative="greater")
print(f"Air puff greater than safe: {len(both)} sessions; median difference {(both.post_minus_pre_airpuff - both.post_minus_pre_safe).median():.4f},"
f" positive in {int((both.post_minus_pre_airpuff > both.post_minus_pre_safe).sum())}; one-tailed Wilcoxon p = {res.pvalue:.2g}")
# Laps of each template (O1 reading c: reward to reward, so each lap includes the time on the
# starting platform after its reward, which the M8 laps leave out).
lap_rows = []
for s in sessions:
if s["session"] not in set(puff_r.session):
continue
row = {"session": s["session"], "rat": s["rat"]}
for kind in ("puff", "safe"):
iv = hpcbla.lap_intervals(s, kind)
d_ = np.diff(iv).ravel()
row.update({f"{kind}_laps": len(iv), f"{kind}_s": round(float(d_.sum())), f"{kind}_lap_median_s": round(float(np.median(d_)), 1),
f"{kind}_lap_max_s": round(float(d_.max()), 1)})
lap_rows.append(row)
lap_table = pd.DataFrame(lap_rows)
print("\nLaps per template:")
print(lap_table.to_string(index=False))
print("Lap duration over all these laps: median", round(float(np.median(np.concatenate([np.diff(hpcbla.lap_intervals(s, k)).ravel()
for s in sessions if s["session"] in set(puff_r.session) for k in ("puff", "safe")]))), 1), "s")
few = lap_table.session[(lap_table.puff_laps < 3) | (lap_table.safe_laps < 3)].tolist()
print("Sessions with fewer than 3 laps in a template:", few)
keep = ~puff_r.session.isin(few)
print("Variant (description only), sessions with at least 3 laps in each template:")
one_tailed(puff_r[keep.to_numpy()], " Air puff template")
one_tailed(safe_r[~safe_r.session.isin(few).to_numpy()], " Safe template")
b3 = both[~both.index.isin(few)]
print(f" Air puff greater than safe: {len(b3)} sessions, positive in {int((b3.post_minus_pre_airpuff > b3.post_minus_pre_safe).sum())};"
f" one-tailed Wilcoxon p = {stats.wilcoxon(b3.post_minus_pre_airpuff, b3.post_minus_pre_safe, alternative='greater').pvalue:.2g}")
direction-all-cell-pairs, its three tests with all cells:
Air puff template, post greater than pre: 14 sessions {'Rat08': 8, 'Rat10': 3, 'Rat11': 3}; post minus pre median 0.0211, positive in 10; one-tailed Wilcoxon p = 0.01
sessions median positive
rat
Rat08 8 0.0274 7
Rat10 3 0.0151 2
Rat11 3 -0.0003 1
Safe template, post greater than pre: 14 sessions {'Rat08': 8, 'Rat10': 3, 'Rat11': 3}; post minus pre median 0.0171, positive in 12; one-tailed Wilcoxon p = 0.0043
sessions median positive
rat
Rat08 8 0.0066 7
Rat10 3 0.0205 3
Rat11 3 0.0347 2
Air puff greater than safe: 14 sessions; median difference 0.0054, positive in 8; one-tailed Wilcoxon p = 0.48
Laps per template:
session rat puff_laps puff_s puff_lap_median_s puff_lap_max_s safe_laps safe_s safe_lap_median_s safe_lap_max_s
Rat08-20130708 Rat08 2 1996 997.9 1976.4 1 23 23.1 23.1
Rat08-20130710 Rat08 6 1438 191.9 504.5 5 1965 465.4 596.8
Rat08-20130711 Rat08 4 2757 660.5 1322.3 4 560 140.1 176.9
Rat08-20130712 Rat08 11 1202 66.7 284.8 11 1897 52.9 1097.1
Rat08-20130713 Rat08 9 1819 115.1 609.5 9 926 79.5 189.5
Rat08-20130715 Rat08 20 1407 56.5 182.1 20 819 29.6 130.4
Rat08-20130716 Rat08 24 1289 46.4 179.8 24 1114 24.2 366.7
Rat08-20130717 Rat08 26 972 28.4 111.9 26 797 18.5 172.6
Rat10-20140624 Rat10 19 826 37.9 97.5 18 1038 51.8 111.9
Rat10-20140626 Rat10 17 625 32.2 75.0 17 721 40.3 89.9
Rat10-20140704 Rat10 25 589 19.2 58.8 25 536 16.7 54.0
Rat11-20150325 Rat11 30 398 10.3 39.8 30 347 9.3 47.0
Rat11-20150326 Rat11 31 328 10.3 18.3 30 547 13.4 77.3
Rat11-20150328 Rat11 32 624 12.7 73.8 31 370 9.6 38.0
Lap duration over all these laps: median 22.6 s
Sessions with fewer than 3 laps in a template: ['Rat08-20130708']
Variant (description only), sessions with at least 3 laps in each template:
Air puff template: 13 sessions {'Rat08': 7, 'Rat10': 3, 'Rat11': 3}; post minus pre median 0.0233, positive in 10; one-tailed Wilcoxon p = 0.0052
sessions median positive
rat
Rat08 7 0.0315 7
Rat10 3 0.0151 2
Rat11 3 -0.0003 1
Safe template: 13 sessions {'Rat08': 7, 'Rat10': 3, 'Rat11': 3}; post minus pre median 0.0205, positive in 11; one-tailed Wilcoxon p = 0.0052
sessions median positive
rat
Rat08 7 0.0080 6
Rat10 3 0.0205 3
Rat11 3 0.0347 2
Air puff greater than safe: 13 sessions, positive in 8; one-tailed Wilcoxon p = 0.37
# Variants (description only) for the platform time and stops that reading c laps hold: laps of
# at most 30 s, and each lap cut to its last 10 s (the run before the arrival reward). Both
# values are the author's, not the claims file's.
LAP_VARIANTS = {"laps of at most 30 s": {"max_s": 30.0}, "last 10 s of each lap": {"last_s": 10.0}}
started = time.time()
var_rows = []
for s in sessions:
if s["session"] not in set(puff_r.session):
continue
hpc, bla = hpcbla.units_of(s, hpcbla.HPC_LABELS), hpcbla.units_of(s, hpcbla.BLA_LABELS)
for variant, kw in LAP_VARIANTS.items():
for template, kind in (("airpuff", "puff"), ("safe", "safe")):
iv = hpcbla.trim_laps(hpcbla.lap_intervals(s, kind), **kw)
if len(iv) == 0:
continue
out = hpcbla.reactivation_strength(s, hpc, bla, iv)
var_rows.append({"variant": variant, "session": s["session"], "rat": s["rat"], "template": template, "laps": len(iv),
"template_s": round(float(np.diff(iv).sum()), 1), "pre_peak": out["pre_peak"], "post_peak": out["post_peak"]})
rv = pd.DataFrame(var_rows)
rv["post_minus_pre"] = rv.post_peak - rv.pre_peak
print(f"computed in {time.time() - started:.0f} s")
print(rv.pivot_table(index="session", columns=["variant", "template"], values="laps").to_string())
for variant in LAP_VARIANTS:
d = rv[rv.variant == variant]
print(f"\n{variant}:")
one_tailed(d[d.template == "airpuff"], " Air puff template")
one_tailed(d[d.template == "safe"], " Safe template")
b = d[d.template == "airpuff"].set_index("session")[["post_minus_pre"]].join(
d[d.template == "safe"].set_index("session")[["post_minus_pre"]], lsuffix="_airpuff", rsuffix="_safe").dropna()
print(f" Air puff greater than safe: {len(b)} sessions, positive in {int((b.post_minus_pre_airpuff > b.post_minus_pre_safe).sum())};"
f" one-tailed Wilcoxon p = {stats.wilcoxon(b.post_minus_pre_airpuff, b.post_minus_pre_safe, alternative='greater').pvalue:.2g}")
computed in 18 s
variant laps of at most 30 s last 10 s of each lap
template airpuff safe airpuff safe
session
Rat08-20130708 1.0 1.0 2.0 1.0
Rat08-20130710 NaN NaN 6.0 5.0
Rat08-20130711 NaN NaN 4.0 4.0
Rat08-20130712 3.0 3.0 11.0 11.0
Rat08-20130713 NaN NaN 9.0 9.0
Rat08-20130715 1.0 10.0 20.0 20.0
Rat08-20130716 7.0 19.0 24.0 24.0
Rat08-20130717 17.0 20.0 26.0 26.0
Rat10-20140624 4.0 1.0 19.0 18.0
Rat10-20140626 8.0 4.0 17.0 17.0
Rat10-20140704 21.0 20.0 25.0 25.0
Rat11-20150325 28.0 29.0 30.0 30.0
Rat11-20150326 31.0 27.0 31.0 30.0
Rat11-20150328 26.0 30.0 32.0 31.0
laps of at most 30 s:
Air puff template: 11 sessions {'Rat08': 5, 'Rat10': 3, 'Rat11': 3}; post minus pre median 0.0103, positive in 8; one-tailed Wilcoxon p = 0.034
sessions median positive
rat
Rat08 5 0.0107 4
Rat10 3 0.0103 3
Rat11 3 -0.0003 1
Safe template: 11 sessions {'Rat08': 5, 'Rat10': 3, 'Rat11': 3}; post minus pre median 0.0302, positive in 8; one-tailed Wilcoxon p = 0.021
sessions median positive
rat
Rat08 5 0.0146 4
Rat10 3 0.0302 2
Rat11 3 0.0517 2
Air puff greater than safe: 11 sessions, positive in 5; one-tailed Wilcoxon p = 0.71
last 10 s of each lap:
Air puff template: 14 sessions {'Rat08': 8, 'Rat10': 3, 'Rat11': 3}; post minus pre median 0.0091, positive in 9; one-tailed Wilcoxon p = 0.034
sessions median positive
rat
Rat08 8 0.0091 5
Rat10 3 0.0186 2
Rat11 3 0.0042 2
Safe template: 14 sessions {'Rat08': 8, 'Rat10': 3, 'Rat11': 3}; post minus pre median 0.0101, positive in 10; one-tailed Wilcoxon p = 0.052
sessions median positive
rat
Rat08 8 0.0023 5
Rat10 3 0.0173 3
Rat11 3 0.0556 2
Air puff greater than safe: 14 sessions, positive in 8; one-tailed Wilcoxon p = 0.36
Result: reproduced. With all cells over 14 sessions, the air puff template is more reactivated in post-NREM than in pre-NREM (10 of 14, one-tailed p = 0.01), so is the safe template (12 of 14, p = 0.0043), and the air puff change does not exceed the safe one (8 of 14, one-tailed p = 0.48). Both within-trajectory tests have p < 0.05 and the between-trajectory test has p >= 0.05, as the criterion requires. By rat, the air puff template is more reactivated in post-NREM in 7 of 8 Rat08, 2 of 3 Rat10 and 1 of 3 Rat11 sessions, and the safe template in 7 of 8, 3 of 3 and 2 of 3, so the air puff result rests mostly on Rat08. Laps are O1 reading c, which as answered decides the claim alone here; the paper had 27 sessions. These laps run from reward to reward, so they include the time on the starting platform after each reward and long stops: the median lap lasts 22.6 s, but in the early Rat08 sessions laps last minutes (Rat08-20130708 has 2 air puff laps of 1996 s in all and one safe lap of 23 s). The templates thus hold awake activity off the track that the M8 laps leave out. Leaving out Rat08-20130708, the only session with fewer than 3 laps in a template, gives the same result (air puff p = 0.0052, safe p = 0.0052, air puff greater than safe p = 0.37, 13 sessions). Long laps are not limited to that session (median air puff laps of 660 s in Rat08-20130711, 192 s in Rat08-20130710 and 46 to 115 s in four other Rat08 sessions), so two further variants, with values of the author's, limit the platform time and stops. With laps of at most 30 s (11 sessions with such laps in both templates) the result is the same: air puff p = 0.034, safe p = 0.021, air puff greater than safe p = 0.71. With each lap cut to its last 10 s (14 sessions) the air puff test holds (9 of 14, p = 0.034) but the safe test does not reach p < 0.05 (10 of 14, p = 0.052), and air puff greater than safe has p = 0.36. The result under the criterion's laps thus depends in part on the time around the run that reading c laps hold. These variants are not the claims file's and decide nothing; the outcome stays as the criterion gives it.
bla-rate-same-across-directions (blocked)¶
Firing rates of BLA cells do not differ significantly between the safe and the air puff trajectories.
The criterion runs two tests, on BLA pyramidal cells and on all BLA cells, under each reading of
the sessions and each reading of the rate (O7): reproduced when neither test has p < 0.05 under
every reading, not reproduced when either test has p < 0.05 under every reading. The pyramidal
test cannot be run (no putative types), but the all-cell test alone could decide "not reproduced"
if it had p < 0.05 under every reading. O7 reading b (the mean of each cell's position
firing-rate curve on each lap type) needs position, which only tracking the raw video would give
(blocked as too large, as the behaviour claims). So the claim is blocked by O7 reading b. The
readings that can be run are shown: per BLA unit, the spike count over the air puff laps and over
the safe laps of the run divided by their durations (O7 reading a), two-sided Wilcoxon signed rank
test, over the cells of the sessions of direction-all-cell-pairs (reading a of the sessions) and
over all BLA cells of sessions with a run and a puff direction (reading b). Laps are O1 reading c,
and so include the time on the starting platform after each reward.
rows = []
for s in sessions:
if s["units"] is None or hpcbla.puff_direction(s)[0] is None:
continue
ivs = {k: hpcbla.lap_intervals(s, k) for k in ("puff", "safe")}
if any(len(v) == 0 for v in ivs.values()):
continue
for i in hpcbla.units_of(s, hpcbla.BLA_LABELS):
t = s["units"]["trains"][i]
rows.append({"session": s["session"], "rat": s["rat"], "unit": i,
**{f"{k}_hz": hpcbla.count_in(t, iv) / np.diff(iv).sum() for k, iv in ivs.items()}})
dir_fr = pd.DataFrame(rows)
for label, d in (("sessions of direction-all-cell-pairs", dir_fr[dir_fr.session.isin(both.index)]), ("all sessions with a puff direction", dir_fr)):
x = d.puff_hz - d.safe_hz
print(f"{label}: {len(d)} BLA cells in {d.session.nunique()} sessions; median air puff minus safe {x.median():.4f} Hz,"
f" Wilcoxon p = {stats.wilcoxon(d.puff_hz, d.safe_hz).pvalue:.2g}")
# Variants (description only): the same rates over the trimmed laps of the cell above.
rows = []
for s in sessions:
if s["units"] is None or hpcbla.puff_direction(s)[0] is None:
continue
for variant, kw in LAP_VARIANTS.items():
ivs = {k: hpcbla.trim_laps(hpcbla.lap_intervals(s, k), **kw) for k in ("puff", "safe")}
if any(len(v) == 0 for v in ivs.values()):
continue
for i in hpcbla.units_of(s, hpcbla.BLA_LABELS):
t = s["units"]["trains"][i]
rows.append({"variant": variant, "session": s["session"], "unit": i,
**{f"{k}_hz": hpcbla.count_in(t, iv) / np.diff(iv).sum() for k, iv in ivs.items()}})
dir_fr_var = pd.DataFrame(rows)
for variant in LAP_VARIANTS:
v = dir_fr_var[dir_fr_var.variant == variant]
for label, d in (("sessions of direction-all-cell-pairs", v[v.session.isin(both.index)]), ("all sessions with a puff direction", v)):
x = d.puff_hz - d.safe_hz
print(f"{variant}, {label}: {len(d)} BLA cells in {d.session.nunique()} sessions; median air puff minus safe {x.median():.4f} Hz,"
f" Wilcoxon p = {stats.wilcoxon(d.puff_hz, d.safe_hz).pvalue:.2g}")
sessions of direction-all-cell-pairs: 739 BLA cells in 14 sessions; median air puff minus safe 0.0036 Hz, Wilcoxon p = 0.022 all sessions with a puff direction: 1033 BLA cells in 21 sessions; median air puff minus safe 0.0032 Hz, Wilcoxon p = 0.0048 laps of at most 30 s, sessions of direction-all-cell-pairs: 466 BLA cells in 11 sessions; median air puff minus safe 0.0050 Hz, Wilcoxon p = 0.016 laps of at most 30 s, all sessions with a puff direction: 760 BLA cells in 18 sessions; median air puff minus safe 0.0030 Hz, Wilcoxon p = 0.0043 last 10 s of each lap, sessions of direction-all-cell-pairs: 739 BLA cells in 14 sessions; median air puff minus safe 0.0000 Hz, Wilcoxon p = 0.75 last 10 s of each lap, all sessions with a puff direction: 1033 BLA cells in 21 sessions; median air puff minus safe 0.0000 Hz, Wilcoxon p = 0.82
Result: blocked (too large). O7 reading b needs position, which only tracking the raw video would give. Under the readings that can be run, the all-cell test finds a difference: the rate on air puff laps exceeds that on safe laps by a median of 0.0036 Hz, p = 0.022 (739 cells, 14 sessions), and by 0.0032 Hz, p = 0.0048 (1033 cells, 21 sessions). That is the direction of "not reproduced", which the criterion gives when either test has p < 0.05 under every reading; only O7 reading b keeps the claim from that outcome. The pyramidal test cannot be run (no putative types), but it is not needed for "not reproduced". The difference may come from the time on the platforms and the stops that these laps hold, not from the trajectory: with laps of at most 30 s the difference remains (0.0050 Hz, p = 0.016, 466 cells in 11 sessions; 0.0030 Hz, p = 0.0043, 760 cells in 18 sessions), but with each lap cut to its last 10 s it is gone (median 0.0000 Hz, p = 0.75 and 0.82). These variants are the author's and decide nothing.
ripple-rate-unchanged¶
The occurrence rate of SPW-Rs does not differ between pre-NREM and post-NREM.
A claim of absence, judged on the paper's test: per session with ripple events, ripple peaks per second of pre-NREM and of post-NREM (RippleRatePrePost.m), two-sided Wilcoxon signed rank test.
rows = []
for s in sessions:
pre, post, run = hpcbla.sleep_epochs(s)
if s["states"] is None or s["ripples"] is None or len(pre) == 0 or len(post) == 0:
continue
rows.append({"session": s["session"], "rat": s["rat"], **hpcbla.ripple_rates(s)})
rip_rate = pd.DataFrame(rows)
rip_rate["post_minus_pre"] = rip_rate.post_rate - rip_rate.pre_rate
print(rip_rate.round(3).to_string(index=False))
res = stats.wilcoxon(rip_rate.pre_rate, rip_rate.post_rate)
print(f"\n{len(rip_rate)} sessions {rip_rate.rat.value_counts().sort_index().to_dict()}; median pre {rip_rate.pre_rate.median():.3f} Hz,"
f" post {rip_rate.post_rate.median():.3f} Hz; median post minus pre {rip_rate.post_minus_pre.median():.4f} Hz,"
f" higher in post in {int((rip_rate.post_minus_pre > 0).sum())}; Wilcoxon two-sided p = {res.pvalue:.2g}")
for rat, g in rip_rate.groupby("rat"):
print(f" {rat}: {len(g)} sessions, median post minus pre {g.post_minus_pre.median():.4f} Hz, higher in post in {int((g.post_minus_pre > 0).sum())}")
session rat pre_nrem_s pre_ripples pre_rate post_nrem_s post_ripples post_rate post_minus_pre
Rat08-20130708 Rat08 4564.0 3581 0.785 4502.000 3993 0.887 0.102
Rat08-20130709 Rat08 4638.0 3952 0.852 5376.000 4919 0.915 0.063
Rat08-20130710 Rat08 3956.0 3021 0.764 4851.000 4222 0.870 0.107
Rat08-20130711 Rat08 3496.0 3152 0.902 3201.000 2987 0.933 0.032
Rat08-20130712 Rat08 3281.0 2753 0.839 3125.000 2241 0.717 -0.122
Rat08-20130713 Rat08 4466.0 2743 0.614 3852.000 2494 0.647 0.033
Rat08-20130715 Rat08 1193.0 1116 0.935 2935.000 3008 1.025 0.089
Rat08-20130716 Rat08 2639.0 2243 0.850 4784.000 4649 0.972 0.122
Rat08-20130717 Rat08 4780.0 3451 0.722 3487.000 3006 0.862 0.140
Rat10-20140624 Rat10 5058.0 3105 0.614 4370.000 2395 0.548 -0.066
Rat10-20140626 Rat10 2374.0 1451 0.611 2870.000 1730 0.603 -0.008
Rat10-20140701 Rat10 6092.0 3882 0.637 4787.000 2810 0.587 -0.050
Rat10-20140702 Rat10 3018.0 1695 0.562 2384.888 1000 0.419 -0.142
Rat10-20140704 Rat10 3818.4 1666 0.436 5514.000 2661 0.483 0.046
Rat10-20140705 Rat10 1973.0 1035 0.525 2722.000 1563 0.574 0.050
Rat10-20140708 Rat10 1732.0 1033 0.596 3382.000 2191 0.648 0.051
Rat11-20150313 Rat11 5000.0 1732 0.346 3076.000 1088 0.354 0.007
Rat11-20150315 Rat11 2837.0 1299 0.458 3325.000 1573 0.473 0.015
Rat11-20150316 Rat11 2337.0 1050 0.449 3098.000 906 0.292 -0.157
Rat11-20150325 Rat11 3317.0 1558 0.470 2853.000 1291 0.453 -0.017
Rat11-20150326 Rat11 4604.0 1526 0.331 4050.000 1306 0.322 -0.009
Rat11-20150327 Rat11 2128.0 868 0.408 5508.000 2388 0.434 0.026
Rat11-20150328 Rat11 2222.0 819 0.369 3341.000 1576 0.472 0.103
23 sessions {'Rat08': 9, 'Rat10': 7, 'Rat11': 7}; median pre 0.611 Hz, post 0.587 Hz; median post minus pre 0.0315 Hz, higher in post in 15; Wilcoxon two-sided p = 0.26
Rat08: 9 sessions, median post minus pre 0.0894 Hz, higher in post in 8
Rat10: 7 sessions, median post minus pre -0.0084 Hz, higher in post in 3
Rat11: 7 sessions, median post minus pre 0.0073 Hz, higher in post in 4
# Context for the claim of absence (decides nothing): the effect size against the paper's, a
# bootstrap interval of the median change, and the change this many sessions would detect.
rng_ripple = np.random.default_rng(6100)
x = rip_rate.post_minus_pre.to_numpy()
n = len(x)
lo, hi = hpcbla.median_ci(x, rng_ripple)
pre_median = rip_rate.pre_rate.median()
print(f"Effect size r = z / sqrt(n): here {abs(hpcbla.effect_r(res.pvalue, n, np.median(x))):.2f} ({n} sessions);"
f" the paper {abs(-1.289) / np.sqrt(41):.2f} (z = -1.289, 41 sessions)")
print(f"Median post minus pre {np.median(x):.4f} Hz, bootstrap 95% interval {lo:.4f} to {hi:.4f} Hz"
f" ({100 * lo / pre_median:.1f}% to {100 * hi / pre_median:.1f}% of the median pre-NREM rate {pre_median:.3f} Hz)")
for n_ in (n, 41):
shift = hpcbla.min_detectable_shift(x, n_, rng_ripple)
print(f"Shift of the median detected with 80% power at {n_} sessions: {shift:.4f} Hz ({100 * shift / pre_median:.1f}% of the median rate);"
f" power for a 10% change: {hpcbla.signed_rank_power(x, n_, 0.1 * pre_median, rng_ripple):.2f}")
print(f"\nNotebook run time: {(time.time() - notebook_started) / 60:.1f} min")
Effect size r = z / sqrt(n): here 0.24 (23 sessions); the paper 0.20 (z = -1.289, 41 sessions) Median post minus pre 0.0315 Hz, bootstrap 95% interval -0.0084 to 0.0514 Hz (-1.4% to 8.4% of the median pre-NREM rate 0.611 Hz)
Shift of the median detected with 80% power at 23 sessions: 0.0679 Hz (11.1% of the median rate); power for a 10% change: 0.75 Shift of the median detected with 80% power at 41 sessions: 0.0503 Hz (8.2% of the median rate); power for a 10% change: 0.95 Notebook run time: 7.3 min
Result: reproduced. Over 23 sessions (the paper had 41), the median ripple rate is 0.611 Hz in pre-NREM and 0.587 Hz in post-NREM; post minus pre has a median of 0.0315 Hz, is positive in 15, and the two-sided Wilcoxon p is 0.26. By rat the median change is 0.089 Hz in Rat08 (higher in post in 8 of 9), -0.008 Hz in Rat10 (3 of 7) and 0.007 Hz in Rat11 (4 of 7). As context, the effect size r is 0.24 here against the paper's 0.20, the bootstrap 95% interval of the median change is -0.008 to 0.051 Hz (-1.4% to 8.4% of the median pre-NREM rate), and with 23 sessions a change of the median of about 0.068 Hz (11% of the rate) would be detected with 80% power (0.050 Hz, 8%, with the paper's 41). The data exclude a change larger than about a tenth of the rate, not a smaller one.
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome | Earlier attempt (-01) |
|---|---|---|
unit-counts |
not reproduced (40 of 61 sessions, 4702 of 7,390 units) | no claim |
rats-learn-puff-location |
blocked: too large (no position) | blocked: too large |
each-rat-learns-puff-location |
blocked: too large | no claim |
dz-speed-slower-post-run |
blocked: too large | no claim |
puff-location-no-speed-bias |
blocked: too large | no claim |
habituation-over-days |
reproduced (O1 reading c) | no claim |
bla-pyr-sleep-firing-increase |
reproduced (152 identified cells, one rat) | no claim |
bla-int-sleep-firing |
not reproduced (REM/wake p = 0.020; NREM p = 0.26) | no claim |
hpc-bla-reactivation-post-sleep |
blocked: metadata insufficient | reproduced (another measure) |
reactivation-all-cell-pairs |
reproduced (p = 0.00031, 15 sessions) | no claim |
reactivation-decay |
blocked: metadata insufficient | no claim |
no-rem-reactivation |
blocked: metadata insufficient | no claim |
piriform-weaker-reactivation |
blocked: metadata insufficient | no claim |
cen-no-reactivation |
reproduced (NREM p = 0.30, REM p = 0.25; weak evidence of absence, 9 and 8 sessions) | no claim |
intra-structure-reactivation |
blocked: metadata insufficient | no claim |
bla-ripple-up-subgroup |
blocked: metadata insufficient | reproduced (another test) |
bla-ripple-down-subgroup |
blocked: metadata insufficient | no claim |
ripple-modulated-cells-larger-reactivation |
not reproduced (all cells alone, rank sum p = 0.27, 12 sessions) | no claim |
run-correlated-pair-fractions |
blocked: metadata insufficient | no claim |
reactivation-involves-ripple-up-cells |
blocked: metadata insufficient | reproduced (another measure) |
pair-group-ev |
blocked: metadata insufficient | no claim |
contributing-up-cells-gain-increase |
blocked: metadata insufficient | no claim |
ev-independent-of-firing-rate |
blocked: metadata insufficient | no claim |
reactivation-peaks-at-ripples |
partially reproduced (all cells alone, p = 0.0042, 15 sessions) | reproduced (another measure) |
airpuff-trajectory-ripple-reactivation |
blocked: metadata insufficient | no claim |
safe-trajectory-no-ripple-reactivation |
blocked: metadata insufficient | no claim |
aversive-direction-stronger |
blocked: metadata insufficient | not reproduced (another measure) |
direction-all-cell-pairs |
reproduced (p = 0.01, 0.0043 and 0.48, 14 sessions; with laps cut to their last 10 s the safe test has p = 0.052) | no claim |
bla-rate-same-across-directions |
blocked: too large (O7 reading b needs position; the runnable readings give p < 0.05, but not with laps cut to their last 10 s) | no claim |
ripple-rate-unchanged |
reproduced (p = 0.26, 23 sessions; bootstrap interval of the change -1.4% to 8.4% of the rate) | no claim |
contributing-pairs-run-sleep-coupling |
blocked: metadata insufficient | no claim |
contributing-quartiles-run-sleep-coupling |
blocked: metadata insufficient | no claim |
How the analysis differs from the paper's. Every method is the claims file's, from the paper's
Methods and code. The differences come from the data: the dandiset holds 40 of the paper's 61
sessions; it has no putative cell types, so the cells of bla-pyr-sleep-firing-increase and
bla-int-sleep-firing are the identified excitatory and inhibitory units of five Rat08 files (O4
reading a, which decides alone as answered), and every claim that needs putative pyramidal
cells is blocked, except that reactivation-peaks-at-ripples and
ripple-modulated-cells-larger-reactivation are judged on their all-cell tests alone, as their
criteria now allow; it has no position, so laps are the intervals between rewards at opposite ends
(O1 reading c, which decides alone as answered), which include the time on the starting platform
and the stops that the M8 laps leave out, and the behaviour claims and
bla-rate-same-across-directions (O7 reading b) are blocked; brain regions are the dandiset's labels read as the paper's structures; units 0 and 1 of
ten files are left out of the analyses because their spike times are mixed. Bins belong to an
interval when their centre does, and ripples to NREM when their peak does (standard choices).
What was found about the dandiset.
- The putative cell types the paper uses are missing; a
cell_typecolumn in five of 33 files with units marks only identified excitatory and inhibitory cells and is described in one line. There are no waveforms. A short-latency cross-correlogram check finds connections of the marked sign for nearly every marked BLA unit. - It holds 40 recording days against the paper's 61, and 4702 units against 7,390.
- Still present from the earlier attempt (the files are unchanged): the epochs table is in
milliseconds although it says seconds; video timestamps are in milliseconds, stated as seconds,
and are not on the session clock; in ten files the units index is wrong at units 0 and 1; seven
files have no units table, sixteen no sleep states and fourteen no ripple events; there is no
position tracking; the unit region labels, the
trialstable and the ripple labels are undescribed. - The archive's asset metadata lists ElectricalSeries and ProcessingModule for exactly the files that lack them, and not for those that have them (all 40 files wrong for both).
- The dandiset links no analysis code; the paper links its code.