NeuroPAL whole-brain activity (Yemini et al. 2021) from the paper's Zenodo record¶
Paper. Yemini et al. (2021), Cell 184:272-288, doi:10.1016/j.cell.2020.12.012; claims written from the accepted author manuscript (PMC10494711) with Methods S1 and S2.
Data. These claims use data the dandiset (DANDI:000715, version 0.241009.1514) does not hold:
the whole-brain calcium activity of 21 heads and 21 tails of strain OH16230, which the paper
deposited on Zenodo, doi:10.5281/zenodo.3906530 (one
published version; files Head_Activity_OH16230.mat and Tail_Activity_OH16230.mat in
NeuroPAL data.zip). They are data: elsewhere claims and are counted apart from the dandiset.
The dandiset's own claims are in reproduce_yemini_2021.ipynb.
| Claim | Status here |
|---|---|
| fig6d-ase-nacl-asymmetry | tested |
| fig6ef-nacl-novel-interneurons | tested |
| fig6de-asj-ric-stimulus-specific | tested |
| tables6-tail-pvq-and-post-stimulus | tested |
| figs1h-avf-cyclic-activity | tested |
| other whole-brain claims (6C, 6D AWC, 6F, 6G-I, 6J-K) | not tested; see the record |
import time
import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
import wholebrain as wb
import neuropal as npl # names the dandiset of this entry; its data are not used here
pd.set_option("display.width", 200)
pd.set_option("display.max_rows", 200)
t0 = time.time()
raw = {part: wb.load(part) for part in ("head", "tail")}
print(f"read the head and tail files from Zenodo {wb.ZENODO_DOI} in {time.time() - t0:.0f} s "
f"(dandiset of this entry: {npl.DANDISET_ID} {npl.DANDISET_VERSION}, not used here)")
read the head and tail files from Zenodo 10.5281/zenodo.3906530 in 143 s (dandiset of this entry: 000715 0.241009.1514, not used here)
Data check¶
D = {part: wb.tidy(m) for part, m in raw.items()}
for part, d in D.items():
lens = [len(x) for x in d["traces"].values()]
vals = np.concatenate([x[~np.isnan(x)] for x in d["traces"].values()])
print(f"{part}: {len(d['files'])} files, {len({k[0] for k in d['traces']})} neuron classes, {len(d['traces'])} traces, "
f"fps {d['fps']}, trace length {min(lens)} to {max(lens)} samples ({min(lens) / d['fps']:.0f} to {max(lens) / d['fps']:.0f} s), "
f"normalised values {vals.min():.2f} to {vals.max():.2f}, NaN samples per trace (median) "
f"{int(np.median([np.isnan(x).sum() for x in d['traces'].values()]))}")
print(f" stimuli {d['stim_names']}, windows {d['stim_times'].tolist()} s, orders per file:",
pd.Series([tuple(r) for r in d["stims"]]).value_counts().to_dict())
print("\nPaper: 21 heads (189 neurons) and 21 tails (42 neurons), 4 min at about 4 Hz, 10 s pulses one minute apart, in rotated order.")
head: 21 files, 109 neuron classes, 3550 traces, fps 4.0, trace length 936 to 960 samples (234 to 240 s), normalised values -12.63 to 34.24, NaN samples per trace (median) 5
stimuli ['butanone', 'pentanedione', 'nacl'], windows [[60.5, 70.5], [120.5, 130.5], [180.5, 190.5]] s, orders per file: {(np.int64(2), np.int64(1), np.int64(3)): 7, (np.int64(3), np.int64(2), np.int64(1)): 7, (np.int64(1), np.int64(3), np.int64(2)): 7}
tail: 21 files, 30 neuron classes, 797 traces, fps 3.87, trace length 899 to 929 samples (232 to 240 s), normalised values -6.40 to 11.78, NaN samples per trace (median) 5
stimuli ['butanone', 'pentanedione', 'nacl'], windows [[60.5, 70.5], [120.5, 130.5], [180.5, 190.5]] s, orders per file: {(np.int64(1), np.int64(3), np.int64(2)): 7, (np.int64(2), np.int64(1), np.int64(3)): 7, (np.int64(3), np.int64(2), np.int64(1)): 7}
Paper: 21 heads (189 neurons) and 21 tails (42 neurons), 4 min at about 4 Hz, 10 s pulses one minute apart, in rotated order.
# The order of stimuli: each row of `stims` is a permutation that is its own inverse, so the two
# readings of the README ("use these indices with both stim_names and stim_times") give the same
# windows. The labels are checked against the published responses the paper used to select the
# worms: ASE excitation to NaCl and AWC inhibition to odors, and against windows shifted in time.
H = D["head"]
print("Every order is its own inverse:", all((np.argsort(r - 1) == r - 1).all() for r in H["stims"]))
rows = []
for stim, cls, sign in (("nacl", "ASE", 1), ("butanone", "AWC", -1), ("pentanedione", "AWC", -1)):
for shift in (0, -30, 30):
Hs = dict(H, stim_times=H["stim_times"] + shift)
ch = {}
for k in H["traces"]:
if k[0] == cls:
pre, dur, _ = wb.window_means(Hs, k, stim)
ch.setdefault(k[2], []).append(dur - pre)
worms = sum(any(sign * c > 0 for c in v) for v in ch.values())
allc = np.concatenate(list(ch.values()))
rows.append(dict(stimulus=stim, neuron=cls, shift_s=shift, worms_with_published_sign=worms, of=len(ch),
mean_change=allc.mean(), paired_t_p=stats.ttest_1samp(allc, 0).pvalue))
print(pd.DataFrame(rows).round(4).to_string(index=False))
Every order is its own inverse: True
stimulus neuron shift_s worms_with_published_sign of mean_change paired_t_p
nacl ASE 0 21 21 0.7145 0.0000
nacl ASE -30 10 21 -0.0597 0.0364
nacl ASE 30 9 21 -0.0472 0.0073
butanone AWC 0 20 21 -0.9978 0.0000
butanone AWC -30 12 21 0.0880 0.3704
butanone AWC 30 17 21 -0.0026 0.9801
pentanedione AWC 0 19 21 -0.8076 0.0000
pentanedione AWC -30 20 21 -0.3587 0.0000
pentanedione AWC 30 11 21 0.2942 0.0041
With the stated windows ASE rises to NaCl and AWC falls to both odors in nearly every worm, and the changes shrink or vanish when the windows are moved 30 s, so the stimulus labels and times are consistent with the responses.
The STAR Methods remove neurons within 5 um whose traces correlate at 0.95 or more; the README says the files keep them. Three readings are run: no removal, removal of both members of each pair (the literal reading of "Neurons within 5 um distance of each other, whose traces showed at least 95% correlation, were deemed too similar and removed"), and removal of one member. The paper supports none over the others. Removing both drops about a fifth of the traces (next cell), so the readings can differ; a claim whose outcome differs between them is blocked.
import copy
readings = {}
for label, kw in (("none", None), ("both removed", dict(both=True)), ("one removed", dict(both=False))):
readings[label] = {}
for part in ("head", "tail"):
d = copy.deepcopy(D[part])
n_removed = len(wb.remove_similar(d, **kw)) if kw else 0
readings[label][part] = d
print(f"{label:13s} {part}: {n_removed} traces removed, {len(d['traces'])} kept")
res = {lab: {part: {s: pd.DataFrame(wb.responses(d, s)).set_index("neuron") for s in d["stim_names"]}
for part, d in r.items()} for lab, r in readings.items()}
for lab in res:
print(lab, "significant responses (q <= 0.05), head:", {s: int((t.q <= 0.05).sum()) for s, t in res[lab]["head"].items()},
"tail:", {s: int((t.q <= 0.05).sum()) for s, t in res[lab]["tail"].items()})
none head: 0 traces removed, 3550 kept none tail: 0 traces removed, 797 kept
both removed head: 647 traces removed, 2903 kept both removed tail: 161 traces removed, 636 kept
one removed head: 353 traces removed, 3197 kept one removed tail: 88 traces removed, 709 kept
none significant responses (q <= 0.05), head: {'butanone': 22, 'pentanedione': 45, 'nacl': 40} tail: {'butanone': 0, 'pentanedione': 1, 'nacl': 0}
both removed significant responses (q <= 0.05), head: {'butanone': 22, 'pentanedione': 28, 'nacl': 43} tail: {'butanone': 0, 'pentanedione': 1, 'nacl': 0}
one removed significant responses (q <= 0.05), head: {'butanone': 24, 'pentanedione': 33, 'nacl': 41} tail: {'butanone': 0, 'pentanedione': 1, 'nacl': 0}
Storey's pi0, the estimated share of neurons that do not respond, scales every q-value. The paper cites Storey (2002) without giving how pi0 was estimated. The main analysis estimates it at lambda = 0.5. Two other estimates are shown: the bootstrap choice of lambda over 0.01 to 0.95 (Storey 2002 section 9), the default of MATLAB's mafdr, and pi0 = 1, the Benjamini-Hochberg bound, which is not Storey's method and gives the largest q-values. The next cell prints pi0 per part and stimulus and the q-values the claims below turn on.
six = ["AIN", "AVF", "AVH", "AVJ", "I2", "MI"]
PI0 = {"lambda 0.5 (main)": 0.5, "bootstrap (mafdr default)": "bootstrap", "pi0 = 1 (BH bound)": 0}
res_pi0 = {e: {lab: {part: {s: (res[lab][part][s] if lam == 0.5 else pd.DataFrame(wb.responses(d, s, lam=lam)).set_index("neuron"))
for s in d["stim_names"]} for part, d in r.items()} for lab, r in readings.items()} for e, lam in PI0.items()}
rows = []
for e, lam in PI0.items():
for lab in readings:
h, t = res_pi0[e][lab]["head"], res_pi0[e][lab]["tail"]
rows.append(dict(pi0_estimate=e, reading=lab,
**{f"pi0_{part}_{s[:4]}": wb.storey_pi0(res_pi0[e][lab][part][s].p.to_numpy(), lam) for part in ("head", "tail") for s in ("nacl", "butanone", "pentanedione")},
six_nacl_max_q=h["nacl"].loc[six].q.max(), six_nacl_max_q_neuron=h["nacl"].loc[six].q.idxmax(),
ASJ_pent_q=h["pentanedione"].loc["ASJ"].q, RIC_nacl_q=h["nacl"].loc["RIC"].q, RIC_pent_q=h["pentanedione"].loc["RIC"].q,
PVQ_pent_q=t["pentanedione"].loc["PVQ"].q,
tail_post_responders=" / ".join(str(int((t[s].q_post <= 0.05).sum())) for s in ("nacl", "butanone", "pentanedione"))))
pi0_tab = pd.DataFrame(rows)
print("pi0 per part and stimulus, and the q-values the claims turn on (tail post-stimulus responders: NaCl / 2-butanone / 2,3-pentanedione):")
print(pi0_tab.iloc[:, :8].round(3).to_string(index=False))
print(pi0_tab.drop(columns=pi0_tab.columns[2:8]).to_string(index=False, float_format=lambda v: f"{v:.4g}"))
pi0 per part and stimulus, and the q-values the claims turn on (tail post-stimulus responders: NaCl / 2-butanone / 2,3-pentanedione):
pi0_estimate reading pi0_head_nacl pi0_head_buta pi0_head_pent pi0_tail_nacl pi0_tail_buta pi0_tail_pent
lambda 0.5 (main) none 0.468 0.574 0.383 0.526 0.632 0.947
lambda 0.5 (main) both removed 0.473 0.452 0.495 0.737 0.947 0.947
lambda 0.5 (main) one removed 0.447 0.532 0.426 0.632 0.842 0.842
bootstrap (mafdr default) none 0.433 0.515 0.365 0.000 0.263 0.831
bootstrap (mafdr default) both removed 0.400 0.452 0.430 0.351 0.000 0.770
bootstrap (mafdr default) one removed 0.411 0.521 0.397 0.658 0.686 0.000
pi0 = 1 (BH bound) none 1.000 1.000 1.000 1.000 1.000 1.000
pi0 = 1 (BH bound) both removed 1.000 1.000 1.000 1.000 1.000 1.000
pi0 = 1 (BH bound) one removed 1.000 1.000 1.000 1.000 1.000 1.000
pi0_estimate reading six_nacl_max_q six_nacl_max_q_neuron ASJ_pent_q RIC_nacl_q RIC_pent_q PVQ_pent_q tail_post_responders
lambda 0.5 (main) none 0.04096 AVF 0.01878 0.06535 0.00443 0.000587 3 / 4 / 11
lambda 0.5 (main) both removed 0.03927 AVJ 0.06288 0.04402 0.01117 0.0007736 0 / 3 / 11
lambda 0.5 (main) one removed 0.03775 AVF 0.03301 0.05376 0.0189 0.0006876 1 / 4 / 11
bootstrap (mafdr default) none 0.03787 AVF 0.01788 0.06042 0.004219 0.000515 14 / 4 / 17
bootstrap (mafdr default) both removed 0.03321 AVJ 0.05468 0.03723 0.009713 0.0006289 0 / 11 / 12
bootstrap (mafdr default) one removed 0.03469 AVF 0.03077 0.0494 0.01762 0 12 / 11 / 11
pi0 = 1 (BH bound) none 0.0875 AVF 0.04903 0.1396 0.01157 0.0006197 3 / 4 / 9
pi0 = 1 (BH bound) both removed 0.08301 AVJ 0.1271 0.09305 0.02258 0.0008165 0 / 1 / 5
pi0 = 1 (BH bound) one removed 0.08449 AVF 0.07757 0.1203 0.04441 0.0008165 1 / 1 / 6
fig6d-ase-nacl-asymmetry¶
NaCl (160 mM) evokes left/right asymmetric responses in the ASE neurons, with ASER increasing its calcium.
Criterion: two-tailed unpaired t test between left and right responses, p of at most 0.05.
rows = []
for lab, r in readings.items():
d = r["head"]
ch = {side: np.array([np.subtract(*wb.window_means(d, k, "nacl")[1::-1]) for k in d["traces"] if k[0] == "ASE" and k[1] == side])
for side in ("L", "R")}
t = stats.ttest_ind(ch["L"], ch["R"])
rows.append(dict(reading=lab, n_L=len(ch["L"]), n_R=len(ch["R"]), mean_change_ASEL=ch["L"].mean(), mean_change_ASER=ch["R"].mean(),
ASER_up_worms=int((ch["R"] > 0).sum()), t=t.statistic, p=t.pvalue,
ASER_vs_0_p=stats.ttest_1samp(ch["R"], 0).pvalue))
print(pd.DataFrame(rows).to_string(index=False))
reading n_L n_R mean_change_ASEL mean_change_ASER ASER_up_worms t p ASER_vs_0_p
none 21 21 1.013928 0.415088 14 3.816853 0.000460 0.003153
both removed 17 17 1.073725 0.461996 11 3.380186 0.001922 0.007121
one removed 21 21 1.013928 0.415088 14 3.816853 0.000460 0.003153
Result: reproduced. ASEL rises by 1.01 and ASER by 0.42 (normalised units) on average over 21 worms each, a left/right difference with t = 3.82, p = 0.00046; ASER rises in 14 of 21 worms (mean against 0, p = 0.0032). With the similar neurons removed (both members) p = 0.0019 and ASER still rises (0.46).
fig6ef-nacl-novel-interneurons¶
NaCl evokes significant responses in six interneurons with no previously known function, AIN, AVF, AVH, AVJ, I2 and MI.
six = ["AIN", "AVF", "AVH", "AVJ", "I2", "MI"]
for lab in res:
t = res[lab]["head"]["nacl"].loc[six, ["n", "mean_change", "p", "q"]]
print(f"{lab}: all six significant (q <= 0.05): {bool((t.q <= 0.05).all())}")
print(t.to_string(float_format=lambda v: f"{v:.4g}"))
none: all six significant (q <= 0.05): True
n mean_change p q
neuron
AIN 40 0.03169 0.004011 0.009288
AVF 31 0.07415 0.03351 0.04096
AVH 42 0.2308 1.72e-06 1.892e-05
AVJ 41 0.1222 0.01705 0.02886
I2 40 0.04212 0.0004223 0.001327
MI 21 0.3128 0.01308 0.02303
both removed: all six significant (q <= 0.05): True
n mean_change p q
neuron
AIN 37 0.0339 0.001304 0.004414
AVF 28 0.08357 0.02905 0.0376
AVH 29 0.2505 3.328e-05 0.0002929
AVJ 32 0.1302 0.03298 0.03927
I2 34 0.04306 0.001625 0.004766
MI 11 0.1596 0.003421 0.008818
one removed: all six significant (q <= 0.05): True
n mean_change p q
neuron
AIN 40 0.03169 0.004011 0.00819
AVF 31 0.07415 0.03351 0.03775
AVH 34 0.2474 9.228e-06 9.689e-05
AVJ 36 0.1227 0.02516 0.03262
I2 40 0.04212 0.0004223 0.00168
MI 11 0.1596 0.003421 0.007975
Result: reproduced. All six neurons respond to NaCl with q at most 0.05 under the three readings (largest q 0.041 for AVF without removal, 0.039 for AVJ with both members removed), and all six rise. Storey's pi0 for the head's NaCl tests is 0.45 to 0.47 at lambda = 0.5 and 0.40 to 0.43 with the bootstrap choice of lambda, which gives a largest q of 0.033 to 0.038; with pi0 = 1, the Benjamini-Hochberg bound and not Storey's method, AVF and AVJ have q of 0.062 to 0.088 under every reading.
fig6de-asj-ric-stimulus-specific¶
The sets of responding neurons overlap but differ by stimulus; ASJ is excited by NaCl and inhibited by both odors, and RIC is excited only by 2,3-pentanedione.
for lab in res:
rows = []
for n in ("ASJ", "RIC"):
for s in ("nacl", "butanone", "pentanedione"):
t = res[lab]["head"][s].loc[n]
rows.append(dict(neuron=n, stimulus=s, n=int(t.n), mean_change=t.mean_change, q=t.q, significant=bool(t.q <= 0.05)))
tab = pd.DataFrame(rows)
asj = tab[(tab.neuron == "ASJ")].set_index("stimulus")
ric = tab[(tab.neuron == "RIC")].set_index("stimulus")
asj_ok = asj.loc["nacl"].significant and asj.loc["nacl"].mean_change > 0 and all(asj.loc[o].significant and asj.loc[o].mean_change < 0 for o in ("butanone", "pentanedione"))
ric_ok = ric.loc["pentanedione"].significant and ric.loc["pentanedione"].mean_change > 0 and not ric.loc["nacl"].significant and not ric.loc["butanone"].significant
sets = {s: set(res[lab]["head"][s].index[res[lab]["head"][s].q <= 0.05]) for s in ("nacl", "butanone", "pentanedione")}
print(f"{lab}: ASJ as stated {asj_ok}; RIC as stated {ric_ok}; responding sets {[len(v) for v in sets.values()]}, "
f"in all three {len(set.intersection(*sets.values()))}, union {len(set.union(*sets.values()))}")
print(tab.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
none: ASJ as stated True; RIC as stated True; responding sets [40, 22, 45], in all three 10, union 65 neuron stimulus n mean_change q significant ASJ nacl 37 0.2248 0.001305 True ASJ butanone 37 -0.1566 0.0009404 True ASJ pentanedione 37 -0.1028 0.01878 True RIC nacl 42 0.01418 0.06535 False RIC butanone 42 -0.0008588 0.5528 False RIC pentanedione 42 0.03472 0.00443 True both removed: ASJ as stated False; RIC as stated False; responding sets [43, 22, 28], in all three 9, union 60 neuron stimulus n mean_change q significant ASJ nacl 32 0.2456 0.002234 True ASJ butanone 32 -0.1628 0.001725 True ASJ pentanedione 32 -0.09073 0.06288 False RIC nacl 37 0.01635 0.04402 True RIC butanone 37 -0.004085 0.366 False RIC pentanedione 37 0.03436 0.01117 True one removed: ASJ as stated True; RIC as stated True; responding sets [41, 24, 33], in all three 10, union 63 neuron stimulus n mean_change q significant ASJ nacl 36 0.2241 0.00168 True ASJ butanone 36 -0.153 0.001427 True ASJ pentanedione 36 -0.0974 0.03301 True RIC nacl 39 0.01545 0.05376 False RIC butanone 39 -0.002711 0.4694 False RIC pentanedione 39 0.02977 0.0189 True
Result: blocked: the outcome depends on a detail the paper leaves open. Without removal of similar neurons, and with one member of each pair removed, both statements hold: ASJ is excited by NaCl and inhibited by both odors (q 0.019 and 0.033 for 2,3-pentanedione), and RIC responds only to 2,3-pentanedione (q 0.0044 and 0.019; NaCl q 0.065 and 0.054). With both members removed, ASJ's fall with 2,3-pentanedione is not significant (q 0.063) and RIC also responds to NaCl (q 0.044). The responding sets overlap and differ under every reading (9 or 10 neurons respond to all three stimuli out of 60 to 65). With the bootstrap choice of Storey's lambda the same pattern holds: both statements hold without removal and with one member removed (RIC to NaCl q 0.060 and 0.049), and with both removed ASJ to 2,3-pentanedione has q 0.055 and RIC to NaCl 0.037.
tables6-tail-pvq-and-post-stimulus¶
In the tail, PVQ is significantly excited by 2,3-pentanedione, and sensory and interneurons show significant post-stimulus responses to all three stimuli.
for lab in res:
t = res[lab]["tail"]
pvq = t["pentanedione"].loc["PVQ"]
print(f"{lab}: PVQ to pentanedione n {int(pvq.n)}, mean change {pvq.mean_change:.4f}, q {pvq.q:.4g}; to butanone q "
f"{t['butanone'].loc['PVQ'].q:.3g}, to NaCl q {t['nacl'].loc['PVQ'].q:.3g}")
for s in ("nacl", "butanone", "pentanedione"):
post = t[s][(t[s].q_post <= 0.05)]
print(f" post-stimulus responders to {s}: " + ", ".join(f"{n} (q {v:.3f})" for n, v in post.q_post.items()))
none: PVQ to pentanedione n 34, mean change 0.0484, q 0.000587; to butanone q 0.369, to NaCl q 0.181 post-stimulus responders to nacl: LUA (q 0.003), PHB (q 0.003), PVW (q 0.003) post-stimulus responders to butanone: LUA (q 0.001), PHA (q 0.024), PHB (q 0.001), PVR (q 0.024) post-stimulus responders to pentanedione: ALN (q 0.015), DVA (q 0.013), DVB (q 0.006), PHA (q 0.006), PHB (q 0.004), PHC (q 0.023), PLM (q 0.049), PLN (q 0.006), PVN (q 0.008), PVQ (q 0.023), PVW (q 0.024) both removed: PVQ to pentanedione n 33, mean change 0.0491, q 0.0007736; to butanone q 0.806, to NaCl q 0.317 post-stimulus responders to nacl: post-stimulus responders to butanone: LUA (q 0.012), PHB (q 0.027), PVR (q 0.027) post-stimulus responders to pentanedione: ALN (q 0.022), DVA (q 0.036), DVB (q 0.022), LUA (q 0.037), PHA (q 0.022), PHB (q 0.036), PHC (q 0.027), PLN (q 0.022), PVN (q 0.022), PVQ (q 0.035), PVW (q 0.035) one removed: PVQ to pentanedione n 33, mean change 0.0491, q 0.0006876; to butanone q 0.586, to NaCl q 0.244 post-stimulus responders to nacl: LUA (q 0.011) post-stimulus responders to butanone: LUA (q 0.001), PHA (q 0.016), PHB (q 0.010), PVR (q 0.010) post-stimulus responders to pentanedione: ALN (q 0.018), DVA (q 0.018), DVB (q 0.018), PHA (q 0.015), PHB (q 0.034), PHC (q 0.023), PLM (q 0.049), PLN (q 0.015), PVN (q 0.017), PVQ (q 0.031), PVW (q 0.031)
Result: blocked: the outcome depends on a detail the paper leaves open. PVQ is excited by 2,3-pentanedione under every reading (mean change 0.048 to 0.049, q 0.0006 to 0.0008) and not by the other two stimuli. Post-stimulus responders exist for all three stimuli without removal (NaCl LUA, PHB and PVW) and with one member removed (NaCl LUA only), but with both members removed none remains for NaCl. The claim holds in full on the criterion under two readings and in part under the third, so the rule stated above, that a claim whose outcome differs between the readings is blocked, gives blocked (methods_underspecified), as for fig6de-asj-ric-stimulus-specific. The first record gave partially reproduced, the outcome every reading supports; the change was proposed after verification round 2 and agreed in round 3. The same holds with the bootstrap choice of Storey's lambda (none for NaCl with both members removed); there pi0 is 0 for some of the tail's small sets of tests, which sets all their q-values to 0, so the bootstrap is not informative for the tail. Which responders are sensory neurons and which interneurons is not tested: the record holds no neuron types and the claims file gives none.
figs1h-avf-cyclic-activity¶
In many animals AVF shows robust cyclical activity of about 0.3 Hz throughout the recording, regardless of stimuli.
Figure S1H counts peaks per 4-minute interval in traces median filtered at 0.5 s (3 samples, 0.75 s, the shortest odd window at 4 Hz), from five worms it does not name. The peak rule is not given; it is run at three prominences, and the frequency of the largest spectral power is shown as a check. All AVF traces are used.
d = D["head"]
keys = [k for k in d["traces"] if k[0] == "AVF"]
rows = []
for k in keys:
x = d["traces"][k]
rows.append(dict(trace=f"{k[1]}{d['files'][k[2]]}", **{f"peaks_per_s_prom{p}": wb.peak_rate(x, d["fps"], p) for p in (0.5, 1.0, 1.5)},
spectral_peak_hz=wb.spectral_peak(x, d["fps"])))
avf = pd.DataFrame(rows)
print(f"{len(avf)} AVF traces from {len({k[2] for k in keys})} worms")
print(avf.describe().loc[["mean", "25%", "50%", "75%", "min", "max"]].round(3).to_string())
for c in avf.columns[1:]:
print(f"{c}: traces within 0.27 to 0.33 Hz {int(avf[c].between(0.27, 0.33).sum())} of {len(avf)}")
31 AVF traces from 18 worms
peaks_per_s_prom0.5 peaks_per_s_prom1.0 peaks_per_s_prom1.5 spectral_peak_hz
mean 0.188 0.117 0.071 0.152
25% 0.155 0.090 0.040 0.062
50% 0.197 0.109 0.071 0.188
75% 0.218 0.151 0.101 0.234
min 0.088 0.038 0.013 0.062
max 0.277 0.235 0.164 0.297
peaks_per_s_prom0.5: traces within 0.27 to 0.33 Hz 1 of 31
peaks_per_s_prom1.0: traces within 0.27 to 0.33 Hz 0 of 31
peaks_per_s_prom1.5: traces within 0.27 to 0.33 Hz 0 of 31
spectral_peak_hz: traces within 0.27 to 0.33 Hz 3 of 31
fig, ax = plt.subplots(figsize=(9, 2.5))
k = keys[0]
x = d["traces"][k]
t = np.arange(len(x)) / d["fps"]
ax.plot(t, x, lw=0.7, color="0.2")
for s, e in d["stim_times"]:
ax.axvspan(s, e, color="C1", alpha=0.3)
ax.set_xlabel("time (s)"); ax.set_ylabel("AVF, normalised")
plt.tight_layout(); plt.show()
print(f"trace {k}: range {np.nanmin(x):.2f} to {np.nanmax(x):.2f}; no point outside the axes")
trace ('AVF', 'L', 0): range -0.59 to 1.39; no point outside the axes
Result: not reproduced, under an assumed peak rule. Over 31 AVF traces from 18 worms the median peak rate is 0.20, 0.11 and 0.07 per second at prominences of 0.5, 1 and 1.5 standard deviations, and the median frequency of the largest spectral power is 0.19 Hz; one trace (peak count) and three traces (spectrum) of 31 fall in 0.27 to 0.33 Hz. The paper's five worms are not named, so a subset at 0.3 Hz cannot be excluded; this analysis pools all 31 traces from 18 worms.
Summary¶
| Claim | Outcome |
|---|---|
| fig6d-ase-nacl-asymmetry | reproduced |
| fig6ef-nacl-novel-interneurons | reproduced |
| fig6de-asj-ric-stimulus-specific | blocked: holds without removal or with one member removed, fails with both removed |
| tables6-tail-pvq-and-post-stimulus | blocked: holds in full without removal or with one member removed, in part with both removed |
| figs1h-avf-cyclic-activity | not reproduced (rests on an assumed peak rule) |
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.
Differences from the paper's method: the removal of similar neurons (5 um, r of 0.95 or more) is run three ways because the text does not say whether both members of a pair were dropped (removing both drops 647 of 3,550 head and 161 of 797 tail traces); the post-stimulus window is taken as 10 s; Storey's pi0 is estimated at lambda = 0.5, with the bootstrap choice of lambda and pi0 = 1 shown; the peak rule of Figure S1H is not given. About the data: the record gives no neuron types and no list of the neurons the paper removed, and the files keep the similar neurons the paper excluded.