Reproducing Bouabid, Vu et al. (2026) from dandiset 001692¶
Paper. Bouabid S, Vu M-AT, Noggle C, Vietti-Michelina S, Brimblecombe K, Platt N, Zhang L, Joshi A, Cragg S, Howe MW (2026). An anatomical hotspot for striatal dopamine-acetylcholine interactions during reward and movement. bioRxiv. doi:10.64898/2026.01.20.700614. This attempt used the full text of the preprint (version 1, CC-BY).
Dandiset. DANDI:001692, version 0.260303.1724. The paper's Key Resources Table (Zenodo doi:10.5281/zenodo.18202873) and its code's README cite this version, which is also the only published one, so it is pinned here.
What the dandiset holds. Multi-fiber photometry of DA (GRAB-rDA3m, red) and ACh (GRAB-ACh3.0, green) in head-fixed mice: spontaneous and unpredicted-reward sessions of 5 mice (UG27 to UG31, Figures 1 to 3, 5 and 6), Pavlovian sessions of 3 mice (AD1 to AD3, Figure 4), sessions of 3 control mice with non-functional sensors (AD4 to AD6, Extended Data Figure 1), and stimulation sessions of 6 DAT-cre mice (ADS*, Figure 7). The ex vivo slice data of Figure 8 are not in it; their source data are on Zenodo and are used for the Figure 8 claims.
Code read. The paper's analysis code, github.com/HoweLab/BouabidVu2026 at commit 833fa3d (MATLAB). Where the Methods leave a choice open the code's choice is followed; where they disagree the Methods are followed and the code's choice is a variant.
| Claim | Tested here |
|---|---|
| fig1e-dominant-polarity: dominant correlation signs and lead (132/127, 42/31 of 174) | yes |
| fig1f-adls-hotspot: negative hotspot in aDLS | yes |
| fig1h-in-vs-out: more negative inside the hotspot | yes |
| fig1f-positive-posterior: positive correlations posterior | yes |
| ed2-hotspot-per-mouse: in-out difference in each mouse | yes |
| ed1-controls-no-organization: no hotspot in control mice | yes |
| fig2a-directionality: DA leads opposite-sign pairs, ACh same-sign pairs | yes |
| fig2-magnitude-correlations: 121, 81, 97, 78 of 174 sites | yes |
| fig2d-magnitude-in-out: opposite-sign magnitude correlations stronger inside | yes |
| fig2b-pairing-overlaps-hotspot: pairing hotspots overlap the correlation hotspot | yes |
| ed3e-transient-rates-in-hotspot: higher transient rates inside | yes |
| fig3b-reward-correlation: 114 negative, 58 positive, latencies | yes |
| fig3d-reward-in-out: more negative inside | yes |
| fig3e-reward-hotspot-overlap: reward hotspot overlaps | yes |
| fig4-cue-sequence: ACh peak, DA peak, ACh dip at 49 of 75 sites | yes |
| fig4c-cue-correlation: 39 negative, 26 positive, latencies | yes |
| fig4h-extinction-correlation: 41 of 74, lag 117.9 ms | yes |
| fig4d-cue-in-out, fig4j-extinction-in-out | yes |
| fig5d-initiation-responses: 131, 119, 73, 32 of 174 | yes |
| fig5-latency-order: DA dip, ACh and DA peaks, ACh dip | yes |
| fig5e-initiation-topography | yes |
| fig6a-initiation-correlation, fig6c-initiation-in-out, fig6g-invigoration-in-out | yes |
| fig7-da-release (179/180), fig7e-ach-suppression (116/180) | yes |
| ed9-cohort-hotspot, fig7f-ach-dip-in-out | yes |
| fig8b, fig8e, fig8fg: slice ANOVAs (Zenodo source data) | yes, data elsewhere |
| fig8c-cocaine-da-similar | blocked (methods underspecified) |
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy import stats
import daach
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
started = time.time()
assets = daach.list_assets()
subject_of = lambda a: a["path"].split("/")[0][4:]
UG = [a for a in assets if subject_of(a).startswith("UG")]
ADS = [a for a in assets if subject_of(a).startswith("ADS")]
AD = [a for a in assets if subject_of(a).startswith("AD") and not subject_of(a).startswith("ADS")]
ug = daach.load_many(UG, series="hp_art")
ad_events = daach.load_many(AD, func=daach.read_events)
ads = daach.load_many(ADS, series="raw")
def light_reward_fraction(e):
ev = e["events"]
if "Light" not in ev:
return np.nan
rewards = ev["Reward"][0] if "Reward" in ev else np.array([])
return np.mean([((rewards > a) & (rewards < b)).any() for a, b in zip(*ev["Light"])])
pav_table = pd.DataFrame([{"path": e["path"], "mouse": e["subject"], "session": e["session_id"],
"indicators": "+".join(e["indicators"]), "events": ",".join(sorted(e["events"])),
"n_light": len(e["events"].get("Light", ([],))[0]),
"n_reward": len(e["events"].get("Reward", ([],))[0]),
"light_rewarded": light_reward_fraction(e)} for e in ad_events])
by_path = {a["path"]: a for a in assets}
selected = []
for m in ("AD1", "AD2", "AD3"):
t = pav_table[pav_table.mouse == m].sort_values("session")
selected += [(p, "learning") for p in t[t.light_rewarded == 1.0].path.iloc[-3:]]
selected += [(p, "extinction") for p in t[np.isclose(t.light_rewarded, 0.2)].path.iloc[-3:]]
pav = daach.load_many([by_path[p] for p, _ in selected], series="hp")
for s, (_, phase) in zip(pav, selected):
s["phase"] = phase
control_paths = pav_table[pav_table.mouse.isin(["AD4", "AD5", "AD6"]) & (pav_table.n_light == 0)
& (pav_table.n_reward > 0)].path
ctrl = daach.load_many([by_path[p] for p in control_paths], series="hp")
print(f"{len(assets)} assets in dandiset {daach.DANDISET_ID} version {daach.DANDISET_VERSION}; streamed "
f"{len(ug)} UG, {len(ad_events)} AD event tables, {len(pav)} Pavlovian, {len(ctrl)} control and "
f"{len(ads)} ADS sessions in {time.time() - started:.0f} s")
339 assets in dandiset 001692 version 0.260303.1724; streamed 105 UG, 211 AD event tables, 18 Pavlovian, 15 control and 23 ADS sessions in 224 s
Data check¶
Each check of the metadata against the values is stated under its output.
files = pd.DataFrame([{"mouse": subject_of(a), "size_GB": a["size"] / 1e9} for a in assets])
summary = files.groupby("mouse").agg(files=("size_GB", "size"), total_GB=("size_GB", "sum")).round(1)
ind = {s["subject"]: "+".join(s["ACh_indicator"] + s["DA_indicator"]) for s in ug + ads if "missing" not in s}
ind.update({e["subject"]: "+".join(e["indicators"]) for e in ad_events})
summary["indicators"] = summary.index.map(ind)
print(summary.to_string())
print("Subjects:", len(summary), "; the Methods list 17 in vivo mice (5 spontaneous, 3 Pavlovian, 3 control, 6 DAT-cre).")
ug_ok = [s for s in ug if "missing" not in s]
print("\nSessions skipped because their high-pass, artifact-removed series is absent (the code skips empty ones):",
[(s["path"], s["missing"]) for s in ug if "missing" in s])
spont = [s for s in ug_ok if "Reward" not in s["events"]]
reward = [s for s in ug_ok if "Reward" in s["events"]]
print(f"UG27-UG31 sessions: {len(ug)} files, {len(spont)} without rewards (spontaneous), {len(reward)} with rewards")
print("Event types in spontaneous sessions:", sorted({k for s in spont for k in s["events"]}))
files total_GB indicators
mouse
AD1 38 519.8 ACh3.0+rDA3m
AD2 26 584.2 ACh3.0+rDA3m
AD3 33 751.6 ACh3.0+rDA3m
AD4 29 488.5 ACh4.3mut+tdTomato
AD5 41 712.7 ACh4.3mut+tdTomato
AD6 44 756.6 ACh4.3mut+tdTomato
ADS12 4 10.5 ACh3.0+rDA3m
ADS13 4 10.3 ACh3.0+rDA3m
ADS16 4 10.8 ACh3.0+rDA3m
ADS17 4 9.7 ACh3.0+rDA3m
ADS19 4 10.7 ACh3.0+rDA3m
ADS6 3 5.5 ACh3.0+rDA3m
UG27 22 150.2 ACh3.0+rDA3m
UG28 20 127.3 ACh3.0+rDA3m
UG29 24 160.7 ACh3.0+rDA3m
UG30 18 108.3 ACh3.0+rDA3m
UG31 21 128.4 ACh3.0+rDA3m
Subjects: 17 ; the Methods list 17 in vivo mice (5 spontaneous, 3 Pavlovian, 3 control, 6 DAT-cre).
Sessions skipped because their high-pass, artifact-removed series is absent (the code skips empty ones): [('sub-UG27/sub-UG27_ses-240313r_ophys.nwb', ['Green'])]
UG27-UG31 sessions: 105 files, 61 without rewards (spontaneous), 43 with rewards
Event types in spontaneous sessions: [np.str_('Lick')]
desc = daach.load_many([by_path["sub-UG27/sub-UG27_ses-240214r_ophys.nwb"], AD[[subject_of(a) for a in AD].index("AD4")]],
func=daach.read_descriptions)
for name, d in zip(("UG27", "AD4"), desc):
print(f"--- {name}")
for k in ("ExcitationSource470", "ExcitationSource570", "ACh3.0", "rDA3m", "ACh4.3mut", "tdTomato"):
if k in d:
print(f"{k}: {d[k]!r}")
print("session_description:", repr(d["session_description"][:420]))
print("pharmacology (control virus):", repr(d["pharmacology"][-260:]))
--- UG27 ExcitationSource470: 'Blue excitation light (470 nm LED, Thorlabs, No. SOLIS-470C) and violet excitation light (for the isosbestic control)\nwere coupled into the optic fiber such that a power of 0.75 mW was delivered to the fiber tip.\nThen, 470 nm and 405 nm excitation were alternated at 100 Hz using a waveform generator,\neach filtered with a corresponding filter.\n' ExcitationSource570: '405 nm, Thorlabs, No. SOLIS-570C' ACh3.0: 'green acetyl choline sensor' rDA3m: 'Red dopamine sensor' session_description: 'This dataset contains fiber photometry recordings from multi-fiber arrays implanted in the striatum in mice running \non a treadmill wheel or treadmill ball. Some sessions may include water reward delivery, where a water spout mounted \non a post delivered water rewards (9 μL) at random time intervals (randomly drawn from a 5-30s uniform distribution) \nthrough a water spout and solenoid valve gated electronically. Som' pharmacology (control virus): ' 30 locations chosen to maximize expression around fiber tips. For control experiments, 500 nl of a 1:67 mix \nof pAAV-CAG-tdTomato (Adgene #59462) and AAV9-hysn-ACh4.3-mut (WZ Biosciences) 60, was injected into the striatum of WT mice \nusing the same strategy.' --- AD4 ExcitationSource470: 'Blue excitation light (470 nm LED, Thorlabs, No. SOLIS-470C) and violet excitation light (for the isosbestic control)\nwere coupled into the optic fiber such that a power of 0.75 mW was delivered to the fiber tip.\nThen, 470 nm and 405 nm excitation were alternated at 100 Hz using a waveform generator,\neach filtered with a corresponding filter.\n' ExcitationSource570: '405 nm, Thorlabs, No. SOLIS-570C' ACh4.3mut: 'mutant green acetyl choline sensor' tdTomato: 'Red static marker' session_description: 'This dataset contains fiber photometry recordings from multi-fiber arrays implanted in the striatum in mice running \non a treadmill wheel or treadmill ball. Some sessions may include water reward delivery, where a water spout mounted \non a post delivered water rewards (9 μL) at random time intervals (randomly drawn from a 5-30s uniform distribution) \nthrough a water spout and solenoid valve gated electronically. Som' pharmacology (control virus): ' 30 locations chosen to maximize expression around fiber tips. For control experiments, 500 nl of a 1:67 mix \nof pAAV-CAG-tdTomato (Adgene #59462) and AAV9-hysn-ACh4.3-mut (WZ Biosciences) 60, was injected into the striatum of WT mice \nusing the same strategy.'
The files of AD4 to AD6 name the mutant ACh sensor and tdTomato; all others name ACh3.0 and rDA3m.
The files call the mutant sensor ACh4.3mut (device and surgery text), where the Methods and the Key
Resources Table name ACh3.0-mut. The description of the 470 nm source describes a different setup
(470 and 405 nm alternated at 100 Hz, with an isosbestic control), and that of the 570 nm source
reads "405 nm"; the Methods alternate 470 and 570 nm at 36 Hz. The session description gives 9 uL
rewards at 5 to 30 s intervals, where the Methods give 5 uL at 10 to 30 s, and its micro sign is
stored mis-encoded ("9 μL").
The cohorts are told apart only by subject name and indicator; no field in the files states the
experiment or the figure a session belongs to. UG29's session 240308 has licks but no rewards and is
treated as spontaneous, as its name has no r.
rows = []
for name, ss in (("UG", ug_ok), ("Pavlovian", pav), ("control", ctrl), ("ADS", ads)):
for m, sm_ in daach.by_mouse(ss).items():
c0 = sm_[0]["coords"]
same = all(s["coords"].shape == c0.shape and np.allclose(s["coords"], c0) for s in sm_)
rows.append({"cohort": name, "mouse": m, "sessions": len(sm_), "fibers_per_channel": sm_[0]["ACh_n_fibers"],
"included_sites": len(c0), "coords_same_in_all_sessions": same,
"AP_range": f"{c0[:, 0].min():.2f} to {c0[:, 0].max():.2f}",
"ML_range": f"{c0[:, 1].min():.2f} to {c0[:, 1].max():.2f}",
"DV_range": f"{c0[:, 2].min():.2f} to {c0[:, 2].max():.2f}"})
sites = pd.DataFrame(rows)
print(sites.to_string(index=False))
print("\nIncluded sites per cohort:", sites.groupby("cohort").included_sites.sum().to_dict())
print("Paper: 174 sites in 5 mice (Figures 1-3, 5, 6); 75 sites in 3 mice (Figure 4); 180 sites in 6 mice (Figure 7).")
cohort mouse sessions fibers_per_channel included_sites coords_same_in_all_sessions AP_range ML_range DV_range
UG UG27 21 55 23 True -0.99 to 0.94 0.94 to 2.76 -4.47 to -2.39
UG UG28 20 54 41 True -1.42 to 1.73 0.62 to 3.03 -4.56 to -2.28
UG UG29 24 58 48 True -1.40 to 1.65 0.76 to 3.47 -4.46 to -2.10
UG UG30 18 46 28 True -1.39 to 1.31 0.82 to 3.09 -4.21 to -2.37
UG UG31 21 54 34 True -1.38 to 1.60 0.63 to 3.08 -4.36 to -2.32
Pavlovian AD1 6 65 17 True -0.66 to 1.17 0.97 to 2.60 -4.84 to -3.19
Pavlovian AD2 6 60 24 True -1.35 to 1.10 0.93 to 3.38 -4.88 to -2.67
Pavlovian AD3 6 60 34 True -1.63 to 1.00 0.50 to 3.39 -4.63 to -2.50
control AD4 5 59 42 True -1.50 to 1.51 0.68 to 3.38 -5.95 to -2.54
control AD5 5 59 41 True -1.68 to 1.41 0.91 to 3.38 -5.01 to -2.49
control AD6 5 54 42 True -2.02 to 1.24 0.83 to 3.56 -4.88 to -2.38
ADS ADS12 4 63 13 True -0.54 to 1.42 0.45 to 2.34 -5.07 to -3.44
ADS ADS13 4 57 23 True -1.28 to 0.74 0.84 to 3.03 -4.55 to -2.34
ADS ADS16 4 58 37 True -1.45 to 1.43 1.12 to 3.52 -4.68 to -3.12
ADS ADS17 4 66 27 True -0.58 to 1.83 0.90 to 2.31 -4.40 to -2.74
ADS ADS19 4 66 27 True -0.86 to 1.50 0.37 to 2.77 -5.12 to -2.98
ADS ADS6 3 67 53 True -1.73 to 1.42 1.08 to 3.74 -4.77 to -2.07
Included sites per cohort: {'ADS': 180, 'Pavlovian': 75, 'UG': 174, 'control': 125}
Paper: 174 sites in 5 mice (Figures 1-3, 5, 6); 75 sites in 3 mice (Figure 4); 180 sites in 6 mice (Figure 7).
The included sites match the paper's counts: 174 for UG27 to UG31, 75 for AD1 to AD3 and 180 for the six ADS mice. The fiber coordinates (mm from bregma, AP, ML and DV) are the same in every session of a mouse, which is consistent with each site being the same fiber across sessions; the analysis treats it so, as the paper does, and the per-mouse breakdowns below are the check on pooling. All sites lie in the right hemisphere (ML 0.37 to 3.74 mm).
rows = []
for s in ug_ok + pav + ctrl + ads:
t = s["t"]
rows.append({"mouse": s["subject"], "n": len(s["DA"]), "stated_rate_green": s["ACh_stated_rate"],
"stated_rate_red": s["DA_stated_rate"], "behavior_dt_ms": np.median(np.diff(t)) * 1000,
"behavior_rate": (len(t) - 1) / (t[-1] - t[0]), "n_behavior": len(t),
"green_start_minus_behavior_ms": (s["ACh_start"] - t[0]) * 1000,
"red_start_minus_green_ms": (s["DA_start"] - s["ACh_start"]) * 1000,
"minutes_at_stated_rate": len(s["DA"]) / s["ACh_stated_rate"] / 60,
"minutes_from_timestamps": (t[-1] - t[0]) / 60})
clock = pd.DataFrame(rows)
print(clock.describe().loc[["min", "50%", "max"]].round(3).to_string())
print("\nSessions whose photometry and behavior sample counts differ:", int((clock.n != clock.n_behavior).sum()))
n stated_rate_green stated_rate_red behavior_dt_ms behavior_rate n_behavior green_start_minus_behavior_ms red_start_minus_green_ms minutes_at_stated_rate minutes_from_timestamps min 3429.0 8.999 8.999 55.5 18.000 3429.0 0.0 -28.0 6.350 3.174 50% 10989.0 9.000 9.000 55.5 18.000 10989.0 0.0 -27.5 20.349 10.174 max 34603.0 9.000 9.000 55.5 18.004 34603.0 0.0 28.0 64.078 32.039 Sessions whose photometry and behavior sample counts differ: 0
Every photometry series states a rate of 9.0 Hz, but each has as many samples as the behavior
series, starts at the behavior's first timestamp, and the behavior timestamps step by 55.6 ms
(18.0 Hz). The Methods give 18 Hz per channel and the code uses sr = 18. At the stated rate the
recordings would last twice as long as their timestamps. The data support 18 Hz; it is used
throughout, and the stated rate is a metadata error. The red series starts 27.5 ms (half a sample)
before the green, as the alternating excitation implies; DA and ACh are paired by sample index, as
in the code.
rows = []
for name, ss in (("UG hp_art", ug_ok), ("Pavlovian hp", pav), ("control hp", ctrl), ("ADS raw", ads)):
for ch in ("DA", "ACh"):
x = np.concatenate([s[ch].ravel() for s in ss])
flat = sum(int((np.nanstd(s[ch], 0) == 0).sum()) for s in ss)
rows.append({"series": name, "channel": ch, "p1": np.nanpercentile(x, 1), "median": np.nanmedian(x),
"p99": np.nanpercentile(x, 99), "NaN_fraction": np.isnan(x).mean(),
"constant_site_sessions": flat})
print("dF/F (dimensionless; the files state the unit 'n.a.'):")
print(pd.DataFrame(rows).round(4).to_string(index=False))
vel = pd.DataFrame([{"mouse": s["subject"], "min": s["velocity"].min(), "median": np.median(s["velocity"]),
"p99": np.percentile(s["velocity"], 99), "max": s["velocity"].max(),
"frac_above_0.0208": (daach.smooth_velocity(s["velocity"]) > 0.0208).mean()} for s in spont])
print("\nWheel velocity (m/s) in spontaneous sessions, per mouse (median over sessions):")
print(vel.groupby("mouse").median().round(4).to_string())
dF/F (dimensionless; the files state the unit 'n.a.'):
series channel p1 median p99 NaN_fraction constant_site_sessions
UG hp_art DA -0.0326 -0.0000 0.0376 0.0100 0
UG hp_art ACh -0.0208 -0.0004 0.0279 0.0001 0
Pavlovian hp DA -0.0280 -0.0005 0.0377 0.0000 0
Pavlovian hp ACh -0.0202 -0.0004 0.0269 0.0000 0
control hp DA -0.0097 0.0000 0.0097 0.0000 19
control hp ACh -0.0088 0.0000 0.0088 0.0000 0
ADS raw DA -0.0218 0.0000 0.0398 0.0000 0
ADS raw ACh -0.0194 0.0000 0.0268 0.0000 0
Wheel velocity (m/s) in spontaneous sessions, per mouse (median over sessions):
min median p99 max frac_above_0.0208
mouse
UG27 -0.1467 0.0136 0.4668 0.5748 0.4870
UG28 -0.1774 0.0000 0.1332 0.2442 0.1133
UG29 -0.1035 0.0000 0.1434 0.2771 0.1385
UG30 -0.1348 0.1952 0.4405 0.4869 0.7895
UG31 -0.1051 0.0000 0.1383 0.3316 0.1209
The dF/F values are small fractions, as dF/F should be (1st to 99th percentiles between -0.033 and 0.040). The artifact-removed series holds NaNs where artifacts were cut. Some control sites have a constant red trace in some sessions; their correlation is undefined there and is averaged over the other sessions. Wheel velocity is in m/s with plausible values (99th percentiles 0.13 to 0.47 m/s). The mice differ in how much they run: the smoothed velocity exceeds the stillness threshold 11% to 79% of the time.
pav_table["has_hp_art"] = [any("HpArt" in k for k in e["series"]) for e in ad_events]
print("AD files holding an artifact-removed (HpArt) series:", pav_table.groupby("mouse").has_hp_art.sum().to_dict(),
"of", pav_table.groupby("mouse").size().to_dict())
print(pav_table.groupby(["mouse", "events", "light_rewarded"], dropna=False).size().rename("sessions").to_string())
print("\nSelected Pavlovian sessions (last 3 per phase, as the code takes):")
print(pd.DataFrame([{"mouse": s["subject"], "session": s["session_id"], "phase": s["phase"],
"light_cues": len(s["events"]["Light"][0]), "rewards": len(s["events"]["Reward"][0])}
for s in pav]).to_string(index=False))
print("\nControl sessions used:", len(ctrl), {m: len(v) for m, v in daach.by_mouse(ctrl).items()})
stim_counts = pd.DataFrame([{"mouse": s["subject"], "session": s["session_id"],
"stimulations": len(s["events"].get("Optogenetic Stimulation", ([],))[0]),
"rewards": len(s["events"].get("Reward", ([],))[0])} for s in ads])
print("\nStimulation sessions:")
print(stim_counts.to_string(index=False))
AD files holding an artifact-removed (HpArt) series: {'AD1': 0, 'AD2': 0, 'AD3': 0, 'AD4': 0, 'AD5': 0, 'AD6': 0} of {'AD1': 38, 'AD2': 26, 'AD3': 33, 'AD4': 29, 'AD5': 41, 'AD6': 44}
mouse events light_rewarded
AD1 Lick,Light,Reward 0.2 14
0.8 12
1.0 12
AD2 Lick,Light,Reward 0.2 6
0.8 8
1.0 12
AD3 Lick,Light,Reward 0.2 11
0.8 10
1.0 12
AD4 Lick,Light,Reward 0.2 3
0.8 6
1.0 15
Lick,Reward NaN 5
AD5 Lick,Light,Reward 0.2 7
0.8 7
1.0 22
Lick,Reward NaN 5
AD6 Lick,Light,Reward 0.2 9
0.8 10
1.0 20
Lick,Reward NaN 5
Selected Pavlovian sessions (last 3 per phase, as the code takes):
mouse session phase light_cues rewards
AD1 240411 learning 30 68
AD1 240412 learning 30 68
AD1 240416 learning 30 68
AD1 240502 extinction 30 38
AD1 240503 extinction 30 38
AD1 240506 extinction 30 38
AD2 240411 learning 30 68
AD2 240412 learning 30 68
AD2 240416 learning 30 68
AD2 240422 extinction 30 38
AD2 240423 extinction 30 38
AD2 240424 extinction 30 38
AD3 240411 learning 30 68
AD3 240412 learning 30 68
AD3 240416 learning 30 68
AD3 240429 extinction 30 38
AD3 240430 extinction 30 38
AD3 240501 extinction 30 38
Control sessions used: 15 {'AD4': 5, 'AD5': 5, 'AD6': 5}
Stimulation sessions:
mouse session stimulations rewards
ADS12 250702 13 16
ADS12 250703 13 16
ADS12 250708 13 16
ADS12 250712 13 16
ADS13 250702 13 16
ADS13 250703 13 16
ADS13 250708 13 16
ADS13 250712 13 16
ADS16 250702 13 16
ADS16 250703 13 16
ADS16 250708 13 16
ADS16 250712 13 16
ADS17 250702 13 16
ADS17 250703 13 16
ADS17 250708 13 16
ADS17 250712 13 16
ADS19 250702 13 16
ADS19 250703 13 16
ADS19 250708 13 16
ADS19 250712 13 16
ADS6 250328 10 10
ADS6 250401 10 10
ADS6 250402 10 10
Pavlovian sessions fall into three groups by the fraction of light cues followed by a reward: 1.0 (learning), 0.2 (light cue extinguished) and 0.8 (tone cue extinguished). Learning sessions have 68 rewards: 30 after light cues, 30 after tone cues and 8 uncued, as the Methods describe. The files hold no tone cue events, so tone trials cannot be analysed; the paper's Figure 4 uses light trials only, as its code shows. Which phase a session belongs to is not stated in the files and is inferred here from the events.
def fingerprint(s):
return (len(s["DA"]), tuple(np.round(np.nan_to_num(s["DA"][:200, :3]).ravel(), 6)))
seen, dups = {}, []
for s in ug_ok + pav + ctrl + ads:
k = fingerprint(s)
if k in seen:
dups.append((seen[k], s["path"]))
seen[k] = s["path"]
print("Sessions with identical length and identical first 200 DA samples at 3 sites:", dups if dups else "none")
same_len = pd.Series([len(s["DA"]) for s in ug_ok]).value_counts()
print("UG sessions sharing a sample count with another:", int((same_len > 1).sum()), "counts;",
"none share traces (above)")
Sessions with identical length and identical first 200 DA samples at 3 sites: none UG sessions sharing a sample count with another: 23 counts; none share traces (above)
def response_contrast(sessions, events_of, shift_s=10.0):
# mean over sites of |mean dF/F 0-1 s minus -0.5-0 s| at the events and at events shifted later
out = []
for m, ss in daach.by_mouse(sessions).items():
vals = {"true": [], "shifted": []}
for s in ss:
e = events_of(s)
for key, ev in (("true", e), ("shifted", e + int(shift_s * daach.SR))):
ev = ev[(ev > 9) & (ev < len(s["DA"]) - 19)]
if not len(ev):
continue
for ch in ("DA", "ACh"):
r = daach.raster(s[ch], ev, -9, 18)
mu = np.nanmean(r, 2)
vals[key].append(np.nanmean(np.abs(mu[9:].mean(0) - mu[:9].mean(0))))
out.append({"mouse": m, "true": np.mean(vals["true"]), "shifted_10s": np.mean(vals["shifted"])})
return pd.DataFrame(out)
checks = {
"reward sessions, first lick": (reward, lambda s: daach.first_licks(s)[1]),
"Pavlovian, light cue": (pav, lambda s: daach.event_index(s, s["events"]["Light"][0])),
"ADS, stimulation": (ads, lambda s: daach.event_index(s, s["events"]["Optogenetic Stimulation"][0])),
"spontaneous, initiation": (spont, lambda s: daach.initiations(daach.smooth_velocity(s["velocity"]))[0]),
}
for name, (ss, f) in checks.items():
t = response_contrast(ss, f)
t["ratio"] = t.true / t.shifted_10s
print(name)
print(t.round(4).to_string(index=False))
reward sessions, first lick mouse true shifted_10s ratio UG27 0.0150 0.0013 11.8339 UG28 0.0124 0.0012 10.0586 UG29 0.0115 0.0013 9.1609 UG30 0.0179 0.0017 10.5749 UG31 0.0146 0.0015 9.9095
Pavlovian, light cue mouse true shifted_10s ratio AD1 0.0049 0.0014 3.4238 AD2 0.0029 0.0013 2.3267 AD3 0.0054 0.0010 5.5923 ADS, stimulation mouse true shifted_10s ratio ADS12 0.0211 0.0022 9.3879 ADS13 0.0084 0.0017 5.0192 ADS16 0.0134 0.0014 9.3411 ADS17 0.0220 0.0014 15.7039 ADS19 0.0118 0.0017 6.7759 ADS6 0.0122 0.0013 9.1231
spontaneous, initiation mouse true shifted_10s ratio UG27 0.0039 0.0013 2.9552 UG28 0.0042 0.0033 1.2796 UG29 0.0030 0.0020 1.4550 UG30 0.0044 0.0026 1.6993 UG31 0.0065 0.0026 2.5080
For every mouse and each kind of event, the photometry response in the second after the event is larger at the event times than at times shifted by 10 s (the ratios above), so the events, the velocity and the photometry share one clock in every cohort. This shows alignment to within the response's time scale (about a second); finer offsets would not show here.
fig1e-dominant-polarity¶
Across 174 sites in 5 mice, the dominant DA-ACh cross-correlation during spontaneous behavior is negative at 132 sites, 127 of them with DA leading, and positive at 42 sites, 31 of them with ACh leading.
xc1 = daach.spontaneous_xcorr(spont)
mice1, coords1 = daach.site_table(spont)
dom1, lag1 = xc1["dominant"], xc1["lag"]
def polarity_counts(dom, lag):
return {"negative": int((dom < 0).sum()), "negative, DA leads": int(((dom < 0) & (lag > 0)).sum()),
"positive": int((dom > 0).sum()), "positive, ACh leads": int(((dom > 0) & (lag < 0)).sum()),
"lag 0": int((lag == 0).sum()), "sites": len(dom)}
expected = {"negative": 132, "negative, DA leads": 127, "positive": 42, "positive, ACh leads": 31}
c1 = polarity_counts(dom1, lag1)
print("Counts:", c1)
print("Within 10% of the paper:", {k: abs(c1[k] - v) <= 0.1 * v for k, v in expected.items()})
per = pd.DataFrame([{"mouse": m, **polarity_counts(dom1[mice1 == m], lag1[mice1 == m])} for m in np.unique(mice1)])
print(per.to_string(index=False))
print("Dominant r range", dom1.min().round(3), "to", dom1.max().round(3), "; latency range (ms)",
lag1.min() * 1000 / 18, "to", lag1.max() * 1000 / 18)
xc1c = daach.spontaneous_xcorr(spont, concatenate=True)
print("Variant, correlation over concatenated sessions:", polarity_counts(xc1c["dominant"], xc1c["lag"]))
Counts: {'negative': 133, 'negative, DA leads': 128, 'positive': 41, 'positive, ACh leads': 32, 'lag 0': 6, 'sites': 174}
Within 10% of the paper: {'negative': True, 'negative, DA leads': True, 'positive': True, 'positive, ACh leads': True}
mouse negative negative, DA leads positive positive, ACh leads lag 0 sites
UG27 19 19 4 4 0 23
UG28 33 33 8 4 3 41
UG29 27 23 21 16 3 48
UG30 25 24 3 3 0 28
UG31 29 29 5 5 0 34
Dominant r range -0.525 to 0.399 ; latency range (ms) -777.7777777777778 to 388.8888888888889
Variant, correlation over concatenated sessions: {'negative': 133, 'negative, DA leads': 128, 'positive': 41, 'positive, ACh leads': 30, 'lag 0': 7, 'sites': 174}
fig, ax = plt.subplots(figsize=(4.5, 3.5))
ax.scatter(lag1 * 1000 / 18, dom1, s=10, c=np.where(dom1 < 0, "tab:blue", "tab:red"))
ax.axhline(0, color="k", lw=0.5)
ax.axvline(0, color="k", lw=0.5)
ax.set_xlabel("latency (ms), ACh leads <- -> DA leads")
ax.set_ylabel("dominant Pearson r")
ax.set_title("Spontaneous DA-ACh, 174 sites")
plt.show()
print("All points are inside the axes (latency -1000 to 1000 ms).")
All points are inside the axes (latency -1000 to 1000 ms).
Result: reproduced. Of 174 sites, the dominant correlation is negative at 133 (128 with DA leading) and positive at 41 (32 with ACh leading), against 132, 127, 42 and 31 in the paper; all four are within 10%. Every mouse has more negative than positive sites. Correlating concatenated sessions gives the same counts except 30 positive sites with ACh leading.
fig1f-adls-hotspot¶
The dominant correlations form a significant hotspot, identified by local Moran's I, in the anterior dorsolateral striatum (aDLS), where the strongest correlations are negative.
The hotspot procedure is fixed in the claims file: linear interpolation onto a 100 um grid, padding, Gaussian smoothing, local Moran's I against 200 shuffles of the values across sites, the 95th percentile, and the connected volume of at least 0.125 mm3 with the largest mean magnitude. One grid is shared by all cohorts so that sites of other cohorts can be placed in this hotspot.
mice_c, coords_c = daach.site_table(ctrl)
mice_p, coords_p = daach.site_table(pav)
mice_s, coords_s = daach.site_table(ads)
GRID = daach.Grid(np.concatenate([coords1, coords_c, coords_p, coords_s]))
interp1 = daach.Interpolator(GRID, coords1)
rng = np.random.default_rng(101)
hs1 = daach.hotspot(interp1, dom1, rng)
in1 = daach.inside(GRID, hs1, coords1)
med = np.median(coords1, 0)
print("Grid (AP, ML, DV voxels):", GRID.shape, "; voxels inside the convex hull of the sites:", int(interp1.inside.sum()))
print("Candidate components:", hs1["n_components"], "; hotspot voxels:", len(hs1["voxels"]),
f"({len(hs1['voxels']) * 1e-3:.3f} mm3)")
print(f"Hotspot mean smoothed r {hs1['mean']:.3f}; striatal (hull) mean {hs1['striatal_mean']:.3f}")
print("Centroid (AP, ML, DV mm):", hs1["centroid"].round(2), "; median site coordinates:", med.round(2))
print("Anterior, lateral, dorsal of the median:", list(hs1["centroid"] > med))
print("Sites inside:", int(in1.sum()), "by mouse:", pd.Series(mice1[in1]).value_counts().to_dict())
for pct, min_vox in ((99, 125), (95, 1000), (90, 125)):
v = daach.hotspot(interp1, dom1, rng, pct=pct, min_vox=min_vox, null=hs1["null"])
print(f"Variant, null percentile {pct}, minimum {min_vox * 1e-3:.3f} mm3: components {v['n_components']}, "
f"mean {v['mean']:.3f}, voxels {0 if v['voxels'] is None else len(v['voxels'])}, centroid {v['centroid'].round(2)}")
shuf1 = pd.DataFrame(daach.hotspot_shuffles(interp1, dom1, hs1["null"], np.random.default_rng(103), n=100))
neg1 = shuf1.found & (shuf1["mean"] < shuf1.striatal_mean)
print(f"Same procedure on dominant r shuffled across sites (100 shuffles, same null): hotspot found in {int(shuf1.found.sum())}, "
f"with a mean below the striatal mean in {int(neg1.sum())}; voxels median {shuf1.voxels[shuf1.found].median():.0f}, "
f"max {shuf1.voxels.max()}; at least as large as the real one ({len(hs1['voxels'])} voxels) in {int((shuf1.voxels >= len(hs1['voxels'])).sum())}; "
f"mean at or below {hs1['mean']:.3f} in {int((shuf1['mean'] <= hs1['mean']).sum())}")
Grid (AP, ML, DV voxels): (51, 46, 51) ; voxels inside the convex hull of the sites: 6155
Candidate components: 2 ; hotspot voxels: 1010 (1.010 mm3)
Hotspot mean smoothed r -0.368; striatal (hull) mean -0.185
Centroid (AP, ML, DV mm): [ 0.63 2.02 -3.14] ; median site coordinates: [ 0.1 2.04 -3.21]
Anterior, lateral, dorsal of the median: [np.True_, np.False_, np.True_]
Sites inside: 28 by mouse: {'UG31': 10, 'UG28': 8, 'UG30': 5, 'UG27': 3, 'UG29': 2}
Variant, null percentile 99, minimum 0.125 mm3: components 1, mean -0.388, voxels 518, centroid [ 0.64 2.09 -3.11]
Variant, null percentile 95, minimum 1.000 mm3: components 1, mean -0.368, voxels 1010, centroid [ 0.63 2.02 -3.14]
Variant, null percentile 90, minimum 0.125 mm3: components 2, mean -0.356, voxels 1313, centroid [ 0.63 1.98 -3.16]
Same procedure on dominant r shuffled across sites (100 shuffles, same null): hotspot found in 67, with a mean below the striatal mean in 55; voxels median 212, max 482; at least as large as the real one (1010 voxels) in 0; mean at or below -0.368 in 0
def show_map(hs, title, vmax=0.5):
vol = hs["smooth"]
mask = np.zeros(GRID.shape, bool)
if hs["voxels"] is not None:
mask.ravel()[hs["voxels"]] = True
fig, axs = plt.subplots(1, 2, figsize=(8, 3.5))
ap, ml, dv = GRID.axes
axs[0].imshow(np.nanmean(vol, 2), origin="lower", cmap="RdBu_r", vmin=-vmax, vmax=vmax, aspect="auto",
extent=[ml[0], ml[-1], ap[0], ap[-1]])
axs[0].contour(ml, ap, mask.any(2), levels=[0.5], colors="k")
axs[0].set_xlabel("ML (mm)"); axs[0].set_ylabel("AP (mm)"); axs[0].set_title("axial, mean projection")
im = axs[1].imshow(np.nanmean(vol, 1).T, origin="lower", cmap="RdBu_r", vmin=-vmax, vmax=vmax, aspect="auto",
extent=[ap[0], ap[-1], dv[0], dv[-1]])
axs[1].contour(ap, dv, mask.any(1).T, levels=[0.5], colors="k")
axs[1].set_xlabel("AP (mm)"); axs[1].set_ylabel("DV (mm)"); axs[1].set_title("sagittal, mean projection")
axs[1].invert_xaxis()
fig.colorbar(im, ax=axs, shrink=0.8)
fig.suptitle(title)
plt.show()
print(f"Smoothed values range {np.nanmin(vol):.3f} to {np.nanmax(vol):.3f}; color scale clipped at +/-{vmax}; "
f"voxels beyond it: {int((np.abs(vol) > vmax).sum())} of {int((~np.isnan(vol)).sum())}")
show_map(hs1, "Dominant spontaneous DA-ACh r; black, hotspot")
Smoothed values range -0.473 to 0.105; color scale clipped at +/-0.5; voxels beyond it: 0 of 6155
Result: partially reproduced. A significant hotspot is found (1.01 mm3) with a negative mean smoothed r (-0.368, against -0.185 over the hull). Its centroid (AP 0.63, ML 2.02, DV -3.14 mm) is anterior and dorsal of the median site (0.10, 2.04, -3.21) but 0.02 mm medial of it, so by the criterion fixed beforehand it is not in the lateral half. With the 99th percentile of the null the hotspot shrinks to 518 voxels and its centroid moves to ML 2.09 mm, lateral of the median. Sites of all five mice fall inside it (28 sites). Finding a hotspot is not evidence on its own: with the values shuffled across sites the procedure finds one in 67 of 100 shuffles. None of those is as large as the real one (at most 482 voxels against 1010) or has a mean as low (-0.368).
fig1h-in-vs-out¶
Sites inside the aDLS hotspot have more negative dominant correlations than sites outside (Wilcoxon rank-sum p = 1.67e-16).
r = daach.in_out(dom1, in1)
print(r)
print("Paper: p = 1.67e-16. The test is partly circular: the hotspot is drawn from the same values.")
{'n_in': 28, 'n_out': 146, 'median_in': np.float64(-0.40661912524303806), 'median_out': np.float64(-0.13754268545237336), 'p': np.float64(1.8707508174182813e-13), 'p_exact': nan}
Paper: p = 1.67e-16. The test is partly circular: the hotspot is drawn from the same values.
Result: reproduced. The 28 sites inside the hotspot have a median dominant r of -0.41 against -0.14 for the 146 outside (rank-sum p = 1.9e-13; the paper 1.67e-16).
fig1f-positive-posterior¶
Dominant positive correlations lie in the posterior striatum.
ap = coords1[:, 0]
t = stats.mannwhitneyu(ap[dom1 > 0], ap[dom1 < 0], alternative="less")
print(f"AP of positive sites: median {np.median(ap[dom1 > 0]):.2f} mm (n {int((dom1 > 0).sum())}); "
f"negative sites: median {np.median(ap[dom1 < 0]):.2f} mm (n {int((dom1 < 0).sum())}); one-sided p = {t.pvalue:.3g}")
for m in np.unique(mice1):
k = mice1 == m
a, b = ap[k & (dom1 > 0)], ap[k & (dom1 < 0)]
p = stats.mannwhitneyu(a, b, alternative="less").pvalue if len(a) and len(b) else np.nan
print(f" {m}: positive n {len(a)}, median AP {np.median(a) if len(a) else np.nan:.2f}; negative n {len(b)}, "
f"median AP {np.median(b):.2f}; p = {p:.3g}")
print("Variant, positive sites against all others (DV and ML medians too):",
{"ML_pos": np.median(coords1[dom1 > 0, 1]).round(2), "ML_neg": np.median(coords1[dom1 < 0, 1]).round(2),
"DV_pos": np.median(coords1[dom1 > 0, 2]).round(2), "DV_neg": np.median(coords1[dom1 < 0, 2]).round(2)})
AP of positive sites: median -0.40 mm (n 41); negative sites: median 0.24 mm (n 133); one-sided p = 4.28e-07
UG27: positive n 4, median AP -0.67; negative n 19, median AP 0.12; p = 0.00599
UG28: positive n 8, median AP -0.17; negative n 33, median AP 0.38; p = 0.178
UG29: positive n 21, median AP -0.44; negative n 27, median AP 0.46; p = 6.28e-05
UG30: positive n 3, median AP 0.10; negative n 25, median AP 0.34; p = 0.117
UG31: positive n 5, median AP -0.94; negative n 29, median AP 0.01; p = 0.0206
Variant, positive sites against all others (DV and ML medians too): {'ML_pos': np.float64(2.03), 'ML_neg': np.float64(2.04), 'DV_pos': np.float64(-3.43), 'DV_neg': np.float64(-3.16)}
Result: reproduced. Sites with a positive dominant correlation lie more posterior (median AP -0.40 mm, 41 sites) than sites with a negative one (0.24 mm, 133 sites), one-sided p = 4.3e-07. All five mice show the same direction, significant in three (UG27, UG29, UG31).
ed2-hotspot-per-mouse¶
The hotspot is consistent across mice; in each mouse, sites inside the group hotspot have more negative dominant correlations than sites outside.
rows = []
for m in np.unique(mice1):
k = mice1 == m
rows.append({"mouse": m, **daach.in_out(dom1[k], in1[k])})
per2 = pd.DataFrame(rows)
print(per2.round(4).to_string(index=False))
ok = per2[(per2.n_in >= 2) & (per2.n_out >= 2)]
print("Mice tested (>= 2 sites each side):", len(ok), "; median inside lower:", int((ok.median_in < ok.median_out).sum()),
"; and p < 0.05:", int(((ok.median_in < ok.median_out) & (ok.p < 0.05)).sum()))
mouse n_in n_out median_in median_out p p_exact UG27 3 20 -0.4527 -0.2647 0.0071 0.0011 UG28 8 33 -0.3836 -0.1124 0.0000 0.0000 UG29 2 46 -0.3172 -0.0288 0.0839 0.0869 UG30 5 23 -0.4201 -0.2564 0.0069 0.0040 UG31 10 24 -0.4299 -0.1477 0.0000 0.0000 Mice tested (>= 2 sites each side): 5 ; median inside lower: 5 ; and p < 0.05: 4
Result: reproduced. In all five mice the median dominant r is lower inside the group hotspot than outside, with p < 0.05 in four (UG27, UG28, UG30, UG31); UG29 has 2 sites inside (p = 0.084, exact p = 0.087).
ed1-controls-no-organization¶
Mice expressing non-functional sensors (n = 3) show only sparse, weak positive red-green correlations with no anatomical organization; no hotspot is found.
xcc = daach.spontaneous_xcorr(ctrl, mask=daach.reward_mask)
domc = xcc["dominant"]
keep = ~np.isnan(domc)
print("Control sites:", len(domc), "; with no defined correlation in any session (dropped):", int((~keep).sum()))
domc, mice_ck, coords_ck = domc[keep], mice_c[keep], coords_c[keep]
print(f"Control dominant r: median {np.median(domc):.3f}, positive at {(domc > 0).mean():.2f} of sites, "
f"median |r| {np.median(np.abs(domc)):.3f}; sensor mice median |r| {np.median(np.abs(dom1)):.3f}")
print(f"|r| > 0.2: control {(np.abs(domc) > 0.2).mean():.2f}, sensor {(np.abs(dom1) > 0.2).mean():.2f}")
print(pd.DataFrame([{"mouse": m, "sites": int((mice_ck == m).sum()), "median_r": np.median(domc[mice_ck == m]),
"frac_positive": (domc[mice_ck == m] > 0).mean()} for m in np.unique(mice_ck)]).round(3).to_string(index=False))
interp_c = daach.Interpolator(GRID, coords_ck)
hsc = daach.hotspot(interp_c, domc, np.random.default_rng(102))
print(f"Hotspot procedure: components {hsc['n_components']}, hotspot mean {hsc['mean']:.3f}, "
f"voxels {0 if hsc['voxels'] is None else len(hsc['voxels'])}, centroid {hsc['centroid'].round(2)}")
xcc2 = daach.spontaneous_xcorr(ctrl)
d2 = xcc2["dominant"][~np.isnan(xcc2["dominant"])]
print(f"Variant, reward periods kept: median r {np.median(d2):.3f}, positive {(d2 > 0).mean():.2f}")
shufc = pd.DataFrame(daach.hotspot_shuffles(interp_c, domc, hsc["null"], np.random.default_rng(103), n=100))
print(f"Control values shuffled across sites (100 shuffles): hotspot found in {int(shufc.found.sum())}; voxels median "
f"{shufc.voxels[shufc.found].median():.0f}; at least {len(hsc['voxels'])} voxels in {int((shufc.voxels >= len(hsc['voxels'])).sum())}; "
f"mean at or above {hsc['mean']:.3f} in {int((shufc['mean'] >= hsc['mean']).sum())}")
Control sites: 125 ; with no defined correlation in any session (dropped): 1 Control dominant r: median 0.202, positive at 0.68 of sites, median |r| 0.238; sensor mice median |r| 0.240 |r| > 0.2: control 0.65, sensor 0.54 mouse sites median_r frac_positive AD4 42 0.245 0.738 AD5 41 0.096 0.561 AD6 41 0.208 0.732
Hotspot procedure: components 2, hotspot mean 0.329, voxels 551, centroid [ 0.68 1.29 -4.14]
Variant, reward periods kept: median r 0.221, positive 0.75
Control values shuffled across sites (100 shuffles): hotspot found in 73; voxels median 208; at least 551 voxels in 0; mean at or above 0.329 in 31
Result: not reproduced. In the control mice the hotspot procedure finds a hotspot (mean smoothed r +0.33, 551 voxels, centroid AP 0.68, ML 1.29, DV -4.14 mm). The procedure finds a hotspot in 73 of 100 shuffles of the control values, but none of those reaches 551 voxels (31 have a mean of at least +0.33), so by the paper's own test the control correlations are organized in space. The correlations are mostly positive (68% of 124 sites, median r 0.20) but not weak: the median |r| is 0.238, as in the sensor mice (0.240), and 65% of control sites exceed |r| = 0.2 against 54% of sensor sites. Two differences from the paper's method could add correlated red-green signal here: the control files hold no artifact-removed series, so the high-pass series is used without artifact removal, and the test used the 15 uncued reward sessions of these mice, not all of their ITI data.
fig2a-directionality¶
For coincident transients of opposite polarity the DA transient leads; for same-polarity pairs the ACh transient leads.
t0 = time.time()
tp = pd.DataFrame(daach.transient_pairs(ug_ok, np.random.default_rng(201)))
print(f"Transients and pairs computed over {len(ug_ok)} sessions in {time.time() - t0:.0f} s; sites {len(tp)}")
names = [daach.pair_name(t) for t in daach.PAIR_TYPES]
rows = []
for k in names:
di = tp[f"di|{k}"].dropna()
rows.append({"pair": k, "sites": len(di), "median_DI": di.median(), "frac_negative": (di < 0).mean(),
"wilcoxon_p": stats.wilcoxon(di).pvalue})
print(pd.DataFrame(rows).to_string(index=False))
print("\nMedian directionality index by mouse (negative: DA leads):")
print(tp.groupby("mouse")[[f"di|{k}" for k in names]].median().round(3).to_string())
Transients and pairs computed over 104 sessions in 58 s; sites 174
pair sites median_DI frac_negative wilcoxon_p
ACh trough / DA peak 174 -0.267096 0.913793 2.508659e-28
ACh peak / DA trough 174 -0.298606 0.936782 4.311079e-29
ACh peak / DA peak 174 0.332209 0.028736 6.832557e-30
ACh trough / DA trough 174 0.119870 0.126437 8.288437e-25
Median directionality index by mouse (negative: DA leads):
di|ACh trough / DA peak di|ACh peak / DA trough di|ACh peak / DA peak di|ACh trough / DA trough
mouse
UG27 -0.316 -0.389 0.311 0.150
UG28 -0.284 -0.284 0.354 0.125
UG29 -0.146 -0.191 0.253 0.066
UG30 -0.306 -0.407 0.393 0.201
UG31 -0.276 -0.315 0.331 0.116
Result: reproduced. The median directionality index is -0.27 (ACh trough with DA peak) and -0.30 (ACh peak with DA trough), DA leading, and +0.33 (peak with peak) and +0.12 (trough with trough), ACh leading, each with Wilcoxon p below 1e-24 over 174 sites. Each mouse has the same four signs.
fig2-magnitude-correlations¶
Magnitudes of paired transients are negatively correlated at many sites for opposite-sign pairs (121 and 81 of 174) and positively for same-sign pairs (97 and 78 of 174), with no site significant in the other direction.
exp = {names[0]: ("negative", 121), names[1]: ("negative", 81), names[2]: ("positive", 97), names[3]: ("positive", 78)}
rows = []
for k in names:
neg, pos = int((tp[f"hi|{k}"] < 0).sum()), int((tp[f"lo|{k}"] > 0).sum())
sign, e = exp[k]
got = neg if sign == "negative" else pos
code_neg = int(((tp[f"code_p|{k}"] < 0.05) & (tp[f"code_r|{k}"] < 0)).sum())
code_pos = int(((tp[f"code_p|{k}"] < 0.05) & (tp[f"code_r|{k}"] > 0)).sum())
rows.append({"pair": k, "sig_negative": neg, "sig_positive": pos, "expected": f"{e} {sign}",
"within_10pct": abs(got - e) <= 0.1 * e, "sites_<100_pairs": int((tp[f"npairs|{k}"] < 100).sum()),
"median_pairs": int(tp[f"npairs|{k}"].median()), "median_mean_r": tp[f"mag|{k}"].median(),
"code_variant_neg": code_neg, "code_variant_pos": code_pos})
print(pd.DataFrame(rows).round(3).to_string(index=False))
print("\nSignificant sites of the stated sign, by mouse:")
by = {}
for k in names:
col = (tp[f"hi|{k}"] < 0) if exp[k][0] == "negative" else (tp[f"lo|{k}"] > 0)
by[k] = col.groupby(tp.mouse).sum()
print(pd.DataFrame(by).to_string())
pair sig_negative sig_positive expected within_10pct sites_<100_pairs median_pairs median_mean_r code_variant_neg code_variant_pos
ACh trough / DA peak 114 0 121 negative True 0 1293 -0.244 149 1
ACh peak / DA trough 69 0 81 negative False 0 1454 -0.170 144 0
ACh peak / DA peak 0 97 97 positive True 0 1151 0.211 0 148
ACh trough / DA trough 0 78 78 positive True 0 1483 0.175 0 142
Significant sites of the stated sign, by mouse:
ACh trough / DA peak ACh peak / DA trough ACh peak / DA peak ACh trough / DA trough
mouse
UG27 5 6 6 5
UG28 32 18 27 26
UG29 23 18 25 20
UG30 28 22 23 22
UG31 26 5 16 5
Result: partially reproduced. Significant sites of the stated sign are 114 (ACh trough with DA peak), 69 (ACh peak with DA trough), 97 (peak with peak) and 78 (trough with trough), against 121, 81, 97 and 78; three are within 10%, and the second is 3.9 below its lower bound (72.9). No site is significant in the other direction for any type, and every site has at least 100 pairs. The counts depend on the random draws by several sites: the first run, with a different random stream, gave 115, 70, 90 and 76. The code's computation, Pearson r over all pairs of the leading order, finds 142 to 149 sites per type.
fig2d-magnitude-in-out¶
Magnitude correlations of opposite-sign pairs, but not same-sign pairs, are stronger inside the aDLS hotspot than outside.
in_tp = in1 # transient_pairs returns sites in the same order as site_table
rows = []
for k in names:
v = np.where(tp[f"npairs|{k}"] >= 100, tp[f"mag|{k}"], np.nan)
rows.append({"pair": k, **daach.in_out(v, in_tp)})
print(pd.DataFrame(rows).round(4).to_string(index=False))
print("Paper: p = 0.017 and 0.006 for the opposite-sign pairs, 0.091 and 0.612 for the same-sign pairs.")
pair n_in n_out median_in median_out p p_exact
ACh trough / DA peak 28 146 -0.3252 -0.2323 0.0001 NaN
ACh peak / DA trough 28 146 -0.1861 -0.1626 0.1134 NaN
ACh peak / DA peak 28 146 0.2450 0.1984 0.1079 NaN
ACh trough / DA trough 28 146 0.1821 0.1725 0.2169 NaN
Paper: p = 0.017 and 0.006 for the opposite-sign pairs, 0.091 and 0.612 for the same-sign pairs.
Result: partially reproduced. For ACh troughs with DA peaks the magnitude correlation is more negative inside the hotspot (median -0.33 against -0.23, p = 0.0001); for ACh peaks with DA troughs the difference is small and not significant (-0.19 against -0.16, p = 0.11). Neither same-sign type differs (p = 0.11 and 0.22), as the paper states.
fig2b-pairing-overlaps-hotspot¶
The frequency of paired transients is highest in aDLS, in a region overlapping the correlation hotspot, for all four pair types (volume overlap test, p < 0.01 for each).
rng = np.random.default_rng(202)
rows = []
for k in names:
occ = tp[f"occ|{k}"].to_numpy()
h = daach.hotspot(interp1, occ, rng, how="max")
ov = daach.overlap_test(hs1, interp1.inside, h, interp1.inside, GRID.shape, rng)
io = daach.in_out(occ, in1)
rows.append({"pair": k, "median_occurrence": np.median(occ), "hotspot_mean": h["mean"],
"striatal_mean": h["striatal_mean"], "hotspot_voxels": ov["n_b"], "overlap_voxels": ov["overlap"],
"null_median_overlap": ov.get("null_median"), "overlap_p": ov["p"],
"centroid_AP_ML_DV": h["centroid"].round(2), "in_median": io["median_in"], "out_median": io["median_out"],
"in_out_p": io["p"]})
print(pd.DataFrame(rows).round(4).to_string(index=False))
print("Paper overlap p: 0.005, 0.005, 0.002, < 0.001 (10000 random volumes; here 200, so the smallest possible p is 0.005).")
print("\nChance-corrected occurrence (actual minus the null of DA and ACh from different sessions, 10 shuffles):")
rng = np.random.default_rng(203)
rows = []
for k in names:
exc = (tp[f"occ|{k}"] - tp[f"null_occ|{k}"]).to_numpy()
h = daach.hotspot(interp1, exc, rng, how="max")
ov = daach.overlap_test(hs1, interp1.inside, h, interp1.inside, GRID.shape, rng)
io = daach.in_out(exc, in1)
rows.append({"pair": k, "median_null": tp[f"null_occ|{k}"].median(), "median_excess": np.median(exc),
"frac_sites_excess>0": (exc > 0).mean(), "hotspot_mean": h["mean"], "striatal_mean": h["striatal_mean"],
"hotspot_voxels": ov["n_b"], "overlap_voxels": ov["overlap"], "null_median_overlap": ov.get("null_median"),
"overlap_p": ov["p"], "in_median": io["median_in"], "out_median": io["median_out"], "in_out_p": io["p"]})
print(pd.DataFrame(rows).round(4).to_string(index=False))
pair median_occurrence hotspot_mean striatal_mean hotspot_voxels overlap_voxels null_median_overlap overlap_p centroid_AP_ML_DV in_median out_median in_out_p
ACh trough / DA peak 0.3589 0.4273 0.3394 975 780 96.0 0.005 [0.59, 2.03, -3.13] 0.4311 0.3312 0.0
ACh peak / DA trough 0.3738 0.4570 0.3632 1133 826 98.5 0.015 [0.64, 1.87, -3.14] 0.4752 0.3403 0.0
ACh peak / DA peak 0.3450 0.4200 0.3325 943 782 88.5 0.005 [0.62, 2.01, -3.14] 0.4289 0.3163 0.0
ACh trough / DA trough 0.3448 0.4703 0.3487 954 815 111.0 0.000 [0.63, 2.05, -3.17] 0.4968 0.3169 0.0
Paper overlap p: 0.005, 0.005, 0.002, < 0.001 (10000 random volumes; here 200, so the smallest possible p is 0.005).
Chance-corrected occurrence (actual minus the null of DA and ACh from different sessions, 10 shuffles):
pair median_null median_excess frac_sites_excess>0 hotspot_mean striatal_mean hotspot_voxels overlap_voxels null_median_overlap overlap_p in_median out_median in_out_p
ACh trough / DA peak 0.2542 0.1082 1.0 0.1487 0.1059 802 606 63.5 0.030 0.1408 0.1036 0.0
ACh peak / DA trough 0.2707 0.1084 1.0 0.1613 0.1116 1122 802 112.5 0.005 0.1724 0.0943 0.0
ACh peak / DA peak 0.2495 0.0916 1.0 0.1494 0.0954 701 566 31.0 0.015 0.1466 0.0815 0.0
ACh trough / DA trough 0.2720 0.0842 1.0 0.1663 0.0983 886 764 67.5 0.000 0.1761 0.0705 0.0
Result: reproduced. For each pair type the occurrence hotspot shares 780 to 826 voxels with the correlation hotspot, against null medians of 88 to 111, with overlap p of 0.005, 0.015, 0.005 and below 0.005 (the paper 0.005, 0.005, 0.002 and below 0.001). The occurrence hotspots have centroids at AP 0.59 to 0.64, ML 1.87 to 2.05 and DV -3.13 to -3.17 mm, and occurrence is higher inside the correlation hotspot than outside for each type (rank-sum p < 0.001). After subtracting the chance occurrence (DA and ACh transients from different sessions of the same site, 10 shuffles), every site keeps an excess for every type (medians 0.08 to 0.11), the excess is higher inside the correlation hotspot (p < 0.001), and its hotspots overlap the correlation hotspot with p of 0.030, 0.005, 0.015 and below 0.005.
ed3e-transient-rates-in-hotspot¶
DA and ACh transient rates are higher at sites inside the global correlation hotspot than outside.
rows = []
for k in ("DA_pos", "DA_neg", "ACh_pos", "ACh_neg"):
rows.append({"transients": k, **daach.in_out(tp[f"rate_{k}"].to_numpy(), in1)})
print("Transients per second (reward periods left out):")
print(pd.DataFrame(rows).round(4).to_string(index=False))
Transients per second (reward periods left out):
transients n_in n_out median_in median_out p p_exact
DA_pos 28 146 0.2261 0.1760 0.0000 NaN
DA_neg 28 146 0.2866 0.2301 0.0000 NaN
ACh_pos 28 146 0.2392 0.2068 0.0013 NaN
ACh_neg 28 146 0.2691 0.2249 0.0000 NaN
Result: reproduced. All four transient types are more frequent inside the hotspot: 0.226 against 0.176 per second for DA peaks, 0.287 against 0.230 for DA troughs, 0.239 against 0.207 for ACh peaks and 0.269 against 0.225 for ACh troughs (medians per second of usable data; p of 0.0013 or less).
fig3b-reward-correlation¶
After unpredicted rewards, trial-by-trial DA peak magnitude correlates with ACh; dominant correlations are negative at 114 of 174 sites with DA leading (177.9 +/- 96.7 ms) and positive at 58 of 174 with ACh leading (-41.2 +/- 218.9 ms).
IDX3 = np.arange(-9, 28)
er3 = daach.event_rasters(reward, lambda s: daach.first_licks(s)[1], -9, 27)
res3 = {m: daach.event_correlation(d["rasters"]["DA"], d["rasters"]["ACh"], 1, IDX3, np.arange(0, 19)) for m, d in er3.items()}
dom3 = daach.stack(res3, lambda v: v["dominant"])
lat3 = daach.stack(res3, lambda v: v["lag"]) * 1000 / 18
mice3, coords3 = daach.site_table(reward)
assert np.allclose(coords3, coords1)
print("Rewards with a consumption lick, per mouse:", {m: len(d["events"]) for m, d in er3.items()})
def corr_summary(dom, lat):
neg, pos = dom < 0, dom > 0
return {"sites": int((~np.isnan(dom)).sum()), "negative": int(neg.sum()), "negative_lag>0": int((neg & (lat > 0)).sum()),
"mean_latency_negative_ms": lat[neg].mean(), "sd_negative": lat[neg].std(ddof=1),
"positive": int(pos.sum()), "positive_lag<0": int((pos & (lat < 0)).sum()),
"mean_latency_positive_ms": lat[pos].mean(), "sd_positive": lat[pos].std(ddof=1)}
s3 = corr_summary(dom3, lat3)
print({k: round(v, 1) if isinstance(v, float) else v for k, v in s3.items()})
checks = [("negative", 114, 11.4), ("positive", 58, 5.8), ("mean_latency_negative_ms", 177.9, 9.1),
("mean_latency_positive_ms", -41.2, 28.7)]
print("Within tolerance:", {k: abs(s3[k] - e) <= t for k, e, t in checks})
print(pd.DataFrame([{"mouse": m, **{k: (round(v, 1) if isinstance(v, float) else v) for k, v in
corr_summary(dom3[mice3 == m], lat3[mice3 == m]).items()}} for m in np.unique(mice3)]).to_string(index=False))
Rewards with a consumption lick, per mouse: {'UG27': 176, 'UG28': 273, 'UG29': 355, 'UG30': 213, 'UG31': 288}
{'sites': 174, 'negative': 122, 'negative_lag>0': 120, 'mean_latency_negative_ms': np.float64(172.1), 'sd_negative': np.float64(78.0), 'positive': 52, 'positive_lag<0': 37, 'mean_latency_positive_ms': np.float64(-59.8), 'sd_positive': np.float64(166.4)}
Within tolerance: {'negative': True, 'positive': False, 'mean_latency_negative_ms': np.True_, 'mean_latency_positive_ms': np.True_}
mouse sites negative negative_lag>0 mean_latency_negative_ms sd_negative positive positive_lag<0 mean_latency_positive_ms sd_positive
UG27 23 21 21 169.3 56.9 2 2 -138.9 117.9
UG28 41 24 22 178.2 75.0 17 12 -0.0 259.1
UG29 48 27 27 174.9 73.4 21 13 -63.5 84.8
UG30 28 22 22 143.9 70.0 6 5 -148.1 75.9
UG31 34 28 28 188.5 100.7 6 5 -101.9 64.9
Result: partially reproduced. The dominant reward correlation is negative at 122 sites with a mean latency of 172.1 ms (the paper 114 and 177.9 ms) and positive at 52 sites with a mean latency of -59.8 ms (the paper 58 and -41.2 ms). Three of the four values are within tolerance; the positive count is 0.2 below its lower bound (52.2). DA leads at 120 of the 122 negative sites and ACh at 37 of the 52 positive ones. Negative sites outnumber positive ones in each mouse.
fig3d-reward-in-out¶
Reward-evoked DA-ACh correlations are significantly more negative inside the global hotspot than outside.
print(daach.in_out(dom3, in1))
print(pd.DataFrame([{"mouse": m, **daach.in_out(dom3[mice3 == m], in1[mice3 == m])} for m in np.unique(mice3)]).round(4).to_string(index=False))
{'n_in': 28, 'n_out': 146, 'median_in': np.float64(-0.5152765765454144), 'median_out': np.float64(-0.30573910454878656), 'p': np.float64(7.250515759083041e-09), 'p_exact': nan}
mouse n_in n_out median_in median_out p p_exact
UG27 3 20 -0.5164 -0.4886 0.3858 0.4043
UG28 8 33 -0.5543 0.1128 0.0004 0.0001
UG29 2 46 -0.4942 -0.1682 0.0597 0.0532
UG30 5 23 -0.6274 -0.4349 0.0164 0.0128
UG31 10 24 -0.4464 -0.3192 0.0006 0.0003
Result: reproduced. Inside the hotspot the median reward r is -0.52 against -0.31 outside (p = 7.3e-09). The difference has the same sign in every mouse and is significant in three (UG28, UG30, UG31).
fig3e-reward-hotspot-overlap¶
Reward anti-correlations form an aDLS hotspot that significantly overlaps the hotspot of spontaneous correlations (volume overlap test, p = 0.008).
rng = np.random.default_rng(301)
hs3 = daach.hotspot(interp1, dom3, rng)
ov3 = daach.overlap_test(hs1, interp1.inside, hs3, interp1.inside, GRID.shape, rng)
print(f"Reward hotspot: components {hs3['n_components']}, mean {hs3['mean']:.3f}, voxels {ov3['n_b']}, "
f"centroid {hs3['centroid'].round(2)}")
print("Overlap with the spontaneous hotspot:", ov3)
show_map(hs3, "Dominant reward DA-ACh r; black, reward hotspot", vmax=0.6)
Reward hotspot: components 2, mean -0.526, voxels 1022, centroid [ 0.54 2.07 -3.26]
Overlap with the spontaneous hotspot: {'overlap': 721, 'p': 0.01, 'n_a': 1010, 'n_b': 1022, 'null_median': 116.0}
Smoothed values range -0.649 to 0.335; color scale clipped at +/-0.6; voxels beyond it: 59 of 6155
Result: reproduced. The reward correlations form a negative hotspot (mean -0.53, 1022 voxels, centroid AP 0.54, ML 2.07, DV -3.26 mm) that shares 721 voxels with the spontaneous hotspot, against a null median of 116 (overlap p = 0.01; the paper 0.008).
fig4-cue-sequence¶
After Pavlovian conditioning, the conditioned light cue evokes an ACh peak, then a DA peak, then an ACh dip at a majority of sites (49 of 75).
W4 = (-9, 54)
IDX4 = np.arange(W4[0], W4[1] + 1)
def light_events(s):
return daach.event_index(s, s["events"]["Light"][0])
def iti_exclusion(s):
cues = light_events(s)
rew = daach.event_index(s, s["events"]["Reward"][0])
return [(i, i + int(9 * daach.SR)) for i in cues] + [(i, i + int(5 * daach.SR)) for i in rew]
learn = [s for s in pav if s["phase"] == "learning"]
ext = [s for s in pav if s["phase"] == "extinction"]
rng = np.random.default_rng(401)
er4 = daach.event_rasters(learn, light_events, *W4, exclude_of=iti_exclusion)
sig4 = {w: daach.triggered_significance(er4, *W4, w, rng) for w in ((9, 27), (9, 36))}
mice4, coords4 = daach.site_table(learn)
def sequence(sig, rows, alpha="methods"):
out = []
for m in sorted(sig):
a, d = sig[m]["ACh"], sig[m]["DA"]
ap, an = a[alpha]
dp, _ = d[alpha]
w = slice(rows[0], rows[1] + 1)
ta, td, tn = np.argmax(a["mean"][w], 0), np.argmax(d["mean"][w], 0), np.argmin(a["mean"][w], 0)
out.append(ap & dp & an & (ta < td) & (td < tn))
return np.concatenate(out)
for label, rows in (("(a) 0 to 1 s", (9, 27)), ("(b) 0 to 1.5 s", (9, 36))):
seq = sequence(sig4[rows], rows)
seq_code = sequence(sig4[rows], rows, "code")
print(f"Reading {label}: sequence at {int(seq.sum())} of {len(seq)} sites; by mouse "
f"{pd.Series(seq).groupby(mice4).sum().to_dict()}; code percentiles {int(seq_code.sum())}")
s = sig4[(9, 27)]
print("Significant features in 0 to 1 s (Methods percentiles):",
{k: int(daach.stack(s, f).sum()) for k, f in (("ACh peak", lambda v: v["ACh"]["methods"][0]),
("ACh dip", lambda v: v["ACh"]["methods"][1]),
("DA peak", lambda v: v["DA"]["methods"][0]))})
print("Light-cue trials per mouse:", {m: len(d["events"]) for m, d in er4.items()})
Reading (a) 0 to 1 s: sequence at 58 of 75 sites; by mouse {'AD1': 12, 'AD2': 19, 'AD3': 27}; code percentiles 55
Reading (b) 0 to 1.5 s: sequence at 58 of 75 sites; by mouse {'AD1': 12, 'AD2': 19, 'AD3': 27}; code percentiles 52
Significant features in 0 to 1 s (Methods percentiles): {'ACh peak': 68, 'ACh dip': 72, 'DA peak': 75}
Light-cue trials per mouse: {'AD1': 90, 'AD2': 90, 'AD3': 90}
Result: not reproduced. The light cue evokes an ACh peak, then a DA peak, then an ACh dip at 58 of 75 sites under both readings, more than the paper's 49 and above the tolerance (44.1 to 53.9). The statement's qualitative part holds: the sequence is at a majority of sites overall and in each mouse (12 of 17, 19 of 24, 27 of 34). With the code's percentiles the count is 55 or 52.
fig4c-cue-correlation¶
For conditioned light cues, ACh correlates with the trial-by-trial DA peak negatively at lagging times at 39 of 75 sites (202.28 +/- 128.4 ms) and positively at leading times at 26 of 75 (-79.1 +/- 224.0 ms).
res4 = {m: daach.event_correlation(d["rasters"]["DA"], d["rasters"]["ACh"], 1, IDX4, np.arange(0, 19)) for m, d in er4.items()}
dom4 = daach.stack(res4, lambda v: v["dominant"])
lat4 = daach.stack(res4, lambda v: v["lag"]) * 1000 / 18
sigDA4 = daach.stack(sig4[(9, 27)], lambda v: v["DA"]["methods"][0])
checks = [("negative", 39, 3.9), ("positive", 26, 2.6), ("mean_latency_negative_ms", 202.28, 20.6),
("mean_latency_positive_ms", -79.1, 43.9)]
for label, sel in (("(a) sites with a significant DA peak", sigDA4), ("(b) all sites", np.ones_like(sigDA4))):
sm_ = corr_summary(dom4[sel], lat4[sel])
print(label, {k: round(v, 1) if isinstance(v, float) else v for k, v in sm_.items()})
print(" within tolerance:", {k: abs(sm_[k] - e) <= t for k, e, t in checks})
print(pd.DataFrame([{"mouse": m, **{k: (round(v, 1) if isinstance(v, float) else v) for k, v in
corr_summary(dom4[mice4 == m], lat4[mice4 == m]).items()}} for m in np.unique(mice4)]).to_string(index=False))
print("Random-time candidates for the cue null, per mouse:", {m: len(d["sampleable"]) for m, d in er4.items()},
"; light-cue trials:", {m: len(d["events"]) for m, d in er4.items()})
sigDA4c = daach.stack(sig4[(9, 27)], lambda v: v["DA"]["code"][0])
sm_ = corr_summary(dom4[sigDA4c], lat4[sigDA4c])
print(f"Sites with a significant cue DA peak under the code's percentiles: {int(sigDA4c.sum())};",
{k: round(v, 1) if isinstance(v, float) else v for k, v in sm_.items()})
(a) sites with a significant DA peak {'sites': 75, 'negative': 41, 'negative_lag>0': 37, 'mean_latency_negative_ms': np.float64(206.0), 'sd_negative': np.float64(174.5), 'positive': 34, 'positive_lag<0': 22, 'mean_latency_positive_ms': np.float64(-83.3), 'sd_positive': np.float64(311.1)}
within tolerance: {'negative': True, 'positive': False, 'mean_latency_negative_ms': np.True_, 'mean_latency_positive_ms': np.True_}
(b) all sites {'sites': 75, 'negative': 41, 'negative_lag>0': 37, 'mean_latency_negative_ms': np.float64(206.0), 'sd_negative': np.float64(174.5), 'positive': 34, 'positive_lag<0': 22, 'mean_latency_positive_ms': np.float64(-83.3), 'sd_positive': np.float64(311.1)}
within tolerance: {'negative': True, 'positive': False, 'mean_latency_negative_ms': np.True_, 'mean_latency_positive_ms': np.True_}
mouse sites negative negative_lag>0 mean_latency_negative_ms sd_negative positive positive_lag<0 mean_latency_positive_ms sd_positive
AD1 17 12 11 175.9 113.2 5 4 -55.6 333.3
AD2 24 9 8 216.0 189.3 15 11 -225.9 276.6
AD3 34 20 18 219.4 201.9 14 7 59.5 288.6
Random-time candidates for the cue null, per mouse: {'AD1': 68868, 'AD2': 69375, 'AD3': 74196} ; light-cue trials: {'AD1': 90, 'AD2': 90, 'AD3': 90}
Sites with a significant cue DA peak under the code's percentiles: 75; {'sites': 75, 'negative': 41, 'negative_lag>0': 37, 'mean_latency_negative_ms': np.float64(206.0), 'sd_negative': np.float64(174.5), 'positive': 34, 'positive_lag<0': 22, 'mean_latency_positive_ms': np.float64(-83.3), 'sd_positive': np.float64(311.1)}
Result: partially reproduced. All 75 sites have a significant cue DA peak, so both readings give the same values: 41 negative sites with a mean latency of 206.0 ms (the paper 39 and 202.28 ms) and 34 positive sites with a mean latency of -83.3 ms (the paper 26 and -79.1 ms). Three values are within tolerance; the positive count is above its range (23.4 to 28.6). ACh lags at 37 of the 41 negative sites and leads at 22 of the 34 positive ones. The code's percentiles also find a cue DA peak at all 75 sites; the null is drawn from about 69000 to 74000 candidate times per mouse, so its size does not explain the paper's 65.
fig4h-extinction-correlation¶
After cue-specific extinction, cue-evoked DA dips correlate negatively with lagging ACh at 41 of 74 sites, mean lag 117.9 +/- 176.2 ms.
rng = np.random.default_rng(402)
er4x = daach.event_rasters(ext, light_events, *W4, exclude_of=iti_exclusion)
sig4x = daach.triggered_significance(er4x, *W4, (9, 36), rng)
mice4x, coords4x = daach.site_table(ext)
assert np.allclose(coords4x, coords4)
dip4x = daach.stack(sig4x, lambda v: v["DA"]["methods"][1])
print("Light-cue trials per mouse:", {m: len(d["events"]) for m, d in er4x.items()})
print("Sites with a significant DA dip in 0 to 1.5 s:", int(dip4x.sum()), "of", len(dip4x),
"; by mouse", pd.Series(dip4x).groupby(mice4x).sum().to_dict(),
"; code percentiles:", int(daach.stack(sig4x, lambda v: v["DA"]["code"][1]).sum()))
for pol, label in ((-1, "DA trough (Methods, main)"), (1, "DA peak (code's polarity, variant)")):
rx = {m: daach.event_correlation(d["rasters"]["DA"], d["rasters"]["ACh"], pol, IDX4, np.arange(0, 28)) for m, d in er4x.items()}
d4x = daach.stack(rx, lambda v: v["dominant"])
l4x = daach.stack(rx, lambda v: v["lag"]) * 1000 / 18
if pol == -1:
dom4x, lat4x = d4x, l4x
for rl, sel in (("(a) significant-dip sites", dip4x), ("(b) all sites", np.ones_like(dip4x))):
neg = sel & (d4x < 0)
print(f"{label}, {rl}: {int(sel.sum())} sites, negative {int(neg.sum())} (lag > 0: {int((neg & (l4x > 0)).sum())}), "
f"mean lag {l4x[neg].mean():.1f} ms (sd {l4x[neg].std(ddof=1):.1f}); within tolerance: count "
f"{abs(neg.sum() - 41) <= 4.1}, lag {abs(l4x[neg].mean() - 117.9) <= 27.5}")
print(pd.DataFrame([{"mouse": m, "dip_sites": int(dip4x[mice4x == m].sum()),
"negative": int((dip4x & (dom4x < 0))[mice4x == m].sum()),
"mean_lag_neg": lat4x[(mice4x == m) & dip4x & (dom4x < 0)].mean()} for m in np.unique(mice4x)]).round(1).to_string(index=False))
Light-cue trials per mouse: {'AD1': 90, 'AD2': 90, 'AD3': 90}
Sites with a significant DA dip in 0 to 1.5 s: 61 of 75 ; by mouse {'AD1': 10, 'AD2': 19, 'AD3': 32} ; code percentiles: 42
DA trough (Methods, main), (a) significant-dip sites: 61 sites, negative 38 (lag > 0: 28), mean lag 118.4 ms (sd 266.0); within tolerance: count True, lag True DA trough (Methods, main), (b) all sites: 75 sites, negative 47 (lag > 0: 34), mean lag 114.7 ms (sd 251.2); within tolerance: count False, lag True
DA peak (code's polarity, variant), (a) significant-dip sites: 61 sites, negative 44 (lag > 0: 38), mean lag 137.6 ms (sd 152.3); within tolerance: count True, lag True DA peak (code's polarity, variant), (b) all sites: 75 sites, negative 50 (lag > 0: 44), mean lag 140.0 ms (sd 144.7); within tolerance: count False, lag True mouse dip_sites negative mean_lag_neg AD1 10 9 234.6 AD2 19 7 -79.4 AD3 32 22 133.8
Result: reproduced. Under reading (a), 38 of the 61 sites with a significant cue DA dip have a dominant negative correlation with a mean lag of 118.4 ms (the paper 41 sites and 117.9 ms); both are within tolerance. Over all 75 sites (reading b) 47 are negative with a mean lag of 114.7 ms. ACh lags at 28 of the 38 negative sites, and in AD2 the mean lag of negative sites is -79.4 ms. The code's polarity (the DA peak) gives 44 negative sites with a mean lag of 137.6 ms.
fig4d-cue-in-out¶
Conditioned-cue DA-ACh correlations are more negative at sites inside the global aDLS hotspot than outside.
in4 = daach.inside(GRID, hs1, coords4)
print("AD1-AD3 sites inside the global hotspot:", int(in4.sum()), "by mouse", pd.Series(in4).groupby(mice4).sum().to_dict())
print("Sites with a significant DA peak:", daach.in_out(np.where(sigDA4, dom4, np.nan), in4))
print("Variant, all sites:", daach.in_out(dom4, in4))
interp4 = daach.Interpolator(GRID, coords4)
rng = np.random.default_rng(403)
hs4 = daach.hotspot(interp4, dom4, rng)
print(f"Variant, cue hotspot: mean {hs4['mean']:.3f}, centroid {hs4['centroid'].round(2)}, overlap",
daach.overlap_test(hs1, interp1.inside, hs4, interp4.inside, GRID.shape, rng))
AD1-AD3 sites inside the global hotspot: 9 by mouse {'AD1': 2, 'AD2': 2, 'AD3': 5}
Sites with a significant DA peak: {'n_in': 9, 'n_out': 66, 'median_in': np.float64(-0.30992225798198275), 'median_out': np.float64(0.01905363306720756), 'p': np.float64(0.02394033082268308), 'p_exact': np.float64(0.022122117824525993)}
Variant, all sites: {'n_in': 9, 'n_out': 66, 'median_in': np.float64(-0.30992225798198275), 'median_out': np.float64(0.01905363306720756), 'p': np.float64(0.02394033082268308), 'p_exact': np.float64(0.022122117824525993)}
Variant, cue hotspot: mean -0.298, centroid [ 0.65 1.66 -3.41], overlap {'overlap': 147, 'p': 0.175, 'n_a': 1010, 'n_b': 471, 'null_median': 3.5}
Result: reproduced. The 9 AD1 to AD3 sites inside the global hotspot have a median cue r of -0.31 against 0.02 for the 66 outside (p = 0.024, exact p = 0.022). A hotspot drawn from the cue correlations themselves is negative, but its overlap with the global hotspot is not significant (147 voxels, p = 0.175).
fig4j-extinction-in-out¶
After extinction, DA dip to ACh correlations are more negative at sites inside the global aDLS hotspot than outside.
print("Sites with a significant DA dip:", daach.in_out(np.where(dip4x, dom4x, np.nan), in4))
print("Variant, all sites:", daach.in_out(dom4x, in4))
rng = np.random.default_rng(404)
hs4x = daach.hotspot(daach.Interpolator(GRID, coords4[dip4x]), dom4x[dip4x], rng)
print(f"Variant, extinction hotspot over dip sites: mean {hs4x['mean']:.3f}, centroid {hs4x['centroid'].round(2)}, overlap",
daach.overlap_test(hs1, interp1.inside, hs4x, daach.Interpolator(GRID, coords4[dip4x]).inside, GRID.shape, rng))
Sites with a significant DA dip: {'n_in': 8, 'n_out': 53, 'median_in': np.float64(-0.3494790084589413), 'median_out': np.float64(-0.21004865793496025), 'p': np.float64(0.03533521336710137), 'p_exact': np.float64(0.03350226630317037)}
Variant, all sites: {'n_in': 9, 'n_out': 66, 'median_in': np.float64(-0.3522801629626127), 'median_out': np.float64(-0.20537897012266856), 'p': np.float64(0.018477459683827645), 'p_exact': np.float64(0.01663639241412914)}
Variant, extinction hotspot over dip sites: mean -0.300, centroid [ 0.81 1.61 -3.66], overlap {'overlap': 128, 'p': 0.17, 'n_a': 1010, 'n_b': 356, 'null_median': 9.0}
Result: reproduced. Among sites with a significant DA dip, the 8 inside the global hotspot have a median r of -0.35 against -0.21 for the 53 outside (p = 0.035, exact p = 0.034); over all sites p = 0.018. The extinction hotspot's own overlap with the global hotspot is not significant (128 voxels, p = 0.17).
fig5d-initiation-responses¶
At locomotion initiations, ACh rises and then falls at 131 of 174 sites; DA dips at 119 of 174 sites, rises at 73, and both rises and dips at 32.
W5 = (-18, 36)
IDX5 = np.arange(W5[0], W5[1] + 1)
SIG5 = (9, 36) # rows for -0.5 to 1 s
def init_rasters(loco):
out = {}
for m, d in loco.items():
samp = daach.sampleable_times(len(d["vel"]), *W5, valid=~np.isnan(d["vel"]))
samp = np.setdiff1d(samp, d["init"])
out[m] = {"cat": d["cat"], "events": d["init"], "sampleable": samp,
"rasters": {k: daach.raster(d["cat"][k], d["init"], *W5) for k in ("DA", "ACh")}}
return out
def init_counts(sig, alpha="methods"):
dp = daach.stack(sig, lambda v: v["DA"][alpha][0]); dn = daach.stack(sig, lambda v: v["DA"][alpha][1])
ap = daach.stack(sig, lambda v: v["ACh"][alpha][0]); an = daach.stack(sig, lambda v: v["ACh"][alpha][1])
ta = daach.stack(sig, lambda v: np.argmax(v["ACh"]["mean"][SIG5[0]:SIG5[1] + 1], 0))
tn = daach.stack(sig, lambda v: np.argmin(v["ACh"]["mean"][SIG5[0]:SIG5[1] + 1], 0))
return {"ACh rise then fall (a)": int((ap & an & (ta < tn)).sum()), "ACh rise (b)": int(ap.sum()),
"DA dip": int(dn.sum()), "DA rise": int(dp.sum()), "DA both": int((dp & dn).sum()),
"ACh dip": int(an.sum()), "ACh both, any order": int((ap & an).sum())}
loco = daach.locomotion_events(spont)
print("Initiations (bout >= 2.27 s) per mouse:", {m: len(d["init"]) for m, d in loco.items()},
"; all clean onsets:", {m: len(d["init_all"]) for m, d in loco.items()})
er5 = init_rasters(loco)
sig5 = daach.triggered_significance(er5, *W5, SIG5, np.random.default_rng(501))
c5 = init_counts(sig5)
print("Counts (Methods percentiles):", c5)
exp5 = {"DA dip": 119, "DA rise": 73, "DA both": 32}
print("Within 10%:", {k: abs(c5[k] - v) <= 0.1 * v for k, v in exp5.items()},
"; ACh (a)", abs(c5["ACh rise then fall (a)"] - 131) <= 13.1, "; ACh (b)", abs(c5["ACh rise (b)"] - 131) <= 13.1)
print("Variant, code percentiles (0.5/99.5):", init_counts(sig5, "code"))
mice5 = np.concatenate([[m] * d["cat"]["DA"].shape[1] for m, d in loco.items()])
per5 = pd.DataFrame([{"mouse": m, **init_counts({m: sig5[m]})} for m in sig5])
print(per5.to_string(index=False))
loco22 = daach.locomotion_events(spont, bout_s=2.2)
sig22 = daach.triggered_significance(init_rasters(loco22), *W5, SIG5, np.random.default_rng(502))
print("Variant, the code's bout of 2.2 s:", {m: len(d["init"]) for m, d in loco22.items()}, init_counts(sig22))
Initiations (bout >= 2.27 s) per mouse: {'UG27': 501, 'UG28': 92, 'UG29': 135, 'UG30': 125, 'UG31': 149} ; all clean onsets: {'UG27': 857, 'UG28': 351, 'UG29': 476, 'UG30': 280, 'UG31': 685}
Counts (Methods percentiles): {'ACh rise then fall (a)': 107, 'ACh rise (b)': 147, 'DA dip': 129, 'DA rise': 77, 'DA both': 46, 'ACh dip': 148, 'ACh both, any order': 135}
Within 10%: {'DA dip': True, 'DA rise': True, 'DA both': False} ; ACh (a) False ; ACh (b) False
Variant, code percentiles (0.5/99.5): {'ACh rise then fall (a)': 97, 'ACh rise (b)': 136, 'DA dip': 114, 'DA rise': 61, 'DA both': 26, 'ACh dip': 129, 'ACh both, any order': 117}
mouse ACh rise then fall (a) ACh rise (b) DA dip DA rise DA both ACh dip ACh both, any order
UG27 19 23 20 6 3 23 23
UG28 30 34 29 16 8 34 31
UG29 12 36 36 14 8 34 29
UG30 19 25 20 10 5 25 24
UG31 27 29 24 31 22 32 28
Variant, the code's bout of 2.2 s: {'UG27': 512, 'UG28': 93, 'UG29': 140, 'UG30': 128, 'UG31': 157} {'ACh rise then fall (a)': 106, 'ACh rise (b)': 147, 'DA dip': 132, 'DA rise': 74, 'DA both': 44, 'ACh dip': 149, 'ACh both, any order': 136}
fig, axs = plt.subplots(1, 2, figsize=(9, 3.2), sharex=True)
t = IDX5 / 18
for ax, ch in zip(axs, ("DA", "ACh")):
mean = np.concatenate([sig5[m][ch]["mean"] for m in sig5], 1)
dn = daach.stack(sig5, lambda v: v[ch]["methods"][1])
dp = daach.stack(sig5, lambda v: v[ch]["methods"][0])
ax.plot(t, mean[:, dn & ~dp].mean(1), label=f"dip only (n {int((dn & ~dp).sum())})")
ax.plot(t, mean[:, dp & ~dn].mean(1), label=f"rise only (n {int((dp & ~dn).sum())})")
ax.plot(t, mean[:, dp & dn].mean(1), label=f"both (n {int((dp & dn).sum())})")
ax.axvline(0, color="k", lw=0.5)
ax.set_title(f"{ch} at initiation"); ax.set_xlabel("s"); ax.legend(fontsize=7)
print(ch, "plotted means range", np.nanmin(mean.mean(1)).round(4), "to", np.nanmax(mean.mean(1)).round(4))
axs[0].set_ylabel("dF/F")
plt.show()
DA plotted means range -0.0031 to 0.0022 ACh plotted means range -0.0032 to 0.0048
Result: partially reproduced. At initiations, 129 sites have a significant DA dip and 77 a DA rise (the paper 119 and 73, both within 10%), but 46 have both (the paper 32). The ACh count is 107 under reading (a), a rise followed by a dip, and 147 under reading (b), a rise; neither is within 117.9 to 144.1. The code's percentiles give 114 dips, 61 rises and 26 sites with both. UG27 contributes 501 initiations, 3.4 to 5.4 times as many as each other mouse (92 to 149).
fig5-latency-order¶
Initiation-related signals occur in sequence; DA dips first (median 111.11 ms), then ACh peaks (222.22 ms) and DA peaks (227.78 ms), then ACh dips (611.11 ms).
def latencies(sig, alpha="methods"):
out = {}
for ch in ("DA", "ACh"):
mean = np.concatenate([sig[m][ch]["mean"] for m in sorted(sig)], 1)[SIG5[0]:SIG5[1] + 1]
t = IDX5[SIG5[0]:SIG5[1] + 1] * 1000 / 18
pk = daach.stack(sig, lambda v: v[ch][alpha][0]); dp = daach.stack(sig, lambda v: v[ch][alpha][1])
out[f"{ch} peak"] = t[np.argmax(mean, 0)][pk]
out[f"{ch} dip"] = t[np.argmin(mean, 0)][dp]
return out
lat5 = latencies(sig5)
exp_l = {"DA dip": 111.11, "ACh peak": 222.22, "DA peak": 227.78, "ACh dip": 611.11}
rows = [{"feature": k, "sites": len(v), "median_ms": np.median(v), "IQR": np.percentile(v, [25, 75]).round(1),
"paper": exp_l[k], "within_55.6ms": abs(np.median(v) - exp_l[k]) <= 55.6} for k, v in lat5.items()]
print(pd.DataFrame(rows).to_string(index=False))
m_ = {k: np.median(v) for k, v in lat5.items()}
print("Order DA dip < ACh peak, DA peak < ACh dip:", m_["DA dip"] < min(m_["ACh peak"], m_["DA peak"]) and
max(m_["ACh peak"], m_["DA peak"]) < m_["ACh dip"])
print("By mouse (median ms):")
print(pd.DataFrame({m: {k: np.median(v) if len(v) else np.nan for k, v in latencies({m: sig5[m]}).items()} for m in sig5}).round(1).to_string())
print("Variant, code percentiles:", {k: round(float(np.median(v)), 1) for k, v in latencies(sig5, "code").items()})
feature sites median_ms IQR paper within_55.6ms
DA peak 77 277.777778 [-444.4, 500.0] 227.78 True
DA dip 129 111.111111 [55.6, 777.8] 111.11 True
ACh peak 147 277.777778 [222.2, 277.8] 222.22 True
ACh dip 148 666.666667 [555.6, 722.2] 611.11 True
Order DA dip < ACh peak, DA peak < ACh dip: True
By mouse (median ms):
UG27 UG28 UG29 UG30 UG31
DA peak 527.8 500.0 388.9 -222.2 -388.9
DA dip 111.1 55.6 111.1 111.1 1000.0
ACh peak 277.8 277.8 277.8 222.2 333.3
ACh dip 722.2 722.2 -111.1 666.7 722.2
Variant, code percentiles: {'DA peak': 277.8, 'DA dip': 138.9, 'ACh peak': 277.8, 'ACh dip': 722.2}
Result: reproduced. Median latencies are 111.1 ms for DA dips, 277.8 ms for ACh peaks and for DA peaks, and 666.7 ms for ACh dips: the paper's order, each within one sample of its value (111.11, 222.22, 227.78 and 611.11 ms). The ACh peak and ACh dip medians are exactly one sample later than the paper's, at the edge of the tolerance. Per mouse the order is less regular: the DA peak median ranges from -388.9 to 527.8 ms, and UG31's DA dip median is 1000 ms.
fig5e-initiation-topography¶
Initiation-related DA dips are concentrated in the aDLS and DA peaks in dorsomedial and posterior striatum; ACh dips are localized to the aDLS.
mean_da = np.concatenate([sig5[m]["DA"]["mean"] for m in sorted(sig5)], 1)[SIG5[0]:SIG5[1] + 1]
mx, mn = mean_da.max(0), mean_da.min(0)
dom_da = np.where(np.abs(mn) > np.abs(mx), mn, mx)
ach_dip = daach.stack(sig5, lambda v: v["ACh"]["methods"][1])
coords5 = coords1
rows = []
for label, a, b in (("negative vs positive dominant DA", dom_da < 0, dom_da > 0), ("ACh dip vs none", ach_dip, ~ach_dip)):
for d, name in ((0, "AP"), (1, "ML"), (2, "DV")):
p = stats.mannwhitneyu(coords5[a, d], coords5[b, d], alternative="greater").pvalue
rows.append({"comparison": label, "axis": name, "n": f"{int(a.sum())} vs {int(b.sum())}",
"median_a": np.median(coords5[a, d]), "median_b": np.median(coords5[b, d]), "one_sided_p_greater": p})
print(pd.DataFrame(rows).round(4).to_string(index=False))
print("Inside the spontaneous hotspot: DA negative", int((in1 & (dom_da < 0)).sum()), "of", int(in1.sum()),
"; ACh dip", int((in1 & ach_dip).sum()), "of", int(in1.sum()))
comparison axis n median_a median_b one_sided_p_greater
negative vs positive dominant DA AP 106 vs 68 0.355 -0.475 0.0000
negative vs positive dominant DA ML 106 vs 68 1.880 2.300 0.9962
negative vs positive dominant DA DV 106 vs 68 -3.480 -2.850 1.0000
ACh dip vs none AP 148 vs 26 0.100 0.055 0.3349
ACh dip vs none ML 148 vs 26 2.045 1.765 0.2470
ACh dip vs none DV 148 vs 26 -3.190 -3.295 0.0887
Inside the spontaneous hotspot: DA negative 28 of 28 ; ACh dip 28 of 28
Result: not reproduced. Sites whose dominant DA response is a dip are more anterior than sites whose dominant response is a rise (median AP 0.36 against -0.48 mm, p < 0.0001), but more medial (ML 1.88 against 2.30 mm; p = 0.996 for lateral) and more ventral, so part (i) fails on ML. Sites with a significant ACh dip (148 of 174) are not more anterior or lateral than the 26 without (p = 0.33 and 0.25). All 28 sites inside the correlation hotspot have a dominant DA dip and a significant ACh dip.
fig6a-initiation-correlation¶
At sites with initiation-related DA dips (119 sites), trial-by-trial DA dip magnitudes correlate negatively with ACh with DA leading (mean latency 193.80 +/- 177.68 ms).
res6 = {m: daach.event_correlation(d["rasters"]["DA"], d["rasters"]["ACh"], -1, IDX5, np.arange(-9, 19)) for m, d in er5.items()}
dom6 = daach.stack(res6, lambda v: v["dominant"])
lat6 = daach.stack(res6, lambda v: v["lag"]) * 1000 / 18
dip5 = daach.stack(sig5, lambda v: v["DA"]["methods"][1])
neg = dip5 & (dom6 < 0)
print(f"Significant-dip sites {int(dip5.sum())}; dominant negative {int(neg.sum())} ({neg.sum() / dip5.sum():.2f}); "
f"negative with lag > 0: {int((neg & (lat6 > 0)).sum())}")
print(f"(a) mean latency over negative sites {lat6[neg].mean():.1f} ms (sd {lat6[neg].std(ddof=1):.1f}); "
f"(b) over all dip sites {lat6[dip5].mean():.1f} ms (sd {lat6[dip5].std(ddof=1):.1f}); paper 193.80 +/- 177.68, tolerance 16.3")
print(pd.DataFrame([{"mouse": m, "dip_sites": int(dip5[mice5 == m].sum()), "negative": int(neg[mice5 == m].sum()),
"mean_latency_neg": lat6[neg & (mice5 == m)].mean()} for m in np.unique(mice5)]).round(1).to_string(index=False))
# Synthetic check of the trough search next to NaNs: a trough one sample after a NaN must be found.
t_rng = np.random.default_rng(0)
L_ = np.arange(-9, 28)
ref_ = 0.01 * t_rng.normal(size=(len(L_), 1, 40))
amp_ = t_rng.uniform(0.5, 2, 40)
k_ = int(np.flatnonzero(L_ == 5)[0])
ref_[k_, 0] -= amp_
ref_[k_ - 1, 0, :10] = np.nan
mag_ = daach.event_correlation(ref_, ref_, -1, L_, np.arange(0, 19))["magnitude"][0]
assert np.allclose(mag_, ref_[k_, 0]), "trough next to a NaN missed"
print("Synthetic check passed: troughs next to a NaN are found (40 of 40).")
n_nan = sum(int(np.isnan(d["rasters"]["DA"][(IDX5 >= -18) & (IDX5 <= 27)]).any(0).sum()) for d in er5.values())
n_all = sum(d["rasters"]["DA"].shape[1] * d["rasters"]["DA"].shape[2] for d in er5.values())
print(f"Site-trials with a NaN in the DA search range (-1 to 1.5 s): {n_nan} of {n_all}")
Significant-dip sites 129; dominant negative 98 (0.76); negative with lag > 0: 88 (a) mean latency over negative sites 188.8 ms (sd 157.3); (b) over all dip sites 98.2 ms (sd 230.9); paper 193.80 +/- 177.68, tolerance 16.3 mouse dip_sites negative mean_latency_neg UG27 20 18 145.1 UG28 29 25 186.7 UG29 36 24 219.9 UG30 20 17 232.0 UG31 24 14 142.9 Synthetic check passed: troughs next to a NaN are found (40 of 40). Site-trials with a NaN in the DA search range (-1 to 1.5 s): 2197 of 30341
Result: reproduced. Of the 129 sites with a significant DA dip, 98 (76%) have a dominant negative correlation, 88 of them with ACh lagging, with a mean latency of 188.8 ms (reading a; the paper 193.80 ms, tolerance 177.5 to 210.1). Over all dip sites (reading b) the mean is 98.2 ms. Most dip sites are negative in every mouse. The trough search now ignores NaNs; before this fix, a trough within 2 samples of a NaN was missed. 2197 of 30341 site-trials have a NaN in the range searched, so at most that many were affected. Before the fix the result was 99 sites with a mean latency of 182.4 ms.
fig6c-initiation-in-out¶
Initiation DA dip to ACh correlations are concentrated in the aDLS and more negative inside the global hotspot than outside.
print(daach.in_out(np.where(dip5, dom6, np.nan), in1))
print("Variant, all sites:", daach.in_out(dom6, in1))
{'n_in': 28, 'n_out': 101, 'median_in': np.float64(-0.3516358833405038), 'median_out': np.float64(-0.21052444299903356), 'p': np.float64(1.7143028453334125e-05), 'p_exact': nan}
Variant, all sites: {'n_in': 28, 'n_out': 146, 'median_in': np.float64(-0.3516358833405038), 'median_out': np.float64(-0.1972676297686128), 'p': np.float64(1.7390734663708561e-06), 'p_exact': nan}
Result: reproduced. Among dip sites the initiation correlation is more negative inside the hotspot (median -0.35, 28 sites) than outside (-0.21, 101 sites), p = 1.7e-05.
fig6g-invigoration-in-out¶
At locomotion invigorations, GLM-estimated DA peak magnitudes correlate negatively with ACh at lagging times, and these correlations are more negative inside the global hotspot (131 sites with invigoration DA peaks).
t0 = time.time()
inv = daach.invigoration_analysis(loco)
print(f"GLM fits in {time.time() - t0:.0f} s; invigorations per mouse:", {m: v["n_invig"] for m, v in inv.items()})
sig_inv = daach.stack(inv, lambda v: v["sig_da_peak"])
res6g = {m: daach.event_correlation(v["rasters"]["DA"], v["rasters"]["ACh"], 1, IDX5, np.arange(-9, 19)) for m, v in inv.items()}
dom6g = daach.stack(res6g, lambda v: v["dominant"])
lat6g = daach.stack(res6g, lambda v: v["lag"]) * 1000 / 18
neg = sig_inv & (dom6g < 0)
print(f"Sites with an invigoration DA peak: {int(sig_inv.sum())} (paper 131); dominant negative {int(neg.sum())}, "
f"lag > 0 among them {int((neg & (lat6g > 0)).sum())}, mean lag {lat6g[neg].mean():.1f} ms")
print("In vs out among selected sites:", daach.in_out(np.where(sig_inv, dom6g, np.nan), in1))
print("Variant, all sites:", daach.in_out(dom6g, in1))
print(pd.DataFrame([{"mouse": m, "selected": int(sig_inv[mice5 == m].sum()), "negative": int(neg[mice5 == m].sum())}
for m in np.unique(mice5)]).to_string(index=False))
GLM fits in 50 s; invigorations per mouse: {'UG27': 501, 'UG28': 92, 'UG29': 135, 'UG30': 125, 'UG31': 149}
Sites with an invigoration DA peak: 160 (paper 131); dominant negative 105, lag > 0 among them 97, mean lag 179.4 ms
In vs out among selected sites: {'n_in': 28, 'n_out': 132, 'median_in': np.float64(-0.33570808109480876), 'median_out': np.float64(-0.20747085269871357), 'p': np.float64(4.507854180118807e-05), 'p_exact': nan}
Variant, all sites: {'n_in': 28, 'n_out': 146, 'median_in': np.float64(-0.33570808109480876), 'median_out': np.float64(-0.19706490015416175), 'p': np.float64(2.3505223167175328e-05), 'p_exact': nan}
mouse selected negative
UG27 23 17
UG28 38 26
UG29 40 21
UG30 27 18
UG31 32 23
Result: reproduced. The GLM finds an invigoration DA peak at 160 sites (the paper 131). Of these, 105 have a dominant negative correlation, 97 of them with ACh lagging (mean lag 179.4 ms), and the correlations are more negative inside the hotspot (median -0.34, 28 sites) than outside (-0.21, 132 sites), p = 4.5e-05.
fig7-da-release¶
Optogenetic stimulation of midbrain DA neurons evokes significant DA release at 179 of 180 sites.
W7 = (-9, 45)
SIG7 = (7, 45) # rows for -2 to 36 samples
def iti_stims(s):
st = daach.event_index(s, s["events"]["Optogenetic Stimulation"][0])
rew = daach.event_index(s, s["events"]["Reward"][0])
return np.array([i for i in st if not (np.abs(rew - i) <= 18).any()], int)
def stim_exclusion(s):
rew, lick = daach.first_licks(s)
st = daach.event_index(s, s["events"]["Optogenetic Stimulation"][0])
return ([(i, i + int(5 * daach.SR) + 1) for i in np.concatenate([rew, lick])] + [(i, i + 19) for i in st])
er7 = daach.event_rasters(ads, iti_stims, *W7, exclude_of=stim_exclusion)
sig7 = daach.triggered_significance(er7, *W7, SIG7, np.random.default_rng(701))
print("ITI stimulation trains per mouse:", {m: len(d["events"]) for m, d in er7.items()},
"; all trains:", {m: int(sum(len(s["events"]["Optogenetic Stimulation"][0]) for s in v)) for m, v in daach.by_mouse(ads).items()})
rows = []
for m in sig7:
rows.append({"mouse": m, "sites": len(sig7[m]["DA"]["methods"][0]), "DA_peak": int(sig7[m]["DA"]["methods"][0].sum()),
"DA_peak_code": int(sig7[m]["DA"]["code"][0].sum()), "ACh_dip": int(sig7[m]["ACh"]["methods"][1].sum()),
"ACh_dip_code": int(sig7[m]["ACh"]["code"][1].sum())})
t7 = pd.DataFrame(rows)
print(t7.to_string(index=False))
print("Totals:", t7.drop(columns="mouse").sum().to_dict())
ITI stimulation trains per mouse: {'ADS12': 20, 'ADS13': 20, 'ADS16': 20, 'ADS17': 20, 'ADS19': 20, 'ADS6': 15} ; all trains: {'ADS12': 52, 'ADS13': 52, 'ADS16': 52, 'ADS17': 52, 'ADS19': 52, 'ADS6': 30}
mouse sites DA_peak DA_peak_code ACh_dip ACh_dip_code
ADS12 13 13 13 10 9
ADS13 23 22 21 21 19
ADS16 37 37 37 28 22
ADS17 27 27 27 27 24
ADS19 27 27 27 10 6
ADS6 53 53 53 51 48
Totals: {'sites': 180, 'DA_peak': 179, 'DA_peak_code': 178, 'ACh_dip': 147, 'ACh_dip_code': 128}
Result: reproduced. Stimulation trains in the ITI (20 per mouse, 15 for ADS6) evoke a significant DA peak at 179 of 180 sites (178 with the code's percentiles), as in the paper.
fig7e-ach-suppression¶
DA stimulation evokes significant ACh suppression at a subset of sites (116 of 180, 64%).
n_dip = int(t7.ACh_dip.sum())
print(f"Sites with a significant ACh dip: {n_dip} of {int(t7.sites.sum())} ({n_dip / t7.sites.sum():.2f}); "
f"code percentiles {int(t7.ACh_dip_code.sum())}; paper 116 (0.64), tolerance 11.6")
fig, ax = plt.subplots(figsize=(4.5, 3))
for ch, c in (("DA", "tab:red"), ("ACh", "tab:green")):
mean = np.concatenate([sig7[m][ch]["mean"] for m in sorted(sig7)], 1)
ax.plot(np.arange(W7[0], W7[1] + 1) / 18, mean.mean(1), color=c, label=ch)
print(ch, "mean over sites, range", mean.mean(1).min().round(4), "to", mean.mean(1).max().round(4))
ax.axvline(0, color="k", lw=0.5); ax.set_xlabel("s from stimulation onset"); ax.set_ylabel("dF/F"); ax.legend()
plt.show()
Sites with a significant ACh dip: 147 of 180 (0.82); code percentiles 128; paper 116 (0.64), tolerance 11.6 DA mean over sites, range -0.0006 to 0.0251 ACh mean over sites, range -0.0051 to 0.0012
Result: partially reproduced. A significant ACh dip follows stimulation at 147 of 180 sites (82%), against 116 (64%) in the paper; with the code's percentiles 128 (71%). The 33 sites without a dip are concentrated in ADS19 (17 of its 27 sites) and ADS16 (9 of 37).
ed9-cohort-hotspot¶
The stimulation cohort replicates the aDLS anti-correlation hotspot in spontaneous release outside stimulation periods.
def ads_mask(mode):
def f(s):
n = len(s["DA"])
rew, lick = daach.first_licks(s)
st = daach.event_index(s, s["events"]["Optogenetic Stimulation"][0])
a, b, c = daach.time_since(n, rew) < 6, daach.time_since(n, lick) < 6, daach.time_since(n, st) < 3
return (a & b & c) if mode == "and" else (a | b | c)
return f
frac = {mode: sum(int(ads_mask(mode)(s).sum()) for s in ads) / sum(len(s["DA"]) for s in ads) for mode in ("and", "or")}
print(f"Fraction of samples left out: AND {frac['and']:.4f}, OR {frac['or']:.4f}")
def post_stim(s):
return daach.time_since(len(s["DA"]), daach.event_index(s, s["events"]["Optogenetic Stimulation"][0])) < 3
kept = {mode: sum(int((post_stim(s) & ~ads_mask(mode)(s)).sum()) for s in ads) / sum(int(post_stim(s).sum()) for s in ads)
for mode in ("and", "or")}
print(f"Fraction of samples within 3 s after stimulation that are kept: AND {kept['and']:.4f}, OR {kept['or']:.4f}")
xc9 = daach.spontaneous_xcorr(ads, mask=ads_mask("and"), concatenate=True)
dom9, lag9 = xc9["dominant"], xc9["lag"]
print("Cohort counts (code's AND exclusion):", polarity_counts(dom9, lag9))
in9 = daach.inside(GRID, hs1, coords_s)
print("Sites inside the original hotspot:", int(in9.sum()), "by mouse", pd.Series(in9).groupby(mice_s).sum().to_dict())
print("(i) in vs out of the original hotspot:", daach.in_out(dom9, in9))
interp9 = daach.Interpolator(GRID, coords_s)
rng = np.random.default_rng(901)
hs9 = daach.hotspot(interp9, dom9, rng)
ov9 = daach.overlap_test(hs1, interp1.inside, hs9, interp9.inside, GRID.shape, rng)
print(f"(ii) cohort hotspot: mean {hs9['mean']:.3f}, voxels {ov9['n_b']}, centroid {hs9['centroid'].round(2)}; overlap", ov9)
xc9o = daach.spontaneous_xcorr(ads, mask=ads_mask("or"), concatenate=True)
print("Variant, OR exclusion:", polarity_counts(xc9o["dominant"], xc9o["lag"]), daach.in_out(xc9o["dominant"], in9))
hs9o = daach.hotspot(interp9, xc9o["dominant"], rng)
print(f"Variant, OR exclusion hotspot: mean {hs9o['mean']:.3f}, centroid {hs9o['centroid'].round(2)}; overlap",
daach.overlap_test(hs1, interp1.inside, hs9o, interp9.inside, GRID.shape, rng))
Fraction of samples left out: AND 0.0649, OR 0.3270 Fraction of samples within 3 s after stimulation that are kept: AND 0.4244, OR 0.0000
Cohort counts (code's AND exclusion): {'negative': 167, 'negative, DA leads': 163, 'positive': 13, 'positive, ACh leads': 10, 'lag 0': 4, 'sites': 180}
Sites inside the original hotspot: 9 by mouse {'ADS12': 0, 'ADS13': 3, 'ADS16': 2, 'ADS17': 1, 'ADS19': 3, 'ADS6': 0}
(i) in vs out of the original hotspot: {'n_in': 9, 'n_out': 171, 'median_in': np.float64(-0.41868141769856504), 'median_out': np.float64(-0.3011621770529894), 'p': np.float64(0.048949152994268996), 'p_exact': np.float64(0.0477655485703631)}
(ii) cohort hotspot: mean -0.439, voxels 1536, centroid [ 0.71 2. -3.45]; overlap {'overlap': 564, 'p': 0.085, 'n_a': 1010, 'n_b': 1536, 'null_median': 65.0}
Variant, OR exclusion: {'negative': 157, 'negative, DA leads': 148, 'positive': 23, 'positive, ACh leads': 19, 'lag 0': 5, 'sites': 180} {'n_in': 9, 'n_out': 171, 'median_in': np.float64(-0.4002437409948406), 'median_out': np.float64(-0.2498742858498871), 'p': np.float64(0.04122775063991309), 'p_exact': np.float64(0.039808044380070136)}
Variant, OR exclusion hotspot: mean -0.414, centroid [ 0.72 2.08 -3.43]; overlap {'overlap': 551, 'p': 0.085, 'n_a': 1010, 'n_b': 1432, 'null_median': 81.0}
Result: partially reproduced. In the stimulation cohort, the 9 sites inside the original hotspot have more negative dominant correlations (median -0.42) than the 171 outside (-0.30), p = 0.049 (exact p = 0.048); the OR exclusion gives p = 0.041. The cohort's own hotspot is negative (mean -0.44, 1536 voxels, centroid AP 0.71, ML 2.00, DV -3.45 mm) and shares 564 voxels with the original against a null median of 65, but the overlap test gives p = 0.085. 167 of the 180 sites have a dominant negative correlation. The code's AND exclusion leaves out 6.5% of samples and the OR exclusion 32.7%. Of the samples within 3 s after a stimulation, AND keeps 42% and OR none; with OR the cohort hotspot is similar (mean -0.41, centroid AP 0.72, ML 2.08, DV -3.43 mm) and its overlap p is also 0.085.
fig7f-ach-dip-in-out¶
Stimulation-evoked ACh dips are larger inside the cohort's correlation hotspot than outside, and their own hotspot overlaps it.
amin = np.concatenate([np.nanmean(np.nanmin(er7[m]["rasters"]["ACh"][SIG7[0]:SIG7[1] + 1], 0), 1) for m in sorted(er7)])
in9c = daach.inside(GRID, hs9, coords_s)
print("Sites inside the cohort hotspot:", int(in9c.sum()))
print("ACh minima, in vs out of the cohort hotspot:", daach.in_out(amin, in9c))
print("Variant, in vs out of the original hotspot:", daach.in_out(amin, in9))
rng = np.random.default_rng(902)
hs7 = daach.hotspot(interp9, amin, rng, how="min")
print(f"Variant, ACh-minimum hotspot: mean {hs7['mean']:.4f}, centroid {hs7['centroid'].round(2)}; overlap with cohort hotspot",
daach.overlap_test(hs9, interp9.inside, hs7, interp9.inside, GRID.shape, rng))
print(pd.DataFrame([{"mouse": m, **daach.in_out(amin[mice_s == m], in9c[mice_s == m])} for m in np.unique(mice_s)]).round(4).to_string(index=False))
in9o = daach.inside(GRID, hs9o, coords_s)
print("Variant, in vs out of the OR-exclusion cohort hotspot:", daach.in_out(amin, in9o))
Sites inside the cohort hotspot: 17
ACh minima, in vs out of the cohort hotspot: {'n_in': 17, 'n_out': 163, 'median_in': np.float32(-0.017874172), 'median_out': np.float32(-0.014852333), 'p': np.float32(0.00054375955), 'p_exact': nan}
Variant, in vs out of the original hotspot: {'n_in': 9, 'n_out': 171, 'median_in': np.float32(-0.017507488), 'median_out': np.float32(-0.015073419), 'p': np.float32(0.09419366), 'p_exact': np.float32(0.09460347)}
Variant, ACh-minimum hotspot: mean -0.0199, centroid [ 0.78 2.03 -3.4 ]; overlap with cohort hotspot {'overlap': 755, 'p': 0.065, 'n_a': 1536, 'n_b': 1357, 'null_median': 45.0}
mouse n_in n_out median_in median_out p p_exact
ADS12 1 12 -0.0203 -0.0146 0.4615 0.4615
ADS13 1 22 -0.0216 -0.0158 0.4975 0.6087
ADS16 3 34 -0.0201 -0.0118 0.0797 0.0808
ADS17 4 23 -0.0168 -0.0137 0.2323 0.2430
ADS19 4 23 -0.0171 -0.0109 0.0105 0.0060
ADS6 4 49 -0.0192 -0.0160 0.1174 0.1202
Variant, in vs out of the OR-exclusion cohort hotspot: {'n_in': 19, 'n_out': 161, 'median_in': np.float32(-0.018816296), 'median_out': np.float32(-0.014852333), 'p': np.float32(0.00026237094), 'p_exact': nan}
Result: reproduced. ACh minima after stimulation are deeper at the 17 sites inside the cohort hotspot (median -0.0179 dF/F) than at the 163 outside (-0.0149), p = 0.0005. The AND exclusion keeps 42% of the samples within 3 s after a stimulation; with the hotspot of the OR exclusion, which leaves them all out, the 19 sites inside also have deeper minima (median -0.0188 against -0.0149, p = 0.0003). Against the original hotspot (9 sites inside) the difference has the same sign but p = 0.094, and the ACh-minimum hotspot's overlap with the cohort hotspot gives p = 0.065, so the claim's second part (the dips' own hotspot overlaps the cohort hotspot) is not supported at p < 0.05; the criterion treats it as a variant. Within each of the six mice the inside sites have lower minima, significantly so in one (ADS19).
fig8b-quinpirole-region¶
Quinpirole suppresses evoked ACh release in DLS and NAcC, more strongly in DLS; region F(1,40) = 105.9, interaction F(4,40) = 3.014, stimulus F(4,40) = 17.80.
The slice data are not in the dandiset. The paper's source data for Figure 8 are on Zenodo (doi:10.5281/zenodo.18202873) and are read from there.
FIG8 = "https://zenodo.org/records/18202873/files/Bouabid_Vu%20%282026%29_Source_Data_Fig8.xlsx?download=1"
fig8 = pd.read_excel(FIG8, sheet_name=None)
print("Sheets:", list(fig8))
def anova(sheet, factor):
d = fig8[sheet].rename(columns={factor: "g"})
# Type III sums of squares with sum-to-zero contrasts, as GraphPad Prism (named in the Methods)
# computes an ordinary two-way ANOVA; the groups are unbalanced in some sheets
fit = smf.ols("Peak_ACh ~ C(g, Sum) * C(Stim, Sum)", data=d).fit()
tab = sm.stats.anova_lm(fit, typ=3)
t2 = sm.stats.anova_lm(smf.ols("Peak_ACh ~ C(g) * C(Stim)", data=d).fit(), typ=2)
print(sheet, "| rows", len(d), "| slices per group", d.groupby("g").Slice_ID.nunique().to_dict(),
"| mouse ids per group", d.groupby("g").Mouse_ID.nunique().to_dict(), "| group means", d.groupby("g").Peak_ACh.mean().round(1).to_dict())
print(tab.round(4).to_string())
print("Variant, Type II F:", t2.F.dropna().round(3).tolist())
return tab
t8b = anova("Figure_8b", "Brain_Region")
print(fig8["Figure_8b"].groupby(["Brain_Region", "Slice_ID"]).Stim.apply(list).to_string())
Sheets: ['Figure_8a_LP', 'Figure_8a_RP', 'Figure_8b', 'Figure_8c_TP', 'Figure_8c_LP', 'Figure_8c_RP', 'Figure_8d_LP', 'Figure_8d_RP', 'Figure_8e', 'Figure_8f', 'Figure_8g']
Figure_8b | rows 50 | slices per group {'DLS': 6, 'NACc': 5} | mouse ids per group {'DLS': 4, 'NACc': 4} | group means {'DLS': 25.1, 'NACc': 68.1}
sum_sq df F PR(>F)
Intercept 108430.7103 1.0 496.4379 0.000
C(g, Sum) 23123.1232 1.0 105.8666 0.000
C(Stim, Sum) 15550.2608 4.0 17.7988 0.000
C(g, Sum):C(Stim, Sum) 2633.2482 4.0 3.0140 0.029
Residual 8736.6984 40.0 NaN NaN
Variant, Type II F: [105.867, 17.799, 3.014]
Brain_Region Slice_ID
DLS slice_A_4 [1, 2, 3, 4, 5]
slice_B_2 [1, 2, 3, 4, 5]
slice_B_3 [1, 2, 3, 4]
slice_B_4 [5]
slice_C_1 [1, 2, 3, 4, 5]
slice_D_3 [1, 2, 3, 4, 5]
NACc slice_E_3 [1, 2, 3, 4, 5]
slice_E_4 [1, 2, 3, 4, 5]
slice_F_1 [1, 2, 3, 4, 5]
slice_G_1 [1, 2, 3, 4, 5]
slice_H_1 [1, 2, 3, 4, 5]
Result: reproduced. The two-way ANOVA on the source data (Type III sums of squares) gives region F = 105.87, stimulus F = 17.80 and interaction F = 3.014 (the paper 105.9, 17.80 and 3.014), with DLS lower (25.1 against 68.1); Type II gives the same values. The DLS group holds 6 slice labels from 4 mouse labels, one of them (slice_B_4) with only stimulus 5, where the legend states 5 slices from 5 mice.
fig8e-cocaine-region¶
Cocaine reduces evoked ACh release more in DLS than in NAcC; interaction F(4,55) = 10.75, stimulus F(4,55) = 36.37, region F(1,55) = 199.9.
t8e = anova("Figure_8e", "Brain_Region")
Figure_8e | rows 65 | slices per group {'DLS': 8, 'NACc': 5} | mouse ids per group {'DLS': 5, 'NACc': 4} | group means {'DLS': 32.8, 'NACc': 74.1}
sum_sq df F PR(>F)
Intercept 175779.0436 1.0 1338.5560 0.0
C(g, Sum) 26245.1661 1.0 199.8567 0.0
C(Stim, Sum) 19106.3355 4.0 36.3736 0.0
C(g, Sum):C(Stim, Sum) 5649.1102 4.0 10.7545 0.0
Residual 7222.5944 55.0 NaN NaN
Variant, Type II F: [199.857, 47.676, 10.754]
Result: reproduced. With Type III sums of squares, region F = 199.86, stimulus F = 36.37 and interaction F = 10.75 match the paper (199.9, 36.37 and 10.75), with DLS lower (32.8 against 74.1). The design is unbalanced (8 DLS and 5 NAcC slices); Type II sums of squares give a stimulus F of 47.68.
fig8fg-d2-block¶
The D2 antagonist attenuates the cocaine-induced suppression in both regions; block effect F(1,55) = 231.9 in DLS and F(1,40) = 32.57 in NAcC.
t8f = anova("Figure_8f", "Condition")
t8g = anova("Figure_8g", "Condition")
Figure_8f | rows 65 | slices per group {'Cocaine': 8, 'Cocaine + L-741.626': 5} | mouse ids per group {'Cocaine': 5, 'Cocaine + L-741.626': 5} | group means {'Cocaine': 32.8, 'Cocaine + L-741.626': 78.6}
sum_sq df F PR(>F)
Intercept 190946.1590 1.0 1370.8162 0.0
C(g, Sum) 32298.3258 1.0 231.8720 0.0
C(Stim, Sum) 19200.6840 4.0 34.4608 0.0
C(g, Sum):C(Stim, Sum) 5784.8634 4.0 10.3825 0.0
Residual 7661.1576 55.0 NaN NaN
Variant, Type II F: [231.872, 45.039, 10.382]
Figure_8g | rows 50 | slices per group {'Cocaine': 5, 'Cocaine + L-741.626': 5} | mouse ids per group {'Cocaine': 4, 'Cocaine + L-741.626': 5} | group means {'Cocaine': 74.1, 'Cocaine + L-741.626': 107.9}
sum_sq df F PR(>F)
Intercept 413811.4457 1.0 946.4767 0.0000
C(g, Sum) 14241.4979 1.0 32.5734 0.0000
C(Stim, Sum) 2648.9180 4.0 1.5147 0.2163
C(g, Sum):C(Stim, Sum) 1522.7768 4.0 0.8707 0.4899
Residual 17488.5003 40.0 NaN NaN
Variant, Type II F: [32.573, 1.515, 0.871]
Result: reproduced. The D2R block effect is F = 231.87 in DLS and 32.57 in NAcC (the paper 231.9 and 32.57), and the antagonist raises ACh in both (DLS 32.8 to 78.6, NAcC 74.1 to 107.9). With Type III sums of squares the DLS stimulus F is 34.46, as in the legend (Type II gives 45.04).
fig8c-cocaine-da-similar (blocked)¶
Cocaine elevates evoked DA release to a similar degree in DLS and NAcC (cocaine by region interaction F(8,36) = 0.20).
print("Figure_8c_TP columns:", fig8["Figure_8c_TP"].columns.tolist(), "| rows", len(fig8["Figure_8c_TP"]))
for sh in ("Figure_8c_LP", "Figure_8c_RP"):
d = fig8[sh]
print(sh, "| mice", d.Mouse_ID.unique().tolist(), "| conditions", d.Condition.unique().tolist(),
"| time points", d.Timepoints.nunique(), d.Timepoints.min(), "to", d.Timepoints.max())
pk = d.groupby(["Mouse_ID", "Condition"]).Mean_DA.max().unstack()
print(" peak DA per mouse:", pk.round(3).to_dict("index"), "| cocaine / control:", (pk.Cocaine / pk.Control).round(2).tolist(),
"| mean ratio", round((pk.Cocaine / pk.Control).mean(), 2))
Figure_8c_TP columns: ['Voltage', 'Current'] | rows 403
Figure_8c_LP | mice ['mouse_A', 'mouse_B', 'mouse_C'] | conditions ['Control', 'Cocaine'] | time points 24 0.625 to 3.5
peak DA per mouse: {'mouse_A': {'Cocaine': 1.354, 'Control': 0.995}, 'mouse_B': {'Cocaine': 1.302, 'Control': 0.992}, 'mouse_C': {'Cocaine': 1.772, 'Control': 1.001}} | cocaine / control: [1.36, 1.31, 1.77] | mean ratio 1.48
Figure_8c_RP | mice ['mouse_D', 'mouse_E', 'mouse_F'] | conditions ['Control', 'Cocaine'] | time points 24 0.625 to 3.5
peak DA per mouse: {'mouse_D': {'Cocaine': 2.047, 'Control': 0.977}, 'mouse_E': {'Cocaine': 1.13, 'Control': 1.03}, 'mouse_F': {'Cocaine': 1.239, 'Control': 0.996}} | cocaine / control: [2.09, 1.1, 1.24] | mean ratio 1.48
The sheet Figure_8c_TP holds a voltammogram (voltage and current); Figure_8c_LP and _RP hold mean DA time courses (24 time points) per mouse and condition, three mice per region. The legend's ANOVA has 8 degrees of freedom for cocaine and 36 residual, which fits neither a two-level cocaine factor nor 24 time points, so the test cannot be repeated. Descriptively, the peak DA with cocaine is 1.48 times the control peak on average in both regions (DLS 1.31 to 1.77, NAcC 1.10 to 2.09 per mouse).
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome |
|---|---|
| fig1e-dominant-polarity | reproduced |
| fig1f-adls-hotspot | partially reproduced (centroid 0.02 mm medial of the median site) |
| fig1h-in-vs-out | reproduced |
| fig1f-positive-posterior | reproduced |
| ed2-hotspot-per-mouse | reproduced |
| ed1-controls-no-organization | not reproduced (control correlations as strong as sensor ones; a hotspot is found) |
| fig2a-directionality | reproduced |
| fig2-magnitude-correlations | partially reproduced |
| fig2d-magnitude-in-out | partially reproduced |
| fig2b-pairing-overlaps-hotspot | reproduced |
| ed3e-transient-rates-in-hotspot | reproduced |
| fig3b-reward-correlation | partially reproduced |
| fig3d-reward-in-out | reproduced |
| fig3e-reward-hotspot-overlap | reproduced |
| fig4-cue-sequence | not reproduced (58 of 75 sites, the paper 49) |
| fig4c-cue-correlation | partially reproduced |
| fig4h-extinction-correlation | reproduced |
| fig4d-cue-in-out | reproduced |
| fig4j-extinction-in-out | reproduced |
| fig5d-initiation-responses | partially reproduced |
| fig5-latency-order | reproduced |
| fig5e-initiation-topography | not reproduced |
| fig6a-initiation-correlation | reproduced |
| fig6c-initiation-in-out | reproduced |
| fig6g-invigoration-in-out | reproduced |
| fig7-da-release | reproduced |
| fig7e-ach-suppression | partially reproduced |
| ed9-cohort-hotspot | partially reproduced |
| fig7f-ach-dip-in-out | reproduced |
| fig8b-quinpirole-region (Zenodo data) | reproduced |
| fig8e-cocaine-region (Zenodo data) | reproduced |
| fig8fg-d2-block (Zenodo data) | reproduced |
| fig8c-cocaine-da-similar (Zenodo data) | blocked, methods underspecified |
How the analysis differs from the paper's. The hotspot procedure is a simplified port of the
paper's MATLAB pipeline: linear in place of natural-neighbor interpolation, 100 um in place of 50 um
voxels, the convex hull of the sites in place of the Allen striatum mask, and 200 in place of 10000
shuffles and random volumes. The chance occurrence of paired transients uses 10 in place of 10000
shuffles of sessions. Bootstrap nulls use 1000 samples in place of 10000. MATLAB's lowpass
is replaced by a zero-phase Butterworth filter. Pavlovian phases are inferred from the files'
events. The control analysis uses the 15 uncued reward sessions of AD4 to AD6. Pairs of transients
are matched by the Methods' rule (the first partner), a simplification of the code's matching.
The Figure 8 ANOVAs use Type III sums of squares, as Prism does.
What was found about the dandiset. The photometry series state a rate of 9 Hz; the timestamps, the Methods and the code give 18 Hz. The files do not say which experiment, cohort or Pavlovian phase a session belongs to; this was inferred from subject names, indicators and events. Tone cue events are absent from the Pavlovian files. One reward session of UG27 (240313r) lacks the artifact-removed green series. The site counts match the paper (174, 75 and 180), and the data support every spontaneous, reward and locomotion analysis tested. The excitation source descriptions differ from the Methods (405 nm isosbestic light at 100 Hz, where the Methods give 470 and 570 nm at 36 Hz), the mutant ACh sensor is named ACh4.3mut where the paper names ACh3.0-mut, and the session description's reward volume and intervals differ from the Methods. The paper's Figure 8 source data are on Zenodo, not in the dandiset.
print(f"Notebook ran in {(time.time() - started) / 60:.1f} min")
Notebook ran in 14.2 min