Reproducing Buccino et al. (2020) from dandiset 000034

Paper. Buccino AP, Hurwitz CL, Garcia S, Magland J, Siegle JH, Hurwitz R, Hennig MH (2020). SpikeInterface, a unified framework for spike sorting. eLife 9:e61834. doi:10.7554/eLife.61834. The claims file (claims/000034/10.7554-elife.61834.yaml) was written from the version of record (full text from Europe PMC, PMC7704107, with the figure-supplement legends, the Methods and the author response) and the paper's notebooks in SpikeInterface/spikeinterface.github.io (last changed at ee08fb6).

Dandiset. DANDI:000034, version 0.211030.0713 (doi:10.48324/dandi.000034/0.211030.0713). The paper cites the draft (dandiset/000034/draft), as it was published before any version existed. This is the only published version; the data check below shows that it holds the same six assets (same asset ids and sizes) as today's draft, so the choice does not change the data.

What the dandiset holds. The four recordings the paper analyses (mouse Neuropixels, simulated MEArec Neuropixels with ground truth, rat cortex Neuropixels, mouse retina Biocam) and the two experts' manual curations of the paper's Kilosort2 output on the mouse recording. It holds no output of any of the six spike sorters the paper compares.

Why most claims are blocked. Every finding except the presence of the recordings and the curator counts of fig4a rests on running spike sorters on these recordings: most need all six, and three (fig2e-ks2-unit-types, fig2e-ks2-curated-unit-types, fig3s2) need Kilosort2 alone. Kilosort2, IronClust and HDSort are MATLAB programs (Kilosort2 also needs a CUDA GPU); they cannot be put in a locked Python environment, and this runner has neither MATLAB nor a GPU. The section "Why the sorter claims are blocked" also measures how fast the recordings can be read: chunk by chunk through h5py it would take hours, but one bulk read of a file takes minutes, so reading time alone would not block a sorter run.

Claim	Tested here	Blocked by
recordings-as-described (link)	yes
fig1b-unit-counts	no	not_runnable
fig1c-low-agreement	no	not_runnable
fig1d-most-units-unmatched	no	not_runnable
fig1s2-agreement-scores-high	no	not_runnable
fig1s3-s4-other-recordings-disagree	no	not_runnable
fig2c-sim-unit-counts	no	not_runnable
fig2s1a-precision-recall-bias	no	not_runnable
fig2s1b-accuracy-snr	no	not_runnable
fig2d-ks2-highest-scores	no	not_runnable
fig2e-ks2-unit-types	no	not_runnable
fig2e-ks2-curated-unit-types	no	not_runnable
fig2e-tdc-few-fp-hds-many-fp	no	not_runnable
fig3a-sim-low-agreement	no	not_runnable
fig3b-sorters-prone-to-low-agreement	no	not_runnable
fig3c-fp-found-by-one-sorter	no	not_runnable
fig3s1-ks2-ic-best-pair	no	not_runnable
fig3s1-pairs-remove-fp-four-keep-tp	no	not_runnable
fig3d-snr-rises-with-agreement	no	not_runnable
fig3s2-snr-does-not-separate-ks2-fp	no	not_runnable
fig4a-curators-agree	in part (the curator counts; the 174 needs the Kilosort2 run)
fig4b-single-sorter-match-curated	no	not_runnable
fig4c-consensus-match-curated	no	not_runnable
fig4d-nonconsensus-rarely-curated	no	not_runnable

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.

import time
import warnings

import numpy as np
import pandas as pd

import si000034 as si

warnings.filterwarnings("ignore")
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
started = time.time()
print(f"dandiset {si.DANDISET_ID} version {si.DANDISET_VERSION}")

dandiset 000034 version 0.211030.0713

Data check

Each check of the metadata against the values is stated under its output.

pinned = pd.DataFrame(si.list_assets(si.DANDISET_VERSION))
draft = pd.DataFrame(si.list_assets("draft"))
print("Assets of version", si.DANDISET_VERSION)
print(pinned.to_string(index=False))
same = pinned.merge(draft, on=["path", "size", "asset_id"], how="outer", indicator=True)
print("\nAssets in both the pinned version and the draft, same id and size:", int((same._merge == "both").sum()),
      "; only in one:", int((same._merge != "both").sum()))
info = si.version_info()
print("Version name:", info["metadata"]["name"])
print("Description:", info["metadata"]["description"][-250:])
print("Species in the archive's asset summary:", [s["name"] for s in info["metadata"].get("assetsSummary", {}).get("species", [])])

Assets of version 0.211030.0713
                                                                         path        size                             asset_id
                      sub-mouse412804/sub-mouse412804_ses-20200824T155542.nwb     9693752 af1745f6-e14a-473a-987d-cef7e29a3956
                      sub-mouse412804/sub-mouse412804_ses-20200824T155543.nwb    11732540 c696fc2b-d2e6-4e27-8775-01657193c4a2
    sub-P29-16-05-14-retina02-left/sub-P29-16-05-14-retina02-left_ecephys.nwb  6470908640 9822a813-0dec-4d07-810b-1c13341c168d
                                  sub-mouse412804/sub-mouse412804_ecephys.nwb 13411878072 3e0d4539-902e-46b4-9cfb-119e0f69fbd0
                                                    sub-c1/sub-c1_ecephys.nwb 24916614016 76c211c2-d588-48a7-a4c8-ea103ff0570a
