Reproducing Boran et al. (2022) from dandiset 000575¶
Paper. Boran E, Hilfiker P, Stieglitz L, Sarnthein J, Klaver P (2022). Persistent neuronal firing
in the medial temporal lobe supports performance and workload of visual working memory in humans.
NeuroImage 254, 119123. doi:10.1016/j.neuroimage.2022.119123.
The claims file was written from the version of record (Methods, Results, figure legends and
Table 1; the paper has no supplement) and is claims/000575/10.1016-j.neuroimage.2022.119123.yaml.
Dandiset. DANDI:000575, version 0.231010.1811 (doi:10.48324/dandi.000575/0.231010.1811). The paper cites no version of the dandiset (its data statement names the GIN repository USZ_NCH/Human_MTL_units_visual_WM, of which the dandiset is a converted form). 0.231010.1811 is the only published version; the draft lists the same 18 assets with the same total size. The files are streamed; nothing is downloaded.
What it holds. 18 NWB files for 16 sessions of 13 participants (two pairs of files of participant 12 are copies): trial tables with set size and the correct answer, the answer per trial, task epochs, 809 sorted units with spike times inside the trial windows and their mean waveforms, and, in 12 files, LFP (not used).
Outside inputs. None for the analysis. The conversion code named in the dandiset's metadata (github.com/janhohenheim/usz-neuro-conversion, commit c31dad9) was read to interpret the response timestamps (data check); no values are taken from it.
| Claim | Tested | Outcome |
|---|---|---|
| behavior-accuracy-by-set-size | yes | reproduced |
| behavior-capacity | yes | reproduced |
| table1-capacity-and-neurons | yes | reproduced |
| behavior-rt-set-size | yes | partially reproduced |
| neurons-recorded | yes | reproduced (assumed) |
| spike-sorting-quality | yes | not reproduced |
| firing-rates-by-period | yes | reproduced |
| fig2c-maintenance-neurons | yes | reproduced (assumed) |
| fig2d-workload-modulated-maintenance-neurons | yes | not reproduced |
| maintenance-neurons-correct-vs-incorrect | yes | not reproduced |
| fig3b-test-neurons | yes | partially reproduced |
| fig4a-dpca-separates-set-sizes | yes | not reproduced |
| fig4b-speed-lowest-in-maintenance | yes | reproduced |
| fig4c-distance-largest-in-maintenance | yes | partially reproduced |
| fig4d-decode-workload-all | yes | not reproduced |
| fig4d-decode-workload-regions | yes | not reproduced |
| decode-performance-not-significant | yes | reproduced |
The notebook takes 16 to 27 minutes, depending on the machine, most of it in the permutation nulls, the dPCA runs and the decoders; reading the files takes about 20 s. Random seeds are fixed, one generator per section.
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from scipy.ndimage import convolve1d
import mtl_wm as m
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
NOTEBOOK_START = time.time()
started = time.time()
files, duplicates, S = m.load_all()
print(f"read {len(files)} files of dandiset {m.DANDISET_ID} version {m.DANDISET_VERSION} in {time.time() - started:.0f} s")
for s in S:
s["counts"] = m.period_counts(s)
UNITS = pd.concat([s["units"].assign(key=s["key"], participant=int(s["subject"]), session=s["session"]) for s in S],
ignore_index=True)
TRIALS = pd.concat([s["trials"].assign(key=s["key"], participant=int(s["subject"])) for s in S], ignore_index=True)
read 18 files of dandiset 000575 version 0.231010.1811 in 15 s
Data check¶
Each check of the metadata against the values is stated under its output.
print("Files in the version:", len(files))
print("Files that repeat another file (same trial starts, spike counts and first spikes), counted once:")
for a, b in duplicates:
print(f" {a}\n = {b}")
info = pd.DataFrame([{
"key": s["key"], "path": s["path"].split("/")[1][:40], "identifier": s["identifier"].split("WM_")[1],
"session_start": s["session_start"][:16], "trials": len(s["trials"]), "correct": int(s["trials"].correct.sum()),
"units": len(s["units"]), "spikes": int(s["units"].n_spikes.sum()), "has_LFP": s["has_lfp"]} for s in S])
print(f"\nSessions kept: {len(S)}; participants: {info.key.str[:3].nunique()}")
print(info.to_string(index=False))
print("\nTrials in all:", len(TRIALS), "; correct:", int(TRIALS.correct.sum()), "; paper: 192 trials per session, 2678 correct trials")
print("Trials per session, min and max:", info.trials.min(), info.trials.max())
Files in the version: 18
Files that repeat another file (same trial starts, spike counts and first spikes), counted once:
sub-12/sub-12_ses-20161217T103800_obj-cpvbt2_behavior+ecephys.nwb
= sub-12/sub-12_ses-20161217T103800_obj-16z48k1_behavior+ecephys.nwb
sub-12/sub-12_ses-20161217T103800_obj-p4v231_behavior+ecephys.nwb
= sub-12/sub-12_ses-20161217T103800_obj-1h9ryn9_behavior+ecephys.nwb
Sessions kept: 16; participants: 13
key path identifier session_start trials correct units spikes has_LFP
P01S1 sub-01_ses-20140827T193800_behavior+ecep subject01_session01 2014-08-27T19:38 189 154 33 55038 True
P02S1 sub-02_ses-20171011T152100_behavior+ecep subject02_session01 2017-10-11T15:21 192 159 36 45296 False
P03S1 sub-03_ses-20170113T105000_behavior+ecep subject03_session01 2017-01-13T10:50 192 159 52 78193 True
P03S2 sub-03_ses-20170113T111800_behavior+ecep subject03_session02 2017-01-13T11:18 191 163 54 97892 False
P04S1 sub-04_ses-20170118T161000_behavior+ecep subject04_session01 2017-01-18T16:10 191 170 85 250782 True
P05S1 sub-05_ses-20180704T132100_behavior+ecep subject05_session01 2018-07-04T13:21 192 166 29 70068 True
P06S1 sub-06_ses-20171027T113600_behavior+ecep subject06_session01 2017-10-27T11:36 192 172 96 204755 True
P07S1 sub-07_ses-20181027T183600_behavior+ecep subject07_session01 2018-10-27T18:36 189 167 48 73676 True
P08S1 sub-08_ses-20190116T110300_behavior+ecep subject08_session01 2019-01-16T11:03 190 170 87 178178 True
P08S2 sub-08_ses-20190116T125000_behavior+ecep subject08_session02 2019-01-16T12:50 191 168 63 133575 True
P09S1 sub-09_ses-20150305T112000_behavior+ecep subject09_session01 2015-03-05T11:20 174 153 23 71066 True
P10S1 sub-10_ses-20160303T114400_behavior+ecep subject10_session01 2016-03-03T11:44 188 169 62 209000 True
P11S1 sub-11_ses-20190515T195800_behavior+ecep subject11_session01 2019-05-15T19:58 192 177 11 21645 False
P12S1 sub-12_ses-20161217T103800_obj-1h9ryn9_b subject12_session01 2016-12-17T10:38 191 178 29 33670 True
P12S2 sub-12_ses-20161217T103800_obj-16z48k1_b subject12_session02 2016-12-17T10:38 192 177 28 35008 False
P13S1 sub-13_ses-20160213T110700_behavior+ecep subject13_session01 2016-02-13T11:07 192 176 73 136837 True
Trials in all: 3038 ; correct: 2678 ; paper: 192 trials per session, 2678 correct trials
Trials per session, min and max: 174 192
The version holds 18 files for 16 sessions: participant 12 has four files with one session start
time, and they form two pairs of identical recordings (the check above compares trial starts,
spike counts and first spikes). The file identifiers (subject12_session01 and session02) name
the two sessions; each is counted once. Both sessions of participant 12 carry the same start time
(10:38), so the start time does not tell them apart. Every session has 174 to 192 trials, where the
paper and the file descriptions state 192 per session; the trial ids run without gaps, so the
missing trials are absent from the files, not marked. The correct trials number 2678, the count
the paper gives. Four files (participants 2, 3 session 2, 11 and 12 session 2) hold no LFP.
# Periods: the epochs table against the windows the helper uses (trial start = fixation onset).
rows = []
for s in S:
ep, tr = s["epochs"], s["trials"]
starts = {tag: ep[ep.tag == tag].start.values for tag in ["Fixation", "Memory", "Retention", "Probe"]}
st = tr.start.values
rows.append({"key": s["key"], "n_fixation_epochs": len(starts["Fixation"]),
"max |fixation - trial start| (s)": np.abs(starts["Fixation"] - st).max(),
"memory - start": np.unique(np.round(starts["Memory"] - st, 4)).tolist(),
"retention - start": np.unique(np.round(starts["Retention"] - st, 4)).tolist(),
"probe - start": np.unique(np.round(starts["Probe"] - st, 4)).tolist(),
"trial length": np.unique(np.round(tr.stop.values - st, 4)).tolist(),
"gap to next trial (s), min": np.round((st[1:] - tr.stop.values[:-1]).min(), 3)})
print(pd.DataFrame(rows).to_string(index=False))
print("\nHelper periods (s from trial start):", m.PERIODS)
inside = []
for s in S:
st, sp = s["trials"].start.values, s["trials"].stop.values
for t in s["spike_times"]:
i = np.searchsorted(st, t, "right") - 1
inside.append(np.mean((i >= 0) & (t <= sp[np.clip(i, 0, None)])))
print(f"Fraction of each unit's spikes inside a trial window: min {min(inside):.4f}, max {max(inside):.4f}")
print("Observation intervals per unit equal the number of trials in every session:",
all(s["n_obs_intervals"] == len(s["units"]) * len(s["trials"]) for s in S))
key n_fixation_epochs max |fixation - trial start| (s) memory - start retention - start probe - start trial length gap to next trial (s), min
P01S1 189 0.0 [0.9] [1.7] [2.6] [4.6] 1.194
P02S1 192 0.0 [0.9] [1.7] [2.6] [4.6] 1.194
P03S1 192 0.0 [0.9] [1.7] [2.6] [4.6] 1.194
P03S2 191 0.0 [0.9] [1.7] [2.6] [4.6] 1.194
P04S1 191 0.0 [0.9] [1.7] [2.6] [4.6] 1.194
P05S1 192 0.0 [0.9] [1.7] [2.6] [4.6] 1.211
P06S1 192 0.0 [0.9] [1.7] [2.6] [4.6] 1.211
P07S1 189 0.0 [0.9] [1.7] [2.6] [4.6] 1.194
P08S1 190 0.0 [0.9] [1.7] [2.6] [4.6] 1.194
P08S2 191 0.0 [0.9] [1.7] [2.6] [4.6] 1.194
P09S1 174 0.0 [0.9] [1.7] [2.6] [4.6] 1.111
P10S1 188 0.0 [0.9] [1.7] [2.6] [4.6] 1.211
P11S1 192 0.0 [0.9] [1.7] [2.6] [4.6] 1.194
P12S1 191 0.0 [0.9] [1.7] [2.6] [4.6] 1.194
P12S2 192 0.0 [0.9] [1.7] [2.6] [4.6] 1.194
P13S1 192 0.0 [0.9] [1.7] [2.6] [4.6] 1.045
Helper periods (s from trial start): {'fixation': (0.0, 0.9), 'encoding': (0.9, 1.7), 'maintenance': (1.7, 2.6), 'test': (2.6, 4.6)}
Fraction of each unit's spikes inside a trial window: min 1.0000, max 1.0000 Observation intervals per unit equal the number of trials in every session: True
In every session the epochs give fixation at trial start, the memory array at 0.9 s, the retention (maintenance) period at 1.7 s and the test array (probe) at 2.6 s, and every trial ends 2.0 s after test onset (4.6 s). These are the periods of Methods 2.2 and 2.3 and of the helper. The fixation period is a fixed 0.9 s in the files, where the task description has a 0.4 s warning signal and a 0.4 to 0.5 s jittered black dot. All spikes lie inside trial windows, and each unit's observation intervals are the trial windows: the files hold spikes only from fixation onset to test offset, not during the inter-trial intervals.
# Units: the microwire label and the location column.
print("Waveforms:", S[0]["waveform_mean"].shape[1], "samples at", S[0]["waveform_rate"], "Hz, unit:", S[0]["waveform_unit"])
wm = np.concatenate([s["waveform_mean"] for s in S]); wsd = np.concatenate([s["waveform_sd"] for s in S])
print(f"Peak |mean waveform| over units: {np.abs(wm).max(1).min() * 1e6:.1f} to {np.abs(wm).max(1).max() * 1e6:.1f} uV; "
f"waveform SD {wsd.min() * 1e6:.1f} to {wsd.max() * 1e6:.1f} uV")
print("\nUnits per microwire bundle (label without the wire number):")
print(UNITS.groupby("bundle").size().to_string())
print("\nBundles not named in the implanted electrodes of the subject description:")
for s in S:
implanted = s["subject_description"].split("Implanted electrodes: ")[1].split("\n")[0].split(",")
extra = sorted(set(s["units"].bundle) - set(implanted))
if extra:
print(f" {s['key']}: units on {extra}; implanted: {implanted}")
Waveforms: 64 samples at 32000.0 Hz, unit: volts Peak |mean waveform| over units: 15.0 to 152.5 uV; waveform SD 1.6 to 30.1 uV Units per microwire bundle (label without the wire number): bundle AHL 159 AHR 97 AL 91 AR 135 ECL 71 ECR 146 HL 2 HR 1 PCL 2 PHL 44 PHR 61 Bundles not named in the implanted electrodes of the subject description: P08S1: units on ['PCL']; implanted: ['AHL', 'AHR', 'AL', 'AR', 'ECL', 'ECR', 'PHL', 'PHR'] P08S2: units on ['PCL']; implanted: ['AHL', 'AHR', 'AL', 'AR', 'ECL', 'ECR', 'PHL', 'PHR'] P11S1: units on ['HL', 'HR']; implanted: ['AHL', 'AHR', 'AL', 'AR']
# Region of each unit. The claims file (neurons-recorded) names two readings: (a) the anatomical
# location of the microwire tips (Figure 1d legend), taken here as the atlas label of the bundle's
# tip contact (contact 1; the microwires protrude 4 mm from the tip, Methods 2.3), and (b) the
# electrode name (AH and PH hippocampus, EC entorhinal cortex, A amygdala, Table 1 legend). The
# `location` column of the unit's own electrode row is shown as well.
import h5py, remfile
def tip_locations(s):
f = h5py.File(remfile.File(m._s3_url(s["asset_id"])), "r")
el = f["general/extracellular_ephys/electrodes"]
lab = [m._text(x) for x in el["label"][:]]
loc = [m._text(x) for x in el["location"][:]]
f.close()
contacts = {}
for l, a in zip(lab, loc):
contacts.setdefault(l.rstrip("0123456789"), {})[int(l[len(l.rstrip("0123456789")):])] = a.split(",")[0]
return {b: [c[i] for i in sorted(c)] for b, c in contacts.items()} # contact 1 (tip) first
for s in S:
s["contacts"] = tip_locations(s)
s["tip"] = {b: c[0] for b, c in s["contacts"].items()}
ATLAS = {"Hipp": "Hipp", "PhG": "Ent", "Amyg": "Amg"}
UNITS["tip_atlas"] = [S[[x["key"] for x in S].index(k)]["tip"].get(b, "?").split(",")[0] for k, b in zip(UNITS.key, UNITS.bundle)]
UNITS["region_tip"] = [ATLAS.get(a, r) for a, r in zip(UNITS.tip_atlas, UNITS.region)]
UNITS["region_location"] = [ATLAS.get(l.split(",")[0], "other") for l in UNITS.location]
# Reading (a) leaves open where the 133 units whose tip contact has no MTL label go. Two choices:
# (a1) the region of the electrode name (column region_tip), (a2) the label of the contact nearest
# the tip that has an MTL label, and the electrode name only where no contact has one (region_tip2).
def nearest_mtl(key, bundle):
labs = S[[x["key"] for x in S].index(key)]["contacts"].get(bundle, [])
return next((ATLAS[a] for a in labs if a in ATLAS), None)
UNITS["region_tip2"] = [ATLAS[a] if a in ATLAS else (nearest_mtl(k, b) or r)
for a, r, k, b in zip(UNITS.tip_atlas, UNITS.region, UNITS.key, UNITS.bundle)]
for s in S:
u = UNITS[UNITS.key == s["key"]]
s["units"]["region_tip"] = u.region_tip.values
s["units"]["region_tip2"] = u.region_tip2.values
print("Region by electrode name (reading b) against the tip contact's atlas label:")
print(pd.crosstab(UNITS.region, UNITS.tip_atlas, margins=True).to_string())
UNITS["region_tip_strict"] = [ATLAS.get(a, "other") for a in UNITS.tip_atlas]
fallback = ~UNITS.tip_atlas.isin(list(ATLAS))
print(f"\nReading (a): {fallback.sum()} units sit on bundles whose tip contact has no hippocampal, parahippocampal or")
print("amygdala label; Figure 1d puts every recording site in one of the three regions, so these take the")
print("region of their electrode name. Their tip labels:", UNITS[fallback].tip_atlas.value_counts().to_dict())
fb = UNITS[fallback].groupby(["key", "bundle"]).agg(units=("label", "size"), electrode_name=("region", "first"),
a2=("region_tip2", "first")).reset_index()
fb["contacts 1 to 4"] = [", ".join(S[[x["key"] for x in S].index(k)]["contacts"][b][:4]) for k, b in zip(fb.key, fb.bundle)]
print("\nThe 133 units by session and bundle, with the atlas labels of the bundle's contacts 1 to 4 (1 is the tip),")
print("the region by electrode name (a1) and by the nearest contact with an MTL label (a2):")
print(fb.to_string(index=False))
print("\nCounts under each reading (paper: Hipp 454, Ent 209, Amg 146):")
print(pd.DataFrame({"(a2) tip, nearest MTL contact": UNITS.region_tip2.value_counts(),
"(a1) tip, electrode-name fallback": UNITS.region_tip.value_counts(),
"(a) tip label only": UNITS.region_tip_strict.value_counts(),
"(b) electrode name": UNITS.region.value_counts(),
"location column": UNITS.region_location.value_counts()}).fillna(0).astype(int).to_string())
ex = UNITS[UNITS.key == "P01S1"][["label", "location"]].head(8)
print("\nExample, P01S1: the location of a unit is the atlas label of the macro contact with its wire's number")
print(ex.to_string(index=False))
Region by electrode name (reading b) against the tip contact's atlas label:
tip_atlas Amyg FuG Hipp PhG unspecific All
region
Amg 102 0 124 0 0 226
Ent 44 0 5 168 0 217
Hipp 0 13 228 0 120 361
other:HL 0 0 2 0 0 2
other:HR 0 0 1 0 0 1
other:PCL 0 0 2 0 0 2
All 146 13 362 168 120 809
Reading (a): 133 units sit on bundles whose tip contact has no hippocampal, parahippocampal or
amygdala label; Figure 1d puts every recording site in one of the three regions, so these take the
region of their electrode name. Their tip labels: {'unspecific': 120, 'FuG': 13}
The 133 units by session and bundle, with the atlas labels of the bundle's contacts 1 to 4 (1 is the tip),
the region by electrode name (a1) and by the nearest contact with an MTL label (a2):
key bundle units electrode_name a2 contacts 1 to 4
P05S1 PHL 1 Hipp Hipp unspecific, unspecific, unspecific, MTG
P06S1 PHR 5 Hipp Hipp unspecific, unspecific, unspecific, MTG
P07S1 AHR 41 Hipp Hipp unspecific, unspecific, BG, INS
P08S1 AHL 27 Hipp Hipp unspecific, Hipp, FuG, ITG
P08S1 PHR 15 Hipp Ent unspecific, PhG, PhG, FuG
P08S2 AHL 18 Hipp Hipp unspecific, Hipp, FuG, ITG
P08S2 PHR 13 Hipp Ent unspecific, PhG, PhG, FuG
P13S1 AHR 13 Hipp Hipp FuG, FuG, unspecific, MTG
Counts under each reading (paper: Hipp 454, Ent 209, Amg 146):
(a2) tip, nearest MTL contact (a1) tip, electrode-name fallback (a) tip label only (b) electrode name location column
Amg 146 146 146 226 46
Ent 196 168 168 217 55
Hipp 467 495 362 361 94
other 0 0 133 0 614
other:HL 0 0 0 2 0
other:HR 0 0 0 1 0
other:PCL 0 0 0 2 0
Example, P01S1: the location of a unit is the atlas label of the macro contact with its wire's number
label location
AHL2 Hipp, Left Hippocampus rHipp, rostral hippocampus
AHL2 Hipp, Left Hippocampus rHipp, rostral hippocampus
AHL2 Hipp, Left Hippocampus rHipp, rostral hippocampus
AHL2 Hipp, Left Hippocampus rHipp, rostral hippocampus
AHL2 Hipp, Left Hippocampus rHipp, rostral hippocampus
AHL2 Hipp, Left Hippocampus rHipp, rostral hippocampus
AHL2 Hipp, Left Hippocampus rHipp, rostral hippocampus
AHL2 Hipp, Left Hippocampus rHipp, rostral hippocampus
The electrodes table has one row per macro contact (AHL1 to AHL8 and so on, MNI coordinates about
5 mm apart, as the contacts of the depth electrodes are), and each unit points to the row whose
number is its microwire's number. The microwires of a bundle all protrude from the tip of the
electrode (Methods 2.3), so the location of a unit is the atlas label of a macro contact, not of
its wire: units on wires 5 to 8 get lateral temporal labels (MTG, STG, ITG). The location column
does not give the units' regions.
The claims file names two readings of the region (neurons-recorded). Reading (a), the anatomical location of the microwire tips as Figure 1d shows them, is taken here as the atlas label of the tip contact of each bundle: Hipp is hippocampus, PhG (parahippocampal gyrus, which holds the entorhinal cortex) entorhinal cortex and Amyg amygdala. The files hold no location for the wires themselves, so this is the nearest the dataset comes to the tip location. 133 units are on bundles whose tip contact has another label (unspecific or FuG, all on hippocampal electrodes); Figure 1d colours every site as one of the three regions, so these keep the region of their electrode name, and the count without them is printed as well.
Where those 133 units go is not fixed by the paper, and it changes the region counts, so both choices are run. (a1) gives them the region of their electrode name. (a2) gives them the label of the contact nearest the tip that has an MTL label, and the electrode name only where no contact of the bundle has one. The crosstab above shows the difference: the 28 units on PHR of participant 8 (contacts 2 and 3 in PhG) go to entorhinal cortex under (a2), the 45 units on AHL of participant 8 (contact 2 in Hipp) stay in hippocampus, and the remaining 60 units (P05 PHL, P06 PHR, P07 AHR, P13 AHR) have no MTL contact and keep their electrode name under both. (a2) gives 467, 196 and 146. It is adopted for reading (a) because it uses the atlas labels of the bundle's own contacts before the electrode name. The grounding is weak: it moves only the 28 PHR units, and it places them by contacts 2 and 3, which lie further from the protruding wire tips than contact 1. It was examined after the verifier's first report pointed at these units, and it is not the paper's assignment (which the files do not hold); the outcomes that depend on it say so. Reading (b) is the electrode name (AH and PH hippocampus, EC entorhinal cortex, A amygdala, Table 1 legend); 5 units are on bundles it does not cover (PCL in participant 8, HL and HR in participant 11, whose subject description lists only AHL, AHR, AL and AR). Reading (a2) gives 467, 196 and 146 units, (a1) 495, 168 and 146, reading (b) 361, 217 and 226, against the paper's 454, 209 and 146. Every per-region result below is given under all three.
# Response time. The files' response timestamps are the conversion's `Response time` + trial start
# + 1.7 s (usz-neuro-conversion, write_behavior). Two readings of what they hold:
rt_stored = TRIALS.response_ts - TRIALS.start - m.RT_OFFSET_S # stored value, read as RT from test onset
rt_as_time = TRIALS.response_ts - (TRIALS.start + m.PERIODS["test"][0]) # timestamp read as the press time
print("Timestamp read as the time of the press, minus test onset (s):")
print(rt_as_time.describe().round(3).to_string())
print(" responses before test onset:", int((rt_as_time < 0).sum()), "of", len(rt_as_time))
print("\nStored Response time (timestamp - trial start - 1.7 s), read as RT from test onset (s):")
print(rt_stored.describe().round(3).to_string())
print(" under 0.2 s:", int((rt_stored < 0.2).sum()), "; over 2 s (after the test array):", int((rt_stored > 2).sum()))
Timestamp read as the time of the press, minus test onset (s): count 3038.000 mean 0.112 std 0.352 min -0.486 25% -0.120 50% 0.026 75% 0.263 max 2.246 responses before test onset: 1385 of 3038 Stored Response time (timestamp - trial start - 1.7 s), read as RT from test onset (s): count 3038.000 mean 1.012 std 0.352 min 0.414 25% 0.780 50% 0.926 75% 1.163 max 3.146 under 0.2 s: 0 ; over 2 s (after the test array): 62
Read as the time of the button press, the response timestamps put 1385 responses (of 3038) before
the test array appears, which cannot be. The conversion code adds the trial start and 1.7 s (its
offset for the event times) to the stored Response time. Taking the stored value as the response
time from test onset gives 0.41 to 3.15 s (median 0.93 s), a plausible range, so the response
timestamps are taken to be shifted by 0.9 s, and RT is computed as timestamp minus trial start
minus 1.7 s. This reading is consistent with the data, not shown by it; it changes only absolute
RTs, not their differences between set sizes. The response description gives no unit of time.
# Within-trial timing: population response around the memory array (z-scored per unit, 50 ms bins).
BIN = 0.05
edges = np.arange(0, m.TRIAL_S, BIN)
psth = []
for s in S:
wc = m.window_counts(s, edges, BIN).mean(1)
sd = wc.std(1, keepdims=True); sd[sd == 0] = 1
psth.append((wc - wc[:, :int(0.9 / BIN)].mean(1, keepdims=True)) / sd)
psth = np.concatenate(psth)
pop = psth.mean(0)
centres = edges + BIN / 2
mem = (centres > 0.9) & (centres < 1.7)
print(f"Population z-score, mean over units: fixation {pop[centres < 0.9].mean():.3f}; peak after memory onset "
f"{pop[mem].max():.3f} at {centres[mem][pop[mem].argmax()] - 0.9:.3f} s after onset; "
f"peak after test onset {pop[centres > 2.6].max():.3f} at {centres[centres > 2.6][pop[centres > 2.6].argmax()] - 2.6:.3f} s")
print("Range of the plotted curve:", pop.min().round(3), pop.max().round(3))
fig, ax = plt.subplots(figsize=(7, 2.5))
ax.plot(centres, pop, color="k")
for x in [0.9, 1.7, 2.6]:
ax.axvline(x, color="grey", lw=0.8)
ax.set_xlabel("time from fixation onset (s)"); ax.set_ylabel("mean z-score")
plt.show()
Population z-score, mean over units: fixation -0.000; peak after memory onset 0.077 at 0.175 s after onset; peak after test onset -0.037 at 0.025 s Range of the plotted curve: -0.252 0.141
# Alignment of spikes with the trial labels: set-size information (sum over units of the squared
# Welch t of low against high workload, all trials, encoding to test) when the spikes of trial
# window i are paired with the label of trial i + k.
K = np.arange(-12, 13)
scan = {}
for s in S:
c = m.period_counts(s, {"em": (0.9, 2.6), "t": (2.6, 4.6)})
ss = s["trials"].set_size.values
n = len(ss)
tot = np.zeros(len(K))
for per in ["em", "t"]:
x = c[per].astype(float)
for j, k in enumerate(K):
xs, lab = (x[:, :n - k], ss[k:]) if k >= 0 else (x[:, -k:], ss[:n + k])
tot[j] += (np.nan_to_num(m.welch_t(xs[:, lab <= 2].T, xs[:, lab >= 4].T)) ** 2).sum()
scan[s["key"]] = tot / len(s["units"]) / 2
scan = pd.DataFrame(scan, index=K).T
pooled = scan.mul([len(s["units"]) for s in S], axis=0).sum() / len(UNITS)
print("Mean squared t per unit and period, pooled over sessions, by shift k (spikes of window i with the label of trial i + k):")
print(pooled.round(3).to_string())
print("\nPer session, the shift with the largest value, and the values at k = 0 and k = 3:")
print(pd.DataFrame({"best k": scan.idxmax(axis=1), "at best": scan.max(axis=1).round(2), "k=0": scan[0].round(2),
"k=3": scan[3].round(2), "median over k": scan.median(axis=1).round(2)}).to_string())
Mean squared t per unit and period, pooled over sessions, by shift k (spikes of window i with the label of trial i + k):
-12 0.986
-11 0.988
-10 0.940
-9 0.961
-8 0.992
-7 0.964
-6 0.992
-5 1.007
-4 0.953
-3 0.943
-2 0.907
-1 0.952
0 0.970
1 0.969
2 1.039
3 1.158
4 1.045
5 0.962
6 0.954
7 0.947
8 0.979
9 0.955
10 0.918
11 0.951
12 0.972
Per session, the shift with the largest value, and the values at k = 0 and k = 3:
best k at best k=0 k=3 median over k
P01S1 6 1.32 0.98 0.87 0.96
P02S1 4 1.33 0.75 0.75 0.99
P03S1 3 1.64 1.16 1.64 0.98
P03S2 4 1.71 1.05 1.07 0.94
P04S1 2 1.34 0.93 1.25 1.02
P05S1 -4 1.39 0.81 1.24 0.95
P06S1 3 1.89 1.13 1.89 1.03
P07S1 -5 1.37 0.73 0.86 0.89
P08S1 0 1.14 1.14 0.98 0.98
P08S2 4 1.21 0.92 0.77 0.93
P09S1 3 1.43 0.78 1.43 0.87
P10S1 -4 1.28 0.91 0.94 0.92
P11S1 -9 1.60 0.61 1.12 0.88
P12S1 -5 1.28 1.21 1.06 1.00
P12S2 9 1.41 0.94 1.07 1.01
P13S1 9 1.28 0.86 1.04 0.90
# The same check by decoding workload (all neurons, correct trials, maintenance windows) after
# moving each window's spikes k trials later; 20 runs and 20 scrambled runs per shift.
centres_mem = np.round(np.arange(0.9, 1.65, 0.1), 2) # window centres, s from memory onset
offs = m.PERIODS["encoding"][0] + centres_mem - 0.5
def workload_inputs(sessions, unit_mask=None):
feats, groups = [], []
for i, s in enumerate(sessions):
t = s["trials"]; ok = np.where(t.correct.values)[0]
wc = m.window_counts(s, offs, 1.0)
if unit_mask is not None:
wc = wc[unit_mask[i]]
if wc.shape[0] == 0:
continue
feats.append(np.moveaxis(wc[:, ok, :], 1, 0))
ss = t.set_size.values[ok]
groups.append([np.where(ss <= 2)[0], np.where(ss >= 4)[0]])
return feats, groups
shift_rows = []
for k in [-3, 0, 2, 3, 4]:
f_, g_ = workload_inputs([m.realign(s, k) for s in S])
r, n = m.decode(f_, g_, 20, 20, 20, "c", True, seed=100 + k)
shift_rows.append({"k": k, "median accuracy, windows 0.9-1.6 s": np.median(np.median(r, 0)).round(3),
"max over windows": np.median(r, 0).max().round(3), "scrambled, median": np.median(n).round(3),
"scrambled, max run": n.max().round(3)})
print(pd.DataFrame(shift_rows).to_string(index=False))
k median accuracy, windows 0.9-1.6 s max over windows scrambled, median scrambled, max run -3 0.492 0.533 0.492 0.667 0 0.458 0.471 0.496 0.642 2 0.498 0.512 0.488 0.633 3 0.598 0.608 0.504 0.650 4 0.535 0.546 0.508 0.625
# Alignment scan with a null: the statistic of the scan above, per session, against 200 random
# orders of the session's set-size labels (which carry no alignment). For each session: the p of the
# value at k = 0 and at k = 3, and the p of the largest value over all shifts against the largest
# value of each shuffled scan (corrected for the 25 shifts). Then the scan pooled over the sessions
# with 192 trials and over those with fewer.
def scan_curve(c, ss):
n = len(ss); out = np.zeros(len(K))
for x in c:
for j, k in enumerate(K):
xs, lab = (x[:, :n - k], ss[k:]) if k >= 0 else (x[:, -k:], ss[:n + k])
out[j] += (np.nan_to_num(m.welch_t(xs[:, lab <= 2].T, xs[:, lab >= 4].T)) ** 2).sum()
return out / x.shape[0] / 2
rng = np.random.default_rng(13)
started = time.time()
rows, obs_c, null_c = [], {}, {}
for s in S:
cc = m.period_counts(s, {"em": (0.9, 2.6), "t": (2.6, 4.6)})
c = [cc["em"].astype(float), cc["t"].astype(float)]
ss = s["trials"].set_size.values
o = scan_curve(c, ss)
nl = np.array([scan_curve(c, rng.permutation(ss)) for _ in range(200)])
obs_c[s["key"]], null_c[s["key"]] = o * len(s["units"]), nl * len(s["units"])
i0, i3 = list(K).index(0), list(K).index(3)
rows.append({"session": s["key"], "trials": len(ss), "units": len(s["units"]), "k=0": o[i0],
"p k=0": (np.sum(nl[:, i0] >= o[i0]) + 1) / 201, "k=3": o[i3], "p k=3": (np.sum(nl[:, i3] >= o[i3]) + 1) / 201,
"best k": K[o.argmax()], "max": o.max(), "p max (corrected)": (np.sum(nl.max(1) >= o.max()) + 1) / 201})
al = pd.DataFrame(rows)
print(al.round(3).to_string(index=False))
print(f"Sessions with p < .05: at k = 0 {int((al['p k=0'] < .05).sum())}, at k = 3 {int((al['p k=3'] < .05).sum())}, "
f"largest over shifts (corrected) {int((al['p max (corrected)'] < .05).sum())} of {len(al)}")
for name, keys in [("192 trials", al.session[al.trials == 192]), ("fewer than 192", al.session[al.trials < 192])]:
o = sum(obs_c[k] for k in keys); nl = sum(null_c[k] for k in keys)
u = al.units[al.session.isin(keys)].sum()
p = [(np.sum(nl[:, j] >= o[j]) + 1) / 201 for j in range(len(K))]
print(f"\nPooled over the {len(keys)} sessions with {name} ({u} units): value per unit by shift, and p against the shuffles")
print(pd.DataFrame({"value": (o / u).round(3), "p": np.round(p, 3)}, index=K).T.to_string())
print(f" largest at k = {K[o.argmax()]}, corrected p = {(np.sum(nl.max(1) >= o.max()) + 1) / 201:.3f}")
print("Columns of the trial table:", list(S[0]["trials"].columns[:5]), "- the ids run 0 to n - 1 in every session:",
all((s["trials"].trial_id.values == np.arange(len(s["trials"]))).all() for s in S))
print(f"({time.time() - started:.0f} s)")
session trials units k=0 p k=0 k=3 p k=3 best k max p max (corrected)
P01S1 189 33 0.977 0.557 0.868 0.751 6 1.318 0.881
P02S1 192 36 0.751 0.900 0.754 0.950 4 1.330 0.687
P03S1 192 52 1.164 0.204 1.639 0.005 3 1.639 0.070
P03S2 191 54 1.052 0.363 1.072 0.318 4 1.715 0.030
P04S1 191 85 0.927 0.706 1.245 0.050 2 1.338 0.313
P05S1 192 29 0.806 0.866 1.242 0.169 -4 1.391 0.731
P06S1 192 96 1.132 0.169 1.893 0.005 3 1.893 0.005
P07S1 189 48 0.727 0.905 0.860 0.687 -5 1.369 0.806
P08S1 190 87 1.144 0.219 0.978 0.587 0 1.144 0.945
P08S2 191 63 0.920 0.687 0.767 0.965 4 1.208 0.896
P09S1 174 23 0.780 0.821 1.429 0.065 3 1.429 0.806
P10S1 188 62 0.907 0.726 0.938 0.632 -4 1.278 0.657
P11S1 192 11 0.611 0.940 1.124 0.333 -9 1.597 0.706
P12S1 191 29 1.211 0.124 1.056 0.398 -5 1.280 0.935
P12S2 192 28 0.936 0.647 1.067 0.328 9 1.409 0.786
P13S1 192 73 0.858 0.801 1.043 0.378 9 1.280 0.692
Sessions with p < .05: at k = 0 0, at k = 3 3, largest over shifts (corrected) 2 of 16
Pooled over the 7 sessions with 192 trials (325 units): value per unit by shift, and p against the shuffles
-12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12
value 1.064 1.086 0.914 0.996 1.011 0.960 1.032 0.957 0.964 0.997 0.901 0.981 0.970 0.920 1.086 1.380 1.035 1.016 0.935 0.914 0.974 1.008 0.895 1.110 1.031
p 0.224 0.109 0.925 0.502 0.473 0.811 0.299 0.741 0.781 0.552 0.950 0.711 0.771 0.896 0.144 0.005 0.303 0.443 0.900 0.955 0.756 0.507 0.930 0.095 0.328
largest at k = 3, corrected p = 0.005
Pooled over the 9 sessions with fewer than 192 (484 units): value per unit by shift, and p against the shuffles
-12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12
value 0.933 0.922 0.957 0.938 0.980 0.967 0.965 1.040 0.946 0.906 0.911 0.932 0.970 1.002 1.008 1.010 1.051 0.925 0.967 0.970 0.982 0.919 0.933 0.844 0.933
p 0.935 0.950 0.846 0.910 0.711 0.801 0.756 0.313 0.891 0.970 0.985 0.900 0.771 0.577 0.488 0.418 0.303 0.955 0.781 0.746 0.672 0.965 0.935 1.000 0.930
largest at k = 4, corrected p = 1.000
Columns of the trial table: ['trial_id', 'start', 'stop', 'set_size', 'solution_same'] - the ids run 0 to n - 1 in every session: True
(18 s)
Within trials the timing is consistent with the epochs: the population rate rises after the memory array (peak 0.175 s after onset, printed above), and many neurons differ between periods (fig2c and fig3b below), so the spikes share the clock of the epochs within each trial. The pairing of spikes with the trial labels does not hold as stored. Workload information in the spikes peaks when the spikes stored in trial window i are paired with the label of trial i + 3, not of trial i (pooled value 1.16 at k = 3 against 0.91 to 1.05 at other shifts), and decoding workload from all neurons gives 0.46 at k = 0, below the scrambled runs (0.50), against 0.60 at k = 3 (printed above).
Against 200 random orders of each session's set-size labels, no session has more set-size information at k = 0 than the shuffles (p of 0.12 to 0.94). The largest value over all shifts is significant after correction for the 25 shifts in 2 of 16 sessions (P06S1 at k = 3, p = 0.005; P03S2 at k = 4, p = 0.03), and 3 sessions have p < .05 at k = 3 (P03S1, P06S1, P04S1). Pooled over the 7 sessions with 192 trials, the scan peaks at k = 3 (1.38 per unit, corrected p = 0.005) and no other shift has p < .05; pooled over the 9 sessions with fewer trials it has no peak (largest at k = 4, corrected p = 1.0). The pooled peak does not hold session by session. Of the 7 complete sessions, only P03S1 and P06S1 have p < .05 at k = 3; P02S1 has p = 0.95 there, and P05S1, P11S1, P12S2 and P13S1 have their largest value at k = -4, -9, 9 and 9. The pooled value is weighted toward P06S1 (96 units) and P03S1 (52 units). Two sessions with fewer trials also point to k = 3 (P04S1, p = 0.050 at k = 3; P09S1, largest at k = 3, p = 0.065 at k = 3). So a shift by three trials is supported in two sessions and in the pooled complete sessions, not in each complete session, and the offset in the other sessions could not be established. The trial table holds no trial numbering from the source (its ids run 0 to n - 1), so where the missing trials were cannot be read from the files. Every claim below is tested on the data as deposited (k = 0); claims that use the trial labels (workload, correct and incorrect, set size) also show, as a variant, the result with every session moved by k = 3. That shift was chosen from these data and is applied to all 16 sessions, where the scan supports it in few of them, so the variant is an indication, not a test.
behavior-accuracy-by-set-size¶
Performance was almost perfect at set sizes 1 (98% correct) and 2 (99%) and fell at 4 (88%) and 6 (73%), a significant effect of set size in a permuted repeated-measures ANOVA over 16 sessions; average performance was 89%.
rng = np.random.default_rng(11)
SS = [1, 2, 4, 6]
def rm_anova_F(Y):
"""One within factor (columns), subjects in rows."""
n, k = Y.shape
gm = Y.mean()
ss_cond = n * ((Y.mean(0) - gm) ** 2).sum()
ss_subj = k * ((Y.mean(1) - gm) ** 2).sum()
ss_err = ((Y - gm) ** 2).sum() - ss_cond - ss_subj
return (ss_cond / (k - 1)) / (ss_err / ((n - 1) * (k - 1)))
def per_session(values_by_session, labels_by_session):
return np.array([[v[l == c].mean() for c in SS] for v, l in zip(values_by_session, labels_by_session)])
vals = [s["trials"].correct.values.astype(float) for s in S]
labs = [s["trials"].set_size.values for s in S]
Y = per_session(vals, labs) * 100
pooled_pc = TRIALS.groupby("set_size").correct.mean() * 100
print("Percent correct per session and set size:")
print(pd.DataFrame(Y, index=[s["key"] for s in S], columns=SS).round(1).to_string())
print("\n(a) trials pooled over sessions:", pooled_pc.round(2).to_dict(), "; all set sizes:", round(TRIALS.correct.mean() * 100, 2))
print("(b) mean of session percentages:", dict(zip(SS, Y.mean(0).round(2))), "; all set sizes:",
round(np.mean([v.mean() for v in vals]) * 100, 2))
F = rm_anova_F(Y)
N_PERM = 4999
null = np.empty(N_PERM)
for i in range(N_PERM):
null[i] = rm_anova_F(per_session(vals, [rng.permutation(l) for l in labs]))
print(f"\nRepeated-measures ANOVA: F(3, 45) = {F:.2f}; permutation p (set-size labels shuffled within session, "
f"{N_PERM} permutations) = {(np.sum(null >= F) + 1) / (N_PERM + 1):.4f}; largest null F {null.max():.2f}")
print(f"Parametric p for comparison: {stats.f.sf(F, 3, 45):.2e}")
null_rows = np.array([rm_anova_F(np.array([rng.permutation(r) for r in Y])) for _ in range(N_PERM)])
print(f"Variant, the four session values permuted within each session: p = {(np.sum(null_rows >= F) + 1) / (N_PERM + 1):.4f}")
print("Order 1 and 2 above 4 above 6, reading (a):", bool(min(pooled_pc[1], pooled_pc[2]) > pooled_pc[4] > pooled_pc[6]),
"; reading (b):", bool(min(Y.mean(0)[0], Y.mean(0)[1]) > Y.mean(0)[2] > Y.mean(0)[3]))
Percent correct per session and set size:
1 2 4 6
P01S1 95.7 91.5 76.6 62.5
P02S1 93.8 93.8 75.0 68.8
P03S1 97.9 97.9 70.8 64.6
P03S2 95.8 100.0 70.8 74.5
P04S1 100.0 100.0 83.3 72.9
P05S1 95.8 100.0 77.1 72.9
P06S1 100.0 97.9 87.5 72.9
P07S1 100.0 97.9 87.5 66.7
P08S1 95.7 100.0 91.7 70.2
P08S2 100.0 97.9 78.7 75.0
P09S1 95.1 100.0 90.7 66.7
P10S1 100.0 97.9 93.6 68.1
P11S1 97.9 100.0 89.6 81.2
P12S1 100.0 100.0 89.4 83.3
P12S2 100.0 100.0 95.8 72.9
P13S1 95.8 91.7 91.7 87.5
(a) trials pooled over sessions: {1: 97.76, 2: 97.9, 4: 84.32, 6: 72.6} ; all set sizes: 88.15
(b) mean of session percentages: {1: np.float64(97.73), 2: np.float64(97.9), 4: np.float64(84.36), 6: np.float64(72.54)} ; all set sizes: 88.15
Repeated-measures ANOVA: F(3, 45) = 89.78; permutation p (set-size labels shuffled within session, 4999 permutations) = 0.0002; largest null F 10.81 Parametric p for comparison: 5.12e-19
Variant, the four session values permuted within each session: p = 0.0002 Order 1 and 2 above 4 above 6, reading (a): True ; reading (b): True
Result: reproduced. Pooled over trials, 97.76, 97.90, 84.32 and 72.60% correct at set sizes 1, 2, 4 and 6 and 88.15% overall; the mean of session percentages gives 97.73, 97.90, 84.36, 72.54 and 88.15%. All are within 10% of 98, 99, 88, 73 and 89%, and the order holds. The repeated-measures ANOVA gives F(3, 45) = 89.78 (paper 89.79) with permutation p = 0.0002, the smallest possible with 4999 permutations.
behavior-capacity¶
Across sessions the median working memory capacity was 3.2 by Pashler's K (range 2.2 to 5) and 3 by Cowan's K (range 1.9 to 4.5), so participants held at least 3 squares and set sizes 4 and 6 exceeded capacity.
rows = []
for s in S:
t = s["trials"]
for n_items in SS:
x = t[t.set_size == n_items]
H = (~x[~x.solution_same].answer_same).mean() # change trials answered "different"
FA = (~x[x.solution_same].answer_same).mean() # no-change trials answered "different"
rows.append({"key": s["key"], "participant": int(s["subject"]), "N": n_items, "H": H, "FA": FA,
"KP": (H - FA) * n_items / (1 - FA), "KC": (H - FA) * n_items})
K_tab = pd.DataFrame(rows)
print(K_tab.pivot(index="key", columns="N", values="KP").round(2).to_string())
sess = K_tab.groupby(["key", "participant"]).agg(Kmax=("KP", "max"), KC_max=("KC", "max")).reset_index()
sess["KC_6"] = K_tab[K_tab.N == 6].set_index("key").KC.reindex(sess.key).values
part = sess.groupby("participant")[["Kmax", "KC_max", "KC_6"]].mean()
def summ(x):
return f"median {np.median(x):.2f}, range {np.min(x):.2f} to {np.max(x):.2f}"
for name, tab in [("(a) per participant (13)", part), ("(b) per session (16)", sess)]:
print(f"\n{name}: KP (Kmax) {summ(tab.Kmax)}; KC (i) max over set sizes {summ(tab.KC_max)}; KC (ii) set size 6 {summ(tab.KC_6)}")
print("\nFalse alarm rates, range over sessions and set sizes:", K_tab.FA.min().round(3), "to", K_tab.FA.max().round(3))
N 1 2 4 6 key P01S1 0.96 1.81 2.41 1.64 P02S1 0.95 1.91 2.18 2.25 P03S1 1.00 1.92 2.50 2.21 P03S2 0.92 2.00 2.22 3.58 P04S1 1.00 2.00 2.91 3.00 P05S1 0.92 2.00 2.17 2.87 P06S1 1.00 1.92 3.43 3.88 P07S1 1.00 1.92 3.00 2.12 P08S1 0.96 2.00 3.33 2.35 P08S2 1.00 1.92 2.26 3.13 P09S1 0.95 2.00 3.58 2.37 P10S1 1.00 1.91 3.81 2.25 P11S1 0.96 2.00 4.00 5.00 P12S1 1.00 2.00 3.17 4.00 P12S2 1.00 2.00 3.67 2.75 P13S1 1.00 1.67 3.33 4.50 (a) per participant (13): KP (Kmax) median 3.23, range 2.25 to 5.00; KC (i) max over set sizes median 3.00, range 2.11 to 4.50; KC (ii) set size 6 median 2.67, range 1.50 to 4.50 (b) per session (16): KP (Kmax) median 3.46, range 2.25 to 5.00; KC (i) max over set sizes median 3.00, range 1.92 to 4.50; KC (ii) set size 6 median 2.75, range 1.50 to 4.50 False alarm rates, range over sessions and set sizes: 0.0 to 0.333
Result: reproduced. Per participant (reading a), Pashler's K has median 3.23 and range 2.25 to 5.00, and Cowan's K, taken as its maximum over set sizes (reading i), median 3.00 and range 2.11 to 4.50; all are within the tolerances of 3.2 (2.2 to 5) and 3 (1.9 to 4.5). Per session (reading b) the medians are 3.46 and 3.00, also within tolerance. Cowan's K at set size 6 (reading ii) has median 2.67 per participant, just outside 3 +- 0.3, and 2.75 per session.
table1-capacity-and-neurons¶
Table 1 gives, per participant, the number of sessions, the memory capacity Kmax (2.41, 2.25, 3.04, 3.00, 2.87, 3.88, 3.00, 3.23, 3.58, 3.81, 5.00, 3.83, 4.50 for participants 1 to 13) and the number of neurons (33, 36, 106, 85, 29, 96, 48, 150, 23, 62, 11, 57, 73).
paper_k = [2.41, 2.25, 3.04, 3.00, 2.87, 3.88, 3.00, 3.23, 3.58, 3.81, 5.00, 3.83, 4.50]
paper_n = [33, 36, 106, 85, 29, 96, 48, 150, 23, 62, 11, 57, 73]
t1 = pd.DataFrame({"sessions": sess.groupby("participant").size(),
"Kmax per session": sess.groupby("participant").Kmax.apply(lambda v: ", ".join(f"{x:.2f}" for x in v)),
"Kmax": part.Kmax.round(2), "paper Kmax": paper_k,
"neurons": UNITS.groupby("participant").size(), "paper neurons": paper_n})
t1["Kmax within 10%"] = (t1.Kmax - t1["paper Kmax"]).abs() <= 0.1 * t1["paper Kmax"] + 1e-9
t1["neurons within 10%"] = (t1.neurons - t1["paper neurons"]).abs() <= 0.1 * t1["paper neurons"]
print(t1.to_string())
print("\nKmax within tolerance:", int(t1["Kmax within 10%"].sum()), "of 13; neuron counts within tolerance:", int(t1["neurons within 10%"].sum()), "of 13")
print("Largest Kmax difference:", (t1.Kmax - t1["paper Kmax"]).abs().max().round(3))
sessions Kmax per session Kmax paper Kmax neurons paper neurons Kmax within 10% neurons within 10% participant 1 1 2.41 2.41 2.41 33 33 True True 2 1 2.25 2.25 2.25 36 36 True True 3 2 2.50, 3.58 3.04 3.04 106 106 True True 4 1 3.00 3.00 3.00 85 85 True True 5 1 2.87 2.87 2.87 29 29 True True 6 1 3.88 3.88 3.88 96 96 True True 7 1 3.00 3.00 3.00 48 48 True True 8 2 3.33, 3.13 3.23 3.23 150 150 True True 9 1 3.58 3.58 3.58 23 23 True True 10 1 3.81 3.81 3.81 62 62 True True 11 1 5.00 5.00 5.00 11 11 True True 12 2 4.00, 3.67 3.83 3.83 57 57 True True 13 1 4.50 4.50 4.50 73 73 True True Kmax within tolerance: 13 of 13; neuron counts within tolerance: 13 of 13 Largest Kmax difference: 0.0
Result: reproduced. All 13 Kmax values and all 13 neuron counts equal those of Table 1 (largest Kmax difference 0.00), including the per-session values of the three two-session participants (2.50 and 3.58, 3.33 and 3.13, 4.00 and 3.67).
behavior-rt-set-size¶
Mean response time on correct trials (2678 trials) increased with set size, by 118 ms per item, a significant effect in a permuted repeated-measures ANOVA over 16 sessions.
rng = np.random.default_rng(12)
cor = [s["trials"][s["trials"].correct] for s in S]
rts = [c.rt.values * 1000 for c in cor]
lab = [c.set_size.values for c in cor]
print("Correct trials:", sum(len(c) for c in cor))
x = np.array(SS, float)
for name, fn in [("(a) mean", np.mean), ("(b) median", np.median)]:
Yr = np.array([[fn(r[l == c]) for c in SS] for r, l in zip(rts, lab)])
F = rm_anova_F(Yr)
null = np.array([rm_anova_F(np.array([[fn(r[pl == c]) for c in SS] for r, pl in
zip(rts, [rng.permutation(l) for l in lab])])) for _ in range(4999)])
p = (np.sum(null >= F) + 1) / 5000
s_i = np.mean([np.polyfit(x, row, 1)[0] for row in Yr])
s_ii = np.polyfit(np.concatenate(lab).astype(float), np.concatenate(rts), 1)[0]
s_iii = np.polyfit(x, Yr.mean(0), 1)[0]
print(f"\n{name} RT per session, ms: means over sessions by set size {dict(zip(SS, Yr.mean(0).round(1)))}")
print(f" F(3, 45) = {F:.2f}, permutation p = {p:.4f} (4999 within-session shuffles of set size)")
print(f" slope (i) per-session fits averaged {s_i:.1f} ms/item; (ii) pooled trials {s_ii:.1f}; (iii) fit of the four means {s_iii:.1f}")
print(f" sessions with a positive slope: {sum(np.polyfit(x, row, 1)[0] > 0 for row in Yr)} of 16")
step = np.arange(1, 5.0)
print(f" variant, not a reading of the claims file (asked in questions/000575.yaml): slope on the condition step "
f"1 to 4 in place of the set size, fit of the four means {np.polyfit(step, Yr.mean(0), 1)[0]:.1f} ms/step, "
f"per-session fits averaged {np.mean([np.polyfit(step, row, 1)[0] for row in Yr]):.1f}")
print("\nThe slopes do not depend on the RT offset assumed in the data check: a constant offset leaves them unchanged.")
Correct trials: 2678
(a) mean RT per session, ms: means over sessions by set size {1: np.float64(812.4), 2: np.float64(900.8), 4: np.float64(1075.8), 6: np.float64(1187.1)}
F(3, 45) = 61.36, permutation p = 0.0002 (4999 within-session shuffles of set size)
slope (i) per-session fits averaged 75.8 ms/item; (ii) pooled trials 76.2; (iii) fit of the four means 75.8
sessions with a positive slope: 16 of 16
variant, not a reading of the claims file (asked in questions/000575.yaml): slope on the condition step 1 to 4 in place of the set size, fit of the four means 129.9 ms/step, per-session fits averaged 129.9
(b) median RT per session, ms: means over sessions by set size {1: np.float64(783.8), 2: np.float64(875.4), 4: np.float64(1005.1), 6: np.float64(1134.1)}
F(3, 45) = 48.47, permutation p = 0.0002 (4999 within-session shuffles of set size)
slope (i) per-session fits averaged 68.8 ms/item; (ii) pooled trials 76.2; (iii) fit of the four means 68.8
sessions with a positive slope: 16 of 16
variant, not a reading of the claims file (asked in questions/000575.yaml): slope on the condition step 1 to 4 in place of the set size, fit of the four means 118.1 ms/step, per-session fits averaged 118.1
The slopes do not depend on the RT offset assumed in the data check: a constant offset leaves them unchanged.
Result: partially reproduced. Response time on the 2678 correct trials rises with set size in all 16 sessions, and the ANOVA on per-session median RT gives F(3, 45) = 48.47, the paper's value (permutation p = 0.0002); on mean RT, F = 61.36. The slope is 68.8 to 76.2 ms/item under the three readings of the claims file, outside 118 +- 11.8 ms/item. The criterion has no partial clause; the outcome is partial because the rise, its test and the F value hold and the slope does not. Fitting the median RTs on the condition step (1 to 4) in place of the set size gives 118.1 ms per step, the paper's number (mean RTs, 129.9). That reading is not in the claims file; it is asked as a question for a reader (questions/000575.yaml), and the outcome is judged on the readings the file names.
neurons-recorded¶
809 neurons were recorded, 454 in hippocampus, 209 in entorhinal cortex and 146 in amygdala, 50 plus or minus 25 (range 11 to 96) per session.
per_sess = UNITS.groupby("key").size()
print("Neurons:", len(UNITS))
print("By region, reading (a2) tip location, nearest MTL contact:", UNITS.region_tip2.value_counts().to_dict())
print("By region, reading (a1) tip location, electrode-name fallback:", UNITS.region_tip.value_counts().to_dict())
print(" tip label only (133 units unassigned):", UNITS.region_tip_strict.value_counts().to_dict())
print("By region, reading (b) electrode name:", UNITS.region.value_counts().to_dict())
paper_r = {"Hipp": 454, "Ent": 209, "Amg": 146}
for col, name in [("region_tip2", "(a2)"), ("region_tip", "(a1)"), ("region", "(b)")]:
vc = UNITS[col].value_counts()
print(f"Reading {name} within 10% of the paper:", {g: bool(abs(vc.get(g, 0) - v) <= 0.1 * v) for g, v in paper_r.items()})
print(f"Per session: mean {per_sess.mean():.1f}, SD {per_sess.std(ddof=1):.1f} (ddof 0: {per_sess.std(ddof=0):.1f}), "
f"range {per_sess.min()} to {per_sess.max()}")
for col, name in [("region_tip2", "(a2) tip location, nearest MTL contact"), ("region_tip", "(a1) tip location, electrode-name fallback"),
("region", "(b) electrode name")]:
print(f"\nPer participant and region, reading {name}:")
print(pd.crosstab(UNITS.participant, UNITS[col], margins=True).to_string())
Neurons: 809
By region, reading (a2) tip location, nearest MTL contact: {'Hipp': 467, 'Ent': 196, 'Amg': 146}
By region, reading (a1) tip location, electrode-name fallback: {'Hipp': 495, 'Ent': 168, 'Amg': 146}
tip label only (133 units unassigned): {'Hipp': 362, 'Ent': 168, 'Amg': 146, 'other': 133}
By region, reading (b) electrode name: {'Hipp': 361, 'Amg': 226, 'Ent': 217, 'other:PCL': 2, 'other:HL': 2, 'other:HR': 1}
Reading (a2) within 10% of the paper: {'Hipp': True, 'Ent': True, 'Amg': True}
Reading (a1) within 10% of the paper: {'Hipp': True, 'Ent': False, 'Amg': True}
Reading (b) within 10% of the paper: {'Hipp': False, 'Ent': True, 'Amg': False}
Per session: mean 50.6, SD 25.4 (ddof 0: 24.6), range 11 to 96
Per participant and region, reading (a2) tip location, nearest MTL contact:
region_tip2 Amg Ent Hipp All
participant
1 1 0 32 33
2 30 0 6 36
3 14 19 73 106
4 8 43 34 85
5 0 0 29 29
6 10 60 26 96
7 0 0 48 48
8 0 28 122 150
9 5 0 18 23
10 23 22 17 62
11 8 0 3 11
12 16 0 41 57
13 31 24 18 73
All 146 196 467 809
Per participant and region, reading (a1) tip location, electrode-name fallback:
region_tip Amg Ent Hipp All
participant
1 1 0 32 33
2 30 0 6 36
3 14 19 73 106
4 8 43 34 85
5 0 0 29 29
6 10 60 26 96
7 0 0 48 48
8 0 0 150 150
9 5 0 18 23
10 23 22 17 62
11 8 0 3 11
12 16 0 41 57
13 31 24 18 73
All 146 168 495 809
Per participant and region, reading (b) electrode name:
region Amg Ent Hipp other:HL other:HR other:PCL All
participant
1 4 5 24 0 0 0 33
2 0 30 6 0 0 0 36
3 60 19 27 0 0 0 106
4 8 43 34 0 0 0 85
5 0 0 29 0 0 0 29
6 10 60 26 0 0 0 96
7 0 0 48 0 0 0 48
8 75 0 73 0 0 2 150
9 5 0 18 0 0 0 23
10 23 22 17 0 0 0 62
11 8 0 0 2 1 0 11
12 2 14 41 0 0 0 57
13 31 24 18 0 0 0 73
All 226 217 361 2 1 2 809
Result: reproduced (resting on an assumption). The dandiset holds 809 neurons, 50.6 +- 25.4 per session (range 11 to 96), all within tolerance. By the tip location (reading a) with the 133 units whose tip contact has no MTL label placed by the nearest contact that has one (a2), there are 467 neurons in hippocampus, 196 in entorhinal cortex and 146 in amygdala, all within tolerance of 454, 209 and 146. Placing those units by their electrode name (a1) gives 495, 168 and 146, the entorhinal count outside tolerance, and the electrode name alone (reading b) gives 361, 217 and 226, only the entorhinal count within tolerance. The outcome rests on (a2), adopted because it uses the bundle's own atlas labels before the electrode name, though it places the 28 PHR units of participant 8 by contacts further from the wire tips than contact 1, and examined after the verifier pointed at these units; the paper's own assignment is not in the files.
spike-sorting-quality¶
Over all neurons, 0.79 plus or minus 1.54% of inter-spike intervals were shorter than 3 ms, most neurons were below 0.5%, and the peak signal-to-noise ratio of the mean waveform was 3.81 plus or minus 0.82.
isi_pct = []
for s in S:
st = s["trials"].start.values
for t in s["spike_times"]:
tr = np.searchsorted(st, t, "right") - 1
d = np.diff(t)[np.diff(tr) == 0] # intervals within one trial window only
isi_pct.append(100 * np.mean(d < 0.003) if len(d) else np.nan)
isi_pct = np.array(isi_pct)
print(f"ISI < 3 ms, reading (ii) spikes within trials: {np.nanmean(isi_pct):.2f} +- {np.nanstd(isi_pct, ddof=1):.2f} %, "
f"median {np.nanmedian(isi_pct):.2f}; units below 0.5%: {np.sum(isi_pct < 0.5)} of {np.sum(~np.isnan(isi_pct))} "
f"({np.mean(isi_pct[~np.isnan(isi_pct)] < 0.5):.1%}); units without a within-trial interval: {np.isnan(isi_pct).sum()}; max {np.nanmax(isi_pct):.2f}%")
print("Reading (i), the whole recording, cannot be computed: the files hold spikes only inside trial windows (data check).")
peak = np.abs(wm).max(1)
snr_a = peak / wsd.mean(1)
snr_peak_sd = peak / wsd[np.arange(len(wsd)), np.abs(wm).argmax(1)]
print(f"\nSNR (a), peak |mean waveform| / mean over samples of the waveform SD: {snr_a.mean():.2f} +- {snr_a.std(ddof=1):.2f} "
f"(median {np.median(snr_a):.2f}, range {snr_a.min():.2f} to {snr_a.max():.2f})")
print(f"Variant, the SD at the peak sample: {snr_peak_sd.mean():.2f} +- {snr_peak_sd.std(ddof=1):.2f}")
snr_first = peak / wsd[:, :5].mean(1)
print(f"Variant, the SD of the first 5 samples (before the spike): {snr_first.mean():.2f} +- {snr_first.std(ddof=1):.2f}")
print("Reading (b): the Units table has no noise level from the spike sorter (columns:", list(S[0]["units"].columns[:3]),
"+ waveform_mean, waveform_sd, offset, obs_intervals).")
ISI < 3 ms, reading (ii) spikes within trials: 0.79 +- 1.77 %, median 0.31; units below 0.5%: 506 of 809 (62.5%); units without a within-trial interval: 0; max 26.46% Reading (i), the whole recording, cannot be computed: the files hold spikes only inside trial windows (data check). SNR (a), peak |mean waveform| / mean over samples of the waveform SD: 5.93 +- 2.31 (median 5.31, range 2.42 to 23.90) Variant, the SD at the peak sample: 5.96 +- 1.43 Variant, the SD of the first 5 samples (before the spike): 5.92 +- 2.78 Reading (b): the Units table has no noise level from the spike sorter (columns: ['unit_id', 'label', 'location'] + waveform_mean, waveform_sd, offset, obs_intervals).
Result: not reproduced. Over spikes within trials (reading ii; reading i cannot be computed, as the files hold no spikes between trials) 0.79 +- 1.77% of intervals are shorter than 3 ms: the mean equals the paper's, the SD is outside 1.54 +- 0.154. 62.5% of units are below 0.5%, a majority. The waveform SNR is 5.93 +- 2.31 under reading (a) and 5.96 +- 1.43 with the SD at the peak sample, against 3.81 +- 0.82; reading (b) cannot be computed, as no noise level is stored. Neither the full nor the partial condition of the criterion is met.
firing-rates-by-period¶
The median (IQR) firing rate over all neurons was 1.09 (0.52 to 2.43) Hz over the trial, 1.12 (0.54 to 2.56) Hz during fixation, 1.10 (0.53 to 2.50) Hz during encoding, 1.09 (0.54 to 2.51) Hz during maintenance and 1.02 (0.49 to 2.45) Hz during test.
paper = {"trial": (1.09, 0.52, 2.43), "fixation": (1.12, 0.54, 2.56), "encoding": (1.10, 0.53, 2.50),
"maintenance": (1.09, 0.54, 2.51), "test": (1.02, 0.49, 2.45)}
for reading, use_correct in [("(a) correct trials", True), ("(b) all trials", False)]:
rows = []
for per in paper:
r = []
for s in S:
ok = s["trials"].correct.values if use_correct else np.ones(len(s["trials"]), bool)
if per == "trial":
c = sum(s["counts"][p] for p in m.PERIODS) / m.TRIAL_S
else:
a, b = m.PERIODS[per]
c = s["counts"][per] / (b - a)
r.append(c[:, ok].mean(1))
r = np.concatenate(r)
q = np.percentile(r, [50, 25, 75])
within = all(abs(v - p) <= 0.1 * p + 1e-9 for v, p in zip(q, paper[per]))
rows.append({"period": per, "median": q[0], "Q1": q[1], "Q3": q[2], "paper": paper[per], "all within 10%": within})
print(reading, "; trial = fixation onset to test offset, reading (i)")
print(pd.DataFrame(rows).round(3).to_string(index=False), "\n")
print("Reading (ii), warning signal to the end of the inter-trial interval, cannot be computed: no spikes are stored between trials.")
(a) correct trials ; trial = fixation onset to test offset, reading (i)
period median Q1 Q3 paper all within 10%
trial 1.079 0.531 2.450 (1.09, 0.52, 2.43) True
fixation 1.111 0.536 2.524 (1.12, 0.54, 2.56) True
encoding 1.089 0.511 2.516 (1.1, 0.53, 2.5) True
maintenance 1.091 0.537 2.508 (1.09, 0.54, 2.51) True
test 1.033 0.509 2.338 (1.02, 0.49, 2.45) True
(b) all trials ; trial = fixation onset to test offset, reading (i)
period median Q1 Q3 paper all within 10%
trial 1.088 0.522 2.427 (1.09, 0.52, 2.43) True
fixation 1.117 0.532 2.516 (1.12, 0.54, 2.56) True
encoding 1.086 0.524 2.520 (1.1, 0.53, 2.5) True
maintenance 1.099 0.541 2.482 (1.09, 0.54, 2.51) True
test 1.044 0.508 2.374 (1.02, 0.49, 2.45) True
Reading (ii), warning signal to the end of the inter-trial interval, cannot be computed: no spikes are stored between trials.
Result: reproduced. On correct trials (reading a) the medians are 1.079 Hz over the trial, 1.111 in fixation, 1.089 in encoding, 1.091 in maintenance and 1.033 in test, and every median and quartile is within 10% of the paper's. All trials (reading b) also match on every value. The trial is fixation onset to test offset (reading i); reading (ii) cannot be computed.
fig2c-maintenance-neurons¶
Maintenance neurons, firing more during maintenance than during fixation, were found in all three regions, 55 of 454 in hippocampus, 21 of 209 in entorhinal cortex and 16 of 146 in amygdala, each a significant proportion.
N_PERM, N_REP = 1000, 499
def rates(s, per, ok):
"""Rates (trials, units) in a period; "encmaint" is encoding and maintenance joined (1.7 s)."""
if per == "encmaint":
return ((s["counts"]["encoding"] + s["counts"]["maintenance"])[:, ok] / 1.7).T.astype(float)
a, b = m.PERIODS[per]
return (s["counts"][per][:, ok] / (b - a)).T.astype(float)
def select(per1, per2, paired, seed, trials="correct", region_col="region", n_rep=N_REP, sessions=None):
"""Units with a higher mean in per1 than per2 at p < .05; null counts per region from n_rep
within-trial shuffles of the two period labels."""
rng = np.random.default_rng(seed)
out, null_counts = [], []
reps = [[] for _ in range(n_rep)]
for s in (sessions or S):
ok = s["trials"].correct.values if trials == "correct" else ~s["trials"].correct.values
x, y = rates(s, per1, ok), rates(s, per2, ok)
null = m.null_t(x, y, paired, N_PERM, rng)
t = m.obs_t(x, y, paired)
p = m.p_from_null(t, null)
out.append(pd.DataFrame({"key": s["key"], "region": s["units"][region_col].values, "t": t, "p": p,
"sel": (p < 0.05) & (t > 0)}))
for r in range(n_rep):
xs, ys = m.swap_within_trials(x, y, rng)
ts = m.obs_t(xs, ys, paired)
reps[r].append(pd.DataFrame({"sel": (m.p_from_null(ts, null) < 0.05) & (ts > 0),
"region": s["units"][region_col].values, "region_tip": s["units"]["region_tip"].values,
"region_tip2": s["units"]["region_tip2"].values}))
res = pd.concat(out, ignore_index=True)
for col in ["region_tip", "region_tip2"]:
res[col] = np.concatenate([s["units"][col].values for s in (sessions or S)])
reps = [pd.concat(r) for r in reps]
nulls = pd.DataFrame([r.groupby("region").sel.sum() for r in reps]).fillna(0)
nulls.attrs["tip"] = pd.DataFrame([r.groupby("region_tip").sel.sum() for r in reps]).fillna(0)
nulls.attrs["tip2"] = pd.DataFrame([r.groupby("region_tip2").sel.sum() for r in reps]).fillna(0)
return res, nulls
REGION_READINGS = [("region_tip2", "(a2) tip location, nearest MTL contact"),
("region_tip", "(a1) tip location, electrode-name fallback"), ("region", "(b) electrode name")]
def summarise(res, nulls, regions=m.REGIONS, col="region"):
if col != "region":
res, nulls = res.assign(region=res[col]), nulls.attrs["tip2" if col == "region_tip2" else "tip"]
rows = []
for g in regions:
obs = int(res[res.region == g].sel.sum())
nl = nulls[g] if g in nulls else pd.Series(np.zeros(len(nulls)))
rows.append({"region": g, "neurons": int((res.region == g).sum()), "selected": obs,
"null median": nl.median(), "null 95th pct": np.percentile(nl, 95),
"p": (np.sum(nl >= obs) + 1) / (len(nl) + 1)})
return pd.DataFrame(rows)
MAINT = {}
for name, paired, seed in [("(a) paired", True, 21), ("(b) unpaired", False, 22)]:
res, nulls = select("maintenance", "fixation", paired, seed)
MAINT[name] = res
print(name, "- all regions:", int(res.sel.sum()), "of", len(res))
for col, rname in REGION_READINGS:
print(" regions by reading", rname + ":")
print(summarise(res, nulls, col=col).round(4).to_string(index=False))
print(" per participant:", res.assign(p=res.key.str[1:3]).groupby("p").sel.sum().to_dict(), "\n")
# The claims file's check at 5000 permutations per neuron (paired test, the observed counts only).
N_PERM = 5000
res5, _ = select("maintenance", "fixation", True, 25, n_rep=1)
N_PERM = 1000
print("Paired test with 5000 permutations per neuron: all", int(res5.sel.sum()),
{rname: res5.groupby(col).sel.sum().reindex(m.REGIONS).to_dict() for col, rname in REGION_READINGS})
# How much the counts at 1000 permutations move with the draw of the permutations (10 draws, paired).
draws = [select("maintenance", "fixation", True, 200 + i, n_rep=1)[0] for i in range(10)]
for col, rname in REGION_READINGS:
cnt = pd.DataFrame([d.groupby(col).sel.sum().reindex(m.REGIONS) for d in draws])
print(f"10 draws of 1000 permutations, paired, {rname}: " +
"; ".join(f"{g} median {cnt[g].median():.0f}, range {cnt[g].min()} to {cnt[g].max()}" for g in m.REGIONS))
(a) paired - all regions: 98 of 809
regions by reading (a2) tip location, nearest MTL contact:
region neurons selected null median null 95th pct p
Hipp 467 61 10.0 17.0 0.002
Ent 196 21 4.0 8.0 0.002
Amg 146 16 3.0 6.0 0.002
regions by reading (a1) tip location, electrode-name fallback:
region neurons selected null median null 95th pct p
Hipp 495 62 11.0 17.1 0.002
Ent 168 20 4.0 7.0 0.002
Amg 146 16 3.0 6.0 0.002
regions by reading (b) electrode name:
region neurons selected null median null 95th pct p
Hipp 361 48 8.0 13.0 0.002
Ent 217 28 5.0 9.0 0.002
Amg 226 22 5.0 10.0 0.002
per participant: {'01': 10, '02': 7, '03': 24, '04': 6, '05': 3, '06': 11, '07': 7, '08': 9, '09': 0, '10': 13, '11': 1, '12': 3, '13': 4}
(b) unpaired - all regions: 87 of 809
regions by reading (a2) tip location, nearest MTL contact:
region neurons selected null median null 95th pct p
Hipp 467 54 9.0 15.0 0.002
Ent 196 20 4.0 8.0 0.002
Amg 146 13 3.0 6.0 0.002
regions by reading (a1) tip location, electrode-name fallback:
region neurons selected null median null 95th pct p
Hipp 495 55 9.0 16.0 0.002
Ent 168 19 4.0 7.0 0.002
Amg 146 13 3.0 6.0 0.002
regions by reading (b) electrode name:
region neurons selected null median null 95th pct p
Hipp 361 42 7.0 12.0 0.002
Ent 217 26 5.0 8.0 0.002
Amg 226 18 4.0 8.0 0.002
per participant: {'01': 9, '02': 6, '03': 24, '04': 5, '05': 3, '06': 9, '07': 5, '08': 8, '09': 0, '10': 11, '11': 1, '12': 3, '13': 3}
Paired test with 5000 permutations per neuron: all 94 {'(a2) tip location, nearest MTL contact': {'Hipp': 58, 'Ent': 21, 'Amg': 15}, '(a1) tip location, electrode-name fallback': {'Hipp': 59, 'Ent': 20, 'Amg': 15}, '(b) electrode name': {'Hipp': 47, 'Ent': 28, 'Amg': 19}}
10 draws of 1000 permutations, paired, (a2) tip location, nearest MTL contact: Hipp median 58, range 55 to 59; Ent median 21, range 20 to 22; Amg median 15, range 14 to 16 10 draws of 1000 permutations, paired, (a1) tip location, electrode-name fallback: Hipp median 59, range 56 to 60; Ent median 20, range 19 to 21; Amg median 15, range 14 to 16 10 draws of 1000 permutations, paired, (b) electrode name: Hipp median 47, range 44 to 48; Ent median 28, range 27 to 29; Amg median 20, range 17 to 20
Result: reproduced (resting on an assumption). Maintenance neurons are found in all three regions under every region reading and both tests, each count significant against within-trial shuffles of the period labels (p = 0.002, 499 shuffles); there are 98 in all with the paired test at 1000 permutations per neuron and 94 at 5000, against the paper's 92. With the paired test (reading a) at 5000 permutations, the claims file's check, the counts are 58, 21 and 15 by tip location (a2), within tolerance of 55, 21 and 16, and 59, 20 and 15 with the electrode-name fallback (a1), also within tolerance. At 1000 permutations the counts move with the draw of the permutations: over 10 draws the (a2) counts are 55 to 59, 20 to 22 and 14 to 16, and the draw used for the significance test gives 61, 21 and 16, the hippocampal count just outside tolerance. The unpaired test (reading b) gives 54, 20 and 13 (a2), the amygdala count outside. By electrode name (region reading b) the counts are 47 to 48, 28 and 19 to 22, none within tolerance. The outcome takes the 5000-permutation counts, as the 1000-permutation count is outside tolerance only in some draws, and rests on the tip-location reading of the regions.
fig2d-workload-modulated-maintenance-neurons¶
Among maintenance neurons, a significant number increased firing with workload in hippocampus (10 of 55, 18%) and entorhinal cortex (5 of 21, 24%), and a significant number decreased firing with workload in hippocampus (6 of 55, 11%); in entorhinal cortex 2 decreased and in amygdala 1 increased, not significant.
def workload_counts(maint, seed, sessions=None, n_rep=N_REP, region_col="region"):
"""(a) unpaired permutation t-test of maintenance rates, high (4, 6) against low (1, 2)
workload, correct trials, within the maintenance neurons; null from n_rep shuffles of the
set-size labels across correct trials within each session."""
rng = np.random.default_rng(seed)
rows, reps = [], [[] for _ in range(n_rep)]
for s in (sessions or S):
sel = maint[maint.key == s["key"]].sel.values
if not sel.any():
continue
t_ = s["trials"]; ok = t_.correct.values
r = rates(s, "maintenance", ok)[:, sel]
hi = t_.set_size.values[ok] >= 4
null = m.null_t(r[hi], r[~hi], False, N_PERM, rng)
t = m.obs_t(r[hi], r[~hi], False)
p = m.p_from_null(t, null)
reg = s["units"][region_col].values[sel]
rows.append(pd.DataFrame({"key": s["key"], "region": reg, "t": t, "p": p}))
for k in range(n_rep):
h = rng.permutation(hi)
ts = m.obs_t(r[h], r[~h], False)
ps = m.p_from_null(ts, null)
reps[k].append(pd.DataFrame({"region": reg, "inc": (ps < .05) & (ts > 0), "dec": (ps < .05) & (ts < 0)}))
res = pd.concat(rows, ignore_index=True)
nul = [pd.concat(r).groupby("region")[["inc", "dec"]].sum() for r in reps]
out = []
for g in m.REGIONS:
for d, sign in [("inc", 1), ("dec", -1)]:
obs = int(((res.region == g) & (res.p < .05) & (np.sign(res.t) == sign)).sum())
nl = np.array([n_.loc[g, d] if g in n_.index else 0 for n_ in nul])
out.append({"region": g, "direction": d, "maintenance neurons": int((res.region == g).sum()), "count": obs,
"null median": np.median(nl), "null 95th pct": np.percentile(nl, 95), "p": (np.sum(nl >= obs) + 1) / (len(nl) + 1)})
return pd.DataFrame(out)
for name in MAINT:
for col, rname in REGION_READINGS:
print("Maintenance neurons of fig2c reading", name, "- per-neuron test reading (a); regions:", rname)
print(workload_counts(MAINT[name], 31, region_col=col).round(4).to_string(index=False), "\n")
Maintenance neurons of fig2c reading (a) paired - per-neuron test reading (a); regions: (a2) tip location, nearest MTL contact
region direction maintenance neurons count null median null 95th pct p Hipp inc 61 2 1.0 4.0 0.430 Hipp dec 61 2 1.0 4.0 0.486 Ent inc 21 0 0.0 2.0 1.000 Ent dec 21 2 0.0 2.0 0.088 Amg inc 16 1 0.0 1.0 0.290 Amg dec 16 1 0.0 2.0 0.354 Maintenance neurons of fig2c reading (a) paired - per-neuron test reading (a); regions: (a1) tip location, electrode-name fallback
region direction maintenance neurons count null median null 95th pct p Hipp inc 62 2 1.0 4.0 0.446 Hipp dec 62 2 1.0 4.0 0.492 Ent inc 20 0 0.0 2.0 1.000 Ent dec 20 2 0.0 2.0 0.082 Amg inc 16 1 0.0 1.0 0.290 Amg dec 16 1 0.0 2.0 0.354 Maintenance neurons of fig2c reading (a) paired - per-neuron test reading (a); regions: (b) electrode name
region direction maintenance neurons count null median null 95th pct p Hipp inc 48 2 1.0 3.0 0.298 Hipp dec 48 2 1.0 3.0 0.376 Ent inc 28 1 1.0 2.0 0.504 Ent dec 28 3 1.0 2.0 0.042 Amg inc 22 0 0.0 2.0 1.000 Amg dec 22 0 0.0 2.0 1.000 Maintenance neurons of fig2c reading (b) unpaired - per-neuron test reading (a); regions: (a2) tip location, nearest MTL contact
region direction maintenance neurons count null median null 95th pct p Hipp inc 54 2 1.0 3.0 0.384 Hipp dec 54 2 1.0 4.0 0.414 Ent inc 20 0 0.0 2.0 1.000 Ent dec 20 2 0.0 2.0 0.086 Amg inc 13 1 0.0 1.0 0.226 Amg dec 13 0 0.0 1.0 1.000 Maintenance neurons of fig2c reading (b) unpaired - per-neuron test reading (a); regions: (a1) tip location, electrode-name fallback
region direction maintenance neurons count null median null 95th pct p Hipp inc 55 2 1.0 3.0 0.400 Hipp dec 55 2 1.0 4.0 0.420 Ent inc 19 0 0.0 2.0 1.000 Ent dec 19 2 0.0 2.0 0.080 Amg inc 13 1 0.0 1.0 0.226 Amg dec 13 0 0.0 1.0 1.000 Maintenance neurons of fig2c reading (b) unpaired - per-neuron test reading (a); regions: (b) electrode name
region direction maintenance neurons count null median null 95th pct p Hipp inc 42 2 1.0 3.0 0.258 Hipp dec 42 2 1.0 3.0 0.312 Ent inc 26 1 0.0 2.0 0.476 Ent dec 26 2 0.0 2.0 0.130 Amg inc 18 0 0.0 2.0 1.000 Amg dec 18 0 0.0 2.0 1.000
# Reading (b): cluster-based permutation test over time (Figure 2a legend), rates in 500 ms windows
# stepped by 50 ms over the trial; a neuron counts when a significant cluster overlaps maintenance.
WIN, STEP = 0.5, 0.05
t_off = np.arange(0, m.TRIAL_S - WIN + 1e-9, STEP)
t_mid = t_off + WIN / 2
in_maint = (t_mid + WIN / 2 > 1.7) & (t_mid - WIN / 2 < 2.6)
def cluster_counts(maint, seed, n_perm=1000, n_rep=N_REP):
rng = np.random.default_rng(seed)
rows, reps = [], [[] for _ in range(n_rep)]
for s in S:
sel = maint[maint.key == s["key"]].sel.values
if not sel.any():
continue
t_ = s["trials"]; ok = t_.correct.values
wc = m.window_counts(s, t_off, WIN)[sel][:, ok, :].astype(float) / WIN # units, trials, time
hi = t_.set_size.values[ok] >= 4
n1, n2 = hi.sum(), (~hi).sum()
thr = stats.t.ppf(0.975, n1 + n2 - 2)
def tcurve(h, u):
a, b = wc[u][h], wc[u][~h]
with np.errstate(divide="ignore", invalid="ignore"):
return np.nan_to_num((a.mean(0) - b.mean(0)) / np.sqrt(a.var(0, ddof=1) / len(a) + b.var(0, ddof=1) / len(b)))
perms = [rng.permutation(hi) for _ in range(n_perm)]
reg = s["units"].region.values[sel]
reg_tip = s["units"].region_tip.values[sel]
reg_tip2 = s["units"].region_tip2.values[sel]
for u in range(wc.shape[0]):
null_max = np.array([max([abs(c[2]) for c in m.cluster_masses(tcurve(h, u), thr)] or [0]) for h in perms])
crit = np.percentile(null_max, 95)
def classify(h):
sig = [c for c in m.cluster_masses(tcurve(h, u), thr) if abs(c[2]) > crit and in_maint[c[0]:c[1]].any()]
if not sig:
return 0
return int(np.sign(max(sig, key=lambda c: abs(c[2]))[2]))
rows.append({"region": reg[u], "region_tip": reg_tip[u], "region_tip2": reg_tip2[u], "dir": classify(hi)})
for k in range(n_rep):
reps[k].append((reg[u], reg_tip[u], reg_tip2[u], classify(rng.permutation(hi))))
res = pd.DataFrame(rows)
out = []
for (col, rname), j in zip(REGION_READINGS, [2, 1, 0]):
for g in m.REGIONS:
for d, sign in [("inc", 1), ("dec", -1)]:
obs = int(((res[col] == g) & (res.dir == sign)).sum())
nl = np.array([sum(1 for x in r if x[j] == g and x[3] == sign) for r in reps])
out.append({"regions": rname.split(",")[0], "region": g, "direction": d,
"count": obs, "null median": np.median(nl), "null 95th pct": np.percentile(nl, 95),
"p": (np.sum(nl >= obs) + 1) / (len(nl) + 1)})
return pd.DataFrame(out)
started = time.time()
print("Maintenance neurons of fig2c reading (a) (selected per neuron, so the same under every region reading),")
print("per-neuron test reading (b), 1000 permutations per neuron, 99 label shuffles:")
print(cluster_counts(MAINT["(a) paired"], 32, n_perm=1000, n_rep=99).round(4).to_string(index=False))
print(f"({time.time() - started:.0f} s)")
Maintenance neurons of fig2c reading (a) (selected per neuron, so the same under every region reading), per-neuron test reading (b), 1000 permutations per neuron, 99 label shuffles:
regions region direction count null median null 95th pct p (a2) tip location Hipp inc 1 1.0 2.0 0.53 (a2) tip location Hipp dec 1 1.0 3.0 0.64 (a2) tip location Ent inc 0 0.0 1.0 1.00 (a2) tip location Ent dec 0 0.0 1.0 1.00 (a2) tip location Amg inc 0 0.0 1.1 1.00 (a2) tip location Amg dec 0 0.0 1.0 1.00 (a1) tip location Hipp inc 1 1.0 2.0 0.55 (a1) tip location Hipp dec 1 1.0 3.0 0.65 (a1) tip location Ent inc 0 0.0 1.0 1.00 (a1) tip location Ent dec 0 0.0 1.0 1.00 (a1) tip location Amg inc 0 0.0 1.1 1.00 (a1) tip location Amg dec 0 0.0 1.0 1.00 (b) electrode name Hipp inc 1 0.0 2.0 0.46 (b) electrode name Hipp dec 1 1.0 2.1 0.56 (b) electrode name Ent inc 0 0.0 1.0 1.00 (b) electrode name Ent dec 0 0.0 1.1 1.00 (b) electrode name Amg inc 0 0.0 1.0 1.00 (b) electrode name Amg dec 0 0.0 1.0 1.00 (14 s)
# Variant: every session realigned by k = 3 (data check); maintenance neurons reselected.
S3 = [m.realign(s, 3) for s in S]
for s in S3:
s["counts"] = m.period_counts(s)
res3, _ = select("maintenance", "fixation", True, 24, sessions=S3, n_rep=1)
print("Realigned k = 3: maintenance neurons (a) by region:", res3.groupby("region").sel.sum().to_dict())
print(workload_counts(res3, 33, sessions=S3, n_rep=99).round(4).to_string(index=False))
Realigned k = 3: maintenance neurons (a) by region: {'Amg': 21, 'Ent': 23, 'Hipp': 43, 'other:HL': 0, 'other:HR': 0, 'other:PCL': 0}
region direction maintenance neurons count null median null 95th pct p Hipp inc 43 1 1.0 3.0 0.60 Hipp dec 43 1 1.0 3.0 0.67 Ent inc 23 1 0.0 1.1 0.43 Ent dec 23 1 0.0 2.0 0.44 Amg inc 21 0 0.0 2.0 1.00 Amg dec 21 0 0.0 2.0 1.00
Result: not reproduced. By tip location (a2), among the 61 hippocampal and 21 entorhinal maintenance neurons (paired selection), 2 hippocampal neurons increase with workload (p = 0.43 against shuffled set-size labels), 2 decrease (p = 0.49) and no entorhinal neuron increases; 2 entorhinal neurons decrease (p = 0.09) and 1 amygdala neuron increases (p = 0.29). With the electrode-name fallback (a1) the counts are the same (2, 0 and 2), and by electrode name (b), among 48 and 28 maintenance neurons, the three named counts are 2 (p = 0.30), 1 (p = 0.50) and 2 (p = 0.38). The unpaired selection gives 2, 0 and 2 by tip location and 2, 1 and 2 by electrode name, with no named count below p = 0.25; the cluster-based test (reading b) gives at most 1 per count under every region reading, and the spikes moved by three trials 1, 1 and 1. The counts expected under shuffled labels are 0 to 1. No count reaches 10, 5 or 6, and none of the three named counts is significant.
maintenance-neurons-correct-vs-incorrect¶
More hippocampal neurons were classified as maintenance neurons on correct trials than on incorrect trials (55 of 454, 12%, against 30 of 454, 7%).
PERMS = {}
def hipp_counts(paired, cor_masks, rng, region_col="region"):
"""Hippocampal maintenance neurons on the trials of each mask (per session). One set of
permutations per sample size is drawn once and reused (shuffling keeps the sizes)."""
n = 0
for s, ok in zip(S, cor_masks):
h = (s["units"][region_col] == "Hipp").values
if not h.any():
continue
x, y = rates(s, "maintenance", ok)[:, h], rates(s, "fixation", ok)[:, h]
size = x.shape[0] if paired else 2 * x.shape[0]
if (paired, size) not in PERMS:
PERMS[(paired, size)] = m.make_perms(size, N_PERM, paired, rng)
t = m.obs_t(x, y, paired)
p = m.p_from_null(t, m.null_t(x, y, paired, N_PERM, rng, perms=PERMS[(paired, size)]))
n += int(((p < .05) & (t > 0)).sum())
return n
for name, paired, seed, col in [("(a) paired, regions (b) electrode name", True, 41, "region"),
("(b) unpaired, regions (b) electrode name", False, 42, "region"),
("(a) paired, regions (a1) tip location, electrode-name fallback", True, 43, "region_tip"),
("(b) unpaired, regions (a1) tip location, electrode-name fallback", False, 44, "region_tip"),
("(a) paired, regions (a2) tip location, nearest MTL contact", True, 45, "region_tip2"),
("(b) unpaired, regions (a2) tip location, nearest MTL contact", False, 46, "region_tip2")]:
rng = np.random.default_rng(seed)
cor = [s["trials"].correct.values for s in S]
n_c = hipp_counts(paired, cor, rng, col)
n_i = hipp_counts(paired, [~c for c in cor], rng, col)
started = time.time()
null = []
for r in range(N_REP):
perm = [rng.permutation(c) for c in cor]
null.append(hipp_counts(paired, perm, rng, col) - hipp_counts(paired, [~p for p in perm], rng, col))
null = np.array(null)
p = (np.sum(null >= n_c - n_i) + 1) / (N_REP + 1)
sub = []
for r in range(100):
masks = []
for c in cor:
keep = np.zeros_like(c)
keep[rng.choice(np.where(c)[0], (~c).sum(), replace=False)] = True
masks.append(keep)
sub.append(hipp_counts(paired, masks, rng, col))
print(f"{name}: hippocampal maintenance neurons, correct {n_c}, incorrect {n_i}; difference {n_c - n_i}, "
f"null (correct/incorrect labels shuffled within session, {N_REP}) median {np.median(null):.0f}, "
f"95th pct {np.percentile(null, 95):.0f}, p = {p:.4f} ({time.time() - started:.0f} s)")
print(f" reading (b) of the claim, correct trials subsampled to the incorrect count (100 subsamples): "
f"mean {np.mean(sub):.1f}, range {min(sub)} to {max(sub)}")
print("Incorrect trials per session:", [int((~c).sum()) for c in cor])
(a) paired, regions (b) electrode name: hippocampal maintenance neurons, correct 47, incorrect 13; difference 34, null (correct/incorrect labels shuffled within session, 499) median 31, 95th pct 39, p = 0.3280 (13 s) reading (b) of the claim, correct trials subsampled to the incorrect count (100 subsamples): mean 12.6, range 5 to 25
(b) unpaired, regions (b) electrode name: hippocampal maintenance neurons, correct 42, incorrect 13; difference 29, null (correct/incorrect labels shuffled within session, 499) median 28, 95th pct 35, p = 0.4980 (35 s) reading (b) of the claim, correct trials subsampled to the incorrect count (100 subsamples): mean 11.5, range 5 to 22
(a) paired, regions (a1) tip location, electrode-name fallback: hippocampal maintenance neurons, correct 57, incorrect 19; difference 38, null (correct/incorrect labels shuffled within session, 499) median 41, 95th pct 51, p = 0.7120 (15 s) reading (b) of the claim, correct trials subsampled to the incorrect count (100 subsamples): mean 16.9, range 7 to 29
(b) unpaired, regions (a1) tip location, electrode-name fallback: hippocampal maintenance neurons, correct 55, incorrect 19; difference 36, null (correct/incorrect labels shuffled within session, 499) median 38, 95th pct 47, p = 0.6620 (39 s) reading (b) of the claim, correct trials subsampled to the incorrect count (100 subsamples): mean 15.1, range 6 to 26
(a) paired, regions (a2) tip location, nearest MTL contact: hippocampal maintenance neurons, correct 56, incorrect 18; difference 38, null (correct/incorrect labels shuffled within session, 499) median 40, 95th pct 50, p = 0.6720 (14 s) reading (b) of the claim, correct trials subsampled to the incorrect count (100 subsamples): mean 15.7, range 8 to 27
(b) unpaired, regions (a2) tip location, nearest MTL contact: hippocampal maintenance neurons, correct 54, incorrect 19; difference 35, null (correct/incorrect labels shuffled within session, 499) median 36, 95th pct 45, p = 0.6220 (40 s) reading (b) of the claim, correct trials subsampled to the incorrect count (100 subsamples): mean 15.1, range 6 to 28 Incorrect trials per session: [35, 33, 33, 28, 21, 26, 20, 22, 20, 23, 21, 19, 15, 13, 15, 16]
Result: not reproduced. By tip location (a2), 56 hippocampal neurons are maintenance neurons on correct trials and 18 on incorrect trials (paired test; unpaired 54 and 19), against 55 and 30: the correct count is within tolerance, the incorrect count is not. The difference of 38 is not larger than with correct and incorrect labels shuffled within sessions (median 40, p = 0.67; unpaired p = 0.62). With the electrode-name fallback (a1) the counts are 57 and 19 (p = 0.71), and by electrode name 47 and 13 (p = 0.33; unpaired 42 and 13, p = 0.50). With correct trials subsampled to the number of incorrect ones (reading b) the correct count averages 15.7 (a2), 16.9 (a1) and 12.6 (electrode name), close to the incorrect counts.
fig3b-test-neurons¶
Test neurons, firing more after the test array than during encoding and maintenance, made up 46 of 809 neurons (6%), half the number of maintenance neurons, and were found in all regions (17 of 454 in hippocampus, 19 of 209 in entorhinal cortex, 10 of 146 in amygdala).
TEST = {}
for paired, tag in [(True, "paired"), (False, "unpaired")]:
ra, na = select("test", "maintenance", paired, 51 + paired)
rc, nc = select("test", "encmaint", paired, 53 + paired)
re_, _ = select("test", "encoding", paired, 55 + paired, n_rep=1)
rb = ra.copy(); rb["sel"] = ra.sel.values & re_.sel.values
# (b): null from within-trial shuffles of the three period labels, 99 repetitions, 200 permutations per neuron
rng = np.random.default_rng(57 + paired)
nb = []
for r in range(99):
cnt = pd.Series(dtype=float)
parts = []
for s in S:
ok = s["trials"].correct.values
P = np.stack([rates(s, "test", ok), rates(s, "encoding", ok), rates(s, "maintenance", ok)])
order = np.array([rng.permutation(3) for _ in range(P.shape[1])])
Q = P[order.T, np.arange(P.shape[1])[None, :], :]
okb = np.ones(P.shape[2], bool)
for j in (1, 2):
t = m.obs_t(Q[0], Q[j], paired)
p = m.p_from_null(t, m.null_t(Q[0], Q[j], paired, 200, rng))
okb &= (p < .05) & (t > 0)
parts.append(pd.DataFrame({"sel": okb, "region": s["units"].region.values, "region_tip": s["units"].region_tip.values,
"region_tip2": s["units"].region_tip2.values}))
nb.append(pd.concat(parts))
nbl = nb
nbt = pd.DataFrame([x.groupby("region_tip").sel.sum() for x in nb]).fillna(0)
nb = pd.DataFrame([x.groupby("region").sel.sum() for x in nb]).fillna(0)
nb.attrs["tip"] = nbt
nb.attrs["tip2"] = pd.DataFrame([x.groupby("region_tip2").sel.sum() for x in nbl]).fillna(0)
for name, res, nl in [("(a) test > maintenance", ra, na), ("(b) test > encoding and test > maintenance", rb, nb),
("(c) test > encoding and maintenance joined", rc, nc)]:
TEST[(name, tag)] = res
print(f"{name}, {tag}: all {int(res.sel.sum())} of {len(res)}")
for col, rname in REGION_READINGS:
print(" regions by reading", rname + ":")
print(summarise(res, nl, col=col).round(4).to_string(index=False))
print()
(a) test > maintenance, paired: all 86 of 809 regions by reading (a2) tip location, nearest MTL contact: region neurons selected null median null 95th pct p Hipp 467 35 11.0 18.0 0.002 Ent 196 34 4.0 9.0 0.002 Amg 146 17 3.0 7.0 0.002 regions by reading (a1) tip location, electrode-name fallback: region neurons selected null median null 95th pct p Hipp 495 38 12.0 19.0 0.002 Ent 168 31 4.0 8.0 0.002 Amg 146 17 3.0 7.0 0.002 regions by reading (b) electrode name: region neurons selected null median null 95th pct p Hipp 361 35 9.0 14.1 0.002 Ent 217 40 5.0 10.0 0.002 Amg 226 11 5.0 10.0 0.046 (b) test > encoding and test > maintenance, paired: all 42 of 809 regions by reading (a2) tip location, nearest MTL contact: region neurons selected null median null 95th pct p Hipp 467 22 1.0 4.0 0.01 Ent 196 11 1.0 2.1 0.01 Amg 146 9 0.0 2.0 0.01 regions by reading (a1) tip location, electrode-name fallback: region neurons selected null median null 95th pct p Hipp 495 23 1.0 4.0 0.01 Ent 168 10 0.0 2.0 0.01 Amg 146 9 0.0 2.0 0.01 regions by reading (b) electrode name: region neurons selected null median null 95th pct p Hipp 361 20 1.0 4.0 0.01 Ent 217 18 1.0 2.1 0.01 Amg 226 4 0.0 3.0 0.04 (c) test > encoding and maintenance joined, paired: all 89 of 809 regions by reading (a2) tip location, nearest MTL contact: region neurons selected null median null 95th pct p Hipp 467 43 11.0 18.0 0.002 Ent 196 24 5.0 9.0 0.002 Amg 146 22 4.0 7.0 0.002 regions by reading (a1) tip location, electrode-name fallback: region neurons selected null median null 95th pct p Hipp 495 45 12.0 19.0 0.002 Ent 168 22 4.0 8.0 0.002 Amg 146 22 4.0 7.0 0.002 regions by reading (b) electrode name: region neurons selected null median null 95th pct p Hipp 361 40 9.0 14.0 0.002 Ent 217 32 5.0 10.0 0.002 Amg 226 17 6.0 11.0 0.004
(a) test > maintenance, unpaired: all 79 of 809 regions by reading (a2) tip location, nearest MTL contact: region neurons selected null median null 95th pct p Hipp 467 31 9.0 15.0 0.002 Ent 196 30 4.0 8.0 0.002 Amg 146 18 3.0 6.0 0.002 regions by reading (a1) tip location, electrode-name fallback: region neurons selected null median null 95th pct p Hipp 495 34 10.0 16.0 0.002 Ent 168 27 3.0 7.0 0.002 Amg 146 18 3.0 6.0 0.002 regions by reading (b) electrode name: region neurons selected null median null 95th pct p Hipp 361 31 7.0 12.0 0.002 Ent 217 37 4.0 8.0 0.002 Amg 226 11 4.0 10.0 0.046 (b) test > encoding and test > maintenance, unpaired: all 36 of 809 regions by reading (a2) tip location, nearest MTL contact: region neurons selected null median null 95th pct p Hipp 467 16 1.0 3.0 0.01 Ent 196 9 0.0 2.0 0.01 Amg 146 11 0.0 1.1 0.01 regions by reading (a1) tip location, electrode-name fallback: region neurons selected null median null 95th pct p Hipp 495 17 1.0 3.1 0.01 Ent 168 8 0.0 2.0 0.01 Amg 146 11 0.0 1.1 0.01 regions by reading (b) electrode name: region neurons selected null median null 95th pct p Hipp 361 14 1.0 3.0 0.01 Ent 217 16 1.0 2.0 0.01 Amg 226 6 0.0 2.0 0.01 (c) test > encoding and maintenance joined, unpaired: all 84 of 809 regions by reading (a2) tip location, nearest MTL contact: region neurons selected null median null 95th pct p Hipp 467 42 8.0 15.0 0.002 Ent 196 23 4.0 7.0 0.002 Amg 146 19 3.0 6.0 0.002 regions by reading (a1) tip location, electrode-name fallback: region neurons selected null median null 95th pct p Hipp 495 44 9.0 16.0 0.002 Ent 168 21 3.0 7.0 0.002 Amg 146 19 3.0 6.0 0.002 regions by reading (b) electrode name: region neurons selected null median null 95th pct p Hipp 361 38 6.0 11.1 0.002 Ent 217 30 4.0 8.0 0.002 Amg 226 16 4.0 9.0 0.004
Result: partially reproduced. Test neurons are found in all three regions under every reading and every region reading, each count significant against shuffled period labels (p = 0.002 to 0.046). The totals do not depend on the region reading: test against maintenance (reading a) gives 86 paired and 79 unpaired, test above both encoding and maintenance (reading b) 42 and 36, and test against both periods joined (reading c) 89 and 84, against 46. By tip location (a2), reading (a) paired gives 35, 34 and 17 by region; reading (b) paired 22, 11 and 9, which puts the total and the amygdala count within tolerance of 46 and 10 but not the hippocampal (17) or entorhinal (19) count; reading (b) unpaired 16, 9 and 11, hippocampus and amygdala within tolerance. With the electrode-name fallback (a1) reading (b) paired gives 23, 10 and 9, and by electrode name 20, 18 and 4, the total and the entorhinal count within tolerance. No combination of readings puts all four counts within tolerance. The criterion has no partial clause; the outcome is partial because test neurons are significant in every region and the total under reading (b) matches, while the region counts do not. Under reading (b), paired, test neurons are 43% of the 98 maintenance neurons, about half as the statement says.
fig4a-dpca-separates-set-sizes¶
Demixed PCA of the population separated set sizes 1, 2, 4 and 6; dPCs 2, 3 and 4 explained 36% of the variance (significant against scrambled labels), the separation appeared early in the trial, and the first 15 dPCs explained 81.4% of the variance.
from dPCA import dPCA
DBIN = 0.01 # 10 ms bins; smoothing is at 200 or 167 ms SD
NB = int(round(m.TRIAL_S / DBIN))
for s in S:
s["bins"] = m.window_counts(s, np.arange(NB) * DBIN, DBIN).astype(np.int16)
def kernel(sd):
half = int(round(0.5 / DBIN)); x = np.arange(-half, half + 1) * DBIN # 1000 ms window
k = np.exp(-0.5 * (x / sd) ** 2); return k / k.sum()
def smooth(r, k):
"""Gaussian smoothing along time; at the trial edges the kernel is renormalised to the samples inside."""
return convolve1d(r, k, axis=-1, mode="constant") / convolve1d(np.ones(r.shape[-1]), k, mode="constant")
def cond_trials(s, shuffle_rng=None):
t = s["trials"]; ok = np.where(t.correct.values)[0]
ss = t.set_size.values[ok]
if shuffle_rng is not None:
ss = shuffle_rng.permutation(ss)
return [ok[ss == c] for c in SS]
def population(k, sessions=S, shuffle_rng=None, trial_sets=None):
"""Trial-averaged rates (units, set sizes, time), smoothed, minus each unit's mean fixation rate."""
X = []
for i, s in enumerate(sessions):
idx = trial_sets[i] if trial_sets is not None else cond_trials(s, shuffle_rng)
X.append(np.stack([s["bins"][:, ix, :].mean(1) / DBIN for ix in idx], 1))
X = smooth(np.concatenate(X), k)
return X - X[:, :, :int(0.9 / DBIN)].mean((1, 2))[:, None, None]
def fit_dpca(X, reg=None, n=15, seed=0):
np.random.seed(seed) # the package's randomized SVD draws from numpy's global generator
d = dPCA.dPCA(labels="st", join={"st": ["s", "st"]}, n_components=n, regularizer=reg)
d.protect = ["t"]
d.fit(X)
Xr = X.reshape(X.shape[0], -1); Xr = Xr - Xr.mean(1, keepdims=True); tot = (Xr ** 2).sum()
comps = []
for key in d.D:
for i in range(d.D[key].shape[1]):
z = d.D[key][:, i] @ Xr
comps.append((1 - ((Xr - np.outer(d.P[key][:, i], z)) ** 2).sum() / tot, key, i))
comps.sort(key=lambda c: -c[0]) # explained variance as in the dPCA toolbox (Kobak et al. 2016)
F = np.stack([d.P[k_][:, i] for _, k_, i in comps[:n]], 1); W = np.stack([d.D[k_][:, i] for _, k_, i in comps[:n]], 1)
cum = 1 - ((Xr - F @ (W.T @ Xr)) ** 2).sum() / tot
Z = (W.T @ Xr).reshape(n, *X.shape[1:])
return {"comps": comps, "cum15": cum, "d234": sum(c[0] for c in comps[1:4]),
"setsize": sum(c[0] for c in comps[:n] if c[1] == "st"), "Z": Z, "W": W}
def maint_distance(Z):
sl = slice(int(1.7 / DBIN), int(2.6 / DBIN))
return np.mean([np.linalg.norm(Z[:, a, sl] - Z[:, b, sl], axis=0).mean() for a in range(4) for b in range(a + 1, 4)])
KERNELS = {"(a) SD 200 ms": 0.2, "(b) SD 167 ms": 1 / 6}
FITS = {}
started = time.time()
for kname, sd in KERNELS.items():
X = population(kernel(sd))
FITS[(kname, "(c) none")] = (X, fit_dpca(X))
print(f"population: {X.shape[0]} units x {X.shape[1]} set sizes x {X.shape[2]} time points ({time.time() - started:.0f} s)")
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/dPCA/dPCA.py:422: SyntaxWarning: invalid escape sequence '\s'
\sum_phi ||X[n] - F_\phi D_phi^{-n} X^{-n}||^2
population: 809 units x 4 set sizes x 460 time points (2 s)
# Regularisation reading (d): the toolbox's cross-validated lambda, from trial-by-trial data.
started = time.time()
k = kernel(0.2)
nmax = max(len(ix) for s in S for ix in cond_trials(s))
blocks = []
for s in S:
arr = np.full((nmax, s["bins"].shape[0], 4, NB), np.nan, np.float32)
for c, ix in enumerate(cond_trials(s)):
arr[:len(ix), :, c, :] = np.moveaxis(smooth(s["bins"][:, ix, :].astype(np.float32) / DBIN, k), 1, 0)
blocks.append(arr)
trialX = np.concatenate(blocks, 1); del blocks
np.random.seed(61)
import contextlib, io
with contextlib.redirect_stdout(io.StringIO()):
dcv = dPCA.dPCA(labels="st", join={"st": ["s", "st"]}, n_components=15, regularizer="auto")
dcv.protect = ["t"]; dcv.n_trials = 3
X_unc = np.concatenate([np.stack([s["bins"][:, ix, :].mean(1) / DBIN for ix in cond_trials(s)], 1) for s in S])
X_unc = smooth(X_unc, k)
fixm = X_unc[:, :, :int(0.9 / DBIN)].mean((1, 2))
trialX -= fixm[None, :, None, None]
dcv.fit(X_unc - fixm[:, None, None], trialX)
LAMBDA = dcv.regularizer
del trialX
for kname, sd in KERNELS.items():
X = FITS[(kname, "(c) none")][0]
FITS[(kname, "(d) cross-validated")] = (X, fit_dpca(X, reg=LAMBDA))
print(f"cross-validated lambda (kernel a): {LAMBDA:.3g} ({time.time() - started:.0f} s)")
rows = []
for (kname, rname), (X, f) in FITS.items():
rows.append({"kernel": kname, "regularisation": rname, "dPCs 2-4 %": 100 * f["d234"], "first 15 %": 100 * f["cum15"],
"set-size dPCs in first 15": sum(c[1] == "st" for c in f["comps"][:15]), "set-size variance %": 100 * f["setsize"],
"maintenance distance": maint_distance(f["Z"])})
print(pd.DataFrame(rows).round(2).to_string(index=False))
f = FITS[("(a) SD 200 ms", "(c) none")][1]
print("\nKernel (a), no regularisation, the first 15 dPCs (explained variance %, marginalisation):")
print([(round(100 * v, 2), key) for v, key, _ in f["comps"][:15]])
cross-validated lambda (kernel a): 0.00063 (78 s)
kernel regularisation dPCs 2-4 % first 15 % set-size dPCs in first 15 set-size variance % maintenance distance
(a) SD 200 ms (c) none 32.88 81.41 9 19.50 10.93
(b) SD 167 ms (c) none 30.97 77.74 8 16.64 10.62
(a) SD 200 ms (d) cross-validated 22.21 48.20 9 5.15 1.97
(b) SD 167 ms (d) cross-validated 19.51 43.43 8 4.03 1.81
Kernel (a), no regularisation, the first 15 dPCs (explained variance %, marginalisation):
[(np.float64(24.06), 't'), (np.float64(20.39), 't'), (np.float64(8.12), 't'), (np.float64(4.38), 't'), (np.float64(3.96), 'st'), (np.float64(3.01), 't'), (np.float64(2.79), 'st'), (np.float64(2.46), 'st'), (np.float64(2.17), 'st'), (np.float64(1.96), 't'), (np.float64(1.87), 'st'), (np.float64(1.8), 'st'), (np.float64(1.62), 'st'), (np.float64(1.47), 'st'), (np.float64(1.36), 'st')]
# Null: set-size labels shuffled across correct trials within each session, dPCA refit (99 repetitions).
rng = np.random.default_rng(62)
started = time.time()
null = []
for r in range(99):
Xs = population(kernel(0.2), shuffle_rng=rng)
fs = fit_dpca(Xs, seed=r)
null.append((fs["d234"], fs["setsize"], maint_distance(fs["Z"]), "".join(c[1][-1] for c in fs["comps"][1:4])))
marg = pd.Series([x[3] for x in null]).value_counts()
null = np.array([x[:3] for x in null])
obs = (f["d234"], f["setsize"], maint_distance(f["Z"]))
for j, name in enumerate(["variance of dPCs 2-4", "set-size variance in the first 15 dPCs", "mean distance between set-size trajectories in maintenance"]):
print(f"{name}: observed {obs[j]:.4f}; shuffled median {np.median(null[:, j]):.4f}, 95th pct {np.percentile(null[:, j], 95):.4f}, "
f"max {null[:, j].max():.4f}; p = {(np.sum(null[:, j] >= obs[j]) + 1) / 100:.2f}")
print("Marginalisation of dPCs 2-4 in the shuffled fits (t = time, s = set size, one letter per dPC):", marg.to_dict())
print(f"({time.time() - started:.0f} s)")
Z = f["Z"]
tt = np.arange(NB) * DBIN
fig, axes = plt.subplots(1, 3, figsize=(11, 2.6))
for j, ax in enumerate(axes):
for c, col in zip(range(4), ["tab:blue", "tab:green", "tab:red", "k"]):
ax.plot(tt, Z[j + 1, c], color=col, lw=1)
for x in [0.9, 1.7, 2.6]:
ax.axvline(x, color="grey", lw=0.6)
ax.set_title(f"dPC {j + 2} ({f['comps'][j + 1][1]})"); ax.set_xlabel("s from fixation onset")
plt.tight_layout(); plt.show()
print("Range of the plotted dPCs 2-4:", Z[1:4].min().round(2), Z[1:4].max().round(2))
variance of dPCs 2-4: observed 0.3288; shuffled median 0.3264, 95th pct 0.3325, max 0.3351; p = 0.27
set-size variance in the first 15 dPCs: observed 0.1950; shuffled median 0.1923, 95th pct 0.2025, max 0.2109; p = 0.35
mean distance between set-size trajectories in maintenance: observed 10.9268; shuffled median 10.9558, 95th pct 12.0871, max 12.5837; p = 0.53
Marginalisation of dPCs 2-4 in the shuffled fits (t = time, s = set size, one letter per dPC): {'ttt': 99}
(102 s)
Range of the plotted dPCs 2-4: -11.88 10.34
# Variant: every session realigned by k = 3 (data check).
for s in S3:
s["bins"] = m.window_counts(s, np.arange(NB) * DBIN, DBIN).astype(np.int16)
f3 = fit_dpca(population(kernel(0.2), sessions=S3))
print(f"Realigned k = 3: dPCs 2-4 {100 * f3['d234']:.2f}%, first 15 {100 * f3['cum15']:.2f}%, set-size variance {100 * f3['setsize']:.2f}%, "
f"maintenance distance {maint_distance(f3['Z']):.3f}")
Realigned k = 3: dPCs 2-4 32.59%, first 15 80.34%, set-size variance 19.22%, maintenance distance 12.168
Result: not reproduced. Without regularisation (reading c) the first 15 dPCs explain 81.41% of the variance (kernel a; 77.74% with kernel b) and dPCs 2 to 4 32.88% (30.97%), both within tolerance of 81.4% and 36%; with the cross-validated lambda (reading d) they explain 48.20% and 22.21%. dPCs 2 to 4 are time components here. The set-size components do not separate the set sizes beyond chance: with set-size labels shuffled within sessions the set-size variance is 19.23% (median of 99 shuffles) against 19.50% observed (p = 0.35), the variance of dPCs 2 to 4 gives p = 0.27, and the mean distance between set-size trajectories during maintenance gives p = 0.53. In all 99 shuffled fits dPCs 2 to 4 are also time components, so their variance cannot detect set-size information; the set-size variance of the first 15 dPCs is the test that can, and it is not above the shuffles. The criterion's condition on the null is not met.
fig4b-speed-lowest-in-maintenance¶
Multidimensional speed in the 15-dPC space was highest during encoding and lowest during maintenance, lower than in encoding and test.
PER_SL = {p: slice(int(a / DBIN), int(b / DBIN)) for p, (a, b) in m.PERIODS.items()}
def speed_distance(Z):
sp = np.linalg.norm(np.diff(Z, axis=2), axis=0) / DBIN # (set sizes, time - 1)
dist = np.mean([np.linalg.norm(Z[:, a] - Z[:, b], axis=0) for a in range(4) for b in range(a + 1, 4)], 0)
return ({p: sp[:, sl][:, :sp.shape[1]].mean() for p, sl in PER_SL.items()},
{p: dist[sl].mean() for p, sl in PER_SL.items()})
def runs(reading, n_runs, seed, k_sd=0.2):
rng = np.random.default_rng(seed)
kk = kernel(k_sd)
X_all = population(kk)
out = []
for r in range(n_runs):
if reading == "neurons":
X = X_all[rng.integers(0, X_all.shape[0], int(round(0.2 * X_all.shape[0])))]
else:
sets = [[rng.choice(ix, max(1, int(round(0.2 * len(ix)))), replace=True) for ix in cond_trials(s)] for s in S]
X = population(kk, trial_sets=sets)
fr = fit_dpca(X, seed=r)
sp, di = speed_distance(fr["Z"])
out.append({**{"speed " + p: v for p, v in sp.items()}, **{"distance " + p: v for p, v in di.items()}})
return pd.DataFrame(out)
def perm_t(a, b, rng, n=1000):
t, p = m.perm_ttest(np.asarray(a)[:, None], np.asarray(b)[:, None], False, n, rng)
return t[0], p[0]
RUNS = {}
for reading, n_runs, seed in [("neurons", 500, 71), ("trials", 100, 72)]:
started = time.time()
RUNS[reading] = runs(reading, n_runs, seed)
print(f"reading ({'a' if reading == 'neurons' else 'b'}) subsample {reading}, {n_runs} runs ({time.time() - started:.0f} s)")
print(RUNS[reading].describe().loc[["mean", "50%", "min", "max"]].round(3).to_string())
reading (a) subsample neurons, 500 runs (163 s)
speed fixation speed encoding speed maintenance speed test distance fixation distance encoding distance maintenance distance test
mean 10.425 18.790 12.965 16.470 5.777 5.013 5.375 5.273
50% 10.378 18.779 12.858 15.965 5.773 4.929 5.346 5.246
min 6.762 11.079 8.369 10.494 4.239 3.490 3.615 3.904
max 14.993 28.619 19.005 26.889 7.509 7.796 7.659 7.300
reading (b) subsample trials, 100 runs (103 s)
speed fixation speed encoding speed maintenance speed test distance fixation distance encoding distance maintenance distance test
mean 34.416 51.352 41.341 46.882 30.657 26.049 27.667 26.891
50% 34.541 51.845 41.287 46.964 30.657 26.032 27.637 26.946
min 27.284 43.390 35.318 42.134 25.834 20.871 23.402 24.337
max 41.412 59.047 51.443 52.325 34.376 30.961 32.530 29.778
rng = np.random.default_rng(73)
for reading, R in RUNS.items():
sp = R[[c for c in R if c.startswith("speed")]]
t1, p1 = perm_t(R["speed maintenance"], R["speed encoding"], rng)
t2, p2 = perm_t(R["speed maintenance"], R["speed test"], rng)
print(f"{reading}: mean speed by period {sp.mean().round(3).to_dict()}; highest {sp.mean().idxmax()}, lowest {sp.mean().idxmin()}")
print(f" maintenance vs encoding t = {t1:.1f}, p = {p1:.4f}; maintenance vs test t = {t2:.1f}, p = {p2:.4f}; "
f"runs with maintenance below both: {np.mean((R['speed maintenance'] < R['speed encoding']) & (R['speed maintenance'] < R['speed test'])):.1%}")
f0 = FITS[("(a) SD 200 ms", "(c) none")][1]
sp0, di0 = speed_distance(f0["Z"])
print("\nAll neurons, one fit: speed", {p: round(v, 3) for p, v in sp0.items()}, "; distance", {p: round(v, 3) for p, v in di0.items()})
neurons: mean speed by period {'speed fixation': 10.425, 'speed encoding': 18.79, 'speed maintenance': 12.965, 'speed test': 16.47}; highest speed encoding, lowest speed fixation
maintenance vs encoding t = -35.9, p = 0.0010; maintenance vs test t = -22.1, p = 0.0010; runs with maintenance below both: 98.2%
trials: mean speed by period {'speed fixation': 34.416, 'speed encoding': 51.352, 'speed maintenance': 41.341, 'speed test': 46.882}; highest speed encoding, lowest speed fixation
maintenance vs encoding t = -23.3, p = 0.0010; maintenance vs test t = -14.2, p = 0.0010; runs with maintenance below both: 90.0%
All neurons, one fit: speed {'fixation': np.float64(16.619), 'encoding': np.float64(38.508), 'maintenance': np.float64(24.119), 'test': np.float64(34.454)} ; distance {'fixation': np.float64(11.966), 'encoding': np.float64(10.067), 'maintenance': np.float64(10.927), 'test': np.float64(11.051)}
Result: reproduced. Over 500 dPCA runs on 20% subsamples of neurons (reading a), speed is highest in encoding (mean 18.79) and lower in maintenance (12.97) than in encoding and test (16.47), p = 0.001 for both comparisons, the smallest possible with 1000 permutations; 98.2% of runs have maintenance below both. Subsampling trials (reading b, 100 runs) gives the same order and p values. Speed in fixation (10.43) is lower than in maintenance, so the statement's "lowest during maintenance" does not hold for fixation, a comparison the criterion leaves out.
fig4c-distance-largest-in-maintenance¶
The pairwise distance between set-size trajectories was larger during maintenance than during encoding, and the largest of all periods, consistent with attractors during maintenance.
rng = np.random.default_rng(81)
for reading, R in RUNS.items():
di = R[[c for c in R if c.startswith("distance")]]
t1, p1 = perm_t(R["distance maintenance"], R["distance encoding"], rng)
t2, p2 = perm_t(R["distance maintenance"], R["distance test"], rng)
print(f"{reading}: mean distance by period {di.mean().round(3).to_dict()}; largest {di.mean().idxmax()}")
print(f" maintenance vs encoding t = {t1:.1f}, p = {p1:.4f}; maintenance vs test t = {t2:.1f}, p = {p2:.4f}")
neurons: mean distance by period {'distance fixation': 5.777, 'distance encoding': 5.013, 'distance maintenance': 5.375, 'distance test': 5.273}; largest distance fixation
maintenance vs encoding t = 8.1, p = 0.0010; maintenance vs test t = 2.6, p = 0.0110
trials: mean distance by period {'distance fixation': 30.657, 'distance encoding': 26.049, 'distance maintenance': 27.667, 'distance test': 26.891}; largest distance fixation
maintenance vs encoding t = 6.7, p = 0.0010; maintenance vs test t = 4.0, p = 0.0010
Result: partially reproduced. In the same runs the mean distance between set-size trajectories is larger in maintenance (5.375) than in encoding (5.013), p = 0.001 (reading a; reading b gives 27.67 against 26.05, p = 0.001), and larger than in test (p = 0.011). It is not the largest of the four periods: fixation, before any stimulus, is larger (5.777; 30.66 in reading b). The fig4a null shows that the maintenance distance with all neurons does not exceed that with set-size labels shuffled (p = 0.53), so the difference between periods reflects how the distance varies with the period, not a separation of set sizes; the partial outcome rests on the first condition of the criterion only and carries little weight.
fig4d-decode-workload-all¶
During maintenance, workload (set sizes 1 and 2 against 4 and 6) was decoded from all neurons pooled with median accuracy 87.5%, above chance; leaving out any one participant gave accuracies with IQR 85.7 to 87.6%.
print("Window centres, s from memory array onset:", centres_mem.tolist(), "(maintenance is 0.8 to 1.7 s)")
IA = [list(centres_mem).index(c) for c in (1.2, 1.3)]
def report(r, n, label):
acc = np.median(r, 0) # median over runs, per window
nulls = np.median(n, 0)
stat_a = [acc[i] for i in IA]
p_a = [(np.sum(n[:, i] >= acc[i]) + 1) / (len(n) + 1) for i in IA]
stat_b = np.median(acc)
null_b = np.median(n, 1) # per scrambled run, the median over windows
p_b = (np.sum(null_b >= stat_b) + 1) / (len(n) + 1)
return {"group": label, "(a) 1.2 s": stat_a[0], "p 1.2": p_a[0], "(a) 1.3 s": stat_a[1], "p 1.3": p_a[1],
"(b) median of windows": stat_b, "p (b)": p_b, "scrambled median": np.median(nulls), "scrambled max run": n.max()}
for reading in ["c", "d"]:
f_, g_ = workload_inputs(S)
_, per_fold = m.draw_pseudo(g_, 20, reading, np.random.default_rng(0))
print(f"pseudo-trials ({reading}): {per_fold} per label and fold, so {2 * 19 * per_fold} training and {2 * per_fold} test "
f"pseudo-trials per split (correct trials per session: {min(len(np.concatenate(g)) for g in g_)} to "
f"{max(len(np.concatenate(g)) for g in g_)}; the paper gives 145 to 169 training trials)")
DEC = {}
started = time.time()
for reading in ["c", "d"]:
for z in [True, False]:
n_runs, n_null = (100, 200) if (reading, z) == ("c", True) else (50, 100)
f_, g_ = workload_inputs(S)
r, n = m.decode(f_, g_, n_runs, n_null, 20, reading, z, seed=90 + ord(reading) + z)
DEC[(reading, z)] = (r, n)
row = report(r, n, "all")
print(f"pseudo-trials ({reading}), features {'(e) z-scored' if z else '(f) raw'}; {n_runs} runs, {n_null} scrambled runs")
print(" " + ", ".join(f"{k}: {v:.3f}" if isinstance(v, float) else f"{k}: {v}" for k, v in row.items()))
print(f"({time.time() - started:.0f} s)")
Window centres, s from memory array onset: [0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6] (maintenance is 0.8 to 1.7 s)
pseudo-trials (c): 3 per label and fold, so 114 training and 6 test pseudo-trials per split (correct trials per session: 153 to 178; the paper gives 145 to 169 training trials)
pseudo-trials (d): 5 per label and fold, so 190 training and 10 test pseudo-trials per split (correct trials per session: 153 to 178; the paper gives 145 to 169 training trials)
pseudo-trials (c), features (e) z-scored; 100 runs, 200 scrambled runs group: all, (a) 1.2 s: 0.442, p 1.2: 0.861, (a) 1.3 s: 0.458, p 1.3: 0.761, (b) median of windows: 0.452, p (b): 0.871, scrambled median: 0.498, scrambled max run: 0.692
pseudo-trials (c), features (f) raw; 50 runs, 100 scrambled runs group: all, (a) 1.2 s: 0.483, p 1.2: 0.644, (a) 1.3 s: 0.487, p 1.3: 0.535, (b) median of windows: 0.483, p (b): 0.594, scrambled median: 0.492, scrambled max run: 0.667
pseudo-trials (d), features (e) z-scored; 50 runs, 100 scrambled runs group: all, (a) 1.2 s: 0.440, p 1.2: 0.891, (a) 1.3 s: 0.465, p 1.3: 0.683, (b) median of windows: 0.462, p (b): 0.812, scrambled median: 0.499, scrambled max run: 0.680
pseudo-trials (d), features (f) raw; 50 runs, 100 scrambled runs group: all, (a) 1.2 s: 0.468, p 1.2: 0.743, (a) 1.3 s: 0.465, p 1.3: 0.812, (b) median of windows: 0.466, p (b): 0.802, scrambled median: 0.505, scrambled max run: 0.675 (193 s)
# Leave one participant out, readings (c) and (e), 30 runs each, window readings (a) and (b).
started = time.time()
loo = []
for p in sorted(UNITS.participant.unique()):
keep = [s for s in S if int(s["subject"]) != p]
f_, g_ = workload_inputs(keep)
r, _ = m.decode(f_, g_, 30, 0, 20, "c", True, seed=200 + p)
acc = np.median(r, 0)
loo.append({"left out": p, "(a) 1.2 s": acc[IA[0]], "(a) 1.3 s": acc[IA[1]], "(b)": np.median(acc)})
loo = pd.DataFrame(loo)
print(loo.round(3).to_string(index=False))
print("Quartiles over the 13:", {c: np.percentile(loo[c], [25, 75]).round(3).tolist() for c in loo.columns[1:]})
print(f"({time.time() - started:.0f} s)")
f_, g_ = workload_inputs(S3)
r3, n3 = m.decode(f_, g_, 50, 100, 20, "c", True, seed=300)
print("\nVariant, realigned k = 3, readings (c), (e):", {k: (round(v, 3) if isinstance(v, float) else v) for k, v in report(r3, n3, "all").items()})
left out (a) 1.2 s (a) 1.3 s (b)
1 0.442 0.458 0.458
2 0.433 0.454 0.456
3 0.442 0.467 0.462
4 0.462 0.462 0.460
5 0.454 0.475 0.452
6 0.433 0.450 0.454
7 0.458 0.467 0.469
8 0.442 0.454 0.442
9 0.433 0.442 0.450
10 0.458 0.496 0.467
11 0.454 0.450 0.456
12 0.467 0.467 0.467
13 0.442 0.454 0.458
Quartiles over the 13: {'(a) 1.2 s': [0.442, 0.458], '(a) 1.3 s': [0.454, 0.467], '(b)': [0.454, 0.462]}
(74 s)
Variant, realigned k = 3, readings (c), (e): {'group': 'all', '(a) 1.2 s': np.float64(0.587), 'p 1.2': np.float64(0.079), '(a) 1.3 s': np.float64(0.583), 'p 1.3': np.float64(0.069), '(b) median of windows': np.float64(0.585), 'p (b)': np.float64(0.04), 'scrambled median': np.float64(0.506), 'scrambled max run': np.float64(0.683)}
Result: not reproduced. Decoding workload from all 809 neurons gives median accuracies of 0.442 and 0.458 in the windows centred at 1.2 and 1.3 s after memory onset (reading a) and 0.452 over the maintenance windows (reading b), with p = 0.76 to 0.87 against scrambled labels; the other pseudo-trial and feature readings give 0.440 to 0.487 and p of 0.54 or more. Leaving out a participant gives quartiles of 0.442 to 0.467. None is near 87.5% or significant. The training sets hold 114 pseudo-trials per split under reading (c) and 190 under reading (d), where the paper gives 145 to 169; neither reading uses each session's whole correct-trial set. With the spikes moved by three trials (variant) accuracy is 0.583 to 0.587 (p = 0.04 to 0.08). The accuracy as deposited is below the scrambled runs (0.50), which fits the misalignment of spikes and labels suggested by the data check.
fig4d-decode-workload-regions¶
Workload was decoded during maintenance from hippocampal neurons (75%) and entorhinal neurons (80%) but not from amygdala neurons.
started = time.time()
rows = []
for g in m.REGIONS:
masks = [(s["units"].region == g).values for s in S]
f_, g_ = workload_inputs(S, masks)
r, n = m.decode(f_, g_, 100, 200, 20, "c", True, seed=400 + m.REGIONS.index(g))
rows.append(report(r, n, f"{g}, electrode name ({int(sum(x.sum() for x in masks))} neurons)"))
for col, tag, seed in [("region_tip", "(a1)", 410), ("region_tip2", "(a2)", 420)]:
masks = [(s["units"][col] == g).values for s in S]
if col == "region_tip2" and all((x == y).all() for x, y in zip(masks, [(s["units"].region_tip == g).values for s in S])):
rows.append({**rows[-1], "group": f"{g}, tip location {tag}: the same neurons as (a1)"})
continue
f_, g_ = workload_inputs(S, masks)
r, n = m.decode(f_, g_, 100, 200, 20, "c", True, seed=seed + m.REGIONS.index(g))
rows.append(report(r, n, f"{g}, tip location {tag} ({int(sum(x.sum() for x in masks))} neurons)"))
f3_, g3_ = workload_inputs(S3, [(s["units"].region == g).values for s in S3])
r3, n3 = m.decode(f3_, g3_, 30, 100, 20, "c", True, seed=500 + m.REGIONS.index(g))
rows.append(report(r3, n3, f"{g}, electrode name, realigned k = 3 (30 runs, 100 scrambled)"))
print(pd.DataFrame(rows).round(3).to_string(index=False))
print("Bonferroni threshold over 4 groups: p < .0125;", f"({time.time() - started:.0f} s)")
group (a) 1.2 s p 1.2 (a) 1.3 s p 1.3 (b) median of windows p (b) scrambled median scrambled max run
Hipp, electrode name (361 neurons) 0.492 0.552 0.529 0.303 0.506 0.468 0.500 0.717
Hipp, tip location (a1) (495 neurons) 0.483 0.697 0.492 0.602 0.487 0.706 0.506 0.717
Hipp, tip location (a2) (467 neurons) 0.467 0.741 0.479 0.557 0.481 0.657 0.500 0.683
Hipp, electrode name, realigned k = 3 (30 runs, 100 scrambled) 0.571 0.168 0.554 0.238 0.558 0.089 0.508 0.667
Ent, electrode name (217 neurons) 0.454 0.786 0.458 0.806 0.456 0.841 0.500 0.667
Ent, tip location (a1) (168 neurons) 0.467 0.726 0.467 0.682 0.467 0.796 0.500 0.683
Ent, tip location (a2) (196 neurons) 0.475 0.682 0.500 0.433 0.487 0.572 0.494 0.725
Ent, electrode name, realigned k = 3 (30 runs, 100 scrambled) 0.571 0.079 0.542 0.228 0.571 0.010 0.500 0.667
Amg, electrode name (226 neurons) 0.467 0.761 0.458 0.816 0.462 0.866 0.500 0.708
Amg, tip location (a1) (146 neurons) 0.467 0.711 0.475 0.682 0.467 0.811 0.506 0.725
Amg, tip location (a2): the same neurons as (a1) 0.467 0.711 0.475 0.682 0.467 0.811 0.506 0.725
Amg, electrode name, realigned k = 3 (30 runs, 100 scrambled) 0.512 0.396 0.517 0.297 0.519 0.267 0.498 0.692
Bonferroni threshold over 4 groups: p < .0125; (375 s)
Result: not reproduced. By tip location (a2), hippocampal neurons (467) give median accuracies of 0.467 to 0.481 and entorhinal neurons (196) 0.475 to 0.500, with p of 0.43 or more; with the electrode-name fallback (a1), 0.483 to 0.492 (495 neurons) and 0.467 (168 neurons), p of 0.60 or more; by electrode name, 0.492 to 0.529 (361 neurons) and 0.454 to 0.458 (217 neurons), p of 0.30 or more. The paper reports 75% and 80% with p < .0125. Amygdala neurons are not significant under any region reading (0.458 to 0.475, p of 0.68 or more), as the paper finds, so one of three regions matches. With the spikes moved by three trials (electrode name, 100 scrambled runs), hippocampus gives 0.554 to 0.571 (p of 0.089 or more) and entorhinal cortex 0.542 to 0.571 (p = 0.010 for the median over windows, 0.079 and 0.23 for the single windows).
decode-performance-not-significant¶
Trial-wise performance (correct against incorrect) could not be decoded from MTL firing in any group (all 55%, hippocampus 55%, entorhinal 55%, amygdala 50%; p from .40 to .56).
def perf_inputs(sessions, unit_mask=None):
feats, groups = [], []
for i, s in enumerate(sessions):
t = s["trials"]
wc = m.window_counts(s, offs, 1.0)
if unit_mask is not None:
wc = wc[unit_mask[i]]
if wc.shape[0] == 0:
continue
feats.append(np.moveaxis(wc, 1, 0))
groups.append([np.where(t.correct.values)[0], np.where(~t.correct.values)[0]])
return feats, groups
k_min = min(int((~s["trials"].correct).sum()) for s in S)
print("Smallest number of incorrect trials in a session:", k_min, "-> trials per label; cross-validation splits:", k_min)
started = time.time()
rows = []
for g, col in [("all", None)] + [(g, c) for c in ["region", "region_tip", "region_tip2"] for g in m.REGIONS]:
masks = None if col is None else [(s["units"][col] == g).values for s in S]
f_, g_ = perf_inputs(S, masks)
r, n = m.decode(f_, g_, 50, 100, k_min, "c", True, seed=600 + len(g) + (col == "region_tip") + 2 * (col == "region_tip2"))
rows.append(report(r, n, g if col is None else g + {"region": ", (b) electrode name", "region_tip": ", (a1) tip location",
"region_tip2": ", (a2) tip location"}[col]))
print(pd.DataFrame(rows).round(3).to_string(index=False))
print(f"({time.time() - started:.0f} s)")
Smallest number of incorrect trials in a session: 13 -> trials per label; cross-validation splits: 13
group (a) 1.2 s p 1.2 (a) 1.3 s p 1.3 (b) median of windows p (b) scrambled median scrambled max run
all 0.462 0.634 0.538 0.436 0.500 0.545 0.500 0.846
Hipp, (b) electrode name 0.500 0.614 0.500 0.634 0.500 0.644 0.500 0.885
Ent, (b) electrode name 0.500 0.495 0.538 0.376 0.538 0.267 0.471 0.846
Amg, (b) electrode name 0.462 0.673 0.462 0.703 0.462 0.713 0.500 0.808
Hipp, (a1) tip location 0.462 0.772 0.462 0.772 0.462 0.762 0.500 0.808
Ent, (a1) tip location 0.500 0.604 0.538 0.426 0.519 0.505 0.510 0.885
Amg, (a1) tip location 0.500 0.535 0.500 0.515 0.462 0.693 0.500 0.885
Hipp, (a2) tip location 0.481 0.634 0.500 0.604 0.500 0.545 0.500 0.885
Ent, (a2) tip location 0.538 0.376 0.538 0.416 0.538 0.366 0.500 0.846
Amg, (a2) tip location 0.462 0.653 0.462 0.653 0.462 0.683 0.500 0.885
(85 s)
Result: reproduced. With 13 trials per label and 13 cross-validation splits, decoding correct against incorrect trials gives 0.462 to 0.538 for all neurons and for each region under every region reading and both window readings, with p from 0.27 to 0.77 against scrambled labels; no group reaches p < .0125. The same pipeline does not decode workload (fig4d), which the paper decodes at 87.5%, and 13 folds of one trial per label give a test of little demonstrated sensitivity, so this outcome carries little weight.
Summary¶
| Claim | Outcome | Main numbers |
|---|---|---|
| behavior-accuracy-by-set-size | reproduced | 97.8, 97.9, 84.3, 72.6% correct, 88.2% overall; F(3, 45) = 89.78, p = .0002 |
| behavior-capacity | reproduced | per participant KP median 3.23 (2.25 to 5.00), KC median 3.00 (2.11 to 4.50) |
| table1-capacity-and-neurons | reproduced | all 13 Kmax values and all 13 neuron counts equal Table 1 |
| behavior-rt-set-size | partially reproduced | 2678 correct trials; F = 48.47 on median RT, p = .0002; slope 69 to 76 ms/item against 118 (118.1 per condition step, a reading asked of a reader) |
| neurons-recorded | reproduced (assumed) | 809 neurons, 50.6 +- 25.4 (11 to 96) per session; Hipp 467, Ent 196, Amg 146 by tip location (a2), 495, 168, 146 with the electrode-name fallback (a1), 361, 217, 226 by electrode name, against 454, 209, 146 |
| spike-sorting-quality | not reproduced | ISI < 3 ms 0.79 +- 1.77%, 62.5% of units below 0.5%; SNR 5.93 +- 2.31 against 3.81 +- 0.82 |
| firing-rates-by-period | reproduced | all medians and quartiles within 10% |
| fig2c-maintenance-neurons | reproduced (assumed) | paired, 5000 permutations: 58, 21, 15 by tip location (a2), 59, 20, 15 (a1), 47, 28, 19 by electrode name (each p = .002), against 55, 21, 16; 61, 21, 16 (a2) in the 1000-permutation draw |
| fig2d-workload-modulated-maintenance-neurons | not reproduced | Hipp increase 2, Ent increase 0 or 1, Hipp decrease 2 under either region reading; none significant |
| maintenance-neurons-correct-vs-incorrect | not reproduced | 56 against 18 hippocampal neurons by tip location (p = .67), 47 against 13 by electrode name (p = .33) |
| fig3b-test-neurons | partially reproduced | test neurons significant in all regions; counts 86, 42 or 89 by reading against 46 |
| fig4a-dpca-separates-set-sizes | not reproduced | dPCs 2 to 4 32.9% and first 15 81.4%, but set-size variance does not exceed shuffled labels (p = .35) |
| fig4b-speed-lowest-in-maintenance | reproduced | speed highest in encoding, maintenance below encoding and test (p = .001) |
| fig4c-distance-largest-in-maintenance | partially reproduced | maintenance distance above encoding (p = .001) but fixation is the largest |
| fig4d-decode-workload-all | not reproduced | median accuracy 44 to 49%, not above scrambled labels |
| fig4d-decode-workload-regions | not reproduced | Hipp 47 to 53%, Ent 45 to 50%, Amg 46 to 48% under every region reading, none significant |
| decode-performance-not-significant | reproduced | 46 to 54%, no group below p = .0125 |
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 paper's analysis code is not public, so every method comes from the Methods as the claims file reads them, with the readings it names run side by side. Every per-region result is given under three region readings. Two follow reading (a) of the claims file, the location of the microwire tips, taken as the atlas label of each bundle's tip contact; for the 133 units whose tip contact has no MTL label, (a2) takes the label of the contact nearest the tip with an MTL label (and the electrode name where none has one), and (a1) the electrode name. The third is reading (b), the electrode name. The outcomes rest on (a2), which was chosen after seeing the counts it gives and differs from (a1) only for 28 units of participant 8, placed by contacts that lie further from the wire tips than the tip contact. The files hold no location for the microwires themselves. The number of label-scrambled repetitions is 499 for the neuron counts as in the claims file, but 99 for the dPCA null, the cluster-based reading and reading (b) of the test neurons, and the decoders use 100 resampling runs and 200 scrambled runs (50 and 100 for the other readings and for performance), where the paper used 500 and 1000. The per-neuron tests use 1000 permutations, and 5000 for the fig2c counts, the claims file's check. Spike counts for dPCA are taken in 10 ms bins before smoothing, in place of 1 ms. The speed and distance of Figure 4b and 4c use the claims file's averaging, which the paper does not define.
What was found about the dataset.
- Behavior matches the paper exactly where it can be compared: 2678 correct trials, F = 89.78 for accuracy (paper 89.79), all Table 1 capacities, and F = 48.47 for response time when the per-session value is the median RT.
- Participant 12's two sessions are each deposited twice (four files), and both carry the same start time.
- Sessions hold 174 to 192 trials where the description states 192; the missing trials are not documented.
- Spikes are stored only inside trial windows (fixation onset to test offset).
- The
locationof each unit is the atlas label of the macro contact whose number is the unit's microwire number, not the location of the microwire. The electrode name does not give the paper's three region counts; the atlas labels of each bundle's contacts, read from the tip, give 467, 196 and 146 against 454, 209 and 146. The files do not hold the paper's assignment of units to regions. Five units are on bundles (PCL, HL, HR) that the subject descriptions do not list. - The response timestamps put 1385 of 3038 responses before the test array. Read as the stored response time from test onset plus a 1.7 s offset the conversion adds, they are plausible.
- The workload information in the spikes peaks when the spikes stored in trial window i are paired with the label of trial i + 3. No session has set-size information at k = 0 beyond shuffled labels; pooled over the 7 sessions with 192 trials the information peaks at k = 3 (corrected p = 0.005), and over the 9 sessions with fewer trials at no shift. Single-neuron and population measures of workload are at chance as deposited (decoding about 45%, below the scrambled 50%), and decoding rises to about 59% with every session moved by three trials. Session by session the shift by three trials is supported only in P03S1 and P06S1 among the complete sessions (and weakly in P04S1 and P09S1); the pooled peak is weighted toward those two sessions. The spikes and the trial table may be misaligned, but a single offset is not established. Moving every session by three trials recovers little: fig2d counts of 1, 1 and 1, dPCA set-size variance at the null median, and decoding at 0.585 (p = .04) against the paper's 87.5%. A misalignment could explain part of why the claims that use trial labels (workload, set size, correct against incorrect) fail while those that do not (period contrasts, rates, behavior) hold, but the variants do not show that it explains the failures.
- The paper's 118 ms/item equals the slope of the median response times on the condition step (1 to 4), not on the set size (68.8 ms/item); asked as a question for a reader.
print(f"Notebook run time: {(time.time() - NOTEBOOK_START) / 60:.1f} min")
Notebook run time: 26.0 min