Reproducing Zheng et al. (2022) from dandiset 000207¶
Paper: Zheng J, Schjetnan AGP, Yebra M, Gomes BA, Mosher CP, Kalia SK, Valiante TA, Mamelak AN, Kreiman G, Rutishauser U (2022). Neurons detect cognitive boundaries to structure episodic memories in humans. Nature Neuroscience 25:358-368. doi:10.1038/s41593-022-01020-w. The journal version is behind a paywall; the full text of the author manuscript in PubMed Central (PMC8966433) was used, with its Methods and figure captions. The supplementary tables and figures were not available.
Data: DANDI:000207, version 0.220216.0323, streamed from the archive. This is the version the paper's Data Availability statement cites. It is not the most recent one; the data check below says why it was chosen.
Authors' code: the dandiset links to the laboratory's repository, rutishauserlab/cogboundary-zheng (commit 60f7ece, 2022-02-15). It was read after a first version of this notebook. It holds four scripts that plot rasters of single cells from the NWB files. It has no code that selects boundary or event cells, computes a chance level, runs a statistical test, or links a tested frame to its clip. One thing in it bears on this analysis: its boundary-aligned plot takes as soft-boundary trials every soft boundary of a session, in SB clips and in HB clips, 75 in all. That set is used below for the cell definitions, and the earlier choice (the first soft boundary of each SB clip) is kept as an alternative reading.
When the criteria were written. The claims file was written after the data check, which
already shows the behavioural means, and after a first count of boundary and event cells. The
criteria of fig5d-state-shift, fig3dh-no-onset-response and fig3j-response-latency, and
all tolerances, were written with those first numbers in view. No criterion or tolerance was
changed after the authors' code was read.
What can and cannot be tested¶
Patients watched 90 silent video clips of about 8 s: 30 with no boundary (NB), 30 with one to three soft boundaries (SB, a cut to another scene of the same film) and 30 with a hard boundary (HB, a cut to a different film, 4 s into the clip). Their memory was then tested with scene recognition (is this frame old or new?) and time discrimination (which of two frames came first?).
The dandiset holds, for 19 patients, the spike times of sorted units in the hippocampus, amygdala and parahippocampal gyrus, and the trial tables of the three tasks. It holds no local field potentials, no units outside the medial temporal lobe, no labels saying which units the paper classified as boundary or event cells, and no stimuli.
| Claim | Paper | Here |
|---|---|---|
fig2a-order-accuracy |
Order memory is at chance across a hard boundary (53%), 73% otherwise | Tested |
fig2b-order-reaction-time |
Order judgements are slower across a hard boundary | Tested |
fig2c-order-confidence |
Order judgements are less confident across a hard boundary | Tested |
edfig2a-recognition-similar |
Scene recognition accuracy is the same for the three kinds of clip | Tested |
fig2def-recognition-vs-boundary-distance |
Recognition falls with time since the last boundary | Blocked: the position of a tested frame in its clip is not in the dandiset |
fig3c-boundary-cells |
42 of 580 neurons (7.2%) are boundary cells | Analysed, then blocked: the procedure behind the paper's counts could not be determined |
fig3g-event-cells |
36 of 580 neurons (6.2%) are event cells | Analysed, then blocked: the procedure behind the paper's counts could not be determined |
fig3-cell-regions |
Boundary cells are commonest in the parahippocampal gyrus, event cells in the hippocampus | Analysed, then blocked: the procedure behind the paper's counts could not be determined |
fig3dh-no-onset-response |
Boundary and event cells do not respond to clip or image onsets | Blocked: the cells are not labelled, and too few are found here to test it |
edfig5-clip-onset-cells |
13.1% of neurons respond to clip onsets, 0.7% to clip offsets | Tested |
fig3j-response-latency |
Boundary cells respond 100 ms before event cells | Blocked: as fig3dh-no-onset-response |
fig3-outside-mtl |
Few boundary or event cells outside the medial temporal lobe | Blocked: those neurons are not in the dandiset |
fig4c-boundary-cells-predict-recognition |
Boundary cells fire more before frames that are later recognised | Blocked: test trials with SB or HB clips cannot be linked to their clips |
fig4h-event-cells-phase-predicts-order |
Theta phase locking of event cells predicts order memory | Blocked: no local field potentials |
fig5d-state-shift |
The population state jumps after soft and hard boundaries | Tested |
fig5ij-state-shift-predicts-memory |
The size of the jump predicts recognition and order memory | Blocked: test trials cannot be linked to clips |
fig6-reinstatement, fig6i-reinstatement-timing, fig6jm-reinstatement-confidence |
The state after a boundary is reinstated at recognition | Blocked: test trials cannot be linked to clips |
fig6no-false-alarm-similarity |
False alarms come with reinstatement when the foil resembles a target | Blocked: no similarity ratings, and test trials cannot be linked to clips |
Nine claims were analysed. Six have a tested outcome. The other three, and the eleven that could not be analysed, are blocked.
import dataclasses
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import cogboundary as cb
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
plt.rcParams.update({"figure.dpi": 100, "axes.spines.top": False, "axes.spines.right": False})
COLOR = {"NB": "#2e8b57", "SB": "#1f77b4", "HB": "#c0392b"}
rng = np.random.default_rng(207)
Load the data and check it against its description¶
cogboundary.load_sessions streams the 19 files. Each gives the spike times of its units, the
region of each unit's electrode, and the three trial tables.
sessions = cb.load_sessions()
rows = []
for s in sessions:
rows.append({
"subject": s.subject,
"units": len(s.units),
**{r: int((s.region == r).sum()) for r in cb.REGIONS},
"spikes": s.sanity["spikes"],
"first spike (s)": f"{s.sanity['first spike (s)']:.4g}",
"last spike (s)": f"{s.sanity['last spike (s)']:.6g}",
"last event (s)": f"{s.sanity['last event (s)']:.6g}",
"clips NB/SB/HB": "/".join(str(int((s.encoding.kind == k).sum())) for k in cb.KINDS),
"recognition trials": len(s.recognition),
"order trials": len(s.order),
"spike clock ok": s.aligned,
})
overview = pd.DataFrame(rows).set_index("subject")
print(overview.to_string())
good = [s for s in sessions if s.aligned]
region_all = np.concatenate([s.region for s in sessions])
region_good = np.concatenate([s.region for s in good])
print(f"\nFiles: {len(sessions)}. Units: {len(region_all)} in all files, "
f"{len(region_good)} in the {len(good)} files whose spike clock matches the trial tables.")
print("Units by region, all files: ", {r: int((region_all == r).sum()) for r in cb.REGIONS})
print("Units by region, usable files:", {r: int((region_good == r).sum()) for r in cb.REGIONS})
print("Paper: 580 medial temporal lobe units in 19 subjects (hippocampus 343, amygdala 169, "
"parahippocampal gyrus 68), and 405 units elsewhere.")
print("Regions named in the files:", sorted({f"{r}_{h}" for s in sessions for r, h in zip(s.region, s.side)}))
print("Spike time unit stated in the files:", {s.sanity["spike time unit"] for s in sessions})
units hippocampus amygdala parahippocampus spikes first spike (s) last spike (s) last event (s) clips NB/SB/HB recognition trials order trials spike clock ok
subject
1 43 12 31 0 603265 0.001647 3088.04 3088.05 30/30/30 180 180 True
2 45 21 14 10 312208 0.001498 2808.29 2808.32 30/30/30 180 180 True
3 31 15 14 2 219145 0.01671 2944.46 2944.5 30/30/30 180 180 True
4 35 25 10 0 214278 0.007858 2821.07 2821.07 30/30/30 180 180 True
5 30 16 14 0 180037 0.01916 2954.55 2954.55 30/30/30 180 180 True
6 19 6 13 0 132389 0.02841 3131.44 3131.45 30/30/30 180 180 True
7 34 20 14 0 260304 0.001017 2968.46 2968.46 30/30/30 180 180 True
8 33 12 7 14 326706 0.01813 2912.4 2912.41 30/30/30 180 180 True
9 31 16 15 0 209769 0.01277 2847.93 2847.93 30/30/30 180 180 True
10 27 27 0 0 140083 0.01169 2884.13 2884.14 30/30/30 180 180 True
11 14 6 5 3 128079 0.1075 3012.87 3012.9 30/30/30 180 180 True
12 28 11 17 0 266463 0.004111 2844.64 2844.64 30/30/30 180 180 True
13 48 15 33 0 229005 0.004601 2697.15 2697.17 30/30/30 180 180 True
14 30 18 5 7 323033 0.01343 2675.31 2675.32 30/30/30 180 180 True
15 31 14 11 6 569548 0.0006405 3073.66 3073.67 30/30/30 180 180 True
16 14 6 5 3 292263 0.01015 3101.95 3102.46 30/30/30 180 180 True
17 27 6 15 6 490333 0.01714 2752.97 2752.98 30/30/30 180 180 True
18 29 21 6 2 392084 0.003138 2846.18 2846.18 30/30/30 180 180 True
19 16 12 4 0 176945 -7.02e+12 -7.02021e+12 2892.19 30/30/30 180 180 False
Files: 19. Units: 565 in all files, 549 in the 18 files whose spike clock matches the trial tables.
Units by region, all files: {'hippocampus': 279, 'amygdala': 233, 'parahippocampus': 53}
Units by region, usable files: {'hippocampus': 267, 'amygdala': 229, 'parahippocampus': 53}
Paper: 580 medial temporal lobe units in 19 subjects (hippocampus 343, amygdala 169, parahippocampal gyrus 68), and 405 units elsewhere.
Regions named in the files: ['amygdala_left', 'amygdala_right', 'hippocampus_left', 'hippocampus_right', 'parahippocampus_left', 'parahippocampus_right']
Spike time unit stated in the files: {'not stated'}
enc = pd.concat([s.encoding.assign(subject=s.subject) for s in sessions], ignore_index=True)
rec = pd.concat([s.recognition.assign(subject=s.subject) for s in sessions], ignore_index=True)
order = pd.concat([s.order.assign(subject=s.subject) for s in sessions], ignore_index=True)
print("ENCODING")
print("Clip duration (s): min %.2f, median %.2f, max %.2f" % tuple(
(enc.stop_time - enc.start_time).quantile([0, 0.5, 1])))
print("Fixation before a clip (s): min %.3f, max %.3f" % tuple(
(enc.start_time - enc.fixcross_time).quantile([0, 1])))
print("stimCategory against the condition in the clip's name:")
print(pd.crosstab(enc.kind, enc.stimCategory).reindex(cb.KINDS).to_string())
print("Its description says:", sessions[0].sanity["descriptions"]["encoding stimCategory"])
offsets = enc.groupby("kind")[["b1", "b2", "b3"]].agg(["count", "min", "max"]).reindex(cb.KINDS).round(3)
print("\nBoundary times relative to clip onset (s); count is over 19 subjects x 30 clips:")
print(offsets.to_string())
same = enc.groupby("clip")[["b1", "b2", "b3"]].std().max().max()
print(f"Largest standard deviation of a clip's boundary offset across subjects: {same:.1e} s")
raw_equal_start = int(((enc.boundary2_time - enc.start_time).abs() < 1e-6).sum())
print(f"Clips whose boundary2_time equals their start_time: {raw_equal_start} of {len(enc)} "
f"(NaN: {int(enc.boundary2_time.isna().sum())}).")
print("Its description says:", sessions[0].sanity["descriptions"]["boundary2_time"][:74])
print("Clips with the same name in every file:", enc.groupby("subject")["clip"].apply(frozenset).nunique() == 1,
"| same order in every file:", enc.groupby("subject")["clip"].apply(tuple).nunique() == 1)
durations = np.concatenate([[u[-1] - u[0] for u in s.units] for s in good])
rates = np.concatenate([[len(u) for u in s.units] for s in good]) / durations
print("\nUNITS (usable files)")
print("Firing rate (Hz): min %.2f, 25%% %.2f, median %.2f, 75%% %.2f, max %.1f" % tuple(
np.quantile(rates, [0, 0.25, 0.5, 0.75, 1])))
print("Units with unsorted spike times in the file:",
{s.subject: s.sanity["units with unsorted spike times"] for s in sessions
if s.sanity["units with unsorted spike times"]})
print("Events with code 0, which the description does not list:",
sum(s.sanity["events with code 0"] for s in sessions),
"| their earliest timestamp: %.4g s" % np.nanmin([s.sanity["earliest code-0 timestamp (s)"] for s in sessions]))
ENCODING
Clip duration (s): min 6.87, median 8.03, max 9.17
Fixation before a clip (s): min 0.900, max 1.100
stimCategory against the condition in the clip's name:
stimCategory 0 1 2
kind
NB 570 0 0
SB 0 570 0
HB 0 0 570
Its description says: 1=no boundary (NB), 2=soft boundary (SB), 3=hard boundary (HB). Note - HB clips can, in addition, also contain SB boundaries.
Boundary times relative to clip onset (s); count is over 19 subjects x 30 clips:
b1 b2 b3
count min max count min max count min max
kind
NB 570 4.000 4.000 0 NaN NaN 0 NaN NaN
SB 570 0.552 5.853 285 1.832 4.896 95 3.587 5.828
HB 570 4.000 4.000 285 1.552 6.082 190 5.869 6.280
Largest standard deviation of a clip's boundary offset across subjects: 1.2e-04 s
Clips whose boundary2_time equals their start_time: 1140 of 1710 (NaN: 0).
Its description says: time of second boundary (0 if it does not exist; absolute timestamp). For
Clips with the same name in every file: True | same order in every file: False
UNITS (usable files)
Firing rate (Hz): min 0.01, 25% 0.46, median 1.90, 75% 6.20, max 22.9
Units with unsorted spike times in the file: {'10': 1}
Events with code 0, which the description does not list: 43 | their earliest timestamp: -7.022e+12 s
print("SCENE RECOGNITION")
print("Trials by kind and target/foil, per subject:",
sorted(int(v) for v in rec.groupby(["subject", "kind", "target"]).size().unique()))
print("accuracy code against stimulus and response (old_new: 2 = old, 1 = new, 0 = no response):")
print(pd.crosstab([rec.target.map({True: "target", False: "foil"}), rec.old_new], rec.accuracy).to_string())
print("Description of the accuracy column:", sessions[0].sanity["descriptions"]["accuracy"])
print("resp_value of trials without a response:", sorted(int(v) for v in rec.resp_value[~rec.responded].unique()),
"| description:", sessions[0].sanity["descriptions"]["resp_value"])
print("\nTIME DISCRIMINATION")
print("Trials by kind, per subject:", sorted(int(v) for v in order.groupby(["subject", "kind"]).size().unique()))
print("accuracy code against whether the response key equals the correct side, and whether there was a response:")
print(pd.crosstab([(order.key == order.leftright).rename("key is correct"), order.responded], order.accuracy).to_string())
gap = order.response_time - order.stop_time - order.RT
print("RT (s): min %.2f, median %.2f, max %.2f; negative in %d of %d trials" % (
order.RT.min(), order.RT.median(), order.RT.max(), (order.RT < 0).sum(), len(order)))
print("response_time - stop_time - RT (s): median %.4f, largest deviation %.4f, in trials with a response" % (
gap[order.responded].median(), gap[order.responded].abs().max()))
print("Frame shown for (s): recognition median %.2f, time discrimination median %.2f" % (
(rec.stop_time - rec.start_time).median(), (order.stop_time - order.start_time).median()))
print("Trials without a response: recognition %d of %d, time discrimination %d of %d" % (
(~rec.responded).sum(), len(rec), (~order.responded).sum(), len(order)))
print(" by subject, recognition:", rec[~rec.responded].groupby("subject").size().to_dict())
print(" by subject, time discrimination:", order[~order.responded].groupby("subject").size().to_dict())
print("\nLINK BETWEEN TEST TRIALS AND CLIPS")
clip_numbers = enc["clip"].str.extract(r"_(\d+)$")[0].astype(int)
for name, table in (("recognition", rec), ("time discrimination", order)):
number = table.frameName.str.extract(r"^[A-Z]+_(\d+)_")[0].astype(int)
for kind in cb.KINDS:
m = table.kind == kind
per_subject = table[m].assign(number=number[m]).groupby("subject").number
print(f"{name:>20}, {kind}: frame names carry {number[m].nunique()} distinct numbers "
f"({number[m].min()} to {number[m].max()}); trials per number in one subject: "
f"{sorted(int(v) for v in per_subject.value_counts().unique())}")
print("Encoding clips carry the numbers", clip_numbers.min(), "to", clip_numbers.max(), "for each kind.")
print("\nCOLUMNS OF THE TEST TABLES, as stored in the file (the analysis adds kind, target, correct, responded)")
print(" recognition_table:", [c for c in sessions[0].recognition.columns if c not in ("kind", "target", "correct", "responded")])
print(" timediscrimination_table:", [c for c in sessions[0].order.columns if c not in ("kind", "correct", "responded")])
print(" encoding_table:", [c for c in sessions[0].encoding.columns if c not in ("kind", "b1", "b2", "b3")])
print("Example frame names, recognition:", sorted(rec.frameName.unique())[:3], "...",
"| time discrimination:", sorted(order.frameName.unique())[:3], "...")
tree = cb.file_tree(cb.asset_urls()["1"])
print(f"\nEvery group and dataset in subject 1's file ({len(tree)} items):")
print("\n".join(" " + line for line in tree))
SCENE RECOGNITION
Trials by kind and target/foil, per subject: [30]
accuracy code against stimulus and response (old_new: 2 = old, 1 = new, 0 = no response):
accuracy 1 2
target old_new
foil 0 27 0
1 0 1281
2 402 0
target 0 30 0
1 459 0
2 0 1221
Description of the accuracy column: correct = 2, incorrect 1, no response = 0
resp_value of trials without a response: [1] | description: range = 2-7, -1=no response
TIME DISCRIMINATION
Trials by kind, per subject: [60]
accuracy code against whether the response key equals the correct side, and whether there was a response:
accuracy 1.0 2.0
key is correct responded
False False 23 0
True 1144 0
True False 0 20
True 0 2233
RT (s): min -3.63, median 1.76, max 6.86; negative in 151 of 3420 trials
response_time - stop_time - RT (s): median 0.0000, largest deviation 0.0009, in trials with a response
Frame shown for (s): recognition median 1.53, time discrimination median 2.03
Trials without a response: recognition 57 of 3420, time discrimination 43 of 3420
by subject, recognition: {'10': 6, '12': 26, '14': 1, '15': 2, '16': 3, '17': 1, '19': 1, '3': 1, '4': 4, '5': 1, '6': 5, '7': 1, '9': 5}
by subject, time discrimination: {'10': 3, '12': 15, '16': 2, '19': 1, '2': 4, '4': 2, '6': 5, '7': 5, '8': 1, '9': 5}
LINK BETWEEN TEST TRIALS AND CLIPS
recognition, NB: frame names carry 30 distinct numbers (1 to 30); trials per number in one subject: [2]
recognition, SB: frame names carry 10 distinct numbers (1 to 30); trials per number in one subject: [2, 4, 6, 10, 12]
recognition, HB: frame names carry 10 distinct numbers (1 to 30); trials per number in one subject: [2, 4, 6, 10, 12]
time discrimination, NB: frame names carry 30 distinct numbers (1 to 30); trials per number in one subject: [2]
time discrimination, SB: frame names carry 60 distinct numbers (1 to 60); trials per number in one subject: [1, 3]
time discrimination, HB: frame names carry 60 distinct numbers (1 to 60); trials per number in one subject: [1]
Encoding clips carry the numbers 1 to 30 for each kind.
COLUMNS OF THE TEST TABLES, as stored in the file (the analysis adds kind, target, correct, responded)
recognition_table: ['ExperimentID', 'RT', 'accuracy', 'boundary_type', 'confidence', 'fixcross_time', 'frameName', 'id', 'old_new', 'resp_value', 'response_time', 'start_time', 'stimuli_type', 'stop_time', 'trial_num']
timediscrimination_table: ['ExperimentID', 'RT', 'accuracy', 'boundary_type', 'confidence', 'fixcross_time', 'frameName', 'id', 'key', 'leftright', 'resp_value', 'response_time', 'start_time', 'stop_time', 'trial_num']
encoding_table: ['clip', 'ExperimentID', 'boundary1_time', 'boundary2_time', 'boundary3_time', 'fixcross_time', 'id', 'start_time', 'stimCategory', 'stop_time']
Example frame names, recognition: ['HB_1_frameA.png', 'HB_1_frameA_foil.png', 'HB_1_frameB.png'] ... | time discrimination: ['HB_10_LongA.png', 'HB_10_ShortA.png', 'HB_11_LongA.png'] ...
Every group and dataset in subject 1's file (145 items): acquisition/ acquisition/events/ acquisition/events/data (1720,) acquisition/events/timestamps (1720,) acquisition/experiment_ids/ acquisition/experiment_ids/data (1720,) acquisition/experiment_ids/timestamps (1720,) analysis/ file_create_date (1,) general/ general/devices/ general/devices/NLX-microwires-177/ general/devices/NLX-microwires-178/ general/devices/NLX-microwires-179/ general/devices/NLX-microwires-180/ general/devices/NLX-microwires-182/ general/devices/NLX-microwires-183/ general/devices/NLX-microwires-184/ general/devices/NLX-microwires-190/ general/devices/NLX-microwires-191/ general/devices/NLX-microwires-192/ general/experiment_description () general/extracellular_ephys/ general/extracellular_ephys/NLX-microwires-177/ general/extracellular_ephys/NLX-microwires-178/ general/extracellular_ephys/NLX-microwires-179/ general/extracellular_ephys/NLX-microwires-180/ general/extracellular_ephys/NLX-microwires-182/ general/extracellular_ephys/NLX-microwires-183/ general/extracellular_ephys/NLX-microwires-184/ general/extracellular_ephys/NLX-microwires-190/ general/extracellular_ephys/NLX-microwires-191/ general/extracellular_ephys/NLX-microwires-192/ general/extracellular_ephys/electrodes/ general/extracellular_ephys/electrodes/filtering (10,) general/extracellular_ephys/electrodes/group (10,) general/extracellular_ephys/electrodes/group_name (10,) general/extracellular_ephys/electrodes/id (10,) general/extracellular_ephys/electrodes/imp (10,) general/extracellular_ephys/electrodes/location (10,) general/extracellular_ephys/electrodes/origChannel (10,) general/extracellular_ephys/electrodes/x (10,) general/extracellular_ephys/electrodes/y (10,) general/extracellular_ephys/electrodes/z (10,) general/intracellular_ephys/ general/lab () general/optogenetics/ general/optophysiology/ general/related_publications (1,) general/session_id () general/source_script () general/subject/ general/subject/age () general/subject/sex () general/subject/species () general/subject/subject_id () identifier () intervals/ intervals/encoding_table/ intervals/encoding_table/Clip_name (90,) intervals/encoding_table/ExperimentID (90,) intervals/encoding_table/boundary1_time (90,) intervals/encoding_table/boundary2_time (90,) intervals/encoding_table/boundary3_time (90,) intervals/encoding_table/fixcross_time (90,) intervals/encoding_table/id (90,) intervals/encoding_table/start_time (90,) intervals/encoding_table/stimCategory (90,) intervals/encoding_table/stop_time (90,) intervals/recognition_table/ intervals/recognition_table/ExperimentID (180,) intervals/recognition_table/RT (180,) intervals/recognition_table/accuracy (180,) intervals/recognition_table/boundary_type (180,) intervals/recognition_table/confidence (180,) intervals/recognition_table/fixcross_time (180,) intervals/recognition_table/frameName (180,) intervals/recognition_table/id (180,) intervals/recognition_table/old_new (180,) intervals/recognition_table/resp_value (180,) intervals/recognition_table/response_time (180,) intervals/recognition_table/start_time (180,) intervals/recognition_table/stimuli_type (180,) intervals/recognition_table/stop_time (180,) intervals/recognition_table/trial_num (180,) intervals/timediscrimination_table/ intervals/timediscrimination_table/ExperimentID (180,) intervals/timediscrimination_table/RT (180,) intervals/timediscrimination_table/accuracy (180,) intervals/timediscrimination_table/boundary_type (180,) intervals/timediscrimination_table/confidence (180,) intervals/timediscrimination_table/fixcross_time (180,) intervals/timediscrimination_table/frameName (180,) intervals/timediscrimination_table/id (180,) intervals/timediscrimination_table/key (180,) intervals/timediscrimination_table/leftright (180,) intervals/timediscrimination_table/resp_value (180,) intervals/timediscrimination_table/response_time (180,) intervals/timediscrimination_table/start_time (180,) intervals/timediscrimination_table/stop_time (180,) intervals/timediscrimination_table/trial_num (180,) processing/ scratch/ session_description () session_start_time () specifications/ specifications/core/ specifications/core/2.4.0/ specifications/core/2.4.0/namespace () specifications/core/2.4.0/nwb.base () specifications/core/2.4.0/nwb.behavior () specifications/core/2.4.0/nwb.device () specifications/core/2.4.0/nwb.ecephys () specifications/core/2.4.0/nwb.epoch () specifications/core/2.4.0/nwb.file () specifications/core/2.4.0/nwb.icephys () specifications/core/2.4.0/nwb.image () specifications/core/2.4.0/nwb.misc () specifications/core/2.4.0/nwb.ogen () specifications/core/2.4.0/nwb.ophys () specifications/core/2.4.0/nwb.retinotopy () specifications/hdmf-common/ specifications/hdmf-common/1.5.0/ specifications/hdmf-common/1.5.0/base () specifications/hdmf-common/1.5.0/namespace () specifications/hdmf-common/1.5.0/sparse () specifications/hdmf-common/1.5.0/table () stimulus/ stimulus/presentation/ stimulus/presentation/StimulusPresentation_encoding/ stimulus/presentation/StimulusPresentation_encoding/data (2, 2, 2) stimulus/presentation/StimulusPresentation_encoding/dimension (3,) stimulus/presentation/StimulusPresentation_encoding/distance () stimulus/presentation/StimulusPresentation_encoding/external_file (90,) stimulus/presentation/StimulusPresentation_encoding/field_of_view (3,) stimulus/presentation/StimulusPresentation_encoding/format () stimulus/presentation/StimulusPresentation_encoding/orientation () stimulus/presentation/StimulusPresentation_encoding/timestamps (90,) stimulus/templates/ timestamps_reference_time () units/ units/electrodes (43,) units/id (43,) units/spike_times (603265,) units/spike_times_index (43,)
Are the spikes on the same clock as the trial tables?¶
Two checks. First, the first and last spike of a file should fall inside the span of its
events (the table above). Second, many units change their firing when a clip starts, so the
change in spiking around clip onsets should be largest when the spike clock is left as it is,
and not when it is shifted. cb.onset_locking measures that change for a range of shifts.
The three files that differ most between the pinned version and the most recent version of the dandiset (subjects 13, 17 and 18) are read from the most recent version too.
lags = np.round(np.arange(-30, 30.01, 0.1), 1)
rows = []
for s in good:
score = cb.onset_locking(s, lags)
rows.append({"subject": s.subject, "units": len(s.units),
"shift with the largest change (s)": lags[np.argmax(score)],
"change there": round(score.max(), 2),
"median over shifts": round(float(np.median(score)), 2),
"peak / median": round(score.max() / np.median(score), 1)})
locking = pd.DataFrame(rows).set_index("subject")
print("Pinned version", cb.VERSION)
print(locking.to_string())
clear = locking[locking["peak / median"] >= 2.5]
print(f"\nSessions with a clear peak (at least 2.5 times the median): {len(clear)} of {len(locking)}; "
f"their peaks lie at shifts of {clear.iloc[:, 1].min()} to {clear.iloc[:, 1].max()} s.")
unclear = locking[locking["peak / median"] < 2.5]
near_boundary = unclear.iloc[:, 1].between(4.0, 5.0)
print(f"Of the other {len(unclear)}, {int(near_boundary.sum())} have their largest change at a shift of "
f"{unclear.iloc[:, 1][near_boundary].min()} to {unclear.iloc[:, 1][near_boundary].max()} s, shortly after the "
"4 s at which hard boundaries occur.")
wide = np.round(np.arange(-250, 250.01, 0.1), 1)
latest = cb.load_sessions(cb.LATEST_VERSION, subjects=("13", "17", "18"))
pinned = {s.subject: s for s in sessions}
rows = []
for s in latest:
old = pinned[s.subject]
score = cb.onset_locking(s, wide)
rows.append({
"subject": s.subject,
"first spike, pinned (s)": round(old.sanity["first spike (s)"], 3),
"first spike, latest (s)": round(s.sanity["first spike (s)"], 3),
"last spike - last event, latest (s)": round(s.sanity["last spike (s)"] - s.sanity["last event (s)"], 1),
"encoding onsets moved by (s)": round(float((s.encoding.start_time - old.encoding.start_time).abs().max()), 3),
"recognition onsets moved by (s)": round(float((s.recognition.start_time - old.recognition.start_time).mean()), 2),
"order onsets moved by (s)": round(float((s.order.start_time - old.order.start_time).mean()), 2),
"shift with the largest change (s)": wide[np.argmax(score)],
"peak / median": round(score.max() / np.median(score), 1),
})
print("\nMost recent version", cb.LATEST_VERSION, "against the pinned version")
print(pd.DataFrame(rows).set_index("subject").T.to_string())
def soft_rows_kept(encoding):
# Rows the authors' plotting code would keep as soft boundaries: stored time above 0.
sb, hb = encoding[encoding.kind == "SB"], encoding[encoding.kind == "HB"]
stored = np.concatenate([sb.boundary1_time, sb.boundary2_time, sb.boundary3_time,
hb.boundary2_time, hb.boundary3_time])
return int((stored > 0).sum()), len(stored)
print("\nSoft-boundary rows with a stored time above 0, which is how the authors' plotting code drops "
"missing boundaries (75 soft boundaries exist per session):")
print(" pinned version, subject 13: %d of %d" % soft_rows_kept(pinned["13"].encoding))
print(" most recent version, subject 13: %d of %d" % soft_rows_kept(latest[0].encoding),
"| NaN in boundary2_time:", int(latest[0].encoding.boundary2_time.isna().sum()))
# A peak 4.3 to 4.7 s after clip onset could be a response to the hard boundaries at 4 s, or a
# spike clock that is off by about 4 s. NB clips have no event at 4 s, so the same scan on NB
# clips alone tells the two apart: an offset clock would put their peak at the same shift.
rows = []
for s in good:
row = {"subject": s.subject}
for kind in cb.KINDS:
only = dataclasses.replace(s, encoding=s.encoding[s.encoding.kind == kind])
score = cb.onset_locking(only, lags)
row[f"{kind}: shift of largest change (s)"] = lags[np.argmax(score)]
row[f"{kind}: peak / median"] = round(score.max() / np.median(score), 1)
rows.append(row)
by_kind = pd.DataFrame(rows).set_index("subject")
print("\nThe same scan on each kind of clip alone:")
print(by_kind.to_string())
late = locking.index[(locking["peak / median"] < 2.5) & locking.iloc[:, 1].between(4.0, 5.0)]
hb_late = by_kind.loc[late, "HB: shift of largest change (s)"].between(4.0, 5.0)
nb_late = by_kind.loc[late, "NB: shift of largest change (s)"].between(4.0, 5.0)
print(f"Of the {len(late)} sessions whose largest change over all clips is 4 to 5 s after onset: "
f"HB clips alone peak 4 to 5 s after onset in {int(hb_late.sum())}, NB clips alone in {int(nb_late.sum())}; "
f"NB peak / median in those sessions: {by_kind.loc[late, 'NB: peak / median'].min()} to "
f"{by_kind.loc[late, 'NB: peak / median'].max()}.")
Pinned version 0.220216.0323
units shift with the largest change (s) change there median over shifts peak / median
subject
1 43 4.5 1.32 0.81 1.6
2 45 0.4 2.09 0.87 2.4
3 31 12.0 1.24 0.83 1.5
4 35 4.3 1.32 0.83 1.6
5 30 4.7 1.22 0.81 1.5
6 19 4.5 1.39 0.83 1.7
7 34 4.4 1.56 0.85 1.8
8 33 4.5 1.49 0.84 1.8
9 31 -13.8 1.47 0.85 1.7
10 27 -5.3 1.53 0.84 1.8
11 14 0.5 3.98 0.92 4.3
12 28 4.6 1.46 0.84 1.7
13 48 4.4 1.60 0.82 2.0
14 30 0.5 2.46 0.90 2.7
15 31 0.5 3.72 0.93 4.0
16 14 0.5 8.31 1.02 8.1
17 27 0.5 4.11 0.95 4.3
18 29 0.5 3.91 0.93 4.2
Sessions with a clear peak (at least 2.5 times the median): 6 of 18; their peaks lie at shifts of 0.5 to 0.5 s.
Of the other 12, 8 have their largest change at a shift of 4.3 to 4.7 s, shortly after the 4 s at which hard boundaries occur.
Most recent version 0.230530.1822 against the pinned version subject 13 17 18 first spike, pinned (s) 0.005 0.017 0.003 first spike, latest (s) 191.536 -16.815 9.315 last spike - last event, latest (s) 132.300 -106.200 -65.100 encoding onsets moved by (s) 0.006 0.009 0.005 recognition onsets moved by (s) 39.440 32.290 -53.530 order onsets moved by (s) 59.260 89.380 74.400 shift with the largest change (s) 195.900 -16.300 9.800 peak / median 1.800 4.600 4.600 Soft-boundary rows with a stored time above 0, which is how the authors' plotting code drops missing boundaries (75 soft boundaries exist per session): pinned version, subject 13: 150 of 150 most recent version, subject 13: 75 of 150 | NaN in boundary2_time: 60
The same scan on each kind of clip alone:
NB: shift of largest change (s) NB: peak / median SB: shift of largest change (s) SB: peak / median HB: shift of largest change (s) HB: peak / median
subject
1 5.1 1.4 3.3 1.4 4.5 1.6
2 0.4 1.7 0.5 1.7 0.4 1.7
3 3.5 1.5 -5.2 1.5 1.7 1.4
4 4.8 1.4 1.3 1.5 4.4 1.8
5 10.2 1.6 -5.5 1.5 4.6 1.8
6 17.5 1.7 17.2 1.5 4.6 1.6
7 -14.2 1.7 -15.9 1.5 4.4 1.8
8 3.5 1.4 -5.1 1.4 4.5 2.2
9 3.5 1.5 16.0 1.5 4.4 2.1
10 -4.9 1.8 22.0 1.7 4.4 1.9
11 0.4 3.1 0.5 3.1 0.4 2.7
12 10.7 1.6 4.5 1.4 4.5 2.4
13 3.5 1.9 1.7 1.6 4.4 2.5
14 0.4 1.8 0.5 2.0 0.5 2.2
15 0.5 2.7 0.5 2.7 0.4 2.8
16 0.5 5.3 0.5 5.0 0.6 5.4
17 0.4 2.8 0.5 3.0 0.5 2.9
18 0.5 3.0 0.5 2.9 0.5 2.7
Of the 8 sessions whose largest change over all clips is 4 to 5 s after onset: HB clips alone peak 4 to 5 s after onset in 8, NB clips alone in 1; NB peak / median in those sessions: 1.4 to 1.9.
Checks against the files' own descriptions and the paper:
- Subjects and units. There are 19 files, one per subject, matching the paper's 19 subjects with usable recordings. They hold 565 units, all in the medial temporal lobe, where the paper analyses 580 there; the split by region also differs (hippocampus 279 against 343, amygdala 233 against 169, parahippocampal gyrus 53 against 68). None of the paper's 405 units outside the medial temporal lobe is in the dandiset, nor is the behaviour of the twentieth subject.
- Spike times of subject 19. They lie near -7.02e12 s, while that file's events run from 0 to 2892 s. There is no way to put them on the clock of the trial tables, so the 16 units of subject 19 are left out of every analysis of neurons. Its behaviour is kept.
- Spike clock of the other 18. The first spike is within 0.11 s of time zero and the last spike within about half a second of the last event. In six sessions, spiking changes clearly around clip onsets, and in each the change is largest 0.5 s after the onset and not at another shift, which is what an aligned clock gives. Eight of the other twelve change most 4.3 to 4.7 s after clip onset, shortly after the time of the hard boundaries. That could also be an onset response displaced by a clock that is off by about 4 s, so the scan was repeated on each kind of clip alone. In all eight, HB clips alone peak 4.4 to 4.6 s after onset, while NB and SB clips show no clear peak at any shift (peak / median 1.4 to 1.9 for NB). A displaced onset response would appear in all three kinds of clip; a peak in HB clips only is a response to the hard boundary at 4 s on an aligned clock. Subjects 9 and 10 show the same HB-only peak at 4.4 s. Subject 3 shows no clear peak for any kind of clip (peak / median 1.4 to 1.5), so its alignment is not established; it is kept, and it has 31 of the 549 units. The scan works in steps of 0.1 s, so offsets smaller than about 0.3 s are not ruled out.
- The most recent version is worse. In version 0.230530.1822 the spike times of subjects 13, 17 and 18 are shifted by 191.5, -16.8 and 9.3 s while the encoding trials moved by less than 0.01 s, and the largest change of spiking is found at a shift of 195.9, -16.3 and 9.8 s instead of 4.4, 0.5 and 0.5 s. Spikes and clips no longer line up there. This notebook therefore uses the version the paper cites, in which they do.
- Units of time. The files do not state a unit for spike times. Read as seconds, sessions last 2675 to 3131 s and the median firing rate is 1.9 Hz (0.01 to 22.9 Hz), which is plausible for sorted human single units, so seconds is right.
- Clips. Clips last 6.9 to 9.2 s (the paper: about 8 s) and follow a fixation of 0.9 to 1.1 s, as in the Methods. The hard boundary is 4.000 s into every HB clip, and NB clips carry a nominal boundary at 4.000 s. SB clips have one to three soft boundaries, the first between 0.55 and 5.85 s. Half of the HB clip presentations (285 of 570) also contain a soft boundary; the column description says so, the paper's Methods do not.
- Condition codes.
stimCategoryholds 0, 1 and 2 where its description says 1, 2 and 3. The condition is taken from the clip's name here; the two agree in every row. - Missing boundaries. A second or third boundary that does not exist is stored as the
clip's own start time (1140 of 1710 rows for
boundary2_time), not as 0 as described. - Accuracy codes. The description gives 0 for "no response", but 0 never occurs. The 57
recognition trials without a response are coded incorrect. Of the 43 time discrimination
trials without a response, 23 are coded incorrect and 20 correct.
resp_valueis 1 for these trials, where its description says -1. - Reaction time.
RTequals the response time minus the time the frames disappeared, to within a millisecond, so it is negative when a subject answered while the frames were still shown (151 of 3420 time discrimination trials). Its description says only "Reaction time". - Events. The
eventsseries holds 43 entries with code 0, which its description does not list, with timestamps as low as -7.0e12 s. They are not used. - Linking test trials to clips. Frame names of NB trials carry the numbers 1 to 30, two trials per number in each test, like the 30 NB clips. Frame names of SB and HB trials do not match the clips: in scene recognition they carry only 10 distinct numbers, with up to 12 trials per number, and in time discrimination they run from 1 to 60, while the clips are numbered 1 to 30. Nothing else in the files says which clip a tested frame came from, or where in the clip it was. Every claim that relates activity during a clip to memory for that clip is therefore blocked for the two conditions it is about.
Behaviour¶
The paper reports each measure as a mean and standard deviation over its 20 subjects and compares the three kinds of clip with a one-way ANOVA. The dandiset has 19 of the 20. The same is done here, with paired tests over subjects added. Trials without a response are kept, as coded in the file, in the main numbers and left out in a second row.
def by_subject(table, column, responded_only=False):
"""Mean of a column for each subject (rows) and kind of clip (columns)."""
t = table[table.responded] if responded_only else table
return t.groupby(["subject", "kind"])[column].mean().unstack()[list(cb.KINDS)]
def describe(per_subject, reported, label):
"""Print the mean and SD over subjects, the paper's ANOVA and paired tests against HB."""
summary = pd.DataFrame({"mean": per_subject.mean(), "SD": per_subject.std(),
"paper": pd.Series(reported)}).reindex(list(cb.KINDS)).round(3)
print(label, f"(n = {len(per_subject)} subjects)")
print(summary.T.to_string())
f = stats.f_oneway(*(per_subject[k] for k in cb.KINDS))
print(f" one-way ANOVA over kinds: F(2, {3 * len(per_subject) - 3}) = {f.statistic:.2f}, p = {f.pvalue:.2g}")
for other in ("NB", "SB"):
t = stats.ttest_rel(per_subject["HB"], per_subject[other])
diff = per_subject["HB"] - per_subject[other]
print(f" HB - {other}: mean {diff.mean():+.3f}, paired t = {t.statistic:.2f}, p = {t.pvalue:.2g}; "
f"HB is lower in {(diff < 0).sum()} and higher in {(diff > 0).sum()} of {len(diff)} subjects")
return summary
def dots(ax, per_subject, ylabel):
for i, k in enumerate(cb.KINDS):
x = i + rng.uniform(-0.12, 0.12, len(per_subject))
ax.plot(x, per_subject[k], "o", color=COLOR[k], alpha=0.6, ms=5)
ax.plot([i - 0.25, i + 0.25], [per_subject[k].mean()] * 2, "k-", lw=2)
ax.set_xticks(range(3), cb.KINDS)
ax.set_ylabel(ylabel)
Claim fig2a-order-accuracy¶
"Subjects correctly identified which frame was shown first in 73 ± 7% and 73 ± 8% of trials when the two frames were separated by no boundary (NB) or a soft boundary (SB) ... In contrast, subjects performed significantly worse when discriminating the order of frames separated by a hard boundary (HB: 53 ± 5%)."
accuracy = by_subject(order, "correct")
describe(accuracy, {"NB": 0.73, "SB": 0.73, "HB": 0.53}, "Time discrimination accuracy, all trials")
for k in cb.KINDS:
t = stats.ttest_1samp(accuracy[k], 0.5)
print(f" {k} against chance (0.5): t({len(accuracy) - 1}) = {t.statistic:.2f}, p = {t.pvalue:.2g}")
print()
describe(by_subject(order, "correct", responded_only=True), {"NB": 0.73, "SB": 0.73, "HB": 0.53},
"Time discrimination accuracy, trials with a response")
print("\nBy subject (all trials):")
print(accuracy.round(2).T.to_string())
fig, ax = plt.subplots(figsize=(3.6, 3.4))
dots(ax, accuracy, "Time discrimination accuracy")
ax.axhline(0.5, color="0.6", ls="--", lw=1)
plt.tight_layout()
print(f"\nPlotted values range from {accuracy.min().min():.2f} to {accuracy.max().max():.2f}.")
Time discrimination accuracy, all trials (n = 19 subjects)
NB SB HB
mean 0.728 0.720 0.528
SD 0.074 0.093 0.051
paper 0.730 0.730 0.530
one-way ANOVA over kinds: F(2, 54) = 43.41, p = 5.8e-12
HB - NB: mean -0.200, paired t = -10.26, p = 6e-09; HB is lower in 18 and higher in 1 of 19 subjects
HB - SB: mean -0.192, paired t = -8.10, p = 2.1e-07; HB is lower in 18 and higher in 1 of 19 subjects
NB against chance (0.5): t(18) = 13.41, p = 8.3e-11
SB against chance (0.5): t(18) = 10.28, p = 5.9e-09
HB against chance (0.5): t(18) = 2.39, p = 0.028
Time discrimination accuracy, trials with a response (n = 19 subjects)
NB SB HB
mean 0.730 0.725 0.529
SD 0.075 0.095 0.050
paper 0.730 0.730 0.530
one-way ANOVA over kinds: F(2, 54) = 43.64, p = 5.3e-12
HB - NB: mean -0.201, paired t = -10.14, p = 7.2e-09; HB is lower in 18 and higher in 1 of 19 subjects
HB - SB: mean -0.196, paired t = -8.30, p = 1.4e-07; HB is lower in 18 and higher in 1 of 19 subjects
By subject (all trials):
subject 1 10 11 12 13 14 15 16 17 18 19 2 3 4 5 6 7 8 9
kind
NB 0.77 0.75 0.68 0.77 0.67 0.67 0.72 0.65 0.55 0.73 0.88 0.82 0.78 0.73 0.78 0.67 0.80 0.72 0.70
SB 0.75 0.73 0.68 0.73 0.65 0.65 0.72 0.65 0.48 0.65 0.92 0.85 0.82 0.75 0.77 0.67 0.80 0.68 0.73
HB 0.57 0.57 0.47 0.58 0.57 0.48 0.53 0.52 0.57 0.47 0.58 0.43 0.57 0.47 0.50 0.55 0.62 0.52 0.48
Plotted values range from 0.43 to 0.92.
Result: reproduced. Accuracy is 0.53 across a hard boundary (paper: 0.53), against 0.73 with no boundary and 0.72 with a soft boundary (paper: 0.73 for both). It is lower across a hard boundary in 18 of 19 subjects (paired p < 1e-6 against both), and only just above chance (p = 0.028; the paper gives p = 0.02). Leaving out trials without a response moves no mean by more than 0.01.
Claim fig2b-order-reaction-time¶
"Subjects also showed longer reaction times (HB: 2.10 ± 0.37 seconds; NB: 1.62 ± 0.28 seconds; SB: 1.59 ± 0.34 seconds)."
The file's RT column is used. The data check showed that it is measured from the moment the
frames disappear, about 2 s after they appear, so it is negative when a subject answers while
the frames are still shown. Trials without a response have an RT of 0 in the file.
rt = by_subject(order, "RT")
describe(rt, {"NB": 1.62, "SB": 1.59, "HB": 2.10}, "Time discrimination RT (s), all trials")
print()
describe(by_subject(order, "RT", responded_only=True), {"NB": 1.62, "SB": 1.59, "HB": 2.10},
"Time discrimination RT (s), trials with a response")
from_onset = by_subject(order.assign(latency=order.response_time - order.start_time), "latency", responded_only=True)
print("\nMeasured from the appearance of the frames instead (s), trials with a response:",
from_onset.mean().round(2).to_dict())
fig, ax = plt.subplots(figsize=(3.6, 3.4))
dots(ax, rt, "Reaction time (s)")
plt.tight_layout()
print(f"Plotted values range from {rt.min().min():.2f} to {rt.max().max():.2f} s.")
Time discrimination RT (s), all trials (n = 19 subjects)
NB SB HB
mean 1.624 1.535 2.054
SD 0.316 0.382 0.384
paper 1.620 1.590 2.100
one-way ANOVA over kinds: F(2, 54) = 11.15, p = 8.8e-05
HB - NB: mean +0.430, paired t = 6.17, p = 8e-06; HB is lower in 1 and higher in 18 of 19 subjects
HB - SB: mean +0.519, paired t = 8.05, p = 2.2e-07; HB is lower in 1 and higher in 18 of 19 subjects
Time discrimination RT (s), trials with a response (n = 19 subjects)
NB SB HB
mean 1.638 1.556 2.081
SD 0.307 0.377 0.386
paper 1.620 1.590 2.100
one-way ANOVA over kinds: F(2, 54) = 11.78, p = 5.7e-05
HB - NB: mean +0.443, paired t = 6.27, p = 6.5e-06; HB is lower in 1 and higher in 18 of 19 subjects
HB - SB: mean +0.525, paired t = 7.92, p = 2.8e-07; HB is lower in 1 and higher in 18 of 19 subjects
Measured from the appearance of the frames instead (s), trials with a response: {'NB': 3.67, 'SB': 3.59, 'HB': 4.11}
Plotted values range from 0.78 to 2.65 s.
Result: reproduced. The mean reaction time is 2.05 s across a hard boundary (paper: 2.10 s), against 1.62 s with no boundary and 1.54 s with a soft boundary (paper: 1.62 and 1.59 s). It is longer across a hard boundary in 18 of 19 subjects (paired p < 1e-5 against both). These times run from the disappearance of the frames; counted from their appearance they are about 2 s longer (3.67, 3.59 and 4.11 s), in the same order.
Claim fig2c-order-confidence¶
"... and lower confidence ratings when discriminating the order of frames separated by a hard boundary (HB: 1.95 ± 0.45; NB: 2.52 ± 0.29; SB: 2.59 ± 0.23)." Confidence is 3 for "sure", 2 for "less sure" and 1 for "very unsure". Trials without a response have a confidence of 0 in the file.
confidence = by_subject(order, "confidence")
describe(confidence, {"NB": 2.52, "SB": 2.59, "HB": 1.95}, "Time discrimination confidence, all trials")
print()
describe(by_subject(order, "confidence", responded_only=True), {"NB": 2.52, "SB": 2.59, "HB": 1.95},
"Time discrimination confidence, trials with a response")
fig, ax = plt.subplots(figsize=(3.6, 3.4))
dots(ax, confidence, "Confidence (1 to 3)")
plt.tight_layout()
print(f"\nPlotted values range from {confidence.min().min():.2f} to {confidence.max().max():.2f}.")
Time discrimination confidence, all trials (n = 19 subjects)
NB SB HB
mean 2.454 2.512 1.975
SD 0.263 0.213 0.377
paper 2.520 2.590 1.950
one-way ANOVA over kinds: F(2, 54) = 19.26, p = 4.9e-07
HB - NB: mean -0.479, paired t = -5.51, p = 3.1e-05; HB is lower in 18 and higher in 1 of 19 subjects
HB - SB: mean -0.537, paired t = -6.68, p = 2.9e-06; HB is lower in 18 and higher in 1 of 19 subjects
Time discrimination confidence, trials with a response (n = 19 subjects)
NB SB HB
mean 2.476 2.549 2.003
SD 0.245 0.199 0.386
paper 2.520 2.590 1.950
one-way ANOVA over kinds: F(2, 54) = 20.17, p = 2.9e-07
HB - NB: mean -0.473, paired t = -5.11, p = 7.3e-05; HB is lower in 17 and higher in 2 of 19 subjects
HB - SB: mean -0.546, paired t = -6.74, p = 2.6e-06; HB is lower in 18 and higher in 1 of 19 subjects
Plotted values range from 1.25 to 2.82.
Result: reproduced. Mean confidence is 1.98 across a hard boundary (paper: 1.95), against 2.45 with no boundary and 2.51 with a soft boundary (paper: 2.52 and 2.59). It is lower across a hard boundary in 18 of 19 subjects (paired p < 1e-4 against both). Without the trials that have no response the means are 2.00, 2.48 and 2.55.
Claim edfig2a-recognition-similar¶
"Across all trials, the ability to recognize a frame as 'old' did not differ significantly between the types of boundaries preceding the frame (NB: 76% ± 10%; SB: 75% ± 9%; HB: 75% ± 8%, F(2, 57) = 0.07, p = 0.94)." The text does not say whether foils are counted; both are shown.
reported = {"NB": 0.76, "SB": 0.75, "HB": 0.75}
recognition = by_subject(rec, "correct")
describe(recognition, reported, "Scene recognition accuracy, targets and foils")
print()
describe(by_subject(rec[rec.target], "correct"), reported, "Scene recognition accuracy, targets only (hit rate)")
print()
describe(by_subject(rec, "correct", responded_only=True), reported,
"Scene recognition accuracy, targets and foils, trials with a response")
fig, ax = plt.subplots(figsize=(3.6, 3.4))
dots(ax, recognition, "Scene recognition accuracy")
ax.axhline(0.5, color="0.6", ls="--", lw=1)
plt.tight_layout()
print(f"\nPlotted values range from {recognition.min().min():.2f} to {recognition.max().max():.2f}.")
Scene recognition accuracy, targets and foils (n = 19 subjects)
NB SB HB
mean 0.738 0.739 0.718
SD 0.097 0.074 0.069
paper 0.760 0.750 0.750
one-way ANOVA over kinds: F(2, 54) = 0.38, p = 0.69
HB - NB: mean -0.019, paired t = -1.27, p = 0.22; HB is lower in 13 and higher in 4 of 19 subjects
HB - SB: mean -0.020, paired t = -1.70, p = 0.11; HB is lower in 11 and higher in 6 of 19 subjects
Scene recognition accuracy, targets only (hit rate) (n = 19 subjects)
NB SB HB
mean 0.737 0.718 0.688
SD 0.147 0.132 0.163
paper 0.760 0.750 0.750
one-way ANOVA over kinds: F(2, 54) = 0.53, p = 0.59
HB - NB: mean -0.049, paired t = -1.78, p = 0.092; HB is lower in 12 and higher in 4 of 19 subjects
HB - SB: mean -0.030, paired t = -1.24, p = 0.23; HB is lower in 11 and higher in 8 of 19 subjects
Scene recognition accuracy, targets and foils, trials with a response (n = 19 subjects)
NB SB HB
mean 0.751 0.748 0.734
SD 0.092 0.076 0.074
paper 0.760 0.750 0.750
one-way ANOVA over kinds: F(2, 54) = 0.24, p = 0.79
HB - NB: mean -0.017, paired t = -1.31, p = 0.21; HB is lower in 15 and higher in 4 of 19 subjects
HB - SB: mean -0.013, paired t = -1.27, p = 0.22; HB is lower in 11 and higher in 7 of 19 subjects
Plotted values range from 0.60 to 0.95.
Result: reproduced. Accuracy over targets and foils is 0.74, 0.74 and 0.72 for NB, SB and HB (paper: 0.76, 0.75 and 0.75), each within 0.04 of the reported value, and the paper's one-way ANOVA finds no difference (F(2, 54) = 0.38, p = 0.69; paper: F(2, 57) = 0.07, p = 0.94). The ANOVA also finds no difference for targets alone (p = 0.59) and without the trials that have no response (p = 0.79). For targets alone, though, accuracy for HB (0.69) is 0.06 below the paper's 0.75, outside the 0.04 of the criterion. The paper does not say whether foils are counted; targets and foils are counted here. Accuracy for HB is the lowest of the three in each version, by 0.01 to 0.05, and no paired test of that difference reaches p < 0.05.
Boundary cells and event cells¶
From the Methods. For each neuron, spikes are counted in the second after and the second before each boundary during encoding; clips with no boundary are aligned to their middle.
- Boundary cell: the count after the boundary differs from the count before for SB and for HB but not for NB (p < 0.05, permutation t-test), and the count after the boundary is greater for SB and for HB than for NB (p < 0.05).
- Event cell: the count differs for HB but not for NB or SB, and the count after the boundary is greater for HB than for NB and for SB.
- Chance: the same procedure after shuffling the condition labels of the clips.
What the Methods leave open, and what is done here:
- The permutation. Before against after: the sign of each trial's difference is permuted. Between conditions: the condition labels are permuted, one-sided. 1000 permutations per test.
- Which soft boundaries. The authors' plotting code takes every soft boundary of a session as an SB trial: those of the SB clips and those inside HB clips, 75 per session against 30 NB and 30 HB trials. That set is used. Two alternatives are shown: the first soft boundary of each SB clip (30 trials, the choice of the first version of this notebook) and all soft boundaries of SB clips (50 trials).
- How "for SB and HB" is tested. Each condition is tested on its own, five tests per class. A
looser reading that pools conditions (three tests per class,
cb.find_cells_pooled) is shown as a variant. - Chance. 300 shuffles of the condition labels over a session's boundaries are used where the paper used 1000.
Only the 18 sessions whose spike clock matches the trial tables are used.
cells = pd.concat([cb.find_cells(s, rng) for s in good], ignore_index=True)
tests = [c for c in cells.columns if " vs " in c or " > " in c]
n_units = len(cells)
print(f"{n_units} units in {len(good)} sessions.")
print("Fraction of units with p < 0.05 in each test (5% is expected of units that do not respond):")
print((cells[tests] < cb.ALPHA).mean().round(3).to_string())
silent = 0
for s in good:
times, _ = cb.boundary_times(s.encoding)
silent += int((cb.spike_counts(s.units, times, -cb.WINDOW, cb.WINDOW).sum(axis=1) == 0).sum())
print(f"Units with no spike in any of the windows, which no test can select: {silent}")
print("Boundaries per session:", {k: int((cb.boundary_times(good[0].encoding)[1] == k).sum()) for k in cb.KINDS})
SHUFFLES = 300
chance = cb.chance_counts(good, rng, SHUFFLES)
chance_first = cb.chance_counts(good, rng, SHUFFLES, soft="first")
# HB clips in which no soft boundary comes before the hard boundary (see the per-clip table below).
clean = [dataclasses.replace(s, encoding=s.encoding[~((s.encoding.kind == "HB") & (s.encoding.b2 < s.encoding.b1))])
for s in good]
print("HB clips without a soft boundary before the hard boundary:",
int((clean[0].encoding.kind == "HB").sum()), "of 30")
def everyone(function, group=good, **options):
return pd.concat([function(s, rng, **options) for s in group], ignore_index=True)
variants = {
"five tests, authors' SB trials (75) [main]": cells,
"five tests, first SB of SB clips (30)": everyone(cb.find_cells, soft="first"),
"five tests, all SBs of SB clips (50)": everyone(cb.find_cells, soft="all"),
"main, alpha = 0.10": everyone(cb.find_cells, alpha=0.10),
"main, HB clips without an earlier SB": everyone(cb.find_cells, group=clean),
"pooled reading (three tests), authors' SB trials": everyone(cb.find_cells_pooled),
"pooled reading (three tests), first SB": everyone(cb.find_cells_pooled, soft="first"),
}
variant_table = pd.DataFrame({name: {"boundary cells": int(v.boundary.sum()), "event cells": int(v.event.sum()),
"both": int((v.boundary & v.event).sum())}
for name, v in variants.items()}).T
repeats = [pd.concat([cb.find_cells(s, rng) for s in good], ignore_index=True) for _ in range(20)]
repeat_counts = pd.DataFrame({"boundary": [int(r.boundary.sum()) for r in repeats],
"event": [int(r.event.sum()) for r in repeats]})
variant_table["boundary %"] = (100 * variant_table["boundary cells"] / n_units).round(2)
variant_table["event %"] = (100 * variant_table["event cells"] / n_units).round(2)
chance_alpha10 = cb.chance_counts(good, rng, SHUFFLES, alpha=0.10)
chance_pooled = cb.chance_counts(good, rng, SHUFFLES, pooled=True)
549 units in 18 sessions.
Fraction of units with p < 0.05 in each test (5% is expected of units that do not respond):
NB post vs pre 0.111
SB post vs pre 0.113
HB post vs pre 0.104
SB > NB 0.186
HB > NB 0.093
HB > SB 0.062
Units with no spike in any of the windows, which no test can select: 0
Boundaries per session: {'NB': 30, 'SB': 75, 'HB': 30}
HB clips without a soft boundary before the hard boundary: 17 of 30
Activity before the hard boundary, clip by clip¶
The test of a response to a hard boundary compares the second after it with the second before.
Half of the HB clips contain soft boundaries too, some of them shortly before the hard boundary.
The table shows, for each HB clip, how far the firing of the whole population departs from the
average clip (cb.clip_population_trace, roughly a z-score) in the half second before the hard
boundary and shortly after it.
edges = np.round(np.arange(0.0, 8.001, 0.1), 1)
by_clip = cb.clip_population_trace(good, edges)
t = by_clip.columns.to_numpy()
pre_window, post_window = (t > 3.5) & (t < 4.0), (t > 4.1) & (t < 4.5)
info = good[0].encoding.set_index("clip")
table = pd.DataFrame({
"soft boundaries (s)": [", ".join(f"{b:.2f}" for b in (info.b2[c], info.b3[c]) if np.isfinite(b)) or "-"
for c in by_clip.index],
"3.5 to 4.0 s": by_clip.loc[:, pre_window].mean(axis=1).round(1),
"4.1 to 4.5 s": by_clip.loc[:, post_window].mean(axis=1).round(1),
"largest bin": by_clip.max(axis=1).round(1),
"at (s)": by_clip.idxmax(axis=1),
}, index=by_clip.index)
hb = table[table.index.str.startswith("HB")]
hb = hb.loc[sorted(hb.index, key=lambda c: int(c.split("_")[1]))]
print("HB clips; the hard boundary is at 4.00 s. Population deviation from the average clip:")
print(hb.to_string())
for kind in cb.KINDS:
rows = table[table.index.str.startswith(kind)]
print(f"{kind} clips: mean deviation 3.5 to 4.0 s {rows['3.5 to 4.0 s'].mean():+.2f} "
f"(largest {rows['3.5 to 4.0 s'].max():+.1f}), 4.1 to 4.5 s {rows['4.1 to 4.5 s'].mean():+.2f} "
f"(largest {rows['4.1 to 4.5 s'].max():+.1f})")
high = hb[hb["3.5 to 4.0 s"] > 2].index.tolist()
rest = hb.drop(index=high)
print(f"HB clips with a deviation above 2 before the hard boundary: {high}")
print(f"The other {len(rest)} HB clips: mean deviation 3.5 to 4.0 s {rest['3.5 to 4.0 s'].mean():+.2f}, "
f"4.1 to 4.5 s {rest['4.1 to 4.5 s'].mean():+.2f}; {int((rest['4.1 to 4.5 s'] > 0).sum())} of {len(rest)} "
"are above the average clip after the boundary.")
fig, ax = plt.subplots(figsize=(6.4, 3.2))
for c in hb.index:
ax.plot(t, by_clip.loc[c], color=COLOR["HB"] if c in high else "0.7", lw=1.2 if c in high else 0.6,
zorder=3 if c in high else 1)
ax.axvline(4.0, color="k", lw=1, ls="--")
ax.set_xlabel("Time from clip onset (s)"); ax.set_ylabel("Population deviation")
ax.set_title("HB clips; red: above 2 before the hard boundary", fontsize=9)
plt.tight_layout()
print(f"Plotted values range from {by_clip.loc[hb.index].min().min():.1f} to {by_clip.loc[hb.index].max().max():.1f}.")
HB clips; the hard boundary is at 4.00 s. Population deviation from the average clip:
soft boundaries (s) 3.5 to 4.0 s 4.1 to 4.5 s largest bin at (s)
HB_1 1.55 -0.3 2.0 6.4 1.35
HB_2 2.36 -1.6 1.6 4.4 4.55
HB_3 2.31 -1.1 1.2 9.2 2.45
HB_4 6.08 -0.0 -0.7 5.2 4.95
HB_5 6.06 -0.2 1.6 2.9 6.35
HB_6 2.25, 6.28 -0.2 2.6 4.2 4.35
HB_7 3.40, 6.07 10.8 7.6 13.1 3.95
HB_8 1.85, 5.87 0.5 1.3 6.2 1.85
HB_9 2.46, 5.87 -0.4 1.6 5.4 2.65
HB_10 1.56, 5.93 -1.0 1.4 2.6 4.25
HB_11 2.88, 6.22 8.6 -0.3 12.8 3.75
HB_12 1.64, 5.93 -0.4 0.3 3.5 6.15
HB_13 2.26, 6.21 -0.3 1.6 3.6 4.25
HB_14 2.75, 5.92 2.7 0.5 7.3 3.65
HB_15 3.23, 6.26 9.3 4.0 10.7 4.05
HB_16 - -0.8 1.7 3.3 4.35
HB_17 - -0.1 0.2 3.6 6.85
HB_18 - -0.4 0.8 2.9 5.25
HB_19 - -0.7 1.6 3.6 4.25
HB_20 - -1.0 1.5 4.9 4.35
HB_21 - 0.3 1.0 6.2 4.75
HB_22 - -0.3 1.1 3.5 4.55
HB_23 - -1.3 0.3 3.9 1.35
HB_24 - 0.0 1.6 4.0 4.25
HB_25 - -0.8 0.7 2.7 6.05
HB_26 - -1.6 -0.7 2.6 6.35
HB_27 - -1.2 0.7 3.2 6.95
HB_28 - -0.1 0.7 2.9 4.25
HB_29 - -0.5 0.6 2.6 6.35
HB_30 - -1.3 0.4 2.4 6.65
NB clips: mean deviation 3.5 to 4.0 s -0.58 (largest +0.9), 4.1 to 4.5 s -0.90 (largest +0.8)
SB clips: mean deviation 3.5 to 4.0 s +0.02 (largest +5.3), 4.1 to 4.5 s -0.37 (largest +1.4)
HB clips: mean deviation 3.5 to 4.0 s +0.55 (largest +10.8), 4.1 to 4.5 s +1.28 (largest +7.6)
HB clips with a deviation above 2 before the hard boundary: ['HB_7', 'HB_11', 'HB_14', 'HB_15']
The other 26 HB clips: mean deviation 3.5 to 4.0 s -0.57, 4.1 to 4.5 s +1.03; 24 of 26 are above the average clip after the boundary.
Plotted values range from -4.0 to 13.1.
Four HB clips (HB_7, HB_11, HB_14 and HB_15) drive the population well above the average clip in the half second before the hard boundary (deviations of 2.7 to 10.8, where no NB clip exceeds 0.9). All four contain a soft boundary between 2.75 and 3.40 s. In the other 26 HB clips the population is below average before the hard boundary (-0.57) and above it afterwards (+1.03), and 24 of the 26 are above average after the boundary. So most HB clips do show a rise after the hard boundary, and a few add strong activity to the baseline window of the before-against-after test. The cell counts below include a variant that leaves out the 13 HB clips with a soft boundary before the hard one.
Claim fig3c-boundary-cells¶
"42/580 MTL neurons (7.24%; expected proportion by chance for all MTL neurons = 2.11%) were classified as boundary cells."
def report(label):
found = int(cells[label].sum())
null = chance[f"{label} all"]
print(f"{label.capitalize()} cells: {found} of {n_units} units ({100 * found / n_units:.2f}%).")
print(f"After shuffling condition labels ({SHUFFLES} shuffles): mean {null.mean():.2f} cells "
f"({100 * null.mean() / n_units:.2f}%), 95th percentile {np.percentile(null, 95):.0f}, "
f"maximum {null.max()}; shuffles with at least {found} cells: {(null >= found).sum()} of {SHUFFLES}.")
print("By subject:", cells.groupby("subject")[label].sum()[lambda x: x > 0].to_dict(),
f"(none in the other {cells.subject.nunique() - (cells.groupby('subject')[label].sum() > 0).sum()} sessions)")
null_first = chance_first[f"{label} all"]
first = int(variants["five tests, first SB of SB clips (30)"][label].sum())
print(f"With the first soft boundary of SB clips only: {first} cells; shuffles give mean {null_first.mean():.2f}, "
f"95th percentile {np.percentile(null_first, 95):.0f}, maximum {null_first.max()}; "
f"shuffles with at least {first} cells: {(null_first >= first).sum()} of {SHUFFLES}.")
print(f"Repeating the permutation tests 20 times with new random permutations: "
f"{repeat_counts[label].min()} to {repeat_counts[label].max()} cells (median {repeat_counts[label].median():.0f}).")
print("\nVariants of the procedure:")
print(variant_table.to_string())
print("\nChance levels of the looser readings (same shuffles of condition labels):")
for name, null, variant in (("main, alpha = 0.10", chance_alpha10, "main, alpha = 0.10"),
("pooled reading (three tests), authors' SB trials", chance_pooled,
"pooled reading (three tests), authors' SB trials")):
found_v = int(variants[variant][label].sum())
values = null[f"{label} all"]
print(f" {name}: {found_v} cells ({100 * found_v / n_units:.2f}%); shuffles give mean {values.mean():.2f} "
f"({100 * values.mean() / n_units:.2f}%), 95th percentile {np.percentile(values, 95):.0f}, "
f"maximum {values.max()}; shuffles with at least {found_v} cells: {(values >= found_v).sum()} of {SHUFFLES}.")
report("boundary")
paper_fraction = 42 / 580
print(f"\nThe paper's proportion applied to {n_units} units would be {paper_fraction * n_units:.0f} cells. "
f"Probability of {int(cells.boundary.sum())} or fewer if the true proportion were {paper_fraction:.4f}: "
f"{stats.binom.cdf(int(cells.boundary.sum()), n_units, paper_fraction):.1e}")
Boundary cells: 1 of 549 units (0.18%).
After shuffling condition labels (300 shuffles): mean 0.14 cells (0.03%), 95th percentile 1, maximum 2; shuffles with at least 1 cells: 42 of 300.
By subject: {'18': 1} (none in the other 17 sessions)
With the first soft boundary of SB clips only: 1 cells; shuffles give mean 0.08, 95th percentile 1, maximum 2; shuffles with at least 1 cells: 23 of 300.
Repeating the permutation tests 20 times with new random permutations: 0 to 1 cells (median 1).
Variants of the procedure:
boundary cells event cells both boundary % event %
five tests, authors' SB trials (75) [main] 1 7 0 0.18 1.28
five tests, first SB of SB clips (30) 1 5 0 0.18 0.91
five tests, all SBs of SB clips (50) 1 6 0 0.18 1.09
main, alpha = 0.10 2 12 0 0.36 2.19
main, HB clips without an earlier SB 1 2 0 0.18 0.36
pooled reading (three tests), authors' SB trials 15 10 0 2.73 1.82
pooled reading (three tests), first SB 10 10 1 1.82 1.82
Chance levels of the looser readings (same shuffles of condition labels):
main, alpha = 0.10: 2 cells (0.36%); shuffles give mean 0.55 (0.10%), 95th percentile 2, maximum 3; shuffles with at least 2 cells: 35 of 300.
pooled reading (three tests), authors' SB trials: 15 cells (2.73%); shuffles give mean 3.01 (0.55%), 95th percentile 6, maximum 9; shuffles with at least 15 cells: 0 of 300.
The paper's proportion applied to 549 units would be 40 cells. Probability of 1 or fewer if the true proportion were 0.0724: 5.2e-17
Result: blocked, methods underspecified. With the definition as the Methods state it and the authors' set of soft-boundary trials, 1 of 549 units (0.18%) is a boundary cell, where the paper reports 42 of 580 (7.24%). One or more cells also appear in 42 of 300 label shuffles. The count is 0 or 1 over 20 repeats of the permutation tests, 1 with either alternative set of soft boundaries, 2 at alpha = 0.10, and 1 without the HB clips that have an earlier soft boundary.
The outcome is blocked and not "not reproduced" because the procedure that produced the paper's number cannot be determined. The paper's chance level is 2.11%. The literal five-test definition gives 0.03% here, alpha = 0.10 gives 0.10% and the pooled three-test reading gives 0.55%, so no reading tried reaches the paper's chance level, and the paper's count must come from a procedure that differs from all of them. Neither the Methods nor the authors' published code states the permutation test or how the conditions enter the tests; the code has no cell selection at all.
Under the pooled reading there are 15 boundary cells (2.73%), more than in any of 300 shuffles (mean 3.01, maximum 9), so the data hold units that respond to boundaries beyond chance. Their proportion is still well below the paper's. The units are also not the paper's set: 549 are used here against 580, with a different split by region, and the dandiset does not label the paper's boundary cells, so the selections cannot be compared cell by cell.
Claim fig3g-event-cells¶
"36/580 MTL neurons (6.20%; proportion expected by chance for all MTL neurons = 2.26%) were classified as event cells."
report("event")
paper_fraction = 36 / 580
print(f"\nThe paper's proportion applied to {n_units} units would be {paper_fraction * n_units:.0f} cells. "
f"Probability of {int(cells.event.sum())} or fewer if the true proportion were {paper_fraction:.4f}: "
f"{stats.binom.cdf(int(cells.event.sum()), n_units, paper_fraction):.1e}")
# Is there any response to boundaries in the population? Average over all units of the firing
# rate around each kind of boundary, each unit divided by its own mean rate.
centres = np.round(np.arange(-1.0, 1.5001, 0.05), 2)
traces = {k: [] for k in cb.KINDS}
for s in good:
times, kinds = cb.boundary_times(s.encoding)
for u in s.units:
mean_rate = len(u) / (u[-1] - u[0])
for k in cb.KINDS:
traces[k].append(cb.rate_trace(u, times[kinds == k], centres) / mean_rate)
fig, ax = plt.subplots(figsize=(5.2, 3.2))
for k in cb.KINDS:
mean = np.mean(traces[k], axis=0)
sem = np.std(traces[k], axis=0, ddof=1) / np.sqrt(len(traces[k]))
ax.plot(centres, mean, color=COLOR[k], label=k)
ax.fill_between(centres, mean - sem, mean + sem, color=COLOR[k], alpha=0.2, lw=0)
ax.axvline(0, color="0.6", lw=1)
ax.set_xlabel("Time from boundary (s)"); ax.set_ylabel("Firing rate / unit's mean rate")
ax.legend(frameon=False)
plt.tight_layout()
after, before = (centres > 0.1) & (centres <= 0.4), (centres >= -0.4) & (centres < -0.1)
print(f"\nPopulation average over {n_units} units (200 ms windows, mean +- s.e.m. over units):")
for k in cb.KINDS:
change = np.array(traces[k])[:, after].mean(axis=1) - np.array(traces[k])[:, before].mean(axis=1)
print(f" {k}: rate 0.1 to 0.4 s after minus 0.4 to 0.1 s before the boundary, as a fraction of the mean "
f"rate: {change.mean():+.3f} +- {change.std(ddof=1) / np.sqrt(len(change)):.3f}")
for k in ("SB", "HB"):
for name, window in (("0.1 to 0.4 s after", after), ("0.4 to 0.1 s before", before)):
gap = np.array(traces[k])[:, window].mean(axis=1) - np.array(traces["NB"])[:, window].mean(axis=1)
print(f" {k} minus NB, {name} the boundary: {gap.mean():+.3f} +- {gap.std(ddof=1) / np.sqrt(len(gap)):.3f}")
print("Plotted means range from %.2f to %.2f." % (
min(np.mean(traces[k], axis=0).min() for k in cb.KINDS), max(np.mean(traces[k], axis=0).max() for k in cb.KINDS)))
Event cells: 7 of 549 units (1.28%).
After shuffling condition labels (300 shuffles): mean 1.22 cells (0.22%), 95th percentile 3, maximum 6; shuffles with at least 7 cells: 0 of 300.
By subject: {'1': 1, '11': 1, '13': 1, '14': 1, '17': 1, '18': 1, '5': 1} (none in the other 11 sessions)
With the first soft boundary of SB clips only: 5 cells; shuffles give mean 0.99, 95th percentile 3, maximum 7; shuffles with at least 5 cells: 3 of 300.
Repeating the permutation tests 20 times with new random permutations: 5 to 8 cells (median 6).
Variants of the procedure:
boundary cells event cells both boundary % event %
five tests, authors' SB trials (75) [main] 1 7 0 0.18 1.28
five tests, first SB of SB clips (30) 1 5 0 0.18 0.91
five tests, all SBs of SB clips (50) 1 6 0 0.18 1.09
main, alpha = 0.10 2 12 0 0.36 2.19
main, HB clips without an earlier SB 1 2 0 0.18 0.36
pooled reading (three tests), authors' SB trials 15 10 0 2.73 1.82
pooled reading (three tests), first SB 10 10 1 1.82 1.82
Chance levels of the looser readings (same shuffles of condition labels):
main, alpha = 0.10: 12 cells (2.19%); shuffles give mean 3.93 (0.72%), 95th percentile 7, maximum 11; shuffles with at least 12 cells: 0 of 300.
pooled reading (three tests), authors' SB trials: 10 cells (1.82%); shuffles give mean 3.56 (0.65%), 95th percentile 7, maximum 11; shuffles with at least 10 cells: 4 of 300.
The paper's proportion applied to 549 units would be 34 cells. Probability of 7 or fewer if the true proportion were 0.0621: 1.0e-08
Population average over 549 units (200 ms windows, mean +- s.e.m. over units): NB: rate 0.1 to 0.4 s after minus 0.4 to 0.1 s before the boundary, as a fraction of the mean rate: -0.058 +- 0.038 SB: rate 0.1 to 0.4 s after minus 0.4 to 0.1 s before the boundary, as a fraction of the mean rate: +0.011 +- 0.029 HB: rate 0.1 to 0.4 s after minus 0.4 to 0.1 s before the boundary, as a fraction of the mean rate: +0.032 +- 0.030 SB minus NB, 0.1 to 0.4 s after the boundary: +0.156 +- 0.031 SB minus NB, 0.4 to 0.1 s before the boundary: +0.087 +- 0.041 HB minus NB, 0.1 to 0.4 s after the boundary: +0.197 +- 0.028 HB minus NB, 0.4 to 0.1 s before the boundary: +0.107 +- 0.041 Plotted means range from 0.90 to 1.24.
Result: blocked, methods underspecified. 7 of 549 units (1.28%) meet the definition of an event cell. That is more than in any of 300 label shuffles (mean 1.22, maximum 6), but about a fifth of the paper's 36 of 580 (6.20%). The count is 5 to 8 over 20 repeats of the permutation tests. At alpha = 0.10 there are 12 (2.19%), more than in any of 300 shuffles (mean 3.93), and under the pooled reading 10 (1.82%; 4 of 300 shuffles give as many). With the first soft boundary of SB clips as the SB trials there are 5, and without the 13 HB clips that have an earlier soft boundary there are 2.
So units that respond to hard boundaries and not to soft ones are found beyond chance under
most readings, at a proportion far below the paper's under every reading. The chance level
here is 0.22% (0.72% at alpha = 0.10, 0.65% pooled) against the paper's 2.26%. For the reasons
given under fig3c-boundary-cells, the paper's procedure cannot be recovered, so the
difference in proportion cannot be attributed to the data.
Across all units, firing 0.1 to 0.4 s after a hard boundary is higher than at the same time in clips with no boundary (by 0.197 ± 0.028 of a unit's mean rate), but it is already higher in the 0.3 s before the boundary (0.107 ± 0.041); the per-clip table above traces that to four clips.
Claim fig3-cell-regions¶
"The regions with the largest proportion of boundary cells were the parahippocampal gyrus (n = 18/68, 26.47%), amygdala (n = 12/169, 7.10%), and the hippocampus (n = 12/343, 3.50%)." "The regions with the largest proportion of event cells were the hippocampus (n = 27/343, 7.87%), amygdala (n = 7/169, 4.27%), and parahippocampal gyrus (n = 2/68, 2.94%)."
paper_regions = {"boundary": {"hippocampus": (12, 343), "amygdala": (12, 169), "parahippocampus": (18, 68)},
"event": {"hippocampus": (27, 343), "amygdala": (7, 169), "parahippocampus": (2, 68)}}
rows = []
for label in ("boundary", "event"):
for region in cb.REGIONS:
in_region = cells.region == region
found, total = int(cells[label][in_region].sum()), int(in_region.sum())
null = chance[f"{label} {region}"]
k, n = paper_regions[label][region]
rows.append({"class": label, "region": region, "units": total, "cells": found,
"%": round(100 * found / total, 2),
"chance mean (cells)": round(null.mean(), 2),
"chance 95th percentile": int(np.percentile(null, 95)),
"paper": f"{k}/{n} ({100 * k / n:.1f}%)"})
regions = pd.DataFrame(rows).set_index(["class", "region"])
print(regions.to_string())
print("\nSame, under the pooled reading:")
pooled = variants["pooled reading (three tests), authors' SB trials"]
print(pooled.groupby("region")[["boundary", "event"]].agg(["sum", "size"]).reindex(cb.REGIONS).to_string())
fig, ax = plt.subplots(figsize=(5.4, 3.2))
x = np.arange(len(regions))
ax.bar(x - 0.2, regions["%"], 0.4, color="0.3", label="here")
ax.bar(x + 0.2, [100 * k / n for label in paper_regions for k, n in paper_regions[label].values()], 0.4,
color="0.75", label="paper")
ax.set_xticks(x, [f"{c}\n{r[:5]}." for c, r in regions.index], fontsize=8)
ax.set_ylabel("% of units in the region"); ax.legend(frameon=False)
plt.tight_layout()
print(f"\nPlotted percentages here range from {regions['%'].min():.2f} to {regions['%'].max():.2f}.")
units cells % chance mean (cells) chance 95th percentile paper
class region
boundary hippocampus 267 1 0.37 0.04 0 12/343 (3.5%)
amygdala 229 0 0.00 0.08 1 12/169 (7.1%)
parahippocampus 53 0 0.00 0.02 0 18/68 (26.5%)
event hippocampus 267 2 0.75 0.51 2 27/343 (7.9%)
amygdala 229 3 1.31 0.52 2 7/169 (4.1%)
parahippocampus 53 2 3.77 0.19 1 2/68 (2.9%)
Same, under the pooled reading:
boundary event
sum size sum size
region
hippocampus 6 267 3 267
amygdala 6 229 4 229
parahippocampus 3 53 3 53
Plotted percentages here range from 0.00 to 3.77.
Result: blocked, methods underspecified. This claim depends on the cells selected under
fig3c-boundary-cells and fig3g-event-cells, whose selection procedure cannot be recovered.
With the literal definition, none of the 53 parahippocampal units is a boundary cell (paper: 18
of 68); the single boundary cell is in the hippocampus. Event cells are 2 of 267 (0.75%) in the
hippocampus, within what label shuffling gives there (95th percentile 2), against 3 of 229 in
the amygdala and 2 of 53 in the parahippocampal gyrus. Under the pooled reading boundary cells are
6 of 267 (2.2%), 6 of 229 (2.6%) and 3 of 53 (5.7%) in hippocampus, amygdala and
parahippocampal gyrus, the paper's ordering, and event cells 3 of 267, 4 of 229 and 3 of 53,
with the parahippocampal gyrus highest where the paper has the hippocampus. No chance level by
region was computed for the pooled reading, and the counts are small. The
dandiset's split of units by region also differs from the paper's.
Claim edfig5-clip-onset-cells¶
"Seventy six of 580 (13.1%, p = 0.01) and four of 580 (0.7%, p = 0.43) MTL cells changed their firing rate in response to clip onsets and clip offsets, respectively, but neither of these cells qualified as boundary cells or event cells." A cell responds if its firing differs between the second before and the second after the onset (or offset) over all 90 clips, "p < 0.05, one-tailed permutation t-test". The direction of the tail is not stated; an increase is tested, and the other choices are shown. The second before a clip's onset is the fixation period.
rows = []
for s in good:
for event, column in (("onset", "start_time"), ("offset", "stop_time")):
t = s.encoding[column].to_numpy()
pre, post = cb.spike_counts(s.units, t, -cb.WINDOW, 0.0), cb.spike_counts(s.units, t, 0.0, cb.WINDOW)
rows.append(pd.DataFrame({
"subject": s.subject, "unit": np.arange(len(s.units)), "region": s.region, "event": event,
"increase": cb.paired_permutation_p(pre, post, rng, tail="greater") < cb.ALPHA,
"decrease": cb.paired_permutation_p(post, pre, rng, tail="greater") < cb.ALPHA,
"two-sided": cb.paired_permutation_p(pre, post, rng) < cb.ALPHA}))
responses = pd.concat(rows, ignore_index=True)
summary = responses.groupby("event")[["increase", "decrease", "two-sided"]].sum().loc[["onset", "offset"]]
print(f"Units responding, of {n_units}:")
print(summary.to_string())
print("As a percentage:")
print((100 * summary / n_units).round(1).to_string())
print("Paper: onset 76 of 580 (13.1%), offset 4 of 580 (0.7%).")
for event in ("onset", "offset"):
k = int(summary.loc[event, "increase"])
print(f" {event}, increase: {k} of {n_units}; probability of at least {k} if 5% respond by chance: "
f"{stats.binom.sf(k - 1, n_units, 0.05):.1e}; probability of at most {k}: {stats.binom.cdf(k, n_units, 0.05):.2f}")
onset = responses[responses.event == "onset"].reset_index(drop=True)
offset = responses[responses.event == "offset"].reset_index(drop=True)
print("\nOnset-responsive units (increase) by region:")
print(onset.groupby("region").increase.agg(["sum", "size", "mean"]).reindex(cb.REGIONS).round(3).to_string())
print("By subject (responsive / units):",
{s: f"{int(g.increase.sum())}/{len(g)}" for s, g in onset.groupby("subject", sort=False)})
share = onset.groupby("subject", sort=False).increase.agg(["sum", "size"])
high = share[share["sum"] / share["size"] >= 0.2]
print(f"Sessions in which at least 20% of units respond: {list(high.index)}. They hold {int(high['sum'].sum())} of the "
f"{int(share['sum'].sum())} onset-responsive units and {int(high['size'].sum())} of the {n_units} units.")
k = int(summary.loc["onset", "increase"])
print(f"Probability of at least {k} of {n_units} if the true proportion were the paper's 0.131: "
f"{stats.binom.sf(k - 1, n_units, 0.131):.1e}")
either = onset.increase | offset.increase
print(f"\nOverlap with the cells found above: {int((either & cells.boundary).sum())} of "
f"{int(cells.boundary.sum())} boundary cells and {int((either & cells.event).sum())} of "
f"{int(cells.event.sum())} event cells respond to clip onsets or offsets.")
Units responding, of 549:
increase decrease two-sided
event
onset 98 17 98
offset 26 28 23
As a percentage:
increase decrease two-sided
event
onset 17.9 3.1 17.9
offset 4.7 5.1 4.2
Paper: onset 76 of 580 (13.1%), offset 4 of 580 (0.7%).
onset, increase: 98 of 549; probability of at least 98 if 5% respond by chance: 1.2e-27; probability of at most 98: 1.00
offset, increase: 26 of 549; probability of at least 26 if 5% respond by chance: 6.4e-01; probability of at most 26: 0.44
Onset-responsive units (increase) by region:
sum size mean
region
hippocampus 57 267 0.213
amygdala 35 229 0.153
parahippocampus 6 53 0.113
By subject (responsive / units): {'1': '3/43', '2': '12/45', '3': '1/31', '4': '2/35', '5': '2/30', '6': '1/19', '7': '0/34', '8': '1/33', '9': '4/31', '10': '1/27', '11': '7/14', '12': '1/28', '13': '0/48', '14': '13/30', '15': '12/31', '16': '12/14', '17': '12/27', '18': '14/29'}
Sessions in which at least 20% of units respond: ['2', '11', '14', '15', '16', '17', '18']. They hold 82 of the 98 onset-responsive units and 190 of the 549 units.
Probability of at least 98 of 549 if the true proportion were the paper's 0.131: 9.4e-04
Overlap with the cells found above: 0 of 1 boundary cells and 0 of 7 event cells respond to clip onsets or offsets.
Result: reproduced. 98 of 549 units (17.9%) fire more in the second after a clip starts than in the second before, against 13.1% in the paper; the gap of 4.8 points is just inside the tolerance of 5 set for this claim, and 98 is more than a true proportion of 13.1% would plausibly give (p = 9.4e-04). 26 units (4.7%) fire more after a clip ends, which is the 5% that chance gives; the paper found 0.7%. The onset responses are concentrated: the seven sessions in which at least 20% of units respond hold 82 of the 98; the other eleven together have 16 of 359 (4.5%), about what chance gives. The overall proportion is therefore an average over sessions that differ widely, and it falls inside the tolerance by 0.2 points. The data check above shows an aligned spike clock for ten of those eleven sessions; for subject 3 (1 of 31 units responsive) it could not be established. None of the 1 boundary cell and 7 event cells found here responds to clip onsets or offsets, which agrees with the paper, but 8 cells are too few to make much of it.
Claim fig5d-state-shift¶
"During NB video clips, the neural state exhibited only slow changes as a function of time. In contrast, the neural state changed abruptly following SBs and HBs ... Plotting MDD as a function of time revealed an abrupt change within ~300 ms after the boundary for SB and HB video clips."
From the Methods: each neuron's spikes in each trial are counted in 10 ms bins, smoothed with a Gaussian kernel of 200 ms standard deviation and z-scored over the trial; all medial temporal lobe neurons of all subjects form a pseudopopulation; principal components are computed, and the multidimensional distance (MDD) is the Euclidean distance, in the space of the components that explain 99% of the variance, between the state at a time and the state at the boundary.
What the Methods leave open, and what is done here:
- Matching trials across subjects. Every subject saw the same 90 clips in a different order. Pseudo-trials are built by clip name, so a pseudo-trial is one clip seen by all neurons.
- What the components are fitted on. The paper needed 55 components for 99% of the variance of 580 neurons, which points to condition averages and not single trials. Components are fitted here on the three condition averages from 1 s before to 1.5 s after the boundary. Single pseudo-trials are then projected on the same components.
- The trial for z-scoring is taken as 1 s before to 8 s after clip onset.
- SB clips are aligned to their first soft boundary, one pseudo-trial per clip. The authors' code does not cover this analysis.
The distance from the state at the boundary grows with time in any trajectory, also before the boundary and also without one. The quantities compared are therefore the distance in the second after the boundary against clips with no boundary, and against the second before the boundary.
rng = np.random.default_rng(5) # its own seed, so these tests do not depend on earlier cells
per_session = [cb.smoothed_trials(s) for s in good]
clips, kinds, data = cb.pseudopopulation(per_session)
shift = cb.distance_from_boundary(data, kinds)
times = cb.state_times()
after, before, early = (times > 0) & (times < 1), times < 0, (times > 0) & (times < 0.3)
print(f"Pseudopopulation: {data.shape[1]} units, {data.shape[0]} pseudo-trials "
f"({', '.join(f'{(kinds == k).sum()} {k}' for k in cb.KINDS)}), {data.shape[2]} bins of {cb.BIN * 1000:.0f} ms.")
print(f"Components for {100 * cb.VARIANCE:.0f}% of the variance of the condition averages: {shift['components']} "
f"(paper: 55 for 580 units); the first three explain {100 * shift['explained by first three']:.1f}% "
f"(paper: 43.6%).")
def shift_table(result, trial_kinds):
rows = {}
post = result["trials"][:, after].mean(axis=1)
pre = result["trials"][:, before].mean(axis=1)
for k in cb.KINDS:
m = trial_kinds == k
average = result["average"][k]
rows[k] = {
"average trajectory: 0 to 0.3 s": average[early].mean(),
"average trajectory: 0 to 1 s": average[after].mean(),
"average trajectory: -1 to 0 s": average[before].mean(),
"single trials: 0 to 1 s": post[m].mean(),
" s.e.m.": post[m].std(ddof=1) / np.sqrt(m.sum()),
"single trials: -1 to 0 s": pre[m].mean(),
"single trials: after minus before": (post - pre)[m].mean(),
" s.e.m.": (post - pre)[m].std(ddof=1) / np.sqrt(m.sum()),
}
return pd.DataFrame(rows).round(2), post, pre
table, post, pre = shift_table(shift, kinds)
print("\nDistance from the state at the boundary (units of z-scored activity):")
print(table.to_string())
print("\nOne-sided permutation tests over pseudo-trials (10000 permutations):")
for k in ("SB", "HB"):
p_post = cb.two_sample_permutation_p(post[None, kinds == k], post[None, kinds == "NB"], rng, 10000)[0]
p_diff = cb.two_sample_permutation_p((post - pre)[None, kinds == k], (post - pre)[None, kinds == "NB"], rng, 10000)[0]
print(f" {k} > NB, distance after the boundary: p = {p_post:.4f}; "
f"distance after minus before: p = {p_diff:.4f}")
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4))
for k in cb.KINDS:
axes[0].plot(times, shift["average"][k], color=COLOR[k], label=k)
trials = shift["trials"][kinds == k]
mean, sem = trials.mean(axis=0), trials.std(axis=0, ddof=1) / np.sqrt(len(trials))
axes[1].plot(times, mean, color=COLOR[k], label=k)
axes[1].fill_between(times, mean - sem, mean + sem, color=COLOR[k], alpha=0.2, lw=0)
axes[0].set_title("Condition-average trajectory", fontsize=10)
axes[1].set_title("Single pseudo-trials, mean ± s.e.m.", fontsize=10)
for ax in axes:
ax.axvline(0, color="0.6", lw=1); ax.set_xlabel("Time from boundary (s)"); ax.set_ylabel("MDD")
axes[0].legend(frameon=False)
plt.tight_layout()
print("\nPlotted distances range from 0 to %.1f (average trajectories) and from %.1f to %.1f (single-trial means)." % (
max(shift["average"][k].max() for k in cb.KINDS),
min(shift["trials"][kinds == k].mean(axis=0).min() for k in cb.KINDS),
max(shift["trials"][kinds == k].mean(axis=0).max() for k in cb.KINDS)))
print("Time of the largest distance of the average trajectory in the 1.5 s after the boundary (s):",
{k: round(float(times[times > 0][np.argmax(shift["average"][k][times > 0])]), 2) for k in cb.KINDS})
Pseudopopulation: 549 units, 90 pseudo-trials (30 NB, 30 SB, 30 HB), 250 bins of 10 ms.
Components for 99% of the variance of the condition averages: 22 (paper: 55 for 580 units); the first three explain 50.4% (paper: 43.6%).
Distance from the state at the boundary (units of z-scored activity):
NB SB HB
average trajectory: 0 to 0.3 s 2.21 2.96 3.81
average trajectory: 0 to 1 s 5.09 5.48 7.02
average trajectory: -1 to 0 s 5.35 5.76 6.06
single trials: 0 to 1 s 7.42 8.15 9.03
s.e.m. 0.17 0.18 0.33
single trials: -1 to 0 s 7.39 8.15 8.19
single trials: after minus before 0.04 0.00 0.84
s.e.m. 0.23 0.21 0.19
One-sided permutation tests over pseudo-trials (10000 permutations):
SB > NB, distance after the boundary: p = 0.0030; distance after minus before: p = 0.5494
HB > NB, distance after the boundary: p = 0.0001; distance after minus before: p = 0.0037
Plotted distances range from 0 to 9.9 (average trajectories) and from 0.0 to 12.4 (single-trial means).
Time of the largest distance of the average trajectory in the 1.5 s after the boundary (s): {'NB': 1.5, 'SB': 0.74, 'HB': 1.5}
How the result depends on the choices, and on the subject. For each subject the same analysis is run on that subject's units alone (sessions with fewer than 10 units are left out), as in the paper's Extended Data Figure 8.
rows = []
for variance in (0.90, 0.95, 0.99):
result = cb.distance_from_boundary(data, kinds, variance)
_, v_post, v_pre = shift_table(result, kinds)
row = {"variance kept": variance, "components": result["components"]}
for k in cb.KINDS:
row[f"{k} after - before"] = round(float((v_post - v_pre)[kinds == k].mean()), 2)
for k in ("SB", "HB"):
row[f"p, {k} > NB"] = round(float(cb.two_sample_permutation_p(
(v_post - v_pre)[None, kinds == k], (v_post - v_pre)[None, kinds == "NB"], rng, 10000)[0]), 4)
rows.append(row)
print("Single-trial distance after minus before the boundary, by the variance kept:")
print(pd.DataFrame(rows).set_index("variance kept").to_string())
neither = ~(cells.boundary | cells.event).to_numpy()
result = cb.distance_from_boundary(data[:, neither], kinds)
_, v_post, v_pre = shift_table(result, kinds)
print(f"\nWithout the {int((~neither).sum())} boundary and event cells found above ({int(neither.sum())} units): "
"after minus before =", {k: round(float((v_post - v_pre)[kinds == k].mean()), 2) for k in cb.KINDS})
rows = []
for s, one in zip(good, per_session):
if len(s.units) < 10:
continue
s_clips, s_kinds, s_data = cb.pseudopopulation([one])
result = cb.distance_from_boundary(s_data, s_kinds)
_, s_post, s_pre = shift_table(result, s_kinds)
row = {"subject": s.subject, "units": len(s.units), "components": result["components"]}
for k in cb.KINDS:
row[f"{k} after - before"] = round(float((s_post - s_pre)[s_kinds == k].mean()), 2)
rows.append(row)
by_session = pd.DataFrame(rows).set_index("subject")
print(f"\nBy subject ({len(by_session)} sessions with at least 10 units; {len(good) - len(by_session)} left out):")
print(by_session.T.to_string())
for k in ("SB", "HB"):
diff = by_session[f"{k} after - before"] - by_session["NB after - before"]
w = stats.wilcoxon(diff)
print(f" {k} minus NB: larger in {(diff > 0).sum()} of {len(diff)} subjects, median {diff.median():+.2f}, "
f"Wilcoxon signed-rank p = {w.pvalue:.3f}")
# SB pseudo-trials are aligned to the clip's first soft boundary. Where that comes less than
# 1 s after clip onset, the second "before" the boundary reaches back before the clip started.
first_sb = good[0].encoding.set_index("clip").b1
same = all(np.allclose(s.encoding.set_index("clip").b1.loc[first_sb.index], first_sb) for s in good)
early_sb = first_sb[(first_sb.index.str[:2] == "SB") & (first_sb < 1.0)]
print(f"\nSB clips whose first soft boundary is less than 1 s after clip onset: {len(early_sb)} of 30 "
f"(times {early_sb.min():.2f} to {early_sb.max():.2f} s; same in every session: {same}).")
keep = ~np.isin(clips, early_sb.index)
result = cb.distance_from_boundary(data[keep], kinds[keep])
_, k_post, k_pre = shift_table(result, kinds[keep])
k_kinds = kinds[keep]
print("Without them: after minus before =",
{k: round(float((k_post - k_pre)[k_kinds == k].mean()), 2) for k in cb.KINDS})
for k in ("SB", "HB"):
p_post = cb.two_sample_permutation_p(k_post[None, k_kinds == k], k_post[None, k_kinds == "NB"], rng, 10000)[0]
p_diff = cb.two_sample_permutation_p((k_post - k_pre)[None, k_kinds == k],
(k_post - k_pre)[None, k_kinds == "NB"], rng, 10000)[0]
print(f" {k} > NB, distance after the boundary: p = {p_post:.4f}; after minus before: p = {p_diff:.4f}")
Single-trial distance after minus before the boundary, by the variance kept:
components NB after - before SB after - before HB after - before p, SB > NB p, HB > NB
variance kept
0.90 13 -0.04 -0.20 0.94 0.6939 0.0008
0.95 16 -0.04 -0.07 0.87 0.5310 0.0021
0.99 22 0.04 0.00 0.84 0.5498 0.0035
Without the 8 boundary and event cells found above (541 units): after minus before = {'NB': 0.03, 'SB': 0.02, 'HB': 0.79}
By subject (18 sessions with at least 10 units; 0 left out):
subject 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
units 43.00 45.00 31.00 35.00 30.00 19.00 34.00 33.00 31.00 27.00 14.00 28.00 48.00 30.00 31.00 14.00 27.00 29.00
components 20.00 21.00 18.00 19.00 18.00 15.00 18.00 19.00 17.00 17.00 10.00 17.00 20.00 17.00 16.00 10.00 16.00 17.00
NB after - before -0.13 0.17 0.17 -0.10 0.05 -0.47 0.05 -0.09 0.31 0.12 -0.15 -0.03 -0.41 0.06 0.05 0.19 0.07 0.02
SB after - before 0.40 0.29 0.32 0.10 -0.13 0.14 -0.20 0.39 0.13 0.26 0.01 0.07 -0.01 -0.12 0.09 -0.07 0.45 0.06
HB after - before 0.11 0.17 0.17 0.11 0.03 -0.11 -0.09 0.26 0.03 0.07 0.13 0.22 -0.30 0.12 0.23 0.33 0.25 0.07
SB minus NB: larger in 13 of 18 subjects, median +0.13, Wilcoxon signed-rank p = 0.144
HB minus NB: larger in 12 of 18 subjects, median +0.12, Wilcoxon signed-rank p = 0.020
SB clips whose first soft boundary is less than 1 s after clip onset: 2 of 30 (times 0.55 to 0.98 s; same in every session: True).
Without them: after minus before = {'NB': 0.04, 'SB': -0.03, 'HB': 0.84}
SB > NB, distance after the boundary: p = 0.0013; after minus before: p = 0.5856
HB > NB, distance after the boundary: p = 0.0001; after minus before: p = 0.0034
Result: partially reproduced. After a hard boundary the population state moves further from its state at the boundary than it does in clips with no boundary (single-trial distance in the second after: 9.03 against 7.42, p = 0.0001), and further than in the second before (after minus before: 0.84 against 0.04, p = 0.004). This holds when 90%, 95% or 99% of the variance is kept (p ≤ 0.004), without the 8 boundary and event cells found above (0.79), and in 12 of 18 subjects analysed alone (Wilcoxon p = 0.020).
After a soft boundary the distance is also larger than in clips with no boundary (8.15 against 7.42, p = 0.003), but it is as large in the second before the boundary (8.15), and after minus before is 0.00, no larger than without a boundary (p = 0.55; larger in 13 of 18 subjects, p = 0.14). So a shift that is specific to the time after the boundary is found for hard boundaries and not for soft ones. Two SB clips have their first soft boundary less than 1 s after clip onset, so that their second before the boundary reaches back before the clip started; without them, after minus before is -0.03 for SB (p = 0.59) and 0.84 for HB (p = 0.003).
That the shift is abrupt and complete within about 300 ms is not shown by these numbers. On the condition-average trajectory the distance in the first 300 ms is 2.21 for NB, 2.96 for SB and 3.81 for HB, and the distance keeps growing to the end of the 1.5 s window for NB and HB, and peaks at 0.74 s for SB.
Summary¶
| Claim | Outcome | Note |
|---|---|---|
fig2a-order-accuracy |
Reproduced | 0.53 across a hard boundary, 0.73 and 0.72 otherwise |
fig2b-order-reaction-time |
Reproduced | 2.05 s across a hard boundary, 1.62 and 1.54 s otherwise |
fig2c-order-confidence |
Reproduced | 1.98 across a hard boundary, 2.45 and 2.51 otherwise |
edfig2a-recognition-similar |
Reproduced | 0.74, 0.74 and 0.72; no difference (p = 0.69) |
fig2def-recognition-vs-boundary-distance |
Blocked: metadata insufficient | The position of a tested frame in its clip is not given |
fig3c-boundary-cells |
Blocked: methods underspecified | 1 of 549 units (0.18%) against 42 of 580 (7.24%); no reading reaches the paper's chance level |
fig3g-event-cells |
Blocked: methods underspecified | 7 of 549 units (1.28%), above chance, against 36 of 580 (6.20%) |
fig3-cell-regions |
Blocked: methods underspecified | Depends on the two cell selections above |
fig3dh-no-onset-response |
Blocked: metadata insufficient | The paper's cells are not labelled; 1 and 7 cells found here, below the 10 needed |
edfig5-clip-onset-cells |
Reproduced | 17.9% respond to clip onsets (paper 13.1%), 4.7% to offsets |
fig3j-response-latency |
Blocked: metadata insufficient | As fig3dh-no-onset-response |
fig3-outside-mtl |
Blocked: data missing | No units outside the medial temporal lobe |
fig4c-boundary-cells-predict-recognition |
Blocked: metadata insufficient | SB and HB test trials cannot be linked to clips |
fig4h-event-cells-phase-predicts-order |
Blocked: data missing | No local field potentials |
fig5d-state-shift |
Partially reproduced | A shift after hard boundaries, none specific to soft boundaries |
fig5ij-state-shift-predicts-memory |
Blocked: metadata insufficient | SB and HB test trials cannot be linked to clips |
fig6-reinstatement |
Blocked: metadata insufficient | SB and HB test trials cannot be linked to clips |
fig6i-reinstatement-timing |
Blocked: metadata insufficient | As above, and the position of the tested frame is not given |
fig6jm-reinstatement-confidence |
Blocked: metadata insufficient | SB and HB test trials cannot be linked to clips |
fig6no-false-alarm-similarity |
Blocked: data missing | No similarity ratings |
Of 20 claims, 6 have a tested outcome: 5 reproduced and 1 partially reproduced. 14 are blocked, three of them after the analysis was run: the cell counts were computed, but the procedure behind the paper's counts could not be determined.
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How this analysis differs from the paper's¶
- Subjects and units. Behaviour is from 19 of the paper's 20 subjects. Neurons are 549 units from 18 sessions, against the paper's 580 from 19; the 16 units of subject 19 are left out because their spike times are not on the clock of the trial tables.
- Text and code. The author manuscript was used. The supplementary tables and figures were not available. The laboratory's published code was read; it only plots rasters, and the one thing taken from it is which soft boundaries count as SB trials.
- Permutation tests. The paper's "permutation t-test" is implemented as a sign-flip test for before against after, and a label-permutation test between conditions, with 1000 permutations.
- Soft boundaries. For the cell definitions, every soft boundary of a session is an SB trial (75), as in the authors' plotting code. For the state space, SB clips are aligned to their first soft boundary.
- Chance level. 300 label shuffles instead of 1000.
- Clip-onset cells. The one-tailed test is taken as testing an increase.
- State space. Pseudo-trials are matched by clip name; principal components are fitted on the three condition averages from 1 s before to 1.5 s after the boundary, which are averages of the same trials that are then compared; a trial is z-scored from 1 s before to 8 s after clip onset. The Methods state none of these.
- Tolerances on reported numbers are this registry's; the paper gives none.
Which outcomes depend on the reading¶
The cell counts. fig3c-boundary-cells gives 0 to 2 cells under every five-test variant and 15
(2.73%), above chance, under the pooled reading. fig3g-event-cells gives 2 to 12 cells across
the variants, always far below the paper's proportion, and above chance in most. The chance
levels computed, for four of the readings, are 0.01% to 0.72%, far below the paper's (2.11% and
2.26%), which is why these claims are blocked. The behavioural claims, the clip-onset cells and the state shift do not use the cell
definitions.
What was found about the dandiset¶
The authors' published code does not contain the paper's analyses: no cell selection, chance level, statistics, state-space or reinstatement analysis, and nothing that links a tested frame to its clip. It reads the NWB files and nothing else.
That code drops missing boundaries by testing for a time above 0. In the pinned version a missing boundary is stored as the clip's start time, so nothing would be dropped (150 of 150 kept, where 75 exist) and the code's fixed row ranges for the three conditions would be wrong. In the most recent version missing boundaries are NaN and the code works.
The Methods describe firing-rate traces with a 200 ms window moved in 2 ms steps; the code uses 10 ms bins smoothed with a Gaussian window of 20 bins.
The most recent version (0.230530.1822) has the spike times of subjects 13, 17 and 18 shifted against the encoding trials by 191.5, -16.8 and 9.3 s. The version the paper cites (0.220216.0323) does not have this problem.
In every version read, the spike times of subject 19 are near -7.02e12 s and cannot be used.
Scene recognition and time discrimination trials with soft or hard boundaries cannot be linked to the clips they test, because their frame names are numbered differently from the clips. The position of a tested frame within its clip is not given either.
The dandiset has no local field potentials, no units outside the medial temporal lobe, no labels for boundary and event cells, no similarity ratings and no stimuli, and it lacks the behaviour of one of the paper's 20 subjects.
It has 565 medial temporal lobe units where the paper has 580, and a different split by region.
Several column descriptions do not match the values: the codes of
stimCategory, the value for a missing boundary, the code for "no response" inaccuracyandresp_value, and whatRTis measured from. Twenty time discrimination trials without a response are coded correct.The
eventsseries holds entries with an undocumented code 0 and meaningless timestamps.