sub-MEAREC-250neuron-Neuropixels/sub-MEAREC-250neuron-Neuropixels_ecephys.nwb 29530187056 6d94dcf4-0b38-4323-8250-04fdc7039a66

Assets in both the pinned version and the draft, same id and size: 6 ; only in one: 0

Version name: SpikeInterface, a unified framework for spike sorting
Description: nalysed in the paper: "SpikeInterface, a unified framework for spike sorting." Alessio P. Buccino, Cole L. Hurwitz, Samuel Garcia, Jeremy Magland, Joshua H. Siegle, Roger Hurwitz, Matthias H. Hennig, eLife - doi: https://doi.org/10.7554/eLife.61834; 
Species in the archive's asset summary: ['House mouse']

The pinned version and the draft hold the same six files. The description says the dandiset holds all recordings analysed in the paper. The archive's species summary lists only house mouse, though the rat recording's file gives its species as a rat (table below).

files = {"curated 1": si.CURATED_1, "curated 2": si.CURATED_2, **si.RECORDINGS}
handles, rows = {}, []
for name, path in files.items():
    f = si.open_nwb(path)
    handles[name] = f
    d = si.describe(f)
    rows.append({"name": name, "path": path.split("/")[-1], **d})
desc = pd.DataFrame(rows).set_index("name")
cols = ["subject_id", "species", "age", "session_start_time", "frames", "channels", "rate_hz", "duration_s",
        "dtype", "chunks", "conversion", "offset", "unit", "starting_time", "units", "unit_columns"]
print(desc[cols].to_string())

                               subject_id        species  age                session_start_time      frames  channels       rate_hz  duration_s    dtype      chunks  conversion  offset   unit  starting_time  units                  unit_columns
name                                                                                                                                                                                                                                               
curated 1                     mouse412804   mus musculus       2020-08-24T15:55:42.004284+02:00         NaN       NaN           NaN         NaN      NaN         NaN         NaN     NaN    NaN            NaN  258.0                 [spike_times]
curated 2                     mouse412804   mus musculus       2020-08-24T15:55:43.012201+02:00         NaN       NaN           NaN         NaN      NaN         NaN         NaN     NaN    NaN            NaN  319.0                 [spike_times]
mouse412804                   mouse412804   mus musculus       2020-08-03T11:57:32.640004+01:00  27000000.0     248.0  30000.000000  900.000000    int16  (26368, 1)    0.000001     NaN  volts            0.0    NaN                           NaN
MEArec       MEAREC_250neuron_Neuropixels                      2020-07-27T09:46:20.619184+01:00  19200000.0     384.0  32000.000000  600.000000  float32  (18750, 1)    0.000001     NaN  volts            0.0  250.0  [spike_times, soma_location]
rat c1                                 c1  rattus rattus       2020-08-07T18:16:07.237608+01:00   8100337.0     384.0  30000.000000  270.011233  float64   (7911, 1)    0.000001     NaN  volts            0.0    NaN                           NaN
retina         P29_16_05_14_retina02_left   mus musculus  P29  2020-08-07T09:42:57.434275+01:00   3148820.0    1024.0  23199.090358  135.730322   uint16  (12301, 1)    0.000001     NaN  volts            0.0    NaN                           NaN

Every recording starts at time 0 and has a regular rate. Each stores its raw data in chunks of one channel by several thousand frames, so a read of a short stretch of all channels fetches one chunk per channel. The curated files hold a Units table with spike times only (no rate, no labels, no link to the Kilosort2 units). The simulated file holds its 250 ground-truth units with their soma locations. No file says which cells are excitatory or inhibitory, or how the recordings were cut or filtered (filtering is "none" and location "unknown" for every channel).

cur = {k: si.load_units(handles[k]) for k in ["curated 1", "curated 2"]}
gt = si.load_units(handles["MEArec"])
for name, units in [*cur.items(), ("MEArec ground truth", gt)]:
    t = np.concatenate(list(units.values()))
    n = np.array([len(v) for v in units.values()])
    print(f"{name}: {len(units)} units, ids {min(units)} to {max(units)}, {len(t)} spikes, times {t.min():.4f} to {t.max():.4f} s, "
          f"spikes per unit {n.min()} to {n.max()} (median {np.median(n):.0f}), sorted within every unit: "
          f"{all(np.all(np.diff(v) >= 0) for v in units.values())}")
for name in cur:
    t = np.concatenate(list(cur[name].values()))
    print(f"{name}: largest distance of t*30000 from a whole frame: {np.abs(t * 30000 - np.round(t * 30000)).max():.1e} frames")
ids1, ids2 = set(cur["curated 1"]), set(cur["curated 2"])
print("Unit ids in both curations:", len(ids1 & ids2), "; spike trains identical under the same id:",
      sum(np.array_equal(cur["curated 1"][u], cur["curated 2"][u]) for u in ids1 & ids2))
