Reproducing Gao, Ramachandran, Yu & He (2025) from dandiset 000946¶
Paper. Gao H, Ramachandran S, Yu K, He B (2025). Transcranial focused ultrasound modulates feedforward
and feedback cortico-thalamo-cortical pathways by selectively activating excitatory neurons. J Neurosci
45(23): e2218242025. doi:10.1523/jneurosci.2218-24.2025.
The claims file (claims/000946/10.1523-jneurosci.2218-24.2025.yaml) was written from the version of
record as served by PubMed Central (PMC12139594): main text, Methods and figure legends; the paper has no
supplement and links no analysis code.
Dandiset. DANDI:000946, version 0.250829.1249. The paper cites the dandiset without a version (Methods, Electrophysiological recordings: "All experimental data has been made public in a repository at the DANDI Archive (https://dandiarchive.org/dandiset/000946)"), so the most recent published version is pinned. The data check below compares its file list with the earliest published version, 0.250711.1640; all six published versions list the same 1920 files.
What the dataset holds. One NWB file per recorded area and session, each with a units table (spike
times and a celltype_label column) and a trials table. There are no spike waveforms in any file, and
celltype_label, described as "label 1 = RSU, label 2 = FSU", is 0 for every unit. The paper splits units
into RSUs and FSUs by k-means on waveform measures (Methods, Neuronal type classification), so every claim
about RSUs or FSUs, and every claim built on them, cannot be tested. No session is marked as sham or as
tFUS aimed at CG1. The unit totals of the two count claims do not need the types, and they are tested on
the totals. The auditory cortex units of Figure 11B are not split by type, and that claim is tested.
| Claim | Tested here |
|---|---|
| fig4-unit-counts | Tested on the area totals (the RSU/FSU split cannot be made) |
| fig4a-s1-rsu-time-locked-3000 | Blocked: data missing (no cell types) |
| fig4b-pom-rsu-time-locked-3000 | Blocked: data missing (no cell types) |
| fig4-postsonication-3000 | Blocked: data missing (no cell types) |
| fig4-30hz-responses | Blocked: data missing (no cell types) |
| fig4-fsu-no-change | Blocked: data missing (no cell types) |
| fig5a-sham-no-change | Blocked: data missing (no cell types, no sham sessions marked) |
| fig5b-cg1-control-no-change | Blocked: data missing (no cell types, no CG1 sessions marked) |
| fig6-unit-counts | Tested on the area totals (the RSU/FSU split cannot be made) |
| fig6a-s1-rsu-time-locked-prf, fig6a-s1-rsu-delayed-prf, fig6b-s1-fsu-prf, fig6c-pom-rsu-prf, fig6d-pom-fsu-prf | Blocked: data missing (no cell types) |
| fig7a to fig7d, fig8a, fig8b, fig9a, fig9b (8 claims) | Blocked: data missing (no cell types) |
| fig10a-dsttc-3000, fig10b-dsttc-30, fig10c to fig10f (6 claims) | Blocked: data missing (no cell types) |
| pathway-follows-local-prf | Blocked: depends on fig6a and fig10c |
| fig11b-ac-responses | Tested under every reading |
| fig11c-ac-to-s1-30hz, fig11d-ac-to-s1-3000 | Blocked: data missing (no S1 cell types) |
| auditory-confound-low-prf | Blocked: depends on fig11c and fig11d |
No step of the analysis is random, so no seed is needed.
import json
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
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()
print("dandiset", tfus.DANDISET_ID, "version", tfus.VERSION)
dandiset 000946 version 0.250829.1249
Data check¶
Every file of the pinned version is opened and its metadata, trial count and per-unit spike counts are read (not the spike times). Each check of the metadata against the values is stated under its output.
assets = tfus.list_assets()
early = tfus.list_assets("0.250711.1640")
same = {(a["path"], a["size"]) for a in assets} == {(a["path"], a["size"]) for a in early}
print(f"{len(assets)} files in {tfus.VERSION}; {len(early)} in 0.250711.1640; same paths and sizes: {same}")
t0 = time.time()
inv = tfus.inventory(assets)
print(f"read the headers of {len(inv)} files in {time.time() - t0:.0f} s; files that failed to open: "
f"{(inv.error != '').sum()}")
inv["kind"] = np.where(inv.path.str.startswith("sub-STTC"), "STTC", "units")
print(inv.kind.value_counts().to_string())
1920 files in 0.250829.1249; 1920 in 0.250711.1640; same paths and sizes: True
read the headers of 1920 files in 72 s; files that failed to open: 0 kind units 1550 STTC 370
u = inv[inv.kind == "units"].copy()
print("Columns of the units tables:")
print(u.unit_columns.value_counts().to_string())
print("\nNon-empty acquisition, analysis or processing groups (where waveforms could be):",
int((u.groups != "").sum()), "of", len(u), "files")
print("\nValues of celltype_label over all files (absent where the column is missing):")
print(u.labels.value_counts().to_string())
n_units_all = int(u.n_units.sum())
print(f"\n{n_units_all} units in {len(u)} unit files, before removing duplicate files")
Columns of the units tables: unit_columns spike_times,celltype_label 1516 spike_times 34 Non-empty acquisition, analysis or processing groups (where waveforms could be): 0 of 1550 files Values of celltype_label over all files (absent where the column is missing): labels [0.0] 1516 absent 34 45124 units in 1550 unit files, before removing duplicate files
The units tables hold only spike times and celltype_label. No file has waveforms, mean waveforms or
anything under acquisition, analysis or processing, and celltype_label is 0 in every file where it is
present (it is missing from the rest), though its description says 1 is RSU and 2 is FSU. The RSU/FSU
split of the paper cannot be made from these data.
kept, dropped = tfus.dedupe(u)
print(f"{len(u)} unit files; {len(dropped)} have the same subject, identifier, area and per-unit spike counts "
f"as another file; {len(kept)} distinct recordings remain")
print("Files per group of duplicates:", (dropped.kept.value_counts() + 1).value_counts().to_dict())
ids = dict(zip(u.path, u.asset_id))
t0 = time.time()
h_kept = dict(zip(dropped.kept.unique(), tfus.map_sessions(tfus.spike_hash, [ids[p] for p in dropped.kept.unique()])))
h_drop = tfus.map_sessions(tfus.spike_hash, [ids[p] for p in dropped.dropped])
ident = sum(h_kept[k] == h for k, h in zip(dropped.kept, h_drop))
print(f"duplicate files whose spike times and trial starts are identical to the kept file: {ident} of {len(dropped)} "
f"({time.time() - t0:.0f} s)")
st = u.set_index("path").session_start
gap = pd.Series([abs(pd.Timestamp(st[a]) - pd.Timestamp(st[b])).total_seconds()
for a, b in zip(dropped.kept, dropped.dropped)])
print("session_start_time difference within a duplicate pair (s): median", gap.median(), "max", gap.max())
print(dropped.head(6).to_string(index=False))
1550 unit files; 734 have the same subject, identifier, area and per-unit spike counts as another file; 816 distinct recordings remain
Files per group of duplicates: {2: 730, 3: 2}
duplicate files whose spike times and trial starts are identical to the kept file: 734 of 734 (69 s)
session_start_time difference within a duplicate pair (s): median 0.0 max 57.0
kept dropped
sub-BH494/sub-BH494_ses-20230817T132800_ecephys.nwb sub-BH494/sub-BH494_ses-20230817T132836_ecephys.nwb
sub-BH494/sub-BH494_ses-20230820T123700_ecephys.nwb sub-BH494/sub-BH494_ses-20230820T123756_ecephys.nwb
sub-BH494/sub-BH494_ses-20230820T125000_ecephys.nwb sub-BH494/sub-BH494_ses-20230820T125042_ecephys.nwb
sub-BH494/sub-BH494_ses-20230820T131000_ecephys.nwb sub-BH494/sub-BH494_ses-20230820T131032_ecephys.nwb
sub-BH494/sub-BH494_ses-20230820T144500_ecephys.nwb sub-BH494/sub-BH494_ses-20230820T144522_ecephys.nwb
sub-BH494/sub-BH494_ses-20230820T145800_ecephys.nwb sub-BH494/sub-BH494_ses-20230820T145857_ecephys.nwb
Most recordings are deposited twice, under two session start times a few seconds to a minute apart, with identical spike times and trial starts. Each recording is counted once from here on.
k = kept.copy()
print("Distinct recordings by recorded area and tFUS target (from general/notes):")
print(pd.crosstab(k.area.replace("", "(empty)"), k.stim.replace("", "(empty)"), margins=True).to_string())
print("\nRats:", k.subject_id.nunique(), sorted(k.subject_id.unique()))
groups = k.groupby("subject_id").area.agg(lambda a: ",".join(sorted(set(a) - {""})))
print("\nAreas recorded per rat:")
print(groups.value_counts().to_string())
nodc = k.prf_dc.isna() | k.pd_dc.isna()
print("\nRecordings whose area is empty, by rat:", k[k.area == ""].subject_id.value_counts().to_dict())
print("Recordings whose general/data_collection is empty:", int((k.data_collection.fillna("").str.strip() == "").sum()),
"; by rat:", k[k.data_collection.fillna("").str.strip() == ""].subject_id.value_counts().to_dict())
print("Recordings without a PRF or PD in general/data_collection:", int(nodc.sum()), "of", len(k), "; by rat:",
k[nodc].subject_id.value_counts().to_dict())
both = ~nodc & k.prf_id.notna()
agree = (k.prf_dc[both] == k.prf_id[both]) & (np.isclose(k.pd_dc[both], k.pd_id[both]))
print("Recordings where data_collection and the identifier both give PRF and PD:", int(both.sum()),
"; they agree in", int(agree.sum()))
print(k[both][~agree][["identifier", "data_collection"]].to_string(index=False))
print("Recordings with no PRF or PD in either field:", int((k.prf.isna() | k.pd_us.isna()).sum()),
k[k.prf.isna() | k.pd_us.isna()].identifier.tolist())
print("\nIdentifier suffixes after <rat>_<PRF>_<PD>_<UD>_<n>:")
print(k.identifier.str.split("_").str[5:].str.join("_").replace("", "(none)").value_counts().to_string())
print("\nAny identifier, note or data_collection naming sham, no tFUS, CG1 or cingulate:",
int(k[["identifier", "notes", "data_collection"]].apply(
lambda c: c.str.contains("(?i)sham|cg1|cingul|no tfus|control")).any(axis=1).sum()))
print("\nPD values (us) of the S1 and POm recordings, by tFUS target; recordings in brackets:")
for area in ["S1", "POM"]:
for stim, g in k[k.area == area].groupby("stim"):
print(f" {area}, tFUS at {stim or '(empty)'}:", g.pd_us.value_counts().sort_index().to_dict())
Distinct recordings by recorded area and tFUS target (from general/notes):
stim M1 POM S1 All
area
(empty) 2 0 3 5
AC 0 0 152 152
AC_1 0 0 1 1
AC_2 0 0 1 1
M1 44 14 44 102
M1_R 0 0 1 1
POM 31 19 102 152
S1 75 32 289 396
S1BF 0 0 3 3
S1_1 0 0 1 1
S1_R 0 0 1 1
VPM 0 0 1 1
All 152 65 599 816
Rats: 20 ['BH494', 'BH495', 'BH499', 'BH500', 'BH501', 'BH510', 'BH511', 'BH513', 'BH528', 'BH530', 'BH532', 'BH535', 'BH537', 'BH538', 'BH539', 'BH547', 'BH551', 'BH553', 'BH558', 'BH594']
Areas recorded per rat:
area
M1,POM,S1 6
AC,S1 6
POM,S1 3
M1,M1_R,S1,S1_R 1
M1,POM,S1,S1BF,VPM 1
M1,S1 1
AC,AC_1,AC_2,S1,S1_1 1
S1 1
Recordings whose area is empty, by rat: {'BH500': 5}
Recordings whose general/data_collection is empty: 33 ; by rat: {'BH494': 25, 'BH499': 6, 'BH500': 2}
Recordings without a PRF or PD in general/data_collection: 33 of 816 ; by rat: {'BH494': 25, 'BH499': 6, 'BH500': 2}
Recordings where data_collection and the identifier both give PRF and PD: 782 ; they agree in 782
Empty DataFrame
Columns: [identifier, data_collection]
Index: []
Recordings with no PRF or PD in either field: 0 []
Identifier suffixes after <rat>_<PRF>_<PD>_<UD>_<n>:
identifier
(none) 504
MOTOR 88
POM_S1 68
POM_MOTOR 62
POM_POM 34
POM 21
2 6
POM_S1BF 4
NEW 4
POM_150-95 4
NEAR 2
NEAR_MOTOR 2
VPM_S1BF 2
POM_POM_2 2
POM_155-95 2
POM__150-95 2
back 2
contra 2
dep05 2
front 2
POM _S1 1
Any identifier, note or data_collection naming sham, no tFUS, CG1 or cingulate: 0
PD values (us) of the S1 and POm recordings, by tFUS target; recordings in brackets:
S1, tFUS at M1: {200.0: 75}
S1, tFUS at POM: {200.0: 32}
S1, tFUS at S1: {1.4: 7, 2.0: 7, 4.0: 9, 14.0: 9, 20.0: 17, 40.0: 9, 67.0: 10, 100.0: 14, 134.0: 8, 200.0: 107, 300.0: 9, 400.0: 8, 600.0: 9, 1000.0: 10, 2000.0: 18, 3000.0: 9, 10000.0: 10, 20000.0: 9, 30000.0: 10}
POM, tFUS at M1: {200.0: 31}
POM, tFUS at POM: {200.0: 19}
POM, tFUS at S1: {1.4: 2, 2.0: 2, 4.0: 3, 14.0: 3, 20.0: 5, 40.0: 3, 67.0: 3, 100.0: 3, 134.0: 3, 200.0: 48, 300.0: 3, 400.0: 3, 600.0: 3, 1000.0: 3, 2000.0: 4, 3000.0: 2, 10000.0: 3, 20000.0: 3, 30000.0: 3}
The paper used 16 rats with S1 and POm probes and 12 with S1 and AC probes (Methods, Animal preparation).
The dataset has 20 rats: 10 with POm recordings and 7 with AC recordings; of the other three, BH494 has S1
and M1, BH500 has S1, M1 and recordings whose area is left empty, and BH594 has a single S1 file. It also
holds recordings of M1, and sessions with tFUS aimed at M1 or at POm, which the paper does not describe.
The PRF and PD are taken from general/data_collection, and from the identifier where that field is empty
(the counts above say where, and whether the two agree where both are given). None is marked as sham or
as aimed at CG1, so the sham and control sessions of Figure 5 cannot be found.
k["n_trials"] = k.n_trials.astype(int)
print("Trials per recording:", k.n_trials.describe()[["min", "50%", "max"]].to_dict(),
"; recordings with 500 trials:", int((k.n_trials == 500).sum()), "of", len(k))
s1 = k[(k.stim == "S1") & (k.pd_us == 200) & k.area.isin(["S1", "POM", "AC"])].copy()
print("\nPD 200 us recordings with tFUS aimed at S1, by area and PRF (Hz):")
print(pd.crosstab(s1.area, s1.prf).to_string())
t0 = time.time()
res = tfus.map_sessions(tfus.session_counts, list(s1.asset_id))
s1["res"] = res
print(f"\nloaded spike times of {len(s1)} recordings in {time.time() - t0:.0f} s")
iti = np.concatenate([np.diff(r["start"]) for r in res])
print("Interval between trial starts (s): min %.3f, median %.3f, max %.3f; paper 2.5 s with 10%% jitter"
% (iti.min(), np.median(iti), iti.max()))
dur = np.concatenate([np.asarray(tfus.load_session(s1.asset_id.iloc[0])["stop"]) - s1.res.iloc[0]["start"]])
print("Trial stop minus start in one recording (s):", np.unique(np.round(dur, 3)))
Trials per recording: {'min': 129.0, '50%': 500.0, 'max': 501.0} ; recordings with 500 trials: 693 of 816
PD 200 us recordings with tFUS aimed at S1, by area and PRF (Hz):
prf 30.0 300.0 1500.0 3000.0 4500.0
area
AC 3 3 4 3 4
POM 10 7 10 11 10
S1 22 19 22 22 22
loaded spike times of 172 recordings in 28 s Interval between trial starts (s): min 2.250, median 2.500, max 2.750; paper 2.5 s with 10% jitter
Trial stop minus start in one recording (s): [2.2]
# Are the two areas of a session recorded together? Compare the trial starts of files that share an identifier.
rows = []
for ident, g in s1.groupby("identifier"):
if len(g) < 2:
continue
starts = [r["start"] for r in g.res]
a, b = starts[0], starts[1]
same_n = len(a) == len(b)
rows.append(dict(identifier=ident, areas="+".join(g.area), n=(len(a), len(b)),
same_times=bool(same_n and np.allclose(a, b, atol=1e-3)),
same_intervals=bool(same_n and np.allclose(np.diff(a), np.diff(b), atol=1e-3))))
pairs = pd.DataFrame(rows)
print(f"{len(pairs)} identifiers with two areas; same trial start times: {pairs.same_times.sum()}; "
f"same intervals between trials: {pairs.same_intervals.sum()}")
dstart = []
for ident, g in s1.groupby("identifier"):
if len(g) == 2:
a, b = (pd.Timestamp(x) for x in g.session_start)
dstart.append(abs(a - b).total_seconds() / 60)
print("session_start_time difference between the two area files of a session (min): median %.1f, max %.1f"
% (np.median(dstart), np.max(dstart)))
print(pairs[~pairs.same_intervals].head(10).to_string(index=False))
60 identifiers with two areas; same trial start times: 60; same intervals between trials: 60 session_start_time difference between the two area files of a session (min): median 15.5, max 3964.0 Empty DataFrame Columns: [identifier, areas, n, same_times, same_intervals] Index: []
Clock check¶
start_time of the trials table is taken as sonication onset. If it is, the units' firing should change
at fixed delays after it; with the trial starts shifted by 1 s (inside the 2.25 s minimum interval, where
no tFUS is delivered), the same measure should show no change. The measure is the mean over units of the
normalized rate (reading N1) in each 67 ms window, for every unit with more than 2000 spikes.
def unit_table(frame, min_spikes=tfus.MIN_SPIKES):
rows = []
for _, r in frame.iterrows():
res = r["res"]
for c, cs, ps in zip(res["true"], res["shifted"], res["psth"]):
if c["n_spikes"] <= min_spikes:
continue
rows.append(dict(subject=r.subject_id, identifier=r.identifier, area=r.area, prf=r.prf,
true=c, shifted=cs, psth=ps))
return pd.DataFrame(rows)
units = unit_table(s1)
print("Units with more than 2000 spikes, PD 200 us, tFUS at S1:")
print(pd.crosstab(units.area, units.prf, margins=True).to_string())
rows = []
for (area, prf), g in units.groupby(["area", "prf"]):
for align in ["true", "shifted"]:
vals = {w: np.nanmean([tfus.normalized_rates(c, w, "N1") for c in g[align]]) for w in tfus.WINDOWS}
rows.append(dict(area=area, prf=prf, onsets=align, n=len(g), **{k: round(v, 1) for k, v in vals.items()}))
clock = pd.DataFrame(rows)
print("\nMean normalized rate (%) by window, at the trial starts and 1 s later:")
print(clock.to_string(index=False))
dev = clock.assign(dev=clock[list(tfus.WINDOWS)].sub(100).abs().max(axis=1)).groupby("onsets").dev.median()
print("\nMedian over area and PRF of the largest departure from 100%:", dev.round(1).to_dict())
Units with more than 2000 spikes, PD 200 us, tFUS at S1: prf 30.0 300.0 1500.0 3000.0 4500.0 All area AC 112 89 138 167 127 633 POM 315 267 246 352 290 1470 S1 530 388 493 572 431 2414 All 957 744 877 1091 848 4517
Mean normalized rate (%) by window, at the trial starts and 1 s later:
area prf onsets n 0-67 67-134 134-201 201-268
AC 30.0 true 112 105.9 115.9 138.0 176.4
AC 30.0 shifted 112 92.1 93.4 96.1 94.5
AC 300.0 true 89 99.7 94.9 97.5 94.3
AC 300.0 shifted 89 106.2 99.5 95.1 93.8
AC 1500.0 true 138 140.3 192.8 149.2 147.9
AC 1500.0 shifted 138 111.5 102.6 103.8 109.1
AC 3000.0 true 167 107.1 108.0 223.7 310.5
AC 3000.0 shifted 167 92.5 91.4 88.2 90.1
AC 4500.0 true 127 108.0 144.2 228.8 210.8
AC 4500.0 shifted 127 95.0 88.9 95.4 84.0
POM 30.0 true 315 97.1 99.7 110.7 110.3
POM 30.0 shifted 315 101.5 102.5 99.3 101.1
POM 300.0 true 267 104.1 105.7 114.8 110.8
POM 300.0 shifted 267 99.6 97.9 98.9 102.3
POM 1500.0 true 246 107.7 119.7 157.2 131.5
POM 1500.0 shifted 246 104.9 107.4 102.5 104.4
POM 3000.0 true 352 105.4 100.5 113.4 121.5
POM 3000.0 shifted 352 101.0 97.3 100.6 98.2
POM 4500.0 true 290 99.5 102.8 117.5 128.9
POM 4500.0 shifted 290 101.7 101.5 102.5 103.1
S1 30.0 true 530 100.1 96.0 101.3 105.8
S1 30.0 shifted 530 99.4 102.6 101.0 100.4
S1 300.0 true 388 107.2 102.4 99.6 106.2
S1 300.0 shifted 388 101.3 100.7 101.7 100.6
S1 1500.0 true 493 126.6 104.5 107.6 111.4
S1 1500.0 shifted 493 95.8 94.2 95.6 98.5
S1 3000.0 true 572 114.0 102.4 119.0 178.2
S1 3000.0 shifted 572 96.5 96.3 94.7 93.7
S1 4500.0 true 431 102.4 103.3 115.9 148.6
S1 4500.0 shifted 431 100.8 101.3 99.8 98.1
Median over area and PRF of the largest departure from 100%: {'shifted': 5.8, 'true': 28.9}
centers = (tfus.PSTH_EDGES[:-1] + tfus.PSTH_EDGES[1:]) / 2 * 1000
fig, axes = plt.subplots(1, 3, figsize=(15, 3.5), sharex=True)
for ax, area in zip(axes, ["S1", "POM", "AC"]):
for prf, g in units[units.area == area].groupby("prf"):
m = np.mean(np.vstack(g.psth), axis=0) / 0.01
ax.plot(centers, m, label=f"{int(prf)} Hz (n={len(g)})")
print(f"{area} {int(prf)} Hz: mean rate {m.min():.1f} to {m.max():.1f} spikes/s; "
f"[-500, 0) ms mean {m[centers < 0].mean():.1f}, [0, 300) ms max {m[(centers > 0) & (centers < 300)].max():.1f}")
ax.axvspan(0, 67, color="0.9")
ax.set_title(f"{area}, PD 200 us")
ax.set_xlabel("time from trial start (ms)")
axes[0].set_ylabel("mean rate (spikes/s)")
axes[-1].legend(fontsize=7)
plt.tight_layout()
plt.show()
S1 30 Hz: mean rate 11.4 to 12.7 spikes/s; [-500, 0) ms mean 12.0, [0, 300) ms max 12.7 S1 300 Hz: mean rate 11.4 to 13.3 spikes/s; [-500, 0) ms mean 11.8, [0, 300) ms max 13.0 S1 1500 Hz: mean rate 11.7 to 15.2 spikes/s; [-500, 0) ms mean 12.0, [0, 300) ms max 15.2 S1 3000 Hz: mean rate 11.8 to 16.9 spikes/s; [-500, 0) ms mean 12.3, [0, 300) ms max 16.9 S1 4500 Hz: mean rate 11.9 to 15.9 spikes/s; [-500, 0) ms mean 12.2, [0, 300) ms max 15.9 POM 30 Hz: mean rate 12.5 to 14.7 spikes/s; [-500, 0) ms mean 13.0, [0, 300) ms max 14.7 POM 300 Hz: mean rate 11.6 to 14.0 spikes/s; [-500, 0) ms mean 12.3, [0, 300) ms max 14.0 POM 1500 Hz: mean rate 11.2 to 15.3 spikes/s; [-500, 0) ms mean 11.8, [0, 300) ms max 15.3 POM 3000 Hz: mean rate 11.5 to 14.0 spikes/s; [-500, 0) ms mean 12.0, [0, 300) ms max 14.0 POM 4500 Hz: mean rate 12.5 to 16.7 spikes/s; [-500, 0) ms mean 13.1, [0, 300) ms max 16.7 AC 30 Hz: mean rate 7.2 to 13.8 spikes/s; [-500, 0) ms mean 7.8, [0, 300) ms max 13.8 AC 300 Hz: mean rate 8.0 to 9.8 spikes/s; [-500, 0) ms mean 8.7, [0, 300) ms max 9.6 AC 1500 Hz: mean rate 6.9 to 11.5 spikes/s; [-500, 0) ms mean 7.7, [0, 300) ms max 11.5 AC 3000 Hz: mean rate 6.1 to 12.9 spikes/s; [-500, 0) ms mean 6.6, [0, 300) ms max 12.9 AC 4500 Hz: mean rate 7.7 to 16.6 spikes/s; [-500, 0) ms mean 8.4, [0, 300) ms max 16.6
The rates change after the trial start and not 1 s later (the table above gives both), so the trial start
times share the clock of the spike times and mark the stimulus, as the dataset's description says
("intervals contains each tFUS trial start and end timestamps"). That start_time is the onset of
sonication itself, and not a trigger at a fixed delay, is consistent with the data (the firing changes
within the first 67 ms window at the high PRFs) but not shown. The intervals between trial starts match
the paper's 2.5 s with 10% jitter. In every session with files from two areas, both files carry the same
trial start times, so the two areas were recorded together on one clock (the description says "Each
recording contains two areas"), though their session_start_time values differ (by the minutes printed
above), so that field does not give the start of the recording.
fig4-unit-counts and fig6-unit-counts¶
In the 3000 Hz PRF (60% DC) sessions 312 RSUs and 54 FSUs were recorded from S1 and 233 RSUs and 17 FSUs from POm; in the 30 Hz PRF (0.6% DC) sessions 148 RSUs and 61 FSUs from S1 and 189 RSUs and 43 FSUs from POm.
The PRF series (PD 200 us) comprises 1050 RSUs and 247 FSUs from S1 and 913 RSUs and 106 FSUs from POm.
Criteria (claims file): units with more than 2000 spikes in the PD 200 us sessions of the S1 and POm rats, typed by k-means. Reproduced when every RSU and FSU count is within its tolerance (10%) under one reading; partially reproduced when the total per area (and condition, for Figure 4) is within 10% but the split is not; not reproduced otherwise.
The split needs the cell types, which the dataset does not hold (data check above). The totals do not: the k-means of the header puts every unit in one of its two clusters, so RSU + FSU is the number of units with more than 2000 spikes under every reading K1, K2, P1, P2. If every count of an area is within 10% of the paper's, so is their sum; so an area total more than 10% from the paper's RSU + FSU rules out both the reproduced and the partially reproduced outcome, whatever the split.
Which sessions count is not fixed beyond "the sessions of the S1 and POm rats". The totals are shown for
every selection the identifiers allow: S1 over all S1 recordings, over the rats with POm recordings, and over
the S1 recordings whose session also has a POm file; each also without the recordings whose identifier has a
suffix other than the area pair (NEAR, NEW, 2, POM_S1BF), which may be repeated or moved placements;
and S1 with the files labelled S1_R and S1BF added. POm is counted over all POm recordings and without
POM_S1BF. Only recordings with tFUS aimed at S1 and PD 200 us are used.
def n_above(row, min_spikes=tfus.MIN_SPIKES):
return int(sum(n > min_spikes for n in json.loads(row)))
k["n_good"] = k.n_spikes.map(n_above)
k["suffix"] = k.identifier.str.split("_").str[5:].str.join("_")
pd200 = k[(k.stim == "S1") & (k.pd_us == 200)]
pom_ids = set(pd200[pd200.area == "POM"].identifier)
pom_rats = set(pd200[pd200.area == "POM"].subject_id)
plain = pd200.suffix.isin(["", "POM_S1", "POM _S1"]) # no suffix, or one that names only the area pair
s1a = pd200.area == "S1"
selections = {
"S1, all S1 recordings": s1a,
"S1, rats with POm recordings": s1a & pd200.subject_id.isin(pom_rats),
"S1, sessions with a POm file": s1a & pd200.identifier.isin(pom_ids),
"S1, all, no extra suffix": s1a & plain,
"S1, rats with POm, no extra suffix": s1a & pd200.subject_id.isin(pom_rats) & plain,
"S1, sessions with POm, no extra suffix": s1a & pd200.identifier.isin(pom_ids) & plain,
"S1, all, with S1_R and S1BF files": pd200.area.isin(["S1", "S1_R", "S1BF"]),
"POm, all POm recordings": pd200.area == "POM",
"POm, no extra suffix": (pd200.area == "POM") & plain,
}
paper = {"S1": {3000: 312 + 54, 30: 148 + 61, "all PRFs": 1050 + 247},
"POm": {3000: 233 + 17, 30: 189 + 43, "all PRFs": 913 + 106}}
rows = []
for name, m in selections.items():
g = pd200[m]
r = dict(selection=name)
for prf in [30, 300, 1500, 3000, 4500]:
r[prf] = int(g[g.prf == prf].n_good.sum())
r["all PRFs"] = sum(r[p] for p in [30, 300, 1500, 3000, 4500])
rows.append(r)
tot = pd.DataFrame(rows).set_index("selection")
print("Units with more than 2000 spikes, tFUS at S1, PD 200 us, by PRF (Hz):")
print(tot.to_string())
rows = []
for name in tot.index:
area = name.split(",")[0]
for cond in [3000, 30, "all PRFs"]:
n, ref = int(tot.loc[name, cond]), paper[area][cond]
rows.append(dict(selection=name, condition=cond, units=n, paper_rsu_plus_fsu=ref,
ratio=round(n / ref, 2), within_10pct=abs(n - ref) <= 0.1 * ref))
judged = pd.DataFrame(rows)
print("\nArea totals against the paper's RSU + FSU (Fig. 4 legend for 3000 and 30 Hz, Fig. 6 legend for all PRFs):")
print(judged.to_string(index=False))
Units with more than 2000 spikes, tFUS at S1, PD 200 us, by PRF (Hz):
30 300 1500 3000 4500 all PRFs
selection
S1, all S1 recordings 530 388 493 572 431 2414
S1, rats with POm recordings 390 293 327 383 333 1726
S1, sessions with a POm file 213 141 174 234 224 986
S1, all, no extra suffix 530 388 440 509 431 2298
S1, rats with POm, no extra suffix 390 293 312 320 333 1648
S1, sessions with POm, no extra suffix 213 141 174 234 224 986
S1, all, with S1_R and S1BF files 624 388 500 591 431 2534
POm, all POm recordings 315 267 246 352 290 1470
POm, no extra suffix 315 267 236 339 290 1447
Area totals against the paper's RSU + FSU (Fig. 4 legend for 3000 and 30 Hz, Fig. 6 legend for all PRFs):
selection condition units paper_rsu_plus_fsu ratio within_10pct
S1, all S1 recordings 3000 572 366 1.56 False
S1, all S1 recordings 30 530 209 2.54 False
S1, all S1 recordings all PRFs 2414 1297 1.86 False
S1, rats with POm recordings 3000 383 366 1.05 True
S1, rats with POm recordings 30 390 209 1.87 False
S1, rats with POm recordings all PRFs 1726 1297 1.33 False
S1, sessions with a POm file 3000 234 366 0.64 False
S1, sessions with a POm file 30 213 209 1.02 True
S1, sessions with a POm file all PRFs 986 1297 0.76 False
S1, all, no extra suffix 3000 509 366 1.39 False
S1, all, no extra suffix 30 530 209 2.54 False
S1, all, no extra suffix all PRFs 2298 1297 1.77 False
S1, rats with POm, no extra suffix 3000 320 366 0.87 False
S1, rats with POm, no extra suffix 30 390 209 1.87 False
S1, rats with POm, no extra suffix all PRFs 1648 1297 1.27 False
S1, sessions with POm, no extra suffix 3000 234 366 0.64 False
S1, sessions with POm, no extra suffix 30 213 209 1.02 True
S1, sessions with POm, no extra suffix all PRFs 986 1297 0.76 False
S1, all, with S1_R and S1BF files 3000 591 366 1.61 False
S1, all, with S1_R and S1BF files 30 624 209 2.99 False
S1, all, with S1_R and S1BF files all PRFs 2534 1297 1.95 False
POm, all POm recordings 3000 352 250 1.41 False
POm, all POm recordings 30 315 232 1.36 False
POm, all POm recordings all PRFs 1470 1019 1.44 False
POm, no extra suffix 3000 339 250 1.36 False
POm, no extra suffix 30 315 232 1.36 False
POm, no extra suffix all PRFs 1447 1019 1.42 False
# By rat and identifier suffix, for the selections above (units with more than 2000 spikes)
for area in ["POM", "S1"]:
g = pd200[pd200.area == area]
print(f"\n{area}: units by rat, identifier suffix and PRF (Hz); recordings in brackets")
n = g.pivot_table(index=["subject_id", "suffix"], columns="prf", values="n_good", aggfunc="sum", margins=True)
c = g.pivot_table(index=["subject_id", "suffix"], columns="prf", values="n_good", aggfunc="count", margins=True)
print((n.fillna(0).astype(int).astype(str) + " (" + c.fillna(0).astype(int).astype(str) + ")").to_string())
print("\nFiles of tFUS-at-S1 recordings at PD 200 us whose area is not S1, POM or AC:")
print(pd200[~pd200.area.isin(["S1", "POM", "AC"])][["subject_id", "identifier", "area", "prf", "n_good"]].to_string(index=False))
POM: units by rat, identifier suffix and PRF (Hz); recordings in brackets
prf 30.0 300.0 1500.0 3000.0 4500.0 All
subject_id suffix
BH495 POM_S1 57 (1) 51 (1) 0 (0) 59 (1) 54 (1) 221 (4)
BH499 POM_S1 18 (1) 18 (1) 10 (1) 18 (1) 16 (1) 80 (5)
POM_S1BF 0 (0) 0 (0) 10 (1) 13 (1) 0 (0) 23 (2)
BH501 POM_S1 32 (1) 34 (1) 47 (1) 28 (1) 33 (1) 174 (5)
BH510 POM_S1 24 (1) 33 (1) 15 (1) 23 (1) 24 (1) 119 (5)
BH511 POM_S1 60 (1) 64 (1) 61 (1) 60 (1) 60 (1) 305 (5)
BH513 POM_S1 14 (1) 6 (1) 17 (1) 65 (1) 11 (1) 113 (5)
BH528 POM_S1 58 (1) 61 (1) 58 (1) 59 (1) 46 (1) 282 (5)
BH551 20 (1) 0 (0) 1 (1) 9 (1) 16 (1) 46 (4)
BH553 27 (1) 0 (0) 13 (1) 12 (1) 17 (1) 69 (4)
BH558 5 (1) 0 (0) 14 (1) 6 (1) 13 (1) 38 (4)
All 315 (10) 267 (7) 246 (10) 352 (11) 290 (10) 1470 (48)
S1: units by rat, identifier suffix and PRF (Hz); recordings in brackets
prf 30.0 300.0 1500.0 3000.0 4500.0 All
subject_id suffix
BH494 40 (1) 37 (1) 34 (1) 49 (1) 37 (1) 197 (5)
NEAR 0 (0) 0 (0) 38 (1) 0 (0) 0 (0) 38 (1)
BH495 17 (1) 11 (1) 2 (1) 15 (1) 11 (1) 56 (5)
POM_S1 12 (1) 6 (1) 12 (1) 13 (1) 17 (1) 60 (5)
BH499 6 (1) 10 (1) 5 (1) 16 (1) 0 (0) 37 (4)
POM_S1 24 (1) 18 (1) 16 (1) 23 (1) 29 (1) 110 (5)
BH500 16 (1) 15 (1) 0 (0) 0 (0) 18 (1) 49 (3)
BH501 12 (1) 14 (1) 9 (1) 0 (0) 10 (1) 45 (4)
NEW 0 (0) 0 (0) 15 (1) 7 (1) 0 (0) 22 (2)
POM_S1 9 (1) 7 (1) 12 (1) 6 (1) 2 (1) 36 (5)
BH510 13 (1) 7 (1) 20 (1) 21 (1) 11 (1) 72 (5)
POM_S1 14 (1) 17 (1) 23 (1) 19 (1) 16 (1) 89 (5)
BH511 27 (1) 10 (1) 23 (1) 23 (1) 4 (1) 87 (5)
POM _S1 53 (1) 0 (0) 0 (0) 0 (0) 0 (0) 53 (1)
POM_S1 0 (0) 56 (1) 60 (1) 55 (1) 57 (1) 228 (4)
BH513 39 (1) 67 (1) 48 (1) 0 (0) 57 (1) 211 (4)
2 0 (0) 0 (0) 0 (0) 56 (1) 0 (0) 56 (1)
POM_S1 35 (1) 20 (1) 20 (1) 24 (1) 9 (1) 108 (5)
BH528 10 (1) 33 (1) 19 (1) 11 (1) 16 (1) 89 (5)
POM_S1 12 (1) 17 (1) 18 (1) 13 (1) 8 (1) 68 (5)
BH530 0 (0) 31 (1) 8 (1) 0 (0) 8 (1) 47 (3)
BH535 49 (1) 0 (0) 42 (1) 41 (1) 0 (0) 132 (3)
BH537 0 (0) 4 (1) 0 (0) 0 (0) 2 (1) 6 (2)
BH538 0 (0) 8 (1) 0 (0) 8 (1) 7 (1) 23 (3)
BH539 9 (1) 0 (0) 44 (1) 32 (1) 0 (0) 85 (3)
BH547 26 (1) 0 (0) 0 (0) 0 (0) 26 (1) 52 (2)
BH551 38 (1) 0 (0) 0 (0) 19 (1) 41 (1) 98 (3)
BH553 49 (1) 0 (0) 8 (1) 46 (1) 16 (1) 119 (4)
BH558 20 (1) 0 (0) 17 (1) 16 (1) 29 (1) 82 (4)
BH594 0 (0) 0 (0) 0 (0) 59 (1) 0 (0) 59 (1)
All 530 (22) 388 (19) 493 (22) 572 (22) 431 (22) 2414 (107)
Files of tFUS-at-S1 recordings at PD 200 us whose area is not S1, POM or AC:
subject_id identifier area prf n_good
BH494 BH494_1500_200_67_5 M1 1500.0 3
BH494 BH494_1500_200_67_5_NEAR M1 1500.0 48
BH494 BH494_3000_200_67_5 M1 3000.0 12
BH494 BH494_4500_200_67_5 M1 4500.0 14
BH494 BH494_30_200_67_5 M1 30.0 3
BH494 BH494_300_200_67_5 M1 300.0 3
BH494 BH494_30_200_67_5 M1_R 30.0 12
BH494 BH494_30_200_67_5 S1_R 30.0 94
BH495 BH495_1500_200_67_5 M1 1500.0 10
BH495 BH495_30_200_67_5 M1 30.0 14
BH495 BH495_4500_200_67_5 M1 4500.0 13
BH495 BH495_300_200_67_5 M1 300.0 17
BH495 BH495_3000_200_67_5 M1 3000.0 17
BH499 BH499_3000_200_67_5 M1 3000.0 15
BH499 BH499_3000_200_67_5_POM_S1BF S1BF 3000.0 12
BH499 BH499_3000_200_67_5_VPM_S1BF VPM 3000.0 30
BH499 BH499_3000_200_67_5_VPM_S1BF S1BF 3000.0 7
BH499 BH499_1500_200_67_5_POM_S1BF S1BF 1500.0 7
BH499 BH499_30_200_67_5 M1 30.0 12
BH499 BH499_300_200_67_5 M1 300.0 9
BH499 BH499_1500_200_67_5 M1 1500.0 9
BH500 BH500_3000_200_67_5 3000.0 5
BH500 BH500_1500_200_201_5 1500.0 16
BH500 BH500_1500_200_67_5 1500.0 12
BH500 BH500_4500_200_67_5 M1 4500.0 17
BH500 BH500_30_200_67_5 M1 30.0 12
BH500 BH500_300_200_67_5 M1 300.0 13
BH501 BH501_1500_200_67_5 M1 1500.0 25
BH501 BH501_1500_200_67_5_NEW M1 1500.0 55
BH501 BH501_30_200_67_5 M1 30.0 50
BH501 BH501_300_200_67_5 M1 300.0 51
BH501 BH501_3000_200_67_5_NEW M1 3000.0 49
BH501 BH501_4500_200_67_5 M1 4500.0 52
BH510 BH510_1500_200_67_5 M1 1500.0 16
BH510 BH510_3000_200_67_5 M1 3000.0 15
BH510 BH510_30_200_67_5 M1 30.0 6
BH510 BH510_300_200_67_5 M1 300.0 5
BH510 BH510_4500_200_67_5 M1 4500.0 5
BH511 BH511_1500_200_67_5 M1 1500.0 31
BH511 BH511_30_200_67_5 M1 30.0 33
BH511 BH511_300_200_67_5 M1 300.0 35
BH511 BH511_3000_200_67_5 M1 3000.0 49
BH511 BH511_4500_200_67_5 M1 4500.0 30
BH513 BH513_30_200_67_5 M1 30.0 9
BH513 BH513_300_200_67_5 M1 300.0 10
BH513 BH513_1500_200_67_5 M1 1500.0 11
BH513 BH513_3000_200_67_5_2 M1 3000.0 12
BH513 BH513_4500_200_67_5 M1 4500.0 6
BH528 BH528_1500_200_67_5 M1 1500.0 36
BH528 BH528_3000_200_67_5 M1 3000.0 23
BH528 BH528_300_200_67_5 M1 300.0 28
BH528 BH528_4500_200_67_5 M1 4500.0 28
BH528 BH528_30_200_67_5 M1 30.0 23
# How the area totals change with the spike-count rule (the paper's is more than 2000 spikes)
rows = []
for m in [1000, 2000, 4000, 8000]:
cnt = pd200.n_spikes.map(lambda r: n_above(r, m))
for name in ["S1, all S1 recordings", "S1, rats with POm recordings", "S1, sessions with a POm file",
"POm, all POm recordings"]:
g = pd200[selections[name]].assign(n=cnt[selections[name]])
area = name.split(",")[0]
r = dict(min_spikes=m, selection=name)
for cond in [3000, 30]:
r[cond] = int(g[g.prf == cond].n.sum())
r["all PRFs"] = int(g.n.sum())
r["all within 10%"] = all(abs(r[c] - paper[area][c]) <= 0.1 * paper[area][c] for c in [3000, 30, "all PRFs"])
rows.append(r)
print("Units above each spike-count rule, tFUS at S1, PD 200 us (paper RSU + FSU: S1 366, 209, 1297; POm 250, 232, 1019):")
print(pd.DataFrame(rows).to_string(index=False))
# POm: which pairs of rats, left out, bring all three POm totals within 10%
from itertools import combinations
pom = pd200[pd200.area == "POM"]
rows = []
for pair in combinations(sorted(pom.subject_id.unique()), 2):
rest = pom[~pom.subject_id.isin(pair)]
r = {3000: int(rest[rest.prf == 3000].n_good.sum()), 30: int(rest[rest.prf == 30].n_good.sum()),
"all PRFs": int(rest.n_good.sum())}
if all(abs(r[c] - paper["POm"][c]) <= 0.1 * paper["POm"][c] for c in r):
rows.append({"left_out": " ".join(pair), **r})
print(f"\nPairs of the {pom.subject_id.nunique()} POm rats whose removal brings the POm totals at 3000 Hz, 30 Hz and "
f"over the five PRFs within 10% of 250, 232 and 1019: {len(rows)} of "
f"{len(list(combinations(pom.subject_id.unique(), 2)))}")
print(pd.DataFrame(rows).to_string(index=False))
Units above each spike-count rule, tFUS at S1, PD 200 us (paper RSU + FSU: S1 366, 209, 1297; POm 250, 232, 1019):
min_spikes selection 3000 30 all PRFs all within 10%
1000 S1, all S1 recordings 639 614 2667 False
1000 S1, rats with POm recordings 419 450 1885 False
1000 S1, sessions with a POm file 258 233 1071 False
1000 POm, all POm recordings 376 344 1599 False
2000 S1, all S1 recordings 572 530 2414 False
2000 S1, rats with POm recordings 383 390 1726 False
2000 S1, sessions with a POm file 234 213 986 False
2000 POm, all POm recordings 352 315 1470 False
4000 S1, all S1 recordings 503 456 2110 False
4000 S1, rats with POm recordings 347 343 1537 False
4000 S1, sessions with a POm file 204 182 857 False
4000 POm, all POm recordings 297 275 1257 False
8000 S1, all S1 recordings 446 411 1860 False
8000 S1, rats with POm recordings 314 308 1371 False
8000 S1, sessions with a POm file 174 162 730 False
8000 POm, all POm recordings 242 230 1031 True
Pairs of the 10 POm rats whose removal brings the POm totals at 3000 Hz, 30 Hz and over the five PRFs within 10% of 250, 232 and 1019: 9 of 45
left_out 3000 30 all PRFs
BH495 BH501 265 226 1075
BH499 BH511 261 237 1062
BH499 BH528 262 239 1085
BH501 BH511 264 223 991
BH501 BH528 265 225 1014
BH510 BH511 269 231 1046
BH510 BH528 270 233 1069
BH511 BH513 227 241 1052
BH513 BH528 228 243 1075
Result: not reproduced (both claims). RSU + FSU equals the number of units with more than 2000 spikes
under every typing reading, and the POm totals are outside 10% of the paper's under both POm selections:
at 3000 Hz 352 (339 without the POM_S1BF recording) against 250, at 30 Hz 315 against 232, and over the
five PRFs 1470 (1447) against 1019, ratios 1.36 to 1.44. So no RSU/FSU split could bring the eight (or four)
counts within tolerance, and the POm totals are not within 10% either, which rules out the partial outcome.
The S1 totals depend on the selection (at 3000 Hz from 234 to 591 against 366, at 30 Hz from 213 to 624
against 209, over the five PRFs from 986 to 2534 against 1297). One selection matches at 3000 Hz (the rats
with POm recordings, 383) and another at 30 Hz (the sessions with a POm file, 213), but none matches both
conditions or the five-PRF total. The POm units come from all 10 POm rats (6 to 65 units each at 3000 Hz,
not counting the POM_S1BF file), so no single rat or suffixed recording explains the excess. Leaving out
two of the 10 POm rats does bring all three POm totals within 10% for 9 of the 45 pairs (for example BH511
and BH513: 227, 241 and 1052), but nothing in the data or the paper marks any rat for exclusion, and the
paper reports more POm rats (16), not fewer. The totals also depend on the unit rule: with more than 8000
spikes in place of 2000, the POm totals are 242, 230 and 1031, all within 10%, while no S1 selection fits all
three at any of the four rules shown. The dataset holds 10 POm rats where the paper reports 16, and more POm
units than the paper, so the units the paper counted are not the units with more than 2000 spikes over the
whole of these recordings; the paper's unit rule or spike counting may differ from that reading. The outcome
follows the criterion as written. The RSU/FSU split itself could not be tested.
The STTC files¶
The 370 files under sub-STTC-* hold STTC values between the units of two areas. They are listed here
because the claims file names them as a variant for Figure 10; they carry no unit identities or cell
types.
sttc = inv[inv.kind == "STTC"]
print(sttc.notes.value_counts().to_string())
import h5py, remfile
f = h5py.File(remfile.File(tfus.asset_url(sttc.asset_id.iloc[0])), "r")
print("\n", sttc.path.iloc[0], f["identifier"][()].decode(), f["general/notes"][()].decode())
print(f["units"].attrs["description"])
for c in ["STTC", "STTCAtoB", "STTCBtoA", "STTCA", "STTCB"]:
print(c, f["units"][c].shape, "|", f["units"][c].attrs["description"])
notes STTC results of:AC vsS1; stim at:S1 126 STTC results of:POM vsS1; stim at:S1 89 STTC results of:M1 vsS1; stim at:S1 46 STTC results of:M1 vsS1; stim at:M1 42 STTC results of:POM vsS1; stim at:M1 30 STTC results of:POM vsS1; stim at:POM 18 STTC results of:M1 vsS1; stim at:POM 13 STTC results of:POM vsS1BF; stim at:S1 2 STTC results of:S1BF vsVPM; stim at:S1 2 STTC results of: vs; stim at:S1 1 STTC results of:AC_1 vsAC_2; stim at:S1 1
sub-STTC-BH494/sub-STTC-BH494_ses-20240118T231600_ecephys.nwb BH494_1500_200_67_5_MOTOR STTC results of:M1 vsS1; stim at:M1
STTC results
{time_range = [-0.067 0;0 0.067;0.067 0.134; 0.134 0.201; 0.201 0.268; 0.268 0.335; 0.335 0.402; 0.402 0.469; 0.469 0.536; 0.536 0.603;];
dt= [0.005,0.01,0.015,0.020,0.025,0.030,0.040,0.050,0.067]},
time_range: column 1 is start and column 2 is the end for each caculation time window;
dt is the quantification of the proportion of total recording time that lies [-dt, +dt] of each spike in each train
and the proportion of spikes in spiketrain_i that lies [-dt, +dt] of any spike in spiketrain_j
STTC (10, 499, 20, 9) | STTC without direction, size: number(dt)*number(spike_target1)*number(spike_target2)*trial_number*number(time_range)
STTCAtoB (10, 499, 20, 9) | STTC from A to B, size: number(dt)*number(spike_target1)*number(spike_target2)*trial_number*number(time_range), A ahead B
STTCBtoA (10, 499, 20, 9) | STTC from B to A, size: number(dt)*number(spike_target1)*number(spike_target2)*trial_number*number(time_range), B ahead A STTCA (10, 499, 9, 9) | STTC within A, size: number(dt)*number(spike_target1)*number(spike_target1)*trial_number*number(time_range) STTCB (10, 499, 20, 20) | STTC within B, size: number(dt)*number(spike_target2)*number(spike_target2)*trial_number*number(time_range)
The STTC arrays do not have the shape their descriptions give: the descriptions list a dimension for the nine values of dt, while the arrays have the 10 time windows, the trials and the units of the two areas, and no dt dimension. Which dt the values were computed with is not stated, and no unit is typed, so the stored values cannot stand in for the STTC between S1 RSUs and POm RSUs of Figure 10.
Untyped responses in S1 and POm (a variant; decides no outcome)¶
The claims on S1 and POm need the RSU/FSU split. For orientation only, the same measure on all units of each area, without typing (reading N1, comparison B1, Bonferroni over the five PRFs):
rows = []
for area in ["S1", "POM"]:
for prf, g in units[units.area == area].groupby("prf"):
base = tfus.baseline_sample(list(g["true"]), "N1", "B1")
r = dict(area=area, prf=int(prf), n=len(g))
for w in tfus.WINDOWS:
x = np.array([tfus.normalized_rates(c, w, "N1") for c in g["true"]])
p_up, _ = tfus.vs_baseline(x, base, 5, "greater")
p_dn, _ = tfus.vs_baseline(x, base, 5, "less")
r[w] = f"{np.nanmean(x):.0f}{'+' if p_up < 0.01 else ('-' if p_dn < 0.01 else '')}"
rows.append(r)
print("Mean normalized rate (%) of all units; + above, - below baseline at corrected p < 0.01:")
print(pd.DataFrame(rows).to_string(index=False))
Mean normalized rate (%) of all units; + above, - below baseline at corrected p < 0.01: area prf n 0-67 67-134 134-201 201-268 S1 30 530 100 96- 101+ 106+ S1 300 388 107+ 102+ 100 106 S1 1500 493 127+ 105 108+ 111+ S1 3000 572 114+ 102+ 119+ 178+ S1 4500 431 102+ 103 116+ 149+ POM 30 315 97- 100 111+ 110+ POM 300 267 104+ 106+ 115+ 111+ POM 1500 246 108+ 120+ 157+ 131+ POM 3000 352 105+ 100 113+ 122+ POM 4500 290 100 103+ 117+ 129+
fig11b-ac-responses¶
AC neurons in the S1-AC rats show time-locked responses to 1500 and 3000 Hz PRF and postsonication responses to 30, 1500, 3000 and 4500 Hz, both PRF dependent, least at 300 Hz.
Criterion (claims file): the header's spiking method on all AC units with more than 2000 spikes in the PRF series with 200 us PD; windows 0 to 67 and 201 to 268 ms (the 67 to 134 and 134 to 201 ms bins shown). Reproduced when, under every method reading, both Kruskal-Wallis tests have p < 0.01, exactly 1500 and 3000 Hz are above baseline during sonication, exactly 30, 1500, 3000 and 4500 Hz are above baseline postsonication, and 300 Hz has the lowest mean normalized rate in at least one window; partially reproduced when the sets differ by one PRF in one window or the 300 Hz part fails. Expected: 453 AC neurons over the five PRFs.
The method readings (claims file header): N1 or N2 for the normalized rate; B1 or B2 for the baseline sample; the Bonferroni factor is the five PRF levels of the panel (C1; C2 applies only to a panel with a single comparison, so it is the same here). B3, a paired signed-rank test, is shown as a variant.
ac = units[units.area == "AC"].copy()
PRFS = [30.0, 300.0, 1500.0, 3000.0, 4500.0]
print("AC units with more than 2000 spikes, by rat and PRF:")
print(pd.crosstab(ac.subject, ac.prf, margins=True).to_string())
print(f"\nTotal {len(ac)} (expected 453 +/- 45.3)")
allac = s1[s1.area == "AC"]
print("AC recordings used:", len(allac), "; units in them before the 2000-spike rule:",
int(sum(r["n_units"] for r in allac.res)))
# Every AC-labelled recording with tFUS at S1, by rat and PD (us): which rats have PD 200 us sessions
acall = k[k.area.str.startswith("AC") & (k.stim == "S1")]
print("\nAC-labelled recordings (any PD) by rat and PD; units with more than 2000 spikes in brackets:")
n = acall.pivot_table(index=["subject_id", "area"], columns="pd_us", values="n_good", aggfunc="sum")
c = acall.pivot_table(index=["subject_id", "area"], columns="pd_us", values="n_good", aggfunc="count")
print((n.fillna(0).astype(int).astype(str) + " (" + c.fillna(0).astype(int).astype(str) + ")").to_string())
print("PRF (Hz) of the AC_1 and AC_2 files:", acall[acall.area.isin(["AC_1", "AC_2"])].prf.tolist())
AC units with more than 2000 spikes, by rat and PRF:
prf 30.0 300.0 1500.0 3000.0 4500.0 All
subject
BH530 0 10 5 0 10 25
BH535 45 0 49 47 0 141
BH537 0 25 24 0 20 69
BH538 0 54 0 61 45 160
BH539 11 0 60 59 0 130
BH547 56 0 0 0 52 108
All 112 89 138 167 127 633
Total 633 (expected 453 +/- 45.3)
AC recordings used: 17 ; units in them before the 2000-spike rule: 667
AC-labelled recordings (any PD) by rat and PD; units with more than 2000 spikes in brackets:
pd_us 1.4 2.0 4.0 14.0 20.0 40.0 67.0 100.0 134.0 200.0 300.0 400.0 600.0 1000.0 2000.0 3000.0 10000.0 20000.0 30000.0
subject_id area
BH530 AC 0 (0) 0 (0) 0 (0) 8 (1) 12 (1) 13 (1) 14 (1) 9 (1) 0 (0) 25 (3) 8 (1) 0 (0) 8 (1) 10 (1) 26 (2) 11 (1) 9 (1) 0 (0) 13 (1)
BH532 AC 11 (1) 17 (1) 11 (1) 19 (1) 24 (2) 6 (1) 10 (1) 53 (5) 8 (1) 0 (0) 7 (1) 12 (1) 17 (1) 13 (1) 22 (2) 18 (1) 7 (1) 15 (1) 18 (1)
BH535 AC 48 (1) 55 (1) 59 (1) 57 (1) 102 (2) 59 (1) 50 (1) 64 (1) 61 (1) 141 (3) 58 (1) 60 (1) 58 (1) 49 (1) 103 (2) 55 (1) 63 (1) 57 (1) 58 (1)
BH537 AC 22 (1) 24 (1) 30 (1) 20 (1) 23 (1) 11 (1) 19 (1) 27 (1) 19 (1) 69 (3) 22 (1) 22 (1) 24 (1) 23 (1) 60 (2) 14 (1) 21 (1) 20 (1) 23 (1)
BH538 AC 67 (1) 60 (1) 60 (1) 49 (1) 116 (2) 0 (0) 50 (1) 57 (1) 25 (1) 160 (3) 58 (1) 64 (1) 58 (1) 55 (1) 116 (2) 50 (1) 62 (1) 62 (1) 62 (1)
AC_1 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 30 (1) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)
AC_2 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 31 (1) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)
BH539 AC 60 (1) 43 (1) 60 (1) 62 (1) 119 (2) 63 (1) 61 (1) 44 (1) 62 (1) 130 (3) 47 (1) 62 (1) 59 (1) 59 (1) 100 (2) 58 (1) 21 (1) 44 (1) 21 (1)
BH547 AC 66 (1) 51 (1) 58 (1) 67 (1) 83 (2) 64 (1) 63 (1) 54 (1) 57 (1) 108 (2) 37 (1) 59 (1) 61 (1) 48 (1) 99 (2) 53 (1) 55 (1) 56 (1) 51 (1)
PRF (Hz) of the AC_1 and AC_2 files: [1500.0, 1500.0]
def evaluate(frame, nreading, breading, min_label=""):
out, sig = [], {}
for w in tfus.WINDOWS:
groups, line = [], {}
for p in PRFS:
cs = list(frame[frame.prf == p]["true"])
x = np.array([tfus.normalized_rates(c, w, nreading) for c in cs])
groups.append(x)
if breading == "B3":
from scipy import stats as st
win = np.array([c[w].mean() for c in cs]); bas = np.array([c["base"].mean() for c in cs])
pr = st.wilcoxon(win, bas, alternative="greater").pvalue
pc = min(1, pr * 5)
else:
base = tfus.baseline_sample(cs, nreading, breading)
pc, pr = tfus.vs_baseline(x, base, 5, "greater")
sem = np.nanstd(x, ddof=1) / np.sqrt(np.isfinite(x).sum())
out.append(dict(N=nreading, B=breading, window=w, prf=int(p), n=int(np.isfinite(x).sum()),
mean=round(np.nanmean(x), 1), sem=round(sem, 1), median=round(np.nanmedian(x), 1),
p_corrected=float(f"{pc:.3g}"), above=pc < 0.01))
kw = tfus.kruskal_levels(groups)
for o in out[-5:]:
o["kruskal_p"] = float(f"{kw:.3g}")
return pd.DataFrame(out)
def judge(t):
def sset(w):
return {r.prf for r in t[(t.window == w) & t.above].itertuples()}
during, post = sset("0-67"), sset("201-268")
kw_ok = (t[t.window == "0-67"].kruskal_p.iloc[0] < 0.01) and (t[t.window == "201-268"].kruskal_p.iloc[0] < 0.01)
d_err = len(during ^ {1500, 3000})
p_err = len(post ^ {30, 1500, 3000, 4500})
low = [t[t.window == w].sort_values("mean").prf.iloc[0] for w in ["0-67", "201-268"]]
low300 = 300 in low
if kw_ok and d_err == 0 and p_err == 0 and low300:
verdict = "reproduced"
elif kw_ok and ((d_err + p_err == 1 and low300) or (d_err == 0 and p_err == 0 and not low300)):
verdict = "partially reproduced"
else:
verdict = "not reproduced"
return dict(kruskal_both=kw_ok, during=sorted(during), post=sorted(post), lowest_0_67=low[0],
lowest_201_268=low[1], verdict=verdict)
tables, verdicts = {}, []
for N in ["N1", "N2"]:
for B in ["B1", "B2"]:
t = evaluate(ac, N, B)
tables[(N, B)] = t
verdicts.append(dict(N=N, B=B, **judge(t)))
t_b3 = evaluate(ac, "N1", "B3")
verdicts.append(dict(N="raw rates", B="B3 (variant)", **judge(t_b3)))
print("Outcome under each reading (sets are the PRFs above baseline at corrected p < 0.01):")
print(pd.DataFrame(verdicts).to_string(index=False))
read = [v for v in verdicts if v["B"] in ("B1", "B2")]
outs = {v["verdict"] for v in read}
combined = outs.pop() if len(outs) == 1 else "blocked (methods_underspecified): the readings give different outcomes"
print("\nVerdicts of the four readings N1/N2 x B1/B2:", [f"{v['N']}/{v['B']}: {v['verdict']}" for v in read])
print("Outcome under the claims file header's rule (the one all readings give; blocked when they differ):", combined)
n1 = {v["verdict"] for v in read if v["N"] == "N1"}
print("Verdicts of N1 with B1 and B2 only:", sorted(n1))
Outcome under each reading (sets are the PRFs above baseline at corrected p < 0.01):
N B kruskal_both during post lowest_0_67 lowest_201_268 verdict
N1 B1 True [1500, 3000, 4500] [30, 1500, 3000, 4500] 300 300 partially reproduced
N1 B2 True [300, 1500, 3000, 4500] [30, 1500, 3000, 4500] 300 300 not reproduced
N2 B1 False [] [] 3000 300 not reproduced
N2 B2 False [1500] [4500] 3000 300 not reproduced
raw rates B3 (variant) True [1500, 3000, 4500] [30, 1500, 3000, 4500] 300 300 partially reproduced
Verdicts of the four readings N1/N2 x B1/B2: ['N1/B1: partially reproduced', 'N1/B2: not reproduced', 'N2/B1: not reproduced', 'N2/B2: not reproduced']
Outcome under the claims file header's rule (the one all readings give; blocked when they differ): blocked (methods_underspecified): the readings give different outcomes
Verdicts of N1 with B1 and B2 only: ['not reproduced', 'partially reproduced']
for (N, B), t in tables.items():
print(f"\nReading {N}, {B}: mean, SEM and median normalized rate (%), corrected p against baseline, Kruskal-Wallis p")
print(t.pivot_table(index="prf", columns="window", values=["mean", "sem", "p_corrected"], sort=False)
.reindex(columns=list(tfus.WINDOWS), level=1).to_string())
print("Kruskal-Wallis p by window:", t.groupby("window", sort=False).kruskal_p.first().to_dict())
print("\nVariant B3 (paired signed-rank of raw rates against [-67, 0) ms), corrected p:")
print(t_b3.pivot_table(index="prf", columns="window", values="p_corrected", sort=False).to_string())
Reading N1, B1: mean, SEM and median normalized rate (%), corrected p against baseline, Kruskal-Wallis p
mean sem p_corrected
window 0-67 67-134 134-201 201-268 0-67 67-134 134-201 201-268 0-67 67-134 134-201 201-268
prf
30 105.9 115.9 138.0 176.4 1.4 1.4 3.3 7.5 1.890000e-02 1.330000e-24 2.380000e-26 3.260000e-26
300 99.7 94.9 97.5 94.3 1.4 1.1 1.3 1.1 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
1500 140.3 192.8 149.2 147.9 8.8 24.6 9.1 10.0 2.460000e-08 4.810000e-22 9.070000e-20 7.660000e-15
3000 107.1 108.0 223.7 310.5 1.2 2.0 17.3 28.0 9.720000e-11 1.050000e-01 1.230000e-23 9.350000e-40
4500 108.0 144.2 228.8 210.8 1.5 6.2 9.9 7.4 7.050000e-05 4.860000e-24 2.250000e-42 3.940000e-45
Kruskal-Wallis p by window: {'0-67': 7.83e-06, '67-134': 1.04e-29, '134-201': 4.42e-42, '201-268': 5.16e-49}
Reading N1, B2: mean, SEM and median normalized rate (%), corrected p against baseline, Kruskal-Wallis p
mean sem p_corrected
window 0-67 67-134 134-201 201-268 0-67 67-134 134-201 201-268 0-67 67-134 134-201 201-268
prf
30 105.9 115.9 138.0 176.4 1.4 1.4 3.3 7.5 1.200000e-01 1.490000e-11 4.820000e-18 3.910000e-18
300 99.7 94.9 97.5 94.3 1.4 1.1 1.3 1.1 6.290000e-03 1.000000e+00 3.680000e-01 1.000000e+00
1500 140.3 192.8 149.2 147.9 8.8 24.6 9.1 10.0 6.480000e-12 7.230000e-21 6.440000e-21 1.270000e-17
3000 107.1 108.0 223.7 310.5 1.2 2.0 17.3 28.0 5.180000e-17 1.850000e-09 7.050000e-31 1.580000e-36
4500 108.0 144.2 228.8 210.8 1.5 6.2 9.9 7.4 9.040000e-03 3.140000e-16 4.440000e-38 5.590000e-41
Kruskal-Wallis p by window: {'0-67': 7.83e-06, '67-134': 1.04e-29, '134-201': 4.42e-42, '201-268': 5.16e-49}
Reading N2, B1: mean, SEM and median normalized rate (%), corrected p against baseline, Kruskal-Wallis p
mean sem p_corrected
window 0-67 67-134 134-201 201-268 0-67 67-134 134-201 201-268 0-67 67-134 134-201 201-268
prf
30 63.6 58.1 59.1 64.2 2.4 2.2 2.9 3.4 1.0 1.0 1.0 1.0
300 65.3 58.0 57.1 55.5 2.3 1.9 2.1 1.9 1.0 1.0 1.0 1.0
1500 66.2 65.3 58.8 56.9 2.3 1.9 2.1 2.1 1.0 1.0 1.0 1.0
3000 59.5 51.4 56.8 66.1 1.4 1.6 2.7 3.8 1.0 1.0 1.0 1.0
4500 64.6 66.7 75.5 74.1 1.8 2.4 3.3 3.3 1.0 1.0 1.0 1.0
Kruskal-Wallis p by window: {'0-67': 0.0213, '67-134': 7.2e-07, '134-201': 6.08e-05, '201-268': 0.00366}
Reading N2, B2: mean, SEM and median normalized rate (%), corrected p against baseline, Kruskal-Wallis p
mean sem p_corrected
window 0-67 67-134 134-201 201-268 0-67 67-134 134-201 201-268 0-67 67-134 134-201 201-268
prf
30 63.6 58.1 59.1 64.2 2.4 2.2 2.9 3.4 1.000000 1.00000 1.00000 1.00000
300 65.3 58.0 57.1 55.5 2.3 1.9 2.1 1.9 0.768000 1.00000 1.00000 1.00000
1500 66.2 65.3 58.8 56.9 2.3 1.9 2.1 2.1 0.000179 0.00028 0.27800 0.93800
3000 59.5 51.4 56.8 66.1 1.4 1.6 2.7 3.8 0.341000 1.00000 1.00000 1.00000
4500 64.6 66.7 75.5 74.1 1.8 2.4 3.3 3.3 0.088900 0.30200 0.00258 0.00197
Kruskal-Wallis p by window: {'0-67': 0.0213, '67-134': 7.2e-07, '134-201': 6.08e-05, '201-268': 0.00366}
Variant B3 (paired signed-rank of raw rates against [-67, 0) ms), corrected p:
window 0-67 67-134 134-201 201-268
prf
30 8.840000e-02 6.560000e-15 9.500000e-16 8.120000e-16
300 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
1500 3.760000e-11 5.010000e-16 1.880000e-16 1.930000e-12
3000 6.340000e-12 3.820000e-02 4.750000e-18 1.960000e-24
4500 1.930000e-05 3.490000e-17 4.350000e-22 3.560000e-22
The two normalizations differ in kind. N2 averages per-trial ratios over the trials with a nonzero baseline count; at these firing rates most 67 ms baselines hold no spike, and keeping only the trials with one makes the baseline larger than average, so N2 sits below 100% even where nothing happens. The next cell measures this at the trial starts shifted by 1 s, where there is no stimulus.
rows = []
for N in ["N1", "N2"]:
for w in tfus.WINDOWS:
x = np.array([tfus.normalized_rates(c, w, N) for c in ac["shifted"]])
rows.append(dict(N=N, window=w, mean_at_shifted_onsets=round(np.nanmean(x), 1),
median=round(np.nanmedian(x), 1)))
print("AC units, normalized rate (%) with trial starts shifted by 1 s (no stimulus):")
print(pd.DataFrame(rows).pivot(index="N", columns="window", values="mean_at_shifted_onsets")
[list(tfus.WINDOWS)].to_string())
nz = np.array([np.mean(c["base"] > 0) for c in ac["true"]])
print(f"\nFraction of trials with a nonzero [-67, 0) ms count, over AC units: median {np.median(nz):.2f}, "
f"range {nz.min():.2f} to {nz.max():.2f}")
AC units, normalized rate (%) with trial starts shifted by 1 s (no stimulus): window 0-67 67-134 134-201 201-268 N N1 99.0 94.8 95.4 94.3 N2 61.7 53.3 51.8 49.8 Fraction of trials with a nonzero [-67, 0) ms count, over AC units: median 0.23, range 0.02 to 0.72
# By rat: mean normalized rate (N1) in the two windows of the criterion
rows = []
for (sub, prf), g in ac.groupby(["subject", "prf"]):
r = dict(rat=sub, prf=int(prf), n=len(g))
for w in ["0-67", "201-268"]:
r[w] = round(np.nanmean([tfus.normalized_rates(c, w, "N1") for c in g["true"]]), 1)
rows.append(r)
print("Mean normalized rate (%) by rat, reading N1:")
print(pd.DataFrame(rows).to_string(index=False))
# The 2000-spike rule: how the count and the result change with it
rows = []
for m in [1000, 2000, 4000, 8000]:
sub = unit_table(s1[s1.area == "AC"], min_spikes=m)
j = judge(evaluate(sub, "N1", "B1"))
rows.append(dict(min_spikes=m, n_units=len(sub), during=j["during"], post=j["post"],
lowest=(j["lowest_0_67"], j["lowest_201_268"]), verdict=j["verdict"]))
print("\nReading N1, B1 under other spike-count rules (the paper's is more than 2000):")
print(pd.DataFrame(rows).to_string(index=False))
Mean normalized rate (%) by rat, reading N1: rat prf n 0-67 201-268 BH530 300 10 105.3 96.3 BH530 1500 5 131.4 95.9 BH530 4500 10 107.1 124.0 BH535 30 45 114.3 113.6 BH535 1500 49 112.1 101.4 BH535 3000 47 113.0 113.8 BH537 300 25 104.7 99.7 BH537 1500 24 302.1 333.1 BH537 4500 20 135.3 297.0 BH538 300 54 96.3 91.4 BH538 3000 61 104.9 654.1 BH538 4500 45 106.1 218.2 BH539 30 11 95.7 99.0 BH539 1500 60 99.3 116.2 BH539 3000 59 104.7 112.0 BH547 30 56 101.1 242.1 BH547 4500 52 99.4 187.9
Reading N1, B1 under other spike-count rules (the paper's is more than 2000):
min_spikes n_units during post lowest verdict
1000 667 [1500, 3000, 4500] [30, 1500, 3000, 4500] (300, 300) partially reproduced
2000 633 [1500, 3000, 4500] [30, 1500, 3000, 4500] (300, 300) partially reproduced
4000 445 [1500, 3000, 4500] [30, 1500, 3000, 4500] (30, 300) partially reproduced
8000 320 [1500, 3000, 4500] [30, 1500, 3000, 4500] (30, 300) partially reproduced
# Leave one rat out: the verdict of readings N1/B1 and N1/B2 without each rat in turn
rows = []
for rat in sorted(ac.subject.unique()):
sub = ac[ac.subject != rat]
for B in ["B1", "B2"]:
t = evaluate(sub, "N1", B)
j = judge(t)
m = t.set_index(["window", "prf"])["mean"]
rows.append(dict(left_out=rat, B=B, n_units=len(sub), kw_0_67=t[t.window == "0-67"].kruskal_p.iloc[0],
kw_201_268=t[t.window == "201-268"].kruskal_p.iloc[0], during=j["during"], post=j["post"],
lowest=(j["lowest_0_67"], j["lowest_201_268"]),
mean_4500_0_67=m[("0-67", 4500)], mean_300_0_67=m[("0-67", 300)], verdict=j["verdict"]))
print("Readings N1/B1 and N1/B2 with each rat left out:")
print(pd.DataFrame(rows).to_string(index=False))
Readings N1/B1 and N1/B2 with each rat left out: left_out B n_units kw_0_67 kw_201_268 during post lowest mean_4500_0_67 mean_300_0_67 verdict BH530 B1 608 1.420000e-05 3.170000e-46 [1500, 3000, 4500] [30, 1500, 3000, 4500] (300, 300) 108.1 99.0 partially reproduced BH530 B2 608 1.420000e-05 3.170000e-46 [1500, 3000] [30, 1500, 3000, 4500] (300, 300) 108.1 99.0 reproduced BH535 B1 492 5.130000e-05 4.290000e-44 [1500, 3000, 4500] [30, 1500, 3000, 4500] (300, 300) 108.0 99.7 partially reproduced BH535 B2 492 5.130000e-05 4.290000e-44 [300, 1500, 3000, 4500] [30, 1500, 3000, 4500] (300, 300) 108.0 99.7 not reproduced BH537 B1 564 1.390000e-04 1.660000e-48 [1500, 3000] [30, 1500, 3000, 4500] (300, 300) 102.9 97.7 reproduced BH537 B2 564 1.390000e-04 1.660000e-48 [1500, 3000] [30, 1500, 3000, 4500] (300, 300) 102.9 97.7 reproduced BH538 B1 473 8.130000e-03 4.110000e-33 [1500, 3000, 4500] [30, 1500, 3000, 4500] (300, 300) 109.1 104.9 partially reproduced BH538 B2 473 8.130000e-03 4.110000e-33 [1500, 3000] [30, 1500, 3000, 4500] (300, 300) 109.1 104.9 reproduced BH539 B1 503 6.570000e-14 2.620000e-43 [30, 1500, 3000, 4500] [30, 1500, 3000, 4500] (300, 300) 108.0 99.7 not reproduced BH539 B2 503 6.570000e-14 2.620000e-43 [300, 1500, 3000, 4500] [30, 1500, 3000, 4500] (300, 300) 108.0 99.7 not reproduced BH547 B1 525 2.930000e-06 2.670000e-44 [30, 1500, 3000, 4500] [30, 1500, 3000, 4500] (300, 300) 114.0 99.7 not reproduced BH547 B2 525 2.930000e-06 2.670000e-44 [300, 1500, 3000, 4500] [1500, 3000, 4500] (300, 300) 114.0 99.7 not reproduced
t = tables[("N1", "B1")]
fig, ax = plt.subplots(figsize=(7, 3.5))
for p in PRFS:
q = t[t.prf == int(p)]
ax.errorbar(range(4), q["mean"], yerr=q["sem"], marker="o", label=f"{int(p)} Hz", capsize=3)
ax.axhline(100, color="0.5", lw=0.8)
ax.set_xticks(range(4))
ax.set_xticklabels([f"{w} ms" for w in tfus.WINDOWS])
ax.set_ylabel("normalized rate (%), mean +/- SEM")
ax.set_title("AC units, PD 200 us, reading N1")
ax.legend(fontsize=7)
plt.tight_layout()
plt.show()
print("Plotted means range from", t["mean"].min(), "to", t["mean"].max(), "%; none outside the axes.")
Plotted means range from 94.3 to 310.5 %; none outside the axes.
Result: blocked (methods_underspecified). The four readings the claims file names give different
outcomes. N1 with B1 is partially reproduced: both Kruskal-Wallis tests have p < 0.01 (7.83e-06 during
sonication, 5.16e-49 at 201 to 268 ms), the postsonication set is exactly 30, 1500, 3000 and 4500 Hz, 300 Hz
has the lowest mean in both windows (99.7% and 94.3%), but during sonication 4500 Hz is also above baseline
(108.0 +/- 1.5%, corrected p = 7.05e-05) besides 1500 Hz (140.3 +/- 8.8%) and 3000 Hz (107.1 +/- 1.2%). N1
with B2 is not reproduced: the during set also gains 300 Hz (99.7%, corrected p = 6.29e-03). N2 with B1 and
with B2 are not reproduced (Kruskal-Wallis p = 0.0213 during sonication). The claims file header blocks a
claim whose readings give different outcomes. No reading is better supported than the others: N2 sits at
49.8% to 61.7% with no stimulus (trial starts shifted by 1 s), where N1 gives 94.3% to 99.0%; under N1 the
B1 sample is 100% for every unit, so the two-sample test the paper names becomes a test against a constant;
and B2 uses a window the paper never names. No reading reproduces the claim. The result also depends on
which rats are pooled: with one rat left out, N1 with B1 gives reproduced (without BH537, the rat whose 4500
Hz rise is largest at 135.3%), partially reproduced (without BH530, BH535 or BH538) or not reproduced
(without BH539 or BH547, where 30 Hz joins the during set), and N1 with B2 gives reproduced without BH530,
BH537 or BH538. In all twelve leave-one-out runs 300 Hz has the lowest mean in both windows. The dataset
gives 633 AC units with more than 2000 spikes (expected 453 +/- 45.3) from the 6 rats with AC recordings at
PD 200 us; the seventh AC rat, BH532, has no PD 200 us recording, and the AC_1 and AC_2 files of BH538
are at 1500 Hz and PD 40 us.
Summary¶
| Claim | Outcome |
|---|---|
| fig4-unit-counts, fig6-unit-counts | Not reproduced (POm totals 36 to 44% above the paper's RSU + FSU under every selection of recordings by identifier) |
| fig11b-ac-responses | Blocked: methods underspecified (the readings give different outcomes; none reproduces the claim) |
| fig4a, fig4b, fig4-postsonication-3000, fig4-30hz-responses, fig4-fsu-no-change | Blocked: data missing (cell types) |
| fig5a-sham-no-change, fig5b-cg1-control-no-change | Blocked: data missing (cell types; no sham or CG1 sessions marked) |
| fig6a (2 claims), fig6b, fig6c, fig6d | Blocked: data missing (cell types) |
| fig7a to fig7d, fig8a, fig8b, fig9a, fig9b | Blocked: data missing (cell types) |
| fig10a to fig10f | Blocked: data missing (cell types) |
| pathway-follows-local-prf | Blocked: depends on blocked claims |
| fig11c-ac-to-s1-30hz, fig11d-ac-to-s1-3000 | Blocked: data missing (S1 cell types) |
| auditory-confound-low-prf | Blocked: depends on blocked claims |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's: the normalization and the baseline sample are not fixed by the paper, and every reading is shown; Figure 11B is blocked because they disagree. The count claims are judged on the area totals, which do not depend on the cell types. Duplicate files are counted once.
What was found about the dataset: no spike waveforms and no cell types (celltype_label is 0 for every
unit, though described as 1 = RSU and 2 = FSU), so 30 of the 33 claims cannot be tested, nor the RSU/FSU split
of the two count claims; most recordings
are deposited twice; the sham and CG1 sessions cannot be identified; the dataset has 10 rats with POm
recordings and 7 with AC recordings, where the paper reports 16 and 12, together with M1 recordings and
tFUS aimed at M1 or POm that the paper does not describe; session_start_time differs between the two area
files of one session; the STTC files lack the dt dimension their descriptions give and identify no units.
print(f"notebook ran in {(time.time() - started) / 60:.1f} min")
notebook ran in 2.9 min