Reproducing Markowitz, Gillis, Jay et al. (2023) from DANDI:000559¶
Paper: Markowitz J.E., Gillis W.F., Jay M. et al., "Spontaneous behaviour is structured by reinforcement without explicit reward", Nature 614, 108-117 (2023), doi:10.1038/s41586-022-05611-2. The full text was used (open access, CC BY 4.0).
Dandiset: DANDI:000559, version 0.260528.0858. The paper cites no dandiset version (its Data availability names a Zenodo record). The archive has two published versions, 0.240502.0456 and 0.260528.0858; they hold the same 4,425 assets (identical asset ids, blobs and sizes, checked below) and differ only in dandiset-level metadata (anatomy terms, study targets, ORCIDs, related resources). The most recent one is pinned.
Every file is streamed from the archive; nothing is downloaded to disk. Claims and criteria are in
claims/000559/10.1038-s41586-022-05611-2.yaml.
Code read before the criteria were written: the paper's analysis code,
dattalab/dopamine-reinforces-spontaneous-behavior
at commit 14ed714 (notebooks for every panel and the rl_analysis package; it reads the paper's
Zenodo parquet files, not NWB), and the conversion code,
catalystneuro/datta-lab-to-nwb at commit
47c681e, which shows how the NWB series were made from those parquet files. Where the Methods
leave a parameter open, the analysis follows the paper's code; spontda.py says which.
| Claim | Tested here |
|---|---|
| fig1e-velocity-timescales | yes |
| fig1h-syllable-dlight-differs | yes |
| fig1l-identity-not-decodable | yes |
| fig2b-dlight-usage-across-syllables | yes |
| fig2i-dlight-entropy-across-syllables | yes |
| fig2e-dlight-predicts-future-usage | yes |
| ed7-dms-fewer-transients | yes |
| ed7-dms-no-usage-prediction | yes |
| fig3e-optoda-learning | yes |
| fig3f-cumulative-counts | yes |
| fig3e-learner-count | yes |
| fig3h-stim-post-correlation | yes |
| fig3l-no-kinematic-effect | run; blocked: metadata insufficient (catch trials not marked; the substitute catch set of the criterion fails in no-opsin controls) |
| fig3m-long-stim-velocity | run; blocked: metadata insufficient (as for fig3l) |
| fig2a-food-reward-transients | blocked: data missing (no food-reward sessions or consumption times) |
| fig1f-syllable-vs-movement-onset | blocked: not attempted here |
| fig2k-dlight-predicts-sequence-variability | blocked: not attempted here |
| fig3i-optoda-entropy | blocked: metadata insufficient (catch trials not marked); not attempted |
| ed10d-optoda-vigour | blocked: not attempted here |
| fig4b-endo-da-influence-consistent | blocked: not attempted here |
| fig4c-endo-da-predicts-learning | blocked: not attempted here |
| fig4f-catch-dlight-predicts-learning | blocked: metadata insufficient (catch trials not marked); not attempted |
| fig4h-decoding-model | blocked: not attempted here |
| fig5-rl-model | blocked: not attempted here |
In [1]:
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import spontda as sd
sd.set_threads(1)
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
T0 = time.time()
print("dandiset", sd.DANDISET_ID, "version", sd.DANDISET_VERSION)
dandiset 000559 version 0.260528.0858
Loading¶
The asset list of the pinned version, then every photometry session (DLS and DMS dLight mice, including the closed-loop sessions of the dLight+ChrimsonR mice) and every closed-loop session of the Opto-DA and control mice. Only the series the analyses use are read: syllable labels and their frame clock, dLight and reference dF/F, 2D speed and stimulation times.
In [2]:
assets = sd.asset_table()
nwb = assets[assets.is_nwb].copy()
# Version check: the earlier published version holds the same assets.
from dandi.dandiapi import DandiAPIClient
with DandiAPIClient() as client:
old = {a.path: (a.identifier, a.size) for a in client.get_dandiset(sd.DANDISET_ID, "0.240502.0456").get_assets()}
new = {p: (i, s) for p, i, s in zip(assets.path, assets.asset_id, assets["size"])}
print("assets in 0.240502.0456:", len(old), " in", sd.DANDISET_VERSION + ":", len(new),
" identical path, asset id and size:", sum(old.get(p) == v for p, v in new.items()))
phot_assets = nwb[nwb.subject.str.contains("dlight") & nwb.kind.isin(["photometry", "reinforcement"])
& nwb.path.str.contains("ophys")]
cl_assets = nwb[(nwb.kind == "reinforcement") & ~nwb.subject.str.startswith("dlight-chrimson")]
# closed-loop files of dLight+ChrimsonR mice without a dLight series (read only for the data check)
chr_extra_assets = nwb[(nwb.kind == "reinforcement") & nwb.subject.str.startswith("dlight-chrimson")
& ~nwb.path.str.contains("ophys")]
t = time.time()
phot = sd.parallel_map(sd.load_photometry_session, phot_assets.asset_id, 12)
print(f"photometry files: {len(phot)} in {time.time() - t:.0f} s")
t = time.time()
cl = sd.parallel_map(sd.load_closed_loop_with_speed, cl_assets.asset_id, 12)
print(f"closed-loop files: {len(cl)} in {time.time() - t:.0f} s")
chr_extra = sd.parallel_map(sd.load_closed_loop_session, chr_extra_assets.asset_id, 12)
assets in 0.240502.0456: 4425 in 0.260528.0858: 4425 identical path, asset id and size: 4425
photometry files: 383 in 129 s
closed-loop files: 1340 in 281 s
Data check¶
In [3]:
print("NWB files by session kind and number of subjects:")
print(nwb.groupby("kind").agg(files=("path", "size"), subjects=("subject", "nunique")).to_string())
print("other assets:", (~assets.is_nwb).sum(), "external video files (",
assets.loc[~assets.is_nwb, "path"].str.extract(r"\.(\w+)$")[0].value_counts().to_dict(), ")")
print("total size of NWB files: %.0f GB, of all assets: %.1f TB" % (nwb["size"].sum() / 1e9, assets["size"].sum() / 1e12))
# --- Photometry sessions -------------------------------------------------------------------
rows = []
for p, d in zip(phot_assets.path, phot):
has_syl = d["t"] is not None and len(d["t"]) > 0
row = dict(path=p, mouse=d["subject"], sex=d["sex"], description=d["description"],
day=d["start"].tz_localize(None).normalize(), has_syl=has_syl, has_ophys=d["has_ophys"],
n_pulses=len(d.get("pulse_on", [])), targets=tuple(d.get("targets", [])),
rotated=d.get("clock_rotated", False), backward_steps=d.get("clock_backward_steps", 0),
jump_s=d.get("clock_jump_s", 0.0), repeated=d.get("clock_repeated", 0),
dt_unique=d.get("clock_dt_unique", np.nan), habituation="habit" in d["description"].lower())
if has_syl and d["has_ophys"]:
row.update(sd.quality(d["signal_dff"], d["reference_dff"]))
row["passes"] = sd.passes_quality(row)
row["ophys_on_clock"], row["speed_on_clock"] = d["ophys_on_clock"], d["speed_on_clock"]
row["dt_median"] = float(np.median(np.diff(d["t"])))
row["duration_s"] = d["t"][-1] - d["t"][0]
rows.append(row)
P = pd.DataFrame(rows)
P["area"] = np.where(P.mouse.str.startswith("dms"), "dms", "dls")
P["stim"] = P.n_pulses > 0
P["chrimson"] = P.mouse.str.startswith("dlight-chrimson")
# post sessions: sessions without stimulation on the day after a stimulation day of the same mouse
stim_days = set(zip(P.loc[P.stim, "mouse"], P.loc[P.stim, "day"]))
P["post"] = [(not s) and ((m, d - pd.Timedelta(days=1)) in stim_days) for m, d, s in zip(P.mouse, P.day, P.stim)]
P["passes"] = P["passes"].fillna(False).astype(bool)
# sessions used by the photometry analyses: syllables and dLight, quality passed, no stimulation, not habituation
P["use"] = P.has_syl & P.passes & ~P.stim & ~P.habituation
print()
print("Photometry files (dLight in file) by area and stimulation:")
print(P.groupby(["area", "stim"]).agg(files=("path", "size"), mice=("mouse", "nunique"),
with_syllables=("has_syl", "sum"), pass_quality=("passes", "sum")).to_string())
print("files whose syllable series are empty:", (~P.has_syl).sum(), "; their targets:",
P.loc[~P.has_syl, "targets"].value_counts().to_dict())
print("habituation sessions among photometry files:", P.groupby("area").habituation.sum().to_dict(),
"; mice:", sorted(P.loc[P.habituation, "mouse"].unique()))
use_P = P[P.use]
print()
print("Sessions without stimulation that pass the quality criteria, per mouse:")
print(use_P.groupby(["area", "mouse"]).agg(sessions=("path", "size"), post=("post", "sum")).to_string())
print("DLS mice:", use_P[use_P.area == "dls"].mouse.nunique(), " DMS mice:", use_P[use_P.area == "dms"].mouse.nunique())
print("DLS sessions:", (use_P.area == "dls").sum(), " of which not post:", ((use_P.area == "dls") & ~use_P.post).sum(),
"; DMS sessions:", (use_P.area == "dms").sum())
print("photometry-kind files without a dLight series:",
int((nwb.kind == "photometry").sum() - (phot_assets.path.str.contains("ses-photometry")).sum()),
"; mice whose photometry files never hold dLight:",
sorted(set(nwb.loc[nwb.kind == "photometry", "subject"]) - set(P.mouse)))
# Files left after each filter, per area and per mouse (the Methods give 379 DLS experiments; 378 in Fig. 2)
P["f_no_stim"] = ~P.stim
P["f_syllables"] = P.f_no_stim & P.has_syl
P["f_not_habituation"] = P.f_syllables & ~P.habituation
P["f_quality"] = P.f_not_habituation & P.passes
flt = P.groupby(["area", "mouse"]).agg(files=("path", "size"), no_stim=("f_no_stim", "sum"), with_syllables=("f_syllables", "sum"),
not_habituation=("f_not_habituation", "sum"), pass_quality=("f_quality", "sum"))
print()
print("Photometry files with dLight left after each filter (cumulative), totals by area:")
print(flt.groupby(level=0).sum().to_string())
print("per mouse:")
print(flt.to_string())
q = P[P.has_syl & P.has_ophys]
print()
print("Clock and units, photometry sessions with syllables and dLight (n = %d):" % len(q))
print(" median frame interval (s): %.5f to %.5f" % (q.dt_median.min(), q.dt_median.max()))
print(" session length (s): %.0f to %.0f, median %.0f" % (q.duration_s.min(), q.duration_s.max(), q.duration_s.median()))
print(" fraction of dLight samples on a syllable frame: min %.4f" % q.ophys_on_clock.min())
print(" fraction of speed samples on a syllable frame: min %.4f" % q.speed_on_clock.min())
print(" syllable clocks with rows out of time order (sorted on load):", int(q.rotated.sum()))
r = q[q.rotated]
ratio = r.jump_s / r.duration_s
print(" files by number of non-ascending steps:", r.backward_steps.value_counts().sort_index().to_dict(),
"; files with repeated timestamps:", int((r.repeated > 0).sum()), r.loc[r.repeated > 0, ["mouse", "repeated"]].values.tolist())
print(" largest backward step / session length: min %.3f, 5th percentile %.3f, median %.3f; at least 0.95 in %d of %d files"
% (ratio.min(), ratio.quantile(0.05), ratio.median(), (ratio >= 0.95).sum(), len(r)))
print(" median interval between distinct timestamps after sorting: %.5f to %.5f s" % (r.dt_unique.min(), r.dt_unique.max()))
print(" 90th percentile of smoothed dF/F (%%): %.2f to %.2f; dLight-reference r: %.2f to %.2f; max dF/F: %.3f to %.1f"
% (q.signal_max.min(), q.signal_max.max(), q.signal_reference_corr.min(), q.signal_reference_corr.max(),
q.raw_dff_max.min(), q.raw_dff_max.max()))
allsig = np.concatenate([d["signal_dff"][np.isfinite(d["signal_dff"])] for d in phot if d.get("signal_dff") is not None])
allspd = np.concatenate([d["speed_mm_per_frame"][np.isfinite(d["speed_mm_per_frame"])] for d in phot if d.get("speed_mm_per_frame") is not None])
print(" dLight dF/F (unit a.u.) 1st, 50th, 99th percentile: %.4f %.4f %.4f" % tuple(np.percentile(allsig, [1, 50, 99])))
print(" 2D speed (stated mm/frame) 50th, 99th, max: %.2f %.2f %.1f mm/frame = %.0f %.0f %.0f mm/s"
% (*np.percentile(allspd, [50, 99]), allspd.max(), *(30 * np.r_[np.percentile(allspd, [50, 99]), allspd.max()])))
lab = phot[0]["offline"]
labels_attr = None
print(" syllable labels: values 0 to", max(int(d["offline"].max()) for d in phot if d["t"] is not None and len(d["t"])),
"; labels 0-36 are named, 37+ are 'Uncommon Syllable (frequency < 1%)'")
del allsig, allspd
NWB files by session kind and number of subjects:
files subjects
kind
exp 3 1
keypoint 5 5
photometry 264 18
reinforcement 1576 41
velocity 94 9
other assets: 2483 external video files ( {'avi': 2483} )
total size of NWB files: 434 GB, of all assets: 14.8 TB
Photometry files (dLight in file) by area and stimulation:
files mice with_syllables pass_quality
area stim
dls False 224 14 210 194
True 84 8 84 83
dms False 75 8 75 64
files whose syllable series are empty: 14 ; their targets: {(3, 1): 6, (1, 35): 2, (8, 1): 2, (3, 23): 2, (43, 35): 2}
habituation sessions among photometry files: {'dls': 17, 'dms': 1} ; mice: ['dlight-chrimson-1', 'dlight-chrimson-2', 'dlight-chrimson-3', 'dlight-chrimson-4', 'dlight-chrimson-6', 'dms-dlight-6']
Sessions without stimulation that pass the quality criteria, per mouse:
sessions post
area mouse
dls dlight-chrimson-1 16 9
dlight-chrimson-2 14 8
dlight-chrimson-3 14 9
dlight-chrimson-4 11 4
dlight-chrimson-5 13 7
dlight-chrimson-6 13 5
dlight-chrimson-8 13 8
dlight-chrimson-9 12 6
dls-dlight-1 12 0
dls-dlight-2 13 0
dls-dlight-3 13 0
dls-dlight-4 13 0
dls-dlight-5 10 0
dls-dlight-7 10 0
dms dms-dlight-1 2 0
dms-dlight-10 12 0
dms-dlight-11 11 0
dms-dlight-13 7 0
dms-dlight-14 6 0
dms-dlight-4 8 0
dms-dlight-6 6 0
dms-dlight-8 11 0
DLS mice: 14 DMS mice: 8
DLS sessions: 177 of which not post: 121 ; DMS sessions: 63
photometry-kind files without a dLight series: 111 ; mice whose photometry files never hold dLight: ['dls-dlight-6', 'dls-dlight-8', 'dms-dlight-7', 'dms-dlight-9']
Photometry files with dLight left after each filter (cumulative), totals by area:
files no_stim with_syllables not_habituation pass_quality
area
dls 308 224 210 193 177
dms 75 75 75 74 63
per mouse:
files no_stim with_syllables not_habituation pass_quality
area mouse
dls dlight-chrimson-1 30 20 19 16 16
dlight-chrimson-2 28 20 18 14 14
dlight-chrimson-3 30 19 18 14 14
dlight-chrimson-4 30 21 19 15 11
dlight-chrimson-5 24 13 13 13 13
dlight-chrimson-6 29 18 16 14 13
dlight-chrimson-8 29 17 15 15 13
dlight-chrimson-9 30 18 14 14 12
dls-dlight-1 12 12 12 12 12
dls-dlight-2 14 14 14 14 13
dls-dlight-3 15 15 15 15 13
dls-dlight-4 15 15 15 15 13
dls-dlight-5 10 10 10 10 10
dls-dlight-7 12 12 12 12 10
dms dms-dlight-1 10 10 10 10 2
dms-dlight-10 12 12 12 12 12
dms-dlight-11 11 11 11 11 11
dms-dlight-13 7 7 7 7 7
dms-dlight-14 6 6 6 6 6
dms-dlight-4 8 8 8 8 8
dms-dlight-6 10 10 10 9 6
dms-dlight-8 11 11 11 11 11
Clock and units, photometry sessions with syllables and dLight (n = 362):
median frame interval (s): 0.03337 to 0.03337
session length (s): 1186 to 1800, median 1797
fraction of dLight samples on a syllable frame: min 1.0000
fraction of speed samples on a syllable frame: min 0.9957
syllable clocks with rows out of time order (sorted on load): 28
files by number of non-ascending steps: {1: 28} ; files with repeated timestamps: 0 []
largest backward step / session length: min 1.000, 5th percentile 1.000, median 1.000; at least 0.95 in 28 of 28 files
median interval between distinct timestamps after sorting: 0.03337 to 0.03337 s
90th percentile of smoothed dF/F (%): 0.82 to 7.36; dLight-reference r: -0.02 to 1.00; max dF/F: 0.024 to 228.0
dLight dF/F (unit a.u.) 1st, 50th, 99th percentile: -0.0094 0.0130 0.0766
2D speed (stated mm/frame) 50th, 99th, max: 1.08 9.63 84.1 mm/frame = 33 289 2522 mm/s syllable labels: values 0 to 57 ; labels 0-36 are named, 37+ are 'Uncommon Syllable (frequency < 1%)'
In [4]:
# --- Closed-loop sessions ------------------------------------------------------------------
rows = []
for p, d in zip(cl_assets.path, cl):
pon, poff = np.asarray(d.get("pulse_on", [])), np.asarray(d.get("pulse_off", []))
t = d["t"]
rows.append(dict(path=p, mouse=d["subject"], sex=d["sex"], description=d["description"],
day=d["start"].tz_localize(None).normalize(), start=d["start"],
group=d.get("ogen_location"), excitation_nm=d.get("excitation_lambda"),
targets=tuple(d.get("targets", [])), n_pulses=len(pon),
pulse_s=round(float(np.median(poff - pon)), 3) if len(pon) else np.nan,
power_w=d.get("power_w"), n_frames=0 if t is None else len(t),
passes=sd.session_passes_filters(t), rotated=d.get("clock_rotated", False),
backward_steps=d.get("clock_backward_steps", 0), jump_s=d.get("clock_jump_s", 0.0),
repeated=d.get("clock_repeated", 0), dt_unique=d.get("clock_dt_unique", np.nan),
duration_s=np.nan if t is None or len(t) < 2 else t[-1] - t[0],
habituation="habit" in d["description"].lower()))
C = pd.DataFrame(rows)
C["stim"] = C.n_pulses > 0
C["opto"] = np.where(C.group == "snc (axon)", "Opto-DA", np.where(C.group == "ctrl", "control", "?"))
print("Closed-loop files of Opto-DA and control mice:", len(C), "; mice:", C.mouse.nunique())
print(C.groupby("opto").agg(mice=("mouse", "nunique"), files=("path", "size"), stim_files=("stim", "sum")).to_string())
print("pulse lengths (s) by group:", C[C.stim].groupby("opto").pulse_s.apply(lambda x: sorted(set(x))).to_dict())
print("stimulation power (W):", sorted(set(C.loc[C.stim, "power_w"].round(4))), "; stated excitation wavelength (nm):",
sorted(set(C.excitation_nm)))
print("mice with 250-ms sessions:", C[C.stim & (C.pulse_s == 0.25)].groupby("opto").mouse.nunique().to_dict(),
"; mice with pulsed (5 ms) sessions:", C[C.stim & (C.pulse_s < 0.05)].groupby("opto").mouse.unique().apply(list).to_dict())
print("sessions failing the truncation filters:", (~C.passes).sum(), "; clocks out of time order (sorted on load):", int(C.rotated.sum()))
r = C[C.rotated]
ratio = r.jump_s / r.duration_s
print(" files by number of non-ascending steps:", r.backward_steps.value_counts().sort_index().to_dict(),
"; files with repeated timestamps:", int((r.repeated > 0).sum()), r.loc[r.repeated > 0, ["mouse", "repeated", "n_frames"]].values.tolist())
print(" largest backward step / session length: min %.3f, 5th percentile %.3f, median %.3f; at least 0.95 in %d of %d files"
% (ratio.min(), ratio.quantile(0.05), ratio.median(), (ratio >= 0.95).sum(), len(r)))
print(" median interval between distinct timestamps after sorting: %.5f to %.5f s" % (r.dt_unique.min(), r.dt_unique.max()))
mismatch = C[C.stim & ~C.description.str.lower().str.contains("stim")]
print("files with stimulation whose description does not say 'stim':", len(mismatch))
print(mismatch[["mouse", "description", "n_pulses", "pulse_s"]].to_string())
print("files listing two target syllables:", int((C.targets.apply(len) > 1).sum() + (P.targets.apply(len) > 1).sum()))
print("sex by group:", C.drop_duplicates("mouse").groupby(["opto", "sex"]).size().to_dict())
# dLight+ChrimsonR closed-loop files (with and without a dLight series)
chr_cl = [(p, d) for p, d in zip(phot_assets.path, phot) if "ses-reinforcement" in p and d["subject"].startswith("dlight-chrimson")]
chr_cl += list(zip(chr_extra_assets.path, chr_extra))
R = pd.DataFrame([dict(mouse=d["subject"], n_pulses=len(d.get("pulse_on", [])), excitation_nm=d.get("excitation_lambda"),
pulse_s=round(float(np.median(np.asarray(d["pulse_off"]) - np.asarray(d["pulse_on"]))), 3)
if len(d.get("pulse_on", [])) else np.nan) for p, d in chr_cl])
print()
print("dLight+ChrimsonR closed-loop files:", len(R), "(%d without a dLight series); mice:" % len(chr_extra), R.mouse.nunique(),
"; with stimulation:", int((R.n_pulses > 0).sum()))
print(" pulse lengths (s) of files with stimulation:", R[R.n_pulses > 0].pulse_s.value_counts().to_dict(),
"; stated excitation wavelength (nm):", R.excitation_nm.value_counts(dropna=False).to_dict())
Closed-loop files of Opto-DA and control mice: 1340 ; mice: 33
mice files stim_files
opto
Opto-DA 20 805 311
control 13 535 205
pulse lengths (s) by group: {'Opto-DA': [0.005, 0.25, 2.0], 'control': [0.005, 0.25, 1.8, 1.967, 2.0]}
stimulation power (W): [0.01] ; stated excitation wavelength (nm): [635.0]
mice with 250-ms sessions: {'Opto-DA': 20, 'control': 12} ; mice with pulsed (5 ms) sessions: {'Opto-DA': ['3440', '3441', '3442'], 'control': ['snc-dls-ctrl-6', 'vta-nacc-ctrl-6']}
sessions failing the truncation filters: 3 ; clocks out of time order (sorted on load): 547
files by number of non-ascending steps: {1: 545, 2: 1, 3: 1} ; files with repeated timestamps: 2 [['3155', 53822, 107644], ['3155', 53862, 107724]]
largest backward step / session length: min 0.635, 5th percentile 1.000, median 1.000; at least 0.95 in 530 of 547 files
median interval between distinct timestamps after sorting: 0.03337 to 0.03337 s
files with stimulation whose description does not say 'stim': 5
mouse description n_pulses pulse_s
298 15847 Recording session 2 121 0.250
943 dls-ai32jr-1 Recording session 2 34 0.250
1233 vta-nacc-ctrl-6 Recording session 1 136 2.000
1278 vta-nacc-ctrl-6 Recording session 1 4575 0.005
1314 vta-nacc-ctrl-7 Recording session 2 132 0.250
files listing two target syllables: 248
sex by group: {('Opto-DA', 'F'): 11, ('Opto-DA', 'M'): 9, ('control', 'F'): 5, ('control', 'M'): 7, ('control', 'U'): 1}
dLight+ChrimsonR closed-loop files: 236 (6 without a dLight series); mice: 8 ; with stimulation: 88
pulse lengths (s) of files with stimulation: {0.25: 88} ; stated excitation wavelength (nm): {635.0: 236}
Clock checks¶
dLight, speed and syllable labels share one frame clock in the files (the timestamps match, see above), but that does not show that the labels are aligned to the signals in time. Two contrasts check it. (1) The dLight waveform averaged around syllable onsets, with the true labels and with the labels shifted by a random 5 to 15 s within each session; the paper describes a dip around onset followed by a peak. (2) The fraction of stimulation pulses with an online onset of the target syllable in the 0.5 s before, with true and with shifted labels.
In [5]:
rng = np.random.default_rng(1)
# Process dLight for every photometry session with syllables and dLight
for d in phot:
if d["t"] is not None and len(d["t"]) and d["has_ophys"]:
d["z"] = sd.rolling_z(sd.rereference(d["signal_dff"], d["reference_smoothed"]))
lags = np.arange(-15, 31) # frames, -0.5 to +1 s
wf_true, wf_shift = [], []
for (_, row), d in zip(P.iterrows(), phot):
if not (row.use and row.area == "dls"):
continue
for lab, out in [(d["offline"], wf_true), (np.roll(d["offline"], int(rng.integers(150, 451))), wf_shift)]:
on = sd.onsets(lab)
on = on[(on > 15) & (on < len(lab) - 31)]
out.append(np.nanmean(d["z"][on[:, None] + lags[None, :]], axis=0))
wf_true, wf_shift = np.array(wf_true), np.array(wf_shift)
depth = lambda w: w[:, 15:].max(1) - w[:, 10:16].min(1) # peak after onset minus trough around onset
print("onset-aligned mean dLight (z), group average of sessions, at -0.5, -0.1, 0, +0.1, +0.2, +0.3, +0.5, +1.0 s:")
print(" true ", np.round(wf_true.mean(0)[[0, 12, 15, 18, 21, 24, 30, 45]], 3))
print(" shifted ", np.round(wf_shift.mean(0)[[0, 12, 15, 18, 21, 24, 30, 45]], 3))
print(" trough at %.2f s and peak at %.2f s after onset (true labels)" % (lags[np.argmin(wf_true.mean(0)[:25])] / 30, lags[15 + np.argmax(wf_true.mean(0)[15:])] / 30))
print(" dip-to-peak depth: true median %.3f z, shifted %.3f z; true larger in %d of %d sessions (Wilcoxon P = %.1e)"
% (np.median(depth(wf_true)), np.median(depth(wf_shift)), (depth(wf_true) > depth(wf_shift)).sum(), len(wf_true),
stats.wilcoxon(depth(wf_true), depth(wf_shift)).pvalue))
frac_true, frac_shift, frac_rot = [], [], []
for (_, row), d in zip(C.iterrows(), cl):
if not (row.stim and row.pulse_s == 0.25 and row.passes):
continue
frac_rot.append(row.rotated)
for lab, out in [(d["online"], frac_true), (np.roll(d["online"], int(rng.integers(150, 451))), frac_shift)]:
on = sd.onsets(lab)
ton = d["t"][on[lab[on] == row.targets[0]]]
p = np.asarray(d["pulse_on"])
k = np.searchsorted(ton, p) - 1
out.append(np.mean((k >= 0) & (p - ton[np.maximum(k, 0)] < 0.5)))
print("250-ms stimulation sessions (n = %d): fraction of pulses with a target onset in the 0.5 s before, "
"true labels median %.3f, shifted labels median %.3f" % (len(frac_true), np.median(frac_true), np.median(frac_shift)))
frac_true, frac_shift, frac_rot = np.array(frac_true), np.array(frac_shift), np.array(frac_rot)
for name, m in [("clocks stored in time order", ~frac_rot), ("rotated clocks (sorted on load)", frac_rot)]:
print(f" {name}: n = {m.sum()}, true labels median {np.median(frac_true[m]):.3f} (min {frac_true[m].min():.3f}), "
f"shifted labels median {np.median(frac_shift[m]):.3f}")
onset-aligned mean dLight (z), group average of sessions, at -0.5, -0.1, 0, +0.1, +0.2, +0.3, +0.5, +1.0 s: true [-0.037 -0.077 -0.086 -0.063 -0.015 0.008 -0.002 -0.001] shifted [ 0.001 0.001 -0.001 0. -0. -0.001 0. 0.002] trough at 0.00 s and peak at 0.30 s after onset (true labels) dip-to-peak depth: true median 0.118 z, shifted 0.035 z; true larger in 175 of 177 sessions (Wilcoxon P = 9.7e-31) 250-ms stimulation sessions (n = 383): fraction of pulses with a target onset in the 0.5 s before, true labels median 0.938, shifted labels median 0.045 clocks stored in time order: n = 226, true labels median 0.938 (min 0.667), shifted labels median 0.045 rotated clocks (sorted on load): n = 157, true labels median 0.936 (min 0.695), shifted labels median 0.046
The files describe the dLight-only mice with the same subject description as the optogenetic mice ("DAT-IRES-Cre and Ai32"), where the Methods give C57BL/6J for the DLS and DMS dLight mice, and every closed-loop file states a 635-nm stimulator at "snc (axon)", where the Methods give 470-nm light delivered over DLS for the Ai32 (ChR2) mice. Neither is used in the analysis. The files do not state the genotype or opsin of a mouse; Opto-DA and control mice are told apart by the stimulation site ("snc (axon)" or "ctrl"), and the dLight+ChrimsonR mice by their subject names.
In [6]:
# One row per syllable instance of every usable photometry session (offline labels).
inst_list = []
for sid, ((_, row), d) in enumerate(zip(P.iterrows(), phot)):
if not row.use:
continue
tab = sd.syllable_table(d["t"], d["offline"], d["z"])
tab["session"], tab["mouse"], tab["area"], tab["post"] = sid, row.mouse, row.area, row.post
tab["seq_pos"] = np.arange(len(tab))
tab["t_rel"] = tab["t"] - d["t"][0]
# the target syllable(s) of a dLight+ChrimsonR post session
tab["is_post_target"] = tab.syllable.isin(row.targets) if (row.post and row.chrimson) else False
inst_list.append(tab)
inst = pd.concat(inst_list, ignore_index=True)
del inst_list
print("instances:", len(inst), "in", inst.session.nunique(), "sessions;",
"median syllable duration %.3f s (sessions with rotated clocks %.3f s)" % (
inst.duration.median(), inst[inst.session.isin(np.flatnonzero(P.rotated.values))].duration.median()))
print("instances without a dLight peak (frames missing):", int(inst.peak.isna().sum()))
instances: 746625 in 240 sessions; median syllable duration 0.400 s (sessions with rotated clocks 0.400 s) instances without a dLight peak (frames missing): 70152
fig1e-velocity-timescales¶
DLS dopamine fluctuations are significantly negatively correlated with translational velocity at short timescales (under 10 s), and the correlation reverses and becomes positive at longer timescales, for both dLight transient rates and binned amplitudes.
Sessions: DLS, no stimulation, passing the quality criteria, post sessions of dLight+ChrimsonR mice left out (as the paper's code does for this panel).
In [7]:
rng = np.random.default_rng(101)
timescales = np.arange(0.25, 60.1, 0.5)
n_shuf = 100
sess_1e = [i for i, (_, r) in enumerate(P.iterrows()) if r.use and r.area == "dls" and not r.post]
obs = {m: np.full((len(sess_1e), len(timescales)), np.nan) for m in ["amplitude", "rate"]}
shf = {m: np.full((len(sess_1e), len(timescales), n_shuf), np.nan) for m in ["amplitude", "rate"]}
for si, i in enumerate(sess_1e):
d = phot[i]
z = d["z"]
z = (z - np.nanmean(z)) / np.nanstd(z)
v = d["speed_mm_per_frame"]
tt = d["t"] - d["t"][0]
cross = np.r_[False, (z[:-1] < 1.96) & (z[1:] > 1.96)].astype(float)
ok = np.isfinite(z) & np.isfinite(v)
for ti, ts in enumerate(timescales):
b = (tt[ok] // ts).astype(int)
nb = b.max() + 1
cnt = np.bincount(b, minlength=nb)
keep = cnt >= np.bincount(cnt[cnt > 0]).argmax() # drop bins with fewer than the modal number of frames
vb = (np.bincount(b, v[ok], nb) / np.maximum(cnt, 1))[keep]
for m, x in [("amplitude", z[ok]), ("rate", cross[ok])]:
xb = (np.bincount(b, x, nb) / np.maximum(cnt, 1))[keep]
if len(xb) < 4 or xb.std() == 0:
continue
xc, vc = xb - xb.mean(), vb - vb.mean()
obs[m][si, ti] = (xc @ vc) / np.sqrt((xc @ xc) * (vc @ vc))
perm = np.argsort(rng.random((n_shuf, len(xc))), axis=1)
xs = xc[perm]
shf[m][si, ti] = (xs @ vc) / np.sqrt((xc @ xc) * (vc @ vc))
print("sessions:", len(sess_1e), "mice:", P.iloc[sess_1e].mouse.nunique())
res_1e = {}
for m in ["amplitude", "rate"]:
mean_r = np.nanmean(obs[m], axis=0)
null = np.nanmean(shf[m], axis=0) # timescales x shuffles
lo, hi = np.percentile(null, 2.5, axis=1), np.percentile(null, 97.5, axis=1)
res_1e[m] = pd.DataFrame(dict(timescale=timescales, mean_r=mean_r, null_lo=lo, null_hi=hi,
sig_neg=mean_r < lo, sig_pos=mean_r > hi,
n_sessions=np.isfinite(obs[m]).sum(0)))
r = res_1e[m]
short = r[r.timescale < 10]
print(f"\n{m}: mean r at 0.25 s = {r.mean_r.iloc[0]:.3f} (null {r.null_lo.iloc[0]:.3f} to {r.null_hi.iloc[0]:.3f})")
print(f" timescales < 10 s significantly negative: {short.sig_neg.sum()} of {len(short)}; "
f"range of mean r below 10 s: {short.mean_r.min():.3f} to {short.mean_r.max():.3f}")
long_ = r[r.timescale >= 10]
print(f" timescales >= 10 s significantly positive: {long_.sig_pos.sum()} of {len(long_)}; "
f"max mean r {long_.mean_r.max():.3f} at {long_.timescale[long_.mean_r.idxmax()]:.2f} s; "
f"first sign change at {r.timescale[(r.mean_r > 0).idxmax()] if (r.mean_r > 0).any() else np.nan:.2f} s")
print(r.iloc[[0, 1, 2, 4, 9, 19, 29, 39, 59, 79, 99, 119]].round(4).to_string(index=False))
sessions: 121 mice: 14
amplitude: mean r at 0.25 s = -0.109 (null -0.002 to 0.002)
timescales < 10 s significantly negative: 20 of 20; range of mean r below 10 s: -0.201 to -0.109
timescales >= 10 s significantly positive: 0 of 100; max mean r 0.005 at 54.25 s; first sign change at 54.25 s
timescale mean_r null_lo null_hi sig_neg sig_pos n_sessions
0.25 -0.1086 -0.0022 0.0021 True False 121
0.75 -0.1321 -0.0038 0.0030 True False 121
1.25 -0.1472 -0.0041 0.0047 True False 121
2.25 -0.1685 -0.0063 0.0086 True False 121
4.75 -0.1932 -0.0129 0.0119 True False 121
9.75 -0.1902 -0.0194 0.0179 True False 121
14.75 -0.1723 -0.0221 0.0228 True False 121
19.75 -0.1253 -0.0263 0.0304 True False 121
29.75 -0.0959 -0.0381 0.0454 True False 121
39.75 -0.1002 -0.0484 0.0374 True False 121
49.75 -0.0575 -0.0530 0.0503 True False 121
59.75 -0.0281 -0.0701 0.0506 False False 121
rate: mean r at 0.25 s = -0.062 (null -0.002 to 0.002)
timescales < 10 s significantly negative: 20 of 20; range of mean r below 10 s: -0.123 to -0.062
timescales >= 10 s significantly positive: 47 of 100; max mean r 0.113 at 58.75 s; first sign change at 18.75 s
timescale mean_r null_lo null_hi sig_neg sig_pos n_sessions
0.25 -0.0616 -0.0021 0.0020 True False 121
0.75 -0.0889 -0.0039 0.0042 True False 121
1.25 -0.1038 -0.0054 0.0044 True False 121
2.25 -0.1173 -0.0061 0.0078 True False 121
4.75 -0.1196 -0.0099 0.0109 True False 121
9.75 -0.0733 -0.0175 0.0172 True False 121
14.75 -0.0254 -0.0209 0.0264 True False 121
19.75 0.0021 -0.0260 0.0268 False False 121
29.75 0.0410 -0.0353 0.0320 False True 121
39.75 0.0812 -0.0383 0.0515 False True 121
49.75 0.0923 -0.0503 0.0484 False True 121
59.75 0.0795 -0.0529 0.0549 False True 121
In [8]:
# Per mouse (amplitude and rate): mean r over the mouse's sessions at 0.25 s and the timescale of maximum group r
mice_1e = P.iloc[sess_1e].mouse.values
for m in ["amplitude", "rate"]:
tmax = res_1e[m].mean_r.idxmax()
pm = pd.DataFrame({"mouse": mice_1e, "r_0.25s": obs[m][:, 0], f"r_{timescales[tmax]:.2f}s": obs[m][:, tmax]}).groupby("mouse").mean()
print(m, "- mice with negative r at 0.25 s:", (pm.iloc[:, 0] < 0).sum(), "of", len(pm),
"; mice with positive r at the timescale of maximum group r:", (pm.iloc[:, 1] > 0).sum())
print(pm.round(3).T.to_string())
fig, ax = plt.subplots(1, 2, figsize=(9, 3), sharex=True)
for a, m in zip(ax, ["rate", "amplitude"]):
r = res_1e[m]
a.fill_between(r.timescale, r.null_lo, r.null_hi, color="0.8", label="shuffle 95%")
a.plot(r.timescale, r.mean_r, color="C0")
a.axhline(0, color="k", lw=0.5)
a.set(title=f"dLight {m} vs 2D speed", xlabel="timescale (s)", ylabel="mean Pearson r")
print(f"{m}: plotted mean r from {r.mean_r.min():.3f} to {r.mean_r.max():.3f}; all points within axes")
ax[0].legend()
plt.tight_layout(); plt.show()
amplitude - mice with negative r at 0.25 s: 14 of 14 ; mice with positive r at the timescale of maximum group r: 8 mouse dlight-chrimson-1 dlight-chrimson-2 dlight-chrimson-3 dlight-chrimson-4 dlight-chrimson-5 dlight-chrimson-6 dlight-chrimson-8 dlight-chrimson-9 dls-dlight-1 dls-dlight-2 dls-dlight-3 dls-dlight-4 dls-dlight-5 dls-dlight-7 r_0.25s -0.083 -0.185 -0.190 -0.198 -0.042 -0.175 -0.134 -0.186 -0.061 -0.091 -0.051 -0.090 -0.107 -0.086 r_54.25s -0.088 0.266 0.295 -0.138 0.010 -0.071 -0.037 0.272 0.025 -0.016 0.013 -0.198 0.044 0.008 rate - mice with negative r at 0.25 s: 14 of 14 ; mice with positive r at the timescale of maximum group r: 12 mouse dlight-chrimson-1 dlight-chrimson-2 dlight-chrimson-3 dlight-chrimson-4 dlight-chrimson-5 dlight-chrimson-6 dlight-chrimson-8 dlight-chrimson-9 dls-dlight-1 dls-dlight-2 dls-dlight-3 dls-dlight-4 dls-dlight-5 dls-dlight-7 r_0.25s -0.071 -0.073 -0.061 -0.094 -0.072 -0.075 -0.061 -0.105 -0.041 -0.056 -0.033 -0.059 -0.066 -0.050 r_58.75s 0.030 0.031 0.356 0.205 -0.116 0.185 0.138 0.107 -0.029 0.082 0.101 0.189 0.152 0.189 rate: plotted mean r from -0.123 to 0.113; all points within axes amplitude: plotted mean r from -0.201 to 0.005; all points within axes
Result: partially reproduced. Over 121 sessions from 14 mice, both measures are significantly negatively correlated with 2D speed at every timescale below 10 s (mean r -0.201 to -0.109 for the binned amplitude and -0.123 to -0.062 for the transient rate), and every one of the 14 mice is negative at 0.25 s for both. The transient rate reverses: it turns positive from 18.75 s and is significantly positive at 47 of the 100 timescales of 10 s or more (largest 0.113 at 58.75 s). The binned amplitude does not: it is still significantly negative at 49.75 s and is significantly positive at none of the 100 longer timescales (largest mean r 0.005). The claim holds for one of the two measures.
fig1h-syllable-dlight-differs¶
Average syllable-associated dLight amplitudes differ across syllables (Kruskal-Wallis H = 209.29, P < 1e-25, n = 518 mouse-syllable pairs).
DLS sessions without stimulation passing the quality criteria; in dLight+ChrimsonR post sessions the target syllable's instances are left out.
In [9]:
d1h = inst[(inst.area == "dls") & (inst.syllable < sd.N_NAMED_SYLLABLES) & ~inst.is_post_target & inst.peak.notna()]
ms = d1h.groupby(["mouse", "syllable"]).peak.agg(["mean", "size"]).reset_index()
groups = [g["mean"].values for _, g in ms.groupby("syllable")]
H, p = stats.kruskal(*groups)
print(f"mouse-syllable pairs: {len(ms)} ({ms.mouse.nunique()} mice x {ms.syllable.nunique()} syllables); "
f"instances per pair: median {ms['size'].median():.0f}, min {ms['size'].min()}")
print(f"Kruskal-Wallis H = {H:.2f}, P = {p:.2e}")
sy = ms.groupby("syllable")["mean"].agg(["mean", "min", "max"])
print(f"per-syllable mean of mouse averages: {sy['mean'].min():.3f} to {sy['mean'].max():.3f} z")
# Variants: per-mouse KW, and pairs with at least 50 instances
per_mouse = {m: stats.kruskal(*[g.peak.values for _, g in gg.groupby("syllable")]) for m, gg in d1h.groupby("mouse")}
print("per-mouse KW on instances: P < 0.05 in", sum(v.pvalue < 0.05 for v in per_mouse.values()), "of", len(per_mouse),
"mice; H from %.0f to %.0f" % (min(v.statistic for v in per_mouse.values()), max(v.statistic for v in per_mouse.values())))
ms50 = ms[ms["size"] >= 50]
H50, p50 = stats.kruskal(*[g["mean"].values for _, g in ms50.groupby("syllable")])
print(f"pairs with >= 50 instances: n = {len(ms50)}, H = {H50:.2f}, P = {p50:.2e}")
# Variant: mouse-syllable averages over all sessions including the post-session targets
ms_all = inst[(inst.area == "dls") & (inst.syllable < sd.N_NAMED_SYLLABLES)].groupby(["mouse", "syllable"]).peak.mean().reset_index()
Ha, pa = stats.kruskal(*[g.peak.values for _, g in ms_all.groupby("syllable")])
print(f"without leaving out post-session targets: n = {len(ms_all)}, H = {Ha:.2f}, P = {pa:.2e}")
# Variants: by session set, and with the sessions that fail the quality criteria added
def kw_pairs(df):
m = df.groupby(["mouse", "syllable"]).peak.mean().reset_index()
h, pv = stats.kruskal(*[g.peak.values for _, g in m.groupby("syllable")])
return len(m), m.mouse.nunique(), h, pv
sess_path = P.path.values
d1h_ph = d1h[[("ses-photometry" in sess_path[s]) for s in d1h.session]]
print("only ses-photometry files (n pairs, mice, H, P): %d, %d, %.2f, %.2e" % kw_pairs(d1h_ph))
print("only ses-reinforcement files without stimulation: %d, %d, %.2f, %.2e" % kw_pairs(d1h[[("ses-reinforcement" in sess_path[s]) for s in d1h.session]]))
extra = []
for (_, row), d in zip(P.iterrows(), phot):
if row.has_syl and row.has_ophys and row.area == "dls" and not row.stim and not row.habituation and not row.passes:
tab = sd.syllable_table(d["t"], d["offline"], d["z"])
tab = tab[(tab.syllable < sd.N_NAMED_SYLLABLES) & tab.peak.notna()]
if row.post and row.chrimson:
tab = tab[~tab.syllable.isin(row.targets)]
tab["mouse"] = row.mouse
extra.append(tab[["mouse", "syllable", "peak"]])
print(f"DLS sessions without stimulation failing the quality criteria: {len(extra)}")
print("with them added (n pairs, mice, H, P): %d, %d, %.2f, %.2e" % kw_pairs(pd.concat([d1h[["mouse", "syllable", "peak"]]] + extra)))
mouse-syllable pairs: 518 (14 mice x 37 syllables); instances per pair: median 698, min 26 Kruskal-Wallis H = 266.40, P = 5.98e-37 per-syllable mean of mouse averages: 0.675 to 1.236 z per-mouse KW on instances: P < 0.05 in 14 of 14 mice; H from 496 to 1158 pairs with >= 50 instances: n = 515, H = 268.58, P = 2.30e-37 without leaving out post-session targets: n = 518, H = 266.99, P = 4.61e-37
only ses-photometry files (n pairs, mice, H, P): 222, 6, 160.91, 1.01e-17 only ses-reinforcement files without stimulation: 296, 8, 162.01, 6.53e-18 DLS sessions without stimulation failing the quality criteria: 9 with them added (n pairs, mice, H, P): 518, 14, 265.62, 8.40e-37
Result: partially reproduced. The 518 mouse-syllable pairs (14 mice x 37 syllables) match the paper's n exactly, and syllable-associated dLight differs across syllables with Kruskal-Wallis P = 5.98e-37, below the paper's 1e-25. H is 266.40, 27 % above the paper's 209.29 and outside the 10 % tolerance. The variants give 266.99 (target instances kept), 268.58 (pairs with at least 50 instances) and 265.62 (the 9 sessions failing the quality criteria added), so neither the quality filter nor the post-session rule accounts for the difference. Per-syllable means of the mouse averages span 0.675 to 1.236 z, and the per-mouse tests are significant in 14 of 14 mice.
fig1l-identity-not-decodable¶
Random forest classifiers cannot predict which syllable is being expressed from the coincident dLight peak amplitude or waveform.
In [10]:
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedKFold, cross_val_score
rng = np.random.default_rng(301)
d1l = inst[(inst.area == "dls") & (inst.syllable < sd.N_NAMED_SYLLABLES) & ~inst.is_post_target & inst.peak.notna()]
W = np.stack(d1l.waveform.values)
okw = np.isfinite(W).all(1)
d1l, W = d1l[okw], W[okw]
pick = np.concatenate([rng.choice(np.flatnonzero(d1l.syllable.values == s), size=min(1000, (d1l.syllable.values == s).sum()), replace=False)
for s in range(sd.N_NAMED_SYLLABLES)])
y = d1l.syllable.values[pick]
X_peak = d1l.peak.values[pick, None]
X_wave = np.hstack([W[pick], np.diff(W[pick], axis=1)])
print("instances:", len(y), "; per syllable min", np.bincount(y).min(), "; chance 1/37 = %.4f, criterion < %.4f" % (1 / 37, 2 / 37))
cv = StratifiedKFold(5, shuffle=True, random_state=301)
res_1l = {}
for name, X in [("peak amplitude", X_peak), ("waveform + derivative", X_wave)]:
clf = RandomForestClassifier(n_estimators=100, min_samples_leaf=5, n_jobs=4, random_state=301)
acc = cross_val_score(clf, X, y, cv=cv, scoring="accuracy")
acc_sh = cross_val_score(clf, X, rng.permutation(y), cv=cv, scoring="accuracy")
res_1l[name] = acc.mean()
print(f"{name}: held-out accuracy {acc.mean():.4f} (folds {acc.min():.4f} to {acc.max():.4f}); "
f"shuffled labels {acc_sh.mean():.4f}; ratio to chance {acc.mean() * 37:.2f}")
instances: 37000 ; per syllable min 1000 ; chance 1/37 = 0.0270, criterion < 0.0541
peak amplitude: held-out accuracy 0.0271 (folds 0.0258 to 0.0284); shuffled labels 0.0278; ratio to chance 1.00
waveform + derivative: held-out accuracy 0.0323 (folds 0.0300 to 0.0349); shuffled labels 0.0279; ratio to chance 1.20
Result: reproduced. Held-out accuracy for the 37 syllables is 0.0271 from the peak amplitude and 0.0323 from the waveform and its derivative, 1.00 and 1.20 times chance (0.0270) and below the criterion of 0.0541; with shuffled labels it is 0.0278 and 0.0279 (one shuffle each).
fig2b-dlight-usage-across-syllables and fig2i-dlight-entropy-across-syllables¶
Syllables associated on average with more DLS dopamine in an experiment tend to be used more.
Mice string together less predictable sequences after syllables associated with more dopamine; average syllable-associated dLight correlates positively with the entropy of the syllable's outgoing transitions.
Per experiment: mean dLight of instances longer than 0.3 s, z-scored across syllables; counts (run-length onsets) and the entropy of outgoing transitions over all labels. In dLight+ChrimsonR post sessions the target syllable is left out.
In [11]:
def per_experiment_features(inst_area):
rows = []
for sid, g in inst_area.groupby("session"):
lab = g.syllable.values
h = sd.outgoing_entropy(lab, 58)
cnt = np.bincount(lab, minlength=58)
long_ = g[(g.duration > 0.3) & g.peak.notna()]
mpk = long_.groupby("syllable").peak.mean()
named = mpk[mpk.index < sd.N_NAMED_SYLLABLES]
zpk = (named - named.mean()) / named.std()
excl = set(g.loc[g.is_post_target, "syllable"])
for s, v in zpk.items():
if s in excl:
continue
rows.append(dict(session=sid, mouse=g.mouse.iat[0], syllable=s, dlight=v, count=cnt[s], entropy=h[s]))
ef = pd.DataFrame(rows)
mf = ef.groupby(["mouse", "syllable"]).agg(dlight=("dlight", "mean"), count=("count", "mean"),
entropy=("entropy", "mean"), total=("count", "sum")).reset_index()
return ef, mf[mf.total >= 100]
def regression_test(mf, y, rng, n_shuf=200):
r_obs = sd.huber_heldout_r(mf.dlight.values, mf[y].values, rng)
null = np.array([sd.huber_heldout_r(rng.permutation(mf.dlight.values), mf[y].values, rng) for _ in range(n_shuf)])
return r_obs, (null >= r_obs).mean(), null
ef_dls, mf_dls = per_experiment_features(inst[inst.area == "dls"])
print("DLS: experiments", ef_dls.session.nunique(), "; mouse-syllable pairs with >= 100 instances:", len(mf_dls),
"in", mf_dls.mouse.nunique(), "mice")
rng = np.random.default_rng(401)
res_2b = {}
for y in ["count", "entropy"]:
r_obs, p, null = regression_test(mf_dls, y, rng)
res_2b[y] = (r_obs, p)
rp = stats.pearsonr(mf_dls.dlight, mf_dls[y])
p_txt = f"P < {1 / len(null):.3f}" if p == 0 else f"P = {p:.3f}"
print(f"{y}: held-out r = {r_obs:.3f}, shuffle {p_txt} over {len(null)} shuffles (null 95th percentile {np.percentile(null, 95):.3f}); "
f"in-sample Pearson r = {rp[0]:.3f}")
DLS: experiments 177 ; mouse-syllable pairs with >= 100 instances: 507 in 14 mice
count: held-out r = 0.165, shuffle P < 0.005 over 200 shuffles (null 95th percentile 0.051); in-sample Pearson r = 0.196
entropy: held-out r = 0.271, shuffle P < 0.005 over 200 shuffles (null 95th percentile 0.042); in-sample Pearson r = 0.284
In [12]:
# Variants: counts and entropy z-scored within mouse; per-mouse correlations; DMS for reference
rng = np.random.default_rng(402)
mfz = mf_dls.copy()
for y in ["count", "entropy"]:
mfz[y] = mfz.groupby("mouse")[y].transform(lambda x: (x - x.mean()) / x.std())
r_obs, p, _ = regression_test(mfz, y, rng, n_shuf=100)
pm = mf_dls.groupby("mouse").apply(lambda g: stats.spearmanr(g.dlight, g[y])[0])
print(f"{y} z-scored within mouse: held-out r = {r_obs:.3f}, shuffle P = {p:.3f} (0.000 means P < 0.01, 100 shuffles); "
f"per-mouse Spearman rho > 0 in {(pm > 0).sum()} of {len(pm)} mice (median {pm.median():.3f})")
ef_dms, mf_dms = per_experiment_features(inst[inst.area == "dms"])
for y in ["count", "entropy"]:
r_obs, p, _ = regression_test(mf_dms, y, rng, n_shuf=100)
print(f"DMS {y}: pairs {len(mf_dms)}, held-out r = {r_obs:.3f}, shuffle P = {p:.3f} (0.000 means P < 0.01, 100 shuffles)")
fig, ax = plt.subplots(1, 2, figsize=(8, 3))
for a, y in zip(ax, ["count", "entropy"]):
a.scatter(mf_dls.dlight, mf_dls[y], s=6)
a.set(xlabel="syllable-associated dLight (z across syllables)", ylabel=y)
print("plotted ranges: dLight %.2f to %.2f, count %.1f to %.1f, entropy %.2f to %.2f; all points within axes"
% (mf_dls.dlight.min(), mf_dls.dlight.max(), mf_dls["count"].min(), mf_dls["count"].max(), mf_dls.entropy.min(), mf_dls.entropy.max()))
plt.tight_layout(); plt.show()
count z-scored within mouse: held-out r = 0.159, shuffle P = 0.000 (0.000 means P < 0.01, 100 shuffles); per-mouse Spearman rho > 0 in 14 of 14 mice (median 0.179)
entropy z-scored within mouse: held-out r = 0.268, shuffle P = 0.000 (0.000 means P < 0.01, 100 shuffles); per-mouse Spearman rho > 0 in 14 of 14 mice (median 0.291)
DMS count: pairs 250, held-out r = 0.026, shuffle P = 0.170 (0.000 means P < 0.01, 100 shuffles)
DMS entropy: pairs 250, held-out r = 0.252, shuffle P = 0.000 (0.000 means P < 0.01, 100 shuffles) plotted ranges: dLight -1.37 to 1.45, count 9.6 to 259.5, entropy 1.22 to 3.14; all points within axes
Result: fig2b reproduced; fig2i reproduced. Over 507 mouse-syllable pairs in 14 mice, the held-out correlation between syllable-associated dLight and syllable count is r = 0.165, and with outgoing transition entropy r = 0.271; none of the 200 shuffles reached either (P < 0.005; 95th percentiles 0.051 and 0.042). Rank correlations within each mouse are positive in 14 of 14 mice for both. In DMS, for reference, the count relation is not significant (r = 0.026, P = 0.170) and the entropy relation is (r = 0.252, P < 0.01 over 100 shuffles).
fig2e-dlight-predicts-future-usage¶
When an instance of a syllable coincides with a higher-amplitude dopamine transient, that syllable is used more over the following minutes.
Sessions: no stimulation, passing the quality criteria, post sessions of dLight+ChrimsonR mice left out.
In [13]:
n_bins = np.arange(10, 400, 20)
def future_usage(inst_area, rng, n_shuf=100, segments=5):
# Mean over mouse-syllable pairs of the correlation between an instance's dLight peak and
# the count of the same syllable among the next n syllables, with a within-group shuffle null.
g = inst_area.sort_values(["session", "seq_pos"])
sess, lab, pk, pos = g.session.values.astype(np.int64), g.syllable.values.astype(np.int64), g.peak.values, g.seq_pos.values.astype(np.int64)
sess_len = g.groupby("session").seq_pos.transform("max").values + 1
seg = np.minimum((g.t_rel.values / (1800.0 / segments)).astype(int), segments - 1)
named = (lab < sd.N_NAMED_SYLLABLES) & np.isfinite(pk)
# occurrences sorted by (session, syllable, position)
key = sess * 100 + lab
big = 10 ** 6
comb = key * big + pos
order = np.argsort(comb, kind="stable")
comb_sorted = comb[order]
n_upto_self = np.searchsorted(comb_sorted, comb, side="right") # occurrences up to and including i
gid = pd.factorize(sess * 1000 + lab * 10 + seg)[0]
mouse_of_sess = g.groupby("session").mouse.first()
# shuffled peaks within (session, syllable, segment) groups; only finite peaks are permuted,
# among the instances that have one, so that the instances without a peak stay the same
fin = np.isfinite(pk)
by_group = np.lexsort((np.arange(len(gid)), ~fin, gid))
shuffles = []
for _ in range(n_shuf):
o = np.lexsort((rng.random(len(gid)), ~fin, gid))
x = np.empty(len(gid)); x[by_group] = pk[o]
assert np.array_equal(np.isfinite(x), fin)
shuffles.append(x)
out = []
for n in n_bins:
idx = np.flatnonzero(named & (pos + n < sess_len))
cnt = (np.searchsorted(comb_sorted, comb[idx] + n, side="right") - n_upto_self[idx]).astype(float)
gg = pd.factorize(gid[idx])[0]
k = gg.max() + 1
first = pd.Series(np.arange(len(gg))).groupby(gg).first().values
g_sess, g_syl = sess[idx][first], lab[idx][first]
nn = np.bincount(gg, minlength=k).astype(float)
sy, syy = np.bincount(gg, cnt, k), np.bincount(gg, cnt * cnt, k)
def summarize(x):
sx, sxx, sxy = np.bincount(gg, x, k), np.bincount(gg, x * x, k), np.bincount(gg, x * cnt, k)
vx, vy = sxx - sx ** 2 / nn, syy - sy ** 2 / nn
ok = (nn >= 3) & (vx > 1e-9 * np.maximum(sxx, 1)) & (vy > 1e-9 * np.maximum(syy, 1))
r = np.full(k, np.nan)
r[ok] = (sxy - sx * sy / nn)[ok] / np.sqrt(vx[ok] * vy[ok])
n_groups.append(int(ok.sum()))
df = pd.DataFrame(dict(session=g_sess, syllable=g_syl, r=r))
df = df.groupby(["session", "syllable"]).r.mean().reset_index()
df["mouse"] = df.session.map(mouse_of_sess)
return df.groupby(["mouse", "syllable"]).r.mean()
n_groups = []
pp = summarize(pk[idx])
null = np.array([summarize(s_[idx]).mean() for s_ in shuffles])
pm = pp.groupby(level=0).mean()
out.append(dict(n=n, mean_r=pp.mean(), n_pairs=int(pp.notna().sum()), groups_obs=n_groups[0],
groups_null_min=min(n_groups[1:]), groups_null_max=max(n_groups[1:]), null_mean=null.mean(),
null_95=np.percentile(null, 95), p=(null >= pp.mean()).mean(),
mice_positive=int((pm > 0).sum()), mice=len(pm)))
return pd.DataFrame(out)
t = time.time()
rng = np.random.default_rng(501)
fu_dls = future_usage(inst[(inst.area == "dls") & ~inst.post], rng)
print(f"computed in {time.time() - t:.0f} s")
print(fu_dls.round(4).to_string(index=False))
sel = fu_dls[fu_dls.n <= 190]
print(f"n from 10 to 190: significant (P < 0.05) at {(sel.p < 0.05).sum()} of {len(sel)}; "
f"mean r {sel.mean_r.min():.4f} to {sel.mean_r.max():.4f}; all n: significant at {(fu_dls.p < 0.05).sum()} of {len(fu_dls)}")
computed in 13 s n mean_r n_pairs groups_obs groups_null_min groups_null_max null_mean null_95 p mice_positive mice 10 0.0251 516 15480 15480 15480 0.0004 0.0076 0.00 11 14 30 0.0221 516 17539 17539 17539 0.0003 0.0068 0.00 12 14 50 0.0182 518 18163 18163 18163 0.0002 0.0075 0.00 11 14 70 0.0135 518 18485 18485 18485 0.0004 0.0072 0.00 12 14 90 0.0100 518 18649 18649 18649 0.0011 0.0075 0.00 10 14 110 0.0071 518 18744 18744 18744 0.0010 0.0087 0.09 10 14 130 0.0044 518 18788 18788 18788 0.0006 0.0083 0.15 8 14 150 0.0044 518 18818 18818 18818 0.0006 0.0072 0.18 8 14 170 0.0072 518 18828 18828 18828 0.0006 0.0065 0.05 9 14 190 0.0099 518 18834 18834 18834 0.0007 0.0068 0.00 11 14 210 0.0099 518 18822 18822 18822 0.0010 0.0066 0.01 9 14 230 0.0103 518 18798 18798 18798 0.0008 0.0065 0.01 11 14 250 0.0095 518 18754 18754 18754 0.0006 0.0066 0.01 10 14 270 0.0090 518 18744 18744 18744 0.0006 0.0066 0.02 11 14 290 0.0107 518 18686 18686 18686 0.0006 0.0061 0.02 10 14 310 0.0118 518 18642 18642 18642 0.0008 0.0066 0.01 10 14 330 0.0096 518 18596 18596 18596 0.0008 0.0069 0.01 8 14 350 0.0065 518 18522 18522 18522 0.0008 0.0072 0.08 8 14 370 0.0047 518 18476 18476 18476 0.0008 0.0069 0.15 8 14 390 0.0043 518 18395 18395 18395 0.0007 0.0066 0.23 9 14 n from 10 to 190: significant (P < 0.05) at 6 of 10; mean r 0.0044 to 0.0251; all n: significant at 13 of 20
Result: partially reproduced. The mean correlation between an instance's dLight and the use of the same syllable over the next n syllables is positive at every n (0.0043 to 0.0251) and above the shuffles from n = 10 to 90 (P < 0.01), with 11 of 14 mice positive at n = 10. It is not significant at n = 110, 130, 150 and 170 (P = 0.09, 0.15, 0.18 and 0.05), so it holds at 6 of the 10 values from 10 to 190 and the criterion's "throughout" fails; it is significant again from 190 to 330. The shuffles permute only the instances that have a dLight peak, and the same number of groups enters the observed value and every shuffle (15,480 at n = 10). The effect is small (r = 0.0251 at n = 10).
ed7-dms-fewer-transients and ed7-dms-no-usage-prediction¶
dLight recordings in dorsomedial striatum (DMS) show less frequent dopamine transients than in DLS.
DMS dopamine transients do not predict future syllable usage.
In [14]:
def transient_rate(z, threshold=1.0, schmitt=3, schmitt_threshold=0.25):
z = (z - np.nanmean(z)) / np.nanstd(z)
zz = np.where(np.isfinite(z), z, -np.inf)
c = np.flatnonzero((zz[:-1] < threshold) & (zz[1:] > threshold))
good = [not (zz[k + 1:k + schmitt + 1] < schmitt_threshold).any() for k in c]
minutes = np.isfinite(z).sum() / sd.FS / 60
return np.sum(good) / minutes
rows = []
for (_, r), d in zip(P.iterrows(), phot):
if r.use:
rows.append(dict(mouse=r.mouse, area=r.area, rate=transient_rate(d["z"]), rate_196=transient_rate(d["z"], 1.96)))
tr = pd.DataFrame(rows)
tm = tr.groupby(["area", "mouse"]).mean().reset_index()
for col in ["rate", "rate_196"]:
a, b = tm.loc[tm.area == "dms", col], tm.loc[tm.area == "dls", col]
u = stats.mannwhitneyu(a, b, alternative="two-sided")
print(f"{col}: DMS median {a.median():.2f}/min ({len(a)} mice, {a.min():.2f} to {a.max():.2f}), "
f"DLS median {b.median():.2f}/min ({len(b)} mice, {b.min():.2f} to {b.max():.2f}); U = {u.statistic:.0f}, P = {u.pvalue:.4f}")
print(tm.round(2).to_string(index=False))
rate: DMS median 56.48/min (8 mice, 48.57 to 59.94), DLS median 62.12/min (14 mice, 56.33 to 73.97); U = 13, P = 0.0022 rate_196: DMS median 22.10/min (8 mice, 20.54 to 24.12), DLS median 25.21/min (14 mice, 21.53 to 28.08); U = 9, P = 0.0006 area mouse rate rate_196 dls dlight-chrimson-1 64.02 25.37 dls dlight-chrimson-2 56.33 21.53 dls dlight-chrimson-3 58.32 24.82 dls dlight-chrimson-4 56.55 21.74 dls dlight-chrimson-5 73.97 24.80 dls dlight-chrimson-6 56.64 24.40 dls dlight-chrimson-8 59.60 24.48 dls dlight-chrimson-9 60.33 25.06 dls dls-dlight-1 67.03 27.78 dls dls-dlight-2 64.29 26.36 dls dls-dlight-3 66.92 28.08 dls dls-dlight-4 62.01 25.62 dls dls-dlight-5 62.22 27.47 dls dls-dlight-7 62.37 26.31 dms dms-dlight-1 59.94 23.29 dms dms-dlight-10 57.14 23.35 dms dms-dlight-11 53.47 20.54 dms dms-dlight-13 56.48 24.12 dms dms-dlight-14 48.57 21.51 dms dms-dlight-4 56.80 21.19 dms dms-dlight-6 56.48 22.51 dms dms-dlight-8 55.05 21.69
In [15]:
t = time.time()
rng = np.random.default_rng(601)
fu_dms = future_usage(inst[inst.area == "dms"], rng)
print(f"computed in {time.time() - t:.0f} s")
print(fu_dms.round(4).to_string(index=False))
sel = fu_dms[fu_dms.n <= 190]
print(f"DMS, n from 10 to 190: significant (P < 0.05) at {(sel.p < 0.05).sum()} of {len(sel)}; mean r {sel.mean_r.min():.4f} to {sel.mean_r.max():.4f}")
sel_dls = fu_dls[fu_dls.n <= 190]
print(f"DLS for comparison: mean r {sel_dls.mean_r.min():.4f} to {sel_dls.mean_r.max():.4f}")
fig, ax = plt.subplots(figsize=(5, 3))
for lbl, f in [("DLS", fu_dls), ("DMS", fu_dms)]:
ax.plot(f.n, f.mean_r, label=lbl)
ax.plot(f.n, f.null_95, ls=":", color=ax.lines[-1].get_color())
ax.axhline(0, color="k", lw=0.5); ax.set(xlabel="next n syllables", ylabel="mean r (dLight, future count)"); ax.legend()
print("plotted: mean r %.4f to %.4f; all points within axes" % (min(fu_dls.mean_r.min(), fu_dms.mean_r.min()), max(fu_dls.mean_r.max(), fu_dms.mean_r.max())))
plt.tight_layout(); plt.show()
computed in 8 s n mean_r n_pairs groups_obs groups_null_min groups_null_max null_mean null_95 p mice_positive mice 10 0.0331 289 7506 7506 7506 0.0008 0.0155 0.00 8 8 30 0.0318 292 8578 8578 8578 0.0007 0.0140 0.00 8 8 50 0.0258 292 8907 8907 8907 0.0010 0.0114 0.00 8 8 70 0.0269 293 9095 9095 9095 0.0009 0.0117 0.00 8 8 90 0.0207 293 9169 9169 9169 0.0003 0.0127 0.00 7 8 110 0.0123 293 9195 9195 9195 0.0001 0.0111 0.04 5 8 130 0.0139 293 9199 9199 9199 0.0003 0.0114 0.01 5 8 150 0.0142 293 9188 9188 9188 0.0006 0.0108 0.01 6 8 170 0.0110 293 9163 9163 9163 0.0008 0.0122 0.07 7 8 190 0.0182 293 9127 9127 9127 0.0006 0.0118 0.00 6 8 210 0.0193 293 9087 9087 9087 0.0002 0.0129 0.00 7 8 230 0.0135 293 9063 9063 9063 0.0010 0.0118 0.04 7 8 250 0.0107 293 9018 9018 9018 0.0015 0.0126 0.12 7 8 270 0.0103 293 8971 8971 8971 0.0016 0.0154 0.14 7 8 290 0.0120 293 8939 8939 8939 0.0017 0.0132 0.07 5 8 310 0.0098 293 8880 8880 8880 0.0014 0.0133 0.14 5 8 330 0.0123 293 8813 8813 8813 0.0011 0.0137 0.07 6 8 350 0.0088 293 8743 8743 8743 0.0015 0.0128 0.16 5 8 370 0.0106 293 8670 8670 8670 0.0015 0.0134 0.11 6 8 390 0.0096 293 8571 8571 8571 0.0011 0.0135 0.14 6 8 DMS, n from 10 to 190: significant (P < 0.05) at 9 of 10; mean r 0.0110 to 0.0331 DLS for comparison: mean r 0.0044 to 0.0251 plotted: mean r 0.0043 to 0.0331; all points within axes
Result: ed7-dms-fewer-transients reproduced; ed7-dms-no-usage-prediction not reproduced. DMS mice have fewer transients than DLS mice: a median 56.48 against 62.12 per minute (U = 13, P = 0.0022), and 22.10 against 25.21 with a 1.96 threshold (P = 0.0006). By the test of fig2e-dlight-predicts-future-usage, DMS dLight does predict the later use of the same syllable: the mean r is 0.0110 to 0.0331 for n from 10 to 190, significant at 9 of those 10 values, positive in 8 of 8 mice at n = 10, and larger than in DLS (0.0044 to 0.0251). The paper's Extended Data Fig. 7 has no lagged future-usage test; it supports the claim with the encoding model (panel e) and the title "DMS dopamine does not correlate with syllable usage or entropy". The across-syllable count test of fig2b, run on DMS above, agrees with the paper (held-out r = 0.026, shuffle P = 0.170). That this attempt's lagged test finds DMS stronger than DLS suggests that it differs from the paper's (equal time segments, no z-scoring to the shuffles); the outcome follows the pre-set criterion.
fig3e-optoda-learning, fig3f-cumulative-counts, fig3e-learner-count, fig3h-stim-post-correlation¶
Pairing a target syllable with optogenetic dopamine release (Opto-DA) increases its use relative to baseline more than in no-opsin controls (log2 fold change, one-sided Mann-Whitney U = 197, P = 0.002, f = 0.82).
The cumulative increase in target counts relative to baseline is larger in Opto-DA mice than in controls (U = 184, P = 0.007, f = 0.77).
Nine Opto-DA mice exceed the largest change in target counts of any control ("learners").
In learners, the change in target usage during stimulation predicts the change in the post-stimulation sessions (r = 0.89, P = 0.005).
In [16]:
counts = {}
for i, ((_, r), d) in enumerate(zip(C.iterrows(), cl)):
if r.passes and d["t"] is not None:
counts[i] = d
rows = []
for (mouse, day), g in C[C.stim & (C.pulse_s == 0.25)].groupby(["mouse", "day"]):
g = g[g.passes]
target = g.targets.iat[0][0] if len(g) else None
if len(g) == 0 or not all(t[0] == target for t in g.targets):
continue
mc = C[(C.mouse == mouse) & ~C.stim & ~C.habituation & C.passes]
bday = sd.monday_baseline(day, mc.day)
if bday is None:
continue
base = mc[mc.day == bday]
post = mc[mc.day == day + pd.Timedelta(days=1)]
cnt = lambda idx: [int(sd.binned_onset_counts(cl[i]["t"], cl[i]["online"], target).sum()) for i in idx]
S, B, Q = cnt(g.index), cnt(base.index), cnt(post.index)
b = np.mean(B)
rows.append(dict(mouse=mouse, opto=g.opto.iat[0], day=day, target=target, n_stim=len(S), stim=sum(S),
baseline_day=bday, baseline_monday=bday.dayofweek == 0, n_base=len(B), baseline_mean=b,
n_post=len(Q), post=sum(Q) if len(Q) else np.nan,
log2fc=np.log2((sum(S) + 2) / (len(S) * b + 2)), change=sum(S) - len(S) * b,
post_change=(sum(Q) - len(Q) * b) if len(Q) else np.nan,
post_log2fc=np.log2((sum(Q) + 2) / (len(Q) * b + 2)) if len(Q) else np.nan))
T3 = pd.DataFrame(rows)
T3 = T3[T3.n_stim >= 2]
keep = T3.groupby("mouse").target.transform("size") >= 2
print("targets dropped for fewer than 2 stimulation sessions or mice with < 2 targets:",
int((~keep).sum()))
T3 = T3[keep]
print("targets:", len(T3), "; per mouse:", T3.groupby("mouse").size().value_counts().to_dict(),
"; baseline day a Monday for", int(T3.baseline_monday.sum()), "; with post sessions:", int(T3.n_post.gt(0).sum()))
M3 = T3.groupby(["opto", "mouse"]).agg(log2fc=("log2fc", "mean"), change=("change", "mean"),
post_change=("post_change", "mean"), n_targets=("target", "size")).reset_index()
print("mice:", M3.groupby("opto").size().to_dict())
res3 = {}
for col, cid in [("log2fc", "fig3e"), ("change", "fig3f")]:
a, b = M3.loc[M3.opto == "Opto-DA", col], M3.loc[M3.opto == "control", col]
u = stats.mannwhitneyu(a, b, alternative="greater")
res3[col] = u
print(f"{cid} {col}: Opto-DA median {a.median():.3f} (n={len(a)}), control median {b.median():.3f} (n={len(b)}); "
f"U = {u.statistic:.0f}, one-sided P = {u.pvalue:.4f}, f = {u.statistic / (len(a) * len(b)):.2f}")
targets dropped for fewer than 2 stimulation sessions or mice with < 2 targets: 0
targets: 191 ; per mouse: {6: 31, 5: 1} ; baseline day a Monday for 167 ; with post sessions: 184
mice: {'Opto-DA': 20, 'control': 12}
fig3e log2fc: Opto-DA median 0.330 (n=20), control median 0.081 (n=12); U = 198, one-sided P = 0.0013, f = 0.82
fig3f change: Opto-DA median 27.583 (n=20), control median -2.333 (n=12); U = 182, one-sided P = 0.0083, f = 0.76
In [17]:
# Variants: per-target values pooled; baseline as the mean of all baseline days in the 5 days before
for col in ["log2fc", "change"]:
a, b = T3.loc[T3.opto == "Opto-DA", col], T3.loc[T3.opto == "control", col]
u = stats.mannwhitneyu(a, b, alternative="greater")
print(f"per target {col}: Opto-DA {len(a)} targets median {a.median():.3f}, control {len(b)} median {b.median():.3f}; "
f"U = {u.statistic:.0f}, P = {u.pvalue:.2e}")
print(M3.sort_values(["opto", "log2fc"]).round(3).to_string(index=False))
fig, ax = plt.subplots(1, 2, figsize=(6, 3))
for a_, col in zip(ax, ["log2fc", "change"]):
for k, grp in enumerate(["Opto-DA", "control"]):
v = M3.loc[M3.opto == grp, col]
a_.scatter(np.full(len(v), k) + np.random.default_rng(1).uniform(-0.1, 0.1, len(v)), v, s=10)
a_.set_xticks([0, 1]); a_.set_xticklabels(["Opto-DA", "control"]); a_.set_ylabel(col); a_.axhline(0, color="k", lw=0.5)
print("plotted ranges: log2fc %.2f to %.2f, change %.1f to %.1f; all points within axes"
% (M3.log2fc.min(), M3.log2fc.max(), M3.change.min(), M3.change.max()))
plt.tight_layout(); plt.show()
per target log2fc: Opto-DA 119 targets median 0.224, control 72 median 0.062; U = 5377, P = 1.59e-03 per target change: Opto-DA 119 targets median 37.000, control 72 median 8.000; U = 5298, P = 3.10e-03 opto mouse log2fc change post_change n_targets Opto-DA 1546 -0.109 -16.000 -32.333 6 Opto-DA 410 0.099 28.000 -7.200 6 Opto-DA 241 0.100 -4.667 -61.667 6 Opto-DA dls-ai32jr-5 0.132 4.000 -59.833 6 Opto-DA 3441 0.139 -9.667 4.500 6 Opto-DA 3440 0.142 1.333 -9.667 6 Opto-DA 239 0.210 27.167 -20.833 6 Opto-DA dls-ai32jr-2 0.217 5.200 -32.000 5 Opto-DA 15847 0.226 27.167 -9.167 6 Opto-DA dls-ai32jr-3 0.316 77.667 29.000 6 Opto-DA dls-ai32jr-4 0.344 1.833 -24.833 6 Opto-DA 242 0.358 6.833 -16.500 6 Opto-DA 1529 0.391 92.667 82.000 6 Opto-DA 3442 0.398 47.333 -9.500 6 Opto-DA dls-ai32jr-1 0.441 41.667 37.500 6 Opto-DA 15823 0.536 68.000 49.667 6 Opto-DA 240 0.586 116.000 44.333 6 Opto-DA 1544 0.660 179.833 145.333 6 Opto-DA 15822 0.678 64.000 33.000 6 Opto-DA 1528 0.828 170.333 111.833 6 control snc-dls-ctrl-6 -0.111 -66.333 -71.417 6 control 209 -0.081 -2.833 -37.600 6 control snc-dls-ctrl-7 -0.033 -3.833 -39.500 6 control 1561 -0.033 -11.000 2.333 6 control 1527 -0.026 -5.333 44.000 6 control 197 0.044 25.833 9.200 6 control 200 0.118 35.333 -36.000 6 control 208 0.153 -1.833 -7.400 6 control 211 0.173 21.833 -16.400 6 control vta-nacc-ctrl-7 0.189 1.000 -6.833 6 control 10 0.194 -2.833 -38.200 6 control vta-nacc-ctrl-6 0.341 54.500 59.000 6 plotted ranges: log2fc -0.11 to 0.83, change -66.3 to 179.8; all points within axes
Result: fig3e-optoda-learning reproduced; fig3f-cumulative-counts reproduced. The files give the paper's 20 Opto-DA and 12 control mice with 250-ms sessions. For the log2 fold change in target counts, U = 198, one-sided P = 0.0013, f = 0.82 (paper 197, 0.002, 0.82); for the change in counts, U = 182, P = 0.0083, f = 0.76 (paper 184, 0.007, 0.77). Both U are within tolerance.
In [18]:
thr = M3.loc[M3.opto == "control", "change"].max()
learners = M3[(M3.opto == "Opto-DA") & (M3.change > thr)]
print(f"largest control change: {thr:.1f}; learners: {len(learners)} of {(M3.opto == 'Opto-DA').sum()} Opto-DA mice")
print("learners:", sorted(learners.mouse))
# The nine learners printed by the paper's code (fig_3e ... optoda learning.ipynb, commit 14ed714)
code_learners = ["1528", "1529", "1544", "15822", "15823", "240", "3442", "dls-ai32jr-1", "dls-ai32jr-3"]
print("in both:", sorted(set(learners.mouse) & set(code_learners)), "; only here:", sorted(set(learners.mouse) - set(code_learners)),
"; only in the code's list:", sorted(set(code_learners) - set(learners.mouse)))
# Threshold variant: learners by log2 fold change above the control maximum
thr2 = M3.loc[M3.opto == "control", "log2fc"].max()
print(f"variant, by log2 fold change (control max {thr2:.3f}): {((M3.opto == 'Opto-DA') & (M3.log2fc > thr2)).sum()} learners")
# Per target, for the two code learners not found here and the control that sets the threshold
print("per target: stimulation counts (two sessions), baseline day and its mean count, change")
cols = ["mouse", "opto", "day", "target", "n_stim", "stim", "baseline_day", "n_base", "baseline_mean", "change", "log2fc"]
print(T3[T3.mouse.isin(["3442", "dls-ai32jr-1", "vta-nacc-ctrl-6"])][cols].round(2).to_string(index=False))
largest control change: 54.5; learners: 7 of 20 Opto-DA mice
learners: ['1528', '1529', '1544', '15822', '15823', '240', 'dls-ai32jr-3']
in both: ['1528', '1529', '1544', '15822', '15823', '240', 'dls-ai32jr-3'] ; only here: [] ; only in the code's list: ['3442', 'dls-ai32jr-1']
variant, by log2 fold change (control max 0.341): 10 learners
per target: stimulation counts (two sessions), baseline day and its mean count, change
mouse opto day target n_stim stim baseline_day n_base baseline_mean change log2fc
3442 Opto-DA 2021-01-19 8 2 1087 2021-01-18 2 529.0 29.0 0.04
3442 Opto-DA 2021-01-21 1 2 350 2021-01-18 2 101.0 148.0 0.79
3442 Opto-DA 2021-01-26 23 2 231 2021-01-25 2 80.0 71.0 0.52
3442 Opto-DA 2021-01-28 3 2 386 2021-01-25 2 217.5 -49.0 -0.17
3442 Opto-DA 2021-02-02 43 2 156 2021-02-01 2 39.5 77.0 0.96
3442 Opto-DA 2021-02-04 35 2 49 2021-02-01 2 20.5 8.0 0.25
dls-ai32jr-1 Opto-DA 2021-09-21 1 2 124 2021-09-20 2 15.5 93.0 1.93
dls-ai32jr-1 Opto-DA 2021-09-23 8 2 1096 2021-09-20 2 509.0 78.0 0.11
dls-ai32jr-1 Opto-DA 2021-09-30 3 2 353 2021-09-29 2 189.5 -26.0 -0.10
dls-ai32jr-1 Opto-DA 2021-10-14 23 2 287 2021-10-13 2 110.5 66.0 0.37
dls-ai32jr-1 Opto-DA 2021-10-19 35 2 128 2021-10-18 2 57.5 13.0 0.15
dls-ai32jr-1 Opto-DA 2021-10-21 43 2 214 2021-10-18 2 94.0 26.0 0.19
vta-nacc-ctrl-6 control 2021-01-19 8 2 841 2021-01-18 2 357.0 127.0 0.24
vta-nacc-ctrl-6 control 2021-01-21 1 2 219 2021-01-18 2 92.5 34.0 0.24
vta-nacc-ctrl-6 control 2021-01-26 23 2 128 2021-01-25 2 49.0 30.0 0.38
vta-nacc-ctrl-6 control 2021-01-28 3 2 176 2021-01-25 2 113.0 -50.0 -0.36
vta-nacc-ctrl-6 control 2021-02-04 35 2 279 2021-02-01 2 63.5 152.0 1.12
vta-nacc-ctrl-6 control 2021-02-09 43 2 132 2021-02-08 2 49.0 34.0 0.42
/tmp/ipykernel_2184/211501021.py:15: UserWarning: obj.round has no effect with datetime, timedelta, or period dtypes. Use obj.dt.round(...) instead. print(T3[T3.mouse.isin(["3442", "dls-ai32jr-1", "vta-nacc-ctrl-6"])][cols].round(2).to_string(index=False))
In [19]:
rng = np.random.default_rng(701)
def r_shuffle(x, y, rng, n=1000):
r = stats.pearsonr(x, y)[0]
null = np.array([stats.pearsonr(x, rng.permutation(y))[0] for _ in range(n)])
return r, (null >= r).mean(), (np.abs(null) >= abs(r)).mean()
for name, sel in [("learners (this attempt)", M3.mouse.isin(learners.mouse)),
("learners (code's list)", M3.mouse.isin(code_learners)),
("controls", M3.opto == "control")]:
s = M3[sel].dropna(subset=["post_change"])
r, p, p2 = r_shuffle(s.post_change.values, s.change.values, rng)
print(f"{name}: n = {len(s)}, r = {r:.3f}, one-sided shuffle P = {p:.3f} (two-sided {p2:.3f}); "
f"stim change {s.change.min():.1f} to {s.change.max():.1f}, post change {s.post_change.min():.1f} to {s.post_change.max():.1f}")
learners (this attempt): n = 7, r = 0.891, one-sided shuffle P = 0.017 (two-sided 0.017); stim change 64.0 to 179.8, post change 29.0 to 145.3
learners (code's list): n = 9, r = 0.888, one-sided shuffle P = 0.001 (two-sided 0.001); stim change 41.7 to 179.8, post change -9.5 to 145.3
controls: n = 12, r = 0.591, one-sided shuffle P = 0.017 (two-sided 0.039); stim change -66.3 to 54.5, post change -71.4 to 59.0
Result: fig3e-learner-count not reproduced; fig3h-stim-post-correlation reproduced. Seven of 20 Opto-DA mice exceed the largest control change (54.5, in vta-nacc-ctrl-6), where the paper has nine. All seven are among the nine the paper's code printed; the other two (3442 at 47.3 and dls-ai32jr-1 at 41.7) fall below that control. Their per-target changes range from -49 to 148 (3442) and -26 to 93 (dls-ai32jr-1), against -50 to 152 for vta-nacc-ctrl-6, with the baseline day 1 to 3 days before the stimulation day in every case. Counting mice above the control maximum of the log2 fold change gives 10. Among the seven learners, the change during stimulation and in the post sessions correlate with r = 0.891 (one-sided shuffle P = 0.017), within tolerance of the paper's 0.89; with the code's nine learners r = 0.888 (P = 0.001). The controls also correlate (r = 0.591, one-sided P = 0.017, two-sided 0.039), where the paper reports P = 0.082.
fig3l-no-kinematic-effect and fig3m-long-stim-velocity¶
Brief (250 ms) Opto-DA does not change movement at the moment of stimulation; no comparison between stimulation and catch trials in any mouse is significant.
Extending stimulation to 3 s (pulsed, 25 Hz) makes mice increase their velocity.
In [20]:
LAT = 0.17
def speed_change(d, align, post_s):
t, v = d["t"], d["speed_mm_per_frame"] * sd.FS
out = []
for a in align:
pre = v[(t >= a - 1) & (t < a)]
post = v[(t >= a) & (t < a + post_s)]
out.append(np.nanmean(post) - np.nanmean(pre) if np.isfinite(pre).any() and np.isfinite(post).any() else np.nan)
return np.array(out)
def trials(r, d, pulsed, post_s, catch_exclusion=None, isolate=None):
# isolate: keep only catch trials with no stimulation pulse (any pulse of a train) within this many s
p_all = np.asarray(d["pulse_on"])
p = p_all[np.r_[True, np.diff(p_all) > 1.0]] if pulsed else p_all
lab = d["online"]
on = sd.onsets(lab)
ton = d["t"][on[lab[on] == r.targets[0]]]
k = np.searchsorted(p, ton)
nxt = np.where(k < len(p), p[np.minimum(k, len(p) - 1)], np.inf)
catch = ton[(nxt - ton) > 0.5]
if catch_exclusion is not None:
prv = np.where(k > 0, p[np.maximum(k - 1, 0)], -np.inf)
catch = catch[(catch - prv[(nxt - ton) > 0.5]) > catch_exclusion]
if isolate is not None:
j = np.searchsorted(p_all, catch)
dprev = np.where(j > 0, catch - p_all[np.maximum(j - 1, 0)], np.inf)
dnext = np.where(j < len(p_all), p_all[np.minimum(j, len(p_all) - 1)] - catch, np.inf)
catch = catch[np.minimum(dprev, dnext) >= isolate]
return speed_change(d, p, post_s), speed_change(d, catch + LAT, post_s)
def per_mouse(pulse_sel, pulsed, post_s, catch_exclusion=None, isolate=None):
rows = []
for mouse, g in C[C.stim & pulse_sel & C.passes].groupby("mouse"):
st, ca, frac = [], [], []
for i, r in g.iterrows():
a, b = trials(r, cl[i], pulsed, post_s, catch_exclusion, isolate)
st.append(a); ca.append(b); frac.append(len(a) / (len(a) + len(b)))
st, ca = np.concatenate(st), np.concatenate(ca)
stim_frac = len(st) / (len(st) + len(ca)) # all detections of the mouse, before dropping trials without speed
st, ca = st[np.isfinite(st)], ca[np.isfinite(ca)]
u = stats.mannwhitneyu(st, ca, alternative="greater")
rows.append(dict(mouse=mouse, opto=g.opto.iat[0], sessions=len(g), n_stim=len(st), n_catch=len(ca),
stim_frac=stim_frac, stim_frac_min=min(frac), stim_frac_max=max(frac),
stim_mean=st.mean(), catch_mean=ca.mean(), diff=st.mean() - ca.mean(), U=u.statistic, p=u.pvalue))
return pd.DataFrame(rows)
from statsmodels.stats.multitest import multipletests
k3l = per_mouse(C.pulse_s == 0.25, False, 1.0)
k3l = k3l[k3l.mouse.isin(M3.mouse)]
for grp in ["Opto-DA", "control"]:
s = k3l[k3l.opto == grp]
holm = multipletests(s.p, method="holm")[1]
print(f"{grp}: {len(s)} mice; P < 0.05 in {(s.p < 0.05).sum()} (Holm-adjusted: {(holm < 0.05).sum()}); "
f"smallest P {s.p.min():.4f}; stim minus catch speed change, median {s['diff'].median():.2f} mm/s "
f"(range {s['diff'].min():.2f} to {s['diff'].max():.2f}); trials per mouse: stim {s.n_stim.min()}-{s.n_stim.max()}, catch {s.n_catch.min()}-{s.n_catch.max()}")
print(k3l.round(4).to_string(index=False))
# Fraction of target detections that were stimulated under the substitute (stimulations / (stimulations + catch)),
# per session; the Methods give stimulation on 75% of detections, so true catch trials are about 25%.
print("stimulated fraction of detections, 250 ms: per mouse median %.3f (range %.3f to %.3f); per session %.3f to %.3f (Methods: 0.75)"
% (k3l.stim_frac.median(), k3l.stim_frac.min(), k3l.stim_frac.max(), k3l.stim_frac_min.min(), k3l.stim_frac_max.max()))
# Variant: catch trials with no stimulation pulse within 2 s either side
k3l_iso = per_mouse(C.pulse_s == 0.25, False, 1.0, isolate=2.0)
k3l_iso = k3l_iso[k3l_iso.mouse.isin(M3.mouse)]
for grp in ["Opto-DA", "control"]:
s = k3l_iso[k3l_iso.opto == grp]
print(f"variant, catch >= 2 s from any pulse, {grp}: P < 0.05 in {(s.p < 0.05).sum()} of {len(s)} mice; "
f"stim minus catch median {s['diff'].median():.2f} mm/s (range {s['diff'].min():.2f} to {s['diff'].max():.2f}); "
f"catch trials per mouse {s.n_catch.min()}-{s.n_catch.max()}")
Opto-DA: 20 mice; P < 0.05 in 9 (Holm-adjusted: 6); smallest P 0.0000; stim minus catch speed change, median -0.07 mm/s (range -15.28 to 7.16); trials per mouse: stim 871-1687, catch 602-1168
control: 12 mice; P < 0.05 in 9 (Holm-adjusted: 7); smallest P 0.0000; stim minus catch speed change, median 3.96 mm/s (range -6.74 to 7.25); trials per mouse: stim 649-1201, catch 556-935
mouse opto sessions n_stim n_catch stim_frac stim_frac_min stim_frac_max stim_mean catch_mean diff U p
10 control 12 1076 777 0.5807 0.3978 0.6621 3.7501 -2.1131 5.8633 468576.0 0.0000
1527 control 12 649 556 0.5386 0.3556 0.6073 2.9069 0.4869 2.4200 186706.0 0.1484
1528 Opto-DA 12 1270 1067 0.5434 0.3833 0.5960 -7.4977 -5.7658 -1.7319 678536.0 0.4757
1529 Opto-DA 12 898 827 0.5206 0.4330 0.6522 2.6387 -2.5233 5.1620 403197.0 0.0010
1544 Opto-DA 12 1687 1168 0.5907 0.4699 0.6529 -6.5669 -3.0629 -3.5040 957771.0 0.8974
1546 Opto-DA 12 871 747 0.5383 0.4822 0.6750 -4.2382 -1.0746 -3.1636 316349.0 0.8308
1561 control 12 702 616 0.5326 0.3621 0.6406 -3.6173 3.1202 -6.7375 205601.0 0.9382
15822 Opto-DA 12 1305 766 0.6301 0.4118 0.7013 0.1909 -5.6262 5.8172 541177.0 0.0008
15823 Opto-DA 12 1244 864 0.5901 0.1111 0.6593 -1.8405 -4.5102 2.6697 560419.0 0.0470
15847 Opto-DA 12 1524 918 0.6241 0.2353 0.6914 -1.6329 -5.2212 3.5883 731411.0 0.0294
197 control 12 1090 810 0.5737 0.4118 0.6667 -3.8856 -8.8879 5.0023 476232.0 0.0016
200 control 12 1159 779 0.5977 0.4595 0.6912 -3.7577 -4.0504 0.2927 465004.0 0.1306
208 control 12 1201 817 0.5951 0.4314 0.6783 2.2388 -3.0937 5.3325 525373.0 0.0034
209 control 12 1140 838 0.5763 0.2000 0.6149 3.7894 0.8666 2.9228 520183.0 0.0004
211 control 12 1169 811 0.5904 0.3636 0.6673 -0.1319 -5.5585 5.4266 517825.0 0.0002
239 Opto-DA 12 1052 697 0.6011 0.3333 0.6588 2.1176 -1.6083 3.7259 399811.0 0.0007
240 Opto-DA 12 1355 960 0.5853 0.4186 0.6291 -9.6742 -8.9168 -0.7574 666535.0 0.1543
241 Opto-DA 12 1245 847 0.5951 0.3000 0.6504 -4.3966 -4.0377 -0.3589 536976.0 0.2368
242 Opto-DA 12 957 699 0.5779 0.4130 0.6364 -3.0819 -4.8628 1.7809 352916.0 0.0275
3440 Opto-DA 12 1012 616 0.6216 0.1429 0.7170 -2.1231 -1.3262 -0.7969 312766.0 0.4537
3441 Opto-DA 12 1116 665 0.6266 0.3333 0.6699 3.8230 3.6054 0.2175 386205.0 0.0747
3442 Opto-DA 12 1513 952 0.6138 0.4737 0.6410 0.7248 -2.3477 3.0725 770967.0 0.0016
410 Opto-DA 12 1405 891 0.6119 0.4444 0.6494 3.0158 -1.5056 4.5214 676229.0 0.0006
dls-ai32jr-1 Opto-DA 12 1496 909 0.6220 0.2759 0.6915 -13.5075 1.7746 -15.2821 585892.0 1.0000
dls-ai32jr-2 Opto-DA 11 1163 777 0.5995 0.3000 0.6913 -0.9196 0.7286 -1.6482 444109.0 0.7384
dls-ai32jr-3 Opto-DA 12 1557 964 0.6176 0.4857 0.7108 -6.5567 -3.5355 -3.0212 741832.0 0.6867
dls-ai32jr-4 Opto-DA 12 1321 602 0.6869 0.2000 0.7760 5.6441 -1.5161 7.1603 461589.0 0.0000
dls-ai32jr-5 Opto-DA 12 1436 947 0.6026 0.3913 0.6604 -12.5504 -0.6035 -11.9468 618160.0 0.9999
snc-dls-ctrl-6 control 12 1176 935 0.5568 0.4706 0.6257 2.3623 -4.4315 6.7938 632573.0 0.0000
snc-dls-ctrl-7 control 12 916 811 0.5301 0.4265 0.6000 -2.4837 -4.7890 2.3053 392006.0 0.0234
vta-nacc-ctrl-6 control 12 1095 813 0.5739 0.4634 0.7121 -5.6527 -6.8225 1.1698 467874.0 0.0279
vta-nacc-ctrl-7 control 12 1189 764 0.6088 0.4211 0.6733 -1.7129 -8.9647 7.2518 503324.0 0.0000
stimulated fraction of detections, 250 ms: per mouse median 0.593 (range 0.521 to 0.687); per session 0.111 to 0.776 (Methods: 0.75)
variant, catch >= 2 s from any pulse, Opto-DA: P < 0.05 in 8 of 20 mice; stim minus catch median 1.34 mm/s (range -16.08 to 7.59); catch trials per mouse 392-733 variant, catch >= 2 s from any pulse, control: P < 0.05 in 9 of 12 mice; stim minus catch median 3.85 mm/s (range -7.69 to 8.79); catch trials per mouse 417-688
In [21]:
k3m = per_mouse(C.pulse_s < 0.05, True, 3.0)
print(k3m.round(4).to_string(index=False))
s = k3m[k3m.opto == "Opto-DA"]
print(f"Opto-DA mice with 3-s pulsed sessions: {len(s)}; significant (P < 0.05): {(s.p < 0.05).sum()}")
# Variant: catch trials at least 4 s after the last stimulation onset; and the 2-s continuous sessions
k3m_b = per_mouse(C.pulse_s < 0.05, True, 3.0, catch_exclusion=4.0)
print("variant, catch >= 4 s after a stimulation:", k3m_b[["mouse", "opto", "n_catch", "diff", "p"]].round(4).to_dict("records"))
k3m_iso = per_mouse(C.pulse_s < 0.05, True, 3.0, isolate=2.0)
print("variant, catch >= 2 s from any pulse:", k3m_iso[["mouse", "opto", "n_stim", "n_catch", "diff", "p"]].round(4).to_dict("records"))
print("stimulated fraction of detections, 3 s pulsed: per mouse %.3f to %.3f; per session %.3f to %.3f (Methods: 0.75)"
% (k3m.stim_frac.min(), k3m.stim_frac.max(), k3m.stim_frac_min.min(), k3m.stim_frac_max.max()))
k2s = per_mouse(C.pulse_s == 2.0, False, 2.0)
print("2-s continuous stimulation:", k2s[["mouse", "opto", "n_stim", "n_catch", "diff", "p"]].round(4).to_dict("records"))
mouse opto sessions n_stim n_catch stim_frac stim_frac_min stim_frac_max stim_mean catch_mean diff U p
3440 Opto-DA 12 606 1159 0.3433 0.2956 0.5208 -5.2571 -8.0433 2.7862 364729.0 0.0913
3441 Opto-DA 12 473 885 0.3483 0.2759 0.5156 -0.3315 -4.9099 4.5784 222898.0 0.0242
3442 Opto-DA 12 485 1120 0.3022 0.2045 0.6000 0.7648 -4.5904 5.3552 285397.0 0.0528
snc-dls-ctrl-6 control 12 528 1052 0.3342 0.2764 0.5294 0.6401 -6.6825 7.3226 307550.0 0.0002
vta-nacc-ctrl-6 control 12 386 828 0.3180 0.2525 0.3776 -2.4424 -8.6218 6.1794 173427.0 0.0083
Opto-DA mice with 3-s pulsed sessions: 3; significant (P < 0.05): 1
variant, catch >= 4 s after a stimulation: [{'mouse': '3440', 'opto': 'Opto-DA', 'n_catch': 798, 'diff': 1.3291, 'p': 0.1495}, {'mouse': '3441', 'opto': 'Opto-DA', 'n_catch': 617, 'diff': 3.121, 'p': 0.096}, {'mouse': '3442', 'opto': 'Opto-DA', 'n_catch': 804, 'diff': 3.9417, 'p': 0.1433}, {'mouse': 'snc-dls-ctrl-6', 'opto': 'control', 'n_catch': 737, 'diff': 7.5441, 'p': 0.0004}, {'mouse': 'vta-nacc-ctrl-6', 'opto': 'control', 'n_catch': 626, 'diff': 5.2267, 'p': 0.0162}]
variant, catch >= 2 s from any pulse: [{'mouse': '3440', 'opto': 'Opto-DA', 'n_stim': 606, 'n_catch': 528, 'diff': 2.3788, 'p': 0.1035}, {'mouse': '3441', 'opto': 'Opto-DA', 'n_stim': 473, 'n_catch': 446, 'diff': 3.9169, 'p': 0.0642}, {'mouse': '3442', 'opto': 'Opto-DA', 'n_stim': 485, 'n_catch': 457, 'diff': 5.0331, 'p': 0.0817}, {'mouse': 'snc-dls-ctrl-6', 'opto': 'control', 'n_stim': 528, 'n_catch': 522, 'diff': 7.826, 'p': 0.0007}, {'mouse': 'vta-nacc-ctrl-6', 'opto': 'control', 'n_stim': 386, 'n_catch': 513, 'diff': 6.0543, 'p': 0.0102}]
stimulated fraction of detections, 3 s pulsed: per mouse 0.302 to 0.348; per session 0.205 to 0.600 (Methods: 0.75)
2-s continuous stimulation: [{'mouse': '3155', 'opto': 'control', 'n_stim': 820, 'n_catch': 900, 'diff': 1.8017, 'p': 0.1979}, {'mouse': '3440', 'opto': 'Opto-DA', 'n_stim': 954, 'n_catch': 1020, 'diff': 1.9131, 'p': 0.0907}, {'mouse': '3441', 'opto': 'Opto-DA', 'n_stim': 831, 'n_catch': 858, 'diff': 12.4525, 'p': 0.0}, {'mouse': '3442', 'opto': 'Opto-DA', 'n_stim': 1103, 'n_catch': 1147, 'diff': 4.7368, 'p': 0.001}, {'mouse': 'snc-dls-ctrl-6', 'opto': 'control', 'n_stim': 877, 'n_catch': 1066, 'diff': 9.5721, 'p': 0.0}, {'mouse': 'vta-nacc-ctrl-6', 'opto': 'control', 'n_stim': 568, 'n_catch': 681, 'diff': 10.8034, 'p': 0.0}]
In [22]:
# How the substitute catch trials behave in controls, and Opto-DA against controls
for name, k in [("250 ms", k3l), ("3 s pulsed", k3m), ("250 ms, catch >= 2 s from any pulse", k3l_iso),
("3 s pulsed, catch >= 2 s from any pulse", k3m_iso)]:
a, b = k[k.opto == "Opto-DA"], k[k.opto == "control"]
print(f"{name}: mice with stimulation above catch at P < 0.05: Opto-DA {(a.p < 0.05).sum()} of {len(a)}, "
f"control {(b.p < 0.05).sum()} of {len(b)}")
print(f" stimulation minus catch (mm/s): Opto-DA median {a['diff'].median():.2f}, control median {b['diff'].median():.2f}; "
f"speed change after stimulation: Opto-DA median {a.stim_mean.median():.2f}, control median {b.stim_mean.median():.2f}")
if len(a) >= 3 and len(b) >= 3:
u = stats.mannwhitneyu(a["diff"], b["diff"], alternative="two-sided")
print(f" Opto-DA against control on stimulation minus catch: U = {u.statistic:.0f}, P = {u.pvalue:.3f}")
250 ms: mice with stimulation above catch at P < 0.05: Opto-DA 9 of 20, control 9 of 12 stimulation minus catch (mm/s): Opto-DA median -0.07, control median 3.96; speed change after stimulation: Opto-DA median -1.98, control median -0.92 Opto-DA against control on stimulation minus catch: U = 68, P = 0.045 3 s pulsed: mice with stimulation above catch at P < 0.05: Opto-DA 1 of 3, control 2 of 2 stimulation minus catch (mm/s): Opto-DA median 4.58, control median 6.75; speed change after stimulation: Opto-DA median -0.33, control median -0.90 250 ms, catch >= 2 s from any pulse: mice with stimulation above catch at P < 0.05: Opto-DA 8 of 20, control 9 of 12 stimulation minus catch (mm/s): Opto-DA median 1.34, control median 3.85; speed change after stimulation: Opto-DA median -1.98, control median -0.92 Opto-DA against control on stimulation minus catch: U = 65, P = 0.034 3 s pulsed, catch >= 2 s from any pulse: mice with stimulation above catch at P < 0.05: Opto-DA 0 of 3, control 2 of 2 stimulation minus catch (mm/s): Opto-DA median 3.92, control median 6.94; speed change after stimulation: Opto-DA median -0.33, control median -0.90
Result: fig3l-no-kinematic-effect and fig3m-long-stim-velocity blocked (metadata insufficient). The files mark stimulation onsets but not catch trials (target detections on which stimulation was withheld). The criterion took catch trials as online target onsets with no pulse in the next 0.5 s. This is not the paper's catch set: the Methods give stimulation on 75 % of detections, while here stimulations are a median 0.593 of detections per mouse in the 250-ms sessions (0.521 to 0.687 over mice, 0.111 to 0.776 over sessions) and 0.302 to 0.348 per mouse in the 3-s sessions (0.205 to 0.600 over sessions). The substitute does not isolate the effect of light: in no-opsin controls, stimulation trials differ from these catch trials at P < 0.05 in 9 of 12 mice with 250-ms pulses and in 2 of 2 with 3-s pulses, with differences at least as large as in Opto-DA mice (median stimulation minus catch 3.96 mm/s for controls and -0.07 mm/s for Opto-DA at 250 ms; 6.75 against 4.58 mm/s at 3 s). Keeping only catch trials with no pulse within 2 s leaves the controls significant (9 of 12 and 2 of 2). Taken literally, the criteria would give fig3l not reproduced (9 of 20 Opto-DA mice significant) and fig3m partially reproduced (1 of 3); with catch trials 2 s from any pulse, 8 of 20 and 0 of 3. Since the controls fail, neither claim can be judged with these data.
Summary¶
| Claim | Outcome |
|---|---|
| fig1e-velocity-timescales | partially reproduced (negative at short timescales for both measures; reversal for transient rate only) |
| fig1h-syllable-dlight-differs | partially reproduced (P = 6.0e-37, n = 518; H 266.40 against 209.29) |
| fig1l-identity-not-decodable | reproduced |
| fig2b-dlight-usage-across-syllables | reproduced |
| fig2i-dlight-entropy-across-syllables | reproduced |
| fig2e-dlight-predicts-future-usage | partially reproduced (significant at 6 of 10 n up to 190) |
| ed7-dms-fewer-transients | reproduced |
| ed7-dms-no-usage-prediction | not reproduced (DMS dLight predicts future use by the same test) |
| fig3e-optoda-learning | reproduced (U = 198, P = 0.0013) |
| fig3f-cumulative-counts | reproduced (U = 182, P = 0.0083) |
| fig3e-learner-count | not reproduced (7 learners, not 9) |
| fig3h-stim-post-correlation | reproduced (r = 0.891) |
| fig3l-no-kinematic-effect | blocked: metadata insufficient (catch trials not marked; substitute fails in controls) |
| fig3m-long-stim-velocity | blocked: metadata insufficient (catch trials not marked; substitute fails in controls) |
| fig2a-food-reward-transients | blocked: data missing |
| fig3i-optoda-entropy, fig4f-catch-dlight-predicts-learning | blocked: metadata insufficient (catch trials not marked); not attempted |
| fig1f, fig2k, ed10d, fig4b, fig4c, fig4h, fig5 | blocked: not attempted here |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's:
- dLight is re-referenced and z-scored here from the stored dF/F series, following the Methods; the paper's analyses used a pre-computed z-scored trace that is not in the files.
- The dandiset holds fewer photometry sessions than the paper's code used: 177 DLS sessions without stimulation are used here (224 files, 210 with syllables, 193 after leaving out habituation sessions, 177 passing the quality criteria), where the Methods give 379 experiments.
- Habituation sessions (17 DLS, 1 DMS) are left out of the photometry analyses, as they are of the closed-loop analyses; the Methods do not say whether the paper's experiments include them.
- Fig. 2e: five equal time segments stand in for the paper's five behavioural segments, and r is compared with shuffles directly instead of being z-scored to them.
- Fig. 2b and 2i: 20 repeats of fivefold cross-validation and 200 shuffles, where the paper used 100 repeats and 1,000 shuffles; a plain Shannon entropy where the code used a Dirichlet estimate.
- Fig. 1l: a 100-tree sklearn random forest (at least 5 instances per leaf) on 1,000 instances per syllable and one label shuffle, where the paper used cuML with 1,000 trees on all instances and 100 shuffles.
- Fig. 3: the baseline and post days are identified from session dates and stimulation, as the files give no session number or session type.
- Fig. 3l and 3m: catch trials are a substitute (target onsets with no pulse in the next 0.5 s), which covers more detections than the paper's 25 % and fails in no-opsin controls.
What was found about the dandiset:
- It holds the 14 DLS dLight mice, 8 DMS dLight mice, 20 Opto-DA mice and 12 no-opsin controls of the paper, and the closed-loop sessions reproduce the paper's Fig. 3e and 3f statistics closely.
- Catch trials of the closed-loop experiments are not marked, nor are the opsin and genotype of each mouse; Opto-DA and control mice can be told apart only by the stimulation site.
- Of the dLight+ChrimsonR mice, only the 250-ms sessions are present, not the 2-s and pulsed sessions the Methods describe; DLS photometry files without stimulation number 224 (177 usable), where the Methods give 379 experiments.
- 14 files hold empty syllable series; 28 photometry and 547 closed-loop files store their frames out of time order (in all but two, one backward step, at least 0.95 of the session length in 558 of 575), and two closed-loop files of mouse 3155 store every frame twice; 5 files with stimulation are described as recording sessions; every closed-loop file states a 635-nm stimulator, including those of Ai32 mice stimulated with 470-nm light; the dLight-only mice carry the subject description of the DAT-IRES-Cre and Ai32 mice.
- No food-reward sessions or consumption times are in the dandiset.
In [23]:
print(f"notebook run time: {(time.time() - T0) / 60:.1f} min")
notebook run time: 12.0 min