soma = handles["MEArec"]["units/soma_location"][:]
print("MEArec soma locations (um): x", soma[:, 0].min().round(1), "to", soma[:, 0].max().round(1),
      "; y", soma[:, 1].min().round(1), "to", soma[:, 1].max().round(1), "; z", soma[:, 2].min().round(1), "to", soma[:, 2].max().round(1))

curated 1: 258 units, ids 0 to 563, 1190752 spikes, times 0.0032 to 899.9996 s, spikes per unit 143 to 46664 (median 3104), sorted within every unit: True
curated 2: 319 units, ids 0 to 573, 1445570 spikes, times 0.0032 to 899.9996 s, spikes per unit 19 to 46664 (median 3081), sorted within every unit: True
MEArec ground truth: 250 units, ids 0 to 249, 1026644 spikes, times 0.0005 to 599.9987 s, spikes per unit 1454 to 12418 (median 3208), sorted within every unit: True
curated 1: largest distance of t*30000 from a whole frame: 1.9e-09 frames
curated 2: largest distance of t*30000 from a whole frame: 1.9e-09 frames
Unit ids in both curations: 192 ; spike trains identical under the same id: 181
MEArec soma locations (um): x 10.0 to 79.9 ; y -53.9 to 53.0 ; z -1932.5 to 1939.8

The curated spike times lie within the 900 s of the mouse recording and fall on whole frames at 30 kHz, which is consistent with spike times taken from that recording. The two curations share many unit ids with identical spike trains; these are the Kilosort2 clusters both experts kept unchanged, as expected of two curations of one sorting. They are two different files, not one recording deposited twice. The ground-truth spikes lie within the simulated recording's 600 s.

for name in ["curated 1", "curated 2", "mouse412804"]:
    print(f"== {name}: {files[name]}")
    print("\n".join(si.hdf5_tree(handles[name])))
    if "units" in handles[name]:
        print("Units attributes:", {k: si._text(v) if not hasattr(v, "tolist") else v.tolist() for k, v in handles[name]["units"].attrs.items()})
    print()
el = handles["mouse412804"]["general/extracellular_ephys/electrodes"]
cols = [si._text(c) for c in el.attrs["colnames"]]
print("mouse412804 electrode columns:", cols)
print("  a gain or channel_conversion column:", [c for c in cols if "gain" in c or "conversion" in c] or "none")
es = handles["mouse412804"]["acquisition/ElectricalSeries"]
print("  ElectricalSeries members:", list(es), "| data attributes:", {k: si._text(v) for k, v in es["data"].attrs.items()})

== curated 1: sub-mouse412804/sub-mouse412804_ses-20200824T155542.nwb
acquisition/
analysis/
file_create_date (1,) object
general/
general/subject/
general/subject/age () object
general/subject/description () object
general/subject/sex () object
general/subject/species () object
general/subject/subject_id () object
identifier () object
processing/
session_description () object
session_start_time () object
stimulus/
stimulus/presentation/
stimulus/templates/
timestamps_reference_time () object
units/
units/id (258,) int32
units/spike_times (1190752,) float64
units/spike_times_index (258,) uint32
Units attributes: {'colnames': ['spike_times'], 'description': 'Autogenerated by NWBFile', 'namespace': 'core', 'neurodata_type': 'Units', 'object_id': '35c9e03f-cdf3-4f15-8dcf-fc3b17cd2c95'}

== curated 2: sub-mouse412804/sub-mouse412804_ses-20200824T155543.nwb
acquisition/
analysis/
file_create_date (1,) object
general/
general/subject/
general/subject/age () object
general/subject/description () object
general/subject/sex () object
general/subject/species () object
general/subject/subject_id () object
identifier () object
processing/
session_description () object
session_start_time () object
stimulus/
stimulus/presentation/
stimulus/templates/
timestamps_reference_time () object
units/
units/id (319,) int32
units/spike_times (1445570,) float64
units/spike_times_index (319,) uint32
Units attributes: {'colnames': ['spike_times'], 'description': 'Autogenerated by NWBFile', 'namespace': 'core', 'neurodata_type': 'Units', 'object_id': '871289c4-97af-4d08-8421-0abba902667e'}

== mouse412804: sub-mouse412804/sub-mouse412804_ecephys.nwb

