Reproducing Zaki et al. (2024) from dandiset 000718

Paper. Zaki Y, Pennington ZT, Morales-Rodriguez D, et al. (2024). Offline ensemble co-reactivation links memories across days. Nature 637, 145-155 (issue of January 2025; published online 6 November 2024). doi:10.1038/s41586-024-08168-4. This attempt used the full text, which is open access under CC BY (PMC11666460).

Dandiset. DANDI:000718, version 0.260825.1902. The paper cites no DANDI version (its data statement says the imaging data are available "upon reasonable request"). This is the only published version, so it is pinned; the draft holds the same 20 assets (same asset ids and sizes).

What the dandiset holds. Two of the nine mice of the calcium imaging and EEG/EMG experiment of Figure 5 and Extended Data Figures 9 and 10: Ca_EEG2-1 (high shock, 1.5 mA) and Ca_EEG3-4 (low shock, 0.25 mA). For each, the Miniscope recordings of encoding and recall sessions with Minian output, freezing from ezTrack, and the first 10 min offline recording after aversive encoding (of the 24 the Methods describe per offline period); for Ca_EEG3-4 also the first after neutral encoding. Ca_EEG2-1 has no neutral or novel recall session. Only Ca_EEG2-1 has cell registration (CellReg tables, in its Week file, with a week of telemetry). The mice of the one-hour offline experiments of Figures 2 to 4 and 6, and of the behaviour-only, inhibitory tagging and no-shock experiments, are not in it.

Outside inputs. One claim (source-data-link) compares the dandiset with the paper's per-mouse source data, from the public repository denisecailab/RetrospectiveMemoryLinkingData_2024 at commit a5081645d3d3e5743ade58ea632c28b003794345 (files ExtendedData9/ExtData9b.csv, ExtData9c.csv and Figure5/Fig_5c.csv). The few reference values used are typed in that section. The analysis code the paper names (github.com/denisecailab/RetrospectiveMemoryLinkingAnalysis_2024) returned 404 and could not be read; the conversion code (catalystneuro/cai-lab-to-nwb at commit a37dadafecaec507669c25baa705333460e0f48d) was read for the meaning of the files.

Claim	Tested here
source-data-link: the deposited mice match the paper's per-mouse freezing	tested (with outside data)
ed9b-aversive-freezing-high-vs-low: high-shock mice freeze more in aversive recall	tested, one mouse per group
offline-bursts-above-shuffle: offline bursts exceed shuffled expectation	tested, three offline recordings
offline-skew-above-shuffle: offline population activity is more skewed than shuffled	tested, three offline recordings
burst-locomotion-pause: bursts coincide with a brief pause in movement	tested, with telemetry activity
fig5c-wake-overlap-neutral-coburst: overlap-neutral co-bursting exceeds overlap-aversive in wake, high shock	tested, one mouse, one recording of 24
Behaviour-only experiments (Figure 1b and 1c, Extended Data Figures 1 and 2), three claims	blocked: data missing
One-hour imaging and single-experience experiments (Figure 1d to 1g, Figures 3, 4 and 6, Extended Data Figures 3, 4, 7), eleven claims	blocked: data missing
Inhibitory tagging and no-shock experiments (Extended Data Figures 5, 6, 8), two claims	blocked: data missing
EEG cohort claims needing the missing mice, sessions or registration (Figure 5c sleep and low shock, Extended Data Figures 9c, 9d to 9g, 10), four claims	blocked: data missing

import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
import time
import warnings
from concurrent.futures import ThreadPoolExecutor

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats

import memlink as ml

pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 30)
warnings.filterwarnings("ignore", category=RuntimeWarning)

started = time.time()
SESSIONS = {
    "Ca_EEG2-1": ["NeutralExposure", "FC", "OfflineDay2Session1", "Recall1"],
    "Ca_EEG3-4": ["NeutralExposure", "OfflineDay1Session1", "FC", "OfflineDay2Session1", "Recall1", "Recall2", "Recall3"],
}
keys = [(m, s) for m, ss in SESSIONS.items() for s in ss]
with ThreadPoolExecutor(6) as ex:
    loaded = list(ex.map(lambda k: ml.load_session(ml.session_path(*k)), keys))
data = dict(zip(keys, loaded))
week = ml.load_week("Ca_EEG2-1")
print("NWB files at the pinned version:", len(ml.asset_urls()))
print(f"loaded {len(data)} session files and the Week file of Ca_EEG2-1 from dandiset {ml.DANDISET_ID} "
      f"version {ml.DANDISET_VERSION} in {time.time() - started:.0f} s")

NWB files at the pinned version: 12
loaded 11 session files and the Week file of Ca_EEG2-1 from dandiset 000718 version 0.260825.1902 in 21 s

Data check

Each check of the metadata against the values is stated under its output.

rows = []
for (m, s), d in data.items():
    rows.append({"mouse": m, "session": s, "session_id": d["session_id"], "start": d["start"][:19],
                 "frames": d["S"].shape[0], "rois": d["S"].shape[1], "t_first": round(d["t"][0], 3),
                 "t_last": round(d["t"][-1], 1), "dt_ms": round(1000 * np.median(np.diff(d["t"])), 1),
                 "max_gap_ms": round(1000 * np.diff(d["t"]).max()), "S_min": round(float(d["S"].min()), 2),
                 "S_max": round(float(d["S"].max())), "frac_zero": round(float((d["S"] == 0).mean()), 3),
                 "silent_rois": int((d["S"].max(0) == 0).sum()),
                 "video_frames": d["ops_n"], "video_t_last": round(d["ops_last"], 1),
                 "video_dt_ms": round(1000 * d["ops_dt"], 1)})
info = pd.DataFrame(rows)
print(info.to_string(index=False))
print("\nSubjects:", {m: (data[(m, SESSIONS[m][0])]["species"], data[(m, SESSIONS[m][0])]["strain"],
                       data[(m, SESSIONS[m][0])]["sex"], data[(m, SESSIONS[m][0])]["age"]) for m in SESSIONS})
print("Deconvolved unit:", {d["S_unit"] for d in data.values()}, "; ROI ids are 0..n-1 in every file:",
      all((d["roi_ids"] == np.arange(len(d["roi_ids"]))).all() for d in data.values()))
print("\nSession descriptions:")
for (m, s), d in data.items():
    print(f"  {m} {s}: {d['description'][:150]}")
print("\nFull description of the offline recordings (all three are the same):",
      {data[k]["description"] for k in [("Ca_EEG2-1", "OfflineDay2Session1"), ("Ca_EEG3-4", "OfflineDay1Session1"), ("Ca_EEG3-4", "OfflineDay2Session1")]})
bad = info[(info.video_t_last < 0.9 * info.t_last)]
print("\nFiles whose OnePhotonSeries timestamps end before 90% of the deconvolved timestamps:")
print(bad[["mouse", "session", "video_frames", "video_t_last", "video_dt_ms", "t_last"]].to_string(index=False))
print("\nTelemetry sampling rate in the offline files (Hz):",
      {" ".join(k): data[k]["activity_rate"] for k in data if "activity" in data[k]}, "; Week file:", week["eeg_rate"])

    mouse             session          session_id               start  frames  rois  t_first  t_last  dt_ms  max_gap_ms  S_min  S_max  frac_zero  silent_rois  video_frames  video_t_last  video_dt_ms
Ca_EEG2-1     NeutralExposure     NeutralExposure 2021-10-12T10:44:18    8880   623    0.335   599.9   67.0         906    0.0   2354      0.985            0         17770         296.1         16.7
Ca_EEG2-1                  FC                  FC 2021-10-14T10:11:24    4076   597    0.336   274.0   67.0          70    0.0   2174      0.981            0          8161         274.0         34.0
Ca_EEG2-1 OfflineDay2Session1 OfflineDay2Session1 2021-10-14T09:24:27    8929   625    0.335   599.8   67.0          69    0.0   2346      0.984            0         17868         599.9         34.0
Ca_EEG2-1             Recall1             Recall1 2021-10-15T10:33:19    4463   336    0.335   300.0   67.0          69    0.0   1806      0.987            0          8935         300.0         34.0
Ca_EEG3-4     NeutralExposure     NeutralExposure 2022-09-17T09:23:03    8931   851    0.335   599.9   67.0          69    0.0   1141      0.976            0         17871         599.9         34.0
Ca_EEG3-4 OfflineDay1Session1 OfflineDay1Session1 2022-09-17T09:40:40    8933   711    0.335   600.1   67.0          69    0.0   1099      0.976            0         17876         297.9         16.7
Ca_EEG3-4                  FC                  FC 2022-09-19T09:18:41    4076   783    0.335   274.0   67.0          69    0.0    850      0.969            0          8161         274.0         34.0
Ca_EEG3-4 OfflineDay2Session1 OfflineDay2Session1 2022-09-19T09:30:07    8933   858    0.335   600.1   67.0          69    0.0    907      0.977            0         17876         600.1         34.0
Ca_EEG3-4             Recall1             Recall1 2022-09-20T09:19:43    4463   758    0.335   299.9   67.0          72    0.0    943      0.973            0          8935         299.9         34.0
Ca_EEG3-4             Recall2             Recall2 2022-09-21T08:41:45    4462   831    0.333   299.9   67.0         101    0.0   1054      0.973            0          8934         148.9         16.7
Ca_EEG3-4             Recall3             Recall3 2022-09-22T08:40:12    4463   809    0.336   299.9   67.0          70    0.0   1029      0.972            0          8935         299.9         34.0

Subjects: {'Ca_EEG2-1': ('Mus musculus', 'C57BL/6J', 'M', 'P12W/P18W'), 'Ca_EEG3-4': ('Mus musculus', 'C57BL/6J', 'M', 'P12W/P18W')}
Deconvolved unit: {'n.a.'} ; ROI ids are 0..n-1 in every file: True

