Reproducing Gerken, Darcher et al. (2025) from dandiset 001616¶
Paper: Gerken F, Darcher A, Gonçalves PJ, Rapp R, Elezi I, Niediek J, Kehl MS, Reber TP, Liebe S, Macke JH, Mormann F, Leal-Taixé L (2025). Decoding movie content from neuronal population activity in the human medial temporal lobe. eLife reviewed preprint 106758, version 1. doi:10.7554/eLife.106758.1. This attempt used the full text, which is open access, and the supplementary file with Tables S1 and S2.
Dandiset: DANDI:001616, version 0.260702.0824 (SUMMER). The paper cites no DANDI version; its data statement says the data will be published separately. This is the most recent published version. The earlier published version, 0.260519.1133, holds the same 29 files (same blobs) and 28 older copies of the same recordings with identical spike times (shown in the data check), so the analysis would read the same data from either; the pinned version holds each recording once.
Code read: the dataset's repository mormannlab/SUMMER at
commit 59fc242 (Zenodo 10.5281/zenodo.21165561 archives it). It holds the NWB conversion, the
responsive-unit and cluster code behind Figure 2, plotting notebooks that hold the paper's LSTM
results as numbers, and a simple linear training baseline. It does not hold the paper's LSTM or its
hyperparameters. The fold layout, the grouping of hippocampal subregions and the onset rule here
follow it, as stated in summer.py.
No data from outside the dandiset is used. Each claim's section has its own random generator, and
the logistic regression is seeded. The notebook takes about 22 minutes, most of it in the logistic
regression fits. Three companion notebooks use the same helper module: top_ranked_units.ipynb
(fig7a and fig7b), and responsive_zeta_cuts.ipynb and responsive_zeta_labels.ipynb, which run
the ZETA test of the dataset's code for every unit as a variant of the responsive-unit claims.
| Claim | Status |
|---|---|
| dataset-composition | tested |
| fig2b-phc-camera-cuts | tested; blocked, the responsive test is not calibrated |
| tables1-camera-cuts-regions | tested; blocked, depends on the responsive test |
| fig2b-phc-scene-cuts | tested; blocked, the responsive test is not calibrated |
| tables1-few-semantic-responses | tested; blocked, the responsive test is not calibrated |
| tables1-no-indoor-response | tested; blocked, the responsive test is not calibrated |
| results-cluster-regions | tested; blocked, depends on the responsive test |
| fig3b-single-neuron-threshold | tested |
| fig3c-logreg-characters-location | tested |
| fig3c-logreg-cuts-lower | tested |
| fig3e-temporal-gap | tested, with the logistic regression in place of the LSTM |
| fig4-patients-below-pooled | tested |
| fig4-named-patients | tested, assuming subject ids are the paper's patient numbers |
| fig5-region-decoding | tested, with the logistic regression in place of the LSTM |
| fig5-faces-presence | tested, with the logistic regression in place of the LSTM |
| fig6-responsive-subsets | tested, with the logistic regression; blocked, depends on the responsive test |
| fig7b-top500-overlap | tested, in top_ranked_units.ipynb |
| fig7a-500-neurons | tested, with the logistic regression in place of the LSTM, in top_ranked_units.ipynb |
| fig3c-lstm-decoding | blocked: the LSTM's code and hyperparameters are not published |
| fig3d-temporal-shuffle | blocked: as above |
import os
for var in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
os.environ.setdefault(var, "1")
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import summer
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)
started_all = started = time.time()
sessions = summer.load_all()
units = summer.unit_table(sessions)
spikes = summer.movie_spikes(sessions)
print(f"loaded {len(sessions)} files, {len(units)} units, from dandiset {summer.DANDISET_ID} "
f"version {summer.DANDISET_VERSION} in {time.time() - started:.0f} s")
loaded 29 files, 2286 units, from dandiset 001616 version 0.260702.0824 in 117 s
Data check¶
Each check of the metadata against the values is stated under its output.
# Versions: what the pinned version holds against the earlier published version.
pinned = summer.asset_urls(summer.DANDISET_VERSION)
earlier = summer.asset_urls("0.260519.1133")
blob = {p: b for p, b, _ in pinned}
same = [p for p, b, _ in earlier if blob.get(p) == b]
extra = [(p, u) for p, b, u in earlier if p not in blob]
print(f"version {summer.DANDISET_VERSION}: {len(pinned)} NWB files; version 0.260519.1133: {len(earlier)} NWB files")
print(f"files of 0.260519.1133 also in the pinned version, same path and blob: {len(same)} of {len(pinned)}")
print(f"files only in 0.260519.1133: {len(extra)}, named like {extra[0][0]}")
from concurrent.futures import ThreadPoolExecutor
with ThreadPoolExecutor(8) as pool:
old = list(pool.map(lambda a: summer.unit_summary(a[1]), extra))
rows = []
for o in old:
s = next(s for s in sessions if s["subject_id"] == o["subject_id"])
new_spikes = np.sort(np.concatenate(s["spikes_ms"]))
same_n = len(new_spikes) == len(o["spikes_ms"])
rows.append({"patient": int(o["subject_id"]), "units_old": o["n_units"], "units_pinned": len(s["units"]),
"spikes_old": len(o["spikes_ms"]), "spikes_pinned": len(new_spikes),
"same_spike_times": bool(same_n and np.allclose(new_spikes, o["spikes_ms"], atol=1e-3)),
"same_regions": sorted(o["regions"]) == sorted(s["units"].brain_region)})
old_cmp = pd.DataFrame(rows).sort_values("patient")
print(old_cmp.to_string(index=False))
print("patients of the pinned version without an older copy:",
sorted(set(int(s["subject_id"]) for s in sessions) - set(old_cmp.patient)))
print("older copies with the same spike times as the pinned file:", int(old_cmp.same_spike_times.sum()), "of", len(old_cmp))
version 0.260702.0824: 29 NWB files; version 0.260519.1133: 57 NWB files files of 0.260519.1133 also in the pinned version, same path and blob: 29 of 29 files only in 0.260519.1133: 28, named like sub-1/sub-1_ses-1_behavior+ecephys.nwb
patient units_old units_pinned spikes_old spikes_pinned same_spike_times same_regions
1 40 40 397655 397655 True True
2 83 83 535482 535482 True True
4 120 120 2087660 2087660 True True
5 79 79 948767 948767 True True
6 99 99 1271620 1271620 True True
7 84 84 1577631 1577631 True True
8 45 45 261971 261971 True True
9 41 41 301121 301121 True True
10 108 108 1809222 1809222 True True
11 31 31 362082 362082 True True
12 52 52 409441 409441 True True
13 137 137 1451503 1451503 True True
14 30 30 148301 148301 True True
15 97 97 674127 674127 True True
16 67 67 528123 528123 True True
17 101 101 967681 967681 True True
18 40 40 245152 245152 True True
19 105 105 729081 729081 True True
20 88 88 1159216 1159216 True True
21 109 109 1291303 1291303 True True
22 79 79 1258690 1258690 True True
23 50 50 303595 303595 True True
24 88 88 403093 403093 True True
25 110 110 515371 515371 True True
26 85 85 420460 420460 True True
27 90 90 910449 910449 True True
28 111 111 1559006 1559006 True True
29 66 66 229656 229656 True True
patients of the pinned version without an older copy: [3]
older copies with the same spike times as the pinned file: 28 of 28
The earlier published version holds the 29 files of the pinned version unchanged and 28 more,
named ses-1_behavior+ecephys, one per patient except patient 3. Each of these holds the same
units, regions and spike times as the pinned file of its patient, so that version deposits 28
recordings twice. The pinned version holds each recording once.
pat = pd.DataFrame([{"patient": int(s["subject_id"]), "sex": s["sex"], "age": s["age"], "units": len(s["units"]),
"single_units": int(s["units"].is_single_unit.sum())} for s in sessions])
ages = pat.age.str.extract(r"P(\d+)Y")[0].astype(int)
print(f"patients {len(pat)}; female {int((pat.sex == 'F').sum())}, male {int((pat.sex == 'M').sum())}; "
f"ages {ages.min()} to {ages.max()} years (median {ages.median():.0f})")
print(f"units {pat.units.sum()} ({pat.single_units.sum()} marked single units); per patient {pat.units.min()} to {pat.units.max()}")
print("\nUnits by `brain_region` as stored:")
print(units.brain_region.value_counts().to_string())
print("\nUnits by region as grouped here (H = AH, MH, PH; FF, LG and PRC are 'other'):")
print(units.region.value_counts().to_string())
print("\nUnits per patient and region:")
print(pd.crosstab(units.patient, units.region, margins=True).to_string())
mismatch = 0
for s in sessions:
region_of = dict(zip(s["electrodes"].csc_nr, s["electrodes"].brain_region))
mismatch += int((s["units"].csc_nr.map(region_of) != s["units"].brain_region).sum())
print("\nunits whose region differs from that of their electrode:", mismatch)
print("units table description:", sessions[0]["units_description"])
patients 29; female 17, male 12; ages 19 to 62 years (median 40) units 2286 (1190 marked single units); per patient 30 to 137 Units by `brain_region` as stored: brain_region A 559 EC 493 AH 418 PH 412 PHC 293 PIC 68 MH 16 LG 15 FF 6 PRC 6 Units by region as grouped here (H = AH, MH, PH; FF, LG and PRC are 'other'): region H 846 A 559 EC 493 PHC 293 PIC 68 other 27 Units per patient and region: region A EC H PHC PIC other All patient 1 8 16 10 6 0 0 40 2 36 12 25 10 0 0 83 3 19 10 18 4 0 0 51 4 27 17 66 10 0 0 120 5 19 37 16 7 0 0 79 6 7 52 20 20 0 0 99 7 20 30 33 1 0 0 84 8 19 12 14 0 0 0 45 9 19 11 0 11 0 0 41 10 19 15 59 15 0 0 108 11 2 0 25 4 0 0 31 12 0 29 19 4 0 0 52 13 34 26 52 19 0 6 137 14 5 4 10 11 0 0 30 15 35 23 20 19 0 0 97 16 14 17 36 0 0 0 67 17 16 32 41 12 0 0 101 18 25 0 15 0 0 0 40 19 19 29 36 6 0 15 105 20 15 10 47 16 0 0 88 21 24 19 51 15 0 0 109 22 30 4 43 2 0 0 79 23 7 2 18 9 14 0 50 24 33 12 35 0 2 6 88 25 18 0 53 24 15 0 110 26 32 22 13 14 4 0 85 27 18 18 26 27 1 0 90 28 18 16 25 27 25 0 111 29 21 18 20 0 7 0 66 All 559 493 846 293 68 27 2286 units whose region differs from that of their electrode: 0 units table description: Spike times (milliseconds), waveforms (uV), and associated spike sorting metrics for 1's movie session. Spike times are reference to the start of the movie, such that a spike time == 0 would indicate a spike at the exact onset of the movie.
29 patients, 17 of them female, and 2286 units, 30 to 137 per patient, as the paper states. The
paper's Results give ages 22 to 63 for these 29 patients; the files give 19 to 62, the range the
Methods give for all 46 patients recorded. The region counts differ from the paper's (A 580, H 794,
EC 440, PHC 373, PIC 99): grouping AH, MH and PH as hippocampus, as the dataset's code does, gives
A 559, H 846, EC 493, PHC 293 and PIC 68, and 27 units lie in regions the paper does not name (FF,
LG, PRC). The totals agree, so the units were assigned to regions differently. Every unit's region
matches that of its electrode. No unit is labelled H; the dataset's training code
(ML_framework/nwb_loading/nwb_loading.py) compares brain_region with H when asked for the
hippocampus, so it would select no unit. The units table states its spike times in milliseconds,
where NWB expects seconds.
rows = []
for s in sessions:
w = s["cleaned_watchlogs"]
raw = s["raw_watchlogs"]
pa = summer.pauses(w)
sp = np.concatenate(s["spikes_ms"])
in_pause = sum(int(((sp > a) & (sp < b)).sum()) for a, b in pa)
rows.append({"patient": int(s["subject_id"]), "frames_logged": len(w),
"analysis_frames_covered": float(np.isin(np.round(s["frame_times_analysis"], 2), np.round(w.pts, 2)).mean()),
"pts_first_s": w.pts.iloc[0], "pts_last_s": w.pts.iloc[-1],
"pauses": len(pa), "pause_min": (pa[:, 1] - pa[:, 0]).sum() / 60000 if len(pa) else 0.0,
"frame_interval_ms": float(np.median(np.diff(w.rec_ms[np.r_[True, np.diff(w.pts) != 0]]))),
"raw_backward_steps": int((np.diff(raw.pts) < 0).sum()),
"spike_min_ms": sp.min(), "spike_max_ms": sp.max(), "rec_last_ms": w.rec_ms.iloc[-1],
"spikes_in_pauses": in_pause, "spikes_outside_log": int(((sp < w.rec_ms.iloc[0]) | (sp > w.rec_ms.iloc[-1])).sum())})
tb = pd.DataFrame(rows)
print(tb.round(3).to_string(index=False))
print("\nmedian frame interval on the recording clock, ms:", tb.frame_interval_ms.min().round(2), "to", tb.frame_interval_ms.max().round(2))
print("analysis frames (credits excluded):", len(sessions[0]["frame_times_analysis"]), "from",
sessions[0]["frame_times_analysis"][0], "to", sessions[0]["frame_times_analysis"][-1], "s; all frames:",
len(sessions[0]["frame_times_base"]))
print("80 ms bin edges identical in all files:", all(np.array_equal(s["edges_s"], sessions[0]["edges_s"]) for s in sessions),
"; bins:", len(sessions[0]["edges_s"]) - 1)
print("annotation and indicator tables identical in all files:",
all(s["annotations_base"].equals(sessions[0]["annotations_base"]) for s in sessions),
all(all(np.array_equal(s["indicator"][k], sessions[0]["indicator"][k]) for k in sessions[0]["indicator"]) for s in sessions))
# Do the patient-aligned annotations (stated in seconds in the column, ms in the table description)
# follow from the movie-time annotations through the watchlog?
dev = []
for s in sessions:
b = s["annotations_base"]
a = s["annotations_patient_aligned"]
for lab in ("summer", "camera-cuts"):
bs = b[(b.label == lab) & (b.value == 1)].start.to_numpy()
as_ = a[(a.label == lab) & (a.value == 1)].start.to_numpy()
w = s["cleaned_watchlogs"]
inside = (bs >= w.pts.iloc[0]) & (bs <= w.pts.iloc[-1])
mapped = summer.movie_to_rec(bs[inside], w)
nearest = np.abs(as_[None, :] - mapped[:, None]).min(1) if len(as_) else np.array([np.nan])
dev.append({"patient": int(s["subject_id"]), "label": lab, "base_onsets": int(inside.sum()),
"aligned_onsets": len(as_), "median_abs_diff_ms": float(np.median(nearest))})
dev = pd.DataFrame(dev)
print("\nmovie-time onsets mapped through the watchlog against `annotations_patient_aligned` (ms):")
print(dev.groupby("label")[["base_onsets", "aligned_onsets", "median_abs_diff_ms"]].describe().round(1).T.to_string())
print("column description:", repr(sessions[0]["annotations_patient_aligned_unit"]))
patient frames_logged analysis_frames_covered pts_first_s pts_last_s pauses pause_min frame_interval_ms raw_backward_steps spike_min_ms spike_max_ms rec_last_ms spikes_in_pauses spikes_outside_log
1 120102 1.0 0.04 4803.84 6 57.339 40.998 39 4.463 8244127.247 8244130.789 157810 0
2 120702 1.0 17.08 4845.00 3 14.487 41.037 29 8.243 5714201.575 5714204.221 81354 2086
3 125746 1.0 0.04 5029.68 5 20.367 41.026 11 -147824.206 6251529.734 6251542.292 100166 2
4 125681 1.0 2.48 5029.68 0 0.000 41.025 140 2.771 5180257.677 5180260.387 0 60519
5 119511 1.0 3.68 4784.08 0 0.000 41.015 0 0.024 6664708.111 4784060.669 0 282760
6 125690 1.0 2.24 5029.64 4 3.277 41.022 1 0.181 5609625.519 5226128.630 37895 84469
7 123441 1.0 19.40 4956.96 1 0.766 41.019 0 23.279 5106875.850 5084884.551 14495 36360
8 120633 1.0 0.04 4825.24 2 0.391 41.006 0 -1135450.700 4848548.072 4848595.609 1418 1
9 125744 1.0 0.04 5029.68 3 4.579 41.019 11 16.254 5304368.213 5304378.473 15939 0
10 122996 1.0 0.04 4919.44 11 17.222 41.018 0 3.036 5952549.468 5952552.932 298658 0
11 122820 1.0 18.16 4930.88 1 0.007 41.010 0 5.393 4997858.356 4997893.238 26 6485
12 121184 1.0 0.04 4847.12 6 2.464 41.013 17 3.277 5015299.911 5015305.729 12292 1363
13 121178 1.0 0.04 4847.08 2 0.679 41.019 0 1.516 4887760.941 4887761.057 11541 0
14 120107 1.0 0.04 4804.24 1 0.224 41.016 2 57.136 4957026.635 4957071.850 373 3696
15 119518 1.0 0.04 4780.60 3 0.411 40.690 0 14.561 4805149.968 4805154.321 3275 0
16 119679 1.0 0.04 4787.04 3 1.544 40.687 0 5.937 4879581.675 4879581.957 10194 0
17 125743 1.0 0.04 5029.68 1 0.448 40.681 0 1.130 5056355.760 5056356.625 5500 0
18 125746 1.0 0.04 5029.68 4 3.713 40.685 0 9.273 5252356.975 5252402.529 10915 0
19 125747 1.0 0.04 5029.68 5 5.287 40.682 0 2.975 5346769.925 5346770.154 41103 0
20 125743 1.0 0.04 5029.68 1 0.764 40.682 0 0.036 5075466.610 5075469.125 11686 0
21 119841 1.0 6.64 4800.24 0 0.000 40.687 29 6.391 5060253.234 4888067.860 0 63303
22 118881 1.0 23.80 4778.88 3 14.520 40.685 1 3.696 5879539.455 5737284.380 178808 49642
23 125663 1.0 3.12 5029.60 0 0.000 40.687 0 8.495 5217031.723 5216948.401 0 10397
24 119925 1.0 0.40 4797.36 0 0.000 41.111 1 12.227 4838102.447 4838105.711 0 3984
25 125199 1.0 0.04 5007.92 3 0.922 41.787 0 4.218 5063058.618 5063068.965 5986 0
26 123858 1.0 0.04 4954.28 1 1.279 41.162 0 7.233 5030976.937 5030993.064 6426 0
27 121717 1.0 0.04 4868.68 1 1.036 41.191 0 0.997 5092099.839 4930779.294 10541 29711
28 120694 1.0 0.04 4827.76 1 0.815 41.094 0 0.482 4878111.985 4878119.780 16115 0
29 119994 1.0 0.04 4799.64 3 0.764 41.103 0 17.722 4845384.709 4845388.835 3206 0
median frame interval on the recording clock, ms: 40.68 to 41.79
analysis frames (credits excluded): 118162 from 36.72 to 4763.16 s; all frames: 125743
80 ms bin edges identical in all files: True ; bins: 59080
annotation and indicator tables identical in all files: True True
movie-time onsets mapped through the watchlog against `annotations_patient_aligned` (ms):
label camera-cuts summer
base_onsets count 29.0 29.0
mean 957.3 141.0
std 1.0 0.0
min 955.0 141.0
25% 956.0 141.0
50% 958.0 141.0
75% 958.0 141.0
max 958.0 141.0
aligned_onsets count 29.0 29.0
mean 975.4 146.6
std 50.8 13.6
min 949.0 141.0
25% 958.0 141.0
50% 958.0 141.0
75% 968.0 142.0
max 1181.0 194.0
median_abs_diff_ms count 29.0 29.0
mean 41.0 41.0
std 0.2 0.2
min 40.7 40.7
25% 40.7 40.7
50% 41.0 41.0
75% 41.0 41.0
max 41.8 41.8
column description: 'Start time of epoch, in seconds'
Every frame of the analysis portion (36.72 to 4763.16 s, credits excluded) is in every patient's cleaned watchlog. Playback was paused in 24 patients, for up to 57 minutes in all (patient 1); spikes in pauses are dropped (up to 298658 in patient 10). Frames advance every 40.68 to 41.79 ms of the recording clock while the film has 40 ms frames, so the recording clock is in milliseconds and playback ran 2 to 4% slower than 25 frames per second; spikes are mapped to movie time through the watchlog, which absorbs this. Spikes span the logged playback in every file, and 15 files also hold spikes outside it. The patient-aligned annotations agree with the movie-time annotations mapped through the watchlog to a median of one frame (40.7 to 41.8 ms). Their start and stop columns are described as seconds but hold milliseconds, as the table description says. In some patients that table holds more onsets than the film has (up to 1181 camera cuts against 958), apparently from replayed portions; the analysis uses the movie-time annotations. All files carry the same annotations, indicator functions and 80 ms bin edges (59080 bins).
s0 = sessions[0]
fa = s0["frame_times_analysis"]
analysis = (s0["frame_times_base"] >= fa[0]) & (s0["frame_times_base"] <= fa[-1])
rows = []
for name, (label, value) in summer.LABELS.items():
v = s0["indicator"][label]
rows.append({"label": name, "dandiset_label": label, "values": str(sorted(np.unique(v).tolist())),
"fraction_of_analysis_frames": float(np.mean(v[analysis] == value)),
"segments": int(((s0["annotations_base"].label == label) & (s0["annotations_base"].value == value)).sum())})
lab = pd.DataFrame(rows)
print(lab.round(3).to_string(index=False))
print("Paper (Results): Summer 50/50, Tom 80/20, Persons 95/5, McKenzie 10/90 of frames.")
cc = s0["indicator"]["camera-cuts"]
sc = s0["indicator"]["scenes"]
movie_s = s0["frame_times_base"][-1]
print(f"\ncamera cuts: {int(cc.sum())} single-frame events, one every {movie_s / cc.sum():.2f} s (supplement: one every 3 s)")
print(f"scene cuts: {int(sc.sum())}, mean scene length {movie_s / sc.sum():.1f} s (supplement: 36.15 s); "
f"scene cut frames that are also camera cut frames: {int(((sc == 1) & (cc == 1)).sum())}")
print("frames of indoor-setting by value (1 indoor, 0 outdoor, 99 neither):",
pd.Series(s0["indicator"]["indoor-setting"][analysis]).value_counts().to_dict())
print("face labels hold only 0 and 1 (the Methods also exclude frames):",
[sorted(np.unique(s0["indicator"][k]).tolist()) for k in ("summer-faces", "tom-faces")])
print("all 53 labels:", ", ".join(sorted(s0["indicator"])))
label dandiset_label values fraction_of_analysis_frames segments
Summer summer [0, 1] 0.527 141
Tom tom [0, 1] 0.792 201
McKenzie mckenzie [0, 1] 0.090 55
Persons persons [0, 1] 0.946 89
Summer Faces summer-faces [0, 1] 0.298 978
Tom Faces tom-faces [0, 1] 0.448 2077
Summer Presence summer-presence [0, 1] 0.643 41
Camera Cuts camera-cuts [0, 1] 0.008 958
Scene Cuts scenes [0, 1] 0.001 133
Indoor indoor-setting [0, 1, 99] 0.672 44
Outdoor indoor-setting [0, 1, 99] 0.321 50
Paper (Results): Summer 50/50, Tom 80/20, Persons 95/5, McKenzie 10/90 of frames.
camera cuts: 958 single-frame events, one every 5.25 s (supplement: one every 3 s)
scene cuts: 133, mean scene length 37.8 s (supplement: 36.15 s); scene cut frames that are also camera cut frames: 121
frames of indoor-setting by value (1 indoor, 0 outdoor, 99 neither): {1: 79390, 0: 37933, 99: 839}
face labels hold only 0 and 1 (the Methods also exclude frames): [[0, 1], [0, 1]]
all 53 labels: alison, autumn, beach, bus, cafe, camera-cuts, car, days-of-summer, douche, elevator, family-home, gallery, ikea, indoor-setting, karaoke-bar, mckenzie, mckenzie-faces, millie, millie-faces, office, other-cafe, park, paul, paul-faces, persons, punch-bar, rachel, rachel-faces, record-store, restaurant, rhoda, scenes, secretary, soccer-field, street, summer, summer-apartment, summer-body-sequence, summer-child-bedroom, summer-faces, summer-presence, summer-speaking, the-graduate, theater, tom, tom-apartment, tom-child-bedroom, tom-faces, tom-speaking, train, vance, vance-faces, wedding-venue
The fractions of frames match the paper's descriptions (Summer 0.53, Tom 0.79, Persons 0.95,
McKenzie 0.09 of the analysis frames). The dandiset has 958 camera cuts, one every 5.25 s, where
the supplement says one every 3 s. Its 133 scene cuts give a mean scene of 37.8 s (supplement
36.15 s), and 121 of them fall on a camera cut frame. Indoor/Outdoor is indoor-setting: 1 indoor,
0 outdoor, and 99 for 839 analysis frames that are neither. The face labels hold only 0 and 1,
while the Methods define face labels with a third state for excluded frames; those frames cannot be
identified from the dandiset.
# Duplicates within the pinned version: identical spike trains across units, identical watchlogs.
keys = {}
dupes = []
for i, t in enumerate(spikes):
if len(t) < 20:
continue
k = tuple(np.round(t[:20], 4))
if k in keys:
dupes.append((keys[k], i))
keys.setdefault(k, i)
print("pairs of units with the same first 20 spike times:", len(dupes))
for a, b in dupes:
print(" ", units.iloc[a][["patient", "unit_id", "brain_region", "n_spikes"]].to_dict(),
units.iloc[b][["patient", "unit_id", "brain_region", "n_spikes"]].to_dict())
w0 = {int(s["subject_id"]): tuple(np.round(s["cleaned_watchlogs"].rec_ms.to_numpy()[100:110], 3)) for s in sessions}
print("patients sharing a watchlog:", len(w0) - len(set(w0.values())))
print("units with no spike during playback:", int(sum(len(t) == 0 for t in spikes)))
pairs of units with the same first 20 spike times: 0 patients sharing a watchlog: 0 units with no spike during playback: 2
No two units share their first 20 spike times and no two patients share a watchlog, so no recording appears twice in the pinned version. Two units have no spike during playback.
# Clock check: does each patient's population change its rate after camera cuts at the annotated
# times, and not at the same number of random times? Index = mean over units with 100 or more
# spikes of |rate 200-600 ms after - rate 400-0 ms before| / unit's mean rate; increases and
# decreases both count, as regions differ in sign (pooled PSTH below).
rng_clock = np.random.default_rng(1)
ann = sessions[0]["annotations_base"]
cuts = summer.onsets(ann, "camera-cuts", 1, min_presence=None)
def cut_index(ts, ev):
vals = []
for t in ts:
if len(t) < 100:
continue
rate = len(t) / (t[-1] - t[0])
after = summer._counts_in(t, ev + 0.2, 0.4).mean() / 0.4
before = summer._counts_in(t, ev - 0.4, 0.4).mean() / 0.4
vals.append(abs(after - before) / rate)
return np.mean(vals)
rows = []
for s in sessions:
p = int(s["subject_id"])
ts = [np.sort(spikes[i]) for i in units.index[units.patient == p]]
lo, hi = s["cleaned_watchlogs"].pts.iloc[0] + 1, s["cleaned_watchlogs"].pts.iloc[-1] - 1
ev = cuts[(cuts > lo) & (cuts < hi)]
true = cut_index(ts, ev)
null = np.array([cut_index(ts, np.sort(rng_clock.uniform(lo, hi, len(ev)))) for _ in range(20)])
rows.append({"patient": p, "cuts": len(ev), "index_true": true, "index_shift_2s": cut_index(ts, ev + 2.0),
"null_mean": null.mean(), "null_sd": null.std(), "z": (true - null.mean()) / null.std()})
clock = pd.DataFrame(rows)
print(clock.round(3).to_string(index=False))
print("patients with z above 3:", int((clock.z > 3).sum()), "of", len(clock), "; not above 3:", clock.patient[clock.z <= 3].tolist())
# Pooled PSTH around camera cuts, by region (rate / unit's mean rate, 100 ms bins from -1 to 1 s).
bins = np.arange(-1.0, 1.0001, 0.1)
ps = {}
for region in ("A", "H", "EC", "PHC"):
vals = []
for i in units.index[units.region == region]:
t = np.sort(spikes[i])
if len(t) < 100:
continue
rate = len(t) / (t[-1] - t[0])
vals.append(summer._counts_in(t, (cuts[:, None] + bins[None, :-1]).ravel(), 0.1).reshape(len(cuts), -1).mean(0) / 0.1 / rate)
ps[region] = np.mean(vals, 0)
print("\npooled rate around camera cuts / mean rate, bins starting at (s):")
print(pd.DataFrame(ps, index=np.round(bins[:-1], 1)).T.round(3).to_string())
patient cuts index_true index_shift_2s null_mean null_sd z
1 884 0.127 0.124 0.112 0.022 0.670
2 883 0.114 0.091 0.087 0.009 2.873
3 884 0.160 0.086 0.087 0.014 5.247
4 883 0.073 0.068 0.066 0.006 1.214
5 883 0.122 0.079 0.091 0.010 2.983
6 883 0.118 0.110 0.099 0.010 1.995
7 883 0.112 0.095 0.091 0.008 2.548
8 884 0.087 0.107 0.103 0.010 -1.628
9 884 0.098 0.088 0.111 0.013 -0.999
10 884 0.110 0.064 0.071 0.006 6.301
11 883 0.097 0.112 0.090 0.020 0.344
12 884 0.112 0.104 0.093 0.009 2.223
13 884 0.127 0.086 0.087 0.009 4.485
14 884 0.148 0.147 0.119 0.020 1.393
15 884 0.152 0.103 0.115 0.009 4.090
16 884 0.144 0.089 0.104 0.011 3.759
17 884 0.125 0.086 0.095 0.009 3.239
18 884 0.192 0.141 0.145 0.025 1.878
19 884 0.159 0.099 0.102 0.008 6.761
20 884 0.226 0.097 0.091 0.007 19.908
21 883 0.109 0.090 0.091 0.008 2.274
22 882 0.126 0.101 0.101 0.008 3.183
23 883 0.171 0.102 0.098 0.016 4.684
24 883 0.157 0.133 0.132 0.011 2.226
25 884 0.154 0.117 0.119 0.012 2.897
26 884 0.138 0.135 0.120 0.013 1.402
27 884 0.086 0.075 0.087 0.009 -0.141
28 884 0.066 0.061 0.066 0.006 0.116
29 884 0.167 0.132 0.110 0.012 4.546
patients with z above 3: 11 of 29 ; not above 3: [1, 2, 4, 5, 6, 7, 8, 9, 11, 12, 14, 18, 21, 24, 25, 26, 27, 28]
pooled rate around camera cuts / mean rate, bins starting at (s):
-1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 -0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
A 1.000 0.995 1.009 1.001 1.007 0.999 1.020 1.021 1.028 1.002 1.008 0.987 0.970 0.951 0.959 0.959 0.974 1.014 1.002 1.031
H 1.004 0.988 1.014 1.008 0.996 0.988 1.013 0.996 0.995 0.994 1.003 0.980 0.986 0.968 0.980 0.993 0.994 0.999 1.011 0.994
EC 0.978 0.994 0.979 0.984 1.004 0.998 0.969 0.999 1.023 1.004 1.004 0.983 0.970 0.947 0.971 0.976 0.989 1.005 1.000 1.030
PHC 1.002 0.995 0.988 1.013 0.988 0.997 0.969 0.975 0.983 0.992 0.973 0.995 1.093 1.078 1.050 1.028 1.026 1.028 1.035 1.002
In 11 of 29 patients the population changes its rate after camera cuts more than after random times (z above 3). In the other 18 (patients 1, 2, 4, 5, 6, 7, 8, 9, 11, 12, 14, 18, 21, 24, 25, 26, 27 and 28) this check does not establish the alignment, though it does not contradict it. Pooled over patients, the PHC rate rises to 1.09 times its mean 200 to 300 ms after a cut, and the A and EC rates fall to 0.95 at 300 to 400 ms. Decoding Summer with its labels shifted by 300 s gives kappa 0.015 against 0.301 with the true labels (fig3c section); per patient, the true labels give the higher kappa in 18 of 29 patients (fig4 section). The spikes and annotations are consistent with one clock: shown for the pooled data, and for 11 patients one by one.
Responsive units¶
Shared by the next six claims. For every unit and the 11 labels of Table S1, the Methods' test as
fixed in the claims file: onsets preceded by at least 1 s without the label and, for labels with a
duration, lasting at least 1 s; the code's onset rule (no minimum duration) is computed as a
variant. Indoor and Outdoor are the values 1 and 0 of indoor-setting.
# Check of the onset rule on a small table: segments at 5-8 s, 8.5-12 s and 20-21 s give onsets at
# 5 s (first segment, nothing before it) and 20 s; 8.5 s follows the first segment by 0.5 s.
toy = pd.DataFrame({"label": "x", "value": 1, "start": [5.0, 8.5, 20.0], "stop": [8.0, 12.0, 21.0]})
assert list(summer.onsets(toy, "x", 1, min_presence=None)) == [5.0, 20.0]
assert list(summer.onsets(toy, "x", 1, min_presence=2.0)) == [5.0]
LABEL_ORDER = ["Tom", "Summer", "McKenzie", "Persons", "Tom Faces", "Summer Faces", "Summer Presence",
"Camera Cuts", "Scene Cuts", "Indoor", "Outdoor"]
events, events_code = {}, {}
for name in LABEL_ORDER:
label, value = summer.LABELS[name]
events[name] = summer.onsets(ann, label, value, min_presence=None if label in summer.CUT_LABELS else 1.0)
events_code[name] = summer.onsets(ann, label, value, min_presence=None)
print(pd.DataFrame({"onsets (Methods rule)": {k: len(v) for k, v in events.items()},
"onsets (code rule)": {k: len(v) for k, v in events_code.items()}}).to_string())
started = time.time()
jobs = []
for s in sessions:
p = int(s["subject_id"])
idx = units.index[units.patient == p].to_numpy()
rng_ = (s["cleaned_watchlogs"].pts.iloc[0], s["cleaned_watchlogs"].pts.iloc[-1])
jobs.append((summer.responsive_table, ([spikes[i] for i in idx], rng_, events)))
jobs.append((summer.responsive_table, ([spikes[i] for i in idx], rng_,
{k: v for k, v in events_code.items() if summer.LABELS[k][0] not in summer.CUT_LABELS})))
out = summer.run_jobs(jobs)
resp_tables = []
for j, s in enumerate(sessions):
idx = units.index[units.patient == int(s["subject_id"])].to_numpy()
for k, rule in ((2 * j, "methods"), (2 * j + 1, "code")):
df = out[k].copy()
df["unit"] = idx[df.unit.to_numpy()]
df["rule"] = rule
resp_tables.append(df)
resp = pd.concat(resp_tables, ignore_index=True)
resp = resp.merge(units[["patient", "region", "brain_region"]], left_on="unit", right_index=True)
print(f"responsive tests done in {time.time() - started:.0f} s; rows {len(resp)}")
P = resp[resp.rule == "methods"].pivot(index="unit", columns="label", values="p")
print("units with p = 1 because they spike in fewer than a third of the pseudo-trials, by label:")
print((resp[resp.rule == "methods"].assign(low=lambda d: d.active < 1 / 3).groupby("label").low.sum()).to_string())
def responsive_counts(pv, alpha):
r = pv <= alpha
by_region = r.groupby(units.region).sum()
return r.sum(), by_region
onsets (Methods rule) onsets (code rule) Tom 152 158 Summer 114 115 McKenzie 50 52 Persons 68 71 Tom Faces 205 403 Summer Faces 138 255 Summer Presence 40 40 Camera Cuts 884 884 Scene Cuts 133 133 Indoor 44 44 Outdoor 48 48
responsive tests done in 122 s; rows 45720 units with p = 1 because they spike in fewer than a third of the pseudo-trials, by label: label Camera Cuts 548 Indoor 553 McKenzie 551 Outdoor 550 Persons 556 Scene Cuts 561 Summer 543 Summer Faces 559 Summer Presence 559 Tom 558 Tom Faces 559
The Methods' test is checked against chance: the same test with random event times (as many as there are camera cut onsets, within each patient's playback), for every unit; then with as many random times as each label has onsets; and on synthetic units firing as stationary Poisson processes with random events.
rng_calib = np.random.default_rng(21)
null_rows = []
for s in sessions:
p = int(s["subject_id"])
lo, hi = s["cleaned_watchlogs"].pts.iloc[0] + 1, s["cleaned_watchlogs"].pts.iloc[-1] - 1
ev = np.sort(rng_calib.uniform(lo, hi, len(events["Camera Cuts"])))
for i in units.index[units.patient == p]:
null_rows.append({"unit": i, "region": units.region[i], "p": summer.wilcoxon_simes(spikes[i], ev)[0]})
null = pd.DataFrame(null_rows)
print("Methods' test at random times, fraction of units 'responsive':")
print(null.groupby("region").p.agg(**{"alpha_0.001": lambda v: np.mean(v <= 0.001), "alpha_0.01": lambda v: np.mean(v <= 0.01),
"units": "size"}).round(3).to_string())
print("all units:", {a: round(float(np.mean(null.p <= a)), 3) for a in summer.ALPHAS})
null_p = null.set_index("unit").p
def null_by_label(spk, lo, hi, counts, seed):
r = np.random.default_rng(seed)
out = {}
for name, k in counts.items():
ev = np.sort(r.uniform(lo, hi, k))
out[name] = [summer.wilcoxon_simes(t, ev)[0] for t in spk]
return out
# Random times, as many as each label's onsets (Methods rule), every unit; one seed per patient.
counts_ = {name: len(events[name]) for name in LABEL_ORDER}
jobs = []
for s in sessions:
p = int(s["subject_id"])
lo, hi = s["cleaned_watchlogs"].pts.iloc[0] + 1, s["cleaned_watchlogs"].pts.iloc[-1] - 1
jobs.append((null_by_label, ([spikes[i] for i in units.index[units.patient == p]], lo, hi, counts_, 2100 + p)))
out = summer.run_jobs(jobs)
null_lab = pd.concat([pd.DataFrame(o, index=units.index[units.patient == int(s["subject_id"])]) for o, s in zip(out, sessions)])
print("\nMethods' test at random times, as many as each label's onsets, fraction of all units 'responsive':")
print(pd.DataFrame({"onsets": counts_, **{f"alpha_{a}": (null_lab <= a).mean().round(4) for a in summer.ALPHAS}}).to_string())
synth = []
for rate in (2.0, 10.0):
for n_ev in (40, 114, 884):
ps_ = []
for _ in range(100):
t = np.sort(rng_calib.uniform(0, 5000, rng_calib.poisson(rate * 5000)))
ps_.append(summer.wilcoxon_simes(t, np.sort(rng_calib.uniform(10, 4990, n_ev)))[0])
synth.append({"rate_hz": rate, "events": n_ev, "alpha_0.001": np.mean(np.array(ps_) <= 0.001),
"alpha_0.01": np.mean(np.array(ps_) <= 0.01)})
print("\nsynthetic Poisson units, 100 each, fraction 'responsive':")
print(pd.DataFrame(synth).to_string(index=False))
Methods' test at random times, fraction of units 'responsive':
alpha_0.001 alpha_0.01 units
region
A 0.710 0.733 559
EC 0.734 0.744 493
H 0.618 0.658 846
PHC 0.724 0.741 293
PIC 0.765 0.794 68
other 0.667 0.704 27
all units: {0.001: 0.684, 0.01: 0.71}
Methods' test at random times, as many as each label's onsets, fraction of all units 'responsive':
onsets alpha_0.001 alpha_0.01
Tom 152 0.5026 0.5914
Summer 114 0.4405 0.5665
McKenzie 50 0.1715 0.3876
Persons 68 0.2598 0.4598
Tom Faces 205 0.5547 0.6238
Summer Faces 138 0.4873 0.5866
Summer Presence 40 0.0954 0.3066
Camera Cuts 884 0.6920 0.7270
Scene Cuts 133 0.4825 0.5731
Indoor 44 0.1146 0.3508
Outdoor 48 0.1479 0.3758
synthetic Poisson units, 100 each, fraction 'responsive':
rate_hz events alpha_0.001 alpha_0.01
2.0 40 0.14 0.46
2.0 114 0.64 0.97
2.0 884 1.00 1.00
10.0 40 0.00 0.01
10.0 114 0.02 0.03
10.0 884 0.37 0.68
fig2b-phc-camera-cuts¶
Over half of the parahippocampal neurons responded to camera cuts (193 of 373, 51.74%).
phc = units.index[units.region == "PHC"]
cc_p = P.loc[phc, "Camera Cuts"]
print(f"PHC units: {len(phc)} (paper: 373)")
for alpha in summer.ALPHAS:
n = int((cc_p <= alpha).sum())
print(f"alpha {alpha}: {n} of {len(phc)} responsive, {100 * n / len(phc):.2f}% (paper: 193 of 373, 51.74%; tolerance 46.57 to 56.91%)")
print("\nThreshold sweep (fraction of PHC units, %):",
{a: round(100 * float((cc_p <= a).mean()), 2) for a in (1e-4, 1e-3, 1e-2, 5e-2)})
print("random event times (calibration above), PHC units responsive (%):",
{a: round(100 * float((null_p[phc] <= a).mean()), 2) for a in summer.ALPHAS})
print("\nPer patient, PHC units responsive to camera cuts at alpha 0.001 and 0.01:")
per = pd.DataFrame({"units": cc_p.groupby(units.patient).size(),
"resp_0.001": (cc_p <= 0.001).groupby(units.patient).sum(),
"resp_0.01": (cc_p <= 0.01).groupby(units.patient).sum()})
per["pct_0.001"] = 100 * per["resp_0.001"] / per.units
print(per.round(1).to_string())
PHC units: 293 (paper: 373)
alpha 0.001: 222 of 293 responsive, 75.77% (paper: 193 of 373, 51.74%; tolerance 46.57 to 56.91%)
alpha 0.01: 228 of 293 responsive, 77.82% (paper: 193 of 373, 51.74%; tolerance 46.57 to 56.91%)
Threshold sweep (fraction of PHC units, %): {0.0001: 72.7, 0.001: 75.77, 0.01: 77.82, 0.05: 79.18}
random event times (calibration above), PHC units responsive (%): {0.001: 72.35, 0.01: 74.06}
Per patient, PHC units responsive to camera cuts at alpha 0.001 and 0.01:
units resp_0.001 resp_0.01 pct_0.001
patient
1 6 6 6 100.0
2 10 5 5 50.0
3 4 1 1 25.0
4 10 8 8 80.0
5 7 5 5 71.4
6 20 19 20 95.0
7 1 0 0 0.0
9 11 6 6 54.5
10 15 13 13 86.7
11 4 0 0 0.0
12 4 4 4 100.0
13 19 19 19 100.0
14 11 5 5 45.5
15 19 9 11 47.4
17 12 11 11 91.7
19 6 4 4 66.7
20 16 14 14 87.5
21 15 14 14 93.3
22 2 0 0 0.0
23 9 8 8 88.9
25 24 17 18 70.8
26 14 12 12 85.7
27 27 20 22 74.1
28 27 22 22 81.5
The ZETA test of the dataset's code, shown here for PHC in the first round, is now run for every
unit and 9 labels of Table S1 in the companion notebooks responsive_zeta_cuts.ipynb and
responsive_zeta_labels.ipynb, each with its own check at random times. Tom Faces and Summer Faces
are not run; no responsive-unit claim uses them.
Result: blocked (methods_underspecified). The Methods' test, as fixed in the claims file,
counts 222 of 293 PHC units (75.77%) responsive to camera cuts at alpha 0.001 and 228 (77.82%) at
0.01, against 51.74% in the paper. The same test counts 72.35% of PHC units, and 68.4% of all
units, responsive at random event times (69.2% with as many random times as there are camera cut
onsets, in the per-label check), and every synthetic Poisson unit at 2 Hz with 884 random events.
Its false-positive rate is far from its alpha, so the procedure here is not the one that produced
the paper's counts, and the claim cannot be judged. The ZETA test of the dataset's code
(responsive_zeta_cuts.ipynb) is not anticonservative at random times (none of 400 units at alpha
0.001, 1.0% at 0.01); it finds 35 of 293 PHC units (11.95%) responsive at 0.001 and 67 (22.87%) at 0.01,
well below the paper's 51.74%, so the paper's count does not come from that test as run here either.
tables1-camera-cuts-regions¶
Camera cuts evoked responses in a large share of neurons in every region (amygdala 46.90%, entorhinal cortex 45.68%, hippocampus 38.92%).
paper_cc = {"A": 46.90, "EC": 45.68, "H": 38.92, "PHC": 51.74}
rows = []
for region in ("A", "EC", "H", "PHC", "PIC", "other"):
idx = units.index[units.region == region]
for alpha in summer.ALPHAS:
n = int((P.loc[idx, "Camera Cuts"] <= alpha).sum())
rows.append({"region": region, "alpha": alpha, "units": len(idx), "responsive": n,
"pct": 100 * n / len(idx), "paper_pct": paper_cc.get(region, np.nan)})
cc_reg = pd.DataFrame(rows)
cc_reg["within_10pct"] = (cc_reg.pct - cc_reg.paper_pct).abs() <= 0.1 * cc_reg.paper_pct
print(cc_reg.round(2).to_string(index=False))
tot = {a: int((P["Camera Cuts"] <= a).sum()) for a in summer.ALPHAS}
print("all units:", {a: f"{n} of {len(P)} ({100 * n / len(P):.2f}%)" for a, n in tot.items()}, "(paper: 1025, 44.84%)")
print("\nPer patient fraction responsive at alpha 0.001 (%), by region (median and range over patients with 5 or more units):")
fr = (P["Camera Cuts"] <= 0.001).groupby([units.region, units.patient]).agg(["mean", "size"])
fr = fr[fr["size"] >= 5]
print((100 * fr["mean"]).groupby(level=0).describe()[["count", "min", "50%", "max"]].round(1).to_string())
region alpha units responsive pct paper_pct within_10pct
A 0.00 559 412 73.70 46.90 False
A 0.01 559 421 75.31 46.90 False
EC 0.00 493 363 73.63 45.68 False
EC 0.01 493 373 75.66 45.68 False
H 0.00 846 546 64.54 38.92 False
H 0.01 846 583 68.91 38.92 False
PHC 0.00 293 222 75.77 51.74 False
PHC 0.01 293 228 77.82 51.74 False
PIC 0.00 68 54 79.41 NaN False
PIC 0.01 68 55 80.88 NaN False
other 0.00 27 19 70.37 NaN False
other 0.01 27 20 74.07 NaN False
all units: {0.001: '1616 of 2286 (70.69%)', 0.01: '1680 of 2286 (73.49%)'} (paper: 1025, 44.84%)
Per patient fraction responsive at alpha 0.001 (%), by region (median and range over patients with 5 or more units):
count min 50% max
region
A 27.0 42.9 71.4 94.7
EC 23.0 36.4 76.5 93.8
H 28.0 40.0 63.3 83.3
PHC 19.0 45.5 81.5 100.0
PIC 4.0 71.4 82.9 92.0
other 3.0 33.3 80.0 83.3
Result: blocked (methods_underspecified), as it rests on the selection of
fig2b-phc-camera-cuts. The miscalibrated test counts 64 to 76% of the units of every region
responsive (A 73.70%, EC 73.63%, H 64.54%, PHC 75.77% at alpha 0.001), against 39 to 52% in the
paper, and 62 to 73% of the units of each region at random times. The ZETA test of the dataset's
code counts A 3.22%, EC 0.61% and H 1.30% (responsive_zeta_cuts.ipynb).
fig2b-phc-scene-cuts¶
Few parahippocampal neurons responded to scene cuts (25 of 373, 6.70%), but the modulation was consistent across the responsive neurons (cluster permutation test, p = 0.018).
rng_scene = np.random.default_rng(23)
psth = {rule: {name: np.vstack(resp[(resp.rule == rule) & (resp.label == name)].set_index("unit").loc[units.index, "psth"])
for name in resp[resp.rule == rule].label.unique()} for rule in ("methods",)}
psth_trial = {name: np.vstack(resp[(resp.rule == "methods") & (resp.label == name)].set_index("unit").loc[units.index, "psth_trial_z"])
for name in LABEL_ORDER}
phc_mask = (units.region == "PHC").to_numpy()
sc_p = P.loc[phc, "Scene Cuts"]
for alpha in summer.ALPHAS:
r = (sc_p <= alpha).to_numpy()
obs, pc = summer.cluster_test(psth["methods"]["Scene Cuts"][phc_mask], r, rng_scene)
obs2, pc2 = summer.cluster_test(psth_trial["Scene Cuts"][phc_mask], r, rng_scene)
print(f"alpha {alpha}: {int(r.sum())} of {len(r)} PHC units responsive ({100 * r.mean():.2f}%; paper 25 of 373, 6.70%, "
f"tolerance 4.16 to 9.24%); largest cluster sum |t| {obs:.1f}, p = {pc:.3f}; "
f"with the code's per-trial z-scoring {obs2:.1f}, p = {pc2:.3f}")
print("threshold sweep (%):", {a: round(100 * float((sc_p <= a).mean()), 2) for a in (1e-4, 1e-3, 1e-2, 5e-2)})
print("onsets used:", len(events["Scene Cuts"]))
alpha 0.001: 138 of 293 PHC units responsive (47.10%; paper 25 of 373, 6.70%, tolerance 4.16 to 9.24%); largest cluster sum |t| 3.0, p = 0.055; with the code's per-trial z-scoring 3.2, p = 0.020
alpha 0.01: 180 of 293 PHC units responsive (61.43%; paper 25 of 373, 6.70%, tolerance 4.16 to 9.24%); largest cluster sum |t| 0.0, p = 1.000; with the code's per-trial z-scoring 0.0, p = 1.000
threshold sweep (%): {0.0001: 34.47, 0.001: 47.1, 0.01: 61.43, 0.05: 68.6}
onsets used: 133
Result: blocked (methods_underspecified). The Methods' test counts 138 of 293 PHC units
(47.10%) responsive to scene cuts at alpha 0.001 and 180 (61.43%) at 0.01, against 6.70%; with as
many random times as there are scene cut onsets it counts 48.25% of all units, so the test is not
calibrated. With this split, the largest cluster has p = 0.055 with the Methods' z-scoring and
p = 0.020 with the code's per-trial z-scoring at alpha 0.001; at 0.01 there is no cluster. The ZETA
test finds 12 of 293 PHC units (4.10%) at alpha 0.001 and 29 (9.90%) at 0.01, either side of the
tolerance of 4.16 to 9.24%; with 1000 resamplings in place of 100 these are 11 and 31
(responsive_zeta_cuts.ipynb).
tables1-few-semantic-responses¶
Only a few neurons responded to the onset of semantic features (Summer 46, Tom 216, Persons 6 and Outdoor 4 of 2286 neurons), in contrast to camera cuts (1025).
table_s1 = {"Tom": 216, "Summer": 46, "McKenzie": 0, "Persons": 6, "Tom Faces": 843, "Summer Faces": 634,
"Summer Presence": 0, "Camera Cuts": 1025, "Scene Cuts": 106, "Indoor": 0, "Outdoor": 4}
Pc = resp[resp.rule == "code"].pivot(index="unit", columns="label", values="p")
rows = []
for name in LABEL_ORDER:
row = {"label": name, "onsets": len(events[name]), "paper_total": table_s1[name]}
for alpha in summer.ALPHAS:
r = P[name] <= alpha
row[f"total_{alpha}"] = int(r.sum())
row[f"pct_{alpha}"] = round(100 * float(r.mean()), 2)
for region in ("A", "EC", "H", "PHC"):
row[f"{region}_{alpha}"] = int(r[units.region == region].sum())
if name in Pc:
row[f"code_rule_total_{alpha}"] = int((Pc[name] <= alpha).sum())
rows.append(row)
s1 = pd.DataFrame(rows)
print(s1.to_string(index=False))
print("\nTolerance: Summer 46 +- 13.2 (32.8 to 59.2), Tom 216 +- 21.6 (194.4 to 237.6).")
semantic = ["Summer", "Tom", "McKenzie", "Persons", "Indoor", "Outdoor", "Summer Presence"]
for alpha in summer.ALPHAS:
print(f"alpha {alpha}: largest semantic fraction {s1.set_index('label').loc[semantic, f'pct_{alpha}'].max():.2f}% "
f"({s1.set_index('label').loc[semantic, f'pct_{alpha}'].idxmax()}); camera cuts {s1.set_index('label').loc['Camera Cuts', f'pct_{alpha}']:.2f}%")
print("\nPer patient, units responsive to Summer and Tom at alpha 0.001:")
print(pd.DataFrame({"units": units.groupby("patient").size(),
"Summer": (P["Summer"] <= 0.001).groupby(units.patient).sum(),
"Tom": (P["Tom"] <= 0.001).groupby(units.patient).sum()}).T.to_string())
label onsets paper_total total_0.001 pct_0.001 A_0.001 EC_0.001 H_0.001 PHC_0.001 code_rule_total_0.001 total_0.01 pct_0.01 A_0.01 EC_0.01 H_0.01 PHC_0.01 code_rule_total_0.01
Tom 152 216 1174 51.36 316 267 394 143 1191.0 1370 59.93 365 309 454 181 1384.0
Summer 114 46 1013 44.31 278 224 343 126 1019.0 1289 56.39 341 290 432 166 1285.0
McKenzie 50 0 439 19.20 128 86 153 56 440.0 923 40.38 249 217 310 111 954.0
Persons 68 6 658 28.78 197 130 212 87 722.0 1090 47.68 294 261 352 138 1101.0
Tom Faces 205 843 1302 56.96 335 304 433 173 1431.0 1449 63.39 370 335 486 193 1542.0
Summer Faces 138 634 1100 48.12 304 243 361 135 1330.0 1321 57.79 350 289 437 179 1455.0
Summer Presence 40 0 227 9.93 67 44 88 22 227.0 734 32.11 200 168 248 87 734.0
Camera Cuts 884 1025 1616 70.69 412 363 546 222 NaN 1680 73.49 421 373 583 228 NaN
Scene Cuts 133 106 1121 49.04 307 251 374 138 NaN 1343 58.75 349 300 453 180 NaN
Indoor 44 0 299 13.08 90 62 101 35 299.0 810 35.43 212 193 271 102 810.0
Outdoor 48 4 375 16.40 114 65 126 52 375.0 873 38.19 228 188 296 117 873.0
Tolerance: Summer 46 +- 13.2 (32.8 to 59.2), Tom 216 +- 21.6 (194.4 to 237.6).
alpha 0.001: largest semantic fraction 51.36% (Tom); camera cuts 70.69%
alpha 0.01: largest semantic fraction 59.93% (Tom); camera cuts 73.49%
Per patient, units responsive to Summer and Tom at alpha 0.001:
patient 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
units 40 83 51 120 79 99 84 45 41 108 31 52 137 30 97 67 101 40 105 88 109 79 50 88 110 85 90 111 66
Summer 17 47 20 51 41 29 34 23 23 49 7 26 60 11 47 29 39 12 39 46 54 28 23 37 50 39 41 54 37
Tom 24 51 22 64 44 50 34 27 22 51 9 32 62 13 48 32 51 16 50 48 61 35 28 48 61 45 48 58 40
Result: blocked (methods_underspecified). The test counts 1013 units responsive to Summer and
1174 to Tom at alpha 0.001 (paper 46 and 216), and 9.93 to 51.36% of all units for each semantic
label. With as many random times as each label has onsets it counts about as many: 44.05% of units
for Summer's 114 (44.31% at the real onsets), 50.26% for Tom's 152 (51.36%), 11.46% for Indoor's 44
(13.08%). The counts are dominated by false positives, so neither the counts nor the contrast with
camera cuts (70.69%) can be judged. The ZETA test of the dataset's code, conservative at random times
(none of 800 units at alpha 0.001, 0.12 to 0.38% at 0.01), finds 11 units for Summer and 13 for Tom
(responsive_zeta_labels.ipynb), below the paper's 46 and 216.
tables1-no-indoor-response¶
No neuron responded to the onset of indoor scenes.
for name in ("Indoor", "McKenzie", "Summer Presence"):
print(name, f"({len(events[name])} onsets):", {a: int((P[name] <= a).sum()) for a in summer.ALPHAS},
"| code onset rule:", {a: int((Pc[name] <= a).sum()) for a in summer.ALPHAS})
print("smallest Indoor p:", f"{P['Indoor'].min():.2g}", "; units testable for Indoor (spiking in a third of the pseudo-trials):",
int((resp[(resp.rule == 'methods') & (resp.label == 'Indoor')].active >= 1 / 3).sum()))
print("reference, all units at random times (calibration above), responsive:", {a: f"{100 * float((null_p <= a).mean()):.2f}%" for a in summer.ALPHAS})
Indoor (44 onsets): {0.001: 299, 0.01: 810} | code onset rule: {0.001: 299, 0.01: 810}
McKenzie (50 onsets): {0.001: 439, 0.01: 923} | code onset rule: {0.001: 440, 0.01: 954}
Summer Presence (40 onsets): {0.001: 227, 0.01: 734} | code onset rule: {0.001: 227, 0.01: 734}
smallest Indoor p: 4.2e-07 ; units testable for Indoor (spiking in a third of the pseudo-trials): 1733
reference, all units at random times (calibration above), responsive: {0.001: '68.42%', 0.01: '71.04%'}
Result: blocked (methods_underspecified). The test as fixed counts 299 units responsive to
Indoor onsets at alpha 0.001 (810 at 0.01), and 439 for McKenzie and 227 for Summer Presence, where
Table S1 has 0 for each. With 44 random times in place of the Indoor onsets it counts 11.46% of all
units. An absence cannot be judged with a test whose
false-positive rate is this far from its alpha. The ZETA test of the dataset's code finds 2 units
for Indoor at alpha 0.001 (responsive_zeta_labels.ipynb); it is not the paper's test.
results-cluster-regions¶
By the cluster permutation test, parahippocampal neurons responded to the largest set of features (Summer, Persons, Camera Cuts, Scene Cuts, Outdoor), compared with the hippocampus (Tom, Camera Cuts, Outdoor), the amygdala (Camera Cuts) and the entorhinal cortex (none).
rng_cluster = np.random.default_rng(24)
paper_sig = {"PHC": {"Summer", "Persons", "Camera Cuts", "Scene Cuts", "Outdoor"},
"H": {"Tom", "Camera Cuts", "Outdoor"}, "A": {"Camera Cuts"}, "EC": set()}
CL = ["Summer", "Tom", "McKenzie", "Persons", "Camera Cuts", "Scene Cuts", "Indoor", "Outdoor"]
rows = []
started = time.time()
for alpha in summer.ALPHAS:
for region in ("A", "H", "EC", "PHC"):
m = (units.region == region).to_numpy()
for name in CL:
r = (P.loc[units.index[m], name] <= alpha).to_numpy()
obs, pc = summer.cluster_test(psth["methods"][name][m], r, rng_cluster)
rows.append({"alpha": alpha, "region": region, "label": name, "responsive": int(r.sum()),
"cluster_sum": obs, "p": pc, "significant": pc < 0.05, "paper": name in paper_sig[region]})
cl = pd.DataFrame(rows)
cl["agree"] = cl.significant == cl.paper
print(f"cluster tests in {time.time() - started:.0f} s")
for alpha in summer.ALPHAS:
t = cl[cl.alpha == alpha]
print(f"\nalpha {alpha}: cells agreeing with the paper {int(t.agree.sum())} of {len(t)}")
print(t.pivot(index="label", columns="region", values="p").round(3).to_string())
print("significant labels per region:", t[t.significant].groupby("region").label.apply(list).to_dict())
print("responsive units per cell:")
print(t.pivot(index="label", columns="region", values="responsive").to_string())
cluster tests in 34 s
alpha 0.001: cells agreeing with the paper 11 of 32
region A EC H PHC
label
Camera Cuts 1.000 1.000 1.000 1.000
Indoor 0.033 1.000 0.000 0.037
McKenzie 0.000 0.042 0.002 0.001
Outdoor 0.000 0.011 0.000 0.001
Persons 0.000 0.009 0.000 0.000
Scene Cuts 0.000 1.000 0.000 0.060
Summer 0.001 0.054 0.001 0.039
Tom 1.000 0.077 0.044 0.025
significant labels per region: {'A': ['Summer', 'McKenzie', 'Persons', 'Scene Cuts', 'Indoor', 'Outdoor'], 'EC': ['McKenzie', 'Persons', 'Outdoor'], 'H': ['Summer', 'Tom', 'McKenzie', 'Persons', 'Scene Cuts', 'Indoor', 'Outdoor'], 'PHC': ['Summer', 'Tom', 'McKenzie', 'Persons', 'Indoor', 'Outdoor']}
responsive units per cell:
region A EC H PHC
label
Camera Cuts 412 363 546 222
Indoor 90 62 101 35
McKenzie 128 86 153 56
Outdoor 114 65 126 52
Persons 197 130 212 87
Scene Cuts 307 251 374 138
Summer 278 224 343 126
Tom 316 267 394 143
alpha 0.01: cells agreeing with the paper 12 of 32
region A EC H PHC
label
Camera Cuts 1.000 1.000 1.000 1.000
Indoor 0.000 1.000 0.000 0.000
McKenzie 0.005 1.000 0.001 0.001
Outdoor 0.022 0.015 0.000 0.000
Persons 0.048 1.000 0.000 0.002
Scene Cuts 0.000 1.000 0.000 1.000
Summer 1.000 1.000 0.033 0.091
Tom 1.000 0.003 0.083 1.000
significant labels per region: {'A': ['McKenzie', 'Persons', 'Scene Cuts', 'Indoor', 'Outdoor'], 'EC': ['Tom', 'Outdoor'], 'H': ['Summer', 'McKenzie', 'Persons', 'Scene Cuts', 'Indoor', 'Outdoor'], 'PHC': ['McKenzie', 'Persons', 'Indoor', 'Outdoor']}
responsive units per cell:
region A EC H PHC
label
Camera Cuts 421 373 583 228
Indoor 212 193 271 102
McKenzie 249 217 310 111
Outdoor 228 188 296 117
Persons 294 261 352 138
Scene Cuts 349 300 453 180
Summer 341 290 432 166
Tom 365 309 454 181
Result: blocked (methods_underspecified), as the split into responsive and non-responsive units comes from the miscalibrated test. With that split, 11 of 32 cells agree with the paper at alpha 0.001 and 12 at 0.01. Camera cuts give no cluster in any region (p = 1.000), and at alpha 0.001 PHC has 6 significant labels against 7 in H, 6 in A and 3 in EC.
Decoding data¶
Samples are the 59080 bins of 80 ms of the analysis portion; each unit's feature is its spike count from 800 ms before to 800 ms after the bin. The labels are those of the decoding claims.
edges = sessions[0]["edges_s"]
counts = summer.binned_counts(spikes, edges)
S = summer.sample_windows(spikes, edges)
n = S.shape[0]
DEC = {"Summer": "summer", "Tom": "tom", "McKenzie": "mckenzie", "Persons": "persons",
"Indoor/Outdoor": "indoor-setting", "Scene Cuts": "scenes", "Camera Cuts": "camera-cuts"}
Y = {name: summer.bin_labels(sessions[0], label) for name, label in DEC.items()}
print("samples x units:", S.shape, "; spikes per unit per bin, median:", float(np.median(counts.mean(1))))
print(pd.DataFrame({name: {"positive": float(np.mean(y == 1)), "negative": float(np.mean(y == 0)),
"excluded": float(np.mean((y != 0) & (y != 1)))} for name, y in Y.items()}).round(3).to_string())
for f in range(5):
tr, va, te = summer.split(n, f)
gap_tr_te = np.min(np.abs(edges[tr][:, None] - edges[te][None, ::50])) if f == 0 else None
print(f"fold {f + 1}: train {len(tr)}, validation {len(va)}, test {len(te)}"
+ (f"; smallest distance from a training to a test sample, every 50th test sample checked: {gap_tr_te:.1f} s" if f == 0 else ""))
samples x units: (59080, 2286) ; spikes per unit per bin, median: 0.07335815842924848
Summer Tom McKenzie Persons Indoor/Outdoor Scene Cuts Camera Cuts
positive 0.527 0.792 0.09 0.946 0.672 0.014 0.105
negative 0.473 0.208 0.91 0.054 0.321 0.986 0.895
excluded 0.000 0.000 0.00 0.000 0.007 0.000 0.000
fold 1: train 30822, validation 7194, test 7194; smallest distance from a training to a test sample, every 50th test sample checked: 34.5 s
fold 2: train 30822, validation 7194, test 7194
fold 3: train 30822, validation 7194, test 7194
fold 4: train 30822, validation 7194, test 7194
fold 5: train 30822, validation 7194, test 7194
fig3b-single-neuron-threshold¶
Thresholding the activity of single neurons does not reliably predict whether Summer is on screen; most neurons perform near chance, below the population decoder.
y = Y["Summer"]
folds = [summer.split(n, f) for f in range(5)]
started = time.time()
single = np.array([[summer.threshold_unit(S[:, i], y, va, te) for (tr, va, te) in folds] for i in range(S.shape[1])])
mean_k = single.mean(1)
print(f"single-unit thresholds in {time.time() - started:.0f} s")
print(f"median unit kappa {np.median(mean_k):.4f}; fraction of units below 0.05: {np.mean(mean_k < 0.05):.3f}; "
f"best unit {mean_k.max():.3f} (patient {units.patient[mean_k.argmax()]}, {units.region[mean_k.argmax()]}); "
f"units at or above 0.31: {int((mean_k >= 0.31).sum())}")
print("quantiles of unit kappa (5, 25, 50, 75, 95, 99%):", np.round(np.quantile(mean_k, [.05, .25, .5, .75, .95, .99]), 3))
print("by region, median and max:")
print(pd.Series(mean_k).groupby(units.region).agg(["median", "max", "size"]).round(3).to_string())
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist(mean_k, bins=60)
ax.axvline(0.31, color="k", ls="--", label="paper LSTM, all units")
ax.set_xlabel("single-unit kappa, Summer (mean of 5 folds)")
ax.set_ylabel("units")
ax.legend()
plt.show()
print(f"plotted range {mean_k.min():.3f} to {mean_k.max():.3f}; all points shown")
single-unit thresholds in 6 s
median unit kappa 0.0076; fraction of units below 0.05: 0.898; best unit 0.276 (patient 7, A); units at or above 0.31: 0
quantiles of unit kappa (5, 25, 50, 75, 95, 99%): [-0.04 -0.01 0.008 0.027 0.068 0.116]
by region, median and max:
median max size
region
A 0.009 0.276 559
EC 0.007 0.120 493
H 0.007 0.194 846
PHC 0.010 0.118 293
PIC 0.006 0.102 68
other 0.017 0.078 27
plotted range -0.125 to 0.276; all points shown
Result: reproduced. The median unit's kappa for Summer is 0.0076, 89.8% of units are below 0.05, and the best unit reaches 0.276 (patient 7, amygdala). No unit reaches the paper's population value of 0.31, nor the pooled logistic regression here (0.301).
fig3c-logreg-characters-location and fig3c-logreg-cuts-lower¶
A logistic regression on spike counts decodes characters and location about as well as the LSTM; pooled over patients it reaches kappa 0.28 for Summer and 0.31 for Indoor/Outdoor.
For scene and camera cuts the logistic regression performs below the LSTM, by 0.1 and 0.07, that is about kappa 0.10 and 0.11.
rng_lr = np.random.default_rng(31)
started = time.time()
jobs = [(summer.fit_logreg, (S, Y[name], *folds[f])) for name in DEC for f in range(5)]
fits = dict(zip([(name, f) for name in DEC for f in range(5)], summer.run_jobs(jobs)))
print(f"{len(jobs)} pooled fits in {time.time() - started:.0f} s")
lstm = {"Summer": (0.3086, 0.0551), "Tom": (0.2300, 0.0604), "McKenzie": (0.3162, 0.0815), "Persons": (0.3599, 0.0477),
"Indoor/Outdoor": (0.3115, 0.0597), "Scene Cuts": (0.2007, 0.0321), "Camera Cuts": (0.1833, 0.0149)}
rows = []
for name in DEC:
ks = np.array([fits[(name, f)]["kappa"] for f in range(5)])
perm_p = []
for f in range(5):
r = fits[(name, f)]
null = summer.shuffled_kappas(r["y_test"], r["pred"], rng_lr)
perm_p.append((1 + np.sum(null >= r["kappa"])) / (1 + len(null)))
rows.append({"label": name, "kappa_mean": ks.mean(), "kappa_sem": ks.std(ddof=1) / np.sqrt(5),
"folds": np.round(ks, 3).tolist(), "C": [fits[(name, f)]["C"] for f in range(5)],
"nonzero_units": [int((fits[(name, f)]["coef"] != 0).sum()) for f in range(5)],
"max_perm_p": max(perm_p), "lstm_paper": lstm[name][0], "lstm_sem": lstm[name][1]})
lr = pd.DataFrame(rows).set_index("label")
print(lr.round(4).to_string())
print("\npaper, logistic regression: Summer 0.28 (tolerance 0.225 to 0.335), Indoor/Outdoor 0.31 (0.25 to 0.37); "
"Scene Cuts 0.10 (0.068 to 0.132), Camera Cuts 0.11 (0.095 to 0.125)")
print("thresholds, LSTM less SEM: Tom 0.17, McKenzie 0.235, Persons 0.312")
print("\ntest kappa at each C, mean over folds:")
print(pd.DataFrame({name: {c: np.mean([fits[(name, f)]["kappa_by_C"][c] for f in range(5)]) for c in summer.C_GRID} for name in DEC}).round(3).to_string())
35 pooled fits in 182 s
kappa_mean kappa_sem folds C nonzero_units max_perm_p lstm_paper lstm_sem
label
Summer 0.3011 0.0560 [0.459, 0.299, 0.316, 0.108, 0.324] [0.003, 0.01, 0.003, 0.03, 0.03] [347, 620, 325, 878, 885] 0.0010 0.3086 0.0551
Tom 0.2184 0.0608 [0.275, 0.31, 0.334, 0.002, 0.172] [0.003, 0.003, 0.01, 0.001, 0.001] [373, 363, 633, 75, 88] 0.4476 0.2300 0.0604
McKenzie 0.3239 0.0732 [0.127, 0.552, 0.322, 0.399, 0.22] [0.003, 0.003, 0.003, 0.003, 0.03] [270, 242, 212, 211, 506] 0.0010 0.3162 0.0815
Persons 0.4100 0.0256 [0.413, 0.426, 0.437, 0.313, 0.461] [0.01, 0.01, 0.003, 0.03, 0.03] [320, 383, 169, 435, 444] 0.0010 0.3599 0.0477
Indoor/Outdoor 0.3592 0.0490 [0.377, 0.413, 0.492, 0.314, 0.201] [0.01, 0.01, 0.001, 0.003, 0.001] [533, 526, 63, 240, 71] 0.0010 0.3115 0.0597
Scene Cuts 0.0771 0.0265 [0.084, 0.048, 0.135, -0.008, 0.127] [0.03, 0.01, 0.003, 0.01, 0.01] [419, 312, 185, 300, 294] 0.8881 0.2007 0.0321
Camera Cuts 0.0961 0.0198 [0.098, 0.05, 0.099, 0.166, 0.067] [0.003, 0.001, 0.003, 0.003, 0.01] [329, 35, 361, 364, 895] 0.0010 0.1833 0.0149
paper, logistic regression: Summer 0.28 (tolerance 0.225 to 0.335), Indoor/Outdoor 0.31 (0.25 to 0.37); Scene Cuts 0.10 (0.068 to 0.132), Camera Cuts 0.11 (0.095 to 0.125)
thresholds, LSTM less SEM: Tom 0.17, McKenzie 0.235, Persons 0.312
test kappa at each C, mean over folds:
Summer Tom McKenzie Persons Indoor/Outdoor Scene Cuts Camera Cuts
0.001 0.266 0.200 0.298 0.330 0.388 0.073 0.084
0.003 0.311 0.236 0.330 0.418 0.393 0.104 0.118
0.010 0.310 0.194 0.310 0.418 0.364 0.098 0.096
0.030 0.304 0.174 0.294 0.408 0.349 0.077 0.070
# Clock check by decoding: the same fits for Summer with the labels shifted by 300 s (circularly).
started = time.time()
shifted = np.roll(Y["Summer"], int(300 / 0.08))
sh = summer.run_jobs([(summer.fit_logreg, (S, shifted, *folds[f])) for f in range(5)])
print(f"Summer, labels shifted by 300 s: kappa {np.mean([r['kappa'] for r in sh]):.3f} "
f"(folds {np.round([r['kappa'] for r in sh], 3).tolist()}); true labels {lr.loc['Summer', 'kappa_mean']:.3f}; "
f"{time.time() - started:.0f} s")
Summer, labels shifted by 300 s: kappa 0.015 (folds [0.01, 0.044, -0.113, 0.055, 0.081]); true labels 0.301; 27 s
Result: reproduced. Pooled over patients, the logistic regression gives kappa 0.301 for Summer (SEM 0.056; paper 0.28) and 0.359 for Indoor/Outdoor (SEM 0.049; paper 0.31), both within tolerance. Tom (0.218), McKenzie (0.324) and Persons (0.410) reach the LSTM values less their SEM (0.17, 0.235, 0.312). Every fold is above chance by the test-label shuffle (p = 0.001, the smallest possible with 1000 shuffles) for Summer, McKenzie, Persons and Indoor/Outdoor; for Tom one fold has kappa 0.002 (p = 0.45). The value of C chosen on validation varies across folds (0.001 to 0.03); over the grid, the test kappa ranges from 0.266 to 0.311 for Summer and from 0.349 to 0.393 for Indoor/Outdoor. With the Summer labels shifted by 300 s the kappa is 0.015.
Result: reproduced. The logistic regression gives kappa 0.0771 for Scene Cuts (SEM 0.027) and 0.0961 for Camera Cuts (SEM 0.020), within the tolerances around 0.10 and 0.11 (Camera Cuts by 0.001) and below the LSTM values of 0.20 and 0.18. Camera Cuts are above chance in every fold (p = 0.001); Scene Cuts are not in one fold (kappa -0.008, p = 0.89).
fig3e-temporal-gap¶
Without temporal gaps between training, validation and test samples, decoding of Summer is almost perfect (kappa 0.99), and it falls to 0.31 as the gap grows to 32 s.
rng_gap = np.random.default_rng(35)
y = Y["Summer"]
sizes = tuple(len(i) for i in folds[0])
started = time.time()
jobs = [(summer.fit_logreg, (S, y, *summer.random_split(n, sizes, rng_gap))) for _ in range(5)]
for gap in (4, 8, 16):
jobs += [(summer.fit_logreg, (S, y, *summer.split(n, f, gap_s=gap))) for f in range(5)]
out = summer.run_jobs(jobs)
gap_rows = {"random": [r["kappa"] for r in out[:5]]}
for j, gap in enumerate((4, 8, 16)):
gap_rows[f"{gap} s"] = [r["kappa"] for r in out[5 + 5 * j: 10 + 5 * j]]
gap_rows["32 s"] = [fits[("Summer", f)]["kappa"] for f in range(5)]
paper_gap = {"random": 0.987, "4 s": 0.584, "8 s": 0.539, "16 s": 0.429, "32 s": 0.309}
gt = pd.DataFrame({k: {"kappa_mean": np.mean(v), "sem": np.std(v, ddof=1) / np.sqrt(5), "min": np.min(v), "max": np.max(v),
"paper_lstm": paper_gap[k]} for k, v in gap_rows.items()}).T
print(f"{len(jobs)} fits in {time.time() - started:.0f} s")
print(gt.round(3).to_string())
print("(paper LSTM values from the dataset's code, buffer_lengths/comparison_performance.ipynb)")
20 fits in 79 s
kappa_mean sem min max paper_lstm
random 0.948 0.002 0.942 0.953 0.987
4 s 0.555 0.012 0.531 0.597 0.584
8 s 0.524 0.028 0.445 0.593 0.539
16 s 0.451 0.019 0.383 0.486 0.429
32 s 0.301 0.056 0.108 0.459 0.309
(paper LSTM values from the dataset's code, buffer_lengths/comparison_performance.ipynb)
Result: reproduced. With samples assigned to training, validation and test at random, the logistic regression decodes Summer with kappa 0.948 (5 draws, 0.942 to 0.953). With gaps of 4, 8, 16 and 32 s it gives 0.555, 0.524, 0.451 and 0.301, close to the LSTM values the dataset's code holds for the same gaps (0.987 random, then 0.584, 0.539, 0.429 and 0.309).
fig4-patients-below-pooled¶
Decoding from a single patient's neurons is lower than from the pooled population; no patient matches the pooled performance, for Summer, Scene Cuts or Indoor/Outdoor.
rng_pat = np.random.default_rng(33)
started = time.time()
PAT_LABELS = ("Summer", "Scene Cuts", "Indoor/Outdoor")
Yp = {k: Y[k] for k in PAT_LABELS} | {"Summer shifted 300 s": shifted}
pats = sorted(units.patient.unique())
out = summer.run_jobs([(summer.fit_folds, (S, units.index[units.patient == p].to_numpy(), Yp, folds)) for p in pats])
pfits = {(p, name, f): r for p, o in zip(pats, out) for (name, f), r in o.items()}
print(f"{len(pfits)} per-patient fits in {time.time() - started:.0f} s")
pt = pd.DataFrame({name: {p: np.mean([pfits[(p, name, f)]["kappa"] for f in range(5)]) for p in pats} for name in Yp})
pt.insert(0, "units", units.groupby("patient").size())
pperm = pd.DataFrame({f"{name} folds p<=0.001": {p: sum(summer.permutation_p(pfits[(p, name, f)], rng_pat) <= 0.001 for f in range(5))
for p in pats} for name in PAT_LABELS})
print(pt.join(pperm).round(3).to_string())
for name in PAT_LABELS:
pooled = lr.loc[name, "kappa_mean"]
print(f"{name}: pooled {pooled:.3f}; best patient {pt[name].idxmax()} at {pt[name].max():.3f}; "
f"patients at or above pooled: {int((pt[name] >= pooled).sum())}; median patient {pt[name].median():.3f}; "
f"patients above chance (p <= 0.001) in all 5 folds: {int((pperm[f'{name} folds p<=0.001'] == 5).sum())}")
print("\nClock check per patient: Summer with true labels against labels shifted by 300 s:")
print(f"patients whose true kappa exceeds the shifted one: {int((pt['Summer'] > pt['Summer shifted 300 s']).sum())} of {len(pats)}; "
f"shifted kappa from {pt['Summer shifted 300 s'].min():.3f} to {pt['Summer shifted 300 s'].max():.3f} "
f"(median {pt['Summer shifted 300 s'].median():.3f}); true from {pt['Summer'].min():.3f} to {pt['Summer'].max():.3f}")
580 per-patient fits in 105 s
units Summer Scene Cuts Indoor/Outdoor Summer shifted 300 s Summer folds p<=0.001 Scene Cuts folds p<=0.001 Indoor/Outdoor folds p<=0.001 1 40 -0.014 -0.002 0.161 -0.024 3 0 5 2 83 0.004 0.034 0.256 -0.033 3 5 5 3 51 -0.017 0.012 0.191 0.024 2 2 5 4 120 -0.028 0.006 0.170 -0.022 2 1 5 5 79 0.061 0.010 0.192 -0.015 2 3 5 6 99 0.005 0.003 0.075 0.082 2 1 3 7 84 0.268 0.007 0.061 0.014 5 2 3 8 45 -0.004 -0.001 0.073 -0.017 2 0 4 9 41 0.082 -0.000 0.058 0.100 5 0 4 10 108 0.227 0.025 0.037 -0.086 5 5 2 11 31 -0.095 0.001 0.137 0.007 1 2 5 12 52 -0.041 0.016 0.099 -0.004 1 4 4 13 137 0.130 0.019 0.210 -0.041 4 3 5 14 30 -0.003 0.009 0.094 0.044 1 2 3 15 97 0.197 0.013 0.072 -0.025 5 3 4 16 67 -0.031 0.011 0.008 -0.105 0 2 2 17 101 0.082 0.005 0.055 -0.010 3 1 3 18 40 0.100 0.009 0.050 -0.096 4 2 4 19 105 0.060 0.036 0.209 0.024 2 5 5 20 88 -0.030 0.043 0.217 -0.025 1 5 5 21 109 0.113 0.005 0.153 -0.045 4 2 4 22 79 0.241 0.026 0.157 0.067 5 5 4 23 50 0.101 0.007 0.062 -0.016 4 1 3 24 88 0.077 0.023 0.109 -0.047 4 4 5 25 110 -0.011 0.011 0.201 0.003 2 1 5 26 85 0.028 0.009 0.072 0.031 2 2 3 27 90 0.013 0.005 0.005 0.099 2 1 2 28 111 0.064 0.023 0.086 0.054 2 4 3 29 66 0.127 0.011 0.139 0.045 4 2 5 Summer: pooled 0.301; best patient 7 at 0.268; patients at or above pooled: 0; median patient 0.060; patients above chance (p <= 0.001) in all 5 folds: 5 Scene Cuts: pooled 0.077; best patient 20 at 0.043; patients at or above pooled: 0; median patient 0.010; patients above chance (p <= 0.001) in all 5 folds: 5 Indoor/Outdoor: pooled 0.359; best patient 2 at 0.256; patients at or above pooled: 0; median patient 0.099; patients above chance (p <= 0.001) in all 5 folds: 12 Clock check per patient: Summer with true labels against labels shifted by 300 s: patients whose true kappa exceeds the shifted one: 18 of 29; shifted kappa from -0.105 to 0.100 (median -0.010); true from -0.095 to 0.268
Result: reproduced. No patient reaches the pooled kappa for any of the three labels. For Summer the best patient is 7 at 0.268 against 0.301 pooled (median patient 0.060); for Scene Cuts patient 20 at 0.043 against 0.077 (median 0.010); for Indoor/Outdoor patient 2 at 0.256 against 0.359 (median 0.099). By the test-label shuffle, 5 patients are above chance (p = 0.001) in all five folds for Summer, 5 for Scene Cuts and 12 for Indoor/Outdoor. With the Summer labels shifted by 300 s, per-patient kappas range from -0.105 to 0.100 (median -0.010).
fig4-named-patients¶
For Summer, patients 7, 10, 15 and 22 exceed kappa 0.2; for Indoor/Outdoor, patients 2, 13, 20 and 25 exceed 0.2.
named = {"Summer": [7, 10, 15, 22], "Indoor/Outdoor": [2, 13, 20, 25], "Scene Cuts": [2, 10, 19, 20]}
hits = 0
for name, ps in named.items():
vals = pt.loc[ps, name]
above = pt.index[pt[name] > 0.2].tolist()
if name != "Scene Cuts":
hits += int((vals > 0.2).sum())
print(f"{name}: named patients {dict(zip(ps, vals.round(3)))}; all patients above 0.2: {above}; "
f"top 4 here: {pt[name].nlargest(4).index.tolist()}")
print("named Summer and Indoor/Outdoor pairs above 0.2:", hits, "of 8")
print("(Scene Cuts: the paper says only that these patients show slightly better results)")
Summer: named patients {7: 0.268, 10: 0.227, 15: 0.197, 22: 0.241}; all patients above 0.2: [7, 10, 22]; top 4 here: [7, 22, 10, 15]
Indoor/Outdoor: named patients {2: 0.256, 13: 0.21, 20: 0.217, 25: 0.201}; all patients above 0.2: [2, 13, 19, 20, 25]; top 4 here: [2, 20, 13, 19]
Scene Cuts: named patients {2: 0.034, 10: 0.025, 19: 0.036, 20: 0.043}; all patients above 0.2: []; top 4 here: [20, 19, 2, 22]
named Summer and Indoor/Outdoor pairs above 0.2: 7 of 8
(Scene Cuts: the paper says only that these patients show slightly better results)
Result: reproduced. Taking subject ids as the paper's patient numbers, 7 of the 8 named pairs exceed 0.2: for Summer patients 7 (0.268), 10 (0.227) and 22 (0.241), with 15 at 0.197; for Indoor/Outdoor patients 2 (0.256), 13 (0.210), 20 (0.217) and 25 (0.201). The four best patients here are 7, 22, 10 and 15 for Summer and 2, 20, 13 and 19 for Indoor/Outdoor. The four Summer patients are above chance in all five folds. For Scene Cuts, three of the patients the paper names (20, 19, 2) are among the four best here; patient 10 (0.025) is fifth, behind patient 22. That the named patients come out on top supports the identity of the numbering, without proving it.
fig5-region-decoding¶
Decoding from single regions depends on the feature; the amygdala is best and the hippocampus worst for Summer, and the parahippocampal cortex is best and the amygdala worst for Indoor/Outdoor.
started = time.time()
REG_LABELS = ("Summer", "Indoor/Outdoor")
REGS = ("A", "H", "EC", "PHC")
out = summer.run_jobs([(summer.fit_folds, (S, units.index[units.region == region].to_numpy(), {k: Y[k]}, folds))
for region in REGS for k in REG_LABELS])
rfits = {(region, name, f): r for (region, k), o in zip([(g, k) for g in REGS for k in REG_LABELS], out) for (name, f), r in o.items()}
print(f"{len(rfits)} region fits in {time.time() - started:.0f} s")
rt = pd.DataFrame({name: {region: np.mean([rfits[(region, name, f)]["kappa"] for f in range(5)]) for region in ("A", "H", "EC", "PHC")}
for name in REG_LABELS})
rs = pd.DataFrame({name: {region: np.std([rfits[(region, name, f)]["kappa"] for f in range(5)], ddof=1) / np.sqrt(5)
for region in ("A", "H", "EC", "PHC")} for name in REG_LABELS})
print("mean kappa:")
print(rt.round(3).to_string())
print("SEM:")
print(rs.round(3).to_string())
rng_reg = np.random.default_rng(34)
print("test-shuffle p per fold:")
for region in REGS:
for name in REG_LABELS:
print(f" {region} {name}: {[round(summer.permutation_p(rfits[(region, name, f)], rng_reg), 4) for f in range(5)]}")
print("units:", units.region.value_counts()[["A", "H", "EC", "PHC"]].to_dict())
print("paper (LSTM): Summer A 0.285, H 0.091, EC 0.191, PHC 0.171 (dataset code); Indoor/Outdoor PHC 0.30, H 0.26, EC 0.15, A 0.10 (Results)")
cond = {"Summer, A highest": rt.Summer.idxmax() == "A", "Summer, H lowest": rt.Summer.idxmin() == "H",
"Indoor/Outdoor, PHC highest": rt["Indoor/Outdoor"].idxmax() == "PHC", "Indoor/Outdoor, A lowest": rt["Indoor/Outdoor"].idxmin() == "A"}
print(cond, "; holding:", sum(cond.values()), "of 4")
40 region fits in 53 s
mean kappa:
Summer Indoor/Outdoor
A 0.291 0.122
H 0.068 0.302
EC 0.149 0.117
PHC 0.096 0.357
SEM:
Summer Indoor/Outdoor
A 0.056 0.041
H 0.041 0.040
EC 0.060 0.052
PHC 0.045 0.043
test-shuffle p per fold:
A Summer: [0.001, 0.001, 0.001, 0.001, 0.001]
A Indoor/Outdoor: [0.001, 0.001, 0.001, 0.001, 0.9011]
H Summer: [0.001, 0.0569, 1.0, 0.001, 0.001]
H Indoor/Outdoor: [0.001, 0.001, 0.001, 0.001, 0.001]
EC Summer: [0.001, 0.003, 0.001, 0.9261, 0.001]
EC Indoor/Outdoor: [0.001, 0.001, 0.001, 0.001, 1.0]
PHC Summer: [0.001, 0.994, 0.001, 0.001, 0.001]
PHC Indoor/Outdoor: [0.001, 0.001, 0.001, 0.001, 0.001]
units: {'A': 559, 'H': 846, 'EC': 493, 'PHC': 293}
paper (LSTM): Summer A 0.285, H 0.091, EC 0.191, PHC 0.171 (dataset code); Indoor/Outdoor PHC 0.30, H 0.26, EC 0.15, A 0.10 (Results)
{'Summer, A highest': True, 'Summer, H lowest': True, 'Indoor/Outdoor, PHC highest': True, 'Indoor/Outdoor, A lowest': False} ; holding: 3 of 4
Result: partially reproduced. With the logistic regression, Summer is decoded best from the amygdala (0.291) and worst from the hippocampus (0.068; EC 0.149, PHC 0.096), and Indoor/Outdoor best from PHC (0.357), then H (0.302). The amygdala is not the lowest for Indoor/Outdoor: A (0.122) and EC (0.117) differ by less than their SEMs (0.041, 0.052). Three of the four orderings hold. Not every region's fold is above chance: by the test-label shuffle, one fold each for A Indoor/Outdoor, EC Summer, EC Indoor/Outdoor and PHC Summer, and two for H Summer, have p above 0.05. The regions here hold different units from the paper's (PHC 293 against 373).
fig5-faces-presence¶
Face labels decode slightly better than character labels (Summer 0.33, Tom 0.27), and Summer Presence (0.27) is decoded best from the amygdala (0.26).
Tested with the logistic regression in place of the LSTM. The dandiset's face labels hold only 0 and 1, so the frames the Methods exclude (a face neither shown nor clearly absent) are negatives here.
rng_face = np.random.default_rng(38)
FACE = {"Summer Faces": "summer-faces", "Tom Faces": "tom-faces", "Summer Presence": "summer-presence"}
for name, label in FACE.items():
Y[name] = summer.bin_labels(sessions[0], label)
print(pd.DataFrame({name: {"positive": float(np.mean(Y[name] == 1)), "negative": float(np.mean(Y[name] == 0))}
for name in ("Summer", "Summer Faces", "Summer Presence", "Tom", "Tom Faces")}).round(3).to_string())
started = time.time()
keys = [(name, f) for name in FACE for f in range(5)]
ffits = dict(zip(keys, summer.run_jobs([(summer.fit_logreg, (S, Y[name], *folds[f])) for name, f in keys])))
out = summer.run_jobs([(summer.fit_folds, (S, units.index[units.region == region].to_numpy(), {k: Y[k]}, folds))
for region in REGS for k in FACE])
frfits = {(region, name, f): r for (region, k), o in zip([(g, k) for g in REGS for k in FACE], out) for (name, f), r in o.items()}
print(f"{len(ffits) + len(frfits)} fits in {time.time() - started:.0f} s")
rows = []
for name in ("Summer", "Summer Faces", "Summer Presence", "Tom", "Tom Faces"):
fl = [fits[(name, f)] if name in DEC else ffits[(name, f)] for f in range(5)]
ks = np.array([r["kappa"] for r in fl])
row = {"label": name, "pooled": ks.mean(), "pooled_sem": ks.std(ddof=1) / np.sqrt(5), "folds": np.round(ks, 3).tolist(),
"max_perm_p": max(summer.permutation_p(r, rng_face) for r in fl)}
for region in REGS:
rk = [(rfits if (region, name, 0) in rfits else frfits).get((region, name, f), {}).get("kappa", np.nan) for f in range(5)]
row[region] = np.mean(rk)
rows.append(row)
ft = pd.DataFrame(rows).set_index("label")
print(ft.round(4).to_string())
print("paper (LSTM): Summer 0.31, Summer Faces 0.33, Tom 0.23, Tom Faces 0.27, Summer Presence 0.27 (amygdala 0.26)")
cond = {"a: Summer Faces above Summer": ft.loc["Summer Faces", "pooled"] > ft.loc["Summer", "pooled"],
"b: Tom Faces above Tom": ft.loc["Tom Faces", "pooled"] > ft.loc["Tom", "pooled"],
"c: Summer Presence highest in A": ft.loc["Summer Presence", list(REGS)].astype(float).idxmax() == "A"}
print(cond, "; holding:", sum(bool(v) for v in cond.values()), "of 3")
Summer Summer Faces Summer Presence Tom Tom Faces positive 0.527 0.298 0.643 0.792 0.448 negative 0.473 0.702 0.357 0.208 0.552
75 fits in 152 s
pooled pooled_sem folds max_perm_p A H EC PHC
label
Summer 0.3011 0.0560 [0.459, 0.299, 0.316, 0.108, 0.324] 0.0010 0.2909 0.0678 0.1487 0.0965
Summer Faces 0.2166 0.0534 [0.335, 0.175, 0.331, 0.051, 0.191] 0.0010 0.1926 0.0659 0.1298 0.1215
Summer Presence 0.2343 0.0732 [0.504, 0.059, 0.209, 0.216, 0.184] 0.0010 0.2917 0.1026 0.1097 0.1323
Tom 0.2184 0.0608 [0.275, 0.31, 0.334, 0.002, 0.172] 0.4396 NaN NaN NaN NaN
Tom Faces 0.1666 0.0373 [0.136, 0.181, 0.305, 0.108, 0.102] 0.0010 0.1209 0.0697 0.1046 0.0837
paper (LSTM): Summer 0.31, Summer Faces 0.33, Tom 0.23, Tom Faces 0.27, Summer Presence 0.27 (amygdala 0.26)
{'a: Summer Faces above Summer': np.False_, 'b: Tom Faces above Tom': np.False_, 'c: Summer Presence highest in A': True} ; holding: 1 of 3
Result: partially reproduced. One of the three conditions holds. Summer Presence is decoded best from the amygdala (0.292, against PHC 0.132, EC 0.110 and H 0.103; pooled 0.234). The face labels decode worse than the character labels here, not better: Summer Faces 0.217 against Summer 0.301, Tom Faces 0.167 against Tom 0.218. Every fold of the three new labels is above chance (p = 0.001). The face labels here keep as negatives the frames the Methods exclude, which may lower their kappa; that cannot be checked, as the dandiset does not mark those frames.
fig6-responsive-subsets¶
Responsive neurons drive decoding of camera cuts but not of Summer; without them, Summer is decoded about as well as from all neurons while camera cuts drop the most, and responsive neurons alone decode camera cuts about as well as all neurons but Summer less well.
rng_sub = np.random.default_rng(36)
started = time.time()
keys, jobs = [], []
sets = {}
for name in ("Summer", "Camera Cuts"):
r = (P[name] <= 0.001).to_numpy()
sets[(name, "responsive only")] = np.where(r)[0]
sets[(name, "non-responsive only")] = np.where(~r)[0]
# Variant: non-responsive units that spike in at least a third of the pseudo-trials (units
# below that get p = 1 without a test).
act = (resp[(resp.rule == "methods") & (resp.label == name)].set_index("unit").active.loc[units.index] >= 1 / 3).to_numpy()
sets[(name, "non-responsive, active only")] = np.where(~r & act)[0]
print(f"{name}: responsive {int(r.sum())}; non-responsive {int((~r).sum())}, of which active {int((~r & act).sum())} "
f"and inactive (p = 1 untested) {int((~r & ~act).sum())}")
if r.sum() <= (~r).sum():
for d in range(3):
sets[(name, f"non-responsive, size-matched {d + 1}")] = np.sort(rng_sub.choice(np.where(~r)[0], r.sum(), replace=False))
else:
print(f"{name}: {int(r.sum())} responsive and {int((~r).sum())} non-responsive units; no size-matched non-responsive set")
out = summer.run_jobs([(summer.fit_folds, (S, cols, {name: Y[name]}, folds)) for (name, kind), cols in sets.items()])
sfits = {(name, kind, f): o[(name, f)] for ((name, kind), cols), o in zip(sets.items(), out) for f in range(5)}
print(f"{len(sfits)} subset fits in {time.time() - started:.0f} s")
rows = []
for name in ("Summer", "Camera Cuts"):
rows.append({"label": name, "set": "complete", "units": S.shape[1], "kappa": lr.loc[name, "kappa_mean"], "sem": lr.loc[name, "kappa_sem"]})
for (nm, kind), cols in sets.items():
if nm == name:
ks = [sfits[(name, kind, f)]["kappa"] for f in range(5)]
rows.append({"label": name, "set": kind, "units": len(cols), "kappa": np.mean(ks), "sem": np.std(ks, ddof=1) / np.sqrt(5)})
sub = pd.DataFrame(rows)
print(sub.round(3).to_string(index=False))
k = sub.set_index(["label", "set"]).kappa
drop_s = k["Summer", "complete"] - k["Summer", "non-responsive only"]
drop_c = k["Camera Cuts", "complete"] - k["Camera Cuts", "non-responsive only"]
cond = {"a: Summer without responsive >= complete - 0.05": k["Summer", "non-responsive only"] >= k["Summer", "complete"] - 0.05,
"b: Summer responsive only < complete - 0.05": k["Summer", "responsive only"] < k["Summer", "complete"] - 0.05,
"c: drop without responsive larger for camera cuts": drop_c > drop_s,
"d: camera cuts responsive only >= complete - 0.05": k["Camera Cuts", "responsive only"] >= k["Camera Cuts", "complete"] - 0.05}
print(f"drop without responsive units: Summer {drop_s:.3f}, Camera Cuts {drop_c:.3f}")
print(f"variant, drop to non-responsive active units: Summer {k['Summer', 'complete'] - k['Summer', 'non-responsive, active only']:.3f}, "
f"Camera Cuts {k['Camera Cuts', 'complete'] - k['Camera Cuts', 'non-responsive, active only']:.3f}")
print(cond, "; holding:", sum(bool(v) for v in cond.values()), "of 4")
Summer: responsive 1013; non-responsive 1273, of which active 730 and inactive (p = 1 untested) 543 Camera Cuts: responsive 1616; non-responsive 670, of which active 122 and inactive (p = 1 untested) 548 Camera Cuts: 1616 responsive and 670 non-responsive units; no size-matched non-responsive set
45 subset fits in 117 s
label set units kappa sem
Summer complete 2286 0.301 0.056
Summer responsive only 1013 0.256 0.061
Summer non-responsive only 1273 0.294 0.062
Summer non-responsive, active only 730 0.241 0.041
Summer non-responsive, size-matched 1 1013 0.278 0.042
Summer non-responsive, size-matched 2 1013 0.262 0.045
Summer non-responsive, size-matched 3 1013 0.300 0.069
Camera Cuts complete 2286 0.096 0.020
Camera Cuts responsive only 1616 0.115 0.015
Camera Cuts non-responsive only 670 0.013 0.009
Camera Cuts non-responsive, active only 122 0.034 0.005
drop without responsive units: Summer 0.007, Camera Cuts 0.083
variant, drop to non-responsive active units: Summer 0.060, Camera Cuts 0.062
{'a: Summer without responsive >= complete - 0.05': np.True_, 'b: Summer responsive only < complete - 0.05': np.False_, 'c: drop without responsive larger for camera cuts': np.True_, 'd: camera cuts responsive only >= complete - 0.05': np.True_} ; holding: 3 of 4
Result: blocked (methods_underspecified), as the responsive sets come from the miscalibrated test (1013 of 2286 units for Summer and 1616 for Camera Cuts at alpha 0.001). With these sets, three of the four conditions of the criterion hold: for Summer, 0.294 without responsive units against 0.301 for all units; for Camera Cuts, 0.115 with responsive units only against 0.096 for all, and 0.013 without them. Condition (b) does not hold: Summer with responsive units only gives 0.256, above the complete value less 0.05 (0.251). Size-matched non-responsive sets give 0.262 to 0.300 for Summer; none can be drawn for Camera Cuts, which has 670 non-responsive units. Of these 670, 548 spike in fewer than a third of the pseudo-trials and get p = 1 untested. Keeping only the non-responsive units that spike, Camera Cuts give 0.034 (122 units) and Summer 0.241 (730 units), drops of 0.062 and 0.060 from all units, so with active units alone the drop is about the same for camera cuts and for Summer.
fig7a-500-neurons and fig7b-top500-overlap¶
Both are in the companion notebook top_ranked_units.ipynb, which fits the 20 sub-splits, ranks the
units and decodes Summer from the top-ranked units.
Blocked claims¶
fig3c-lstm-decoding and fig3d-temporal-shuffle rest on the paper's LSTM. Its code is not in the
paper's materials, the dataset's GitHub repository or its Zenodo archive (which hold a two-layer
linear baseline), and the Methods say that the number of linear layers, hidden sizes, batch size,
learning rate, weight decay and dropout were chosen for each label without giving them. These claims
are blocked with methods_underspecified. Training an LSTM for 700 epochs on 2286 inputs for each
label and fold would also not fit the time limit on a CPU; this was not tried. fig3d shuffles the
order of the 20 time bins of each input, which the logistic regression on window sums cannot see,
so no stand-in tests it. fig5-faces-presence and fig7a-500-neurons, also LSTM results in the paper,
are tested above and in top_ranked_units.ipynb with the logistic regression in their place.
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 |
|---|---|
| dataset-composition | partially reproduced |
| fig2b-phc-camera-cuts | blocked (methods underspecified) |
| tables1-camera-cuts-regions | blocked (methods underspecified) |
| fig2b-phc-scene-cuts | blocked (methods underspecified) |
| tables1-few-semantic-responses | blocked (methods underspecified) |
| tables1-no-indoor-response | blocked (methods underspecified) |
| results-cluster-regions | blocked (methods underspecified) |
| fig3b-single-neuron-threshold | reproduced |
| fig3c-logreg-characters-location | reproduced |
| fig3c-logreg-cuts-lower | reproduced |
| fig3e-temporal-gap | reproduced |
| fig4-patients-below-pooled | reproduced |
| fig4-named-patients | reproduced |
| fig5-region-decoding | partially reproduced |
| fig5-faces-presence | partially reproduced |
| fig6-responsive-subsets | blocked (methods underspecified) |
| fig7b-top500-overlap | partially reproduced (top_ranked_units.ipynb) |
| fig7a-500-neurons | partially reproduced (top_ranked_units.ipynb) |
| fig3c-lstm-decoding | blocked (methods underspecified) |
| fig3d-temporal-shuffle | blocked (methods underspecified) |
How the analysis differs from the paper's:
- The paper's LSTM is not run. Its logistic regression, which the paper reports for the patient-wise analysis and Figure 3c, is run as the Methods describe, with class weights in place of oversampling and C chosen from four values. It stands in for the LSTM in fig3e, fig5 (regions, faces and presence), fig6 and fig7a.
- The responsive-unit test follows the Methods as written. On this data it is not calibrated (68.4% of units responsive at random times at alpha 0.001 with as many events as camera cuts, 9.5 to 55.5% with as many as each other label's onsets), so the claims that rest on it are blocked. The dataset code's ZETA test is conservative at random times and gives counts far below Table S1 (camera cuts 79 units against 1025, Summer 11 against 46), so neither test as run here gives the paper's counts.
- Onsets of a label include its first segment, and the decoding windows at the ends of the analysis portion count the spikes just outside it; both were changed after the first verification round.
- Regions are the dandiset's
brain_region, grouped as in the dataset's code; their counts differ from the paper's. - Cut labels are extended to 520 ms after the cut for scene cuts too; the Methods state it for camera cuts.
What was found about the dandiset:
- It holds the 2286 units of the 29 patients, the movie annotations and the watchlogs needed for the analyses, and the logistic regression results of the paper are recovered from it.
- The units' regions give counts that differ from the paper's (PHC 293 against 373, EC 493 against 440, PIC 68 against 99, 27 units in other regions), with no mapping between the two.
- Spike times are stored in milliseconds; the patient-aligned annotation columns are described as seconds but hold milliseconds.
- The face labels lack the excluded state the Methods use; the face decoding here keeps those frames as negatives.
- The earlier published version 0.260519.1133 holds 28 of the recordings twice.
- The dataset's code identifies responsive units with the ZETA test, not the test of the Methods, and its training code selects no hippocampal unit.
print(f"total run time {(time.time() - started_all) / 60:.1f} min")
total run time 19.6 min