Reproducing Kyzar et al. (2024) from dandiset 000469¶
Paper: Kyzar, Kamiński, Brzezicka, Reed, Chung, Mamelak & Rutishauser (2024), "Dataset of human-single neuron activity during a Sternberg working memory task", Sci Data 11, 89, doi:10.1038/s41597-024-02943-8.
Dataset: DANDI:000469, version 0.240123.1806 (doi:10.48324/dandi.000469/0.240123.1806). The paper cites the dandiset without a version (ref. 46, "Kyzar M, 2023, DANDI Archive"), so the most recent published version is pinned. The four published versions (0.231213.2047, 0.240122.1800, 0.240122.1812, 0.240123.1806) and the draft list the same 41 assets with the same asset identifiers and sizes, so the choice does not change the data read; the cell below checks this against the archive.
Text the claims file was written from: the version of record (Europe PMC, PMC10796636), with
Methods, Technical Validation, legends and Tables 1 to 3; the paper has no supplement. The
claims file also read the paper's code, github.com/rutishauserlab/workingmem-release-NWB at
commit f8594e0. Where a method detail below is taken from that code, the helper module
sternberg.py cites the file.
| Claim | Tested here |
|---|---|
| subjects-sessions | tested |
| sternberg-unit-yield | tested |
| screening-unit-yield | tested |
| behavior-accuracy | tested |
| behavior-rt-load | tested |
| behavior-accuracy-load | tested |
| isi-violations | tested |
| firing-rate | tested |
| projection-distance | tested |
| waveform-snr | tested |
| cv2 | tested |
| isolation-distance | tested |
| screening-concept-mtl | tested |
| screening-concept-mfc | tested |
| fig4-screening-examples | tested |
| sternberg-concept-mtl | tested |
| persistent-activity | tested |
| fig5-sternberg-examples | tested |
| maintenance-cells-mfc | tested |
| probe-cells-mfc | tested; the readings the criterion names are all run (see its section) |
import time
import numpy as np
import pandas as pd
import requests
from scipy import stats
import matplotlib.pyplot as plt
import sternberg as sb
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
started = time.time()
# Versions of the dandiset: the assets of each published version, compared with the pinned one.
api = f"https://api.dandiarchive.org/api/dandisets/{sb.DANDISET_ID}/versions/"
versions = [v["version"] for v in requests.get(api, timeout=60).json()["results"]]
def assets(version):
r = requests.get(f"{api}{version}/assets/", params={"page_size": 200}, timeout=60).json()["results"]
return sorted((a["path"], a["asset_id"], a["size"]) for a in r)
pinned = assets(sb.DANDISET_VERSION)
for v in versions:
print(f"{v:>14}: {len(assets(v))} assets, identical to {sb.DANDISET_VERSION}: {assets(v) == pinned}")
draft: 41 assets, identical to 0.240123.1806: True
0.231213.2047: 41 assets, identical to 0.240123.1806: True
0.240122.1800: 41 assets, identical to 0.240123.1806: True
0.240122.1812: 41 assets, identical to 0.240123.1806: True
0.240123.1806: 41 assets, identical to 0.240123.1806: True
sessions = sb.load_all()
sc = [s for s in sessions if s["session"] == 1]
st = [s for s in sessions if s["session"] == 2]
units_sc, units_st = sb.units_table(sc), sb.units_table(st)
print(len(sessions), "files;", len(sc), "screening,", len(st), "Sternberg")
41 files; 20 screening, 21 Sternberg
Data check¶
What was loaded, the ranges and units of every variable used, and the checks of the description against the data.
# Files, subjects and sessions
files = pd.DataFrame([{"subject": s["subject"], "task": {1: "screening", 2: "Sternberg"}[s["session"]],
"identifier": s["identifier"], "age": s["age"], "sex": s["sex"],
"year": s["session_start_time"][:4], "units": len(s["units"]),
"trials": len(s["trials"]), "images": len(set(s["stim_ids"]))}
for s in sessions]).sort_values(["subject", "task"])
print(files.to_string(index=False))
native = files.assign(native=files.identifier.str.split("_", n=2).str[2])
print("\nsubject IDs:", files.subject.nunique(), "| native IDs:", native.native.nunique(),
"| native IDs without the '_2' suffix:", native.native.str.replace("_2", "").nunique())
ages = files.drop_duplicates("subject").age.str.strip("PY").astype(int)
print("ages in the files:", ages.min(), "to", ages.max(), "years")
# Units: ids are 0..n-1 in every file (the code's 1-based unit index is id + 1)
print("unit ids 0..n-1 in every file:",
all((s["units"].unit_id.to_numpy() == np.arange(len(s["units"]))).all() for s in sessions))
print("\nelectrode locations:", sorted(set(units_st.location) | set(units_sc.location)))
print("units by area and task:")
print(pd.DataFrame({"screening": units_sc.area.value_counts(), "Sternberg": units_st.area.value_counts()}))
# Spike times
for name, u in (("screening", units_sc), ("Sternberg", units_st)):
allspk = np.concatenate(u.spike_times.to_list())
unsorted = [(r.subject, r.unit_id, int(np.sum(np.diff(r.spike_times) < 0)),
float(np.diff(r.spike_times).min())) for r in u.itertuples() if np.any(np.diff(r.spike_times) < 0)]
print(f"{name}: {len(allspk)} spikes, times {allspk.min():.1f} to {allspk.max():.1f} s; "
f"spikes per unit {u.spike_times.map(len).min()} to {u.spike_times.map(len).max()}; "
f"units with unsorted spike times (subject, unit, decreasing steps, most negative step s): {unsorted}")
# Identical spike trains deposited twice (first 20 spike times of every unit, all files)
seen = {}
for s in sessions:
for r in s["units"].itertuples():
seen.setdefault(tuple(np.round(np.sort(r.spike_times)[:20], 6)), []).append((s["path"], r.unit_id))
print("spike trains found in more than one place:", [v for v in seen.values() if len(v) > 1])
# Trials: Sternberg
t_all = pd.concat([s["trials"].assign(subject=s["subject"]) for s in st], ignore_index=True)
print("\nSternberg trials:", len(t_all), "| per session:", sorted(t_all.groupby("subject").size().unique()))
print("loads:", t_all.loads.value_counts().sort_index().to_dict(),
"| Enc1 image IDs:", sorted(t_all.loadsEnc1_PicIDs.unique()),
"| Enc2:", sorted(t_all.loadsEnc2_PicIDs.unique()), "| Enc3:", sorted(t_all.loadsEnc3_PicIDs.unique()),
"| probe:", sorted(t_all.loadsProbe_PicIDs.unique()))
print("response_accuracy values:", sorted(t_all.response_accuracy.unique()),
"| probe_in_out values:", sorted(t_all.probe_in_out.unique()))
print("probe_in_out agrees with the probe being among the encoded images:",
((t_all.loadsProbe_PicIDs == t_all.loadsEnc1_PicIDs) | (t_all.loadsProbe_PicIDs == t_all.loadsEnc2_PicIDs)
| (t_all.loadsProbe_PicIDs == t_all.loadsEnc3_PicIDs)).astype(int).eq(t_all.probe_in_out).mean())
print("Enc2/Enc3 IDs are 0 exactly when the load leaves them out:",
((t_all.loadsEnc2_PicIDs == 0) == (t_all.loads < 2)).all(), ((t_all.loadsEnc3_PicIDs == 0) == (t_all.loads < 3)).all())
dur = pd.DataFrame({
"fixation to Enc1 (s)": t_all.timestamps_Encoding1 - t_all.timestamps_FixationCross,
"Enc1 shown (s)": t_all.timestamps_Encoding1_end - t_all.timestamps_Encoding1,
"maintenance to probe (s)": t_all.timestamps_Probe - t_all.timestamps_Maintenance,
"reaction time (s)": t_all.timestamps_Response - t_all.timestamps_Probe})
print(dur.describe(percentiles=[0.01, 0.5, 0.99]).round(3).to_string())
print("start_time equals timestamps_FixationCross:", (t_all.start_time == t_all.timestamps_FixationCross).all(),
"| stop_time equals timestamps_Response:", (t_all.stop_time == t_all.timestamps_Response).all())
for k in (2, 3):
on = t_all[f"timestamps_Encoding{k}"]
print(f"Enc{k} timestamp is 0 or NaN exactly when the load leaves it out:",
((on.isna() | (on == 0)) == (t_all.loads < k)).all(), "| values when left out:",
sorted(on[t_all.loads < k].fillna(-1).unique()), "| smallest when shown: %.1f s" % on[t_all.loads >= k].min())
odd = t_all[(dur["maintenance to probe (s)"] < 2) | (dur["fixation to Enc1 (s)"] > 2) | (dur["fixation to Enc1 (s)"] < 0.5)]
print("trials with maintenance under 2 s, or fixation under 0.5 s or over 2 s:")
print(odd[["subject", "loads", "timestamps_FixationCross", "timestamps_Encoding1", "timestamps_Encoding1_end",
"timestamps_Maintenance", "timestamps_Probe"]].to_string())
# Trials: screening
for s in sc:
assert np.allclose(s["trials"].start_time.to_numpy(), s["stim_times"])
counts = {s["subject"]: sorted(set(np.bincount(s["stim_ids"])[np.unique(s["stim_ids"])])) for s in sc}
print("\nscreening: start_time equals the StimulusPresentation timestamps in every file; presentations of each image:",
sorted({c for v in counts.values() for c in v}), "| images per session:",
sorted({len(set(s["stim_ids"])) for s in sc}))
dur_sc = np.concatenate([np.diff(s["trials"].start_time) for s in sc])
print("screening onset to next onset (s): median %.3f, 1st to 99th percentile %.3f to %.3f" %
(np.median(dur_sc), *np.percentile(dur_sc, [1, 99])))
# Quality metrics stored per unit
print("\nstored quality metrics, Sternberg units:")
print(units_st[["isolation_distance", "proj_dist", "mean_snr", "peak_snr"]].describe().round(3).to_string())
print("NaN isolation distance:", units_st.isolation_distance.isna().sum(), "| zero projection distance:",
(units_st.proj_dist == 0).sum())
subject task identifier age sex year units trials images
1 Sternberg SBID_1_P42CS P25Y F 2016 40 135 6
1 screening SCID_1_P42CS P25Y F 2016 40 378 63
2 Sternberg SBID_2_P43CS P42Y F 2016 32 135 6
2 screening SCID_2_P43CS P42Y F 2016 22 378 63
3 Sternberg SBID_3_P44CS P53Y F 2016 28 135 6
3 screening SCID_3_P44CS P53Y F 2016 40 378 63
4 Sternberg SBID_4_P47CS P32Y M 2017 76 135 6
4 screening SCID_4_P47CS P32Y M 2017 63 378 63
5 Sternberg SBID_5_P48CS P32Y F 2017 50 135 6
5 screening SCID_5_P48CS P32Y F 2017 36 378 63
6 Sternberg SBID_6_P31CS P32Y M 2014 35 108 6
6 screening SCID_6_P31CS P32Y M 2014 28 324 54
7 Sternberg SBID_7_P32CS P19Y M 2014 22 108 6
7 screening SCID_7_P32CS P19Y M 2014 34 324 54
8 Sternberg SBID_8_P33CS P44Y F 2014 74 108 6
8 screening SCID_8_P33CS P44Y F 2014 77 324 54
9 Sternberg SBID_9_P34CS P70Y M 2014 43 108 6
9 screening SCID_9_P34CS P70Y M 2014 53 324 54
10 Sternberg SBID_10_P35CS P63Y M 2014 28 135 6
10 screening SCID_10_P35CS P63Y M 2014 74 378 63
11 Sternberg SBID_11_P36CS P45Y M 2015 74 135 6
11 screening SCID_11_P36CS P45Y M 2015 94 378 63
12 Sternberg SBID_12_P37CS P33Y F 2015 22 135 6
12 screening SCID_12_P37CS P33Y F 2015 18 378 63
13 Sternberg SBID_13_P39CS P26Y M 2015 58 135 6
13 screening SCID_13_P39CS P26Y M 2015 103 378 63
14 Sternberg SBID_14_P39CS_2 P26Y M 2015 65 135 6
14 screening SCID_14_P39CS_2 P26Y M 2015 73 378 63
15 Sternberg SBID_15_P40CS P25Y M 2016 55 135 6
15 screening SCID_15_P40CS P25Y M 2016 32 378 63
16 Sternberg SBID_16_P48HMH P54Y M 2014 33 108 6
16 screening SCID_16_P48HMH P54Y M 2014 34 324 54
17 Sternberg SBID_17_P49HMH P54Y F 2015 13 135 6
17 screening SCID_17_P49HMH P54Y F 2015 8 378 63
18 Sternberg SBID_18_P51HMH P24Y M 2016 118 135 6
18 screening SCID_18_P51HMH P24Y M 2016 53 378 63
19 Sternberg SBID_19_P47HMH P20Y M 2014 11 108 6
20 Sternberg SBID_20_P49CS P24Y F 2017 5 135 6
20 screening SCID_20_P49CS P24Y F 2017 7 378 63
21 Sternberg SBID_21_P51CS P17Y M 2017 20 135 6
21 screening SCID_21_P51CS P17Y M 2017 18 378 63
subject IDs: 21 | native IDs: 21 | native IDs without the '_2' suffix: 20
ages in the files: 17 to 70 years
unit ids 0..n-1 in every file: True
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 task:
screening Sternberg
area
Amg 262 259
Hipp 159 190
dACC 221 171
pre-SMA 236 250
vmPFC 29 32
screening: 1833840 spikes, times 0.0 to 871.3 s; spikes per unit 55 to 19686; units with unsorted spike times (subject, unit, decreasing steps, most negative step s): [(11, 26, 1, -0.00975000000000037), (11, 32, 1, -0.012375249999999838), (5, 18, 1, -0.0009064999999850443)]
Sternberg: 2682324 spikes, times 0.0 to 1593.2 s; spikes per unit 88 to 43427; units with unsorted spike times (subject, unit, decreasing steps, most negative step s): [(2, 22, 1, -0.0003440000000409782)]
spike trains found in more than one place: []
Sternberg trials: 2673 | per session: [np.int64(108), np.int64(135)]
loads: {1: 891, 2: 891, 3: 891} | Enc1 image IDs: [np.uint8(1), np.uint8(2), np.uint8(3), np.uint8(4), np.uint8(5)] | Enc2: [np.uint8(0), np.uint8(1), np.uint8(2), np.uint8(3), np.uint8(4), np.uint8(5)] | Enc3: [np.uint8(0), np.uint8(1), np.uint8(2), np.uint8(3), np.uint8(4), np.uint8(5)] | probe: [np.uint8(1), np.uint8(2), np.uint8(3), np.uint8(4), np.uint8(5)]
response_accuracy values: [np.uint8(0), np.uint8(1)] | probe_in_out values: [np.uint8(0), np.uint8(1)]
probe_in_out agrees with the probe being among the encoded images: 1.0
Enc2/Enc3 IDs are 0 exactly when the load leaves them out: True True
fixation to Enc1 (s) Enc1 shown (s) maintenance to probe (s) reaction time (s)
count 2673.000 2673.000 2673.000 2673.000
mean 1.059 1.236 2.679 1.416
std 0.261 0.415 0.081 1.291
min 0.000 1.002 0.000 0.337
1% 0.928 1.002 2.556 0.490
50% 1.054 1.017 2.684 1.124
99% 1.200 2.507 2.813 5.850
max 14.070 2.509 2.836 30.829
start_time equals timestamps_FixationCross: True | stop_time equals timestamps_Response: True
Enc2 timestamp is 0 or NaN exactly when the load leaves it out: True | values when left out: [np.float64(0.0)] | smallest when shown: 2.0 s
Enc3 timestamp is 0 or NaN exactly when the load leaves it out: True | values when left out: [np.float64(0.0)] | smallest when shown: 12.4 s
trials with maintenance under 2 s, or fixation under 0.5 s or over 2 s:
subject loads timestamps_FixationCross timestamps_Encoding1 timestamps_Encoding1_end timestamps_Maintenance timestamps_Probe
1769 21 2 484.259399 498.329706 499.342518 500.509236 503.162704
2593 9 1 525.715789 526.783195 529.288256 529.288256 529.288256
2653 9 3 1189.469398 1189.469429 1191.975647 1197.326207 1199.983518
screening: start_time equals the StimulusPresentation timestamps in every file; presentations of each image: [np.int64(6)] | images per session: [54, 63]
screening onset to next onset (s): median 1.500, 1st to 99th percentile 1.134 to 5.571
stored quality metrics, Sternberg units:
isolation_distance proj_dist mean_snr peak_snr
count 676.000 902.000 902.000 902.000
mean 211.683 11.209 1.647 6.625
std 623.354 8.511 0.954 3.632
min 4.151 0.000 0.417 2.650
25% 15.839 5.593 1.017 4.477
50% 31.138 10.428 1.338 5.301
75% 99.408 16.435 2.034 7.546
max 7928.563 54.125 7.895 39.319
NaN isolation distance: 226 | zero projection distance: 168
The trial table timestamps are checked against the TTL stream in acquisition/events (codes from
Tables 2 and 3 of the paper). The spikes are checked against the trial clock with a contrast
between conditions: the fraction of units whose response 200 to 1000 ms after image onset differs
between images, with the onsets shifted by -2 to +2 s. On the right clock it peaks at 0 s; at
±1 s the window falls on the next or previous image (screening, and Encoding 2 in Sternberg).
ttl = pd.DataFrame([dict(file=s["path"].split("/")[1], **sb.ttl_mismatch(s)) for s in sessions]).set_index("file")
print("trials whose timestamp has no TTL of the expected code within 1 ms, summed over files:")
print(ttl.sum().to_string())
shifts = (-2.0, -1.0, 0.0, 1.0, 2.0)
sel = pd.DataFrame([{"file": s["path"].split("/")[1], "task": s["session"], "units": len(s["units"]),
**sb.selectivity_by_shift(s, shifts)} for s in sessions]).set_index("file")
for task in (1, 2):
sub = sel[sel.task == task]
w = sub[list(shifts)].mul(sub.units, axis=0).sum() / sub.units.sum()
print(f"\n{'screening' if task == 1 else 'Sternberg'}: fraction of selective units by shift (s), pooled:",
w.round(3).to_dict())
not_peak = sel[sel[list(shifts)].idxmax(axis=1) != 0.0]
print("\nfiles where shift 0 is not the highest (clock not established by this check):")
print(not_peak.round(3).to_string())
trials whose timestamp has no TTL of the expected code within 1 ms, summed over files: start_time 0.0 timestamps_FixationCross 0.0 timestamps_Encoding1 0.0 timestamps_Maintenance 0.0 timestamps_Probe 0.0 timestamps_Response 0.0
screening: fraction of selective units by shift (s), pooled: {-2.0: 0.034, -1.0: 0.055, 0.0: 0.181, 1.0: 0.058, 2.0: 0.037}
Sternberg: fraction of selective units by shift (s), pooled: {-2.0: 0.053, -1.0: 0.039, 0.0: 0.142, 1.0: 0.073, 2.0: 0.049}
files where shift 0 is not the highest (clock not established by this check):
task units -2.0 -1.0 0.0 1.0 2.0
file
sub-1_ses-1_ecephys+image.nwb 1 40 0.025 0.075 0.075 0.050 0.025
sub-1_ses-2_ecephys+image.nwb 2 40 0.075 0.075 0.050 0.000 0.025
sub-10_ses-1_ecephys+image.nwb 1 74 0.054 0.068 0.068 0.014 0.027
sub-10_ses-2_ecephys+image.nwb 2 28 0.071 0.107 0.036 0.036 0.000
sub-15_ses-1_ecephys+image.nwb 1 32 0.062 0.000 0.031 0.062 0.094
sub-15_ses-2_ecephys+image.nwb 2 55 0.018 0.036 0.036 0.036 0.091
sub-17_ses-1_ecephys+image.nwb 1 8 0.000 0.000 0.000 0.000 0.000
sub-17_ses-2_ecephys+image.nwb 2 13 0.154 0.000 0.077 0.000 0.000
sub-19_ses-2_ecephys+image.nwb 2 11 0.000 0.000 0.091 0.182 0.000
sub-2_ses-1_ecephys+image.nwb 1 22 0.045 0.091 0.045 0.000 0.091
sub-20_ses-1_ecephys+image.nwb 1 7 0.000 0.143 0.000 0.000 0.000
sub-20_ses-2_ecephys+image.nwb 2 5 0.000 0.000 0.000 0.000 0.000
sub-3_ses-1_ecephys+image.nwb 1 40 0.025 0.025 0.050 0.075 0.025
sub-5_ses-2_ecephys+image.nwb 2 50 0.100 0.060 0.080 0.080 0.020
sub-6_ses-1_ecephys+image.nwb 1 28 0.071 0.000 0.000 0.036 0.000
sub-6_ses-2_ecephys+image.nwb 2 35 0.000 0.057 0.000 0.114 0.057
sub-8_ses-2_ecephys+image.nwb 2 74 0.027 0.041 0.068 0.108 0.041
sub-9_ses-1_ecephys+image.nwb 1 53 0.019 0.075 0.075 0.038 0.094
# Column descriptions stored in the files, for the columns used here and the electrodes table
d = pd.DataFrame([{"file": s["path"].split("/")[1], **s["descriptions"]} for s in sessions]).set_index("file")
for col in ("units/waveforms_isolation_distance", "units/waveforms_mean_proj_dist", "trials/loadsEnc1_PicIDs",
"trials/loadsEnc2_PicIDs", "trials/loadsEnc3_PicIDs"):
print(f"{col}:", d[col].dropna().value_counts().to_dict())
print("electrodes table, descriptions per column over all files:")
for col in [c for c in d.columns if c.startswith("electrodes/")]:
print(f" {col}:", d[col].dropna().value_counts().to_dict())
units/waveforms_isolation_distance: {'Cluster Isolation distance, computed using all waveforms in the cluster.': 41}
units/waveforms_mean_proj_dist: {'Cluster Isolation distance, computed using all waveforms in the cluster.': 41}
trials/loadsEnc1_PicIDs: {'Picture ID for Enc1 loads.': 21}
trials/loadsEnc2_PicIDs: {'Picture ID for Enc2 loads.': 21}
trials/loadsEnc3_PicIDs: {'Picture ID for Enc1 loads.': 21}
electrodes table, descriptions per column over all files:
electrodes/filtering: {'my description': 41}
electrodes/group: {'my description': 41}
electrodes/group_name: {'my description': 41}
electrodes/id: {'': 41}
electrodes/location: {'my description': 41}
electrodes/origChannel: {'my description': 41}
electrodes/x: {'my description': 41}
electrodes/y: {'my description': 41}
electrodes/z: {'my description': 41}
Unit and session counts¶
subjects-sessions¶
The dataset holds 41 recording sessions from 21 patients, every patient with a Sternberg session and all but one (P47HMH) with a screening session.
by_subject = files.pivot_table(index="subject", columns="task", values="units", aggfunc="count").fillna(0).astype(int)
print("subjects:", files.subject.nunique(), "| sessions:", len(files), "| Sternberg sessions:",
(files.task == "Sternberg").sum(), "| screening sessions:", (files.task == "screening").sum())
missing = by_subject[by_subject.screening == 0].index.tolist()
print("subjects without a screening session:", missing,
files[files.subject.isin(missing)].identifier.tolist())
print("subjects 13 and 14:", files[files.subject.isin([13, 14])][["subject", "identifier", "age", "sex", "year"]]
.drop_duplicates("subject").to_string(index=False))
for s_ in sessions:
if s_["subject"] in (13, 14):
e_ = s_["electrodes"]
print(f"subject {s_['subject']}, {s_['path'].split('/')[1].split('_')[1]}: {len(s_['units'])} units, {len(e_)} electrodes,",
e_.location.value_counts().sort_index().to_dict())
print("ages:", ages.min(), "to", ages.max(), "(text 17 to 71; Table 1 17 to 70)")
subjects: 21 | sessions: 41 | Sternberg sessions: 21 | screening sessions: 20
subjects without a screening session: [19] ['SBID_19_P47HMH']
subjects 13 and 14: subject identifier age sex year
13 SBID_13_P39CS P26Y M 2015
14 SBID_14_P39CS_2 P26Y M 2015
subject 13, ses-1: 103 units, 48 electrodes, {'amygdala_left': 8, 'amygdala_right': 8, 'dorsal_anterior_cingulate_cortex_left': 2, 'dorsal_anterior_cingulate_cortex_right': 4, 'hippocampus_left': 7, 'hippocampus_right': 4, 'pre_supplementary_motor_area_left': 8, 'pre_supplementary_motor_area_right': 7}
subject 13, ses-2: 58 units, 33 electrodes, {'amygdala_left': 2, 'amygdala_right': 8, 'dorsal_anterior_cingulate_cortex_left': 1, 'dorsal_anterior_cingulate_cortex_right': 2, 'hippocampus_left': 7, 'hippocampus_right': 1, 'pre_supplementary_motor_area_left': 7, 'pre_supplementary_motor_area_right': 5}
subject 14, ses-1: 73 units, 34 electrodes, {'amygdala_left': 1, 'amygdala_right': 8, 'dorsal_anterior_cingulate_cortex_right': 5, 'hippocampus_left': 6, 'hippocampus_right': 2, 'pre_supplementary_motor_area_left': 4, 'pre_supplementary_motor_area_right': 8}
subject 14, ses-2: 65 units, 34 electrodes, {'amygdala_right': 7, 'dorsal_anterior_cingulate_cortex_left': 5, 'dorsal_anterior_cingulate_cortex_right': 3, 'hippocampus_left': 6, 'hippocampus_right': 2, 'pre_supplementary_motor_area_left': 5, 'pre_supplementary_motor_area_right': 6}
ages: 17 to 70 (text 17 to 71; Table 1 17 to 70)
Result: reproduced. The files hold 21 subject IDs, 41 sessions and 21 Sternberg sessions; subject 19 (P47HMH) is the only one without a screening session. Subjects 13 and 14 carry the native IDs P39CS and P39CS_2, with the same age, sex and recording year, so the 21 subjects may be 20 people; their files hold different units and electrodes (subject 13: 48 and 33 electrodes, subject 14: 34 and 34). Ages in the files run from 17 to 70, as in Table 1.
sternberg-unit-yield¶
In the Sternberg task 902 units were recorded, 190 in the hippocampus, 259 in the amygdala, 171 in the dACC, 250 in the pre-SMA and 32 in the vmPFC.
areas = ["Hipp", "Amg", "dACC", "pre-SMA", "vmPFC"]
print("Sternberg units:", len(units_st), units_st.area.value_counts().reindex(areas).to_dict())
print("MTL:", units_st.area.isin(sb.MTL).sum(), "| MFC:", units_st.area.isin(sb.MFC).sum())
print(pd.crosstab(units_st.subject, units_st.area).reindex(columns=areas, fill_value=0).T.to_string())
Sternberg units: 902 {'Hipp': 190, 'Amg': 259, 'dACC': 171, 'pre-SMA': 250, 'vmPFC': 32}
MTL: 449 | MFC: 453
subject 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
area
Hipp 2 15 4 0 13 1 0 25 6 8 20 0 11 17 15 6 0 33 6 0 8
Amg 6 3 1 42 6 0 16 2 34 4 6 4 25 19 20 27 6 31 5 1 1
dACC 11 0 9 1 20 10 0 33 0 0 20 11 6 10 8 0 3 19 0 1 9
pre-SMA 21 0 14 18 8 24 6 14 3 16 28 7 16 19 12 0 4 35 0 3 2
vmPFC 0 14 0 15 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Result: reproduced. 902 Sternberg units: 190 hippocampus, 259 amygdala, 171 dACC, 250 pre-SMA and 32 vmPFC, every count equal to the paper's.
screening-unit-yield¶
In the screening task 907 units were recorded, 449 of them in the MTL.
n_mtl_sc = units_sc.area.isin(sb.MTL).sum()
print("screening units:", len(units_sc), units_sc.area.value_counts().reindex(areas).to_dict())
print("MTL:", n_mtl_sc, "| MFC:", units_sc.area.isin(sb.MFC).sum())
for reading in (449, 421):
print(f"MTL {n_mtl_sc} within 10% of {reading}: {abs(n_mtl_sc - reading) <= 0.1 * reading}")
print("Sternberg MTL for comparison:", units_st.area.isin(sb.MTL).sum())
screening units: 907 {'Hipp': 159, 'Amg': 262, 'dACC': 221, 'pre-SMA': 236, 'vmPFC': 29}
MTL: 421 | MFC: 486
MTL 421 within 10% of 449: True
MTL 421 within 10% of 421: True
Sternberg MTL for comparison: 449
Result: reproduced. 907 screening units, as stated. The MTL holds 421 of them (159 hippocampus, 262 amygdala), which is within 10% of 449 and equal to the 421 implied by the screening percentages. 449 is the Sternberg MTL count, so the text appears to give the Sternberg number for the screening task.
Behaviour¶
behavior-accuracy¶
Subjects answered on average 88.73% of Sternberg trials correctly.
acc = pd.Series({s["subject"]: 100 * s["trials"].response_accuracy.mean() for s in st}).sort_index()
print("percent correct per session:", acc.round(2).to_dict())
print("mean over sessions: %.2f%%, s.e.m. %.2f, range %.2f to %.2f" % (acc.mean(), acc.std() / np.sqrt(len(acc)), acc.min(), acc.max()))
percent correct per session: {1: 76.3, 2: 91.11, 3: 91.85, 4: 94.07, 5: 85.93, 6: 75.0, 7: 95.37, 8: 88.89, 9: 87.04, 10: 88.89, 11: 90.37, 12: 91.11, 13: 84.44, 14: 91.11, 15: 83.7, 16: 94.44, 17: 95.56, 18: 93.33, 19: 94.44, 20: 85.19, 21: 85.19}
mean over sessions: 88.73%, s.e.m. 1.25, range 75.00 to 95.56
Result: reproduced. Mean accuracy over the 21 sessions is 88.73%, the paper's value.
behavior-rt-load¶
Reaction times grew with memory load; a repeated-measures ANOVA of per-session median reaction times found a load effect.
for label, correct in (("correct trials (main)", True), ("all trials (variant)", False)):
m = sb.rt_medians(st, correct_only=correct)
f, d1, d2, p = sb.rm_anova_load(m)
print(f"{label}: mean of session medians by load (s) {m.mean(0).round(3).tolist()}; "
f"F({d1},{d2}) = {f:.2f}, p = {p:.4f}; sessions with load 3 slower than load 1: {(m[:, 2] > m[:, 0]).sum()}/{len(m)}")
m = sb.rt_medians(st)
print("median RT range over sessions and loads: %.3f to %.3f s" % (m.min(), m.max()))
fig, ax = plt.subplots(figsize=(3.5, 3))
ax.plot([1, 2, 3], m.T, color="0.7", lw=0.8)
ax.plot([1, 2, 3], m.mean(0), color="C0", lw=2)
ax.set(xlabel="load", ylabel="median RT, correct trials (s)", xticks=[1, 2, 3])
plt.show()
correct trials (main): mean of session medians by load (s) [1.055, 1.13, 1.198]; F(2,40) = 7.92, p = 0.0013; sessions with load 3 slower than load 1: 17/21 all trials (variant): mean of session medians by load (s) [1.098, 1.154, 1.224]; F(2,40) = 4.08, p = 0.0244; sessions with load 3 slower than load 1: 17/21 median RT range over sessions and loads: 0.677 to 1.841 s
Result: reproduced. On correct trials the mean of the session medians rises from 1.055 s (load 1) to 1.130 and 1.198 s, and the load effect is F(2,40) = 7.92, p = 0.0013, the paper's values. On all trials the effect is weaker (F(2,40) = 4.08, p = 0.024).
behavior-accuracy-load¶
Accuracy was lower at higher memory loads.
acc_load = np.array([[100 * s["trials"][s["trials"].loads == L].response_accuracy.mean() for L in (1, 2, 3)] for s in st])
print("mean percent correct by load:", acc_load.mean(0).round(2).tolist())
print("per subject:")
print(pd.DataFrame(acc_load, index=[s["subject"] for s in st], columns=["load 1", "load 2", "load 3"]).round(1).sort_index().T.to_string())
mean percent correct by load: [91.22, 88.62, 86.35]
per subject:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
load 1 68.9 91.1 93.3 100.0 95.6 83.3 97.2 94.4 88.9 97.8 93.3 95.6 88.9 95.6 82.2 97.2 100.0 95.6 94.4 77.8 84.4
load 2 77.8 91.1 93.3 88.9 77.8 69.4 94.4 91.7 83.3 84.4 86.7 95.6 86.7 91.1 88.9 91.7 91.1 95.6 97.2 95.6 88.9
load 3 82.2 91.1 88.9 93.3 84.4 72.2 94.4 80.6 88.9 84.4 91.1 82.2 77.8 86.7 80.0 94.4 95.6 88.9 91.7 82.2 82.2
Result: reproduced. Mean accuracy falls from 91.22% at load 1 to 88.62% at load 2 and 86.35% at load 3.
Spike sorting quality (Sternberg units)¶
All six measures follow NWB_QA_graphs.m; spike times are used as stored, as the code does.
isi = np.array([sb.isi_violation_pct(t) for t in units_st.spike_times])
rate = np.array([sb.mean_rate(t) for t in units_st.spike_times])
cv2s = np.array([sb.cv2(t) for t in units_st.spike_times])
print("units:", len(units_st))
units: 902
isi-violations¶
Across Sternberg units, refractory period violations (ISIs shorter than 3 ms) were 0.41% of ISIs on average, with a standard deviation of 0.50%.
print("all units (main reading): mean %.3f%%, s.d. %.3f%%, n = %d" % (isi.mean(), isi.std(ddof=1), len(isi)))
keep = isi < 4
print("units below 4%% (code's thresh_isi): mean %.3f%%, s.d. %.3f%%, n = %d; left out: %s" %
(isi[keep].mean(), isi[keep].std(ddof=1), keep.sum(), np.round(isi[~keep], 2).tolist()))
print("range: %.3f to %.3f%%" % (isi.min(), isi.max()))
all units (main reading): mean 0.428%, s.d. 0.625%, n = 902 units below 4% (code's thresh_isi): mean 0.416%, s.d. 0.498%, n = 901; left out: [11.74] range: 0.000 to 11.739%
Result: reproduced. Over all 902 units the mean is 0.428% (within tolerance) and the s.d. 0.625% (outside 0.50 ± 0.05). With the code's filter, which leaves out one unit, the values are 0.416% and 0.498%, both within tolerance; the claim holds on that reading. The one unit left out has 11.74% of its ISIs below 3 ms.
firing-rate¶
Sternberg units fired at 2.55 Hz on average, ranging from 0.11 to 33.82 Hz.
print("mean %.3f Hz, min %.3f Hz, max %.3f Hz, median %.3f Hz" % (rate.mean(), rate.min(), rate.max(), np.median(rate)))
mean 2.546 Hz, min 0.108 Hz, max 33.824 Hz, median 1.387 Hz
Result: reproduced. Mean 2.546 Hz, range 0.108 to 33.824 Hz, against 2.55 and 0.11 to 33.82 Hz.
projection-distance¶
The pairwise projection distance between units isolated on the same wire was 13.77 plus or minus 7.32 s.d. of the signal.
pdist = units_st.proj_dist.to_numpy()
nz = pdist[pdist != 0]
print("non-zero values: %d of %d; mean %.3f, s.d. %.3f, range %.2f to %.2f" % (len(nz), len(pdist), nz.mean(), nz.std(ddof=1), nz.min(), nz.max()))
non-zero values: 734 of 902; mean 13.774, s.d. 7.325, range 3.72 to 54.13
Result: reproduced. Over the 734 non-zero stored values, 13.774 ± 7.325.
waveform-snr¶
The ratio between the s.d. of the noise and the peak amplitude of the mean waveform (an SNR) was 6.62 plus or minus 3.63 across Sternberg units.
for col in ("peak_snr", "mean_snr"):
v = units_st[col]
print("%s: mean %.3f, s.d. %.3f, range %.2f to %.2f, n = %d" % (col, v.mean(), v.std(), v.min(), v.max(), v.notna().sum()))
peak_snr: mean 6.625, s.d. 3.632, range 2.65 to 39.32, n = 902 mean_snr: mean 1.647, s.d. 0.954, range 0.42 to 7.90, n = 902
Result: reproduced. The stored peak SNR (reading a) is 6.625 ± 3.632; the mean SNR (reading b) is 1.647 ± 0.954 and does not match.
cv2¶
The modified coefficient of variation of the ISIs (CV2) was 0.94 plus or minus 0.14 across Sternberg units.
print("CV2: mean %.3f, s.d. %.3f, range %.3f to %.3f, undefined %d" % (np.nanmean(cv2s), np.nanstd(cv2s, ddof=1), np.nanmin(cv2s), np.nanmax(cv2s), np.isnan(cv2s).sum()))
CV2: mean 0.938, s.d. 0.138, range 0.253 to 1.924, undefined 0
Result: reproduced. CV2 is 0.938 ± 0.138.
isolation-distance¶
The median isolation distance of Sternberg units was 31.14.
iso = units_st.isolation_distance.dropna()
print("defined for %d of %d units; median %.3f, range %.2f to %.2f" % (len(iso), len(units_st), iso.median(), iso.min(), iso.max()))
defined for 676 of 902 units; median 31.138, range 4.15 to 7928.56
Result: reproduced. Median 31.14 over the 676 units with a defined value.
Selective cells¶
Tests follow the Methods ("Cell selection"): permutation tests with 2000 permutations. The code's parametric tests are shown as a variant. Each claim's section has its own random generator.
screening-concept-mtl¶
In the screening task, 138 MTL units (32.78%) qualified as concept cells.
rng_sc = np.random.default_rng(469_01)
conc_sc = sb.screening_concept(sc, rng_sc)
conc_sc["MTL"] = conc_sc.area.isin(sb.MTL)
conc_sc["MFC"] = conc_sc.area.isin(sb.MFC)
def summary(df, group):
g = df[df[group]]
return {"units": len(g), **{c: (int(g[c].sum()), round(100 * g[c].mean(), 2))
for c in ("concept_two", "concept_one", "parametric")}}
print("columns: (count, percent); concept_two = permutation tests, second test two-sided; "
"concept_one = one-sided, larger; parametric = the code's anovan tests")
print("MTL:", summary(conc_sc, "MTL"))
print("MTL count as a percent of 449 (the text's MTL count):", {c: round(100 * conc_sc[conc_sc.MTL][c].sum() / 449, 2)
for c in ("concept_two", "concept_one", "parametric")})
print("units passing the permuted ANOVA:", int((conc_sc[conc_sc.MTL].p_anova < 0.05).sum()), "MTL")
print("per subject (MTL units, concept_two):")
print(conc_sc[conc_sc.MTL].groupby("subject").agg(units=("concept_two", "size"), concept=("concept_two", "sum")).T.to_string())
columns: (count, percent); concept_two = permutation tests, second test two-sided; concept_one = one-sided, larger; parametric = the code's anovan tests
MTL: {'units': 421, 'concept_two': (140, np.float64(33.25)), 'concept_one': (140, np.float64(33.25)), 'parametric': (138, np.float64(32.78))}
MTL count as a percent of 449 (the text's MTL count): {'concept_two': np.float64(31.18), 'concept_one': np.float64(31.18), 'parametric': np.float64(30.73)}
units passing the permuted ANOVA: 140 MTL
per subject (MTL units, concept_two):
subject 1 2 3 4 5 7 8 9 10 11 12 13 14 15 16 17 18 20 21
units 9 11 7 43 14 21 25 15 36 47 4 59 41 15 34 2 27 1 10
concept 3 1 1 13 2 7 5 0 6 14 2 32 24 1 9 0 15 0 5
Result: reproduced. Of 421 MTL screening units, 140 (33.25%) are concept cells under both readings of the second test (two-sided and one-sided give the same cells); both values are within tolerance of 138 and 32.78%. The code's parametric tests give 138 (32.78%), the paper's numbers exactly.
screening-concept-mfc¶
In the screening task only 26 MFC units (5.35%) were concept cells, far fewer than in the MTL, so that concept cells are largely found in the MTL.
print("MFC:", summary(conc_sc, "MFC"))
print("MTL percent (two-sided) %.2f%% against MFC %.2f%%" % (100 * conc_sc[conc_sc.MTL].concept_two.mean(),
100 * conc_sc[conc_sc.MFC].concept_two.mean()))
# The MFC count is near the upper bound of the tolerance (28.6 units, 5.885%): repeat with three other seeds.
mfc_sessions = [dict(s, units=s["units"][s["units"].area.isin(sb.MFC)]) for s in sc]
for seed in (469_11, 469_12, 469_13):
rep = sb.screening_concept(mfc_sessions, np.random.default_rng(seed))
print(f"seed {seed}: MFC concept cells two-sided {rep.concept_two.sum()} ({100 * rep.concept_two.mean():.2f}%), "
f"one-sided {rep.concept_one.sum()} ({100 * rep.concept_one.mean():.2f}%)")
# Variant: units whose p lies between 0.02 and 0.10 for either test re-tested with 20000 permutations,
# to reduce the Monte Carlo error of 2000 permutations (same test, finer resolution).
refined = sb.screening_concept_refined(sc, conc_sc.drop(columns=["MTL", "MFC"]), np.random.default_rng(469_14))
for group, areas_ in (("MTL", sb.MTL), ("MFC", sb.MFC)):
g = refined[refined.area.isin(areas_)]
print(f"20000 permutations near 0.05 ({int(g.refined.sum())} {group} units re-tested): {group} concept cells "
f"two-sided {int(g.concept_two.sum())} ({100 * g.concept_two.mean():.2f}%), "
f"one-sided {int(g.concept_one.sum())} ({100 * g.concept_one.mean():.2f}%)")
print("per area (two-sided):", conc_sc.groupby("area").concept_two.agg(["sum", "size"]).T.to_dict())
MFC: {'units': 486, 'concept_two': (27, np.float64(5.56)), 'concept_one': (27, np.float64(5.56)), 'parametric': (26, np.float64(5.35))}
MTL percent (two-sided) 33.25% against MFC 5.56%
seed 46911: MFC concept cells two-sided 29 (5.97%), one-sided 29 (5.97%)
seed 46912: MFC concept cells two-sided 28 (5.76%), one-sided 28 (5.76%)
seed 46913: MFC concept cells two-sided 27 (5.56%), one-sided 27 (5.56%)
20000 permutations near 0.05 (132 MTL units re-tested): MTL concept cells two-sided 139 (33.02%), one-sided 139 (33.02%)
20000 permutations near 0.05 (142 MFC units re-tested): MFC concept cells two-sided 26 (5.35%), one-sided 26 (5.35%)
per area (two-sided): {'Amg': {'sum': 102, 'size': 262}, 'Hipp': {'sum': 38, 'size': 159}, 'dACC': {'sum': 15, 'size': 221}, 'pre-SMA': {'sum': 9, 'size': 236}, 'vmPFC': {'sum': 3, 'size': 29}}
Result: reproduced. Of 486 MFC screening units, 27 (5.56%) are concept cells under both readings, within tolerance of 26 and 5.35%, and far below the MTL's 33.25%. The count sits near the upper bound of the tolerance (28.6): with three other seeds it is 29, 28 and 27, and one of them (29, 5.97%) is outside. Re-testing the 142 MFC units with a p-value between 0.02 and 0.10 with 20000 permutations, which reduces the Monte Carlo error of the same test, gives 26 (5.35%), and 139 (33.02%) in the MTL. The code's parametric tests give 26 (5.35%), the paper's numbers.
fig4-screening-examples¶
Five example screening units (Fig. 4a to e) are concept cells.
examples_sc = [(4, 16), (7, 24), (15, 3), (16, 22), (21, 15)] # NWB_SB_plotCell_Screening.m, 1-based
rows = [conc_sc[(conc_sc.subject == s_) & (conc_sc.unit_id == u - 1)] for s_, u in examples_sc]
print(pd.concat(rows)[["subject", "unit_id", "area", "pref", "p_anova", "p_t_two", "p_t_one",
"concept_two", "concept_one", "parametric"]].to_string(index=False))
subject unit_id area pref p_anova p_t_two p_t_one concept_two concept_one parametric
4 15 Amg 37 0.004998 0.0005 0.0005 True True True
7 23 Amg 30 0.000500 0.0005 0.0005 True True True
15 2 Amg 53 0.027986 0.0010 0.0010 True True True
16 21 Amg 41 0.000500 0.0005 0.0005 True True True
21 14 Hipp 41 0.005497 0.0005 0.0005 True True True
Result: reproduced. All five units pass the concept cell test under both readings (permuted ANOVA p from 0.0005 to 0.028, second test p at most 0.001), and under the code's tests.
sternberg-concept-mtl¶
In the Sternberg task, 94 MTL units (20.94%) qualified as concept cells during the encoding period.
rng_st = np.random.default_rng(469_02)
conc_st = sb.sternberg_concept(st, rng_st)
conc_st["MTL"] = conc_st.area.isin(sb.MTL)
conc_st["MFC"] = conc_st.area.isin(sb.MFC)
print("MTL:", summary(conc_st, "MTL"))
print("all areas:", summary(conc_st.assign(all=True), "all"))
print("per subject (MTL units, concept_two):")
print(conc_st[conc_st.MTL].groupby("subject").agg(units=("concept_two", "size"), concept=("concept_two", "sum")).T.to_string())# Repeat with other seeds: the one-sided count is near the upper bound of the tolerance (103.4 units).
for seed in (469_31, 469_32, 469_33):
rep = sb.sternberg_concept(st, np.random.default_rng(seed))
rep_mtl = rep[rep.area.isin(sb.MTL)]
print(f"seed {seed}:", {k: (int(rep_mtl[k].sum()), round(100 * rep_mtl[k].mean(), 2),
abs(rep_mtl[k].sum() - 94) <= 9.4 and abs(100 * rep_mtl[k].mean() - 20.94) <= 2.094)
for k in ("concept_two", "concept_one")})
MTL: {'units': 449, 'concept_two': (96, np.float64(21.38)), 'concept_one': (100, np.float64(22.27)), 'parametric': (94, np.float64(20.94))}
all areas: {'units': 902, 'concept_two': (123, np.float64(13.64)), 'concept_one': (128, np.float64(14.19)), 'parametric': (119, np.float64(13.19))}
per subject (MTL units, concept_two):
subject 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
units 8 18 5 42 19 1 16 27 40 12 26 4 36 36 35 33 6 64 11 1 9
concept 0 3 1 8 0 0 6 4 7 0 6 3 17 20 2 10 0 6 1 0 2
seed 46931: {'concept_two': (95, np.float64(21.16), np.True_), 'concept_one': (99, np.float64(22.05), np.True_)}
seed 46932: {'concept_two': (97, np.float64(21.6), np.True_), 'concept_one': (100, np.float64(22.27), np.True_)}
seed 46933: {'concept_two': (96, np.float64(21.38), np.True_), 'concept_one': (99, np.float64(22.05), np.True_)}
Result: reproduced. Of 449 MTL Sternberg units, 96 (21.38%) are concept cells with a two-sided second test and 100 (22.27%) with a one-sided one; both are within tolerance of 94 and 20.94%. The code's parametric tests give 94 (20.94%), the paper's numbers. With three other seeds the counts are 95 to 97 (two-sided) and 99 to 100 (one-sided), all within tolerance.
persistent-activity¶
As a group, concept cells fired more during the maintenance period when their preferred image was held in memory than when it was not (4.76 against 2.67 Hz).
Readings, from the criterion: cells (a) MTL concept cells or (b) concept cells in any area; preferred trials (i) the preferred image shown in any encoding period or (ii) as the Encoding 1 image. The concept cells are those of sternberg-concept-mtl, under both readings of its second test. Mean and s.d. over cells of the maintenance rate (0 to 2500 ms), and a paired one-sided t-test. The variant is the same comparison with the Encoding 1 window rates (200 to 1000 ms), as the code computes it with the intended mask.
expected = (4.76, 4.71, 2.67, 3.66)
tol = (0.476, 0.471, 0.267, 0.366)
def persistent_table(conc):
rows = []
for test in ("concept_two", "concept_one", "parametric"):
for cells in ("MTL", "all"):
g = conc[conc[test] & (conc.MTL if cells == "MTL" else True)]
for window in ("maint", "enc1"):
for trials in ("any", "enc1"):
a, b = g[f"{window}_pref_{trials}"], g[f"{window}_nonpref_{trials}"]
vals = (a.mean(), a.std(), b.mean(), b.std())
rows.append({"test": test, "cells": cells, "window": window, "pref trials": trials, "n": len(g),
"pref mean": vals[0], "pref sd": vals[1], "nonpref mean": vals[2], "nonpref sd": vals[3],
"p (paired, one-sided)": stats.ttest_rel(a, b, alternative="greater").pvalue,
"cells pref > nonpref": int((a > b).sum()),
"all four within tolerance": all(abs(v - e) <= t for v, e, t in zip(vals, expected, tol))})
return pd.DataFrame(rows)
pa = persistent_table(conc_st)
print("window maint = 0-2500 ms after maintenance onset (the criterion); enc1 = 200-1000 ms after Encoding 1 onset (variant)")
print(pa.round({c: 3 for c in pa.columns if c.startswith(("pref ", "nonpref "))}).to_string(index=False, formatters={"p (paired, one-sided)": "{:.1e}".format}))
window maint = 0-2500 ms after maintenance onset (the criterion); enc1 = 200-1000 ms after Encoding 1 onset (variant)
test cells window pref trials n pref mean pref sd nonpref mean nonpref sd p (paired, one-sided) cells pref > nonpref all four within tolerance
concept_two MTL maint any 96 2.567 3.768 2.116 3.671 3.9e-08 83 False
concept_two MTL maint enc1 96 2.607 3.812 2.219 3.670 2.8e-08 78 False
concept_two MTL enc1 any 96 3.223 4.129 2.173 3.660 3.7e-16 95 False
concept_two MTL enc1 enc1 96 4.312 4.869 2.166 3.635 2.9e-16 96 False
concept_two all maint any 123 2.943 3.659 2.566 3.619 2.7e-08 98 False
concept_two all maint enc1 123 2.991 3.694 2.649 3.608 3.8e-09 97 False
concept_two all enc1 any 123 3.609 4.013 2.585 3.593 1.5e-21 122 False
concept_two all enc1 enc1 123 4.631 4.644 2.588 3.578 8.6e-22 123 True
concept_one MTL maint any 100 2.743 4.064 2.324 4.054 1.7e-07 84 False
concept_one MTL maint enc1 100 2.774 4.062 2.421 4.045 6.3e-07 80 False
concept_one MTL enc1 any 100 3.381 4.285 2.313 3.746 3.6e-16 99 False
concept_one MTL enc1 enc1 100 4.460 4.989 2.312 3.745 5.7e-17 100 False
concept_one all maint any 128 3.083 3.880 2.734 3.903 1.3e-07 99 False
concept_one all maint enc1 128 3.126 3.878 2.811 3.886 7.5e-08 100 False
concept_one all enc1 any 128 3.735 4.121 2.701 3.645 1.8e-21 127 False
concept_one all enc1 enc1 128 4.748 4.726 2.710 3.648 1.3e-22 128 True
parametric MTL maint any 94 2.657 3.816 2.204 3.732 6.2e-08 80 False
parametric MTL maint enc1 94 2.698 3.860 2.307 3.728 4.4e-08 75 False
parametric MTL enc1 any 94 3.359 4.200 2.279 3.740 2.0e-16 93 False
parametric MTL enc1 enc1 94 4.466 4.924 2.274 3.718 2.4e-16 94 False
parametric all maint any 119 2.997 3.710 2.622 3.682 7.9e-08 93 False
parametric all maint enc1 119 3.047 3.746 2.704 3.667 9.8e-09 92 False
parametric all enc1 any 119 3.722 4.087 2.668 3.669 1.4e-21 118 False
parametric all enc1 enc1 119 4.764 4.709 2.673 3.656 1.4e-21 119 True
# Variant: the three trials flagged in the data check left out, concept test and rates recomputed (same seed)
def flagged(t):
fix = t.timestamps_Encoding1 - t.timestamps_FixationCross
return ((t.timestamps_Probe - t.timestamps_Maintenance) < 2) | (fix > 2) | (fix < 0.5)
st_clean = [dict(s, trials=s["trials"][~flagged(s["trials"])].reset_index(drop=True)) for s in st]
print("trials left out:", sum(int(flagged(s["trials"]).sum()) for s in st), "of", sum(len(s["trials"]) for s in st))
conc_clean = sb.sternberg_concept(st_clean, np.random.default_rng(469_02))
conc_clean["MTL"] = conc_clean.area.isin(sb.MTL)
print("MTL concept cells without them:", summary(conc_clean, "MTL"), "| all areas:", summary(conc_clean.assign(all=True), "all"))
pa_clean = persistent_table(conc_clean)
print(pa_clean.round({c: 3 for c in pa_clean.columns if c.startswith(("pref ", "nonpref "))})
.to_string(index=False, formatters={"p (paired, one-sided)": "{:.1e}".format}))
trials left out: 3 of 2673
MTL concept cells without them: {'units': 449, 'concept_two': (97, np.float64(21.6)), 'concept_one': (102, np.float64(22.72)), 'parametric': (94, np.float64(20.94))} | all areas: {'units': 902, 'concept_two': (124, np.float64(13.75)), 'concept_one': (130, np.float64(14.41)), 'parametric': (119, np.float64(13.19))}
test cells window pref trials n pref mean pref sd nonpref mean nonpref sd p (paired, one-sided) cells pref > nonpref all four within tolerance
concept_two MTL maint any 97 2.567 3.748 2.125 3.653 5.4e-08 83 False
concept_two MTL maint enc1 97 2.608 3.793 2.226 3.652 3.4e-08 78 False
concept_two MTL enc1 any 97 3.227 4.108 2.176 3.641 2.0e-16 96 False
concept_two MTL enc1 enc1 97 4.311 4.847 2.170 3.616 2.0e-16 97 False
concept_two all maint any 124 2.940 3.645 2.570 3.605 3.7e-08 98 False
concept_two all maint enc1 124 2.988 3.679 2.651 3.594 4.4e-09 97 False
concept_two all enc1 any 124 3.609 3.998 2.585 3.579 8.4e-22 123 False
concept_two all enc1 enc1 124 4.627 4.629 2.589 3.564 6.0e-22 124 True
concept_one MTL maint any 102 2.724 4.028 2.319 4.016 2.7e-07 84 False
concept_one MTL maint enc1 102 2.753 4.027 2.413 4.007 1.0e-06 80 False
concept_one MTL enc1 any 102 3.362 4.249 2.303 3.710 2.1e-16 101 False
concept_one MTL enc1 enc1 102 4.430 4.951 2.303 3.709 3.7e-17 102 False
concept_one all maint any 130 3.063 3.854 2.724 3.875 2.0e-07 99 False
concept_one all maint enc1 130 3.104 3.854 2.799 3.859 1.2e-07 100 False
concept_one all enc1 any 130 3.715 4.095 2.688 3.620 1.1e-21 129 False
concept_one all enc1 enc1 130 4.719 4.701 2.697 3.622 8.5e-23 130 True
parametric MTL maint any 94 2.656 3.816 2.204 3.732 6.6e-08 80 False
parametric MTL maint enc1 94 2.695 3.860 2.308 3.728 5.5e-08 74 False
parametric MTL enc1 any 94 3.361 4.200 2.280 3.741 2.0e-16 93 False
parametric MTL enc1 enc1 94 4.471 4.927 2.275 3.718 2.4e-16 94 False
parametric all maint any 119 2.996 3.710 2.623 3.682 8.3e-08 93 False
parametric all maint enc1 119 3.045 3.746 2.704 3.667 1.2e-08 91 False
parametric all enc1 any 119 3.724 4.088 2.668 3.669 1.4e-21 118 False
parametric all enc1 enc1 119 4.768 4.712 2.674 3.656 1.4e-21 119 True
Result: not reproduced. In the maintenance window, under every cell and trial reading, concept cells fire more on preferred trials (MTL cells, two-sided test, 2.57 against 2.12 Hz with the preferred image in any encoding period, and 2.61 against 2.22 Hz with it as the Encoding 1 image; all cells 2.94 against 2.57 and 2.99 against 2.65 Hz; every p below 1e-5), but the preferred means are 2.57 to 3.13 Hz, far below the paper's 4.76 Hz, and no reading puts the four values within tolerance. The paper's values come out exactly from the variant: all 119 concept cells of the code's parametric test, Encoding 1 window, preferred as the Encoding 1 image, give 4.76 ± 4.71 against 2.67 ± 3.66 Hz. With the permutation concept cells the same variant gives 4.63 ± 4.64 against 2.59 ± 3.58 Hz (123 cells) and 4.75 ± 4.73 against 2.71 ± 3.65 Hz (128 cells), also within tolerance. The paper's numbers therefore appear to describe the encoding response, not maintenance activity. Leaving out the three trials flagged in the data check and recomputing the concept cells (97 and 102 MTL cells, 94 parametric) changes none of this: the maintenance-window preferred means are 2.57 to 3.10 Hz, no maintenance reading is within tolerance, and the parametric Encoding 1 variant gives 4.77 ± 4.71 against 2.67 ± 3.66 Hz.
fig5-sternberg-examples¶
Two example amygdala units (Fig. 5a,b) are concept cells in the Sternberg task that also show persistent activity during maintenance.
examples_st = [(7, 9), (16, 16)] # NWB_SB_plotCell_Sternberg.m, 1-based
rows = [conc_st[(conc_st.subject == s_) & (conc_st.unit_id == u - 1)] for s_, u in examples_st]
print(pd.concat(rows)[["subject", "unit_id", "area", "pref", "p_anova", "p_t_two", "p_t_one", "concept_two",
"concept_one", "parametric", "maint_pref_any", "maint_nonpref_any"]].round(4).to_string(index=False))
subject unit_id area pref p_anova p_t_two p_t_one concept_two concept_one parametric maint_pref_any maint_nonpref_any
7 8 Amg 5 0.0005 0.0005 0.0005 True True True 5.0512 1.4277
16 15 Amg 5 0.0005 0.0005 0.0005 True True True 5.7333 4.1697
Result: reproduced. Both units are amygdala concept cells under both readings (and the code's tests), and both fire more during maintenance on trials that showed the preferred image (5.05 against 1.43 Hz, and 5.73 against 4.17 Hz).
maintenance-cells-mfc¶
In the MFC, 144 units (31.79%) qualified as maintenance cells.
Readings, from the criterion: permutation by (a) label shuffling or (b) sign flipping; for concept cells, non-preferred trials (i) Encoding 1 image not preferred or (ii) preferred image not shown in any encoding period. Concept cells by the two-sided and one-sided tests of sternberg-concept-mtl. Variant: the code's paired one-sided t-tests.
rng_maint = np.random.default_rng(469_03)
maint = sb.sternberg_maintenance(st, conc_st, rng_maint)
maint = maint.merge(conc_st[["subject", "unit_id", "concept_two", "concept_one", "parametric"]], on=["subject", "unit_id"])
mfc = maint.area.isin(sb.MFC)
# The same permutations with ties counted as exceedances (>=), a variant
maint_ge = sb.sternberg_maintenance(st, conc_st, np.random.default_rng(469_03), ties=True)
maint_ge = maint_ge.merge(conc_st[["subject", "unit_id", "concept_two", "concept_one", "parametric"]], on=["subject", "unit_id"])
rows = []
for perm in ("shuffle", "flip", "t", "shuffle_ge", "flip_ge"):
for nonpref in ("enc1", "any"):
for test in ("concept_two", "concept_one", "parametric"):
df, pm = (maint_ge, perm[:-3]) if perm.endswith("_ge") else (maint, perm)
is_m = (df[f"p_{pm}"] < 0.05) & (~df[test] | (df[f"p2_{pm}_{nonpref}"] < 0.05))
n = int(is_m[mfc].sum())
rows.append({"permutation": perm, "non-preferred": nonpref, "concept test": test, "MFC maintenance": n,
"percent": round(100 * n / mfc.sum(), 2),
"within tolerance": abs(n - 144) <= 14.4 and abs(100 * n / mfc.sum() - 31.79) <= 3.179})
print("MFC units:", int(mfc.sum()), "(t = the code's paired t-tests, a variant; _ge = ties counted, a variant)")
print(pd.DataFrame(rows).to_string(index=False))
main = (maint.p_shuffle < 0.05) & (~maint.concept_two | (maint.p2_shuffle_enc1 < 0.05))
print("per area, label shuffling, Encoding 1 non-preferred:", maint[main].area.value_counts().to_dict())
print("per subject (MFC):", maint[mfc & main].groupby("subject").size().to_dict())
MFC units: 453 (t = the code's paired t-tests, a variant; _ge = ties counted, a variant)
permutation non-preferred concept test MFC maintenance percent within tolerance
shuffle enc1 concept_two 140 30.91 True
shuffle enc1 concept_one 140 30.91 True
shuffle enc1 parametric 140 30.91 True
shuffle any concept_two 139 30.68 True
shuffle any concept_one 139 30.68 True
shuffle any parametric 139 30.68 True
flip enc1 concept_two 144 31.79 True
flip enc1 concept_one 144 31.79 True
flip enc1 parametric 144 31.79 True
flip any concept_two 143 31.57 True
flip any concept_one 143 31.57 True
flip any parametric 143 31.57 True
t enc1 concept_two 144 31.79 True
t enc1 concept_one 144 31.79 True
t enc1 parametric 144 31.79 True
t any concept_two 143 31.57 True
t any concept_one 143 31.57 True
t any parametric 143 31.57 True
shuffle_ge enc1 concept_two 140 30.91 True
shuffle_ge enc1 concept_one 140 30.91 True
shuffle_ge enc1 parametric 140 30.91 True
shuffle_ge any concept_two 139 30.68 True
shuffle_ge any concept_one 139 30.68 True
shuffle_ge any parametric 139 30.68 True
flip_ge enc1 concept_two 142 31.35 True
flip_ge enc1 concept_one 142 31.35 True
flip_ge enc1 parametric 142 31.35 True
flip_ge any concept_two 141 31.13 True
flip_ge any concept_one 141 31.13 True
flip_ge any parametric 141 31.13 True
per area, label shuffling, Encoding 1 non-preferred: {'pre-SMA': 82, 'Amg': 64, 'Hipp': 58, 'dACC': 46, 'vmPFC': 12}
per subject (MFC): {1: 13, 2: 4, 3: 13, 4: 9, 5: 4, 6: 11, 7: 1, 8: 22, 9: 1, 10: 3, 11: 23, 12: 4, 13: 9, 14: 8, 15: 6, 18: 6, 21: 3}
Result: reproduced. Of 453 MFC units, 139 to 144 (30.68% to 31.79%) are maintenance cells
under the four readings and every concept cell set, all within tolerance of 144 and 31.79%. The
code's paired t-tests give 144 (31.79%), the paper's numbers. Counting permuted
values equal to the observed one as exceedances (>=, a variant on the same permutations)
leaves label shuffling at 140 and 139 and lowers sign flipping to 142 and 141, all within
tolerance.
probe-cells-mfc¶
In the MFC, 122 units (26.93%) qualified as probe cells.
Readings, from the criterion: encoding rate (i) pooled over the encoding periods shown or (ii) Encoding 1 only; permutation (a) label shuffling or (b) sign flipping. Variant: the code's paired one-sided t-tests, with its Encoding 3 window at Encoding 2 onset.
def probe_counts(df):
out = {}
for col in ("p_shuffle_pooled", "p_shuffle_enc1", "p_flip_pooled", "p_flip_enc1", "p_code"):
n = int((df[col] < 0.05).sum())
out[col] = (n, round(100 * n / len(df), 2), abs(n - 122) <= 12.2 and abs(100 * n / len(df) - 26.93) <= 2.693)
return out
rng_probe = np.random.default_rng(469_04)
probe = sb.sternberg_probe(st, rng_probe, areas=sb.MFC)
print("MFC units:", len(probe), "| (count, percent, within tolerance) per reading:")
for k, v in probe_counts(probe).items():
print(f" {k}: {v}")
# The label-shuffling count is close to the lower bound of the tolerance (109.8 units): repeat with other seeds.
for seed in (469_21, 469_22, 469_23):
rep = sb.sternberg_probe(st, np.random.default_rng(seed), areas=sb.MFC)
print(f"seed {seed}:", {k: v for k, v in probe_counts(rep).items() if k != "p_code"})
# The same permutations as the main seed with ties counted as exceedances (>=), a variant
probe_ge = sb.sternberg_probe(st, np.random.default_rng(469_04), areas=sb.MFC, ties=True)
print("ties counted:", {k: v for k, v in probe_counts(probe_ge).items() if k != "p_code"})
print("per area (label shuffling, pooled):", probe[probe.p_shuffle_pooled < 0.05].area.value_counts().to_dict())
MFC units: 453 | (count, percent, within tolerance) per reading: p_shuffle_pooled: (109, 24.06, False) p_shuffle_enc1: (104, 22.96, False) p_flip_pooled: (118, 26.05, True) p_flip_enc1: (111, 24.5, True) p_code: (122, 26.93, True)
seed 46921: {'p_shuffle_pooled': (111, 24.5, True), 'p_shuffle_enc1': (105, 23.18, False), 'p_flip_pooled': (120, 26.49, True), 'p_flip_enc1': (114, 25.17, True)}
seed 46922: {'p_shuffle_pooled': (110, 24.28, True), 'p_shuffle_enc1': (105, 23.18, False), 'p_flip_pooled': (124, 27.37, True), 'p_flip_enc1': (110, 24.28, True)}
seed 46923: {'p_shuffle_pooled': (110, 24.28, True), 'p_shuffle_enc1': (104, 22.96, False), 'p_flip_pooled': (122, 26.93, True), 'p_flip_enc1': (113, 24.94, True)}
ties counted: {'p_shuffle_pooled': (108, 23.84, False), 'p_shuffle_enc1': (103, 22.74, False), 'p_flip_pooled': (118, 26.05, True), 'p_flip_enc1': (111, 24.5, True)}
per area (label shuffling, pooled): {'pre-SMA': 90, 'dACC': 18, 'vmPFC': 1}
Result: blocked (methods underspecified). The readings disagree. With sign flipping, 118
(26.05%, pooled encoding) and 111 (24.50%, Encoding 1) units are probe cells, within tolerance.
With label shuffling, 109 (24.06%) and 104 (22.96%) are, outside it. With three other seeds,
label shuffling with pooled encoding gives 110 or 111, just inside the bound of 109.8, while
label shuffling with Encoding 1 gives 104 or 105 every time, outside it; sign flipping stays
within tolerance (110 to 124). Counting ties as exceedances (>=, same permutations as the main
seed) gives 108 and 103 for label shuffling, outside, and 118 and 111 for sign flipping, within
tolerance, so the readings still disagree. The code's paired
t-tests give 122 (26.93%), the paper's numbers. The Methods name "permutation t-tests" without
saying whether the trials are paired, and that choice largely decides the outcome. Under label
shuffling the encoding rate also matters: the pooled rate is within tolerance on 3 of the 4
seeds and Encoding 1 on none.
Summary¶
| Claim | Outcome |
|---|---|
| subjects-sessions | reproduced |
| sternberg-unit-yield | reproduced |
| screening-unit-yield | reproduced (MTL count is 421, not 449) |
| behavior-accuracy | reproduced |
| behavior-rt-load | reproduced |
| behavior-accuracy-load | reproduced |
| isi-violations | reproduced (on the code's filter) |
| firing-rate | reproduced |
| projection-distance | reproduced |
| waveform-snr | reproduced (peak SNR) |
| cv2 | reproduced |
| isolation-distance | reproduced |
| screening-concept-mtl | reproduced |
| screening-concept-mfc | reproduced |
| fig4-screening-examples | reproduced |
| sternberg-concept-mtl | reproduced |
| persistent-activity | not reproduced |
| fig5-sternberg-examples | reproduced |
| maintenance-cells-mfc | reproduced |
| probe-cells-mfc | blocked, methods underspecified |
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 the analysis differs from the paper's: the selection tests are the permutation tests the Methods name, with 2000 permutations, where the paper's code runs parametric tests; the code's tests are shown as variants and give the paper's counts exactly. The number of permutations and the form of the permutation for paired rates are not given in the paper.
What was found about the dataset: it holds every session, unit and trial the paper describes,
and the paper's behaviour numbers, unit counts by area and quality measures are recovered from
it, with three exceptions: the screening MTL count (449) is the Sternberg count, and the files
give 421; the ISI violation s.d. matches only with the code's filter (0.625% over all units);
and the SNR matches the stored peak SNR, not the mean SNR. The paper's
persistent activity values are Encoding 1 window rates, not maintenance rates. Subjects 13 and 14
are P39CS and P39CS_2. Four units have spike times out of order. In one Sternberg trial (subject 9)
Encoding 1 lasts 2.5 s and the maintenance and probe onsets coincide, in another (subject 9)
the fixation cross and Encoding 1 onsets are 0.03 ms apart, so the baseline window falls on the
image, and in one (subject 21) the fixation cross lasts 14.07 s. These trials are kept, as the
code keeps them; leaving them out does not change the persistent activity result. Two column descriptions are wrong
(waveforms_mean_proj_dist is described as the isolation distance, loadsEnc3_PicIDs as Enc1),
and every column of the electrodes table but id is described as "my description".
print("run time of this notebook: %.1f min" % ((time.time() - started) / 60))
run time of this notebook: 4.9 min