Session descriptions:
  Ca_EEG2-1 NeutralExposure: Neutral Exposure session: mouse was exposed to a neutral context for 10 min to explore. Context: overhead external light, external fan off, box fan on
  Ca_EEG2-1 FC: Fear Conditioning session: after a baseline period of 2 min, mouse received three 2s foot shocks of 1.5, with an intershock interval of 1 min. Then, 3
  Ca_EEG2-1 OfflineDay2Session1: After Neutral Exposure and Fear Conditioning sessions, mice were taken out of the testing chambers and immediately placed in their homecage (scope was
  Ca_EEG2-1 Recall1: First Recall session: mouse was placed in shock context for 5 min. Context: overhead external light, external fan at medium level, box fan on, even gr
  Ca_EEG3-4 NeutralExposure: Neutral Exposure session: mouse was exposed to a neutral context for 10 min to explore. Context: overhead external light, external fan off, box fan on
  Ca_EEG3-4 OfflineDay1Session1: After Neutral Exposure and Fear Conditioning sessions, mice were taken out of the testing chambers and immediately placed in their homecage (scope was
  Ca_EEG3-4 FC: Fear Conditioning session: after a baseline period of 2 min, mouse received three 2s foot shocks of 0.25, with an intershock interval of 1 min. Then, 
  Ca_EEG3-4 OfflineDay2Session1: After Neutral Exposure and Fear Conditioning sessions, mice were taken out of the testing chambers and immediately placed in their homecage (scope was
  Ca_EEG3-4 Recall1: First Recall session: mouse was placed in shock context for 5 min. Context: overhead external light, external fan at medium level, box fan on, even gr
  Ca_EEG3-4 Recall2: Second Recall session: mouse was placed in Novel context for 5 min. Context: overhead external light, external fan off, box fan on, bath mat floor, A 
  Ca_EEG3-4 Recall3: Third Recall session: mouse was placed in Neutral context for 5 min. Context: overhead external light, external fan off, box fan on, smooth floor, whi

Full description of the offline recordings (all three are the same): {'After Neutral Exposure and Fear Conditioning sessions, mice were taken out of the testing chambers and immediately placed in their homecage (scope was not removed).The homecage was placed in a dark grey storage bin with a webcam on top of the bin, taped to a wooden plank, looking down into the homecage. Mouse behavior and calcium were recorded for an hour.'}

Files whose OnePhotonSeries timestamps end before 90% of the deconvolved timestamps:
    mouse             session  video_frames  video_t_last  video_dt_ms  t_last
Ca_EEG2-1     NeutralExposure         17770         296.1         16.7   599.9
Ca_EEG3-4 OfflineDay1Session1         17876         297.9         16.7   600.1
Ca_EEG3-4             Recall2          8934         148.9         16.7   299.9

Telemetry sampling rate in the offline files (Hz): {'Ca_EEG2-1 OfflineDay2Session1': 500.0, 'Ca_EEG3-4 OfflineDay1Session1': 500.0, 'Ca_EEG3-4 OfflineDay2Session1': 500.0} ; Week file: 500.0

Deconvolved activity is sampled at 15 Hz (67 ms), as the Methods state for the EEG/EMG experiment (videos downsampled by 2 from 30 Hz), with one gap of 0.9 s in Ca_EEG2-1 NeutralExposure. Its unit is "n.a." and its values are non-negative with about 97 to 99% zeros, as expected of deconvolved activity. Sessions last 10 min (neutral exposure and offline), 5 min (recall) and 274 s (FC: 2 min baseline, three shocks 1 min apart, 30 s after the last), as the Methods describe. In three files the OnePhotonSeries timestamps run at 16.7 ms and end at half the recording's length, while the frame count fits 30 Hz over the full length; these raw-video timestamps are wrong (a finding). They are not used here. Subject metadata give age as a range (P12W/P18W), as the Methods' 12 to 18 weeks. The telemetry (EEG, EMG, activity) is sampled at 500 Hz in every offline file and in the Week file, while the Methods say the HD-X02 implants recorded at 100 Hz. The files are trusted here, and the difference is recorded as a finding.

# Session order and clocks. Behaviour sessions take their start time from the lab's session-time
# sheet; offline recordings from the Miniscope metadata (conversion code). The Week file of
# Ca_EEG2-1 places every session on its telemetry clock.
ep = pd.DataFrame([(k, a, b) for k, (a, b) in week["epochs"].items()], columns=["session", "start_s", "stop_s"])
ep = ep.sort_values("start_s")
print("Week file start:", week["start"], "; telemetry", week["eeg_n"], "samples at", week["eeg_rate"], "Hz =",
      round(week["eeg_n"] / week["eeg_rate"] / 86400, 2), "days")
print("Epoch table of the Week file, all", len(ep), "rows, in order of start time:")
print(ep.to_string(index=False))

# Where is the offline file's EEG in the week recording? Exact match of its first samples.
off = data[("Ca_EEG2-1", "OfflineDay2Session1")]
a0 = week["epochs"]["OfflineDay2Session1"][0]
seg = ml.week_signal(a0 - 5, a0 + 15)
head = off["eeg_head"][:2000]
hits = [i for i in range(len(seg) - len(head)) if np.allclose(seg[i:i + len(head)], head)]
print("\nOffline-day-2 session 1 EEG found in the Week EEG at", [round(a0 - 5 + i / week["eeg_rate"], 3) for i in hits],
      "s; epoch start", round(a0, 3), "s")

# The telemetry drops to zero while the mouse is out of its home cage (in a behaviour chamber).
def dropouts(t0, t1, step=10.0):
    x = ml.week_signal(t0, t1)
    n = int(step * week["eeg_rate"])
    z = (x[: len(x) // n * n].reshape(-1, n) == 0).mean(1) > 0.9
    edges = np.diff(np.concatenate([[0], z.astype(int), [0]]))
    return [(t0 + i * step, t0 + j * step) for i, j in zip(np.where(edges == 1)[0], np.where(edges == -1)[0])]

for name in ["NeutralExposure", "FC", "Recall1"]:
    a, b = week["epochs"][name]
    drops = [(round(x), round(y)) for x, y in dropouts(a - 7200, a + 3600) if y - x >= 60]
    print(f"{name}: epoch {a:.0f} to {b:.0f} s; telemetry absent (>= 60 s, within -2 h to +1 h) at {drops}; "
          f"epoch moved 1 h earlier lies inside an absence: {any(x <= a - 3600 and b - 3600 <= y for x, y in drops)}")

# Ca_EEG3-4: gap between the end of each encoding session and the start of the offline recording after it.
def start_dt(m, s):
    return pd.Timestamp(data[(m, s)]["start"])
for enc, offl in [("NeutralExposure", "OfflineDay1Session1"), ("FC", "OfflineDay2Session1")]:
    gap = (start_dt("Ca_EEG3-4", offl) - start_dt("Ca_EEG3-4", enc)).total_seconds() - data[("Ca_EEG3-4", enc)]["t"][-1]
    print(f"Ca_EEG3-4: {offl} starts {gap / 60:.1f} min after the end of {enc}")
gap = (start_dt("Ca_EEG2-1", "OfflineDay2Session1") - start_dt("Ca_EEG2-1", "FC")).total_seconds()
print(f"Ca_EEG2-1: OfflineDay2Session1 starts {gap / 60:.1f} min after the start of FC, by the files' start times")

Week file start: 2021-10-08T07:00:00+00:00 ; telemetry 404395000 samples at 500.0 Hz = 9.36 days
Epoch table of the Week file, all 27 rows, in order of start time:
             session    start_s     stop_s
 OfflineDay1Session1 356757.595 356813.764
 OfflineDay1Session2 356897.828 357005.289
 OfflineDay1Session3 357044.154 357312.538
 OfflineDay1Session4 357350.761 357417.826
     NeutralExposure 359058.744 359658.717
 OfflineDay1Session5 364070.226 364370.844
 OfflineDay2Session1 527067.677 527667.538
 OfflineDay2Session2 528867.463 529467.359
                  FC 529884.777 530158.743
 OfflineDay2Session3 530667.282 531267.178
 OfflineDay2Session4 532467.102 533066.997
 OfflineDay2Session5 534266.922 534866.817
 OfflineDay2Session6 536066.738 536666.632
 OfflineDay2Session7 537866.553 538466.447
 OfflineDay2Session8 539666.372 540266.265
 OfflineDay2Session9 541466.187 542066.080
OfflineDay2Session10 543266.005 543865.900
OfflineDay2Session11 545065.825 545665.720
OfflineDay2Session12 546865.644 547465.538
OfflineDay2Session13 548665.430 549265.358
OfflineDay2Session14 550465.248 551065.177
OfflineDay2Session15 552265.067 552864.995
OfflineDay2Session16 554064.885 554664.780
OfflineDay2Session17 555864.703 556464.597
OfflineDay2Session18 557664.522 558264.416
OfflineDay2Session19 559464.340 560064.236
             Recall1 617599.999 617899.952

Offline-day-2 session 1 EEG found in the Week EEG at [np.float64(527067.679)] s; epoch start 527067.677 s

NeutralExposure: epoch 359059 to 359659 s; telemetry absent (>= 60 s, within -2 h to +1 h) at [(355279, 356499)]; epoch moved 1 h earlier lies inside an absence: True

FC: epoch 529885 to 530159 s; telemetry absent (>= 60 s, within -2 h to +1 h) at [(526175, 526955)]; epoch moved 1 h earlier lies inside an absence: True

Recall1: epoch 617600 to 617900 s; telemetry absent (>= 60 s, within -2 h to +1 h) at [(613870, 614750)]; epoch moved 1 h earlier lies inside an absence: True
Ca_EEG3-4: OfflineDay1Session1 starts 7.6 min after the end of NeutralExposure
Ca_EEG3-4: OfflineDay2Session1 starts 6.9 min after the end of FC
Ca_EEG2-1: OfflineDay2Session1 starts -47.0 min after the start of FC, by the files' start times

The Week file's epoch table has 27 rows: neutral exposure, FC, Recall1, five offline-day-1 recordings (the first four last 56 to 268 s) and 19 offline-day-2 recordings, against the 24 per offline period of the Methods. It puts the first offline recording after aversive encoding (OfflineDay2Session1, starting 527068 s) 47 min before the FC session (529885 s), and the session start times in the files say the same (09:24 against 10:11). The telemetry shows otherwise: the EEG of the offline file sits exactly at the epoch's start (527067.679 against 527067.677 s), and the telemetry is absent, as when the mouse is out of its home cage, from 526175 to 526955 s, which contains the FC epoch moved 1 h earlier and ends 113 s before that offline recording. The next offline recording starts 1800 s after it, and the others follow every 1800 s (10 min of recording and a 20 min break, as the Methods describe), so no offline-day-2 recording lies between the FC session, read 1 h earlier, and OfflineDay2Session1. The output above gives the same one-hour lead for the neutral exposure and Recall1 epochs. The data are consistent with the behaviour-session times of Ca_EEG2-1 being about 1 h late, and with the offline recording being the first after FC, as its name says; it is treated so here, and the start times are recorded as a finding. If this order were wrong, the fig5c analysis below would not test the claim. For Ca_EEG3-4 the offline recordings start 7.6 and 6.9 min after the end of the session they follow.

# Behaviour: freezing intervals against motion, and shocks against motion.
rows = []
for (m, s), d in data.items():
    if "motion" not in d:
        continue
    inside = ml.in_intervals(d["motion_t"], d["freeze"])
    rows.append({"mouse": m, "session": s, "motion_frames": len(d["motion"]), "motion_s": round(d["motion_t"][-1], 1),
                 "motion_range": f"{d['motion'].min():.0f}..{d['motion'].max():.0f}", "intervals": len(d["freeze"]),
                 "min_interval_s": round(float(np.min(d['freeze'][:, 1] - d['freeze'][:, 0])), 2),
                 "median_motion_in_freeze": np.median(d["motion"][inside]), "median_motion_outside": np.median(d["motion"][~inside])})
print(pd.DataFrame(rows).to_string(index=False))
print("\nFreezing interval parameters (Ca_EEG2-1 FC):", " ".join(data[("Ca_EEG2-1", "FC")]["freeze_description"].split())[:300])
print("\nShocks (s) and the motion response: mean motion in the 2 s shock against the 10 s before, and against 2 s windows at 1000 random times")
rng_check = np.random.default_rng(1)
for m in SESSIONS:
    d = data[(m, "FC")]
    shock = [d["motion"][(d["motion_t"] >= a) & (d["motion_t"] < b)].mean() for a, b in d["shocks"]]
    pre = [d["motion"][(d["motion_t"] >= a - 10) & (d["motion_t"] < a)].mean() for a, b in d["shocks"]]
    rnd = [d["motion"][(d["motion_t"] >= r) & (d["motion_t"] < r + 2)].mean() for r in rng_check.uniform(0, 270, 1000)]
    print(f"  {m}: shocks {d['shocks'][:, 0].tolist()} amplitude {d['shock_amplitude'].tolist()} mA; "
          f"motion during shocks {np.round(shock)}, before {np.round(pre)}, random windows 95th percentile {np.percentile(rnd, 95):.0f}")

    mouse         session  motion_frames  motion_s motion_range  intervals  min_interval_s  median_motion_in_freeze  median_motion_outside
Ca_EEG2-1 NeutralExposure          17996     599.8      0..4568         43             0.5                     32.0                  821.0
Ca_EEG2-1              FC           8217     273.9      0..7547         22             0.5                     24.0                  756.0
Ca_EEG2-1         Recall1           8998     299.9      0..1530         58             0.5                      0.0                  257.0
Ca_EEG3-4 NeutralExposure          17999     599.9      0..4880         56             0.5                      9.0                  714.0
Ca_EEG3-4              FC           8217     273.9      0..4021         21             0.5                      6.0                  780.0
Ca_EEG3-4         Recall1           8999     299.9      0..3711         33             0.5                     12.5                  693.0
Ca_EEG3-4         Recall2           8997     299.9      0..4426         41             0.5                      3.0                  536.0
Ca_EEG3-4         Recall3           8998     299.9      0..4709         45             0.5                     10.0                  524.0

Freezing interval parameters (Ca_EEG2-1 FC): Freezing behavior intervals generated using EzTrack software for file Ca_EEG2-1_FC.wmv. Parameters used include a motion cutoff of 10.0, freeze threshold of 200.0, and a minimum freeze duration of 15.0. - Freeze threshold: The maximum amount of motion the animal can display while still being classif

Shocks (s) and the motion response: mean motion in the 2 s shock against the 10 s before, and against 2 s windows at 1000 random times
  Ca_EEG2-1: shocks [120.0, 180.0, 240.0] amplitude [1.5, 1.5, 1.5] mA; motion during shocks [3387. 3644. 3848.], before [562. 645. 259.], random windows 95th percentile 1878
  Ca_EEG3-4: shocks [120.0, 180.0, 240.0] amplitude [0.25, 0.25, 0.25] mA; motion during shocks [733. 366. 841.], before [538. 570. 590.], random windows 95th percentile 1537

Motion inside the freezing intervals is near zero and outside them in the hundreds, so the intervals agree with the motion they were computed from. The shock table gives 1.5 mA for Ca_EEG2-1 and 0.25 mA for Ca_EEG3-4 at 120, 180 and 240 s, as the session descriptions and the Methods say. Motion during the high shocks of Ca_EEG2-1 (3387 to 3848) is above the 95th percentile of random 2 s windows (1878), consistent with the shock times sharing the behaviour clock. The low shocks of Ca_EEG3-4 give motion within the range of random windows (366 to 841 against 1537), so for that session the shock times could not be checked this way.

# Cell registration (Ca_EEG2-1 only).
reg = week["registration"]
print("Registration tables in the Week file:", len(week["tables"]), "; table used:", ml.REGISTRATION_TABLE,
      "; rows", len(reg["FC"]))
size = {"FC": data[("Ca_EEG2-1", "FC")]["S"].shape[1], "NeutralExposure": data[("Ca_EEG2-1", "NeutralExposure")]["S"].shape[1],
        "OfflineDay2Session1": off["S"].shape[1], "Recall1": data[("Ca_EEG2-1", "Recall1")]["S"].shape[1]}
for col, v in reg.items():
    ok = v[v >= 0]
    print(f"  {col}: values {v.min()}..{v.max()} (-9999 = no match), matched {len(ok)} of {size[col]} ROIs, "
          f"each ROI once: {len(np.unique(ok)) == len(ok)}, every ROI present: {len(np.unique(ok)) == size[col]}")
rng_check = np.random.default_rng(2)
oc = reg["OfflineDay2Session1"]
for col in ["NeutralExposure", "FC", "Recall1"]:
    msk = (oc >= 0) & (reg[col] >= 0)
    dist = np.linalg.norm(off["centroids"][oc[msk]] - data[("Ca_EEG2-1", col)]["centroids"][reg[col][msk]], axis=1)
    rnd = np.linalg.norm(off["centroids"][rng_check.choice(oc[oc >= 0], msk.sum())] - data[("Ca_EEG2-1", col)]["centroids"][reg[col][msk]], axis=1)
    print(f"  offline-{col} matches: {msk.sum()}; centroid distance median {np.median(dist):.1f} px "
          f"(95th pct {np.percentile(dist, 95):.1f}); random pairs median {np.median(rnd):.0f} px")
labels = ml.ensembles(reg, "OfflineDay2Session1", off["S"].shape[1])
print("\nEnsembles of the offline recording:", pd.Series(labels).value_counts().to_dict())
print("Ca_EEG3-4 has no Week file and no registration table:",
      not any("EEG3-4" in p and "Week" in p for p in ml.asset_urls()))

Registration tables in the Week file: 19 ; table used: OfflineDay2Session1vsConditioningSessions ; rows 1395
  FC: values -9999..596 (-9999 = no match), matched 597 of 597 ROIs, each ROI once: True, every ROI present: True
  NeutralExposure: values -9999..622 (-9999 = no match), matched 623 of 623 ROIs, each ROI once: True, every ROI present: True
  OfflineDay2Session1: values -9999..624 (-9999 = no match), matched 625 of 625 ROIs, each ROI once: True, every ROI present: True
  Recall1: values -9999..335 (-9999 = no match), matched 336 of 336 ROIs, each ROI once: True, every ROI present: True
  offline-NeutralExposure matches: 263; centroid distance median 3.2 px (95th pct 5.8); random pairs median 196 px
  offline-FC matches: 246; centroid distance median 2.0 px (95th pct 4.4); random pairs median 185 px
  offline-Recall1 matches: 143; centroid distance median 3.2 px (95th pct 6.3); random pairs median 169 px

Ensembles of the offline recording: {'remaining': 282, 'overlap': 166, 'neutral': 97, 'aversive': 80}
Ca_EEG3-4 has no Week file and no registration table: True

Every ROI of each of the four sessions appears once in the table, and matched ROIs lie 2 to 3 px apart against 169 to 196 px for random pairs, so the table is consistent with the footprints it was computed from. The offline recording of Ca_EEG2-1 splits into the counts of neutral, aversive, overlap and remaining cells printed above. The dandiset holds 19 tables for offline day 2 but only the first of those recordings, and no registration for Ca_EEG3-4.

# Shared clock between calcium and the other streams: decode each stream from the binned activity
# with the stream shifted in time; a peak at zero shift means a shared clock.
shifts = list(range(-20, 21, 2))
rows = []
for (m, s), d in data.items():
    tests = {}
    if "motion" in d:
        tests["motion"] = (d["motion"].astype(float), d["motion_t"])
    if "activity" in d:
        mt, mv = ml.movement(d)
        tests["movement"] = (mv, mt)
        st = ml.state_of_frames(d["t"], d["sleep"])
        tests["wake_label"] = ((st == "wake").astype(float), d["t"])
    for name, (y, ty) in tests.items():
        r = ml.shift_decoding(d["S"], d["t"], y, ty, shifts)
        best = max(r, key=r.get)
        rows.append({"mouse": m, "session": s, "stream": name, "r_at_0": round(r[0], 3), "best_shift_s": best,
                     "r_best": round(r[best], 3), "max_r_beyond_4s": round(max(v for k, v in r.items() if abs(k) > 4), 3)})
clock = pd.DataFrame(rows)
print(clock.to_string(index=False))

    mouse             session     stream  r_at_0  best_shift_s  r_best  max_r_beyond_4s
Ca_EEG2-1     NeutralExposure     motion   0.308             0   0.308           -0.014
Ca_EEG2-1                  FC     motion   0.166             0   0.166            0.077
Ca_EEG2-1 OfflineDay2Session1   movement   0.228             0   0.228            0.122
Ca_EEG2-1 OfflineDay2Session1 wake_label   0.241             2   0.242            0.103
Ca_EEG2-1             Recall1     motion   0.449             0   0.449           -0.068
Ca_EEG3-4     NeutralExposure     motion   0.508             0   0.508            0.126
Ca_EEG3-4 OfflineDay1Session1   movement   0.063            -2   0.159            0.132
Ca_EEG3-4 OfflineDay1Session1 wake_label   0.182             0   0.182            0.148
Ca_EEG3-4                  FC     motion   0.471             0   0.471            0.029
Ca_EEG3-4 OfflineDay2Session1   movement   0.244             2   0.255            0.104
Ca_EEG3-4 OfflineDay2Session1 wake_label   0.187             0   0.187            0.119
Ca_EEG3-4             Recall1     motion   0.469             0   0.469            0.265
Ca_EEG3-4             Recall2     motion   0.526             0   0.526            0.069
Ca_EEG3-4             Recall3     motion   0.427            -2   0.427            0.034

For every behaviour session, motion decoded from the calcium activity is best at zero shift, or equally good at the neighbouring 2 s step (Ca_EEG3-4 Recall3), and lower beyond 4 s, so imaging and behaviour are consistent with a shared clock. Telemetry movement peaks at 0 or 2 s in the offline recordings of Ca_EEG2-1 and of Ca_EEG3-4 day 2. In Ca_EEG3-4 day 1 it is weak (r = 0.06 at zero, 0.16 at -2 s, 0.13 beyond 4 s), so a shared clock of its telemetry could not be established this way. The wake label (6 s epochs alternating with quiet wake) peaks within 2 s of zero in all three offline recordings, with a small margin over other shifts in Ca_EEG3-4 day 1 (0.18 against 0.15).

# Recordings deposited more than once: frame counts, maximum projections, telemetry.
keys_ = list(data)
pairs = []
for i in range(len(keys_)):
    for j in range(i + 1, len(keys_)):
        a, b = data[keys_[i]], data[keys_[j]]
        same_shape = a["S"].shape == b["S"].shape
        r = np.corrcoef(a["max_proj"].ravel(), b["max_proj"].ravel())[0, 1] if a["max_proj"].shape == b["max_proj"].shape else np.nan
        same_eeg = ("eeg_head" in a and "eeg_head" in b and np.array_equal(a["eeg_head"], b["eeg_head"]))
        pairs.append({"a": " ".join(keys_[i]), "b": " ".join(keys_[j]), "same_S_shape": same_shape,
                      "max_proj_r": round(r, 3), "same_eeg": same_eeg})
pairs = pd.DataFrame(pairs)
print("Pairs with the same trace shape or identical EEG:", int((pairs.same_S_shape | pairs.same_eeg).sum()))
print("Highest maximum-projection correlations:")
print(pairs.sort_values("max_proj_r", ascending=False).head(6).to_string(index=False))

Pairs with the same trace shape or identical EEG: 0
Highest maximum-projection correlations:
                        a                             b  same_S_shape  max_proj_r  same_eeg
Ca_EEG3-4 NeutralExposure Ca_EEG3-4 OfflineDay1Session1         False       0.863     False
             Ca_EEG3-4 FC Ca_EEG3-4 OfflineDay2Session1         False       0.853     False
        Ca_EEG3-4 Recall1             Ca_EEG3-4 Recall2         False       0.845     False
             Ca_EEG2-1 FC Ca_EEG2-1 OfflineDay2Session1         False       0.833     False
Ca_EEG3-4 NeutralExposure Ca_EEG3-4 OfflineDay2Session1         False       0.830     False
Ca_EEG2-1 NeutralExposure Ca_EEG2-1 OfflineDay2Session1         False       0.825     False

No two sessions share a trace shape or telemetry, so no recording is deposited twice. Maximum projections of one mouse correlate across sessions (up to 0.86), as expected of one field of view imaged through a fixed baseplate; the ROI counts differ, so these are separate recordings.

# Sleep states in the offline recordings.
for m, s in [("Ca_EEG2-1", "OfflineDay2Session1"), ("Ca_EEG3-4", "OfflineDay1Session1"), ("Ca_EEG3-4", "OfflineDay2Session1")]:
    st, sp, lab = data[(m, s)]["sleep"]
    dur = pd.Series(sp - st).groupby(lab).sum().round(1).to_dict()
    frames = pd.Series(ml.state_of_frames(data[(m, s)]["t"], data[(m, s)]["sleep"])).value_counts().to_dict()
    print(f"{m} {s}: intervals {len(st)}, {st.min():.1f} to {sp.max():.1f} s; seconds per state {dur}; calcium frames per state {frames}")
print("\nSleepIntervals description:", ml.open_nwb(ml.session_path("Ca_EEG2-1", "OfflineDay2Session1"))["processing/sleep/SleepIntervals"].attrs["description"])
print("sleep_state column description:", " ".join(ml.open_nwb(ml.session_path("Ca_EEG2-1", "OfflineDay2Session1"))["processing/sleep/SleepIntervals/sleep_state"].attrs["description"].split()))

Ca_EEG2-1 OfflineDay2Session1: intervals 26, 0.0 to 599.3 s; seconds per state {'quiet wake': 482.9, 'wake': 114.7}; calcium frames per state {np.str_('quiet wake'): 7203, np.str_('wake'): 1719, 'none': 7}
Ca_EEG3-4 OfflineDay1Session1: intervals 14, 0.0 to 599.3 s; seconds per state {'quiet wake': 480.3, 'wake': 118.1}; calcium frames per state {np.str_('quiet wake'): 7155, np.str_('wake'): 1767, 'none': 11}
Ca_EEG3-4 OfflineDay2Session1: intervals 36, 0.0 to 599.3 s; seconds per state {'quiet wake': 348.2, 'wake': 248.7}; calcium frames per state {np.str_('quiet wake'): 5205, np.str_('wake'): 3717, 'none': 11}

SleepIntervals description: Sleep states classified with custom algorithm using the data from the HD-X02 sensor (EEG, EMG, temperature, etc.).

sleep_state column description: Sleep State Classification, it can be one of the following: - 'quiet wake': The animal is awake and alert but at rest, without engaging in any voluntary movement or significant cognitive activity - 'rem': Sleep stage characterized by Rapid Eye Movements (REM) beneath closed eyelids. - 'sws': Slow-Wave Sleep (SWS) is a stage of non-REM (NREM) sleep characterized by slow, high-amplitude brain waves - 'wake': State of full consciousness when the animal is alert, responsive to the environment, and capable of voluntary movement

The three offline recordings hold only wake and quiet wake, in 6 s epochs (the Methods' scoring epoch), with quiet wake in 58 to 81% of the time. The file describes the states as classified "with custom algorithm" and defines quiet wake as awake and at rest. The paper's states are wake, NREM and REM, and its scoring has no quiet-wake class, so the criterion of fig5c fixed two readings before the analysis. The next cell applies the Methods' own scoring to the EMG of Ca_EEG2-1 to see where quiet wake falls. No NREM or REM falls in these first recordings, consistent with the mice staying awake for about 30 min after encoding (Extended Data Figure 10c).

# The Methods' sleep scoring applied to Ca_EEG2-1: 6 s epochs, two-state Gaussian mixture on EMG,
# fitted over the 12 h offline period that starts with OfflineDay2Session1 (telemetry clock of the
# Week file). Which class do the deposited `wake` and `quiet wake` epochs fall in?
a_off = week["epochs"]["OfflineDay2Session1"][0]
emg12 = ml.week_signal(a_off, a_off + 12 * 3600, "EMGSignal")
score = {"log_rms_12h": ml.emg_wake_gmm(emg12, week["eeg_rate"]),
         "rms_12h": ml.emg_wake_gmm(emg12, week["eeg_rate"], log=False),
         "log_rms_first_6h": ml.emg_wake_gmm(emg12[: len(emg12) // 2], week["eeg_rate"])}
t_off = data[("Ca_EEG2-1", "OfflineDay2Session1")]["t"]
lab_off = ml.state_of_frames(t_off, data[("Ca_EEG2-1", "OfflineDay2Session1")]["sleep"])
rows = []
for name, sc in score.items():
    gmm_frames = sc[np.minimum((t_off // ml.SCORE_EPOCH_S).astype(int), len(sc) - 1)]
    for lab_ in ("wake", "quiet wake"):
        sel = lab_off == lab_
        rows.append({"scoring": name, "deposited_label": lab_, "frames": int(sel.sum()),
                     "pct_scored_wake": round(100 * (gmm_frames[sel] == 1).mean(), 1),
                     "pct_scored_sleep": round(100 * (gmm_frames[sel] == 0).mean(), 1)})
    print(f"{name}: {(sc == 1).sum()} of {len(sc)} epochs ({100 * (sc == 1).mean():.1f}%) scored wake; epochs without telemetry: {(sc == -1).sum()}")
print(pd.DataFrame(rows).to_string(index=False))
gmm_wake_frames = score["log_rms_12h"][np.minimum((t_off // ml.SCORE_EPOCH_S).astype(int), len(score["log_rms_12h"]) - 1)] == 1

# ActivitySignal: description, distribution, relation to EMG and to the deposited labels (6 s epochs).
print("\nActivitySignal description:", " ".join(ml.open_nwb(ml.session_path("Ca_EEG2-1", "OfflineDay2Session1"))["acquisition/ActivitySignal"].attrs["description"].split()))
for k in [("Ca_EEG2-1", "OfflineDay2Session1"), ("Ca_EEG3-4", "OfflineDay1Session1"), ("Ca_EEG3-4", "OfflineDay2Session1")]:
    d = data[k]
    f = ml.open_nwb(d["path"])
    emg = f["acquisition/EMGSignal/data"][:]
    act = d["activity"]
    n = int(ml.SCORE_EPOCH_S * d["activity_rate"])
    m = min(len(act), len(emg)) // n
    act_e = np.abs(act[: m * n]).reshape(m, n).mean(1)
    emg_e = emg[: m * n].reshape(m, n).std(1)
    mid = (np.arange(m) + 0.5) * ml.SCORE_EPOCH_S + d["activity_t0"]
    lab_e = ml.state_of_frames(mid, d["sleep"])
    q = np.percentile(act, [0, 1, 50, 99, 100])
    print(f"{' '.join(k)}: ActivitySignal unit {f['acquisition/ActivitySignal/data'].attrs['unit']}, "
          f"percentiles 0/1/50/99/100: {np.round(q, 3).tolist()}; "
          f"Spearman r of |activity| with EMG SD over {m} epochs {stats.spearmanr(act_e, emg_e).statistic:.2f}; "
          f"mean |activity| wake {act_e[lab_e == 'wake'].mean():.3f} (n={int((lab_e == 'wake').sum())}), "
          f"quiet wake {act_e[lab_e == 'quiet wake'].mean():.3f} (n={int((lab_e == 'quiet wake').sum())}); "
          f"EMG SD wake {emg_e[lab_e == 'wake'].mean():.2e}, quiet wake {emg_e[lab_e == 'quiet wake'].mean():.2e} V")

log_rms_12h: 2637 of 7200 epochs (36.6%) scored wake; epochs without telemetry: 0
rms_12h: 0 of 7200 epochs (0.0%) scored wake; epochs without telemetry: 0
log_rms_first_6h: 894 of 3600 epochs (24.8%) scored wake; epochs without telemetry: 0
         scoring deposited_label  frames  pct_scored_wake  pct_scored_sleep
     log_rms_12h            wake    1719             98.5               1.5
     log_rms_12h      quiet wake    7203             94.1               5.9
         rms_12h            wake    1719              0.0             100.0
         rms_12h      quiet wake    7203              0.0             100.0
log_rms_first_6h            wake    1719             98.5               1.5
log_rms_first_6h      quiet wake    7203             97.9               2.1

ActivitySignal description: Activity signal recorder with HD-X02 wireless telemetry probe. It refers to the motion of the probe relative to the receiver and it can be used as a proxy for locomotion.

Ca_EEG2-1 OfflineDay2Session1: ActivitySignal unit n.a., percentiles 0/1/50/99/100: [-0.913, -0.269, 0.002, 1.918, 4.001]; Spearman r of |activity| with EMG SD over 99 epochs 0.26; mean |activity| wake 0.322 (n=19), quiet wake 0.080 (n=80); EMG SD wake 5.63e-05, quiet wake 5.03e-05 V

Ca_EEG3-4 OfflineDay1Session1: ActivitySignal unit n.a., percentiles 0/1/50/99/100: [-0.695, -0.285, 0.0, 1.183, 2.145]; Spearman r of |activity| with EMG SD over 100 epochs 0.55; mean |activity| wake 0.349 (n=20), quiet wake 0.023 (n=80); EMG SD wake 7.77e-05, quiet wake 5.80e-05 V

Ca_EEG3-4 OfflineDay2Session1: ActivitySignal unit n.a., percentiles 0/1/50/99/100: [-0.935, -0.349, 0.005, 2.03, 5.004]; Spearman r of |activity| with EMG SD over 100 epochs 0.11; mean |activity| wake 0.419 (n=42), quiet wake 0.055 (n=58); EMG SD wake 7.55e-05, quiet wake 6.89e-05 V

Scored as the Methods describe (6 s epochs, a two-state Gaussian mixture on EMG, the higher state wake), with the log EMG RMS of each epoch over the 12 h after aversive encoding, 36.6% of epochs are wake, and 94.1% of the frames labelled quiet wake in the offline recording fall in the wake state (98.5% of those labelled wake). Fitted on the first 6 h only, 97.9% of quiet-wake frames are wake. Fitted on the RMS without the log, the mixture puts every epoch of the 12 h in one state, so it does not separate sleep from wake. Under the Methods' scoring, the deposited quiet wake is wake. The data support this mapping, but it rests on an EMG measure (log RMS) chosen here after seeing the data, as the paper does not name its EMG feature; it is consistent with the data, not shown.

The ActivitySignal is described in the files as the motion of the probe relative to the receiver, "a proxy for locomotion"; the Methods call the same channel locomotion. Its mean magnitude is 4 to 15 times higher in wake epochs than in quiet wake epochs in all three recordings (0.322 against 0.080, 0.349 against 0.023, 0.419 against 0.055), and it correlates with the epoch's EMG SD at Spearman r = 0.26, 0.55 and 0.11. This is consistent with it tracking movement, at the coarse scale of 6 s epochs; it does not show that it resolves the brief slowing of Figure 2d.

source-data-link

The two deposited mice are the mice of the same names in the paper's per-mouse source data, with the same shock group, and the freezing computed from the dandiset's freezing intervals agrees with the per-mouse freezing the paper reports for them.

Reference values typed from the source-data repository at commit a508164 (ExtData9b.csv, ExtData9c.csv, Fig_5c.csv). The Figure 5c values are over all 24 offline recordings and are shown for context only.

SOURCE_GROUP = {"Ca_EEG2-1": 1.5, "Ca_EEG3-4": 0.25}            # ExtData9b/9c, Amplitude column
SOURCE_FREEZING = [("Ca_EEG2-1", "Recall1", "aversive", 63.86974883307402),   # ExtData9b.csv
                   ("Ca_EEG3-4", "Recall1", "aversive", 8.456495166129569),   # ExtData9b.csv
                   ("Ca_EEG3-4", "Recall3", "neutral", 16.903756390308956),   # ExtData9c.csv
                   ("Ca_EEG3-4", "Recall2", "novel", 14.60486828942981)]      # ExtData9c.csv
for m in SESSIONS:
    sid = {data[(m, s)]["subject_id"] for s in SESSIONS[m]}
    amp = data[(m, "FC")]["shock_amplitude"]
    print(f"{m}: subject_id in files {sid}; shock amplitude {sorted(set(amp.tolist()))} mA; source data group {SOURCE_GROUP[m]} mA")
rows = []
for m, s, ctx, ref in SOURCE_FREEZING:
    d = data[(m, s)]
    a, b = ml.freezing_percent(d)
    rows.append({"mouse": m, "session": s, "context_in_description": d["description"].split(":")[1].split(".")[0].strip()[:45],
                 "paper_context": ctx, "paper": round(ref, 2), "tolerance": round(0.1 * ref, 2),
                 "reading_a": round(a, 2), "reading_b": round(b, 2),
                 "a_within": abs(a - ref) <= 0.1 * ref, "b_within": abs(b - ref) <= 0.1 * ref})
link = pd.DataFrame(rows)
print(link.to_string(index=False))
print("Within tolerance: reading (a)", int(link.a_within.sum()), "of 4; reading (b)", int(link.b_within.sum()), "of 4")

Ca_EEG2-1: subject_id in files {'Ca_EEG2-1'}; shock amplitude [1.5] mA; source data group 1.5 mA
Ca_EEG3-4: subject_id in files {'Ca_EEG3-4'}; shock amplitude [0.25] mA; source data group 0.25 mA
    mouse session                        context_in_description paper_context  paper  tolerance  reading_a  reading_b  a_within  b_within
Ca_EEG2-1 Recall1   mouse was placed in shock context for 5 min      aversive  63.87       6.39      64.09      63.45      True      True
Ca_EEG3-4 Recall1   mouse was placed in shock context for 5 min      aversive   8.46       0.85       8.82       8.46      True      True
Ca_EEG3-4 Recall3 mouse was placed in Neutral context for 5 min       neutral  16.90       1.69      17.40      16.91      True      True
Ca_EEG3-4 Recall2   mouse was placed in Novel context for 5 min         novel  14.60       1.46      15.06      14.61      True      True
Within tolerance: reading (a) 4 of 4; reading (b) 4 of 4

# Other freezing values of the deposited sessions, for reference (no paper value to compare).
for (m, s), d in data.items():
    if "freeze" in d and (m, s) not in [(x[0], x[1]) for x in SOURCE_FREEZING]:
        a, b = ml.freezing_percent(d)
        print(f"{m} {s}: freezing (a) {a:.2f}%, (b) {b:.2f}%")
print("\nContext only, Fig_5c.csv wake co-bursting over all 24 recordings (%):",
      {"Ca_EEG2-1": {"OxN": 10.37, "OxA": 8.39}, "Ca_EEG3-4": {"OxN": 14.35, "OxA": 18.99}})

Ca_EEG2-1 NeutralExposure: freezing (a) 6.18%, (b) 5.95%
Ca_EEG2-1 FC: freezing (a) 11.99%, (b) 11.72%
Ca_EEG3-4 NeutralExposure: freezing (a) 10.41%, (b) 10.10%
Ca_EEG3-4 FC: freezing (a) 6.26%, (b) 6.00%

Context only, Fig_5c.csv wake co-bursting over all 24 recordings (%): {'Ca_EEG2-1': {'OxN': 10.37, 'OxA': 8.39}, 'Ca_EEG3-4': {'OxN': 14.35, 'OxA': 18.99}}

Result: reproduced. Both deposited mice appear in the source data under the subject ids of the files, with the same shock amplitude (1.5 mA for Ca_EEG2-1, 0.25 mA for Ca_EEG3-4). All four freezing values are within 10% of the paper's per-mouse values under both readings: reading (a) gives 64.09, 8.82, 17.40 and 15.06% against 63.87, 8.46, 16.90 and 14.60%, and reading (b) gives 63.45, 8.46, 16.91 and 14.61%. The session descriptions of Recall2 and Recall3 of Ca_EEG3-4 name the novel and neutral contexts, the assignment the source data use.

ed9b-aversive-freezing-high-vs-low

In the calcium imaging and EEG/EMG experiment, high-shock mice froze more during aversive recall than low-shock mice (t7 = 8.99, p = 0.020).

One mouse per group is deposited, so the comparison is between the two mice, with 30 s bins of the 5 min Recall1 session as the units of a within-session test.

fr = {m: ml.freezing_percent(data[(m, "Recall1")])[0] for m in SESSIONS}
bins = {m: ml.binned_freezing(data[(m, "Recall1")]) for m in SESSIONS}
for m in SESSIONS:
    print(f"{m} ({SOURCE_GROUP[m]} mA): Recall1 freezing {fr[m]:.2f}%; 30 s bins n={len(bins[m])}, "
          f"range {bins[m].min():.1f} to {bins[m].max():.1f}%, median {np.median(bins[m]):.1f}%")
u = stats.mannwhitneyu(bins["Ca_EEG2-1"], bins["Ca_EEG3-4"], alternative="two-sided")
print(f"Difference high minus low: {fr['Ca_EEG2-1'] - fr['Ca_EEG3-4']:.2f} percentage points; "
      f"Mann-Whitney on bins U = {u.statistic:.0f}, p = {u.pvalue:.2g}")
for w in (15.0, 60.0):
    b1, b2 = ml.binned_freezing(data[("Ca_EEG2-1", "Recall1")], w), ml.binned_freezing(data[("Ca_EEG3-4", "Recall1")], w)
    print(f"Variant, {w:.0f} s bins: n = {len(b1)} and {len(b2)}, Mann-Whitney p = {stats.mannwhitneyu(b1, b2).pvalue:.2g}")

Ca_EEG2-1 (1.5 mA): Recall1 freezing 64.09%; 30 s bins n=10, range 32.0 to 87.0%, median 68.7%
Ca_EEG3-4 (0.25 mA): Recall1 freezing 8.82%; 30 s bins n=10, range 2.2 to 13.4%, median 10.2%
Difference high minus low: 55.27 percentage points; Mann-Whitney on bins U = 100, p = 0.00018
Variant, 15 s bins: n = 20 and 20, Mann-Whitney p = 8.9e-07
Variant, 60 s bins: n = 5 and 5, Mann-Whitney p = 0.0079

Result: reproduced. The high-shock mouse froze during 64.09% of aversive recall and the low-shock mouse during 8.82%, a difference of 55.27 percentage points. The 30 s bins of the two sessions differ (Mann-Whitney p = 0.00018), as do 15 s and 60 s bins (p = 8.9e-07 and 0.0079). This compares one mouse per group; the paper's t-test over 5 and 4 mice cannot be run on the dandiset. The 30 s bins are consecutive parts of one session of each mouse, so they are not independent and the p-values overstate the evidence; the difference in freezing itself is large.

offline-bursts-above-shuffle

During offline periods, hippocampal population bursts occur more often than expected from shuffled neuronal activity, so CA1 neurons fire in a more coordinated way than chance.

The paper shows this for the one-hour offline recordings of Extended Data Figure 3; the dandiset holds three 10 min offline recordings of the EEG/EMG experiment, and the same test is run on each.

OFFLINE = [("Ca_EEG2-1", "OfflineDay2Session1"), ("Ca_EEG3-4", "OfflineDay1Session1"), ("Ca_EEG3-4", "OfflineDay2Session1")]
prep = {}
for k in OFFLINE:
    d = data[k]
    Z, keep = ml.zscore_cells(d["S"])
    pop = ml.population(Z)
    prep[k] = {"Z": Z, "keep": keep, "pop": pop, "minutes": (d["t"][-1] - d["t"][0]) / 60}
    ev = ml.burst_events(pop)
    lens = (ev[:, 1] - ev[:, 0]) / 15
    print(f"{' '.join(k)}: cells {Z.shape[1]} (left out {int((~keep).sum())}), {prep[k]['minutes']:.2f} min, "
          f"population z range {pop.min():.2f} to {pop.max():.2f}, bursts {len(ev)} "
          f"({len(ev) / prep[k]['minutes']:.2f}/min), duration {lens.min():.2f} to {lens.max():.2f} s (median {np.median(lens):.2f})")

rng_bursts = np.random.default_rng(7180)
rows = []
for k in OFFLINE:
    p = prep[k]
    for thr in (1.5, 2.0, 2.5):
        true = len(ml.burst_events(p["pop"], thr)) / p["minutes"]
        freq, _ = ml.shuffle_surrogates(p["Z"], ml.N_SHUFFLES, rng_bursts, p["minutes"], thr)
        rows.append({"recording": " ".join(k), "threshold_z": thr, "true_per_min": round(true, 2),
                     "shuffle_mean": round(freq.mean(), 2), "shuffle_p95": round(np.percentile(freq, 95), 2),
                     "shuffle_max": round(freq.max(), 2), "frac_shuffles_ge_true": (freq >= true).mean(),
                     "ratio_true_to_mean": round(true / freq.mean(), 2)})
burst_shuffle = pd.DataFrame(rows)
print(burst_shuffle.to_string(index=False))
main = burst_shuffle[burst_shuffle.threshold_z == 2.0]
print("\nAt z = 2: true above the shuffle 95th percentile in", int((main.true_per_min > main.shuffle_p95).sum()),
      "of 3; above the shuffle mean in", int((main.true_per_min > main.shuffle_mean).sum()), "of 3")

# The paper's chance level for this procedure (Extended Data Figure 3m, source data ExtData3m.csv,
# per-mouse range over 8 mice of another cohort, one-hour recordings at 30 Hz), against ours.
PAPER_ED3M = {"true_per_min": (4.51, 4.92), "shuffle_per_min": (3.38, 4.39), "shuffle_over_true": (0.72, 0.96)}
print("Paper (ED3m):", PAPER_ED3M)
print("Here, z = 2:", {"true_per_min": (main.true_per_min.min(), main.true_per_min.max()),
                       "shuffle_per_min": (main.shuffle_mean.min(), main.shuffle_mean.max()),
                       "shuffle_over_true": (round((main.shuffle_mean / main.true_per_min).min(), 2),
                                             round((main.shuffle_mean / main.true_per_min).max(), 2))})

Ca_EEG2-1 OfflineDay2Session1: cells 625 (left out 0), 9.99 min, population z range -1.06 to 6.39, bursts 87 (8.71/min), duration 0.07 to 1.27 s (median 0.27)
Ca_EEG3-4 OfflineDay1Session1: cells 711 (left out 0), 10.00 min, population z range -1.30 to 5.38, bursts 64 (6.40/min), duration 0.07 to 1.33 s (median 0.33)
Ca_EEG3-4 OfflineDay2Session1: cells 858 (left out 0), 10.00 min, population z range -1.23 to 7.88, bursts 53 (5.30/min), duration 0.07 to 1.80 s (median 0.40)

                    recording  threshold_z  true_per_min  shuffle_mean  shuffle_p95  shuffle_max  frac_shuffles_ge_true  ratio_true_to_mean
Ca_EEG2-1 OfflineDay2Session1          1.5         13.21         20.57        22.32        24.02                  1.000                0.64
Ca_EEG2-1 OfflineDay2Session1          2.0          8.71         10.94        12.11        13.01                  0.999                0.80
Ca_EEG2-1 OfflineDay2Session1          2.5          6.01          5.23         6.11         7.01                  0.085                1.15
Ca_EEG3-4 OfflineDay1Session1          1.5         11.20         17.69        19.41        20.81                  1.000                0.63
Ca_EEG3-4 OfflineDay1Session1          2.0          6.40          8.79         9.90        11.40                  1.000                0.73
Ca_EEG3-4 OfflineDay1Session1          2.5          4.30          3.87         4.70         5.70                  0.216                1.11
Ca_EEG3-4 OfflineDay2Session1          1.5         10.30         17.27        18.91        20.81                  1.000                0.60
Ca_EEG3-4 OfflineDay2Session1          2.0          5.30          8.48         9.61        10.70                  1.000                0.63
Ca_EEG3-4 OfflineDay2Session1          2.5          3.30          3.65         4.40         5.60                  0.795                0.90

At z = 2: true above the shuffle 95th percentile in 0 of 3; above the shuffle mean in 0 of 3
Paper (ED3m): {'true_per_min': (4.51, 4.92), 'shuffle_per_min': (3.38, 4.39), 'shuffle_over_true': (0.72, 0.96)}
Here, z = 2: {'true_per_min': (np.float64(5.3), np.float64(8.71)), 'shuffle_per_min': (np.float64(8.48), np.float64(10.94)), 'shuffle_over_true': (np.float64(1.26), np.float64(1.6))}

# A variant: count surrogate bursts against the true trace's own mean + 2 SD (no re-z-scoring).
rng_var = np.random.default_rng(7181)
for k in OFFLINE:
    p = prep[k]
    raw = p["Z"].mean(1)
    thr_abs = raw.mean() + 2 * raw.std()
    freq, _ = ml.shuffle_surrogates(p["Z"], ml.N_SHUFFLES, rng_var, p["minutes"], fixed=thr_abs)
    true = len(ml.burst_events(raw, thr_abs)) / p["minutes"]
    print(f"{' '.join(k)}: threshold fixed at the true trace's mean + 2 SD, true {true:.2f}/min, "
          f"surrogates mean {freq.mean():.2f}, 95th percentile {np.percentile(freq, 95):.2f}, max {freq.max():.2f}")
rng_merge = np.random.default_rng(7184)
for k in OFFLINE:
    p = prep[k]
    true = len(ml.burst_events(p["pop"], 2.0, merge=15)) / p["minutes"]
    freq, _ = ml.shuffle_surrogates(p["Z"], ml.N_SHUFFLES, rng_merge, p["minutes"], merge=15)
    print(f"{' '.join(k)}: events closer than 1 s merged, true {true:.2f}/min, surrogates mean {freq.mean():.2f}, "
          f"95th percentile {np.percentile(freq, 95):.2f}, fraction of surrogates >= true {(freq >= true).mean():.3f}")
fig, axes = plt.subplots(1, 3, figsize=(12, 2.6), sharey=True)
for ax, k in zip(axes, OFFLINE):
    t = data[k]["t"]
    ax.plot(t, prep[k]["pop"], lw=0.4, color="0.2")
    ax.axhline(2, color="tab:red", lw=0.8)
    ax.set_title(" ".join(k), fontsize=9)
    ax.set_xlabel("s")
axes[0].set_ylabel("population z")
plt.tight_layout()
plt.show()
print("Plotted population z ranges:", {" ".join(k): (round(prep[k]["pop"].min(), 2), round(prep[k]["pop"].max(), 2)) for k in OFFLINE},
      "; all points within the axes")

Ca_EEG2-1 OfflineDay2Session1: threshold fixed at the true trace's mean + 2 SD, true 8.71/min, surrogates mean 0.10, 95th percentile 0.30, max 0.50

Ca_EEG3-4 OfflineDay1Session1: threshold fixed at the true trace's mean + 2 SD, true 6.40/min, surrogates mean 0.01, 95th percentile 0.10, max 0.20

Ca_EEG3-4 OfflineDay2Session1: threshold fixed at the true trace's mean + 2 SD, true 5.30/min, surrogates mean 0.00, 95th percentile 0.00, max 0.10

Ca_EEG2-1 OfflineDay2Session1: events closer than 1 s merged, true 7.31/min, surrogates mean 8.52, 95th percentile 9.41, fraction of surrogates >= true 0.988

Ca_EEG3-4 OfflineDay1Session1: events closer than 1 s merged, true 6.00/min, surrogates mean 6.87, 95th percentile 7.70, fraction of surrogates >= true 0.962

Ca_EEG3-4 OfflineDay2Session1: events closer than 1 s merged, true 4.40/min, surrogates mean 6.56, 95th percentile 7.40, fraction of surrogates >= true 1.000

Plotted population z ranges: {'Ca_EEG2-1 OfflineDay2Session1': (np.float32(-1.06), np.float32(6.39)), 'Ca_EEG3-4 OfflineDay1Session1': (np.float32(-1.3), np.float32(5.38)), 'Ca_EEG3-4 OfflineDay2Session1': (np.float32(-1.23), np.float32(7.88))} ; all points within the axes

# Positive and negative control of the shuffle tests on synthetic data of the size of an offline
# recording (10 min at 15 Hz, 400 cells, each active independently in 1% of frames). "Coordinated"
# adds 6 events per minute in which a random 30% of cells are active together for 2 frames.
rng_ctrl = np.random.default_rng(7185)
T_s, N_s = 9000, 400
min_s = T_s / 15 / 60
rows = []
for coordinated in (True, False):
    S_s = (rng_ctrl.random((T_s, N_s)) < 0.01) * rng_ctrl.exponential(1, (T_s, N_s))
    if coordinated:
        for t0 in rng_ctrl.choice(T_s - 2, int(6 * min_s), replace=False):
            on = rng_ctrl.random(N_s) < 0.3
            S_s[t0:t0 + 2, on] += rng_ctrl.exponential(1, (2, on.sum()))
    Z_s, _ = ml.zscore_cells(S_s)
    pop_s = ml.population(Z_s)
    raw_s = Z_s.mean(1)
    thr_s = raw_s.mean() + 2 * raw_s.std()
    fq, sk = ml.shuffle_surrogates(Z_s, ml.N_SHUFFLES, rng_ctrl, min_s)
    fq_fix, _ = ml.shuffle_surrogates(Z_s, ml.N_SHUFFLES, rng_ctrl, min_s, fixed=thr_s)
    rows.append({"synthetic": "coordinated" if coordinated else "independent",
                 "criterion_true_per_min": round(len(ml.burst_events(pop_s)) / min_s, 2),
                 "criterion_shuffle_mean": round(fq.mean(), 2), "criterion_shuffle_p95": round(np.percentile(fq, 95), 2),
                 "fixed_true_per_min": round(len(ml.burst_events(raw_s, thr_s)) / min_s, 2),
                 "fixed_shuffle_mean": round(fq_fix.mean(), 2), "fixed_shuffle_p95": round(np.percentile(fq_fix, 95), 2),
                 "true_skew": round(stats.skew(pop_s), 2), "shuffle_skew_p95": round(np.percentile(sk, 95), 2)})
print(pd.DataFrame(rows).to_string(index=False))

  synthetic  criterion_true_per_min  criterion_shuffle_mean  criterion_shuffle_p95  fixed_true_per_min  fixed_shuffle_mean  fixed_shuffle_p95  true_skew  shuffle_skew_p95
coordinated                     5.6                   34.98                   36.9                 5.6                0.00                0.0       8.54              0.94
independent                    34.8                   37.36                   39.3                34.8               38.37               40.7       0.98              1.09

Result: blocked (methods underspecified). With surrogates z-scored and thresholded at z = 2 like the data, as the criterion fixed, the true burst rate is below the surrogate mean in all three recordings (5.30 to 8.71 against 8.48 to 10.94 per min; true over surrogate 0.63 to 0.80), which by the criterion would be not reproduced. The paper's chance level for this test, from its source data for Extended Data Figure 3m, is 3.38 to 4.39 surrogate bursts per min against 4.51 to 4.92 true (surrogate over true 0.72 to 0.96); this procedure gives 1.26 to 1.60. The procedures differ, and the Methods do not say how bursts are counted in a surrogate. The variants disagree with each other: merging events closer than 1 s leaves the true rate below the surrogates (4.40 to 7.31 against 6.56 to 8.52 per min), while thresholding surrogates at the true trace's own mean + 2 SD puts the true rate far above them (5.30 to 8.71 against at most 0.50 per min). At z = 2.5 the true rate is above the surrogate mean in two of three recordings and above the 95th percentile in none.

The controls on synthetic data show why. With 6 coordinated events per minute planted, the criterion's test (surrogates z-scored again and thresholded at z = 2) gives 5.6 true bursts per min against a surrogate mean of 34.98, and with independent cells 34.8 against 37.36. Re-z-scoring each surrogate makes it cross z = 2 as often as any trace with the same spread, so the test cannot detect coordination and its "not reproduced" says nothing about the claim. The fixed-threshold variant separates the two cases (5.6 against 0.00 per min with coordination, 34.8 against 38.37 without), but its chance level, surrogate means of at most 0.10 bursts per min against 5.30 to 8.71 true, is far below the paper's (surrogate over true 0.72 to 0.96). Neither matches the paper's procedure.

offline-skew-above-shuffle

The mean population activity during offline periods is more right-skewed than in shuffled data, showing brief burst-like activations that need coordinated activity of groups of neurons.

rng_skew = np.random.default_rng(7182)
rows = []
for k in OFFLINE:
    p = prep[k]
    true = stats.skew(p["pop"])
    _, sk = ml.shuffle_surrogates(p["Z"], ml.N_SHUFFLES, rng_skew, p["minutes"])
    rows.append({"recording": " ".join(k), "true_skew": round(true, 2), "shuffle_mean": round(sk.mean(), 2),
                 "shuffle_p95": round(np.percentile(sk, 95), 2), "shuffle_max": round(sk.max(), 2),
                 "frac_shuffles_ge_true": (sk >= true).mean(), "ratio": round(true / sk.mean(), 2)})
skew_shuffle = pd.DataFrame(rows)
print(skew_shuffle.to_string(index=False))
print("True skew above the shuffle 95th percentile in", int((skew_shuffle.true_skew > skew_shuffle.shuffle_p95).sum()), "of 3")
PAPER_ED3O = {"true": (1.76, 3.65), "shuffle": (0.47, 1.91)}   # ExtData3o.csv, per-mouse range, other cohort
print("Paper (ED3o):", PAPER_ED3O)

                    recording  true_skew  shuffle_mean  shuffle_p95  shuffle_max  frac_shuffles_ge_true  ratio
Ca_EEG2-1 OfflineDay2Session1       1.65          0.60         0.72         0.82                    0.0   2.74
Ca_EEG3-4 OfflineDay1Session1       1.34          0.40         0.51         0.63                    0.0   3.32
Ca_EEG3-4 OfflineDay2Session1       1.72          0.37         0.48         0.61                    0.0   4.69
True skew above the shuffle 95th percentile in 3 of 3
Paper (ED3o): {'true': (1.76, 3.65), 'shuffle': (0.47, 1.91)}

Result: reproduced. The skewness of the population trace is 1.65, 1.34 and 1.72 in the three offline recordings, above all 1000 surrogates in each (surrogate means 0.60, 0.40 and 0.37; 95th percentiles 0.72, 0.51 and 0.48; true over surrogate mean 2.74 to 4.69). Both arms are lower than in the paper's Extended Data Figure 3o source data, from another cohort recorded for 1 h at 30 Hz: the true skews here are all below the paper's lowest (1.76 to 3.65), and two of the three surrogate means are below the paper's lowest shuffled value (0.47 to 1.91). The procedure's chance level is therefore not shown to match the paper's; the outcome rests on the true skew exceeding every surrogate. In the synthetic controls the skew test separates coordinated activity (8.54 against a surrogate 95th percentile of 0.94) from independent cells (0.98 against 1.09), so this test, unlike the burst count, detects coordination.

burst-locomotion-pause

Offline population bursts coincide with the mouse briefly slowing down about 1 s before the burst and resuming its movement about 1 to 2 s later.

Movement is the magnitude of the telemetry activity signal in 100 ms bins (the paper tracked the mouse with a webcam, not deposited).

rng_move = np.random.default_rng(7183)
rows, traces = [], {}
lags = np.arange(-5, 5, ml.MOVE_BIN_S)
for k in OFFLINE:
    d, p = data[k], prep[k]
    mt, mv = ml.movement(d)
    ev = ml.burst_events(p["pop"])
    zero = {"a_peak": d["t"][ml.event_peaks(p["pop"], ev)], "b_onset": d["t"][ev[:, 0]],
            "random": np.sort(rng_move.uniform(d["t"][0] + 5, d["t"][-1] - 5, 300))}
    for reading, c in zero.items():
        pre = ml.window_mean(mt, mv, c, -1.0, 0.0)
        base = ml.window_mean(mt, mv, c, -5.0, -3.0)
        post = ml.window_mean(mt, mv, c, 1.0, 2.0)
        ok = ~np.isnan(pre) & ~np.isnan(base)
        w = stats.wilcoxon(pre[ok], base[ok], alternative="less")
        rows.append({"recording": " ".join(k), "time_zero": reading, "n": int(ok.sum()), "base_-5_-3": round(np.mean(base[ok]), 4),
                     "pre_-1_0": round(np.mean(pre[ok]), 4), "post_1_2": round(np.nanmean(post), 4),
                     "median_pre_minus_base": round(np.median(pre[ok] - base[ok]), 4),
                     "frac_lower": round(np.mean(pre[ok] < base[ok]), 2), "wilcoxon_p": w.pvalue})
        traces[(k, reading)] = np.array([np.nanmean(ml.window_mean(mt, mv, c, l, l + ml.MOVE_BIN_S)) for l in lags])
pause = pd.DataFrame(rows)
print(pause.to_string(index=False))
for reading in ("a_peak", "b_onset"):
    sub = pause[pause.time_zero == reading]
    print(f"Reading {reading}: p < 0.05 in {int((sub.wilcoxon_p < 0.05).sum())} of 3 recordings")

                    recording time_zero   n  base_-5_-3  pre_-1_0  post_1_2  median_pre_minus_base  frac_lower  wilcoxon_p
Ca_EEG2-1 OfflineDay2Session1    a_peak  85      0.2120    0.2247    0.1550                -0.0001        0.51    0.763133
Ca_EEG2-1 OfflineDay2Session1   b_onset  85      0.2094    0.2203    0.1670                 0.0004        0.48    0.884622
Ca_EEG2-1 OfflineDay2Session1    random 300      0.1091    0.1139    0.1052                -0.0004        0.63    0.015938
Ca_EEG3-4 OfflineDay1Session1    a_peak  63      0.1237    0.1278    0.1403                 0.0002        0.32    0.894867
Ca_EEG3-4 OfflineDay1Session1   b_onset  63      0.1300    0.1422    0.1494                 0.0003        0.30    0.939916
Ca_EEG3-4 OfflineDay1Session1    random 300      0.0921    0.1083    0.0896                 0.0002        0.27    0.999937
Ca_EEG3-4 OfflineDay2Session1    a_peak  53      0.2825    0.2574    0.1287                -0.0010        0.58    0.393577
Ca_EEG3-4 OfflineDay2Session1   b_onset  53      0.2777    0.2856    0.1294                -0.0011        0.60    0.414146
Ca_EEG3-4 OfflineDay2Session1    random 300      0.2257    0.2274    0.2274                 0.0026        0.47    0.919418
Reading a_peak: p < 0.05 in 0 of 3 recordings
Reading b_onset: p < 0.05 in 0 of 3 recordings

fig, axes = plt.subplots(1, 3, figsize=(12, 2.8), sharey=True)
for ax, k in zip(axes, OFFLINE):
    for reading, col in [("a_peak", "tab:blue"), ("b_onset", "tab:orange"), ("random", "0.6")]:
        ax.plot(lags + ml.MOVE_BIN_S / 2, traces[(k, reading)], color=col, label=reading)
    ax.axvline(0, color="k", lw=0.5)
    ax.set_title(" ".join(k), fontsize=9)
    ax.set_xlabel("s from time zero")
axes[0].set_ylabel("|activity| (telemetry, a.u.)")
axes[0].legend(fontsize=7)
plt.tight_layout()
plt.show()
print("Plotted ranges:", {f"{' '.join(k)} {r}": (round(np.nanmin(v), 4), round(np.nanmax(v), 4)) for (k, r), v in traces.items()})
print("Movement in 1 s bins against the population trace, Spearman r per recording:")
for k in OFFLINE:
    mt, mv = ml.movement(data[k])
    edges = np.arange(data[k]["t"][0], data[k]["t"][-1], 1.0)
    pm = ml.bin_matrix(prep[k]["pop"][:, None], data[k]["t"], edges)[:, 0]
    mm = ml.bin_matrix(mv[:, None], mt, edges)[:, 0]
    r = stats.spearmanr(pm, mm)
    print(f"  {' '.join(k)}: r = {r.statistic:.3f}, p = {r.pvalue:.2g}")

Plotted ranges: {'Ca_EEG2-1 OfflineDay2Session1 a_peak': (np.float64(0.101), np.float64(0.2713)), 'Ca_EEG2-1 OfflineDay2Session1 b_onset': (np.float64(0.0886), np.float64(0.2875)), 'Ca_EEG2-1 OfflineDay2Session1 random': (np.float64(0.0849), np.float64(0.1411)), 'Ca_EEG3-4 OfflineDay1Session1 a_peak': (np.float64(0.0925), np.float64(0.2199)), 'Ca_EEG3-4 OfflineDay1Session1 b_onset': (np.float64(0.1067), np.float64(0.2328)), 'Ca_EEG3-4 OfflineDay1Session1 random': (np.float64(0.0805), np.float64(0.1182)), 'Ca_EEG3-4 OfflineDay2Session1 a_peak': (np.float64(0.0936), np.float64(0.3676)), 'Ca_EEG3-4 OfflineDay2Session1 b_onset': (np.float64(0.0936), np.float64(0.348)), 'Ca_EEG3-4 OfflineDay2Session1 random': (np.float64(0.1961), np.float64(0.2554))}
Movement in 1 s bins against the population trace, Spearman r per recording:
  Ca_EEG2-1 OfflineDay2Session1: r = 0.086, p = 0.035
  Ca_EEG3-4 OfflineDay1Session1: r = 0.324, p = 3.8e-16
  Ca_EEG3-4 OfflineDay2Session1: r = 0.069, p = 0.09

Result: not reproduced. Telemetry movement in the second before a burst is not lower than 3 to 5 s before it in any recording, under either reading of time zero (peak: one-sided Wilcoxon p = 0.76, 0.89 and 0.39; onset: p = 0.88, 0.94 and 0.41), and the share of bursts with lower movement before them is 0.30 to 0.60. Movement 1 to 2 s after the burst peak is lower than in the second before it in two of three recordings (0.155 against 0.225, and 0.129 against 0.257), which is not a resumption of movement. The test is not calibrated: at 300 random times it gives p = 0.016 in Ca_EEG2-1 (0.9999 and 0.92 in the others), as overlapping windows of a slowly varying signal are not independent. The test rests on the telemetry ActivitySignal, which the files describe as motion of the probe relative to the receiver and a proxy for locomotion; it is higher in wake than in quiet wake (data check), but its relation to the calcium activity is weak (movement decoded at r = 0.06 to 0.24 at zero shift; Spearman r with the population trace 0.07 to 0.32 in 1 s bins). The paper used webcam tracking, which is not deposited, in another cohort. The result says that this proxy does not show the pause, and is weak evidence that the pause is absent.

fig5c-wake-overlap-neutral-coburst

After a strong (high-shock) aversive experience, during wake in the offline period the overlap ensemble co-bursts more with the neutral ensemble than with the aversive ensemble (t4 = 4.94, P = 0.016).

Ca_EEG2-1, first offline recording after aversive encoding (10 min of the 4 h the paper used).

k = ("Ca_EEG2-1", "OfflineDay2Session1")
d, p = data[k], prep[k]
lab = labels[p["keep"]]
ens_names = ["neutral", "aversive", "overlap", "remaining"]
ens_z = np.column_stack([ml.zscore(p["Z"][:, lab == e].mean(1)) for e in ens_names])
state = ml.state_of_frames(d["t"], d["sleep"])


def coburst_table(thr=2.0):
    ev = ml.burst_events(p["pop"], thr)
    pk = ml.event_peaks(p["pop"], ev)
    part = ml.participates(ens_z, ev, thr)
    N, A, O = part[:, 0], part[:, 1], part[:, 2]
    rows = []
    for wake_reading, states in [("a_wake+quiet", {"wake", "quiet wake"}), ("b_wake_only", {"wake"}), ("all_bursts", None),
                                 ("methods_emg_wake", "gmm")]:
        if states is None:
            sel = np.ones(len(ev), bool)
        elif states == "gmm":
            sel = gmm_wake_frames[pk]
        else:
            sel = np.isin(state[pk], list(states))
        for co_reading in ("i_inclusive", "ii_exclusive"):
            oxn = O & N if co_reading == "i_inclusive" else O & N & ~A
            oxa = O & A if co_reading == "i_inclusive" else O & A & ~N
            b, c = int((oxn & ~oxa & sel).sum()), int((oxa & ~oxn & sel).sum())
            mc = stats.binomtest(min(b, c), b + c, 0.5).pvalue if b + c else np.nan
            rows.append({"threshold": thr, "wake": wake_reading, "co": co_reading, "bursts": int(sel.sum()),
                         "OxN_pct": round(100 * oxn[sel].mean(), 2), "OxA_pct": round(100 * oxa[sel].mean(), 2),
                         "OxN_only": b, "OxA_only": c, "mcnemar_p": mc})
    return pd.DataFrame(rows), ev, pk, part


cob, ev, pk, part = coburst_table(2.0)
print("Cells per ensemble (active cells):", {e: int((lab == e).sum()) for e in ens_names})
print("Bursts:", len(ev), "; by state of the peak frame:", pd.Series(state[pk]).value_counts().to_dict())
print("Independent and joint participation of the ensembles, % of all bursts:",
      {e: round(100 * part[:, i].mean(), 1) for i, e in enumerate(ens_names)})
print(cob.to_string(index=False))
main = cob[cob.wake.isin(["a_wake+quiet", "b_wake_only"])]
print("\nThe criterion's four readings (a, b by i, ii); methods_emg_wake is the Methods' EMG scoring, a variant.")
print("Readings with OxN > OxA:", int((main.OxN_pct > main.OxA_pct).sum()), "of 4; of these with McNemar p < 0.05:",
      int(((main.OxN_pct > main.OxA_pct) & (main.mcnemar_p < 0.05)).sum()))

Cells per ensemble (active cells): {'neutral': 97, 'aversive': 80, 'overlap': 166, 'remaining': 282}
Bursts: 87 ; by state of the peak frame: {np.str_('quiet wake'): 62, np.str_('wake'): 25}
Independent and joint participation of the ensembles, % of all bursts: {'neutral': np.float64(42.5), 'aversive': np.float64(50.6), 'overlap': np.float64(63.2), 'remaining': np.float64(74.7)}
 threshold             wake           co  bursts  OxN_pct  OxA_pct  OxN_only  OxA_only  mcnemar_p
       2.0     a_wake+quiet  i_inclusive      87    26.44    35.63         7        15   0.133801
       2.0     a_wake+quiet ii_exclusive      87     8.05    17.24         7        15   0.133801
       2.0      b_wake_only  i_inclusive      25    32.00    28.00         2         1   1.000000
       2.0      b_wake_only ii_exclusive      25     8.00     4.00         2         1   1.000000
       2.0       all_bursts  i_inclusive      87    26.44    35.63         7        15   0.133801
       2.0       all_bursts ii_exclusive      87     8.05    17.24         7        15   0.133801
       2.0 methods_emg_wake  i_inclusive      87    26.44    35.63         7        15   0.133801
       2.0 methods_emg_wake ii_exclusive      87     8.05    17.24         7        15   0.133801

The criterion's four readings (a, b by i, ii); methods_emg_wake is the Methods' EMG scoring, a variant.
Readings with OxN > OxA: 2 of 4; of these with McNemar p < 0.05: 0

# Variants: burst threshold, and co-participation outside bursts (frames not in any burst, in 0.5 s windows).
print(pd.concat([coburst_table(1.5)[0], coburst_table(2.5)[0]]).to_string(index=False))
inburst = np.zeros(len(p["pop"]), bool)
for a, b in ev:
    inburst[a:b] = True
win = 7
rows = []
for wake_reading, states in [("a_wake+quiet", {"wake", "quiet wake"}), ("b_wake_only", {"wake"})]:
    starts = [i for i in range(0, len(p["pop"]) - win, win)
              if not inburst[i:i + win].any() and all(x in states for x in state[i:i + win])]
    pa = np.array([(ens_z[i:i + win] > 2).any(0) for i in starts])
    Nn, An, On = pa[:, 0], pa[:, 1], pa[:, 2]
    rows.append({"wake": wake_reading, "windows": len(starts), "OxN_pct": round(100 * (On & Nn).mean(), 2),
                 "OxA_pct": round(100 * (On & An).mean(), 2)})
print("\nOutside bursts (inclusive reading; windows with every frame in the state):")
print(pd.DataFrame(rows).to_string(index=False))
# Context: participation of each ensemble in bursts (fraction of its cells above z = 2 per burst), as in Figure 3d-f.
cellpart = ml.participates(p["Z"], ev)
frac = {e: cellpart[:, lab == e].mean(1) for e in ens_names}
print("\nMean fraction of each ensemble's cells participating per burst:", {e: round(float(v.mean()), 3) for e, v in frac.items()})
for e in ["neutral", "aversive", "overlap"]:
    w = stats.wilcoxon(frac[e], frac["remaining"])
    print(f"  {e} vs remaining across {len(ev)} bursts: median difference {np.median(frac[e] - frac['remaining']):.3f}, Wilcoxon p = {w.pvalue:.2g}")

 threshold             wake           co  bursts  OxN_pct  OxA_pct  OxN_only  OxA_only  mcnemar_p
       1.5     a_wake+quiet  i_inclusive     132    35.61    40.91        16        23   0.336784
       1.5     a_wake+quiet ii_exclusive     132    12.12    17.42        16        23   0.336784
       1.5      b_wake_only  i_inclusive      31    38.71    45.16         3         5   0.726562
       1.5      b_wake_only ii_exclusive      31     9.68    16.13         3         5   0.726562
       1.5       all_bursts  i_inclusive     132    35.61    40.91        16        23   0.336784
       1.5       all_bursts ii_exclusive     132    12.12    17.42        16        23   0.336784
       1.5 methods_emg_wake  i_inclusive     129    34.88    41.09        15        23   0.255875
       1.5 methods_emg_wake ii_exclusive     129    11.63    17.83        15        23   0.255875
       2.5     a_wake+quiet  i_inclusive      60    21.67    30.00         4         9   0.266846
       2.5     a_wake+quiet ii_exclusive      60     6.67    15.00         4         9   0.266846
       2.5      b_wake_only  i_inclusive      20    20.00    15.00         1         0   1.000000
       2.5      b_wake_only ii_exclusive      20     5.00     0.00         1         0   1.000000
       2.5       all_bursts  i_inclusive      60    21.67    30.00         4         9   0.266846
       2.5       all_bursts ii_exclusive      60     6.67    15.00         4         9   0.266846
       2.5 methods_emg_wake  i_inclusive      60    21.67    30.00         4         9   0.266846
       2.5 methods_emg_wake ii_exclusive      60     6.67    15.00         4         9   0.266846

Outside bursts (inclusive reading; windows with every frame in the state):
        wake  windows  OxN_pct  OxA_pct
a_wake+quiet     1131     0.44     0.88
 b_wake_only      196     0.51     1.02

Mean fraction of each ensemble's cells participating per burst: {'neutral': 0.053, 'aversive': 0.059, 'overlap': 0.057, 'remaining': 0.047}
  neutral vs remaining across 87 bursts: median difference 0.005, Wilcoxon p = 0.041
  aversive vs remaining across 87 bursts: median difference 0.009, Wilcoxon p = 0.0005
  overlap vs remaining across 87 bursts: median difference 0.010, Wilcoxon p = 0.0003

Result: not reproduced. In the one deposited high-shock recording, with wake read as wake and quiet wake together, overlap-neutral co-bursting (OxN) is below overlap-aversive co-bursting (OxA): 26.4 against 35.6% of all 87 bursts with both ensembles counted whatever the third does, 8.0 against 17.2% with the third excluded; 7 against 15 discordant bursts, McNemar p = 0.13). This is the reading the Methods imply: their sleep scoring, applied to this mouse's EMG in the data check, puts 94.1% of quiet-wake frames in wake, and all 87 bursts fall in its wake (the methods_emg_wake rows give the same numbers). With wake read as the wake label alone, the other reading the criterion fixed, OxN exceeds OxA (32.0 against 28.0% of 25 bursts, and 8.0 against 4.0%), on 2 against 1 discordant bursts (McNemar p = 1.0); that reading picks out epochs the Methods' scoring does not separate, so the outcome is judged on the first. The order under the Methods' reading holds at thresholds 1.5 (34.9 against 41.1%) and 2.5 (21.7 against 30.0%), and no difference reaches p < 0.05. The evidence is weak: this is the first of the 24 offline recordings of one of the five high-shock mice; the ensemble traces are z-scored within this 10 min recording, while the Methods z-score them over all 24 concatenated, so the levels are not comparable with the paper's values for this mouse over all 24 recordings (10.37% OxN and 8.39% OxA in wake). With both ensembles counted whatever the third does, co-participation here (26 to 36%) is 2.5 to 4 times those values; with the third excluded (8.0 and 17.2% under the Methods' reading) it is on their scale, but OxN is the lower; and the paper's test across mice cannot be run. Outside bursts both co-participations are at most 1.02% of 0.5 s windows. The neutral, aversive and overlap ensembles each have a larger fraction of their cells in bursts than the remaining cells (median differences 0.005 to 0.010, Wilcoxon p = 0.041 to 0.0003).

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
source-data-link: the deposited mice match the paper's per-mouse freezing (outside data)	reproduced
ed9b-aversive-freezing-high-vs-low: high-shock mouse freezes more in aversive recall	reproduced (one mouse per group)
offline-bursts-above-shuffle: bursts exceed shuffled expectation	blocked: methods underspecified (chance level differs from the paper's)
offline-skew-above-shuffle: population activity more skewed than shuffled	reproduced
burst-locomotion-pause: movement pauses around bursts	not reproduced (telemetry proxy, weak evidence)
fig5c-wake-overlap-neutral-coburst: overlap-neutral exceeds overlap-aversive co-bursting in wake, high shock	not reproduced (one recording of one mouse; changed from partially reproduced after the first verification)
Behaviour-only experiments, three claims	blocked: data missing
One-hour imaging and single-experience experiments, eleven claims	blocked: data missing
Inhibitory tagging and no-shock experiments, two claims	blocked: data missing
EEG cohort claims needing missing mice, sessions or registration, four claims	blocked: data missing

How the analysis differs from the paper's: The paper's analysis code could not be read (its repository returns 404), so bursts, participation and co-participation follow the Methods, with burst events taken as maximal runs above threshold. The EEG cohort claims are tested on one or two mice of nine and on the first offline recording of 24. The burst and skew tests are run on the EEG cohort's 10 min recordings at 15 Hz, while the paper shows them for another cohort with one-hour recordings at 30 Hz. Movement around bursts uses the telemetry activity signal in place of webcam tracking. Quiet wake is not a state of the paper; the criterion fixed two mappings, and the Methods' EMG scoring, applied to Ca_EEG2-1, puts quiet wake in wake, so fig5c is judged on that mapping. In fig5c the ensemble traces are z-scored within one 10 min recording, not over 24 concatenated recordings as in the Methods. The burst shuffle test of the criterion cannot detect coordinated firing (synthetic controls), and no variant matches the paper's chance level.

What was found about the dandiset:

• It holds two of the nine mice of the EEG/EMG experiment and, of each offline period, only the first 10 min recording of 24. Ca_EEG2-1 lacks its offline-day-1 recordings and its neutral and novel recall sessions. The mice of the one-hour offline, behaviour-only, inhibitory tagging and no-shock experiments are absent, so most of the paper's findings cannot be tested.
• Ca_EEG3-4 has no cell registration, so its ensembles cannot be formed. The Week file of Ca_EEG2-1 holds 19 registration tables, for recordings of which only one is deposited.
• The session start times of the Ca_EEG2-1 behaviour sessions, and the Week epochs built from them, are about 1 h later than the telemetry and Miniscope clock; as deposited, the offline recording after aversive encoding starts 47 min before the FC session.
• The OnePhotonSeries timestamps of three files run at 60 Hz and cover half the recording.
• The sleep states include quiet wake, which the paper does not use or map to its wake, NREM and REM; the three offline recordings hold no NREM or REM.
• The telemetry is sampled at 500 Hz; the Methods say 100 Hz.
• The offline session description says the recordings lasted an hour with a webcam; the deposited offline recordings last 10 min and hold telemetry, not video.
• The freezing intervals reproduce the paper's per-mouse freezing for both mice to within 5% (1% with freezing as summed interval time), and the registration tables agree with the ROI centroids.

print(f"Notebook ran in {(time.time() - started) / 60:.1f} min")

Notebook ran in 1.2 min