acquisition/
acquisition/ElectricalSeries/
acquisition/ElectricalSeries/data (27000000, 248) int16
acquisition/ElectricalSeries/electrodes (248,) int64
acquisition/ElectricalSeries/starting_time () float64
analysis/
file_create_date (1,) object
general/
general/devices/
general/devices/Device/
general/extracellular_ephys/
general/extracellular_ephys/0/
general/extracellular_ephys/electrodes/
general/extracellular_ephys/electrodes/filtering (248,) object
general/extracellular_ephys/electrodes/group (248,) object
general/extracellular_ephys/electrodes/group_name (248,) object
general/extracellular_ephys/electrodes/id (248,) int64
general/extracellular_ephys/electrodes/imp (248,) float64
general/extracellular_ephys/electrodes/location (248,) object
general/extracellular_ephys/electrodes/rel_x (248,) float64
general/extracellular_ephys/electrodes/rel_y (248,) float64
general/extracellular_ephys/electrodes/x (248,) float64
general/extracellular_ephys/electrodes/y (248,) float64
general/extracellular_ephys/electrodes/z (248,) float64
general/subject/
general/subject/age () object
general/subject/description () object
general/subject/sex () object
general/subject/species () object
general/subject/subject_id () object
identifier () object
processing/
session_description () object
session_start_time () object
stimulus/
stimulus/presentation/
stimulus/templates/
timestamps_reference_time () object

mouse412804 electrode columns: ['x', 'y', 'z', 'imp', 'location', 'filtering', 'group', 'group_name', 'rel_x', 'rel_y']
  a gain or channel_conversion column: none
  ElectricalSeries members: ['data', 'electrodes', 'starting_time'] | data attributes: {'conversion': '1e-06', 'resolution': '-1.0', 'unit': 'volts'}

The curated files hold nothing besides a Units table of unit ids and spike times: no acquisition, no processing module, no column or attribute naming the Kilosort2 clusters, and the Units description is the default "Autogenerated by NWBFile". The 174 Kilosort2 units both curators discarded cannot be counted from them. The mouse recording stores int16 values with a conversion of 1e-6 V and gives no gain or channel_conversion per channel, so the scale of its values cannot be checked against the file itself (see below).

env = si.environment()
print("Software the paper's sorters need, on this runner:", env)

Software the paper's sorters need, on this runner: {'matlab': None, 'nvcc': None, 'nvidia-smi': None, 'gpu': None, 'mpiexec': None}

Clock and units of the raw data

The curated spike times and the ground-truth spike times are checked against the raw traces they belong to: the mean trace at each unit's spike frames, shifted by -1.5 to +1.5 ms, should have its deepest trough at lag 0 if the spike times share the recording's clock, and random frames give the reference. Only the first chunk of each channel is read (1.76 s of the mouse recording, 0.59 s of the simulation), because streaming is slow; the time of each read is printed and used below.

rng_clock = np.random.default_rng(34001)
reads = []
lags = np.arange(-45, 46)

mouse = handles["mouse412804"]
n_mouse = mouse["acquisition/ElectricalSeries/data"].chunks[0] * 2
x, secs = si.read_block(mouse, n_mouse)
reads.append({"recording": "mouse412804", "frames": n_mouse, "channels": x.shape[1], "seconds_of_data": n_mouse / 30000, "read_s": secs})
print(f"mouse412804: {x.shape} {x.dtype}, values {x.min()} to {x.max()} (stated conversion 1e-6 V per unit, scale not checked), read in {secs:.1f} s")
y = si.bandpass_cmr(x, 30000.0)
print(f"  band-passed: median channel SD {np.median(y.std(0)):.1f} (value x 1e-6 V, read as uV)")
for name in cur:
    trains = {u: si.to_frames(t) for u, t in cur[name].items()}
    r = si.lag_scan(y, trains, lags, min_spikes=5, rng=rng_clock)
    print(f"  {name}: {len(r['unit'])} units with >= 5 spikes in the block; lag of the deepest trough (frames): median "
          f"{np.median(r['best_lag']):.0f}, 10th to 90th percentile {np.percentile(r['best_lag'], 10):.0f} to "
          f"{np.percentile(r['best_lag'], 90):.0f}, within 3 frames: {np.mean(np.abs(r['best_lag']) <= 3):.2f}; "
          f"trough at lag 0 median {np.median(r['amp_at_0']):.0f} uV, at random frames {np.median(r['amp_random']):.0f} uV")
del x, y

mouse412804: (52736, 248) int16, values -3441 to 1528 (stated conversion 1e-6 V per unit, scale not checked), read in 55.2 s

  band-passed: median channel SD 49.9 (value x 1e-6 V, read as uV)
  curated 1: 84 units with >= 5 spikes in the block; lag of the deepest trough (frames): median 0, 10th to 90th percentile -2 to 3, within 3 frames: 0.88; trough at lag 0 median 455 uV, at random frames 55 uV

  curated 2: 111 units with >= 5 spikes in the block; lag of the deepest trough (frames): median 0, 10th to 90th percentile -2 to 2, within 3 frames: 0.86; trough at lag 0 median 378 uV, at random frames 63 uV

mea = handles["MEArec"]
n_mea = mea["acquisition/ElectricalSeries/data"].chunks[0]
x, secs = si.read_block(mea, n_mea)
reads.append({"recording": "MEArec", "frames": n_mea, "channels": x.shape[1], "seconds_of_data": n_mea / 32000, "read_s": secs})
print(f"MEArec: {x.shape} {x.dtype}, values {x.min():.1f} to {x.max():.1f} uV, read in {secs:.1f} s")
mad_sd = np.median(np.abs(x - np.median(x, 0)), 0) / 0.6745
print(f"  noise SD per channel from the median absolute deviation: median {np.median(mad_sd):.2f} uV, range "
      f"{mad_sd.min():.2f} to {mad_sd.max():.2f} (paper: additive Gaussian noise of 10 uV SD)")
