Reproducing Chandravadia et al. (2020) from dandiset 000004¶
Paper. Chandravadia N, Liang D, Schjetnan AGP, Carlson A, Faraut M, Chung JM, Reed CM, Dichter B, Maoz U, Kalia SK, Valiante TA, Mamelak AN, Rutishauser U (2020). A NWB-based dataset and processing pipeline of human single-neuron activity during a declarative memory task. Scientific Data 7, 78. doi:10.1038/s41597-020-0415-9. The article is open access (CC BY 4.0); this attempt used its full text, including Tables 1 and 2 and Online-only Table 1.
Dandiset. DANDI:000004, version 0.220126.1852. The paper cites no DANDI version: its data statement points to the Open Science Framework (doi:10.17605/OSF.IO/HV7JA), and the dandiset was created on 2020-03-16, after the paper appeared. The archive has two published versions, 0.220126.1851 and 0.220126.1852, published a minute apart with the same 87 assets (checked in the first data cell) and metadata that differs only in version fields. The most recent is pinned. The draft holds the same files.
What the dandiset holds. One NWB file per session: 87 sessions of 59 subjects, with the trials table, the TTL event stream, the stimulus images, the spike-sorted units with mean waveforms, and the electrodes table. There are no raw broadband recordings.
Code read. rutishauserlab/recogmem-release-NWB
at commit 1d95162 (2025-01-14): the Python and MATLAB export and analysis code, the session list
defineNOsessions_release.ini, and plotFigure5.ipynb. The session list and the demo script
python/export/demo/export.py are fetched from that commit where the notebook relies on them, and
the lines used are printed. The code settles what the Methods leave
open (the AUC, the z-ROC, the high/low confidence split, the label coding). Where it disagrees with
the Methods (the analysis window, and correct trials only for memory selectivity), the Methods give
the main result and the code's choice is shown as a variant. plotFigure5.ipynb writes the Figure
5f percentages into the plot as constants, so the code does not show how they were computed.
Outside inputs. Online-only Table 1 of the paper is copied into the helper module as text, to compare with the files' subject metadata. Table 1 of Faraut et al. 2018 (Scientific Data 5, 180010, doi:10.1038/sdata.2018.10), the paper describing the earlier release, is copied the same way; it lists that release's sessions and is used only to mark which sessions of the dandiset are new (previous-release-link and previous-release-additions).
| Claim | Status |
|---|---|
| counts-neurons-subjects: 1863 neurons, 59 subjects | tested |
| counts-sessions: 87 sessions of Online-only Table 1 | tested |
| table1-demographics: age and sex of every subject | tested |
| task-structure: 100 encoding and 100 recognition trials, 5 categories | tested |
| new-old-label-coding: label 0 is old, 1 is new | tested |
| ttl-codes: event codes of Table 1 | tested |
| experiment-ids: block ids of Table 2 | tested |
| session-dates: day set to the first of the month | tested |
| units-table: columns, 100 kHz waveforms, one electrode per unit | tested |
| electrodes-table: area, coordinates, impedance, filter, channel | tested |
| isi-criterion: fewer than 3% of ISIs below 3 ms | tested |
| fig5a-mtl-locations: all units in amygdala or hippocampus | tested |
| fig5c-auc-all: AUC 0.74 ± 0.10 | tested |
| fig5c-auc-sites: AUC 0.74 ± 0.09 (CS) and 0.69 ± 0.13 (TWH) | tested |
| fig5b-roc-asymmetric: asymmetric ROC | tested |
| fig5d-confidence-accuracy: high confidence more accurate | tested |
| vs-cells-proportion: 343 VS cells (18.4%) | tested |
| ms-cells-proportion: 146 MS cells (7.8%) | tested |
| dual-cells-proportion: 26 dual cells (1.4%) | tested |
| fig5f-site-proportions: comparable proportions at CS and TWH | tested |
| twh-sessions: 9 TWH sessions | tested |
| fig4-example-vs-neuron: neuron 132-177-4-RA | tested |
| previous-release-link: earlier release's sessions found here | tested (with outside input) |
| previous-release-additions: 23 sessions, 17 subjects, 288 neurons added | tested (with outside input) |
| fig3-example-vs-neuron: neuron 144-140-3-RH | blocked (metadata insufficient) |
In [1]:
import time
import warnings
from itertools import combinations
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import recogmem as rm
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
sessions = rm.load_all()
by_id = {s["session"]: s for s in sessions}
print(f"loaded {len(sessions)} files from dandiset {rm.DANDISET_ID} version {rm.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
loaded 87 files from dandiset 000004 version 0.220126.1852 in 72 s
Data check¶
Each check of the metadata against the values is stated under its output.
In [2]:
info = pd.DataFrame([{
"session": s["session"], "label": s["label"], "subject": s["subject_id"], "site": s["site"],
"institution": s["institution"], "start": s["session_start_time"][:10], "age": s["age"], "sex": s["sex"],
"trials": len(s["trials"]), "learn": int((s["trials"].stim_phase == "learn").sum()),
"recog": int((s["trials"].stim_phase == "recog").sum()), "units": len(s["spikes"]),
"electrodes": len(s["electrodes"]), "categories": s["trials"].stimCategory.nunique(),
"experiment_ids": tuple(int(x) for x in np.unique(s["experiment_ids"])),
"data_collection": s["data_collection"]} for s in sessions])
print("Files:", len(info), " subjects:", info.subject.nunique(), " species:", sorted({s["species"] for s in sessions}))
print("\nBy site (CS = subjects labelled H or C, TWH = labelled T):")
print(info.groupby("site").agg(sessions=("session", "size"), subjects=("subject", "nunique"), units=("units", "sum")).to_string())
print("\nInstitution strings:", info.institution.value_counts().to_dict())
print("Devices:", pd.Series([d for s in sessions for d in s["devices"]]).value_counts().to_dict())
print("\nTrials per file (learn, recog):", info.groupby(["learn", "recog"]).size().to_dict())
print("Units per file: min", info.units.min(), "median", info.units.median(), "max", info.units.max(), "total", info.units.sum())
print("Is there a session 144 (Figure 3 example)?", 144 in by_id)
client = rm.DandiAPIClient()
versions = {}
for v in ("0.220126.1851", "0.220126.1852", "draft"):
versions[v] = {a.path: (a.identifier, a.size) for a in client.get_dandiset(rm.DANDISET_ID, v).get_assets()}
print("Published versions 0.220126.1851 and 0.220126.1852 hold the same assets (ids and sizes):", versions["0.220126.1851"] == versions["0.220126.1852"])
print("The draft holds the same paths and sizes as the pinned version:",
{p: s for p, (_, s) in versions["draft"].items()} == {p: s for p, (_, s) in versions["0.220126.1852"].items()})
print()
print(info.drop(columns=["institution", "data_collection"]).to_string(index=False))
Files: 87 subjects: 59 species: ['Human']
By site (CS = subjects labelled H or C, TWH = labelled T):
sessions subjects units
site
CS 78 51 1787
TWH 9 8 77
Institution strings: {'Cedars-Sinai Medical Center': 44, 'Hunigton Memorial Hospital': 28, 'Toronto Western Hospital': 9, 'Hunington Memorial Hospital': 6}
Devices: {'Neuralynx-Atlas': 53, 'Neuralynx-cheetah': 29, 'Neuralynx-Cheetah': 5}
Trials per file (learn, recog): {(50, 100): 7, (97, 100): 1, (100, 100): 79}
Units per file: min 1 median 18.0 max 64 total 1864
Is there a session 144 (Figure 3 example)? False
Published versions 0.220126.1851 and 0.220126.1852 hold the same assets (ids and sizes): True
The draft holds the same paths and sizes as the pinned version: True
session label subject site start age sex trials learn recog units electrodes categories experiment_ids
5 H09 P9HMH CS 2006-03-01 55 M 200 100 100 24 14 5 (80, 81)
6 H09 P9HMH CS 2006-03-01 55 M 200 100 100 19 13 5 (83, 84)
7 H10 P10HMH CS 2006-09-01 37 M 200 100 100 38 18 5 (80, 81)
9 H11 P11HMH CS 2006-11-01 16 M 200 100 100 25 9 5 (80, 81)
17 H14 P14HMH CS 2007-06-01 31 M 200 100 100 23 9 5 (83, 84)
18 H14 P14HMH CS 2007-06-01 31 M 150 50 100 33 10 5 (80, 81)
20 H15 P15HMH CS 2007-09-01 45 M 150 50 100 20 11 5 (80, 81)
21 H15 P15HMH CS 2007-09-01 45 M 150 50 100 14 7 5 (83, 84)
23 H16 P16HMH CS 2007-10-01 34 F 150 50 100 18 9 5 (83, 84)
24 H16 P16HMH CS 2007-10-01 34 F 200 100 100 21 18 5 (80, 81)
26 H18 P18HMH CS 2008-06-01 40 M 150 50 100 33 16 5 (80, 81)
27 H19 P19HMH CS 2008-06-01 34 M 150 50 100 16 6 5 (83, 84)
28 H19 P19HMH CS 2008-06-01 34 M 200 100 100 15 7 5 (80, 81)
32 H17 P17HMH CS 2008-05-01 19 M 200 100 100 43 16 5 (80, 81)
38 H21 P21HMH CS 2009-01-01 20 M 200 100 100 20 8 5 (83, 84)
39 H21 P21HMH CS 2009-01-01 20 M 200 100 100 18 8 5 (80, 81)
41 H21 P21HMH CS 2009-01-01 20 M 200 100 100 15 6 5 (88, 89)
43 H23 P23HMH CS 2009-02-01 40 M 200 100 100 23 8 5 (88, 89)
44 H23 P23HMH CS 2009-02-01 40 M 150 50 100 23 8 5 (83, 84)
47 H23 P23HMH CS 2009-02-01 40 M 200 100 100 17 7 5 (80, 81)
48 H28 P28HMH CS 2010-03-13 22 M 200 100 100 45 19 5 (88, 89)
49 H29 P29HMH CS 2010-12-01 17 F 200 100 100 27 10 5 (88, 89)
50 H31 P31HMH CS 2011-04-01 30 M 200 100 100 14 7 5 (83, 84)
52 H33 P33HMH CS 2011-07-01 30 M 200 100 100 18 8 5 (80, 81)
54 H42 P42HMH CS 2013-03-01 29 M 200 100 100 32 13 5 (80, 81)
55 H42 P42HMH CS 2013-03-01 29 M 200 100 100 32 13 5 (80, 82)
56 H43 P43HMH CS 2013-04-01 27 F 200 100 100 15 6 5 (80, 81)
58 H27 P27HMH CS 2010-02-01 40 M 200 100 100 49 14 5 (83, 84)
59 CS24 P24CS CS 2012-09-01 57 F 200 100 100 21 7 5 (80, 81)
60 CS24 P24CS CS 2012-09-01 57 F 200 100 100 22 9 5 (83, 84)
61 CS25 P25CS CS 2012-09-01 36 F 200 100 100 8 6 5 (80, 81)
63 CS25 P25CS CS 2012-09-01 36 F 200 100 100 14 11 5 (83, 84)
64 CS26 P26CS CS 2012-12-01 56 F 200 100 100 33 9 5 (80, 81)
66 CS26 P26CS CS 2012-12-01 56 F 200 100 100 11 6 5 (83, 84)
67 CS27 P27CS CS 2013-01-01 44 M 200 100 100 15 8 5 (80, 81)
68 H44 P44HMH CS 2000-01-01 57 F 200 100 100 19 7 5 (80, 81)
69 CS29 P29CS CS 2013-10-01 19 M 200 100 100 64 17 5 (80, 81)
70 CS29 P29CS CS 2013-10-01 19 M 200 100 100 33 13 5 (80, 81)
73 CS31 P31CS CS 2014-02-01 32 M 200 100 100 4 2 5 (80, 81)
74 CS32 P32CS CS 2014-02-01 19 M 200 100 100 35 9 5 (80, 81)
76 CS33 P33CS CS 2014-03-01 44 F 200 100 100 53 14 5 (83, 84)
77 CS33 P33CS CS 2014-03-01 44 F 200 100 100 30 11 5 (81, 82)
78 CS33 P33CS CS 2014-03-01 44 F 200 100 100 13 9 5 (80, 81)
85 CS34 P34CS CS 2014-12-01 70 M 200 100 100 6 4 5 (83, 84)
92 H47 P47HMH CS 2014-06-01 20 M 200 100 100 7 3 5 (80, 81)
93 CS39 P39CS CS 2015-08-01 26 M 200 100 100 35 13 5 (80, 81)
96 CS37 P37CS CS 2015-03-01 33 F 200 100 100 13 8 5 (80, 81)
97 H47 P47HMH CS 2014-06-01 20 M 200 100 100 10 7 5 (80, 81)
99 H48 P48HMH CS 2014-12-01 32 F 200 100 100 30 20 5 (80, 81)
100 CS39 P39CS CS 2015-08-01 26 M 200 100 100 42 16 5 (82, 83)
101 CS40 P40CS CS 2016-01-01 25 M 200 100 100 14 10 5 (80, 81)
102 CS38 P38CS CS 2015-07-01 63 F 200 100 100 44 17 5 (80, 81)
104 H51 P51HMH CS 2016-02-01 24 M 200 100 100 42 19 5 (80, 81)
105 H51 P51HMH CS 2016-02-01 24 M 200 100 100 43 19 5 (82, 83)
111 CS42 P42CS CS 2016-08-01 25 F 200 100 100 14 10 5 (80, 81)
112 CS42 P42CS CS 2016-08-01 25 F 200 100 100 12 7 5 (80, 81)
113 CS43 P43CS CS 2016-11-01 42 F 200 100 100 21 14 5 (80, 81)
114 CS44 P44CS CS 2016-09-01 53 F 200 100 100 15 9 5 (80, 81)
115 CS47 P47CS CS 2017-02-01 32 M 200 100 100 35 17 5 (80, 81)
116 CS48 P48CS CS 2017-03-01 32 F 200 100 100 36 14 5 (80, 81)
117 CS49 P49CS CS 2017-05-01 24 F 200 100 100 11 7 5 (80, 81)
118 CS49 P49CS CS 2017-05-01 24 F 200 100 100 13 11 5 (80, 81)
119 CS51 P51CS CS 2017-07-01 17 M 200 100 100 16 10 5 (80, 81)
120 CS51 P51CS CS 2017-07-01 17 M 200 100 100 15 8 5 (80, 81)
122 CS53 P53CS CS 2017-11-01 60 M 200 100 100 27 15 5 (80, 81)
123 CS53 P53CS CS 2017-11-01 60 M 197 97 100 6 4 5 (80, 81)
124 CS54 P54CS CS 2018-01-01 59 F 200 100 100 43 17 5 (80, 81)
125 CS54 P54CS CS 2018-01-01 59 F 200 100 100 46 24 5 (80, 81)
126 CS54 P54CS CS 2018-01-01 59 F 200 100 100 34 15 5 (80, 81)
127 CS55 P55CS CS 2018-03-01 43 F 200 100 100 8 6 5 (80, 81)
128 CS56 P56CS CS 2018-04-01 48 M 200 100 100 7 7 5 (80, 81)
130 CS57 P57CS CS 2018-06-01 46 M 200 100 100 16 11 5 (80, 81)
131 CS57 P57CS CS 2018-06-01 46 M 200 100 100 11 9 5 (80, 81)
132 CS58 P58CS CS 2018-07-01 32 F 200 100 100 8 8 5 (80, 81)
133 CS58 P58CS CS 2018-07-01 32 F 200 100 100 2 2 5 (80, 81)
134 CS60 P60CS CS 2018-10-01 67 M 200 100 100 14 8 5 (80, 81)
135 CS61 P61CS CS 2019-02-01 52 F 200 100 100 13 12 5 (80, 81)
136 CS62 P62CS CS 2019-04-01 25 F 200 100 100 13 9 5 (80, 81)
2001 T88 TWH088 TWH 2018-03-01 26 M 200 100 100 9 7 5 (80, 81)
2002 T88 TWH088 TWH 2018-03-01 26 M 200 100 100 4 4 5 (80, 81)
2003 T90 TWH090 TWH 2018-04-01 20 M 200 100 100 5 4 5 (80, 81)
2004 T89 TWH089 TWH 2018-03-01 45 F 200 100 100 31 21 5 (80, 81)
2005 T98 TWH098 TWH 2018-09-01 46 F 200 100 100 13 11 5 (80, 81)
2007 T107 TWH107 TWH 2018-12-01 64 F 200 100 100 3 3 5 (80, 81)
2008 T100 TWH100 TWH 2018-10-01 19 F 200 100 100 1 1 5 (80, 81)
2009 T103 TWH103 TWH 2018-11-01 49 M 200 100 100 3 2 5 (80, 81)
2011 T101 TWH101 TWH 2018-10-01 25 F 200 100 100 8 8 5 (80, 81)
The archive's summary gives 59 subjects and 87 files, which the files confirm, and the two published versions hold the same assets. There are 1864 unit rows, one more than the paper's 1863. Eight files have fewer than 100 learning trials (seven with 50 and one with 97). The institution is written "Hunigton Memorial Hospital" in 28 files and "Hunington Memorial Hospital" in 6, and the device both "Neuralynx-cheetah" and "Neuralynx-Cheetah". There is no session 144, the session of the example neuron of Figure 3.
In [3]:
# Units and coding of the variables used: times, responses, labels, categories.
rows = []
for s in sessions:
t = s["trials"]
r = rm.recognition(s)
ev, et = s["events"], s["events_t"]
stim_on = et[ev == 1]
on_idx = np.flatnonzero(ev == 1)
resp_ev = np.array([ev[i + 3] if i + 3 < len(ev) else np.nan for i in on_idx])
sp = np.concatenate(s["spikes"])
rows.append({
"session": s["session"],
"stim_s_median": np.median(t.stim_off_time - t.stim_on_time),
"start_is_stim_on": bool(np.allclose(t.start_time, t.stim_on_time)),
"onsets_match_events": bool(len(stim_on) == len(t) and np.allclose(np.sort(stim_on), np.sort(t.stim_on_time.values))),
"resp_mismatch": int(np.sum(resp_ev != np.where(t.stim_phase == "recog", t.response_value, t.response_value + 20)))
if len(resp_ev) == len(t) else -1,
"unlisted_after_question": sorted({int(x) for x in resp_ev} - set(range(20, 22)) - set(range(31, 37))),
"resp_min": int(r.response.min()), "resp_max": int(r.response.max()),
"labels": tuple(sorted(r.label.unique())),
"spikes_from_s": sp.min() - t.start_time.min(), "spikes_to_s": sp.max() - t.stop_time.max(),
"trial_s_median": np.median(t.stop_time - t.start_time)})
units = pd.DataFrame(rows)
print("Stimulus duration (s), median per file: range", units.stim_s_median.min().round(3), "to", units.stim_s_median.max().round(3))
print("Trial duration, stimulus onset to end of trial (s), median per file:", units.trial_s_median.min().round(2), "to", units.trial_s_median.max().round(2))
print("start_time equals stim_on_time in every file:", units.start_is_stim_on.all())
print("Trial onsets equal the stimulus-onset events (TTL 1) in every file:", units.onsets_match_events.all())
print("Per trial, the event three after stimulus onset (the response the paper's MATLAB code reads) against response_value")
print("(recognition: 31-36; learning: response_value + 20 against 20/21): files with every trial equal:", int((units.resp_mismatch == 0).sum()), "of", len(units))
print(units[units.resp_mismatch != 0][["session", "resp_mismatch", "unlisted_after_question"]].to_string(index=False))
print("Recognition responses (response_value - 30) span", units.resp_min.min(), "to", units.resp_max.max(), "; labels:", sorted({x for l in units.labels for x in l}))
print("First spike relative to first trial onset (s): range", units.spikes_from_s.min().round(1), "to", units.spikes_from_s.max().round(1))
print("Last spike relative to last trial end (s): range", units.spikes_to_s.min().round(1), "to", units.spikes_to_s.max().round(1))
print("Files whose first spike is more than 60 s before the first trial:", len(units[units.spikes_from_s < -60]), units[units.spikes_from_s < -60].session.tolist())
print("\nDescription of new_old_labels_recog in the files:", repr(sessions[0]["trial_descriptions"]["new_old_labels_recog"]))
print("Learning response_value values:", sorted(set(np.concatenate([s["trials"][s["trials"].stim_phase == "learn"].response_value.unique() for s in sessions]))))
Stimulus duration (s), median per file: range 0.998 to 2.016
Trial duration, stimulus onset to end of trial (s), median per file: 2.87 to 5.81
start_time equals stim_on_time in every file: True
Trial onsets equal the stimulus-onset events (TTL 1) in every file: True
Per trial, the event three after stimulus onset (the response the paper's MATLAB code reads) against response_value
(recognition: 31-36; learning: response_value + 20 against 20/21): files with every trial equal: 81 of 87
session resp_mismatch unlisted_after_question
58 1 [23]
123 18 []
127 1 [3]
2001 6 [3]
2002 9 [3]
2003 4 [3]
Recognition responses (response_value - 30) span 1 to 6 ; labels: [np.int64(0), np.int64(1)]
First spike relative to first trial onset (s): range -2462.3 to 17.4
Last spike relative to last trial end (s): range -46.6 to 2055.1
Files whose first spike is more than 60 s before the first trial: 15 [6, 7, 9, 17, 18, 24, 26, 28, 39, 41, 43, 47, 52, 63, 85]
Description of new_old_labels_recog in the files: "The Ground truth Labels for New or Old Stimulus. 0 == Old Stimuli \n (presented during the learning phase), 1 = New Stimuli (not seen )'during learning phase"
Learning response_value values: [np.float64(0.0), np.float64(1.0)]
Times are in seconds on one clock per file: trial start times equal the stimulus-onset events and stimuli last a median of 1.0 to 2.0 s per file (the Methods give no duration). In 81 of 87 files the trials table's response_value equals, trial by trial, the response event the paper's MATLAB code reads; in the six others 1 to 18 trials differ, where extra codes (3 or 23) or a shifted sequence follow the question. The analysis takes responses from response_value minus 30, which spans 1 to 6. In learning trials the table stores 0 and 1 for the events 20 and 21. The column description of new_old_labels_recog says 0 is old; this is tested under new-old-label-coding. In 15 files spikes start more than 60 s before the first trial (up to 2462 s), so the units cover more than the task.
In [4]:
# Recordings deposited more than once: identical spike trains, trial times or event sequences.
pairs = []
for a, b in combinations(sessions, 2):
same_spikes = len(a["spikes"]) == len(b["spikes"]) and all(
len(x) == len(y) and np.array_equal(x, y) for x, y in zip(a["spikes"], b["spikes"]))
ta, tb = a["trials"], b["trials"]
n = min(len(ta), len(tb))
shared_trials = int(np.sum(np.isin(np.round(ta.start_time.values, 4), np.round(tb.start_time.values, 4))))
if same_spikes or shared_trials >= 10:
pairs.append({"a": a["session"], "b": b["session"], "units": len(a["spikes"]), "identical_spike_trains": same_spikes,
"shared_trial_onsets": shared_trials, "trials_a": len(ta), "trials_b": len(tb),
"exp_ids_a": tuple(int(x) for x in np.unique(a["experiment_ids"])),
"exp_ids_b": tuple(int(x) for x in np.unique(b["experiment_ids"])),
"a_trials_s": (round(ta.start_time.min()), round(ta.stop_time.max())),
"b_trials_s": (round(tb.start_time.min()), round(tb.stop_time.max())),
"spikes_s": (round(np.concatenate(a["spikes"]).min()), round(np.concatenate(a["spikes"]).max()))})
pairs = pd.DataFrame(pairs)
print(pairs.to_string(index=False))
DUPLICATE_SECOND = sorted(pairs[pairs.identical_spike_trains].b.tolist())
print("\nSessions whose units repeat another session's (dropped when each recording is counted once):", DUPLICATE_SECOND)
print("Units in them:", sum(len(by_id[k]["spikes"]) for k in DUPLICATE_SECOND))
same_pair = pairs[(pairs.a == 54) & (pairs.b == 55)]
r54, r55 = rm.recognition(by_id[54]), rm.recognition(by_id[55])
l54 = by_id[54]["trials"][by_id[54]["trials"].stim_phase == "learn"]
l55 = by_id[55]["trials"][by_id[55]["trials"].stim_phase == "learn"]
print("Sessions 54 and 55: learning trials identical:", l54.reset_index(drop=True).equals(l55.reset_index(drop=True)),
"; recognition onsets shared:", int(np.isin(r54.start_time.round(4), r55.start_time.round(4)).sum()),
"; recognition images shared:", len(set(r54.external_image_file) & set(r55.external_image_file)))
a b units identical_spike_trains shared_trial_onsets trials_a trials_b exp_ids_a exp_ids_b a_trials_s b_trials_s spikes_s 43 44 23 True 0 200 150 (88, 89) (83, 84) (3292, 4478) (1495, 3035) (1494, 4478) 54 55 32 True 100 200 200 (80, 81) (80, 82) (1611, 3626) (1611, 4525) (1610, 4525) Sessions whose units repeat another session's (dropped when each recording is counted once): [44, 55] Units in them: 55 Sessions 54 and 55: learning trials identical: True ; recognition onsets shared: 0 ; recognition images shared: 0
Two pairs of files hold identical spike trains for every unit. Sessions 43 and 44 (subject H23) are two different task blocks of one continuous recording: the spike trains span both blocks, session 44's trials come first and session 43's later. Sessions 54 and 55 (subject H42) share the same learning block, trial for trial, and differ in their recognition blocks (experiment ids 81 and 82), so session 55 is a second recognition test of the same encoding block, deposited with a copy of the learning trials and the same spike trains. In both pairs the units are the same neurons. The paper counts neurons per session, so its 1863 includes them twice; results are also given with sessions 44 and 55 dropped.
In [5]:
# Clock check: category selectivity of the units in 1 s windows shifted around stimulus onset, with
# the true labels and with the labels of the previous trial. A response to the stimulus can only
# appear after onset and with the true labels; an offset between spikes and trials would move it.
offsets = np.arange(-1.6, 1.61, 0.2)
frac_true, frac_shift = [], []
per_unit = []
for s in sessions:
r = rm.recognition(s)
on, cat = r.start_time.values, r.category.values
for k, sp in enumerate(s["spikes"]):
p_true = [rm.anova_p(rm.window_rates(sp, on, (o, o + 1.0)), cat) for o in offsets]
p_shift = [rm.anova_p(rm.window_rates(sp, on, (o, o + 1.0)), np.roll(cat, 1)) for o in offsets]
per_unit.append({"session": s["session"], "p_true": p_true, "p_shift": p_shift,
"p_pre": rm.anova_p(rm.window_rates(sp, on, (-1.0, 0.0)), cat),
"p_post": rm.anova_p(rm.window_rates(sp, on, rm.WINDOW_METHODS), cat)})
P = np.array([u["p_true"] for u in per_unit]); PS = np.array([u["p_shift"] for u in per_unit])
clock = pd.DataFrame({"window_start_s": offsets.round(1), "frac_p<0.05_true_labels": (P < 0.05).mean(0).round(3),
"frac_p<0.05_previous_trial_labels": (PS < 0.05).mean(0).round(3)})
print(clock.to_string(index=False))
pu = pd.DataFrame(per_unit)
ses = pu.groupby("session").agg(units=("p_post", "size"), sig_pre=("p_pre", lambda p: int((p < 0.05).sum())),
sig_post=("p_post", lambda p: int((p < 0.05).sum())))
ses["binom_p_post_vs_5pct"] = [stats.binomtest(k, n, 0.05, alternative="greater").pvalue for k, n in zip(ses.sig_post, ses.units)]
ses["established"] = (ses.binom_p_post_vs_5pct < 0.05) & (ses.sig_post > ses.sig_pre)
print("\nSessions where the stimulus response (more category-selective units after onset than expected by chance,")
print("and more than before onset) establishes the shared clock:", int(ses.established.sum()), "of", len(ses))
print("Not established (too few units or no excess):", ses[~ses.established].index.tolist())
print("Units in sessions not established:", int(ses[~ses.established].units.sum()))
fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(offsets, (P < 0.05).mean(0), "o-", label="true labels")
ax.plot(offsets, (PS < 0.05).mean(0), "s-", label="previous trial's labels")
ax.axhline(0.05, color="k", lw=0.8)
ax.set_xlabel("window start relative to stimulus onset (s), 1 s windows")
ax.set_ylabel("fraction of units p < 0.05")
ax.legend(); plt.tight_layout(); plt.show()
lo_, hi_ = min((P < 0.05).mean(0).min(), (PS < 0.05).mean(0).min()), max((P < 0.05).mean(0).max(), (PS < 0.05).mean(0).max())
print(f"Plotted fractions range from {lo_:.3f} to {hi_:.3f}; none outside the axes.")
window_start_s frac_p<0.05_true_labels frac_p<0.05_previous_trial_labels
-1.6 0.047 0.055
-1.4 0.044 0.051
-1.2 0.041 0.055
-1.0 0.049 0.060
-0.8 0.056 0.056
-0.6 0.064 0.055
-0.4 0.101 0.048
-0.2 0.127 0.047
-0.0 0.164 0.056
0.2 0.178 0.047
0.4 0.167 0.054
0.6 0.141 0.048
0.8 0.119 0.056
1.0 0.106 0.054
1.2 0.085 0.050
1.4 0.079 0.047
1.6 0.070 0.039
Sessions where the stimulus response (more category-selective units after onset than expected by chance,
and more than before onset) establishes the shared clock: 38 of 87
Not established (too few units or no excess): [5, 6, 9, 17, 18, 20, 21, 24, 26, 27, 38, 39, 47, 50, 52, 55, 56, 61, 63, 64, 66, 67, 68, 78, 85, 92, 96, 97, 101, 111, 113, 114, 116, 117, 118, 120, 122, 125, 127, 128, 133, 135, 136, 2001, 2002, 2003, 2005, 2007, 2008]
Units in sessions not established: 801
Plotted fractions range from 0.039 to 0.178; none outside the axes.
With the true labels, the fraction of category-selective units is 4.1% to 4.9% for windows that end at or before stimulus onset (start at -1.0 s or earlier), rises once the window reaches past 200 ms after onset, peaks at 17.8% for the window starting at 0.2 s and declines after. With the labels of the previous trial it stays between 3.9% and 6.0% at every offset. Spikes and trials therefore share a clock overall. Per session this is shown for 38 of 87 sessions; the 49 others, listed above with 801 units, have too few units or too few selective ones to tell, and for them the shared clock is consistent with the data but assumed.
counts-neurons-subjects¶
The dataset holds 1863 single neurons recorded from the medial temporal lobes of 59 subjects.
In [6]:
n_units = int(info.units.sum())
n_units_once = int(info[~info.session.isin(DUPLICATE_SECOND)].units.sum())
n_subjects = info.subject.nunique()
print(f"(a) every unit row: {n_units} neurons (paper 1863, tolerance 1676.7 to 2049.3; difference {n_units - 1863:+d})")
print(f"(b) duplicated spike trains counted once: {n_units_once} neurons (difference {n_units_once - 1863:+d})")
print(f"subjects: {n_subjects} (paper 59)")
print("units by site:", info.groupby("site").units.sum().to_dict())
print("\nunits per site:", info.groupby("site").units.sum().to_dict())
print("units per session:")
print(info.set_index("session").units.to_frame().T.to_string())
fewest = sorted(((len(sp), s["session"], k) for s in sessions for k, sp in enumerate(s["spikes"])))[:5]
print("five units with the fewest spikes (spikes, session, unit):", fewest)
(a) every unit row: 1864 neurons (paper 1863, tolerance 1676.7 to 2049.3; difference +1)
(b) duplicated spike trains counted once: 1809 neurons (difference -54)
subjects: 59 (paper 59)
units by site: {'CS': 1787, 'TWH': 77}
units per site: {'CS': 1787, 'TWH': 77}
units per session:
session 5 6 7 9 17 18 20 21 23 24 26 27 28 32 38 39 41 43 44 47 48 49 50 52 54 55 56 58 59 60 61 63 64 66 67 68 69 70 73 74 76 77 78 85 92 93 96 97 99 100 101 102 104 105 111 112 113 114 115 116 117 118 119 120 122 123 124 125 126 127 128 130 131 132 133 134 135 136 2001 2002 2003 2004 2005 2007 2008 2009 2011
units 24 19 38 25 23 33 20 14 18 21 33 16 15 43 20 18 15 23 23 17 45 27 14 18 32 32 15 49 21 22 8 14 33 11 15 19 64 33 4 35 53 30 13 6 7 35 13 10 30 42 14 44 42 43 14 12 21 15 35 36 11 13 16 15 27 6 43 46 34 8 7 16 11 8 2 14 13 13 9 4 5 31 13 3 1 3 8
five units with the fewest spikes (spikes, session, unit): [(22, 28, 7), (30, 69, 12), (41, 23, 15), (46, 20, 0), (51, 113, 1)]
Result: reproduced. The files hold 1864 unit rows, one more than the paper's 1863, and 59 subjects, both within tolerance. Counting the 55 units of sessions 44 and 55 once, as their spike trains repeat those of sessions 43 and 54, gives 1809, also within tolerance. The files hold 1787 units in CS sessions and 77 in TWH sessions. The paper gives no count per session or site, so the extra unit cannot be identified; no unit is empty (the fewest spikes is 22).
counts-sessions¶
The patients performed 87 sessions, listed by subject and session id in Online-only Table 1 (the Methods give 89).
In [7]:
table = rm.online_table_1()
files_set = set(zip(info.label.str.replace("CS", "C").str.replace(r"^([A-Z])0*", r"\1", regex=True), info.session))
table_set = set(zip(table.table_subject.str.replace(r"^([A-Z])0*", r"\1", regex=True), table.session))
print("sessions in the files:", len(info), "; in Online-only Table 1:", len(table), "; Methods: 89")
print("(subject, session) pairs only in the files:", sorted(files_set - table_set))
print("(subject, session) pairs only in the table:", sorted(table_set - files_set))
print("Task variant against the stimulus set of the file:")
var = table.set_index("session").variant
print(pd.crosstab(info.session.map(var).rename("variant (table)"),
pd.Series([sorted(s["trials"].category_name.unique())[0] for s in sessions], name="first category name")).to_string())
sessions in the files: 87 ; in Online-only Table 1: 87 ; Methods: 89 (subject, session) pairs only in the files: [] (subject, session) pairs only in the table: [] Task variant against the stimulus set of the file: first category name 1cars fruit houses variant (table) 1 0 0 42 2 0 20 0 3 25 0 0
Result: reproduced. The 87 files match the 87 (subject, session) pairs of Online-only Table 1 exactly, with none missing or extra; the Methods' 89 sessions disagree with the paper's own table. Each table variant corresponds to one stimulus set in the files.
table1-demographics¶
Online-only Table 1 gives the age and sex of every subject; there are 25 female and 34 male subjects.
In [8]:
subj = info.groupby("subject").agg(label=("label", "first"), age=("age", "first"), sex=("sex", "first"),
ages=("age", "nunique"), sexes=("sex", "nunique")).reset_index()
subj["session"] = info.groupby("subject").session.first().values
m = subj.merge(table[["session", "table_subject", "age", "sex", "diagnosis"]], on="session", suffixes=("_file", "_table"))
m["age_file"] = m.age_file.astype(int)
m["match"] = (m.age_file == m.age_table) & (m.sex_file == m.sex_table)
diag = {s["subject_id"]: s["diagnosis"] for s in sessions}
m["diagnosis_file"] = m.subject.map(diag)
print("Subjects with one age and one sex across their files:", int(((subj.ages == 1) & (subj.sexes == 1)).sum()), "of", len(subj))
print("Subjects matching the table in age and sex:", int(m.match.sum()), "of", len(m))
print(m[~m.match][["subject", "table_subject", "age_file", "age_table", "sex_file", "sex_table"]].to_string(index=False))
print("\nFemale / male in the files:", subj.sex.value_counts().to_dict(), "; in the table: F 25, M 34")
print("\nDiagnosis, files against table, where the text differs:")
diag_diff = m[m.diagnosis_file.str.lower() != m.diagnosis.str.lower()]
print("subjects whose diagnosis text differs:", len(diag_diff), "of", len(m))
print(diag_diff[["table_subject", "diagnosis_file", "diagnosis"]].to_string(index=False))
code_list = rm.code_session_list()
print(f"\nThe code's session list (defineNOsessions_release.ini at {rm.CODE_COMMIT[:7]}): {int(code_list.active.sum())} active sections,",
f"commented out: {code_list[~code_list.active].session.tolist()}")
cl = code_list[code_list.active].assign(age=lambda d: d.age.astype(int))
fs = info.merge(cl[["session", "ID", "age", "sex"]], on="session", how="left", suffixes=("_file", "_code"))
print("files with a section in the code's list:", int(fs.ID.notna().sum()), "of", len(fs),
"; with the same age and sex as the file:", int(((fs.age_file.astype(int) == fs.age_code) & (fs.sex_file == fs.sex_code)).sum()))
print("the three subjects that differ from Online-only Table 1, in the code's list:")
print(fs[fs.session.isin(m[~m.match].session)][["session", "label", "ID", "age_file", "age_code", "sex_file", "sex_code"]].to_string(index=False))
Subjects with one age and one sex across their files: 59 of 59
Subjects matching the table in age and sex: 56 of 59
subject table_subject age_file age_table sex_file sex_table
P24CS C24 57 47 F F
P33HMH H33 30 29 M M
P48HMH H48 32 54 F M
Female / male in the files: {'M': 33, 'F': 26} ; in the table: F 25, M 34
Diagnosis, files against table, where the text differs:
subjects whose diagnosis text differs: 19 of 59
table_subject diagnosis_file diagnosis
H10 Left Lateral Frontal Left Frontal
H16 Right Lateral Frontal Right Frontal
H17 Left Lateral Frontal Left Inferior Frontal
H19 Left Lateral Frontal Left Frontal
C25 Bilateral Independent Temporal Bilateral Independent Mesial Temporal
H29 Left Other Left Deep Insula
C32 Not Localized Not Localized (Generalized)
C34 Bilateral Independent Temporal Bilateral Mesial Temporal
C39 Right Other Right Insula
C40 Right Other Right Motor Cortex
C51 Not Localized Not Localized (No Seizures)
H51 Bilateral Independent Frontal Bilateral Frontal and Temporal
C53 Bilateral Independent Temporal Bilateral Independent Mesial Temporal Lobe
C56 Bilateral Independent Temporal Multifocal, Bitemporal + Orbitofrontal
C57 Left Other Left Neocortical Parietal
C58 Right Lateral Frontal Right Frontal Neocortical
C61 Left Lateral Frontal Left Mesial Temporal
T89 Right Lateral Frontal Right Frontal
T90 Bilateral Occipital Cortex Occipital Cortex
The code's session list (defineNOsessions_release.ini at 1d95162): 87 active sections,
commented out: [98]
files with a section in the code's list: 87 of 87 ; with the same age and sex as the file: 87
the three subjects that differ from Online-only Table 1, in the code's list:
session label ID age_file age_code sex_file sex_code
52 H33 H33 30 30 M M
59 CS24 CS24 57 57 F F
99 H48 H48 32 32 F F
Result: partially reproduced. Age and sex match the table for 56 of 59 subjects. C24 is 57 in the files and 47 in the table, H33 is 30 and 29, and H48 is a 32-year-old woman in the files and a 54-year-old man in the table. The files hold 26 female and 33 male subjects, within tolerance of the table's 25 and 34. The diagnosis text differs for 19 of 59 subjects. For eight of them (C25, H29, C32, C39, C40, C51, C53, C57) the files give a broader category, such as Right Other for Right Insula. For six (H10, H16, H17, H19, C58, T89) the files give Lateral Frontal, which adds Lateral to the table's Frontal, except for H17, whose table entry is Left Inferior Frontal. The five others differ in substance: C61 (Left Lateral Frontal against Left Mesial Temporal), C34 (Bilateral Independent Temporal against Bilateral Mesial Temporal), C56 and H51, where the files leave out part of the table's diagnosis (orbitofrontal, temporal), and T90, where they add "Bilateral". The code's session list has a section for each of the 87 files with the same age and sex as the file, so for C24, H33 and H48 the files follow the code and the table differs from both.
task-structure¶
Each session has 100 encoding trials showing 100 novel images and 100 recognition trials, of which 50 show new and 50 old images, and the images come from 5 visual categories.
In [9]:
rows = []
for s in sessions:
t = s["trials"]; learn = t[t.stim_phase == "learn"]; r = rm.recognition(s)
(shape, img_t), = s["stimulus"].values()
rows.append({"session": s["session"], "learn": len(learn), "learn_distinct_images": learn.external_image_file.nunique(),
"recog": len(r), "label_counts": tuple(r.label.value_counts().sort_index()), "categories": t.stimCategory.nunique(),
"category_names": t.category_name.nunique(), "images_stored": shape[0],
"image_times_are_onsets": bool(np.allclose(np.sort(img_t), np.sort(t.start_time.values)))})
task = pd.DataFrame(rows)
ok_learn = (task.learn == 100) & (task.learn_distinct_images == 100)
ok_recog = (task.recog == 100) & (task.label_counts == (50, 50))
ok_cat = (task.categories == 5) & (task.category_names == 5)
print("100 encoding trials with 100 distinct images:", int(ok_learn.sum()), "of", len(task))
print("100 recognition trials, 50 per label:", int(ok_recog.sum()), "of", len(task))
print("5 categories (ids and names):", int(ok_cat.sum()), "of", len(task))
print("One stored image per trial, timestamped at trial onset:", int(((task.images_stored == task.learn + task.recog) & task.image_times_are_onsets).sum()), "of", len(task))
print("\nSessions with other encoding blocks:")
print(task[~ok_learn].to_string(index=False))
print("\nStimulus storage: one ImageSeries per file,", sorted({k for s in sessions for k in s["stimulus"]}),
"(the paper describes stimuli_learn_x and stimuli_recog_x entries)")
100 encoding trials with 100 distinct images: 79 of 87
100 recognition trials, 50 per label: 87 of 87
5 categories (ids and names): 87 of 87
One stored image per trial, timestamped at trial onset: 87 of 87
Sessions with other encoding blocks:
session learn learn_distinct_images recog label_counts categories category_names images_stored image_times_are_onsets
18 50 50 100 (50, 50) 5 5 150 True
20 50 50 100 (50, 50) 5 5 150 True
21 50 50 100 (50, 50) 5 5 150 True
23 50 50 100 (50, 50) 5 5 150 True
26 50 50 100 (50, 50) 5 5 150 True
27 50 50 100 (50, 50) 5 5 150 True
44 50 50 100 (50, 50) 5 5 150 True
123 97 97 100 (50, 50) 5 5 197 True
Stimulus storage: one ImageSeries per file, ['StimulusPresentation'] (the paper describes stimuli_learn_x and stimuli_recog_x entries)
Result: partially reproduced. All 87 sessions have 100 recognition trials with 50 of each label and 5 categories, and one stored image per trial. The encoding block has 100 trials with 100 distinct images in 79 sessions; seven sessions have 50 encoding trials and session 123 has 97. The images are stored as one image series per file, not as the stimuli_learn_x and stimuli_recog_x entries the paper describes. In the eight sessions with shorter encoding blocks, 2% to 48% of the images shown as old in recognition have no encoding trial in the file (shown under new-old-label-coding), so the file does not hold the encoding of every old image.
new-old-label-coding¶
In the trials table, new_old_labels_recog is 0 for old stimuli (seen during encoding) and 1 for new stimuli.
In [10]:
rows = []
for s in sessions:
t = s["trials"]; r = rm.recognition(s)
learned = set(t[t.stim_phase == "learn"].external_image_file)
seen = r.external_image_file.isin(learned)
rows.append({"session": s["session"],
"label0_seen_in_encoding": seen[r.label == 0].mean(), "label1_seen_in_encoding": seen[r.label == 1].mean(),
"label0_old_responses": (r.response[r.label == 0] >= 4).mean(), "label1_old_responses": (r.response[r.label == 1] >= 4).mean()})
lab = pd.DataFrame(rows)
print("Fraction of recognition trials whose image file was shown during encoding, median over sessions:")
print(" label 0:", lab.label0_seen_in_encoding.median(), " label 1:", lab.label1_seen_in_encoding.median())
print("Sessions where every label-1 image and no label-0 image was shown in encoding:",
int(((lab.label1_seen_in_encoding == 1) & (lab.label0_seen_in_encoding == 0)).sum()), "of", len(lab))
print("Sessions where label 0 matches the stated coding (all label-0 images seen, no label-1 image seen):",
int(((lab.label0_seen_in_encoding == 1) & (lab.label1_seen_in_encoding == 0)).sum()))
print("\nFraction of 'old' responses (4-6), median over sessions: label 0", round(lab.label0_old_responses.median(), 3),
" label 1", round(lab.label1_old_responses.median(), 3))
print("Sessions with more 'old' responses to label 1 than to label 0:", int((lab.label1_old_responses > lab.label0_old_responses).sum()), "of", len(lab))
print("Sign test, label 1 above label 0: p =", stats.binomtest(int((lab.label1_old_responses > lab.label0_old_responses).sum()),
int((lab.label1_old_responses != lab.label0_old_responses).sum())).pvalue)
print("Sessions with imperfect image match:")
print(lab[(lab.label1_seen_in_encoding < 1) | (lab.label0_seen_in_encoding > 0)].round(3).to_string(index=False))
Fraction of recognition trials whose image file was shown during encoding, median over sessions:
label 0: 0.0 label 1: 1.0
Sessions where every label-1 image and no label-0 image was shown in encoding: 79 of 87
Sessions where label 0 matches the stated coding (all label-0 images seen, no label-1 image seen): 0
Fraction of 'old' responses (4-6), median over sessions: label 0 0.24 label 1 0.7
Sessions with more 'old' responses to label 1 than to label 0: 86 of 87
Sign test, label 1 above label 0: p = 1.137373342260413e-24
Sessions with imperfect image match:
session label0_seen_in_encoding label1_seen_in_encoding label0_old_responses label1_old_responses
18 0.0 0.58 0.02 0.44
20 0.0 0.58 0.76 0.88
21 0.0 0.52 0.70 0.80
23 0.0 0.52 0.06 0.58
26 0.0 0.58 0.14 0.44
27 0.0 0.52 0.30 0.82
44 0.0 0.52 0.06 0.42
123 0.0 0.98 0.14 0.60
Result: not reproduced. The coding is the reverse of the stated one. In 79 of 87 sessions every label-1 image and no label-0 image was shown during encoding, and in no session does label 0 mark the encoded images. Label-1 trials draw "old" responses more often than label-0 trials in 86 of 87 sessions (median 0.70 against 0.24, sign test p = 1.1e-24). Label 1 is old, as the paper's code reads it; the paper's text and the files' column description say 0. In the eight sessions with shortened encoding blocks only 52% to 98% of the label-1 images are among the stored encoding trials, so encoding trials are missing from those files.
ttl-codes¶
The events stream holds the TTL markers 55, 1, 2, 3, 20 or 21, 31 to 36, 6 and 66 with the meanings of Table 1, and the events and experiment_ids streams have one entry each, in the same order.
In [11]:
LISTED = {55, 1, 2, 3, 20, 21, 31, 32, 33, 34, 35, 36, 6, 66}
rows = []
for s in sessions:
ev, et = s["events"], s["events_t"]
same = len(ev) == len(s["experiment_ids"]) and np.allclose(et, s["experiment_ids_t"])
other = sorted(set(ev.astype(int)) - LISTED)
on = np.flatnonzero(ev == 1)
t = s["trials"]
phase = t.stim_phase.values[np.abs(t.start_time.values[None, :] - et[on][:, None]).argmin(1)]
seq_ok = phase_ok = 0
for i, ph in zip(on, phase):
nxt = ev[i + 1:i + 5]
if len(nxt) == 4 and nxt[0] == 2 and nxt[1] == 3 and (nxt[2] in (20, 21) or 31 <= nxt[2] <= 36) and nxt[3] == 6:
seq_ok += 1
phase_ok += (nxt[2] in (20, 21)) if ph == "learn" else (31 <= nxt[2] <= 36)
rows.append({"session": s["session"], "events": len(ev), "same_length_and_times": same, "unlisted_codes": tuple(other),
"trials": len(on), "trials_in_order": seq_ok, "response_kind_of_phase": phase_ok,
"starts_55": int((ev == 55).sum()), "ends_66": int((ev == 66).sum())})
ttl = pd.DataFrame(rows)
print("events and experiment_ids of equal length and timestamps:", int(ttl.same_length_and_times.sum()), "of", len(ttl))
print("files with codes not in Table 1:", ttl[ttl.unlisted_codes != ()][["session", "unlisted_codes"]].to_string(index=False))
print("trials showing 1, 2, 3, response, 6 in order:", int(ttl.trials_in_order.sum()), "of", int(ttl.trials.sum()))
print("trials in order whose response is of the phase's kind (20 or 21 in learning, 31 to 36 in recognition):",
int(ttl.response_kind_of_phase.sum()), "of", int(ttl.trials_in_order.sum()))
print("(start 55, end 66) markers per file: number of files", ttl.groupby(["starts_55", "ends_66"]).size().to_dict())
print("files without a 55:", ttl[ttl.starts_55 == 0].session.tolist(), "; files with other than two 66:", ttl[ttl.ends_66 != 2].session.tolist())
ttl["all_ok"] = (ttl.same_length_and_times & (ttl.unlisted_codes == ()) & (ttl.trials_in_order == ttl.trials)
& (ttl.response_kind_of_phase == ttl.trials_in_order))
print("files meeting every condition:", int(ttl.all_ok.sum()), "of", len(ttl))
print(ttl[~ttl.all_ok].to_string(index=False))
events and experiment_ids of equal length and timestamps: 87 of 87
files with codes not in Table 1: session unlisted_codes
43 (98, 99)
58 (23,)
93 (99,)
123 (0,)
127 (0,)
128 (4,)
2001 (22, 23)
2002 (23,)
2003 (23,)
trials showing 1, 2, 3, response, 6 in order: 17019 of 17047
trials in order whose response is of the phase's kind (20 or 21 in learning, 31 to 36 in recognition): 17019 of 17019
(start 55, end 66) markers per file: number of files {(0, 2): 7, (1, 2): 1, (2, 2): 78, (2, 3): 1}
files without a 55: [92, 93, 96, 97, 99, 100, 101] ; files with other than two 66: [2001]
files meeting every condition: 78 of 87
session events same_length_and_times unlisted_codes trials trials_in_order response_kind_of_phase starts_55 ends_66 all_ok
43 1006 True (98, 99) 200 200 200 2 2 False
58 1005 True (23,) 200 199 199 2 2 False
93 1004 True (99,) 200 200 200 0 2 False
123 991 True (0,) 197 196 196 1 2 False
127 1005 True (0,) 200 199 199 2 2 False
128 1005 True (4,) 200 199 199 2 2 False
2001 1013 True (22, 23) 200 193 193 2 3 False
2002 1014 True (23,) 200 190 190 2 2 False
2003 1011 True (23,) 200 193 193 2 2 False
Result: partially reproduced. In all 87 files the events and experiment_ids streams have the same length and timestamps, and 17019 of 17047 trials show onset, offset, question, response and end in order. 78 of 87 files meet every condition; the nine others hold codes Table 1 does not list (0, 4, 22, 23, 98 or 99). The start and end markers are not part of the condition: seven files (sessions 92, 93, 96, 97, 99, 100 and 101) have no start-of-experiment code 55, six of which meet every condition, and session 2001 has three end codes 66. In every trial that shows the sequence, the response is of its phase's kind (20 or 21 in learning, 31 to 36 in recognition).
experiment-ids¶
The learning block is labelled with one of 80, 83 or 88 and the recognition block with one of 81, 84 or 89.
In [12]:
rows = []
for s in sessions:
ids = s["experiment_ids"]; ev = s["events"]
first_resp = {int(e): int(ev[(ids == e) & (((ev >= 20) & (ev <= 21)) | ((ev >= 31) & (ev <= 36)))][0]) for e in np.unique(ids)}
learn_id = [e for e, r in first_resp.items() if r in (20, 21)]
recog_id = [e for e, r in first_resp.items() if 31 <= r <= 36]
rows.append({"session": s["session"], "learn_id": learn_id[0] if learn_id else None,
"recog_id": recog_id[0] if recog_id else None, "data_collection": s["data_collection"]})
eid = pd.DataFrame(rows)
eid["ok"] = eid.learn_id.isin([80, 83, 88]) & eid.recog_id.isin([81, 84, 89])
print("Block identified by its responses (20/21 learning, 31-36 recognition).")
print("files whose ids are in the stated sets:", int(eid.ok.sum()), "of", len(eid))
print(eid[~eid.ok].to_string(index=False))
print("\nid pairs:", eid.groupby(["learn_id", "recog_id"]).size().to_dict())
Block identified by its responses (20/21 learning, 31-36 recognition).
files whose ids are in the stated sets: 83 of 87
session learn_id recog_id data_collection ok
55 80 82 learning: 80, recognition: 82 False
77 81 82 learning: 81, recognition: 82 False
100 82 83 learning: 82, recognition: 83 False
105 82 83 learning: 82, recognition: 83 False
id pairs: {(80, 81): 65, (80, 82): 1, (81, 82): 1, (82, 83): 2, (83, 84): 14, (88, 89): 4}
Result: partially reproduced. 83 of 87 files use the ids of Table 2. Session 55 labels its recognition block 82, session 77 uses 81 for learning and 82 for recognition, and sessions 100 and 105 use 82 for learning and 83 for recognition.
session-dates¶
session_start_time keeps only the year and month of the experiment, with the day set to the first of the month for all sessions.
In [13]:
dates = pd.DataFrame({"session": [s["session"] for s in sessions], "start": [s["session_start_time"] for s in sessions]})
dates["day"] = dates.start.str[8:10]; dates["year"] = dates.start.str[:4].astype(int)
print("files on the first of a month:", int((dates.day == "01").sum()), "of", len(dates))
print(dates[dates.day != "01"].to_string(index=False))
print("years:", dates.year.min(), "to", dates.year.max(), "; by year:", dates.year.value_counts().sort_index().to_dict())
files on the first of a month: 86 of 87
session start day year
48 2010-03-13T00:00:00-08:00 13 2010
years: 2000 to 2019 ; by year: {2000: 1, 2006: 4, 2007: 6, 2008: 4, 2009: 6, 2010: 3, 2011: 2, 2012: 6, 2013: 6, 2014: 9, 2015: 4, 2016: 7, 2017: 8, 2018: 19, 2019: 2}
Result: partially reproduced. 86 of 87 files start on the first of a month. Session 48 starts on 2010-03-13, a day that was not set to the first. Session 68 is dated 2000-01-01, outside the span of the other recordings (2006 to 2019).
units-table¶
The units table holds spike_times, waveform_mean_encoding, waveform_mean_recognition (sampled at 100 kHz), origClusterID, SNR, IsolDist and electrodes, with one electrode per neuron.
In [14]:
needed = {"spike_times", "waveform_mean_encoding", "waveform_mean_recognition", "origClusterID", "SNR", "IsolationDist", "electrodes"}
cols_ok = sum(needed <= set(s["unit_columns"]) for s in sessions)
one_elec = sum(bool(np.array_equal(s["unit_electrode_index"], np.arange(1, len(s["spikes"]) + 1))) for s in sessions)
rates = np.unique(np.concatenate([s["waveform_rate"] for s in sessions]))
snr = np.concatenate([s["snr"] for s in sessions]); iso = np.concatenate([s["isolation_distance"] for s in sessions])
print("files with every column (IsolDist is named IsolationDist):", cols_ok, "of", len(sessions))
print("columns in every file:", sorted(set.intersection(*[set(s["unit_columns"]) for s in sessions])))
print("files where every unit links to exactly one electrode:", one_elec, "of", len(sessions))
print("stated waveform sampling rate (Hz):", rates, "; waveform samples:", sorted({s["waveform_shape"][1] for s in sessions}))
print(f"SNR: {np.isfinite(snr).sum()} of {len(snr)} finite; {np.isinf(snr).sum()} infinite; finite range {np.nanmin(snr[np.isfinite(snr)]):.2f} to {np.nanmax(snr[np.isfinite(snr)]):.2f}")
print(f"IsolationDist: {np.isfinite(iso).sum()} of {len(iso)} finite; {np.isnan(iso).sum()} NaN; finite range {np.nanmin(iso):.3g} to {np.nanmax(iso[np.isfinite(iso)]):.3g}; above 1e4: {(iso > 1e4).sum()}")
unsorted = [(s["session"], k) for s in sessions for k, sp in enumerate(s["spikes"]) if np.any(np.diff(sp) < 0)]
print("units whose spike times are not in increasing order:", len(unsorted), "in sessions", sorted({a for a, _ in unsorted}))
for a, k in unsorted:
sp = by_id[a]["spikes"][k]
n_dup = int(np.sum(np.diff(np.sort(sp)) == 0))
if n_dup:
print(f" session {a} unit {k}: {len(sp)} spike times, {n_dup} of them repeat another exactly")
files with every column (IsolDist is named IsolationDist): 87 of 87 columns in every file: ['IsolationDist', 'SNR', 'electrodes', 'electrodes_index', 'id', 'origClusterID', 'spike_times', 'spike_times_index', 'waveform_mean_encoding', 'waveform_mean_recognition', 'waveform_mean_sampling_rate'] files where every unit links to exactly one electrode: 87 of 87 stated waveform sampling rate (Hz): [98400.] ; waveform samples: [256] SNR: 1778 of 1864 finite; 11 infinite; finite range 0.38 to 11.89 IsolationDist: 1371 of 1864 finite; 493 NaN; finite range 4.51 to 8.08e+32; above 1e4: 4
units whose spike times are not in increasing order: 13 in sessions [21, 58, 68, 70, 73, 77, 128] session 58 unit 35: 3268 spike times, 1634 of them repeat another exactly session 58 unit 37: 824 spike times, 412 of them repeat another exactly
Result: partially reproduced. Every file has every column (IsolDist is named IsolationDist), and every unit links to one electrode. The stated waveform sampling rate is 98.4 kHz, within the tolerance of the paper's 100 kHz. The quality metrics are not finite for every unit: SNR is finite for 1778 of 1864 units, 11 being infinite, and the isolation distance is NaN for 493 units, with 4 values above 1e4 (up to 8.08e+32). The spike times of 13 units in 7 sessions are not in increasing order, and two units of session 58 hold every spike time twice.
electrodes-table¶
The electrodes table gives each electrode's brain area, MNI coordinates, impedance, the filter applied before spike sorting (300 to 3000 Hz) and the original hardware channel.
In [15]:
E = pd.concat([s["electrodes"].assign(session=s["session"]) for s in sessions], ignore_index=True)
print("electrodes:", len(E), "; columns in every file:", sorted(set.intersection(*[set(s["electrodes"].columns) for s in sessions])))
for c in ["location", "x", "y", "z", "imp", "filtering", "origChannel"]:
v = E[c]
empty = v.isna() | (v.astype(str).str.strip() == "")
print(f"{c:12s} empty for {int(empty.sum())} of {len(E)}; files with any empty: {E[empty].session.nunique()}")
print("filtering values:", E.filtering.value_counts().to_dict())
print("sessions with missing coordinates:", sorted(E[E.x.isna()].session.unique()))
print("coordinate ranges (mm): x", E.x.min(), "to", E.x.max(), "; y", E.y.min(), "to", E.y.max(), "; z", E.z.min(), "to", E.z.max())
side = E.dropna(subset=["x"])
print("left-named electrodes with x < 0:", int((side.location.str.startswith("Left") & (side.x < 0)).sum()), "of", int(side.location.str.startswith("Left").sum()),
"; right-named with x > 0:", int((side.location.str.startswith("Right") & (side.x > 0)).sum()), "of", int(side.location.str.startswith("Right").sum()))
print("imp: dtype", E.imp.dtype, "; distinct stored values:", [repr(v) for v in pd.unique(E.imp)][:10])
electrodes: 887 ; columns in every file: ['filtering', 'group_name', 'id', 'imp', 'location', 'origChannel', 'x', 'y', 'z']
location empty for 0 of 887; files with any empty: 0
x empty for 99 of 887; files with any empty: 8
y empty for 99 of 887; files with any empty: 8
z empty for 99 of 887; files with any empty: 8
imp empty for 887 of 887; files with any empty: 87
filtering empty for 0 of 887; files with any empty: 0
origChannel empty for 0 of 887; files with any empty: 0
filtering values: {'300-3000Hz': 887}
sessions with missing coordinates: [np.int64(7), np.int64(20), np.int64(21), np.int64(23), np.int64(24), np.int64(50), np.int64(52), np.int64(2004)]
coordinate ranges (mm): x -33.0 to 34.97 ; y -26.5 to 6.88 ; z -35.9736 to -7.74
left-named electrodes with x < 0: 364 of 364 ; right-named with x > 0: 424 of 424
imp: dtype float64 ; distinct stored values: ['np.float64(nan)']
Result: not reproduced. Location, filter (300-3000Hz) and original channel are filled for all 887 electrodes, and the coordinates match the hemisphere of the location wherever they are given. The impedance is stored as NaN (float64) for all 887 electrodes, and the coordinates are missing for 99 electrodes in 8 sessions.
isi-criterion¶
Clusters were accepted as units only when fewer than 3% of their interspike intervals were shorter than 3 ms.
In [16]:
isi = pd.DataFrame([{"session": s["session"], "unit": k, "spikes": len(sp), "isi_sorted": rm.isi_fraction(sp),
"isi_as_stored": rm.isi_fraction(sp, sort=False)} for s in sessions for k, sp in enumerate(s["spikes"])])
print("units:", len(isi))
for c in ["isi_sorted", "isi_as_stored"]:
print(f"{c}: below 3% {int((isi[c] < 0.03).sum())} ({100 * (isi[c] < 0.03).mean():.1f}%); median {100 * isi[c].median():.2f}%; max {100 * isi[c].max():.1f}%")
for thr in (0.03, 0.05, 0.10):
print(f"units at or above {thr:.0%} (sorted): {int((isi.isi_sorted >= thr).sum())}")
print(isi[isi.isi_sorted >= 0.03].round(4).to_string(index=False))
print("by site, fraction below 3%:", isi.assign(site=isi.session.map({s['session']: s['site'] for s in sessions})).groupby("site").isi_sorted.apply(lambda v: round((v < 0.03).mean(), 3)).to_dict())
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist(100 * isi.isi_sorted.clip(upper=0.2), bins=50)
ax.axvline(3, color="k"); ax.set_xlabel("% of ISIs < 3 ms (clipped at 20%)"); ax.set_ylabel("units")
plt.tight_layout(); plt.show()
print("units above 20% (clipped into the last bin):", int((isi.isi_sorted > 0.2).sum()))
units: 1864
isi_sorted: below 3% 1858 (99.7%); median 0.19%; max 50.1%
isi_as_stored: below 3% 1860 (99.8%); median 0.19%; max 9.0%
units at or above 3% (sorted): 6
units at or above 5% (sorted): 4
units at or above 10% (sorted): 2
session unit spikes isi_sorted isi_as_stored
58 35 3268 0.5002 0.0003
58 37 824 0.5006 0.0012
66 2 1176 0.0647 0.0647
124 2 4511 0.0302 0.0302
126 6 5530 0.0901 0.0901
126 7 2018 0.0377 0.0377
by site, fraction below 3%: {'CS': 0.997, 'TWH': 1.0}
units above 20% (clipped into the last bin): 2
Result: partially reproduced. 1858 of 1864 units (99.7%) have fewer than 3% of their interspike intervals below 3 ms, above the 95% the criterion sets for a partial result. Six units are at or above 3%, two of them above 20%. With the spike times as stored (unsorted in 13 units), 1860 units pass.
fig5a-mtl-locations¶
All neurons were recorded from the medial temporal lobe, in the amygdala and the hippocampus.
In [17]:
loc = pd.Series([s["electrodes"].location.values[e] for s in sessions for e in s["unit_electrode"]])
print("unit locations:", loc.value_counts().to_dict())
print("units in amygdala or hippocampus:", int(loc.str.contains("Amygdala|Hippocampus").sum()), "of", len(loc))
unit locations:
{'Right Amygdala': 594, 'Left Amygdala': 498, 'Right Hippocampus': 407, 'Left Hippocampus': 365}
units in amygdala or hippocampus: 1864 of 1864
Result: reproduced. All 1864 units are on electrodes located in the left or right amygdala or hippocampus.
fig5c-auc-all¶
Subjects distinguished new from old images; the behavioural AUC averaged over all sessions is 0.74 ± 0.10.
In [18]:
beh = []
for s in sessions:
r = rm.recognition(s)
lab_, resp = r.label.values, r.response.values
hi_f, lo_f = rm.confidence_accuracy(lab_, resp, "fixed")
hi_d, lo_d = rm.confidence_accuracy(lab_, resp, "dynamic")
beh.append({"session": s["session"], "site": s["site"], "auc": rm.auc(lab_, resp), "zroc_slope": rm.zroc_slope(lab_, resp), "zroc_degenerate": rm.zroc_degenerate(lab_, resp),
"acc_high_fixed": hi_f, "acc_low_fixed": lo_f, "acc_high_dynamic": hi_d, "acc_low_dynamic": lo_d,
"hit_rate": np.mean(resp[lab_ == rm.OLD] >= 4), "fa_rate": np.mean(resp[lab_ == rm.NEW] >= 4)})
beh = pd.DataFrame(beh)
print(f"AUC over {len(beh)} sessions: mean {beh.auc.mean():.3f}, SD {beh.auc.std(ddof=1):.3f} (paper 0.74 ± 0.10)")
print(f"range {beh.auc.min():.3f} to {beh.auc.max():.3f}; sessions below 0.6: {int((beh.auc < 0.6).sum())}; at or below 0.5: {int((beh.auc <= 0.5).sum())}")
print(f"chance AUC is 0.5; one-sample t-test of AUC against 0.5: t = {stats.ttest_1samp(beh.auc, 0.5).statistic:.1f}, p = {stats.ttest_1samp(beh.auc, 0.5).pvalue:.2g}")
print(f"hit rate {beh.hit_rate.mean():.3f}, false alarm rate {beh.fa_rate.mean():.3f} (means over sessions)")
once = beh[~beh.session.isin(DUPLICATE_SECOND)]
print(f"variant, sessions 44 and 55 dropped: mean {once.auc.mean():.3f}, SD {once.auc.std(ddof=1):.3f}")
beh["auc_stated_coding"] = [rm.auc(1 - rm.recognition(s).label.values, rm.recognition(s).response.values) for s in sessions]
print(f"variant, label coding as the paper states (0 = old), recomputed: mean AUC {beh.auc_stated_coding.mean():.3f}, "
f"range {beh.auc_stated_coding.min():.3f} to {beh.auc_stated_coding.max():.3f}")
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist(beh.auc, bins=20); ax.set_xlabel("AUC"); ax.set_ylabel("sessions"); plt.tight_layout(); plt.show()
AUC over 87 sessions: mean 0.753, SD 0.094 (paper 0.74 ± 0.10)
range 0.445 to 0.919; sessions below 0.6: 6; at or below 0.5: 1 chance AUC is 0.5; one-sample t-test of AUC against 0.5: t = 25.1, p = 2.2e-41 hit rate 0.693, false alarm rate 0.269 (means over sessions) variant, sessions 44 and 55 dropped: mean 0.754, SD 0.095
variant, label coding as the paper states (0 = old), recomputed: mean AUC 0.247, range 0.081 to 0.555
Result: reproduced. The AUC over 87 sessions is 0.753 with a standard deviation of 0.094, both within tolerance of 0.74 ± 0.10, and well above the chance level of 0.5 (p = 2.2e-41). Recomputed with the label coding the paper states, the mean AUC is 0.247 (0.081 to 0.555), below chance.
fig5c-auc-sites¶
Performance was similar at the two study sites, with an average AUC of 0.74 ± 0.09 for CS sessions and 0.69 ± 0.13 for TWH sessions.
In [19]:
g = beh.groupby("site").auc.agg(["size", "mean", lambda v: v.std(ddof=1)]).rename(columns={"<lambda_0>": "sd"})
print(g.round(3).to_string())
print("paper: CS 0.74 ± 0.09, TWH 0.69 ± 0.13")
w = stats.ttest_ind(beh[beh.site == "CS"].auc, beh[beh.site == "TWH"].auc, equal_var=False)
print(f"CS against TWH, Welch's t-test: t = {w.statistic:.2f}, p = {w.pvalue:.3f}; difference of means {g.loc['CS', 'mean'] - g.loc['TWH', 'mean']:.3f}")
print("TWH sessions:", beh[beh.site == "TWH"][["session", "auc"]].round(3).to_dict("records"))
size mean sd
site
CS 78 0.760 0.085
TWH 9 0.688 0.138
paper: CS 0.74 ± 0.09, TWH 0.69 ± 0.13
CS against TWH, Welch's t-test: t = 1.53, p = 0.161; difference of means 0.072
TWH sessions: [{'session': 2001, 'auc': 0.798}, {'session': 2002, 'auc': 0.787}, {'session': 2003, 'auc': 0.628}, {'session': 2004, 'auc': 0.544}, {'session': 2005, 'auc': 0.445}, {'session': 2007, 'auc': 0.774}, {'session': 2008, 'auc': 0.738}, {'session': 2009, 'auc': 0.616}, {'session': 2011, 'auc': 0.865}]
Result: reproduced. CS sessions have an AUC of 0.760 ± 0.085 (78 sessions) and TWH sessions 0.688 ± 0.138 (9 sessions), all four values within tolerance. The sites do not differ significantly (Welch's t-test p = 0.161).
fig5b-roc-asymmetric¶
The ROC curve is asymmetric (p < 0.05), as expected for declarative memory.
In [20]:
# Sessions whose five interior z(FA) or z(H) points are all equal have no real z-ROC slope: with
# equal z(FA) the regression has no unique fit (it can return a negative slope), with equal z(H)
# it returns 0. They are left out of the main result and kept in a variant.
syn_lab = np.r_[np.ones(50), np.zeros(50)].astype(int)
syn_resp = np.ones(100, int); syn_resp[:36] = 6; syn_resp[50:54] = 6
print(f"synthetic session answering only 1 and 6 (36 of 50 old, 4 of 50 new at 6): slope {rm.zroc_slope(syn_lab, syn_resp):.4f}, "
f"degenerate: {rm.zroc_degenerate(syn_lab, syn_resp)}")
print("sessions with a degenerate z-ROC (left out):")
print(beh[beh.zroc_degenerate][["session", "site", "auc", "zroc_slope", "hit_rate", "fa_rate"]].round(3).to_string(index=False))
sl = beh[~beh.zroc_degenerate].zroc_slope
t1 = stats.ttest_1samp(sl, 1.0)
wx = stats.wilcoxon(sl - 1.0)
print(f"\nmain, z-ROC slope over {len(sl)} sessions: mean {sl.mean():.3f}, median {sl.median():.3f}, SD {sl.std(ddof=1):.3f}, range {sl.min():.2f} to {sl.max():.2f}")
print(f"one-sample t-test against 1: t = {t1.statistic:.2f}, p = {t1.pvalue:.4g}; sessions with slope below 1: {int((sl < 1).sum())} of {len(sl)}")
print(f"variant, Wilcoxon signed-rank against 1: W = {wx.statistic:.0f}, p = {wx.pvalue:.4g}")
sl_all = beh.zroc_slope
print(f"variant, all {len(sl_all)} sessions with the degenerate ones: mean {sl_all.mean():.3f}, below 1 in {int((sl_all < 1).sum())}, "
f"t-test p = {stats.ttest_1samp(sl_all, 1.0).pvalue:.4g}, Wilcoxon p = {stats.wilcoxon(sl_all - 1.0).pvalue:.4g}")
pooled_lab = np.concatenate([rm.recognition(s).label.values for s in sessions])
pooled_resp = np.concatenate([rm.recognition(s).response.values for s in sessions])
print(f"variant, slope of the pooled ROC (all trials of all sessions): {rm.zroc_slope(pooled_lab, pooled_resp):.3f}")
once = beh[~beh.session.isin(DUPLICATE_SECOND) & ~beh.zroc_degenerate].zroc_slope
print(f"variant, sessions 44 and 55 also dropped: {len(once)} sessions, mean {once.mean():.3f}, p = {stats.ttest_1samp(once, 1.0).pvalue:.2g}")
print("\nper session: AUC, z-ROC slope, hit and false alarm rates")
print(beh[["session", "site", "auc", "zroc_slope", "zroc_degenerate", "hit_rate", "fa_rate"]].round(3).to_string(index=False))
synthetic session answering only 1 and 6 (36 of 50 old, 4 of 50 new at 6): slope -0.2074, degenerate: True
sessions with a degenerate z-ROC (left out):
session site auc zroc_slope hit_rate fa_rate
59 CS 0.845 -0.000 0.88 0.20
60 CS 0.764 -0.247 0.70 0.16
93 CS 0.510 0.415 0.40 0.38
99 CS 0.817 -0.207 0.72 0.08
119 CS 0.845 -0.535 0.86 0.16
2008 TWH 0.738 -0.413 0.72 0.24
main, z-ROC slope over 81 sessions: mean 0.792, median 0.795, SD 0.289, range 0.04 to 1.67
one-sample t-test against 1: t = -6.49, p = 6.587e-09; sessions with slope below 1: 64 of 81
variant, Wilcoxon signed-rank against 1: W = 487, p = 3.292e-08
variant, all 87 sessions with the degenerate ones: mean 0.726, below 1 in 70, t-test p = 1.631e-09, Wilcoxon p = 1.586e-09
variant, slope of the pooled ROC (all trials of all sessions): 0.816
variant, sessions 44 and 55 also dropped: 79 sessions, mean 0.795, p = 2.1e-08
per session: AUC, z-ROC slope, hit and false alarm rates
session site auc zroc_slope zroc_degenerate hit_rate fa_rate
5 CS 0.767 0.775 False 0.70 0.26
6 CS 0.720 0.629 False 0.72 0.32
7 CS 0.658 0.952 False 0.66 0.46
9 CS 0.800 0.882 False 0.70 0.24
17 CS 0.729 1.065 False 0.66 0.30
18 CS 0.812 0.915 False 0.44 0.02
20 CS 0.604 0.505 False 0.88 0.76
21 CS 0.571 1.144 False 0.80 0.70
23 CS 0.768 0.102 False 0.58 0.06
24 CS 0.676 1.669 False 0.56 0.20
26 CS 0.681 1.080 False 0.44 0.14
27 CS 0.830 0.895 False 0.82 0.30
28 CS 0.749 0.642 False 0.80 0.48
32 CS 0.794 0.795 False 0.70 0.24
38 CS 0.736 0.736 False 0.92 0.60
39 CS 0.798 0.816 False 0.84 0.46
41 CS 0.721 0.922 False 0.80 0.36
43 CS 0.789 0.536 False 0.60 0.10
44 CS 0.683 0.585 False 0.42 0.06
47 CS 0.734 0.822 False 0.68 0.34
48 CS 0.779 0.959 False 0.62 0.18
49 CS 0.858 0.689 False 0.66 0.16
50 CS 0.703 1.297 False 0.70 0.38
52 CS 0.890 0.959 False 0.86 0.26
54 CS 0.852 0.630 False 0.70 0.14
55 CS 0.741 0.712 False 0.70 0.32
56 CS 0.697 0.831 False 0.58 0.10
58 CS 0.886 1.074 False 0.80 0.16
59 CS 0.845 -0.000 True 0.88 0.20
60 CS 0.764 -0.247 True 0.70 0.16
61 CS 0.805 1.014 False 0.82 0.28
63 CS 0.681 0.787 False 0.54 0.22
64 CS 0.669 0.533 False 0.70 0.36
66 CS 0.654 1.360 False 0.80 0.50
67 CS 0.574 0.621 False 0.50 0.34
68 CS 0.631 0.953 False 0.78 0.58
69 CS 0.794 1.110 False 0.80 0.32
70 CS 0.785 0.626 False 0.70 0.24
73 CS 0.760 0.208 False 0.58 0.02
74 CS 0.785 0.558 False 0.68 0.28
76 CS 0.792 0.545 False 0.58 0.06
77 CS 0.754 0.796 False 0.52 0.16
78 CS 0.864 0.653 False 0.80 0.22
85 CS 0.624 1.083 False 0.52 0.36
92 CS 0.885 0.569 False 0.86 0.20
93 CS 0.510 0.415 True 0.40 0.38
96 CS 0.809 0.920 False 0.68 0.20
97 CS 0.856 0.719 False 0.72 0.22
99 CS 0.817 -0.207 True 0.72 0.08
100 CS 0.712 0.042 False 0.46 0.02
101 CS 0.823 0.573 False 0.68 0.18
102 CS 0.836 1.385 False 0.98 0.40
104 CS 0.849 0.640 False 0.92 0.40
105 CS 0.753 0.833 False 0.84 0.48
111 CS 0.838 0.560 False 0.64 0.08
112 CS 0.783 0.635 False 0.74 0.24
113 CS 0.768 0.560 False 0.48 0.08
114 CS 0.769 0.983 False 0.64 0.22
115 CS 0.919 0.806 False 0.94 0.18
116 CS 0.876 0.778 False 0.84 0.26
117 CS 0.797 0.894 False 0.64 0.16
118 CS 0.872 0.430 False 0.84 0.30
119 CS 0.845 -0.535 True 0.86 0.16
120 CS 0.772 0.471 False 0.68 0.16
122 CS 0.828 1.048 False 0.78 0.26
123 CS 0.797 0.732 False 0.60 0.14
124 CS 0.853 0.609 False 0.64 0.04
125 CS 0.757 0.660 False 0.56 0.10
126 CS 0.740 0.821 False 0.66 0.18
127 CS 0.637 1.055 False 0.82 0.48
128 CS 0.597 1.411 False 0.78 0.58
130 CS 0.658 1.197 False 0.90 0.64
131 CS 0.633 0.567 False 0.78 0.54
132 CS 0.825 0.755 False 0.62 0.08
133 CS 0.744 0.650 False 0.50 0.14
134 CS 0.730 0.289 False 0.68 0.18
135 CS 0.741 0.294 False 0.56 0.02
136 CS 0.852 0.880 False 0.74 0.24
2001 TWH 0.798 1.257 False 0.80 0.32
2002 TWH 0.787 0.822 False 0.64 0.24
2003 TWH 0.628 0.804 False 0.26 0.08
2004 TWH 0.544 0.485 False 0.18 0.08
2005 TWH 0.445 0.870 False 0.68 0.72
2007 TWH 0.774 0.741 False 0.90 0.44
2008 TWH 0.738 -0.413 True 0.72 0.24
2009 TWH 0.616 1.110 False 0.86 0.68
2011 TWH 0.865 0.799 False 0.86 0.14
Result: reproduced. Six sessions (59, 60, 93, 99, 119 and 2008) answer with too few confidence levels for a z-ROC: their five interior z(FA) or z(H) points are all equal, so the fitted slope is not an estimate (it is negative in four of them and 0 in one). Without them, the z-ROC slope is below 1 in 64 of 81 sessions, with a mean of 0.792 (t-test against 1, p = 6.6e-09; Wilcoxon p = 3.3e-08). With all 87 sessions the mean is 0.726, below 1 in 70 (t-test p = 1.631e-09, Wilcoxon p = 1.586e-09). The pooled ROC has a slope of 0.816.
fig5d-confidence-accuracy¶
Accuracy is significantly higher in high than in low confidence trials (P < 0.001, paired t test) across all sessions.
In [21]:
for split in ("fixed", "dynamic"):
h, l = beh[f"acc_high_{split}"], beh[f"acc_low_{split}"]
ok = h.notna() & l.notna()
tr = stats.ttest_rel(h[ok], l[ok])
tu = stats.ttest_ind(h[ok], l[ok])
print(f"({'a' if split == 'fixed' else 'b'}) {split} split: {int(ok.sum())} sessions, {int((~ok).sum())} left out (no responses in a group): "
f"{beh[~ok].session.tolist()}")
print(f" high {h[ok].mean():.1f}% against low {l[ok].mean():.1f}% (means); high > low in {int((h[ok] > l[ok]).sum())} sessions")
print(f" paired t-test: t = {tr.statistic:.2f}, p = {tr.pvalue:.2g}; variant, unpaired (the code's ttest2): p = {tu.pvalue:.2g}")
print(f" mean difference {np.mean(h[ok] - l[ok]):.1f} points, 2.5th to 97.5th percentile of the per-session differences {np.percentile(h[ok] - l[ok], [2.5, 97.5]).round(1)}")
(a) fixed split: 84 sessions, 3 left out (no responses in a group): [93, 99, 2008]
high 80.1% against low 60.9% (means); high > low in 73 sessions
paired t-test: t = 10.16, p = 3.2e-16; variant, unpaired (the code's ttest2): p = 3.4e-17
mean difference 19.2 points, 2.5th to 97.5th percentile of the per-session differences [-15.1 47.2]
(b) dynamic split: 84 sessions, 3 left out (no responses in a group): [93, 99, 2008]
high 78.3% against low 60.4% (means); high > low in 69 sessions
paired t-test: t = 8.87, p = 1.2e-13; variant, unpaired (the code's ttest2): p = 1.3e-15
mean difference 17.9 points, 2.5th to 97.5th percentile of the per-session differences [-20.3 47.2]
Result: reproduced. Under the fixed split, accuracy is 80.1% for high and 60.9% for low confidence over 84 sessions (paired t-test p = 3.2e-16); under the dynamic split 78.3% and 60.4% (p = 1.2e-13). Three sessions without responses in one group are left out. The unpaired test the code runs also gives p below 0.001.
Single-unit selectivity¶
Every unit is tested on the recognition trials of its file. The firing rate in the window is the spike count divided by the window length. Results for VS, MS and dual cells come from this one computation; each claim's section prints its own numbers.
In [22]:
rng_vs = np.random.default_rng(20200304)
rng_ms = np.random.default_rng(78)
rng_ms2 = np.random.default_rng(79) # a second seed, to show how much the bootstrap's noise moves the counts
rng_null = np.random.default_rng(415)
rng_null_code = np.random.default_rng(416) # null of the code's choices, kept apart so the other nulls do not move
N_NULL = 10
rows = []
started = time.time()
for s in sessions:
r = rm.recognition(s)
on, cat, lab_, resp = r.start_time.values, r.category.values, r.label.values, r.response.values
correct = ((resp >= 4) & (lab_ == rm.OLD)) | ((resp <= 3) & (lab_ == rm.NEW))
names = r.groupby("category").category_name.first()
for k, sp in enumerate(s["spikes"]):
row = {"session": s["session"], "site": s["site"], "unit": k,
"channel": int(s["electrodes"].origChannel.values[s["unit_electrode"][k]]),
"electrode_row": int(s["unit_electrode"][k]),
"location": s["electrodes"].location.values[s["unit_electrode"][k]]}
for wname, w in (("methods", rm.WINDOW_METHODS), ("code", rm.WINDOW_CODE)):
x = rm.window_rates(sp, on, w)
row[f"p_vs_{wname}"] = rm.anova_p(x, cat)
row[f"p_ms_{wname}_all"] = rm.bootstrap_p(x[lab_ == rm.OLD], x[lab_ == rm.NEW], rng_ms)
row[f"p_ms_{wname}_correct"] = rm.bootstrap_p(x[(lab_ == rm.OLD) & correct], x[(lab_ == rm.NEW) & correct], rng_ms)
row[f"p_ms_{wname}_all_seed2"] = rm.bootstrap_p(x[lab_ == rm.OLD], x[lab_ == rm.NEW], rng_ms2)
row[f"p_ms_{wname}_correct_seed2"] = rm.bootstrap_p(x[(lab_ == rm.OLD) & correct], x[(lab_ == rm.NEW) & correct], rng_ms2)
row[f"best_{wname}"] = names[pd.Series(x).groupby(cat).mean().idxmax()]
x = rm.window_rates(sp, on, rm.WINDOW_METHODS)
row["null_vs"] = np.mean([rm.anova_p(x, rng_null.permutation(cat)) < 0.05 for _ in range(N_NULL)])
shuffled = [rng_null.permutation(lab_) for _ in range(N_NULL)]
row["null_ms"] = np.mean([rm.bootstrap_p(x[L == rm.OLD], x[L == rm.NEW], rng_null) < 0.05 for L in shuffled])
# null of the code's choices: labels shuffled among the correct trials, code's window
xc, lc = rm.window_rates(sp, on, rm.WINDOW_CODE)[correct], lab_[correct]
row["null_ms_code_correct"] = np.mean([rm.bootstrap_p(xc[L == rm.OLD], xc[L == rm.NEW], rng_null_code) < 0.05
for L in (rng_null_code.permutation(lc) for _ in range(N_NULL))])
row["rate_hz"] = len(sp) / (s["trials"].stop_time.max() - s["trials"].start_time.min())
rows.append(row)
cells = pd.DataFrame(rows)
cells["vs"] = cells.p_vs_methods < 0.05
cells["ms"] = cells.p_ms_methods_all < 0.05
cells["once"] = ~cells.session.isin(DUPLICATE_SECOND)
print(f"{len(cells)} units tested in {time.time() - started:.0f} s; none dropped. Recognition trials per unit: 100 in every file.")
print(f"units with no spike in any recognition window (Methods' window): {int((cells.p_vs_methods == 1).sum())}")
print("correct trials per session (old, new):", sorted({(int((((rm.recognition(s).response >= 4) & (rm.recognition(s).label == 1))).sum()),
int((((rm.recognition(s).response <= 3) & (rm.recognition(s).label == 0))).sum())) for s in sessions})[:5], "... (first five)")
1864 units tested in 21 s; none dropped. Recognition trials per unit: 100 in every file. units with no spike in any recognition window (Methods' window): 6
correct trials per session (old, new): [(9, 46), (13, 46), (20, 31), (21, 47), (22, 43)] ... (first five)
vs-cells-proportion¶
343 of the 1863 cells (18.4%) are visually selective, by a one-way ANOVA (P < 0.05) of the firing rate in recognition trials, in a 1 s window starting 200 ms after stimulus onset, across the 5 visual categories.
In [23]:
def share(mask, base=None):
base = cells if base is None else base
return f"{int(mask.sum())} of {len(base)} ({100 * mask.mean():.1f}%)"
print("main (Methods' window, every unit):", share(cells.vs), " paper: 343 of 1863 (18.4%)")
print("variant, code's window 200-1700 ms:", share(cells.p_vs_code < 0.05))
print("variant, duplicated spike trains counted once:", share(cells[cells.once].vs, cells[cells.once]))
print(f"null, categories shuffled within session ({N_NULL} shuffles per unit): {100 * cells.null_vs.mean():.1f}% (expected 5%)")
for thr in (0.01, 0.05, 0.10):
print(f"threshold p < {thr}: {100 * (cells.p_vs_methods < thr).mean():.1f}%")
print("\nby session (main), sessions with at least 20 units:")
bs = cells.groupby("session").agg(units=("vs", "size"), vs=("vs", "sum"))
print(bs[bs.units >= 20].assign(pct=lambda d: (100 * d.vs / d.units).round(1)).T.to_string())
print("range of per-session percentages (sessions with at least 20 units):", (100 * bs[bs.units >= 20].vs / bs[bs.units >= 20].units).min().round(1),
"to", (100 * bs[bs.units >= 20].vs / bs[bs.units >= 20].units).max().round(1))
print("by area:", cells.groupby("location").vs.apply(lambda v: f"{int(v.sum())}/{len(v)}").to_dict())
main (Methods' window, every unit): 332 of 1864 (17.8%) paper: 343 of 1863 (18.4%)
variant, code's window 200-1700 ms: 343 of 1864 (18.4%)
variant, duplicated spike trains counted once: 326 of 1809 (18.0%)
null, categories shuffled within session (10 shuffles per unit): 4.7% (expected 5%)
threshold p < 0.01: 8.6%
threshold p < 0.05: 17.8%
threshold p < 0.1: 25.6%
by session (main), sessions with at least 20 units:
session 5 7 9 17 18 20 24 26 32 38 43 44 48 49 54 55 58 59 60 64 69 70 74 76 77 93 99 100 102 104 105 113 115 116 122 124 125 126 2004
units 24.0 38.0 25.0 23.0 33.0 20.0 21.0 33.0 43.0 20.0 23.0 23.0 45.0 27.0 32.0 32.0 49.0 21.0 22.0 33.0 64.0 33.0 35.0 53.0 30.0 35.0 30.0 42.0 44.0 42.0 43.0 21.0 35.0 36.0 27.0 43.0 46.0 34.0 31.0
vs 2.0 8.0 3.0 2.0 3.0 0.0 3.0 0.0 13.0 0.0 7.0 5.0 8.0 8.0 6.0 1.0 17.0 5.0 5.0 1.0 9.0 5.0 10.0 11.0 8.0 9.0 8.0 12.0 8.0 15.0 10.0 2.0 10.0 3.0 2.0 7.0 3.0 6.0 6.0
pct 8.3 21.1 12.0 8.7 9.1 0.0 14.3 0.0 30.2 0.0 30.4 21.7 17.8 29.6 18.8 3.1 34.7 23.8 22.7 3.0 14.1 15.2 28.6 20.8 26.7 25.7 26.7 28.6 18.2 35.7 23.3 9.5 28.6 8.3 7.4 16.3 6.5 17.6 19.4
range of per-session percentages (sessions with at least 20 units): 0.0 to 35.7
by area: {'Left Amygdala': '81/498', 'Left Hippocampus': '60/365', 'Right Amygdala': '116/594', 'Right Hippocampus': '75/407'}
Result: reproduced. With the Methods' window, 332 of 1864 units (17.8%) are visually selective, within tolerance of 343 (18.4%), against 4.7% with shuffled categories. The code's window gives 343 of 1864 (18.4%), the paper's numbers. Counting duplicated spike trains once gives 18.0%.
ms-cells-proportion¶
146 of the 1863 cells (7.8%) are memory selective, by a two-tailed bootstrap comparison of mean firing rates (1,000 runs, P < 0.05) between the 50 novel and 50 familiar recognition trials in a 1 s window starting 200 ms after stimulus onset.
In [24]:
print("main (Methods' window, all 50 and 50 trials):", share(cells.ms), " paper: 146 of 1863 (7.8%)")
print("variant, Methods' window, correct trials only:", share(cells.p_ms_methods_correct < 0.05))
print("variant, code's window, all trials:", share(cells.p_ms_code_all < 0.05))
print("variant, code's choices (code's window, correct trials only):", share(cells.p_ms_code_correct < 0.05))
print("variant, duplicated spike trains counted once:", share(cells[cells.once].ms, cells[cells.once]))
print(f"null, labels shuffled within session ({N_NULL} shuffles per unit): {100 * cells.null_ms.mean():.1f}% (expected 5%)")
for thr in (0.01, 0.05, 0.10):
print(f"threshold p < {thr}: {100 * (cells.p_ms_methods_all < thr).mean():.1f}%")
print("by site:", cells.groupby("site").ms.apply(lambda v: f"{int(v.sum())}/{len(v)}").to_dict())
print("by area:", cells.groupby("location").ms.apply(lambda v: f"{int(v.sum())}/{len(v)}").to_dict())
print("\nsecond bootstrap seed:")
for col, name in (("p_ms_methods_all", "main"), ("p_ms_methods_correct", "Methods' window, correct trials"),
("p_ms_code_all", "code's window, all trials"), ("p_ms_code_correct", "code's choices")):
a, b = cells[col] < 0.05, cells[col + "_seed2"] < 0.05
print(f" {name}: {int(a.sum())} and {int(b.sum())}; units selected under one seed only: {int((a != b).sum())}")
n, k, f0 = len(cells), int(cells.ms.sum()), cells.null_ms.mean()
lo0, hi0 = stats.binomtest(int(round(f0 * n * N_NULL)), n * N_NULL).proportion_ci()
print(f"\nnull fraction {100 * f0:.2f}% from {n * N_NULL} shuffled tests, 95% CI {100 * lo0:.2f}% to {100 * hi0:.2f}%")
print(f"main count {k} against the null: expected {n * f0:.1f}, binomial p (one-sided) = {stats.binomtest(k, n, f0, alternative='greater').pvalue:.3g}; "
f"against 5%: expected {n * 0.05:.1f}, p = {stats.binomtest(k, n, 0.05, alternative='greater').pvalue:.3g}")
k2 = int((cells.p_ms_code_correct < 0.05).sum())
print(f"code's choices count {k2} against 5%: p = {stats.binomtest(k2, n, 0.05, alternative='greater').pvalue:.3g}")
f2 = cells.null_ms_code_correct.mean()
lo2, hi2 = stats.binomtest(int(round(f2 * n * N_NULL)), n * N_NULL).proportion_ci()
print(f"null of the code's choices (labels shuffled among correct trials, code's window, {N_NULL} shuffles per unit): "
f"{100 * f2:.2f}%, 95% CI {100 * lo2:.2f}% to {100 * hi2:.2f}%")
print(f"code's choices count {k2} against this null: expected {n * f2:.1f}, binomial p (one-sided) = "
f"{stats.binomtest(k2, n, f2, alternative='greater').pvalue:.3g}")
ms_cols = ["p_ms_methods_all", "p_ms_methods_all_seed2", "p_ms_methods_correct", "p_ms_methods_correct_seed2",
"p_ms_code_all", "p_ms_code_all_seed2", "p_ms_code_correct", "p_ms_code_correct_seed2"]
any_ms = (cells[ms_cols] < 0.05).any(axis=1)
print(f"\nper-unit p-values of every unit memory selective under some reading or seed ({int(any_ms.sum())} units):")
print(cells[any_ms][["session", "unit", "p_vs_methods", "p_vs_code"] + ms_cols].round(3).to_string(index=False))
main (Methods' window, all 50 and 50 trials): 110 of 1864 (5.9%) paper: 146 of 1863 (7.8%)
variant, Methods' window, correct trials only: 143 of 1864 (7.7%)
variant, code's window, all trials: 108 of 1864 (5.8%)
variant, code's choices (code's window, correct trials only): 147 of 1864 (7.9%)
variant, duplicated spike trains counted once: 105 of 1809 (5.8%)
null, labels shuffled within session (10 shuffles per unit): 5.4% (expected 5%)
threshold p < 0.01: 1.1%
threshold p < 0.05: 5.9%
threshold p < 0.1: 11.6%
by site: {'CS': '103/1787', 'TWH': '7/77'}
by area: {'Left Amygdala': '25/498', 'Left Hippocampus': '26/365', 'Right Amygdala': '36/594', 'Right Hippocampus': '23/407'}
second bootstrap seed:
main: 110 and 118; units selected under one seed only: 22
Methods' window, correct trials: 143 and 140; units selected under one seed only: 17
code's window, all trials: 108 and 109; units selected under one seed only: 15
code's choices: 147 and 149; units selected under one seed only: 10
null fraction 5.39% from 18640 shuffled tests, 95% CI 5.07% to 5.73%
main count 110 against the null: expected 100.5, binomial p (one-sided) = 0.177; against 5%: expected 93.2, p = 0.0443
code's choices count 147 against 5%: p = 7.12e-08
null of the code's choices (labels shuffled among correct trials, code's window, 10 shuffles per unit): 5.58%, 95% CI 5.26% to 5.92%
code's choices count 147 against this null: expected 104.1, binomial p (one-sided) = 2.49e-05
per-unit p-values of every unit memory selective under some reading or seed (288 units):
session unit p_vs_methods p_vs_code p_ms_methods_all p_ms_methods_all_seed2 p_ms_methods_correct p_ms_methods_correct_seed2 p_ms_code_all p_ms_code_all_seed2 p_ms_code_correct p_ms_code_correct_seed2
5 6 0.817 0.860 0.047 0.045 0.287 0.272 0.569 0.557 0.857 0.841
5 12 0.475 0.046 0.031 0.031 0.033 0.043 0.036 0.042 0.022 0.027
5 14 0.392 0.333 0.043 0.050 0.374 0.356 0.112 0.129 0.488 0.480
5 17 0.515 0.658 0.027 0.021 0.198 0.202 0.038 0.028 0.271 0.292
6 5 0.206 0.552 0.142 0.153 0.030 0.027 0.023 0.029 0.003 0.005
6 7 0.444 0.725 0.017 0.015 0.061 0.046 0.101 0.088 0.211 0.217
6 10 0.178 0.716 0.047 0.054 0.001 0.001 0.126 0.112 0.020 0.023
6 14 0.314 0.370 0.123 0.112 0.043 0.045 0.541 0.531 0.309 0.313
7 1 0.955 0.960 0.020 0.012 0.111 0.110 0.026 0.032 0.127 0.141
7 13 0.672 0.842 0.010 0.004 0.030 0.028 0.004 0.005 0.048 0.047
7 16 0.985 0.886 0.003 0.001 0.031 0.018 0.069 0.067 0.164 0.180
7 35 0.818 0.925 0.335 0.332 0.029 0.027 0.651 0.651 0.214 0.197
7 37 0.000 0.000 0.018 0.013 0.060 0.068 0.165 0.182 0.533 0.532
9 0 0.462 0.452 0.156 0.152 0.041 0.036 0.444 0.436 0.402 0.420
9 3 0.641 0.273 0.472 0.451 0.048 0.041 0.568 0.571 0.056 0.062
9 21 0.444 0.160 0.027 0.031 0.090 0.067 0.007 0.010 0.039 0.030
17 15 0.341 0.446 0.030 0.020 0.044 0.054 0.160 0.159 0.110 0.127
17 16 0.721 0.561 0.038 0.033 0.227 0.239 0.109 0.092 0.217 0.229
17 21 0.491 0.124 0.539 0.553 0.050 0.043 0.598 0.608 0.109 0.088
18 3 0.673 0.297 0.072 0.067 0.084 0.091 0.070 0.070 0.039 0.032
18 7 0.564 0.917 0.000 0.001 0.000 0.000 0.119 0.129 0.145 0.147
18 8 0.451 0.190 0.058 0.071 0.005 0.007 0.121 0.143 0.004 0.011
18 9 0.256 0.038 0.593 0.614 0.029 0.024 0.192 0.210 0.028 0.028
18 11 0.092 0.190 0.365 0.386 0.042 0.047 0.363 0.335 0.922 0.924
18 12 0.000 0.010 0.140 0.159 0.001 0.001 0.346 0.337 0.034 0.026
18 15 0.306 0.106 0.145 0.133 0.044 0.048 0.022 0.025 0.020 0.016
18 20 0.083 0.039 0.728 0.727 0.022 0.021 0.533 0.526 0.000 0.002
18 21 0.669 0.041 0.650 0.578 0.181 0.166 0.123 0.122 0.005 0.005
18 23 0.510 0.842 0.531 0.544 0.226 0.257 0.210 0.234 0.034 0.049
18 25 0.163 0.073 0.067 0.060 0.039 0.058 0.009 0.012 0.000 0.002
18 27 0.102 0.307 0.316 0.323 0.025 0.031 0.115 0.109 0.007 0.007
18 29 0.810 0.555 0.002 0.000 0.019 0.023 0.000 0.002 0.007 0.013
18 30 0.294 0.050 0.452 0.434 0.158 0.182 0.083 0.070 0.045 0.047
18 31 0.553 0.313 0.014 0.018 0.007 0.005 0.006 0.002 0.006 0.008
20 2 0.478 0.521 0.073 0.088 0.075 0.083 0.052 0.054 0.015 0.015
20 11 0.584 0.721 0.451 0.444 0.531 0.540 0.677 0.669 0.035 0.046
20 14 0.560 0.460 0.015 0.018 0.513 0.492 0.031 0.024 0.924 0.920
21 0 0.535 0.845 0.152 0.155 0.194 0.190 0.062 0.035 0.051 0.055
21 6 0.654 0.721 0.009 0.011 0.353 0.345 0.007 0.010 0.132 0.136
23 1 0.436 0.231 0.197 0.187 0.125 0.100 0.111 0.101 0.027 0.021
23 4 0.019 0.072 0.117 0.126 0.081 0.092 0.010 0.013 0.017 0.020
23 11 0.423 0.215 0.178 0.151 0.118 0.125 0.029 0.032 0.063 0.062
24 5 0.149 0.618 0.089 0.112 0.019 0.027 0.272 0.278 0.064 0.071
24 9 0.825 0.850 0.212 0.220 0.312 0.335 0.024 0.021 0.078 0.061
26 9 0.379 0.805 0.010 0.007 0.062 0.069 0.019 0.011 0.069 0.057
26 18 0.716 0.974 0.053 0.044 0.119 0.115 0.270 0.226 0.168 0.159
26 25 0.672 0.544 0.065 0.065 0.036 0.035 0.273 0.267 0.163 0.164
26 28 0.157 0.486 0.002 0.001 0.000 0.000 0.332 0.325 0.035 0.030
26 32 0.220 0.023 0.520 0.536 0.112 0.098 0.463 0.455 0.000 0.000
27 8 0.821 0.639 0.024 0.025 0.047 0.053 0.144 0.173 0.313 0.316
27 9 0.876 0.897 0.081 0.063 0.355 0.385 0.044 0.047 0.259 0.281
32 1 0.608 0.897 0.176 0.189 0.022 0.025 0.098 0.091 0.015 0.007
32 11 0.471 0.515 0.019 0.017 0.024 0.024 0.114 0.109 0.172 0.172
32 18 0.968 0.779 0.592 0.630 0.171 0.202 0.152 0.128 0.023 0.025
32 33 0.689 0.947 0.166 0.186 0.115 0.114 0.111 0.106 0.052 0.045
32 37 0.099 0.115 0.222 0.254 0.146 0.124 0.049 0.055 0.028 0.028
38 1 0.475 0.673 0.144 0.142 0.342 0.356 0.010 0.019 0.049 0.048
38 5 0.560 0.560 0.165 0.177 0.044 0.046 0.238 0.235 0.186 0.189
38 8 0.736 0.909 0.070 0.076 0.030 0.030 0.284 0.309 0.835 0.836
38 10 0.689 0.753 0.018 0.023 0.038 0.032 0.037 0.041 0.039 0.044
38 17 0.123 0.312 0.313 0.293 0.372 0.369 0.008 0.009 0.023 0.030
39 8 0.328 0.121 0.014 0.015 0.013 0.016 0.022 0.023 0.019 0.013
39 9 0.628 0.649 0.086 0.068 0.117 0.087 0.048 0.036 0.109 0.103
41 1 0.000 0.000 0.103 0.090 0.025 0.032 0.160 0.170 0.059 0.067
41 2 0.369 0.507 0.078 0.060 0.192 0.220 0.030 0.034 0.159 0.181
41 10 0.202 0.833 0.412 0.426 0.040 0.031 0.223 0.244 0.019 0.021
41 11 0.562 0.626 0.765 0.723 0.738 0.751 0.020 0.023 0.004 0.003
41 13 0.919 0.868 0.039 0.024 0.006 0.006 0.117 0.112 0.022 0.018
43 1 0.288 0.471 0.019 0.014 0.109 0.115 0.046 0.040 0.100 0.098
43 10 0.455 0.293 0.055 0.062 0.077 0.066 0.042 0.035 0.036 0.035
43 12 0.668 0.641 0.362 0.366 0.105 0.085 0.048 0.057 0.017 0.046
44 0 0.009 0.218 0.335 0.325 0.334 0.323 0.092 0.109 0.040 0.046
44 3 0.051 0.128 0.082 0.097 0.305 0.305 0.023 0.020 0.041 0.031
44 10 0.072 0.286 0.414 0.419 0.109 0.094 0.060 0.035 0.035 0.048
44 13 0.254 0.404 0.146 0.131 0.054 0.037 0.055 0.065 0.061 0.055
44 14 0.295 0.061 0.158 0.149 0.170 0.189 0.123 0.132 0.043 0.037
47 0 0.672 0.333 0.042 0.039 0.098 0.137 0.368 0.371 0.596 0.599
47 2 0.209 0.254 0.269 0.261 0.006 0.013 0.214 0.215 0.069 0.065
47 5 0.645 0.422 0.030 0.028 0.217 0.259 0.269 0.244 0.974 0.967
47 6 0.221 0.096 0.016 0.013 0.000 0.000 0.175 0.170 0.053 0.043
47 12 0.526 0.185 0.087 0.100 0.067 0.063 0.087 0.090 0.052 0.039
48 2 0.304 0.351 0.767 0.776 0.228 0.247 0.479 0.491 0.021 0.032
48 15 0.611 0.374 0.053 0.044 0.028 0.029 0.061 0.059 0.053 0.041
48 17 0.088 0.127 0.461 0.420 0.240 0.258 0.082 0.090 0.022 0.023
48 26 0.470 0.352 0.180 0.167 0.229 0.210 0.047 0.062 0.112 0.086
48 30 0.064 0.022 0.063 0.052 0.047 0.041 0.933 0.935 0.796 0.807
48 33 0.056 0.226 0.167 0.180 0.396 0.423 0.036 0.038 0.118 0.122
48 38 0.125 0.558 0.561 0.529 0.121 0.112 0.543 0.550 0.034 0.037
49 16 0.066 0.197 0.042 0.042 0.077 0.072 0.151 0.167 0.336 0.305
49 18 0.246 0.560 0.049 0.050 0.116 0.111 0.138 0.133 0.175 0.189
52 1 0.781 0.776 0.214 0.276 0.120 0.116 0.173 0.141 0.040 0.029
52 4 0.540 0.582 0.369 0.384 0.007 0.011 0.138 0.168 0.000 0.001
52 5 0.676 0.598 0.071 0.083 0.027 0.025 0.021 0.025 0.005 0.009
52 13 0.483 0.346 0.072 0.073 0.038 0.032 0.170 0.211 0.169 0.164
52 17 0.741 0.244 0.053 0.043 0.046 0.060 0.011 0.008 0.043 0.037
54 0 0.012 0.072 0.084 0.075 0.082 0.085 0.058 0.057 0.040 0.043
54 7 0.002 0.003 0.032 0.032 0.000 0.001 0.145 0.132 0.000 0.003
54 8 0.071 0.020 0.046 0.055 0.039 0.047 0.309 0.298 0.271 0.266
54 12 0.002 0.006 0.413 0.422 0.028 0.027 0.680 0.666 0.089 0.090
54 15 0.468 0.485 0.039 0.035 0.067 0.071 0.476 0.491 0.529 0.536
54 16 0.779 0.560 0.022 0.019 0.059 0.058 0.016 0.018 0.022 0.021
55 3 0.729 0.696 0.202 0.208 0.039 0.053 0.055 0.076 0.013 0.008
55 6 0.027 0.303 0.026 0.028 0.016 0.008 0.039 0.039 0.054 0.067
55 19 0.165 0.432 0.039 0.032 0.064 0.067 0.073 0.069 0.088 0.086
55 20 0.821 0.966 0.036 0.032 0.024 0.023 0.039 0.026 0.044 0.053
55 21 0.468 0.185 0.031 0.037 0.092 0.079 1.000 1.000 0.513 0.558
55 24 0.846 0.897 0.039 0.049 0.040 0.043 0.027 0.038 0.057 0.053
55 29 0.900 0.935 0.057 0.054 0.075 0.083 0.035 0.030 0.074 0.063
56 4 0.002 0.021 0.826 0.828 0.031 0.032 0.846 0.834 0.058 0.070
58 9 0.758 0.827 0.232 0.249 0.007 0.005 0.492 0.491 0.234 0.236
58 21 0.394 0.659 0.006 0.003 0.002 0.006 0.002 0.000 0.001 0.001
58 24 0.444 0.741 0.047 0.045 0.121 0.123 0.129 0.138 0.445 0.459
58 28 0.000 0.000 0.215 0.218 0.047 0.049 0.214 0.213 0.146 0.131
58 30 0.527 0.424 0.022 0.029 0.020 0.014 0.127 0.126 0.164 0.138
59 17 0.110 0.031 0.114 0.089 0.029 0.035 0.092 0.095 0.023 0.034
60 0 0.263 0.272 0.040 0.038 0.062 0.069 0.062 0.054 0.078 0.085
60 3 0.512 0.829 0.051 0.049 0.147 0.136 0.025 0.019 0.080 0.065
60 11 0.163 0.029 0.084 0.123 0.047 0.060 0.116 0.122 0.043 0.033
60 15 0.169 0.091 0.528 0.515 0.916 0.921 0.037 0.046 0.166 0.155
61 2 0.101 0.215 0.035 0.040 0.081 0.086 0.043 0.042 0.053 0.062
61 3 0.327 0.851 0.083 0.106 0.075 0.054 0.020 0.016 0.022 0.015
63 2 0.764 0.883 0.263 0.244 0.071 0.065 0.059 0.055 0.010 0.008
63 9 0.149 0.175 0.046 0.041 0.187 0.179 0.068 0.069 0.155 0.151
63 11 0.999 0.761 0.020 0.017 0.011 0.010 0.062 0.052 0.066 0.078
64 10 0.686 0.501 0.868 0.838 0.134 0.132 0.473 0.502 0.028 0.026
64 19 0.036 0.218 0.065 0.048 0.249 0.260 0.186 0.155 0.634 0.619
64 21 0.877 0.836 0.007 0.007 0.066 0.050 0.013 0.013 0.094 0.087
64 32 0.126 0.135 0.047 0.045 0.247 0.241 0.060 0.072 0.279 0.298
66 2 0.798 0.462 0.016 0.012 0.044 0.040 0.004 0.004 0.114 0.111
67 13 0.663 0.422 0.035 0.060 0.004 0.006 0.179 0.170 0.007 0.010
67 14 0.645 0.169 0.061 0.049 0.446 0.457 0.312 0.294 0.891 0.879
68 6 0.445 0.335 0.756 0.776 0.739 0.743 0.045 0.049 0.091 0.094
68 13 0.284 0.539 0.104 0.091 0.003 0.007 0.255 0.258 0.043 0.036
69 4 0.007 0.009 0.156 0.178 0.079 0.101 0.064 0.088 0.040 0.047
69 14 0.524 0.940 0.841 0.850 0.043 0.049 0.714 0.705 0.109 0.095
69 15 0.991 0.887 0.088 0.081 0.043 0.037 0.100 0.100 0.016 0.010
69 22 0.516 0.360 0.051 0.042 0.201 0.192 0.110 0.121 0.388 0.366
69 30 0.462 0.396 0.030 0.035 0.058 0.055 0.147 0.145 0.356 0.351
69 51 0.567 0.277 0.020 0.014 0.013 0.005 0.003 0.003 0.004 0.005
69 56 0.197 0.212 0.317 0.320 0.125 0.108 0.032 0.025 0.028 0.030
69 61 0.833 0.207 0.569 0.583 0.029 0.033 0.825 0.829 0.536 0.515
70 12 0.041 0.129 0.195 0.204 0.056 0.052 0.369 0.350 0.048 0.046
70 16 0.320 0.182 0.025 0.035 0.055 0.069 0.019 0.024 0.031 0.034
70 26 0.150 0.310 0.469 0.449 0.033 0.033 0.219 0.226 0.008 0.014
74 4 0.131 0.168 0.092 0.084 0.108 0.146 0.042 0.057 0.311 0.336
74 6 0.414 0.567 0.138 0.123 0.106 0.120 0.029 0.033 0.022 0.019
74 22 0.020 0.075 0.078 0.066 0.009 0.008 0.092 0.086 0.023 0.018
74 28 0.226 0.037 0.017 0.022 0.103 0.079 0.201 0.201 0.111 0.100
74 34 0.433 0.644 0.104 0.103 0.055 0.063 0.055 0.072 0.045 0.033
76 3 0.030 0.022 0.230 0.210 0.036 0.025 0.186 0.173 0.009 0.016
76 5 0.120 0.144 0.338 0.371 0.063 0.052 0.057 0.057 0.007 0.003
76 12 0.035 0.093 0.075 0.077 0.145 0.130 0.020 0.011 0.021 0.024
76 15 0.074 0.046 0.141 0.163 0.063 0.065 0.090 0.104 0.022 0.025
76 24 0.025 0.032 0.247 0.248 0.500 0.488 0.050 0.045 0.110 0.085
76 26 0.524 0.133 0.227 0.238 0.010 0.023 0.321 0.297 0.085 0.089
76 30 0.093 0.292 0.061 0.064 0.000 0.000 0.099 0.093 0.007 0.005
76 31 0.309 0.247 0.044 0.039 0.154 0.146 0.198 0.194 0.376 0.375
76 36 0.001 0.004 0.478 0.471 0.027 0.024 0.468 0.478 0.039 0.027
76 37 0.056 0.099 0.664 0.666 0.059 0.075 0.205 0.172 0.000 0.000
76 39 0.081 0.190 1.000 1.000 0.045 0.048 0.863 0.879 0.076 0.062
76 44 0.496 0.600 0.302 0.330 0.044 0.049 0.074 0.061 0.027 0.020
77 2 0.006 0.000 0.800 0.807 0.004 0.002 0.655 0.635 0.002 0.005
77 3 0.666 0.808 0.438 0.437 0.274 0.288 0.088 0.090 0.024 0.020
77 6 0.353 0.145 0.024 0.024 0.071 0.080 0.016 0.034 0.167 0.153
77 7 0.468 0.039 0.037 0.030 0.057 0.061 0.073 0.082 0.155 0.153
77 25 0.185 0.295 0.459 0.469 0.045 0.038 0.195 0.209 0.019 0.030
78 2 0.002 0.000 0.377 0.402 0.145 0.118 0.104 0.091 0.010 0.007
78 10 0.160 0.058 0.367 0.393 0.192 0.192 0.334 0.343 0.047 0.043
78 11 0.222 0.122 0.024 0.024 0.029 0.027 0.084 0.069 0.117 0.133
85 3 0.124 0.119 0.236 0.234 0.110 0.116 0.192 0.194 0.040 0.054
92 6 0.360 0.495 0.036 0.033 0.045 0.046 0.036 0.024 0.037 0.029
93 3 0.429 0.684 0.165 0.190 0.065 0.053 0.394 0.405 0.035 0.037
93 5 0.892 0.871 0.203 0.206 0.038 0.039 0.096 0.073 0.005 0.006
93 8 0.622 0.776 0.006 0.003 0.028 0.036 0.024 0.021 0.186 0.174
93 10 0.038 0.018 0.261 0.239 0.485 0.480 0.029 0.020 0.481 0.488
93 21 0.030 0.098 0.011 0.006 0.055 0.049 0.118 0.136 0.068 0.059
96 7 0.000 0.000 0.006 0.011 0.005 0.002 0.025 0.015 0.014 0.020
99 5 0.353 0.317 0.014 0.019 0.080 0.085 0.012 0.006 0.043 0.051
99 6 0.614 0.442 0.014 0.017 0.037 0.034 0.010 0.007 0.032 0.027
99 14 0.004 0.001 0.012 0.011 0.028 0.039 0.030 0.029 0.132 0.129
99 20 0.073 0.069 0.011 0.014 0.001 0.003 0.032 0.019 0.023 0.018
99 21 0.368 0.305 0.094 0.103 0.021 0.021 0.173 0.194 0.041 0.038
99 22 0.002 0.009 0.002 0.002 0.005 0.007 0.017 0.019 0.020 0.037
100 1 0.823 0.408 0.021 0.028 0.001 0.004 0.043 0.039 0.080 0.087
100 2 0.058 0.074 0.101 0.105 0.013 0.010 0.032 0.026 0.000 0.000
100 10 0.013 0.011 0.794 0.782 0.047 0.061 0.971 0.963 0.083 0.077
100 18 0.137 0.069 0.354 0.321 0.390 0.404 0.018 0.023 0.022 0.025
100 40 0.064 0.033 0.021 0.028 0.012 0.021 0.184 0.190 0.126 0.141
100 41 0.137 0.339 0.382 0.414 0.242 0.218 0.029 0.019 0.088 0.083
101 6 0.725 0.927 0.104 0.102 0.046 0.034 0.009 0.014 0.009 0.009
101 12 0.344 0.180 0.001 0.005 0.001 0.001 0.000 0.003 0.001 0.002
102 3 0.623 0.540 0.076 0.071 0.035 0.045 0.033 0.028 0.030 0.044
102 4 0.126 0.287 0.280 0.296 0.135 0.154 0.058 0.069 0.017 0.012
102 22 0.642 0.315 0.351 0.363 0.664 0.652 0.369 0.340 0.023 0.030
102 38 0.320 0.167 0.096 0.089 0.015 0.015 0.052 0.038 0.007 0.009
102 40 0.171 0.326 0.013 0.006 0.066 0.061 0.025 0.021 0.160 0.144
102 41 0.354 0.665 0.061 0.055 0.186 0.190 0.047 0.046 0.097 0.069
104 2 0.036 0.122 0.081 0.065 0.005 0.005 0.272 0.254 0.020 0.020
104 3 0.085 0.048 0.374 0.375 0.328 0.341 0.028 0.030 0.081 0.078
104 15 0.079 0.009 0.060 0.038 0.126 0.128 0.057 0.035 0.106 0.099
104 16 0.612 0.415 0.169 0.186 0.163 0.174 0.033 0.026 0.034 0.038
104 20 0.105 0.103 0.633 0.641 0.030 0.037 0.562 0.569 0.033 0.032
104 22 0.006 0.007 0.041 0.045 0.066 0.048 0.006 0.007 0.018 0.021
104 23 0.348 0.271 0.074 0.088 0.177 0.172 0.026 0.019 0.033 0.028
104 24 0.661 0.190 0.032 0.042 0.120 0.132 0.765 0.778 0.771 0.756
104 36 0.001 0.002 0.043 0.043 0.060 0.060 0.086 0.090 0.189 0.166
105 7 0.000 0.001 0.057 0.071 0.003 0.001 0.063 0.087 0.000 0.001
105 8 0.025 0.048 0.219 0.202 0.124 0.131 0.144 0.156 0.048 0.047
105 9 0.000 0.002 0.086 0.093 0.061 0.055 0.080 0.073 0.028 0.023
105 10 0.758 0.471 0.028 0.025 0.602 0.590 0.042 0.046 0.904 0.910
105 12 0.976 0.861 0.007 0.010 0.026 0.033 0.007 0.004 0.009 0.012
105 18 0.966 0.996 0.013 0.011 0.028 0.027 0.133 0.139 0.168 0.162
105 19 0.487 0.883 0.674 0.701 0.032 0.042 0.863 0.864 0.282 0.274
105 25 0.044 0.162 0.000 0.002 0.018 0.007 0.000 0.002 0.003 0.005
105 30 0.677 0.296 0.219 0.187 0.070 0.066 0.199 0.252 0.044 0.041
105 37 0.246 0.450 0.048 0.051 0.082 0.079 0.049 0.050 0.181 0.178
105 41 0.170 0.382 0.017 0.011 0.042 0.048 0.066 0.068 0.225 0.224
111 3 0.518 0.674 0.058 0.049 0.218 0.226 0.076 0.079 0.499 0.504
111 9 0.685 0.130 0.005 0.004 0.000 0.001 0.030 0.029 0.014 0.026
111 12 0.539 0.966 0.078 0.075 0.501 0.506 0.018 0.022 0.364 0.347
112 6 0.779 0.794 0.345 0.354 0.372 0.359 0.014 0.021 0.023 0.036
113 5 0.879 0.974 0.167 0.149 0.268 0.272 0.070 0.061 0.022 0.018
113 10 0.634 0.674 0.119 0.126 0.200 0.178 0.035 0.041 0.155 0.134
114 0 0.000 0.000 0.080 0.079 0.009 0.015 0.154 0.133 0.026 0.019
114 5 0.082 0.009 0.386 0.335 0.016 0.019 0.189 0.192 0.045 0.046
114 6 0.671 0.448 0.627 0.623 0.106 0.093 0.137 0.150 0.005 0.007
115 5 0.374 0.132 0.225 0.220 0.156 0.155 0.053 0.049 0.021 0.022
115 17 0.997 0.660 0.077 0.085 0.056 0.042 0.070 0.086 0.035 0.037
115 20 0.071 0.171 0.027 0.026 0.016 0.021 0.090 0.082 0.070 0.073
115 22 0.097 0.323 0.005 0.006 0.010 0.018 0.009 0.006 0.007 0.004
115 32 0.042 0.079 0.002 0.006 0.005 0.006 0.007 0.004 0.007 0.013
116 3 0.787 0.135 0.156 0.150 0.046 0.041 0.022 0.024 0.003 0.011
116 8 0.778 0.438 0.013 0.020 0.138 0.143 0.014 0.019 0.096 0.093
116 14 0.523 0.229 0.004 0.004 0.027 0.023 0.423 0.431 0.784 0.782
116 16 0.777 0.357 0.505 0.522 0.170 0.155 0.119 0.107 0.021 0.030
116 28 0.544 0.584 0.059 0.043 0.049 0.064 0.230 0.221 0.283 0.317
116 29 0.156 0.407 0.060 0.040 0.177 0.168 0.183 0.208 0.404 0.423
117 2 0.233 0.219 0.015 0.012 0.001 0.000 0.006 0.003 0.000 0.000
117 3 0.237 0.124 0.010 0.019 0.004 0.002 0.003 0.006 0.001 0.000
117 8 0.428 0.131 0.008 0.017 0.001 0.000 0.004 0.003 0.000 0.000
117 10 0.056 0.052 0.169 0.179 0.107 0.106 0.017 0.022 0.036 0.043
118 5 0.389 0.541 0.063 0.049 0.003 0.001 0.675 0.686 0.353 0.352
118 10 0.741 0.530 0.351 0.359 0.095 0.118 0.055 0.062 0.039 0.039
119 8 0.507 0.236 0.026 0.024 0.024 0.014 0.033 0.024 0.035 0.041
119 9 0.087 0.105 0.032 0.049 0.078 0.099 0.147 0.143 0.230 0.236
122 26 0.716 0.955 0.186 0.207 0.388 0.421 0.057 0.048 0.186 0.188
123 3 0.540 0.764 0.045 0.053 0.382 0.375 0.029 0.034 0.462 0.472
124 1 0.291 0.181 0.140 0.147 0.042 0.045 0.221 0.236 0.336 0.327
124 12 0.159 0.027 0.010 0.016 0.037 0.028 0.041 0.037 0.040 0.043
124 21 0.215 0.379 0.035 0.033 0.007 0.008 0.390 0.389 0.295 0.241
124 23 0.518 0.533 0.515 0.491 0.893 0.908 0.047 0.054 0.137 0.111
124 32 0.131 0.013 0.285 0.328 0.039 0.049 0.203 0.201 0.035 0.030
125 2 0.472 0.444 0.072 0.044 0.163 0.182 0.049 0.059 0.066 0.049
125 9 0.657 0.571 0.010 0.015 0.020 0.024 0.015 0.014 0.029 0.038
125 20 0.668 0.775 0.069 0.069 0.001 0.005 0.224 0.242 0.066 0.048
125 41 0.301 0.184 0.297 0.334 0.010 0.008 0.313 0.306 0.047 0.052
125 43 0.403 0.508 0.912 0.899 0.021 0.038 0.702 0.679 0.252 0.223
126 6 0.139 0.036 0.016 0.015 0.176 0.170 0.027 0.025 0.221 0.260
126 8 0.038 0.019 0.016 0.009 0.020 0.023 0.126 0.137 0.182 0.170
126 10 0.346 0.559 0.247 0.253 0.060 0.063 0.131 0.115 0.037 0.037
126 11 0.269 0.719 0.039 0.035 0.108 0.105 0.023 0.028 0.085 0.072
126 29 0.546 0.406 0.027 0.045 0.018 0.014 0.107 0.085 0.089 0.090
126 31 0.359 0.665 0.018 0.019 0.035 0.046 0.258 0.276 0.483 0.516
126 32 0.969 0.649 0.257 0.248 0.011 0.023 0.185 0.192 0.067 0.065
127 5 0.955 0.979 0.000 0.002 0.041 0.029 0.071 0.076 0.358 0.371
128 0 0.018 0.129 0.222 0.206 0.005 0.007 0.135 0.127 0.015 0.029
130 5 0.733 0.326 0.099 0.081 0.056 0.057 0.059 0.064 0.028 0.040
130 10 0.163 0.044 0.471 0.474 0.021 0.018 0.753 0.750 0.322 0.314
131 5 0.032 0.401 0.164 0.142 0.027 0.032 0.154 0.143 0.023 0.026
133 1 0.487 0.589 0.021 0.021 0.036 0.027 0.045 0.042 0.020 0.030
134 5 0.018 0.048 0.087 0.085 0.218 0.257 0.042 0.045 0.108 0.085
134 11 0.133 0.605 0.439 0.449 0.445 0.436 0.055 0.039 0.085 0.073
136 1 0.524 0.809 0.591 0.599 0.041 0.040 0.642 0.615 0.107 0.093
2001 2 0.668 0.935 0.040 0.047 0.038 0.032 0.112 0.114 0.084 0.066
2001 6 0.457 0.333 0.050 0.044 0.199 0.187 0.311 0.281 0.525 0.488
2001 7 0.622 0.655 0.099 0.082 0.015 0.014 0.856 0.860 0.248 0.288
2003 4 0.933 0.719 1.000 1.000 0.046 0.059 0.395 0.377 0.021 0.012
2004 2 0.307 0.329 0.041 0.032 0.586 0.584 0.117 0.128 0.496 0.482
2004 4 0.006 0.008 0.070 0.036 0.040 0.055 0.218 0.202 0.145 0.134
2004 11 0.318 0.154 0.160 0.138 0.019 0.015 0.008 0.009 0.227 0.192
2004 12 0.427 0.653 0.009 0.008 0.000 0.002 0.002 0.002 0.000 0.000
2004 13 0.063 0.090 0.030 0.033 0.000 0.000 0.015 0.018 0.000 0.000
2004 14 0.794 0.787 0.011 0.004 0.000 0.000 0.001 0.000 0.000 0.000
2004 17 0.403 0.043 0.791 0.801 0.051 0.036 0.712 0.696 0.017 0.021
2004 19 0.613 0.774 0.031 0.033 0.000 0.000 0.000 0.000 0.000 0.000
2004 21 0.016 0.016 0.117 0.134 0.397 0.401 0.044 0.045 0.287 0.326
2004 30 0.056 0.205 0.200 0.180 0.004 0.006 0.136 0.132 0.005 0.010
2005 2 0.099 0.088 0.011 0.016 0.036 0.035 0.000 0.001 0.008 0.009
Result: not reproduced. The test the Methods describe, all 50 new and 50 old trials in a 1 s window, finds 110 of 1864 units (5.9%), outside the tolerance of 146 (7.8%). This is not distinguishable from the 5.4% found with shuffled labels (one-sided binomial p = 0.18) and only just above 5% (p = 0.044). A second bootstrap seed gives 118 (6.3%), also outside the tolerance. The code's choices, correct trials only in a 1.5 s window, give 147 and 149 (7.9% and 8.0%) under the two seeds, and correct trials in the Methods' window 143 and 140, all matching the paper. The code's choices also exceed their own null, labels shuffled among the correct trials in the code's window (5.58%, 147 against 104.1 expected, one-sided binomial p = 2.5e-05). The outcome comes from the Methods disagreeing with the paper's code on the trials used, not from the data: the paper's number follows its code.
dual-cells-proportion¶
26 cells (1.4%) are both visually and memory selective.
In [25]:
dual = cells.vs & cells.ms
print("main:", share(dual), " paper: 26 of 1863 (1.4%)")
print("expected if VS and MS were independent:", f"{len(cells) * cells.vs.mean() * cells.ms.mean():.1f} units")
print("variant, code's choices:", share((cells.p_vs_code < 0.05) & (cells.p_ms_code_correct < 0.05)))
print("variant, duplicated spike trains counted once:", share(dual[cells.once], cells[cells.once]))
print("seed 2, main:", share(cells.vs & (cells.p_ms_methods_all_seed2 < 0.05)),
"; code's choices:", share((cells.p_vs_code < 0.05) & (cells.p_ms_code_correct_seed2 < 0.05)))
main: 12 of 1864 (0.6%) paper: 26 of 1863 (1.4%) expected if VS and MS were independent: 19.6 units variant, code's choices: 26 of 1864 (1.4%) variant, duplicated spike trains counted once: 11 of 1809 (0.6%) seed 2, main: 14 of 1864 (0.8%) ; code's choices: 26 of 1864 (1.4%)
Result: not reproduced. Under the Methods' tests 12 units (0.6%) are both VS and MS (14 with the second bootstrap seed), outside the tolerance of 26 (1.4%) and below the 19.6 expected if the two were independent. The code's choices give 26 (1.4%) under both seeds, the paper's numbers. As for ms-cells-proportion, the outcome comes from the Methods disagreeing with the code, not from the data.
fig5f-site-proportions¶
The percentages of VS, MS and dual cells are comparable between CS sessions (18.4%, 7.7% and 1.4%) and the new TWH sessions (19.4%, 10.4% and 1.3%).
In [26]:
cells["dual"] = cells.vs & cells.ms
PAPER = {("vs", "CS"): 18.4, ("ms", "CS"): 7.7, ("dual", "CS"): 1.4, ("vs", "TWH"): 19.4, ("ms", "TWH"): 10.4, ("dual", "TWH"): 1.3}
within = 0
for kind in ("vs", "ms", "dual"):
a = cells[cells.site == "CS"][kind]; b = cells[cells.site == "TWH"][kind]
tab2 = [[int(a.sum()), int((~a).sum())], [int(b.sum()), int((~b).sum())]]
p = stats.fisher_exact(tab2).pvalue
for site, v in (("CS", a), ("TWH", b)):
ok = abs(100 * v.mean() - PAPER[(kind, site)]) <= 0.1 * PAPER[(kind, site)]
within += ok
print(f"{kind:4s} {site:3s}: {int(v.sum())} of {len(v)} = {100 * v.mean():.1f}% (paper {PAPER[(kind, site)]}%, within 10%: {ok})")
print(f" Fisher's exact test CS against TWH: p = {p:.3f}")
print("values within tolerance:", within, "of 6")
code_sel = {"vs": cells.p_vs_code < 0.05, "ms": cells.p_ms_code_correct < 0.05}
code_sel["dual"] = code_sel["vs"] & code_sel["ms"]
for kind, sel in code_sel.items():
print(f"variant, code's choices, {kind}:", cells.assign(x=sel).groupby("site").x.apply(lambda v: f"{int(v.sum())} of {len(v)} = {100 * v.mean():.1f}%").to_dict())
print("\n95% Clopper-Pearson intervals of the main percentages:")
for kind in ("vs", "ms", "dual"):
for site in ("CS", "TWH"):
v = cells[cells.site == site][kind]
ci = stats.binomtest(int(v.sum()), len(v)).proportion_ci()
print(f" {kind:4s} {site:3s}: {100 * v.mean():.1f}% ({100 * ci.low:.1f} to {100 * ci.high:.1f})")
vs CS : 315 of 1787 = 17.6% (paper 18.4%, within 10%: True)
vs TWH: 17 of 77 = 22.1% (paper 19.4%, within 10%: False)
Fisher's exact test CS against TWH: p = 0.360
ms CS : 103 of 1787 = 5.8% (paper 7.7%, within 10%: False)
ms TWH: 7 of 77 = 9.1% (paper 10.4%, within 10%: False)
Fisher's exact test CS against TWH: p = 0.215
dual CS : 12 of 1787 = 0.7% (paper 1.4%, within 10%: False)
dual TWH: 0 of 77 = 0.0% (paper 1.3%, within 10%: False)
Fisher's exact test CS against TWH: p = 1.000
values within tolerance: 1 of 6
variant, code's choices, vs: {'CS': '328 of 1787 = 18.4%', 'TWH': '15 of 77 = 19.5%'}
variant, code's choices, ms: {'CS': '139 of 1787 = 7.8%', 'TWH': '8 of 77 = 10.4%'}
variant, code's choices, dual: {'CS': '25 of 1787 = 1.4%', 'TWH': '1 of 77 = 1.3%'}
95% Clopper-Pearson intervals of the main percentages:
vs CS : 17.6% (15.9 to 19.5)
vs TWH: 22.1% (13.4 to 33.0)
ms CS : 5.8% (4.7 to 6.9)
ms TWH: 9.1% (3.7 to 17.8)
dual CS : 0.7% (0.3 to 1.2)
dual TWH: 0.0% (0.0 to 4.7)
Result: partially reproduced. The sites do not differ by Fisher's exact test (VS p = 0.360, MS p = 0.215, dual p = 1.000), but only 1 of the 6 site values is within tolerance under the Methods' tests (CS VS 17.6%). The code's choices give CS 18.4%, 7.8% and 1.4% and TWH 19.5%, 10.4% and 1.3%, the values of Figure 5f to within 0.1 point, so the partial outcome comes from the Methods disagreeing with the code. TWH has 77 units, so its percentages rest on 17, 7 and 0 cells under the main tests, with 95% intervals of 13.4% to 33.0%, 3.7% to 17.8% and 0.0% to 4.7%; the tests have little power, and "no difference" carries little weight.
twh-sessions¶
The dataset includes 9 sessions from a new institution, Toronto Western Hospital.
In [27]:
twh = info[info.site == "TWH"]
print("TWH sessions:", len(twh), twh.session.tolist(), "; subjects:", twh.subject.nunique(), "; units:", int(twh.units.sum()))
print("institution of these files:", twh.institution.unique().tolist())
TWH sessions: 9 [2001, 2002, 2003, 2004, 2005, 2007, 2008, 2009, 2011] ; subjects: 8 ; units: 77 institution of these files: ['Toronto Western Hospital']
Result: reproduced. Nine sessions (2001 to 2011) of 8 subjects come from Toronto Western Hospital, with 77 units.
fig4-example-vs-neuron¶
Neuron 132-177-4-RA (session 132, channel 177, cell 4, right amygdala) is visually selective, responding most to images of animals.
In [28]:
ex = cells[(cells.session == 132) & (cells.channel == 177)]
print(ex[["session", "unit", "channel", "electrode_row", "location", "p_vs_methods", "p_vs_code", "best_methods", "best_code"]].to_string(index=False))
s = by_id[132]; r = rm.recognition(s); k = int(ex.unit.iloc[0])
for wname, w in (("Methods' window", rm.WINDOW_METHODS), ("code's window", rm.WINDOW_CODE)):
x = rm.window_rates(s["spikes"][k], r.start_time.values, w)
print(wname, "mean rate (Hz) by category:", pd.Series(x).groupby(r.category_name.values).mean().round(2).to_dict())
print("units on channel 177 in session 132:", len(ex))
print("units of session 132 (0-based unit index, electrode row, channel):",
cells[cells.session == 132][["unit", "electrode_row", "channel"]].values.tolist())
session unit channel electrode_row location p_vs_methods p_vs_code best_methods best_code
132 4 177 4 Right Amygdala 1.197068e-22 7.885374e-19 5animals 5animals
Methods' window mean rate (Hz) by category: {'1cars': 0.45, '2food': 0.5, '3people': 0.85, '4spatial': 0.25, '5animals': 6.75}
code's window mean rate (Hz) by category: {'1cars': 0.4, '2food': 0.37, '3people': 0.7, '4spatial': 0.23, '5animals': 5.6}
units on channel 177 in session 132: 1
units of session 132 (0-based unit index, electrode row, channel): [[0, 0, 146], [1, 1, 148], [2, 2, 149], [3, 3, 150], [4, 4, 177], [5, 5, 181], [6, 6, 182], [7, 7, 183]]
Result: reproduced. The single unit on channel 177 of session 132 sits in the right amygdala at electrode row 4, and it is the fifth of the file's eight units (0-based index 4), so the label's cell number equals both its electrode row and its 0-based unit index. It is category selective (ANOVA p = 1.2e-22 in the Methods' window, 7.9e-19 in the code's) and fires most to animals (6.75 Hz against 0.25 to 0.85 Hz for the other categories).
previous-release-link¶
The sessions of the earlier release (Faraut et al. 2018) are found in this dandiset under the same subject, session id and task variant.
The earlier release's sessions come from Table 1 of Faraut et al. 2018 (doi:10.1038/sdata.2018.10), copied into the helper module. The task variant of a file is read from its stimulus set, through the variant Online-only Table 1 gives each stimulus set (shown under counts-sessions).
In [29]:
norm = lambda v: v.str.replace("CS", "C").str.replace(r"^([A-Z])0*", r"\1", regex=True)
early = rm.faraut_2018_table_1()
early["subject_norm"] = norm(early.table_subject)
first_cat = pd.Series({s["session"]: sorted(s["trials"].category_name.unique())[0] for s in sessions})
cat_to_variant = pd.crosstab(first_cat.rename("cat"), pd.Series(first_cat.index.map(var), index=first_cat.index, name="variant")).idxmax(axis=1)
files = info[["session", "label", "age", "sex"]].assign(subject_norm=lambda d: norm(d.label),
variant_file=lambda d: d.session.map(first_cat).map(cat_to_variant).values)
lk = early.merge(files, on="session", how="left", suffixes=("_2018", "_file"))
lk["found"] = lk.label.notna()
lk["agrees"] = lk.found & (lk.subject_norm_2018 == lk.subject_norm_file) & (lk.variant == lk.variant_file)
print("sessions in Faraut et al. 2018 Table 1:", len(early), "; subjects:", early.table_subject.nunique())
print("found in the dandiset by session id:", int(lk.found.sum()), "of", len(lk), f"({100 * lk.found.mean():.1f}%)")
print("found and agreeing on subject and task variant:", int(lk.agrees.sum()), "of", len(lk), f"({100 * lk.agrees.mean():.1f}%)")
print("\nnot found or not agreeing:")
print(lk[~lk.agrees][["table_subject", "session", "variant", "label", "variant_file"]].to_string(index=False))
subj18 = lk[lk.found].groupby("table_subject").agg(age_2018=("age_2018", "first"), sex_2018=("sex_2018", "first"),
age_file=("age_file", "first"), sex_file=("sex_file", "first"))
subj18["age_file"] = subj18.age_file.astype(int)
same = (subj18.age_2018 == subj18.age_file) & (subj18.sex_2018 == subj18.sex_file)
print("\nmatched subjects with the same age and sex in the files as in the 2018 table:", int(same.sum()), "of", len(subj18))
print(subj18[~same].to_string())
print("\nsession 98 in the code's session list:")
print(code_list[code_list.session == 98].to_string(index=False))
sessions in Faraut et al. 2018 Table 1: 65 ; subjects: 42
found in the dandiset by session id: 64 of 65 (98.5%)
found and agreeing on subject and task variant: 64 of 65 (98.5%)
not found or not agreeing:
table_subject session variant label variant_file
H48 98 2 NaN NaN
matched subjects with the same age and sex in the files as in the 2018 table: 39 of 42
age_2018 sex_2018 age_file sex_file
table_subject
C24 47 F 57 F
H33 29 M 30 M
H48 54 M 32 F
session 98 in the code's session list:
session active name ID age sex filename
98 False P48HMH_NOID98 H48 32 F P48HMH_NOID98.nwb
Result: partially reproduced. 64 of the 65 sessions of Faraut et al. 2018 (98.5%) are in the dandiset under the same subject and task variant. Session 98 of subject H48 is not; its section in the code's session list is commented out. Age and sex agree for 39 of the 42 matched subjects; C24, H33 and H48 differ, the same three subjects that differ from Online-only Table 1, and in each the two papers' tables agree with each other.
previous-release-additions¶
Relative to the earlier release of this task's data (Faraut et al. 2018), this release adds 23 sessions, 17 subjects and 288 neurons.
In [30]:
new_ses = info[~info.session.isin(early.session)]
old_subj = set(norm(info[info.session.isin(early.session)].label))
new_subj = sorted(set(norm(new_ses.label)) - old_subj)
print("new sessions:", len(new_ses), new_ses.session.tolist(), " paper: 23")
print("new subjects:", len(new_subj), new_subj, " paper: 17")
print("units in the new sessions:", int(new_ses.units.sum()), " paper: 288")
print(" of which CS:", int(new_ses[new_ses.site == "CS"].units.sum()), "in", int((new_ses.site == "CS").sum()), "sessions;",
"TWH:", int(new_ses[new_ses.site == "TWH"].units.sum()), "in", int((new_ses.site == "TWH").sum()), "sessions")
print("subjects of the new sessions that also have earlier sessions:", sorted(set(norm(new_ses.label)) & old_subj))
print("\nvariant, differences of the totals the two papers state (2020 minus 2018): sessions 87 - 65 =", 87 - 65,
"; subjects 59 - 42 =", 59 - 42, "; neurons 1863 - 1576 =", 1863 - 1576)
print("units in the dandiset's sessions that are in the 2018 table:", int(info[info.session.isin(early.session)].units.sum()),
" (Faraut et al. 2018: 1576 neurons in 65 sessions)")
print("Faraut et al. 2018 neurons plus the stated additions: 1576 + 288 =", 1576 + 288, "; unit rows in the dandiset:", int(info.units.sum()))
new sessions: 23 [122, 123, 124, 125, 126, 127, 128, 130, 131, 132, 133, 134, 135, 136, 2001, 2002, 2003, 2004, 2005, 2007, 2008, 2009, 2011] paper: 23 new subjects: 17 ['C53', 'C54', 'C55', 'C56', 'C57', 'C58', 'C60', 'C61', 'C62', 'T100', 'T101', 'T103', 'T107', 'T88', 'T89', 'T90', 'T98'] paper: 17 units in the new sessions: 325 paper: 288 of which CS: 248 in 14 sessions; TWH: 77 in 9 sessions subjects of the new sessions that also have earlier sessions: [] variant, differences of the totals the two papers state (2020 minus 2018): sessions 87 - 65 = 22 ; subjects 59 - 42 = 17 ; neurons 1863 - 1576 = 287 units in the dandiset's sessions that are in the 2018 table: 1539 (Faraut et al. 2018: 1576 neurons in 65 sessions) Faraut et al. 2018 neurons plus the stated additions: 1576 + 288 = 1864 ; unit rows in the dandiset: 1864
Result: partially reproduced. Against the session list of Faraut et al. 2018, the dandiset has 23 new sessions (14 CS, 9 TWH) from 17 new subjects, the paper's numbers. The new sessions hold 325 units, outside the tolerance of 288; the sessions also in the earlier release hold 1539, where Faraut et al. 2018 gives 1576 neurons for its 65 sessions. The paper's 288 equals 1864 minus 1576, which suggests it is a difference of totals, not a count of the new sessions. The missing session 98 (previous-release-link) does not change the new sessions or subjects; it may account for part of the 37 units fewer in the earlier sessions.
fig3-example-vs-neuron (blocked)¶
Neuron 144-140-3-RH (session 144, channel 140, cell 3, right hippocampus) is visually selective for phones.
Session 144 is not in the dandiset (data check) or in Online-only Table 1, so the claim is blocked.
As a weak check of whether the neuron is in the dandiset under another session id, every unit on
original channel 140 in the right hippocampus is listed with its category ANOVA and preferred
category. This does not decide the claim. The paper's code is read for the same purpose, not run: the
lines of its demo script python/export/demo/export.py that name session 144, and the code's
session list for subject T103, are printed. The cell number of the label is tested against the
Figure 4 neuron, whose label 132-177-4-RA has cell number 4.
In [31]:
cand = cells[(cells.channel == 140) & (cells.location == "Right Hippocampus")]
print("units on channel 140 in the right hippocampus:", len(cand))
cand = cand.assign(cluster=[int(by_id[a]["cluster_id"][k]) for a, k in zip(cand.session, cand.unit)])
print(cand[["session", "unit", "electrode_row", "cluster", "p_vs_methods", "best_methods", "p_vs_code", "best_code"]].to_string(index=False))
print("of these, VS (Methods' window) and preferring a phone category:",
int((cand.vs & cand.best_methods.str.contains("hone", case=False)).sum()))
print("category names containing 'phone':", sorted({n for s in sessions for n in s["trials"].category_name.unique() if "hone" in n.lower()}))
print("subject of session 2009:", by_id[2009]["label"], "; sessions of this subject:", info[info.label == by_id[2009]["label"]].session.tolist())
print("\nexport.py at the pinned commit, lines naming the session:",
[l for l in rm.code_file("python/export/demo/export.py").splitlines() if "144" in l or "TWH103" in l])
print("the code's session list, sections of subject T103:")
print(code_list[code_list.ID == "T103"].to_string(index=False))
print("\ncell number of the label against the Figure 4 neuron (label cell number 4):")
f4 = cells[(cells.session == 132) & (cells.channel == 177)].iloc[0]
print(f" Figure 4 unit: 0-based unit index {f4.unit}, 1-based {f4.unit + 1}, electrode row {f4.electrode_row}")
for _, c in cand.iterrows():
print(f" session {c.session} candidate: 0-based unit index {c.unit}, 1-based {c.unit + 1}, electrode row {c.electrode_row}")
print("units of session 2009 (0-based unit index, electrode row, channel):",
cells[cells.session == 2009][["unit", "electrode_row", "channel"]].values.tolist())
units on channel 140 in the right hippocampus: 2
session unit electrode_row cluster p_vs_methods best_methods p_vs_code best_code
2009 1 1 352 2.668321e-02 smallAnimal 4.351638e-03 smallAnimal
2009 2 1 414 1.344250e-09 phones 5.022522e-08 phones
of these, VS (Methods' window) and preferring a phone category: 1
category names containing 'phone': ['phones']
subject of session 2009: T103 ; sessions of this subject: [2009]
export.py at the pinned commit, lines naming the session: ["#session: 'TWH103_111918'", 'NOID = 144', "io = NWBHDF5IO('TWH103_111918.nwb', mode='w')"]
the code's session list, sections of subject T103:
session active name ID age sex filename
2009 True TWH103_NOID2009 T103 49 M TWH103_NOID2009.nwb
cell number of the label against the Figure 4 neuron (label cell number 4):
Figure 4 unit: 0-based unit index 4, 1-based 5, electrode row 4
session 2009 candidate: 0-based unit index 1, 1-based 2, electrode row 1
session 2009 candidate: 0-based unit index 2, 1-based 3, electrode row 1
units of session 2009 (0-based unit index, electrode row, channel): [[0, 0, 137], [1, 1, 140], [2, 1, 140]]
Result: blocked (metadata insufficient). No file has session 144. The paper's demo script sets NOID 144 for the session it names TWH103_111918, and the code's session list gives subject T103 one section, session 2009. In session 2009 two units are on channel 140 in the right hippocampus, both at electrode row 1; unit 2 is visually selective (p = 1.3e-09) and fires most to phones. The Figure 4 neuron sets the convention for the cell number: its cell number 4 is its electrode row and its 0-based unit index, and its 1-based index is 5. Under that convention the phone-selective unit is cell 2 (unit index) or 1 (electrode row), not 3; only its 1-based index is 3. The candidate matches the label's channel, area and category, but neither its cell number nor a recorded mapping of 144 to 2009 identifies it as the neuron of Figure 3.
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome |
|---|---|
| counts-neurons-subjects | reproduced |
| counts-sessions | reproduced |
| table1-demographics | partially reproduced |
| task-structure | partially reproduced |
| new-old-label-coding | not reproduced |
| ttl-codes | partially reproduced |
| experiment-ids | partially reproduced |
| session-dates | partially reproduced |
| units-table | partially reproduced |
| electrodes-table | not reproduced |
| isi-criterion | partially reproduced |
| fig5a-mtl-locations | reproduced |
| fig5c-auc-all | reproduced |
| fig5c-auc-sites | reproduced |
| fig5b-roc-asymmetric | reproduced |
| fig5d-confidence-accuracy | reproduced |
| vs-cells-proportion | reproduced |
| ms-cells-proportion | not reproduced |
| dual-cells-proportion | not reproduced |
| fig5f-site-proportions | partially reproduced |
| twh-sessions | reproduced |
| fig4-example-vs-neuron | reproduced |
| previous-release-link | partially reproduced (with outside input) |
| previous-release-additions | partially reproduced (with outside input) |
| fig3-example-vs-neuron | blocked: session 144 is not in the dandiset, and nothing maps it to a session there |
How the analysis differs from the paper's.
- Selectivity uses the Methods' window (200 to 1200 ms) and, for memory selectivity, all 50 new and 50 old trials. The paper's code uses 200 to 1700 ms and correct trials only; with these choices the counts are 343 VS, 147 MS and 26 dual cells, against the paper's 343, 146 and 26. The MS, dual and Figure 5f outcomes rest on the Methods' description, and come from its disagreement with the code.
- The z-ROC slope is averaged over the 81 sessions whose interior points allow a fit; the six sessions left out are shown as a variant.
- The bootstrap is the code's algorithm with a fixed seed; a second seed moves the main MS count from 110 to 118 and the code's choices from 147 to 149.
- The paper names no test for the ROC asymmetry or for "comparable" proportions; a t-test of the z-ROC slope against 1 and Fisher's exact test were chosen.
- Accuracy by confidence uses a paired test, as the Results say; the code runs an unpaired one, shown as a variant.
- Units are counted per file, as in the paper; counts with the duplicated spike trains of sessions 44 and 55 dropped are shown.
- The earlier release's sessions are taken from Table 1 of Faraut et al. 2018, not from its Dryad archive, which a script could not fetch.
What was found about the dandiset.
- The new/old label is coded 1 = old, the reverse of the paper's text and the files' own column description. A user following the description would compute an AUC of 0.25.
- Sessions 43 and 44, and 54 and 55, hold identical spike trains: the same neurons are deposited under two sessions, so 55 of the 1864 units are repeats. Session 55 also repeats session 54's learning block.
- Seven sessions have 50 encoding trials and one has 97, where the paper describes 100; part of the images later shown as old are missing from their encoding blocks.
- Electrode impedance is empty for every electrode, coordinates are missing for 99 electrodes, the isolation distance is NaN for 493 units, and SNR is not finite for 86.
- Spike times are unsorted in 13 units, and two units of session 58 hold each spike time twice.
- Age or sex differ from Online-only Table 1 for three subjects. One session date keeps its day, and another is set to 2000-01-01.
- Event codes outside Table 1 occur in nine files and experiment ids outside Table 2 in four.
- The Methods give 89 sessions where the paper's table and the files have 87, and the unit count is 1864 where the paper gives 1863.
- The paper's Figure 3 example neuron is from session 144, which is not in the dandiset. The paper's code names NOID 144 as a recording of subject T103, whose one session here is 2009; one unit of session 2009 on channel 140 in the right hippocampus prefers phones, but no record maps 144 to 2009.
- Session 98 of the earlier release (subject H48) is not in the dandiset; its section in the code's session list is commented out. The new sessions hold 325 units, not the 288 the paper gives.