Daume et al. (2024) on dandiset 000673: behaviour and single neurons¶
Paper. Daume J, Kaminski J, Schjetnan AGP, Salimpour Y, Khan U, Kyzar M, Reed CM, Anderson WS, Valiante TA, Mamelak AN, Rutishauser U (2024). Control of working memory by phase-amplitude coupling of human hippocampal neurons. Nature 629, 393-401. doi:10.1038/s41586-024-07309-z.
Dandiset. DANDI:000673, version 0.240118.2135. The paper's reference 59 cites this version ("v.0.240118.2135"), so it is pinned, as the protocol asks. Four later versions exist (0.240524.1751, 0.240524.1758, 0.250122.0110 and the draft); in 0.250122.0110 every one of the 44 files has a new asset with a size a few hundred kilobytes smaller, with the same paths. The analysis here does not use the later versions.
Text the claims were written from. The version of record (PMC11078732), with the Methods, figure legends, Extended Data legends and Supplementary Information, and the paper's code (github.com/rutishauserlab/SBCAT-release-NWB at commit a25e8cc3).
This notebook tests the claims on behaviour and single-neuron firing. The LFP claims are in
reproduce_daume_2024_lfp.ipynb and reproduce_daume_2024_lfp_controls.ipynb. Claims tested here:
| Claim | Tested here |
|---|---|
| cohort-counts | yes |
| fig1d-rt-load | yes |
| fig1d-accuracy-load | yes |
| behaviour-overall-accuracy | yes |
| category-cell-counts | yes |
| fig3b-persistent-activity | yes |
| edf4b-persistent-activity-by-area | yes |
| fig3c-category-load | yes |
| fig3d-category-correct-incorrect | yes |
| category-no-rt-effect | yes |
| edf8fgh-frontal-controls | panel h only; blocked (panels f and g need SFC, not run) |
| PAC neuron claims | in reproduce_daume_2024_pac_neurons.ipynb |
| fig3h-pac-category-firing | in reproduce_daume_2024_lfp.ipynb |
| SFC, noise correlation and decoding claims | blocked (see the attempt record) |
In [1]:
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import sbcat
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
assets = sbcat.list_assets()
S = sbcat.load_all(assets)
print(f"loaded {len(S)} files of dandiset {sbcat.DANDISET_ID} version {sbcat.DANDISET_VERSION} in {time.time() - started:.0f} s")
loaded 44 files of dandiset 000673 version 0.240118.2135 in 133 s
Data check¶
In [2]:
info = pd.DataFrame([{
"path": s["path"], "subject": s["subject"], "session": s["session"], "age": s["age"], "sex": s["sex"],
"trials": len(s["trials"]), "load1": int((s["trials"].load == 1).sum()), "load3": int((s["trials"].load == 3).sum()),
"units": len(s["units"]), "lfp_channels": len(s["lfp_electrodes"]), "electrodes": len(s["electrodes"]),
"lfp_rate": s["lfp_info"]["rate"] if s["lfp_info"] else np.nan,
"lfp_unit": s["lfp_info"]["unit"] if s["lfp_info"] else None} for s in S])
print(info.to_string(index=False))
print("\nsessions:", len(info), "; patients:", info.subject.nunique(),
"; sessions per patient:", info.subject.value_counts().value_counts().sort_index().to_dict())
pat = info.drop_duplicates("subject")
ages = pat.age.str.strip("PY").astype(int)
print("patients by sex:", pat.sex.value_counts().to_dict(), f"; age {ages.mean():.2f} +/- {ages.std(ddof=1):.2f} years")
print("trials per session:", info.trials.value_counts().to_dict())
print("files without LFP:", info.path[info.lfp_channels == 0].tolist())
print("LFP description:", {s["lfp_info"]["description"] for s in S if s["lfp_info"]})
path subject session age sex trials load1 load3 units lfp_channels electrodes lfp_rate lfp_unit
sub-1/sub-1_ses-1_ecephys+image.nwb 1 1 P43Y F 140 70 70 46 70 74 400.0 microvolts
sub-1/sub-1_ses-2_ecephys+image.nwb 1 2 P43Y F 140 70 70 41 69 73 400.0 microvolts
sub-10/sub-10_ses-1_ecephys+image.nwb 10 1 P54Y M 140 70 70 36 44 54 400.0 microvolts
sub-11/sub-11_ses-1_ecephys+image.nwb 11 1 P41Y F 139 70 69 36 65 67 400.0 microvolts
sub-12/sub-12_ses-1_ecephys+image.nwb 12 1 P30Y F 140 70 70 24 55 64 400.0 microvolts
sub-12/sub-12_ses-2_ecephys+image.nwb 12 2 P30Y F 140 70 70 12 53 56 400.0 microvolts
sub-13/sub-13_ses-1_ecephys+image.nwb 13 1 P40Y M 140 70 70 27 70 70 400.0 microvolts
sub-14/sub-14_ses-1_ecephys+image.nwb 14 1 P25Y F 140 70 70 13 49 52 400.0 microvolts
sub-15/sub-15_ses-1_ecephys+image.nwb 15 1 P58Y F 140 70 70 41 71 71 400.0 microvolts
sub-16/sub-16_ses-1_ecephys+image.nwb 16 1 P24Y F 140 70 70 63 64 64 400.0 microvolts
sub-17/sub-17_ses-1_ecephys+image.nwb 17 1 P46Y F 140 70 70 80 57 60 400.0 microvolts
sub-18/sub-18_ses-1_ecephys+image.nwb 18 1 P54Y F 140 70 70 27 14 14 400.0 microvolts
sub-19/sub-19_ses-1_ecephys+image.nwb 19 1 P42Y F 140 70 70 95 59 72 400.0 microvolts
sub-19/sub-19_ses-2_ecephys+image.nwb 19 2 P42Y F 140 70 70 72 58 71 400.0 microvolts
sub-2/sub-2_ses-1_ecephys+image.nwb 2 1 P48Y M 140 70 70 42 52 54 400.0 microvolts
sub-20/sub-20_ses-1_ecephys+image.nwb 20 1 P26Y M 140 70 70 28 0 14 NaN NaN
sub-21/sub-21_ses-1_ecephys+image.nwb 21 1 P45Y F 140 70 70 40 56 59 400.0 microvolts
sub-22/sub-22_ses-1_ecephys+image.nwb 22 1 P20Y M 128 63 65 42 24 33 400.0 microvolts
sub-22/sub-22_ses-2_ecephys+image.nwb 22 2 P20Y M 140 70 70 13 38 38 400.0 microvolts
sub-22/sub-22_ses-3_ecephys+image.nwb 22 3 P20Y M 140 70 70 13 38 40 400.0 microvolts
sub-23/sub-23_ses-1_ecephys+image.nwb 23 1 P59Y F 140 70 70 18 33 33 400.0 microvolts
sub-24/sub-24_ses-1_ecephys+image.nwb 24 1 P23Y M 140 70 70 7 14 21 400.0 microvolts
sub-25/sub-25_ses-1_ecephys+image.nwb 25 1 P58Y F 140 70 70 11 14 18 400.0 microvolts
sub-26/sub-26_ses-1_ecephys+image.nwb 26 1 P25Y F 140 70 70 21 36 37 400.0 microvolts
sub-26/sub-26_ses-2_ecephys+image.nwb 26 2 P25Y F 140 70 70 14 35 35 400.0 microvolts
sub-27/sub-27_ses-1_ecephys+image.nwb 27 1 P49Y M 140 70 70 16 29 29 400.0 microvolts
sub-28/sub-28_ses-1_ecephys+image.nwb 28 1 P26Y M 140 70 70 21 28 36 400.0 microvolts
sub-29/sub-29_ses-1_ecephys+image.nwb 29 1 P28Y M 140 70 70 31 62 62 400.0 microvolts
sub-3/sub-3_ses-1_ecephys+image.nwb 3 1 P32Y F 140 70 70 53 69 73 400.0 microvolts
sub-30/sub-30_ses-1_ecephys+image.nwb 30 1 P38Y M 140 70 70 16 33 34 400.0 microvolts
sub-31/sub-31_ses-1_ecephys+image.nwb 31 1 P36Y M 140 70 70 32 42 42 400.0 microvolts
sub-31/sub-31_ses-2_ecephys+image.nwb 31 2 P36Y M 140 70 70 18 35 35 400.0 microvolts
sub-32/sub-32_ses-1_ecephys+image.nwb 32 1 P28Y M 140 70 70 36 47 49 400.0 microvolts
sub-33/sub-33_ses-1_ecephys+image.nwb 33 1 P25Y F 140 70 70 30 42 42 400.0 microvolts
sub-34/sub-34_ses-1_ecephys+image.nwb 34 1 P62Y M 140 70 70 12 16 16 400.0 microvolts
sub-35/sub-35_ses-1_ecephys+image.nwb 35 1 P45Y M 140 70 70 1 8 8 400.0 microvolts
sub-36/sub-36_ses-1_ecephys+image.nwb 36 1 P27Y F 140 70 70 10 8 9 400.0 microvolts
sub-4/sub-4_ses-1_ecephys+image.nwb 4 1 P67Y M 140 70 70 51 69 72 400.0 microvolts
sub-5/sub-5_ses-1_ecephys+image.nwb 5 1 P52Y F 140 70 70 39 63 63 400.0 microvolts
sub-5/sub-5_ses-2_ecephys+image.nwb 5 2 P52Y F 140 70 70 36 53 54 400.0 microvolts
sub-6/sub-6_ses-1_ecephys+image.nwb 6 1 P25Y F 140 70 70 24 27 27 400.0 microvolts
sub-7/sub-7_ses-1_ecephys+image.nwb 7 1 P63Y F 140 70 70 112 52 53 400.0 microvolts
sub-8/sub-8_ses-1_ecephys+image.nwb 8 1 P55Y F 140 70 70 33 42 57 400.0 microvolts
sub-9/sub-9_ses-1_ecephys+image.nwb 9 1 P38Y F 140 70 70 21 59 63 400.0 microvolts
sessions: 44 ; patients: 36 ; sessions per patient: {1: 29, 2: 6, 3: 1}
patients by sex: {'F': 21, 'M': 15} ; age 40.47 +/- 13.76 years
trials per session: {140: 42, 139: 1, 128: 1}
files without LFP: ['sub-20/sub-20_ses-1_ecephys+image.nwb']
LFP description: {'These are LFP recordings that have spike potentials removed and is downsampled to 400Hz'}
In [3]:
units = pd.concat([s["units"].assign(path=s["path"], subject=s["subject"]) for s in S], ignore_index=True)
lfp = pd.concat([s["electrodes"].iloc[s["lfp_electrodes"]].assign(path=s["path"]) for s in S], ignore_index=True)
print("electrode locations:", sorted(set(units.location) | set(lfp.location)))
print("\nunits by area and hemisphere:\n", units.pivot_table(index="area", columns="hemisphere", values="unit", aggfunc="count").to_string())
print("\nLFP channels by area and hemisphere:\n", lfp.pivot_table(index="area", columns="hemisphere", values="x", aggfunc="count").to_string())
# Each unit's electrode, and each LFP channel, is a row of the electrodes table; electrodes of one
# area and hemisphere in a session share coordinates (one bundle of 8 wires).
bundles = pd.concat([s["electrodes"].assign(path=s["path"]) for s in S]).groupby(["path", "location"]).agg(
wires=("x", "size"), distinct_xyz=("x", lambda v: len(set(v))))
print("\nwires per bundle (area and hemisphere of a session):", bundles.wires.value_counts().sort_index().to_dict(),
"; bundles with more than one coordinate:", int((bundles.distinct_xyz > 1).sum()), "of", len(bundles))
# Columns of the units table (for the waveform measure of Extended Data Fig. 8d,e)
ucols = pd.Series({s["path"]: tuple(sorted(s["units_columns"])) for s in S})
print("\nunits table columns, by the number of files that have them:", ucols.value_counts().to_dict())
wf = [s["units_columns"].get("waveform_mean") for s in S]
print("waveform_mean shapes (first 3 files):", wf[:3], "; files without it:", sum(w is None for w in wf))
electrode locations: ['amygdala_left', 'amygdala_right', 'dorsal_anterior_cingulate_cortex_left', 'dorsal_anterior_cingulate_cortex_right', 'hippocampus_left', 'hippocampus_right', 'pre_supplementary_motor_area_left', 'pre_supplementary_motor_area_right', 'ventral_medial_prefrontal_cortex_left', 'ventral_medial_prefrontal_cortex_right']
units by area and hemisphere:
hemisphere left right
area
amygdala 233 263
dACC 109 79
hippocampus 203 157
preSMA 110 94
vmPFC 96 110
LFP channels by area and hemisphere:
hemisphere left right
area
amygdala 194 227
dACC 166 159
hippocampus 290 296
preSMA 144 139
vmPFC 153 154
wires per bundle (area and hemisphere of a session): {1: 8, 2: 3, 3: 2, 4: 9, 5: 5, 6: 17, 7: 180, 8: 36, 9: 1, 10: 1, 11: 1, 13: 4, 14: 15, 15: 3} ; bundles with more than one coordinate: 26 of 285
units table columns, by the number of files that have them: {('clusterID_orig', 'electrodes', 'id', 'spike_times', 'spike_times_index', 'waveform_mean', 'waveform_sd', 'waveforms', 'waveforms_index', 'waveforms_index_index', 'waveforms_isolation_distance', 'waveforms_mean_proj_dist', 'waveforms_mean_snr', 'waveforms_peak_snr'): 43, ('clusterID_orig', 'electrodes', 'id', 'spike_times', 'spike_times_index', 'waveforms', 'waveforms_index', 'waveforms_index_index', 'waveforms_isolation_distance', 'waveforms_mean_proj_dist', 'waveforms_mean_snr', 'waveforms_peak_snr'): 1}
waveform_mean shapes (first 3 files): [(46, 256), (41, 256), (36, 256)] ; files without it: 1
In [4]:
tr_all = pd.concat([s["trials"].assign(path=s["path"]) for s in S], ignore_index=True)
d = tr_all
print("trial columns:", [c for c in S[0]["trials"].columns])
print("loads:", d.load.value_counts().to_dict(), "; response_accuracy:", d.response_accuracy.value_counts().to_dict(),
"; probe_in_out:", d.probe_in_out.value_counts().to_dict())
dur = pd.DataFrame({
"fixation_s": d.timestamps_Encoding1 - d.timestamps_FixationCross,
"picture1_s": d.timestamps_Encoding1_end - d.timestamps_Encoding1,
"maintenance_s": d.timestamps_Probe - d.timestamps_Maintenance,
"rt_s": d.rt})
print("\n", dur.describe().round(3).to_string())
l3 = d[d.load == 3]
gap = pd.concat([l3.timestamps_Encoding2 - l3.timestamps_Encoding1_end, l3.timestamps_Encoding3 - l3.timestamps_Encoding2_end])
print("\nload 3 blank between pictures (s): min %.3f, max %.3f (Methods: 0.017 to 0.200)" % (gap.min(), gap.max()))
print("maintenance onset equals the end of the last picture:",
bool(np.allclose(np.where(d.load == 1, d.timestamps_Encoding1_end, d.timestamps_Encoding3_end), d.timestamps_Maintenance)))
cats = pd.concat([d[f"cat_{k}"] for k in ["Encoding1", "Encoding2", "Encoding3", "Probe"]])
print("category codes (first digit of the picture ID):", cats.value_counts().sort_index().to_dict())
l3c = l3[["cat_Encoding1", "cat_Encoding2", "cat_Encoding3"]].to_numpy()
print("load 3 trials with three different categories:", int((np.array([len(set(r)) for r in l3c]) == 3).sum()), "of", len(l3c))
match = d.PicIDs_Probe.to_numpy()[:, None] == d[["PicIDs_Encoding1", "PicIDs_Encoding2", "PicIDs_Encoding3"]].to_numpy()
print("probe_in_out equals 'probe is one of the encoded pictures':", bool(np.all(match.any(1) == (d.probe_in_out == 1))))
print("RT above 10 s:", int((d.rt > 10).sum()), "trials; max %.1f s" % d.rt.max())
trial columns: ['PicIDs_Encoding1', 'PicIDs_Encoding2', 'PicIDs_Encoding3', 'PicIDs_Probe', 'loads', 'probe_in_out', 'response_accuracy', 'start_time', 'stop_time', 'timestamps_Encoding1', 'timestamps_Encoding1_end', 'timestamps_Encoding2', 'timestamps_Encoding2_end', 'timestamps_Encoding3', 'timestamps_Encoding3_end', 'timestamps_FixationCross', 'timestamps_Maintenance', 'timestamps_Probe', 'timestamps_Response', 'cat_Encoding1', 'cat_Encoding2', 'cat_Encoding3', 'cat_Probe', 'rt', 'correct', 'load']
loads: {3: 3074, 1: 3073} ; response_accuracy: {1: 5722, 0: 425} ; probe_in_out: {1: 3074, 0: 3073}
fixation_s picture1_s maintenance_s rt_s
count 6147.000 6147.000 6147.000 6147.000
mean 1.079 2.012 2.692 1.687
std 0.068 0.005 0.067 2.226
min 0.922 2.004 2.530 0.149
25% 1.033 2.007 2.635 0.932
50% 1.084 2.016 2.690 1.231
75% 1.124 2.016 2.745 1.736
max 1.431 2.036 3.203 74.214
load 3 blank between pictures (s): min 0.011, max 0.341 (Methods: 0.017 to 0.200)
maintenance onset equals the end of the last picture: True
category codes (first digit of the picture ID): {0: 6146, 1: 3822, 2: 3729, 3: 3688, 4: 3560, 5: 3643}
load 3 trials with three different categories: 3074 of 3074
probe_in_out equals 'probe is one of the encoded pictures': True
RT above 10 s: 69 trials; max 74.2 s
In [5]:
# The TTL events against the trials table: every trial timestamp should be an event with the
# code the file's description gives (11 fixation, 1-3 pictures, 6 maintenance, 7 probe, 8 response).
codes = {"timestamps_FixationCross": 11, "timestamps_Encoding1": 1, "timestamps_Encoding2": 2, "timestamps_Encoding3": 3,
"timestamps_Maintenance": 6, "timestamps_Probe": 7, "timestamps_Response": 8}
rows = []
for s in S:
ev = s["events"]
for col, code in codes.items():
t = s["trials"][col].to_numpy()
t = t[t != 0]
evt = np.sort(ev.t[ev.code == code].to_numpy())
i = np.clip(np.searchsorted(evt, t), 1, len(evt) - 1) if len(evt) > 1 else np.zeros(len(t), int)
dist = np.minimum(np.abs(evt[i] - t), np.abs(evt[i - 1] - t)) if len(evt) > 1 else np.full(len(t), np.inf)
rows.append({"path": s["path"], "column": col, "n": len(t), "max_offset_ms": 1000 * dist.max()})
evchk = pd.DataFrame(rows)
print(evchk.groupby("column").max_offset_ms.describe().round(3).to_string())
print("sessions with a trial timestamp more than 1 ms from its TTL event:",
sorted(set(evchk.path[evchk.max_offset_ms > 1])))
count mean std min 25% 50% 75% max column timestamps_Encoding1 44.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 timestamps_Encoding2 44.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 timestamps_Encoding3 44.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 timestamps_FixationCross 44.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 timestamps_Maintenance 44.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 timestamps_Probe 44.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 timestamps_Response 44.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 sessions with a trial timestamp more than 1 ms from its TTL event: []
In [6]:
# Spikes against the trial clock: spike times span the task, and spikes respond to picture onset.
rows = []
for s in S:
allsp = np.concatenate([sp for sp in s["spikes"] if len(sp)]) if len(s["spikes"]) else np.array([])
tr = s["trials"]
rows.append({"path": s["path"], "spike_min": allsp.min(), "spike_max": allsp.max(),
"task_start": tr.start_time.min(), "task_end": tr.stop_time.max(),
"units_with_unsorted_spike_times": len(s["unsorted_units"])})
spk = pd.DataFrame(rows)
print(spk.round(1).to_string(index=False))
print("units whose stored spike times are not sorted (sorted on loading), with the largest backward step (s) and the number of backward steps:",
{s["path"]: s["unsorted_units"] for s in S if s["unsorted_units"]})
rates_all = np.array([len(sp) / (s["trials"].stop_time.max() - s["trials"].start_time.min()) for s in S for sp in s["spikes"]])
print("mean rate over the task per unit (Hz): median %.2f, range %.3f to %.1f" % (np.median(rates_all), rates_all.min(), rates_all.max()))
path spike_min spike_max task_start task_end units_with_unsorted_spike_times
sub-1/sub-1_ses-1_ecephys+image.nwb 0.0 1375.7 1.0 1364.1 0
sub-1/sub-1_ses-2_ecephys+image.nwb 0.0 1303.1 0.2 1284.3 0
sub-10/sub-10_ses-1_ecephys+image.nwb 0.0 1423.4 0.0 1390.0 0
sub-11/sub-11_ses-1_ecephys+image.nwb 0.0 1436.9 0.0 1428.0 0
sub-12/sub-12_ses-1_ecephys+image.nwb 0.0 1391.8 0.0 1375.2 0
sub-12/sub-12_ses-2_ecephys+image.nwb 0.1 1363.0 0.0 1357.3 0
sub-13/sub-13_ses-1_ecephys+image.nwb 0.0 1340.3 0.1 1328.4 3
sub-14/sub-14_ses-1_ecephys+image.nwb 0.1 1447.6 0.2 1435.9 0
sub-15/sub-15_ses-1_ecephys+image.nwb 0.0 1557.7 0.0 1536.3 0
sub-16/sub-16_ses-1_ecephys+image.nwb 0.0 1428.3 0.0 1416.5 0
sub-17/sub-17_ses-1_ecephys+image.nwb 0.0 1328.4 0.2 1322.4 0
sub-18/sub-18_ses-1_ecephys+image.nwb 0.0 1452.2 0.0 1395.1 0
sub-19/sub-19_ses-1_ecephys+image.nwb 0.0 1447.3 0.0 1402.8 0
sub-19/sub-19_ses-2_ecephys+image.nwb 0.0 1392.2 0.0 1375.5 0
sub-2/sub-2_ses-1_ecephys+image.nwb 0.0 1545.1 0.0 1499.0 0
sub-20/sub-20_ses-1_ecephys+image.nwb 0.0 1236.9 0.0 1226.6 0
sub-21/sub-21_ses-1_ecephys+image.nwb 0.0 1425.4 0.7 1413.6 0
sub-22/sub-22_ses-1_ecephys+image.nwb 0.0 1220.8 0.6 1219.7 0
sub-22/sub-22_ses-2_ecephys+image.nwb 0.1 1452.6 0.7 1438.4 1
sub-22/sub-22_ses-3_ecephys+image.nwb 0.1 1332.8 0.7 1322.2 0
sub-23/sub-23_ses-1_ecephys+image.nwb 0.0 1414.7 0.1 1404.8 0
sub-24/sub-24_ses-1_ecephys+image.nwb 0.1 1738.6 0.6 1720.9 0
sub-25/sub-25_ses-1_ecephys+image.nwb 0.0 1526.0 0.0 1504.2 0
sub-26/sub-26_ses-1_ecephys+image.nwb 0.1 1543.5 0.0 1331.8 0
sub-26/sub-26_ses-2_ecephys+image.nwb 0.0 1326.0 0.1 1300.4 0
sub-27/sub-27_ses-1_ecephys+image.nwb 0.1 1482.8 0.0 1466.8 0
sub-28/sub-28_ses-1_ecephys+image.nwb 0.0 1690.9 0.7 1671.2 0
sub-29/sub-29_ses-1_ecephys+image.nwb 0.0 1516.2 0.0 1504.6 0
sub-3/sub-3_ses-1_ecephys+image.nwb 0.0 1422.7 0.0 1412.2 0
sub-30/sub-30_ses-1_ecephys+image.nwb 0.1 1805.0 0.1 1669.9 0
sub-31/sub-31_ses-1_ecephys+image.nwb 0.0 1324.2 0.1 1309.5 0
sub-31/sub-31_ses-2_ecephys+image.nwb 0.1 1304.2 0.2 1274.2 0
sub-32/sub-32_ses-1_ecephys+image.nwb 0.1 1358.8 0.0 1349.1 0
sub-33/sub-33_ses-1_ecephys+image.nwb 0.0 1336.2 0.0 1330.0 0
sub-34/sub-34_ses-1_ecephys+image.nwb 0.0 1315.7 0.0 1305.3 0
sub-35/sub-35_ses-1_ecephys+image.nwb 7.9 1517.7 0.2 1494.6 0
sub-36/sub-36_ses-1_ecephys+image.nwb 0.1 1433.6 0.0 1423.3 0
sub-4/sub-4_ses-1_ecephys+image.nwb 0.0 1369.2 0.0 1355.2 0
sub-5/sub-5_ses-1_ecephys+image.nwb 0.0 1334.5 0.0 1327.0 0
sub-5/sub-5_ses-2_ecephys+image.nwb 0.0 1470.8 0.0 1400.1 0
sub-6/sub-6_ses-1_ecephys+image.nwb 0.0 1573.3 0.0 1555.3 0
sub-7/sub-7_ses-1_ecephys+image.nwb 0.0 1899.2 0.0 1882.4 0
sub-8/sub-8_ses-1_ecephys+image.nwb 0.1 1342.9 0.0 1334.7 0
sub-9/sub-9_ses-1_ecephys+image.nwb 0.0 1288.3 0.0 1278.3 0
units whose stored spike times are not sorted (sorted on loading), with the largest backward step (s) and the number of backward steps: {'sub-13/sub-13_ses-1_ecephys+image.nwb': {6: (-1318.5681, 1), 14: (-1338.6463, 3), 15: (-1329.3442, 1)}, 'sub-22/sub-22_ses-2_ecephys+image.nwb': {10: (-0.0033, 24)}}
mean rate over the task per unit (Hz): median 1.57, range 0.001 to 22.9
In [7]:
# Clock check by contrast: across MTL units, the share with a category effect (one-way ANOVA on
# the 0.2-1.0 s window after picture onset, p < 0.05) at the true onsets and with onsets shifted.
shifts = [-3.0, -1.5, -0.6, 0.0, 0.6, 1.5, 3.0]
res = {sh: [] for sh in shifts}
for s in S:
ts, cats_ = sbcat.presentations(s)
for u in np.where(s["units"].area.isin(sbcat.MTL))[0]:
for sh in shifts:
r = sbcat.count_in(s["spikes"][u], ts + sh + 0.2, ts + sh + 1.0)
groups = [r[cats_ == c] for c in np.unique(cats_)]
res[sh].append(stats.f_oneway(*groups).pvalue if r.std() > 0 else 1.0)
print("share of MTL units with a category effect, by shift of the picture onsets (s):")
print({sh: round(float(np.mean(np.array(v) < 0.05)), 3) for sh, v in res.items()})
share of MTL units with a category effect, by shift of the picture onsets (s):
{-3.0: 0.049, -1.5: 0.053, -0.6: 0.079, 0.0: 0.331, 0.6: 0.186, 1.5: 0.091, 3.0: 0.039}
In [8]:
# Recordings deposited twice: compare spike trains and trial tables across files.
sig = []
for s in S:
first = tuple(np.round(np.concatenate([sp[:5] for sp in s["spikes"]]), 4)[:50])
sig.append((s["path"], tuple(np.round(s["trials"].rt.to_numpy()[:20], 4)), first))
dups = [(a[0], b[0]) for i, a in enumerate(sig) for b in sig[i + 1:] if a[1] == b[1] or a[2] == b[2]]
print("pairs of files with identical RT sequences or identical first spike times:", dups if dups else "none")
pairs of files with identical RT sequences or identical first spike times: none
All counts match the paper (checked under cohort-counts below). The trials table agrees with the
TTL events, the maintenance period starts at the end of the last picture and lasts 2.53 to 3.20 s
(the Methods say 2.55 to 2.85 s), load 3 trials hold three categories, and probe_in_out is
whether the probe was one of the encoded pictures. The share of MTL units with a category effect
peaks at the true picture onsets and falls with the shift, so the spike times share the trial
clock (consistent with the data; the shift contrast shows a response locked to the onsets).
cohort-counts¶
The study recorded 36 patients in 44 sessions, with 1,454 single neurons (360 hippocampus, 496 amygdala, 204 pre-SMA, 188 dACC, 206 vmPFC) and 1,922 microwire LFP channels (586, 421, 283, 325 and 307).
In [9]:
paper_units = {"hippocampus": 360, "amygdala": 496, "preSMA": 204, "dACC": 188, "vmPFC": 206}
paper_lfp = {"hippocampus": 586, "amygdala": 421, "preSMA": 283, "dACC": 325, "vmPFC": 307}
cc = pd.DataFrame({"units": units.area.value_counts(), "paper_units": pd.Series(paper_units),
"lfp": lfp.area.value_counts(), "paper_lfp": pd.Series(paper_lfp)}).loc[sbcat.AREAS]
print(cc.to_string())
print("totals: sessions %d (44), patients %d (36), units %d (1454), LFP channels %d (1922)" % (
len(info), info.subject.nunique(), len(units), len(lfp)))
print("units per patient: min %d, median %d, max %d" % tuple(units.groupby("subject").size().agg(["min", "median", "max"])))
units paper_units lfp paper_lfp hippocampus 360 360 586 586 amygdala 496 496 421 421 preSMA 204 204 283 283 dACC 188 188 325 325 vmPFC 206 206 307 307 totals: sessions 44 (44), patients 36 (36), units 1454 (1454), LFP channels 1922 (1922) units per patient: min 1, median 34, max 167
Result: reproduced. 44 sessions, 36 patients, 1,454 units and 1,922 LFP channels, and every per-area count equals the paper's.
fig1d-rt-load and fig1d-accuracy-load¶
Patients responded more slowly in load 3 than in load 1 trials (1.46 s against 1.33 s).
Patients were less accurate in load 3 than in load 1 trials (91.60% against 95.71%).
In [10]:
rng = np.random.default_rng(101)
beh = []
for s in S:
tr = s["trials"]
row = {"path": s["path"], "subject": s["subject"]}
for L in (1, 3):
t = tr[tr.load == L]
row[f"acc{L}"] = 100 * t.correct.mean()
row[f"rt_median{L}"] = t.rt[t.correct].median()
row[f"rt_mean{L}"] = t.rt[t.correct].mean()
beh.append(row)
beh = pd.DataFrame(beh)
for kind in ["rt_median", "rt_mean"]:
r = sbcat.perm_ttest_paired(beh[f"{kind}3"], beh[f"{kind}1"], rng)
print(f"{kind}: load 1 {beh[f'{kind}1'].mean():.3f} s, load 3 {beh[f'{kind}3'].mean():.3f} s; "
f"t{r['df']} = {r['t']:.2f}, permutation p = {r['p']:.4f}; sessions with load 3 slower: "
f"{int((beh[f'{kind}3'] > beh[f'{kind}1']).sum())} of {len(beh)}")
print("range of session median RT (s): %.2f to %.2f" % (beh[["rt_median1", "rt_median3"]].min().min(), beh[["rt_median1", "rt_median3"]].max().max()))
by_pat = beh.groupby("subject")[["rt_median1", "rt_median3"]].mean()
print("patients with load 3 slower (median RT, sessions averaged):", int((by_pat.rt_median3 > by_pat.rt_median1).sum()), "of", len(by_pat))
rt_median: load 1 1.233 s, load 3 1.360 s; t43 = 5.33, permutation p = 0.0001; sessions with load 3 slower: 36 of 44 rt_mean: load 1 1.633 s, load 3 1.596 s; t43 = -0.48, permutation p = 0.6628; sessions with load 3 slower: 28 of 44 range of session median RT (s): 0.55 to 3.11 patients with load 3 slower (median RT, sessions averaged): 29 of 36
In [11]:
rng = np.random.default_rng(102)
r = sbcat.perm_ttest_paired(beh.acc3, beh.acc1, rng)
print(f"accuracy: load 1 {beh.acc1.mean():.2f}%, load 3 {beh.acc3.mean():.2f}%; t{r['df']} = {r['t']:.2f}, permutation p = {r['p']:.4f}")
print("sessions with load 3 less accurate:", int((beh.acc3 < beh.acc1).sum()), "; equal:", int((beh.acc3 == beh.acc1).sum()), "of", len(beh))
by_pat = beh.groupby("subject")[["acc1", "acc3"]].mean()
print("patients with load 3 less accurate:", int((by_pat.acc3 < by_pat.acc1).sum()), "of", len(by_pat))
fig, ax = plt.subplots(1, 2, figsize=(7, 3))
for a, k, lab in [(ax[0], "rt_median", "median RT (s)"), (ax[1], "acc", "accuracy (%)")]:
a.plot([1, 3], beh[[f"{k}1", f"{k}3"]].T, color="0.7", lw=0.8)
a.plot([1, 3], beh[[f"{k}1", f"{k}3"]].mean(), "k-o")
a.set_xticks([1, 3]); a.set_xlabel("load"); a.set_ylabel(lab)
plt.tight_layout(); plt.show()
accuracy: load 1 95.10%, load 3 91.02%; t43 = -5.24, permutation p = 0.0001 sessions with load 3 less accurate: 32 ; equal: 5 of 44 patients with load 3 less accurate: 27 of 36
Result: reproduced (RT). The mean over sessions of the median RT of correct trials is 1.233 s in load 1 and 1.360 s in load 3, both within tolerance (1.33 and 1.46 s, 10%), and the paired permutation test gives t43 = 5.33, p = 0.0001 with load 3 slower. With the mean RT per session (variant) load 1 is 1.633 s and load 3 1.596 s (p = 0.66); RTs reach 74.2 s, with 69 trials above 10 s.
Result: reproduced (accuracy). 95.10% in load 1 and 91.02% in load 3, within tolerance of 95.71% and 91.60%, with t43 = -5.24, p = 0.0001 and load 3 less accurate.
behaviour-overall-accuracy¶
Mean accuracy was 93.66 +/- 7.04% (s.d. over sessions), and 78.34 +/- 20.09% of errors were false negatives.
In [12]:
acc = np.array([100 * s["trials"].correct.mean() for s in S])
fn = []
for s in S:
tr = s["trials"]
err = tr[~tr.correct]
if len(err):
fn.append(100 * (err.probe_in_out == 1).mean())
fn = np.array(fn)
print(f"accuracy over sessions: {acc.mean():.2f} +/- {acc.std(ddof=1):.2f}% (range {acc.min():.1f} to {acc.max():.1f})")
print(f"false negatives as a share of errors: {fn.mean():.2f} +/- {fn.std(ddof=1):.2f}% over {len(fn)} sessions with errors")
allerr = tr_all[~tr_all.correct]
print(f"pooled over all errors: {100 * (allerr.probe_in_out == 1).mean():.2f}% of {len(allerr)} errors")
accuracy over sessions: 93.05 +/- 7.22% (range 74.8 to 100.0) false negatives as a share of errors: 76.08 +/- 20.11% over 41 sessions with errors pooled over all errors: 75.29% of 425 errors
Result: reproduced. Mean accuracy is 93.05 +/- 7.22% (paper 93.66 +/- 7.04%), and false negatives are 76.08 +/- 20.11% of errors over the 41 sessions with errors (paper 78.34 +/- 20.09%); both are within tolerance.
category-cell-counts¶
Category neurons were found above chance in the hippocampus (89, 24.72%), amygdala (181, 36.49%) and vmPFC (37, 17.96%), but not in the pre-SMA and dACC.
In [13]:
t0 = time.time()
cat_table, cat_null = sbcat.category_table(S, seed=103)
print(f"selection and 500 label shuffles in {time.time() - t0:.0f} s")
expected = {"hippocampus": 89, "amygdala": 181, "preSMA": 18, "dACC": 14, "vmPFC": 37}
rows = []
for a in sbcat.AREAS:
m = (cat_table.area == a).to_numpy()
obs = int(cat_table.is_category[m].sum())
null = cat_null[m].sum(0)
rows.append({"area": a, "neurons": int(m.sum()), "category": obs, "percent": 100 * obs / m.sum(),
"paper": expected[a], "within_10pct": abs(obs - expected[a]) <= 0.1 * expected[a],
"null_mean": null.mean(), "null_p99": np.percentile(null, 99),
"p_shuffle": (np.sum(null >= obs) + 1) / (len(null) + 1), "above_p99": obs > np.percentile(null, 99)})
cat_counts = pd.DataFrame(rows)
print(cat_counts.round(3).to_string(index=False))
print("\ncategory neurons per patient (MTL):")
print(cat_table[cat_table.area.isin(sbcat.MTL)].groupby("subject").is_category.agg(["sum", "size"]).T.to_string())
print("\npreferred categories of MTL category neurons:",
cat_table[cat_table.is_category & cat_table.area.isin(sbcat.MTL)].preferred.value_counts().sort_index().to_dict())
# How the count depends on the post hoc threshold (one-sided 0.05 as two-sided 0.1, the code's)
for thr in [0.05, 0.1]:
c = cat_table[(cat_table.p_anova < 0.05) & (cat_table.p_posthoc < thr)].area.value_counts().reindex(sbcat.AREAS)
print(f"post hoc two-sided p < {thr}:", c.to_dict())
selection and 500 label shuffles in 117 s
area neurons category percent paper within_10pct null_mean null_p99 p_shuffle above_p99
hippocampus 360 94 26.111 89 True 16.794 26.01 0.002 True
amygdala 496 181 36.492 181 True 23.458 37.00 0.002 True
preSMA 204 18 8.824 18 True 9.450 17.00 0.010 True
dACC 188 13 6.915 14 True 8.300 16.00 0.088 False
vmPFC 206 40 19.417 37 True 9.522 18.00 0.002 True
category neurons per patient (MTL):
subject 1 10 11 12 13 14 15 16 17 18 19 2 20 21 22 23 24 25 26 27 28 29 3 30 31 32 33 34 35 36 4 5 6 7 8 9
sum 22 3 6 7 0 2 7 17 14 1 28 9 4 21 6 1 1 4 21 3 9 3 8 1 17 8 7 4 0 1 9 11 4 3 9 4
size 48 23 15 21 3 4 22 30 45 27 77 25 28 40 68 6 7 11 35 12 21 22 13 16 50 27 30 12 1 10 21 22 12 17 23 12
preferred categories of MTL category neurons: {1: 88, 2: 56, 3: 44, 4: 48, 5: 39}
post hoc two-sided p < 0.05: {'hippocampus': 87, 'amygdala': 174, 'preSMA': 18, 'dACC': 14, 'vmPFC': 41}
post hoc two-sided p < 0.1: {'hippocampus': 94, 'amygdala': 182, 'preSMA': 19, 'dACC': 15, 'vmPFC': 43}
Result: partially reproduced. Every count is within 10% of the paper's (hippocampus 94, amygdala 181, vmPFC 40, pre-SMA 18, dACC 13), and the hippocampus, amygdala and vmPFC exceed the 99th percentile of the shuffle null (p = 0.002 each). The dACC does not (p = 0.088), but the pre-SMA's 18 also exceeds its 99th percentile of 17 (p = 0.010), where the paper reports no excess. The pre-SMA sits at the boundary: 18 is one above the 99th percentile, and p = 0.010 equals the 0.01 threshold, so with 500 shuffles the two forms of the rule do not separate it from chance clearly.
Maintenance rates of category neurons¶
Per neuron, the mean rate during maintenance (0 to 2.5 s) over correct preferred trials (preferred category among the encoded pictures, G4) and correct non-preferred trials, and the baseline (-0.9 to -0.3 s before picture 1) over all correct trials. Reading (a) of O6 uses Hz; reading (b) the percentage change from each neuron's baseline (neurons with a baseline of 0 Hz are left out under (b)).
In [14]:
cat_units = cat_table[cat_table.is_category].reset_index(drop=True)
cat_rates = sbcat.trial_rates(S, cat_units)
def neuron_means(rates, cond_fn, reading):
"""Mean maintenance value per neuron over the trials cond_fn selects, in Hz or % of baseline."""
out = []
for r in rates:
c = cond_fn(r)
base = r["base"][r["correct"]].mean()
v = r["maint"][c].mean() if c.sum() else np.nan
if reading == "b":
v = 100 * (v - base) / base if base > 0 else np.nan
out.append(v)
return np.array(out)
def base_means(rates, reading):
if reading == "b":
return np.array([0.0 if r["base"][r["correct"]].mean() > 0 else np.nan for r in rates])
return np.array([r["base"][r["correct"]].mean() for r in rates])
def persistence_tests(sel, rng):
rows = []
for reading in ["a", "b"]:
rs = [cat_rates[i] for i in np.where(sel)[0]]
pref = neuron_means(rs, lambda r: r["correct"] & r["pref"], reading)
npref = neuron_means(rs, lambda r: r["correct"] & ~r["pref"], reading)
base = base_means(rs, reading)
for name, x, y in [("preferred vs baseline", pref, base), ("non-preferred vs baseline", npref, base),
("preferred vs non-preferred", pref, npref)]:
t = sbcat.perm_ttest_paired(x, y, rng)
rows.append({"reading": reading, "test": name, "n": t["n"], "mean_x": t["mean_x"], "mean_y": t["mean_y"],
"t": t["t"], "p": t["p"]})
return pd.DataFrame(rows)
fig3b-persistent-activity¶
MTL category neurons stayed above baseline during maintenance in both preferred and non-preferred trials, and fired more when their preferred category was held in memory.
In [15]:
rng = np.random.default_rng(104)
mtl = cat_units.area.isin(sbcat.MTL).to_numpy()
fig3b = persistence_tests(mtl, rng)
print(f"MTL category neurons: {mtl.sum()} (paper 270)")
print(fig3b.round(4).to_string(index=False))
print("neurons left out under (b), baseline 0 Hz:", int(np.isnan(base_means([cat_rates[i] for i in np.where(mtl)[0]], 'b')).sum()))
d_pp = neuron_means([cat_rates[i] for i in np.where(mtl)[0]], lambda r: r["correct"] & r["pref"], "a") - \
neuron_means([cat_rates[i] for i in np.where(mtl)[0]], lambda r: r["correct"] & ~r["pref"], "a")
per_pat = pd.Series(d_pp).groupby(cat_units.subject[mtl].to_numpy()).mean()
print("patients with preferred above non-preferred (mean over their neurons, Hz):", int((per_pat > 0).sum()), "of", len(per_pat))
print("range of the per-neuron preferred minus non-preferred difference (Hz): %.2f to %.2f" % (np.nanmin(d_pp), np.nanmax(d_pp)))
MTL category neurons: 275 (paper 270)
reading test n mean_x mean_y t p
a preferred vs baseline 275 2.8119 2.5599 5.7494 0.0001
a non-preferred vs baseline 275 2.7568 2.5599 3.7324 0.0001
a preferred vs non-preferred 275 2.8119 2.7568 1.7381 0.0724
b preferred vs baseline 275 32.0290 0.0000 3.9289 0.0001
b non-preferred vs baseline 275 20.4872 0.0000 5.1801 0.0001
b preferred vs non-preferred 275 32.0290 20.4872 2.2613 0.0013
neurons left out under (b), baseline 0 Hz: 0
patients with preferred above non-preferred (mean over their neurons, Hz): 24 of 34
range of the per-neuron preferred minus non-preferred difference (Hz): -6.31 to 1.33
Result: partially reproduced. Over 275 MTL category neurons (paper 270), maintenance rates exceed baseline in preferred and non-preferred trials under both readings (p = 0.0001 each). The preferred against non-preferred test gives p = 0.0013 in % of baseline (b) but p = 0.072 in Hz (a), so the criterion (all three tests under both readings) holds only in part. Under reading (b) alone, the reading adopted for the claims that leave the choice open, all three tests hold.
edf4b-persistent-activity-by-area¶
Persistent and selective maintenance activity of category neurons held in the hippocampus and amygdala separately, but vmPFC category neurons returned to baseline and did not differ between preferred and non-preferred trials.
In [16]:
rng = np.random.default_rng(105)
tabs = []
for a in ["hippocampus", "amygdala", "vmPFC"]:
sel = (cat_units.area == a).to_numpy()
t = persistence_tests(sel, rng)
t.insert(0, "area", a)
tabs.append(t)
edf4b = pd.concat(tabs, ignore_index=True)
print(edf4b.round(4).to_string(index=False))
area reading test n mean_x mean_y t p
hippocampus a preferred vs baseline 94 2.7287 2.4642 3.7924 0.0005
hippocampus a non-preferred vs baseline 94 2.6777 2.4642 3.0754 0.0010
hippocampus a preferred vs non-preferred 94 2.7287 2.6777 1.1859 0.2429
hippocampus b preferred vs baseline 94 55.7764 0.0000 2.5197 0.0001
hippocampus b non-preferred vs baseline 94 33.0136 0.0000 3.4020 0.0001
hippocampus b preferred vs non-preferred 94 55.7764 33.0136 1.6387 0.0188
amygdala a preferred vs baseline 181 2.8551 2.6096 4.3833 0.0001
amygdala a non-preferred vs baseline 181 2.7978 2.6096 2.6257 0.0024
amygdala a preferred vs non-preferred 181 2.8551 2.7978 1.3386 0.1834
amygdala b preferred vs baseline 181 19.6961 0.0000 4.4327 0.0001
amygdala b non-preferred vs baseline 181 13.9818 0.0000 4.3734 0.0001
amygdala b preferred vs non-preferred 181 19.6961 13.9818 2.0297 0.0385
vmPFC a preferred vs baseline 40 3.1048 2.9583 1.4370 0.1605
vmPFC a non-preferred vs baseline 40 3.0451 2.9583 0.8592 0.3934
vmPFC a preferred vs non-preferred 40 3.1048 3.0451 1.3684 0.1803
vmPFC b preferred vs baseline 40 10.7547 0.0000 1.8902 0.0644
vmPFC b non-preferred vs baseline 40 6.7503 0.0000 1.1899 0.2356
vmPFC b preferred vs non-preferred 40 10.7547 6.7503 1.5805 0.1232
Result: partially reproduced. In % of baseline (b) the hippocampus (94 neurons) and amygdala (181) pass all three tests (p <= 0.0385) and the vmPFC (40) passes none (p >= 0.064). In Hz (a) the vmPFC again passes none, but the preferred against non-preferred test fails in the hippocampus (p = 0.24) and amygdala (p = 0.18).
fig3c-category-load¶
In preferred trials, MTL category neurons fired more during maintenance in load 1 than in load 3, while no load effect appeared in non-preferred trials.
In [17]:
rng = np.random.default_rng(106)
rs = [cat_rates[i] for i in np.where(mtl)[0]]
rows = []
for reading in ["a", "b"]:
for name, pref in [("preferred", True), ("non-preferred", False)]:
x1 = neuron_means(rs, lambda r: r["correct"] & (r["pref"] == pref) & (r["load"] == 1), reading)
x3 = neuron_means(rs, lambda r: r["correct"] & (r["pref"] == pref) & (r["load"] == 3), reading)
t = sbcat.perm_ttest_paired(x1, x3, rng)
rows.append({"reading": reading, "trials": name, "n": t["n"], "load1": t["mean_x"], "load3": t["mean_y"], "t": t["t"], "p": t["p"]})
fig3c = pd.DataFrame(rows)
print(fig3c.round(4).to_string(index=False))
reading trials n load1 load3 t p
a preferred 275 2.8995 2.7804 2.2006 0.0223
a non-preferred 275 2.7759 2.7149 1.3315 0.1957
b preferred 275 35.2870 30.8506 0.4728 0.7731
b non-preferred 275 19.4940 22.5588 -1.1490 0.2583
Result: not reproduced, resting on an assumption. In % of baseline (b), the reading adopted (the legends show rates as % change from baseline and give the P values with those panels), preferred trials do not differ between loads (35.29 against 30.85%, t = 0.47, p = 0.77). In Hz (a) load 1 is higher in preferred trials (2.90 against 2.78 Hz, t = 2.20, p = 0.022) and non-preferred trials do not differ (p = 0.20), which would reproduce the claim.
fig3d-category-correct-incorrect¶
MTL category neurons fired more during maintenance in correct than in incorrect trials.
Readings: trials (a) all, (b) preferred only; minimum incorrect trials per neuron 1, 2 or 3 (O7); rates in Hz or % of baseline (O6).
In [18]:
rng = np.random.default_rng(107)
rows = []
for trials in ["all", "preferred"]:
for min_inc in [1, 2, 3]:
for reading in ["a", "b"]:
keep = []
for r in rs:
base = (r["pref"] if trials == "preferred" else np.ones(len(r["pref"]), bool))
keep.append(((~r["correct"]) & base).sum() >= min_inc)
keep = np.array(keep)
rk = [r for r, k in zip(rs, keep) if k]
tsel = (lambda r: r["pref"]) if trials == "preferred" else (lambda r: np.ones(len(r["pref"]), bool))
xc = neuron_means(rk, lambda r: r["correct"] & tsel(r), reading)
xi = neuron_means(rk, lambda r: ~r["correct"] & tsel(r), reading)
t = sbcat.perm_ttest_paired(xc, xi, rng)
rows.append({"trials": trials, "min_incorrect": min_inc, "reading": reading, "n": t["n"],
"correct": t["mean_x"], "incorrect": t["mean_y"], "t": t["t"], "p": t["p"]})
fig3d = pd.DataFrame(rows)
print(fig3d.round(4).to_string(index=False))
print("paper: n = 246, t245 = 2.43, p = 0.02")
trials min_incorrect reading n correct incorrect t p
all 1 a 251 2.8263 2.7646 0.8447 0.3985
all 1 b 251 26.5556 9.2744 2.5027 0.0035
all 2 a 229 2.8454 2.7420 1.5902 0.1176
all 2 b 229 26.5791 10.8279 2.1649 0.0100
all 3 a 209 2.9656 2.8599 1.5142 0.1315
all 3 b 209 24.1645 11.4863 1.6721 0.0545
preferred 1 a 217 3.0318 2.9993 0.3607 0.7319
preferred 1 b 217 30.7961 15.4536 1.4366 0.1042
preferred 2 a 168 3.0211 3.0265 -0.0607 0.9504
preferred 2 b 168 34.7653 16.4492 1.3662 0.1125
preferred 3 a 130 3.0210 3.0343 -0.1414 0.8944
preferred 3 b 130 27.3016 19.1330 1.4689 0.1486
paper: n = 246, t245 = 2.43, p = 0.02
Result: reproduced, resting on an assumption. In % of baseline (the adopted reading), over all trials, the 251 neurons with at least one incorrect trial (paper 246) fire more in correct than incorrect trials (26.56 against 9.27%, t = 2.50, p = 0.0035); with at least two incorrect trials p = 0.010, with three p = 0.055. In Hz the test fails (p >= 0.12), as it does on preferred trials only (p >= 0.10).
category-no-rt-effect¶
Category neuron maintenance rates did not differ between fast and slow RT trials in the hippocampus or amygdala.
In [19]:
rng = np.random.default_rng(108)
rows = []
for a in sbcat.MTL:
sel = np.where((cat_units.area == a).to_numpy())[0]
for reading in ["a", "b"]:
fast_v, slow_v = [], []
for i in sel:
r = cat_rates[i]
base = r["base"][r["correct"]].mean()
use = r["correct"] & r["pref"]
fast, slow = sbcat.rt_split(r["load"], r["rt"], use)
fv = np.nanmean([r["maint"][fast & (r["load"] == L)].mean() for L in (1, 3)])
sv = np.nanmean([r["maint"][slow & (r["load"] == L)].mean() for L in (1, 3)])
if reading == "b":
fv, sv = (100 * (fv - base) / base, 100 * (sv - base) / base) if base > 0 else (np.nan, np.nan)
fast_v.append(fv); slow_v.append(sv)
t = sbcat.perm_ttest_paired(fast_v, slow_v, rng)
rows.append({"area": a, "reading": reading, "n": t["n"], "fast": t["mean_x"], "slow": t["mean_y"], "t": t["t"], "p": t["p"]})
rt_tab = pd.DataFrame(rows)
print(rt_tab.round(4).to_string(index=False))
print("paper: hippocampus t88 = 0.73, p = 0.47; amygdala t180 = 1.38, p = 0.18")
area reading n fast slow t p hippocampus a 94 2.7559 2.7670 -0.1873 0.8529 hippocampus b 94 58.6398 49.5278 1.2597 0.2172 amygdala a 181 2.9248 2.8348 1.6628 0.0976 amygdala b 181 23.1773 20.9565 0.6869 0.5083 paper: hippocampus t88 = 0.73, p = 0.47; amygdala t180 = 1.38, p = 0.18
Result: reproduced. No area and reading gives p < 0.05 (hippocampus p = 0.85 and 0.22, amygdala p = 0.098 and 0.51).
edf8fgh-frontal-controls (panel h only)¶
Across all vmPFC neurons, neither within-region SFC nor SFC to hippocampal LFPs differed between loads, and maintenance firing rates of all neurons in the pre-SMA, dACC and vmPFC did not differ between correct and incorrect trials.
Panels f and g need spike-field coherence, which this attempt did not run, so the claim is blocked; panel h is shown here for the record.
In [20]:
rng = np.random.default_rng(109)
front = cat_table[cat_table.area.isin(["preSMA", "dACC", "vmPFC"])].reset_index(drop=True)
front_rates = sbcat.trial_rates(S, front)
rows = []
for a, n_paper in [("preSMA", 201), ("dACC", 180), ("vmPFC", 201)]:
sel = np.where((front.area == a).to_numpy())[0]
for min_inc in [1, 2, 3]:
for reading in ["a", "b"]:
rk = [front_rates[i] for i in sel if (~front_rates[i]["correct"]).sum() >= min_inc]
xc = neuron_means(rk, lambda r: r["correct"], reading)
xi = neuron_means(rk, lambda r: ~r["correct"], reading)
t = sbcat.perm_ttest_paired(xc, xi, rng)
rows.append({"area": a, "paper_n": n_paper, "min_incorrect": min_inc, "reading": reading, "n": t["n"],
"correct": t["mean_x"], "incorrect": t["mean_y"], "t": t["t"], "p": t["p"]})
print(pd.DataFrame(rows).round(4).to_string(index=False))
area paper_n min_incorrect reading n correct incorrect t p preSMA 201 1 a 204 2.9053 3.0681 -2.0700 0.0393 preSMA 201 1 b 201 28.4704 32.5767 -0.6259 0.5433 preSMA 201 2 a 197 2.9560 3.0857 -1.6883 0.0890 preSMA 201 2 b 194 30.1583 31.7510 -0.2458 0.8152 preSMA 201 3 a 188 2.9859 3.1153 -1.6178 0.1089 preSMA 201 3 b 185 31.9476 33.1020 -0.1705 0.8698 dACC 180 1 a 183 2.6865 2.7003 -0.2229 0.8293 dACC 180 1 b 180 13.8052 16.9608 -0.3807 0.7105 dACC 180 2 a 167 2.6841 2.6811 0.0603 0.9515 dACC 180 2 b 164 14.0666 19.3424 -0.5911 0.5692 dACC 180 3 a 167 2.6841 2.6811 0.0603 0.9550 dACC 180 3 b 164 14.0666 19.3424 -0.5911 0.5775 vmPFC 201 1 a 202 2.6177 2.5844 0.5273 0.5999 vmPFC 201 1 b 201 5.3232 5.0051 0.0609 0.9700 vmPFC 201 2 a 192 2.5091 2.4649 0.7756 0.4278 vmPFC 201 2 b 191 4.9283 5.5121 -0.1086 0.9498 vmPFC 201 3 a 192 2.5091 2.4649 0.7756 0.4316 vmPFC 201 3 b 191 4.9283 5.5121 -0.1086 0.9478
Result: blocked (panels f and g need spike-field coherence, not run). Panel h: in % of baseline (the adopted reading of O6) with at least one incorrect trial, the n are the paper's exactly (pre-SMA 201, dACC 180, vmPFC 201) and no area differs (p = 0.54, 0.71, 0.97); in Hz with at least one incorrect trial the pre-SMA differs (incorrect 3.07 against correct 2.91 Hz, p = 0.039), with two or three incorrect trials it does not (p = 0.089, 0.11). That the % of baseline reading gives the paper's n in all three areas (neurons with no baseline spike drop out) supports it as the reading of the paper's firing rate tests.
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome |
|---|---|
| cohort-counts | reproduced |
| fig1d-rt-load | reproduced |
| fig1d-accuracy-load | reproduced |
| behaviour-overall-accuracy | reproduced |
| category-cell-counts | partially reproduced (pre-SMA also above chance) |
| fig3b-persistent-activity | partially reproduced (selectivity only in % of baseline) |
| edf4b-persistent-activity-by-area | partially reproduced |
| fig3c-category-load | not reproduced, resting on an assumption (% of baseline) |
| fig3d-category-correct-incorrect | reproduced, resting on an assumption (% of baseline) |
| category-no-rt-effect | reproduced |
| edf8fgh-frontal-controls | blocked (SFC not run; panel h holds in % of baseline) |
The LFP claims are in reproduce_daume_2024_lfp.ipynb and reproduce_daume_2024_lfp_controls.ipynb, and the PAC neuron claims in
reproduce_daume_2024_pac_neurons.ipynb; the SFC, noise correlation and decoding claims are
blocked (attempt record).
How the analysis differs from the paper's: the category selection follows the paper's code, with the shuffled-label null drawn by shuffling the labels of a session's presentations together; the permutation tests flip the signs of paired differences (10,000 draws); whether rates are compared in Hz or in % of baseline is not stated, both are shown, and % of baseline is adopted where they disagree; units recorded in several sessions of one patient are counted once per session, as the paper's totals do.
What was found about the dataset: every count of the paper is in it; the trials table agrees with
the TTL events; spike times are unsorted in four units of two files (three in sub-13 jump back by
about 1,320 s, as if two spike trains were concatenated); the maintenance period lasts 2.53 to
3.20 s, against 2.55 to 2.85 s in the Methods; and RTs reach 74 s. The units table holds waveforms (waveforms, and waveform_mean of 256
samples in 43 of the 44 files).
In [21]:
print(f"notebook ran in {(time.time() - started) / 60:.1f} min")
notebook ran in 4.3 min