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. The article is not open access (it is not in PubMed Central, and its HAL record holds metadata only), so this attempt used the abstract alone (PubMed 28892057). Every claim is a sentence of the abstract, none has a reported number, and every method choice was fixed in the claims file before the analysis.
Dandiset. DANDI:000061, version
0.250624.0426. The abstract cites no DANDI version. There are two published versions,
0.210812.1517 and 0.250624.0426; they hold the same 40 files (same blobs and sizes) and differ only
in metadata (contributor affiliations and ORCIDs, and a relatedResource link to the paper, added
in 2025), none of which the analysis uses. The most recent published version is pinned.
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. There is no position tracking; the files hold the raw video frames, from which position is not derived.
Code. The dandiset links no analysis code, and the abstract names none. The parameters are
common choices, listed in the helper module hpcbla.py and the claims file.
| Claim | Tested here |
|---|---|
hpc-bla-reactivation-post-sleep: joint hippocampus and BLA patterns of the run reactivate in post-sleep Non-REM |
tested |
reactivation-peaks-at-ripples: the reactivations peak during sharp wave-ripples |
tested |
bla-ripple-up-subgroup: a subgroup of BLA cells is positively modulated during ripples |
tested |
reactivation-involves-ripple-up-cells: the reactivations involve those BLA cells |
tested |
aversive-direction-stronger: reactivation is stronger for the air-puff direction than the safe one |
tested |
rats-learn-puff-location: rats learned the location of the air puff |
blocked (too large: position could only come from tracking the raw video frames, see the data check) |
In [1]:
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 253 s
Data check¶
Each check of the metadata against the values is stated under its output.
In [2]:
from dandi.dandiapi import DandiAPIClient
with DandiAPIClient() as client:
meta = client.get_dandiset(hpcbla.DANDISET_ID, hpcbla.DANDISET_VERSION).get_raw_metadata()
summary = meta["assetsSummary"]
print("Archive summary: files", summary["numberOfFiles"], "; subjects", summary["numberOfSubjects"],
"; species", [s["name"] for s in summary["species"]], "; measured", summary["variableMeasured"])
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)),
"bla_complex": int(sum(region.get(k, 0) for k in hpcbla.BLA_COMPLEX_LABELS)),
"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()))
Archive summary: files 40 ; subjects 4 ; species ['Brown rat'] ; measured ['Units', 'LFP', 'ElectricalSeries', 'ProcessingModule']
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 bla_complex states ripples puffs rewards pre_sleep_s post_sleep_s run_s lfp raw videos trials
Rat08-20130708 Rat08 166 9 15 60 True 9473 3 9 9890 10664 3603 True True 2 215.0
Rat08-20130709 Rat08 201 34 78 99 True 10713 2 22 9725 11239 3621 True True 3 185.0
Rat08-20130710 Rat08 190 31 87 87 True 9251 9 24 10098 10224 3607 True True 3 113.0
Rat08-20130711 Rat08 201 24 92 92 True 8193 5 36 8903 8425 3612 True True 3 289.0
Rat08-20130712 Rat08 199 59 40 40 True 6118 14 69 10318 9222 3328 True True 3 266.0
Rat08-20130713 Rat08 252 60 94 94 True 6835 20 44 9199 9766 3244 True True 3 183.0
Rat08-20130715 Rat08 166 20 39 39 True 6945 34 82 9266 9334 2730 True False 3 NaN
Rat08-20130716 Rat08 187 19 34 73 True 8934 30 81 10663 11500 2486 True False 3 NaN
Rat08-20130717 Rat08 197 42 17 76 True 7857 31 99 9671 10985 1797 True False 3 NaN
Rat08-20130718 Rat08 219 50 0 51 False 0 49 146 11154 11121 2425 True False 3 NaN
Rat08-20130719 Rat08 157 0 0 48 False 0 29 83 7054 10068 2126 True False 3 NaN
Rat08-20130720 Rat08 182 0 0 7 False 0 38 92 10801 9430 1819 True False 3 NaN
Rat08-20130722 Rat08 0 0 0 0 False 0 40 108 8573 11377 1856 True False 3 NaN
Rat09-20140324 Rat09 136 0 35 66 True 0 20 37 9887 10655 3065 True True 3 137.0
Rat09-20140325 Rat09 97 0 49 55 False 0 37 62 6871 9671 3111 True True 3 167.0
Rat09-20140326 Rat09 140 0 108 108 False 0 78 162 12038 8735 2708 True True 3 276.0
Rat09-20140328 Rat09 75 0 0 0 False 0 54 129 0 9763 2027 True True 2 69.0
Rat09-20140329 Rat09 67 0 67 67 False 0 50 137 8932 8349 1935 True True 3 219.0
Rat09-20140402 Rat09 0 0 0 0 False 0 59 158 9800 10256 2226 True True 3 271.0
Rat09-20140403 Rat09 0 0 0 0 False 0 40 98 8984 7443 1299 True True 3 190.0
Rat09-20140404 Rat09 0 0 0 0 False 0 52 123 8982 10688 1347 True True 3 240.0
Rat09-20140405 Rat09 0 0 0 0 False 0 44 116 5513 8597 908 True True 3 185.0
Rat09-20140407 Rat09 0 0 0 0 False 0 0 28 8248 5062 621 True True 1 71.0
Rat10-20140624 Rat10 101 29 37 40 True 7114 26 79 8141 9106 1890 True False 3 NaN
Rat10-20140626 Rat10 76 15 33 33 True 9302 20 66 8126 7554 1364 True False 3 NaN
Rat10-20140701 Rat10 40 12 6 23 True 9275 14 72 9766 12766 2618 True False 3 NaN
Rat10-20140702 Rat10 61 18 2 36 True 4534 38 162 8551 9870 1333 True False 3 NaN
Rat10-20140704 Rat10 57 17 9 9 True 5503 9 111 8103 9854 1388 True False 3 NaN
Rat10-20140705 Rat10 49 18 0 6 True 5717 53 158 10876 9559 1492 True False 3 NaN
Rat10-20140708 Rat10 56 7 0 3 True 5934 24 152 7925 9390 1353 True False 3 NaN
Rat11-20150313 Rat11 136 32 0 16 True 3724 19 56 11423 9948 1461 True False 3 NaN
Rat11-20150315 Rat11 127 22 0 39 True 4942 30 102 13153 9933 1067 True False 3 NaN
Rat11-20150316 Rat11 190 57 0 90 True 3503 46 90 11336 12113 1246 True False 3 NaN
Rat11-20150325 Rat11 161 8 86 153 True 4219 32 89 14305 10695 812 True False 3 NaN
Rat11-20150326 Rat11 202 27 74 158 True 3823 32 85 11065 10850 903 True False 3 NaN
Rat11-20150327 Rat11 214 68 0 104 True 4446 30 82 11965 11606 1042 True False 3 NaN
Rat11-20150328 Rat11 167 50 84 117 True 4130 37 111 11595 11582 1023 True False 3 NaN
Rat11-20150331 Rat11 156 83 32 68 False 5155 21 71 13324 7614 1025 True False 3 NaN
Rat11-20150401 Rat11 77 18 0 20 False 4384 28 78 10826 10173 941 True False 3 NaN
Rat11-20150403 Rat11 0 0 0 0 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. Seven files have no units table, sixteen have no sleep states, and the Rat09 files have no ripple events and no hippocampal units, so Rat09 cannot enter any analysis that needs both regions. Raw voltage is in 16 of the 40 files; the LFP and video are in all of them.
In [3]:
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("\nUsed as hippocampus:", sorted(hpcbla.HPC_LABELS), "; as BLA:", sorted(hpcbla.BLA_LABELS),
"; wider BLA 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']
Used as hippocampus: ['hip', 'hpc'] ; as BLA: ['bla'] ; wider BLA variant: ['bla', 'blv', 'bmp', 'ladl'] (case ignored)
The unit labels are short codes with no description in the files or the dandiset metadata. Some
appear in two cases (bla and BLA, hpc and HIP, pir and PIR); in Rat08-20130717 the
upper-case labels are on shanks 13 to 20 and the lower-case ones on shanks 1 to 12, so the case
marks a probe, not a duplicate (no two units of that file share a spike count and first spike;
see the duplicate check below). Others (ladl, blv, bmp, den, ven, cen, pir,
plc01, PLCo1, cpu, no label) are not defined anywhere. The electrode table gives every
channel the location unknown. Rat09-20140328 has no location column at all.
In [4]:
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": (min(a for _, a, _ in s["states"]), max(b for _, _, b 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"]]})
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
Rat08-20130708 24759200.0 24759.2 1000.0 24759.2 1250.0 (5.0, 24754.0)
Rat08-20130709 25510500.0 25510.5 1000.0 25510.5 1250.0 (5.0, 25505.0)
Rat08-20130710 25195700.0 25195.7 1000.0 25195.7 1250.0 (5.0, 25190.0)
Rat08-20130711 21867700.0 21867.7 1000.0 21867.7 1250.0 (5.0, 21862.0)
Rat08-20130712 23980425.0 23980.4 1000.0 23980.4 1250.0 (5.0, 23975.0)
Rat08-20130713 22680000.0 22680.0 1000.0 22680.0 1250.0 (5.0, 22675.0)
Rat08-20130715 21827500.0 21827.5 1000.0 21827.5 1250.0 (5.0, 21822.0)
Rat08-20130716 25112700.0 25112.7 1000.0 25112.7 1250.0 (5.0, 25107.0)
Rat08-20130717 23243600.0 23243.6 1000.0 23243.6 1250.0 (5.0, 23238.0)
Rat08-20130718 25651100.0 25651.1 1000.0 25651.1 1250.0 None
Rat08-20130719 19752575.0 19752.6 1000.0 19752.6 1250.0 None
Rat08-20130720 22498700.0 22498.7 1000.0 22498.7 1250.0 None
Rat08-20130722 22295200.0 22295.2 1000.0 NaN 1250.0 None
Rat09-20140324 23930200.0 23930.2 1000.0 23930.2 1250.0 (5.0, 23925.0)
Rat09-20140325 20033500.0 20033.5 1000.0 20033.5 1250.0 None
Rat09-20140326 24306500.0 24306.5 1000.0 24306.5 1250.0 None
Rat09-20140328 12660300.0 12660.3 1000.0 12660.3 1250.0 None
Rat09-20140329 20108100.0 20108.1 1000.0 20108.1 1250.0 None
Rat09-20140402 23074300.0 23074.3 1000.0 NaN 1250.0 None
Rat09-20140403 18526600.0 18526.6 1000.0 NaN 1250.0 None
Rat09-20140404 21850700.0 21850.7 1000.0 NaN 1250.0 None
Rat09-20140405 15987200.0 15987.2 1000.0 NaN 1250.0 None
Rat09-20140407 13930600.0 13930.6 1000.0 NaN 1250.0 None
Rat10-20140624 20018200.0 20018.2 1000.0 20018.2 1250.0 (5.0, 20013.0)
Rat10-20140626 17839200.0 17839.2 1000.0 17839.2 1250.0 (5.0, 17834.0)
Rat10-20140701 26196300.0 26196.3 1000.0 26196.3 1250.0 (5.0, 26191.0)
Rat10-20140702 21061687.5 21061.7 1000.0 21061.7 1250.0 (5.0, 21056.0)
Rat10-20140704 20183000.0 20183.0 1000.0 20183.0 1250.0 (5.0, 20178.0)
Rat10-20140705 22905400.0 22905.4 1000.0 22905.4 1250.0 (5.0, 22900.0)
Rat10-20140708 19519200.0 19519.2 1000.0 19519.1 1250.0 (5.0, 19514.0)
Rat11-20150313 23213900.0 23213.9 1000.0 23213.9 1250.0 (5.0, 23208.0)
Rat11-20150315 24989700.0 24989.7 1000.0 24989.7 1250.0 (5.0, 24984.0)
Rat11-20150316 25062400.0 25062.4 1000.0 25062.4 1250.0 (5.0, 25057.0)
Rat11-20150325 26214300.0 26214.3 1000.0 26214.3 1250.0 (5.0, 26209.0)
Rat11-20150326 23205600.0 23205.6 1000.0 23205.6 1250.0 (5.0, 23200.0)
Rat11-20150327 25025187.5 25025.2 1000.0 25025.2 1250.0 (5.0, 25020.0)
Rat11-20150328 25004100.0 25004.1 1000.0 25004.1 1250.0 (5.0, 24999.0)
Rat11-20150331 22777000.0 22777.0 1000.0 22777.0 1250.0 None
Rat11-20150401 22477800.0 22477.8 1000.0 22477.8 1250.0 None
Rat11-20150403 9613500.0 9613.5 1000.0 NaN 1250.0 None
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')]}
The epoch table says its times are in seconds, but its last stop time is 1000 times the LFP
duration (and the last spike time) 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. Three files tag their sleep epochs only as sleep; there the sleep before the
first run is taken as pre-sleep and the sleep after it as post-sleep. Of these only
Rat08-20130708 enters the analysis (the other two have no units).
In [5]:
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 run videos of the early
Rat08 files span about half of their run epoch). The videos are not used here. 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.
It is not used.
In [6]:
# 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")
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.6 s (0.28 s per frame of 691200 bytes) At that rate the run frames alone would take about 149 hours to stream
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 printed above and running a tracker, which does not fit the time limit; their timestamps would also have to be aligned to the session clock first (see above). This is why the behaviour claim is blocked as too large.
In [7]:
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.
In ten files the spike times of unit 1 are not sorted: they fall back once, always after 65536, 131072, 196608, 327680 or 458752 spikes (multiples of 65536), 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.
In [8]:
# 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 = [i for i in u["trains"] if u["region"][i] in 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 at 0.33 to 1.0 per second of Non-REM and at most 0.19 per second of REM (median 0.022); with the states shifted by 600 s the two rates move toward each other (0.37 and 0.32). 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.
In [9]:
# (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 = [i for i in u["trains"] if u["region"][i] in hpcbla.HPC_LABELS]
laps = hpcbla.run_laps(s)
if len(ids) < 5 or len(laps) < 6:
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-20130709 7 1.000 0.714 Rat08-20130710 13 1.000 0.692 Rat08-20130711 9 0.889 0.222 Rat08-20130712 23 1.000 0.783 Rat08-20130713 20 1.000 0.800 Rat08-20130715 41 0.976 0.756 Rat08-20130716 48 1.000 0.833 Rat08-20130717 54 1.000 0.944 Rat08-20130718 89 0.989 0.876 Rat10-20140624 39 0.949 0.744 Rat10-20140626 34 1.000 0.853 Rat10-20140701 28 0.857 0.821 Rat10-20140702 71 0.986 0.845 Rat10-20140704 52 1.000 0.962 Rat10-20140705 94 0.947 0.883 Rat10-20140708 82 0.988 0.866 Rat11-20150313 39 0.949 0.872 Rat11-20150315 54 0.963 0.852 Rat11-20150316 66 0.955 0.939 Rat11-20150325 60 1.000 0.850 Rat11-20150326 63 1.000 0.968 Rat11-20150327 59 0.983 0.898 Rat11-20150328 65 1.000 0.908 Rat11-20150331 32 1.000 0.938 Rat11-20150401 54 0.926 0.870 Median accuracy 0.989 ; with laps shifted by half a lap 0.853 ; chance 0.5
Hippocampal rates tell the lap direction apart in nearly every lap with the 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.
In [10]:
# (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.
rows = []
for s in sessions:
u = s["units"]
if not u or "puf" not in s["events"] or len(s["events"]["puf"][1]) < 5:
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 = [i for i in u["trains"] if u["region"][i] in labels_]
if len(ids) < 3:
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()))
print("Puffs inside laps of the session's air-puff direction, in these files:", int(sum(puff[d].where(puff.direction == d, 0).sum() for d in ("LR", "RL"))),
"of", int(puff.puffs.sum()))
session puffs direction LR RL none bla_change bla_change_shift30 hpc_change hpc_change_shift30 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 Rat09-20140324 20 RL 1 19 0 0.141 0.084 NaN NaN Rat09-20140325 37 LR 36 1 0 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 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 Puffs inside laps of the session's air-puff direction, in these files: 927 of 957
In the 31 files above (units and at least 5 puffs), nearly all puffs (927 of 957) fall inside laps of one direction in each session, which gives the air-puff direction; in the 15 sessions used for the direction claim, 295 of 304 do (see that section). The puff side changes between days in the same rat. The pooled response to puffs is small in most files (median change 0.04 in BLA and in hippocampus); four files show a BLA increase above 0.3, against one with the shifted times. This is consistent with the puffs sharing the spike clock but does not show it for every file.
In [11]:
# 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; each session is analysed on its own and results are given per session and per rat.
In [12]:
rows = []
for s in sessions:
if s["units"] is None:
continue
hpc, bla = hpcbla.unit_sets(s)
hpc_w, bla_w = hpcbla.unit_sets(s, hpcbla.BLA_COMPLEX_LABELS)
rows.append({"session": s["session"], "rat": s["rat"], "hpc_used": len(hpc), "bla_used": len(bla),
"bla_complex_used": len(bla_w), "usable": hpcbla.usable(s),
"usable_wider_bla": hpcbla.usable(s, hpcbla.BLA_COMPLEX_LABELS),
"puff_direction": hpcbla.puff_direction(s)[0]})
use = pd.DataFrame(rows)
print(use.to_string(index=False))
used = [s for s in sessions if hpcbla.usable(s)]
print("\nSessions used:", len(used), {r: sum(s["rat"] == r for s in used) for r in sorted({s["rat"] for s in used})})
print("Sessions with the wider BLA labels:", int(use.usable_wider_bla.sum()))
print("Units used (hippocampus, BLA):", int(use.hpc_used[use.usable].sum()), int(use.bla_used[use.usable].sum()))
session rat hpc_used bla_used bla_complex_used usable usable_wider_bla puff_direction
Rat08-20130708 Rat08 7 10 39 True True LR
Rat08-20130709 Rat08 26 45 62 True True LR
Rat08-20130710 Rat08 26 45 45 True True RL
Rat08-20130711 Rat08 18 43 43 True True RL
Rat08-20130712 Rat08 39 23 23 True True LR
Rat08-20130713 Rat08 44 63 63 True True LR
Rat08-20130715 Rat08 15 33 33 True True RL
Rat08-20130716 Rat08 15 26 63 True True LR
Rat08-20130717 Rat08 36 9 63 True True RL
Rat08-20130718 Rat08 23 0 50 False False RL
Rat08-20130719 Rat08 0 0 45 False False LR
Rat08-20130720 Rat08 0 0 7 False False LR
Rat09-20140324 Rat09 0 23 48 False False RL
Rat09-20140325 Rat09 0 28 32 False False LR
Rat09-20140326 Rat09 0 57 57 False False RL
Rat09-20140328 Rat09 0 0 0 False False RL
Rat09-20140329 Rat09 0 46 46 False False RL
Rat10-20140624 Rat10 21 19 22 True True RL
Rat10-20140626 Rat10 14 16 16 True True RL
Rat10-20140701 Rat10 11 2 16 False True RL
Rat10-20140702 Rat10 17 0 28 False True LR
Rat10-20140704 Rat10 13 6 6 True True LR
Rat10-20140705 Rat10 16 0 4 False False RL
Rat10-20140708 Rat10 5 0 3 False False RL
Rat11-20150313 Rat11 27 0 9 False True RL
Rat11-20150315 Rat11 20 0 20 False True RL
Rat11-20150316 Rat11 43 0 59 False True RL
Rat11-20150325 Rat11 8 44 79 True True RL
Rat11-20150326 Rat11 20 47 96 True True LR
Rat11-20150327 Rat11 49 0 73 False True RL
Rat11-20150328 Rat11 39 54 80 True True RL
Rat11-20150331 Rat11 70 28 57 False False RL
Rat11-20150401 Rat11 14 0 19 False False LR
Sessions used: 15 {'Rat08': 9, 'Rat10': 3, 'Rat11': 3}
Sessions with the wider BLA labels: 21
Units used (hippocampus, BLA): 341 483
Under the claims file's rule (sleep states, ripples, a sleep epoch before and after the run, and at least five hippocampal and five BLA units at 0.1 Hz or more in the run), 15 sessions from three rats are used: nine of Rat08, three of Rat10 and three of Rat11, with 341 hippocampal and 483 BLA units. With the wider BLA labels, 21 sessions qualify.
In [13]:
# Epoch boundaries (s) of the sessions used: every run epoch, the one used, and the sleep epochs.
rows = []
for s in used:
pre, post, run = hpcbla.sleep_epochs(s)
rows.append({"session": s["session"], "run_epochs": [(round(a), round(b)) for t, a, b in s["epochs"] if t == "run"],
"run_used": (round(run[0]), round(run[1])), "pre_sleep": [(round(a), round(b)) for a, b in pre],
"post_sleep": [(round(a), round(b)) for a, b in post]})
print(pd.DataFrame(rows).to_string(index=False))
session run_epochs run_used pre_sleep post_sleep Rat08-20130708 [(9890, 13493), (24157, 24759)] (9890, 13493) [(0, 9890)] [(13493, 24157)] Rat08-20130709 [(10042, 13663)] (10042, 13663) [(317, 10042)] [(13663, 19954), (19954, 24902)] Rat08-20130710 [(10414, 14021)] (10414, 14021) [(315, 10414)] [(14021, 24245)] Rat08-20130711 [(9212, 12824)] (9212, 12824) [(308, 9212)] [(12824, 21249)] Rat08-20130712 [(10613, 13941)] (10613, 13941) [(295, 10613)] [(13941, 17743), (17743, 23163)] Rat08-20130713 [(9406, 12649)] (9406, 12649) [(206, 9406)] [(12649, 22415)] Rat08-20130715 [(9390, 12120)] (9390, 12120) [(124, 9390)] [(12120, 21454)] Rat08-20130716 [(10813, 13299)] (10813, 13299) [(150, 10813)] [(13299, 24799)] Rat08-20130717 [(9861, 11658)] (9861, 11658) [(190, 9861)] [(11658, 22642)] Rat10-20140624 [(8372, 10262)] (8372, 10262) [(231, 8372)] [(10262, 19368)] Rat10-20140626 [(8317, 9681)] (8317, 9681) [(191, 8317)] [(9681, 15423), (15423, 17236)] Rat10-20140704 [(8289, 9677)] (8289, 9677) [(186, 8289)] [(9677, 19531)] Rat11-20150325 [(14520, 15332)] (14520, 15332) [(215, 14520)] [(15332, 26027)] Rat11-20150326 [(11255, 12158)] (11255, 12158) [(190, 11255)] [(12158, 23008)] Rat11-20150328 [(11782, 12805)] (11782, 12805) [(187, 11782)] [(12805, 24387)]
Each session used has one run epoch, except Rat08-20130708, which has a second run of about 600 s after the post-sleep; there the first run is used, and its post-sleep is the sleep between the two runs. Pre- and post-sleep are the epochs next to the run used.
hpc-bla-reactivation-post-sleep¶
Coordinated hippocampus and BLA activity patterns of the training run are reactivated during non-REM sleep after training, more than during non-REM sleep before it.
Joint patterns are the singular vectors of the hippocampus by BLA cross-correlation matrix of the run (50 ms bins); a pattern is significant when its singular value exceeds the 95th percentile of the largest one over 200 shuffles of the BLA trains. Reactivation strength is (z_hpc . u)(z_bla . v) in each Non-REM bin.
In [14]:
rng = np.random.default_rng(61001)
started = time.time()
main = {s["session"]: hpcbla.reactivation(s, rng) for s in used}
patterns = pd.DataFrame([r for out in main.values() for r in out["rows"]])
patterns["change"] = patterns.post - patterns.pre
sig = patterns[patterns.significant].copy()
print(f"computed in {time.time() - started:.0f} s")
print("Significant patterns:", len(sig), "of", len(patterns), "; sessions with at least one:",
sig.session.nunique(), "of", len(used))
per_session = sig.groupby(["rat", "session"]).agg(patterns=("change", "size"), pre=("pre", "median"),
post=("post", "median"), change=("change", "median"),
positive=("change", lambda x: int((x > 0).sum())))
print(per_session.round(4).to_string())
def wilcoxon(x):
x = np.asarray(x)
return stats.wilcoxon(x).pvalue if len(x) > 0 and np.any(x != 0) else np.nan
print("\nPooled: median post minus pre", round(sig.change.median(), 4), "; positive", int((sig.change > 0).sum()),
"of", len(sig), "; Wilcoxon p =", f"{wilcoxon(sig.change):.2g}")
by_rat = sig.groupby("rat").change.agg(["size", "median", lambda x: int((x > 0).sum()), wilcoxon])
by_rat.columns = ["patterns", "median_change", "positive", "wilcoxon_p"]
print(by_rat.to_string())
print("\nReference: non-significant patterns, median post minus pre",
round(patterns[~patterns.significant].change.median(), 4), "over", int((~patterns.significant).sum()))
print("Session medians above zero:", int((per_session.change > 0).sum()), "of", len(per_session),
"; Wilcoxon over session medians p =", f"{wilcoxon(per_session.change):.2g}")
computed in 26 s
Significant patterns: 36 of 289 ; sessions with at least one: 15 of 15
patterns pre post change positive
rat session
Rat08 Rat08-20130708 2 0.0031 0.0085 0.0054 1
Rat08-20130709 2 0.0554 0.0825 0.0271 2
Rat08-20130710 2 0.0573 0.0626 0.0053 1
Rat08-20130711 3 0.0143 0.0364 0.0305 3
Rat08-20130712 1 0.0124 0.0322 0.0198 1
Rat08-20130713 4 0.0225 0.0291 0.0142 3
Rat08-20130715 3 0.0026 0.0308 0.0188 3
Rat08-20130716 2 0.0239 0.0438 0.0199 2
Rat08-20130717 2 0.0234 0.0642 0.0409 2
Rat10 Rat10-20140624 3 0.0066 0.0036 0.0032 2
Rat10-20140626 4 0.0093 0.0081 -0.0011 2
Rat10-20140704 1 0.0147 0.0241 0.0093 1
Rat11 Rat11-20150325 2 0.0326 0.0560 0.0234 1
Rat11-20150326 2 0.0365 0.0500 0.0135 2
Rat11-20150328 3 0.0165 0.0355 0.0190 2
Pooled: median post minus pre 0.0119 ; positive 28 of 36 ; Wilcoxon p = 3.6e-05
patterns median_change positive wilcoxon_p
rat
Rat08 21 0.017012 18 0.000131
Rat10 8 0.002416 5 0.742188
Rat11 7 0.019026 5 0.109375
Reference: non-significant patterns, median post minus pre 0.0001 over 253
Session medians above zero: 14 of 15 ; Wilcoxon over session medians p = 0.00012
In [15]:
# Variants: 25 and 100 ms bins, the wider BLA labels, and the first pattern of every session.
rng_variants = np.random.default_rng(61002)
rows = []
for name, kwargs, sess in (("25 ms bins", {"bin_s": 0.025}, used), ("100 ms bins", {"bin_s": 0.1}, used),
("wider BLA labels", {"bla_labels": hpcbla.BLA_COMPLEX_LABELS},
[s for s in sessions if s["units"] and hpcbla.usable(s, hpcbla.BLA_COMPLEX_LABELS)])):
df = pd.DataFrame([r for s in sess for r in hpcbla.reactivation(s, rng_variants, keep_traces=False, **kwargs)["rows"]])
df["change"] = df.post - df.pre
for subset, d in (("significant", df[df.significant]), ("first pattern", df[df.pattern == 0])):
rows.append({"variant": name, "patterns": subset, "n": len(d), "sessions": d.session.nunique(),
"median_change": round(d.change.median(), 4), "positive": int((d.change > 0).sum()),
"p": wilcoxon(d.change),
"rats_median_above_0": int((d.groupby("rat").change.median() > 0).sum()), "rats": d.rat.nunique()})
first = patterns[patterns.pattern == 0]
rows.append({"variant": "main (50 ms)", "patterns": "first pattern", "n": len(first), "sessions": first.session.nunique(),
"median_change": round(first.change.median(), 4), "positive": int((first.change > 0).sum()),
"p": wilcoxon(first.change), "rats_median_above_0": int((first.groupby("rat").change.median() > 0).sum()),
"rats": first.rat.nunique()})
print(pd.DataFrame(rows).to_string(index=False))
variant patterns n sessions median_change positive p rats_median_above_0 rats
25 ms bins significant 38 15 0.0062 29 0.000009 2 3
25 ms bins first pattern 15 15 0.0064 12 0.003357 3 3
100 ms bins significant 36 15 0.0213 29 0.000011 3 3
100 ms bins first pattern 15 15 0.0223 13 0.001526 3 3
wider BLA labels significant 57 21 0.0101 46 0.000001 3 3
wider BLA labels first pattern 21 21 0.0156 18 0.002482 3 3
main (50 ms) first pattern 15 15 0.0123 12 0.008362 2 3
In [16]:
fig, ax = plt.subplots(figsize=(5, 4))
for rat, d in sig.groupby("rat"):
ax.scatter(d.pre, d.post, s=14, label=rat)
lim = [min(sig.pre.min(), sig.post.min()), max(sig.pre.max(), sig.post.max())]
ax.plot(lim, lim, "k--", lw=0.8)
ax.set_xlabel("pre-sleep Non-REM strength")
ax.set_ylabel("post-sleep Non-REM strength")
ax.legend(fontsize=8)
plt.show()
print("Plotted ranges: pre", round(sig.pre.min(), 4), "to", round(sig.pre.max(), 4), "; post",
round(sig.post.min(), 4), "to", round(sig.post.max(), 4), "; no point outside the axes")
Plotted ranges: pre -0.0145 to 0.1687 ; post -0.0246 to 0.1617 ; no point outside the axes
Result: reproduced. All 15 sessions have at least one significant joint pattern (36 in all). Over these patterns the post-sleep minus pre-sleep Non-REM strength has a median of 0.0119, is positive for 28 of 36, and the Wilcoxon p is 3.6e-05. The median is above zero in all three rats (Rat08 0.0170, 18 of 21; Rat10 0.0024, 5 of 8; Rat11 0.0190, 5 of 7), though in Rat10 and Rat11 alone the test is not significant (p = 0.74 and 0.11). Non-significant patterns change by a median of 0.0001. Every variant (25 and 100 ms bins, the wider BLA labels, the first pattern of each session) gives a positive median with p of 0.0084 or less.
reactivation-peaks-at-ripples¶
The hippocampus and BLA reactivations peak during hippocampal sharp wave-ripples.
The ripple-peak-triggered mean of each significant pattern's strength in post-sleep Non-REM, from -1 to +1 s in 50 ms bins, and reading (b): the post minus pre change in bins that overlap a ripple against bins outside ripples.
In [17]:
traces = pd.DataFrame([t for out in main.values() for t in out["traces"]])
lags_s = (np.arange(-20, 21)) * hpcbla.BIN_S
curves = np.vstack(traces.curve)
curves_pre = np.vstack(traces.curve_pre)
mean_curve = curves.mean(0)
print("Patterns:", len(traces), "; post-sleep Non-REM ripples per session: median", int(traces.groupby("session").n_ripples.first().median()))
print("Reading (a): the mean curve is highest at lag", round(lags_s[np.argmax(mean_curve)], 3), "s")
print(pd.DataFrame({"lag_s": lags_s, "post": mean_curve, "pre": curves_pre.mean(0)}).iloc[14:27].round(4).to_string(index=False))
flank = (np.abs(lags_s) >= 0.5 - 1e-9)
peak_minus_flank = curves[:, 20] - curves[:, flank].mean(1)
print("Peak bin minus flanks (0.5 to 1 s): median", round(float(np.median(peak_minus_flank)), 4), "; positive",
int((peak_minus_flank > 0).sum()), "of", len(peak_minus_flank), "; Wilcoxon p =", f"{wilcoxon(peak_minus_flank):.2g}")
print("Ratio of the mean curve at the peak bin to its flanks:", round(float(mean_curve[20] / mean_curve[flank].mean()), 2))
traces["peak_minus_flank"] = peak_minus_flank
traces["argmax_lag_s"] = lags_s[np.argmax(curves, 1)]
print("Patterns whose own curve peaks within 100 ms of the ripple:", int((np.abs(traces.argmax_lag_s) <= 0.1 + 1e-9).sum()), "of", len(traces))
print(traces.groupby("rat").agg(patterns=("pattern", "size"), median_peak_minus_flank=("peak_minus_flank", "median"),
positive=("peak_minus_flank", lambda x: int((x > 0).sum()))).round(4).to_string())
pre_pmf = curves_pre[:, 20] - curves_pre[:, flank].mean(1)
print("Reference, pre-sleep: peak bin minus flanks median", round(float(np.median(pre_pmf)), 4),
"; post minus pre of this difference: median", round(float(np.median(peak_minus_flank - pre_pmf)), 4),
", positive", int((peak_minus_flank - pre_pmf > 0).sum()), "of", len(traces), ", Wilcoxon p =",
f"{wilcoxon(peak_minus_flank - pre_pmf):.2g}")
Patterns: 36 ; post-sleep Non-REM ripples per session: median 2661
Reading (a): the mean curve is highest at lag 0.0 s
lag_s post pre
-0.30 0.0559 0.0278
-0.25 0.0573 0.0253
-0.20 0.0683 0.0283
-0.15 0.0736 0.0351
-0.10 0.0779 0.0326
-0.05 0.0917 0.0562
0.00 0.2234 0.1362
0.05 0.0892 0.0481
0.10 0.0879 0.0428
0.15 0.0940 0.0501
0.20 0.0748 0.0339
0.25 0.0686 0.0399
0.30 0.0794 0.0416
Peak bin minus flanks (0.5 to 1 s): median 0.1032 ; positive 30 of 36 ; Wilcoxon p = 6e-06
Ratio of the mean curve at the peak bin to its flanks: 4.36
Patterns whose own curve peaks within 100 ms of the ripple: 22 of 36
patterns median_peak_minus_flank positive
rat
Rat08 21 0.1067 16
Rat10 8 0.0647 8
Rat11 7 0.1182 6
Reference, pre-sleep: peak bin minus flanks median 0.0419 ; post minus pre of this difference: median 0.0611 , positive 23 of 36 , Wilcoxon p = 0.022
In [18]:
traces["change_rip"] = traces.post_rip - traces.pre_rip
traces["change_out"] = traces.post_out - traces.pre_out
d = traces.change_rip - traces.change_out
print("Reading (b): post minus pre in ripple bins, median", round(traces.change_rip.median(), 4),
"; outside ripples", round(traces.change_out.median(), 4))
print("Difference (ripple minus outside): median", round(float(d.median()), 4), "; positive", int((d > 0).sum()), "of", len(d),
"; Wilcoxon p =", f"{wilcoxon(d):.2g}")
print(traces.assign(d=d).groupby("rat").d.agg(["size", "median", lambda x: int((x > 0).sum())]).round(4).to_string())
fig, ax = plt.subplots(figsize=(6, 3.5))
ax.plot(lags_s, mean_curve, label="post-sleep")
ax.plot(lags_s, curves_pre.mean(0), label="pre-sleep")
ax.axvline(0, color="k", lw=0.5)
ax.set_xlabel("time from ripple peak (s)")
ax.set_ylabel("mean reactivation strength")
ax.legend()
plt.show()
print("Plotted range:", round(float(min(mean_curve.min(), curves_pre.mean(0).min())), 4), "to",
round(float(max(mean_curve.max(), curves_pre.mean(0).max())), 4))
Reading (b): post minus pre in ripple bins, median 0.0654 ; outside ripples 0.0115
Difference (ripple minus outside): median 0.0537 ; positive 25 of 36 ; Wilcoxon p = 0.0033
size median <lambda_0>
rat
Rat08 21 0.0731 15
Rat10 8 0.0257 5
Rat11 7 0.0509 5
Plotted range: 0.0163 to 0.2234
In [19]:
# Session level: the per-pattern measures above averaged within each session, tested over sessions.
traces["d_rip"] = traces.change_rip - traces.change_out
traces["pmf_post_minus_pre"] = traces.peak_minus_flank - pre_pmf
by_session = traces.groupby(["rat", "session"])[["peak_minus_flank", "pmf_post_minus_pre", "d_rip"]].mean()
print(by_session.round(4).to_string())
for col, label in (("peak_minus_flank", "reading (a), peak minus flanks"),
("pmf_post_minus_pre", "reading (a), post minus pre of peak minus flanks"),
("d_rip", "reading (b), ripple minus outside change")):
x = by_session[col]
print(f"{label}: median over {len(x)} sessions {x.median():.4f}, positive {int((x > 0).sum())}, Wilcoxon p = {wilcoxon(x):.2g}")
peak_minus_flank pmf_post_minus_pre d_rip
rat session
Rat08 Rat08-20130708 -0.0364 -0.0543 -0.0303
Rat08-20130709 0.3269 0.0677 0.0713
Rat08-20130710 0.1195 -0.0354 -0.0052
Rat08-20130711 0.1400 0.0662 0.0574
Rat08-20130712 0.1155 0.0819 0.1167
Rat08-20130713 -0.0367 -0.0132 -0.0160
Rat08-20130715 0.1078 0.0981 0.1109
Rat08-20130716 0.1447 0.0658 0.0677
Rat08-20130717 0.5242 0.4777 0.3814
Rat10 Rat10-20140624 0.0403 -0.0275 -0.0075
Rat10-20140626 0.0756 0.0440 0.0349
Rat10-20140704 0.2362 0.4545 0.1815
Rat11 Rat11-20150325 0.2357 0.1476 0.1495
Rat11-20150326 0.3735 0.1452 0.0785
Rat11-20150328 0.4831 -0.1272 -0.1206
reading (a), peak minus flanks: median over 15 sessions 0.1400, positive 13, Wilcoxon p = 0.00031
reading (a), post minus pre of peak minus flanks: median over 15 sessions 0.0662, positive 10, Wilcoxon p = 0.03
reading (b), ripple minus outside change: median over 15 sessions 0.0677, positive 10, Wilcoxon p = 0.03
Result: reproduced. Under reading (a), the mean ripple-triggered strength in post-sleep is highest in the bin of the ripple peak, 4.36 times its mean 0.5 to 1 s away; the peak bin exceeds the flanks in 30 of 36 patterns (median 0.103, p = 6e-06), with a positive median in every rat, and in 13 of 15 sessions at session level (p = 0.00031). Much of this peak is already present in pre-sleep (0.136 at the ripple, against 0.223 in post-sleep), so reading (a) alone does not separate reactivation of the run from co-activation at ripples in general. The part that does, the peak over flanks being larger in post-sleep than in pre-sleep, holds in 23 of 36 patterns (p = 0.022) and 10 of 15 sessions (p = 0.03). Under reading (b), the post minus pre change is larger in ripple bins than outside them (median difference 0.0537, 25 of 36 patterns, p = 0.0033; 10 of 15 sessions, p = 0.03).
bla-ripple-up-subgroup¶
A subgroup of BLA cells is positively modulated during hippocampal sharp wave-ripples.
Spikes of each BLA unit in the windows of the Non-REM ripples of pre- and post-sleep, against 500 shifts of all windows together by 0.5 to 5 s. The chance level is the fraction selected with the real windows replaced by the windows shifted by 7 s.
In [20]:
rng_rip = np.random.default_rng(61003)
started = time.time()
rows = []
for s in sessions:
u = s["units"]
if s["states"] is None or s["ripples"] is None or u is None:
continue
pre, post, _ = hpcbla.sleep_epochs(s)
nrem = hpcbla.intersect(hpcbla.intervals_of(s, "Non-REM"), np.vstack([pre, post]))
r = s["ripples"]
keep = hpcbla.in_intervals(r["peak"], nrem)
start, stop = r["start"][keep], r["stop"][keep]
end = max(b for _, _, b in s["epochs"])
used_bla = set(hpcbla.unit_sets(s)[1])
for i in u["trains"]:
if u["region"][i] not in hpcbla.BLA_LABELS:
continue
real, chance = hpcbla.ripple_modulation(u["trains"][i], [(start, stop), (start + hpcbla.CHANCE_SHIFT_S, stop + hpcbla.CHANCE_SHIFT_S)], end, rng_rip)
rows.append({"session": s["session"], "rat": s["rat"], "unit": i, "ripples": int(keep.sum()),
"in_reactivation_set": i in used_bla and hpcbla.usable(s),
"observed": real["observed"], "expected": real["null_mean"], "p_up": real["p_up"], "p_down": real["p_down"],
"chance_p_up": chance["p_up"], "chance_p_down": chance["p_down"]})
mod = pd.DataFrame(rows)
mod["up"] = mod.p_up <= hpcbla.RIPPLE_ALPHA
mod["down"] = mod.p_down <= hpcbla.RIPPLE_ALPHA
mod["chance_up"] = mod.chance_p_up <= hpcbla.RIPPLE_ALPHA
print(f"computed in {time.time() - started:.0f} s; BLA units", len(mod), "in", mod.session.nunique(), "sessions;",
"units with no spike in any window", int((mod.observed == 0).sum()))
chance = mod.chance_up.mean()
n_up = int(mod.up.sum())
print("Positively modulated:", n_up, "of", len(mod), f"({mod.up.mean():.3f}); negatively: {int(mod.down.sum())}",
f"({mod.down.mean():.3f}); chance level from the 7 s shifted windows: {chance:.3f} ({int(mod.chance_up.sum())} units)")
print("Binomial test of the positive fraction against the chance level: p =",
f"{stats.binomtest(n_up, len(mod), max(chance, 1 / len(mod)), alternative='greater').pvalue:.2g}")
by = mod.groupby("rat").agg(units=("up", "size"), up=("up", "sum"), down=("down", "sum"), chance_up=("chance_up", "sum"))
by["up_fraction"] = by.up / by.units
by["chance_fraction"] = by.chance_up / by.units
by["p"] = [stats.binomtest(int(r.up), int(r.units), max(chance, 1 / len(mod)), alternative="greater").pvalue for r in by.itertuples()]
# Against each rat's own chance level (at least one unit, so that a chance of zero still gives a test).
by["p_own_chance"] = [stats.binomtest(int(r.up), int(r.units), max(r.chance_fraction, 1 / r.units), alternative="greater").pvalue
for r in by.itertuples()]
print(by.round(4).to_string())
print(mod.groupby("session").agg(units=("up", "size"), up=("up", "sum"), down=("down", "sum")).T.to_string())
print("\nThreshold variants (fraction up, chance fraction):")
for alpha in (0.05, 0.01, 0.002):
print(f" p <= {alpha}: {np.mean(mod.p_up <= alpha):.3f}, chance {np.mean(mod.chance_p_up <= alpha):.3f}")
computed in 92 s; BLA units 822 in 17 sessions; units with no spike in any window 0
Positively modulated: 103 of 822 (0.125); negatively: 104 (0.127); chance level from the 7 s shifted windows: 0.007 (6 units)
Binomial test of the positive fraction against the chance level: p = 9.8e-90
units up down chance_up up_fraction chance_fraction p p_own_chance
rat
Rat08 496 51 78 3 0.1028 0.0060 0.0000 0.0000
Rat10 84 4 14 0 0.0476 0.0000 0.0034 0.0183
Rat11 242 48 12 3 0.1983 0.0124 0.0000 0.0000
session Rat08-20130708 Rat08-20130709 Rat08-20130710 Rat08-20130711 Rat08-20130712 Rat08-20130713 Rat08-20130715 Rat08-20130716 Rat08-20130717 Rat10-20140624 Rat10-20140626 Rat10-20140701 Rat10-20140702 Rat10-20140704 Rat11-20150325 Rat11-20150326 Rat11-20150328
units 15 78 87 92 40 94 39 34 17 37 33 4 1 9 86 72 84
up 1 11 5 8 3 6 5 8 4 4 0 0 0 0 14 16 18
down 2 22 16 19 7 8 1 3 0 4 7 2 0 1 3 6 3
Threshold variants (fraction up, chance fraction):
p <= 0.05: 0.163, chance 0.056
p <= 0.01: 0.125, chance 0.007
p <= 0.002: 0.108, chance 0.002
In [21]:
# Where the shifted windows land: the fraction of shifted ripple peaks that fall outside the
# Non-REM of pre- and post-sleep, for the 500 random offsets and for the 7 s shift.
rng_null = np.random.default_rng(61004)
lo, hi = (int(x / hpcbla.RIPPLE_RES_S) for x in hpcbla.RIPPLE_SHIFT_RANGE_S)
rows = []
for s in sessions:
if s["states"] is None or s["ripples"] is None or s["units"] is None or s["session"] not in set(mod.session):
continue
pre, post, _ = hpcbla.sleep_epochs(s)
nrem = hpcbla.intersect(hpcbla.intervals_of(s, "Non-REM"), np.vstack([pre, post]))
peaks = s["ripples"]["peak"][hpcbla.in_intervals(s["ripples"]["peak"], nrem)]
offsets = rng_null.integers(lo, hi + 1, hpcbla.RIPPLE_SHIFTS) * rng_null.choice([-1, 1], hpcbla.RIPPLE_SHIFTS) * hpcbla.RIPPLE_RES_S
out_random = np.mean([np.mean(~hpcbla.in_intervals(peaks + o, nrem)) for o in offsets])
rows.append({"session": s["session"], "ripples": len(peaks), "outside_nrem_random_shifts": out_random,
"outside_nrem_7s": np.mean(~hpcbla.in_intervals(peaks + hpcbla.CHANCE_SHIFT_S, nrem))})
null_place = pd.DataFrame(rows)
print(null_place.round(3).to_string(index=False))
print("Median over sessions: outside Non-REM with random shifts", round(float(null_place.outside_nrem_random_shifts.median()), 3),
"; with the 7 s shift", round(float(null_place.outside_nrem_7s.median()), 3))
session ripples outside_nrem_random_shifts outside_nrem_7s Rat08-20130708 7574 0.016 0.037 Rat08-20130709 8871 0.018 0.050 Rat08-20130710 7243 0.011 0.028 Rat08-20130711 6139 0.032 0.098 Rat08-20130712 4994 0.021 0.059 Rat08-20130713 5237 0.009 0.015 Rat08-20130715 4124 0.009 0.014 Rat08-20130716 6892 0.010 0.021 Rat08-20130717 6457 0.005 0.006 Rat10-20140624 5500 0.006 0.010 Rat10-20140626 3181 0.017 0.029 Rat10-20140701 6692 0.007 0.009 Rat10-20140702 2695 0.007 0.009 Rat10-20140704 4327 0.007 0.011 Rat11-20150325 2849 0.010 0.011 Rat11-20150326 2832 0.008 0.008 Rat11-20150328 2395 0.005 0.008 Median over sessions: outside Non-REM with random shifts 0.009 ; with the 7 s shift 0.014
Result: reproduced. 103 of 822 BLA units (12.5%) fire more in Non-REM ripples than in the shifted windows, against a chance level of 0.7% from the procedure run on windows shifted by 7 s (binomial p = 9.8e-90), and the fraction is above chance in each rat (Rat08 10.3%, Rat10 4.8%, Rat11 19.8%), against the pooled chance level and against each rat's own (Rat10 p = 0.0034 and 0.018). It is a subgroup, since the fraction is well below 1. A similar number of units (104, 12.7%) fire less in ripples. The shifted windows stay mostly in Non-REM (a median of 0.9% of the randomly shifted ripple peaks and 1.4% of the 7 s shifted ones fall outside it). Units are counted once per session, so a cell recorded on several days may be counted more than once.
reactivation-involves-ripple-up-cells¶
The hippocampus and BLA reactivations involve the subgroup of BLA cells that are positively modulated during sharp wave-ripples.
For each significant pattern of the run, the mean squared BLA weight of the positively modulated units minus that of the other BLA units used.
In [22]:
up_units = {(r.session, r.unit) for r in mod[mod.up].itertuples()}
rows = []
for r in sig.itertuples():
is_up = np.array([(r.session, i) in up_units for i in r.bla_units])
if is_up.any() and (~is_up).any():
w2 = np.asarray(r.bla_weights) ** 2
rows.append({"session": r.session, "rat": r.rat, "pattern": r.pattern, "change": r.change,
"n_up": int(is_up.sum()), "n_other": int((~is_up).sum()),
"w2_up": w2[is_up].mean(), "w2_other": w2[~is_up].mean()})
inv = pd.DataFrame(rows)
inv["diff"] = inv.w2_up - inv.w2_other
print("Patterns with both kinds of BLA unit:", len(inv), "of", len(sig), "in", inv.session.nunique(), "sessions")
print("Mean squared weight, up units: median", round(inv.w2_up.median(), 4), "; other units:", round(inv.w2_other.median(), 4),
"; equal weights would give", round(float((1 / (inv.n_up + inv.n_other)).median()), 4))
print("Difference: median", round(inv["diff"].median(), 4), "; positive", int((inv["diff"] > 0).sum()), "of", len(inv),
"; Wilcoxon p =", f"{wilcoxon(inv['diff']):.2g}")
print(inv.groupby("rat")["diff"].agg(["size", "median", lambda x: int((x > 0).sum())]).round(4).to_string())
reac = inv[inv.change > 0]
print("Variant, patterns with post minus pre above zero:", len(reac), "; median difference", round(reac["diff"].median(), 4),
"; positive", int((reac["diff"] > 0).sum()), "; Wilcoxon p =", f"{wilcoxon(reac['diff']):.2g}")
print("Fraction of the BLA weight (sum of squares) on up units: median", round(float((inv.w2_up * inv.n_up / (inv.w2_up * inv.n_up + inv.w2_other * inv.n_other)).median()), 3),
"; fraction of BLA units that are up: median", round(float((inv.n_up / (inv.n_up + inv.n_other)).median()), 3))
Patterns with both kinds of BLA unit: 31 of 36 in 13 sessions
Mean squared weight, up units: median 0.0505 ; other units: 0.0181 ; equal weights would give 0.0233
Difference: median 0.0339 ; positive 24 of 31 ; Wilcoxon p = 0.0013
size median <lambda_0>
rat
Rat08 21 0.0519 17
Rat10 3 -0.0413 1
Rat11 7 0.0147 6
Variant, patterns with post minus pre above zero: 25 ; median difference 0.0339 ; positive 19 ; Wilcoxon p = 0.0018
Fraction of the BLA weight (sum of squares) on up units: median 0.342 ; fraction of BLA units that are up: median 0.2
In [23]:
# Firing rate. Run-epoch rate of each BLA unit used; within each pattern, each up unit is paired
# with the other BLA unit of the nearest log rate (with replacement), and the mean squared weight
# of the up units is compared with that of their matched units.
run_rate = {}
for s in used:
_, _, run = hpcbla.sleep_epochs(s)
for i in hpcbla.unit_sets(s)[1]:
t = s["units"]["trains"][i]
run_rate[(s["session"], i)] = np.count_nonzero((t >= run[0]) & (t < run[1])) / (run[1] - run[0])
unit_rates = pd.DataFrame([{"session": k[0], "unit": k[1], "rate": v, "up": k in up_units} for k, v in run_rate.items()])
print("Run-epoch rate of BLA units used (Hz): up units median", round(unit_rates[unit_rates.up].rate.median(), 2),
f"(n = {int(unit_rates.up.sum())}); other units median", round(unit_rates[~unit_rates.up].rate.median(), 2),
f"(n = {int((~unit_rates.up).sum())}); Mann-Whitney p = {stats.mannwhitneyu(unit_rates[unit_rates.up].rate, unit_rates[~unit_rates.up].rate).pvalue:.2g}")
rows = []
for r in sig.itertuples():
is_up = np.array([(r.session, i) in up_units for i in r.bla_units])
if not (is_up.any() and (~is_up).any()):
continue
w2 = np.asarray(r.bla_weights) ** 2
log_rate = np.log(np.array([run_rate[(r.session, i)] for i in r.bla_units]))
others = np.flatnonzero(~is_up)
match = np.array([others[np.argmin(np.abs(log_rate[others] - log_rate[j]))] for j in np.flatnonzero(is_up)])
rows.append({"session": r.session, "rat": r.rat, "pattern": r.pattern,
"rho_rate_w2": stats.spearmanr(log_rate, w2)[0],
"matched_diff": w2[is_up].mean() - w2[match].mean(),
"matched_log_rate_gap": float(np.median(np.abs(log_rate[is_up] - log_rate[match])))})
rate_ctl = pd.DataFrame(rows)
print("Spearman correlation of run rate and squared weight within a pattern: median", round(rate_ctl.rho_rate_w2.median(), 3),
"; positive in", int((rate_ctl.rho_rate_w2 > 0).sum()), "of", len(rate_ctl))
print("Rate-matched difference (up minus matched other): median", round(rate_ctl.matched_diff.median(), 4), "; positive",
int((rate_ctl.matched_diff > 0).sum()), "of", len(rate_ctl), f"; Wilcoxon p = {wilcoxon(rate_ctl.matched_diff):.2g}",
"; median ratio of matched rates", round(float(np.exp(rate_ctl.matched_log_rate_gap.median())), 2))
inv_session = inv.groupby(["rat", "session"])["diff"].mean()
print("\nSession level: mean difference per session, median", round(inv_session.median(), 4), "; positive",
int((inv_session > 0).sum()), "of", len(inv_session), f"; Wilcoxon p = {wilcoxon(inv_session):.2g}")
matched_session = rate_ctl.groupby("session").matched_diff.mean()
print("Session level, rate-matched: median", round(matched_session.median(), 4), "; positive", int((matched_session > 0).sum()),
"of", len(matched_session), f"; Wilcoxon p = {wilcoxon(matched_session):.2g}")
Run-epoch rate of BLA units used (Hz): up units median 1.66 (n = 84); other units median 0.41 (n = 399); Mann-Whitney p = 7.2e-11 Spearman correlation of run rate and squared weight within a pattern: median 0.463 ; positive in 30 of 31 Rate-matched difference (up minus matched other): median 0.0355 ; positive 24 of 31 ; Wilcoxon p = 0.012 ; median ratio of matched rates 1.09 Session level: mean difference per session, median 0.0269 ; positive 10 of 13 ; Wilcoxon p = 0.027 Session level, rate-matched: median 0.0219 ; positive 9 of 13 ; Wilcoxon p = 0.22
Result: reproduced. In 31 significant patterns that have both kinds of BLA unit, the mean squared weight of the ripple-up units is above that of the other BLA units (medians 0.0505 and 0.0181; difference positive in 24 of 31, p = 0.0013). The ripple-up units are a median 20% of the BLA units of a pattern and carry 34% of its BLA weight. It holds in Rat08 (17 of 21) and Rat11 (6 of 7) and not in Rat10 (1 of 3), and in the patterns that reactivate (p = 0.0018). Firing rate accounts for part of it: up units fire faster in the run (median 1.66 against 0.41 Hz), and squared weight rises with rate within a pattern (median Spearman 0.46). Against other units of matched rate the up units still carry more weight (24 of 31 patterns, p = 0.012). At session level the difference holds (10 of 13 sessions, p = 0.027), but the rate-matched difference does not (9 of 13, p = 0.22).
aversive-direction-stronger¶
Reactivation is stronger for the hippocampus and BLA correlation patterns of the run direction that held the air puff than for those of the safe direction.
Templates are built from the run bins inside laps of each direction (laps from the reward times; the air-puff direction holds most puffs).
In [24]:
rng_dir = np.random.default_rng(61005)
started = time.time()
directional = {s["session"]: hpcbla.reactivation(s, rng_dir, templates=("puff", "safe"), keep_traces=False) for s in used}
equal = {s["session"]: hpcbla.reactivation(s, rng_dir, templates=("puff", "safe"), equal_bins=True, keep_traces=False) for s in used}
print(f"computed in {time.time() - started:.0f} s")
def direction_table(results):
df = pd.DataFrame([r for out in results.values() for r in out["rows"]])
df["change"] = df.post - df.pre
return df
dpat = direction_table(directional)
print("Template bins per session (puff, safe):", {k: (v["template_bins"].get("puff"), v["template_bins"].get("safe")) for k, v in directional.items()})
def compare(df, value="change", subset="significant"):
d = df[df.significant] if subset == "significant" else df[df.pattern == 0]
per = d.groupby(["rat", "session", "template"])[value].mean().unstack("template")
per = per.dropna(subset=["puff", "safe"]) if {"puff", "safe"} <= set(per.columns) else per.iloc[0:0]
per["diff"] = per.puff - per.safe
return d, per
d, per = compare(dpat)
print("\nSignificant patterns: puff", int((d.template == "puff").sum()), "; safe", int((d.template == "safe").sum()))
print(per.round(4).to_string())
print("\nSessions with both:", len(per), "; median puff minus safe", round(per["diff"].median(), 4), "; positive",
int((per["diff"] > 0).sum()), "of", len(per), "; Wilcoxon p =", f"{wilcoxon(per['diff']):.2g}")
print(per.groupby("rat")["diff"].agg(["size", "median", lambda x: int((x > 0).sum())]).round(4).to_string())
computed in 48 s
Template bins per session (puff, safe): {'Rat08-20130708': (39915, 31622), 'Rat08-20130709': (66804, 5610), 'Rat08-20130710': (28766, 43378), 'Rat08-20130711': (55135, 11400), 'Rat08-20130712': (24040, 38048), 'Rat08-20130713': (46133, 18737), 'Rat08-20130715': (28145, 16782), 'Rat08-20130716': (25790, 22282), 'Rat08-20130717': (19838, 16102), 'Rat10-20140624': (16528, 21264), 'Rat10-20140626': (12509, 14409), 'Rat10-20140704': (16634, 11126), 'Rat11-20150325': (7960, 6947), 'Rat11-20150326': (6567, 11501), 'Rat11-20150328': (12479, 7979)}
Significant patterns: puff 35 ; safe 31
template puff safe diff
rat session
Rat08 Rat08-20130708 0.0047 0.0171 -0.0124
Rat08-20130709 0.0282 0.0061 0.0221
Rat08-20130710 0.0177 -0.0058 0.0235
Rat08-20130711 0.0249 -0.0044 0.0293
Rat08-20130712 0.0045 0.0151 -0.0106
Rat08-20130713 0.0129 -0.0013 0.0142
Rat08-20130715 0.0112 -0.0022 0.0133
Rat08-20130716 0.0201 0.0185 0.0016
Rat08-20130717 0.0300 0.0372 -0.0073
Rat10 Rat10-20140624 -0.0040 0.0116 -0.0156
Rat10-20140626 0.0150 -0.0051 0.0200
Rat10-20140704 0.0119 -0.0032 0.0151
Rat11 Rat11-20150325 0.0161 0.0209 -0.0048
Rat11-20150326 0.0184 0.0040 0.0144
Rat11-20150328 0.0045 0.0120 -0.0075
Sessions with both: 15 ; median puff minus safe 0.0133 ; positive 9 of 15 ; Wilcoxon p = 0.11
size median <lambda_0>
rat
Rat08 9 0.0133 6
Rat10 3 0.0151 2
Rat11 3 -0.0048 1
In [25]:
rows = []
pooled = d.groupby("template").change
mw = stats.mannwhitneyu(d[d.template == "puff"].change, d[d.template == "safe"].change)
rows.append({"variant": "pooled per pattern (Mann-Whitney)", "n": f"{int((d.template == 'puff').sum())}/{int((d.template == 'safe').sum())}",
"median_puff": round(pooled.median()["puff"], 4), "median_safe": round(pooled.median()["safe"], 4), "p": mw.pvalue})
for name, df, subset, value in (("first pattern of each direction", dpat, "first", "change"),
("equal numbers of bins", direction_table(equal), "significant", "change"),
("post-sleep strength alone", dpat, "significant", "post")):
_, p = compare(df, value, subset)
rows.append({"variant": name, "n": len(p), "median_puff": round(p.puff.median(), 4), "median_safe": round(p.safe.median(), 4),
"median_diff": round(p["diff"].median(), 4), "positive": int((p["diff"] > 0).sum()), "p": wilcoxon(p["diff"])})
print(pd.DataFrame(rows).to_string(index=False))
variant n median_puff median_safe p median_diff positive
pooled per pattern (Mann-Whitney) 35/31 0.0131 0.0054 0.082831 NaN NaN
first pattern of each direction 15 0.0240 0.0004 0.047913 0.0227 10.0
equal numbers of bins 14 0.0155 0.0081 0.325806 0.0024 8.0
post-sleep strength alone 15 0.0338 0.0236 0.454285 0.0014 9.0
In [26]:
# Laps of the run epoch. A lap runs from one reward to the next at the other end, so it holds the
# time spent at the starting end as well as the run. The shortest laps of a session (its 10th
# percentile of lap duration) are taken as the time needed to travel the track; the variant below
# builds the templates from only that last part of each lap, before the reward that ends it.
rows, lap_tail = [], {}
for s in used:
direction, counts = hpcbla.puff_direction(s)
laps_ = hpcbla.run_laps(s)
lap_tail[s["session"]] = float(np.percentile([b - a for a, b, _ in laps_], 10))
row = {"session": s["session"], "puffs": int(sum(counts.values())), "puffs_in_puff_laps": counts[direction], "travel_s": round(lap_tail[s["session"]], 1)}
for name in ("puff", "safe"):
dur = np.array([b - a for a, b, dd in laps_ if (dd == direction) == (name == "puff")])
row[f"{name}_laps"] = len(dur)
row[f"{name}_median_s"] = round(float(np.median(dur)), 1)
row[f"{name}_max_s"] = round(float(dur.max()), 1)
row[f"{name}_bins"] = directional[s["session"]]["template_bins"][name]
row[f"{name}_travel_fraction"] = round(float(np.minimum(dur, lap_tail[s["session"]]).sum() / dur.sum()), 3)
rows.append(row)
lap_table = pd.DataFrame(rows)
print(lap_table.to_string(index=False))
print("Puffs inside laps of the air-puff direction, in these sessions:", int(lap_table.puffs_in_puff_laps.sum()), "of", int(lap_table.puffs.sum()))
print("Fraction of template time inside the travel part of laps: median puff", lap_table.puff_travel_fraction.median(),
"; safe", lap_table.safe_travel_fraction.median())
rng_tail = np.random.default_rng(61006)
tail = {s["session"]: hpcbla.reactivation(s, rng_tail, templates=("puff", "safe"), lap_tail_s=lap_tail[s["session"]],
keep_traces=False) for s in used}
rows = []
few = set(lap_table.session[lap_table.puffs < 5])
for name, df, keep in (("main", dpat, None), ("sessions with at least 5 puffs", dpat, few),
("travel part of each lap only", direction_table(tail), None)):
_, p = compare(df)
if keep:
p = p[~p.index.get_level_values("session").isin(keep)]
rows.append({"variant": name, "sessions": len(p), "median_diff": round(p["diff"].median(), 4),
"positive": int((p["diff"] > 0).sum()), "p": wilcoxon(p["diff"])})
print("\nSessions left out for fewer than 5 puffs:", sorted(few))
print(pd.DataFrame(rows).to_string(index=False))
session puffs puffs_in_puff_laps travel_s puff_laps puff_median_s puff_max_s puff_bins puff_travel_fraction safe_laps safe_median_s safe_max_s safe_bins safe_travel_fraction Rat08-20130708 3 3 20.5 2 997.9 1976.4 39915 0.020 2 790.5 1558.0 31622 0.026 Rat08-20130709 2 2 27.5 4 500.2 2309.9 66804 0.033 3 87.7 169.0 5610 0.281 Rat08-20130710 9 8 74.5 6 191.9 504.5 28766 0.284 7 238.4 596.8 43378 0.240 Rat08-20130711 5 4 84.0 4 660.5 1322.3 55135 0.122 5 103.5 176.9 11400 0.608 Rat08-20130712 14 14 20.5 11 66.7 284.8 24040 0.187 12 47.2 1097.1 38048 0.119 Rat08-20130713 20 19 55.0 10 121.9 609.5 46133 0.238 10 79.0 189.5 18737 0.540 Rat08-20130715 34 33 23.2 20 56.5 182.1 28145 0.330 21 29.1 130.4 16782 0.572 Rat08-20130716 30 28 17.0 24 46.4 179.8 25790 0.314 24 24.2 366.7 22282 0.361 Rat08-20130717 31 31 15.0 27 28.2 111.9 19838 0.405 27 18.3 172.6 16102 0.493 Rat10-20140624 26 24 21.4 19 37.9 97.5 16528 0.489 20 51.0 111.9 21264 0.380 Rat10-20140626 20 20 18.1 17 32.2 75.0 12509 0.489 17 40.3 89.9 14409 0.423 Rat10-20140704 9 8 12.5 26 19.3 242.5 16634 0.391 26 16.7 54.0 11126 0.567 Rat11-20150325 32 32 7.7 30 10.3 39.8 7960 0.563 30 9.3 47.0 6947 0.666 Rat11-20150326 32 32 7.9 31 10.3 18.3 6567 0.742 32 13.9 77.3 11501 0.430 Rat11-20150328 37 37 7.6 32 12.7 73.8 12479 0.389 33 9.6 38.0 7979 0.627 Puffs inside laps of the air-puff direction, in these sessions: 295 of 304 Fraction of template time inside the travel part of laps: median puff 0.33 ; safe 0.43
Sessions left out for fewer than 5 puffs: ['Rat08-20130708', 'Rat08-20130709']
variant sessions median_diff positive p
main 15 0.0133 9 0.106995
sessions with at least 5 puffs 13 0.0133 8 0.146484
travel part of each lap only 14 0.0052 9 0.193726
Result: not reproduced. Across the 15 sessions the post minus pre reactivation of the air-puff templates exceeds that of the safe templates by a median of 0.0133, in 9 of 15 sessions, with Wilcoxon p = 0.11, above the 0.05 of the criterion and the 0.1 of its partial outcome. By rat the difference is positive in 6 of 9 Rat08 sessions, 2 of 3 Rat10 and 1 of 3 Rat11. The variants lean the same way with mixed strength: the first pattern of each direction gives p = 0.048 (10 of 15), the pooled per-pattern test p = 0.083, equal numbers of bins p = 0.33 (8 of 14), post-sleep strength alone p = 0.45, leaving out the two sessions with fewer than 5 puffs p = 0.15 (8 of 13), and templates from the travel part of each lap only p = 0.19 (9 of 14). The laps from reward times include the time spent at the end where a lap starts: a median of 33% of the air-puff and 43% of the safe template time lies in the travel part, and in the early Rat08 sessions laps last up to 2310 s, so the directional templates are not templates of running alone.
rats-learn-puff-location (blocked)¶
Rats learned the location of an aversive air puff on a linear track.
Blocked as too large: the dandiset holds no position, and tracking the rat in the raw video frames does not fit the time limit (see the data check). No criterion was written for this claim. The lap durations below, from the reward times, are shown as a description only: they say how long the rat took from one end to the other, not where on the track it slowed or stopped, and do not decide an outcome.
In [27]:
# Laps (reward to reward at the other end) in the pre-run, run and post-run epochs, by direction,
# in every file with rewards and a puff direction. A lap belongs to the epoch holding its midpoint.
rows = []
for s in sessions:
if not all(k in s["events"] for k in ("lrw", "rrw", "puf")):
continue
direction, _ = hpcbla.puff_direction(s)
if direction is None:
continue
for a, b, dd in hpcbla.laps(s):
tag = [t for t, x, y in s["epochs"] if x <= (a + b) / 2 < y and t in ("prerun", "run", "postrun")]
rows.append({"session": s["session"], "rat": s["rat"], "epoch": tag[0] if tag else "other",
"direction": "puff" if dd == direction else "safe", "duration": b - a})
all_laps = pd.DataFrame(rows)
print("Laps by epoch:", all_laps.epoch.value_counts().to_dict(), "(laps outside these epochs are left out)")
all_laps = all_laps[all_laps.epoch != "other"]
table = all_laps.groupby(["rat", "session", "epoch", "direction"]).duration.agg(["size", "median"]).unstack(["epoch", "direction"])
print(table.round(1).to_string())
ratio = all_laps.groupby(["session", "epoch", "direction"]).duration.median().unstack("direction")
ratio = (ratio.puff / ratio.safe).unstack("epoch")[["prerun", "run", "postrun"]]
print("\nMedian lap duration, puff direction over safe direction (reference 1):")
for epoch in ratio:
x = ratio[epoch].dropna()
print(f" {epoch}: median ratio {x.median():.2f} over {len(x)} sessions, above 1 in {int((x > 1).sum())},"
f" Wilcoxon of log ratio p = {wilcoxon(np.log(x)):.2g}")
# Air puffs by epoch, in every file with puff events.
puff_epochs = pd.DataFrame([{"session": s["session"], "epoch": next((e for e, x, y in s["epochs"] if x <= t < y), "none")}
for s in sessions if "puf" in s["events"] for t in s["events"]["puf"][1]])
print("\nAir puffs by epoch, over", puff_epochs.session.nunique(), "files:", puff_epochs.epoch.value_counts().to_dict())
print(f"Notebook run time: {(time.time() - notebook_started) / 60:.1f} min")
Laps by epoch: {np.str_('run'): 2008, np.str_('postrun'): 964, np.str_('prerun'): 387, 'other': 50} (laps outside these epochs are left out)
size median
epoch run postrun prerun run postrun prerun
direction puff safe puff safe puff safe puff safe puff safe puff safe
rat session
Rat08 Rat08-20130708 4.0 3.0 NaN NaN NaN NaN 216.5 28.4 NaN NaN NaN NaN
Rat08-20130709 2.0 3.0 2.0 3.0 4.0 5.0 500.2 87.7 124.9 26.9 37.0 24.9
Rat08-20130710 6.0 5.0 NaN NaN 5.0 5.0 191.9 465.4 NaN NaN 17.4 17.0
Rat08-20130711 4.0 4.0 9.0 9.0 4.0 3.0 660.5 140.1 29.2 18.1 31.8 35.4
Rat08-20130712 11.0 11.0 16.0 16.0 6.0 6.0 66.7 52.9 14.7 16.5 16.3 20.1
Rat08-20130713 9.0 9.0 7.0 7.0 5.0 4.0 115.1 79.5 15.5 15.4 17.3 15.8
Rat08-20130715 20.0 21.0 14.0 14.0 5.0 4.0 56.5 29.1 11.8 12.3 10.3 12.3
Rat08-20130716 24.0 24.0 11.0 10.0 4.0 5.0 46.4 24.2 11.8 12.2 15.1 13.4
Rat08-20130717 26.0 26.0 17.0 17.0 5.0 5.0 28.4 18.5 15.9 13.7 14.4 17.0
Rat08-20130718 43.0 44.0 16.0 16.0 11.0 11.0 21.6 14.8 16.0 11.4 13.3 12.8
Rat08-20130719 27.0 28.0 6.0 5.0 6.0 7.0 20.6 26.1 16.1 22.1 10.3 15.6
Rat08-20130720 35.0 35.0 4.0 4.0 6.0 5.0 17.6 21.8 18.5 23.9 10.8 20.4
Rat08-20130722 40.0 40.0 8.0 8.0 4.0 5.0 19.7 15.0 17.7 13.7 27.5 7.1
Rat09 Rat09-20140324 11.0 10.0 4.0 3.0 3.0 3.0 104.7 49.8 24.8 24.8 15.7 16.2
Rat09-20140325 22.0 21.0 4.0 3.0 5.0 4.0 52.2 37.4 25.9 19.1 18.0 18.9
Rat09-20140326 54.0 54.0 23.0 23.0 2.0 2.0 16.6 16.9 13.0 13.7 60.6 26.8
Rat09-20140328 40.0 41.0 16.0 16.0 7.0 6.0 16.1 15.4 13.2 16.2 13.6 11.4
Rat09-20140329 39.0 39.0 21.0 21.0 7.0 7.0 17.2 15.4 11.7 16.2 12.8 15.0
Rat09-20140402 50.0 51.0 20.0 21.0 7.0 6.0 18.0 14.3 12.9 9.9 11.5 12.8
Rat09-20140403 31.0 32.0 9.0 10.0 7.0 6.0 16.3 13.7 8.4 10.3 11.3 14.3
Rat09-20140404 42.0 42.0 11.0 10.0 7.0 7.0 11.9 13.1 15.7 21.0 12.4 10.4
Rat09-20140405 27.0 27.0 24.0 24.0 6.0 5.0 16.4 10.5 13.1 14.0 12.2 13.0
Rat10 Rat10-20140624 19.0 18.0 14.0 13.0 6.0 6.0 37.9 51.8 20.5 20.6 18.3 16.2
Rat10-20140626 17.0 17.0 13.0 12.0 2.0 2.0 32.2 40.3 20.7 19.2 28.2 37.3
Rat10-20140701 13.0 13.0 20.0 21.0 1.0 1.0 54.4 80.3 15.1 15.1 19.7 34.4
Rat10-20140702 34.0 35.0 40.0 40.0 5.0 5.0 15.2 13.2 11.3 10.6 15.3 13.9
Rat10-20140704 25.0 25.0 23.0 23.0 6.0 6.0 19.2 16.7 13.1 11.0 15.5 12.9
Rat10-20140705 46.0 46.0 25.0 25.0 6.0 7.0 13.2 13.4 11.0 14.1 12.6 14.9
Rat10-20140708 41.0 40.0 29.0 28.0 5.0 6.0 13.5 12.1 9.6 11.8 12.5 17.9
Rat11 Rat11-20150313 19.0 19.0 3.0 3.0 5.0 4.0 24.6 16.2 41.3 15.9 18.3 18.6
Rat11-20150315 26.0 27.0 16.0 16.0 6.0 6.0 15.1 14.8 17.7 14.6 14.2 14.6
Rat11-20150316 32.0 33.0 5.0 5.0 6.0 6.0 14.3 11.9 15.0 14.2 13.1 12.8
Rat11-20150325 30.0 30.0 8.0 7.0 5.0 6.0 10.3 9.3 12.4 11.4 16.0 16.1
Rat11-20150326 31.0 30.0 5.0 5.0 6.0 5.0 10.3 13.4 15.4 16.1 14.6 15.3
Rat11-20150327 29.0 28.0 7.0 6.0 5.0 4.0 12.5 9.4 16.7 13.8 16.6 15.4
Rat11-20150328 32.0 31.0 17.0 17.0 5.0 5.0 12.7 9.6 16.0 12.2 14.0 15.8
Rat11-20150331 16.0 16.0 13.0 13.0 5.0 5.0 19.7 17.0 15.9 21.4 15.8 18.4
Rat11-20150401 27.0 26.0 5.0 5.0 6.0 6.0 11.7 16.8 15.5 17.7 18.5 26.1
Median lap duration, puff direction over safe direction (reference 1):
prerun: median ratio 0.95 over 37 sessions, above 1 in 14, Wilcoxon of log ratio p = 0.19
run: median ratio 1.18 over 38 sessions, above 1 in 27, Wilcoxon of log ratio p = 0.0063
postrun: median ratio 1.00 over 36 sessions, above 1 in 19, Wilcoxon of log ratio p = 0.62
Air puffs by epoch, over 38 files: {np.str_('run'): 1195, np.str_('presleep'): 2}
Notebook run time: 9.4 min
In the run, laps toward the air-puff direction take longer than laps the other way (median ratio 1.18 over 38 sessions, above 1 in 27, p = 0.0063); in the pre-run and the post-run, when no puffs are given (1195 of the 1197 puffs in 38 files fall in the run, 2 in the pre-sleep), the two directions take about as long (0.95 and 1.00). Slower laps where puffs are given are what both learning and the puffs themselves would produce, so this does not show that the rats learned where the puff was.
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 |
|---|---|
hpc-bla-reactivation-post-sleep |
reproduced |
reactivation-peaks-at-ripples |
reproduced |
bla-ripple-up-subgroup |
reproduced |
reactivation-involves-ripple-up-cells |
reproduced |
aversive-direction-stronger |
not reproduced (p = 0.11, in the claimed direction) |
rats-learn-puff-location |
blocked: too large (position only in the raw video frames) |
How the analysis differs from the paper's. The full text was not available, so the methods
are common choices fixed before the analysis: hippocampus and BLA as the hpc/HIP and
bla/BLA labels, 50 ms bins, joint patterns from the singular vectors of the run's
hippocampus by BLA cross-correlation matrix with a shuffle threshold, reactivation strength as the
product of the two projections, ripple windows from the rip annotations, laps and the air-puff
direction from the reward and puff times (there is no position), so a lap also holds the time
spent at the end where it starts. The paper's session selection,
cell classification, bin size and statistics are unknown, and 15 sessions from three rats are
used; Rat09 has no hippocampal units or ripples here.
What was found about the dandiset.
- The epochs table is in milliseconds although it says seconds; video timestamps are also in milliseconds, stated as seconds, and are not on the session clock.
- In ten files the units index is wrong at the first unit: unit 1 begins with the end of unit 0 (a fall back in time after a multiple of 65536 spikes).
- Seven files have no units table, sixteen no sleep states and fourteen no ripple events; the Rat09 files have no hippocampal units.
- The unit region labels are undescribed short codes, in two cases for the same region in some
files, and one file has no labels; the electrode table gives every channel as
unknown. - There is no position tracking, only raw video frames, so the behaviour on the track cannot be
measured here; the
trialstable is undescribed.