print(f"  channel means: median {np.median(x.mean(0)):.2f} uV, largest in size {np.abs(x.mean(0)).max():.2f} uV")
from scipy.signal import welch
freq, psd = welch(x.astype("float64"), fs=32000.0, nperseg=4096, axis=0)
psd = np.median(psd, axis=1)
bands = [(0, 100), (100, 300), (300, 3000), (3000, 6000), (6000, 16000)]
print("  median power per Hz over channels (uV^2/Hz) by band:",
      "; ".join(f"{lo}-{hi} Hz {psd[(freq >= lo) & (freq < hi)].mean():.4f}" for lo, hi in bands))
y = si.bandpass_cmr(x, 32000.0)
mad_bp = np.median(np.abs(y - np.median(y, 0)), 0) / 0.6745
print(f"  noise SD from the median absolute deviation after the 300 to 6000 Hz band-pass: median {np.median(mad_bp):.2f} uV")
trains = {u: si.to_frames(t, 32000.0) for u, t in gt.items()}
r = si.lag_scan(y, trains, lags, min_spikes=2, rng=rng_clock)
print(f"  ground truth: {len(r['unit'])} units with >= 2 spikes in the block; lag of the deepest trough: median "
      f"{np.median(r['best_lag']):.0f}, within 3 frames: {np.mean(np.abs(r['best_lag']) <= 3):.2f}; trough at lag 0 median "
      f"{np.median(r['amp_at_0']):.0f} uV, at random frames {np.median(r['amp_random']):.0f} uV")
del x, y

MEArec: (18750, 384) float32, values -620.8 to 216.9 uV, read in 74.3 s

  noise SD per channel from the median absolute deviation: median 15.91 uV, range 14.63 to 17.60 (paper: additive Gaussian noise of 10 uV SD)
  channel means: median 0.54 uV, largest in size 2.17 uV

  median power per Hz over channels (uV^2/Hz) by band: 0-100 Hz 0.0973; 100-300 Hz 0.1503; 300-3000 Hz 0.0546; 3000-6000 Hz 0.0210; 6000-16000 Hz 0.0067

  noise SD from the median absolute deviation after the 300 to 6000 Hz band-pass: median 11.98 uV
  ground truth: 213 units with >= 2 spikes in the block; lag of the deepest trough: median 0, within 3 frames: 0.63; trough at lag 0 median 93 uV, at random frames 30 uV

rat = handles["rat c1"]
n_rat = rat["acquisition/ElectricalSeries/data"].chunks[0]
x, secs = si.read_block(rat, n_rat, slice(0, 128))
reads.append({"recording": "rat c1", "frames": n_rat, "channels": x.shape[1], "seconds_of_data": n_rat / 30000, "read_s": secs})
print(f"rat c1, channels 0 to 127: {x.dtype}, values {x.min():.1f} to {x.max():.1f}, per-channel mean "
      f"{x.mean(0).min():.1f} to {x.mean(0).max():.1f}, median channel SD {np.median(x.std(0)):.1f}, read in {secs:.1f} s")
ret = handles["retina"]
n_ret = ret["acquisition/ElectricalSeries/data"].chunks[0]
x, secs = si.read_block(ret, n_ret, slice(0, 256))
reads.append({"recording": "retina", "frames": n_ret, "channels": x.shape[1], "seconds_of_data": n_ret / desc.loc["retina", "rate_hz"], "read_s": secs})
vals, counts = np.unique(x, return_counts=True)
print(f"retina, channels 0 to 255: {x.dtype}, values {x.min()} to {x.max()}, median {np.median(x):.0f}, median channel SD "
      f"{np.median(x.std(0)):.1f}, channels constant over the block: {(x.std(0) == 0).sum()}, samples at {x.max()}: "
      f"{counts[-1]}, read in {secs:.1f} s")
print("retina stated conversion", desc.loc["retina", "conversion"], desc.loc["retina", "unit"], "; offset attribute:", desc.loc["retina", "offset"])
del x

rat c1, channels 0 to 127: float64, values -58.5 to 821.3, per-channel mean 15.3 to 762.3, median channel SD 15.6, read in 26.7 s

retina, channels 0 to 255: uint16, values 1812 to 4096, median 2049, median channel SD 9.7, channels constant over the block: 0, samples at 4096: 157260, read in 36.6 s
retina stated conversion 1e-06 volts ; offset attribute: nan

The mouse curations and the simulated ground truth share the clocks of their recordings: for most units the trough of the mean trace is at lag 0 or within a few frames, and much deeper than at random frames (numbers above). This is shown on the first seconds of each recording only.

The scale of the mouse values is not checked. Read with the stated conversion (1e-6 V per unit), the band-passed noise SD is about 50 uV and the median curated trough about 455 uV (printed above), several times what Neuropixels recordings usually show. Raw Neuropixels 1.0 action-potential counts at about 0.195 uV per bit would give about 10 uV of noise, so the values are consistent with raw counts stored with a conversion of 1e-6, but the file holds no gain to show it (cell above). No claim tested here uses amplitudes; the clock check uses the lag of the trough, not its size.

