Working Memory of Multi-Object Scenes in Primate Frontal Cortex (Watters, Gabel, Tenenbaum & Jazayeri 2026): reproduction on DANDI 000620¶
Paper. N. Watters, J. Gabel, J. Tenenbaum, M. Jazayeri, "Working Memory of Multi-Object Scenes in Primate Frontal Cortex", bioRxiv preprint version 1, doi:10.64898/2026.01.27.702062. The claims file (claims/000620/10.64898-2026.01.27.702062.yaml) was written from the bioRxiv full text (main text, Methods, legends of Extended Data Figures 1 to 11) and the paper's code, https://github.com/jazlab/multi_object_memory_2025 at commit 6e4b333b78cd4461e6373a42378535e8fd392cdb, which this attempt also reads for method details (cited by file in wm.py and below).
Dataset. DANDI 000620, "Multi-object working memory in macaque frontal cortex", version 0.261019.2000 (doi:10.48324/dandi.000620/0.261019.2000), streamed with remfile and h5py. Why this version: the paper cites the dandiset without a version, so the most recent published version is pinned. The four published versions (0.260127.2208, 0.261014.1424, 0.261019.2000 and the draft) have the same 142 assets and the same total size (5,490,107,572,879 bytes), and the release notes of 0.261019.2000 say "No changes to data", so no earlier version matches the claims file's description better. Each session has three files: raw ecephys (5.5 TB in all, not used), spikesorting (units with spike times, quality, probe, stable intervals and per-trial stability) and behavior+task (the trials table, eye position and reward line). This notebook uses the trials tables and the units tables.
What is tested where. This notebook tests the behaviour, unit-count, selectivity and correlation claims. The decoding claims (ext3ab-decoding, ext3c-decoder-orthogonality) are in decoding_watters_2026.ipynb, and variants of its decoders that are not the code's in decoding_variants_watters_2026.ipynb. The claims that need the paper's encoding models fitted (13 claims) are blocked as too_large; the last section times the computation.
| claim | where | status |
|---|---|---|
| fig1b-performance-load | this notebook | tested |
| fig1c-reaction-times | this notebook | tested |
| ext1-misbinding-errors | this notebook | tested |
| triangle-unit-counts | this notebook | tested |
| ring-unit-counts | this notebook | tested |
| ext2-stimulus-delay-activity | this notebook | tested |
| ext3ab-decoding | decoding notebook | tested |
| ext3c-decoder-orthogonality | decoding notebook | tested |
| fig2b-selectivity-range | this notebook | tested |
| fig2d-signal-noise-triangle | this notebook | tested |
| fig5ab-ring-behavior | this notebook | tested |
| fig5cd-ring-tuning | this notebook | tested |
| fig5e-signal-noise-ring | this notebook | tested |
| fig4a-gain-best-triangle | this notebook, timing only | blocked, too_large |
| ext6-gain-best-areas | timing only | blocked, too_large |
| ext5-general-slot | timing only | blocked, too_large |
| fig4b-identifiability | timing only | blocked, too_large |
| nslots-two-best | timing only | blocked, too_large |
| ext11-kappa-tuning | timing only | blocked, too_large |
| ext7-gains-distributed | timing only | blocked, too_large |
| fig4c-gain-performance-triangle | timing only | blocked, too_large |
| fig4d-gain-rt-triangle | timing only | blocked, too_large |
| ext8-slot-switching-behavior | timing only | blocked, too_large |
| fig5f-gain-best-ring | timing only | blocked, too_large |
| fig5g-gain-performance-ring | timing only | blocked, too_large |
| fig5h-gain-rt-ring | timing only | blocked, too_large |
Shared definitions T1 to T8 are those of the claims file's opening comment; wm.py says where each parameter comes from.
import warnings
warnings.filterwarnings("ignore")
import time
import numpy as np, pandas as pd
import wm
t_start = time.time()
assets = wm.list_assets()
print(f"Dandiset {wm.DANDISET} version {wm.VERSION}: {len(assets)} assets,", assets.kind.value_counts().to_dict())
print("Raw ecephys files (not used):", f"{assets[assets.kind == 'ecephys']['size'].sum() / 1e12:.2f} TB")
keys = sorted(set(zip(assets.subject, assets.session)))
S = [wm.load_session(assets, subj, sess) for subj, sess in keys]
print(f"Loaded {len(S)} sessions in {time.time() - t_start:.0f} s")
Dandiset 000620 version 0.261019.2000: 142 assets, {'behavior+task': 48, 'ecephys': 47, 'spikesorting': 47}
Raw ecephys files (not used): 5.47 TB
Loaded 48 sessions in 264 s
import numpy as np, pandas as pd, scipy.stats as st
import matplotlib
import matplotlib.pyplot as plt
pd.set_option("display.width", 200); pd.set_option("display.max_columns", 30)
MONKEY = {"Perle": "P", "Elgar": "E"}
for x in S:
x["trials"]["subject"] = x["subject"]; x["trials"]["session"] = x["session"]
x["trials"]["task"] = x["task"]; x["trials"]["trial"] = np.arange(len(x["trials"]))
# Physiology sessions of each dataset: sessions with spike sorting. Ring sessions are those
# with 1- and 2-object trials (Methods, Behavioral Task); Perle 2022-06-06 has 2-object
# trials only and is kept apart (see the data check).
tri = [x for x in S if x["task"] == "triangle" and x["units"] is not None]
ring = [x for x in S if x["task"] == "ring" and x["units"] is not None
and (x["trials"].n_obj == 1).any()]
other = [x for x in S if x not in tri and x not in ring]
T = pd.concat([x["trials"] for x in S], ignore_index=True)
T["physiology"] = T.session.isin([x["session"] for x in tri + ring])
print("Triangle sessions:", {m: sum(x["subject"] == m for x in tri) for m in MONKEY})
print("Ring sessions:", {m: sum(x["subject"] == m for x in ring) for m in MONKEY})
print("Other sessions:", [(x["subject"], x["session"], x["task"], x["units"] is not None,
sorted(x["trials"].n_obj.unique())) for x in other])
Triangle sessions: {'Perle': 10, 'Elgar': 14}
Ring sessions: {'Perle': 10, 'Elgar': 12}
Other sessions: [('Perle', '2022-06-06', 'ring', True, [np.int64(2)]), ('Perle', '2022-06-08', 'ring', False, [np.int64(1), np.int64(2)])]
Data check¶
Printed below: sessions, trials and units by subject, session and condition; ranges and units of every variable used; the reward against the response error (which correctness radius the task used); and the clock check (location information appears after stimulus onset and in the delay, and not before onset or with shifted labels).
rows = []
for x in S:
tr, u = x["trials"], x["units"]
c = tr[tr.completed]
rows.append(dict(
monkey=MONKEY[x["subject"]], session=x["session"], task=x["task"],
trials=len(tr), completed=len(c),
n1=(tr.n_obj == 1).sum(), n2=(tr.n_obj == 2).sum(), n3=(tr.n_obj == 3).sum(),
on_tri=round(tr.on_triangle.mean(), 3),
units=np.nan if u is None else len(u),
good=np.nan if u is None else (u.quality == "good").sum(),
DMFC=np.nan if u is None else (u.area == "DMFC").sum(),
FEF=np.nan if u is None else (u.area == "FEF").sum()))
inv = pd.DataFrame(rows)
print(inv.to_string(index=False))
print("\nQuality labels:", sorted(set(np.concatenate([x["units"].quality.unique() for x in S if x["units"] is not None]))))
print("Probes:", sorted(set(np.concatenate([x["units"].probe.unique() for x in S if x["units"] is not None]))))
c = T[T.completed]
def rng_(v): return f"min {np.nanmin(v):.4f}, median {np.nanmedian(v):.4f}, max {np.nanmax(v):.4f}"
print("\nPhase durations of completed trials (s):")
print(" stimulus (stimulus to delay onset):", rng_(c.phase_delay_time - c.phase_stimulus_time))
print(" delay (delay to cue onset): ", rng_(c.phase_cue_time - c.phase_delay_time))
print(" cue (cue to response-phase onset): ", rng_(c.phase_response_time - c.phase_cue_time))
print(" reaction time (response_time - cue onset):", rng_(c.reaction_time))
print(" share of completed trials with reaction time < 0.18 s:", round(float((c.reaction_time < 0.18).mean()), 4))
print(" share with reaction time < cue duration (0.25 s):", round(float((c.reaction_time < 0.25).mean()), 4))
r_obj = np.hypot(T[["x0", "x1", "x2"]].values, T[["y0", "y1", "y2"]].values)
print("\nObject eccentricity (arena widths):", rng_(r_obj), "; expected 0.35 = 8.4 deg / 24 deg")
print("Response position x, y (arena widths from centre):", rng_(c.resp_x), "|", rng_(c.resp_y))
print("Identities:", sorted(set(T.id0.dropna()) | set(T.id1.dropna()) | set(T.id2.dropna())))
print("delay_object_blanks true:", int(T.delay_object_blanks.sum()), "; nonzero object velocities:", int((T.max_speed > 0).sum()),
"; trials without exactly one target:", int((T.n_target_flags != 1).sum()))
tt = T[(T.task == "triangle") & T.on_triangle]
print("Distinct triangle object locations on on-triangle trials:", sorted(set(tt[["loc0","loc1","loc2"]].values.ravel()) - {-9}))
print("Units: spikes unsorted in", sum(x["units"].attrs["n_unsorted_units"] for x in S if x["units"] is not None),
"units; spikes outside the stable intervals:", sum(x["units"].attrs["n_spikes_outside_intervals"] for x in S if x["units"] is not None))
allrates = np.concatenate([np.nanmean(x["rates"]["delay"], 1) for x in tri + ring])
print("Unit mean delay-phase rate (Hz):", rng_(allrates), "; units with no completed stable trial:", int(np.isnan(allrates).sum()))
print("closed-loop minus post-hoc response time (s):", rng_(c.closed_loop_response_time - c.response_time))
# Clock check: for each physiology session, the share of units whose rate differs across the
# three object locations on 1-object trials (Kruskal-Wallis p < 0.01), before stimulus onset
# (-300 to 0 ms), early in the stimulus (50 to 350 ms) and in the delay, and the last two with
# the labels rolled by one trial (a control that keeps the label counts).
clk = []
for x in tri + ring:
tr = x["trials"]
m = (tr.n_obj == 1).values & tr.completed.values
if x["task"] == "triangle":
m &= tr.on_triangle.values
lab = tr.loc0.values
else:
lab = wm.ring_sector(tr.theta0.values)
res = {}
for w, shift in [("pre", 0), ("early", 0), ("early", 1), ("delay", 0), ("delay", 1)]:
R = x["rates"][w][:, m]
L = np.roll(lab[m], shift)
sig = []
for r in R:
ok = np.isfinite(r)
g = [r[ok & (L == k)] for k in range(3)]
if min(len(v) for v in g) < 5 or np.ptp(r[ok]) == 0:
continue
sig.append(st.kruskal(*g).pvalue < 0.01)
res[f"{w}{'_shifted' if shift else ''}"] = np.mean(sig) if sig else np.nan
clk.append(dict(monkey=MONKEY[x["subject"]], session=x["session"], task=x["task"], **res))
clk = pd.DataFrame(clk)
print("\nClock check, share of units location-modulated (Kruskal-Wallis p < 0.01):")
print(clk.round(3).to_string(index=False))
print("Sessions where early modulation is not above both controls:",
clk[(clk.early <= clk.pre) | (clk.early <= clk.early_shifted)].session.tolist())
print("Sessions where delay modulation is not above the shifted control:",
clk[clk.delay <= clk.delay_shifted].session.tolist())
# Correctness radius (T2): which radius did the task reward? Reward against the closed-loop
# response error (the position the task used for reward), completed trials of all sessions.
c = T[T.completed]
bins = [0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 10]
tab = c.groupby(pd.cut(c.target_err_closed_loop, bins, right=False), observed=True).rewarded.agg(["size", "mean"])
tab.columns = ["completed trials", "share rewarded"]
print("\nShare of completed trials rewarded, by closed-loop distance to the target (arena widths; 0.1 = 2.4 deg, 0.2 = 4.8 deg):")
print(tab.round(4).to_string())
print("\nElectrode locations named in the spikesorting files:",
sorted(set(l for x in S if x["units"] is not None for l in x["units"].attrs["electrode_locations"])),
"; the area of a unit is taken from its probe (s0 = DMFC Neuropixels, vprobe0/1 = FEF V-probes; claims file, code unit_filters.py).")
# Recordings deposited twice: compare trial timing (relative to the first trial) and the
# sorted spike counts of the units between every pair of sessions.
import hashlib
fp = {}
for x in S:
tr = x["trials"]
k = ("trials", hashlib.md5(np.round(tr.phase_stimulus_time.values - tr.phase_stimulus_time.values[0], 4).tobytes()).hexdigest())
fp.setdefault(k, []).append(x["session"])
if x["units"] is not None:
k = ("units", hashlib.md5(np.sort(x["units"].n_spikes.values).tobytes()).hexdigest())
fp.setdefault(k, []).append(x["session"])
print("Sessions sharing trial timing or sorted spike counts:", [v for v in fp.values() if len(v) > 1])
monkey session task trials completed n1 n2 n3 on_tri units good DMFC FEF
E 2022-08-19 triangle 2847 1253 234 928 1685 0.900 323.0 123.0 249.0 74.0
E 2022-08-20 triangle 2619 1361 251 991 1377 0.908 414.0 99.0 248.0 166.0
E 2022-08-21 triangle 2928 1517 289 992 1647 0.913 229.0 63.0 156.0 73.0
E 2022-08-22 triangle 1964 1157 187 722 1055 0.916 234.0 118.0 164.0 70.0
E 2022-08-23 triangle 1816 1299 149 631 1036 0.888 217.0 55.0 176.0 41.0
E 2022-08-24 triangle 1401 849 130 501 770 0.897 283.0 59.0 205.0 78.0
E 2022-08-25 triangle 2083 1353 205 768 1110 0.916 164.0 48.0 137.0 27.0
E 2022-08-26 triangle 2132 1514 240 845 1047 0.916 141.0 50.0 94.0 47.0
E 2022-08-31 triangle 1578 1238 154 595 829 0.911 206.0 83.0 151.0 55.0
E 2022-09-01 triangle 1947 1373 189 720 1038 0.907 225.0 65.0 147.0 78.0
E 2022-09-02 triangle 1512 1193 148 574 790 0.921 151.0 50.0 151.0 0.0
E 2022-09-03 triangle 2267 1735 240 931 1096 0.921 203.0 37.0 148.0 55.0
E 2022-09-04 triangle 1954 1336 200 751 1003 0.908 290.0 85.0 160.0 130.0
E 2022-09-05 triangle 1598 1238 172 621 805 0.913 181.0 61.0 77.0 104.0
E 2022-09-06 ring 2080 1335 389 1691 0 0.003 197.0 44.0 118.0 79.0
E 2022-09-07 ring 2307 1297 437 1870 0 0.003 78.0 13.0 78.0 0.0
E 2022-09-09 ring 2186 1656 481 1705 0 0.004 127.0 34.0 127.0 0.0
E 2022-09-10 ring 1828 1448 376 1452 0 0.003 165.0 42.0 59.0 106.0
E 2022-09-11 ring 2435 1722 488 1947 0 0.004 148.0 70.0 112.0 36.0
E 2022-09-12 ring 2008 1402 428 1580 0 0.003 134.0 38.0 84.0 50.0
E 2022-09-13 ring 1922 1361 387 1535 0 0.002 71.0 19.0 0.0 71.0
E 2022-09-14 ring 2019 1350 414 1605 0 0.000 98.0 30.0 0.0 98.0
E 2022-09-15 ring 1849 1339 384 1465 0 0.004 93.0 19.0 0.0 93.0
E 2022-09-18 ring 1722 1323 315 1407 0 0.001 77.0 11.0 62.0 15.0
E 2022-09-19 ring 1598 1214 322 1276 0 0.003 114.0 43.0 114.0 0.0
E 2022-09-20 ring 1680 1242 306 1374 0 0.002 63.0 14.0 39.0 24.0
P 2022-05-26 triangle 1756 1002 241 587 928 1.000 101.0 42.0 75.0 26.0
P 2022-05-27 triangle 1629 1239 200 602 827 1.000 111.0 35.0 30.0 81.0
P 2022-05-28 triangle 1840 1412 198 711 931 1.000 137.0 54.0 73.0 64.0
P 2022-05-29 triangle 2040 1326 275 768 997 0.904 129.0 53.0 28.0 101.0
P 2022-05-30 triangle 1555 992 223 545 787 0.920 185.0 76.0 85.0 100.0
P 2022-05-31 triangle 1481 957 172 555 754 0.911 200.0 112.0 113.0 87.0
P 2022-06-01 triangle 1621 1097 209 648 764 0.915 226.0 121.0 154.0 72.0
P 2022-06-03 triangle 2091 1335 270 791 1030 0.905 198.0 80.0 105.0 93.0
P 2022-06-04 triangle 1862 1123 234 687 941 0.924 243.0 80.0 86.0 157.0
P 2022-06-05 triangle 1338 896 193 507 638 0.911 108.0 27.0 18.0 90.0
P 2022-06-06 ring 2378 1559 0 2378 0 0.000 77.0 24.0 6.0 71.0
P 2022-06-07 ring 2660 1481 538 2122 0 0.002 14.0 7.0 14.0 0.0
P 2022-06-08 ring 1517 1001 281 1236 0 0.001 NaN NaN NaN NaN
P 2022-06-09 ring 2417 1416 488 1929 0 0.002 193.0 54.0 75.0 118.0
P 2022-06-10 ring 1476 1206 295 1181 0 0.002 216.0 107.0 110.0 106.0
P 2022-06-11 ring 1474 1170 280 1194 0 0.001 207.0 110.0 107.0 100.0
P 2022-06-12 ring 1724 1275 330 1394 0 0.003 231.0 106.0 124.0 107.0
P 2022-06-13 ring 1365 968 257 1108 0 0.001 232.0 104.0 148.0 84.0
P 2022-06-14 ring 1762 965 347 1415 0 0.002 182.0 65.0 104.0 78.0
P 2022-06-15 ring 2291 1516 465 1826 0 0.001 136.0 72.0 100.0 36.0
P 2022-06-16 ring 1372 968 283 1089 0 0.002 217.0 104.0 90.0 127.0
P 2022-06-19 ring 3678 1645 748 2930 0 0.004 193.0 75.0 97.0 96.0
Quality labels: ['good', 'mua']
Probes: ['s0', 'vprobe0', 'vprobe1']
Phase durations of completed trials (s):
stimulus (stimulus to delay onset): min 0.9997, median 1.0000, max 1.4167
delay (delay to cue onset): min 0.9998, median 1.0000, max 1.3333
cue (cue to response-phase onset): min 0.2498, median 0.2500, max 0.3167
reaction time (response_time - cue onset): min -0.0125, median 0.2531, max 1.6594
share of completed trials with reaction time < 0.18 s: 0.0109
share with reaction time < cue duration (0.25 s): 0.4786
Object eccentricity (arena widths): min 0.3500, median 0.3500, max 0.3500 ; expected 0.35 = 8.4 deg / 24 deg
Response position x, y (arena widths from centre): min -1.6099, median -0.1016, max 0.7682 | min -0.9328, median 0.0185, max 0.5415
Identities: ['a', 'b', 'c']
delay_object_blanks true: 0 ; nonzero object velocities: 0 ; trials without exactly one target: 0
Distinct triangle object locations on on-triangle trials: [np.int64(0), np.int64(1), np.int64(2)]
Units: spikes unsorted in 0 units; spikes outside the stable intervals: 0
Unit mean delay-phase rate (Hz): min 0.0093, median 5.9615, max 131.7202 ; units with no completed stable trial: 0 closed-loop minus post-hoc response time (s): min -0.2832, median 0.1341, max 1.5242
Clock check, share of units location-modulated (Kruskal-Wallis p < 0.01):
monkey session task pre early early_shifted delay delay_shifted
E 2022-08-19 triangle 0.003 0.035 0.009 0.076 0.013
E 2022-08-20 triangle 0.010 0.051 0.017 0.059 0.002
E 2022-08-21 triangle 0.004 0.048 0.017 0.044 0.013
E 2022-08-22 triangle 0.004 0.095 0.009 0.078 0.017
E 2022-08-23 triangle 0.023 0.019 0.019 0.051 0.023
E 2022-08-24 triangle 0.018 0.067 0.014 0.057 0.007
E 2022-08-25 triangle 0.006 0.085 0.006 0.037 0.018
E 2022-08-26 triangle 0.021 0.086 0.014 0.071 0.014
E 2022-08-31 triangle 0.005 0.176 0.015 0.190 0.024
E 2022-09-01 triangle 0.009 0.049 0.018 0.098 0.009
E 2022-09-02 triangle 0.007 0.013 0.000 0.013 0.013
E 2022-09-03 triangle 0.010 0.064 0.025 0.153 0.000
E 2022-09-04 triangle 0.003 0.107 0.014 0.103 0.014
E 2022-09-05 triangle 0.000 0.144 0.033 0.105 0.000
P 2022-05-26 triangle 0.010 0.277 0.010 0.317 0.000
P 2022-05-27 triangle 0.027 0.145 0.018 0.245 0.018
P 2022-05-28 triangle 0.000 0.244 0.015 0.296 0.015
P 2022-05-29 triangle 0.016 0.262 0.008 0.214 0.016
P 2022-05-30 triangle 0.011 0.239 0.005 0.402 0.027
P 2022-05-31 triangle 0.015 0.337 0.010 0.452 0.000
P 2022-06-01 triangle 0.004 0.426 0.027 0.635 0.000
P 2022-06-03 triangle 0.005 0.311 0.020 0.347 0.005
P 2022-06-04 triangle 0.008 0.229 0.008 0.196 0.000
P 2022-06-05 triangle 0.009 0.204 0.028 0.185 0.019
E 2022-09-06 ring 0.005 0.046 0.005 0.066 0.000
E 2022-09-07 ring 0.000 0.051 0.000 0.064 0.013
E 2022-09-09 ring 0.000 0.165 0.016 0.299 0.000
E 2022-09-10 ring 0.006 0.170 0.000 0.067 0.006
E 2022-09-11 ring 0.014 0.041 0.014 0.082 0.007
E 2022-09-12 ring 0.007 0.075 0.000 0.104 0.007
E 2022-09-13 ring 0.000 0.143 0.000 0.183 0.014
E 2022-09-14 ring 0.010 0.214 0.000 0.071 0.000
E 2022-09-15 ring 0.000 0.269 0.022 0.129 0.000
E 2022-09-18 ring 0.000 0.052 0.026 0.013 0.000
E 2022-09-19 ring 0.000 0.079 0.000 0.018 0.000
E 2022-09-20 ring 0.016 0.048 0.016 0.063 0.016
P 2022-06-07 ring 0.000 0.214 0.000 0.429 0.000
P 2022-06-09 ring 0.005 0.233 0.005 0.280 0.021
P 2022-06-10 ring 0.019 0.370 0.009 0.435 0.023
P 2022-06-11 ring 0.005 0.556 0.014 0.570 0.000
P 2022-06-12 ring 0.017 0.372 0.000 0.303 0.004
P 2022-06-13 ring 0.013 0.194 0.004 0.108 0.013
P 2022-06-14 ring 0.016 0.192 0.000 0.170 0.005
P 2022-06-15 ring 0.000 0.360 0.015 0.324 0.000
P 2022-06-16 ring 0.014 0.323 0.000 0.313 0.009
P 2022-06-19 ring 0.005 0.430 0.000 0.539 0.016
Sessions where early modulation is not above both controls: ['2022-08-23']
Sessions where delay modulation is not above the shifted control: []
Share of completed trials rewarded, by closed-loop distance to the target (arena widths; 0.1 = 2.4 deg, 0.2 = 4.8 deg):
completed trials share rewarded
target_err_closed_loop
[0.0, 0.05) 27581 1.0
[0.05, 0.1) 18064 1.0
[0.1, 0.15) 2902 1.0
[0.15, 0.2) 638 1.0
[0.2, 0.25) 244 0.0
[0.25, 0.3) 95 0.0
[0.3, 10.0) 12130 0.0
Electrode locations named in the spikesorting files: ['unknown'] ; the area of a unit is taken from its probe (s0 = DMFC Neuropixels, vprobe0/1 = FEF V-probes; claims file, code unit_filters.py).
Sessions sharing trial timing or sorted spike counts: []
Checks of the metadata against the data. Sessions: 24 Triangle physiology sessions (P 10, E 14) and 22 Ring sessions (P 10, E 12), as the Methods state; two more Perle sessions are in the dandiset, 2022-06-06 (spike sorting, 2-object trials only, so not a Ring session by the Methods' design) and 2022-06-08 (behaviour only); both are left out. No two sessions share their trial timing or their sorted spike counts. Objects lie at eccentricity 0.35 arena widths (8.4 degrees), with three identities and, on the triangle, three locations, as described. The stimulus and delay phases last 1.000 s (median), as described, but the cue phase (cue onset to response-phase onset) lasts 0.250 s, where the Methods say 180 ms; 47.9% of completed trials have a post-hoc response time earlier than 0.25 s after cue onset, 1.1% earlier than 0.18 s, and the earliest is 12.5 ms before cue onset, so response_time reads as saccade onset (the closed-loop time is 134 ms later at the median). Every completed trial with a closed-loop error under 0.2 arena widths (4.8 degrees) was rewarded and none at 0.2 or more, so the task scored correctness at the code's radius, not the Methods' 2.4 degrees. The electrodes tables name no location ("unknown"); areas come from the probe. Spike times are sorted and all lie inside the units' stable intervals. Clock: in every session the share of location-modulated units in the delay exceeds the share with labels shifted by one trial; early in the stimulus it exceeds both controls in every session except Elgar 2022-08-23 (0.019 against 0.023 before onset), whose delay-phase modulation is above control (0.051 against 0.023). The trial table and the spike times are consistent with a shared clock in every session; for Elgar 2022-08-23 this rests on the delay window only.
Behaviour¶
fig1b-performance-load¶
In both monkeys accuracy falls as the number of objects grows, yet stays far above the chance of saccading to a random stimulus object for two- and three-object stimuli.
Source: Results, Task, Behavior, and Hypotheses, paragraph 1; Figure 1b. Paper: "Animals learned the task and exhibited variable reaction times, with error rates increasing as the number of objects grew ... However, performance remained far above chance for two- and three-object stimuli"
Criterion (claims file): Per monkey, the fraction correct (T2) of completed trials (T1) for 1, 2 and 3 objects; chance is 1/n, as the Fig. 1b legend defines it. Reproduced when, for both monkeys, the fraction falls from 1 to 2 to 3 objects and the 2- and 3-object fractions lie above 1/2 and 1/3. The paper gives no test; "far above" is not given a number, so the margin is reported and not judged. Which trials form the figure is open: (a) Triangle trials on the triangle only, as the paper's figures/figure_1/performance/performance.ipynb plots; (b) Triangle and Ring trials pooled, as the legend says (the Ring has no 3-object trials). Each reading is run under both correctness radii of T2. For reference, that notebook printed, for reading (a) with the code's radius, P 0.953, 0.93, 0.692 and E 0.928, 0.759, 0.561 for 1, 2, 3 objects; these are reported against, not used to judge.
# fig1b-performance-load. Completed trials (T1) of physiology sessions; correctness under both
# radii (T2). Reading (a): Triangle trials on the triangle; (b): those pooled with all Ring trials.
Tp = T[T.completed & T.physiology]
sel_a = Tp[(Tp.task == "triangle") & Tp.on_triangle]
sel_b = pd.concat([sel_a, Tp[Tp.task == "ring"]])
perf = []
for reading, d in [("a", sel_a), ("b", sel_b)]:
for rad in wm.RADII:
g = d.groupby(["subject", "n_obj"])[f"correct_{rad}"].agg(["mean", "size"]).reset_index()
for _, r in g.iterrows():
perf.append(dict(reading=reading, radius=rad, monkey=MONKEY[r.subject], n_obj=r.n_obj,
frac_correct=r["mean"], chance=1 / r.n_obj, n_trials=r["size"]))
perf = pd.DataFrame(perf)
perf["margin"] = perf.frac_correct - perf.chance
print(perf.round(3).to_string(index=False))
ok = {}
for (reading, rad), g in perf.groupby(["reading", "radius"]):
res = []
for m, gm in g.groupby("monkey"):
f = gm.set_index("n_obj").frac_correct
res.append(bool(f[1] > f[2] > f[3] and f[2] > 0.5 and f[3] > 1 / 3))
ok[(reading, rad)] = all(res)
print("\nFalls 1 > 2 > 3 and above chance at 2 and 3, both monkeys:", ok)
fig, axes = plt.subplots(1, 2, figsize=(7, 3), sharey=True)
for ax, m in zip(axes, ["P", "E"]):
for (reading, rad), g in perf[perf.monkey == m].groupby(["reading", "radius"]):
ax.plot(g.n_obj, g.frac_correct, "o-", label=f"{reading}, {rad}")
ax.plot([1, 2, 3], [1, 1/2, 1/3], "k--", label="chance")
ax.set_title(f"Monkey {m}"); ax.set_xlabel("objects"); ax.set_ylim(0, 1)
axes[0].set_ylabel("fraction correct"); axes[1].legend(fontsize=7); plt.show()
print("All plotted values lie in [0, 1].")
reading radius monkey n_obj frac_correct chance n_trials margin
a methods_2.4deg E 1 0.903 1.000 1713 -0.097
a methods_2.4deg E 2 0.749 0.500 6380 0.249
a methods_2.4deg E 3 0.555 0.333 8697 0.222
a methods_2.4deg P 1 0.901 1.000 1266 -0.099
a methods_2.4deg P 2 0.887 0.500 4071 0.387
a methods_2.4deg P 3 0.661 0.333 5376 0.327
a code_4.8deg E 1 0.928 1.000 1713 -0.072
a code_4.8deg E 2 0.759 0.500 6380 0.259
a code_4.8deg E 3 0.561 0.333 8697 0.228
a code_4.8deg P 1 0.953 1.000 1266 -0.047
a code_4.8deg P 2 0.930 0.500 4071 0.430
a code_4.8deg P 3 0.692 0.333 5376 0.359
b methods_2.4deg E 1 0.901 1.000 4693 -0.099
b methods_2.4deg E 2 0.730 0.500 20089 0.230
b methods_2.4deg E 3 0.555 0.333 8697 0.222
b methods_2.4deg P 1 0.837 1.000 3832 -0.163
b methods_2.4deg P 2 0.815 0.500 14115 0.315
b methods_2.4deg P 3 0.661 0.333 5376 0.327
b code_4.8deg E 1 0.939 1.000 4693 -0.061
b code_4.8deg E 2 0.757 0.500 20089 0.257
b code_4.8deg E 3 0.561 0.333 8697 0.228
b code_4.8deg P 1 0.958 1.000 3832 -0.042
b code_4.8deg P 2 0.911 0.500 14115 0.411
b code_4.8deg P 3 0.692 0.333 5376 0.359
Falls 1 > 2 > 3 and above chance at 2 and 3, both monkeys: {('a', 'code_4.8deg'): True, ('a', 'methods_2.4deg'): True, ('b', 'code_4.8deg'): True, ('b', 'methods_2.4deg'): True}
All plotted values lie in [0, 1].
Result: reproduced. In both monkeys and under every reading and radius the fraction correct falls from 1 to 2 to 3 objects and the 2- and 3-object fractions are above 1/2 and 1/3. Reading (a) with the code's radius gives P 0.953, 0.930, 0.692 and E 0.928, 0.759, 0.561, the values the paper's notebook printed. The smallest margin above chance is 0.222 (E, 3 objects, Methods radius).
fig1c-reaction-times¶
Reaction times, from cue onset to the response saccade, vary from trial to trial, with trials faster than 180 ms counted as broken fixation.
Source: Results, Task, Behavior, and Hypotheses, paragraph 1; Figure 1c legend. Paper: "Animals learned the task and exhibited variable reaction times"
Criterion (claims file): Per monkey, reaction times (T3) of completed trials, correct and incorrect apart, as in figures/figure_1/reaction_times/reaction_times.ipynb, which bins them from 100 to 600 ms and shows 150 to 500 ms. Trials with a saccade earlier than 180 ms after cue onset are left out as broken fixation (Fig. 1c legend). Which trials form the figure is open, as for fig1b-performance-load: (a) Triangle trials on the triangle, as the notebook plots; (b) Triangle and Ring trials pooled, as the legend says; both are run. The paper fixes no measure of variability. Reproduced when the distribution of each monkey is not degenerate (interquartile range above one 10 ms bin) and its median lies inside the plotted 150 to 500 ms range; these cut-offs are the agent's reading of "variable" and are recorded as such. The histograms, quartiles and the share of trials under 180 ms before the cut are reported.
# fig1c-reaction-times. Reaction time (T3) of completed trials, correct (code radius) and
# incorrect apart; trials under 180 ms are left out (Fig. 1c legend).
rt_rows = []
fig, axes = plt.subplots(2, 2, figsize=(8, 5), sharex=True)
for i, (reading, d) in enumerate([("a", sel_a), ("b", sel_b)]):
for j, m in enumerate(["Perle", "Elgar"]):
dm = d[d.subject == m]
under = float((dm.reaction_time < 0.18).mean())
dk = dm[dm.reaction_time >= 0.18]
for lab, dd in [("all", dk), ("correct", dk[dk["correct_code_4.8deg"]]),
("incorrect", dk[~dk["correct_code_4.8deg"]])]:
q = np.percentile(dd.reaction_time, [25, 50, 75])
rt_rows.append(dict(reading=reading, monkey=MONKEY[m], trials=lab, n=len(dd),
q25_ms=1000 * q[0], median_ms=1000 * q[1], q75_ms=1000 * q[2],
iqr_ms=1000 * (q[2] - q[0]), share_under_180ms_before_cut=under,
outside_150_500ms=float(((dd.reaction_time < 0.15) | (dd.reaction_time > 0.5)).mean())))
ax = axes[i, j]
bins = np.arange(0.1, 0.61, 0.01)
ax.hist(dk[dk["correct_code_4.8deg"]].reaction_time, bins, alpha=0.6, label="correct")
ax.hist(dk[~dk["correct_code_4.8deg"]].reaction_time, bins, alpha=0.6, label="incorrect")
ax.set_xlim(0.15, 0.5); ax.set_title(f"reading {reading}, Monkey {MONKEY[m]}", fontsize=9)
axes[0, 0].legend(); axes[1, 0].set_xlabel("reaction time (s)"); plt.tight_layout(); plt.show()
rt = pd.DataFrame(rt_rows)
print(rt.round(3).to_string(index=False))
print("(Histograms show 150 to 500 ms; the share of trials outside that range is in column outside_150_500ms.)")
a = rt[rt.trials == "all"]
print("\nIQR above 10 ms and median within 150-500 ms, per reading and monkey:")
print(a.assign(holds=(a.iqr_ms > 10) & a.median_ms.between(150, 500))[["reading", "monkey", "iqr_ms", "median_ms", "holds"]].round(1).to_string(index=False))
reading monkey trials n q25_ms median_ms q75_ms iqr_ms share_under_180ms_before_cut outside_150_500ms
a P all 10532 211.192 229.981 263.232 52.041 0.017 0.004
a P correct 8611 212.931 232.195 266.848 53.917 0.017 0.004
a P incorrect 1921 202.383 220.697 246.986 44.603 0.017 0.004
a E all 16726 238.651 268.522 308.429 69.778 0.004 0.002
a E correct 11277 242.095 272.420 313.883 71.788 0.004 0.002
a E incorrect 5449 231.969 260.669 295.528 63.559 0.004 0.002
b P all 22786 208.883 227.356 259.476 50.593 0.023 0.003
b P correct 19882 210.368 228.837 261.305 50.937 0.023 0.003
b P incorrect 2904 200.484 218.281 246.281 45.797 0.023 0.005
b E all 33377 242.411 273.421 315.508 73.097 0.003 0.003
b E correct 24434 246.736 278.258 322.101 75.365 0.003 0.003
b E incorrect 8943 232.595 261.110 297.199 64.605 0.003 0.003
(Histograms show 150 to 500 ms; the share of trials outside that range is in column outside_150_500ms.)
IQR above 10 ms and median within 150-500 ms, per reading and monkey:
reading monkey iqr_ms median_ms holds
a P 52.0 230.0 True
a E 69.8 268.5 True
b P 50.6 227.4 True
b E 73.1 273.4 True
Result: reproduced. Under both readings each monkey's distribution is broad (interquartile range 50.6 to 73.1 ms) with its median inside 150 to 500 ms (227 to 273 ms). Before the 180 ms cut, 1.7% (P) and 0.4% (E) of reading (a) trials fall under 180 ms. The cut-offs are the claims file's reading of "variable".
ext1-misbinding-errors¶
A common error is a saccade to the location of a non-target stimulus object, reporting the cued identity at the wrong location.
Source: Discussion, Position-identity binding, paragraph 1; Extended Data Fig. 1. Paper: "A common error in our task occurs when the monkey reports the identity of an object at the wrong location ( Extended Data Fig. 1 ), consistent with misbinding in working memory"
Criterion (claims file): Completed Triangle trials on the triangle with 2 or 3 objects, per monkey. Each response is assigned to the triangle vertex within the correctness radius (T2), as in load_triangle_data.py (response_location), and counted by target location and response location as in Extended Data Fig. 1. Among incorrect trials, the share whose response lies at a non-target object's vertex. The paper gives no number for "common"; two readings are run, (a) these responses are the largest class of errors (against responses at an empty vertex or at no vertex), and (b) they are more than half of the errors. Reproduced when both readings hold for both monkeys and both set sizes, not reproduced when neither does; when the readings give different outcomes the claim is methods_underspecified, with the shares reported.
# ext1-misbinding-errors. Completed Triangle trials on the triangle with 2 or 3 objects; each
# response assigned to the vertex within the radius (load_triangle_data.py, response_location).
d = sel_a[sel_a.n_obj >= 2].copy()
ext1 = []
for rad, radius in wm.RADII.items():
d["resp_loc"] = wm.response_location(d, radius)
locs = d[["loc0", "loc1", "loc2"]].values
tloc = locs[np.arange(len(d)), d.target.values]
occupied = (locs == d.resp_loc.values[:, None]).any(1) & (d.resp_loc.values >= 0)
cls = np.where(d[f"correct_{rad}"], "correct",
np.where(d.resp_loc.values < 0, "no vertex",
np.where(d.resp_loc.values == tloc, "target vertex, outside radius",
np.where(occupied, "non-target object", "empty vertex"))))
d["cls"] = cls
for (m, n), g in d.groupby(["subject", "n_obj"]):
e = g[g.cls != "correct"].cls.value_counts()
share = e.get("non-target object", 0) / e.sum()
largest = e.idxmax() == "non-target object"
ext1.append(dict(radius=rad, monkey=MONKEY[m], n_obj=n, n_trials=len(g), n_errors=int(e.sum()),
non_target=int(e.get("non-target object", 0)), empty_vertex=int(e.get("empty vertex", 0)),
no_vertex=int(e.get("no vertex", 0)), target_outside=int(e.get("target vertex, outside radius", 0)),
share_non_target=share, a_largest=bool(largest), b_over_half=bool(share > 0.5)))
ext1 = pd.DataFrame(ext1)
print(ext1.round(3).to_string(index=False))
for rad in wm.RADII:
e = ext1[ext1.radius == rad]
print(f"{rad}: reading (a) holds in all four: {e.a_largest.all()}; reading (b): {e.b_over_half.all()}")
# Counts by target and response location, as Extended Data Fig. 1 plots them (code radius)
d["resp_loc"] = wm.response_location(d, wm.RADII["code_4.8deg"])
d["target_loc"] = d[["loc0", "loc1", "loc2"]].values[np.arange(len(d)), d.target.values]
print("\nTrials by target location (rows) and response location (columns, -1 = none), code radius:")
for (m, n), g in d.groupby(["subject", "n_obj"]):
print(f"Monkey {MONKEY[m]}, {n} objects"); print(pd.crosstab(g.target_loc, g.resp_loc).to_string())
radius monkey n_obj n_trials n_errors non_target empty_vertex no_vertex target_outside share_non_target a_largest b_over_half methods_2.4deg E 2 6380 1604 1388 73 143 0 0.865 True True methods_2.4deg E 3 8697 3866 3740 0 126 0 0.967 True True methods_2.4deg P 2 4071 461 207 5 249 0 0.449 False False methods_2.4deg P 3 5376 1824 1495 0 329 0 0.820 True True code_4.8deg E 2 6380 1537 1418 83 36 0 0.923 True True code_4.8deg E 3 8697 3818 3797 0 21 0 0.994 True True code_4.8deg P 2 4071 283 239 8 36 0 0.845 True True code_4.8deg P 3 5376 1655 1597 0 58 0 0.965 True True methods_2.4deg: reading (a) holds in all four: False; reading (b): False code_4.8deg: reading (a) holds in all four: True; reading (b): True Trials by target location (rows) and response location (columns, -1 = none), code radius: Monkey E, 2 objects resp_loc -1 0 1 2 target_loc 0 6 1840 250 107 1 17 200 1491 272 2 13 270 402 1512 Monkey E, 3 objects resp_loc -1 0 1 2 target_loc 0 11 2139 472 338 1 5 723 1204 953 2 5 557 754 1536 Monkey P, 2 objects resp_loc -1 0 1 2 target_loc 0 8 1360 8 13 1 8 33 1187 110 2 20 10 73 1241 Monkey P, 3 objects resp_loc -1 0 1 2 target_loc 0 8 1555 115 113 1 22 184 879 731 2 28 30 424 1287
Result: reproduced, resting on an assumption. With the code's radius (0.2 arena widths), responses at a non-target object's vertex are the largest class of errors and more than half of them in both monkeys and set sizes (84.5% to 99.4%). With the Methods' radius (2.4 degrees) both readings fail for Monkey P with 2 objects (207 such errors against 249 at no vertex, 44.9%) and hold in the other three. The two radii give different outcomes; the code's radius is adopted because it is the one the task rewarded (data check: every completed trial under 0.2 rewarded, none above) and the one the paper's code uses to assign responses to vertices for Extended Data Fig. 1.
fig5ab-ring-behavior¶
Both monkeys perform well on the Ring dataset, with reaction-time variability like that on the Triangle dataset.
Source: Results, Replication of the Results with the Ring Dataset, paragraph 1; Figure 5a and 5b; Extended Data Fig. 9a and 9b (the Results cite "Extended Data Fig. 8a-b", but the per-monkey Ring figure is the ninth Extended Data legend). Paper: "Both monkeys performed well on the Ring dataset and exhibited reaction-time variability similar to the Triangle dataset"
Criterion (claims file): Completed Ring trials (T1), per monkey and pooled, fraction correct (T2) for 1 and 2 objects and the reaction-time distribution (T3), as figures/figure_5/performance_reaction_time/. The paper gives no number. Reproduced when the 2-object fraction correct is above chance (1/2) for both monkeys and each monkey's Ring reaction-time interquartile range is within a factor of 1.5 of its Triangle one; the factor is the agent's reading of "similar" and the ranges are reported. For reference, the paper's notebook printed P 0.958 and 0.903, E 0.945 and 0.756 for 1 and 2 objects with the code's radius.
# fig5ab-ring-behavior. Completed Ring trials (T1), all conditions (Methods); the code's
# exclusion of the left/right trials shown as a variant.
R5 = Tp[Tp.task == "ring"].copy()
lr = (R5.n_obj == 2) & (np.isclose(R5.theta0, 0, atol=1e-3) | np.isclose(R5.theta0, 2*np.pi, atol=1e-3)) & np.isclose(R5.theta1, np.pi, atol=1e-3)
lr |= (R5.n_obj == 2) & np.isclose(R5.theta0, np.pi, atol=1e-3) & (np.isclose(R5.theta1, 0, atol=1e-3) | np.isclose(R5.theta1, 2*np.pi, atol=1e-3))
R5["left_right"] = lr
print("Left/right 2-object trials among completed Ring trials:", int(lr.sum()), "of", len(R5), f"({lr.mean():.3f})")
rows = []
for variant, d in [("all conditions", R5), ("without left/right", R5[~R5.left_right])]:
for rad in wm.RADII:
for m, dm in list(d.groupby("subject")) + [("pooled", d)]:
for n in [1, 2]:
g = dm[dm.n_obj == n]
rows.append(dict(variant=variant, radius=rad, monkey=MONKEY.get(m, m), n_obj=n,
n_trials=len(g), frac_correct=g[f"correct_{rad}"].mean()))
p5 = pd.DataFrame(rows)
print(p5.round(3).to_string(index=False))
rt5 = []
for m in ["Perle", "Elgar"]:
q_r = np.percentile(R5[(R5.subject == m) & (R5.reaction_time >= 0.18)].reaction_time, [25, 50, 75])
q_t = np.percentile(sel_a[(sel_a.subject == m) & (sel_a.reaction_time >= 0.18)].reaction_time, [25, 50, 75])
rt5.append(dict(monkey=MONKEY[m], ring_median_ms=1000*q_r[1], ring_iqr_ms=1000*(q_r[2]-q_r[0]),
triangle_median_ms=1000*q_t[1], triangle_iqr_ms=1000*(q_t[2]-q_t[0]),
ratio=(q_r[2]-q_r[0])/(q_t[2]-q_t[0])))
rt5 = pd.DataFrame(rt5)
print("\nReaction times (trials >= 180 ms), Ring against Triangle (on-triangle trials):")
print(rt5.round(3).to_string(index=False))
for rad in wm.RADII:
q = p5[(p5.variant == "all conditions") & (p5.radius == rad) & (p5.n_obj == 2) & p5.monkey.isin(["P", "E"])]
print(f"{rad}: 2-object above 1/2 for both monkeys: {bool((q.frac_correct > 0.5).all())}")
print("IQR ratio within a factor of 1.5 for both monkeys:", bool(rt5.ratio.between(1/1.5, 1.5).all()))
fig, ax = plt.subplots(figsize=(4, 3))
for m in ["Perle", "Elgar"]:
ax.hist(R5[(R5.subject == m) & (R5.reaction_time >= 0.18)].reaction_time, np.arange(0.1, 0.61, 0.01), alpha=0.5, label=MONKEY[m])
ax.set_xlim(0.15, 0.5); ax.set_xlabel("Ring reaction time (s)"); ax.legend(); plt.show()
Left/right 2-object trials among completed Ring trials: 1525 of 29299 (0.052)
variant radius monkey n_obj n_trials frac_correct
all conditions methods_2.4deg E 1 2980 0.900
all conditions methods_2.4deg E 2 13709 0.721
all conditions methods_2.4deg P 1 2566 0.806
all conditions methods_2.4deg P 2 10044 0.786
all conditions methods_2.4deg pooled 1 5546 0.856
all conditions methods_2.4deg pooled 2 23753 0.748
all conditions code_4.8deg E 1 2980 0.945
all conditions code_4.8deg E 2 13709 0.756
all conditions code_4.8deg P 1 2566 0.960
all conditions code_4.8deg P 2 10044 0.903
all conditions code_4.8deg pooled 1 5546 0.952
all conditions code_4.8deg pooled 2 23753 0.818
without left/right methods_2.4deg E 1 2980 0.900
without left/right methods_2.4deg E 2 12765 0.712
without left/right methods_2.4deg P 1 2566 0.806
without left/right methods_2.4deg P 2 9463 0.776
without left/right methods_2.4deg pooled 1 5546 0.856
without left/right methods_2.4deg pooled 2 22228 0.740
without left/right code_4.8deg E 1 2980 0.945
without left/right code_4.8deg E 2 12765 0.749
without left/right code_4.8deg P 1 2566 0.960
without left/right code_4.8deg P 2 9463 0.899
without left/right code_4.8deg pooled 1 5546 0.952
without left/right code_4.8deg pooled 2 22228 0.813
Reaction times (trials >= 180 ms), Ring against Triangle (on-triangle trials):
monkey ring_median_ms ring_iqr_ms triangle_median_ms triangle_iqr_ms ratio
P 225.041 48.649 229.981 52.041 0.935
E 278.771 76.601 268.522 69.778 1.098
methods_2.4deg: 2-object above 1/2 for both monkeys: True
code_4.8deg: 2-object above 1/2 for both monkeys: True
IQR ratio within a factor of 1.5 for both monkeys: True
Result: reproduced. The 2-object fraction correct on Ring trials is above 1/2 for both monkeys under both radii (code radius P 0.903, E 0.756; Methods radius P 0.786, E 0.721). The Ring interquartile range of reaction time is 0.935 (P) and 1.098 (E) times the Triangle one, inside the factor of 1.5. With the code's radius the 1-object values are P 0.960 and E 0.945 against the paper's notebook's 0.958 and 0.945, and the 2-object values match its 0.903 and 0.756.
Units¶
triangle-unit-counts¶
The Triangle dataset holds 4,688 neurons, 2,937 in DMFC and 1,751 in FEF, from 24 sessions, 10 from Monkey P and 14 from Monkey E.
Source: Results, The Triangle Dataset, paragraph 1; Methods, Neurophysiology, last paragraph. Paper: "We recorded 4,688 well-isolated neurons (2,937 from DMFC and 1,751 from FEF) in the Triangle dataset across 24 sessions"
Criterion (claims file): Sessions are the counts stated in Methods, Neurophysiology, so no tolerance. Unit counts take the default 10%. Which units are counted is open: (a) good and mua units with delay-phase rate above 1 Hz on Triangle trials, which is how figures/figure_2/selectivity/selectivity.ipynb printed 4688, 2937 and 1751; (b) as (a) and also stable for at least 500 trials (T5, Methods); (c) good units only, reading "well-isolated" literally. Reproduced when all five values match under one reading. DMFC is probe s0, FEF the V-probes.
ring-unit-counts¶
The Ring dataset holds 3,051 neurons, 1,712 in DMFC and 1,339 in FEF, from 22 sessions, 10 from Monkey P and 12 from Monkey E.
Source: Results, Replication of the Results with the Ring Dataset, paragraph 1; Methods, Neurophysiology, last paragraph. Paper: "We recorded 3,051 well-isolated neurons (1,712 from DMFC and 1,339 from FEF) across 22 neurophysiology sessions"
Criterion (claims file): As triangle-unit-counts, with the same three readings, on Ring trials. The paper's Ring notebooks list 20 sessions (10 per monkey); the session dates found are reported against both counts.
# triangle-unit-counts and ring-unit-counts. A unit's delay-phase mean rate is the mean over its
# completed stable trials of the session (T4, T5; run_cache_selectivity.py mean_firing_rate).
U = []
for x in S:
if x["units"] is None:
continue
u = x["units"].copy()
u["subject"], u["session"], u["task"] = x["subject"], x["session"], x["task"]
u["delay_rate"] = np.nanmean(x["rates"]["delay"], 1)
u["stim_rate"] = np.nanmean(x["rates"]["stimulus"], 1)
u["n_completed_stable"] = np.isfinite(x["rates"]["delay"]).sum(1)
u["dataset"] = "triangle" if x in tri else ("ring" if x in ring else "other")
U.append(u)
U = pd.concat(U, ignore_index=True)
readings = {
"a: good+mua, rate > 1 Hz": U.delay_rate > wm.MIN_RATE_HZ,
"b: (a) and stable >= 500 trials": (U.delay_rate > wm.MIN_RATE_HZ) & (U.n_obs_trials >= wm.MIN_STABLE_TRIALS),
"c: good only, rate > 1 Hz": (U.delay_rate > wm.MIN_RATE_HZ) & (U.quality == "good"),
}
cnt = []
for ds in ["triangle", "ring", "other"]:
for name, m in readings.items():
d = U[m & (U.dataset == ds)]
cnt.append(dict(dataset=ds, reading=name, total=len(d), DMFC=(d.area == "DMFC").sum(), FEF=(d.area == "FEF").sum(),
P=(d.subject == "Perle").sum(), E=(d.subject == "Elgar").sum(),
sessions_P=d[d.subject == "Perle"].session.nunique(), sessions_E=d[d.subject == "Elgar"].session.nunique()))
d = U[U.dataset == ds]
cnt.append(dict(dataset=ds, reading="all sorted units", total=len(d), DMFC=(d.area == "DMFC").sum(), FEF=(d.area == "FEF").sum(),
P=(d.subject == "Perle").sum(), E=(d.subject == "Elgar").sum(),
sessions_P=d[d.subject == "Perle"].session.nunique(), sessions_E=d[d.subject == "Elgar"].session.nunique()))
cnt = pd.DataFrame(cnt)
print(cnt.to_string(index=False))
exp = {"triangle": (4688, 2937, 1751, 10, 14), "ring": (3051, 1712, 1339, 10, 12)}
for ds, (tot, dm, fe, sp, se) in exp.items():
for _, r in cnt[(cnt.dataset == ds) & cnt.reading.str.match("[abc]:")].iterrows():
ok = (abs(r.total - tot) <= round(0.1 * tot) and abs(r.DMFC - dm) <= round(0.1 * dm)
and abs(r.FEF - fe) <= round(0.1 * fe) and r.sessions_P == sp and r.sessions_E == se)
print(f"{ds}, {r.reading}: all five values match: {ok}")
print("Ring session dates:", {m: sorted(U[(U.dataset == 'ring') & (U.subject == m)].session.unique()) for m in MONKEY})
d6 = U[(U.dataset == "other") & (U.delay_rate > 1)]
print("Ring with Perle 2022-06-06 added, reading a:", int(cnt.query("dataset=='ring'").iloc[0].total) + len(d6), "units, P sessions 11")
code_ring = ["2022-06-07","2022-06-09","2022-06-10","2022-06-11","2022-06-12","2022-06-13","2022-06-14","2022-06-15","2022-06-16","2022-06-19",
"2022-09-06","2022-09-07","2022-09-09","2022-09-10","2022-09-11","2022-09-12","2022-09-13","2022-09-14","2022-09-15","2022-09-20"]
d20 = U[(U.dataset == "ring") & U.session.isin(code_ring) & (U.delay_rate > 1)]
print("Ring, the 20 sessions of the code's model comparison, reading a:", len(d20), "units; DMFC", (d20.area == "DMFC").sum(), "FEF", (d20.area == "FEF").sum())
dataset reading total DMFC FEF P E sessions_P sessions_E
triangle a: good+mua, rate > 1 Hz 4688 2937 1751 1536 3152 10 14
triangle b: (a) and stable >= 500 trials 4473 2857 1616 1374 3099 10 14
triangle c: good only, rate > 1 Hz 1563 1110 453 614 949 10 14
triangle all sorted units 4899 3030 1869 1638 3261 10 14
ring a: good+mua, rate > 1 Hz 3052 1712 1340 1729 1323 10 12
ring b: (a) and stable >= 500 trials 2936 1655 1281 1648 1288 10 12
ring c: good only, rate > 1 Hz 1123 841 282 762 361 10 12
ring all sorted units 3186 1762 1424 1821 1365 10 12
other a: good+mua, rate > 1 Hz 69 6 63 69 0 1 0
other b: (a) and stable >= 500 trials 68 6 62 68 0 1 0
other c: good only, rate > 1 Hz 18 5 13 18 0 1 0
other all sorted units 77 6 71 77 0 1 0
triangle, a: good+mua, rate > 1 Hz: all five values match: True
triangle, b: (a) and stable >= 500 trials: all five values match: True
triangle, c: good only, rate > 1 Hz: all five values match: False
ring, a: good+mua, rate > 1 Hz: all five values match: True
ring, b: (a) and stable >= 500 trials: all five values match: True
ring, c: good only, rate > 1 Hz: all five values match: False
Ring session dates: {'Perle': ['2022-06-07', '2022-06-09', '2022-06-10', '2022-06-11', '2022-06-12', '2022-06-13', '2022-06-14', '2022-06-15', '2022-06-16', '2022-06-19'], 'Elgar': ['2022-09-06', '2022-09-07', '2022-09-09', '2022-09-10', '2022-09-11', '2022-09-12', '2022-09-13', '2022-09-14', '2022-09-15', '2022-09-18', '2022-09-19', '2022-09-20']}
Ring with Perle 2022-06-06 added, reading a: 3121 units, P sessions 11
Ring, the 20 sessions of the code's model comparison, reading a: 2862 units; DMFC 1537 FEF 1325
Result: reproduced. Reading (a), good and multi-units above 1 Hz in the delay phase, gives 4,688 units, 2,937 in DMFC and 1,751 in FEF, from 10 Perle and 14 Elgar sessions, the paper's numbers exactly. Reading (b), adding the Methods' 500-trial stability rule, gives 4,473, 2,857 and 1,616, also within tolerance. Reading (c), good units only, gives 1,563 and does not match: only 1,563 of the 4,688 "well-isolated neurons" are labelled single units.
Result: reproduced. Reading (a) gives 3,052 Ring units, 1,712 in DMFC and 1,340 in FEF, from 10 Perle and 12 Elgar sessions, against 3,051, 1,712 and 1,339; reading (b) gives 2,936, 1,655 and 1,281, also within tolerance; reading (c) gives 1,123 and does not match. The 22 sessions include Elgar 2022-09-18 and 2022-09-19, which the code's Ring model comparison leaves out (20 sessions, 2,862 units).
ext2-stimulus-delay-activity¶
Neurons are typically active in both the stimulus and the delay phase, and a unit's spatial preference is similar across the two phases.
Source: Results, The Triangle Dataset, paragraph 1; Extended Data Fig. 2a and 2b. Paper: "We found that neurons were typically active in both stimulus and delay phases ( Extended Data Fig. 2 )"
Criterion (claims file): Triangle units (T5), per monkey. Panel a, the stimulus-phase (0 to 1000 ms after stimulus onset) against the delay-phase mean rate (T4) per unit; panel b, per unit and triangle location, the 1-object mean rate at that location divided by the unit's mean over locations, stimulus phase against delay phase, for DMFC and FEF apart, as the legend defines it. The paper gives no test. Readings of "typically active in both": (a) the Pearson and (b) the Spearman correlation across units of panel a's rates is positive at p < 0.05, for both monkeys; (c) more than half of the units, all above 1 Hz in the delay phase by T5, are also above 1 Hz in the stimulus phase (the Methods' activity threshold applied to the other phase). Panel b is judged by readings (a) and (b) per monkey and area. Reproduced when all three readings hold for panel a in both monkeys and (a) and (b) hold for at least three of the four panel b plots; not reproduced when none hold; methods_underspecified when they differ. The second half of the statement is read from the title of the panel b legend, "Correlation in spatial selectivity across stimulus and delay phases"; the three-of-four rule is the agent's reading. The correlations and shares are reported.
# ext2-stimulus-delay-activity. Triangle units above 1 Hz (T5). Panel a: stimulus against delay
# mean rate per unit; panel b: per unit and vertex, the 1-object rate at the vertex divided by
# the unit's mean over the three vertices, stimulus against delay (Extended Data Fig. 2 legend).
UT = U[(U.dataset == "triangle") & (U.delay_rate > wm.MIN_RATE_HZ)]
ra = []
for m, g in UT.groupby("subject"):
pr, sp = st.pearsonr(g.stim_rate, g.delay_rate), st.spearmanr(g.stim_rate, g.delay_rate)
ra.append(dict(monkey=MONKEY[m], n_units=len(g), pearson_r=pr[0], pearson_p=pr[1], spearman_rho=sp[0], spearman_p=sp[1],
share_stim_above_1Hz=float((g.stim_rate > 1).mean())))
ra = pd.DataFrame(ra)
print("Panel a:"); print(ra.to_string(index=False))
a_ok = bool(((ra.pearson_r > 0) & (ra.pearson_p < 0.05)).all() and ((ra.spearman_rho > 0) & (ra.spearman_p < 0.05)).all()
and (ra.share_stim_above_1Hz > 0.5).all())
print("All three readings hold for panel a in both monkeys:", a_ok)
pb = []
for x in tri:
tr = x["trials"]
m = (tr.n_obj == 1).values & tr.on_triangle.values & tr.completed.values
keep = np.nanmean(x["rates"]["delay"], 1) > wm.MIN_RATE_HZ
for ph in ["stimulus", "delay"]:
R = x["rates"][ph][keep][:, m]
loc = tr.loc0.values[m]
mu = np.stack([np.nanmean(R[:, loc == k], 1) for k in range(3)], 1)
x[f"norm_{ph}"] = mu / mu.mean(1, keepdims=True)
areas = x["units"].area.values[keep]
for k in range(3):
for i in range(keep.sum()):
pb.append(dict(monkey=MONKEY[x["subject"]], area=areas[i], loc=k,
stim=x["norm_stimulus"][i, k], delay=x["norm_delay"][i, k]))
pb = pd.DataFrame(pb).replace([np.inf, -np.inf], np.nan).dropna()
rb = []
for (m, ar), g in pb.groupby(["monkey", "area"]):
pr, sp = st.pearsonr(g.stim, g.delay), st.spearmanr(g.stim, g.delay)
rb.append(dict(monkey=m, area=ar, n_points=len(g), pearson_r=pr[0], pearson_p=pr[1], spearman_rho=sp[0], spearman_p=sp[1]))
rb = pd.DataFrame(rb)
print("\nPanel b:"); print(rb.to_string(index=False))
nb_ok = int((((rb.pearson_r > 0) & (rb.pearson_p < 0.05)) & ((rb.spearman_rho > 0) & (rb.spearman_p < 0.05))).sum())
print(f"Panel b plots where readings (a) and (b) both hold: {nb_ok} of 4")
print("Points dropped from panel b (a unit with zero mean over vertices in a phase):",
int(sum(3 * (np.nanmean(x['rates']['delay'], 1) > 1).sum() for x in tri) - len(pb)))
fig, axes = plt.subplots(1, 2, figsize=(8, 3.5))
for ax, m in zip(axes, ["Perle", "Elgar"]):
g = UT[UT.subject == m]
for ar, col in [("DMFC", "tab:orange"), ("FEF", "tab:blue")]:
h = g[g.area == ar]; ax.scatter(h.stim_rate, h.delay_rate, s=3, c=col, label=ar)
ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlabel("stimulus rate (Hz)"); ax.set_ylabel("delay rate (Hz)")
ax.set_title(f"Monkey {MONKEY[m]}"); ax.legend()
plt.show()
print("Ranges: stimulus", UT.stim_rate.min().round(3), "to", UT.stim_rate.max().round(1), "Hz; delay",
UT.delay_rate.min().round(3), "to", UT.delay_rate.max().round(1), "Hz; units with a zero stimulus rate (not on log axes):", int((UT.stim_rate <= 0).sum()))
Panel a:
monkey n_units pearson_r pearson_p spearman_rho spearman_p share_stim_above_1Hz
E 3152 0.971079 0.0 0.963348 0.0 0.991434
P 1536 0.898321 0.0 0.893823 0.0 0.974609
All three readings hold for panel a in both monkeys: True
Panel b:
monkey area n_points pearson_r pearson_p spearman_rho spearman_p
E DMFC 6636 0.397997 9.203846e-251 0.371406 3.934155e-216
E FEF 2820 0.439628 1.291489e-133 0.384681 3.934749e-100
P DMFC 2175 0.599897 9.135321e-213 0.544998 1.579157e-168
P FEF 2433 0.487605 1.599058e-145 0.391610 5.421155e-90
Panel b plots where readings (a) and (b) both hold: 4 of 4
Points dropped from panel b (a unit with zero mean over vertices in a phase): 0
Ranges: stimulus 0.092 to 133.9 Hz; delay 1.0 to 115.3 Hz; units with a zero stimulus rate (not on log axes): 0
Result: reproduced. Stimulus- and delay-phase rates correlate across units (Pearson r 0.898 for P and 0.971 for E, Spearman 0.894 and 0.963, p below 1e-300), and 97.5% (P) and 99.1% (E) of units above 1 Hz in the delay are also above 1 Hz in the stimulus phase. The normalised location preferences of the two phases correlate positively in all four monkey and area plots (Pearson 0.398 to 0.600, Spearman 0.371 to 0.545, all p below 1e-89).
Selectivity¶
AUCs of every unit above 1 Hz, used by fig2b-selectivity-range and fig5cd-ring-tuning, and the count of units entering the models (T7) used to estimate the model fits.
# AUC (T6) of every unit above 1 Hz (T5) on 1-object trials: Triangle trials on the triangle
# labelled by vertex, Ring trials by nearest vertex (three 120-degree sectors); identity by
# fruit. Both p-value methods of T6. Own generator.
rng_auc = np.random.default_rng(620)
ID_NUM = {"a": 0, "b": 1, "c": 2}
auc_rows = []
for x in tri + ring:
tr = x["trials"]; R = x["rates"]["delay"]
m = (tr.n_obj == 1).values & tr.completed.values
if x["task"] == "triangle":
m &= tr.on_triangle.values
pos = tr.loc0.values
else:
pos = wm.ring_sector(tr.theta0.values)
idn = tr.id0.map(ID_NUM).fillna(-1).astype(int).values
rate = np.nanmean(R, 1)
for k in np.where(rate > wm.MIN_RATE_HZ)[0]:
r = R[k]; ok = m & np.isfinite(r)
row = dict(subject=x["subject"], session=x["session"], dataset=x["task"], unit_index=k,
area=x["units"].area.iloc[k], quality=x["units"].quality.iloc[k],
n_obs_trials=x["units"].n_obs_trials.iloc[k], delay_rate=rate[k], n_1obj=int(ok.sum()))
for meth in ["methods", "code"]:
row[f"pos_auc_{meth}"], row[f"pos_p_{meth}"] = wm.unit_auc(r[ok], pos[ok], rng_auc, meth)
row[f"id_auc_{meth}"], row[f"id_p_{meth}"] = wm.unit_auc(r[ok], idn[ok], rng_auc, meth)
auc_rows.append(row)
AUC = pd.DataFrame(auc_rows)
for meth in ["methods", "code"]:
AUC[f"modeled_{meth}"] = (((AUC[f"pos_auc_{meth}"] > 0.65) & (AUC[f"pos_p_{meth}"] < 0.05))
| ((AUC[f"id_auc_{meth}"] > 0.65) & (AUC[f"id_p_{meth}"] < 0.05)))
print("Units with AUCs:", AUC.groupby(["dataset", "subject"]).size().to_dict())
print("Fewest 1-object trials of a unit:", AUC.n_1obj.min())
nm = AUC.groupby(["dataset", "subject", "session"])[["modeled_methods", "modeled_code"]].sum().reset_index()
print("\nUnits entering the models per session (T7: AUC > 0.65 with p < 0.05, rate > 1 Hz):")
print(nm.to_string(index=False))
Units with AUCs: {('ring', 'Elgar'): 1323, ('ring', 'Perle'): 1729, ('triangle', 'Elgar'): 3152, ('triangle', 'Perle'): 1536}
Fewest 1-object trials of a unit: 12
Units entering the models per session (T7: AUC > 0.65 with p < 0.05, rate > 1 Hz):
dataset subject session modeled_methods modeled_code
ring Elgar 2022-09-06 13 29
ring Elgar 2022-09-07 8 15
ring Elgar 2022-09-09 27 42
ring Elgar 2022-09-10 10 20
ring Elgar 2022-09-11 8 23
ring Elgar 2022-09-12 10 15
ring Elgar 2022-09-13 22 29
ring Elgar 2022-09-14 26 29
ring Elgar 2022-09-15 22 27
ring Elgar 2022-09-18 2 6
ring Elgar 2022-09-19 5 19
ring Elgar 2022-09-20 11 13
ring Perle 2022-06-07 7 8
ring Perle 2022-06-09 89 104
ring Perle 2022-06-10 112 119
ring Perle 2022-06-11 114 127
ring Perle 2022-06-12 69 82
ring Perle 2022-06-13 53 78
ring Perle 2022-06-14 40 58
ring Perle 2022-06-15 29 39
ring Perle 2022-06-16 75 91
ring Perle 2022-06-19 79 87
triangle Elgar 2022-08-19 87 130
triangle Elgar 2022-08-20 98 172
triangle Elgar 2022-08-21 28 64
triangle Elgar 2022-08-22 78 115
triangle Elgar 2022-08-23 44 82
triangle Elgar 2022-08-24 87 141
triangle Elgar 2022-08-25 29 54
triangle Elgar 2022-08-26 27 48
triangle Elgar 2022-08-31 84 110
triangle Elgar 2022-09-01 55 89
triangle Elgar 2022-09-02 26 60
triangle Elgar 2022-09-03 54 80
triangle Elgar 2022-09-04 77 120
triangle Elgar 2022-09-05 54 77
triangle Perle 2022-05-26 59 65
triangle Perle 2022-05-27 56 70
triangle Perle 2022-05-28 60 73
triangle Perle 2022-05-29 48 64
triangle Perle 2022-05-30 109 124
triangle Perle 2022-05-31 119 132
triangle Perle 2022-06-01 170 176
triangle Perle 2022-06-03 97 111
triangle Perle 2022-06-04 112 142
triangle Perle 2022-06-05 47 59
fig2b-selectivity-range¶
Delay-phase position and identity AUCs on 1-object trials span a wide range across units.
Source: Results, The Triangle Dataset, paragraph 1; Figure 2b. Paper: "A Receiver Operating Characteristic analysis of firing rates in the delay phase revealed a wide range of selectivities for position and identity across neurons ( Fig. 2b )"
Criterion (claims file): Position and identity AUC (T6) of every Triangle unit (T5) per monkey, with significance at p < 0.05 as the legend shades it. The paper gives no number for "wide". Reproduced when, for both monkeys and both features, the AUCs run from near 0.5 (lowest decile below 0.6) to high values (highest decile above 0.7), with both significant and non-significant units present; the deciles are reported under both p-value methods of T6. These cut-offs are the agent's reading of "wide range" and are recorded as such.
# fig2b-selectivity-range: deciles of position and identity AUC per monkey, Triangle units.
def auc_range(A, label):
rows = []
for meth in ["methods", "code"]:
for m, g in list(A.groupby("subject")) + ([("pooled", A)] if label == "ring" else []):
sig = (g[f"pos_p_{meth}"] < 0.05) | (g[f"id_p_{meth}"] < 0.05)
for feat in ["pos", "id"]:
v = g[f"{feat}_auc_{meth}"].dropna()
rows.append(dict(p_method=meth, monkey=MONKEY.get(m, m), feature=feat, n_units=len(v),
min=v.min(), p10=np.percentile(v, 10), median=np.median(v), p90=np.percentile(v, 90), max=v.max(),
n_significant=int((g[f"{feat}_p_{meth}"] < 0.05).sum()),
n_shaded_significant=int(sig.sum()), n_shaded_not=int((~sig).sum())))
r = pd.DataFrame(rows)
r["holds"] = (r.p10 < 0.6) & (r.p90 > 0.7) & (r.n_shaded_significant > 0) & (r.n_shaded_not > 0)
return r
for variant, A in [("units above 1 Hz", AUC[AUC.dataset == "triangle"]),
("also stable >= 500 trials", AUC[(AUC.dataset == "triangle") & (AUC.n_obs_trials >= 500)])]:
r = auc_range(A, "triangle")
print(f"\nTriangle, {variant}:"); print(r.round(3).to_string(index=False))
print("Holds for both monkeys and both features:", r.groupby("p_method").holds.all().to_dict())
fig, axes = plt.subplots(1, 2, figsize=(8, 3.5))
for ax, m in zip(axes, ["Perle", "Elgar"]):
g = AUC[(AUC.dataset == "triangle") & (AUC.subject == m)]
sig = (g.pos_p_methods < 0.05) | (g.id_p_methods < 0.05)
ax.scatter(g.pos_auc_methods[~sig], g.id_auc_methods[~sig], s=4, c="0.8")
ax.scatter(g.pos_auc_methods[sig], g.id_auc_methods[sig], s=4, c="k")
ax.set_xlim(0.49, 1.01); ax.set_ylim(0.49, 1.01); ax.set_title(f"Monkey {MONKEY[m]}, N={len(g)}")
ax.set_xlabel("position AUC"); ax.set_ylabel("identity AUC")
plt.show()
print("All AUCs lie in [0.5, 1]:", bool(AUC.filter(like="_auc_").stack().between(0.5, 1).all()))
Triangle, units above 1 Hz:
p_method monkey feature n_units min p10 median p90 max n_significant n_shaded_significant n_shaded_not holds
methods E pos 3152 0.504 0.538 0.600 0.706 1.000 849 1156 1996 True
methods E id 3152 0.503 0.533 0.582 0.654 1.000 406 1156 1996 False
methods P pos 1536 0.508 0.565 0.682 0.888 1.000 861 955 581 True
methods P id 1536 0.504 0.539 0.600 0.700 0.969 269 955 581 False
code E pos 3152 0.505 0.547 0.615 0.722 1.000 1473 2001 1151 True
code E id 3152 0.500 0.539 0.597 0.681 1.000 1093 2001 1151 False
code P pos 1536 0.505 0.574 0.691 0.897 1.000 1029 1169 367 True
code P id 1536 0.504 0.543 0.607 0.715 0.982 507 1169 367 True
Holds for both monkeys and both features: {'code': False, 'methods': False}
Triangle, also stable >= 500 trials:
p_method monkey feature n_units min p10 median p90 max n_significant n_shaded_significant n_shaded_not holds
methods E pos 3099 0.504 0.537 0.599 0.704 1.000 842 1145 1954 True
methods E id 3099 0.503 0.533 0.582 0.651 0.829 401 1145 1954 False
methods P pos 1374 0.508 0.563 0.678 0.888 1.000 803 884 490 True
methods P id 1374 0.504 0.538 0.595 0.684 0.836 247 884 490 False
code E pos 3099 0.505 0.547 0.614 0.721 1.000 1462 1982 1117 True
code E id 3099 0.504 0.539 0.597 0.678 0.829 1082 1982 1117 False
code P pos 1374 0.505 0.572 0.687 0.897 1.000 950 1068 306 True
code P id 1374 0.504 0.540 0.601 0.697 0.982 468 1068 306 False
Holds for both monkeys and both features: {'code': False, 'methods': False}
All AUCs lie in [0.5, 1]: True
Result: partially reproduced. Position AUCs span a wide range in both monkeys under both p-value methods (lowest decile 0.538 to 0.574, highest decile 0.706 to 0.897), with significant and non-significant units present. Identity AUCs do not reach the claims file's upper cut-off: their highest decile is 0.654 (E) and 0.700 (P) with the Methods p-value, and 0.681 (E) and 0.715 (P) with the code's, so the rule holds for identity only for Monkey P under the code's p-value. Adding the 500-trial rule changes no conclusion. Identity AUCs still run from 0.50 to 0.97 or more in both monkeys. The criterion defines only the reproduced case; partially reproduced is this attempt's label for a result that holds in part.
fig5cd-ring-tuning¶
Ring neurons in DMFC and FEF show a variety of spatial receptive fields and diverse delay-phase position and identity selectivities, as in the Triangle dataset.
Source: Results, Replication of the Results with the Ring Dataset, paragraph 1; Figure 5c and 5d; Extended Data Fig. 9c. Paper: "These neurons showed a variety of spatial receptive fields in both DMFC and FEF ( Fig. 5c ; Extended Data Fig. 9c ) and diverse selectivities to position and identity in delay-phase firing rates ( Fig. 5d ), similar to the Triangle dataset"
Criterion (claims file): Position and identity AUC (T6) of Ring units, judged as fig2b-selectivity-range. Tuning curves as Methods, Ring Dataset Tuning Curves (256 points on the ring, 100 bootstraps of 100 1-object trials, Gaussian kernel of 2.4 degrees); "variety" is judged by the preferred angles of units with significant position AUC, reproduced when in each area they fall in all three 120-degree sectors.
# fig5cd-ring-tuning. Ring AUC judged as fig2b, per monkey and pooled; tuning curves as Methods,
# Ring Dataset Tuning Curves (256 points on the ring, 100 bootstraps of 100 1-object trials,
# Gaussian kernel sigma 0.1 arena widths = 2.4 deg). Preferred angle = peak of the mean curve.
r5 = auc_range(AUC[AUC.dataset == "ring"], "ring")
print("Ring AUC deciles:"); print(r5.round(3).to_string(index=False))
print("Holds for both monkeys and both features:", r5[r5.monkey != "pooled"].groupby("p_method").holds.all().to_dict())
rng_tc = np.random.default_rng(5)
ang = np.linspace(0, 2 * np.pi, 256, endpoint=False)
P = wm.RADIUS * np.stack([np.cos(ang), np.sin(ang)], 1)
pref = []
for x in ring:
tr = x["trials"]
m = (tr.n_obj == 1).values & tr.completed.values
keep = np.where(np.nanmean(x["rates"]["delay"], 1) > wm.MIN_RATE_HZ)[0]
xy = np.stack([tr.x0.values, tr.y0.values], 1)
for k in keep:
r = x["rates"]["delay"][k]; ok = m & np.isfinite(r)
rr, pp = r[ok], xy[ok]
K = np.exp(-((P[:, None, :] - pp[None]) ** 2).sum(-1) / (2 * 0.1 ** 2)) # 256 x n
idx = rng_tc.integers(0, len(rr), size=(100, 100))
curves = np.stack([(K[:, i] @ rr[i]) / K[:, i].sum(1) for i in idx])
mc = curves.mean(0)
pref.append(dict(subject=x["subject"], session=x["session"], unit_index=k, area=x["units"].area.iloc[k],
pref_angle=ang[np.argmax(mc)], depth_of_tuning=(mc.max() - mc.min()) / mc.mean()))
pref = pd.DataFrame(pref).merge(AUC[AUC.dataset == "ring"][["session", "unit_index", "pos_p_methods", "pos_p_code"]], on=["session", "unit_index"])
pref["sector"] = wm.ring_sector(pref.pref_angle.values)
for meth in ["methods", "code"]:
sg = pref[pref[f"pos_p_{meth}"] < 0.05]
tab = pd.crosstab([sg.subject, sg.area], sg.sector)
print(f"\nPreferred-angle sector (0 right, 1 top-left, 2 bottom-left) of units with significant position AUC, {meth} p:")
print(tab.to_string())
tab2 = pd.crosstab(sg.area, sg.sector)
print("Every area has units in all three sectors (pooled over monkeys):", bool((tab2 > 0).all().all() and tab2.shape[1] == 3),
"; in each monkey and area:", bool((tab > 0).all().all() and tab.shape[1] == 3))
fig, axes = plt.subplots(1, 2, figsize=(7, 3), subplot_kw=dict(projection="polar"))
for ax, ar in zip(axes, ["DMFC", "FEF"]):
sg = pref[(pref.pos_p_methods < 0.05) & (pref.area == ar)]
ax.hist(sg.pref_angle, bins=24); ax.set_title(ar)
plt.show()
Ring AUC deciles:
p_method monkey feature n_units min p10 median p90 max n_significant n_shaded_significant n_shaded_not holds
methods E pos 1323 0.502 0.528 0.573 0.658 0.934 391 507 816 False
methods E id 1323 0.502 0.522 0.558 0.611 0.901 172 507 816 False
methods P pos 1729 0.506 0.541 0.618 0.796 0.992 957 1064 665 True
methods P id 1729 0.502 0.523 0.563 0.630 0.854 283 1064 665 False
methods pooled pos 3052 0.502 0.533 0.593 0.740 0.992 1348 1571 1481 True
methods pooled id 3052 0.502 0.523 0.561 0.623 0.901 455 1571 1481 False
code E pos 1323 0.505 0.535 0.584 0.674 0.932 655 873 450 False
code E id 1323 0.504 0.527 0.569 0.638 0.914 468 873 450 False
code P pos 1729 0.504 0.545 0.628 0.805 0.992 1159 1325 404 True
code P id 1729 0.503 0.527 0.573 0.653 0.939 624 1325 404 False
code pooled pos 3052 0.504 0.539 0.605 0.751 0.992 1814 2198 854 True
code pooled id 3052 0.503 0.527 0.571 0.646 0.939 1092 2198 854 False
Holds for both monkeys and both features: {'code': False, 'methods': False}
Preferred-angle sector (0 right, 1 top-left, 2 bottom-left) of units with significant position AUC, methods p:
sector 0 1 2
subject area
Elgar DMFC 52 79 69
FEF 91 48 52
Perle DMFC 149 214 213
FEF 138 130 113
Every area has units in all three sectors (pooled over monkeys): True ; in each monkey and area: True
Preferred-angle sector (0 right, 1 top-left, 2 bottom-left) of units with significant position AUC, code p:
sector 0 1 2
subject area
Elgar DMFC 92 140 159
FEF 111 72 81
Perle DMFC 167 249 266
FEF 164 160 153
Every area has units in all three sectors (pooled over monkeys): True ; in each monkey and area: True
Result: partially reproduced. Units with significant position AUC have preferred angles in all three sectors in each area and each monkey under both p-value methods (fewest per cell 48), so the variety of receptive fields holds. The AUC part judged as fig2b fails: identity AUCs have a highest decile of 0.611 to 0.653, and Elgar's position AUCs 0.658 (Methods p) and 0.674 (code p), below 0.7; Perle's position AUCs (0.796, 0.805) and the pooled ones (0.740, 0.751) pass. The criterion defines only the reproduced case; partially reproduced is this attempt's label for a result that holds in part.
Noise correlations¶
fig2d-signal-noise-triangle¶
In the Triangle dataset, pairs of neurons with similar 1-object tuning have positive noise correlations on 2-object trials, and dissimilarly tuned pairs negative ones.
Source: Results, The Triangle Dataset, paragraph 2; Figure 2d. Paper: "Pairs of neurons with similar selectivities on 1-object trials tended to have positive noise correlations on 2-object trials ... Conversely, dissimilarly tuned pairs of neurons tended to have negative noise correlations."
Criterion (claims file): Pairs of simultaneously recorded Triangle units, per monkey, delay-phase rates (T4). Main method, Methods, Noise Correlations: signal correlation is the correlation across the 9 1-object conditions of the two units' mean rates, pairs dropped unless each condition has at least 3 trials; noise correlation is the mean over the 18 2-object conditions (identity included) with at least 3 shared trials of the across-trial correlation. Variant, the paper's figures/figure_2/noise_correlations/run_compute_triangle.py: per pair and per triple of two 1-object conditions and their 2-object composite, the signal correlation from the two 1-object conditions with trials shuffled between units, the noise correlation on the composite, units above 1 Hz, rows kept with more than 5 trials in each condition. As in plot_noise_vs_signal_correlation.ipynb, signal correlation is cut into 6 bins with edges at linspace(-0.8, 0.8, 7), ends open, and the mean noise correlation per bin is shown with its 95% confidence interval. The paper gives no test. Reproduced when, for both monkeys, the bin of highest signal correlation has a mean noise correlation above zero and the lowest bin one below zero, each with a 95% interval excluding zero; partially reproduced when only the rise from the lowest to the highest bin holds (non-overlapping intervals). Both methods are run; when they differ the claim is methods_underspecified.
# fig2d-signal-noise-triangle. Main method (Methods, Noise Correlations) and the code's variant;
# 6 bins with edges at linspace(-0.8, 0.8, 7), ends open; mean noise correlation per bin with
# a 95% interval, per monkey.
rng_nc = np.random.default_rng(2)
edges6 = np.linspace(-0.8, 0.8, 7)
NCm, NCc = [], []
for x in tri:
d, ncond = wm.triangle_noise_methods(x)
d["subject"], d["session"] = x["subject"], x["session"]; NCm.append(d)
d = wm.triangle_noise_code(x, rng_nc)
d["subject"], d["session"] = x["subject"], x["session"]; NCc.append(d)
NCm, NCc = pd.concat(NCm, ignore_index=True), pd.concat(NCc, ignore_index=True)
print("Pairs (Methods):", NCm.groupby("subject").size().to_dict(), "; rows (code):", NCc.groupby("subject").size().to_dict())
print("Pairs with a finite noise correlation (Methods):", int(np.isfinite(NCm.noise).sum()), "of", len(NCm))
def judge_bins(b):
hi, lo = b.iloc[-1], b.iloc[0]
full = hi.mean_noise > 0 and hi.ci_lo > 0 and lo.mean_noise < 0 and lo.ci_hi < 0
rise = hi.ci_lo > lo.ci_hi
return full, rise
res2d = {}
fig, axes = plt.subplots(1, 2, figsize=(8, 3), sharey=True)
for name, D, col in [("Methods", NCm, "noise"), ("code", NCc, "noise")]:
for ax, m in zip(axes, ["Perle", "Elgar"]):
b = wm.bin_means(D[D.subject == m].signal, D[D.subject == m][col], edges6)
res2d[(name, MONKEY[m])] = judge_bins(b)
print(f"\n{name}, Monkey {MONKEY[m]}:"); print(b.round(4).to_string(index=False))
ax.errorbar(np.arange(6), b.mean_noise, yerr=[b.mean_noise - b.ci_lo, b.ci_hi - b.mean_noise], marker="o", label=name)
ax.axhline(0, c="k", lw=0.5); ax.set_title(f"Monkey {MONKEY[m]}"); ax.set_xlabel("signal-correlation bin")
axes[0].set_ylabel("mean noise correlation"); axes[0].legend(); plt.show()
print("\n(highest bin > 0 and lowest < 0 with intervals excluding 0, rise with non-overlapping intervals):")
for k, v in res2d.items():
print(k, v)
# per session, Methods: lowest- and highest-bin means
ps = []
for (m, s), g in NCm.groupby(["subject", "session"]):
b = wm.bin_means(g.signal, g.noise, edges6)
ps.append(dict(monkey=MONKEY[m], session=s, pairs=len(g), lowest_bin=b.mean_noise.iloc[0], highest_bin=b.mean_noise.iloc[-1],
pearson_signal_noise=st.pearsonr(*g[["signal", "noise"]].dropna().values.T)[0]))
print(pd.DataFrame(ps).round(3).to_string(index=False))
# Code variant with intervals clustered by pair: a pair contributes up to 18 rows (one per
# triple of conditions), so rows are not independent; the cluster-robust interval counts each
# pair once (wm.bin_means_clustered).
NCc["pair"] = NCc.session + ":" + NCc.unit_0.astype(str) + "-" + NCc.unit_1.astype(str)
print("\nCode variant, rows per pair: median", int(NCc.groupby("pair").size().median()), ", max", int(NCc.groupby("pair").size().max()))
res2d_cl = {}
for m in ["Perle", "Elgar"]:
g = NCc[NCc.subject == m]
b = wm.bin_means_clustered(g.signal, g.noise, g.pair, edges6)
res2d_cl[MONKEY[m]] = judge_bins(b)
print(f"Code variant, Monkey {MONKEY[m]}, intervals clustered by pair:"); print(b.round(4).to_string(index=False))
print("(full, rise) with intervals clustered by pair, code variant:", res2d_cl)
Pairs (Methods): {'Elgar': 328273, 'Perle': 55766} ; rows (code): {'Elgar': 5087445, 'Perle': 1015982}
Pairs with a finite noise correlation (Methods): 380273 of 384039
Methods, Monkey P:
bin lo_edge hi_edge n_pairs mean_noise ci_lo ci_hi
0 -inf -0.5333 6684 -0.0041 -0.0071 -0.0011
1 -0.5333 -0.2667 8988 0.0061 0.0035 0.0088
2 -0.2667 0.0000 10945 0.0135 0.0111 0.0159
3 0.0000 0.2667 10840 0.0204 0.0180 0.0228
4 0.2667 0.5333 9321 0.0306 0.0279 0.0332
5 0.5333 inf 8502 0.0574 0.0545 0.0602
Methods, Monkey E:
bin lo_edge hi_edge n_pairs mean_noise ci_lo ci_hi
0 -inf -0.5333 25927 -0.0230 -0.0244 -0.0216
1 -0.5333 -0.2667 56490 -0.0085 -0.0095 -0.0076
2 -0.2667 0.0000 78217 0.0035 0.0027 0.0043
3 0.0000 0.2667 77179 0.0149 0.0141 0.0157
4 0.2667 0.5333 57778 0.0289 0.0279 0.0298
5 0.5333 inf 29402 0.0490 0.0475 0.0504
code, Monkey P: bin lo_edge hi_edge n_pairs mean_noise ci_lo ci_hi 0 -inf -0.5333 3389 -0.0189 -0.0292 -0.0085 1 -0.5333 -0.2667 39186 -0.0008 -0.0037 0.0022 2 -0.2667 0.0000 458200 0.0178 0.0169 0.0186 3 0.0000 0.2667 462930 0.0242 0.0234 0.0251 4 0.2667 0.5333 46041 0.0416 0.0388 0.0444 5 0.5333 inf 5265 0.0769 0.0685 0.0853
code, Monkey E: bin lo_edge hi_edge n_pairs mean_noise ci_lo ci_hi 0 -inf -0.5333 124 0.0245 -0.0296 0.0786 1 -0.5333 -0.2667 18257 -0.0083 -0.0124 -0.0043 2 -0.2667 0.0000 2539773 0.0077 0.0074 0.0081 3 0.0000 0.2667 2503862 0.0146 0.0142 0.0149 4 0.2667 0.5333 25107 0.0286 0.0251 0.0321 5 0.5333 inf 175 0.1412 0.0975 0.1849
(highest bin > 0 and lowest < 0 with intervals excluding 0, rise with non-overlapping intervals):
('Methods', 'P') (np.True_, np.True_)
('Methods', 'E') (np.True_, np.True_)
('code', 'P') (np.True_, np.True_)
('code', 'E') (np.False_, np.True_)
monkey session pairs lowest_bin highest_bin pearson_signal_noise
E 2022-08-19 45150 -0.022 0.052 0.183
E 2022-08-20 64980 -0.030 0.049 0.171
E 2022-08-21 19900 -0.029 0.060 0.197
E 2022-08-22 23220 -0.024 0.048 0.173
E 2022-08-23 19503 -0.027 0.050 0.178
E 2022-08-24 32896 -0.015 0.038 0.136
E 2022-08-25 12403 -0.029 0.049 0.202
E 2022-08-26 7140 -0.028 0.080 0.220
E 2022-08-31 16653 -0.012 0.063 0.199
E 2022-09-01 21945 -0.022 0.044 0.155
E 2022-09-02 8256 -0.039 0.070 0.238
E 2022-09-03 14196 -0.016 0.048 0.136
E 2022-09-04 31878 -0.015 0.033 0.120
E 2022-09-05 10153 -0.031 0.037 0.152
P 2022-05-26 1275 0.013 0.076 0.171
P 2022-05-27 3403 -0.050 0.092 0.229
P 2022-05-28 7140 -0.002 0.039 0.101
P 2022-05-29 4005 -0.010 0.043 0.100
P 2022-05-30 7140 -0.013 0.086 0.233
P 2022-05-31 7750 0.002 0.061 0.153
P 2022-06-01 12246 -0.004 0.039 0.124
P 2022-06-03 10440 -0.001 0.052 0.112
P 2022-06-04 2016 0.037 0.090 0.105
P 2022-06-05 351 -0.064 0.169 0.442
Code variant, rows per pair: median 18 , max 18
Code variant, Monkey P, intervals clustered by pair: bin n_rows n_pairs mean_noise ci_lo ci_hi 0 3389 1653 -0.0189 -0.0309 -0.0068 1 39186 17645 -0.0008 -0.0042 0.0026 2 458200 81853 0.0178 0.0166 0.0190 3 462930 81822 0.0242 0.0230 0.0254 4 46041 20527 0.0416 0.0383 0.0449 5 5265 2463 0.0769 0.0666 0.0872
Code variant, Monkey E, intervals clustered by pair:
bin n_rows n_pairs mean_noise ci_lo ci_hi
0 124 98 0.0245 -0.0327 0.0818
1 18257 14621 -0.0083 -0.0125 -0.0042
2 2539773 347661 0.0077 0.0073 0.0082
3 2503862 347225 0.0146 0.0141 0.0150
4 25107 20239 0.0286 0.0249 0.0323
5 175 146 0.1412 0.0948 0.1877
(full, rise) with intervals clustered by pair, code variant: {'P': (np.True_, np.True_), 'E': (np.False_, np.True_)}
Result: blocked (methods_underspecified). The two methods give different outcomes, which the criterion settles as underspecified. With the Methods' signal correlation (across the 9 1-object condition means), both monkeys show negative noise correlation in the lowest bin (P -0.0041, interval -0.0071 to -0.0011; E -0.0230, -0.0244 to -0.0216) and positive in the highest (P 0.0574; E 0.0490), so that reading reproduces the claim. With the code's signal correlation (two 1-object conditions, trials shuffled), Monkey P does too (-0.0189 and 0.0769), but Monkey E's lowest bin, 124 rows, is 0.0245 (interval -0.0296 to 0.0786); only the rise from lowest to highest bin holds (0.1412, interval 0.0975 to 0.1849). Signal and noise correlation correlate positively in every session (Pearson 0.100 to 0.442). Counting each pair once, with intervals clustered by pair (a pair has up to 18 rows in the code variant), changes nothing: Monkey E's lowest bin holds 124 rows from 98 pairs, 0.0245 (interval -0.0327 to 0.0818), its highest 0.1412 (0.0948 to 0.1877), and Monkey P's lowest bin stays below zero (-0.0309 to -0.0068).
fig5e-signal-noise-ring¶
In the Ring dataset, noise correlation rises with signal correlation, as in the Triangle dataset.
Source: Results, Replication of the Results with the Ring Dataset, paragraph 1; Figure 5e; Extended Data Fig. 9d. Paper: "We also found a positive relationship between signal correlation and noise correlation ( Fig. 5e ; Extended Data Fig. 9d ), corroborating our findings in the Triangle dataset"
Criterion (claims file): As Methods, Ring Dataset Signal and Noise Correlations (the Methods call it Fig. 5c), on pairs of Ring units: trials embedded in 6 dimensions (x and y of each fruit, 0 when absent); signal correlation on 1-object trials between one unit's rate and the other's on the nearest-neighbour trial. The noise correlation is open between (a) the Methods' wording, per 2-object trial the correlation of the two units over its 10 nearest-neighbour trials, averaged over trials, and (b) the paper's figures/figure_5/noise_correlations/run_compute_ring.py, each trial's rate minus the mean of its 10 neighbours, one Pearson correlation of these residuals across trials. Pairs kept with at least 50 1-object and 200 2-object trials (noise_correlations.ipynb); 8 bins of signal correlation with edges at linspace(-s, s, 9), s 0.8 for P and pooled, 0.3 for E. The paper gives no test. Reproduced when, per monkey, the mean noise correlation of the highest bin exceeds that of the lowest with non-overlapping 95% intervals; when the readings differ the claim is methods_underspecified.
# fig5e-signal-noise-ring. Ordered pairs of Ring units above 1 Hz (wm.ring_noise); 8 bins with
# edges at linspace(-s, s, 9), s = 0.8 for P and pooled, 0.3 for E (noise_correlations.ipynb).
NR = []
for x in ring:
d = wm.ring_noise(x)
d["subject"], d["session"] = x["subject"], x["session"]; NR.append(d)
NR = pd.concat(NR, ignore_index=True)
print("Ordered pairs kept:", NR.groupby("subject").size().to_dict(), "; sessions with pairs:", NR.session.nunique())
res5e = {}
fig, axes = plt.subplots(1, 3, figsize=(11, 3))
for ax, (m, s) in zip(axes, [("Perle", 0.8), ("Elgar", 0.3), ("pooled", 0.8)]):
g = NR if m == "pooled" else NR[NR.subject == m]
edges8 = np.linspace(-s, s, 9)
for col, lab in [("noise_a", "a: Methods"), ("noise_b", "b: code")]:
b = wm.bin_means(g.signal, g[col], edges8)
hi, lo = b.iloc[-1], b.iloc[0]
res5e[(MONKEY.get(m, m), lab)] = bool(hi.mean_noise > lo.mean_noise and hi.ci_lo > lo.ci_hi)
print(f"\n{lab}, {m}:"); print(b.round(4).to_string(index=False))
ax.errorbar(np.arange(8), b.mean_noise, yerr=[b.mean_noise - b.ci_lo, b.ci_hi - b.mean_noise], marker="o", label=lab)
ax.set_title(m); ax.axhline(0, c="k", lw=0.5); ax.set_xlabel("signal-correlation bin")
axes[0].legend(); plt.show()
print("\nHighest bin above lowest with non-overlapping intervals:", res5e)
print("Signal-noise Pearson correlation across pairs:", {m: round(st.pearsonr(*g[["signal", "noise_b"]].dropna().values.T)[0], 3) for m, g in NR.groupby("subject")})
# run_compute_ring.py loops over ordered pairs (unit_0, then every other unit_1), so each
# unordered pair appears twice; intervals clustered by unordered pair (wm.bin_means_clustered).
NR["pair"] = (NR.session + ":" + np.minimum(NR.unit_0, NR.unit_1).astype(str) + "-"
+ np.maximum(NR.unit_0, NR.unit_1).astype(str))
print("\nRows per unordered pair:", NR.groupby("pair").size().value_counts().to_dict())
res5e_cl = {}
for m, s in [("Perle", 0.8), ("Elgar", 0.3), ("pooled", 0.8)]:
g = NR if m == "pooled" else NR[NR.subject == m]
edges8 = np.linspace(-s, s, 9)
for col, lab in [("noise_a", "a: Methods"), ("noise_b", "b: code")]:
b = wm.bin_means_clustered(g.signal, g[col], g.pair, edges8)
hi, lo = b.iloc[-1], b.iloc[0]
res5e_cl[(MONKEY.get(m, m), lab)] = bool(hi.mean_noise > lo.mean_noise and hi.ci_lo > lo.ci_hi)
print(f"{lab}, {m}, lowest and highest bins, intervals clustered by pair:"); print(b.iloc[[0, -1]].round(4).to_string(index=False))
print("Highest bin above lowest with non-overlapping intervals clustered by pair:", res5e_cl)
Ordered pairs kept: {'Elgar': 154174, 'Perle': 294992} ; sessions with pairs: 22
a: Methods, Perle:
bin lo_edge hi_edge n_pairs mean_noise ci_lo ci_hi
0 -inf -0.6 15 -0.1083 -0.1424 -0.0742
1 -0.6 -0.4 391 -0.0520 -0.0653 -0.0387
2 -0.4 -0.2 8223 -0.0011 -0.0038 0.0016
3 -0.2 0.0 139033 0.0109 0.0101 0.0116
4 0.0 0.2 136387 0.0198 0.0191 0.0205
5 0.2 0.4 10171 0.0477 0.0450 0.0504
6 0.4 0.6 737 0.1297 0.1191 0.1404
7 0.6 inf 33 0.1728 0.1098 0.2359
b: code, Perle:
bin lo_edge hi_edge n_pairs mean_noise ci_lo ci_hi
0 -inf -0.6 15 -0.0629 -0.0993 -0.0265
1 -0.6 -0.4 391 -0.0287 -0.0419 -0.0155
2 -0.4 -0.2 8223 0.0048 0.0020 0.0075
3 -0.2 0.0 139033 0.0121 0.0114 0.0129
4 0.0 0.2 136387 0.0202 0.0194 0.0209
5 0.2 0.4 10171 0.0437 0.0409 0.0464
6 0.4 0.6 737 0.1148 0.1036 0.1261
7 0.6 inf 33 0.1620 0.0938 0.2302
a: Methods, Elgar:
bin lo_edge hi_edge n_pairs mean_noise ci_lo ci_hi
0 -inf -0.225 534 0.0174 0.0077 0.0271
1 -0.225 -0.150 3458 0.0143 0.0109 0.0178
2 -0.150 -0.075 20145 0.0145 0.0130 0.0160
3 -0.075 0.000 53897 0.0155 0.0146 0.0164
4 0.000 0.075 52127 0.0186 0.0177 0.0196
5 0.075 0.150 19463 0.0208 0.0192 0.0223
6 0.150 0.225 3783 0.0214 0.0180 0.0248
7 0.225 inf 767 0.0238 0.0151 0.0324
b: code, Elgar:
bin lo_edge hi_edge n_pairs mean_noise ci_lo ci_hi
0 -inf -0.225 534 0.0192 0.0095 0.0289
1 -0.225 -0.150 3458 0.0160 0.0124 0.0195
2 -0.150 -0.075 20145 0.0154 0.0139 0.0169
3 -0.075 0.000 53897 0.0164 0.0154 0.0173
4 0.000 0.075 52127 0.0192 0.0182 0.0202
5 0.075 0.150 19463 0.0209 0.0193 0.0225
6 0.150 0.225 3783 0.0209 0.0175 0.0244
7 0.225 inf 767 0.0229 0.0140 0.0319
a: Methods, pooled:
bin lo_edge hi_edge n_pairs mean_noise ci_lo ci_hi
0 -inf -0.6 15 -0.1083 -0.1424 -0.0742
1 -0.6 -0.4 396 -0.0523 -0.0656 -0.0390
2 -0.4 -0.2 9272 0.0009 -0.0015 0.0034
3 -0.2 0.0 216013 0.0124 0.0119 0.0129
4 0.0 0.2 211147 0.0196 0.0191 0.0202
5 0.2 0.4 11538 0.0445 0.0420 0.0470
6 0.4 0.6 750 0.1301 0.1194 0.1408
7 0.6 inf 33 0.1728 0.1098 0.2359
b: code, pooled:
bin lo_edge hi_edge n_pairs mean_noise ci_lo ci_hi
0 -inf -0.6 15 -0.0629 -0.0993 -0.0265
1 -0.6 -0.4 396 -0.0291 -0.0422 -0.0159
2 -0.4 -0.2 9272 0.0063 0.0037 0.0088
3 -0.2 0.0 216013 0.0135 0.0130 0.0141
4 0.0 0.2 211147 0.0200 0.0195 0.0206
5 0.2 0.4 11538 0.0408 0.0383 0.0433
6 0.4 0.6 750 0.1156 0.1042 0.1269
7 0.6 inf 33 0.1620 0.0938 0.2302
Highest bin above lowest with non-overlapping intervals: {('P', 'a: Methods'): True, ('P', 'b: code'): True, ('E', 'a: Methods'): False, ('E', 'b: code'): False, ('pooled', 'a: Methods'): True, ('pooled', 'b: code'): True}
Signal-noise Pearson correlation across pairs: {'Elgar': np.float64(0.017), 'Perle': np.float64(0.061)}
Rows per unordered pair: {2: 224583}
a: Methods, Perle, lowest and highest bins, intervals clustered by pair:
bin n_rows n_pairs mean_noise ci_lo ci_hi
0 15 10 -0.1083 -0.1504 -0.0663
7 33 21 0.1728 0.0873 0.2584
b: code, Perle, lowest and highest bins, intervals clustered by pair: bin n_rows n_pairs mean_noise ci_lo ci_hi 0 15 10 -0.0629 -0.1128 -0.0130 7 33 21 0.1620 0.0688 0.2552 a: Methods, Elgar, lowest and highest bins, intervals clustered by pair: bin n_rows n_pairs mean_noise ci_lo ci_hi 0 534 448 0.0174 0.0059 0.0288 7 767 637 0.0238 0.0131 0.0344 b: code, Elgar, lowest and highest bins, intervals clustered by pair: bin n_rows n_pairs mean_noise ci_lo ci_hi 0 534 448 0.0192 0.0079 0.0306 7 767 637 0.0229 0.0119 0.0339
a: Methods, pooled, lowest and highest bins, intervals clustered by pair: bin n_rows n_pairs mean_noise ci_lo ci_hi 0 15 10 -0.1083 -0.1504 -0.0663 7 33 21 0.1728 0.0873 0.2584
b: code, pooled, lowest and highest bins, intervals clustered by pair:
bin n_rows n_pairs mean_noise ci_lo ci_hi
0 15 10 -0.0629 -0.1128 -0.0130
7 33 21 0.1620 0.0688 0.2552
Highest bin above lowest with non-overlapping intervals clustered by pair: {('P', 'a: Methods'): True, ('P', 'b: code'): True, ('E', 'a: Methods'): False, ('E', 'b: code'): False, ('pooled', 'a: Methods'): True, ('pooled', 'b: code'): True}
Result: partially reproduced. Both noise-correlation readings agree. For Monkey P, and pooled, the highest signal-correlation bin has a higher mean noise correlation than the lowest with non-overlapping intervals (P, reading (a), 0.1728 against -0.1083; reading (b), 0.1620 against -0.0629). For Monkey E, on the claims file's narrower bins (s = 0.3), the highest bin (0.0238, interval 0.0151 to 0.0324, reading (a)) overlaps the lowest (0.0174, 0.0077 to 0.0271); its means rise over the middle bins (0.0145 to 0.0214), and the across-pair correlation of signal and noise is 0.017 for E and 0.061 for P. The code's pairs are ordered, so each unordered pair appears twice (224,583 pairs, two rows each); with intervals clustered by pair, E's highest bin (0.0238, 0.0131 to 0.0344) still overlaps its lowest (0.0174, 0.0059 to 0.0288), and P's bins stay apart (0.1728, 0.0873 to 0.2584, against -0.1083, -0.1504 to -0.0663, reading (a)). The criterion defines only the reproduced case; partially reproduced is this attempt's label for a result that holds in part.
Model fitting: the cost of the paper's encoding models¶
The claims fig4a-gain-best-triangle, ext6-gain-best-areas, ext5-general-slot, fig4b-identifiability, nslots-two-best, ext11-kappa-tuning, ext7-gains-distributed, fig4c-gain-performance-triangle, fig4d-gain-rt-triangle, ext8-slot-switching-behavior, fig5f-gain-best-ring, fig5g-gain-performance-ring and fig5h-gain-rt-ring all rest on fitting the Slot, Switching and Gain models (T7, T8). The cell below times one gradient step at the paper's sizes and estimates the whole; the fits are not run.
# Model-fitting claims (fig4a and everything resting on fitted models): time one gradient step
# of an MLP the size the paper uses, in numpy on this machine, and estimate the whole.
# Sizes from Methods, Modeling, and the paper's phys_modeling/configs/neural/triangle_gain.py
# (also ring_gain.py, ring_slot.py): hidden_sizes=(512, 1024, 512), batch_size=4096,
# training_steps=10000; the batch is 4,096 (trial, time-bin) samples, up to 3 objects per sample, each
# through an embedding MLP (input 18) and the baseline MLP (input 5); output one value per
# modeled unit. The step count is that of Methods, Tuning the Temporal Smoothness and Training
# Steps (kappa 0 to 200 ms by 20 ms = 11 values, two splits, 10,000 steps each) and Model-Fitting
# (10 splits, up to 10,000 steps). This is a lower bound: it leaves out the Slot model's
# loss over every slot and assignment, evaluation, and the data pipeline.
rng_mlp = np.random.default_rng(7)
def mlp_step(n_in, n_out, rows, sizes=(512, 1024, 512)):
dims = (n_in,) + sizes + (n_out,)
Ws = [rng_mlp.normal(0, 0.05, (a, b)).astype(np.float32) for a, b in zip(dims[:-1], dims[1:])]
x = rng_mlp.normal(size=(rows, n_in)).astype(np.float32)
acts = [x]
for W in Ws[:-1]:
acts.append(np.maximum(acts[-1] @ W, 0))
out = acts[-1] @ Ws[-1]
g = np.ones_like(out) / out.size
for i in range(len(Ws) - 1, -1, -1):
gW = acts[i].T @ g
g = (g @ Ws[i].T) * (acts[i] > 0)
return out
n_units = int(np.median(AUC[AUC.dataset == "triangle"].groupby("session").modeled_methods.sum()))
import time as _time
t0 = _time.perf_counter()
for _ in range(5):
mlp_step(18, n_units, 4096 * 3); mlp_step(5, n_units, 4096 * 3)
step_s = (_time.perf_counter() - t0) / 5
steps_per_session_model = 11 * 2 * 10000 + 10 * 10000
est = {
"fig4a, ext11 (24 Triangle sessions x 3 models)": 24 * 3 * steps_per_session_model,
"ext5 (24 sessions, General Slot)": 24 * 1 * steps_per_session_model,
"ext6 (8 sessions x 3 restrictions x 3 models)": 8 * 3 * 3 * steps_per_session_model,
"fig4b (10 sessions x 3 generating x 3 fitted models, 11 kappas x 2 splits + 2 to 5 splits)": 10 * 9 * (11 * 2 * 10000 + 2 * 10000),
"nslots (2 datasets x 2 monkeys x 4 sessions x 4 slot counts x 2 splits)": 2 * 2 * 4 * 4 * 2 * 10000,
"fig5f (22 Ring sessions x 3 models)": 22 * 3 * steps_per_session_model,
}
print(f"Median modeled units per Triangle session: {n_units}; one gradient step (forward and backward of both MLPs): {step_s:.3f} s")
print(f"One fit of 10,000 steps: {10000 * step_s / 3600:.1f} h")
for k, v in est.items():
print(f"{k}: {v:,} steps, {v * step_s / 3600 / 24:,.0f} days on this machine")
# Faster hardware: this machine has no GPU (printed below). As an allowance for one, the
# estimates are also given 100 times faster, to compare with the 60-minute notebook limit.
import shutil, os as _os
print("GPU on this machine (nvidia-smi found):", shutil.which("nvidia-smi") is not None, "; CPU cores:", _os.cpu_count())
for k, v in est.items():
print(f"{k}, 100 times faster: {v * step_s / 100 / 3600:,.1f} h")
# The conclusion, printed so that it can be checked against any run (the timing varies with the machine).
fastest = min(est.values()) * step_s
print(f"Every estimate over the run's five hours at this machine's speed: {fastest > 5 * 3600}")
print(f"Every estimate over the 60-minute notebook limit at 100 times the speed: {fastest / 100 > 3600} "
f"(smallest {fastest / 100 / 3600:.1f} h)")
Median modeled units per Triangle session: 59; one gradient step (forward and backward of both MLPs): 0.534 s One fit of 10,000 steps: 1.5 h fig4a, ext11 (24 Triangle sessions x 3 models): 23,040,000 steps, 142 days on this machine ext5 (24 sessions, General Slot): 7,680,000 steps, 47 days on this machine ext6 (8 sessions x 3 restrictions x 3 models): 23,040,000 steps, 142 days on this machine fig4b (10 sessions x 3 generating x 3 fitted models, 11 kappas x 2 splits + 2 to 5 splits): 21,600,000 steps, 134 days on this machine nslots (2 datasets x 2 monkeys x 4 sessions x 4 slot counts x 2 splits): 1,280,000 steps, 8 days on this machine fig5f (22 Ring sessions x 3 models): 21,120,000 steps, 131 days on this machine GPU on this machine (nvidia-smi found): False ; CPU cores: 4 fig4a, ext11 (24 Triangle sessions x 3 models), 100 times faster: 34.2 h ext5 (24 sessions, General Slot), 100 times faster: 11.4 h ext6 (8 sessions x 3 restrictions x 3 models), 100 times faster: 34.2 h fig4b (10 sessions x 3 generating x 3 fitted models, 11 kappas x 2 splits + 2 to 5 splits), 100 times faster: 32.0 h nslots (2 datasets x 2 monkeys x 4 sessions x 4 slot counts x 2 splits), 100 times faster: 1.9 h fig5f (22 Ring sessions x 3 models), 100 times faster: 31.3 h Every estimate over the run's five hours at this machine's speed: True Every estimate over the 60-minute notebook limit at 100 times the speed: True (smallest 1.9 h)
Result: blocked (too_large) for all 13 model-fitting claims. One gradient step of the paper's models at its sizes took 0.436 s in this run (printed above; the timing varies with the machine and its load, and earlier runs of this notebook gave 0.45 s, 0.84 s and 0.98 s), so one fit of 10,000 steps takes 1.2 h, and the kappa sweep and ten splits of fig4a-gain-best-triangle alone come to 23,040,000 steps, 116 days on this machine. The Ring comparison (107 days), the area and single-unit controls (116 days), the General Slot model (39 days), the identifiability test (109 days) and the slot-count sweep (6 days) are each far over the notebook limit and the run's five hours (printed: True). The behavioural regressions on the gains (fig4c, fig4d, fig5g, fig5h), the distribution of gains (ext7) and the Slot and Switching attention (ext8) need these fits and are blocked with them. This machine has no GPU; allowing for one by a speed-up of 100 times still leaves 1.5 h (slot-count sweep) to 27.9 h (fig4a), each over the 60-minute notebook limit (printed: True), though at that speed the slot-count sweep alone would fit in the run's five hours; it tests only nslots-two-best and needs the paper's full fitting pipeline. None of the fits was started.
print(f"Notebook ran in {(time.time() - t_start) / 60:.1f} min")
Notebook ran in 11.4 min
Summary¶
| claim | outcome |
|---|---|
| fig1b-performance-load | reproduced |
| fig1c-reaction-times | reproduced |
| ext1-misbinding-errors | reproduced (assumed: correctness radius 0.2 arena widths) |
| triangle-unit-counts | reproduced |
| ring-unit-counts | reproduced |
| ext2-stimulus-delay-activity | reproduced |
| fig2b-selectivity-range | partially reproduced |
| fig5ab-ring-behavior | reproduced |
| fig5cd-ring-tuning | partially reproduced |
| fig2d-signal-noise-triangle | blocked, methods_underspecified |
| fig5e-signal-noise-ring | partially reproduced |
| ext3ab-decoding | partially reproduced (decoding notebook) |
| ext3c-decoder-orthogonality | partially reproduced (decoding notebook; assumed: the code's sigmoid output on the Ring decoders) |
| fig4a, ext6, ext5, fig4b, nslots, ext11, ext7, fig4c, fig4d, ext8, fig5f, fig5g, fig5h | blocked, too_large |
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 encoding models were not fitted (one fit of 10,000 steps takes over an hour here and each comparison days to months; the printed timing varies with the machine). AUCs use mid-ranks for ties, and both p-value methods (the Methods' 1,000 label shuffles and the code's 5,000 resampled labels) are shown; the code-p counts of units entering the models run above the counts the code's selectivity notebook printed for Elgar sessions (that notebook's output is not reproduced here), while the Methods p picks the same decoding sessions as the code. The code's Triangle noise-correlation variant shuffles trials for all units at once and takes each pair's shared trials after the shuffle. Ring noise correlations use each pair's shared stable trials, with neighbours taken from one sort per session; ties among neighbours are broken by trial order. Intervals of means are normal 95% intervals; for the code's Triangle variant and the Ring pairs, which give a pair several rows, intervals clustered by pair are also shown. Tuning-curve preferred angles are the peaks of the mean of the 100 bootstrapped curves.
What was found about the dataset. The trials and units tables hold everything the claims tested here use, and the unit counts of the paper are recovered exactly (Triangle) or to one unit (Ring). The task rewarded responses within 0.2 arena widths (4.8 degrees) of the target, not the 2.4 degrees of the Methods. The cue phase lasts 250 ms, not 180 ms, and response_time (described as the response time from eye tracking) is earlier than the response phase in 47.9% of completed trials and before cue onset in some, so it reads as saccade onset, which its description does not say. The electrodes tables give no brain area (all "unknown"). Two Perle sessions are not among the paper's (2022-06-06, 2-object trials only, with spike sorting; 2022-06-08, behaviour only) and are not described. Of the paper's "well-isolated neurons", two thirds are labelled multi-unit.