Reproducing Ramachandran et al. (2025) from dandiset 000945¶
Paper. Ramachandran S, Gao H, Yttri E, Yu K, He B (2025). Parameter-dependent cell-type specific effects of transcranial focused ultrasound stimulation in an awake head-fixed rodent model. J Neural Eng 22(2) 026022. doi:10.1088/1741-2552/adbb1f. The claims file was written from the version of record as served by Europe PMC (PMC11920894): main text, Methods, figure legends and Discussion; the supplementary materials were not supplied, and no claim rests on them.
Dandiset. DANDI:000945, version 0.250213.2212. This is the version the paper cites (reference 50, doi:10.48324/dandi.000945/0.250213.2212). The four later published versions (0.250711.1746, 0.250829.0425, 0.250829.1250) and the draft hold the same 75 files with the same blobs and sizes, as the archive's asset lists show, so the choice does not change the data analysed.
What the dandiset holds. Sorted units (spike times; a cell-type label in 60 of 75 files) and
the trial table of each recording. There are no waveforms, no raw or filtered voltage, no
electrode positions and no sham recordings. The condition of each recording is coded only in the
file's identifier (<rat>_<PRF>_<TBD>[_2][_anes]).
| Claim | Tested here |
|---|---|
| fig3-cell-type-sorting | Blocked: the k-means features need waveforms, which the files lack; the counts of the dataset's labels are shown |
| fig4-rsu-time-locked-both-states | Blocked: the criterion compares with sham, and there is no sham |
| fig4hi-fsu-delayed-latency | Tested |
| fig4h-awake-latency-wider | Tested; blocked, as the claims file's readings D1 and D2 give different results |
| fig5a-rsu-time-locked-prf-awake | Tested; the sham parts cannot be |
| fig5c-rsu-time-locked-prf-anesth | Tested (two labelled rats); the sham parts cannot be |
| rsu-time-locked-anesthesia-factor | Tested |
| fig5bd-fsu-time-locked-no-prf | Tested |
| fig5e-rsu-delayed-prf-awake | Tested; the sham part cannot be |
| fig5f-fsu-delayed-prf-awake | Tested; the sham part cannot be |
| fig5gh-delayed-prf-anesth | Tested |
| delayed-anesthesia-factor | Tested |
| fig6ad-time-adjacent-batch | Blocked: the second batch of five rats is not in the dandiset |
| fig6ef-rsu-ispta-rise-then-fall | Tested; blocked, as the claims file's readings W1 and W2 give different results |
| fig6fh-fsu-ispta | Tested |
| fig5a-1500-highest-time-locked-mean | Tested |
| fig7a-rsu-time-locked-dc | Tested; the sham part cannot be |
| fig7a-60-highest-time-locked-mean | Tested |
| fig7b-fsu-time-locked-dc | Tested |
| fig7c-rsu-delayed-dc | Blocked: the paper's own 30% DC recordings are not in the dandiset, and the substitutes decide the outcome |
| fig7d-fsu-delayed-dc | Tested; the sham part cannot be |
| sham-no-auditory-response | Blocked: no sham recordings, and no test named |
import time
import warnings
import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import remfile
import statsmodels.api as sm
from scipy import stats
from statsmodels.formula.api import ols
import tfus
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
assets = tfus.asset_list()
sessions = tfus.load_all(assets)
print(f"loaded {len(sessions)} files from dandiset {tfus.DANDISET_ID} version {tfus.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
loaded 75 files from dandiset 000945 version 0.250213.2212 in 45 s
Data check¶
Each check of the metadata against the values is stated under its output.
inv = pd.DataFrame([{
"path": s["path"], "subject_id": s["subject_id"], "identifier": s["identifier"], "rat": s["rat"],
"state": s["state"], "prf": s["prf"], "tbd": s["tbd"], "repeat": s["repeat"],
"description": s["session_description"], "age": s["age"], "start": s["session_start"][:19],
"units": s["n_units"], "labelled": int((~np.isnan(s["label"])).sum()), "trials": s["n_trials"],
"units_cols": " ".join(s["units_colnames"]), "trial_cols": " ".join(s["trials_colnames"]),
"acq_proc": " ".join(s["acquisition"] + s["processing"]) or "-"} for s in sessions])
print("Files:", len(inv), "; subject folders:", inv.subject_id.nunique(), sorted(inv.subject_id.unique()))
print("Rats named in identifiers:", sorted(inv.rat.unique()))
print("Files whose subject_id differs from the rat in the identifier:")
print(inv[inv.subject_id != inv.rat][["path", "identifier"]].to_string(index=False))
print("\nUnits table columns:", inv.units_cols.value_counts().to_dict())
print("Trials table columns:", inv.trial_cols.value_counts().to_dict())
print("Acquisition or processing modules:", inv.acq_proc.value_counts().to_dict())
print("Ages:", inv.age.value_counts().to_dict())
print("Session descriptions by state in the identifier:")
print(pd.crosstab(inv.description, inv.state).to_string())
print("\nTrials per file:", inv.trials.value_counts().to_dict())
print(inv[inv.trials != 500][["identifier", "trials"]].to_string(index=False))
print("\nLabels: files with celltype_label", int((inv.labelled > 0).sum()), "of", len(inv),
"; units in files without it:", int(inv[inv.labelled == 0].units.sum()))
print("\nDandiset description at this version:")
print(tfus.dandiset_description())
print("\nEvery file, by the rat its identifier names and its session start:")
print(inv.sort_values(["rat", "start"])[["rat", "identifier", "subject_id", "start", "path"]].to_string(index=False))
print("Every identifier parses as <rat>_<PRF>_<TBD>[_2][_anes|_anes2]: True (parse_identifier raised no error)")
Files: 75 ; subject folders: 10 ['BH497', 'BH498', 'BH506', 'BH508', 'BH512', 'BH514', 'BH519', 'BH521', 'BH526', 'BH527']
Rats named in identifiers: ['BH498', 'BH506', 'BH508', 'BH512', 'BH519', 'BH521', 'BH526', 'BH527']
Files whose subject_id differs from the rat in the identifier:
path identifier
sub-BH497/sub-BH497_ses-20240310T143729_ecephys.nwb BH498_3000_200_anes
sub-BH497/sub-BH497_ses-20240310T145053_ecephys.nwb BH498_1500_200_anes
sub-BH497/sub-BH497_ses-20240310T145814_ecephys.nwb BH498_300_200_anes
sub-BH497/sub-BH497_ses-20240310T150748_ecephys.nwb BH498_30_200_anes
sub-BH497/sub-BH497_ses-20240310T152652_ecephys.nwb BH498_4500_200_anes
sub-BH514/sub-BH514_ses-20240309T183140_ecephys.nwb BH512_1500_200_anes
sub-BH514/sub-BH514_ses-20240309T183856_ecephys.nwb BH512_300_200_anes
sub-BH514/sub-BH514_ses-20240309T184852_ecephys.nwb BH512_4500_200_anes
sub-BH514/sub-BH514_ses-20240309T185729_ecephys.nwb BH512_3000_200_anes
Units table columns: {'spike_times celltype_label': 60, 'spike_times': 15}
Trials table columns: {'start_time stop_time': 75}
Acquisition or processing modules: {'-': 75}
Ages: {'P24W': 75}
Session descriptions by state in the identifier:
state anesthetized awake
description
Awake Model Anes S1 Stimulation by tFUS 15 0
Awake S1 Stimulation by tFUS 9 51
Trials per file: {500: 70, 203: 2, 238: 1, 328: 1, 339: 1}
identifier trials
BH498_30_200 203
BH508_1500_200_2 238
BH508_3000_200 328
BH508_1500_200 339
BH508_30_200_2 203
Labels: files with celltype_label 60 of 75 ; units in files without it: 960
Dandiset description at this version:
In these recordings, we tested different PRFs (pulse repetition frequencies) of ultrasound stimulation using a 128-element random array ultrasound transducer to stimulate the somatosensory cortex of awake head-fixed rats. Chronic electrophysiological recordings were acquired using 32-channel NeuroNexus electrodes (Model: A1x32-Poly3-10mm-50-177-Z32), chronically implanted into the rat brain. Ultrasound stimulation is delivered every 2.5 seconds with a 10% jitter, and each recording has 500 trials. The PRFs tested were 30 Hz, 300 Hz, 1500 Hz, 3000 Hz, and 4500 Hz, each with a 200 microsecond pulse duration and a 67 ms ultrasound duration. Anesthetized files were performed under 2% isoflurane anesthesia for comparison. File Information All 10 subjects were male rats, implanted with their chronic electrode at 6 months of age and then recordings taken first at 8-10 months, and then some repeats taken at 12 months. Within each subject's folder are recordings for the different PRFs. Most subjects have 5 recordings within, one for each PRF. Some subjects have duplicate recordings taken a few months after the original ones. A few recordings were not included due to excessive noise in the recordings. Files are named in the format SubjectName_PRF_PulseDuration. Each file contains spike time data with the cell type labels included for each neurons, as well as time series data for the onset of each trial of ultrasound stimulation. Every file, by the rat its identifier names and its session start: rat identifier subject_id start path BH498 BH498_3000_200 BH498 2023-09-06T18:09:13 sub-BH498/sub-BH498_ses-20230906T180913_ecephys.nwb BH498 BH498_4500_200_2 BH498 2024-03-05T19:18:33 sub-BH498/sub-BH498_ses-20240305T191833_ecephys.nwb BH498 BH498_1500_200_2 BH498 2024-03-05T20:06:06 sub-BH498/sub-BH498_ses-20240305T200606_ecephys.nwb BH498 BH498_3000_200_2 BH498 2024-03-05T20:45:52 sub-BH498/sub-BH498_ses-20240305T204552_ecephys.nwb BH498 BH498_1500_200 BH498 2024-03-06T12:42:27 sub-BH498/sub-BH498_ses-20240306T124227_ecephys.nwb BH498 BH498_30_200 BH498 2024-03-06T17:58:05 sub-BH498/sub-BH498_ses-20240306T175805_ecephys.nwb BH498 BH498_300_200 BH498 2024-03-06T18:45:36 sub-BH498/sub-BH498_ses-20240306T184536_ecephys.nwb BH498 BH498_4500_200 BH498 2024-03-06T23:13:17 sub-BH498/sub-BH498_ses-20240306T231317_ecephys.nwb BH498 BH498_3000_200_anes BH497 2024-03-10T14:37:29 sub-BH497/sub-BH497_ses-20240310T143729_ecephys.nwb BH498 BH498_3000_200_anes BH498 2024-03-10T14:37:29 sub-BH498/sub-BH498_ses-20240310T143729_ecephys.nwb BH498 BH498_1500_200_anes BH497 2024-03-10T14:50:53 sub-BH497/sub-BH497_ses-20240310T145053_ecephys.nwb BH498 BH498_1500_200_anes BH498 2024-03-10T14:50:53 sub-BH498/sub-BH498_ses-20240310T145053_ecephys.nwb BH498 BH498_300_200_anes BH497 2024-03-10T14:58:14 sub-BH497/sub-BH497_ses-20240310T145814_ecephys.nwb BH498 BH498_300_200_anes BH498 2024-03-10T14:58:14 sub-BH498/sub-BH498_ses-20240310T145814_ecephys.nwb BH498 BH498_30_200_anes BH497 2024-03-10T15:07:48 sub-BH497/sub-BH497_ses-20240310T150748_ecephys.nwb BH498 BH498_30_200_anes BH498 2024-03-10T15:07:48 sub-BH498/sub-BH498_ses-20240310T150748_ecephys.nwb BH498 BH498_4500_200_anes BH497 2024-03-10T15:26:52 sub-BH497/sub-BH497_ses-20240310T152652_ecephys.nwb BH498 BH498_4500_200_anes BH498 2024-03-10T15:26:52 sub-BH498/sub-BH498_ses-20240310T152652_ecephys.nwb BH506 BH506_1500_200 BH506 2023-10-18T18:16:09 sub-BH506/sub-BH506_ses-20231018T181609_ecephys.nwb BH506 BH506_4500_200 BH506 2023-10-26T15:51:06 sub-BH506/sub-BH506_ses-20231026T155106_ecephys.nwb BH506 BH506_3000_200 BH506 2023-10-26T17:27:06 sub-BH506/sub-BH506_ses-20231026T172706_ecephys.nwb BH506 BH506_300_200 BH506 2023-10-26T22:20:22 sub-BH506/sub-BH506_ses-20231026T222022_ecephys.nwb BH506 BH506_30_200 BH506 2023-10-26T23:11:03 sub-BH506/sub-BH506_ses-20231026T231103_ecephys.nwb BH506 BH506_3000_200_2 BH506 2024-03-08T15:06:12 sub-BH506/sub-BH506_ses-20240308T150612_ecephys.nwb BH506 BH506_4500_200_2 BH506 2024-03-08T20:36:35 sub-BH506/sub-BH506_ses-20240308T203635_ecephys.nwb BH506 BH506_1500_40 BH506 2024-03-09T15:43:17 sub-BH506/sub-BH506_ses-20240309T154317_ecephys.nwb BH506 BH506_1500_400 BH506 2024-03-09T16:25:19 sub-BH506/sub-BH506_ses-20240309T162519_ecephys.nwb BH508 BH508_1500_200_2 BH508 2023-08-18T23:26:09 sub-BH508/sub-BH508_ses-20230818T232609_ecephys.nwb BH508 BH508_3000_200 BH508 2023-08-27T14:08:43 sub-BH508/sub-BH508_ses-20230827T140843_ecephys.nwb BH508 BH508_4500_200 BH508 2023-08-27T17:22:57 sub-BH508/sub-BH508_ses-20230827T172257_ecephys.nwb BH508 BH508_300_200 BH508 2023-08-27T23:54:25 sub-BH508/sub-BH508_ses-20230827T235425_ecephys.nwb BH508 BH508_30_200 BH508 2023-10-12T20:07:51 sub-BH508/sub-BH508_ses-20231012T200751_ecephys.nwb BH508 BH508_1500_200 BH508 2023-10-22T18:52:10 sub-BH508/sub-BH508_ses-20231022T185210_ecephys.nwb BH508 BH508_30_200_2 BH508 2024-03-07T00:23:54 sub-BH508/sub-BH508_ses-20240307T002354_ecephys.nwb BH508 BH508_300_200_2 BH508 2024-03-07T13:19:48 sub-BH508/sub-BH508_ses-20240307T131948_ecephys.nwb BH508 BH508_4500_200_2 BH508 2024-03-07T15:19:10 sub-BH508/sub-BH508_ses-20240307T151910_ecephys.nwb BH508 BH508_3000_200_2 BH508 2024-03-07T16:10:26 sub-BH508/sub-BH508_ses-20240307T161026_ecephys.nwb BH512 BH512_1500_200 BH512 2023-10-13T12:59:04 sub-BH512/sub-BH512_ses-20231013T125904_ecephys.nwb BH512 BH512_3000_200 BH512 2023-10-13T14:49:25 sub-BH512/sub-BH512_ses-20231013T144925_ecephys.nwb BH512 BH512_4500_200 BH512 2023-10-17T11:12:24 sub-BH512/sub-BH512_ses-20231017T111224_ecephys.nwb BH512 BH512_300_200 BH512 2023-10-17T13:39:12 sub-BH512/sub-BH512_ses-20231017T133912_ecephys.nwb BH512 BH512_30_200 BH512 2023-10-17T15:00:24 sub-BH512/sub-BH512_ses-20231017T150024_ecephys.nwb BH512 BH512_1500_400_anes BH512 2024-03-09T17:34:51 sub-BH512/sub-BH512_ses-20240309T173451_ecephys.nwb BH512 BH512_1500_40_anes BH512 2024-03-09T17:42:22 sub-BH512/sub-BH512_ses-20240309T174222_ecephys.nwb BH512 BH512_1500_200_anes BH512 2024-03-09T18:31:40 sub-BH512/sub-BH512_ses-20240309T183140_ecephys.nwb BH512 BH512_1500_200_anes BH514 2024-03-09T18:31:40 sub-BH514/sub-BH514_ses-20240309T183140_ecephys.nwb BH512 BH512_300_200_anes BH512 2024-03-09T18:38:56 sub-BH512/sub-BH512_ses-20240309T183856_ecephys.nwb BH512 BH512_300_200_anes BH514 2024-03-09T18:38:56 sub-BH514/sub-BH514_ses-20240309T183856_ecephys.nwb BH512 BH512_4500_200_anes BH512 2024-03-09T18:48:52 sub-BH512/sub-BH512_ses-20240309T184852_ecephys.nwb BH512 BH512_4500_200_anes BH514 2024-03-09T18:48:52 sub-BH514/sub-BH514_ses-20240309T184852_ecephys.nwb BH512 BH512_3000_200_anes BH512 2024-03-09T18:57:29 sub-BH512/sub-BH512_ses-20240309T185729_ecephys.nwb BH512 BH512_3000_200_anes BH514 2024-03-09T18:57:29 sub-BH514/sub-BH514_ses-20240309T185729_ecephys.nwb BH512 BH512_1500_200_anes2 BH512 2024-05-03T22:40:59 sub-BH512/sub-BH512_ses-20240503T224059_ecephys.nwb BH519 BH519_3000_200 BH519 2023-10-13T12:16:04 sub-BH519/sub-BH519_ses-20231013T121604_ecephys.nwb BH519 BH519_4500_200 BH519 2023-10-18T18:48:27 sub-BH519/sub-BH519_ses-20231018T184827_ecephys.nwb BH519 BH519_1500_200 BH519 2023-10-22T16:39:21 sub-BH519/sub-BH519_ses-20231022T163921_ecephys.nwb BH519 BH519_300_200 BH519 2023-10-26T14:46:09 sub-BH519/sub-BH519_ses-20231026T144609_ecephys.nwb BH519 BH519_30_200 BH519 2023-10-26T15:26:35 sub-BH519/sub-BH519_ses-20231026T152635_ecephys.nwb BH521 BH521_3000_200_anes BH521 2024-02-06T14:57:44 sub-BH521/sub-BH521_ses-20240206T145744_ecephys.nwb BH521 BH521_300_200_anes BH521 2024-02-06T16:26:53 sub-BH521/sub-BH521_ses-20240206T162653_ecephys.nwb BH521 BH521_4500_200_anes BH521 2024-02-06T22:28:55 sub-BH521/sub-BH521_ses-20240206T222855_ecephys.nwb BH521 BH521_1500_400 BH521 2024-03-09T14:29:51 sub-BH521/sub-BH521_ses-20240309T142951_ecephys.nwb BH521 BH521_1500_40 BH521 2024-03-09T14:56:24 sub-BH521/sub-BH521_ses-20240309T145624_ecephys.nwb BH526 BH526_1500_200 BH526 2023-10-29T17:24:56 sub-BH526/sub-BH526_ses-20231029T172456_ecephys.nwb BH526 BH526_30_200 BH526 2023-10-29T18:22:26 sub-BH526/sub-BH526_ses-20231029T182226_ecephys.nwb BH526 BH526_300_200 BH526 2023-10-29T18:44:25 sub-BH526/sub-BH526_ses-20231029T184425_ecephys.nwb BH526 BH526_3000_200 BH526 2023-10-29T19:10:50 sub-BH526/sub-BH526_ses-20231029T191050_ecephys.nwb BH526 BH526_4500_200 BH526 2023-10-29T19:25:29 sub-BH526/sub-BH526_ses-20231029T192529_ecephys.nwb BH526 BH526_1500_400 BH526 2024-01-30T18:43:18 sub-BH526/sub-BH526_ses-20240130T184318_ecephys.nwb BH526 BH526_1500_40 BH526 2024-01-30T19:52:06 sub-BH526/sub-BH526_ses-20240130T195206_ecephys.nwb BH527 BH527_4500_200 BH527 2023-10-29T23:33:34 sub-BH527/sub-BH527_ses-20231029T233334_ecephys.nwb BH527 BH527_3000_200 BH527 2023-10-30T00:03:25 sub-BH527/sub-BH527_ses-20231030T000325_ecephys.nwb BH527 BH527_1500_200 BH527 2023-12-01T16:28:33 sub-BH527/sub-BH527_ses-20231201T162833_ecephys.nwb BH527 BH527_300_200 BH527 2024-03-04T22:28:14 sub-BH527/sub-BH527_ses-20240304T222814_ecephys.nwb BH527 BH527_30_200 BH527 2024-03-05T16:58:12 sub-BH527/sub-BH527_ses-20240305T165812_ecephys.nwb Every identifier parses as <rat>_<PRF>_<TBD>[_2][_anes|_anes2]: True (parse_identifier raised no error)
There are no waveforms, voltage traces, electrode positions or processing modules: units hold only spike times and, in 60 files, a label ("1 = RSU, 2 = FSU"). Every file says age P24W (about 5.5 months), where the Methods give 7, 9 and 10 months for the three series. The 15 files without labels are all anesthetized recordings, and they are the only files whose description says "Awake Model Anes"; the labelled anesthetized copies (subject folders BH497 and BH514) are described as "Awake S1 Stimulation". The identifiers of those copies name rats BH498 and BH512. Five files hold fewer than the 500 trials the Methods describe.
The dandiset's description explains the _2 files: recordings were "taken first at 8-10 months,
and then some repeats taken at 12 months", and "Some subjects have duplicate recordings taken a
few months after the original ones". The session start times (table above) do not always agree
with the suffix; see the check after the next one.
dup = tfus.find_duplicates(sessions)
print("Groups of files with identical trial times:")
print(dup[dup.group_size > 1].sort_values(["duplicate_of", "path"], na_position="first")
[["path", "identifier", "kept", "duplicate_of"]].to_string(index=False))
by_path = {s["path"]: s for s in sessions}
rows = []
for r in dup[~dup.kept].itertuples():
a, b = by_path[r.path], by_path[r.duplicate_of]
same_counts = a["n_units"] == b["n_units"] and np.array_equal(a["n_spikes"], b["n_spikes"])
same_first = np.array_equal(a["first_spikes"], b["first_spikes"])
rows.append({"file": r.identifier + " (" + r.subject_id + ")", "kept": b["identifier"] + " (" + b["subject_id"] + ")",
"units": f'{a["n_units"]}/{b["n_units"]}', "spikes": f'{a["n_spikes"].sum()}/{b["n_spikes"].sum()}',
"same spike count per unit": same_counts, "same first 5000 spike times": same_first})
print(pd.DataFrame(rows).to_string(index=False))
def shared_spikes(p1, p2):
out = []
for p in (p1, p2):
url = next(a["url"] for a in assets if a["path"] == p)
with h5py.File(remfile.File(url), "r") as f:
out.append(np.unique(np.round(f["units/spike_times"][:], 5)))
return len(out[0]), len(out[1]), len(np.intersect1d(out[0], out[1]))
n1, n2, n12 = shared_spikes("sub-BH498/sub-BH498_ses-20240306T175805_ecephys.nwb",
"sub-BH508/sub-BH508_ses-20240307T002354_ecephys.nwb")
print(f"\nBH498_30_200 and BH508_30_200_2: distinct spike times {n1} and {n2}, shared to 10 us: {n12}")
for p in ("sub-BH498/sub-BH498_ses-20240306T175805_ecephys.nwb", "sub-BH508/sub-BH508_ses-20240307T002354_ecephys.nwb"):
print(f' {by_path[p]["identifier"]} ({by_path[p]["subject_id"]}): session start {by_path[p]["session_start"][:19]}')
kept_paths = set(dup[dup.kept].path)
kept = [s for s in sessions if s["path"] in kept_paths]
print("Recordings counted once:", len(kept), "of", len(sessions), "files")
Groups of files with identical trial times:
path identifier kept duplicate_of
sub-BH497/sub-BH497_ses-20240310T143729_ecephys.nwb BH498_3000_200_anes True NaN
sub-BH497/sub-BH497_ses-20240310T145053_ecephys.nwb BH498_1500_200_anes True NaN
sub-BH497/sub-BH497_ses-20240310T145814_ecephys.nwb BH498_300_200_anes True NaN
sub-BH497/sub-BH497_ses-20240310T150748_ecephys.nwb BH498_30_200_anes True NaN
sub-BH497/sub-BH497_ses-20240310T152652_ecephys.nwb BH498_4500_200_anes True NaN
sub-BH498/sub-BH498_ses-20240306T175805_ecephys.nwb BH498_30_200 True NaN
sub-BH514/sub-BH514_ses-20240309T183140_ecephys.nwb BH512_1500_200_anes True NaN
sub-BH514/sub-BH514_ses-20240309T183856_ecephys.nwb BH512_300_200_anes True NaN
sub-BH514/sub-BH514_ses-20240309T184852_ecephys.nwb BH512_4500_200_anes True NaN
sub-BH514/sub-BH514_ses-20240309T185729_ecephys.nwb BH512_3000_200_anes True NaN
sub-BH498/sub-BH498_ses-20240310T143729_ecephys.nwb BH498_3000_200_anes False sub-BH497/sub-BH497_ses-20240310T143729_ecephys.nwb
sub-BH498/sub-BH498_ses-20240310T145053_ecephys.nwb BH498_1500_200_anes False sub-BH497/sub-BH497_ses-20240310T145053_ecephys.nwb
sub-BH498/sub-BH498_ses-20240310T145814_ecephys.nwb BH498_300_200_anes False sub-BH497/sub-BH497_ses-20240310T145814_ecephys.nwb
sub-BH498/sub-BH498_ses-20240310T150748_ecephys.nwb BH498_30_200_anes False sub-BH497/sub-BH497_ses-20240310T150748_ecephys.nwb
sub-BH498/sub-BH498_ses-20240310T152652_ecephys.nwb BH498_4500_200_anes False sub-BH497/sub-BH497_ses-20240310T152652_ecephys.nwb
sub-BH508/sub-BH508_ses-20240307T002354_ecephys.nwb BH508_30_200_2 False sub-BH498/sub-BH498_ses-20240306T175805_ecephys.nwb
sub-BH512/sub-BH512_ses-20240309T183140_ecephys.nwb BH512_1500_200_anes False sub-BH514/sub-BH514_ses-20240309T183140_ecephys.nwb
sub-BH512/sub-BH512_ses-20240503T224059_ecephys.nwb BH512_1500_200_anes2 False sub-BH514/sub-BH514_ses-20240309T183140_ecephys.nwb
sub-BH512/sub-BH512_ses-20240309T183856_ecephys.nwb BH512_300_200_anes False sub-BH514/sub-BH514_ses-20240309T183856_ecephys.nwb
sub-BH512/sub-BH512_ses-20240309T184852_ecephys.nwb BH512_4500_200_anes False sub-BH514/sub-BH514_ses-20240309T184852_ecephys.nwb
sub-BH512/sub-BH512_ses-20240309T185729_ecephys.nwb BH512_3000_200_anes False sub-BH514/sub-BH514_ses-20240309T185729_ecephys.nwb
file kept units spikes same spike count per unit same first 5000 spike times
BH498_3000_200_anes (BH498) BH498_3000_200_anes (BH497) 64/64 1200943/1200943 True True
BH498_1500_200_anes (BH498) BH498_1500_200_anes (BH497) 64/64 1145377/1145377 True True
BH498_300_200_anes (BH498) BH498_300_200_anes (BH497) 64/64 1463898/1463898 True True
BH498_30_200_anes (BH498) BH498_30_200_anes (BH497) 64/64 1582396/1582396 True True
BH498_4500_200_anes (BH498) BH498_4500_200_anes (BH497) 64/64 1406942/1406942 True True
BH508_30_200_2 (BH508) BH498_30_200 (BH498) 26/31 32176/35497 False False
BH512_1500_200_anes (BH512) BH512_1500_200_anes (BH514) 64/64 926883/926883 True True
BH512_1500_200_anes2 (BH512) BH512_1500_200_anes (BH514) 64/64 905554/926883 False True
BH512_300_200_anes (BH512) BH512_300_200_anes (BH514) 64/64 1121052/1121052 True True
BH512_4500_200_anes (BH512) BH512_4500_200_anes (BH514) 64/64 1000775/1000775 True True
BH512_3000_200_anes (BH512) BH512_3000_200_anes (BH514) 64/64 904524/904524 True True
BH498_30_200 and BH508_30_200_2: distinct spike times 30188 and 27391, shared to 10 us: 18860 BH498_30_200 (BH498): session start 2024-03-06T17:58:05 BH508_30_200_2 (BH508): session start 2024-03-07T00:23:54 Recordings counted once: 64 of 75 files
Eleven files are second copies and are not kept. The five anesthetized BH498 recordings are each deposited twice,
under BH498 without labels and under BH497 with labels; the BH512 anesthetized recordings at
300, 1500, 3000 and 4500 Hz are deposited under BH512 without labels and under BH514 with labels,
and BH512_1500_200_anes2 is a third copy of the 1500 Hz recording with the same trial times and
a different sort (21,329 fewer spikes). In each of the nine pairs the trial times, the number of spikes of
every unit and the first spike times are identical; the third copy shares the trial times and the
first spike times, but not the spike counts (905,554 against 926,883). The labelled copy is kept, and the rat is
taken from the identifier (BH498, BH512), which the unlabelled copies' subject folders confirm.
BH498_30_200 and BH508_30_200_2 are filed under two rats but share all 203 trial times and
18,860 spike times; they are one recording sorted twice, and only BH498_30_200 is kept (the
other is a _2 repeat). 64 recordings remain.
units = tfus.unit_table(kept)
start = {s["identifier"]: s["session_start"][:19] for s in kept}
units["start"] = units.session.map(start)
# `later`: not the earliest of the rat's recordings of the same state, PRF and TBD (reading Rep3)
units["later"] = units.start != units.groupby(["rat", "state", "prf", "tbd"]).start.transform("min")
cond = units.groupby(["state", "prf", "tbd", "repeat", "rat"]).session.first().reset_index()
print("Recordings kept, by state, PRF, TBD and rat (x = one recording; R = a `_2` repeat):")
cond["mark"] = np.where(cond.repeat, "R", "x")
print(cond.pivot_table(index=["state", "prf", "tbd"], columns="rat", values="mark",
aggfunc=lambda v: "+".join(sorted(v))).fillna("").to_string())
iti = np.concatenate([np.diff(s["trial_times"]) for s in kept])
tl = np.concatenate([s["trial_len"] for s in kept])
print(f"\nInter-trial interval: {iti.min():.3f} to {iti.max():.3f} s, median {np.median(iti):.3f} "
f"(Methods: 2.5 s with 10% jitter); {int((iti > 2.75).sum())} intervals above 2.75 s")
print(f"stop_time - start_time: {tl.min():.3f} to {tl.max():.3f} s")
cover = [(s["identifier"], s["trial_times"][0], s["trial_times"][-1], *s["spike_range"]) for s in kept]
bad = [c for c in cover if c[1] - 0.5 < c[3] - 1 or c[2] + 1.5 > c[4] + 1]
print("Recordings whose spikes do not span the trials from -0.5 to +1.5 s (1 s slack):", bad or "none")
print("Spike times sorted within every unit:", all(s["spikes_sorted"] for s in kept))
rec = units.groupby(["rat", "state", "prf", "tbd", "session"]).agg(
start=("start", "first"), repeat=("repeat", "first"), later=("later", "first"), units=("unit", "size")).reset_index()
multi = rec.groupby(["rat", "state", "prf", "tbd"]).session.transform("size") > 1
print("\nConditions recorded more than once in a rat (repeat = `_2` suffix; later = not the rat's earliest recording of the condition):")
print(rec[multi].to_string(index=False))
print("`_2` recordings that are the earlier of their pair:", rec[multi & rec.repeat & ~rec.later].session.tolist())
print("\nAwake 1500 Hz recordings of the rats with DC recordings (TBD 40 and 400 us = 6% and 60%; 200 us = the PRF series):")
print(rec[(rec.state == "awake") & (rec.prf == 1500) & rec.rat.isin(rec[rec.tbd != 200].rat)].to_string(index=False))
Recordings kept, by state, PRF, TBD and rat (x = one recording; R = a `_2` repeat):
rat BH498 BH506 BH508 BH512 BH519 BH521 BH526 BH527
state prf tbd
anesthetized 30 200 x
300 200 x x x
1500 40 x
200 x x
400 x
3000 200 x x x
4500 200 x x x
awake 30 200 x x x x x x x
300 200 x x R+x x x x x
1500 40 x x x
200 R+x x R+x x x x x
400 x x x
3000 200 R+x R+x R+x x x x x
4500 200 R+x R+x R+x x x x x
Inter-trial interval: 2.250 to 25.336 s, median 2.499 (Methods: 2.5 s with 10% jitter); 3 intervals above 2.75 s
stop_time - start_time: 2.200 to 2.200 s
Recordings whose spikes do not span the trials from -0.5 to +1.5 s (1 s slack): [('BH508_1500_200', np.float64(0.27619328), np.float64(843.23606528), 4.43793408, 843.28824832)]
Spike times sorted within every unit: True
Conditions recorded more than once in a rat (repeat = `_2` suffix; later = not the rat's earliest recording of the condition):
rat state prf tbd session start repeat later units
BH498 awake 1500 200 BH498_1500_200 2024-03-06T12:42:27 False True 37
BH498 awake 1500 200 BH498_1500_200_2 2024-03-05T20:06:06 True False 64
BH498 awake 3000 200 BH498_3000_200 2023-09-06T18:09:13 False False 72
BH498 awake 3000 200 BH498_3000_200_2 2024-03-05T20:45:52 True True 64
BH498 awake 4500 200 BH498_4500_200 2024-03-06T23:13:17 False True 47
BH498 awake 4500 200 BH498_4500_200_2 2024-03-05T19:18:33 True False 64
BH506 awake 3000 200 BH506_3000_200 2023-10-26T17:27:06 False False 71
BH506 awake 3000 200 BH506_3000_200_2 2024-03-08T15:06:12 True True 39
BH506 awake 4500 200 BH506_4500_200 2023-10-26T15:51:06 False False 64
BH506 awake 4500 200 BH506_4500_200_2 2024-03-08T20:36:35 True True 53
BH508 awake 300 200 BH508_300_200 2023-08-27T23:54:25 False False 127
BH508 awake 300 200 BH508_300_200_2 2024-03-07T13:19:48 True True 38
BH508 awake 1500 200 BH508_1500_200 2023-10-22T18:52:10 False True 64
BH508 awake 1500 200 BH508_1500_200_2 2023-08-18T23:26:09 True False 68
BH508 awake 3000 200 BH508_3000_200 2023-08-27T14:08:43 False False 73
BH508 awake 3000 200 BH508_3000_200_2 2024-03-07T16:10:26 True True 33
BH508 awake 4500 200 BH508_4500_200 2023-08-27T17:22:57 False False 71
BH508 awake 4500 200 BH508_4500_200_2 2024-03-07T15:19:10 True True 58
`_2` recordings that are the earlier of their pair: ['BH498_1500_200_2', 'BH498_4500_200_2', 'BH508_1500_200_2']
Awake 1500 Hz recordings of the rats with DC recordings (TBD 40 and 400 us = 6% and 60%; 200 us = the PRF series):
rat state prf tbd session start repeat later units
BH506 awake 1500 40 BH506_1500_40 2024-03-09T15:43:17 False False 56
BH506 awake 1500 200 BH506_1500_200 2023-10-18T18:16:09 False False 71
BH506 awake 1500 400 BH506_1500_400 2024-03-09T16:25:19 False False 60
BH512 awake 1500 200 BH512_1500_200 2023-10-13T12:59:04 False False 78
BH521 awake 1500 40 BH521_1500_40 2024-03-09T14:56:24 False False 63
BH521 awake 1500 400 BH521_1500_400 2024-03-09T14:29:51 False False 83
BH526 awake 1500 40 BH526_1500_40 2024-01-30T19:52:06 False False 34
BH526 awake 1500 200 BH526_1500_200 2023-10-29T17:24:56 False False 68
BH526 awake 1500 400 BH526_1500_400 2024-01-30T18:43:18 False False 43
The dandiset holds the awake PRF series of seven rats (BH498, BH506, BH508, BH512, BH519, BH526,
BH527), ten _2 repeats of awake PRF conditions in BH498, BH506 and BH508, the awake 6% and 60%
DC conditions (TBD 40 and 400 us) of three rats (BH506, BH521, BH526), anesthetized PRF series
for BH498 (all five PRFs), BH512 (no 30 Hz) and BH521 (300, 3000 and 4500 Hz, unlabelled), and
anesthetized 6% and 60% DC recordings of BH512 (unlabelled). The paper describes ten rats in the
first batch, five more in a second batch with a sham condition, and a DC series with its own 30%
recordings at 10 months; none of the second batch, no sham and no separate 30% DC recording are
here. Intervals between trials fit the Methods' 2.5 s with 10% jitter, apart from one long gap.
The trial table's stop_time is always 2.2 s after start_time; the Methods give no 2.2 s
period, and start_time is taken as the trigger, which the next check supports.
The session start times show that the _2 suffix does not always mark the later recording.
Seven of the ten _2 recordings are four to seven months after the recording without the suffix
(BH498 3000 Hz; BH506 3000 and 4500 Hz; BH508 300, 3000 and 4500 Hz). Three are earlier:
BH498_1500_200_2 and BH498_4500_200_2 one day before the unsuffixed recordings (5 and 6 March
2024), and BH508_1500_200_2 two months before (August against October 2023). So which recording
of a pair is the "repeat" the description mentions is not settled by the names, and the readings
below include one based on the dates. The DC recordings (TBD 40 and 400 us) were made in January
2024 (BH526) and March 2024 (BH506, BH521), three to five months after the 1500 Hz, 200 us
recordings of the same rats (October 2023); BH521 has no awake 200 us recording. The other rats'
awake 1500 Hz, 200 us recordings date from August 2023 (BH508) to 6 March 2024 (BH498, 5 and 6
March, three to four days before the DC recordings of 9 March). The 30% group that the duty cycle
claims borrow therefore comes from other sessions, most of them months earlier.
lab_units = units[units.celltype.notna()]
print("Units in kept recordings:", len(units), "; labelled:", len(lab_units),
"; unlabelled (BH512 DC and BH521 anesthetized):", int(units.celltype.isna().sum()))
print("Units with no baseline spike (left out):", int((lab_units.base_hz == 0).sum()))
print("\nBaseline rate (Hz) by label, over units: median, IQR")
print(lab_units.groupby("celltype").base_hz.describe(percentiles=[.25, .5, .75])[["count", "25%", "50%", "75%", "max"]].round(2).to_string())
per_ses = lab_units.groupby(["session", "celltype"]).base_hz.median().unstack()
print(f"Recordings where the FSU median rate exceeds the RSU median: {(per_ses.FSU > per_ses.RSU).sum()} of {len(per_ses.dropna())}")
print("\nLabelled units per state and type (kept recordings):")
print(pd.crosstab([lab_units.state, lab_units.tbd], lab_units.celltype).to_string())
Units in kept recordings: 4202 ; labelled: 3882 ; unlabelled (BH512 DC and BH521 anesthetized): 320
Units with no baseline spike (left out): 0
Baseline rate (Hz) by label, over units: median, IQR
count 25% 50% 75% max
celltype
FSU 1959.0 5.10 8.60 12.49 106.94
RSU 1923.0 2.99 5.84 10.56 91.30
Recordings where the FSU median rate exceeds the RSU median: 46 of 59
Labelled units per state and type (kept recordings):
celltype FSU RSU
state tbd
anesthetized 200 285 291
awake 40 63 90
200 1536 1431
400 75 111
The labels split the units about evenly. FSU, the putative fast-spiking interneurons, are expected to fire faster; the comparison of baseline rates above says whether the labels follow that (it is a consistency check, not a test of the labels).
# Clock check: the population spike count of each recording in 0-67 ms after each trial start
# against the 0.5 s before it, at the true trial times and 1.25 s later (mid-interval).
rows = []
for s in kept:
c = s["counts"].sum(axis=0) # trials x windows, summed over units
i = {k: tfus.window_index(k) for k in ("base", "tl_w1", "shift_base", "shift_tl")}
true = c[:, i["tl_w1"]] / 0.067 - c[:, i["base"]] / 0.5
shift = c[:, i["shift_tl"]] / 0.067 - c[:, i["shift_base"]] / 0.5
rows.append({"session": s["identifier"], "state": s["state"], "prf": s["prf"], "tbd": s["tbd"],
"t_true": stats.ttest_1samp(true, 0).statistic, "t_shifted": stats.ttest_1samp(shift, 0).statistic})
clock = pd.DataFrame(rows)
print(clock.round(2).sort_values(["state", "tbd", "prf"]).to_string(index=False))
print(f"\nRecordings with t > 3 at the true times: {(clock.t_true > 3).sum()} of {len(clock)}; "
f"at shifted times: {(clock.t_shifted > 3).sum()}")
print(f"Median t: true {clock.t_true.median():.2f}, shifted {clock.t_shifted.median():.2f}")
fig, axes = plt.subplots(1, 3, figsize=(13, 3.2), sharey=True)
t = tfus.FINE_EDGES[:-1] + tfus.FINE / 2
for ax, (st, prf, tbd) in zip(axes, [("awake", 1500, 400), ("awake", 3000, 200), ("anesthetized", 1500, 200)]):
for code, name in [(1.0, "RSU"), (2.0, "FSU")]:
prof = []
for s in kept:
if (s["state"], s["prf"], s["tbd"]) != (st, prf, tbd):
continue
p = s["psth5"][s["label"] == code]
b = p[:, :100].mean(1)
prof.append(p[b > 0] / b[b > 0, None])
m = 100 * np.concatenate(prof).mean(0)
seg = {f"{a:+.3f} to {b:+.3f} s": m[(t > a) & (t < b)].mean() for a, b in ((-0.075, 0), (0, 0.075), (0.075, 0.15))}
print(f"{st} {prf} Hz {tbd} us {name}: mean % of baseline", {k: round(float(v), 1) for k, v in seg.items()},
f"; peak {m.max():.0f}% at {t[np.argmax(m)]:.4f} s")
ax.plot(t, m, label=f"{name} (n={len(np.concatenate(prof))}, range {m.min():.0f}-{m.max():.0f}%)")
ax.axvspan(0, 0.067, color="0.85")
ax.set_title(f"{st} {prf} Hz, TBD {tbd} us")
ax.set_xlabel("time from trial start (s)")
ax.legend(fontsize=7)
axes[0].set_ylabel("rate, % of baseline (5 ms bins)")
plt.tight_layout()
plt.show()
session state prf tbd t_true t_shifted
BH512_1500_40_anes anesthetized 1500 40 4.00 0.04
BH498_30_200_anes anesthetized 30 200 -0.62 0.59
BH498_300_200_anes anesthetized 300 200 0.96 -0.06
BH512_300_200_anes anesthetized 300 200 0.38 -0.21
BH521_300_200_anes anesthetized 300 200 1.60 -1.63
BH498_1500_200_anes anesthetized 1500 200 2.33 -0.57
BH512_1500_200_anes anesthetized 1500 200 9.96 1.04
BH498_3000_200_anes anesthetized 3000 200 -0.68 -0.02
BH512_3000_200_anes anesthetized 3000 200 6.57 0.75
BH521_3000_200_anes anesthetized 3000 200 0.78 1.05
BH498_4500_200_anes anesthetized 4500 200 -0.17 -1.79
BH512_4500_200_anes anesthetized 4500 200 0.70 -1.44
BH521_4500_200_anes anesthetized 4500 200 -0.93 1.14
BH512_1500_400_anes anesthetized 1500 400 30.03 1.83
BH506_1500_40 awake 1500 40 40.75 1.34
BH521_1500_40 awake 1500 40 3.97 2.49
BH526_1500_40 awake 1500 40 1.79 -0.05
BH498_30_200 awake 30 200 -0.85 -0.57
BH506_30_200 awake 30 200 0.74 -0.80
BH508_30_200 awake 30 200 -0.18 -2.87
BH512_30_200 awake 30 200 0.78 0.42
BH519_30_200 awake 30 200 -0.57 -0.02
BH526_30_200 awake 30 200 -1.55 -1.48
BH527_30_200 awake 30 200 -0.85 0.11
BH498_300_200 awake 300 200 0.68 1.05
BH506_300_200 awake 300 200 0.15 -0.31
BH508_300_200 awake 300 200 0.06 -0.03
BH508_300_200_2 awake 300 200 -0.76 -0.47
BH512_300_200 awake 300 200 -0.21 -1.12
BH519_300_200 awake 300 200 0.71 1.00
BH526_300_200 awake 300 200 -1.89 2.06
BH527_300_200 awake 300 200 1.76 -0.76
BH498_1500_200_2 awake 1500 200 -0.41 -1.14
BH498_1500_200 awake 1500 200 3.62 0.48
BH506_1500_200 awake 1500 200 2.73 0.74
BH508_1500_200_2 awake 1500 200 0.18 -0.27
BH508_1500_200 awake 1500 200 -0.85 -0.55
BH512_1500_200 awake 1500 200 0.29 -0.29
BH519_1500_200 awake 1500 200 0.08 0.49
BH526_1500_200 awake 1500 200 -2.58 0.00
BH527_1500_200 awake 1500 200 0.33 -0.88
BH498_3000_200 awake 3000 200 1.94 1.78
BH498_3000_200_2 awake 3000 200 1.47 -0.31
BH506_3000_200 awake 3000 200 12.96 1.15
BH506_3000_200_2 awake 3000 200 -1.38 -0.98
BH508_3000_200 awake 3000 200 2.64 -0.19
BH508_3000_200_2 awake 3000 200 -0.19 -1.63
BH512_3000_200 awake 3000 200 -0.56 -1.64
BH519_3000_200 awake 3000 200 0.24 -1.40
BH526_3000_200 awake 3000 200 0.26 -1.96
BH527_3000_200 awake 3000 200 0.70 -1.21
BH498_4500_200_2 awake 4500 200 0.99 0.12
BH498_4500_200 awake 4500 200 1.57 -1.05
BH506_4500_200 awake 4500 200 0.79 -0.07
BH506_4500_200_2 awake 4500 200 0.78 -2.06
BH508_4500_200 awake 4500 200 0.26 1.94
BH508_4500_200_2 awake 4500 200 -0.17 -0.70
BH512_4500_200 awake 4500 200 -1.09 -1.17
BH519_4500_200 awake 4500 200 -1.56 0.84
BH526_4500_200 awake 4500 200 -1.03 1.31
BH527_4500_200 awake 4500 200 -0.05 1.91
BH506_1500_400 awake 1500 400 34.02 0.40
BH521_1500_400 awake 1500 400 4.87 0.59
BH526_1500_400 awake 1500 400 1.45 -1.64
Recordings with t > 3 at the true times: 10 of 64; at shifted times: 0
Median t: true 0.35, shifted -0.06
awake 1500 Hz 400 us RSU: mean % of baseline {'-0.075 to +0.000 s': 103.6, '+0.000 to +0.075 s': 146.3, '+0.075 to +0.150 s': 112.4} ; peak 159% at 0.0275 s
awake 1500 Hz 400 us FSU: mean % of baseline {'-0.075 to +0.000 s': 101.1, '+0.000 to +0.075 s': 196.9, '+0.075 to +0.150 s': 115.1} ; peak 228% at 0.0575 s
awake 3000 Hz 200 us RSU: mean % of baseline {'-0.075 to +0.000 s': 101.2, '+0.000 to +0.075 s': 107.6, '+0.075 to +0.150 s': 100.3} ; peak 117% at 0.0225 s
awake 3000 Hz 200 us FSU: mean % of baseline {'-0.075 to +0.000 s': 100.6, '+0.000 to +0.075 s': 106.5, '+0.075 to +0.150 s': 101.1} ; peak 110% at 0.0575 s
anesthetized 1500 Hz 200 us RSU: mean % of baseline {'-0.075 to +0.000 s': 99.8, '+0.000 to +0.075 s': 106.2, '+0.075 to +0.150 s': 100.2} ; peak 113% at 0.0325 s
anesthetized 1500 Hz 200 us FSU: mean % of baseline {'-0.075 to +0.000 s': 100.2, '+0.000 to +0.075 s': 103.7, '+0.075 to +0.150 s': 102.2} ; peak 111% at 0.5875 s
Ten of the 64 recordings show a population rise in the 67 ms after start_time (t > 3) and none
shows one at the shifted times. In the conditions plotted, the mean rate is near baseline in the
75 ms before start_time (99.8 to 103.6%), rises in the 75 ms after it (103.7 to 196.9%), the
span of the 67 ms ultrasound, and falls back in the next 75 ms (100.2 to 115.1%). So
start_time is the trigger, and spikes and trials share a clock in those recordings, which
include awake and anesthetized ones and four of the eight rats (BH498, BH506, BH512, BH521). In
the other 54 recordings the check cannot tell a shared clock from no response.
This also limits every time-locked claim. Among the awake 200 us PRF-series recordings, only BH498_1500_200 (t = 3.62) and BH506_3000_200 (t = 12.96) show a population rise in 0 to 67 ms, and among the anesthetized 200 us ones only BH512's 1500 and 3000 Hz (t = 9.96 and 6.57). The 40 and 400 us recordings respond most clearly (6 of the 8 at t > 3). Group differences in the time-locked window at 200 us therefore rest on small changes in most recordings.
Readings and selections used in every claim¶
The claims file's header fixes the method; the dataset decides which of its readings can run.
- Cell type. K0: the dataset carries labels, so they are used. K1 and K2 (k-means on the initial and after-hyperpolarization phase durations) cannot run: the files hold no waveforms. That the labels are the paper's k-means output is an assumption; the baseline rates above are consistent with it and do not prove it.
- Trials. The files hold 500 trials, the 505 delivered less the first five (T1, Methods 2.4); T2 (all 505) cannot run. That the 500 are trials 6 to 505 cannot be checked: the dandiset's description says only that "each recording has 500 trials". Five recordings hold 203 to 339 trials and are used as they are.
- Exclusion. E1 = E2: the dataset marks no excluded unit or interval.
- Windows. W1 (exact window) and W2 (overlapping PSTH bins) are both run.
- Groups. G1 only: there are no sham recordings, so G2 cannot run.
- Repeats (not named by the claims file; the dandiset adds them). Ten awake recordings repeat
a condition in the same rat. The paper mentions no repeats, and Methods 2.3 fits one recording per
parameter set per rat without stating it ("Recordings for a single rat for each set of
parameters were completed on the same day" says when a rat's recordings were made), and the dandiset's description calls them repeats
taken "a few months after the original ones". Three readings are run: Rep1 counts every
recording once, repeats included; Rep2 leaves out the
_2recordings; Rep3 keeps the earliest recording of each rat and condition, following the description's timing where the suffix and the dates disagree (three pairs, see above). Rep2 and Rep3 both give one recording per rat and condition. Where Rep1 and those two give different outcomes, one recording per rat and condition is adopted as an assumption (Methods 2.3 and the dandiset's description); where Rep2 and Rep3 differ from each other, neither is preferred, as the BH498 pairs are a day apart. - Anesthetized data. Only two rats' anesthetized recordings carry labels (BH498 at all five PRFs, BH512 without 30 Hz); BH521's three are unlabelled and cannot enter a cell-type test. Rat and PRF are therefore confounded under anesthesia: the 30 Hz group is BH498 alone.
- Sham. The dandiset holds no sham recording. Parts of a criterion that compare with sham are reported as not tested. A claim with an untested part is not judged reproduced: its best outcome is partially reproduced. A claim is settled as not reproduced only when the parts tested rule out the partial outcome whatever the untested parts would show; a claim that could still be partially reproduced if an untested part held is blocked as data missing (fig4 below). The claims file's header provides this for an anesthetized sham; the same is done for the awake sham.
- Readings of the claims file that disagree. Where the readings the claims file names (W1/W2, S1/S2, D1/D2, L1/L2, A1/A2) give different outcomes, the claim is blocked as methods_underspecified, as the header says, with every result shown.
- The 30% duty cycle. The paper's DC series had its own 30% recordings, which are not in the dandiset. The substitutes (C1, C2, defined before the duty cycle claims) and C0, without a 30% group, are shown. A part of a criterion that needs the 30% group is treated as data missing, and a claim whose outcome depends on that part is blocked as data missing.
L = units[units.celltype.notna() & (units.base_hz > 0)].copy()
PRFS = [30, 300, 1500, 3000, 4500]
REPS = ("Rep1", "Rep2", "Rep3")
READINGS = [(w, r) for w in ("W1", "W2") for r in REPS]
def col(kind, state, w):
if kind == "tl":
return "tl_w1" if w == "W1" else "tl_w2"
return tfus.delayed_window(state, w)
def select(state, ct, rep, tbd=200, prf=None):
d = L[(L.state == state) & (L.celltype == ct) & (L.tbd == tbd)]
if prf is not None:
d = d[d.prf == prf]
if rep == "Rep2":
d = d[~d.repeat]
if rep == "Rep3":
d = d[~d.later]
return d
def prf_test(state, ct, kind, w, rep):
d = select(state, ct, rep)
c = col(kind, state, w)
return tfus.group_test({p: d[d.prf == p][c].values for p in PRFS})
def sig(post, a, b):
return bool(tfus.pair(post, a, b).p_bonf < 0.05)
def lower(post, a, b):
"""a significantly below b."""
r = tfus.pair(post, a, b)
return bool(r.p_bonf < 0.05 and r.mean_a < r.mean_b)
def show(o, post, label):
print(f"{label}: Kruskal-Wallis H = {o['H']:.1f}, p = {o['p']:.2g}; n = {o['n']}")
print(" means:", {k: round(v, 1) for k, v in o["mean"].items()},
" medians:", {k: round(v, 1) for k, v in o["median"].items()})
print(post[["a", "b", "mean_a", "mean_b", "p", "p_bonf", "sig"]].round(4).to_string(index=False))
def by_rat(state, ct, kind, w, rep, tbd=200):
d = select(state, ct, rep, tbd=tbd)
c = col(kind, state, w)
return d.pivot_table(index="rat", columns="prf" if tbd == 200 else "dc", values=c, aggfunc="mean").round(1)
print("Labelled units with a baseline, by state and type:", L.groupby(["state", "celltype"]).size().to_dict())
Labelled units with a baseline, by state and type: {('anesthetized', 'FSU'): 285, ('anesthetized', 'RSU'): 291, ('awake', 'FSU'): 1674, ('awake', 'RSU'): 1632}
fig3-cell-type-sorting¶
Units are sorted into fast spiking units (FSU) and regular spiking units (RSU) by k-means clustering of the initial phase and after-hyperpolarization phase durations of their action potentials; in one animal's recordings this gave 28 FSU and 32 RSU, 60 in total.
cnt = units[units.celltype.notna()].groupby(["rat", "session", "state"]).celltype.value_counts().unstack(fill_value=0)
cnt["total"] = cnt.sum(axis=1)
ok = lambda c: (abs(c.FSU - 28) <= 2.8) & (abs(c.RSU - 32) <= 3.2) & (abs(c.total - 60) <= 6)
print("R1, one recording: recordings with all three counts in tolerance (FSU 28 +- 2.8, RSU 32 +- 3.2, total 60 +- 6)")
print(cnt[ok(cnt)].to_string())
print(f"({len(cnt)} labelled recordings; totals range {cnt.total.min()} to {cnt.total.max()}, FSU share "
f"{(cnt.FSU / cnt.total).min():.2f} to {(cnt.FSU / cnt.total).max():.2f})")
lu = units[units.celltype.notna()].copy()
lu["series"] = np.where(lu.tbd == 200, lu.state + " PRF", lu.state + " DC")
cnt2 = lu.groupby(["rat", "series"]).celltype.value_counts().unstack(fill_value=0)
cnt2["total"] = cnt2.sum(axis=1)
print("\nR2, all of one rat's recordings in one series:")
print(cnt2.to_string())
print("R2 matches:", cnt2[ok(cnt2)].index.tolist() or "none")
rep_ses, later_ses = set(units[units.repeat].session), set(units[units.later].session)
for rp, drop in (("Rep2", rep_ses), ("Rep3", later_ses)):
c = cnt[~cnt.index.get_level_values("session").isin(drop)]
print(f"R1 under {rp}: {len(c)} recordings, matches:", c[ok(c)].index.get_level_values("session").tolist())
band = (cnt.total - 60).abs() <= 6
fsu_share = cnt.FSU.sum() / cnt.total.sum()
p_match = [stats.binom.pmf(k, n, fsu_share) for n in cnt.total
for k in range(n + 1) if abs(k - 28) <= 2.8 and abs(n - k - 32) <= 3.2]
print(f"\nBase rate: {int(band.sum())} of {len(cnt)} recordings have a total within 60 +- 6; of these {int(ok(cnt).sum())} "
f"also match on FSU and RSU. If each recording's FSU count were binomial with the pooled FSU share "
f"({fsu_share:.3f}), {sum(p_match):.2f} recordings would be expected to match all three counts.")
print("\nThe clusters' separation on the initial phase and after-hyperpolarization durations cannot be "
"computed: the files hold no waveforms (units columns:", sorted(set(inv.units_cols)), ")")
R1, one recording: recordings with all three counts in tolerance (FSU 28 +- 2.8, RSU 32 +- 3.2, total 60 +- 6)
celltype FSU RSU total
rat session state
BH508 BH508_4500_200_2 awake 28 30 58
BH527 BH527_30_200 awake 28 34 62
(59 labelled recordings; totals range 31 to 127, FSU share 0.16 to 0.88)
R2, all of one rat's recordings in one series:
celltype FSU RSU total
rat series
BH498 anesthetized PRF 163 157 320
awake PRF 202 225 427
BH506 awake DC 35 81 116
awake PRF 234 201 435
BH508 awake PRF 266 350 616
BH512 anesthetized PRF 122 134 256
awake PRF 269 110 379
BH519 awake PRF 202 194 396
BH521 awake DC 61 85 146
BH526 awake DC 42 35 77
awake PRF 193 182 375
BH527 awake PRF 170 169 339
R2 matches: none
R1 under Rep2: 50 recordings, matches: ['BH527_30_200']
R1 under Rep3: 50 recordings, matches: ['BH527_30_200']
Base rate: 22 of 59 recordings have a total within 60 +- 6; of these 2 also match on FSU and RSU. If each recording's FSU count were binomial with the pooled FSU share (0.505), 3.84 recordings would be expected to match all three counts.
The clusters' separation on the initial phase and after-hyperpolarization durations cannot be computed: the files hold no waveforms (units columns: ['spike_times', 'spike_times celltype_label'] )
Result: blocked (data missing). Both the reproduced and the partially reproduced outcome of the criterion require the clusters to separate on the initial phase and after-hyperpolarization durations, and the files hold no waveforms, so neither can be reached. With the dataset's labels (K0), two single recordings give all three counts within tolerance (R1): BH527_30_200 with 28 FSU, 34 RSU and 62 units, and BH508_4500_200_2 with 28 FSU, 30 RSU and 58 units; under Rep2 and Rep3 only BH527_30_200 does, as BH508_4500_200_2 is a later repeat. No rat's series matches (R2; totals 77 to 616). The match is no more than chance: 22 of 59 recordings have a total within 60 +- 6, and if each recording's FSU count were binomial with the pooled FSU share (0.505), 3.84 recordings would be expected to match all three counts, against 2 observed. The labels follow the expected rate difference in 46 of 59 recordings (FSU median baseline rate above the RSU median; 8.60 against 5.84 Hz over units). The other claims use the labels directly, as the claims file's header says for K0, and do not need the waveform features.
fig4-rsu-time-locked-both-states (blocked)¶
RSUs show time-locked responses to tFUS, an increase in firing rate during the 67 ms of stimulation, in both awake and anesthetized recordings.
The criterion compares the 1500 Hz, 200 us RSU responses with the sham sessions of each state. The dandiset holds no sham recording, so the claim is blocked. The responses the criterion would have compared are shown, without a test.
for st in ("awake", "anesthetized"):
for w in ("W1", "W2"):
for rep in ("Rep1", "Rep2"):
d = select(st, "RSU", rep, prf=1500)[col("tl", st, w)]
print(f"{st:12s} {w} {rep}: RSU at 1500 Hz, n = {len(d)}, mean {d.mean():.1f}%, median {d.median():.1f}%, "
f"units above 100%: {(d > 100).sum()}")
awake W1 Rep1: RSU at 1500 Hz, n = 303, mean 98.5%, median 100.7%, units above 100%: 157 awake W1 Rep2: RSU at 1500 Hz, n = 227, mean 98.1%, median 101.9%, units above 100%: 120 awake W2 Rep1: RSU at 1500 Hz, n = 303, mean 98.7%, median 100.9%, units above 100%: 163 awake W2 Rep2: RSU at 1500 Hz, n = 227, mean 98.1%, median 101.2%, units above 100%: 124 anesthetized W1 Rep1: RSU at 1500 Hz, n = 78, mean 107.0%, median 106.4%, units above 100%: 65 anesthetized W1 Rep2: RSU at 1500 Hz, n = 78, mean 107.0%, median 106.4%, units above 100%: 65 anesthetized W2 Rep1: RSU at 1500 Hz, n = 78, mean 106.4%, median 105.5%, units above 100%: 64 anesthetized W2 Rep2: RSU at 1500 Hz, n = 78, mean 106.4%, median 105.5%, units above 100%: 64
Result: blocked (data missing). No sham recording exists. The awake half already fails one of the criterion's conditions: the awake 1500 Hz RSU mean is 98.1 to 98.7% of baseline under every reading, not above 100%. So the claim could not be reproduced whatever a sham would show; at best it would be partially reproduced, on the anesthetized half (106.4 to 107.0%, two rats).
fig4hi-fsu-delayed-latency¶
FSU delayed responses to tFUS at 1500 Hz PRF begin later in awake rats (mean onset 0.51 s, spread 0.14 s) than under anesthesia (0.24 s, spread 0.08 s), and the onset distributions differ significantly by a Wilcoxon rank-sum test.
Sessions: FSU of the 1500 Hz, 200 us recordings. Q1 and Q2 select the same recordings here: the awake 30% DC condition has no recording of its own (it is the 1500 Hz, 200 us recording), and batch 2 is not in the dandiset. Every combination of S1/S2, D1/D2, L1/L2 and Rep1/Rep2 is run.
lat_rows = []
for s in kept:
if s["prf"] != 1500 or s["tbd"] != 200:
continue
for u in range(s["n_units"]):
if s["label"][u] != 2.0:
continue
r = {"session": s["identifier"], "rat": s["rat"], "state": s["state"], "repeat": s["repeat"]}
for sem in ("S1", "S2"):
for first in ("D1", "D2"):
for timing in ("L1", "L2"):
r[f"{sem}{first}{timing}"] = tfus.latency(s["psth"][u], s["counts"][u, :, 0], sem, first, timing)
lat_rows.append(r)
lat = pd.DataFrame(lat_rows)
combos = [c for c in lat.columns if c.startswith("S")]
res = []
for c in combos:
for rep in REPS:
d = {"Rep1": lat, "Rep2": lat[~lat.repeat], "Rep3": lat[~lat.session.isin(later_ses)]}[rep]
aw, an = d[d.state == "awake"][c].dropna(), d[d.state == "anesthetized"][c].dropna()
n_aw, n_an = (d.state == "awake").sum(), (d.state == "anesthetized").sum()
p = stats.ranksums(aw, an).pvalue
res.append({"reading": c + " " + rep, "awake n": f"{len(aw)}/{n_aw}", "anes n": f"{len(an)}/{n_an}",
"awake mean": aw.mean(), "awake SD": aw.std(), "anes mean": an.mean(), "anes SD": an.std(),
"p": p, "anes in tol": abs(an.mean() - 0.24) <= 0.08, "awake in tol": abs(aw.mean() - 0.51) <= 0.14,
"awake later, p<0.001": bool(p < 0.001 and aw.median() > an.median())})
lat_res = pd.DataFrame(res)
print("FSU at 1500 Hz, 200 us; n = units with a significant delayed response / FSU")
print(lat_res.round(4).to_string(index=False))
main = lat[lat.state == "awake"]
print("\nBy rat, S1D1L1 Rep1: units with a delayed response and their mean latency (s)")
print(lat.groupby(["state", "rat"]).S1D1L1.agg(["size", "count", "mean"]).round(3).to_string())
for st in ("awake", "anesthetized"):
d = lat[lat.state == st]
v = d.S1D1L1.dropna()
print(f"{st}: latencies range {v.min():.3f} to {v.max():.3f} s; first allowed bin (D1) holds {np.mean(v < 0.08):.2f} of them")
FSU at 1500 Hz, 200 us; n = units with a significant delayed response / FSU
reading awake n anes n awake mean awake SD anes mean anes SD p anes in tol awake in tol awake later, p<0.001
S1D1L1 Rep1 262/306 39/50 0.4466 0.3844 0.5055 0.3724 0.2332 False True False
S1D1L1 Rep2 217/250 39/50 0.4998 0.3883 0.5055 0.3724 0.8491 False True False
S1D1L1 Rep3 231/272 39/50 0.4561 0.3956 0.5055 0.3724 0.2503 False True False
S1D1L2 Rep1 262/306 39/50 0.4644 0.3844 0.5234 0.3724 0.2332 False True False
S1D1L2 Rep2 217/250 39/50 0.5177 0.3883 0.5234 0.3724 0.8491 False True False
S1D1L2 Rep3 231/272 39/50 0.4739 0.3956 0.5234 0.3724 0.2503 False True False
S1D2L1 Rep1 265/306 39/50 0.4031 0.3812 0.4066 0.4055 0.7328 False True False
S1D2L1 Rep2 219/250 39/50 0.4529 0.3876 0.4066 0.4055 0.2721 False True False
S1D2L1 Rep3 234/272 39/50 0.4174 0.3932 0.4066 0.4055 0.6613 False True False
S1D2L2 Rep1 265/306 39/50 0.4210 0.3812 0.4245 0.4055 0.7328 False True False
S1D2L2 Rep2 219/250 39/50 0.4707 0.3876 0.4245 0.4055 0.2721 False True False
S1D2L2 Rep3 234/272 39/50 0.4353 0.3932 0.4245 0.4055 0.6613 False True False
S2D1L1 Rep1 193/306 34/50 0.5067 0.3875 0.4842 0.3668 0.8307 False True False
S2D1L1 Rep2 167/250 34/50 0.5616 0.3846 0.4842 0.3668 0.2645 False True False
S2D1L1 Rep3 169/272 34/50 0.5036 0.3879 0.4842 0.3668 0.9210 False True False
S2D1L2 Rep1 193/306 34/50 0.5245 0.3875 0.5021 0.3668 0.8307 False True False
S2D1L2 Rep2 167/250 34/50 0.5794 0.3846 0.5021 0.3668 0.2645 False True False
S2D1L2 Rep3 169/272 34/50 0.5214 0.3879 0.5021 0.3668 0.9210 False True False
S2D2L1 Rep1 195/306 34/50 0.4745 0.3907 0.3792 0.3926 0.0739 False True False
S2D2L1 Rep2 169/250 34/50 0.5245 0.3926 0.3792 0.3926 0.0180 False True False
S2D2L1 Rep3 171/272 34/50 0.4731 0.3889 0.3792 0.3926 0.0845 False True False
S2D2L2 Rep1 195/306 34/50 0.4924 0.3907 0.3971 0.3926 0.0739 False True False
S2D2L2 Rep2 169/250 34/50 0.5424 0.3926 0.3971 0.3926 0.0180 False True False
S2D2L2 Rep3 171/272 34/50 0.4909 0.3889 0.3971 0.3926 0.0845 False True False
By rat, S1D1L1 Rep1: units with a delayed response and their mean latency (s)
size count mean
state rat
anesthetized BH498 33 28 0.519
BH512 17 11 0.471
awake BH498 44 43 0.252
BH506 34 31 0.341
BH508 46 33 0.282
BH512 59 44 0.677
BH519 35 35 0.409
BH526 33 25 0.946
BH527 55 51 0.363
awake: latencies range 0.071 to 1.393 s; first allowed bin (D1) holds 0.23 of them
anesthetized: latencies range 0.071 to 1.286 s; first allowed bin (D1) holds 0.13 of them
Result: not reproduced. Under none of the 24 readings (S1/S2, D1/D2, L1/L2, Rep1/Rep2/Rep3) do awake FSU onsets come later than anesthetized ones at p < 0.001: rank-sum p ranges from 0.018 to 0.92, and in 8 readings the anesthetized mean is the later one. The awake means (0.40 to 0.58 s) are within 0.51 +- 0.14 s, but the anesthetized means (0.38 to 0.52 s) are all outside 0.24 +- 0.08 s. Most FSU qualify as having a delayed response (193 to 265 of 306 awake under Rep1, 262 of 306, or 86%, under S1 D1; 34 or 39 of 50 anesthetized), and latencies spread over 0.07 to 1.39 s awake and 0.07 to 1.29 s anesthetized, so the one-SEM rule of the Methods, as the claims file reads it, flags most units and the onsets look largely noise-driven. Rat means range from 0.25 s (BH498) to 0.95 s (BH526) awake. The anesthetized FSU come from two rats.
fig4h-awake-latency-wider¶
Awake FSU delayed-response onset latencies have a wider spread than anesthetized ones and appear to have two peaks.
print(lat_res[["reading", "awake SD", "anes SD"]].assign(wider=lambda d: d["awake SD"] > d["anes SD"]).round(3).to_string(index=False))
fig, ax = plt.subplots(figsize=(6, 3))
bins = tfus.EDGES[tfus.N_BASE_BINS:]
for st in ("anesthetized", "awake"):
v = lat[lat.state == st].S1D1L1.dropna()
h, _ = np.histogram(v, bins)
peaks = [round(float(bins[i]), 3) for i in range(len(h)) if h[i] > 0 and (i == 0 or h[i] > h[i - 1]) and (i == len(h) - 1 or h[i] >= h[i + 1])]
print(f"{st}: histogram (35.7 ms bins from 0 s) counts {h.tolist()}; local maxima at {peaks} s")
ax.hist(v, bins, alpha=0.6, label=f"{st} (n={len(v)}, {v.min():.2f}-{v.max():.2f} s)")
ax.set_xlabel("onset latency (s), S1 D1 L1 Rep1")
ax.set_ylabel("FSU")
ax.legend(fontsize=8)
plt.show()
reading awake SD anes SD wider S1D1L1 Rep1 0.384 0.372 True S1D1L1 Rep2 0.388 0.372 True S1D1L1 Rep3 0.396 0.372 True S1D1L2 Rep1 0.384 0.372 True S1D1L2 Rep2 0.388 0.372 True S1D1L2 Rep3 0.396 0.372 True S1D2L1 Rep1 0.381 0.406 False S1D2L1 Rep2 0.388 0.406 False S1D2L1 Rep3 0.393 0.406 False S1D2L2 Rep1 0.381 0.406 False S1D2L2 Rep2 0.388 0.406 False S1D2L2 Rep3 0.393 0.406 False S2D1L1 Rep1 0.387 0.367 True S2D1L1 Rep2 0.385 0.367 True S2D1L1 Rep3 0.388 0.367 True S2D1L2 Rep1 0.387 0.367 True S2D1L2 Rep2 0.385 0.367 True S2D1L2 Rep3 0.388 0.367 True S2D2L1 Rep1 0.391 0.393 False S2D2L1 Rep2 0.393 0.393 True S2D2L1 Rep3 0.389 0.393 False S2D2L2 Rep1 0.391 0.393 False S2D2L2 Rep2 0.393 0.393 True S2D2L2 Rep3 0.389 0.393 False anesthetized: histogram (35.7 ms bins from 0 s) counts [0, 0, 5, 1, 1, 4, 1, 1, 3, 2, 3, 0, 1, 0, 0, 0, 2, 1, 0, 1, 2, 1, 0, 2, 0, 1, 0, 0, 1, 0, 3, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0]; local maxima at [0.071, 0.179, 0.286, 0.357, 0.429, 0.571, 0.714, 0.821, 0.893, 1.0, 1.071, 1.179, 1.286] s awake: histogram (35.7 ms bins from 0 s) counts [0, 0, 61, 12, 9, 9, 12, 11, 11, 10, 8, 8, 10, 10, 10, 2, 7, 1, 1, 6, 5, 1, 4, 3, 3, 1, 3, 5, 4, 8, 5, 3, 3, 1, 0, 3, 3, 2, 3, 4, 0, 0]; local maxima at [0.071, 0.214, 0.429, 0.571, 0.679, 0.786, 0.964, 1.036, 1.25, 1.393] s
Result: blocked (methods underspecified). The claims file's readings give different results: the awake SD exceeds the anesthetized SD under every D1 reading (0.384 to 0.396 s against 0.367 to 0.372 s), but not under D2, where it is smaller in 10 of 12 readings (0.381 to 0.393 s against 0.393 to 0.406 s) and larger by less than 0.001 s in the other two (0.393 against 0.393 s). The header blocks a claim whose readings disagree. Under every reading the SDs are within 0.03 s of each other, with no contrast like the paper's 0.14 against 0.08 s. Under S1 D1 L1 Rep1 the awake histogram has its largest count in the first allowed bin (61 of 262 at 0.071 s) and several smaller local maxima; the anesthetized histogram (39 units) is flat with many single-count maxima.
fig5a-rsu-time-locked-prf-awake¶
In awake rats PRF is a significant factor for RSU time-locked responses; 30 Hz and 300 Hz PRF give lower responses than every higher PRF (so 30 Hz is also below 300 Hz), 1500, 3000 and 4500 Hz do not differ, and sham differs from every PRF but 30 Hz.
Parts (5) and (6) compare with sham and cannot be tested, so the claim can at best be partially reproduced.
rows = []
for w, rep in READINGS:
o, post = prf_test("awake", "RSU", "tl", w, rep)
if (w, rep) == ("W1", "Rep1"):
show(o, post, "awake RSU 0-67 ms, W1 Rep1")
rows.append({"reading": f"{w} {rep}", "(1) KW": o["p"] < 0.05,
"(2) 30 below each": all(lower(post, 30, b) for b in (300, 1500, 3000, 4500)),
"(3) 300 below each": all(lower(post, 300, b) for b in (1500, 3000, 4500)),
"(4) no pair of 1500/3000/4500": not any(sig(post, a, b) for a, b in ((1500, 3000), (1500, 4500), (3000, 4500))),
"highest mean": max(o["mean"], key=o["mean"].get)})
print(pd.DataFrame(rows).to_string(index=False))
print("\nMean response by rat (W1 Rep1):")
print(by_rat("awake", "RSU", "tl", "W1", "Rep1").to_string())
awake RSU 0-67 ms, W1 Rep1: Kruskal-Wallis H = 78.7, p = 3.3e-16; n = {30: 255, 300: 269, 1500: 303, 3000: 334, 4500: 270}
means: {30: 96.8, 300: 97.6, 1500: 98.5, 3000: 108.2, 4500: 99.2} medians: {30: 96.5, 300: 99.5, 1500: 100.7, 3000: 104.5, 4500: 99.9}
a b mean_a mean_b p p_bonf sig
30 300 96.8249 97.6406 0.0096 0.0955 False
30 1500 96.8249 98.4572 0.0124 0.1241 False
30 3000 96.8249 108.2052 0.0000 0.0000 True
30 4500 96.8249 99.1909 0.0028 0.0277 True
300 1500 97.6406 98.4572 0.6262 1.0000 False
300 3000 97.6406 108.2052 0.0000 0.0000 True
300 4500 97.6406 99.1909 0.8846 1.0000 False
1500 3000 98.4572 108.2052 0.0000 0.0000 True
1500 4500 98.4572 99.1909 0.7678 1.0000 False
3000 4500 108.2052 99.1909 0.0000 0.0000 True
reading (1) KW (2) 30 below each (3) 300 below each (4) no pair of 1500/3000/4500 highest mean
W1 Rep1 True False False False 3000
W1 Rep2 True False False False 3000
W1 Rep3 True False False False 3000
W2 Rep1 True True False False 3000
W2 Rep2 True False False False 3000
W2 Rep3 True True False False 3000
Mean response by rat (W1 Rep1):
prf 30 300 1500 3000 4500
rat
BH498 87.4 104.2 103.9 113.5 103.8
BH506 106.4 100.1 108.9 130.1 103.8
BH508 98.7 95.6 96.1 100.8 97.3
BH512 105.6 95.8 110.7 97.9 90.7
BH519 98.4 108.6 102.5 105.9 87.9
BH526 94.3 64.0 71.1 100.6 96.0
BH527 89.6 110.1 100.1 102.0 99.8
Result: not reproduced. PRF is a significant factor (Kruskal-Wallis p = 3.3e-16, W1 Rep1), but the post hoc pattern differs. The highest mean is at 3000 Hz (108.2%; 1500 Hz 98.5%), and 3000 Hz differs from 1500 and 4500 Hz (corrected p < 1e-4), so part (4) fails under every reading. 300 Hz is not below 1500 or 4500 Hz (part 3 fails under every reading), and 30 Hz is below every higher PRF only under W2 Rep1 and W2 Rep3 (part 2). Parts (5) and (6) need the sham and are not tested; even if both held, (1) with at most three other parts falls short of the four the partial outcome needs. Means by rat vary widely (BH526 at 64 to 101%, BH506 at 100 to 130%).
fig5c-rsu-time-locked-prf-anesth¶
Under anesthesia RSU time-locked responses depend on PRF as in the awake model, the only difference being that 30 Hz and 300 Hz do not differ significantly.
Two labelled rats: BH498 (all five PRFs) and BH512 (no 30 Hz). The repeats reading does not apply (no anesthetized repeats remain). Parts (5) and (6) compare with sham and cannot be tested.
rows = []
for w in ("W1", "W2"):
o, post = prf_test("anesthetized", "RSU", "tl", w, "Rep1")
show(o, post, f"anesthetized RSU 0-67 ms, {w}")
rows.append({"reading": w, "(1) KW": o["p"] < 0.05,
"(2) 30 below 1500/3000/4500": all(lower(post, 30, b) for b in (1500, 3000, 4500)),
"(3) 300 below each": all(lower(post, 300, b) for b in (1500, 3000, 4500)),
"(4) no pair of 1500/3000/4500": not any(sig(post, a, b) for a, b in ((1500, 3000), (1500, 4500), (3000, 4500))),
"(7) 30 = 300": not sig(post, 30, 300)})
print(pd.DataFrame(rows).to_string(index=False))
print("\nMean response by rat (W1):")
print(by_rat("anesthetized", "RSU", "tl", "W1", "Rep1").to_string())
anesthetized RSU 0-67 ms, W1: Kruskal-Wallis H = 38.9, p = 7.3e-08; n = {30: 32, 300: 59, 1500: 78, 3000: 60, 4500: 62}
means: {30: 100.5, 300: 100.8, 1500: 107.0, 3000: 103.1, 4500: 100.7} medians: {30: 100.8, 300: 99.7, 1500: 106.4, 3000: 103.2, 4500: 101.3}
a b mean_a mean_b p p_bonf sig
30 300 100.4574 100.7609 0.9470 1.0000 False
30 1500 100.4574 106.9661 0.0000 0.0001 True
30 3000 100.4574 103.1425 0.0530 0.5303 False
30 4500 100.4574 100.6983 0.8544 1.0000 False
300 1500 100.7609 106.9661 0.0000 0.0000 True
300 3000 100.7609 103.1425 0.0297 0.2972 False
300 4500 100.7609 100.6983 0.8115 1.0000 False
1500 3000 106.9661 103.1425 0.0024 0.0238 True
1500 4500 106.9661 100.6983 0.0000 0.0000 True
3000 4500 103.1425 100.6983 0.0609 0.6088 False
anesthetized RSU 0-67 ms, W2: Kruskal-Wallis H = 33.9, p = 7.7e-07; n = {30: 32, 300: 59, 1500: 78, 3000: 60, 4500: 62}
means: {30: 100.4, 300: 100.7, 1500: 106.4, 3000: 103.0, 4500: 100.6} medians: {30: 100.3, 300: 99.6, 1500: 105.5, 3000: 103.6, 4500: 100.7}
a b mean_a mean_b p p_bonf sig
30 300 100.4218 100.6581 0.8811 1.0000 False
30 1500 100.4218 106.3717 0.0001 0.0006 True
30 3000 100.4218 103.0139 0.0464 0.4636 False
30 4500 100.4218 100.5692 0.8294 1.0000 False
300 1500 100.6581 106.3717 0.0000 0.0000 True
300 3000 100.6581 103.0139 0.0211 0.2107 False
300 4500 100.6581 100.5692 0.6783 1.0000 False
1500 3000 106.3717 103.0139 0.0180 0.1796 False
1500 4500 106.3717 100.5692 0.0000 0.0002 True
3000 4500 103.0139 100.5692 0.0574 0.5744 False
reading (1) KW (2) 30 below 1500/3000/4500 (3) 300 below each (4) no pair of 1500/3000/4500 (7) 30 = 300
W1 True False False False True
W2 True False False False True
Mean response by rat (W1):
prf 30 300 1500 3000 4500
rat
BH498 100.5 101.0 101.7 100.6 100.0
BH512 NaN 100.5 110.4 106.0 101.4
Result: not reproduced. Under anesthesia (two labelled rats) the Kruskal-Wallis is significant (p = 7.3e-08 W1, 7.7e-07 W2) and 30 Hz and 300 Hz do not differ (part 7 holds), but parts (2), (3) and (4) do not: 1500 Hz is the peak (107.0%), 30 Hz is not below 3000 or 4500 Hz, 300 Hz is not below 3000 or 4500 Hz, and 1500 Hz differs from 4500 Hz (and from 3000 Hz under W1). Even if the two sham parts held, (1) with three other parts falls short of the five the partial outcome needs. With two rats, PRF and rat are confounded: the 30 Hz group is BH498 alone, and the 1500 Hz peak comes from BH512 (110.4% against 101.7% in BH498; BH498 ranges 100.0 to 101.7% over all five PRFs).
rsu-time-locked-anesthesia-factor¶
Anesthesia is a significant factor in the RSU time-locked response by two-way ANOVA, and it lowers the average magnitude of RSU responses.
def anova(ct, kind, w, rep, model, drop=()):
d = pd.concat([select(st, ct, rep) for st in ("awake", "anesthetized")])
d = d[~d.session.isin(drop)]
d = d.assign(y=[r[col(kind, r.state, w)] for _, r in d.iterrows()])
form = "y ~ C(state) * C(prf)" if model == "A1" else "y ~ C(state) + C(prf)"
a = sm.stats.anova_lm(ols(form, d).fit(), typ=2)
return {"reading": f"{model} {w} {rep}", "n": len(d), "F state": a.loc["C(state)", "F"], "p state": a.loc["C(state)", "PR(>F)"],
"awake mean": d[d.state == "awake"].y.mean(), "anes mean": d[d.state == "anesthetized"].y.mean()}
def anova_table(ct, kind):
return pd.DataFrame([anova(ct, kind, w, rep, m) for m in ("A1", "A2") for w, rep in READINGS])
t_rsu = anova_table("RSU", "tl")
print(t_rsu.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
print("\nPer-PRF means (W1 Rep1):")
print(pd.concat({st: select(st, "RSU", "Rep1").groupby("prf").tl_w1.mean() for st in ("awake", "anesthetized")}, axis=1).round(1).to_string())
print("\nVariant, the two rats with labelled anesthetized data only (BH498, BH512), W1 Rep1, A1:")
d = pd.concat([select(st, "RSU", "Rep1") for st in ("awake", "anesthetized")])
d = d[d.rat.isin(["BH498", "BH512"])].assign(y=lambda x: x.tl_w1)
a = sm.stats.anova_lm(ols("y ~ C(state) * C(prf)", d).fit(), typ=2)
print(f" F = {a.loc['C(state)', 'F']:.2f}, p = {a.loc['C(state)', 'PR(>F)']:.3g}; awake mean {d[d.state == 'awake'].y.mean():.1f}, "
f"anesthetized {d[d.state == 'anesthetized'].y.mean():.1f}")
reading n F state p state awake mean anes mean
A1 W1 Rep1 1722 6.473 0.01104 100.4 102.9
A1 W1 Rep2 1464 6.142 0.01332 100.3 102.9
A1 W1 Rep3 1485 6.977 0.008342 100.3 102.9
A1 W2 Rep1 1722 5.871 0.0155 100.3 102.6
A1 W2 Rep2 1464 5.667 0.01742 100.1 102.6
A1 W2 Rep3 1485 6.273 0.01236 100.1 102.6
A2 W1 Rep1 1722 6.404 0.01148 100.4 102.9
A2 W1 Rep2 1464 6.063 0.01392 100.3 102.9
A2 W1 Rep3 1485 6.877 0.008822 100.3 102.9
A2 W2 Rep1 1722 5.816 0.01598 100.3 102.6
A2 W2 Rep2 1464 5.603 0.01806 100.1 102.6
A2 W2 Rep3 1485 6.196 0.01291 100.1 102.6
Per-PRF means (W1 Rep1):
awake anesthetized
prf
30 96.8 100.5
300 97.6 100.8
1500 98.5 107.0
3000 108.2 103.1
4500 99.2 100.7
Variant, the two rats with labelled anesthetized data only (BH498, BH512), W1 Rep1, A1:
F = 0.64, p = 0.425; awake mean 103.4, anesthetized 102.9
Result: not reproduced. Anesthesia is a factor at p = 0.0083 to 0.018 over the twelve readings (below 0.01 only under W1 Rep3), but in the reverse direction: the anesthetized mean (102.6 to 102.9%) is above the awake mean (100.1 to 100.4%) under every reading. The anesthetized units come from two rats and the awake units from seven, so the state effect is partly a rat effect: with only the two rats that have labelled anesthetized data, the effect disappears (p = 0.425; awake 103.4%, anesthetized 102.9%).
fig5bd-fsu-time-locked-no-prf¶
FSU time-locked responses show no significant differences across PRF levels, in either the awake or the anesthetized model.
A claim of absence, judged on the Kruskal-Wallis across the five PRFs (G1) at 0.05 and the post hoc tests at Bonferroni-corrected 0.05.
rows = []
for st in ("awake", "anesthetized"):
for w, rep in READINGS:
if st == "anesthetized" and rep != "Rep1":
continue
o, post = prf_test(st, "FSU", "tl", w, rep)
if w == "W1" and rep == "Rep1":
show(o, post, f"{st} FSU 0-67 ms, W1 Rep1")
med = o["median"]
rows.append({"state": st, "reading": f"{w} {rep}", "KW p": o["p"], "post hoc pairs significant": int(post.sig.sum()),
"median range": f"{min(med.values()):.1f}-{max(med.values()):.1f}",
"mean range": f"{min(o['mean'].values()):.1f}-{max(o['mean'].values()):.1f}"})
print(pd.DataFrame(rows).round(6).to_string(index=False))
print("\nMean response by rat, awake (W1 Rep1):")
print(by_rat("awake", "FSU", "tl", "W1", "Rep1").to_string())
awake FSU 0-67 ms, W1 Rep1: Kruskal-Wallis H = 61.0, p = 1.8e-12; n = {30: 232, 300: 286, 1500: 306, 3000: 319, 4500: 393}
means: {30: 99.5, 300: 96.7, 1500: 101.8, 3000: 106.8, 4500: 99.5} medians: {30: 99.7, 300: 98.7, 1500: 100.8, 3000: 103.4, 4500: 99.9}
a b mean_a mean_b p p_bonf sig
30 300 99.5112 96.7398 0.0647 0.6465 False
30 1500 99.5112 101.7893 0.0768 0.7679 False
30 3000 99.5112 106.7696 0.0000 0.0000 True
30 4500 99.5112 99.4510 0.8754 1.0000 False
300 1500 96.7398 101.7893 0.0003 0.0035 True
300 3000 96.7398 106.7696 0.0000 0.0000 True
300 4500 96.7398 99.4510 0.0247 0.2467 False
1500 3000 101.7893 106.7696 0.0005 0.0049 True
1500 4500 101.7893 99.4510 0.0729 0.7288 False
3000 4500 106.7696 99.4510 0.0000 0.0000 True
anesthetized FSU 0-67 ms, W1 Rep1: Kruskal-Wallis H = 19.8, p = 0.00055; n = {30: 32, 300: 69, 1500: 50, 3000: 68, 4500: 66}
means: {30: 100.3, 300: 100.3, 1500: 104.4, 3000: 102.9, 4500: 99.3} medians: {30: 100.0, 300: 100.0, 1500: 103.9, 3000: 103.3, 4500: 99.8}
a b mean_a mean_b p p_bonf sig
30 300 100.2636 100.3250 0.9709 1.0000 False
30 1500 100.2636 104.4308 0.0058 0.0584 False
30 3000 100.2636 102.9424 0.0400 0.3996 False
30 4500 100.2636 99.2805 0.7561 1.0000 False
300 1500 100.3250 104.4308 0.0012 0.0121 True
300 3000 100.3250 102.9424 0.0165 0.1649 False
300 4500 100.3250 99.2805 0.7050 1.0000 False
1500 3000 104.4308 102.9424 0.5240 1.0000 False
1500 4500 104.4308 99.2805 0.0004 0.0043 True
3000 4500 102.9424 99.2805 0.0057 0.0572 False
state reading KW p post hoc pairs significant median range mean range
awake W1 Rep1 0.000000 5 98.7-103.4 96.7-106.8
awake W1 Rep2 0.000000 6 98.7-104.1 96.7-108.7
awake W1 Rep3 0.000000 4 98.7-104.1 96.7-108.7
awake W2 Rep1 0.000000 5 99.3-103.2 97.3-106.4
awake W2 Rep2 0.000000 6 98.8-103.8 97.3-108.3
awake W2 Rep3 0.000000 4 99.2-103.8 97.3-108.3
anesthetized W1 Rep1 0.000553 2 99.8-103.9 99.3-104.4
anesthetized W2 Rep1 0.000778 3 99.8-104.1 99.4-104.2
Mean response by rat, awake (W1 Rep1):
prf 30 300 1500 3000 4500
rat
BH498 98.7 100.4 105.8 105.0 107.8
BH506 101.2 100.6 115.3 124.2 101.8
BH508 99.2 93.8 100.3 103.6 101.0
BH512 102.2 99.6 100.2 98.2 91.0
BH519 95.0 101.7 99.2 101.8 96.6
BH526 97.7 91.6 91.4 100.6 98.7
BH527 99.5 102.8 101.1 105.1 99.6
Result: not reproduced. By the paper's test FSU time-locked responses do differ across PRF: Kruskal-Wallis p = 1.8e-12 awake and 5.5e-04 anesthetized (W1 Rep1), with 4 to 6 significant post hoc pairs awake and 2 or 3 anesthetized under every reading. The effects are small: group medians span 98.7 to 103.4% awake (W1 Rep1) and 99.8 to 103.9% anesthetized; the awake 3000 Hz mean is 106.8%. The anesthetized groups come from two rats (30 Hz from BH498 alone).
fig5e-rsu-delayed-prf-awake¶
Awake RSU delayed responses (0.5 to 1 s) depend on PRF, are highest at 1500 Hz and 3000 Hz with lower responses at both lower and higher PRFs, and every PRF is significantly higher than sham.
Part (3) compares with sham and cannot be tested. The 0.5 to 1 s window is on bin edges, so W1 and W2 are the same.
rows = []
for rep in REPS:
o, post = prf_test("awake", "RSU", "del", "W1", rep)
show(o, post, f"awake RSU 0.5-1 s, {rep}")
m = o["mean"]
top2 = sorted(m, key=m.get)[-2:]
rows.append({"reading": rep, "(1) KW": o["p"] < 0.05,
"(2) R1: 1500 and 3000 above 30, 300, 4500": all(lower(post, a, b) for b in (1500, 3000) for a in (30, 300, 4500)),
"(2) R2: top two means are 1500 and 3000": sorted(top2) == [1500, 3000],
"(2) R2: both above 30 and 4500": all(lower(post, a, b) for b in (1500, 3000) for a in (30, 4500))})
print(pd.DataFrame(rows).to_string(index=False))
print("\nMean response by rat (Rep1):")
print(by_rat("awake", "RSU", "del", "W1", "Rep1").to_string())
awake RSU 0.5-1 s, Rep1: Kruskal-Wallis H = 96.3, p = 6e-20; n = {30: 255, 300: 269, 1500: 303, 3000: 334, 4500: 270}
means: {30: 106.6, 300: 101.2, 1500: 106.7, 3000: 102.3, 4500: 100.0} medians: {30: 105.5, 300: 99.9, 1500: 103.7, 3000: 100.0, 4500: 99.5}
a b mean_a mean_b p p_bonf sig
30 300 106.6032 101.2027 0.0000 0.0000 True
30 1500 106.6032 106.7424 0.0286 0.2863 False
30 3000 106.6032 102.3071 0.0000 0.0000 True
30 4500 106.6032 99.9992 0.0000 0.0000 True
300 1500 101.2027 106.7424 0.0029 0.0292 True
300 3000 101.2027 102.3071 0.6732 1.0000 False
300 4500 101.2027 99.9992 0.1044 1.0000 False
1500 3000 106.7424 102.3071 0.0005 0.0053 True
1500 4500 106.7424 99.9992 0.0000 0.0001 True
3000 4500 102.3071 99.9992 0.1202 1.0000 False
awake RSU 0.5-1 s, Rep2: Kruskal-Wallis H = 73.3, p = 4.5e-15; n = {30: 255, 300: 246, 1500: 227, 3000: 268, 4500: 177}
means: {30: 106.6, 300: 102.8, 1500: 108.8, 3000: 103.3, 4500: 100.3} medians: {30: 105.5, 300: 101.0, 1500: 105.8, 3000: 100.1, 4500: 100.0}
a b mean_a mean_b p p_bonf sig
30 300 106.6032 102.8418 0.0000 0.0000 True
30 1500 106.6032 108.8265 0.9397 1.0000 False
30 3000 106.6032 103.3249 0.0000 0.0000 True
30 4500 106.6032 100.3441 0.0000 0.0000 True
300 1500 102.8418 108.8265 0.0087 0.0871 False
300 3000 102.8418 103.3249 0.1428 1.0000 False
300 4500 102.8418 100.3441 0.0133 0.1332 False
1500 3000 108.8265 103.3249 0.0011 0.0107 True
1500 4500 108.8265 100.3441 0.0000 0.0002 True
3000 4500 103.3249 100.3441 0.2069 1.0000 False
awake RSU 0.5-1 s, Rep3: Kruskal-Wallis H = 73.1, p = 5.1e-15; n = {30: 255, 300: 246, 1500: 236, 3000: 268, 4500: 189}
means: {30: 106.6, 300: 102.8, 1500: 107.3, 3000: 103.3, 4500: 100.5} medians: {30: 105.5, 300: 101.0, 1500: 103.9, 3000: 100.1, 4500: 99.9}
a b mean_a mean_b p p_bonf sig
30 300 106.6032 102.8418 0.0000 0.0000 True
30 1500 106.6032 107.3363 0.0569 0.5689 False
30 3000 106.6032 103.3249 0.0000 0.0000 True
30 4500 106.6032 100.5445 0.0000 0.0000 True
300 1500 102.8418 107.3363 0.0974 0.9738 False
300 3000 102.8418 103.3249 0.1428 1.0000 False
300 4500 102.8418 100.5445 0.0191 0.1906 False
1500 3000 107.3363 103.3249 0.0067 0.0668 False
1500 4500 107.3363 100.5445 0.0002 0.0018 True
3000 4500 103.3249 100.5445 0.2408 1.0000 False
reading (1) KW (2) R1: 1500 and 3000 above 30, 300, 4500 (2) R2: top two means are 1500 and 3000 (2) R2: both above 30 and 4500
Rep1 True False False False
Rep2 True False False False
Rep3 True False False False
Mean response by rat (Rep1):
prf 30 300 1500 3000 4500
rat
BH498 112.0 99.4 101.6 100.7 97.8
BH506 97.7 100.0 116.6 99.0 102.2
BH508 108.0 90.8 103.0 97.1 98.5
BH512 109.6 100.4 85.2 91.5 121.0
BH519 99.4 114.6 128.4 127.0 94.2
BH526 103.8 106.4 92.5 101.0 100.5
BH527 114.6 101.8 113.0 96.1 86.8
Result: blocked (data missing). PRF is a significant factor (p = 6e-20 Rep1, 4.5e-15 Rep2, 5.1e-15 Rep3), but the two highest means are at 1500 Hz and 30 Hz (106.7 and 106.6% Rep1; 108.8 and 106.6% Rep2; 107.3 and 106.6% Rep3), not 1500 and 3000 Hz: 3000 Hz (102.3 to 103.3%) is below 30 Hz (corrected p < 1e-4) and not above 300 Hz. Part (2) fails under both of its readings and every repeat reading.
Part (3), every PRF above sham, needs the sham and cannot be tested. The criterion gives partially reproduced when (1) and (3) hold without (2), so if a sham held the outcome would be partially reproduced, and if it did not, not reproduced. The data present do not decide between the two, so by the rule stated under Readings (as for fig4-rsu-time-locked-both-states) the claim is blocked as data missing. The second verifier pointed this out, and the third agreed. Adding a sham group (G2) would not rescue part (2), as it enlarges the Bonferroni family.
fig5f-fsu-delayed-prf-awake¶
Awake FSU delayed responses depend on PRF, are strongest at 4500 Hz, differ between 30 Hz and the three highest PRFs and between 300 Hz and the two highest, and at 30 Hz and 300 Hz are not significantly higher than sham.
Part (5) compares with sham and cannot be tested.
rows = []
for rep in REPS:
o, post = prf_test("awake", "FSU", "del", "W1", rep)
show(o, post, f"awake FSU 0.5-1 s, {rep}")
rows.append({"reading": rep, "(1) KW": o["p"] < 0.05, "(2) 4500 highest": max(o["mean"], key=o["mean"].get) == 4500,
"(3) 30 differs from 1500/3000/4500": all(sig(post, 30, b) for b in (1500, 3000, 4500)),
"(4) 300 differs from 3000/4500": all(sig(post, 300, b) for b in (3000, 4500))})
print(pd.DataFrame(rows).to_string(index=False))
print("\nMean response by rat (Rep1):")
print(by_rat("awake", "FSU", "del", "W1", "Rep1").to_string())
awake FSU 0.5-1 s, Rep1: Kruskal-Wallis H = 59.0, p = 4.6e-12; n = {30: 232, 300: 286, 1500: 306, 3000: 319, 4500: 393}
means: {30: 103.5, 300: 98.7, 1500: 101.8, 3000: 100.6, 4500: 100.3} medians: {30: 102.0, 300: 99.4, 1500: 100.8, 3000: 100.9, 4500: 99.7}
a b mean_a mean_b p p_bonf sig
30 300 103.4912 98.7202 0.0000 0.0000 True
30 1500 103.4912 101.8436 0.0009 0.0086 True
30 3000 103.4912 100.6370 0.0004 0.0036 True
30 4500 103.4912 100.2888 0.0000 0.0000 True
300 1500 98.7202 101.8436 0.0003 0.0027 True
300 3000 98.7202 100.6370 0.0028 0.0276 True
300 4500 98.7202 100.2888 0.3284 1.0000 False
1500 3000 101.8436 100.6370 0.6806 1.0000 False
1500 4500 101.8436 100.2888 0.0004 0.0042 True
3000 4500 100.6370 100.2888 0.0030 0.0301 True
awake FSU 0.5-1 s, Rep2: Kruskal-Wallis H = 53.8, p = 5.7e-11; n = {30: 232, 300: 271, 1500: 250, 3000: 249, 4500: 311}
means: {30: 103.5, 300: 99.1, 1500: 101.7, 3000: 102.1, 4500: 100.4} medians: {30: 102.0, 300: 99.7, 1500: 100.7, 3000: 100.9, 4500: 99.6}
a b mean_a mean_b p p_bonf sig
30 300 103.4912 99.1447 0.0000 0.0000 True
30 1500 103.4912 101.6684 0.0008 0.0081 True
30 3000 103.4912 102.1280 0.0046 0.0463 True
30 4500 103.4912 100.4356 0.0000 0.0000 True
300 1500 99.1447 101.6684 0.0074 0.0737 False
300 3000 99.1447 102.1280 0.0009 0.0093 True
300 4500 99.1447 100.4356 0.7428 1.0000 False
1500 3000 101.6684 102.1280 0.4950 1.0000 False
1500 4500 101.6684 100.4356 0.0027 0.0274 True
3000 4500 102.1280 100.4356 0.0001 0.0011 True
awake FSU 0.5-1 s, Rep3: Kruskal-Wallis H = 56.5, p = 1.6e-11; n = {30: 232, 300: 271, 1500: 272, 3000: 249, 4500: 316}
means: {30: 103.5, 300: 99.1, 1500: 101.5, 3000: 102.1, 4500: 100.4} medians: {30: 102.0, 300: 99.7, 1500: 100.7, 3000: 100.9, 4500: 99.5}
a b mean_a mean_b p p_bonf sig
30 300 103.4912 99.1447 0.0000 0.0000 True
30 1500 103.4912 101.5019 0.0003 0.0025 True
30 3000 103.4912 102.1280 0.0046 0.0463 True
30 4500 103.4912 100.3760 0.0000 0.0000 True
300 1500 99.1447 101.5019 0.0066 0.0661 False
300 3000 99.1447 102.1280 0.0009 0.0093 True
300 4500 99.1447 100.3760 0.8180 1.0000 False
1500 3000 101.5019 102.1280 0.3665 1.0000 False
1500 4500 101.5019 100.3760 0.0019 0.0187 True
3000 4500 102.1280 100.3760 0.0001 0.0005 True
reading (1) KW (2) 4500 highest (3) 30 differs from 1500/3000/4500 (4) 300 differs from 3000/4500
Rep1 True False True False
Rep2 True False True False
Rep3 True False True False
Mean response by rat (Rep1):
prf 30 300 1500 3000 4500
rat
BH498 126.4 100.2 104.9 103.8 100.2
BH506 100.0 100.7 101.5 100.6 101.0
BH508 104.1 91.2 101.9 89.1 97.4
BH512 105.4 101.0 99.0 98.4 108.9
BH519 98.5 109.8 104.4 110.1 98.0
BH526 100.4 100.6 96.8 101.4 100.4
BH527 111.3 103.5 104.0 100.7 95.6
Result: not reproduced. PRF is a significant factor (p = 4.6e-12 Rep1) and 30 Hz differs from 1500, 3000 and 4500 Hz (part 3), but the strongest mean is at 30 Hz (103.5%), not 4500 Hz (100.3 to 100.4%), and 300 Hz does not differ from 4500 Hz (corrected p = 1.0), under every repeat reading. Part (3) is tested without direction, as the criterion says "differs", but 30 Hz is above the higher PRFs, the opposite of the paper's pattern. Part (5) needs the sham; even if it held, (1) with two other parts falls short of the three the partial outcome needs.
fig5gh-delayed-prf-anesth¶
Under anesthesia PRF is a significant factor for the delayed responses (0.2 to 0.7 s) of both RSUs and FSUs, which broadly prefer higher PRFs.
rows = []
for ct in ("RSU", "FSU"):
for w in ("W1", "W2"):
o, post = prf_test("anesthetized", ct, "del", w, "Rep1")
d = select("anesthetized", ct, "Rep1")[[col("del", "anesthetized", w), "prf"]]
hi = d[d.prf >= 1500].iloc[:, 0].mean()
lo = d[d.prf <= 300].iloc[:, 0].mean()
rows.append({"type": ct, "reading": w, "KW H": o["H"], "KW p": o["p"], "n": o["n"],
"mean 1500-4500": hi, "mean 30-300": lo, "higher PRFs higher": hi > lo})
print(pd.DataFrame(rows).round(4).to_string(index=False))
print("\nMeans by rat and PRF, RSU W1:")
print(by_rat("anesthetized", "RSU", "del", "W1", "Rep1").to_string())
print("FSU W1:")
print(by_rat("anesthetized", "FSU", "del", "W1", "Rep1").to_string())
type reading KW H KW p n mean 1500-4500 mean 30-300 higher PRFs higher
RSU W1 6.5559 0.1613 {30: 32, 300: 59, 1500: 78, 3000: 60, 4500: 62} 99.9680 99.1495 True
RSU W2 7.2563 0.1229 {30: 32, 300: 59, 1500: 78, 3000: 60, 4500: 62} 99.9851 99.1774 True
FSU W1 6.6081 0.1581 {30: 32, 300: 69, 1500: 50, 3000: 68, 4500: 66} 99.8696 99.3486 True
FSU W2 6.2279 0.1828 {30: 32, 300: 69, 1500: 50, 3000: 68, 4500: 66} 99.8810 99.3662 True
Means by rat and PRF, RSU W1:
prf 30 300 1500 3000 4500
rat
BH498 99.3 99.6 99.5 99.7 99.7
BH512 NaN 98.5 100.3 100.4 100.0
FSU W1:
prf 30 300 1500 3000 4500
rat
BH498 99.0 99.6 99.3 99.8 99.4
BH512 NaN 99.5 99.4 100.9 100.1
Result: not reproduced. Under anesthesia PRF is not a significant factor for delayed responses: Kruskal-Wallis p = 0.16 (W1) and 0.12 (W2) for RSU, 0.16 and 0.18 for FSU. The pooled 1500 to 4500 Hz means are above the 30 and 300 Hz means by less than one percentage point (RSU 99.97 against 99.15%, FSU 99.87 against 99.35%, W1). The data are two rats, and the 30 Hz group is BH498 alone; every rat and PRF mean is within 98.5 to 100.9%.
delayed-anesthesia-factor¶
Anesthesia is a significant factor in the delayed response for both RSUs and FSUs by two-way ANOVA, and it lowers the average response magnitude.
tabs = {ct: anova_table(ct, "del") for ct in ("RSU", "FSU")}
for ct, t in tabs.items():
print(ct)
print(t.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
chk = pd.DataFrame({"reading": tabs["RSU"].reading,
"RSU p<0.001": tabs["RSU"]["p state"] < 0.001, "FSU p<0.01": tabs["FSU"]["p state"] < 0.01,
"both p<0.05": (tabs["RSU"]["p state"] < 0.05) & (tabs["FSU"]["p state"] < 0.05),
"anes lower, both types": (tabs["RSU"]["anes mean"] < tabs["RSU"]["awake mean"]) & (tabs["FSU"]["anes mean"] < tabs["FSU"]["awake mean"])})
print(chk.to_string(index=False))
print("\nFSU per-rat means, awake recordings with and without the repeats (0.5-1 s):")
print(select("awake", "FSU", "Rep1").groupby(["rat", "repeat"]).del_awake.agg(["size", "mean"]).round(1).to_string())
print("\nVariant without BH498_30_200, the recording also filed as BH508_30_200_2 (if it is BH508's, it is a later repeat):")
var = pd.DataFrame([{"type": ct, **anova(ct, "del", w, rp, m, drop=("BH498_30_200",))}
for ct in ("RSU", "FSU") for m in ("A1", "A2") for w, rp in READINGS])
print(var.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
RSU
reading n F state p state awake mean anes mean
A1 W1 Rep1 1722 21.63 3.556e-06 103.4 99.71
A1 W1 Rep2 1464 34.93 4.243e-09 104.6 99.71
A1 W1 Rep3 1485 32.45 1.471e-08 104.3 99.71
A1 W2 Rep1 1722 21.38 4.049e-06 103.4 99.73
A1 W2 Rep2 1464 34.63 4.951e-09 104.6 99.73
A1 W2 Rep3 1485 32.15 1.717e-08 104.3 99.73
A2 W1 Rep1 1722 21.52 3.769e-06 103.4 99.71
A2 W1 Rep2 1464 34.68 4.818e-09 104.6 99.71
A2 W1 Rep3 1485 32.3 1.585e-08 104.3 99.71
A2 W2 Rep1 1722 21.27 4.284e-06 103.4 99.73
A2 W2 Rep2 1464 34.38 5.61e-09 104.6 99.73
A2 W2 Rep3 1485 32 1.847e-08 104.3 99.73
FSU
reading n F state p state awake mean anes mean
A1 W1 Rep1 1821 4.316 0.03789 100.9 99.68
A1 W1 Rep2 1598 10.38 0.001302 101.3 99.68
A1 W1 Rep3 1625 10.54 0.001194 101.2 99.68
A1 W2 Rep1 1821 4.202 0.04051 100.9 99.7
A1 W2 Rep2 1598 10.19 0.001438 101.3 99.7
A1 W2 Rep3 1625 10.35 0.001324 101.2 99.7
A2 W1 Rep1 1821 4.296 0.03834 100.9 99.68
A2 W1 Rep2 1598 10.33 0.001335 101.3 99.68
A2 W1 Rep3 1625 10.49 0.001228 101.2 99.68
A2 W2 Rep1 1821 4.182 0.04101 100.9 99.7
A2 W2 Rep2 1598 10.14 0.001477 101.3 99.7
A2 W2 Rep3 1625 10.29 0.001363 101.2 99.7
reading RSU p<0.001 FSU p<0.01 both p<0.05 anes lower, both types
A1 W1 Rep1 True False True True
A1 W1 Rep2 True True True True
A1 W1 Rep3 True True True True
A1 W2 Rep1 True False True True
A1 W2 Rep2 True True True True
A1 W2 Rep3 True True True True
A2 W1 Rep1 True False True True
A2 W1 Rep2 True True True True
A2 W1 Rep3 True True True True
A2 W2 Rep1 True False True True
A2 W2 Rep2 True True True True
A2 W2 Rep3 True True True True
FSU per-rat means, awake recordings with and without the repeats (0.5-1 s):
size mean
rat repeat
BH498 False 110 102.0
True 92 104.4
BH506 False 188 101.1
True 46 99.1
BH508 False 181 97.6
True 85 91.8
BH512 False 269 102.4
BH519 False 202 102.8
BH526 False 193 100.0
BH527 False 170 102.6
Variant without BH498_30_200, the recording also filed as BH508_30_200_2 (if it is BH508's, it is a later repeat):
type reading n F state p state awake mean anes mean RSU A1 W1 Rep1 1696 20.89 5.218e-06 103.2 99.71 RSU A1 W1 Rep2 1438 34.1 6.468e-09 104.4 99.71 RSU A1 W1 Rep3 1459 31.66 2.203e-08 104.1 99.71 RSU A1 W2 Rep1 1696 20.64 5.929e-06 103.2 99.73 RSU A1 W2 Rep2 1438 33.8 7.536e-09 104.4 99.73 RSU A1 W2 Rep3 1459 31.35 2.568e-08 104.1 99.73 RSU A2 W1 Rep1 1696 20.79 5.497e-06 103.2 99.71 RSU A2 W1 Rep2 1438 33.86 7.305e-09 104.4 99.71 RSU A2 W1 Rep3 1459 31.52 2.359e-08 104.1 99.71 RSU A2 W2 Rep1 1696 20.54 6.238e-06 103.2 99.73 RSU A2 W2 Rep2 1438 33.56 8.494e-09 104.4 99.73 RSU A2 W2 Rep3 1459 31.22 2.745e-08 104.1 99.73 FSU A1 W1 Rep1 1816 3.915 0.04802 100.8 99.68 FSU A1 W1 Rep2 1593 9.878 0.001703 101.2 99.68 FSU A1 W1 Rep3 1620 10.03 0.001568 101.1 99.68 FSU A1 W2 Rep1 1816 3.804 0.05127 100.8 99.7 FSU A1 W2 Rep2 1593 9.695 0.001881 101.2 99.7 FSU A1 W2 Rep3 1620 9.84 0.001739 101.1 99.7 FSU A2 W1 Rep1 1816 3.899 0.04846 100.8 99.68 FSU A2 W1 Rep2 1593 9.842 0.001737 101.2 99.68 FSU A2 W1 Rep3 1620 9.991 0.001602 101.1 99.68 FSU A2 W2 Rep1 1816 3.789 0.05176 100.8 99.7 FSU A2 W2 Rep2 1593 9.656 0.001921 101.2 99.7 FSU A2 W2 Rep3 1620 9.797 0.00178 101.1 99.7
Result: reproduced, resting on an assumption. With one recording per rat and condition, the
reading adopted, anesthesia is a significant factor for both types under all four model and
window readings, whether the _2 recordings are left out (Rep2) or the earliest recording of each
pair is kept (Rep3): RSU p = 4.2e-09 to 1.8e-08 and FSU p = 0.0012 to 0.0015, with the
anesthetized mean lower in both (RSU 99.7 against 104.3 to 104.6%, FSU 99.7 against 101.2 to
101.3%). With the repeats (Rep1) the RSU part holds (p < 1e-5) but the FSU effect is p = 0.038 to
0.041, which misses the 0.01 bound: partially reproduced under that reading. Leaving out the
recording filed under both BH498 and BH508 does not change the adopted outcome (FSU p = 0.0016 to
0.0019 under Rep2 and Rep3). Under Rep1 it raises the FSU p to 0.048 to 0.052, above 0.05 under
W2; the RSU part alone still holds, which the criterion counts as partially reproduced ("when
one cell type holds"). So the reproduced outcome rests on the one-recording assumption: with every
recording the claim is partial at best. The anesthetized units come from two rats, so the state effect may also be a rat
effect.
fig6ad-time-adjacent-batch (blocked)¶
In the five rats recorded awake and anesthetized in the same week, time-locked responses to the low PRF and 1500 Hz differ significantly in both states, with a reduced response under anesthesia.
The ten rats in the dandiset are the first batch's (each has awake PRF recordings), and no file is from a sham or a two-PRF protocol; the five rats of the second batch are not in it. Blocked as data missing; no section is run.
fig6ef-rsu-ispta-rise-then-fall¶
Awake RSU mean responses, time-locked and delayed, rise with ISPTA up to around 60 mW/cm2 and then decline.
def ispta_table(ct):
rows = []
for w, rep in READINGS:
for kind in ("tl", "del"):
if kind == "del" and w == "W2":
continue
d = select("awake", ct, rep)
c = col(kind, "awake", w)
g = d.groupby("prf")[c]
m, se = g.mean(), g.sem()
rows.append({"window": kind, "reading": f"{w} {rep}",
**{f"{p} Hz ({tfus.PRF_ISPTA[p]})": f"{m[p]:.1f} +- {se[p]:.1f}" for p in PRFS},
"_m": m})
return rows
rows = ispta_table("RSU")
out = []
for r in rows:
m = r.pop("_m")
peak = m.idxmax()
rise = all(m[PRFS[i]] <= m[PRFS[i + 1]] for i in range(PRFS.index(peak))) if peak in (1500, 3000) else False
out.append({**r, "peak": peak, "peak at 1500/3000": peak in (1500, 3000), "rises 30 to peak": rise,
"4500 below peak": m[4500] < m[peak]})
print("Mean +- SEM (%) by PRF (ISPTA mW/cm2):")
print(pd.DataFrame(out).to_string(index=False))
Mean +- SEM (%) by PRF (ISPTA mW/cm2):
window reading 30 Hz (0.6) 300 Hz (6.03) 1500 Hz (30.13) 3000 Hz (60.26) 4500 Hz (90.39) peak peak at 1500/3000 rises 30 to peak 4500 below peak
tl W1 Rep1 96.8 +- 0.7 97.6 +- 1.0 98.5 +- 1.0 108.2 +- 1.1 99.2 +- 0.7 3000 True True True
del W1 Rep1 106.6 +- 0.5 101.2 +- 0.6 106.7 +- 1.0 102.3 +- 0.8 100.0 +- 0.6 1500 True False True
tl W1 Rep2 96.8 +- 0.7 98.2 +- 1.1 98.1 +- 1.3 108.7 +- 1.3 98.1 +- 0.9 3000 True False True
del W1 Rep2 106.6 +- 0.5 102.8 +- 0.6 108.8 +- 1.3 103.3 +- 0.9 100.3 +- 0.9 1500 True False True
tl W1 Rep3 96.8 +- 0.7 98.2 +- 1.1 98.3 +- 1.1 108.7 +- 1.3 98.0 +- 0.7 3000 True True True
del W1 Rep3 106.6 +- 0.5 102.8 +- 0.6 107.3 +- 1.2 103.3 +- 0.9 100.5 +- 0.8 1500 True False True
tl W2 Rep1 96.4 +- 0.7 97.4 +- 1.0 98.7 +- 1.0 107.7 +- 1.1 99.5 +- 0.7 3000 True True True
tl W2 Rep2 96.4 +- 0.7 98.0 +- 1.1 98.1 +- 1.3 108.2 +- 1.3 98.7 +- 0.9 3000 True True True
tl W2 Rep3 96.4 +- 0.7 98.0 +- 1.1 98.5 +- 1.0 108.2 +- 1.3 98.5 +- 0.7 3000 True True True
Result: blocked (methods underspecified). The delayed means do not rise to their peak under any reading (106.6% at 30 Hz, 101.2 to 102.8% at 300 Hz, 106.7 to 108.8% at 1500 Hz), so at best the claim is partially reproduced, on the time-locked window. That window holds under W2 with every repeat reading and under W1 with Rep1 and Rep3: means rise from 30 Hz (96.4 to 96.8%) to a peak at 3000 Hz (60.26 mW/cm2; 107.7 to 108.7%) and fall at 4500 Hz (98.0 to 99.5%). Under W1 Rep2 it fails: the 300 Hz mean (98.2%) is 0.1 point above the 1500 Hz mean (98.1%), against SEMs of 1.1 and 1.3. With one recording per rat and condition (the reading adopted) W1 and W2 thus disagree under Rep2, and the claims file's header blocks a claim whose readings disagree.
fig6fh-fsu-ispta¶
Awake FSUs show no time-locked response, and their delayed response increases almost linearly with ISPTA.
The time-locked part takes the awake outcome of fig5bd-fsu-time-locked-no-prf.
rows = []
for rep in REPS:
m = select("awake", "FSU", rep).groupby("prf").del_awake.mean()
x = [tfus.PRF_ISPTA[p] for p in PRFS]
y = [m[p] for p in PRFS]
rows.append({"reading": rep, **{f"{p} Hz": round(m[p], 1) for p in PRFS},
"Spearman rho": stats.spearmanr(x, y)[0], "Pearson r": stats.pearsonr(x, y)[0]})
print(pd.DataFrame(rows).round(3).to_string(index=False))
o, post = prf_test("awake", "FSU", "tl", "W1", "Rep1")
print(f"\nTime-locked part (from fig5bd, awake, W1 Rep1): Kruskal-Wallis p = {o['p']:.2g}, "
f"{int(post.sig.sum())} significant post hoc pairs")
reading 30 Hz 300 Hz 1500 Hz 3000 Hz 4500 Hz Spearman rho Pearson r Rep1 103.5 98.7 101.8 100.6 100.3 -0.4 -0.268 Rep2 103.5 99.1 101.7 102.1 100.4 -0.3 -0.170 Rep3 103.5 99.1 101.5 102.1 100.4 -0.3 -0.177 Time-locked part (from fig5bd, awake, W1 Rep1): Kruskal-Wallis p = 1.8e-12, 5 significant post hoc pairs
Result: not reproduced. The time-locked part fails (fig5bd: awake Kruskal-Wallis p = 1.8e-12). The delayed FSU means do not increase with ISPTA: 103.5, 98.7, 101.8, 100.6 and 100.3% from 30 to 4500 Hz (Spearman rho = -0.4, Pearson r = -0.27, Rep1; rho = -0.3 Rep2 and Rep3).
fig5a-1500-highest-time-locked-mean¶
Among PRFs, 1500 Hz gives the highest mean RSU time-locked response in awake rats, although not significantly above 3000 and 4500 Hz.
rows = []
for w, rep in READINGS:
o, post = prf_test("awake", "RSU", "tl", w, rep)
rows.append({"reading": f"{w} {rep}", **{f"{p} Hz": round(o["mean"][p], 1) for p in PRFS},
"highest": max(o["mean"], key=o["mean"].get),
"p_bonf 1500 vs 3000": tfus.pair(post, 1500, 3000).p_bonf, "p_bonf 1500 vs 4500": tfus.pair(post, 1500, 4500).p_bonf})
print(pd.DataFrame(rows).round(4).to_string(index=False))
reading 30 Hz 300 Hz 1500 Hz 3000 Hz 4500 Hz highest p_bonf 1500 vs 3000 p_bonf 1500 vs 4500 W1 Rep1 96.8 97.6 98.5 108.2 99.2 3000 0.0000 1.0 W1 Rep2 96.8 98.2 98.1 108.7 98.1 3000 0.0029 1.0 W1 Rep3 96.8 98.2 98.3 108.7 98.0 3000 0.0000 1.0 W2 Rep1 96.4 97.4 98.7 107.7 99.5 3000 0.0000 1.0 W2 Rep2 96.4 98.0 98.1 108.2 98.7 3000 0.0040 1.0 W2 Rep3 96.4 98.0 98.5 108.2 98.5 3000 0.0002 1.0
Result: not reproduced. The highest mean RSU time-locked response is at 3000 Hz under every reading (107.7 to 108.7%), and it is significantly above 1500 Hz (corrected p = 0.0040 or less); 1500 Hz gives 98.1 to 98.7%.
Duty cycle claims (fig7a to fig7d)¶
The DC series is 1500 Hz with TBD 40, 200 and 400 us (6, 30 and 60%). The dandiset holds 6% and 60% recordings for three rats (BH506, BH521, BH526) and no 30% recording of its own: the paper's 30% recordings, made with the DC series at 10 months, are missing. The only 1500 Hz, 200 us recordings are those of the PRF series: the same rats' were made three to five months before the DC recordings, and the other rats' from August 2023 to 6 March 2024, three to four days before the March DC recordings (data check). The claims file names no reading for this; three come from what the dandiset holds: C0, no 30% group (the tests run on 6% and 60% alone, Bonferroni over one pair); C1, the PRF series' 1500 Hz, 200 us recordings of the same rats (BH506 and BH526; BH521 has none); and C2, every awake 1500 Hz, 200 us recording, run with each repeat reading. C1 and C2 are substitutes: a part of a criterion that needs the 30% group is treated as data missing, and does not decide an outcome. Parts that compare with sham cannot be tested.
DCS = [6, 30, 60]
def dc_test(ct, kind, w, reading):
c = col(kind, "awake", w)
d = L[(L.state == "awake") & (L.prf == 1500) & (L.celltype == ct)]
dc_rats = set(d[d.tbd != 200].rat)
if reading == "C1":
d = d[d.rat.isin(dc_rats)]
if reading == "C0":
d = d[d.dc != 30]
if reading == "C2 Rep2":
d = d[~d.repeat]
if reading == "C2 Rep3":
d = d[~d.later]
return tfus.group_test({dc: d[d.dc == dc][c].values for dc in DCS})
C_READINGS = ("C0", "C1", "C2 Rep1", "C2 Rep2", "C2 Rep3")
DC_READINGS = [(w, r) for w in ("W1", "W2") for r in C_READINGS]
def pbonf(post, a, b):
"""Corrected p of a pair, or nan when a group is absent (C0 has no 30% group)."""
return tfus.pair(post, a, b).p_bonf if {a, b} <= set(post.a) | set(post.b) else np.nan
print("Recordings in the DC comparison:")
d = L[(L.state == "awake") & (L.prf == 1500)]
print(d.groupby(["dc", "rat"]).session.unique().to_string())
Recordings in the DC comparison:
dc rat
6 BH506 [BH506_1500_40]
BH521 [BH521_1500_40]
BH526 [BH526_1500_40]
30 BH498 [BH498_1500_200_2, BH498_1500_200]
BH506 [BH506_1500_200]
BH508 [BH508_1500_200_2, BH508_1500_200]
BH512 [BH512_1500_200]
BH519 [BH519_1500_200]
BH526 [BH526_1500_200]
BH527 [BH527_1500_200]
60 BH506 [BH506_1500_400]
BH521 [BH521_1500_400]
BH526 [BH526_1500_400]
fig7a-rsu-time-locked-dc¶
Awake RSU time-locked responses at 1500 Hz PRF are significantly lower at 6% DC than at 60% DC (and, by the Figure 7A legend, than at 30% DC), and every DC is significantly higher than sham.
rows = []
for w, r in DC_READINGS:
o, post = dc_test("RSU", "tl", w, r)
if (w, r) == ("W1", "C1"):
show(o, post, "awake RSU 0-67 ms by DC, W1 C1")
rows.append({"reading": f"{w} {r}", "n": o["n"], **{f"mean {k}%": round(v, 1) for k, v in o["mean"].items()},
"KW p": o["p"], "(1) KW": o["p"] < 0.05, "(2) 6 below 60": lower(post, 6, 60), "(4) 6 below 30": lower(post, 6, 30) if r != "C0" else None,
"p_bonf 6-60": pbonf(post, 6, 60), "p_bonf 6-30": pbonf(post, 6, 30)})
print(pd.DataFrame(rows).to_string(index=False))
print("\nMean by rat and DC (W1):")
print(by_rat("awake", "RSU", "tl", "W1", "Rep1", tbd=40).join(by_rat("awake", "RSU", "tl", "W1", "Rep1", tbd=400)).join(
select("awake", "RSU", "Rep1", prf=1500).groupby("rat").tl_w1.mean().round(1).rename("30 (1500 Hz, 200 us)"), how="outer").to_string())
awake RSU 0-67 ms by DC, W1 C1: Kruskal-Wallis H = 100.5, p = 1.5e-22; n = {6: 90, 30: 72, 60: 111}
means: {6: 180.8, 30: 90.5, 60: 147.5} medians: {6: 122.0, 30: 87.1, 60: 123.5}
a b mean_a mean_b p p_bonf sig
6 30 180.7961 90.5398 0.0000 0.0 True
6 60 180.7961 147.5271 0.9059 1.0 False
30 60 90.5398 147.5271 0.0000 0.0 True
reading n mean 6% mean 60% KW p (1) KW (2) 6 below 60 (4) 6 below 30 p_bonf 6-60 p_bonf 6-30 mean 30%
W1 C0 {6: 90, 60: 111} 180.8 147.5 9.058538e-01 False False None 0.905854 NaN NaN
W1 C1 {6: 90, 30: 72, 60: 111} 180.8 147.5 1.519712e-22 True False False 1.000000 1.558560e-17 90.5
W1 C2 Rep1 {6: 90, 30: 303, 60: 111} 180.8 147.5 1.181829e-47 True False False 1.000000 1.101738e-27 98.5
W1 C2 Rep2 {6: 90, 30: 227, 60: 111} 180.8 147.5 8.654375e-38 True False False 1.000000 2.823132e-23 98.1
W1 C2 Rep3 {6: 90, 30: 236, 60: 111} 180.8 147.5 1.056297e-45 True False False 1.000000 9.851551e-28 98.3
W2 C0 {6: 90, 60: 111} 176.3 147.1 6.660148e-01 False False None 0.666015 NaN NaN
W2 C1 {6: 90, 30: 72, 60: 111} 176.3 147.1 3.599515e-24 True False False 1.000000 3.617519e-18 90.4
W2 C2 Rep1 {6: 90, 30: 303, 60: 111} 176.3 147.1 8.132655e-51 True False False 1.000000 1.664563e-28 98.7
W2 C2 Rep2 {6: 90, 30: 227, 60: 111} 176.3 147.1 3.888363e-41 True False False 1.000000 2.952743e-24 98.1
W2 C2 Rep3 {6: 90, 30: 236, 60: 111} 176.3 147.1 1.321556e-48 True False False 1.000000 1.545196e-28 98.5
Mean by rat and DC (W1):
6 60 30 (1500 Hz, 200 us)
rat
BH498 NaN NaN 103.9
BH506 249.0 209.9 108.9
BH508 NaN NaN 96.1
BH512 NaN NaN 110.7
BH519 NaN NaN 102.5
BH521 113.1 117.3 NaN
BH526 120.8 113.7 71.1
BH527 NaN NaN 100.1
Result: blocked (data missing). Part (2) does not need the 30% group and fails under every reading: 6% is the highest group, not the lowest, at 180.8% against 147.5% at 60% (W1), and 6% does not differ from 60% (corrected p = 1.0; 0.67 to 0.91 under C0). Parts (1) and (4) need the 30% group and part (3) needs the sham, so those three are not tested, and whatever they would show, the claim cannot be reproduced. With the substitutes the 30% group is far below both (90.4 to 98.7%). The 6% and 60% effects come mostly from BH506 (249.0 and 209.9%); BH521 and BH526 give 113 to 121%.
Why it is blocked: the criterion gives partially reproduced when (1) and at least one of (2) to (4) hold. Part (1) is the Kruskal-Wallis over 6, 30 and 60%, which needs the 30% group: without it (C0) the test is not significant (p = 0.91 W1, 0.67 W2), with the substitutes it is (p = 1.5e-22 or smaller), so the substitutes alone decide it. Part (4), 6% below 30%, also needs the 30% group, and part (3) needs the sham. Partial is therefore reachable if (1) and (3) or (4) held, and the data present do not decide it. By the rule for the 30% group stated under Readings (as for fig7c-rsu-delayed-dc) the claim is blocked as data missing. The second verifier pointed this out, and the third agreed.
fig7a-60-highest-time-locked-mean¶
By mean values 60% DC gives the strongest RSU time-locked response, though it is only significantly higher than 6%.
rows = []
for w, r in DC_READINGS:
o, post = dc_test("RSU", "tl", w, r)
rows.append({"reading": f"{w} {r}", **{f"mean {k}%": round(v, 1) for k, v in o["mean"].items()},
"highest": max(o["mean"], key=o["mean"].get), "p_bonf 60 vs 30": pbonf(post, 60, 30)})
print(pd.DataFrame(rows).to_string(index=False))
reading mean 6% mean 60% highest p_bonf 60 vs 30 mean 30%
W1 C0 180.8 147.5 6 NaN NaN
W1 C1 180.8 147.5 6 2.058986e-19 90.5
W1 C2 Rep1 180.8 147.5 6 1.389623e-33 98.5
W1 C2 Rep2 180.8 147.5 6 9.803830e-28 98.1
W1 C2 Rep3 180.8 147.5 6 2.691738e-33 98.3
W2 C0 176.3 147.1 6 NaN NaN
W2 C1 176.3 147.1 6 2.960234e-21 90.4
W2 C2 Rep1 176.3 147.1 6 1.182691e-36 98.7
W2 C2 Rep2 176.3 147.1 6 4.661095e-31 98.1
W2 C2 Rep3 176.3 147.1 6 3.441677e-36 98.5
Result: not reproduced. The highest mean is at 6% DC under every reading (176.3 to 180.8%), above 60% (147.1 to 147.5%); this does not need the 30% group. With the substitutes, 60% is significantly above 30% (corrected p < 1e-16), where the claim says only 6% differs from it.
fig7b-fsu-time-locked-dc¶
FSUs show no time-locked response across DC levels; no significant differences are observed.
A claim of absence, judged on the Kruskal-Wallis across 6, 30 and 60% (G1) at 0.05 and the post hoc tests at Bonferroni-corrected 0.05.
rows = []
for w, r in DC_READINGS:
o, post = dc_test("FSU", "tl", w, r)
if (w, r) == ("W1", "C1"):
show(o, post, "awake FSU 0-67 ms by DC, W1 C1")
rows.append({"reading": f"{w} {r}", "n": o["n"], **{f"median {k}%": round(v, 1) for k, v in o["median"].items()},
**{f"mean {k}%": round(v, 1) for k, v in o["mean"].items()},
"KW p": o["p"], "significant pairs": int(post.sig.sum())})
print(pd.DataFrame(rows).to_string(index=False))
print("\nMean by rat and DC (W1):")
print(by_rat("awake", "FSU", "tl", "W1", "Rep1", tbd=40).join(by_rat("awake", "FSU", "tl", "W1", "Rep1", tbd=400)).to_string())
awake FSU 0-67 ms by DC, W1 C1: Kruskal-Wallis H = 52.8, p = 3.5e-12; n = {6: 63, 30: 67, 60: 75}
means: {6: 156.8, 30: 103.5, 60: 205.3} medians: {6: 115.0, 30: 100.8, 60: 128.6}
a b mean_a mean_b p p_bonf sig
6 30 156.8051 103.5128 0.0000 0.0000 True
6 60 156.8051 205.3385 0.0054 0.0163 True
30 60 103.5128 205.3385 0.0000 0.0000 True
reading n median 6% median 60% mean 6% mean 60% KW p significant pairs median 30% mean 30%
W1 C0 {6: 63, 60: 75} 115.0 128.6 156.8 205.3 5.426807e-03 1 NaN NaN
W1 C1 {6: 63, 30: 67, 60: 75} 115.0 128.6 156.8 205.3 3.475143e-12 3 100.8 103.5
W1 C2 Rep1 {6: 63, 30: 306, 60: 75} 115.0 128.6 156.8 205.3 5.995601e-32 3 100.8 101.8
W1 C2 Rep2 {6: 63, 30: 250, 60: 75} 115.0 128.6 156.8 205.3 2.792563e-28 3 101.4 102.5
W1 C2 Rep3 {6: 63, 30: 272, 60: 75} 115.0 128.6 156.8 205.3 6.800539e-35 3 100.4 100.8
W2 C0 {6: 63, 60: 75} 114.7 127.2 154.5 200.4 5.426807e-03 1 NaN NaN
W2 C1 {6: 63, 30: 67, 60: 75} 114.7 127.2 154.5 200.4 7.468969e-13 3 100.0 103.0
W2 C2 Rep1 {6: 63, 30: 306, 60: 75} 114.7 127.2 154.5 200.4 4.443666e-33 3 100.4 101.7
W2 C2 Rep2 {6: 63, 30: 250, 60: 75} 114.7 127.2 154.5 200.4 1.508681e-29 3 100.9 102.2
W2 C2 Rep3 {6: 63, 30: 272, 60: 75} 114.7 127.2 154.5 200.4 4.627784e-36 3 100.1 100.7
Mean by rat and DC (W1):
dc 6 60
rat
BH506 338.6 403.3
BH521 110.9 125.6
BH526 119.4 107.9
Result: not reproduced. FSU time-locked responses differ across DC under every reading, and the difference between 6% and 60% does not need the 30% group: under C0 the Kruskal-Wallis gives p = 0.0054 (one pair, significant), and with the substitutes p = 3.5e-12 (W1 C1) or smaller, with all three pairs significant. FSU respond strongly at 6% and 60% (means 156.8 and 205.3%, medians 115.0 and 128.6%) and not at the substitute 30% (median 100.8%, C1); the largest responses are in BH506 (338.6 and 403.3%).
fig7c-rsu-delayed-dc¶
Awake RSU delayed responses prefer 30% DC over 6% and 60%, all three DCs differ (Figure 7C legend), and at 60% DC they do not differ from sham.
Part (3) compares with sham and cannot be tested. The window is on bin edges (W1 = W2).
rows = []
for r in C_READINGS:
o, post = dc_test("RSU", "del", "W1", r)
show(o, post, f"awake RSU 0.5-1 s by DC, {r}")
rows.append({"reading": r, "(1) KW": o["p"] < 0.05, "(2) 30 above 6 and 60": (lower(post, 6, 30) and lower(post, 60, 30)) if r != "C0" else None,
"(4) 6 and 60 differ": sig(post, 6, 60)})
print(pd.DataFrame(rows).to_string(index=False))
awake RSU 0.5-1 s by DC, C0: Kruskal-Wallis H = 5.1, p = 0.024; n = {6: 90, 60: 111}
means: {6: 99.3, 60: 100.6} medians: {6: 99.0, 60: 100.1}
a b mean_a mean_b p p_bonf sig
6 60 99.277 100.6097 0.0236 0.0236 True
awake RSU 0.5-1 s by DC, C1: Kruskal-Wallis H = 3.2, p = 0.21; n = {6: 90, 30: 72, 60: 111}
means: {6: 99.3, 30: 104.9, 60: 100.6} medians: {6: 99.0, 30: 98.7, 60: 100.1}
a b mean_a mean_b p p_bonf sig
6 30 99.2770 104.8808 0.8187 1.0000 False
6 60 99.2770 100.6097 0.0236 0.0709 False
30 60 104.8808 100.6097 0.8808 1.0000 False
awake RSU 0.5-1 s by DC, C2 Rep1: Kruskal-Wallis H = 16.8, p = 0.00022; n = {6: 90, 30: 303, 60: 111}
means: {6: 99.3, 30: 106.7, 60: 100.6} medians: {6: 99.0, 30: 103.7, 60: 100.1}
a b mean_a mean_b p p_bonf sig
6 30 99.2770 106.7424 0.0006 0.0017 True
6 60 99.2770 100.6097 0.0236 0.0709 False
30 60 106.7424 100.6097 0.0069 0.0206 True
awake RSU 0.5-1 s by DC, C2 Rep2: Kruskal-Wallis H = 20.4, p = 3.7e-05; n = {6: 90, 30: 227, 60: 111}
means: {6: 99.3, 30: 108.8, 60: 100.6} medians: {6: 99.0, 30: 105.8, 60: 100.1}
a b mean_a mean_b p p_bonf sig
6 30 99.2770 108.8265 0.0002 0.0006 True
6 60 99.2770 100.6097 0.0236 0.0709 False
30 60 108.8265 100.6097 0.0010 0.0031 True
awake RSU 0.5-1 s by DC, C2 Rep3: Kruskal-Wallis H = 17.1, p = 0.00019; n = {6: 90, 30: 236, 60: 111}
means: {6: 99.3, 30: 107.3, 60: 100.6} medians: {6: 99.0, 30: 103.9, 60: 100.1}
a b mean_a mean_b p p_bonf sig
6 30 99.2770 107.3363 0.0005 0.0014 True
6 60 99.2770 100.6097 0.0236 0.0709 False
30 60 107.3363 100.6097 0.0054 0.0161 True
reading (1) KW (2) 30 above 6 and 60 (4) 6 and 60 differ
C0 True None True
C1 False False False
C2 Rep1 True True False
C2 Rep2 True True False
C2 Rep3 True True False
Result: blocked (data missing). Part (2), 30% above 6% and 60%, needs the paper's 30% recordings, which are not in the dandiset, and the substitutes decide it: with the same rats' 1500 Hz, 200 us recordings (C1) the Kruskal-Wallis is not significant (p = 0.21) and no pair differs; with every such recording (C2) the effect is significant (p = 2.2e-04 Rep1, 3.7e-05 Rep2, 1.9e-04 Rep3) and 30% is above 6% and 60%. Part (4), 6% against 60%, also depends on the 30% group through the Bonferroni family: 6% is below 60% at p = 0.024 alone (C0), and at corrected p = 0.071 over three pairs. Part (3) needs the sham.
fig7d-fsu-delayed-dc¶
Awake FSU delayed responses are significantly higher at 60% DC than at 6%, and at 6% DC they do not differ from sham.
Part (3) compares with sham and cannot be tested.
rows = []
for r in C_READINGS:
o, post = dc_test("FSU", "del", "W1", r)
show(o, post, f"awake FSU 0.5-1 s by DC, {r}")
rows.append({"reading": r, "(1) KW": o["p"] < 0.05, "(2) 60 above 6": lower(post, 6, 60)})
print(pd.DataFrame(rows).to_string(index=False))
awake FSU 0.5-1 s by DC, C0: Kruskal-Wallis H = 7.9, p = 0.0051; n = {6: 63, 60: 75}
means: {6: 98.2, 60: 99.6} medians: {6: 97.5, 60: 99.7}
a b mean_a mean_b p p_bonf sig
6 60 98.2485 99.625 0.0051 0.0051 True
awake FSU 0.5-1 s by DC, C1: Kruskal-Wallis H = 8.9, p = 0.011; n = {6: 63, 30: 67, 60: 75}
means: {6: 98.2, 30: 99.2, 60: 99.6} medians: {6: 97.5, 30: 98.4, 60: 99.7}
a b mean_a mean_b p p_bonf sig
6 30 98.2485 99.220 0.0258 0.0774 False
6 60 98.2485 99.625 0.0051 0.0152 True
30 60 99.2200 99.625 0.3942 1.0000 False
awake FSU 0.5-1 s by DC, C2 Rep1: Kruskal-Wallis H = 29.2, p = 4.5e-07; n = {6: 63, 30: 306, 60: 75}
means: {6: 98.2, 30: 101.8, 60: 99.6} medians: {6: 97.5, 30: 100.8, 60: 99.7}
a b mean_a mean_b p p_bonf sig
6 30 98.2485 101.8436 0.0000 0.0000 True
6 60 98.2485 99.6250 0.0051 0.0152 True
30 60 101.8436 99.6250 0.0141 0.0422 True
awake FSU 0.5-1 s by DC, C2 Rep2: Kruskal-Wallis H = 26.9, p = 1.4e-06; n = {6: 63, 30: 250, 60: 75}
means: {6: 98.2, 30: 101.7, 60: 99.6} medians: {6: 97.5, 30: 100.7, 60: 99.7}
a b mean_a mean_b p p_bonf sig
6 30 98.2485 101.6684 0.0000 0.0000 True
6 60 98.2485 99.6250 0.0051 0.0152 True
30 60 101.6684 99.6250 0.0194 0.0581 False
awake FSU 0.5-1 s by DC, C2 Rep3: Kruskal-Wallis H = 28.2, p = 7.6e-07; n = {6: 63, 30: 272, 60: 75}
means: {6: 98.2, 30: 101.5, 60: 99.6} medians: {6: 97.5, 30: 100.7, 60: 99.7}
a b mean_a mean_b p p_bonf sig
6 30 98.2485 101.5019 0.0000 0.0000 True
6 60 98.2485 99.6250 0.0051 0.0152 True
30 60 101.5019 99.6250 0.0221 0.0664 False
reading (1) KW (2) 60 above 6
C0 True True
C1 True True
C2 Rep1 True True
C2 Rep2 True True
C2 Rep3 True True
Result: partially reproduced. Parts (1) and (2) hold under every reading, C0 included, so they do not depend on the substitute 30% group: the Kruskal-Wallis is significant (p = 0.0051 C0, 0.011 C1, 4.5e-07 to 1.4e-06 C2) and 60% is above 6% (p = 0.0051 alone, corrected p = 0.015 over three pairs). Part (3), 6% against sham, cannot be tested, so the claim is not judged reproduced. The difference is small and both groups are below baseline (means 99.6% at 60% and 98.2% at 6%).
sham-no-auditory-response (blocked)¶
The off-target sham at 1500 Hz PRF produced no significant response, so the activation is not attributed to an auditory or skull-conduction confound.
The criterion blocks this claim (the paper names no test of the sham), and the dandiset holds no sham recording, so the sham's responses cannot be shown either.
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 |
|---|---|
| fig3-cell-type-sorting | Blocked: data missing (no waveforms for the k-means features); label counts match in one or two recordings, no more than chance |
| fig4-rsu-time-locked-both-states | Blocked: data missing (no sham); the awake mean is below 100% |
| fig4hi-fsu-delayed-latency | Not reproduced |
| fig4h-awake-latency-wider | Blocked: methods underspecified (D1 and D2 disagree) |
| fig5a-rsu-time-locked-prf-awake | Not reproduced |
| fig5c-rsu-time-locked-prf-anesth | Not reproduced |
| rsu-time-locked-anesthesia-factor | Not reproduced |
| fig5bd-fsu-time-locked-no-prf | Not reproduced |
| fig5e-rsu-delayed-prf-awake | Blocked: data missing (no sham; partial or not reproduced depends on it) |
| fig5f-fsu-delayed-prf-awake | Not reproduced |
| fig5gh-delayed-prf-anesth | Not reproduced |
| delayed-anesthesia-factor | Reproduced (assumption: one recording per rat and condition, Rep2 or Rep3) |
| fig6ad-time-adjacent-batch | Blocked: data missing (second batch) |
| fig6ef-rsu-ispta-rise-then-fall | Blocked: methods underspecified (W1 and W2 disagree under Rep2; at best partial) |
| fig6fh-fsu-ispta | Not reproduced |
| fig5a-1500-highest-time-locked-mean | Not reproduced |
| fig7a-rsu-time-locked-dc | Blocked: data missing (partial or not reproduced depends on the 30% group and the sham) |
| fig7a-60-highest-time-locked-mean | Not reproduced |
| fig7b-fsu-time-locked-dc | Not reproduced |
| fig7c-rsu-delayed-dc | Blocked: data missing (the paper's 30% DC recordings) |
| fig7d-fsu-delayed-dc | Partially reproduced (the sham part cannot be tested) |
| sham-no-auditory-response | Blocked: methods underspecified, and no sham data |
How the analysis differs from the paper's. Cell types are the dataset's labels; the k-means
on waveform features cannot be run, as there are no waveforms. The files hold 500 trials (the
paper's 505 less five), and five recordings hold 203 to 339. No sham recording exists, so every
comparison with sham is left out, and no claim with a sham part is judged reproduced. The
anesthetized PRF series has labelled units from two rats (the paper: ten), and BH512 has no 30 Hz
recording, so rat and PRF are confounded under anesthesia. The DC series has 6% and 60% recordings
from three rats and no 30% recording of its own; the PRF series' 1500 Hz, 200 us recordings stand
in for it, without deciding any outcome. Ten awake recordings repeat a condition in the same rat;
results are shown with all of them (Rep1), without the _2 recordings (Rep2) and with the
earliest recording of each pair (Rep3). No unit was excluded for a fluctuating baseline or for
movement, as the dataset marks none. Most 200 us recordings show no time-locked population
response, which limits every time-locked claim.
What was found about the dataset. Eleven of the 75 files are second copies of other files:
the anesthetized recordings of BH498 and BH512 are each deposited twice (with labels under
subject folders BH497 and BH514, without labels under BH498 and BH512), BH512's 1500 Hz
anesthetized recording a third time, and one 30 Hz recording is filed under both BH498 and
BH508, with start times 6.5 hours apart. The condition of a recording is coded only in its
identifier. The dandiset's description says the repeats were taken months after the originals,
but three of the ten _2 recordings are earlier than their unsuffixed pair. The 2.2 s
stop_time is not described. Fifteen anesthetized files have no cell-type labels. Every
file gives age P24W, where the paper gives 7 to 10 months and the description 8 to 12, and the
labelled anesthetized copies are described as awake. The second batch of five rats, the sham
condition and the 30% DC recordings are not in the dandiset. In most 200 us recordings no
time-locked population response can be detected (10 of 64 recordings show one at t > 3, among
them 4 of the 6 awake 40 and 400 us recordings), while the 6% and 60% DC recordings respond
strongly, FSU included.
print(f"notebook ran in {time.time() - started:.0f} s")
notebook ran in 50 s