The simulated traces are not band-passed at 300 Hz: their power per Hz below 300 Hz is above that from 300 to 3000 Hz (printed above), though the channel means are near 0 uV. The noise level estimated from the traces is above the paper's 10 uV before and after a 300 to 6000 Hz band-pass, but both estimates include the overlapping spikes of 250 cells, so they do not show what noise was added; this is left open. The rat recording carries a DC offset that differs between channels. The retina recording is stored as unsigned integers centred near 2048 and reaching 4096, which is consistent with raw ADC counts of a 12-bit converter, not with microvolts; read with the stated conversion (1e-6 V per unit, no offset) every channel would sit near +2 mV, and samples at 4096 (printed above) look like saturation, which the paper's notebook removes with blank_saturation. The gain and offset that turn these counts into volts are not in the file. No claim tested here uses the retina's values.

recordings-as-described (link claim)

The dandiset holds the recordings the paper analyses and the files its notebooks name, with the properties the paper and the notebooks give.

Criterion: each of the five named files is present, and a recording of mouse 412804 with 30 kHz, 15 min (27,000,000 frames) and 246 (paper) or 248 (notebook) channels; the simulated recording lasts 10 min and its ground truth has 250 units; the rat recording has 384 channels at 30 kHz; the Biocam recording has 1024 channels at about 23199 Hz.

named = ["sub-mouse412804_ses-20200824T155542.nwb", "sub-mouse412804_ses-20200824T155543.nwb",
         "sub-MEAREC-250neuron-Neuropixels_ecephys.nwb", "sub-c1_ecephys.nwb", "sub-P29-16-05-14-retina02-left_ecephys.nwb"]
present = set(pinned.path.str.split("/").str[-1])
checks = [(f"file {n} present", n in present, "") for n in named]
m, s, r, b = (desc.loc[k] for k in ["mouse412804", "MEArec", "rat c1", "retina"])
checks += [
    ("mouse 412804 recording present", m.subject_id == "mouse412804", f"{m.path}, subject {m.subject_id}"),
    ("mouse: 30 kHz", m.rate_hz == 30000, m.rate_hz),
    ("mouse: 27,000,000 frames (15 min)", m.frames == 27_000_000, f"{m.frames} frames, {m.duration_s:.1f} s"),
    ("mouse: 246 or 248 channels", m.channels in (246, 248), m.channels),
    ("simulated: 10 min", abs(s.duration_s - 600) < 1, f"{s.duration_s:.2f} s at {s.rate_hz:.0f} Hz"),
    ("simulated: 250 ground-truth units", s.units == 250, s.units),
    ("rat: 384 channels", r.channels == 384, r.channels),
    ("rat: 30 kHz", r.rate_hz == 30000, f"{r.rate_hz} Hz, {r.duration_s:.2f} s (notebook printed 270.0 s)"),
    ("retina: 1024 channels", b.channels == 1024, b.channels),
    ("retina: about 23199 Hz", abs(b.rate_hz - 23199) < 1, f"{b.rate_hz:.3f} Hz, {b.duration_s:.2f} s (notebook printed 135.7 s)"),
]
link = pd.DataFrame(checks, columns=["check", "holds", "value"])
print(link.to_string(index=False))
print("\nAll hold:", bool(link.holds.all()), f"({int(link.holds.sum())} of {len(link)})")

                                                    check  holds                                             value
     file sub-mouse412804_ses-20200824T155542.nwb present   True                                                  
     file sub-mouse412804_ses-20200824T155543.nwb present   True                                                  
file sub-MEAREC-250neuron-Neuropixels_ecephys.nwb present   True                                                  
                          file sub-c1_ecephys.nwb present   True                                                  
  file sub-P29-16-05-14-retina02-left_ecephys.nwb present   True                                                  
                           mouse 412804 recording present   True  sub-mouse412804_ecephys.nwb, subject mouse412804
                                            mouse: 30 kHz   True                                           30000.0
                        mouse: 27,000,000 frames (15 min)   True                        27000000.0 frames, 900.0 s
                               mouse: 246 or 248 channels   True                                             248.0
                                        simulated: 10 min   True                              600.00 s at 32000 Hz
                        simulated: 250 ground-truth units   True                                             250.0
                                        rat: 384 channels   True                                             384.0
                                              rat: 30 kHz   True   30000.0 Hz, 270.01 s (notebook printed 270.0 s)
                                    retina: 1024 channels   True                                            1024.0
                                   retina: about 23199 Hz   True 23199.090 Hz, 135.73 s (notebook printed 135.7 s)

All hold: True (15 of 15)

Result: reproduced. All 15 checks above hold: the five named files are present, the mouse 412804 recording has 248 channels (the notebook's count, not the paper's 246) at 30 kHz for 27,000,000 frames, the simulation lasts 600 s with 250 ground-truth units, the rat recording has 384 channels at 30 kHz, and the retina recording 1024 channels at 23199.09 Hz.

