Reproducing Jeong, Taylor, Floeder et al. (2022) from dandiset 000351¶
Paper. Jeong H, Taylor A, Floeder JR, Lohmann M, Mihalas S, Wu B, Zhou M, Burke DA,
Namboodiri VMK (2022). Mesolimbic dopamine release conveys causal associations. Science 378,
eabq6740. doi:10.1126/science.abq6740. The claims file
(claims/000351/10.1126-science.abq6740.yaml) was written from the accepted author manuscript in
PubMed Central (PMC9910357), its figure legends, and its supplement (Materials and Methods,
Supplementary Notes, Table S1). Its criteria cite the Methods and the paper's code,
github.com/namboodirilab/ANCCR at commit cff44144e60a4313a30ffbb4bb8a1f3dd398b525.
Dandiset. DANDI:000351, version 0.260828.1905. The paper cites the dandiset without a version (its code reads the draft). This is the only published version, so it is pinned. The draft has the same 428 assets and size; it is not used, since it can change.
What the dandiset holds. Photometry (dF/F of dLight1.3b in the nucleus accumbens core) and the event log (licks, rewards, cues, laser) of the 8 mice of Experiments 1 to 6, the 14 mice of Experiment 7, and the event logs (no photometry) of the 6 mice of Fig S8L. The simulations of the paper (Figs 2, 3C and F, 5B, 6H, S5, S6, S13) use no recorded data.
Outside inputs. The day lists (startday, endday, learnedday) that the paper's Matlab scripts
use to select sessions are copied into the helper jeong.py from the ANCCR repository at the
commit above, and the notebook reads the scripts at that commit to check the copies; they name sessions by the Day number in the dandiset's file names
(claim session-day-link). The spreadsheet that the dandiset's description names ("Subject and
session information", a Google spreadsheet) is not used: it has no fixed version to pin.
| Claim | Tested here |
|---|---|
| session-day-link, methods-session-counts | yes |
| exp1-excluded-rewards, fig3e, fig3g, figs8k, figs8b, figs8d, figs8g, figs8h, figs8l | yes |
| figs8e | tested; blocked, methods_underspecified (the readings differ) |
| fig4ef, figs10c, fig4hi (two), figs10b, fig4kl, fig4jl, fig4p, fig6b | yes |
| fig5d (two), figs11de | yes |
| fig6d (two), fig6g, fig6i (three), figs12f | yes |
| fig2, fig3cf, fig5b, fig6h, figs6, figs5, figs13 (simulations) | blocked, not_runnable |
import time
import warnings
import numpy as np
import pandas as pd
import pingouin as pg
import statsmodels.api as sm
from scipy import stats
import jeong
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 400)
started = time.time()
assets = jeong.list_assets()
photometry_assets = [a for a in assets if not a["subject"].startswith("DB-")]
lick_only_assets = [a for a in assets if a["subject"].startswith("DB-")]
sessions = jeong.load_many(photometry_assets, dff=True)
sessions += jeong.load_many(lick_only_assets, dff=False)
print(f"{len(assets)} assets in dandiset {jeong.DANDISET_ID} version {jeong.DANDISET_VERSION}; "
f"loaded {len(sessions)} sessions in {time.time() - started:.0f} s")
def subject_sessions(subject, task=None):
return sorted([s for s in sessions if s["subject"] == subject
and (task is None or s["description"] == task)], key=lambda s: s["day"])
def ols_slope(x, y):
fit = sm.OLS(np.asarray(y, float), sm.add_constant(np.asarray(x, float))).fit()
return fit.params[1], fit.tvalues[1], fit.pvalues[1], int(fit.nobs)
def show_t(label, res):
print(f"{label}: t = {res.statistic:.3f}, p = {res.pvalue:.4g}, df = {res.df:.0f}")
428 assets in dandiset 000351 version 0.260828.1905; loaded 428 sessions in 30 s
Data check¶
Each check of the metadata against the values is stated under its output.
info = pd.DataFrame([{
"subject": s["subject"], "day": s["day"], "session": s["session"], "task": s["description"],
"genotype": s["genotype"], "sex": s["sex"], "subject_id": s["subject_id"], "dff": s["has_dff"],
"event_unit": s["event_unit"], "dff_unit": s.get("dff_unit"), "size_MB": s["size"] / 1e6,
"last_event_s": s["time"].max()} for s in sessions])
print("Sessions by subject and task (session_description):")
print(pd.crosstab(info.subject, info.task).to_string())
print("\nGenotype, sex and subject_id per subject:")
print(info.groupby("subject")[["genotype", "sex", "subject_id"]].first().to_string())
print("\nSubjects:", info.subject.nunique(), "; files with dF/F:", int(info.dff.sum()),
"; files without:", int((~info.dff).sum()), "(subjects", sorted(info[~info.dff].subject.unique()), ")")
print("Stated units: event times", info.event_unit.unique().tolist(), "; dF/F", info.dff_unit.dropna().unique().tolist())
Sessions by subject and task (session_description):
task DB_longITI_C1_F1 DB_longITI_C1_F2 DB_longITI_C1_M1 DB_longITI_C1_M2 DB_longITI_C1_M3 DB_longITI_C1_M4 Pavlovian Poisson RandomRewards
subject
DB-longITI-C1-F1 6 0 0 0 0 0 0 0 0
DB-longITI-C1-F2 0 6 0 0 0 0 0 0 0
DB-longITI-C1-M1 0 0 6 0 0 0 0 0 0
DB-longITI-C1-M2 0 0 0 6 0 0 0 0 0
DB-longITI-C1-M3 0 0 0 0 6 0 0 0 0
DB-longITI-C1-M4 0 0 0 0 0 6 0 0 0
HJ-FP-F1 0 0 0 0 0 0 26 4 5
HJ-FP-F2 0 0 0 0 0 0 25 5 8
HJ-FP-M2 0 0 0 0 0 0 33 4 11
HJ-FP-M3 0 0 0 0 0 0 30 2 7
HJ-FP-M4 0 0 0 0 0 0 37 3 8
HJ-FP-M6 0 0 0 0 0 0 30 3 6
HJ-FP-M7 0 0 0 0 0 0 29 4 4
HJ-FP-M8 0 0 0 0 0 0 0 0 7
HJ-FP-WT-stGtACR-F1 0 0 0 0 0 0 5 0 1
HJ-FP-WT-stGtACR-F2 0 0 0 0 0 0 5 0 1
HJ-FP-WT-stGtACR-F3 0 0 0 0 0 0 6 0 1
HJ-FP-WT-stGtACR-M1 0 0 0 0 0 0 5 0 1
HJ-FP-WT-stGtACR-M2 0 0 0 0 0 0 5 0 1
HJ-FP-WT-stGtACR-M3 0 0 0 0 0 0 5 0 1
HJ-FP-datHT-stGtACR-M2 0 0 0 0 0 0 5 0 1
HJ-FP-datHT-stGtACR-M3 0 0 0 0 0 0 6 0 1
HJ-FP-datHT-stGtACR-M4 0 0 0 0 0 0 6 0 1
HJ-FP-datHT-stGtACR-M5 0 0 0 0 0 0 8 0 1
HJ-FP-datHT-stGtACR-M6 0 0 0 0 0 0 7 0 1
HJ-FP-datHT-stGtACR-M7 0 0 0 0 0 0 12 0 1
HJ-FP-datHT-stGtACR-M8 0 0 0 0 0 0 5 0 1
HJ-FP-datWT-stGtACR-F1 0 0 0 0 0 0 7 0 1
Genotype, sex and subject_id per subject:
genotype sex subject_id
subject
DB-longITI-C1-F1 WT F DB_longITI_C1_F1
DB-longITI-C1-F2 WT F DB_longITI_C1_F2
DB-longITI-C1-M1 WT M DB_longITI_C1_M1
DB-longITI-C1-M2 WT M DB_longITI_C1_M2
DB-longITI-C1-M3 WT M DB_longITI_C1_M3
DB-longITI-C1-M4 WT M DB_longITI_C1_M4
HJ-FP-F1 WT F HJ_FP_F1
HJ-FP-F2 WT F HJ_FP_F2
HJ-FP-M2 WT M HJ_FP_M2
HJ-FP-M3 WT M HJ_FP_M3
HJ-FP-M4 WT M HJ_FP_M4
HJ-FP-M6 WT M HJ_FP_M6
HJ-FP-M7 WT M HJ_FP_M7
HJ-FP-M8 WT M HJ_FP_M8
HJ-FP-WT-stGtACR-F1 wt F HJ_FP_WT_stGtACR_F1
HJ-FP-WT-stGtACR-F2 wt F HJ_FP_WT_stGtACR_F2
HJ-FP-WT-stGtACR-F3 wt F HJ_FP_WT_stGtACR_F3
HJ-FP-WT-stGtACR-M1 wt M HJ_FP_WT_stGtACR_M1
HJ-FP-WT-stGtACR-M2 wt M HJ_FP_WT_stGtACR_M2
HJ-FP-WT-stGtACR-M3 wt M HJ_FP_WT_stGtACR_M3
HJ-FP-datHT-stGtACR-M2 DAT-ht M HJ_FP_datHT_stGtACR_M2
HJ-FP-datHT-stGtACR-M3 DAT-ht M HJ_FP_datHT_stGtACR_M3
HJ-FP-datHT-stGtACR-M4 DAT-ht M HJ_FP_datHT_stGtACR_M4
HJ-FP-datHT-stGtACR-M5 DAT-ht M HJ_FP_datHT_stGtACR_M5
HJ-FP-datHT-stGtACR-M6 DAT-ht M HJ_FP_datHT_stGtACR_M6
HJ-FP-datHT-stGtACR-M7 wt M HJ_FP_datHT_stGtACR_M7
HJ-FP-datHT-stGtACR-M8 DAT-ht M HJ_FP_datHT_stGtACR_M8
HJ-FP-datWT-stGtACR-F1 DAT-wt F HJ_FP_datWT_stGtACR_F1
Subjects: 28 ; files with dF/F: 392 ; files without: 36 (subjects ['DB-longITI-C1-F1', 'DB-longITI-C1-F2', 'DB-longITI-C1-M1', 'DB-longITI-C1-M2', 'DB-longITI-C1-M3', 'DB-longITI-C1-M4'] )
Stated units: event times ['milliseconds'] ; dF/F ['%']
The archive counts 28 subjects; the files hold 28 (8 of Experiments 1 to 6, 14 of Experiment 7,
6 of Fig S8L), with the session tasks the paper describes. The Fig S8L files hold no photometry,
which that panel does not use. The Experiment 7 subject HJ-FP-datHT-stGtACR-M7 is named as a
DAT-Cre heterozygote but its genotype field says wt; the code's lists treat it as DAT-Cre
(dathetlist), and so does this analysis. The sessions of datHT-M6 and datHT-M7 lack the task in
their file names, but their session_description gives it.
# Event times: the files say milliseconds. The random reward task has a mean inter-reward
# interval of 12 s (Methods); the long-ITI group has 6, 12 and 30 s by name.
rows = []
for s in sessions:
rw = jeong.events(s, jeong.BGD_REWARD)
if s["description"] in ("RandomRewards",) or s["subject"].startswith("DB-"):
if len(rw) > 10:
label = s["session"].split("-")[1] if s["subject"].startswith("DB-") else "12 (Methods)"
rows.append({"group": "Fig S8L" if s["subject"].startswith("DB-") else "random rewards",
"nominal_IRI": label, "mean_IRI_in_file_units": np.mean(np.diff(rw))})
iri = pd.DataFrame(rows)
print(iri.groupby(["group", "nominal_IRI"]).mean_IRI_in_file_units.describe()[["count", "mean", "min", "max"]].round(2).to_string())
print("\nSession length (last event, file units): ", info.last_event_s.describe()[["min", "50%", "max"]].round(0).to_dict())
count mean min max
group nominal_IRI
Fig S8L 12 12.0 12.11 10.13 14.02
30 12.0 30.24 26.50 36.67
6 12.0 5.98 5.16 6.63
random rewards 12 (Methods) 70.0 11.11 5.23 14.26
Session length (last event, file units): {'min': 481.0, '50%': 2916.0, 'max': 4810.0}
The mean interval between random rewards is about 12 in the file's units, and 6, 12 and 30 in the sessions named for those intervals, and sessions last 481 to 4810 units: event times are in seconds, not the milliseconds the files state. The analysis treats them as seconds (the paper's Python scripts do the same). This is recorded as a finding.
# Photometry sampling and coverage.
rows = []
for s in sessions:
if not s["has_dff"]:
continue
t = s["t"]
rows.append({"subject": s["subject"], "session": s["session"],
"experiment": "7" if "stGtACR" in s["subject"] else "1-6",
"samples": len(t), "dt_ms": round(float(np.nanmedian(np.diff(t))) * 1000, 2),
"nan_time": int(np.isnan(t).sum()), "nan_dff": int(np.isnan(s["dff"]).sum()),
"dff_start": float(np.nanmin(t)), "dff_end": float(np.nanmax(t)),
"first_event": float(s["time"].min()), "last_event": float(s["time"].max()),
"dff_p1": float(np.nanpercentile(s["dff"], 1)), "dff_p99": float(np.nanpercentile(s["dff"], 99)),
"raw_samples": s.get("raw_n"), "raw_t_last": s.get("raw_t_last")})
ph = pd.DataFrame(rows)
print("Median dF/F sample interval (ms) by experiment:")
print(ph.groupby(["experiment", "dt_ms"]).size().to_string())
print("\nNaN timestamps and dF/F values per session:")
print(ph.groupby("experiment")[["nan_time", "nan_dff"]].describe().round(0).to_string())
print("\ndF/F 1st and 99th percentiles over sessions (stated unit %):",
ph.dff_p1.median().round(2), ph.dff_p99.median().round(2))
cover = ph[(ph.dff_start > ph.first_event + 1) | (ph.dff_end < ph.last_event - 1)]
print(f"\nSessions whose dF/F starts after the first event or ends before the last one (by more than 1 s): {len(cover)}")
print(cover[["subject", "session", "dff_start", "dff_end", "first_event", "last_event"]].round(1).to_string(index=False))
raw7 = ph[ph.experiment == "7"]
rate = raw7.raw_samples / raw7.raw_t_last
big = raw7[rate > 1000]
print(f"\nExperiment 7 raw photometry (acquisition/photometry): {len(raw7)} files; {len(big)} have more than 1000 samples "
f"per stated second (median {rate[rate > 1000].median():.0f}), their last raw timestamp {big.raw_t_last.min():.2f} to "
f"{big.raw_t_last.max():.2f} against dF/F ending at {big.dff_end.min():.0f} to {big.dff_end.max():.0f} s; "
f"the other {len(raw7) - len(big)} run at {rate[rate <= 1000].median():.0f} samples per second")
Median dF/F sample interval (ms) by experiment:
experiment dt_ms
1-6 8.0 247
9.0 44
7 10.0 101
NaN timestamps and dF/F values per session:
nan_time nan_dff
count mean std min 25% 50% 75% max count mean std min 25% 50% 75% max
experiment
1-6 291.0 17.0 7.0 0.0 20.0 20.0 20.0 33.0 291.0 349.0 312.0 300.0 310.0 314.0 323.0 4953.0
7 101.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 101.0 13.0 11.0 0.0 0.0 20.0 20.0 54.0
dF/F 1st and 99th percentiles over sessions (stated unit %): -1.43 3.2
Sessions whose dF/F starts after the first event or ends before the last one (by more than 1 s): 13
subject session dff_start dff_end first_event last_event
HJ-FP-F1 Day26 -1.8 1452.1 0.6 1453.5
HJ-FP-F2 Day13 0.1 1952.0 1.2 2045.8
HJ-FP-F2 Day17-learned 0.1 2899.7 1.2 2919.8
HJ-FP-F2 Day25 0.1 2319.6 37.0 2337.6
HJ-FP-M2 Day8 0.1 1094.2 3.1 1154.8
HJ-FP-M2 Day33-8s-bgd6s -7.6 1882.6 6.8 1930.1
HJ-FP-M4 Day3 0.1 1094.2 3.1 1154.8
HJ-FP-M4 Day71 0.1 2368.8 21.6 2370.9
HJ-FP-M6 Day22 -2.1 1376.8 1.2 1563.5
HJ-FP-M6 Day48 -0.0 2363.7 13.0 2407.1
HJ-FP-M7 Day7 0.1 1641.7 4.9 1668.9
HJ-FP-M7 Day43-extinction -0.0 4012.3 6.2 4187.0
HJ-FP-WT-stGtACR-M2 Day6-Pavlovian 1117.3 3239.9 1.6 3208.5
Experiment 7 raw photometry (acquisition/photometry): 101 files; 65 have more than 1000 samples per stated second (median 12047869), their last raw timestamp 0.43 to 4.18 against dF/F ending at 596 to 4133 s; the other 36 run at 130 samples per second
dF/F is sampled every 8 or 9 ms in Experiments 1 to 6 (120 Hz Doric demodulated output in the
Methods) and every 10 ms in Experiment 7, where dF/F starts before the event log (by about 15 s
in most sessions) and has no NaN timestamps. In Experiments 1 to 6 the first 20 or so timestamps
and about 300 dF/F values per session are NaN, the samples before the photometry starts; these are
left out by the measures (as nanmean and min do in the code). The sessions listed lose the end
or the start of their photometry (WT-M2 Day6 has none for its first 1117 s); the measures treat trials
there as the paper's code does. aucsignal2event.m (Experiment 1) sums with sum, so a reward
with a NaN sample in its window gets no value and is dropped. aucsignal.m (Experiments 2 to 6)
sums with nansum, so a trial outside the coverage gets an AUC of 0, and one at its edge, or near
an isolated NaN sample spread by the smoothing, the sum of its other samples. These trials are
counted under each claim, with the result when they are left out. In Experiment 7,
calculate_auc sums the samples present, and a window with none gets no value. In 65 of the 101 Experiment 7 files
the raw photometry holds about 12 million samples per stated second and ends at a timestamp of
0.43 to 4.18, where dF/F runs to 596 to 4133 s: the timestamps of the raw series are in thousands
of seconds (12 kHz, the Doric detection rate of the Methods). The raw series is not used here; this is recorded as a finding.
# dF/F scale: the files say %. The Methods define dF/F as a fraction, (470 - fitted 405) / fitted
# 405. Compare the stored dF/F with that fraction computed from the raw channels of one session
# per photometry system (the paper's baseline correction is not repeated, so only the scale and
# the correlation are checked).
import h5py, remfile
check = [a for a in assets if a["path"].endswith("sub-HJ-FP-M2_ses-Day1.nwb")
or a["path"].endswith("sub-HJ-FP-WT-stGtACR-F3_ses-Day2-Pavlovian.nwb")]
for a in check:
with h5py.File(remfile.File(a["url"]), "r") as f:
raw = f["acquisition/photometry/data"][:]
rt = f["acquisition/photometry/timestamps"][:]
d = f["processing/photometry/dff/data"][:]
dt = f["processing/photometry/dff/timestamps"][:]
good = np.isfinite(rt) & np.isfinite(raw).all(1)
best = None
for sig, ref in [(0, 1), (1, 0)]:
c = np.polyfit(raw[good, ref], raw[good, sig], 1)
fit = np.polyval(c, raw[good, ref])
frac = (raw[good, sig] - fit) / fit
# dF/F timestamps are on the behaviour clock; the raw series starts at its own 0. Map the
# raw samples onto dF/F time by their position in the recording.
td = dt[np.isfinite(dt)]
pos = (rt[good] - rt[good][0]) / (rt[good][-1] - rt[good][0])
stored = np.interp(td[0] + pos * (td[-1] - td[0]), td, d[np.isfinite(dt)])
ok = np.isfinite(stored)
r = np.corrcoef(frac[ok], stored[ok])[0, 1]
slope = np.polyfit(frac[ok], stored[ok], 1)[0]
if best is None or r > best[1]:
best = (sig, r, slope)
print(f"{a['path']}: channel {best[0]} as the signal, r = {best[1]:.2f} with stored dF/F, "
f"stored dF/F = {best[2]:.1f} x the fraction")
sub-HJ-FP-M2/sub-HJ-FP-M2_ses-Day1.nwb: channel 1 as the signal, r = 0.96 with stored dF/F, stored dF/F = 87.8 x the fraction
sub-HJ-FP-WT-stGtACR-F3/sub-HJ-FP-WT-stGtACR-F3_ses-Day2-Pavlovian.nwb: channel 1 as the signal, r = 1.00 with stored dF/F, stored dF/F = 98.8 x the fraction
The stored dF/F follows the fraction computed from the raw channels and is about 100 times it, so
it is in percent, as stated. The paper's Matlab scripts give their bounds for fitting dopamine
responses (claim fig4hi-delay-change-dopamine-unchanged) in units of a fraction times
milliseconds (aucsignal.m multiplies by the bin in ms); the helper divides dF/F by 100 for the
Matlab-derived measures, so that their values are on the scale the Methods define (a fraction).
Whether the code's bounds were set for this scale is not known: they bind on it as well as on
percent (fig4hi-delay-change-dopamine-unchanged), and the analysis assumes the paper's fit was
bounded as here. Ratios, correlations and t-tests do not depend on it. For the
pyPhotometry session the raw samples are placed on the dF/F clock by their position in the
recording, since the dF/F timestamps start about 15 s before the raw series' zero.
# Recordings deposited more than once: identical event logs or dF/F traces across files.
import hashlib
def digest(a):
return hashlib.md5(np.ascontiguousarray(a).tobytes()).hexdigest()
by_events, by_dff, by_licks, by_schedule = {}, {}, {}, {}
for s in sessions:
by_events.setdefault(digest(s["time"]), []).append(s["path"])
by_licks.setdefault(digest(jeong.events(s, jeong.LICK)), []).append(s["path"])
by_schedule.setdefault(digest(np.round(s["time"][np.isin(s["code"], [7, 10, 15, 16])], 2)), []).append(s["path"])
if s["has_dff"]:
by_dff.setdefault(digest(np.nan_to_num(s["dff"])), []).append(s["path"])
for name, d in [("event logs", by_events), ("lick times", by_licks), ("dF/F traces", by_dff),
("reward and cue schedules", by_schedule)]:
dup = [v for v in d.values() if len(v) > 1]
print(f"Files sharing identical {name}: {len(dup)} groups")
for v in dup:
print(" ", " = ".join(p.split("/")[1] for p in v))
Files sharing identical event logs: 0 groups
Files sharing identical lick times: 0 groups
Files sharing identical dF/F traces: 0 groups
Files sharing identical reward and cue schedules: 8 groups
sub-HJ-FP-F1_ses-Day24.nwb = sub-HJ-FP-M3_ses-Day32.nwb
sub-HJ-FP-F1_ses-Day42-extinction.nwb = sub-HJ-FP-M7_ses-Day42.nwb
sub-HJ-FP-F2_ses-Day21.nwb = sub-HJ-FP-M7_ses-Day15.nwb
sub-HJ-FP-M2_ses-Day1.nwb = sub-HJ-FP-M4_ses-Day6.nwb
sub-HJ-FP-M4_ses-Day19.nwb = sub-HJ-FP-M4_ses-Day24.nwb
sub-HJ-FP-WT-stGtACR-F3_ses-Day2-Pavlovian.nwb = sub-HJ-FP-WT-stGtACR-F3_ses-Day7-Pavlovian.nwb
sub-HJ-FP-datHT-stGtACR-M2_ses-Day1-RandomRewards.nwb = sub-HJ-FP-datHT-stGtACR-M8_ses-Day1-RandomRewards.nwb
sub-HJ-FP-datHT-stGtACR-M2_ses-Day7-Pavlovian.nwb = sub-HJ-FP-datWT-stGtACR-F1_ses-Day4-Pavlovian.nwb
No two files share an event log, lick times or a dF/F trace, so no recording is deposited twice. Eight pairs of sessions share the same reward and cue times: the task schedules were reused (across animals, or across days of one animal), while the licks and dF/F differ. Each file is counted as one session.
# Shared clock of dF/F and the event log. Per session, the response to unpredicted events (the
# first lick after a random reward, cue onsets) or the contrast of CS+ against CS- onsets
# (Experiments 2 to 5), with the event times shifted by lags from -2 to 2 s, and by 10 s. If the
# two streams share a clock, the response peaks near lag 0 (within the response latency).
lags = np.round(np.arange(-2, 2.01, 0.25), 2) # the rewards come 3 s after cues; larger lags reach them
csplus_ids = {}
for m in jeong.PAV_MICE:
first = [s for s in subject_sessions(m, "Pavlovian") if s["day"] == jeong.CUMSUM_START[m][0]][0]
csplus_ids[m] = jeong.cs_plus_identity(first)
rows = []
for s in sessions:
if not s["has_dff"]:
continue
prof, n, ref = jeong.response_profile(s, lags, csplus_ids.get(s["subject"]))
group = "7" if "stGtACR" in s["subject"] else {"RandomRewards": "1", "Pavlovian": "2-5", "Poisson": "6"}[s["description"]]
rows.append({"experiment": group, "reference": ref, "subject": s["subject"], "session": s["session"], "n": n,
"peak_lag": lags[int(np.nanargmax(prof))], "at_0": prof[lags == 0][0],
"shift_10s": jeong.response_profile(s, [10.0], csplus_ids.get(s["subject"]))[0][0]})
clock = pd.DataFrame(rows)
print("Lag (s) of the largest response, per session:")
print(clock.groupby(["experiment", "reference"]).peak_lag.describe()[["count", "min", "25%", "50%", "75%", "max"]].to_string())
print("\nSessions with a positive response at lag 0, and with the events shifted by 10 s:")
print(clock.groupby(["experiment", "reference"]).apply(lambda g: pd.Series({
"sessions": len(g), "positive_at_0": int((g.at_0 > 0).sum()), "positive_shift_10s": int((g.shift_10s > 0).sum()),
"median_at_0": g.at_0.median(), "median_shift_10s": g.shift_10s.median()})).round(3).to_string())
far = clock[(clock.peak_lag < -1) | (clock.peak_lag > 1)]
print(f"\nSessions whose peak lag is beyond 1 s from 0: {len(far)}")
print(far[["subject", "session", "reference", "n", "peak_lag", "at_0"]].round(3).to_string(index=False))
Lag (s) of the largest response, per session:
count min 25% 50% 75% max
experiment reference
1 first lick after reward 56.0 -1.00 -0.3125 -0.25 0.000 2.0
2-5 CS+ minus CS- onset 210.0 -2.00 0.0000 0.00 0.000 2.0
6 cue onset 25.0 -2.00 0.0000 0.00 0.000 0.0
7 CS1 onset 87.0 -2.00 0.0000 0.00 0.125 2.0
first lick after reward 14.0 -0.25 0.0000 0.00 0.000 0.0
Sessions with a positive response at lag 0, and with the events shifted by 10 s:
sessions positive_at_0 positive_shift_10s median_at_0 median_shift_10s
experiment reference
1 first lick after reward 56.0 50.0 25.0 1.004 -0.010
2-5 CS+ minus CS- onset 210.0 198.0 62.0 1.528 -0.131
6 cue onset 25.0 24.0 14.0 1.098 0.011
7 CS1 onset 87.0 79.0 29.0 0.737 -0.019
first lick after reward 14.0 14.0 5.0 1.965 -0.008
Sessions whose peak lag is beyond 1 s from 0: 27
subject session reference n peak_lag at_0
HJ-FP-F1 Day26 CS+ minus CS- onset 17 -1.75 -0.087
HJ-FP-F2 Day50-Poisson-cue cue onset 100 -2.00 -0.217
HJ-FP-M2 Day1 first lick after reward 97 2.00 -0.020
HJ-FP-M2 Day37 CS+ minus CS- onset 20 2.00 0.113
HJ-FP-M2 Day38 CS+ minus CS- onset 20 1.75 -0.088
HJ-FP-M2 Day39 CS+ minus CS- onset 20 -1.25 -0.788
HJ-FP-M2 Day40 CS+ minus CS- onset 20 1.75 -0.268
HJ-FP-M2 Day41 CS+ minus CS- onset 20 1.50 -0.462
HJ-FP-M2 Day42 CS+ minus CS- onset 20 2.00 -0.982
HJ-FP-M2 Day43 CS+ minus CS- onset 20 1.75 -0.290
HJ-FP-M3 Day37 CS+ minus CS- onset 20 -2.00 0.087
HJ-FP-M3 Day38 CS+ minus CS- onset 20 -1.75 0.248
HJ-FP-M4 Day44 CS+ minus CS- onset 20 1.50 0.130
HJ-FP-M4 Day46 CS+ minus CS- onset 20 -2.00 0.131
HJ-FP-M4 Day47 CS+ minus CS- onset 20 -1.50 0.233
HJ-FP-M6 Day7 CS+ minus CS- onset 50 1.75 -0.092
HJ-FP-M7 Day29 CS+ minus CS- onset 20 -1.25 0.027
HJ-FP-WT-stGtACR-M1 Day2-Pavlovian CS1 onset 50 2.00 0.179
HJ-FP-WT-stGtACR-M3 Day2-Pavlovian CS1 onset 50 1.50 0.200
HJ-FP-WT-stGtACR-M3 Day3-Pavlovian CS1 onset 50 1.50 0.470
HJ-FP-datHT-stGtACR-M3 Day3-Pavlovian CS1 onset 14 1.75 -0.016
HJ-FP-datHT-stGtACR-M4 Day3-Pavlovian CS1 onset 50 1.25 0.031
HJ-FP-datHT-stGtACR-M4 Day7-Pavlovian CS1 onset 50 1.25 0.107
HJ-FP-datHT-stGtACR-M5 Day3-Pavlovian CS1 onset 50 1.25 -0.007
HJ-FP-datHT-stGtACR-M5 Day6-Pavlovian CS1 onset 50 2.00 -0.010
HJ-FP-datHT-stGtACR-M6 Day8 CS1 onset 50 1.50 0.029
HJ-FP-datHT-stGtACR-M8 Day6-Pavlovian CS1 onset 50 -2.00 -0.003
Across sessions the median lag of the largest response is 0 or -0.25 s in every group, and the quartiles lie within 0.31 s of 0; the response is positive at lag 0 in 50 of 56, 198 of 210, 24 of 25, 79 of 87 and 14 of 14 sessions, and shifting the events by 10 s removes it (medians -0.13 to 0.01). This is consistent with dF/F and the event log sharing one clock in every experiment, including Experiment 7, whose dF/F timestamps start before the event log. It is not established for the 27 sessions listed, whose peak lies more than 1 s from 0: they are sessions with no clear response to the reference event (the first random reward session of M2, the first sessions of conditioning, the background reward sessions where CS+ and CS- responses become similar, DAT-Cre sessions under inhibition), where the peak is not defined by a response.
session-day-link¶
The subjects and "Day" numbers that the paper's analysis code names for each experiment are sessions of the dandiset, so that the code's selection of sessions can be applied to it.
Subjects are matched with _ read as -, and sessions by the number after "Day" in the file name
(load_dandi_url). For each script's day list, the listed days (startday, endday, learnedday and
the single days of omission.m) are checked to be sessions of the task, and every session in each
range is checked for the condition: 2 s cue (3 s cue-reward delay) for Experiment 2, 8 s cue (9 s
delay) for Experiment 3, background rewards for Experiment 5, no delivered reward for Experiment 4,
the trial-less task (Poisson) for Experiment 6. The criterion's "every Day number from startday
to endday" is also counted literally, every integer day in each range.
cond = {(s["subject"], s["day"]): jeong.session_condition(s) for s in sessions if s["subject"] in jeong.EXP1_MICE}
present = {(s["subject"], s["day"]): s["description"] for s in sessions}
print("Subjects of the code's lists missing from the dandiset:",
[m for m in jeong.EXP1_MICE + jeong.DAT_MICE + [jeong.DAT_EXCLUDED] + jeong.CONTROL_MICE + jeong.LONGITI_MICE
if m not in set(info.subject)])
def expect(c, kind):
if c is None:
return False
d = c.get("delays", [])
return {"2s": c["task"] == "Pavlovian" and d == [3.0] and c["n_bgd"] == 0 and c["delivered"] > 0,
"8s": c["task"] == "Pavlovian" and 9.0 in d and c["n_bgd"] == 0 and c["delivered"] > 0,
"8s_or_2s": c["task"] == "Pavlovian" and c["n_bgd"] == 0 and c["delivered"] > 0,
"bgd": c["task"] == "Pavlovian" and c["n_bgd"] > 0,
"ext": c["task"] == "Pavlovian" and c.get("delivered", 1) == 0,
"poisson": c["task"] == "Poisson"}[kind]
checks = [] # (script, subject, kind, day, listed day?, present, condition ok)
def add_range(script, m, kind, a, b, listed=()):
for d in range(a, b + 1):
c = cond.get((m, d))
checks.append({"script": script, "subject": m, "condition": kind, "day": d,
"listed": d in listed, "session": c is not None,
"condition_ok": expect(c, kind) if c is not None else np.nan})
for m in jeong.PAV_MICE:
st, en = jeong.CUMSUM_START[m], jeong.CUMSUM_END[m]
add_range("cumsum_pavlovian C1 (Exp 2)", m, "2s", st[0], en[0], {st[0], en[0], jeong.LEARNED_DAY[m]})
add_range("cumsum_pavlovian C2 pre (Exp 2)", m, "2s", st[1] - 3, st[1] - 1)
add_range("cumsum_pavlovian C2 (Exp 3)", m, "8s", st[1], en[1], {st[1], en[1]})
add_range("cumsum_pavlovian C3 pre (Exp 3)", m, "8s", st[2] - 1, st[2] - 1)
add_range("cumsum_pavlovian C3 (Exp 4)", m, "ext", st[2], en[2], {st[2], en[2]})
add_range("cumsum_pavlovian C4 (Exp 5)", m, "bgd", st[3], en[3], {st[3], en[3]})
bs, be = jeong.BGD_START[m], jeong.BGD_END[m]
add_range("cumsum_bgdrw_vs_extinction (Exp 5)", m, "bgd", bs[0], be[0], {bs[0], be[0]})
add_range("cumsum_bgdrw_vs_extinction (Exp 4)", m, "ext", bs[1], be[1], {bs[1], be[1]})
add_range("intermediate_rsp_trialless (Exp 6)", m, "poisson", jeong.TRIALLESS_LEARNED[m], jeong.TRIALLESS_END[m],
{jeong.TRIALLESS_LEARNED[m], jeong.TRIALLESS_END[m]})
o3, ox = jeong.OMISSION_DAYS[m]
add_range("omission (Exp 3 first)", m, "8s", o3, o3, {o3})
add_range("omission (Exp 4 first)", m, "ext", ox, ox, {ox})
ck = pd.DataFrame(checks)
listed = ck[ck.listed]
print(f"Listed days (start, end, learned, single days): {len(listed)}; sessions: {int(listed.session.sum())}; "
f"with the condition: {int(listed.condition_ok.fillna(False).astype(bool).sum())}")
print(listed[~listed.session | (listed.condition_ok == False)].to_string(index=False))
sess = ck[ck.session]
print(f"\nSessions in the ranges: {len(sess)}; with the condition: {int(sess.condition_ok.astype(bool).sum())}")
print(sess[sess.condition_ok == False].to_string(index=False))
print(f"\nEvery integer day of the ranges: {len(ck)}; of them sessions: {int(ck.session.sum())} "
f"({ck.session.mean():.1%})")
print("Integer days in the ranges with no session:")
print(ck[~ck.session].groupby(["subject"]).day.apply(lambda d: sorted(set(d))).to_string())
Subjects of the code's lists missing from the dandiset: []
Listed days (start, end, learned, single days): 119; sessions: 117; with the condition: 117
script subject condition day listed session condition_ok
cumsum_pavlovian C4 (Exp 5) HJ-FP-M4 bgd 38 True False NaN
cumsum_bgdrw_vs_extinction (Exp 5) HJ-FP-M4 bgd 38 True False NaN
Sessions in the ranges: 320; with the condition: 320
Empty DataFrame
Columns: [script, subject, condition, day, listed, session, condition_ok]
Index: []
Every integer day of the ranges: 331; of them sessions: 320 (96.7%)
Integer days in the ranges with no session:
subject
HJ-FP-F2 [16]
HJ-FP-M2 [23]
HJ-FP-M3 [16, 31]
HJ-FP-M4 [30, 38]
HJ-FP-M6 [49]
# The mouse and day lists as the scripts define them, read from the repository at the pinned
# commit, against the copies in jeong.py used above.
code_lists = jeong.fetch_code_lists()
col = lambda d, c: [d[m][c] for m in jeong.PAV_MICE]
mine = {
"analysis/fig4/cumsum_pavlovian.m": {
"startday": [jeong.CUMSUM_START[m] for m in jeong.PAV_MICE], "endday": [jeong.CUMSUM_END[m] for m in jeong.PAV_MICE],
"learnedday": [[jeong.LEARNED_DAY[m] for m in jeong.PAV_MICE]], "ndaybeforestart": [jeong.CUMSUM_BEFORE],
"ratio": [jeong.CUMSUM_RATIO]},
"analysis/fig4/cumsum_bgdrw_vs_extinction.m": {
"startday": [jeong.BGD_START[m] for m in jeong.PAV_MICE], "endday": [jeong.BGD_END[m] for m in jeong.PAV_MICE]},
"analysis/fig5/intermediate_rsp_trialless.m": {
"learnedday": [[jeong.TRIALLESS_LEARNED[m] for m in jeong.PAV_MICE]], "endday": [[jeong.TRIALLESS_END[m] for m in jeong.PAV_MICE]]},
"analysis/figS10/omission.m": {"startday": [col(jeong.OMISSION_DAYS, 0), col(jeong.OMISSION_DAYS, 1)]},
"analysis/figS10/cue_acquisition.m": {
"startday": [col(jeong.CUMSUM_START, 0)], "learnedday": [[jeong.LEARNED_DAY[m] for m in jeong.PAV_MICE]],
"endday": [col(jeong.CUMSUM_END, 0)]},
"analysis/fig6/backpropagation_pavlovian.m": {"startday": [col(jeong.CUMSUM_START, 0)], "endday": [col(jeong.CUMSUM_END, 0)]},
}
all_match = True
for path, found in code_lists.items():
print(path)
for name, (line, value) in found.items():
if name == "mouseList":
ok = value == jeong.PAV_MICE
else:
ok = np.array_equal(value, np.array(mine[path][name], float))
all_match &= bool(ok)
print(f" {'matches' if ok else 'DIFFERS'}: {line}")
missing = set(mine[path]) - set(found)
if missing:
all_match = False
print(" not found in the script:", missing)
print("Every list in jeong.py equals the script's:", all_match)
analysis/fig4/cumsum_pavlovian.m
matches: mouseList = {'HJ_FP_M2';'HJ_FP_M3';'HJ_FP_M4';'HJ_FP_F1';'HJ_FP_F2';'HJ_FP_M6';'HJ_FP_M7'};
matches: startday = [16,29,70,38; 12,24,61,28; 13,32,66,38; 7,19,42,23; 13,24,46,28; 7,19,47,23; 6,20,43,24]; % first day or each condition
matches: endday = [28,32,72,43; 23,27,63,38; 31,37,67,47; 18,22,43,27; 23,27,47,33; 18,22,50,30; 19,23,44,29]; % last day of each condition
matches: learnedday = [22,17,25,12,17,12,14]; % first day w/ anticipatory behavior
matches: ndaybeforestart = [0,3,1,1]; % number of pre-condition sessions
matches: ratio = [1.5,1.9,2.5,2]; % ratio b/w pre- and post-condition trials - for population analysis
analysis/fig4/cumsum_bgdrw_vs_extinction.m
matches: mouseList = {'HJ_FP_M2';'HJ_FP_M3';'HJ_FP_M4';'HJ_FP_F1';'HJ_FP_F2';'HJ_FP_M6';'HJ_FP_M7'};
matches: startday = [33,70; 28,61; 38,66; 23,42; 28,46; 23,47; 24,43]; % first day or each condition
matches: endday = [43,72; 38,63; 47,67; 27,43; 33,47; 30,50; 29,44]; % last day of each condition
analysis/fig5/intermediate_rsp_trialless.m
matches: mouseList = {'HJ_FP_M2';'HJ_FP_M3';'HJ_FP_M4';'HJ_FP_F1';'HJ_FP_F2';'HJ_FP_M6';'HJ_FP_M7'};
matches: learnedday = [77,67,70,46,51,53,48]; % first day w/ anticipatory lick
matches: endday = [78,68,72,49,54,55,50];
analysis/figS10/omission.m
matches: mouseList = {'HJ_FP_M2';'HJ_FP_M3';'HJ_FP_M4';'HJ_FP_F1';'HJ_FP_F2';'HJ_FP_M6';'HJ_FP_M7'};
matches: startday = [29,24,32,19,24,19,20; 70,61,66,42,46,47,43];
analysis/figS10/cue_acquisition.m
matches: mouseList = {'HJ_FP_M2';'HJ_FP_M3';'HJ_FP_M4';'HJ_FP_F1';'HJ_FP_F2';'HJ_FP_M6';'HJ_FP_M7'};
matches: startday = [16,12,13,7,13,7,6]; % first session
matches: learnedday = [22,17,25,12,17,12,14]; % first session w/ anticipatory lick
matches: endday = [28,23,31,18,23,18,19]; % last session
analysis/fig6/backpropagation_pavlovian.m
matches: mouseList = {'HJ_FP_M2';'HJ_FP_M3';'HJ_FP_M4';'HJ_FP_F1';'HJ_FP_F2';'HJ_FP_M6';'HJ_FP_M7'};
matches: startday = [16,12,13,7,13,7,6]; % first session
matches: endday = [28,23,31,18,23,18,19]; % last session
Every list in jeong.py equals the script's: True
# The 8 s condition's first session in two animals holds both delays (the cue lengthened within
# the session); the condition of each range's sessions, by delays:
for m in jeong.PAV_MICE:
st, en = jeong.CUMSUM_START[m][1], jeong.CUMSUM_END[m][1]
print(m, {d: cond[(m, d)].get("delays") for d in range(st, en + 1) if (m, d) in cond})
HJ-FP-M2 {29: [9.0], 30: [9.0], 31: [9.0], 32: [9.0]}
HJ-FP-M3 {24: [9.0], 25: [9.0], 26: [9.0], 27: [9.0]}
HJ-FP-M4 {32: [9.0], 33: [9.0], 34: [9.0], 35: [9.0], 36: [9.0], 37: [9.0]}
HJ-FP-F1 {19: [9.0], 20: [9.0], 21: [9.0], 22: [9.0]}
HJ-FP-F2 {24: [9.0], 25: [9.0], 26: [9.0], 27: [9.0]}
HJ-FP-M6 {19: [9.0], 20: [9.0], 21: [9.0], 22: [9.0]}
HJ-FP-M7 {20: [9.0], 21: [9.0], 22: [9.0], 23: [9.0]}
Result: partially reproduced. All 28 subjects of the code's lists are in the dandiset. Of the
119 listed days, 117 are sessions with the condition; the exception is HJ-FP-M4 Day38 (listed twice), the
first day of Experiment 5 in both background reward scripts, which is not in the dandiset (its
next session, Day39, has background rewards). All 320 sessions inside the ranges have the condition the
script names. Counted literally, every integer day of the ranges, 320 of 331
(96.7%) are sessions; the 7 missing days within the Pavlovian ranges (M2 Day23, M3 Day16 and Day31, M4 Day30
and Day38, F2 Day16, M6 Day49) fall where the paper's session counts imply a session (see
methods-session-counts). The subjects match and at least 90% of the listed days do, so the link
holds in part. The claims that depend on it are tested on the sessions that are present. The mouse and day lists
in jeong.py equal those the six scripts define at the pinned commit (second output); the
start, end and learned days of cue_acquisition.m and backpropagation_pavlovian.m equal the
initial learning columns of cumsum_pavlovian.m.
methods-session-counts¶
Animals were trained for 7 +/- 2.1 sessions of random rewards, 13.3 +/- 2.7 sessions of initial conditioning, 4.3 +/- 0.8 sessions with the 8 s cue, 8.1 +/- 2.5 sessions with background rewards and 2.6 +/- 0.8 sessions of extinction; in sequential conditioning, 5.6 +/- 2.4 (wild type) and 7.1 +/- 0.7 (DAT-Cre) sessions (mean +/- SD).
Sessions are counted from the dandiset per animal. Experiment 1, the random reward sessions of the 8 mice. Experiments 2, 3 and 4, the sessions in the code's ranges (cumsum_pavlovian.m, startday to endday of conditions 1 to 3). Experiment 5, the code gives two ranges, condition 4 of cumsum_pavlovian.m and cumsum_bgdrw_vs_extinction.m; both are counted, and also the sessions with background rewards in the event log, which is the count from the data. Experiment 7, the conditioning sessions after the random reward session, all 7 mice of each group as in the Methods. The counts of integer days in each range are shown beside them.
counts = []
for m in jeong.EXP1_MICE:
counts.append({"quantity": "Exp 1", "subject": m, "sessions": len(subject_sessions(m, "RandomRewards"))})
for m in jeong.PAV_MICE:
days = {s["day"] for s in subject_sessions(m, "Pavlovian")}
st, en = jeong.CUMSUM_START[m], jeong.CUMSUM_END[m]
for q, c in [("Exp 2", 0), ("Exp 3", 1), ("Exp 4", 2), ("Exp 5 (cumsum_pavlovian C4)", 3)]:
rng_ = range(st[c], en[c] + 1)
counts.append({"quantity": q, "subject": m, "sessions": len(days & set(rng_)), "integer_days": len(rng_)})
rng_ = range(jeong.BGD_START[m][0], jeong.BGD_END[m][0] + 1)
counts.append({"quantity": "Exp 5 (cumsum_bgdrw range)", "subject": m, "sessions": len(days & set(rng_)), "integer_days": len(rng_)})
counts.append({"quantity": "Exp 5 (sessions with background rewards)", "subject": m,
"sessions": sum(cond[(m, d)]["n_bgd"] > 0 for d in days)})
for grp, mice in [("Exp 7 WT", jeong.CONTROL_MICE), ("Exp 7 DAT-Cre", jeong.DAT_MICE + [jeong.DAT_EXCLUDED])]:
for m in mice:
ss = subject_sessions(m)
assert ss[0]["description"] == "RandomRewards" and all(s["description"] == "Pavlovian" for s in ss[1:])
counts.append({"quantity": grp, "subject": m, "sessions": len(ss) - 1})
counts = pd.DataFrame(counts)
summary = counts.groupby("quantity", sort=False).agg(animals=("subject", "size"), mean=("sessions", "mean"),
sd=("sessions", "std"), min=("sessions", "min"),
max=("sessions", "max"), integer_day_mean=("integer_days", "mean"))
paper = {"Exp 1": (7, 2.1), "Exp 2": (13.3, 2.7), "Exp 3": (4.3, 0.8), "Exp 5 (cumsum_bgdrw range)": (8.1, 2.5),
"Exp 5 (sessions with background rewards)": (8.1, 2.5), "Exp 5 (cumsum_pavlovian C4)": (8.1, 2.5),
"Exp 4": (2.6, 0.8), "Exp 7 WT": (5.6, 2.4), "Exp 7 DAT-Cre": (7.1, 0.7)}
summary["paper"] = [paper[q][0] for q in summary.index]
summary["tolerance"] = [paper[q][1] for q in summary.index]
summary["within"] = (summary["mean"] - summary.paper).abs() <= summary.tolerance
print(summary.round(2).to_string())
animals mean sd min max integer_day_mean paper tolerance within quantity Exp 1 8 7.00 2.14 4 11 NaN 7.0 2.1 True Exp 2 7 12.71 2.63 10 18 13.29 13.3 2.7 True Exp 3 7 4.29 0.76 4 6 4.29 4.3 0.8 True Exp 4 7 2.43 0.53 2 3 2.57 2.6 0.8 True Exp 5 (cumsum_pavlovian C4) 7 7.14 1.86 5 10 7.43 8.1 2.5 True Exp 5 (cumsum_bgdrw range) 7 7.86 2.27 5 11 8.14 8.1 2.5 True Exp 5 (sessions with background rewards) 7 7.86 2.27 5 11 NaN 8.1 2.5 True Exp 7 WT 7 5.43 0.79 5 7 NaN 5.6 2.4 True Exp 7 DAT-Cre 7 7.00 2.45 5 12 NaN 7.1 0.7 True
Result: reproduced. Counted from the dandiset, the means are 7.0 random reward sessions (4 to 11, SD 2.14), 12.7 of Experiment 2, 4.3 of Experiment 3, 7.86 of Experiment 5 (sessions with background rewards and the cumsum_bgdrw range; 7.14 with condition 4 of cumsum_pavlovian.m), 2.43 of extinction, and 5.43 and 7.0 sequential conditioning sessions for wild type and DAT-Cre. All seven are within the stated SDs. The counts of integer days in the code's ranges, 13.3, 8.1 (cumsum_bgdrw) and 2.6, equal the paper's means, which suggests that the paper counted sessions that are not in the dandiset (session-day-link). The SDs of the Experiment 7 counts, 0.79 for wild type and 2.45 for DAT-Cre (from datHT-M7's 12 sessions), are close to the paper's 2.4 and 0.7 taken in the other order, so the paper may have swapped them; the criterion tests only the means.
Experiment 1: random rewards¶
The 12 s random reward sessions of the 8 mice, in Day order (the dandiset holds no session named
"30s", so the code's cut leaves all of them). The measures follow corr_randomrewards.m as the
criteria give it; see jeong.exp1_rewards.
exp1 = {}
for m in jeong.EXP1_MICE:
exp1[m] = [(s, jeong.exp1_rewards(s)) for s in subject_sessions(m, "RandomRewards")]
rows = []
for m, lst in exp1.items():
for k, (s, r) in enumerate(lst, start=1):
keep = ~r["out"]
ok = keep & np.isfinite(r["iri"]) & np.isfinite(r["response"])
rows.append({"subject": m, "session": k, "day": s["day"], "rewards": len(r["reward"]),
"short": r["short"].mean(), "nolick": r["nolick"].mean(),
"kept": int(keep.sum()), "kept_with_dff": int(np.isfinite(r["response"][keep]).sum()),
"r_iri": stats.pearsonr(r["iri"][ok], r["response"][ok])[0]})
exp1_sessions = pd.DataFrame(rows)
print(exp1_sessions.groupby("subject")[["rewards", "kept", "kept_with_dff"]].sum().to_string())
rewards kept kept_with_dff subject HJ-FP-F1 500 353 350 HJ-FP-F2 800 614 604 HJ-FP-M2 1100 776 766 HJ-FP-M3 700 517 511 HJ-FP-M4 800 568 557 HJ-FP-M6 600 419 416 HJ-FP-M7 400 300 298 HJ-FP-M8 700 517 511
exp1-excluded-rewards¶
In Experiment 1, 22.6 +/- 2.4% of rewards came less than 3 s after the previous reward and 5.7 +/- 0.1% had no lick before the next reward; both were excluded.
per_animal = exp1_sessions.groupby("subject")[["short", "nolick"]].mean() * 100
print("Per animal (mean over sessions, %):")
print(per_animal.round(2).to_string())
a = per_animal.mean()
b = exp1_sessions[["short", "nolick"]].mean() * 100
print(f"\nReading (a), mean over animals: short IRI {a['short']:.2f}% (SD over animals {per_animal['short'].std():.2f}), "
f"no lick {a['nolick']:.2f}% (SD {per_animal['nolick'].std():.2f})")
print(f"Reading (b), mean over {len(exp1_sessions)} sessions: short IRI {b['short']:.2f}% "
f"(SD {exp1_sessions['short'].std() * 100:.2f}), no lick {b['nolick']:.2f}% (SD {exp1_sessions['nolick'].std() * 100:.2f})")
print("Paper: 22.6 +/- 2.4 and 5.7 +/- 0.1")
Per animal (mean over sessions, %):
short nolick
subject
HJ-FP-F1 25.20 6.20
HJ-FP-F2 20.25 4.00
HJ-FP-M2 23.91 5.82
HJ-FP-M3 20.29 6.14
HJ-FP-M4 22.50 7.38
HJ-FP-M6 26.50 4.50
HJ-FP-M7 20.50 5.50
HJ-FP-M8 21.57 5.71
Reading (a), mean over animals: short IRI 22.59% (SD over animals 2.39), no lick 5.66% (SD 1.04)
Reading (b), mean over 56 sessions: short IRI 22.59% (SD 7.04), no lick 5.68% (SD 4.27)
Paper: 22.6 +/- 2.4 and 5.7 +/- 0.1
Result: reproduced. Under reading (a) 22.59% of rewards follow the previous one by less than 3 s and 5.66% have no lick before the next reward; under reading (b) 22.59% and 5.68%. Both values are within tolerance under both readings. The spreads printed (SD over animals 2.39 and 1.04) suggest the paper's 2.4 is the SD over animals of the first; its 0.1 matches neither.
fig3e-reward-response-increases¶
In experimentally naive mice given unpredicted sucrose, the dopamine response to reward increases with the number of rewards experienced, in every animal (Test 1).
r_trial = {}
for m, lst in exp1.items():
n, resp, num = 0, [], []
for s, r in lst:
keep = ~r["out"]
resp.append(r["response"][keep])
num.append(np.flatnonzero(keep) + n)
n += len(r["reward"])
R, N = np.concatenate(resp), np.concatenate(num)
ok = np.isfinite(R)
r_trial[m] = stats.pearsonr(N[ok], R[ok])[0]
r_trial = pd.Series(r_trial)
print("Per animal Pearson r, response against reward number:")
print(r_trial.round(3).to_string())
print(f"mean r {r_trial.mean():.3f}; animals with r > 0: {(r_trial > 0).sum()} of {len(r_trial)}")
show_t("one-sample t-test against 0", stats.ttest_1samp(r_trial, 0))
print("Paper: t(7) = 4.40, p = 0.0031")
Per animal Pearson r, response against reward number: HJ-FP-M2 0.355 HJ-FP-M3 0.677 HJ-FP-M4 0.068 HJ-FP-F1 0.115 HJ-FP-F2 0.157 HJ-FP-M6 0.424 HJ-FP-M7 0.294 HJ-FP-M8 0.369 mean r 0.308; animals with r > 0: 8 of 8 one-sample t-test against 0: t = 4.405, p = 0.003138, df = 7 Paper: t(7) = 4.40, p = 0.0031
Result: reproduced. The correlation between the reward response and the reward number is positive in all 8 animals (0.068 to 0.677, mean 0.308), t(7) = 4.41, p = 0.0031, the paper's values.
fig3g-response-iri-positive¶
The dopamine response to an unpredicted reward is positively correlated with the interval since the previous reward (Test 2).
r_iri = exp1_sessions.groupby("subject").r_iri.mean()
print("Per animal mean over sessions of the per-session r (response against previous IRI):")
print(r_iri.round(3).to_string())
show_t("one-sample t-test against 0", stats.ttest_1samp(r_iri, 0))
print("Paper: t(7) = 5.95, p = 5.7e-4")
Per animal mean over sessions of the per-session r (response against previous IRI): subject HJ-FP-F1 0.239 HJ-FP-F2 0.176 HJ-FP-M2 0.254 HJ-FP-M3 0.173 HJ-FP-M4 0.405 HJ-FP-M6 0.029 HJ-FP-M7 0.269 HJ-FP-M8 0.247 one-sample t-test against 0: t = 5.948, p = 0.0005715, df = 7 Paper: t(7) = 5.95, p = 5.7e-4
Result: reproduced. The mean per-session correlation with the previous IRI is positive in all 8 animals (0.029 to 0.405), t(7) = 5.95, p = 5.7e-4, the paper's values.
figs8k-iri-correlation-stable¶
The positive correlation between reward response and previous IRI holds across sessions and does not change over the first 7 sessions.
first7 = exp1_sessions[exp1_sessions.session <= 7]
slope, t, p, n = ols_slope(first7.session, first7.r_iri)
print(f"Linear regression of per-session r on session number (1 to 7), n = {n} sessions: slope {slope:.4f}, t = {t:.2f}, p = {p:.4f}")
print("Paper: t = -0.04, p = 0.9707, n = 50")
print("\nMean r over animals per session (and animals):")
print(first7.groupby("session").r_iri.agg(["mean", "size"]).round(3).to_string())
Linear regression of per-session r on session number (1 to 7), n = 50 sessions: slope -0.0005, t = -0.04, p = 0.9686
Paper: t = -0.04, p = 0.9707, n = 50
Mean r over animals per session (and animals):
mean size
session
1 0.188 8
2 0.269 8
3 0.208 8
4 0.275 8
5 0.193 7
6 0.134 6
7 0.288 5
Result: reproduced. Over the 50 sessions the slope of r on session number is -0.0005 per session, t = -0.04, p = 0.97 (paper t = -0.04, p = 0.9707), and the mean r over animals is positive in each of sessions 1 to 7 (0.134 to 0.288).
figs8b-lick-rates-constant¶
Total, consummatory and non-consummatory lick rates do not change across the first 7 random reward sessions.
lick_rows, lick_rw = [], {}
for m, lst in exp1.items():
L, R = [], []
for k, (s, r) in enumerate(lst, start=1):
lk = jeong.exp1_licks(s)
lick_rows.append({"subject": m, "session": k, **lk["rates"]})
L.append(lk["lick_rate_rw"]); R.append(lk["response_rw"])
exp1[m][k - 1] = (s, r, lk)
lick_rw[m] = (np.concatenate(L), np.concatenate(R))
licks = pd.DataFrame(lick_rows)
l7 = licks[licks.session <= 7]
for c in ["total", "consummatory", "nonconsummatory"]:
slope, t, p, n = ols_slope(l7.session, l7[c])
print(f"{c:16s} rate (licks/s): mean {l7[c].mean():.2f}, slope {slope:.4f} per session, t = {t:.2f}, p = {p:.4f}, n = {n}")
print("Paper: t = -0.79, -0.20, -0.58; p = 0.4389, 0.8397, 0.5618")
total rate (licks/s): mean 1.86, slope -0.0382 per session, t = -0.78, p = 0.4389, n = 50 consummatory rate (licks/s): mean 5.31, slope -0.0082 per session, t = -0.20, p = 0.8397, n = 50 nonconsummatory rate (licks/s): mean 0.93, slope -0.0271 per session, t = -0.58, p = 0.5618, n = 50 Paper: t = -0.79, -0.20, -0.58; p = 0.4389, 0.8397, 0.5618
Result: reproduced. None of the three lick rates changes with session over the 50 sessions: t = -0.78 (p = 0.4389), -0.20 (p = 0.8397) and -0.58 (p = 0.5618), matching Table S1.
figs8d-consummatory-lick-uncorrelated¶
The consummatory lick rate after a reward is not correlated with the dopamine response to that reward.
Result: reproduced. The mean correlation is 0.010, t(7) = 0.23, p = 0.82 (paper t = 0.23, p = 0.8245): no correlation is found.
figs8e-lick-bout-responses¶
Dopamine rises at the onset of consummatory lick bouts but falls slightly at the onset of non-consummatory (dry) lick bouts, so the rising reward response is not a general response to licking.
windows = {"a: 0 to 1 s": (0.0, 1.0), "b: -0.5 to 1 s": (-0.5, 1.0), "c: 0 to 0.5 s": (0.0, 0.5)}
zrows = []
for m, lst in exp1.items():
for k, (s, r, lk) in enumerate(lst, start=1):
z = jeong.bout_onset_z(s, lk, windows)
zrows.append({"subject": m, "session": k, **{f"{a}|{b}": v for (a, b), v in z.items()}})
zdf = pd.DataFrame(zrows)
za = zdf.groupby("subject").mean(numeric_only=True)
print("Onsets per session (median): consummatory", zdf["consummatory|n"].median(), "; non-consummatory", zdf["nonconsummatory|n"].median())
for w in windows:
c, n = za[f"consummatory|{w}"], za[f"nonconsummatory|{w}"]
print(f"\nReading {w}: consummatory mean z {c.mean():.3f}, non-consummatory {n.mean():.3f} "
f"(non-consummatory below 0 in {(n < 0).sum()} of 8 animals)")
show_t(" paired, consummatory against non-consummatory", stats.ttest_rel(c, n))
show_t(" consummatory against 0", stats.ttest_1samp(c, 0))
show_t(" non-consummatory against 0", stats.ttest_1samp(n, 0))
print("\nPaper: paired t = 3.28, p = 0.0134; consummatory t = 3.03, p = 0.0191; non-consummatory t = -2.55, p = 0.0379")
Onsets per session (median): consummatory 95.0 ; non-consummatory 124.0 Reading a: 0 to 1 s: consummatory mean z 2.149, non-consummatory -0.147 (non-consummatory below 0 in 6 of 8 animals) paired, consummatory against non-consummatory: t = 7.165, p = 0.000183, df = 7 consummatory against 0: t = 6.300, p = 0.000404, df = 7 non-consummatory against 0: t = -1.652, p = 0.1426, df = 7 Reading b: -0.5 to 1 s: consummatory mean z 1.755, non-consummatory -0.023 (non-consummatory below 0 in 5 of 8 animals) paired, consummatory against non-consummatory: t = 7.978, p = 9.274e-05, df = 7 consummatory against 0: t = 6.451, p = 0.0003498, df = 7 non-consummatory against 0: t = -0.268, p = 0.7968, df = 7 Reading c: 0 to 0.5 s: consummatory mean z 3.066, non-consummatory -0.347 (non-consummatory below 0 in 6 of 8 animals) paired, consummatory against non-consummatory: t = 8.580, p = 5.816e-05, df = 7 consummatory against 0: t = 6.492, p = 0.0003366, df = 7 non-consummatory against 0: t = -2.559, p = 0.03758, df = 7 Paper: paired t = 3.28, p = 0.0134; consummatory t = 3.03, p = 0.0191; non-consummatory t = -2.55, p = 0.0379
Result: blocked (methods_underspecified). Under every reading the consummatory onset response is positive (mean z 1.76 to 3.07, p <= 4.0e-4) and above the non-consummatory one (p <= 1.8e-4). The non-consummatory response is negative on average under all three readings, but significantly so only from 0 to 0.5 s (reading c, z = -0.35, t = -2.56, p = 0.0376, close to the paper's -2.55, p = 0.0379); from 0 to 1 s and from -0.5 to 1 s it is -0.15 (p = 0.14) and -0.02 (p = 0.80). The window is fixed by neither the Methods, the legend nor the code, and the readings give different outcomes, so under the claims file's rule for readings of an open method detail the claim is blocked. The consummatory t values against 0 (6.30 to 6.49) are about twice the paper's 3.03, so the paper's measure of the response differs from these readings in a way not found here.
figs8g-baseline-unchanged¶
Spontaneous baseline dopamine events, away from rewards and licks, do not change in size across random reward sessions.
brows = []
for m, lst in exp1.items():
thr = jeong.peak_threshold(lst[0][0])
for k, (s, r, lk) in enumerate(lst, start=1):
bp = jeong.baseline_peaks(s, thr)
brows.append({"subject": m, "session": k, "threshold": thr,
"reward_response": np.nanmean(r["response"][~r["short"]]), **bp})
base = pd.DataFrame(brows)
print("Peaks per session (all above threshold, kept by the Methods' rule, kept by the code's rule):")
print(base[["n_peaks", "n_methods", "n_code"]].describe().loc[["min", "50%", "max"]].to_string())
bcorr = base.groupby("subject").apply(lambda g: pd.Series({
"r_methods": stats.pearsonr(g.session, g.methods)[0], "r_code": stats.pearsonr(g.session, g["code"])[0],
"r_reward": stats.pearsonr(g.session, g.reward_response)[0], "sessions": len(g)}))
print(bcorr.round(3).to_string())
for c in ["r_methods", "r_code"]:
print(f"{c}: mean r {bcorr[c].mean():.3f}"); show_t(" one-sample t-test against 0", stats.ttest_1samp(bcorr[c], 0))
print("Paper: t = 0.11, p = 0.9190")
Peaks per session (all above threshold, kept by the Methods' rule, kept by the code's rule):
n_peaks n_methods n_code
min 153.0 5.0 4.0
50% 562.5 194.0 246.0
max 1023.0 516.0 597.0
r_methods r_code r_reward sessions
subject
HJ-FP-F1 0.917 0.954 0.546 5.0
HJ-FP-F2 -0.143 -0.142 0.549 8.0
HJ-FP-M2 -0.860 -0.755 0.746 11.0
HJ-FP-M3 0.965 0.954 0.881 7.0
HJ-FP-M4 -0.219 -0.181 0.234 8.0
HJ-FP-M6 0.576 0.576 0.774 6.0
HJ-FP-M7 -0.978 -0.955 0.922 4.0
HJ-FP-M8 -0.532 -0.304 0.818 7.0
r_methods: mean r -0.034
one-sample t-test against 0: t = -0.126, p = 0.9033, df = 7
r_code: mean r 0.018
one-sample t-test against 0: t = 0.071, p = 0.9454, df = 7
Paper: t = 0.11, p = 0.9190
Result: reproduced. With peaks kept by the Methods' rule, the mean correlation of peak size with session is -0.034, t(7) = -0.13, p = 0.90; with the code's rule, t(7) = 0.07, p = 0.95 (paper t = 0.11, p = 0.919). No change is found. Per animal the correlations range from -0.98 to 0.97.
figs8h-baseline-uncorrelated-with-reward-change¶
Across animals, the change of baseline dopamine activity over sessions is not correlated with the change of the reward response over sessions.
Result: reproduced. Across the 8 animals, the session correlation of baseline peak size is not correlated with that of the reward response (Methods' rule r = -0.10, p = 0.81; code's rule r = -0.09, p = 0.83), while the reward response rises with session in every animal (r 0.23 to 0.92).
figs8l-lick-rate-tracks-iri¶
In a separate group of 6 mice, the non-consummatory lick rate depends on the session's mean IRI (6, 12 or 30 s) while the consummatory lick rate does not, so animals learn the average reward rate within at most two sessions.
The IRI of each session is read from its name (Day1-30 is 30 s), and checked against the
intervals in the data check above.
lrows = []
for m in jeong.LONGITI_MICE:
for s in subject_sessions(m):
lrows.append({"subject": m, "day": s["day"], "iri": int(s["session"].split("-")[1]), **jeong.longiti_rates(s)})
lr = pd.DataFrame(lrows)
print(pd.crosstab(lr.subject, lr.iri).to_string())
g = lr.groupby(["subject", "iri"])[["consummatory", "nonconsummatory"]].mean().reset_index()
g["rank"] = g.iri.map({6: 1, 12: 2, 30: 3})
print(g.groupby("iri")[["consummatory", "nonconsummatory"]].mean().round(3).to_string())
for reg in ["iri", "rank"]:
for c in ["nonconsummatory", "consummatory"]:
slope, t, p, n = ols_slope(g[reg], g[c])
print(f"reading {'(a) IRI in s' if reg == 'iri' else '(b) rank'}: {c} slope {slope:.4f}, t = {t:.2f}, p = {p:.4f}, n = {n}")
print("Paper: non-consummatory t = -2.73, p = 0.0147; consummatory t = 1.00, p = 0.3338")
iri 6 12 30
subject
DB-longITI-C1-F1 2 2 2
DB-longITI-C1-F2 2 2 2
DB-longITI-C1-M1 2 2 2
DB-longITI-C1-M2 2 2 2
DB-longITI-C1-M3 2 2 2
DB-longITI-C1-M4 2 2 2
consummatory nonconsummatory
iri
6 5.392 0.851
12 5.843 0.721
30 5.794 0.379
reading (a) IRI in s: nonconsummatory slope -0.0195, t = -2.73, p = 0.0147, n = 18
reading (a) IRI in s: consummatory slope 0.0123, t = 1.00, p = 0.3338, n = 18
reading (b) rank: nonconsummatory slope -0.2356, t = -2.61, p = 0.0189, n = 18
reading (b) rank: consummatory slope 0.2013, t = 1.34, p = 0.1992, n = 18
Paper: non-consummatory t = -2.73, p = 0.0147; consummatory t = 1.00, p = 0.3338
Result: reproduced. With the IRI in seconds the non-consummatory rate falls with the IRI (t = -2.73, p = 0.0147) and the consummatory rate does not change (t = 1.00, p = 0.334), the paper's values; with the rank of the condition, t = -2.61 (p = 0.019) and t = 1.34 (p = 0.199).
Experiments 2 to 5: Pavlovian conditioning¶
Per animal, the CS+ is the cue whose reward is delivered in the first session of Experiment 2.
Per CS+ trial jeong.cs_plus_trials computes the dopamine and lick measures of the criteria.
pav = {}
csplus = {}
for m in jeong.PAV_MICE:
ss = subject_sessions(m, "Pavlovian")
first = [s for s in ss if s["day"] == jeong.CUMSUM_START[m][0]][0]
csplus[m] = jeong.cs_plus_identity(first)
pav[m] = {s["day"]: jeong.cs_plus_trials(s, csplus[m]) for s in ss}
unpaired = [(s["subject"], s["day"]) for s in sessions if s["subject"] in jeong.PAV_MICE and s["description"] == "Pavlovian"
and not jeong.pavlovian_events(s)["paired"]]
print("CS+ identity (1 = code 15, 2 = code 16):", csplus)
print("Sessions where cues and reward events do not pair one to one:", unpaired)
nan_trials = sum(int(np.isnan(v["da"]).sum()) for m in pav for v in pav[m].values())
inc_trials = sum(int(v["da_incomplete"].sum()) for m in pav for v in pav[m].values())
print("CS+ trials:", sum(len(v["da"]) for m in pav for v in pav[m].values()), "; without a dopamine value:", nan_trials,
"; with a NaN smoothed dF/F sample from -1 to 2 s (summed by nansum, as aucsignal.m does):", inc_trials)
R = jeong.matlab_round
def cumsum_condition(m, c, drop_incomplete=False):
# cumsum_pavlovian.m for condition c (0 initial learning, 1 cue duration, 2 extinction).
# With drop_incomplete, trials with a NaN smoothed dF/F sample from -1 to 2 s are left out of
# the series (a variant; the code keeps them with their nansum AUC).
st, en = jeong.CUMSUM_START[m][c], jeong.CUMSUM_END[m][c]
days = [d for d in sorted(pav[m]) if st - jeong.CUMSUM_BEFORE[c] <= d <= en]
cut = jeong.LEARNED_DAY[m] if c == 0 else st
out = {"days": days}
keep = {d: np.ones(len(pav[m][d]["da"]), bool) for d in days}
if c == 1:
keep = {d: np.isfinite(pav[m][d]["lick_median_ms"]) for d in days}
out["dropped_no_lick"] = sum(int((~keep[d]).sum()) for d in days)
if drop_incomplete:
keep = {d: keep[d] & ~pav[m][d]["da_incomplete"] for d in days}
cat = lambda f: np.concatenate([f(pav[m][d])[keep[d]] for d in days])
nb = sum(int(keep[d].sum()) for d in days if d < cut)
da = cat(lambda v: v["da"])
if c == 1:
out["lick_series"] = {"median lick time": cat(lambda v: v["lick_median_ms"])}
else:
out["lick_series"] = {"Methods 0-2 s": cat(lambda v: v["lick_methods"]),
"code cue+trace": cat(lambda v: v["lick_code_cs"] - v["lick_code_base"])}
nt = min(R(nb * jeong.CUMSUM_RATIO[c]), len(da)) # the cap binds only when trials are left out
out.update(nb=nb, nt=nt, da=da, n_incomplete=int(cat(lambda v: v["da_incomplete"])[:nt].sum()))
out["da_abrupt"], out["da_change"] = jeong.abruptness(np.nancumsum(da), nt)
out["lick"] = {k: jeong.abruptness(np.cumsum(v), nt) for k, v in out["lick_series"].items()}
return out
def abrupt_table(c, drop_incomplete=False):
rows = []
for m in jeong.PAV_MICE:
o = cumsum_condition(m, c, drop_incomplete)
for k, (ab, ch) in o["lick"].items():
rows.append({"subject": m, "lick_measure": k, "trials_before": o["nb"], "trials_used": o["nt"],
"incomplete": o["n_incomplete"],
"lick_abrupt": ab, "da_abrupt": o["da_abrupt"], "lick_change": ch, "da_change": o["da_change"]})
return pd.DataFrame(rows)
def report_abrupt(tab, brief=False):
for k, g in tab.groupby("lick_measure", sort=False):
print(f"\nLick measure: {k}")
if not brief:
print(g.drop(columns="lick_measure").round(3).to_string(index=False))
print(f"abruptness lick {g.lick_abrupt.mean():.3f} against dopamine {g.da_abrupt.mean():.3f}; "
f"change trial lick {g.lick_change.mean():.1f} against dopamine {g.da_change.mean():.1f}")
show_t(" paired t, abruptness lick minus dopamine", stats.ttest_rel(g.lick_abrupt, g.da_abrupt))
show_t(" paired t, change trial lick minus dopamine", stats.ttest_rel(g.lick_change, g.da_change))
CS+ identity (1 = code 15, 2 = code 16): {'HJ-FP-M2': 1, 'HJ-FP-M3': 1, 'HJ-FP-M4': 2, 'HJ-FP-F1': 2, 'HJ-FP-F2': 1, 'HJ-FP-M6': 1, 'HJ-FP-M7': 2}
Sessions where cues and reward events do not pair one to one: []
CS+ trials: 8721 ; without a dopamine value: 0 ; with a NaN smoothed dF/F sample from -1 to 2 s (summed by nansum, as aucsignal.m does): 37
fig4ef-dopamine-leads-behavior¶
During initial cue-reward learning, the dopamine response to CS+ develops before anticipatory licking, and behavioral learning is much more abrupt than dopaminergic learning (Test 3).
tab_c1 = abrupt_table(0)
report_abrupt(tab_c1)
print("\nVariant: trials with a NaN dF/F sample from -1 to 2 s left out")
report_abrupt(abrupt_table(0, drop_incomplete=True), brief=True)
print("\nPaper: abruptness t(6) = 9.06, p = 1.0e-4; change trial t(6) = -2.93, p = 0.0263")
Lick measure: Methods 0-2 s subject trials_before trials_used incomplete lick_abrupt da_abrupt lick_change da_change HJ-FP-M2 300 450 3 0.778 0.177 349 131 HJ-FP-M3 200 300 1 1.452 0.270 230 146 HJ-FP-M4 600 900 0 0.701 0.182 608 238 HJ-FP-F1 250 375 1 0.572 0.298 347 167 HJ-FP-F2 127 191 3 0.863 0.305 121 121 HJ-FP-M6 250 375 0 0.683 0.403 273 232 HJ-FP-M7 373 560 3 0.898 0.215 430 202 abruptness lick 0.850 against dopamine 0.264; change trial lick 336.9 against dopamine 176.7 paired t, abruptness lick minus dopamine: t = 5.063, p = 0.002304, df = 6 paired t, change trial lick minus dopamine: t = 3.315, p = 0.01611, df = 6 Lick measure: code cue+trace subject trials_before trials_used incomplete lick_abrupt da_abrupt lick_change da_change HJ-FP-M2 300 450 3 0.699 0.177 333 131 HJ-FP-M3 200 300 1 0.870 0.270 212 146 HJ-FP-M4 600 900 0 0.687 0.182 602 238 HJ-FP-F1 250 375 1 0.570 0.298 270 167 HJ-FP-F2 127 191 3 0.782 0.305 121 121 HJ-FP-M6 250 375 0 0.663 0.403 275 232 HJ-FP-M7 373 560 3 0.719 0.215 380 202 abruptness lick 0.713 against dopamine 0.264; change trial lick 313.3 against dopamine 176.7 paired t, abruptness lick minus dopamine: t = 9.088, p = 9.962e-05, df = 6 paired t, change trial lick minus dopamine: t = 2.930, p = 0.02628, df = 6 Variant: trials with a NaN dF/F sample from -1 to 2 s left out Lick measure: Methods 0-2 s abruptness lick 0.856 against dopamine 0.265; change trial lick 335.6 against dopamine 175.9 paired t, abruptness lick minus dopamine: t = 5.152, p = 0.00211, df = 6 paired t, change trial lick minus dopamine: t = 3.310, p = 0.01619, df = 6 Lick measure: code cue+trace abruptness lick 0.714 against dopamine 0.265; change trial lick 312.1 against dopamine 175.9 paired t, abruptness lick minus dopamine: t = 9.145, p = 9.619e-05, df = 6 paired t, change trial lick minus dopamine: t = 2.928, p = 0.02635, df = 6 Paper: abruptness t(6) = 9.06, p = 1.0e-4; change trial t(6) = -2.93, p = 0.0263
Result: reproduced. With the Methods' lick measure, lick abruptness exceeds dopamine abruptness (paired t = 5.06, p = 0.0023) and the dopamine change trial comes earlier (lick minus dopamine t = 3.31, p = 0.016). With the code's lick measure, t = 9.09 (p = 1.0e-4) and t = 2.93 (p = 0.026), the paper's values up to the sign of the change-trial t, which the paper prints as dopamine minus lick. Of the trials used, 0 to 3 per animal have a NaN smoothed dF/F sample from -1 to 2 s (kept with their nansum AUC, as aucsignal.m does); without them the four t values are 5.15, 3.31, 9.15 and 2.93.
figs10c-cue-response-at-acquisition¶
By the session at which behavior was acquired, the dopamine cue response had reached 93% of its value on the final session (about 100% relative to the mean of the last two sessions).
rows = []
for m in jeong.PAV_MICE:
st, en, ln = jeong.CUMSUM_START[m][0], jeong.CUMSUM_END[m][0], jeong.LEARNED_DAY[m]
days = [d for d in sorted(pav[m]) if st <= d <= en]
cr = np.array([np.nanmean(pav[m][d]["early"]) for d in days])
cs_ = np.array([np.nanmean(pav[m][d]["early"][~pav[m][d]["unsmoothed_incomplete"]]) for d in days])
x = days.index(ln)
rows.append({"subject": m, "sessions": len(days), "learned_index": x + 1, "x": cr[x], "n-1": cr[-2], "n": cr[-1],
"x_over_n": cr[x] / cr[-1], "x_over_last2": cr[x] / np.mean(cr[-2:]),
"incomplete": int(sum(pav[m][d]["unsmoothed_incomplete"][:].sum() for d in (days[x], days[-2], days[-1]))),
"strict_x_over_n": cs_[x] / cs_[-1], "strict_x_over_last2": cs_[x] / np.mean(cs_[-2:])})
acq = pd.DataFrame(rows)
print(acq.round(3).to_string(index=False))
print(f"Reading (a), mean of ratios: x/n {acq.x_over_n.mean():.3f}, x/mean(n-1, n) {acq.x_over_last2.mean():.3f}")
print(f"Reading (b), ratio of means: x/n {acq.x.mean() / acq.n.mean():.3f}, x/mean(n-1, n) {acq.x.mean() / ((acq.n + acq['n-1']) / 2).mean():.3f}")
# Variant: the Fig 4 baseline-subtracted 0 to 2 s measure.
var = []
for m in jeong.PAV_MICE:
st, en, ln = jeong.CUMSUM_START[m][0], jeong.CUMSUM_END[m][0], jeong.LEARNED_DAY[m]
days = [d for d in sorted(pav[m]) if st <= d <= en]
cr = np.array([np.nanmean(pav[m][d]["da"]) for d in days]); x = days.index(ln)
var.append((cr[x] / cr[-1], cr[x] / np.mean(cr[-2:])))
var = np.array(var)
print(f"Variant, Fig 4 measure: mean x/n {var[:, 0].mean():.3f}, x/mean(n-1, n) {var[:, 1].mean():.3f}")
print(f"Without trials with a NaN unsmoothed dF/F sample from -1 to 3 s, reading (a): "
f"x/n {acq.strict_x_over_n.mean():.3f}, x/mean(n-1, n) {acq.strict_x_over_last2.mean():.3f}")
print("Paper: 0.93 and about 1.0")
subject sessions learned_index x n-1 n x_over_n x_over_last2 incomplete strict_x_over_n strict_x_over_last2 HJ-FP-M2 12 7 24.618 28.930 27.334 0.901 0.875 0 0.901 0.875 HJ-FP-M3 11 5 42.167 52.098 71.358 0.591 0.683 0 0.591 0.683 HJ-FP-M4 18 13 44.884 29.275 33.766 1.329 1.424 0 1.329 1.424 HJ-FP-F1 12 6 6.918 11.333 8.751 0.791 0.689 0 0.791 0.689 HJ-FP-F2 10 4 54.039 49.881 45.755 1.181 1.130 1 1.207 1.155 HJ-FP-M6 12 6 17.245 19.885 21.109 0.817 0.841 0 0.817 0.841 HJ-FP-M7 14 9 11.391 13.149 12.278 0.928 0.896 0 0.928 0.896 Reading (a), mean of ratios: x/n 0.934, x/mean(n-1, n) 0.934 Reading (b), ratio of means: x/n 0.913, x/mean(n-1, n) 0.947 Variant, Fig 4 measure: mean x/n 0.985, x/mean(n-1, n) 0.955 Without trials with a NaN unsmoothed dF/F sample from -1 to 3 s, reading (a): x/n 0.938, x/mean(n-1, n) 0.938 Paper: 0.93 and about 1.0
Result: reproduced. Under reading (a), the code's, the cue response on the learned session is 0.934 of the last session's and 0.934 of the mean of the last two, both within tolerance (0.93 +/- 0.093 and 1.0 +/- 0.1); under reading (b), 0.913 and 0.947. Per animal the ratio to the last session ranges from 0.59 to 1.33. One trial of the sessions compared (F2) has a NaN dF/F sample from -1 to 3 s; without it reading (a) gives 0.938 and 0.938.
fig4hi-delay-change-behavior-changes¶
When the cue is lengthened from 2 s to 8 s (cue-reward delay 3 s to 9 s), animals change the timing of anticipatory licking abruptly, and this change is better fit by a changing model than by a constant one.
The changing model is fitted by bounded least squares from 20 random start points (Matlab's fit
starts custom models at random points; the best fit is kept), with this section's own generator.
rng_fig4hi = np.random.default_rng(7)
tab_c2 = abrupt_table(1)
report_abrupt(tab_c2)
fits = []
for m in jeong.PAV_MICE:
o = cumsum_condition(m, 1)
nb, nt = o["nb"], o["nt"]
ly, dy = o["lick_series"]["median lick time"][:nt], o["da"][:nt]
rss_l = jeong.fit_weibull(ly, [0, nb - 10, 0, 0], [9000, nt, 100, 4000], rng_fig4hi)
rss_d = jeong.fit_weibull(dy, [-2, nb - 10, 0, -10], [50, nt, 10, 10], rng_fig4hi)
rc_l, rc_d = np.sum((ly - ly.mean()) ** 2), np.sum((dy - dy.mean()) ** 2)
fits.append({"subject": m, "trials": nt, "dropped_no_lick": o["dropped_no_lick"],
"dAIC_lick_k4": jeong.aic(rc_l, 1, nt) - jeong.aic(rss_l, 4, nt),
"dAIC_da_k4": jeong.aic(rc_d, 1, nt) - jeong.aic(rss_d, 4, nt)})
o = cumsum_condition(m, 1, drop_incomplete=True)
ly, dy = o["lick_series"]["median lick time"][:o["nt"]], o["da"][:o["nt"]]
fits[-1]["incomplete"] = cumsum_condition(m, 1)["n_incomplete"]
fits[-1]["dAIC_lick_strict"] = jeong.aic(np.sum((ly - ly.mean()) ** 2), 1, o["nt"]) - jeong.aic(
jeong.fit_weibull(ly, [0, o["nb"] - 10, 0, 0], [9000, o["nt"], 100, 4000], rng_fig4hi), 4, o["nt"])
fits[-1]["dAIC_da_strict"] = jeong.aic(np.sum((dy - dy.mean()) ** 2), 1, o["nt"]) - jeong.aic(
jeong.fit_weibull(dy, [-2, o["nb"] - 10, 0, -10], [50, o["nt"], 10, 10], rng_fig4hi), 4, o["nt"])
fits = pd.DataFrame(fits)
fits["dAIC_lick_k3"] = fits.dAIC_lick_k4 + 2
fits["dAIC_da_k3"] = fits.dAIC_da_k4 + 2
print("\n", fits.round(2).to_string(index=False))
for c in ["dAIC_lick_k4", "dAIC_lick_k3"]:
print(f"{c}: mean {fits[c].mean():.1f}"); show_t(" one-sample t against 0", stats.ttest_1samp(fits[c], 0))
show_t("dAIC_lick_strict (trials with a NaN dF/F sample from -1 to 2 s left out)", stats.ttest_1samp(fits.dAIC_lick_strict, 0))
print("\nVariant: trials with a NaN dF/F sample from -1 to 2 s left out")
report_abrupt(abrupt_table(1, drop_incomplete=True), brief=True)
print("Paper: abruptness t(6) = 22.92, p = 4.52e-7; lick delta AIC t(6) = 7.49")
Lick measure: median lick time subject trials_before trials_used incomplete lick_abrupt da_abrupt lick_change da_change HJ-FP-M2 146 277 0 0.232 0.063 152 139 HJ-FP-M3 149 283 2 0.250 0.038 152 103 HJ-FP-M4 88 167 0 0.205 0.084 86 97 HJ-FP-F1 146 277 1 0.249 0.094 150 106 HJ-FP-F2 150 285 2 0.233 0.052 148 101 HJ-FP-M6 149 283 0 0.242 0.053 149 175 HJ-FP-M7 143 272 1 0.219 0.064 157 182 abruptness lick 0.233 against dopamine 0.064; change trial lick 142.0 against dopamine 129.0 paired t, abruptness lick minus dopamine: t = 15.337, p = 4.854e-06, df = 6 paired t, change trial lick minus dopamine: t = 1.010, p = 0.3513, df = 6
subject trials dropped_no_lick dAIC_lick_k4 dAIC_da_k4 incomplete dAIC_lick_strict dAIC_da_strict dAIC_lick_k3 dAIC_da_k3 HJ-FP-M2 277 21 485.94 -6.00 0 485.94 -6.00 487.94 -4.00 HJ-FP-M3 283 16 456.44 -168.49 2 452.38 -167.46 458.44 -166.49 HJ-FP-M4 167 19 118.12 32.74 0 118.12 32.74 120.12 34.74 HJ-FP-F1 277 59 471.22 8.85 1 470.76 9.33 473.22 10.85 HJ-FP-F2 285 12 300.43 -6.00 2 296.99 -6.00 302.43 -4.00 HJ-FP-M6 283 12 476.45 -6.00 0 476.45 -6.00 478.45 -4.00 HJ-FP-M7 272 77 370.87 14.33 1 367.96 13.16 372.87 16.33 dAIC_lick_k4: mean 382.8 one-sample t against 0: t = 7.492, p = 0.0002923, df = 6 dAIC_lick_k3: mean 384.8 one-sample t against 0: t = 7.531, p = 0.000284, df = 6 dAIC_lick_strict (trials with a NaN dF/F sample from -1 to 2 s left out): t = 7.462, p = 0.0002987, df = 6 Variant: trials with a NaN dF/F sample from -1 to 2 s left out Lick measure: median lick time abruptness lick 0.233 against dopamine 0.064; change trial lick 141.7 against dopamine 128.7 paired t, abruptness lick minus dopamine: t = 15.400, p = 4.739e-06, df = 6 paired t, change trial lick minus dopamine: t = 1.010, p = 0.3513, df = 6 Paper: abruptness t(6) = 22.92, p = 4.52e-7; lick delta AIC t(6) = 7.49
Result: reproduced. Lick timing changes more abruptly than the dopamine cue response (paired t = 15.34, p = 4.9e-6; paper t = 22.92), and the changing model fits lick timing better than the constant one in every animal (mean delta AIC 382.8, t = 7.49, p = 2.9e-4 with k = 4, the paper's t; t = 7.53 with k = 3). Of the trials used, 0 to 2 per animal have a NaN smoothed dF/F sample from -1 to 2 s; without them the abruptness t is 15.40 and the lick delta AIC t is 7.46. The abruptness t is smaller than the paper's while the delta AIC t equals it; the cause is not found here. One candidate is HJ-FP-M4, whose Day30, inside the three sessions before the delay change, is not in the dandiset, so that it has 88 trials before the change where the others have 143 to 150.
fig4hi-delay-change-dopamine-unchanged¶
After the cue-reward delay is lengthened, the dopamine cue response shows no significant change (Test 4).
for c in ["dAIC_da_k4", "dAIC_da_k3"]:
print(f"{c}: mean {fits[c].mean():.2f}, per animal {fits[c].round(1).tolist()}")
show_t(" one-sample t against 0", stats.ttest_1samp(fits[c], 0))
print("Paper: t(6) = -0.86, p = 0.42")
# The code's bounds on the dopamine fit hold on the scale assumed here (dF/F as a fraction).
# Variants: no bounds (L and S kept above 0), the bounds applied to dF/F in percent (the stored
# scale, 100 times larger values), and without HJ-FP-M3, whose changing model fits worse than
# the constant one.
rng_fig4hi_da = np.random.default_rng(8)
da_var = []
for m in jeong.PAV_MICE:
o = cumsum_condition(m, 1)
nb, nt = o["nb"], o["nt"]
dy = o["da"][:nt]
row = {"subject": m, "mean_response": dy.mean(), "sd_response": dy.std()}
for name, y, bounded in [("unbounded", dy, False), ("percent_scale", dy * 100, True)]:
rss = jeong.fit_weibull(y, [-2, nb - 10, 0, -10], [50, nt, 10, 10], rng_fig4hi_da, bounded=bounded)
row[name] = jeong.aic(np.sum((y - y.mean()) ** 2), 1, nt) - jeong.aic(rss, 4, nt)
da_var.append(row)
da_var = pd.DataFrame(da_var)
print(da_var.round(2).to_string(index=False))
for c in ["unbounded", "percent_scale"]:
show_t(f"{c}: delta AIC (k = 4) against 0", stats.ttest_1samp(da_var[c], 0))
show_t("fraction scale, bounded, without HJ-FP-M3", stats.ttest_1samp(fits[fits.subject != "HJ-FP-M3"].dAIC_da_k4, 0))
show_t("trials with a NaN dF/F sample from -1 to 2 s left out", stats.ttest_1samp(fits.dAIC_da_strict, 0))
dAIC_da_k4: mean -18.65, per animal [-6.0, -168.5, 32.7, 8.8, -6.0, -6.0, 14.3] one-sample t against 0: t = -0.730, p = 0.4929, df = 6 dAIC_da_k3: mean -16.65, per animal [-4.0, -166.5, 34.7, 10.8, -4.0, -4.0, 16.3] one-sample t against 0: t = -0.652, p = 0.5387, df = 6 Paper: t(6) = -0.86, p = 0.42
subject mean_response sd_response unbounded percent_scale HJ-FP-M2 17.80 6.65 34.59 -577.84 HJ-FP-M3 44.55 10.85 19.31 -817.40 HJ-FP-M4 24.25 8.97 37.16 -354.66 HJ-FP-F1 5.37 3.96 14.06 -271.66 HJ-FP-F2 27.75 9.82 35.13 -625.18 HJ-FP-M6 13.22 5.46 37.61 -538.16 HJ-FP-M7 9.32 4.06 22.86 -488.08 unbounded: delta AIC (k = 4) against 0: t = 7.826, p = 0.0002298, df = 6 percent_scale: delta AIC (k = 4) against 0: t = -7.739, p = 0.0002444, df = 6 fraction scale, bounded, without HJ-FP-M3: t = 0.989, p = 0.3679, df = 5 trials with a NaN dF/F sample from -1 to 2 s left out: t = -0.733, p = 0.4912, df = 6
Result: reproduced, resting on an assumption. For the dopamine cue response the mean delta AIC is -18.7 (constant model better on average), t(6) = -0.73, p = 0.49 with k = 4 (paper t = -0.86, p = 0.42), and t = -0.65, p = 0.54 with k = 3. No change is found; without HJ-FP-M3, t(5) = 0.99, p = 0.37. With the code's bounds, which the criterion fixes, this test cannot detect a change in these data. On dF/F as a fraction, the scale used here, b is at most 10 where the mean response of 5 of the 7 animals is 13 to 45; in three animals delta AIC is exactly -6.0, the changing model having collapsed to the constant one, and in HJ-FP-M3 the changing model fits worse than the constant one (-168.5). On dF/F in percent, the stored scale, the bounded model cannot fit the data at all (delta AIC -272 to -817 in every animal), so that variant does not test the scale. Fitted without bounds, the changing model fits better in every animal (delta AIC 14 to 38, t = 7.84, p = 2.3e-4). The outcome rests on the assumption that the paper's fit was bounded as here; the paper's t of -0.86, close to -0.73, suggests it was. The scale the code's bounds were set for is asked of the authors (finding).
figs10b-no-omission-response¶
After the delay change, there is no suppression of dopamine at the old reward time (3 s) in any animal, whereas omission of the reward in extinction does suppress dopamine at the reward time (Test 5).
rows = []
for m in jeong.PAV_MICE:
d3, dx = jeong.OMISSION_DAYS[m]
rows.append({"subject": m, "exp3_methods": np.nansum(pav[m][d3]["omit3_base"]), "ext_methods": np.nansum(pav[m][dx]["omit9_base"]),
"exp3_code": np.nansum(pav[m][d3]["omit3_code"]), "ext_code": np.nansum(pav[m][dx]["omit9_code"]),
"trials": (len(pav[m][d3]["da"]), len(pav[m][dx]["da"])),
"incomplete": (int(pav[m][d3]["omit3_incomplete"].sum()), int(pav[m][dx]["omit9_incomplete"].sum())),
"exp3_strict": np.nansum(pav[m][d3]["omit3_base"][~pav[m][d3]["omit3_incomplete"]]),
"ext_strict": np.nansum(pav[m][dx]["omit9_base"][~pav[m][dx]["omit9_incomplete"]])})
om = pd.DataFrame(rows)
print(om.round(2).to_string(index=False))
for v in ["methods", "code", "strict"]:
e3, ex = om[f"exp3_{v}"], om[f"ext_{v}"]
print(f"\n{v}: Exp 3 at or above 0 in {(e3 >= 0).sum()} of 7; extinction below 0 in {(ex < 0).sum()} of 7, mean {ex.mean():.2f}")
show_t(" paired t, Exp 3 against extinction", stats.ttest_rel(e3, ex))
print("strict: the Methods' measure without trials with a NaN dF/F sample from -1 s to the end of the window")
print("Paper: t = 3.37, p = 0.0150")
subject exp3_methods ext_methods exp3_code ext_code trials incomplete exp3_strict ext_strict HJ-FP-M2 54.12 -213.63 109.75 -21.06 (50, 50) (1, 0) 56.51 -213.63 HJ-FP-M3 488.14 -49.21 54.95 -158.18 (50, 50) (0, 1) 488.14 -33.73 HJ-FP-M4 254.28 -61.07 33.11 -34.96 (50, 50) (0, 0) 254.28 -61.07 HJ-FP-F1 1.96 31.05 18.08 9.77 (50, 50) (0, 0) 1.96 31.05 HJ-FP-F2 330.01 -276.82 118.54 -114.99 (50, 50) (0, 0) 330.01 -276.82 HJ-FP-M6 128.86 -23.33 57.87 -3.58 (50, 50) (0, 0) 128.86 -23.33 HJ-FP-M7 35.05 -17.72 59.09 2.77 (50, 50) (0, 4) 35.05 -7.54 methods: Exp 3 at or above 0 in 7 of 7; extinction below 0 in 6 of 7, mean -87.25 paired t, Exp 3 against extinction: t = 3.033, p = 0.02301, df = 6 code: Exp 3 at or above 0 in 7 of 7; extinction below 0 in 5 of 7, mean -45.75 paired t, Exp 3 against extinction: t = 3.419, p = 0.01417, df = 6 strict: Exp 3 at or above 0 in 7 of 7; extinction below 0 in 6 of 7, mean -83.58 paired t, Exp 3 against extinction: t = 3.012, p = 0.02365, df = 6 strict: the Methods' measure without trials with a NaN dF/F sample from -1 s to the end of the window Paper: t = 3.37, p = 0.0150
Result: reproduced. The cumulative response at the old reward time on the first 8 s session is at or above 0 in all 7 animals, and the extinction omission response is below 0 on average (Methods' measure, mean -87.25, negative in 6 of 7; code's, -45.75, negative in 5 of 7). The paired test gives t = 3.03, p = 0.023 (Methods) and t = 3.42, p = 0.014 (code; paper t = 3.37, p = 0.0150). The extinction response is not negative in every animal, which the criterion's reading (b) allows. Leaving out the 0 to 4 trials per session with a NaN dF/F sample in the windows gives t = 3.01, p = 0.024, with the same signs.
fig4kl-extinction-abrupt-behavior¶
Behavioral learning of extinction is abrupt, much more so than the change of the dopamine cue response.
tab_c3 = abrupt_table(2)
report_abrupt(tab_c3)
print("\nVariant: trials with a NaN dF/F sample from -1 to 2 s left out")
report_abrupt(abrupt_table(2, drop_incomplete=True), brief=True)
print("\nPaper: abruptness t(6) = 5.67, p = 0.0013; change trial t(6) = -2.40, p = 0.0531")
Lick measure: Methods 0-2 s subject trials_before trials_used incomplete lick_abrupt da_abrupt lick_change da_change HJ-FP-M2 50 125 0 0.669 0.104 25 47 HJ-FP-M3 50 125 2 0.135 0.153 18 46 HJ-FP-M4 50 125 0 0.907 0.086 10 44 HJ-FP-F1 50 125 0 2.528 0.130 9 67 HJ-FP-F2 50 125 1 1.657 0.201 12 56 HJ-FP-M6 50 125 1 0.680 0.111 93 70 HJ-FP-M7 50 125 4 1.023 0.087 41 77 abruptness lick 1.086 against dopamine 0.125; change trial lick 29.7 against dopamine 58.1 paired t, abruptness lick minus dopamine: t = 3.283, p = 0.01676, df = 6 paired t, change trial lick minus dopamine: t = -2.953, p = 0.0255, df = 6 Lick measure: code cue+trace subject trials_before trials_used incomplete lick_abrupt da_abrupt lick_change da_change HJ-FP-M2 50 125 0 0.543 0.104 50 47 HJ-FP-M3 50 125 2 0.392 0.153 46 46 HJ-FP-M4 50 125 0 0.428 0.086 46 44 HJ-FP-F1 50 125 0 0.476 0.130 47 67 HJ-FP-F2 50 125 1 0.664 0.201 29 56 HJ-FP-M6 50 125 1 0.322 0.111 30 70 HJ-FP-M7 50 125 4 0.474 0.087 49 77 abruptness lick 0.471 against dopamine 0.125; change trial lick 42.4 against dopamine 58.1 paired t, abruptness lick minus dopamine: t = 9.726, p = 6.786e-05, df = 6 paired t, change trial lick minus dopamine: t = -2.402, p = 0.05314, df = 6 Variant: trials with a NaN dF/F sample from -1 to 2 s left out Lick measure: Methods 0-2 s abruptness lick 1.092 against dopamine 0.121; change trial lick 29.4 against dopamine 58.0 paired t, abruptness lick minus dopamine: t = 3.431, p = 0.01395, df = 6 paired t, change trial lick minus dopamine: t = -2.973, p = 0.02487, df = 6 Lick measure: code cue+trace abruptness lick 0.478 against dopamine 0.121; change trial lick 42.0 against dopamine 58.0 paired t, abruptness lick minus dopamine: t = 9.387, p = 8.299e-05, df = 6 paired t, change trial lick minus dopamine: t = -2.524, p = 0.04505, df = 6 Paper: abruptness t(6) = 5.67, p = 0.0013; change trial t(6) = -2.40, p = 0.0531
Result: reproduced. Lick abruptness exceeds dopamine abruptness with both lick measures (Methods t = 3.28, p = 0.017; code t = 9.73, p = 6.8e-5; paper t = 5.67). The change trial comparison, not part of the finding, gives t = -2.95 (p = 0.026) and t = -2.40 (p = 0.0531, the paper's values). Without the 0 to 4 trials per animal with a NaN dF/F sample from -1 to 2 s, the abruptness t values are 3.43 and 9.39. The code variant's abruptness t (9.7) is larger than the paper's 5.67 while its change-trial t equals the paper's; the cause is not found here.
fig4jl-extinction-cue-response-positive¶
The dopamine cue response remains significantly positive well after animals stop licking to the cue in extinction (Test 6).
For reading (b) the trials after the lick change trial are those from the change trial (the change-trial index plus 4, the first point of the distance) to the end of the extinction sessions, with the code's lick measure, whose change trial matches the paper's test.
rows = []
for m in jeong.PAV_MICE:
o = cumsum_condition(m, 2)
st, en = jeong.CUMSUM_START[m][2], jeong.CUMSUM_END[m][2]
ext_days = [d for d in o["days"] if st <= d <= en]
last = pav[m][ext_days[-1]]["da"]
ch = o["lick"]["code cue+trace"][1] + 4
ch_m = o["lick"]["Methods 0-2 s"][1] + 4
after = o["da"][ch:]
last_inc = pav[m][ext_days[-1]]["da_incomplete"]
tl = stats.ttest_1samp(last[np.isfinite(last)], 0)
rows.append({"subject": m, "extinction_days": ext_days, "last_mean": np.nanmean(last), "last_t": tl.statistic,
"last_p": tl.pvalue, "last_incomplete": int(last_inc.sum()),
"last_mean_strict": np.nanmean(last[~last_inc]), "change_trial": ch, "after_trials": int(np.isfinite(after).sum()),
"after_mean": np.nanmean(after), "change_trial_methods": ch_m, "after_mean_methods": np.nanmean(o["da"][ch_m:])})
ej = pd.DataFrame(rows)
print(ej.round(4).to_string(index=False))
show_t("(a) last extinction session, mean over animals against 0", stats.ttest_1samp(ej.last_mean, 0))
show_t("(b) trials after the lick change trial (code's lick measure), mean over animals against 0", stats.ttest_1samp(ej.after_mean, 0))
show_t("(b) trials after the lick change trial (Methods' lick measure), mean over animals against 0", stats.ttest_1samp(ej.after_mean_methods, 0))
show_t("(a) without trials with a NaN dF/F sample from -1 to 2 s", stats.ttest_1samp(ej.last_mean_strict, 0))
print(f"(c) animals positive at p < 0.05 in the last session: {((ej.last_t > 0) & (ej.last_p < 0.05)).sum()} of 7")
subject extinction_days last_mean last_t last_p last_incomplete last_mean_strict change_trial after_trials after_mean change_trial_methods after_mean_methods HJ-FP-M2 [70, 71, 72] 8.2013 4.6300 0.0000 1 8.3686 54 146 13.2583 29 15.5985 HJ-FP-M3 [61, 62, 63] 12.8625 9.9335 0.0000 0 12.8625 50 150 25.1021 22 30.6391 HJ-FP-M4 [66, 67] 16.0718 15.4202 0.0000 1 16.0301 50 100 16.2842 14 18.6410 HJ-FP-F1 [42, 43] 1.7793 4.2947 0.0001 0 1.7793 51 99 2.1717 13 2.5920 HJ-FP-F2 [46, 47] 14.2016 8.3019 0.0000 0 14.2016 33 117 19.8187 16 22.5886 HJ-FP-M6 [47, 48, 50] 2.0106 4.9324 0.0000 0 2.0106 34 144 3.8696 97 2.7080 HJ-FP-M7 [43, 44] 9.7136 12.4761 0.0000 0 9.7136 53 97 10.4772 45 10.5124 (a) last extinction session, mean over animals against 0: t = 4.314, p = 0.005016, df = 6 (b) trials after the lick change trial (code's lick measure), mean over animals against 0: t = 4.156, p = 0.005971, df = 6 (b) trials after the lick change trial (Methods' lick measure), mean over animals against 0: t = 3.778, p = 0.009202, df = 6 (a) without trials with a NaN dF/F sample from -1 to 2 s: t = 4.332, p = 0.004916, df = 6 (c) animals positive at p < 0.05 in the last session: 7 of 7
Result: reproduced. The CS+ response in the last extinction session is positive across animals (t(6) = 4.31, p = 0.0050; 4.33 without the 2 trials with a NaN dF/F sample), after the lick change trial (reading b, all 7 animal means positive, t(6) = 4.16, p = 0.006 with the code's lick measure and t(6) = 3.78, p = 0.0092 with the Methods'), and within each animal's last session (7 of 7 at p < 0.05).
fig4p-background-faster-decay¶
When unpredicted rewards are added during the ITI, the dopamine cue response stays positive but decreases faster over trials than during extinction (Test 7).
def bgd_vs_ext(start, end, before, P_key="da", strict=False):
# cumsum_bgdrw_vs_extinction.m for one animal and condition: 140 trials from 50 before the
# condition starts; readings (a) to (c) and the post-change mean.
rows = []
for m in jeong.PAV_MICE:
row = {"subject": m}
for c, name in enumerate(["bgd", "ext"]):
st, en = start[m][c], end[m][c]
days = [d for d in sorted(pav[m]) if st - before <= d <= en]
n = sum(len(pav[m][d]["da"]) for d in days[:before])
sel = np.arange(140) + n - 50
data = np.concatenate([pav[m][d]["da"] for d in days])[sel]
inc = np.concatenate([pav[m][d]["da_incomplete"] for d in days])[sel]
if strict:
data = np.where(inc, np.nan, data)
nc = np.nancumsum(data) / np.nansum(data[:50])
row.update({f"{name}_days": f"{days[0]}-{days[-1]} ({len(days)})", f"{name}_incomplete": int(inc.sum()),
f"{name}_pre_sum": np.nansum(data[:50]), f"{name}_a": nc[-1], f"{name}_b": nc[50:].mean(),
f"{name}_c": np.nanmean(data[50:]) / np.nanmean(data[:50]), f"{name}_post": np.nanmean(data[50:])})
rows.append(row)
return pd.DataFrame(rows)
def report_bgd(bg):
for r in "abc":
res = stats.ttest_rel(bg[f"bgd_{r}"], bg[f"ext_{r}"])
print(f" reading ({r}): background below extinction in {(bg[f'bgd_{r}'] < bg[f'ext_{r}']).sum()} of 7; "
f"paired t = {res.statistic:.3f}, p = {res.pvalue:.4g}")
bg = bgd_vs_ext(jeong.BGD_START, jeong.BGD_END, 2)
print(bg.round(3).to_string(index=False))
report_bgd(bg)
show_t("background post-change response against 0", stats.ttest_1samp(bg.bgd_post, 0))
no_m4 = bg[bg.subject != "HJ-FP-M4"]
show_t("variant without HJ-FP-M4 (first background session missing), reading (a)", stats.ttest_rel(no_m4.bgd_a, no_m4.ext_a))
# Readings (a) and (c) are one affine function of the other when no trial is NaN:
# (a) = 1 + (90 / 50) x (c), so they give the same t.
print("largest |(a) - (1 + 1.8 (c))| over animals and conditions:",
float(np.max(np.abs(np.r_[bg.bgd_a - 1 - 1.8 * bg.bgd_c, bg.ext_a - 1 - 1.8 * bg.ext_c]))))
print("\nVariant: trials with a NaN dF/F sample from -1 to 2 s left out")
report_bgd(bgd_vs_ext(jeong.BGD_START, jeong.BGD_END, 2, strict=True))
# Variant of the sessions: those of cumsum_pavlovian.m, conditions 4 (background) and 3
# (extinction), with one session before each; HJ-FP-M2's background range starts at Day38 there
# against Day33 in cumsum_bgdrw_vs_extinction.m.
start_cp = {m: [jeong.CUMSUM_START[m][3], jeong.CUMSUM_START[m][2]] for m in jeong.PAV_MICE}
end_cp = {m: [jeong.CUMSUM_END[m][3], jeong.CUMSUM_END[m][2]] for m in jeong.PAV_MICE}
bg_cp = bgd_vs_ext(start_cp, end_cp, 1)
print("\nVariant: sessions of cumsum_pavlovian.m, one session before each condition")
print(bg_cp[["subject", "bgd_days", "ext_days", "bgd_pre_sum", "ext_pre_sum", "bgd_a", "ext_a", "bgd_b", "ext_b"]].round(3).to_string(index=False))
report_bgd(bg_cp)
print("Paper: t(6) = -3.51, p = 0.0126")
subject bgd_days bgd_incomplete bgd_pre_sum bgd_a bgd_b bgd_c bgd_post ext_days ext_incomplete ext_pre_sum ext_a ext_b ext_c ext_post HJ-FP-M2 31-43 (13) 0 704.902 1.476 1.444 0.264 3.728 68-72 (5) 0 1428.380 2.055 1.590 0.586 16.751 HJ-FP-M3 26-38 (12) 1 2213.951 1.367 1.256 0.204 9.023 59-63 (5) 2 3200.482 1.914 1.509 0.508 32.507 HJ-FP-M4 36-47 (11) 1 1109.466 1.146 1.112 0.081 1.806 64-67 (4) 1 1185.709 2.290 1.632 0.716 16.990 HJ-FP-F1 21-27 (7) 1 119.383 0.942 1.079 -0.032 -0.077 40-43 (4) 0 167.120 2.209 1.712 0.672 2.245 HJ-FP-F2 26-33 (8) 1 1130.466 1.875 1.517 0.486 10.989 44-47 (4) 1 2000.160 1.742 1.437 0.412 16.494 HJ-FP-M6 21-30 (10) 1 288.050 2.445 1.826 0.803 4.625 45-50 (5) 1 330.093 2.158 1.617 0.643 4.248 HJ-FP-M7 22-29 (8) 3 457.992 1.925 1.502 0.514 4.708 41-44 (4) 4 642.000 2.506 1.815 0.837 10.745 reading (a): background below extinction in 5 of 7; paired t = -2.403, p = 0.05307 reading (b): background below extinction in 5 of 7; paired t = -1.968, p = 0.09668 reading (c): background below extinction in 5 of 7; paired t = -2.403, p = 0.05307 background post-change response against 0: t = 3.399, p = 0.01452, df = 6 variant without HJ-FP-M4 (first background session missing), reading (a): t = -1.850, p = 0.1235, df = 5 largest |(a) - (1 + 1.8 (c))| over animals and conditions: 1.887379141862766e-15 Variant: trials with a NaN dF/F sample from -1 to 2 s left out reading (a): background below extinction in 5 of 7; paired t = -2.475, p = 0.04812 reading (b): background below extinction in 5 of 7; paired t = -2.037, p = 0.08778 reading (c): background below extinction in 5 of 7; paired t = -2.493, p = 0.04694 Variant: sessions of cumsum_pavlovian.m, one session before each condition subject bgd_days ext_days bgd_pre_sum ext_pre_sum bgd_a ext_a bgd_b ext_b HJ-FP-M2 37-43 (7) 69-72 (4) -492.936 1428.380 2.746 2.055 1.816 1.590 HJ-FP-M3 27-38 (11) 60-63 (4) 2213.951 3200.482 1.367 1.914 1.256 1.509 HJ-FP-M4 37-47 (10) 65-67 (3) 1109.466 1185.709 1.146 2.290 1.112 1.632 HJ-FP-F1 22-27 (6) 41-43 (3) 119.383 167.120 0.942 2.209 1.079 1.712 HJ-FP-F2 27-33 (7) 45-47 (3) 1130.466 2000.160 1.875 1.742 1.517 1.437 HJ-FP-M6 22-30 (9) 46-50 (4) 37.318 330.093 12.154 2.158 7.372 1.617 HJ-FP-M7 23-29 (7) 42-44 (3) 457.992 642.000 1.925 2.506 1.502 1.815 reading (a): background below extinction in 4 of 7; paired t = 0.687, p = 0.518 reading (b): background below extinction in 4 of 7; paired t = 0.718, p = 0.4995 reading (c): background below extinction in 4 of 7; paired t = 0.687, p = 0.518 Paper: t(6) = -3.51, p = 0.0126
Result: not reproduced. The background reward cue response stays positive after the change (mean over the 90 trials, t(6) = 3.40, p = 0.015). It decays further than in extinction on average, but not significantly under any reading: paired t = -2.40, p = 0.053 under readings (a) and (c), t = -1.97, p = 0.097 under (b) (paper t = -3.51, p = 0.0126); background is below extinction in 5 of 7 animals. Readings (a) and (c) are one affine function of the other, (a) = 1 + 1.8 (c), so only two independent readings were tested. The criterion asks for background below extinction at p < 0.05 under every reading, which no reading meets. HJ-FP-M4's first background reward session (Day38) is not in the dandiset, so its 90 post-change trials start a session later; without that animal reading (a) gives t(5) = -1.85, p = 0.12. Leaving out the 0 to 4 trials per animal and condition with a NaN dF/F sample, (a) and (c) reach p = 0.048 and 0.047, but (b) stays at p = 0.088. With the sessions of cumsum_pavlovian.m (conditions 4 and 3, one session before each), no reading shows a difference (t = 0.69 and 0.72). There HJ-FP-M6's background sum over the 50 trials before the change is 37.3 (288.1 with the sessions of cumsum_bgdrw_vs_extinction.m), so its reading (a) is 12.15, and HJ-FP-M2's is -492.9, as its range starts at Day38 against Day33.
fig6b-no-backpropagation¶
During initial cue-reward learning, dopamine activity does not shift backward from just before reward to the cue; the late (pre-reward) response does not differ from the early (cue) response in the first 100 trials (Test 9).
rows = []
for m in jeong.PAV_MICE:
st, en = jeong.CUMSUM_START[m][0], jeong.CUMSUM_END[m][0]
days = [d for d in sorted(pav[m]) if st <= d <= en]
e = np.concatenate([pav[m][d]["early"] for d in days])
l = np.concatenate([pav[m][d]["late"] for d in days])
p = np.concatenate([pav[m][d]["pre"] for d in days])
inc = np.concatenate([pav[m][d]["unsmoothed_incomplete"] for d in days])
eb, lb = e - p, l - p
s_ = lambda v: np.where(inc, np.nan, v)
rows.append({"subject": m, "trials": len(e), "incomplete": int(inc[:100].sum() + inc[350:400].sum()),
"methods": np.nanmean(lb[:100] - eb[:100]) / np.nanmean(eb[350:400]),
"code": np.nanmean(l[:100] - e[:100]) / np.nanmean(e[350:400]),
"methods_strict": np.nanmean(s_(lb)[:100] - s_(eb)[:100]) / np.nanmean(s_(eb)[350:400])})
bp = pd.DataFrame(rows)
print(bp.round(3).to_string(index=False))
for v in ["methods", "code", "methods_strict"]:
print(f"{v}: mean {bp[v].mean():.3f}"); show_t(" one-sample t of late minus early against 0", stats.ttest_1samp(bp[v], 0))
print("methods_strict leaves out trials with a NaN unsmoothed dF/F sample from -1 to 3 s")
print("Paper: t = -0.16, p = 0.8771")
subject trials incomplete methods code methods_strict HJ-FP-M2 600 0 0.169 0.191 0.169 HJ-FP-M3 550 1 -0.045 -0.047 -0.043 HJ-FP-M4 900 0 -0.088 -0.088 -0.088 HJ-FP-F1 600 0 -0.046 -0.048 -0.046 HJ-FP-F2 473 3 -0.022 -0.022 -0.023 HJ-FP-M6 600 0 0.040 0.043 0.040 HJ-FP-M7 673 4 -0.061 -0.072 -0.059 methods: mean -0.008 one-sample t of late minus early against 0: t = -0.229, p = 0.8266, df = 6 code: mean -0.006 one-sample t of late minus early against 0: t = -0.163, p = 0.876, df = 6 methods_strict: mean -0.007 one-sample t of late minus early against 0: t = -0.216, p = 0.8359, df = 6 methods_strict leaves out trials with a NaN unsmoothed dF/F sample from -1 to 3 s Paper: t = -0.16, p = 0.8771
Result: reproduced. Over the first 100 CS+ trials the normalized late minus early response is -0.008 with the Methods' baseline (t = -0.23, p = 0.83) and -0.006 without it (t = -0.16, p = 0.876, the paper's values). No difference is found. Leaving out the 0 to 4 trials per animal with a NaN dF/F sample from -1 to 3 s gives t = -0.22, p = 0.84.
Experiment 6: trial-less conditioning¶
Sessions from learnedday to endday of intermediate_rsp_trialless.m; jeong.trialless_pairs
follows that script and the Methods.
peak_readings = {(0.5, "last"): 1, (0.5, "mean"): 1, (1.0, "last"): 1, (1.0, "mean"): 1}
rows = []
for m in jeong.PAV_MICE:
for s in subject_sessions(m, "Poisson"):
if jeong.TRIALLESS_LEARNED[m] <= s["day"] <= jeong.TRIALLESS_END[m]:
rows.append({"subject": m, "day": s["day"], **jeong.trialless_pairs(s, peak_readings)})
tl = pd.DataFrame(rows)
print(tl[["subject", "day", "n_cues", "n_pairs", "delay", "cue_prev", "cue_inter", "rew_prev", "rew_inter"]].round(3).to_string(index=False))
tla = tl.groupby("subject").mean(numeric_only=True)
print(f"Cues of pairs with a NaN dF/F sample from -0.5 to 0.5 s: {tl.n_incomplete_cue.sum()} of {2 * tl.n_pairs.sum()}; "
f"their rewards with one from 0 to 0.5 s after the first lick: {tl.n_incomplete_rew.sum()}")
# Variant: those cues and rewards left out of the means.
tl_strict = pd.DataFrame([{"subject": m, **jeong.trialless_pairs(s, peak_readings, drop_incomplete=True)}
for m in jeong.PAV_MICE for s in subject_sessions(m, "Poisson")
if jeong.TRIALLESS_LEARNED[m] <= s["day"] <= jeong.TRIALLESS_END[m]])
tla_strict = tl_strict.groupby("subject").mean(numeric_only=True)
subject day n_cues n_pairs delay cue_prev cue_inter rew_prev rew_inter HJ-FP-M2 77 100 4 3.0 24.242 14.701 5.812 19.980 HJ-FP-M2 78 100 8 3.0 15.669 4.446 3.142 9.095 HJ-FP-M3 67 100 3 3.0 42.848 22.302 28.470 36.174 HJ-FP-M3 68 100 4 3.0 45.361 21.746 6.140 23.780 HJ-FP-M4 70 100 8 3.0 7.447 3.824 23.779 23.713 HJ-FP-M4 71 81 3 3.0 10.644 2.929 27.489 28.158 HJ-FP-M4 72 100 4 3.0 14.185 7.393 15.985 27.952 HJ-FP-F1 46 100 5 3.0 0.825 -1.162 7.089 6.079 HJ-FP-F1 47 100 4 3.0 4.435 1.591 5.620 8.447 HJ-FP-F1 48 100 5 3.0 4.465 1.402 6.276 9.651 HJ-FP-F1 49 100 4 3.0 7.830 0.275 4.585 11.938 HJ-FP-F2 51 100 4 3.0 19.556 9.422 32.733 34.930 HJ-FP-F2 52 100 4 3.0 33.415 17.898 23.596 42.298 HJ-FP-F2 53 100 4 3.0 37.792 18.147 24.590 41.904 HJ-FP-F2 54 100 4 3.0 29.811 14.588 12.551 43.367 HJ-FP-M6 53 100 4 3.0 7.766 4.296 6.844 7.427 HJ-FP-M6 54 100 6 3.0 4.718 3.115 3.715 4.071 HJ-FP-M6 55 100 4 3.0 4.826 4.316 3.441 3.551 HJ-FP-M7 48 100 4 3.0 1.875 0.267 0.788 2.059 HJ-FP-M7 49 100 4 3.0 2.996 1.236 -0.095 1.148 HJ-FP-M7 50 100 5 3.0 1.729 0.676 1.497 2.086 Cues of pairs with a NaN dF/F sample from -0.5 to 0.5 s: 0 of 190; their rewards with one from 0 to 0.5 s after the first lick: 0
fig5d-intermediate-cue-positive¶
In the trial-less task, the dopamine response to an intermediate cue (one arriving before the previous cue's reward) is positive but weaker than the response to the previous cue (Test 8).
ratio_cue = tla.cue_inter / tla.cue_prev
print("Per animal ratio intermediate over previous cue response:", ratio_cue.round(3).to_dict())
print(f"mean ratio {ratio_cue.mean():.3f}; pairs per animal {tl.groupby('subject').n_pairs.sum().to_dict()}")
show_t("one-sample t against 0", stats.ttest_1samp(ratio_cue, 0))
show_t("one-sample t against 1", stats.ttest_1samp(ratio_cue, 1))
print("Paper: t(6) = 6.64, p = 5.6e-4")
r_strict = tla_strict.cue_inter / tla_strict.cue_prev
print(f"Without cues with a NaN dF/F sample in their windows: mean ratio {r_strict.mean():.3f}")
show_t(" one-sample t against 0", stats.ttest_1samp(r_strict, 0))
Per animal ratio intermediate over previous cue response: {'HJ-FP-F1': 0.12, 'HJ-FP-F2': 0.498, 'HJ-FP-M2': 0.48, 'HJ-FP-M3': 0.499, 'HJ-FP-M4': 0.438, 'HJ-FP-M6': 0.677, 'HJ-FP-M7': 0.33}
mean ratio 0.435; pairs per animal {'HJ-FP-F1': 18, 'HJ-FP-F2': 16, 'HJ-FP-M2': 12, 'HJ-FP-M3': 7, 'HJ-FP-M4': 15, 'HJ-FP-M6': 14, 'HJ-FP-M7': 13}
one-sample t against 0: t = 6.656, p = 0.0005558, df = 6
one-sample t against 1: t = -8.655, p = 0.0001311, df = 6
Paper: t(6) = 6.64, p = 5.6e-4
Without cues with a NaN dF/F sample in their windows: mean ratio 0.435
one-sample t against 0: t = 6.656, p = 0.0005558, df = 6
Result: reproduced. The intermediate cue response is positive relative to the previous cue's baseline, at 0.435 of the previous cue's response (0.12 to 0.68 per animal), t(6) = 6.66, p = 5.6e-4 against 0 (paper t = 6.64), and below 1. No cue or reward of the pairs has a NaN dF/F sample in its windows (0 of 190 cues), so leaving such values out changes nothing.
fig5d-intermediate-reward-larger¶
The dopamine response to the reward paired with an intermediate cue is larger than the response to the reward of the previous cue (the ratio is above 1).
ratio_rw = tla.rew_inter / tla.rew_prev
print("Per animal ratio intermediate over previous reward response:", ratio_rw.round(3).to_dict())
print(f"mean ratio {ratio_rw.mean():.3f}")
show_t("one-sample t against 1", stats.ttest_1samp(ratio_rw, 1))
print("Paper: t(6) = 2.95, p = 0.0256")
r_strict = tla_strict.rew_inter / tla_strict.rew_prev
print(f"Without rewards with a NaN dF/F sample in their windows: mean ratio {r_strict.mean():.3f}")
show_t(" one-sample t against 1", stats.ttest_1samp(r_strict, 1))
Per animal ratio intermediate over previous reward response: {'HJ-FP-F1': 1.532, 'HJ-FP-F2': 1.739, 'HJ-FP-M2': 3.247, 'HJ-FP-M3': 1.732, 'HJ-FP-M4': 1.187, 'HJ-FP-M6': 1.075, 'HJ-FP-M7': 2.418}
mean ratio 1.847
one-sample t against 1: t = 2.959, p = 0.02531, df = 6
Paper: t(6) = 2.95, p = 0.0256
Without rewards with a NaN dF/F sample in their windows: mean ratio 1.847
one-sample t against 1: t = 2.959, p = 0.02531, df = 6
Result: reproduced. The mean ratio is 1.85 (1.08 to 3.25), t(6) = 2.96, p = 0.025 against 1 (paper t = 2.95, p = 0.0256). No cue or reward of the pairs has a NaN dF/F sample in its windows (0 of 190 cues), so leaving such values out changes nothing.
figs11de-intermediate-cue-local-baseline¶
With the baseline taken immediately before the intermediate cue and the response measured as a peak, the intermediate cue response is still positive.
for pw, bw in peak_readings:
r = tla[f"peak_{pw}_{bw}_inter"] / tla[f"peak_{pw}_{bw}_prev"]
print(f"peak over 0 to {pw} s, baseline {'last sample' if bw == 'last' else 'mean of 0.5 s'} before the cue: "
f"mean ratio {r.mean():.3f}, min {r.min():.3f}")
show_t(" one-sample t against 0", stats.ttest_1samp(r, 0))
print("Paper: t = 10.42, p = 4.6e-5")
print("Without cues with a NaN dF/F sample in their windows:")
for pw, bw in peak_readings:
r = tla_strict[f"peak_{pw}_{bw}_inter"] / tla_strict[f"peak_{pw}_{bw}_prev"]
res = stats.ttest_1samp(r, 0)
print(f" peak 0 to {pw} s, baseline {bw}: mean ratio {r.mean():.3f}, t = {res.statistic:.3f}, p = {res.pvalue:.3g}")
peak over 0 to 0.5 s, baseline last sample before the cue: mean ratio 0.383, min 0.276 one-sample t against 0: t = 8.355, p = 0.0001598, df = 6 peak over 0 to 0.5 s, baseline mean of 0.5 s before the cue: mean ratio 0.253, min 0.135 one-sample t against 0: t = 6.682, p = 0.0005442, df = 6 peak over 0 to 1.0 s, baseline last sample before the cue: mean ratio 0.448, min 0.282 one-sample t against 0: t = 5.851, p = 0.0011, df = 6 peak over 0 to 1.0 s, baseline mean of 0.5 s before the cue: mean ratio 0.313, min 0.170 one-sample t against 0: t = 5.215, p = 0.001986, df = 6 Paper: t = 10.42, p = 4.6e-5 Without cues with a NaN dF/F sample in their windows: peak 0 to 0.5 s, baseline last: mean ratio 0.383, t = 8.355, p = 0.00016 peak 0 to 0.5 s, baseline mean: mean ratio 0.253, t = 6.682, p = 0.000544 peak 0 to 1.0 s, baseline last: mean ratio 0.448, t = 5.851, p = 0.0011 peak 0 to 1.0 s, baseline mean: mean ratio 0.313, t = 5.215, p = 0.00199
Result: reproduced. Under all four readings the intermediate cue's peak response over the previous cue's is positive (mean ratio 0.25 to 0.45, t = 5.21 to 8.35, p <= 0.002; paper t = 10.42). No cue or reward of the pairs has a NaN dF/F sample in its windows (0 of 190 cues), so leaving such values out changes nothing.
Experiment 7: sequential conditioning with inhibition¶
Per animal the first session in Day order is the random reward session and the rest are conditioning sessions (checked in methods-session-counts). DAT-Cre are datHT-M2 to M7 (M8 excluded, Fig S12D) and controls datWT-F1 and WT-F1, F2, F3, M1, M2, M3.
e7 = {}
for m in jeong.DAT_MICE + [jeong.DAT_EXCLUDED] + jeong.CONTROL_MICE:
ss = subject_sessions(m)
e7[m] = {"random": jeong.random_session_reward(ss[0]), "cond": [jeong.seq_session(s) for s in ss[1:]]}
cues = [(len(c["cs1"]), len(c["cs2"])) for c in e7[m]["cond"]]
print(m, "conditioning sessions", len(cues), "CS1/CS2 per session", cues)
# auc_python gives NaN when a window holds no dF/F sample; the session means below use nanmean,
# so such trials drop out. Count them per measure, for sessions where any are NaN.
nan_rows = []
for m in e7:
for i, c in enumerate(e7[m]["cond"]):
row = {"subject": m, "conditioning_session": i + 1, "of": len(e7[m]["cond"]), "cs1": len(c["cs1"])}
row.update({k: int(np.isnan(c[k]).sum()) for k in ["cs1_15", "early", "late", "rw_15"]})
if any(row[k] for k in ["cs1_15", "early", "late", "rw_15"]):
nan_rows.append(row)
print("\nConditioning sessions with trials without a dF/F value (NaN AUC):")
print(pd.DataFrame(nan_rows).to_string(index=False))
print("HJ-FP-WT-stGtACR-M2, last session, CS1 trials (1-based) without a value:",
(np.flatnonzero(np.isnan(e7["HJ-FP-WT-stGtACR-M2"]["cond"][-1]["cs1_15"])) + 1).tolist())
print("Random reward sessions with a NaN reward response:", [m for m in e7 if not np.isfinite(e7[m]["random"]["rw_15"])])
HJ-FP-datHT-stGtACR-M2 conditioning sessions 5 CS1/CS2 per session [(50, 50), (50, 50), (50, 50), (50, 50), (50, 50)]
HJ-FP-datHT-stGtACR-M3 conditioning sessions 6 CS1/CS2 per session [(50, 50), (14, 14), (50, 50), (50, 50), (50, 50), (50, 50)]
HJ-FP-datHT-stGtACR-M4 conditioning sessions 6 CS1/CS2 per session [(50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50)]
HJ-FP-datHT-stGtACR-M5 conditioning sessions 8 CS1/CS2 per session [(50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50)]
HJ-FP-datHT-stGtACR-M6 conditioning sessions 7 CS1/CS2 per session [(50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50)]
HJ-FP-datHT-stGtACR-M7 conditioning sessions 12 CS1/CS2 per session [(50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50)] HJ-FP-datHT-stGtACR-M8 conditioning sessions 5 CS1/CS2 per session [(50, 50), (50, 50), (50, 50), (50, 50), (50, 50)]
HJ-FP-datWT-stGtACR-F1 conditioning sessions 7 CS1/CS2 per session [(50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50)] HJ-FP-WT-stGtACR-F1 conditioning sessions 5 CS1/CS2 per session [(50, 50), (50, 50), (50, 50), (50, 50), (50, 50)]
HJ-FP-WT-stGtACR-F2 conditioning sessions 5 CS1/CS2 per session [(50, 50), (50, 50), (50, 50), (50, 50), (50, 50)] HJ-FP-WT-stGtACR-F3 conditioning sessions 6 CS1/CS2 per session [(50, 50), (50, 50), (50, 50), (50, 50), (50, 50), (50, 50)]
HJ-FP-WT-stGtACR-M1 conditioning sessions 5 CS1/CS2 per session [(50, 50), (50, 50), (50, 50), (50, 50), (50, 50)] HJ-FP-WT-stGtACR-M2 conditioning sessions 5 CS1/CS2 per session [(50, 50), (49, 49), (50, 50), (50, 50), (50, 50)]
HJ-FP-WT-stGtACR-M3 conditioning sessions 5 CS1/CS2 per session [(50, 50), (50, 50), (50, 50), (50, 50), (25, 25)]
Conditioning sessions with trials without a dF/F value (NaN AUC):
subject conditioning_session of cs1 cs1_15 early late rw_15
HJ-FP-datHT-stGtACR-M3 4 6 50 0 0 0 1
HJ-FP-datHT-stGtACR-M5 6 8 50 0 0 0 12
HJ-FP-WT-stGtACR-M1 1 5 50 0 0 0 2
HJ-FP-WT-stGtACR-M2 5 5 50 17 17 17 18
HJ-FP-WT-stGtACR-M2, last session, CS1 trials (1-based) without a value: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]
Random reward sessions with a NaN reward response: []
In Experiment 7 a trial whose window holds no dF/F sample has no value (auc_python, as
calculate_auc), and the session means leave it out. HJ-FP-WT-stGtACR-M2's last conditioning
session (Day6) has no dF/F for its first 1117 s, so its CS1 and CS2 trials 1 to 17 and 18 of its
reward responses have no value. Elsewhere only reward responses are missing: 1, 12 and 2 in one
session each of datHT-M3, datHT-M5 and WT-M1 (in WT-M1's first session, so its day 1 reward
response, which normalises Fig 6I, is a mean over 48 trials). Of the trials Fig 6D uses (1 to 50 and 151 to
200), only trial 200 of WT-M2 falls in that session (its sessions hold 50, 49, 50, 50 and 50
trials), and it has no value. fig6i-learning-despite-inhibition
and fig6i-lower-cs1-asymptote take that session's mean over its other 33 trials, and the latter
shows the result without the session.
fig6d-cues-increase-together¶
In sequential conditioning without inhibition, the dopamine responses to CS2 and CS1 increase together early in learning (no difference in the first 50 trials) (Test 10).
rows = []
for m in jeong.CONTROL_MICE:
c = e7[m]["cond"]
e = np.concatenate([x["early"] for x in c]); l = np.concatenate([x["late"] for x in c])
p1 = np.concatenate([x["pre1"] for x in c]); p15 = np.concatenate([x["pre15"] for x in c])
row = {"subject": m, "trials": len(e)}
for name, b in [("code (no baseline)", 0 * e), ("Methods, 1 s baseline", p1), ("Methods, 1.5 s baseline", p15)]:
E, L = e - b, l - b
row[name] = np.nanmean(L[:50] - E[:50]) / np.nanmean(E[150:200])
row[name + " late ratio"] = np.nanmean(L[150:200]) / np.nanmean(E[150:200])
rows.append(row)
sq = pd.DataFrame(rows)
print(sq.round(3).to_string(index=False))
for name in ["code (no baseline)", "Methods, 1 s baseline", "Methods, 1.5 s baseline"]:
print(f"{name}: mean CS2 minus CS1 over the first 50 trials {sq[name].mean():.3f}")
show_t(" one-sample (paired) t against 0", stats.ttest_1samp(sq[name], 0))
print("Paper: t = -0.03, p = 0.9752")
subject trials code (no baseline) code (no baseline) late ratio Methods, 1 s baseline Methods, 1 s baseline late ratio Methods, 1.5 s baseline Methods, 1.5 s baseline late ratio HJ-FP-datWT-stGtACR-F1 350 -0.315 0.303 -0.283 0.372 -0.285 0.370 HJ-FP-WT-stGtACR-F1 250 -0.023 0.589 -0.023 0.593 -0.023 0.588 HJ-FP-WT-stGtACR-F2 250 0.024 0.702 0.023 0.711 0.023 0.709 HJ-FP-WT-stGtACR-F3 300 -0.102 0.407 -0.094 0.453 -0.095 0.450 HJ-FP-WT-stGtACR-M1 250 0.014 0.705 0.014 0.702 0.013 0.711 HJ-FP-WT-stGtACR-M2 249 -0.134 0.534 -0.129 0.552 -0.127 0.561 HJ-FP-WT-stGtACR-M3 225 0.515 0.604 0.480 0.631 0.486 0.626 code (no baseline): mean CS2 minus CS1 over the first 50 trials -0.003 one-sample (paired) t against 0: t = -0.032, p = 0.9752, df = 6 Methods, 1 s baseline: mean CS2 minus CS1 over the first 50 trials -0.002 one-sample (paired) t against 0: t = -0.021, p = 0.9841, df = 6 Methods, 1.5 s baseline: mean CS2 minus CS1 over the first 50 trials -0.001 one-sample (paired) t against 0: t = -0.010, p = 0.9921, df = 6 Paper: t = -0.03, p = 0.9752
Result: reproduced. Over the first 50 trials the normalized CS2 minus CS1 response does not differ from 0 under any reading: t = -0.03, p = 0.975 without a baseline (the paper's values), and p = 0.98 and 0.99 with the 1 s and 1.5 s baselines (means -0.003 to -0.001).
fig6d-cs2-below-cs1-late¶
Later in sequential conditioning, the CS2 response falls below the CS1 response (the two separate).
for name in ["code (no baseline)", "Methods, 1 s baseline", "Methods, 1.5 s baseline"]:
v = sq[name + " late ratio"]
print(f"{name}: CS2/CS1 over trials 151 to 200, mean {v.mean():.3f} (range {v.min():.3f} to {v.max():.3f})")
show_t(" one-sample t against 1", stats.ttest_1samp(v, 1))
code (no baseline): CS2/CS1 over trials 151 to 200, mean 0.549 (range 0.303 to 0.705) one-sample t against 1: t = -8.006, p = 0.0002026, df = 6 Methods, 1 s baseline: CS2/CS1 over trials 151 to 200, mean 0.573 (range 0.372 to 0.711) one-sample t against 1: t = -8.994, p = 0.0001056, df = 6 Methods, 1.5 s baseline: CS2/CS1 over trials 151 to 200, mean 0.574 (range 0.370 to 0.711) one-sample t against 1: t = -8.876, p = 0.0001138, df = 6
Result: reproduced. Over trials 151 to 200 the CS2 response is 0.55 of the CS1 response (0.30 to 0.71 per animal), t(6) = -8.01, p = 2.0e-4 against 1, and below 1 under the baseline readings too.
fig6g-inhibition-magnitude¶
Optogenetic inhibition of VTA dopamine neurons strongly suppresses NAcc dopamine during the CS2-to-reward delay in DAT-Cre mice (to about -0.6 times the reward response on day 1) and reduces their reward response, compared with controls.
The baseline is the mean dF/F over the 0.5 s before CS1 (claim), or before every code-15 event,
CS1 and CS2 (the code's reftimes, a variant). For a window of another length the baseline AUC is
scaled to that length, so that each response subtracts a baseline over an equal window (claims
header); this is a choice made here, and it only matters for the 1 s reward window.
rows = []
for m in jeong.DAT_MICE + jeong.CONTROL_MICE + [jeong.DAT_EXCLUDED]:
r, c = e7[m]["random"], e7[m]["cond"][0]
row = {"subject": m, "group": "DAT-Cre" if "datHT" in m else "control"}
for bname, b in [("cs1", c["base05_cs1"]), ("all15", c["base05_all15"])]:
inh = c["inh_bins"][1.0] - b
rw1, rw05 = c["rw1"] - 2 * b, c["rw05"] - b
row[f"{bname}: (a) AUC ratio, 1 s reward"] = inh / rw1
row[f"{bname}: (b) AUC ratio, 0.5 s reward"] = inh / rw05
row[f"{bname}: (c) level ratio, 1 s reward"] = (inh / 0.5) / rw1
row[f"{bname}: day 1 reward, 1 s"] = rw1 / (r["rw1"] - 2 * r["base05"])
row[f"{bname}: day 1 reward, 0.5 s"] = rw05 / (r["rw05"] - r["base05"])
rows.append(row)
inh = pd.DataFrame(rows)
print(inh.round(3).T.to_string())
main = inh[inh.subject != jeong.DAT_EXCLUDED]
for col in [c for c in inh.columns if ":" in c]:
d, w = main[main.group == "DAT-Cre"][col], main[main.group == "control"][col]
res = stats.ttest_ind(d, w)
print(f"{col:42s} DAT-Cre mean {d.mean():7.3f} control {w.mean():7.3f} t = {res.statistic:6.2f}, p = {res.pvalue:.3g}")
print("Paper: inhibition t = -5.74, p = 1.3e-4; day 1 reward t = -3.99, p = 0.0021; DAT-Cre about -0.6")
0 1 2 3 4 5 6 7 8 9 10 11 12 13 subject HJ-FP-datHT-stGtACR-M2 HJ-FP-datHT-stGtACR-M3 HJ-FP-datHT-stGtACR-M4 HJ-FP-datHT-stGtACR-M5 HJ-FP-datHT-stGtACR-M6 HJ-FP-datHT-stGtACR-M7 HJ-FP-datWT-stGtACR-F1 HJ-FP-WT-stGtACR-F1 HJ-FP-WT-stGtACR-F2 HJ-FP-WT-stGtACR-F3 HJ-FP-WT-stGtACR-M1 HJ-FP-WT-stGtACR-M2 HJ-FP-WT-stGtACR-M3 HJ-FP-datHT-stGtACR-M8 group DAT-Cre DAT-Cre DAT-Cre DAT-Cre DAT-Cre DAT-Cre control control control control control control control DAT-Cre cs1: (a) AUC ratio, 1 s reward -0.24 -0.239 -0.309 -0.302 -0.117 -0.275 0.012 0.016 -0.025 0.006 0.084 -0.02 0.051 2.072 cs1: (b) AUC ratio, 0.5 s reward -0.515 -0.61 -0.929 -0.503 -0.24 -0.594 0.017 0.023 -0.028 0.011 0.116 -0.028 0.074 3.351 cs1: (c) level ratio, 1 s reward -0.481 -0.477 -0.618 -0.604 -0.234 -0.551 0.025 0.032 -0.049 0.013 0.168 -0.041 0.102 4.144 cs1: day 1 reward, 1 s 0.478 0.472 0.57 0.647 0.392 0.646 2.478 1.296 1.97 1.25 1.193 0.885 1.271 -0.426 cs1: day 1 reward, 0.5 s 0.48 0.448 0.423 0.71 0.404 0.558 2.801 1.313 1.711 1.229 0.951 1.023 1.408 -0.471 all15: (a) AUC ratio, 1 s reward -0.328 -0.216 -0.304 -0.259 -0.107 -0.443 -0.029 0.007 -0.032 -0.007 0.106 -0.033 0.029 3.08 all15: (b) AUC ratio, 0.5 s reward -0.708 -0.547 -0.911 -0.436 -0.218 -0.972 -0.037 0.01 -0.036 -0.012 0.149 -0.045 0.041 4.436 all15: (c) level ratio, 1 s reward -0.655 -0.432 -0.607 -0.519 -0.213 -0.885 -0.057 0.014 -0.064 -0.015 0.212 -0.066 0.057 6.16 all15: day 1 reward, 1 s 0.428 0.487 0.574 0.683 0.399 0.532 2.286 1.273 1.942 1.216 1.259 0.865 1.212 -0.26 all15: day 1 reward, 0.5 s 0.426 0.466 0.427 0.743 0.411 0.452 2.656 1.296 1.698 1.201 0.987 1.007 1.36 -0.323 cs1: (a) AUC ratio, 1 s reward DAT-Cre mean -0.247 control 0.018 t = -8.62, p = 3.2e-06 cs1: (b) AUC ratio, 0.5 s reward DAT-Cre mean -0.565 control 0.026 t = -6.86, p = 2.72e-05 cs1: (c) level ratio, 1 s reward DAT-Cre mean -0.494 control 0.036 t = -8.62, p = 3.2e-06 cs1: day 1 reward, 1 s DAT-Cre mean 0.534 control 1.478 t = -4.13, p = 0.00168 cs1: day 1 reward, 0.5 s DAT-Cre mean 0.504 control 1.491 t = -3.76, p = 0.00316 all15: (a) AUC ratio, 1 s reward DAT-Cre mean -0.276 control 0.006 t = -5.99, p = 9.06e-05 all15: (b) AUC ratio, 0.5 s reward DAT-Cre mean -0.632 control 0.010 t = -5.74, p = 0.000129 all15: (c) level ratio, 1 s reward DAT-Cre mean -0.552 control 0.012 t = -5.99, p = 9.06e-05 all15: day 1 reward, 1 s DAT-Cre mean 0.517 control 1.436 t = -4.45, p = 0.000977 all15: day 1 reward, 0.5 s DAT-Cre mean 0.487 control 1.458 t = -3.99, p = 0.00211 Paper: inhibition t = -5.74, p = 1.3e-4; day 1 reward t = -3.99, p = 0.0021; DAT-Cre about -0.6
Result: reproduced. With the baseline before CS1, dopamine 1 to 1.5 s after CS2 is lower in the 6 DAT-Cre mice than in the 7 controls (t = -8.62, p = 3.2e-6 under readings a and c; t = -6.86, p = 2.7e-5 under b), and the day 1 reward response relative to the random reward session is lower (1 s window t = -4.13, p = 0.0017). The DAT-Cre mean is -0.565 under reading (b), within tolerance of -0.6, and -0.247 and -0.494 under (a) and (c). With the code's baseline over all code-15 events the tests give t = -5.74 (p = 1.3e-4) and, with the 0.5 s reward window, t = -3.99 (p = 0.0021), the paper's values, and the DAT-Cre mean is -0.632 under (b).
fig6i-learning-despite-inhibition¶
Despite inhibition of dopamine from CS2 to reward, every DAT-Cre animal learns anticipatory licking and acquires a positive dopamine response to CS1, both increasing over sessions as in controls (Test 11).
lc_rows, reg = [], []
for grp, mice in [("control", jeong.CONTROL_MICE), ("DAT-Cre", jeong.DAT_MICE)]:
X, Y, L = [], [], []
for m in mice:
c = e7[m]["cond"]
cs = [np.nanmean(x["cs1_15"]) for x in c]
rw = np.nanmean(c[0]["rw_15"])
lk = [np.nanmean(x["lick_cs1_rw"]) for x in c]
v = np.array(cs[:3] + cs[-2:]) / rw
l = np.array(lk[:3] + lk[-2:])
X += list(range(5)); Y += list(v); L += list(l)
lc_rows.append({"group": grp, "subject": m, "sessions": len(c), **{f"cs1_{i}": v[i] for i in range(5)},
**{f"lick_{i}": l[i] for i in range(5)},
"s12_first": np.nanmean(c[0]["lick_cs1_cs2"]), "s12_last": np.nanmean(c[-1]["lick_cs1_cs2"]),
"s12_code_first": np.nanmean(c[0]["lick_cs1_cs2_code"]), "s12_code_last": np.nanmean(c[-1]["lick_cs1_cs2_code"])})
for name, y in [("CS1 response", Y), ("anticipatory licks", L)]:
slope, t, p, n = ols_slope(X, y)
reg.append({"group": grp, "measure": name, "slope": slope, "t": t, "p": p, "n": n})
lc = pd.DataFrame(lc_rows)
print(lc.round(3).to_string(index=False))
print(pd.DataFrame(reg).round(6).to_string(index=False))
dat = lc[lc.group == "DAT-Cre"]
print(f"DAT-Cre animals with last-session CS1 response > 0: {(dat.cs1_4 > 0).sum()} of {len(dat)}; licks > 0: {(dat.lick_4 > 0).sum()} of {len(dat)}")
print("Paper: CS1 t = 5.49 (WT), 2.55 (DAT-Cre); licks t = 3.61 (WT), 4.92 (DAT-Cre)")
group subject sessions cs1_0 cs1_1 cs1_2 cs1_3 cs1_4 lick_0 lick_1 lick_2 lick_3 lick_4 s12_first s12_last s12_code_first s12_code_last control HJ-FP-datWT-stGtACR-F1 7 0.227 0.373 0.411 0.641 0.386 3.66 8.020 13.36 7.86 7.76 1.02 1.18 1.08 1.44 control HJ-FP-WT-stGtACR-F1 5 0.169 0.477 0.603 0.820 0.763 1.56 7.320 11.04 15.96 8.08 0.74 1.20 0.76 0.86 control HJ-FP-WT-stGtACR-F2 5 0.283 0.818 0.943 1.043 0.997 0.28 0.440 1.62 7.36 0.76 0.06 0.22 -0.36 -0.08 control HJ-FP-WT-stGtACR-F3 6 0.202 0.650 0.731 0.783 0.765 1.18 5.620 6.56 5.54 5.88 0.52 1.22 1.00 1.10 control HJ-FP-WT-stGtACR-M1 5 -0.010 0.334 0.897 0.999 1.383 0.54 1.580 7.24 5.64 9.06 0.72 3.06 0.38 3.26 control HJ-FP-WT-stGtACR-M2 5 0.134 0.572 0.621 0.462 0.459 -0.18 2.347 4.84 1.44 8.06 0.10 1.44 0.22 1.86 control HJ-FP-WT-stGtACR-M3 5 0.088 0.166 0.391 0.555 0.077 0.24 1.580 5.60 2.96 11.16 0.18 4.08 -0.02 4.28 DAT-Cre HJ-FP-datHT-stGtACR-M2 5 0.023 0.246 0.320 0.377 0.612 -0.04 1.960 1.52 5.26 10.78 0.06 4.38 -0.16 3.82 DAT-Cre HJ-FP-datHT-stGtACR-M3 6 0.009 -0.063 0.355 0.486 0.649 -0.16 0.643 1.30 3.04 1.32 -0.12 -0.06 0.10 -0.08 DAT-Cre HJ-FP-datHT-stGtACR-M4 6 0.044 0.137 0.206 0.311 0.275 0.24 0.960 1.60 7.64 9.14 -0.18 3.10 0.28 2.52 DAT-Cre HJ-FP-datHT-stGtACR-M5 8 -0.093 -0.117 -0.049 0.113 0.119 0.18 0.620 0.76 5.82 4.04 -0.14 0.36 0.00 0.68 DAT-Cre HJ-FP-datHT-stGtACR-M6 7 0.096 -0.016 0.048 0.368 0.221 0.86 1.380 3.00 8.18 3.24 0.34 0.58 1.14 0.50 DAT-Cre HJ-FP-datHT-stGtACR-M7 12 0.371 0.619 0.710 0.315 0.282 2.94 3.680 1.26 6.30 0.52 1.00 0.26 1.48 0.44 group measure slope t p n control CS1 response 0.134134 4.172377 0.000206 35 control anticipatory licks 1.525901 3.674701 0.000838 35 DAT-Cre CS1 response 0.076320 2.873745 0.007658 30 DAT-Cre anticipatory licks 1.283952 4.218591 0.000233 30 DAT-Cre animals with last-session CS1 response > 0: 6 of 6; licks > 0: 6 of 6 Paper: CS1 t = 5.49 (WT), 2.55 (DAT-Cre); licks t = 3.61 (WT), 4.92 (DAT-Cre)
Result: reproduced. Over the first three and last two sessions, the normalized CS1 response rises in controls (t = 4.17, p = 2.1e-4, 35 values) and in DAT-Cre mice (t = 2.87, p = 0.0077, 30 values), and so do anticipatory licks (t = 3.67, p = 8.4e-4; t = 4.22, p = 2.3e-4). In the last session every DAT-Cre animal has a positive CS1 response and positive anticipatory licking.
fig6i-lower-cs1-asymptote¶
With inhibition, the CS1 dopamine response reaches a lower level in the last session than in controls.
w, d = lc[lc.group == "control"].cs1_4, lc[lc.group == "DAT-Cre"].cs1_4
print(f"last-session normalized CS1 response: control mean {w.mean():.3f} (n {len(w)}), DAT-Cre {d.mean():.3f} (n {len(d)})")
show_t("two-sample t, control against DAT-Cre", stats.ttest_ind(w, d))
print("Paper: t = 2.80, p = 0.0171")
# Variants of the session and normalisation choices: the last session's CS1 response normalised
# by the random reward session's reward response (in place of the first conditioning session's),
# not normalised, with HJ-FP-WT-stGtACR-M3's last session (the shortest) replaced by the one
# before it, and with HJ-FP-WT-stGtACR-M2's last session (Day6, no dF/F for its first 1117 s)
# replaced by the one before it.
var = []
for grp, mice in [("control", jeong.CONTROL_MICE), ("DAT-Cre", jeong.DAT_MICE)]:
for m in mice:
c = e7[m]["cond"]
cs = [np.nanmean(x["cs1_15"]) for x in c]
last = -2 if m == "HJ-FP-WT-stGtACR-M3" else -1
last_m2 = -2 if m == "HJ-FP-WT-stGtACR-M2" else -1
var.append({"group": grp, "subject": m, "last_trials": len(c[-1]["cs1"]),
"last_nan": int(np.isnan(c[-1]["cs1_15"]).sum()),
"first_rw": np.nanmean(c[0]["rw_15"]), "random_rw": e7[m]["random"]["rw_15"],
"raw_last": cs[-1], "claim": cs[-1] / np.nanmean(c[0]["rw_15"]),
"by_random_rw": cs[-1] / e7[m]["random"]["rw_15"],
"m3_previous": cs[last] / np.nanmean(c[0]["rw_15"]),
"m2_previous": cs[last_m2] / np.nanmean(c[0]["rw_15"])})
var = pd.DataFrame(var)
print(var.round(3).to_string(index=False))
for col in ["claim", "by_random_rw", "raw_last", "m3_previous", "m2_previous"]:
w, d = var[var.group == "control"][col], var[var.group == "DAT-Cre"][col]
print(f"{col}: control mean {w.mean():.3f}, DAT-Cre {d.mean():.3f}"); show_t(" two-sample t", stats.ttest_ind(w, d))
last-session normalized CS1 response: control mean 0.690 (n 7), DAT-Cre 0.359 (n 6) two-sample t, control against DAT-Cre: t = 1.700, p = 0.1172, df = 11 Paper: t = 2.80, p = 0.0171 group subject last_trials last_nan first_rw random_rw raw_last claim by_random_rw m3_previous m2_previous control HJ-FP-datWT-stGtACR-F1 50 0 3.640 1.534 1.405 0.386 0.916 0.386 0.386 control HJ-FP-WT-stGtACR-F1 50 0 4.285 3.189 3.267 0.763 1.025 0.763 0.763 control HJ-FP-WT-stGtACR-F2 50 0 8.126 3.456 8.105 0.997 2.345 0.997 0.997 control HJ-FP-WT-stGtACR-F3 50 0 13.362 10.407 10.220 0.765 0.982 0.765 0.765 control HJ-FP-WT-stGtACR-M1 50 0 3.409 3.081 4.715 1.383 1.530 1.383 1.383 control HJ-FP-WT-stGtACR-M2 50 17 8.449 8.697 3.881 0.459 0.446 0.459 0.462 control HJ-FP-WT-stGtACR-M3 25 0 3.242 2.887 0.251 0.077 0.087 0.555 0.077 DAT-Cre HJ-FP-datHT-stGtACR-M2 50 0 0.486 1.096 0.297 0.612 0.271 0.612 0.612 DAT-Cre HJ-FP-datHT-stGtACR-M3 50 0 1.049 2.144 0.680 0.649 0.317 0.649 0.649 DAT-Cre HJ-FP-datHT-stGtACR-M4 50 0 0.655 1.097 0.180 0.275 0.164 0.275 0.275 DAT-Cre HJ-FP-datHT-stGtACR-M5 50 0 0.412 0.607 0.049 0.119 0.081 0.119 0.119 DAT-Cre HJ-FP-datHT-stGtACR-M6 50 0 0.201 0.499 0.044 0.221 0.089 0.221 0.221 DAT-Cre HJ-FP-datHT-stGtACR-M7 50 0 3.703 5.208 1.043 0.282 0.200 0.282 0.282 claim: control mean 0.690, DAT-Cre 0.359 two-sample t: t = 1.700, p = 0.1172, df = 11 by_random_rw: control mean 1.047, DAT-Cre 0.187 two-sample t: t = 2.836, p = 0.01619, df = 11 raw_last: control mean 4.549, DAT-Cre 0.382 two-sample t: t = 2.848, p = 0.01586, df = 11 m3_previous: control mean 0.758, DAT-Cre 0.359 two-sample t: t = 2.434, p = 0.03317, df = 11 m2_previous: control mean 0.690, DAT-Cre 0.359 two-sample t: t = 1.703, p = 0.1167, df = 11
Result: not reproduced. The last-session CS1 response is lower in DAT-Cre mice on average (0.359 against 0.690 of the first conditioning session's reward response) but not significantly, t(11) = 1.70, p = 0.117 (paper t = 2.80, p = 0.0171). One control (WT-M3, 0.077) is below every DAT-Cre value; its last session has 25 trials, where every other animal's has 50. The outcome turns on these choices: with WT-M3's previous session in place of its last, t = 2.43 (p = 0.033); normalised by the random reward session's reward response, t = 2.84 (p = 0.016); not normalised, t = 2.85 (p = 0.016). WT-M2's last session has no dF/F for its first 1117 s, so 17 of its 50 CS1 trials have no value and are left out of its mean; with its previous session in place of it, its value is 0.462 against 0.459 and t = 1.70 (p = 0.117), so these trials do not bear on the outcome. The normalisation of the criterion is the one Supplementary Note 10 and learnign_curve_sequential.py name ("normalized by reward response on day 1 of conditioning"), and no passage excludes a short session, so the outcome is kept on the criterion's measure.
fig6i-last-session-licking-similar¶
In the last session, anticipatory licking does not differ between DAT-Cre and control animals.
w, d = lc[lc.group == "control"].lick_4, lc[lc.group == "DAT-Cre"].lick_4
print(f"last-session licks (CS1 to reward minus baseline): control mean {w.mean():.2f}, DAT-Cre {d.mean():.2f}, difference {w.mean() - d.mean():.2f}")
show_t("two-sample t, control against DAT-Cre", stats.ttest_ind(w, d))
print("Paper: t = 0.88, p = 0.3999")
last-session licks (CS1 to reward minus baseline): control mean 7.25, DAT-Cre 4.84, difference 2.41 two-sample t, control against DAT-Cre: t = 1.165, p = 0.2686, df = 11 Paper: t = 0.88, p = 0.3999
Result: reproduced. Last-session anticipatory licking is 7.25 licks in controls and 4.84 in DAT-Cre mice, t(11) = 1.17, p = 0.27 (paper t = 0.88, p = 0.40): no difference is found.
figs12f-pre-cs2-licking¶
Anticipatory licking between CS1 and CS2, before inhibition starts, increases from the first to the last session (main effect of session), with no effect of genotype and no genotype by session interaction.
for label, a, b in [("claim (baseline 1.5 s before CS1)", "s12_first", "s12_last"),
("code variant (baseline 1500 s before CS1)", "s12_code_first", "s12_code_last")]:
long = pd.concat([lc[["subject", "group", a]].rename(columns={a: "licks"}).assign(session="first"),
lc[["subject", "group", b]].rename(columns={b: "licks"}).assign(session="last")])
print(f"\n{label}: mean first {lc[a].mean():.2f}, last {lc[b].mean():.2f}")
print(pg.mixed_anova(long, dv="licks", within="session", between="group", subject="subject").round(4).to_string(index=False))
print("\nPaper: session F = 9.07, p = 0.0118; genotype F = 0.25, p = 0.6241; interaction F = 0.07, p = 0.7915")
claim (baseline 1.5 s before CS1): mean first 0.33, last 1.62
Source SS DF1 DF2 MS F p_unc np2 eps
group 0.6865 1 11 0.6865 0.5681 0.4668 0.0491 NaN
session 10.7522 1 11 10.7522 7.3903 0.0200 0.4019 1.0
Interaction 0.0005 1 11 0.0005 0.0003 0.9855 0.0000 1.0
code variant (baseline 1500 s before CS1): mean first 0.45, last 1.58
Source SS DF1 DF2 MS F p_unc np2 eps
group 0.3532 1 11 0.3532 0.3431 0.5699 0.0302 NaN
session 8.3112 1 11 8.3112 5.2269 0.0431 0.3221 1.0
Interaction 0.4710 1 11 0.4710 0.2962 0.5971 0.0262 1.0
Paper: session F = 9.07, p = 0.0118; genotype F = 0.25, p = 0.6241; interaction F = 0.07, p = 0.7915
Result: reproduced. Pre-CS2 licking rises from the first to the last session (from 0.33 to 1.62 licks; session F(1, 11) = 7.39, p = 0.020), with no genotype effect (F = 0.57, p = 0.47) and no interaction (F = 0.0003, p = 0.99). With the code's baseline window the session effect is F = 5.23, p = 0.043, and the others stay non-significant.
Simulation claims (blocked)¶
fig2-anccr-classic-results, fig3cf-simulated-unpredicted-rewards, fig5b-simulated-intermediate-cue, fig6h-simulated-inhibition, figs6-learns-causal-structure, figs5-timescale-invariance and figs13-ramping are simulations of the ANCCR and TDRL models. They use no recorded data, and the paper's simulation code is Matlab (simulation/ and functions/calculateANCCR.m of the ANCCR repository). Matlab cannot be installed in the entry's locked Python environment, and GNU Octave, which might run the scripts, is not on the runner and cannot be installed from the lock file either. A port of the models to Python was not attempted. They are blocked as not_runnable, the protocol's reason for a method whose software cannot be run here; nothing is computed for them.
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Data | Outcome |
|---|---|---|
| session-day-link | with labels | partially reproduced |
| methods-session-counts | with labels | reproduced |
| exp1-excluded-rewards | dataset | reproduced |
| fig3e-reward-response-increases | dataset | reproduced |
| fig3g-response-iri-positive | dataset | reproduced |
| figs8k-iri-correlation-stable | dataset | reproduced |
| figs8b-lick-rates-constant | dataset | reproduced |
| figs8d-consummatory-lick-uncorrelated | dataset | reproduced |
| figs8e-lick-bout-responses | dataset | blocked: methods underspecified |
| figs8g-baseline-unchanged | dataset | reproduced |
| figs8h-baseline-uncorrelated-with-reward-change | dataset | reproduced |
| figs8l-lick-rate-tracks-iri | dataset | reproduced |
| fig4ef-dopamine-leads-behavior | with labels | reproduced |
| figs10c-cue-response-at-acquisition | with labels | reproduced |
| fig4hi-delay-change-behavior-changes | with labels | reproduced |
| fig4hi-delay-change-dopamine-unchanged | with labels | reproduced |
| figs10b-no-omission-response | with labels | reproduced |
| fig4kl-extinction-abrupt-behavior | with labels | reproduced |
| fig4jl-extinction-cue-response-positive | with labels | reproduced |
| fig4p-background-faster-decay | with labels | not reproduced |
| fig5d-intermediate-cue-positive | with labels | reproduced |
| fig5d-intermediate-reward-larger | with labels | reproduced |
| figs11de-intermediate-cue-local-baseline | with labels | reproduced |
| fig6b-no-backpropagation | with labels | reproduced |
| fig6d-cues-increase-together | dataset | reproduced |
| fig6d-cs2-below-cs1-late | dataset | reproduced |
| fig6g-inhibition-magnitude | dataset | reproduced |
| fig6i-learning-despite-inhibition | dataset | reproduced |
| fig6i-lower-cs1-asymptote | dataset | not reproduced |
| fig6i-last-session-licking-similar | dataset | reproduced |
| figs12f-pre-cs2-licking | dataset | reproduced |
| seven simulation claims | dataset | blocked: not runnable |
How the analysis differs from the paper's. The measures follow the paper's code at the commit the claims file cites, ported to Python; where the Methods and the code differ, both are shown. The Matlab curve fits start at random points; here the best of 20 seeded starts is kept. dF/F is divided by 100 for the Matlab-derived measures, as the files store percent and the Methods define dF/F as a fraction. The code's bounds on the dopamine model bind on either scale, so the dopamine model comparison (fig4hi-delay-change-dopamine-unchanged) cannot detect a change here; without the bounds it finds one in every animal. That outcome rests on the assumption that the paper's fit was bounded as here. As in the code, trials with NaN dF/F samples in their windows keep the nansum of the others (0 outside the coverage); they are 37 of 8721 CS+ trials, and leaving them out changes no outcome. In Experiment 7, where a trial with no dF/F sample in its window has no value and drops out of the session mean, 17 of the 50 CS1 trials of HJ-FP-WT-stGtACR-M2's last session are such trials; replacing that session changes no outcome. Seven integer days inside the code's day ranges have no session in the dandiset, and the analysis runs on the sessions present.
What was found about the dandiset.
- Event times are in seconds, though the files state milliseconds.
- Seven integer days within the code's day ranges have no session in the dandiset (M2 Day23, M3 Day16 and Day31, M4 Day30 and Day38, F2 Day16, M6 Day49). M4 Day38 is a listed start day, its first background reward session, so that session is missing; the paper's session counts of Experiments 2, 4 and 5 equal the counts with all seven days, which suggests the others are missing sessions too.
- In the large Experiment 7 files the raw photometry timestamps are in thousands of seconds, though stated in seconds.
- The DAT-Cre mouse HJ-FP-datHT-stGtACR-M7 has genotype
wtin its file. - Some sessions lose the start or the end of their photometry relative to the event log.
- No recording is deposited twice; dF/F and the event log are consistent with one clock, except in 27 sessions with no clear response to the reference event, where the check cannot show it.
print(f"Notebook run time: {(time.time() - started) / 60:.1f} min")
Notebook run time: 1.4 min