Reproducing Hadjas et al. (2025) from dandiset 001528¶
Paper. Hadjas LC, Kollman GJ, Linderhof L, Xia M, Mansur S, Saint-Pierre M, Lim B, Lee EB, Cicchetti F, Awatramani R, Hollon NG, Hnasko TS (2025). Parkinson's Disease-vulnerable and -resilient dopamine neurons display opposite responses to excitatory input. bioRxiv, doi:10.1101/2025.06.03.657460. This attempt used the full text of the preprint (version 1, CC BY), which has no data or code availability statement.
Dandiset. DANDI:001528, version draft.
The dandiset has no published version, and the paper cites none, so the draft is the only version
there is. The draft read here was last modified on 2026-09-01T15:27:47Z (290 assets); the first
cell below checks the draft's modification time against that and prints whether it has changed.
A draft can change, so the numbers below hold for that state of it.
What the dandiset holds. Fiber photometry with optogenetic stimulation of STN or PPN inputs to the substantia nigra in 120 mice: SNr GABA neurons, SNc pan-DA neurons, GRAB-DA in DLS and TS, DA terminals in DLS and TS, Vglut2+ DA neurons and Anxa1+ DA neurons (Figures 2 to 5, S3, S4). Each file holds the lab's dF/F, the 465 and 405 nm signals and the stimulation epochs. The human tissue counts of Figure 1 and Figure S1 are not in it.
Code read. The dandiset links its conversion code,
catalystneuro/hnasko-lab-to-nwb at commit
deaefa7b1776a5ffb98fd50c5ff4d7b899bdaf59 (the lotfi_2025 conversion, its metadata, and a
Python re-implementation of the lab's MATLAB pipeline with a description of the MATLAB scripts),
and an example notebook in dandi/example-notebooks
at commit 77a123e54e0d185e37768603da1faff79b62421f (001528/HnaskoLab/Lotfi_2025), which shows
how to load the files. The lab's MATLAB scripts are not public.
Outside inputs. None. Groups, inputs and sites come from the files' own metadata.
No random numbers are drawn anywhere in this notebook.
| Claim | Data | Status |
|---|---|---|
| cohort-sizes | dandiset | tested |
| fig2-gaba-excited | dandiset | tested |
| fig2-gaba-frequency | dandiset | tested |
| fig3b-da-stn-250ms-peak | dandiset | tested |
| fig3c-da-stn-duration | dandiset | tested |
| fig3d-da-stn-frequency-ns | dandiset | tested |
| fig3g-da-ppn-duration | dandiset | tested |
| fig3h-da-ppn-frequency-ns | dandiset | tested |
| fig4cf-ts-dls-opposite | dandiset | tested |
| fig4dg-ts-frequency | dandiset | tested |
| fig4eh-dls-frequency | dandiset | tested |
| fig4eh-dls-5hz-increase | dandiset | tested |
| fig4cf-dls-rebound | dandiset | tested |
| figs3-terminals-opposite | dandiset | tested |
| fig5-vglut2-excited | dandiset | tested |
| fig5-vglut2-frequency | dandiset | tested |
| fig5-anxa1-dls-inhibited | dandiset | blocked: data missing (analysis shown) |
| fig5-anxa1-dls-frequency | dandiset | tested |
| fig5-anxa1-dls-rebound | dandiset | tested |
| figs4-anxa1-soma-inhibited | dandiset | tested |
| fig5-opposite-subtypes | dandiset | tested |
| fig1de-th-anxa1-reduced, fig1f-anxa1-de-enriched, figs1-calb1-aldh1a1-counts, figs1-if-aldh1a1 (human tissue) | dandiset | blocked: data missing |
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import photometry as ph
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()
version = ph.check_version()
print("Dandiset", ph.DANDISET_ID, "version", ph.DANDISET_VERSION, "| draft modified", version["modified"],
"| assets", version["asset_count"], "| same as the pinned draft (", ph.DRAFT_MODIFIED, "):", version["matches_pinned"])
recs = ph.load_all()
print(f"loaded {len(recs)} NWB files in {time.time() - started:.0f} s")
Dandiset 001528 version draft | draft modified 2026-09-01T15:27:47.076206+00:00 | assets 290 | same as the pinned draft ( 2026-09-01T15:27:47.076206Z ): True
loaded 258 NWB files in 56 s
Data check¶
Each check of the metadata against the values is stated under its output.
info = pd.DataFrame([{k: r[k] for k in ["path", "subject", "subject_id_raw", "group", "site", "input", "protocol",
"session_id", "sex", "genotype", "raw_rate", "session_start"]}
| {"kind": ph.recording_kind(r), "fp_table": "/".join(r["fp_table_locations"])} for r in recs])
print("NWB files:", len(info), "| mice:", info.subject.nunique(), "| groups per mouse (max):", info.groupby("subject").group.nunique().max())
print("\nMice per group, input, recorded site and protocol:")
print(info.groupby(["group", "input", "site", "protocol"]).subject.nunique().unstack().to_string())
print("\nGenotype per group:")
print(info.groupby("group").genotype.unique().to_string())
print("\nSex per group (mice):")
print(info.drop_duplicates("subject").groupby(["group", "sex"]).size().unstack(fill_value=0).to_string())
print("\nRaw sampling rate (Hz) per group:")
print(info.groupby(["group", "raw_rate"]).subject.nunique().unstack(fill_value=0).to_string())
print("\nSubject ids written in lower case in the files:", sorted(info[info.subject_id_raw != info.subject].subject_id_raw.unique()))
print("Session start time zone offsets:", info.session_start.str[-6:].value_counts().to_dict())
NWB files: 258 | mice: 120 | groups per mouse (max): 1
Mice per group, input, recorded site and protocol:
protocol durations frequencies
group input site
Anxa1 PPN DLS 4 4
SNc 10 10
STN DLS 5 5
SNc 8 8
DAterm STN DLS 8 8
TS 9 9
GABA PPN SNr 8 8
STN SNr 7 7
GRABDA PPN DLS 9 9
TS 9 9
STN DLS 7 7
TS 6 6
Vglut2 PPN SNc 10 10
STN SNc 9 9
panDA PPN SNc 9 9
STN SNc 11 11
Genotype per group:
group
Anxa1 [Vglut2-phiC X Anxa1-cre]
DAterm [VGlut2-cre x DAT-flp]
GABA [VGlut2-cre x VGAT-flp]
GRABDA [VGlut2-cre]
Vglut2 [VGlut2-cre X DAT-flp]
panDA [VGLUT2-Cre x DAT-Flp mice]
Sex per group (mice):
sex F M
group
Anxa1 7 11
DAterm 7 10
GABA 8 7
GRABDA 18 13
Vglut2 12 7
panDA 14 6
Raw sampling rate (Hz) per group:
raw_rate 6103.515625 24414.062500 24414.062500
group
Anxa1 18 0 0
DAterm 0 17 17
GABA 8 7 7
GRABDA 31 0 0
Vglut2 16 3 3
panDA 14 6 6
Subject ids written in lower case in the files: ['c3015', 'c3016', 'c4379', 'c4708', 'c4709', 'c4977', 'c4978', 'c5113']
Session start time zone offsets: {'-00:01': 258}
The files hold 120 mice, each in one group, which fits the dandiset's summary (120 subjects). The
group of each file is read from its experiment description, the site from the name of its dF/F
series, and the input from the name of its opsin (ChrimsonR in STN or PPN). For the Anxa1 mice,
nine have separate DLS and SNc files (session ids ending in _DLS and _SNc), and nine have one
pair of files without a site suffix, whose dF/F series and raw signal are named SNc; the fiber
photometry table of every Anxa1 file lists both an SNc and a DLS fiber. Every session start time
has the offset -00:01, the local mean time of London that pytz gives when a time zone is
attached with replace; the true offset is lost (a metadata error, it does not affect the
analysis). Eight subject ids are written in lower case (c3015 and others) while the others are
upper case.
rows = []
for r in recs:
dt = np.diff(r["dff_t"])
rows.append({"path": r["path"], "protocol": r["protocol"], "dff_unit": r["dff_unit"], "gc_unit": r["gc_unit"],
"n": len(r["dff"]), "t0": r["dff_t"][0], "dt_median_ms": 1000 * np.median(dt),
"regular": bool(np.all(np.abs(dt - 0.01) < 1e-6) or (np.abs(dt[dt < 1] - 0.01) < 1e-6).all()),
"longest_gap_s": dt.max(), "same_times": np.array_equal(r["dff_t"], r["gc_t"]) and np.array_equal(r["dff_t"], r["af_t"]),
"timestamps": r["dff_has_timestamps"], "nan_frac": np.isnan(r["dff"]).mean(),
"dff_p01": np.nanpercentile(r["dff"], 0.1), "dff_p999": np.nanpercentile(r["dff"], 99.9),
"gc_median": np.nanmedian(r["gc"]), "af_median": np.nanmedian(r["af"])})
series = pd.DataFrame(rows)
print("Stated units of dF/F and of the 465 nm signal:", series.dff_unit.unique().tolist(), series.gc_unit.unique().tolist())
print("All 100 Hz series start at 0 s:", bool((series.t0 == 0).all()), "| median sample interval (ms):", series.dt_median_ms.round(4).unique().tolist())
print("dF/F, 465 and 405 nm series share one time base in every file:", bool(series.same_times.all()))
print("Files whose dF/F has explicit timestamps:", int(series.timestamps.sum()), "of", len(series))
print("Files with a gap in time above 1 s:", int((series.longest_gap_s > 1).sum()), "; longest gaps (s):", series.longest_gap_s.nlargest(3).round(1).tolist())
print("Files with missing dF/F values:", int((series.nan_frac > 0).sum()), "; by protocol:", series[series.nan_frac > 0].protocol.value_counts().to_dict(),
"; largest fraction:", round(series.nan_frac.max(), 3))
print("dF/F 0.1th percentile across files: %.1f to %.1f; 99.9th percentile: %.1f to %.1f" %
(series.dff_p01.min(), series.dff_p01.max(), series.dff_p999.min(), series.dff_p999.max()))
print("465 nm median across files: %.4f to %.4f; 405 nm: %.4f to %.4f (a.u.)" %
(series.gc_median.min(), series.gc_median.max(), series.af_median.min(), series.af_median.max()))
Stated units of dF/F and of the 465 nm signal: ['a.u.'] ['a.u.']
All 100 Hz series start at 0 s: True | median sample interval (ms): [10.0]
dF/F, 465 and 405 nm series share one time base in every file: True
Files whose dF/F has explicit timestamps: 89 of 258
Files with a gap in time above 1 s: 34 ; longest gaps (s): [171981.0, 171929.9, 171904.8]
Files with missing dF/F values: 88 ; by protocol: {'frequencies': 88} ; largest fraction: 0.548
dF/F 0.1th percentile across files: -89.4 to -1.5; 99.9th percentile: 1.7 to 249.4
465 nm median across files: 0.0131 to 0.6156; 405 nm: 0.0137 to 0.3375 (a.u.)
All 100 Hz series start at 0 s and are sampled every 10 ms; the dF/F and the two signals share
one time base. Some files of the frequencies protocol join separate recordings: their timestamps
jump (up to about two days, below) and the joins are filled with missing values. The dF/F is
stated in "a.u.", and the series description gives (Gc - fitted_af) / fitted_af; its values (for
example 99.9th percentiles of 1.7 to 249) are in percent, as in the MATLAB pipeline (dF = 100 * (Gc - fit) / fit), and the recomputation below confirms this. The unit is a metadata error; the
analysis uses the values as percent.
eps = []
for r in recs:
e = ph.label_epochs(r).assign(path=r["path"], group=r["group"], protocol=r["protocol"])
eps.append(e)
eps = pd.concat(eps, ignore_index=True)
print("Stimulation power (mW) and wavelength (nm) in the epoch tables:", sorted(eps.power_mW.unique()), sorted(eps.wavelength_nm.unique()))
print("\nEpochs per protocol, train duration and frequency (files, epochs, dropped as not a train, outside the signal, not contiguous):")
tab = eps.groupby(["protocol", "duration", "frequency"]).agg(
files=("path", "nunique"), epochs=("start", "size"), not_train=("is_train", lambda x: int((~x).sum())),
outside=("in_signal", lambda x: int((~x).sum())), not_contiguous=("contiguous", lambda x: int((~x).sum())),
usable=("usable", "sum"), with_missing=("nan_in_window", "sum"))
print(tab.to_string())
per_file = eps[eps.usable].groupby(["path", "protocol", "duration", "frequency"]).size().rename("trials").reset_index()
print("\nTrials per file and condition:", per_file.trials.value_counts().sort_index().to_dict())
ten = sorted(per_file[per_file.trials >= 10].path.str.split("/").str[0].unique())
print("Files with 10 trials per condition belong to", len(ten), "mice:", ten)
sixteen = eps[eps.duration == 16.0].path.str.split("/").str[0].unique()
print("Mice with 16 s trains:", len(sixteen), "| groups:", eps[eps.duration == 16.0].group.unique().tolist())
four = eps[(eps.protocol == "frequencies") & (eps.duration == 4.0)].path.nunique()
print("Frequencies-protocol files with 4 s trains:", four, "of", eps[eps.protocol == "frequencies"].path.nunique())
print("Number of pulses against frequency x duration (median ratio):",
round(float(np.median(eps[eps.usable].pulses / (eps[eps.usable].frequency * eps[eps.usable].duration))), 3))
no4 = sorted(set(eps[eps.protocol == "frequencies"].path) - set(eps[(eps.protocol == "frequencies") & (eps.duration == 4.0)].path))
print("Frequencies-protocol files without 4 s trains:", [p.split("/")[0] for p in no4])
gaps = []
for r in recs:
e = ph.label_epochs(r)
e = e[e.usable].sort_values("start")
gaps.append(pd.DataFrame({"protocol": r["protocol"], "onset_to_onset": np.diff(e.start.values), "offset_to_onset": e.start.values[1:] - e.stop.values[:-1]}))
gaps = pd.concat(gaps)
gaps = gaps[gaps.onset_to_onset < 100]
print("Median interval between consecutive trains (s), onset to onset and offset to next onset:")
print(gaps.groupby("protocol").median().round(2).to_string())
print("Train lengths within 15 ms of the labelled duration:", round(float((np.abs(eps[eps.usable].length - eps[eps.usable].duration) < 0.015).mean()), 3))
Stimulation power (mW) and wavelength (nm) in the epoch tables: [np.float64(5.0)] [np.float64(635.0)]
Epochs per protocol, train duration and frequency (files, epochs, dropped as not a train, outside the signal, not contiguous):
files epochs not_train outside not_contiguous usable with_missing
protocol duration frequency
durations 0.25 40.0 129 710 0 1 1 709 0
1.00 40.0 129 710 0 1 1 709 0
4.00 40.0 129 710 0 4 4 706 0
16.00 40.0 30 215 0 0 0 215 0
frequencies 0.25 5.0 126 699 3 0 0 696 0
10.0 129 714 3 0 0 711 0
20.0 129 714 3 0 0 711 0
40.0 129 711 0 0 0 711 0
1.00 5.0 126 695 0 0 0 695 0
10.0 129 710 0 0 0 710 0
20.0 129 710 0 0 0 710 0
40.0 129 710 0 0 0 710 0
4.00 5.0 112 560 0 0 0 560 0
10.0 115 575 0 0 0 575 0
20.0 115 575 0 0 0 575 0
40.0 115 575 0 0 0 575 0
Trials per file and condition: {4: 6, 5: 1734, 6: 4, 10: 156}
Files with 10 trials per condition belong to 13 mice: ['sub-B7514', 'sub-B7515', 'sub-B7532', 'sub-B7616', 'sub-B8364', 'sub-B8369', 'sub-B8370', 'sub-B8427', 'sub-B8429', 'sub-B8569', 'sub-B8603', 'sub-B8612', 'sub-B8677']
Mice with 16 s trains: 30 | groups: ['GABA', 'panDA', 'DAterm']
Frequencies-protocol files with 4 s trains: 115 of 129
Number of pulses against frequency x duration (median ratio): 1.0
Frequencies-protocol files without 4 s trains: ['sub-B7514', 'sub-B7515', 'sub-B7532', 'sub-B7616', 'sub-B8364', 'sub-B8369', 'sub-B8370', 'sub-B8427', 'sub-B8429', 'sub-B8569', 'sub-B8603', 'sub-B8612', 'sub-B8677', 'sub-C4608']
Median interval between consecutive trains (s), onset to onset and offset to next onset:
onset_to_onset offset_to_onset
protocol
durations 11.05 10.05
frequencies 21.00 20.00
Train lengths within 15 ms of the labelled duration: 0.995
The epoch tables give 5 mW at 635 nm for every train, as the Methods state, and the number of
pulses in each train fits its frequency and length. The intervals between trains are 10 s in the
durations sessions and 20 s in the frequencies sessions (offset to next onset), as the Methods
state; the files' session description of the frequencies protocol says 10 s, which the data do
not support. Three things differ from the Methods. First, 13 mice (SNr GABA and SNc pan-DA mice
with ids starting with B) have 10 trials per condition, where the Methods state 5. Second, 30
mice (SNr GABA, pan-DA and DA terminal mice recorded on the rig with a 24 kHz raw rate) also have
16 s trains, which the paper does not mention; they are not used. Third, 14 frequencies sessions
have no 4 s trains (the 13 mice above, and C4608, whose 4 s block the conversion notes describe
as a separate recording that is not in the dandiset). Fifteen epochs are dropped: nine
zero-length entries and six trains whose window runs past the end of the recording or across a
gap. The conversion code (convert_session.py) names the 1 s block of the STN to SNr GABA
frequencies sessions s1s10mW, which suggests 10 mW, while the epoch tables and the Methods give
5 mW; the data cannot settle this.
# Recordings deposited twice: identical dF/F, or near-identical traces of the same length.
md5 = pd.Series({r["path"]: r["dff_md5"] for r in recs})
print("Files sharing an identical dF/F array:", int(md5.duplicated(keep=False).sum()))
by_len = {}
for r in recs:
by_len.setdefault(len(r["dff"]), []).append(r)
pairs = []
for n, rs in by_len.items():
for i in range(len(rs)):
for j in range(i + 1, len(rs)):
a, b = rs[i]["dff"], rs[j]["dff"]
ok = np.isfinite(a) & np.isfinite(b)
if ok.sum() > 100 and np.corrcoef(a[ok], b[ok])[0, 1] > 0.99:
pairs.append((rs[i]["path"], rs[j]["path"]))
print("Pairs of same-length files with dF/F correlation above 0.99:", pairs)
starts = info.groupby(["subject", "session_start"]).size()
print("Mice with two files starting at the same time:", starts[starts > 1].index.tolist())
Files sharing an identical dF/F array: 0 Pairs of same-length files with dF/F correlation above 0.99: [] Mice with two files starting at the same time: []
No recording is deposited twice: no two files share a dF/F array, no two files of the same length have nearly the same trace, and no mouse has two files with the same start time. Each recording is counted once.
t0 = time.time()
rec_rows = []
for r in recs:
x = ph.recompute_dff(r["gc"], r["af"], r["dff_t"])
ok = np.isfinite(x) & np.isfinite(r["dff"])
rec_rows.append({"path": r["path"], "protocol": r["protocol"], "r": np.corrcoef(x[ok], r["dff"][ok])[0, 1],
"median_abs_diff": np.median(np.abs(x[ok] - r["dff"][ok])), "dff_sd": np.nanstd(r["dff"])})
recomp = pd.DataFrame(rec_rows)
print(f"recomputed dF/F for {len(recomp)} files in {time.time() - t0:.0f} s")
print(recomp.groupby("protocol").agg(files=("r", "size"), r_median=("r", "median"), r_min=("r", "min"),
above_0_99=("r", lambda x: int((x > 0.99).sum())),
abs_diff_median=("median_abs_diff", "median")).round(4).to_string())
print("\nLowest correlations:")
print(recomp.nsmallest(6, "r").round(3).to_string(index=False))
recomputed dF/F for 258 files in 370 s
files r_median r_min above_0_99 abs_diff_median
protocol
durations 129 0.9999 0.3364 113 0.0711
frequencies 129 0.9988 0.3091 104 0.0511
Lowest correlations:
path protocol r median_abs_diff dff_sd
sub-B7532/sub-B7532_ses-varying-frequencies.nwb frequencies 0.309 1.754 19.032
sub-C2618/sub-C2618_ses-varying-durations.nwb durations 0.336 0.970 7.457
sub-B7514/sub-B7514_ses-varying-frequencies.nwb frequencies 0.441 0.243 3.793
sub-C5906/sub-C5906_ses-varying-frequencies.nwb frequencies 0.560 0.033 1.027
sub-C3798/sub-C3798_ses-varying-durations.nwb durations 0.679 0.455 5.906
sub-C4561/sub-C4561_ses-varying-durations_image.nwb durations 0.691 1.022 3.371
The stored dF/F agrees with a recomputation from the stored 465 and 405 nm signals following the Methods (biexponential fit of the 405 nm signal, scaled to the 465 nm signal by robust regression, in percent): median correlation 0.9999 in durations sessions and 0.9988 in frequencies sessions, above 0.99 in 217 of 258 files, with a median absolute difference below 0.1 % dF/F. This is consistent with the stored series being the lab's dF/F in percent. In a few files the two differ more (correlation down to 0.31), where the fit of a segment differs; the lab's exact fitting options are not public. The analysis uses the stored dF/F.
# Clock check: responses at the true onsets against onsets shifted by -3 s and +3 s (within the
# 10 to 20 s between trains). Per recording, the mean over conditions of |during - pre|.
def abs_response(r, shift):
ct = ph.condition_traces(r, shift=shift)
vals = [abs(ph.window_mean(c.trace, c.t, 0, c.duration) - ph.window_mean(c.trace, c.t, -c.duration, 0))
for _, c in ct.iterrows()]
return float(np.mean(vals)) if vals else np.nan
clock = pd.DataFrame([{"path": r["path"], "subject": r["subject"], "group": r["group"], "input": r["input"], "site": r["site"], "protocol": r["protocol"],
**{f"shift_{s:+.0f}": abs_response(r, s) for s in (0.0, -3.0, 3.0)}} for r in recs])
clock["ratio"] = clock["shift_+0"] / clock[["shift_-3", "shift_+3"]].max(axis=1)
print("Median |during - pre| (% dF/F) at true and shifted onsets, and the median ratio true / largest shifted:")
print(clock.groupby(["group", "protocol"])[["shift_+0", "shift_-3", "shift_+3", "ratio"]].median().round(2).to_string())
print("\nRecordings whose true-onset response exceeds both shifted ones:", int((clock.ratio > 1).sum()), "of", len(clock))
print("Recordings where it does not (the clock could not be shown from the response):")
print(clock[clock.ratio <= 1][["path", "group", "input", "site", "shift_+0", "shift_-3", "shift_+3"]].round(2).to_string(index=False))
print("\nRecordings where it does not, by group, input and protocol:")
print(clock[clock.ratio <= 1].groupby(["group", "input", "protocol"]).size().to_string())
Median |during - pre| (% dF/F) at true and shifted onsets, and the median ratio true / largest shifted:
shift_+0 shift_-3 shift_+3 ratio
group protocol
Anxa1 durations 0.54 0.17 0.29 2.11
frequencies 0.35 0.14 0.18 1.80
DAterm durations 3.30 1.19 1.50 2.75
frequencies 2.62 0.62 1.05 1.79
GABA durations 9.03 1.82 0.89 4.40
frequencies 3.87 0.43 0.71 6.27
GRABDA durations 6.38 1.73 3.27 2.70
frequencies 3.67 1.21 1.88 2.03
Vglut2 durations 3.83 0.79 1.41 3.14
frequencies 3.12 0.87 1.44 1.92
panDA durations 2.13 1.08 1.18 1.19
frequencies 1.86 0.73 0.93 1.56
Recordings whose true-onset response exceeds both shifted ones: 235 of 258
Recordings where it does not (the clock could not be shown from the response):
path group input site shift_+0 shift_-3 shift_+3
sub-C3164/sub-C3164_ses-varying-durations.nwb panDA PPN SNc 0.84 1.74 0.28
sub-C3164/sub-C3164_ses-varying-frequencies.nwb panDA PPN SNc 0.75 0.94 0.76
sub-C3163/sub-C3163_ses-varying-frequencies.nwb panDA PPN SNc 0.14 0.09 0.16
sub-C4217/sub-C4217_ses-varying-durations.nwb panDA STN SNc 1.09 1.14 1.95
sub-C4218/sub-C4218_ses-varying-durations.nwb panDA STN SNc 2.63 1.01 2.70
sub-C4217/sub-C4217_ses-varying-frequencies.nwb panDA STN SNc 0.55 0.32 0.97
sub-C4279/sub-C4279_ses-varying-durations.nwb panDA STN SNc 2.59 2.05 5.20
sub-C4407/sub-C4407_ses-varying-durations.nwb panDA STN SNc 2.41 2.87 5.47
sub-C4279/sub-C4279_ses-varying-frequencies.nwb panDA STN SNc 1.36 1.09 1.95
sub-C4407/sub-C4407_ses-varying-frequencies.nwb panDA STN SNc 2.19 1.53 2.69
sub-C4254/sub-C4254_ses-varying-frequencies.nwb GRABDA PPN DLS 1.09 1.39 1.30
sub-C5007/sub-C5007_ses-varying-durations.nwb panDA PPN SNc 0.27 0.25 0.38
sub-C5008/sub-C5008_ses-varying-durations.nwb panDA PPN SNc 0.30 0.70 2.02
sub-B8427/sub-B8427_ses-varying-frequencies.nwb panDA STN SNc 0.51 0.23 0.56
sub-C0176/sub-C0176_ses-varying-durations.nwb DAterm STN DLS 1.67 2.05 1.84
sub-C0177/sub-C0177_ses-varying-durations.nwb DAterm STN DLS 3.30 1.57 4.35
sub-C1498/sub-C1498_ses-varying-frequencies.nwb Vglut2 STN SNc 3.75 2.33 4.49
sub-C5903/sub-C5903_ses-varying-frequencies.nwb Anxa1 STN SNc 0.16 0.14 0.17
sub-C5992/sub-C5992_ses-varying-durations.nwb Anxa1 PPN SNc 0.20 0.10 0.21
sub-C6900/sub-C6900_ses-varying-frequencies.nwb Anxa1 STN SNc 0.16 0.13 0.25
sub-C7244/sub-C7244_ses-varying-durations.nwb Anxa1 PPN SNc 0.12 0.08 0.14
sub-C7248/sub-C7248_ses-varying-frequencies.nwb Anxa1 PPN SNc 0.13 0.14 0.17
sub-C7249/sub-C7249_ses-varying-frequencies.nwb Anxa1 PPN SNc 0.13 0.15 0.12
Recordings where it does not, by group, input and protocol:
group input protocol
Anxa1 PPN durations 2
frequencies 2
STN frequencies 2
DAterm STN durations 2
GRABDA PPN frequencies 1
Vglut2 STN frequencies 1
panDA PPN durations 3
frequencies 2
STN durations 4
frequencies 4
In every group the median response at the true onsets is larger than at onsets shifted by 3 s either way (median ratio 1.2 to 6.3 by group and protocol), and in 235 of 258 recordings the response at the true onsets exceeds both shifted ones. The epoch tables share the clock of the dF/F. The 23 recordings listed last, mostly SNc pan-DA recordings with small and variable responses and unsuffixed Anxa1 recordings, are those where the clock could not be shown from the response; they are used as they are. By input, they include three PPN pan-DA durations recordings (C3164, C5007, C5008), four STN pan-DA durations recordings, two DA terminal DLS recordings and one frequencies recording each of GRAB-DA DLS and Vglut2 STN.
# The Anxa1 recordings: the mice with DLS and SNc files, and those with one unsuffixed file.
anx = []
for r in recs:
if r["group"] != "Anxa1":
continue
ct = ph.condition_traces(r)
c = ct[(ct.duration == 1.0) & (ct.frequency == 40.0)].iloc[0]
m = ph.response_measures(c.trace, c.t, 1.0)
anx.append({"subject": r["subject"], "input": r["input"], "protocol": r["protocol"], "kind": ph.recording_kind(r),
"site": r["site"], "label": r["site"] if ph.recording_kind(r) == "dual" else "unsuffixed (SNc)",
"dff_sd": np.nanstd(r["dff"]), "during_1s": m["during"], "post_1s": m["post"], "start": r["session_start"][:10]})
anx = pd.DataFrame(anx)
print("Mice per input and kind of Anxa1 recording:")
print(anx[anx.protocol == "durations"].groupby(["input", "label"]).subject.nunique().unstack().to_string())
print("\nMedian over mice (durations protocol, 1 s 40 Hz trains, % dF/F):")
print(anx[anx.protocol == "durations"].groupby("label")[["dff_sd", "during_1s", "post_1s"]].median().round(2).to_string())
Mice per input and kind of Anxa1 recording:
label DLS SNc unsuffixed (SNc)
input
PPN 4 4 6
STN 5 5 3
Median over mice (durations protocol, 1 s 40 Hz trains, % dF/F):
dff_sd during_1s post_1s
label
DLS 1.34 -1.31 0.97
SNc 0.75 -0.40 0.07
unsuffixed (SNc) 0.69 -0.27 -0.06
The paper reports Anxa1+ DLS terminals in 8 (STN) and 10 (PPN) mice and SNc cell bodies in 5 and 4. The dandiset holds DLS recordings for 5 and 4 mice, the mice that also have an SNc file, and SNc-labeled recordings for 8 and 10. The paper's DLS counts equal the dandiset's mice with a DLS file plus those with an unsuffixed file (5 + 3 and 4 + 6), and its SNc counts equal the mice with both files. Two explanations fit the counts: the unsuffixed recordings are DLS recordings labelled SNc, or the dandiset lacks the DLS recordings of nine mice and holds SNc recordings the paper did not analyse ("cell body signals were only detected in a subset of mice"). The signals favour the second: the unsuffixed recordings resemble the SNc recordings of the dual mice (small inhibition, no rebound, similar variability), not their DLS recordings (larger inhibition and a rebound). The tests of figs4-anxa1-soma-inhibited and fig5-anxa1-dls-inhibited below bear on this too: the paper's cell-body statistics are matched by the dual mice's SNc recordings alone, and its terminal statistics are not matched by adding the unsuffixed recordings to DLS. This is consistent with, but does not prove, that the labels are right and that the DLS recordings of nine mice are missing. The analysis keeps the dandiset's labels and shows the other reading as a variant.
R = ph.response_table(recs)
print("Response table:", len(R), "recording x condition rows;", R.subject.nunique(), "mice")
print("Trials per mouse and condition:", R.n_trials.value_counts().sort_index().to_dict(), "| dropped trials:", int(R.n_dropped.sum()))
for col in ["during", "pre", "post"]:
print(f"{col}: {R[col].min():.2f} to {R[col].max():.2f} % dF/F")
print("Peak latency range (s):", R.peak_latency.min(), "to", R.peak_latency.max())
one = R[np.isclose(R.duration, 1.0)]
print(f"1 s trains: largest |pre| {one.pre.abs().max():.2g} % (the pre-stimulus window is the baseline window); "
f"other durations: |pre| up to {R[~np.isclose(R.duration, 1.0)].pre.abs().max():.2f} %")
Response table: 1900 recording x condition rows; 120 mice
Trials per mouse and condition: {4: 6, 5: 1734, 6: 4, 10: 156} | dropped trials: 15
during: -25.37 to 89.23 % dF/F
pre: -5.00 to 4.74 % dF/F
post: -13.64 to 59.19 % dF/F
Peak latency range (s): 0.0 to 2.24
1 s trains: largest |pre| 2.5e-15 % (the pre-stimulus window is the baseline window); other durations: |pre| up to 5.00 %
The response measures, per mouse and condition, are computed once here and used by every claim.
For 1 s trains the pre-stimulus window (1 s before onset) is the baseline window, so pre is 0
by construction (printed above). Every paired test of a 1 s train against its pre-stimulus window
is therefore a one-sample test of the response against 0, and in the two-way mixed ANOVAs on 1 s
trains the site F equals the site by time interaction F. The paper's two-way ANOVA values show
the same pattern (site F 40.9 and interaction of the same size).
def cohort(group, inp, site=None, protocol="durations", duration=1.0, frequency=40.0, kind=None):
"""One row per mouse: the response measures for one condition."""
q = (R.group == group) & (R.input == inp) & (R.protocol == protocol) & np.isclose(R.duration, duration) & np.isclose(R.frequency, frequency)
if site is not None:
q &= R.site == site
if kind is not None:
q &= R.kind == kind
out = R[q]
assert not out.subject.duplicated().any(), "a mouse appears twice"
return out.set_index("subject")
def freq_matrix(group, inp, duration, site=None, kind=None, measure="during"):
"""Mice x frequency (5, 10, 20, 40 Hz) matrix of one measure, frequencies protocol."""
q = (R.group == group) & (R.input == inp) & (R.protocol == "frequencies") & np.isclose(R.duration, duration)
if site is not None:
q &= R.site == site
if kind is not None:
q &= R.kind == kind
m = R[q].pivot_table(index="subject", columns="frequency", values=measure)
return m.reindex(columns=list(ph.FREQUENCIES))
def fmt(res):
keep = {k: (round(v, 4) if isinstance(v, float) else v) for k, v in res.items()}
return ", ".join(f"{k} {v}" for k, v in keep.items())
def plot_traces(ax, group, inp, site=None, protocol="durations", duration=1.0, frequency=40.0, kind=None, label=None):
q = (R.group == group) & (R.input == inp) & (R.protocol == protocol) & np.isclose(R.duration, duration) & np.isclose(R.frequency, frequency)
if site is not None:
q &= R.site == site
if kind is not None:
q &= R.kind == kind
paths = R[q].path.tolist()
traces = []
for r in recs:
if r["path"] in paths:
ct = ph.condition_traces(r)
c = ct[np.isclose(ct.duration, duration) & np.isclose(ct.frequency, frequency)].iloc[0]
traces.append(c.trace)
t = c.t
m = np.vstack(traces)
mean, sem = np.nanmean(m, axis=0), np.nanstd(m, axis=0) / np.sqrt(len(m))
ax.plot(t, mean, label=f"{label or site} (n={len(m)})")
ax.fill_between(t, mean - sem, mean + sem, alpha=0.3)
ax.axvspan(0, duration, color="red", alpha=0.08)
ax.set_xlim(-3, 8)
return float(np.nanmin(mean)), float(np.nanmax(mean))
cohort-sizes¶
The number of mice per experimental group is as given in the figure legends.
paper_n = {("GABA", "STN", "SNr"): 7, ("GABA", "PPN", "SNr"): 8, ("panDA", "STN", "SNc"): 11, ("panDA", "PPN", "SNc"): 9,
("GRABDA", "STN", "TS"): 6, ("GRABDA", "STN", "DLS"): 7, ("GRABDA", "PPN", "TS"): 9, ("GRABDA", "PPN", "DLS"): 9,
("DAterm", "STN", "TS"): 9, ("DAterm", "STN", "DLS"): 9, ("Vglut2", "STN", "SNc"): 10, ("Vglut2", "PPN", "SNc"): 10,
("Anxa1", "STN", "DLS"): 8, ("Anxa1", "PPN", "DLS"): 10, ("Anxa1", "STN", "SNc"): 5, ("Anxa1", "PPN", "SNc"): 4}
dur_n = R[R.protocol == "durations"].groupby(["group", "input", "site"]).subject.nunique()
freq_n = R[R.protocol == "frequencies"].groupby(["group", "input", "site"]).subject.nunique()
alt = R[(R.protocol == "durations") & (R.group == "Anxa1")].assign(
site_alt=lambda d: np.where(d.kind == "single", "DLS", d.site)).groupby(["input", "site_alt"]).subject.nunique()
ctab = pd.DataFrame([{"group": g, "input": i, "site": s, "paper": n, "dandiset_durations": int(dur_n.get((g, i, s), 0)),
"dandiset_frequencies": int(freq_n.get((g, i, s), 0)),
"anxa1_unsuffixed_as_DLS": int(alt.get((i, s), 0)) if g == "Anxa1" else np.nan} for (g, i, s), n in paper_n.items()])
ctab["match"] = ctab.paper == ctab.dandiset_durations
print(ctab.to_string(index=False))
print("\nCounts matching:", int(ctab.match.sum()), "of", len(ctab))
a = ctab[ctab.group == "Anxa1"]
print("Anxa1 counts matching if the unsuffixed recordings were DLS:", int((a.paper == a.anxa1_unsuffixed_as_DLS).sum()), "of", len(a))
print("Mice by sex across the dandiset:", info.drop_duplicates("subject").sex.value_counts().to_dict())
dur_ids = R[R.protocol == "durations"].groupby(["group", "input", "site"]).subject.unique()
for key in [("DAterm", "STN", "TS"), ("DAterm", "STN", "DLS"), ("Vglut2", "STN", "SNc"), ("Vglut2", "PPN", "SNc")]:
print(f"{key}: {len(dur_ids[key])} mice", sorted(dur_ids[key]))
both = set(dur_ids[("DAterm", "STN", "TS")]) & set(dur_ids[("DAterm", "STN", "DLS")])
print("DA terminal mice recorded in both TS and DLS:", len(both))
print("\nVglut2 STN: exact two-sided Wilcoxon p of the paper's W (paper p 0.012, 0.008, 0.82):")
for w in (41, 43, -5):
print(f" W = {w}: n 9 -> p {ph.wilcoxon_exact_p(w, 9):.4f}; n 10 -> p {ph.wilcoxon_exact_p(w, 10):.4f}")
group input site paper dandiset_durations dandiset_frequencies anxa1_unsuffixed_as_DLS match
GABA STN SNr 7 7 7 NaN True
GABA PPN SNr 8 8 8 NaN True
panDA STN SNc 11 11 11 NaN True
panDA PPN SNc 9 9 9 NaN True
GRABDA STN TS 6 6 6 NaN True
GRABDA STN DLS 7 7 7 NaN True
GRABDA PPN TS 9 9 9 NaN True
GRABDA PPN DLS 9 9 9 NaN True
DAterm STN TS 9 9 9 NaN True
DAterm STN DLS 9 8 8 NaN False
Vglut2 STN SNc 10 9 9 NaN False
Vglut2 PPN SNc 10 10 10 NaN True
Anxa1 STN DLS 8 5 5 8.0 False
Anxa1 PPN DLS 10 4 4 10.0 False
Anxa1 STN SNc 5 8 8 5.0 False
Anxa1 PPN SNc 4 10 10 4.0 False
Counts matching: 10 of 16
Anxa1 counts matching if the unsuffixed recordings were DLS: 4 of 4
Mice by sex across the dandiset: {'F': 66, 'M': 54}
('DAterm', 'STN', 'TS'): 9 mice ['B8821', 'B9071', 'B9262', 'B9266', 'B9962', 'C0679', 'C0680', 'C0683', 'C0932']
('DAterm', 'STN', 'DLS'): 8 mice ['B8627', 'B8630', 'C0176', 'C0177', 'C0178', 'C0180', 'C0925', 'C0926']
('Vglut2', 'STN', 'SNc'): 9 mice ['C1496', 'C1498', 'C1499', 'C2450', 'C2454', 'C2733', 'C4871', 'C4875', 'C4891']
('Vglut2', 'PPN', 'SNc'): 10 mice ['C3096', 'C3412', 'C3798', 'C4707', 'C5109', 'C5112', 'C5367', 'C5368', 'C6454', 'C6455']
DA terminal mice recorded in both TS and DLS: 0
Vglut2 STN: exact two-sided Wilcoxon p of the paper's W (paper p 0.012, 0.008, 0.82):
W = 41: n 9 -> p 0.0117; n 10 -> p 0.0371
W = 43: n 9 -> p 0.0078; n 10 -> p 0.0273
W = -5: n 9 -> p 0.8203; n 10 -> p 0.8457
Result: not reproduced. 10 of 16 counts match. The SNr GABA, pan-DA and GRAB-DA groups match exactly, as do DA terminals in TS (9) and Vglut2 PPN (10). DA terminals in DLS have 8 mice (paper 9) and Vglut2 STN 9 (paper 10). The four Anxa1 counts differ: DLS 5 and 4 against 8 and 10, SNc 8 and 10 against 5 and 4. All four Anxa1 counts would match if the unsuffixed recordings were DLS recordings; the data check and the Anxa1 claims below favour the reading that the paper's DLS recordings of those nine mice are not in the dandiset. For Vglut2 STN the paper's own Wilcoxon statistics fit 9 mice: W = 41, 43 and -5 give exact p = 0.0117, 0.0078 and 0.82 with n = 9 (the paper's 0.012, 0.008 and 0.82) and 0.037, 0.027 and 0.85 with n = 10, so the legend's 10 is likely an error in the paper, not a missing mouse. Counted that way 11 of 16 match, still below 12. The DA terminal mice recorded in TS (9) and DLS (8) are different mice, so the data do not show which DLS mouse is absent.
fig2-gaba-excited¶
Optogenetic stimulation of STN or PPN inputs increases SNr GABA neuron activity (paired t-test on the 1 s train, STN t = 3.8, p = 0.008; PPN t = 9.2, p < 0.0001).
fig, axes = plt.subplots(1, 2, figsize=(10, 3), sharey=True)
for ax, inp in zip(axes, ["STN", "PPN"]):
lo, hi = plot_traces(ax, "GABA", inp, "SNr", duration=1.0, label=f"SNr GABA, {inp}")
ax.set_title(f"{inp}, 1 s 40 Hz"); ax.set_xlabel("s from onset"); ax.legend()
print(f"{inp}: plotted group mean ranges {lo:.2f} to {hi:.2f} % dF/F (all within the axes)")
axes[0].set_ylabel("dF/F (%)")
plt.show()
for inp, paper in [("STN", "t = 3.8, p = 0.008"), ("PPN", "t = 9.2, p < 0.0001")]:
for dur in [0.25, 1.0, 4.0]:
c = cohort("GABA", inp, "SNr", duration=dur)
res = ph.paired_t(c.during, c.pre)
tag = f"(paper, 1 s: {paper})" if dur == 1.0 else ""
print(f"{inp} {dur} s: n {res['n']}, mean during {c.during.mean():.2f}, pre {c.pre.mean():.2f} % ; t {res['stat']:.2f}, p {res['p']:.2g} {tag}")
c = cohort("GABA", "STN", "SNr", duration=1.0)
print("\nSTN per mouse, during 1 s train (%):", c.during.round(2).to_dict())
print("Variant, mice of the 24 kHz rig (10 trials) vs others, STN 1 s:", c.groupby(c.n_trials).during.agg(["size", "mean"]).round(2).to_dict())
STN: plotted group mean ranges -0.53 to 28.71 % dF/F (all within the axes) PPN: plotted group mean ranges -1.36 to 8.19 % dF/F (all within the axes)
STN 0.25 s: n 7, mean during 5.55, pre -0.01 % ; t 4.20, p 0.0057
STN 1.0 s: n 7, mean during 18.89, pre -0.00 % ; t 3.86, p 0.0083 (paper, 1 s: t = 3.8, p = 0.008)
STN 4.0 s: n 7, mean during 35.28, pre -0.08 % ; t 4.06, p 0.0066
PPN 0.25 s: n 8, mean during 3.12, pre 0.07 % ; t 10.72, p 1.4e-05
PPN 1.0 s: n 8, mean during 6.45, pre 0.00 % ; t 9.19, p 3.7e-05 (paper, 1 s: t = 9.2, p < 0.0001)
PPN 4.0 s: n 8, mean during 9.08, pre -0.29 % ; t 6.56, p 0.00031
STN per mouse, during 1 s train (%): {'B7532': 11.45, 'B7515': 9.68, 'B7514': 5.53, 'B7616': 9.75, 'B8569': 25.32, 'B8603': 31.67, 'B8677': 38.86}
Variant, mice of the 24 kHz rig (10 trials) vs others, STN 1 s: {'size': {10: 7}, 'mean': {10: 18.89}}
Result: reproduced. During 1 s 40 Hz trains SNr GABA activity rises above the pre-stimulus window for STN (mean 18.89 %, t = 3.86, p = 0.0083, n = 7; paper t = 3.8, p = 0.008) and PPN (6.45 %, t = 9.19, p = 3.7e-05, n = 8; paper t = 9.2). The 0.25 s and 4 s trains also give significant increases (p from 1.4e-05 to 0.0066). All seven STN mice were recorded with 10 trials per condition.
fig2-gaba-frequency¶
The SNr GABA response scales with stimulation frequency from 5 to 40 Hz (Friedman test on 1 s trains, STN F = 21.0, p = 0.0001; PPN F = 18.2, p = 0.0004).
for inp, paper in [("STN", "F = 21.0, p = 0.0001"), ("PPN", "F = 18.2, p = 0.0004")]:
m = freq_matrix("GABA", inp, 1.0, "SNr")
res = ph.friedman(m.values)
print(f"{inp}: mean during (%) by frequency {m.mean().round(2).to_dict()}; mice {len(m)}, complete {res['n']}")
print(f" Friedman chi2 {res['stat']:.2f}, p {res['p']:.2g} (paper: {paper}); mice with 40 Hz > 5 Hz: {int((m[40.0] > m[5.0]).sum())} of {len(m)}")
print(f" RM ANOVA variant: {fmt(ph.rm_anova(m.values))}")
for dur in [0.25, 4.0]:
mv = freq_matrix("GABA", inp, dur, "SNr")
if mv.dropna().shape[0] >= 3:
rv = ph.friedman(mv.values)
print(f" variant {dur} s trains: mice {rv['n']}, means {mv.mean().round(2).to_dict()}, Friedman {rv['stat']:.2f}, p {rv['p']:.2g}")
else:
print(f" variant {dur} s trains: {mv.dropna().shape[0]} mice with all four frequencies")
STN: mean during (%) by frequency {5.0: 4.07, 10.0: 9.6, 20.0: 22.61, 40.0: 32.66}; mice 7, complete 7
Friedman chi2 21.00, p 0.00011 (paper: F = 21.0, p = 0.0001); mice with 40 Hz > 5 Hz: 7 of 7
RM ANOVA variant: test RM ANOVA (GG), n 7, stat 11.5227, df1 1.0164, df2 6.0982, epsilon 0.3388, p 0.0141, p_uncorrected 0.0002
variant 0.25 s trains: mice 7, means {5.0: 1.78, 10.0: 3.04, 20.0: 6.36, 40.0: 9.39}, Friedman 21.00, p 0.00011
variant 4.0 s trains: 0 mice with all four frequencies
PPN: mean during (%) by frequency {5.0: 1.57, 10.0: 2.77, 20.0: 3.9, 40.0: 6.1}; mice 8, complete 8
Friedman chi2 18.15, p 0.00041 (paper: F = 18.2, p = 0.0004); mice with 40 Hz > 5 Hz: 8 of 8
RM ANOVA variant: test RM ANOVA (GG), n 8, stat 24.9185, df1 1.6251, df2 11.376, epsilon 0.5417, p 0.0001, p_uncorrected 0.0
variant 0.25 s trains: mice 8, means {5.0: 0.69, 10.0: 1.12, 20.0: 2.38, 40.0: 3.35}, Friedman 18.15, p 0.00041
variant 4.0 s trains: mice 7, means {5.0: 1.86, 10.0: 3.08, 20.0: 5.39, 40.0: 8.73}, Friedman 21.00, p 0.00011
Result: reproduced. For 1 s trains the response rises from 5 to 40 Hz for STN (4.07 to 32.66 %, Friedman 21.00, p = 0.00011, 40 Hz above 5 Hz in 7 of 7 mice) and PPN (1.57 to 6.10 %, 18.15, p = 0.00041, 8 of 8), the paper's values to the first decimal. The repeated-measures ANOVA variant gives p = 0.014 and 0.0001, and 0.25 s trains give the same Friedman values.
fig3b-da-stn-250ms-peak¶
A 250 ms 40 Hz train of STN inputs elicits a net increase in SNc pan-DA neuron activity that peaks after stimulus offset (latency to peak 0.50 +/- 0.03 s, mean +/- SEM).
c = cohort("panDA", "STN", "SNc", duration=0.25)
lat = c.peak_latency
print("Per mouse latency to peak (s, 0 to 2.25 s after onset):", lat.round(2).to_dict())
print(f"Mean {lat.mean():.3f} s, SEM {lat.sem():.3f} s, n {len(lat)} (paper 0.50 +/- 0.03 s); mice peaking after offset (> 0.25 s): {int((lat > 0.25).sum())}")
net, short = [], []
for r in recs:
if r["path"] in c.path.values:
ct = ph.condition_traces(r)
cc = ct[np.isclose(ct.duration, 0.25)].iloc[0]
net.append(ph.window_mean(cc.trace, cc.t, 0.0, 1.25))
w = (cc.t >= 0) & (cc.t < 1.25)
short.append(cc.t[w][np.nanargmax(cc.trace[w])])
net, short = np.array(net), np.array(short)
res = stats.ttest_1samp(net, 0)
print(f"Net response over 0 to 1.25 s: mean {net.mean():.2f} % (range {net.min():.2f} to {net.max():.2f}), t {res.statistic:.2f}, p {res.pvalue:.2g}")
print(f"Variant, peak window 0 to 1.25 s: mean latency {short.mean():.3f} s, SEM {short.std(ddof=1) / np.sqrt(len(short)):.3f}")
fig, ax = plt.subplots(figsize=(5, 3))
lo, hi = plot_traces(ax, "panDA", "STN", "SNc", duration=0.25, label="SNc pan-DA, STN 0.25 s")
ax.set_xlabel("s from onset"); ax.set_ylabel("dF/F (%)"); ax.legend(); plt.show()
print(f"plotted mean ranges {lo:.2f} to {hi:.2f} %")
# The mice whose peak lies after 1.25 s: their traces from -1 to 5 s.
late = lat[lat > 1.25].index.tolist()
fig, ax = plt.subplots(figsize=(6, 3))
for r in recs:
if r["path"] in c.path.values and r["subject"] in late:
ct = ph.condition_traces(r); cc = ct[np.isclose(ct.duration, 0.25)].iloc[0]
w = (cc.t >= -1) & (cc.t < 5)
ax.plot(cc.t[w], cc.trace[w], label=r["subject"])
e = (cc.t >= 0) & (cc.t < 1.25); l = (cc.t >= 1.25) & (cc.t < 2.25); x = (cc.t >= 0) & (cc.t < 5)
print(f"{r['subject']}: n trials {cc.n_trials}; max over 0 to 1.25 s {np.nanmax(cc.trace[e]):.2f} % at {cc.t[e][np.nanargmax(cc.trace[e])]:.2f} s; "
f"max over 1.25 to 2.25 s {np.nanmax(cc.trace[l]):.2f} % at {cc.t[l][np.nanargmax(cc.trace[l])]:.2f} s; "
f"max over 0 to 5 s at {cc.t[x][np.nanargmax(cc.trace[x])]:.2f} s; mean over 2 to 5 s {ph.window_mean(cc.trace, cc.t, 2, 5):.2f} %; "
f"range -1 to 5 s {np.nanmin(cc.trace[w]):.2f} to {np.nanmax(cc.trace[w]):.2f} %")
ax.axvspan(0, 0.25, color="red", alpha=0.08); ax.axvline(2.25, color="k", ls=":")
ax.set_xlabel("s from onset"); ax.set_ylabel("dF/F (%)"); ax.set_title("Mice peaking after 1.25 s, 0.25 s trains"); ax.legend(); plt.show()
Per mouse latency to peak (s, 0 to 2.25 s after onset): {'C4217': 2.03, 'C4218': 0.47, 'C4279': 0.46, 'C4280': 0.58, 'C4407': 0.57, 'B8370': 0.23, 'B8364': 0.51, 'B8369': 0.6, 'B8427': 0.47, 'B8429': 0.49, 'B8612': 2.24}
Mean 0.786 s, SEM 0.204 s, n 11 (paper 0.50 +/- 0.03 s); mice peaking after offset (> 0.25 s): 10
Net response over 0 to 1.25 s: mean 2.89 % (range -1.35 to 6.78), t 3.28, p 0.0083
Variant, peak window 0 to 1.25 s: mean latency 0.492 s, SEM 0.033
plotted mean ranges -0.44 to 5.68 % C4217: n trials 5; max over 0 to 1.25 s 1.28 % at 0.41 s; max over 1.25 to 2.25 s 1.39 % at 2.03 s; max over 0 to 5 s at 2.97 s; mean over 2 to 5 s 0.69 %; range -1 to 5 s -1.44 to 2.02 % B8612: n trials 10; max over 0 to 1.25 s 0.30 % at 0.62 s; max over 1.25 to 2.25 s 0.54 % at 2.24 s; max over 0 to 5 s at 2.78 s; mean over 2 to 5 s 0.55 %; range -1 to 5 s -0.92 to 1.15 %
Result: partially reproduced. With the fixed window (peak from onset to 2.25 s) the mean latency is 0.786 s (SEM 0.204, n = 11), outside 0.50 +/- 0.03 s, because two mice peak at the window's end (C4217 at 2.03 s and B8612 at 2.24 s). Ten of 11 mice peak after offset and the net response over 0 to 1.25 s is positive (2.89 %, t = 3.28, p = 0.0083). With a window ending 1.25 s after onset (a variant, which does not change the outcome) the mean latency is 0.492 s, SEM 0.033, the paper's value. The two late-peaking mice have small responses (ranges -1.44 to 2.02 % and -0.92 to 1.15 % from -1 to 5 s), an early local maximum at 0.41 and 0.62 s (1.28 and 0.30 %), and a slow rise that keeps going past the window (maximum over 0 to 5 s at 2.97 and 2.78 s). This looks like a late rebound or slow drift, not the early peak, and C4217 is one of the recordings whose clock could not be shown in the data check.
fig3c-da-stn-duration¶
For SNc pan-DA neurons the STN-evoked response varies with train duration; with 1 and 4 s trains an initial increase is followed by inhibition (Friedman test, F = 13.0, p = 0.005).
def duration_tests(group, inp, site):
rows = {}
for d in [0.25, 1.0, 4.0]:
c = cohort(group, inp, site, duration=d)
rows[d] = c
during = pd.DataFrame({d: rows[d].during for d in rows})
auc = pd.DataFrame({d: rows[d].auc for d in rows})
print("Mean during (%) per duration:", during.mean().round(2).to_dict(), "| mice:", len(during))
out = {}
out["a"] = ph.friedman(during.values)
out["b"] = ph.friedman(np.column_stack([rows[1.0].pre.reindex(during.index).values, during.values]))
out["c"] = ph.friedman(auc.values)
for k, v in out.items():
print(f" reading ({k}): Friedman chi2 {v['stat']:.2f}, df {3 if k == 'b' else 2}, p {v['p']:.2g}, n {v['n']}")
return rows, out
rows, out = duration_tests("panDA", "STN", "SNc")
early, late = [], []
for r in recs:
if r["path"] in rows[4.0].path.values:
ct = ph.condition_traces(r); cc = ct[np.isclose(ct.duration, 4.0) & np.isclose(ct.frequency, 40.0)].iloc[0]
early.append(ph.window_mean(cc.trace, cc.t, 0.0, 0.5)); late.append(ph.window_mean(cc.trace, cc.t, 3.0, 4.0))
early, late = np.array(early), np.array(late)
res = stats.ttest_rel(early, late)
print(f"4 s trains: mean over 0 to 0.5 s {early.mean():.2f} %, over 3 to 4 s {late.mean():.2f} %; paired t {res.statistic:.2f}, p {res.pvalue:.2g}, n {len(early)}")
print(f"Mice with late below early: {int((late < early).sum())} of {len(early)}; with late below 0: {int((late < 0).sum())}")
fig, ax = plt.subplots(figsize=(6, 3))
for d in [0.25, 1.0, 4.0]:
lo, hi = plot_traces(ax, "panDA", "STN", "SNc", duration=d, label=f"{d} s")
print(f"{d} s: plotted mean ranges {lo:.2f} to {hi:.2f} %")
ax.set_xlabel("s from onset"); ax.set_ylabel("dF/F (%)"); ax.legend(); plt.show()
Mean during (%) per duration: {0.25: 0.52, 1.0: 0.97, 4.0: -3.13} | mice: 11
reading (a): Friedman chi2 13.64, df 2, p 0.0011, n 11
reading (b): Friedman chi2 12.93, df 3, p 0.0048, n 11
reading (c): Friedman chi2 5.64, df 2, p 0.06, n 11
4 s trains: mean over 0 to 0.5 s 0.71 %, over 3 to 4 s -5.05 %; paired t 4.81, p 0.00071, n 11 Mice with late below early: 11 of 11; with late below 0: 11 0.25 s: plotted mean ranges -0.44 to 5.68 %
1.0 s: plotted mean ranges -0.27 to 5.56 %
4.0 s: plotted mean ranges -5.28 to 8.80 %
Result: reproduced. The STN-evoked pan-DA response differs between train durations (mean during 0.52, 0.97 and -3.13 % for 0.25, 1 and 4 s): reading (a), three durations, Friedman 13.64, p = 0.0011; reading (b), with the pre-stimulus window as a fourth group, 12.93, p = 0.0048, which matches the paper's F = 13.0, p = 0.005; reading (c), area under the curve, p = 0.06. In 4 s trains the mean over the first 0.5 s (0.71 %) is above that over the last 1 s (-5.05 %) in all 11 mice (paired t = 4.81, p = 0.00071).
fig3d-da-stn-frequency-ns¶
The amplitude of the SNc pan-DA response during STN stimulation does not scale with frequency (Friedman test, F = 4.2, p = 0.24).
for dur in [1.0, 0.25, 4.0]:
m = freq_matrix("panDA", "STN", dur, "SNc")
if m.dropna().shape[0] < 3:
print(f"{dur} s: {m.dropna().shape[0]} complete mice"); continue
res = ph.friedman(m.values)
tag = "main (paper F = 4.2, p = 0.24)" if dur == 1.0 else "variant"
print(f"{dur} s [{tag}]: means {m.mean().round(2).to_dict()}, Friedman chi2 {res['stat']:.2f}, p {res['p']:.2g}, n {res['n']}; "
f"40 minus 5 Hz {np.nanmean(m[40.0] - m[5.0]):.2f} %")
1.0 s [main (paper F = 4.2, p = 0.24)]: means {5.0: 2.26, 10.0: 1.57, 20.0: 1.13, 40.0: 2.47}, Friedman chi2 4.20, p 0.24, n 11; 40 minus 5 Hz 0.21 %
0.25 s [variant]: means {5.0: 0.46, 10.0: 0.57, 20.0: 0.35, 40.0: 0.93}, Friedman chi2 3.55, p 0.31, n 11; 40 minus 5 Hz 0.47 %
4.0 s [variant]: means {5.0: 0.39, 10.0: -1.06, 20.0: -0.5, 40.0: 0.39}, Friedman chi2 4.92, p 0.18, n 5; 40 minus 5 Hz -0.00 %
Result: reproduced. The paper's test finds no effect of frequency on the STN-evoked pan-DA response to 1 s trains: Friedman 4.20, p = 0.24 (paper F = 4.2, p = 0.24), 40 minus 5 Hz 0.21 %. The 0.25 s (p = 0.31) and 4 s (p = 0.18, 5 complete mice) variants also find none.
fig3g-da-ppn-duration¶
PPN stimulation also evokes SNc pan-DA responses that vary with train duration (Friedman test, F = 11.5, p = 0.009).
rows_p, out_p = duration_tests("panDA", "PPN", "SNc")
fig, ax = plt.subplots(figsize=(6, 3))
for d in [0.25, 1.0, 4.0]:
lo, hi = plot_traces(ax, "panDA", "PPN", "SNc", duration=d, label=f"{d} s")
print(f"{d} s: plotted mean ranges {lo:.2f} to {hi:.2f} %")
ax.set_xlabel("s from onset"); ax.set_ylabel("dF/F (%)"); ax.legend(); plt.show()
unclear = clock[(clock.ratio <= 1) & (clock.group == "panDA") & (clock.protocol == "durations")]
print("pan-DA durations recordings whose clock could not be shown from the response:", unclear.groupby("input").subject.apply(sorted).to_dict())
ppn_bad = set(unclear[unclear.input == "PPN"].subject)
during_p = pd.DataFrame({d: rows_p[d].during for d in rows_p})
keep = during_p.drop(index=[s for s in during_p.index if s in ppn_bad])
print(f"Variant without those PPN mice: n {len(keep)}, means {keep.mean().round(2).to_dict()}")
if len(keep) >= 3:
rk = ph.friedman(keep.values)
pre_k = rows_p[1.0].pre.reindex(keep.index).values
rb = ph.friedman(np.column_stack([pre_k, keep.values]))
print(f" reading (a) Friedman chi2 {rk['stat']:.2f}, p {rk['p']:.2g}; reading (b) chi2 {rb['stat']:.2f}, p {rb['p']:.2g}")
Mean during (%) per duration: {0.25: -0.6, 1.0: 1.03, 4.0: 0.7} | mice: 9
reading (a): Friedman chi2 2.67, df 2, p 0.26, n 9
reading (b): Friedman chi2 4.07, df 3, p 0.25, n 9
reading (c): Friedman chi2 1.56, df 2, p 0.46, n 9
0.25 s: plotted mean ranges -1.57 to 2.49 % 1.0 s: plotted mean ranges -1.60 to 2.14 %
4.0 s: plotted mean ranges -0.51 to 2.04 %
pan-DA durations recordings whose clock could not be shown from the response: {'PPN': ['C3164', 'C5007', 'C5008'], 'STN': ['C4217', 'C4218', 'C4279', 'C4407']}
Variant without those PPN mice: n 6, means {0.25: -0.73, 1.0: 1.63, 4.0: 1.08}
reading (a) Friedman chi2 4.33, p 0.11; reading (b) chi2 5.40, p 0.14
Result: not reproduced. For PPN stimulation the pan-DA response does not differ between train durations under any reading (means -0.60, 1.03 and 0.70 %; Friedman p = 0.26, 0.25 and 0.46 for readings (a) to (c), n = 9), where the paper reports F = 11.5, p = 0.009. Three of the nine PPN durations recordings (C3164, C5007, C5008) are among those whose clock could not be shown from the response; without them (n = 6, a variant) the test still finds no effect (reading (a) p = 0.11, reading (b) p = 0.14), so the failure does not come from those recordings alone.
fig3h-da-ppn-frequency-ns¶
The SNc pan-DA response to PPN stimulation does not vary with frequency (repeated-measures one-way ANOVA, F2,15 = 1.35, p = 0.28).
for dur in [1.0, 0.25, 4.0]:
m = freq_matrix("panDA", "PPN", dur, "SNc")
if m.dropna().shape[0] < 3:
print(f"{dur} s: {m.dropna().shape[0]} complete mice"); continue
res = ph.rm_anova(m.values)
tag = "main (paper F2,15 = 1.35, p = 0.28)" if dur == 1.0 else "variant"
print(f"{dur} s [{tag}]: means {m.mean().round(2).to_dict()}; {fmt(res)}; 40 minus 5 Hz {np.nanmean(m[40.0] - m[5.0]):.2f} %")
1.0 s [main (paper F2,15 = 1.35, p = 0.28)]: means {5.0: -1.27, 10.0: -1.5, 20.0: -1.38, 40.0: 0.03}; test RM ANOVA (GG), n 9, stat 1.3506, df1 1.849, df2 14.7924, epsilon 0.6163, p 0.2872, p_uncorrected 0.2816; 40 minus 5 Hz 1.30 %
0.25 s [variant]: means {5.0: -0.5, 10.0: -0.58, 20.0: -0.83, 40.0: -0.74}; test RM ANOVA (GG), n 9, stat 0.1266, df1 1.5701, df2 12.5608, epsilon 0.5234, p 0.8346, p_uncorrected 0.9434; 40 minus 5 Hz -0.24 %
4.0 s [variant]: means {5.0: -2.78, 10.0: -2.93, 20.0: -2.39, 40.0: -1.41}; test RM ANOVA (GG), n 9, stat 1.217, df1 2.5573, df2 20.4582, epsilon 0.8524, p 0.3247, p_uncorrected 0.325; 40 minus 5 Hz 1.38 %
Result: reproduced. The paper's test finds no effect of frequency on the PPN-evoked pan-DA response to 1 s trains: F = 1.35 on Greenhouse-Geisser corrected df 1.85 and 14.8, p = 0.29 (paper F2,15 = 1.35, p = 0.28); 40 minus 5 Hz is 1.30 %. The 0.25 s and 4 s variants give p = 0.83 and 0.32.
fig4cf-ts-dls-opposite¶
STN or PPN stimulation evokes opposite DA release in TS (increase) and DLS (decrease) (two-way ANOVA, striatal site, STN F1,11 = 40.9, PPN F1,16 = 36.4, p < 0.0001; site by time interaction of the same size).
def site_contrast(group, inp, duration=1.0):
ts, dls = cohort(group, inp, "TS", duration=duration), cohort(group, inp, "DLS", duration=duration)
res = ph.mixed_2x2(ts.pre, ts.during, dls.pre, dls.during)
print(f"{group} {inp} {duration} s: TS n {len(ts)} mean during {ts.during.mean():.2f} % (mice above 0: {int((ts.during > 0).sum())}); "
f"DLS n {len(dls)} mean during {dls.during.mean():.2f} % (mice below 0: {int((dls.during < 0).sum())})")
print(f" site F{res['df']} {res['F_site']:.1f}, p {res['p_site']:.2g}; interaction F {res['F_interaction']:.1f}, p {res['p_interaction']:.2g}")
return ts, dls, res
for inp, paper in [("STN", "F1,11 = 40.9"), ("PPN", "F1,16 = 36.4")]:
print(f"paper: {paper}")
site_contrast("GRABDA", inp)
for d in [0.25, 4.0]:
site_contrast("GRABDA", inp, d)
fig, axes = plt.subplots(1, 2, figsize=(10, 3), sharey=True)
for ax, inp in zip(axes, ["STN", "PPN"]):
for s in ["TS", "DLS"]:
lo, hi = plot_traces(ax, "GRABDA", inp, s, duration=1.0)
print(f"{inp} {s}: plotted mean ranges {lo:.2f} to {hi:.2f} %")
ax.set_title(f"GRAB-DA, {inp}, 1 s 40 Hz"); ax.set_xlabel("s from onset"); ax.legend()
axes[0].set_ylabel("dF/F (%)"); plt.show()
paper: F1,11 = 40.9 GRABDA STN 1.0 s: TS n 6 mean during 12.45 % (mice above 0: 6); DLS n 7 mean during -8.91 % (mice below 0: 7) site F(1, 11) 40.9, p 5.1e-05; interaction F 40.9, p 5.1e-05 GRABDA STN 0.25 s: TS n 6 mean during 5.37 % (mice above 0: 6); DLS n 7 mean during -4.28 % (mice below 0: 7) site F(1, 11) 25.4, p 0.00038; interaction F 28.7, p 0.00023 GRABDA STN 4.0 s: TS n 6 mean during 8.74 % (mice above 0: 6); DLS n 7 mean during -12.55 % (mice below 0: 7) site F(1, 11) 54.2, p 1.4e-05; interaction F 72.7, p 3.5e-06 paper: F1,16 = 36.4 GRABDA PPN 1.0 s: TS n 9 mean during 11.39 % (mice above 0: 9); DLS n 9 mean during -5.75 % (mice below 0: 9) site F(1, 16) 36.4, p 1.7e-05; interaction F 36.4, p 1.7e-05 GRABDA PPN 0.25 s: TS n 9 mean during 6.04 % (mice above 0: 9); DLS n 9 mean during -3.19 % (mice below 0: 9) site F(1, 16) 27.7, p 7.7e-05; interaction F 24.4, p 0.00015 GRABDA PPN 4.0 s: TS n 9 mean during 5.49 % (mice above 0: 7); DLS n 9 mean during -8.05 % (mice below 0: 9) site F(1, 16) 27.7, p 7.8e-05; interaction F 36.0, p 1.9e-05
STN TS: plotted mean ranges -3.28 to 16.29 % STN DLS: plotted mean ranges -13.22 to 5.40 %
PPN TS: plotted mean ranges -1.27 to 13.92 %
PPN DLS: plotted mean ranges -6.95 to 10.33 %
Result: reproduced. For 1 s 40 Hz trains GRAB-DA rises in TS and falls in DLS with both inputs: STN, TS 12.45 % (6 of 6 mice above 0), DLS -8.91 % (7 of 7 below 0), site by time interaction F(1, 11) = 40.9, p = 5.1e-05; PPN, TS 11.39 % (9 of 9), DLS -5.75 % (9 of 9), F(1, 16) = 36.4, p = 1.7e-05. These are the paper's values. The 0.25 s and 4 s trains give the same pattern (interaction p from 3.5e-06 to 0.00023). For 1 s trains the pre-stimulus window is the baseline, so the site F and the interaction F are equal by construction, as in the paper.
fig4dg-ts-frequency¶
In TS, the DA increase evoked by either input scales with frequency, for all train durations (STN, repeated-measures ANOVA, 0.25 s F2,8 = 7.5, 1 s F2,7 = 13.6, 4 s F2,9 = 16.9; PPN, Friedman, F = 20.0, 24.0, 23.2).
def frequency_tests(group, site, plan, direction, kind=None):
"""plan: list of (input, duration, test name, paper text). direction +1: 40 Hz above 5 Hz."""
out = []
for inp, dur, test, paper in plan:
m = freq_matrix(group, inp, dur, site, kind=kind)
res = ph.friedman(m.values) if test == "Friedman" else ph.rm_anova(m.values)
other = ph.rm_anova(m.values) if test == "Friedman" else ph.friedman(m.values)
diff = np.nanmean(m[40.0] - m[5.0])
ok = res["p"] < 0.05 and np.sign(diff) == direction
incomplete = m[~np.isfinite(m.values).all(axis=1)]
dropped = "; ".join(f"{s} lacks {'/'.join(f'{f:g}' for f in row.index[row.isna()])} Hz" for s, row in incomplete.iterrows())
out.append({"input": inp, "duration": dur, "test": test, "mice": len(m), "n": res["n"], "stat": round(res["stat"], 2),
"p": res["p"], "mean_5": m[5.0].mean(), "mean_10": m[10.0].mean(), "mean_20": m[20.0].mean(),
"mean_40": m[40.0].mean(), "holds": ok, "paper": paper, "other_test_p": other["p"], "max_friedman": 3 * res["n"], "dropped": dropped or "-"})
t = pd.DataFrame(out)
print(t.round(4).to_string(index=False))
print("Tests that hold:", int(t.holds.sum()), "of", len(t))
return t
plan = [("STN", 0.25, "RM", "F2,8 = 7.5, p = 0.02"), ("STN", 1.0, "RM", "F2,7 = 13.6, p = 0.004"), ("STN", 4.0, "RM", "F2,9 = 16.9, p = 0.001"),
("PPN", 0.25, "Friedman", "F = 20.0, p = 0.0002"), ("PPN", 1.0, "Friedman", "F = 24.0, p < 0.0001"), ("PPN", 4.0, "Friedman", "F = 23.2, p < 0.0001")]
ts_freq = frequency_tests("GRABDA", "TS", plan, +1)
input duration test mice n stat p mean_5 mean_10 mean_20 mean_40 holds paper other_test_p max_friedman dropped STN 0.25 RM 6 6 7.62 0.0191 2.6026 3.8561 5.3315 7.1496 True F2,8 = 7.5, p = 0.02 0.0081 18 - STN 1.00 RM 6 6 13.76 0.0044 3.3894 7.1187 9.4613 12.8385 True F2,7 = 13.6, p = 0.004 0.0018 18 - STN 4.00 RM 6 6 17.31 0.0011 2.0716 6.0386 8.5342 9.3275 True F2,9 = 16.9, p = 0.001 0.0010 18 - PPN 0.25 Friedman 9 9 19.80 0.0002 3.0627 3.3258 5.5891 9.2448 True F = 20.0, p = 0.0002 0.0070 27 - PPN 1.00 Friedman 9 9 24.33 0.0000 2.7028 3.0110 8.0197 12.9861 True F = 24.0, p < 0.0001 0.0014 27 - PPN 4.00 Friedman 9 9 23.27 0.0000 0.6376 0.6997 4.9207 7.9162 True F = 23.2, p < 0.0001 0.0034 27 - Tests that hold: 6 of 6
Result: reproduced. In TS all six tests find an effect of frequency with the 40 Hz mean above the 5 Hz mean: STN repeated-measures ANOVA F = 7.62, 13.76 and 17.31 (p = 0.019, 0.0044, 0.0011; paper 7.5, 13.6, 16.9) and PPN Friedman 19.80, 24.33 and 23.27 (p = 0.0002 or below; paper 20.0, 24.0, 23.2).
fig4eh-dls-frequency¶
In DLS, the DA inhibition scales with frequency and is most prominent at higher frequencies (STN, repeated-measures ANOVA, 0.25 s F = 7.5, 1 s F = 12.8, 4 s F = 46.8; PPN, 0.25 s Friedman F = 12.3, 1 s and 4 s repeated-measures ANOVA F = 12.8 and 10.2).
plan = [("STN", 0.25, "RM", "F2,9 = 7.5, p = 0.015"), ("STN", 1.0, "RM", "F1,11 = 12.8, p = 0.0005"), ("STN", 4.0, "RM", "F2,13 = 46.8, p < 0.0001"),
("PPN", 0.25, "Friedman", "F = 12.3, p = 0.006"), ("PPN", 1.0, "RM", "F1,11 = 12.8, p = 0.002"), ("PPN", 4.0, "RM", "F2,13 = 10.2, p = 0.004")]
dls_freq = frequency_tests("GRABDA", "DLS", plan, -1)
print("\nThe paper gives F1,11 = 12.8 for both STN and PPN 1 s trains, with p = 0.0005 and 0.002; F(1, 11) = 12.8 gives p =",
round(stats.f.sf(12.8, 1, 11), 4))
input duration test mice n stat p mean_5 mean_10 mean_20 mean_40 holds paper other_test_p max_friedman dropped STN 0.25 RM 7 7 7.89 0.0131 -0.3924 -2.4723 -3.0830 -2.9360 True F2,9 = 7.5, p = 0.015 0.0106 21 - STN 1.00 RM 7 7 26.17 0.0000 1.5052 -0.6185 -5.3302 -5.7811 True F1,11 = 12.8, p = 0.0005 0.0005 21 - STN 4.00 RM 7 7 47.23 0.0000 5.4077 -1.6904 -6.5861 -8.2738 True F2,13 = 46.8, p < 0.0001 0.0003 21 - PPN 0.25 Friedman 9 9 11.53 0.0092 0.1900 -1.9987 -2.4408 -1.4803 True F = 12.3, p = 0.006 0.0131 27 - PPN 1.00 RM 9 9 12.98 0.0023 1.4044 -1.0476 -3.7206 -4.7846 True F1,11 = 12.8, p = 0.002 0.0010 27 - PPN 4.00 RM 9 9 15.92 0.0004 2.0076 -2.0989 -5.0765 -5.5154 True F2,13 = 10.2, p = 0.004 0.0003 27 - Tests that hold: 6 of 6 The paper gives F1,11 = 12.8 for both STN and PPN 1 s trains, with p = 0.0005 and 0.002; F(1, 11) = 12.8 gives p = 0.0043
Result: reproduced. In DLS all six tests find an effect of frequency with the 40 Hz mean below the 5 Hz mean (p from below 0.0001 to 0.013). STN 0.25 s and 4 s give F = 7.89 and 47.23 (paper 7.5 and 46.8), PPN 0.25 s Friedman 11.53 (paper 12.3) and PPN 1 s F = 12.98 (paper 12.8). STN 1 s gives F = 26.17, where the paper reports F1,11 = 12.8 with p = 0.0005, the same F and df as it reports for PPN 1 s with p = 0.002; F(1, 11) = 12.8 gives p = 0.0043, so one of the two entries in the paper is likely a copy.
fig4eh-dls-5hz-increase¶
Low-frequency (5 Hz) stimulation of either STN or PPN increases DLS dopamine.
for inp in ["STN", "PPN"]:
for dur in [1.0, 0.25, 4.0]:
c = cohort("GRABDA", inp, "DLS", protocol="frequencies", duration=dur, frequency=5.0)
if not len(c):
continue
res = stats.ttest_1samp(c.during, 0)
tag = {1.0: "main, Figure 4E/4H", 4.0: "variant, Figure S3F/S3L, cited for this sentence", 0.25: "variant, Figure S3E/S3K, not cited"}[dur]
print(f"{inp} {dur} s 5 Hz [{tag}]: n {len(c)}, mean during {c.during.mean():.2f} % (range {c.during.min():.2f} to {c.during.max():.2f}), "
f"mice above 0 {int((c.during > 0).sum())}, t {res.statistic:.2f}, p {res.pvalue:.2g}; mean in 1 s after offset {c.post.mean():.2f} %")
four = {inp: stats.ttest_1samp(cohort("GRABDA", inp, "DLS", protocol="frequencies", duration=4.0, frequency=5.0).during, 0) for inp in ["STN", "PPN"]}
print("Reading by the 4 s panels the sentence also cites: inputs with a mean above 0 at p < 0.05:", sum(r.pvalue < 0.05 and r.statistic > 0 for r in four.values()), "of 2")
STN 1.0 s 5 Hz [main, Figure 4E/4H]: n 7, mean during 1.51 % (range -2.65 to 5.08), mice above 0 4, t 1.39, p 0.22; mean in 1 s after offset 4.29 % STN 0.25 s 5 Hz [variant, Figure S3E/S3K, not cited]: n 7, mean during -0.39 % (range -2.03 to 1.16), mice above 0 4, t -0.74, p 0.49; mean in 1 s after offset 3.41 % STN 4.0 s 5 Hz [variant, Figure S3F/S3L, cited for this sentence]: n 7, mean during 5.41 % (range 2.90 to 8.59), mice above 0 7, t 6.52, p 0.00062; mean in 1 s after offset 6.76 % PPN 1.0 s 5 Hz [main, Figure 4E/4H]: n 9, mean during 1.40 % (range -2.16 to 5.34), mice above 0 6, t 1.67, p 0.13; mean in 1 s after offset 0.96 % PPN 0.25 s 5 Hz [variant, Figure S3E/S3K, not cited]: n 9, mean during 0.19 % (range -4.10 to 1.48), mice above 0 6, t 0.32, p 0.76; mean in 1 s after offset 2.65 % PPN 4.0 s 5 Hz [variant, Figure S3F/S3L, cited for this sentence]: n 9, mean during 2.01 % (range -2.25 to 5.69), mice above 0 6, t 2.09, p 0.07; mean in 1 s after offset 4.00 % Reading by the 4 s panels the sentence also cites: inputs with a mean above 0 at p < 0.05: 1 of 2
Result: not reproduced. For 1 s trains at 5 Hz (Figure 4E and 4H, the reading fixed in the criterion) the mean GRAB-DA signal in DLS during the train is positive but not significantly so for either input: STN 1.51 % (4 of 7 mice above 0, p = 0.22), PPN 1.40 % (6 of 9, p = 0.13). The sentence also cites Figure S3F and S3L, the 4 s DLS frequency panels. In 4 s trains at 5 Hz the signal is above 0 for STN (5.41 %, 7 of 7, p = 0.00062) and not significantly for PPN (2.01 %, 6 of 9, p = 0.07), 1 of 2 inputs, which would be partially reproduced under that reading; it does not change the outcome. The 0.25 s trains (Figure S3E and S3K, not cited) show no increase (p = 0.49 and 0.76). The 1 s after offset is positive in every condition (0.96 to 6.76 %), as the paper's account of summed rebounds would predict.
fig4cf-dls-rebound¶
In DLS the DA inhibition is followed by a post-stimulus rebound above baseline.
for inp in ["STN", "PPN"]:
for dur in [1.0, 0.25, 4.0]:
c = cohort("GRABDA", inp, "DLS", duration=dur)
res = stats.ttest_1samp(c.post, 0)
tag = "main" if dur == 1.0 else "variant"
print(f"{inp} {dur} s [{tag}]: n {len(c)}, mean during {c.during.mean():.2f} %, mean in 1 s after offset {c.post.mean():.2f} % "
f"(range {c.post.min():.2f} to {c.post.max():.2f}), mice above 0 {int((c.post > 0).sum())}, t {res.statistic:.2f}, p {res.pvalue:.2g}")
STN 1.0 s [main]: n 7, mean during -8.91 %, mean in 1 s after offset 0.87 % (range -2.04 to 3.93), mice above 0 5, t 1.17, p 0.29 STN 0.25 s [variant]: n 7, mean during -4.28 %, mean in 1 s after offset 2.75 % (range 1.09 to 5.58), mice above 0 7, t 4.19, p 0.0058 STN 4.0 s [variant]: n 7, mean during -12.55 %, mean in 1 s after offset 4.09 % (range -0.17 to 7.51), mice above 0 6, t 3.61, p 0.011 PPN 1.0 s [main]: n 9, mean during -5.75 %, mean in 1 s after offset 4.53 % (range -1.71 to 12.29), mice above 0 8, t 2.99, p 0.017 PPN 0.25 s [variant]: n 9, mean during -3.19 %, mean in 1 s after offset 4.36 % (range 0.66 to 8.76), mice above 0 9, t 4.00, p 0.0039 PPN 4.0 s [variant]: n 9, mean during -8.05 %, mean in 1 s after offset 3.56 % (range -2.49 to 7.48), mice above 0 7, t 3.10, p 0.015
Result: partially reproduced. After 1 s 40 Hz trains GRAB-DA in DLS rises above baseline for PPN (4.53 %, 8 of 9 mice, p = 0.017) but not significantly for STN (0.87 %, 5 of 7, p = 0.29). After 0.25 s and 4 s trains the rebound is significant for both inputs (p from 0.0039 to 0.015).
figs3-terminals-opposite¶
Calcium in DA terminals in TS and DLS shows opposite responses to STN stimulation, as the GRAB-DA recordings do.
ts_t, dls_t, res_t = site_contrast("DAterm", "STN")
for d in [0.25, 4.0]:
site_contrast("DAterm", "STN", d)
print(f"DLS terminals, mean in 1 s after offset of 1 s trains: {dls_t.post.mean():.2f} % (mice above 0: {int((dls_t.post > 0).sum())} of {len(dls_t)})")
fig, ax = plt.subplots(figsize=(5, 3))
for s in ["TS", "DLS"]:
lo, hi = plot_traces(ax, "DAterm", "STN", s, duration=1.0)
print(f"{s}: plotted mean ranges {lo:.2f} to {hi:.2f} %")
ax.set_title("DA terminal GCaMP, STN, 1 s 40 Hz"); ax.set_xlabel("s from onset"); ax.set_ylabel("dF/F (%)"); ax.legend(); plt.show()
DAterm STN 1.0 s: TS n 9 mean during 8.96 % (mice above 0: 9); DLS n 8 mean during -4.41 % (mice below 0: 7) site F(1, 15) 26.2, p 0.00013; interaction F 26.2, p 0.00013 DAterm STN 0.25 s: TS n 9 mean during 5.60 % (mice above 0: 9); DLS n 8 mean during -3.50 % (mice below 0: 8) site F(1, 15) 23.7, p 0.0002; interaction F 29.6, p 6.9e-05 DAterm STN 4.0 s: TS n 9 mean during 2.73 % (mice above 0: 9); DLS n 8 mean during -2.77 % (mice below 0: 7) site F(1, 15) 12.2, p 0.0032; interaction F 25.0, p 0.00016 DLS terminals, mean in 1 s after offset of 1 s trains: 6.73 % (mice above 0: 8 of 8)
TS: plotted mean ranges -1.83 to 13.97 %
DLS: plotted mean ranges -5.44 to 16.81 %
Result: reproduced. STN stimulation raises DA terminal calcium in TS (8.96 %, 9 of 9 mice above 0) and lowers it in DLS (-4.41 %, 7 of 8 below 0) during 1 s trains; site by time interaction F(1, 15) = 26.2, p = 0.00013. The 0.25 s and 4 s trains give the same pattern, and the DLS terminals rebound after offset (6.73 %, 8 of 8 above 0).
fig5-vglut2-excited¶
Stimulation of STN or PPN inputs increases the activity of Vglut2+ DA neurons in SNc (STN, Wilcoxon, 0.25 s W = 41, p = 0.012, 1 s W = 43, p = 0.008, 4 s W = -5, p = 0.82; PPN, paired t-test, 0.25 s t = 2.9, p = 0.02, 1 s t = 3.4, p = 0.008, 4 s t = 3.0, p = 0.01).
def pre_during_tests(group, site, plan, kind=None, label=""):
out = []
for inp, dur, test, paper in plan:
c = cohort(group, inp, site, duration=dur, kind=kind)
res = ph.wilcoxon_w(c.during, c.pre) if test == "Wilcoxon" else ph.paired_t(c.during, c.pre)
other = ph.paired_t(c.during, c.pre) if test == "Wilcoxon" else ph.wilcoxon_w(c.during, c.pre)
out.append({"input": inp, "duration": dur, "test": res["test"], "n": res["n"], "stat": round(res["stat"], 2), "p": res["p"],
"mean_during": c.during.mean(), "mean_pre": c.pre.mean(), "mice_up": int((c.during > c.pre).sum()),
"mean_post": c.post.mean(), "paper": paper, "other_test_p": other["p"]})
t = pd.DataFrame(out)
print(label); print(t.round(4).to_string(index=False))
return t
plan = [("STN", 0.25, "Wilcoxon", "W = 41, p = 0.012"), ("STN", 1.0, "Wilcoxon", "W = 43, p = 0.008"), ("STN", 4.0, "Wilcoxon", "W = -5, p = 0.82"),
("PPN", 0.25, "t", "t = 2.9, p = 0.02"), ("PPN", 1.0, "t", "t = 3.4, p = 0.008"), ("PPN", 4.0, "t", "t = 3.0, p = 0.01")]
vg = pre_during_tests("Vglut2", "SNc", plan)
core = vg[vg.duration < 4]
print("\nTests of 0.25 and 1 s trains with a significant increase:", int(((core.p < 0.05) & (core.mean_during > core.mean_pre)).sum()), "of 4")
fig, axes = plt.subplots(1, 2, figsize=(10, 3), sharey=True)
for ax, inp in zip(axes, ["STN", "PPN"]):
lo, hi = plot_traces(ax, "Vglut2", inp, "SNc", duration=1.0, label="Vglut2 SNc")
ax.set_title(f"{inp}, 1 s 40 Hz"); ax.legend(); ax.set_xlabel("s from onset")
print(f"{inp}: plotted mean ranges {lo:.2f} to {hi:.2f} %")
axes[0].set_ylabel("dF/F (%)"); plt.show()
input duration test n stat p mean_during mean_pre mice_up mean_post paper other_test_p STN 0.25 Wilcoxon 9 43.00 0.0078 3.6763 -0.2581 8 2.5570 W = 41, p = 0.012 0.0157 STN 1.00 Wilcoxon 9 43.00 0.0078 6.5357 0.0000 8 1.6618 W = 43, p = 0.008 0.0209 STN 4.00 Wilcoxon 9 -5.00 0.8203 0.8797 0.5611 3 -1.2967 W = -5, p = 0.82 0.6995 PPN 0.25 paired t 10 2.89 0.0178 3.0364 -0.0900 8 3.4144 t = 2.9, p = 0.02 0.0098 PPN 1.00 paired t 10 3.42 0.0076 9.3509 -0.0000 9 5.4061 t = 3.4, p = 0.008 0.0098 PPN 4.00 paired t 10 2.56 0.0305 6.3781 -0.3895 9 1.5030 t = 3.0, p = 0.01 0.0273 Tests of 0.25 and 1 s trains with a significant increase: 4 of 4
STN: plotted mean ranges -0.86 to 8.47 % PPN: plotted mean ranges -2.07 to 12.22 %
Result: reproduced. Vglut2+ DA neurons increase their activity during 0.25 s and 1 s trains with both inputs: STN Wilcoxon W = 43 and 43 (p = 0.0078 each, n = 9; paper 41 and 43), PPN paired t = 2.89 and 3.42 (p = 0.018 and 0.0076, n = 10; paper 2.9 and 3.4). For 4 s trains STN gives W = -5, p = 0.82 (the paper's values) and PPN t = 2.56, p = 0.031 (paper t = 3.0, p = 0.01).
fig5-vglut2-frequency¶
The Vglut2+ DA neuron increase scales with frequency for short STN trains and for all PPN durations (Friedman, STN 0.25 s F = 15.5, p = 0.001, 1 s F = 5.9, p = 0.11, 4 s F = 7.0, p = 0.07; PPN F = 18.8, 24.8, 20.2, p <= 0.0003).
plan = [("STN", 0.25, "Friedman", "F = 15.5, p = 0.001"), ("STN", 1.0, "Friedman", "F = 5.9, p = 0.11"), ("STN", 4.0, "Friedman", "F = 7.0, p = 0.07"),
("PPN", 0.25, "Friedman", "F = 18.8, p = 0.0003"), ("PPN", 1.0, "Friedman", "F = 24.8, p < 0.0001"), ("PPN", 4.0, "Friedman", "F = 20.2, p = 0.0001")]
vg_freq = frequency_tests("Vglut2", "SNc", plan, +1)
core = vg_freq[~((vg_freq.input == "STN") & (vg_freq.duration > 0.25))]
print("Tests the claim rests on (STN 0.25 s, PPN all) that hold:", int(core.holds.sum()), "of", len(core))
print("The Friedman statistic over four conditions is at most 3n; the paper's PPN 1 s value 24.8 needs n >=", int(np.ceil(24.8 / 3)),
"and 18.8 and 20.2 need n >=", int(np.ceil(18.8 / 3)), "and", int(np.ceil(20.2 / 3)))
input duration test mice n stat p mean_5 mean_10 mean_20 mean_40 holds paper other_test_p max_friedman dropped STN 0.25 Friedman 9 9 13.53 0.0036 1.6972 1.4616 2.7505 3.4994 True F = 15.5, p = 0.001 0.0962 27 - STN 1.00 Friedman 9 9 2.73 0.4346 3.2617 4.0913 3.2141 5.2406 False F = 5.9, p = 0.11 0.2114 27 - STN 4.00 Friedman 9 9 5.13 0.1623 1.9980 -0.0191 0.5809 -0.0084 False F = 7.0, p = 0.07 0.4047 27 - PPN 0.25 Friedman 10 7 11.23 0.0106 0.2882 1.0644 2.8083 3.9936 True F = 18.8, p = 0.0003 0.0316 21 C3412 lacks 5 Hz; C3798 lacks 5 Hz; C4707 lacks 5 Hz PPN 1.00 Friedman 10 7 17.91 0.0005 -0.5696 1.5787 5.0631 10.5561 True F = 24.8, p < 0.0001 0.0037 21 C3412 lacks 5 Hz; C3798 lacks 5 Hz; C4707 lacks 5 Hz PPN 4.00 Friedman 10 7 16.20 0.0010 -2.2717 -1.0520 2.8456 7.1468 True F = 20.2, p = 0.0001 0.0085 21 C3412 lacks 5 Hz; C3798 lacks 5 Hz; C4707 lacks 5 Hz Tests that hold: 4 of 6 Tests the claim rests on (STN 0.25 s, PPN all) that hold: 4 of 4 The Friedman statistic over four conditions is at most 3n; the paper's PPN 1 s value 24.8 needs n >= 9 and 18.8 and 20.2 need n >= 7 and 7
Result: reproduced. The four tests the claim rests on find an effect of frequency with the 40 Hz mean above the 5 Hz mean: STN 0.25 s (Friedman 13.53, p = 0.0036) and PPN 0.25, 1 and 4 s (11.23, 17.91 and 16.20; p = 0.011, 0.0005 and 0.0010). The PPN tests use 7 of 10 mice: C3412, C3798 and C4707 have no 5 Hz trains in the dandiset. The Friedman statistic over four conditions is at most 3n, and the paper's 24.8 for PPN 1 s needs at least 9 mice, so the paper had 5 Hz trains for mice whose 5 Hz trains are not deposited. As in the paper, STN 1 s and 4 s are not significant (p = 0.43 and 0.16; paper 0.11 and 0.07).
fig5-anxa1-dls-inhibited¶
In DLS terminals of SNc Anxa1+ neurons, STN and PPN stimulation produce inhibition (STN, 0.25 s t = 6.4, p = 0.0004, 1 s W = -36, p = 0.005, 4 s W = -36, p = 0.008; PPN, 0.25 s t = 3.9, p = 0.004, 1 s W = -55, p = 0.002, 4 s t = 3.8, p = 0.004).
plan = [("STN", 0.25, "t", "t = 6.4, p = 0.0004"), ("STN", 1.0, "Wilcoxon", "W = -36, p = 0.005"), ("STN", 4.0, "Wilcoxon", "W = -36, p = 0.008"),
("PPN", 0.25, "t", "t = 3.9, p = 0.004"), ("PPN", 1.0, "Wilcoxon", "W = -55, p = 0.002"), ("PPN", 4.0, "t", "t = 3.8, p = 0.004")]
ax_dls = pre_during_tests("Anxa1", "DLS", plan, label="Recordings labelled DLS (main):")
print("Significant decreases:", int(((ax_dls.p < 0.05) & (ax_dls.mean_during < ax_dls.mean_pre)).sum()), "of 6")
# Variant: the unsuffixed (SNc-labelled) recordings counted as DLS.
R_alt = R.copy()
R_alt.loc[(R_alt.group == "Anxa1") & (R_alt.kind == "single"), "site"] = "DLS"
R_main = R
R = R_alt
ax_alt = pre_during_tests("Anxa1", "DLS", plan, label="\nVariant, unsuffixed recordings added to DLS:")
R = R_main
print("Variant significant decreases:", int(((ax_alt.p < 0.05) & (ax_alt.mean_during < ax_alt.mean_pre)).sum()), "of 6")
print("\nSmallest attainable two-sided Wilcoxon p: n 5 ->", ph.wilcoxon_exact_p(15, 5), "; n 4 ->", ph.wilcoxon_exact_p(10, 4))
for w, paper_p in [(-36, "0.005 and 0.008"), (-55, "0.002")]:
n_all = [n for n in range(3, 15) if n * (n + 1) // 2 == abs(w)]
print(f"Paper W = {w} equals -n(n+1)/2 (every mouse down) for n = {n_all}; exact p {ph.wilcoxon_exact_p(w, n_all[0]):.4f} (paper {paper_p})")
for t_val, p_val in [(6.4, 0.0004), (3.9, 0.004), (3.8, 0.004)]:
print(f"Paper t = {t_val}, p = {p_val}: two-sided p with df", {df: round(2 * stats.t.sf(t_val, df), 4) for df in (3, 4, 7, 9)})
Recordings labelled DLS (main):
input duration test n stat p mean_during mean_pre mice_up mean_post paper other_test_p
STN 0.25 paired t 5 -6.14 0.0036 -1.6481 -0.1752 0 0.9765 t = 6.4, p = 0.0004 0.0625
STN 1.00 Wilcoxon 5 -15.00 0.0625 -1.7306 -0.0000 0 1.0182 W = -36, p = 0.005 0.0189
STN 4.00 Wilcoxon 5 -15.00 0.0625 -1.4565 -0.0399 0 1.4806 W = -36, p = 0.008 0.0181
PPN 0.25 paired t 4 -3.46 0.0405 -1.7299 -0.0643 0 0.5551 t = 3.9, p = 0.004 0.1250
PPN 1.00 Wilcoxon 4 -10.00 0.1250 -1.4833 -0.0000 0 0.5686 W = -55, p = 0.002 0.0520
PPN 4.00 paired t 4 -2.09 0.1276 -0.7699 0.0327 1 1.0790 t = 3.8, p = 0.004 0.2500
Significant decreases: 2 of 6
Variant, unsuffixed recordings added to DLS:
input duration test n stat p mean_during mean_pre mice_up mean_post paper other_test_p
STN 0.25 paired t 8 -2.65 0.0328 -1.0140 -0.1273 2 0.6213 t = 6.4, p = 0.0004 0.0391
STN 1.00 Wilcoxon 8 -34.00 0.0156 -1.1472 -0.0000 1 0.6539 W = -36, p = 0.005 0.0236
STN 4.00 Wilcoxon 8 -36.00 0.0078 -1.1160 -0.0331 0 0.7424 W = -36, p = 0.008 0.0061
PPN 0.25 paired t 10 -2.47 0.0354 -0.8254 -0.0584 1 0.2544 t = 3.9, p = 0.004 0.0137
PPN 1.00 Wilcoxon 10 -51.00 0.0059 -0.7327 -0.0000 1 0.1513 W = -55, p = 0.002 0.0258
PPN 4.00 paired t 10 -2.98 0.0153 -0.5625 -0.0532 1 0.1957 t = 3.8, p = 0.004 0.0098
Variant significant decreases: 6 of 6
Smallest attainable two-sided Wilcoxon p: n 5 -> 0.0625 ; n 4 -> 0.125
Paper W = -36 equals -n(n+1)/2 (every mouse down) for n = [8]; exact p 0.0078 (paper 0.005 and 0.008)
Paper W = -55 equals -n(n+1)/2 (every mouse down) for n = [10]; exact p 0.0020 (paper 0.002)
Paper t = 6.4, p = 0.0004: two-sided p with df {3: np.float64(0.0077), 4: np.float64(0.0031), 7: np.float64(0.0004), 9: np.float64(0.0001)}
Paper t = 3.9, p = 0.004: two-sided p with df {3: np.float64(0.0299), 4: np.float64(0.0175), 7: np.float64(0.0059), 9: np.float64(0.0036)}
Paper t = 3.8, p = 0.004: two-sided p with df {3: np.float64(0.032), 4: np.float64(0.0191), 7: np.float64(0.0067), 9: np.float64(0.0042)}
Result: blocked (data missing). Anxa1 DLS terminals fall during the train: no STN mouse (of 5) goes up in any condition, and no PPN mouse (of 4) except one in 4 s trains. But the dandiset holds DLS recordings for 5 (STN) and 4 (PPN) mice, where the paper has 8 and 10. At these numbers the Wilcoxon tests the paper uses for three of the six conditions cannot reach p < 0.05 (smallest attainable p 0.0625 for 5 mice and 0.125 for 4), so "reproduced" (6 of 6) cannot be reached with the deposited mice. The three t-tests can reach it, so "partially reproduced" (3 to 5 of 6) could be reached; 2 of 6 tests are significant, the two t-tests other than PPN 4 s (p = 0.13). The paper's statistics fix its n: W = -36 and -55 are -n(n+1)/2 for n = 8 and 10 (every mouse down; exact p 0.0078 and 0.0020), and t = 6.4 with p = 0.0004 needs 7 df (8 mice), where 4 df give p = 0.0031. With the unsuffixed recordings added to DLS (a variant), all 6 tests are significant. Its STN 4 s Wilcoxon gives the paper's W = -36, but the other two W (-34 and -51 against -36 and -55) and the three t values (2.65, 2.47 and 2.98 against 6.4, 3.9 and 3.8) are not the paper's, so those recordings are not the paper's DLS data. The DLS recordings of nine of the paper's mice are not in the dandiset, and the outcome is blocked for that reason.
fig5-anxa1-dls-frequency¶
The Anxa1+ terminal inhibition scales with frequency (STN, 0.25 s F = 18.0, 1 s F = 15.5, 4 s F = 15.3; PPN, 0.25 s F = 11.5, 1 s F2,15 = 11.0, 4 s F2,18 = 7.2).
plan = [("STN", 0.25, "Friedman", "F = 18.0, p = 0.0003"), ("STN", 1.0, "Friedman", "F = 15.5, p = 0.015"), ("STN", 4.0, "Friedman", "F = 15.3, p = 0.001"),
("PPN", 0.25, "Friedman", "F = 11.5, p = 0.0002"), ("PPN", 1.0, "RM", "F2,15 = 11.0, p = 0.006"), ("PPN", 4.0, "RM", "F2,18 = 7.2, p = 0.002")]
ax_freq = frequency_tests("Anxa1", "DLS", plan, -1)
R = R_alt
print("\nVariant, unsuffixed recordings added to DLS:")
ax_freq_alt = frequency_tests("Anxa1", "DLS", plan, -1)
R = R_main
print("\nChi-square p for the paper's Friedman values with 3 df: F = 11.5 ->", round(stats.chi2.sf(11.5, 3), 4),
"(paper 0.0002); F = 15.5 ->", round(stats.chi2.sf(15.5, 3), 4), "(paper 0.015)")
for lab, d1, d2 in [("PPN 1 s, F2,15", 2, 15), ("PPN 4 s, F2,18", 2, 18)]:
lo, hi = (d2 - 0.5) / (d1 + 0.5) + 1, (d2 + 0.5) / (d1 - 0.5) + 1
print(f"{lab}: Greenhouse-Geisser df2/df1 = n - 1, so with the dfs rounded n lies between {lo:.1f} and {hi:.1f} ({d2 / d1 + 1:.1f} at the stated values); here n =",
int(np.isfinite(freq_matrix("Anxa1", "PPN", 1.0 if "1 s" in lab else 4.0, "DLS").values).all(axis=1).sum()))
input duration test mice n stat p mean_5 mean_10 mean_20 mean_40 holds paper other_test_p max_friedman dropped STN 0.25 Friedman 5 5 10.68 0.0136 -0.0479 -1.3536 -1.9146 -1.7067 True F = 18.0, p = 0.0003 0.0006 15 - STN 1.00 Friedman 5 5 10.20 0.0169 0.5834 -0.9902 -1.4720 -1.0517 True F = 15.5, p = 0.015 0.0490 15 - STN 4.00 Friedman 5 5 11.16 0.0109 0.3309 -0.4909 -0.9318 -0.8393 True F = 15.3, p = 0.001 0.0087 15 - PPN 0.25 Friedman 4 4 9.30 0.0256 -0.4329 -1.3009 -1.5115 -1.7255 True F = 11.5, p = 0.0002 0.0328 12 - PPN 1.00 RM 4 4 15.82 0.0043 -0.1726 -0.6284 -1.4385 -1.0717 True F2,15 = 11.0, p = 0.006 0.0112 12 - PPN 4.00 RM 4 4 1.95 0.2447 -0.3125 -0.3430 -0.6810 -0.5362 False F2,18 = 7.2, p = 0.002 0.2725 12 - Tests that hold: 5 of 6 Variant, unsuffixed recordings added to DLS: input duration test mice n stat p mean_5 mean_10 mean_20 mean_40 holds paper other_test_p max_friedman dropped STN 0.25 Friedman 8 8 13.05 0.0045 -0.0071 -0.8222 -1.2517 -1.1307 True F = 18.0, p = 0.0003 0.0051 24 - STN 1.00 Friedman 8 8 5.25 0.1544 0.3251 -0.6036 -0.9916 -0.6594 False F = 15.5, p = 0.015 0.0666 24 - STN 4.00 Friedman 8 8 13.95 0.0030 0.1744 -0.3331 -0.7693 -0.6600 True F = 15.3, p = 0.001 0.0048 24 - PPN 0.25 Friedman 10 10 4.08 0.2530 -0.2523 -0.6348 -0.7226 -0.7940 False F = 11.5, p = 0.0002 0.0470 30 - PPN 1.00 RM 10 10 4.50 0.0475 -0.1448 -0.3388 -0.6842 -0.4745 True F2,15 = 11.0, p = 0.006 0.1689 30 - PPN 4.00 RM 10 10 3.15 0.0602 -0.1882 -0.1641 -0.3996 -0.3372 False F2,18 = 7.2, p = 0.002 0.1447 30 - Tests that hold: 3 of 6 Chi-square p for the paper's Friedman values with 3 df: F = 11.5 -> 0.0093 (paper 0.0002); F = 15.5 -> 0.0014 (paper 0.015) PPN 1 s, F2,15: Greenhouse-Geisser df2/df1 = n - 1, so with the dfs rounded n lies between 6.8 and 11.3 (8.5 at the stated values); here n = 4 PPN 4 s, F2,18: Greenhouse-Geisser df2/df1 = n - 1, so with the dfs rounded n lies between 8.0 and 13.3 (10.0 at the stated values); here n = 4
Result: partially reproduced. Among Anxa1 DLS recordings, 5 of 6 tests find an effect of frequency with the 40 Hz mean below the 5 Hz mean (p from 0.0043 to 0.026, with 5 STN and 4 PPN mice); PPN 4 s does not (F = 1.95, p = 0.24). With the unsuffixed recordings added, 3 of 6 hold. The samples are small for the same reason as in fig5-anxa1-dls-inhibited: the paper's F2,15 and F2,18 for PPN 1 and 4 s imply about 8.5 and 10 mice (between 6.8 and 13.3 with the dfs rounded), against 4 here. Two of the paper's Friedman entries do not fit their p-values (F = 11.5 gives p = 0.0093 against the stated 0.0002; F = 15.5 gives 0.0014 against 0.015), and the Figure 5 legend says the STN panels (N, O) show no significant effect of frequency while the Results report one for all three STN durations.
fig5-anxa1-dls-rebound¶
Anxa1+ terminals in DLS show a transient increase in activity at stimulus offset, for both inputs.
for label, RR in [("main, labelled DLS", R_main), ("variant, unsuffixed added", R_alt)]:
R = RR
for inp in ["STN", "PPN"]:
for dur in [1.0, 0.25, 4.0]:
c = cohort("Anxa1", inp, "DLS", duration=dur)
res = stats.ttest_1samp(c.post, 0)
print(f"[{label}] {inp} {dur} s: n {len(c)}, mean in 1 s after offset {c.post.mean():.2f} % (range {c.post.min():.2f} to {c.post.max():.2f}), "
f"mice above 0 {int((c.post > 0).sum())}, t {res.statistic:.2f}, p {res.pvalue:.2g}")
R = R_main
fig, axes = plt.subplots(1, 2, figsize=(10, 3), sharey=True)
for ax, inp in zip(axes, ["STN", "PPN"]):
for s, kind, lab in [("DLS", None, "DLS"), ("SNc", "dual", "SNc, dual mice"), ("SNc", "single", "unsuffixed")]:
lo, hi = plot_traces(ax, "Anxa1", inp, s, duration=1.0, kind=kind, label=lab)
print(f"{inp} {lab}: plotted mean ranges {lo:.2f} to {hi:.2f} %")
ax.set_title(f"Anxa1, {inp}, 1 s 40 Hz"); ax.legend(); ax.set_xlabel("s from onset")
axes[0].set_ylabel("dF/F (%)"); plt.show()
[main, labelled DLS] STN 1.0 s: n 5, mean in 1 s after offset 1.02 % (range 0.64 to 1.65), mice above 0 5, t 5.79, p 0.0044 [main, labelled DLS] STN 0.25 s: n 5, mean in 1 s after offset 0.98 % (range 0.36 to 2.11), mice above 0 5, t 3.11, p 0.036 [main, labelled DLS] STN 4.0 s: n 5, mean in 1 s after offset 1.48 % (range 0.54 to 2.63), mice above 0 5, t 3.46, p 0.026
[main, labelled DLS] PPN 1.0 s: n 4, mean in 1 s after offset 0.57 % (range -0.20 to 1.03), mice above 0 3, t 1.99, p 0.14 [main, labelled DLS] PPN 0.25 s: n 4, mean in 1 s after offset 0.56 % (range 0.04 to 1.55), mice above 0 4, t 1.64, p 0.2 [main, labelled DLS] PPN 4.0 s: n 4, mean in 1 s after offset 1.08 % (range -0.01 to 2.35), mice above 0 3, t 2.22, p 0.11 [variant, unsuffixed added] STN 1.0 s: n 8, mean in 1 s after offset 0.65 % (range -0.48 to 1.65), mice above 0 7, t 2.90, p 0.023 [variant, unsuffixed added] STN 0.25 s: n 8, mean in 1 s after offset 0.62 % (range -0.04 to 2.11), mice above 0 6, t 2.43, p 0.046 [variant, unsuffixed added] STN 4.0 s: n 8, mean in 1 s after offset 0.74 % (range -1.45 to 2.63), mice above 0 6, t 1.58, p 0.16 [variant, unsuffixed added] PPN 1.0 s: n 10, mean in 1 s after offset 0.15 % (range -0.62 to 1.03), mice above 0 4, t 0.89, p 0.4 [variant, unsuffixed added] PPN 0.25 s: n 10, mean in 1 s after offset 0.25 % (range -0.19 to 1.55), mice above 0 7, t 1.65, p 0.13 [variant, unsuffixed added] PPN 4.0 s: n 10, mean in 1 s after offset 0.20 % (range -0.72 to 2.35), mice above 0 3, t 0.64, p 0.54
STN DLS: plotted mean ranges -2.51 to 2.28 % STN SNc, dual mice: plotted mean ranges -0.87 to 0.60 % STN unsuffixed: plotted mean ranges -0.39 to 0.60 %
PPN DLS: plotted mean ranges -2.61 to 2.22 % PPN SNc, dual mice: plotted mean ranges -0.58 to 0.54 %
PPN unsuffixed: plotted mean ranges -0.42 to 0.26 %
Result: partially reproduced. After 1 s trains Anxa1 DLS terminals rise above baseline for STN (1.02 %, 5 of 5 mice, p = 0.0044) but not significantly for PPN (0.57 %, 3 of 4, p = 0.14, n = 4). The same holds for 0.25 s and 4 s trains (STN p = 0.036 and 0.026; PPN p = 0.20 and 0.11).
figs4-anxa1-soma-inhibited¶
Anxa1+ SNc cell bodies show stimulation-induced inhibition (paired t-test, STN 0.25 s t = 3.1, 1 s t = 4.0, 4 s t = 7.2, n = 5; PPN t = 4.8, 8.4, 3.3, n = 4) followed by post-stimulus rebound activation.
plan = [("STN", 0.25, "t", "t = 3.1, p = 0.035"), ("STN", 1.0, "t", "t = 4.0, p = 0.016"), ("STN", 4.0, "t", "t = 7.2, p = 0.002"),
("PPN", 0.25, "t", "t = 4.8, p = 0.018"), ("PPN", 1.0, "t", "t = 8.4, p = 0.004"), ("PPN", 4.0, "t", "t = 3.3, p = 0.044")]
soma = pre_during_tests("Anxa1", "SNc", plan, label="Recordings labelled SNc (main):")
print("Significant decreases:", int(((soma.p < 0.05) & (soma.mean_during < soma.mean_pre)).sum()), "of 6")
soma_d = pre_during_tests("Anxa1", "SNc", plan, kind="dual", label="\nVariant, mice with both a DLS and an SNc file (the paper's n):")
print("Variant significant decreases:", int(((soma_d.p < 0.05) & (soma_d.mean_during < soma_d.mean_pre)).sum()), "of 6")
for kind in [None, "dual"]:
for inp in ["STN", "PPN"]:
c = cohort("Anxa1", inp, "SNc", duration=1.0, kind=kind)
res = stats.ttest_1samp(c.post, 0)
print(f"Rebound, {kind or 'all SNc'} {inp} 1 s: mean in 1 s after offset {c.post.mean():.2f} %, mice above 0 {int((c.post > 0).sum())} of {len(c)}, t {res.statistic:.2f}, p {res.pvalue:.2g}")
Recordings labelled SNc (main): input duration test n stat p mean_during mean_pre mice_up mean_post paper other_test_p STN 0.25 paired t 8 -1.56 0.1617 -0.2756 -0.0131 2 0.1376 t = 3.1, p = 0.035 0.1953 STN 1.00 paired t 8 -3.33 0.0126 -0.3260 0.0000 1 -0.0771 t = 4.0, p = 0.016 0.0234 STN 4.00 paired t 8 -6.45 0.0003 -0.7679 0.0064 0 -0.1692 t = 7.2, p = 0.002 0.0078 PPN 0.25 paired t 10 -3.03 0.0143 -0.3023 0.0153 1 0.1142 t = 4.8, p = 0.018 0.0137 PPN 1.00 paired t 10 -4.11 0.0027 -0.2757 -0.0000 1 0.0288 t = 8.4, p = 0.004 0.0059 PPN 4.00 paired t 10 -4.61 0.0013 -0.4542 -0.0634 0 -0.2156 t = 3.3, p = 0.044 0.0020 Significant decreases: 5 of 6 Variant, mice with both a DLS and an SNc file (the paper's n): input duration test n stat p mean_during mean_pre mice_up mean_post paper other_test_p STN 0.25 paired t 5 -3.13 0.0353 -0.4667 0.0074 0 0.2025 t = 3.1, p = 0.035 0.0625 STN 1.00 paired t 5 -3.98 0.0164 -0.4167 0.0000 0 -0.1514 t = 4.0, p = 0.016 0.0625 STN 4.00 paired t 5 -7.26 0.0019 -0.8996 0.0233 0 0.0220 t = 7.2, p = 0.002 0.0625 PPN 0.25 paired t 4 -4.72 0.0180 -0.4221 0.1199 0 0.2046 t = 4.8, p = 0.018 0.1250 PPN 1.00 paired t 4 -8.35 0.0036 -0.3407 -0.0000 0 0.2623 t = 8.4, p = 0.004 0.1250 PPN 4.00 paired t 4 -3.33 0.0447 -0.4991 0.0073 0 0.0507 t = 3.3, p = 0.044 0.1250 Variant significant decreases: 6 of 6 Rebound, all SNc STN 1 s: mean in 1 s after offset -0.08 %, mice above 0 3 of 8, t -0.47, p 0.65 Rebound, all SNc PPN 1 s: mean in 1 s after offset 0.03 %, mice above 0 5 of 10, t 0.26, p 0.8 Rebound, dual STN 1 s: mean in 1 s after offset -0.15 %, mice above 0 1 of 5, t -0.68, p 0.53 Rebound, dual PPN 1 s: mean in 1 s after offset 0.26 %, mice above 0 4 of 4, t 1.88, p 0.16
Result: partially reproduced. Among the recordings labelled SNc (8 STN and 10 PPN mice), 5 of 6 tests show a significant decrease during the train (p from 0.0003 to 0.014); STN 0.25 s does not (t = -1.56, p = 0.16). Restricted to the mice with both a DLS and an SNc file (5 and 4, the paper's n; a variant), all 6 are significant with the paper's statistics (t = 3.13, 3.98, 7.26, 4.72, 8.35, 3.33 against 3.1, 4.0, 7.2, 4.8, 8.4, 3.3). No post-stimulus rebound is seen in the cell bodies after 1 s trains (all SNc, STN -0.08 % and PPN 0.03 %, p = 0.65 and 0.80; dual mice -0.15 % and 0.26 %, p = 0.53 and 0.16).
fig5-opposite-subtypes¶
Excitatory inputs from STN and PPN activate PD-resilient Vglut2+ DA neurons but inhibit PD-vulnerable Anxa1+ DA neurons.
for label, RR in [("main", R_main), ("variant, unsuffixed Anxa1 added to DLS", R_alt)]:
R = RR
for inp in ["STN", "PPN"]:
v, a = cohort("Vglut2", inp, "SNc", duration=1.0), cohort("Anxa1", inp, "DLS", duration=1.0)
res = stats.ttest_ind(v.during, a.during, equal_var=False)
print(f"[{label}] {inp}: Vglut2 n {len(v)} mean during {v.during.mean():.2f} % ; Anxa1 DLS n {len(a)} mean {a.during.mean():.2f} % ; "
f"Welch t {res.statistic:.2f}, p {res.pvalue:.2g}")
R = R_main
for inp in ["STN", "PPN"]:
v, a = cohort("Vglut2", inp, "SNc", duration=1.0), cohort("Anxa1", inp, "SNc", duration=1.0)
res = stats.ttest_ind(v.during, a.during, equal_var=False)
print(f"[variant, Anxa1 SNc] {inp}: Anxa1 SNc n {len(a)} mean {a.during.mean():.2f} % ; Welch t {res.statistic:.2f}, p {res.pvalue:.2g}")
[main] STN: Vglut2 n 9 mean during 6.54 % ; Anxa1 DLS n 5 mean -1.73 % ; Welch t 3.56, p 0.0066 [main] PPN: Vglut2 n 10 mean during 9.35 % ; Anxa1 DLS n 4 mean -1.48 % ; Welch t 3.90, p 0.0032 [variant, unsuffixed Anxa1 added to DLS] STN: Vglut2 n 9 mean during 6.54 % ; Anxa1 DLS n 8 mean -1.15 % ; Welch t 3.32, p 0.0097 [variant, unsuffixed Anxa1 added to DLS] PPN: Vglut2 n 10 mean during 9.35 % ; Anxa1 DLS n 10 mean -0.73 % ; Welch t 3.67, p 0.005 [variant, Anxa1 SNc] STN: Anxa1 SNc n 8 mean -0.33 % ; Welch t 3.01, p 0.017 [variant, Anxa1 SNc] PPN: Anxa1 SNc n 10 mean -0.28 % ; Welch t 3.52, p 0.0065
Result: reproduced. During 1 s 40 Hz trains Vglut2+ SNc neurons rise and Anxa1+ DLS terminals fall with both inputs: STN 6.54 % against -1.73 % (Welch t = 3.56, p = 0.0066), PPN 9.35 % against -1.48 % (t = 3.90, p = 0.0032). The difference holds with the unsuffixed Anxa1 recordings added to DLS (p = 0.0097 and 0.005) and against the Anxa1 SNc recordings (p = 0.017 and 0.0065).
print(f"Notebook ran in {(time.time() - started) / 60:.1f} min")
Notebook ran in 7.4 min
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 | Data | Outcome |
|---|---|---|
| cohort-sizes | dandiset | not reproduced |
| fig2-gaba-excited | dandiset | reproduced |
| fig2-gaba-frequency | dandiset | reproduced |
| fig3b-da-stn-250ms-peak | dandiset | partially reproduced |
| fig3c-da-stn-duration | dandiset | reproduced |
| fig3d-da-stn-frequency-ns | dandiset | reproduced |
| fig3g-da-ppn-duration | dandiset | not reproduced |
| fig3h-da-ppn-frequency-ns | dandiset | reproduced |
| fig4cf-ts-dls-opposite | dandiset | reproduced |
| fig4dg-ts-frequency | dandiset | reproduced |
| fig4eh-dls-frequency | dandiset | reproduced |
| fig4eh-dls-5hz-increase | dandiset | not reproduced |
| fig4cf-dls-rebound | dandiset | partially reproduced |
| figs3-terminals-opposite | dandiset | reproduced |
| fig5-vglut2-excited | dandiset | reproduced |
| fig5-vglut2-frequency | dandiset | reproduced |
| fig5-anxa1-dls-inhibited | dandiset | blocked: data missing |
| fig5-anxa1-dls-frequency | dandiset | partially reproduced |
| fig5-anxa1-dls-rebound | dandiset | partially reproduced |
| figs4-anxa1-soma-inhibited | dandiset | partially reproduced |
| fig5-opposite-subtypes | dandiset | reproduced |
| fig1de, fig1f, figs1 counts, figs1 IF (human tissue) | dandiset | blocked: data missing |
How the analysis differs from the paper's. It uses the lab's dF/F stored in the files, the Methods' 1 s baseline, and windows the Methods do not give: the mean during the train, a pre-stimulus window of the train's length, and the 1 s after offset. Where a panel shows one train duration without saying which, the 1 s train is used. The paper chose between parametric and non-parametric tests by a Shapiro-Wilk test; here each comparison uses the test the paper reports for it. Repeated-measures ANOVAs are Greenhouse-Geisser corrected, as the paper's fractional degrees of freedom imply. For many comparisons this reproduces the paper's statistic to the first decimal, which suggests the windows match the lab's. The Anxa1 groups follow the dandiset's site labels.
For 1 s trains the pre-stimulus window is the baseline window, so tests of a 1 s train against its pre-stimulus window are tests against 0.
What was found about the dandiset.
- It holds the photometry of Figures 2 to 5, S3 and S4 for 120 mice; the human tissue counts of Figure 1 and S1 and the histology of Figure S2 are not in it.
- For the Anxa1 mice it holds DLS terminal recordings for 5 (STN) and 4 (PPN) mice where the paper has 8 and 10, and SNc recordings for 8 and 10 where the paper has 5 and 4. The nine extra SNc recordings carry no site suffix in their file names; their signals resemble the SNc recordings, and the paper's cell-body statistics are matched by the other SNc recordings alone. The DLS recordings of nine mice appear to be missing.
- DA terminals in DLS have 8 mice (paper 9). Vglut2+ neurons with STN stimulation have 9 mice where the legend gives 10, but the paper's Wilcoxon statistics fit 9, so the legend is likely wrong.
- Three Vglut2 PPN mice (C3412, C3798, C4707) have no 5 Hz trains, while the paper's Friedman statistic needs at least 9 mice with all four frequencies.
- The dF/F is stated in "a.u." but is in percent.
- Thirteen mice have 10 trials per condition (Methods: 5), 30 mice have 16 s trains the paper does not mention, and 14 frequencies sessions have no 4 s trains, one of them (C4608) because a separate recording was not converted.
- Every session start time carries the offset -00:01, an artefact of attaching the London time
zone with
replace, so the true offset is lost; eight subject ids are lower case. - The files' description of the frequencies protocol gives an interval of 10 s between trains; the epochs and the Methods give 20 s.
- The conversion code names the 1 s block of the STN to SNr GABA frequencies sessions as 10 mW, while the files and the Methods give 5 mW.
- The epoch tables share the clock of the dF/F (the responses lock to the true onsets in 235 of 258 recordings), and no recording is deposited twice.
- The dandiset has no published version, so results hold for the draft as modified on 2026-09-01.