Why the sorter claims are blocked

Every other claim needs the units that the six sorters (HerdingSpikes2, Kilosort2, IronClust, Tridesclous, SpyKING Circus, HDSort) return on these recordings; the dandiset holds none of them, and the claims file fixes the method as running those sorters. Below: what the runner offers, and how long reading each recording whole would take at the streaming rate measured above.

print("MATLAB:", env["matlab"], "| CUDA compiler:", env["nvcc"], "| GPU:", env["gpu"], "| MPI:", env["mpiexec"])
rt = pd.DataFrame(reads)
chans = desc.channels.reindex(rt.recording).values
frames_per_chunk = np.array([desc.loc[k, "chunks"][0] for k in rt.recording])
rt["chunks_read"] = rt.channels * np.ceil(rt.frames / frames_per_chunk)
rt["s_per_chunk"] = rt.read_s / rt.chunks_read
rt["chunks_whole_file"] = np.ceil(desc.frames.reindex(rt.recording).values / frames_per_chunk) * chans
rt["hours_whole_file"] = rt.chunks_whole_file * rt.s_per_chunk / 3600
print(rt.round(3).to_string(index=False))
print(f"\nShortest estimated time to read a whole recording chunk by chunk through h5py: {rt.hours_whole_file.min():.1f} h; notebook limit 0.5 h.")
bulk = si.bulk_read_rate(si.MOUSE, 512 * 2**20)
print(f"\nOne ranged GET of the first {bulk['bytes'] / 2**20:.0f} MB of the mouse recording: {bulk['seconds']:.1f} s, "
      f"{bulk['mb_per_s']:.0f} MB/s (held in memory, discarded)")
sizes = pd.DataFrame(si.list_assets()).set_index("path")["size"]
for name, path in si.RECORDINGS.items():
    print(f"  {name}: {sizes[path] / 2**30:.2f} GiB, read whole at that rate in {sizes[path] / 2**20 / bulk['mb_per_s'] / 60:.1f} min")

MATLAB: None | CUDA compiler: None | GPU: None | MPI: None
  recording  frames  channels  seconds_of_data  read_s  chunks_read  s_per_chunk  chunks_whole_file  hours_whole_file
mouse412804   52736       248            1.758  55.188        496.0        0.111           253952.0             7.849
     MEArec   18750       384            0.586  74.321        384.0        0.194           393216.0            21.140
     rat c1    7911       128            0.264  26.697        128.0        0.209           393216.0            22.782
     retina   12301       256            0.530  36.617        256.0        0.143           262144.0            10.415

Shortest estimated time to read a whole recording chunk by chunk through h5py: 7.8 h; notebook limit 0.5 h.

One ranged GET of the first 512 MB of the mouse recording: 9.5 s, 54 MB/s (held in memory, discarded)

  mouse412804: 12.49 GiB, read whole at that rate in 4.0 min
  MEArec: 27.50 GiB, read whole at that rate in 8.7 min
  rat c1: 23.21 GiB, read whole at that rate in 7.4 min
  retina: 6.03 GiB, read whole at that rate in 1.9 min

Kilosort2, IronClust and HDSort need MATLAB, and Kilosort2 a CUDA GPU; none is on this runner and none can be installed from a lock file. While preparing (not run here), Tridesclous 1.6.0 and SpyKING Circus 0.9.7 installed with pip on Python 3.13, HerdingSpikes 0.3.7 failed to build, and SpyKING Circus has no MPI to run with. Most sorter claims need the output of all six sorters (or, for the retina, of four, three of them in MATLAB); fig2e-ks2-unit-types, fig2e-ks2-curated-unit-types and fig3s2-snr-does-not-separate-ks2-fp need Kilosort2 alone, which needs MATLAB and a GPU. So these claims are blocked as not_runnable whatever the Python sorters could do.

Reading time does not add to the block. Read chunk by chunk through h5py, as above, a whole recording would take hours, but one bulk ranged read of the mouse file runs at the rate printed above, so each recording could be read whole in minutes (printed above). The block rests on the missing MATLAB and GPU alone.

fig4a-curators-agree (partially reproduced; the 174 untested)

The two experts' manual curations of the Kilosort2 output agree on most curated units (226 of 351, 64.2%), and 174 Kilosort2 units are discarded by both.

The claim's count of 174 needs the Kilosort2 run, which is blocked (fig1b-unit-counts); the curated files hold no link to the Kilosort2 clusters (cell above), so it cannot be had from them. Its other numbers come from the two curated files alone under the criterion's second reading: compare_two_sorters(curated 1, curated 2) with defaults, 226 = matched pairs, 351 = the two counts summed less the matched pairs. The criterion gives reproduced when all four values are within tolerance under one reading, and partially reproduced when 226, 351 and 64.2% are and 174 is not. The first reading (compare_multiple_sorters with the Kilosort2 run) cannot be run.

added = si.legacy_numpy_aliases()
import spikecomparison, spikeextractors
print("spikecomparison", spikecomparison.__version__, "| spikeextractors", spikeextractors.__version__, "| NumPy aliases restored:", added)
tr1 = {u: si.to_frames(t) for u, t in cur["curated 1"].items()}
tr2 = {u: si.to_frames(t) for u, t in cur["curated 2"].items()}
comp = si.compare_two(tr1, tr2)
match = comp.hungarian_match_12
n1, n2, k = len(tr1), len(tr2), int((match != -1).sum())
print(f"curated 1: {n1} units; curated 2: {n2} units")
print(f"matched pairs (both curators): {k}; units of either curator: {n1 + n2 - k}; ratio {100 * k / (n1 + n2 - k):.2f}%")
print("expected: 226 +- 22.6; 351 +- 35.1; 64.2 +- 6.42%")
scores = np.array([comp.agreement_scores.loc[a, b] for a, b in match.items() if b != -1])
print(f"agreement scores of matched pairs: min {scores.min():.3f}, median {np.median(scores):.3f}, equal to 1: {(scores == 1).sum()}")
print("\nDependence on the minimum score (criterion: the default, 0.5):")
for ms in [0.3, 0.4, 0.5, 0.6, 0.7, 0.8]:
    c = si.compare_two(tr1, tr2, match_score=ms)
    kk = int((c.hungarian_match_12 != -1).sum())
    print(f"  match_score {ms}: both {kk}, either {n1 + n2 - kk}, {100 * kk / (n1 + n2 - kk):.2f}%")
print(f"\nrun time so far {(time.time() - started) / 60:.1f} min")

spikecomparison 0.2.6 | spikeextractors 0.8.4 | NumPy aliases restored: ['float', 'int', 'Inf', 'product', 'in1d']

curated 1: 258 units; curated 2: 319 units
matched pairs (both curators): 226; units of either curator: 351; ratio 64.39%
expected: 226 +- 22.6; 351 +- 35.1; 64.2 +- 6.42%
agreement scores of matched pairs: min 0.510, median 1.000, equal to 1: 147

Dependence on the minimum score (criterion: the default, 0.5):

  match_score 0.3: both 231, either 346, 66.76%

  match_score 0.4: both 229, either 348, 65.80%

  match_score 0.5: both 226, either 351, 64.39%

  match_score 0.6: both 222, either 355, 62.54%

  match_score 0.7: both 216, either 361, 59.83%

  match_score 0.8: both 209, either 368, 56.79%

run time so far 13.7 min

Result: partially reproduced. The curators share 226 matched units of 351 (64.39%). These are the paper's 226 and 351 exactly, and 64.39% is within 6.42 of the paper's 64.2%. With a minimum score of 0.3 to 0.8 the shared count stays inside its tolerance (printed above). The ratio stays inside its tolerance up to 0.7 but at 0.8 falls below the lower bound of 57.78%; the criterion fixes the default of 0.5, so this does not decide the outcome. The 174 Kilosort2 units both curators discarded could not be counted, as that needs the Kilosort2 run. Under the criterion the claim is reproduced if the 174 holds and partially reproduced if it does not, so with the 174 untested it is partially reproduced: at least this, and reproduced only if a Kilosort2 run, which this attempt could not make, confirmed the 174.

Summary

Claim	Outcome
recordings-as-described (link)	reproduced
fig4a-curators-agree	partially reproduced (226 of 351 curated units shared, as in the paper; the 174 needs the Kilosort2 run and is untested)
the 22 other claims (fig1b to fig4d)	blocked, not_runnable

How the analysis differs from the paper's. No sorter was run. The paper ran six sorters, three of them in MATLAB, on the four recordings; this attempt could only check the recordings and compare the two curations with each other. The curations were compared with spikecomparison 0.2.6, the version the paper's notebooks ran, on NumPy 2 with the removed NumPy aliases restored.

What was found about the dataset.

• The dandiset holds the four recordings and the two curations the paper names, with the stated rates, channel counts, durations and number of ground-truth units. It holds no sorter output, so the paper's comparisons can be repeated only by running the sorters again.
• The mouse recording has 248 channels; the paper says 246 active channels. Which two are inactive is not marked.
• The curated sortings carry no rate, no unit labels and no link to the Kilosort2 units they were curated from; their spike times are consistent with the mouse recording's clock (trough at lag 0).
• The retina recording is stored as unsigned counts centred near 2048 with a stated conversion of 1e-6 V and no offset, which would give a DC level of about 2 mV; the gain and offset of the Biocam's converter are not in the file.
• Every raw series is chunked one channel per chunk, so reading a stretch of all channels through h5py is very slow (table above); one bulk read of a whole file is fast.
• The mouse recording's int16 values, read with its conversion of 1e-6 V, give a band-passed noise SD of about 50 uV, which is consistent with raw Neuropixels counts; the file gives no gain to show whether the values are counts or microvolts.
• The curated files hold unit ids and spike times only; nothing in them links a unit to a Kilosort2 cluster.
• The archive's species summary lists only house mouse, though one recording is from a rat.
• The paper's 64.2% for 226 of 351 is 64.39%.

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.