Reproducing Vu et al. (2024) from dandiset 001084¶
Paper. Vu MT, Brown EH, Wen MJ, Noggle CA, Zhang Z, Monk KJ, Bouabid S, Mroz L, Graham BM, Zhuo Y, Li Y, Otchy TM, Tian L, Davison IG, Boas DA, Howe MW (2024). Targeted micro-fiber arrays for measuring and manipulating localized multi-scale neural dynamics over large, deep brain volumes during behavior. Neuron 112, 909-923. doi:10.1016/j.neuron.2023.12.011. This attempt used the full text of the NIH author manuscript in PubMed Central (PMC10957316). The supplementary PDF could not be retrieved, so the supplementary figure legends were not read.
Dandiset. DANDI:001084, version 0.241023.2011. The paper cites no DANDI version (its data statement points to Zenodo). This is the only published version, so it is pinned; the draft has the same 180 assets and size.
Code read. The paper's analysis code on Zenodo (doi:10.5281/zenodo.10272874), which settles several choices the Methods leave open (the 11-frame tone delay, the spline basis, the latency rule, how stimulation times are found), and the conversion code catalystneuro/howe-lab-to-nwb at commit 3c53055.
Outside inputs. One claim (fig4-link-published-results) reads the paper's per-fiber results
(glm_spline_results.mat, Zenodo Fig04.zip) to check that the dandiset holds the data behind them.
Every other claim uses the dandiset alone.
Random numbers: each section that draws them has its own generator with a fixed seed.
| Claim | Data | Status |
|---|---|---|
| fig1c-isosbestic-control | dandiset | tested |
| figs2d-correlation-distance | dandiset | tested |
| figs2ef-no-bundle-crosstalk | dandiset | tested |
| fig1d-chronic-reward-stability | dandiset | tested |
| fig4-link-published-results | with inputs | tested |
| fig4e-light-increases-widespread | dandiset | tested |
| fig4e-light-anterior-posterior-gradient | dandiset | tested |
| fig4f-tone-tail-increase-anterior-decrease | dandiset | tested |
| fig4f-decreases-lateral-ventral | dandiset | tested |
| figs4ac-patterns-consistent-across-mice | dandiset | tested |
| figs4eh-patterns-not-movement | dandiset | tested |
| figs5-isosbestic-no-pattern | dandiset | tested |
| fig5ab-response-latencies | dandiset | tested |
| fig5g-light-latency-gradient | dandiset | tested |
| fig5h-tone-latency-gradient | dandiset | tested |
| figs6cd-gradients-individual-mice | dandiset | tested |
| fig6-inclusion-counts | dandiset | tested |
| fig6gh-release-decays-with-distance | dandiset | tested |
| fig6g-dorsal-spread | dandiset | tested |
| fig6-tdtomato-control-no-response | dandiset | tested |
| figs7a-photoswitch-confined | dandiset | tested |
| figs1-no-tissue-damage | Zenodo | blocked: not attempted (data only on Zenodo) |
| figs3-targeting-errors | dandiset | blocked: data_missing (planned targets not public) |
| fig3-freely-moving-signals | Zenodo | blocked: not attempted (data only on Zenodo) |
| fig7d-ventral-not-dorsal-jaw | Zenodo | blocked: not attempted (data only on Zenodo) |
| fig7ef-single-ventral-fiber | Zenodo | blocked: not attempted (data only on Zenodo) |
import os
for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
os.environ.setdefault(_v, "1") # many small least-squares fits
import io
import time
import warnings
from concurrent.futures import ThreadPoolExecutor
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
import scipy.io as sio
import statsmodels.api as sm
from scipy import stats
import vu2024 as V
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
started = time.time()
print("dandiset", V.DANDISET_ID, "version", V.DANDISET_VERSION)
dandiset 001084 version 0.241023.2011
Data check¶
Each check of the description against the values is stated under its output.
r = requests.get(f"https://api.dandiarchive.org/api/dandisets/{V.DANDISET_ID}/versions/{V.DANDISET_VERSION}/")
r.raise_for_status()
meta = r.json()
summary = meta["assetsSummary"]
print("Archive summary: subjects", summary["numberOfSubjects"], "; files", summary["numberOfFiles"],
"; species", [s["name"] for s in summary["species"]])
print("approach:", [a["name"] for a in summary["approach"]])
print("measurementTechnique:", [a["name"] for a in summary["measurementTechnique"]])
print("variableMeasured:", summary["variableMeasured"])
paths = V.list_nwb_paths()
print("\nNWB files:", len(paths))
def survey(path):
f = V.open_nwb(path)
sub = f["general/subject"]
row = {"path": path, "subject": V._text(sub["subject_id"][()]), "genotype": V._text(sub["genotype"][()]),
"sex": V._text(sub["sex"][()]), "session": V._text(f["general/session_id"][()]),
"start": V._text(f["session_start_time"][()])[:19]}
series = [k for k in f["acquisition"] if k.startswith("FiberPhotometryResponseSeries")]
row["channels"] = ",".join(s.replace("FiberPhotometryResponseSeries", "") for s in series)
g = f["acquisition"][series[0]]
row["frames"] = g["data"].shape[0]
row["rate_hz"] = round(float(g["starting_time"].attrs["rate"]), 2)
row["fibers"] = g["data"].shape[1]
tab = V.fiber_table(f, series[0])
row["included"] = int(tab.included.sum())
row["fiber table columns"] = ",".join(sorted(f["general/FiberPhotometry/FiberPhotometryTable"].keys()))
ip = f["general/optophysiology"]
row["indicator"] = ",".join(sorted({V._text(ip[k]["indicator"][()]) for k in ip}))
ev = V.read_events(f)
for kind in ("Light", "Tone", "Reward", "Lick"):
row[kind] = int((ev.event_type == kind).sum())
row["videos"] = sum("Video" in k for k in f["acquisition"])
return row
with ThreadPoolExecutor(16) as ex:
files = pd.DataFrame(list(ex.map(survey, paths)))
files = files.sort_values(["subject", "session"]).reset_index(drop=True)
print("Distinct subjects in the files:", files.subject.nunique())
per_subject = files.groupby("subject").agg(
sessions=("session", "size"), genotype=("genotype", "first"), sex=("sex", "first"),
indicator=("indicator", "first"), channels=("channels", "first"), rate_hz=("rate_hz", "median"),
fibers=("fibers", "first"), included=("included", "first"), light=("Light", "sum"), tone=("Tone", "sum"),
reward=("Reward", "sum"), first=("start", "min"), last=("start", "max"))
print(per_subject.to_string())
Archive summary: subjects 16 ; files 180 ; species ['Mus musculus - House mouse'] approach: ['microscopy approach; cell population imaging'] measurementTechnique: ['analytical technique', 'surgical technique'] variableMeasured: ['OpticalChannel', 'ImagingPlane', 'PlaneSegmentation', 'ProcessingModule']
NWB files: 116
Distinct subjects in the files: 16
sessions genotype sex indicator channels rate_hz fibers included light tone reward first last
subject
842 2 D1-Cre M GCaMP7f Green,GreenIsosbestic 18.000 82 61 0 0 95 2022-11-18T16:10:03 2022-11-23T17:44:45
DL15 46 wildtype F dLight1.3b Green 29.995 64 44 0 0 1923 2021-08-19T11:44:21 2021-10-15T07:56:38
DL18 1 wildtype F dLight1.3b Green 29.990 103 19 14 44 0 2021-11-10T16:12:41 2021-11-10T16:12:41
DL20 1 wildtype F dLight1.3b Green 29.990 84 53 20 20 0 2022-01-11T11:41:45 2022-01-11T11:41:45
DL21 1 wildtype F dLight1.3b Green 29.990 74 51 18 18 0 2022-01-12T11:11:14 2022-01-12T11:11:14
DL23 1 wildtype M dLight1.3b Green 30.000 74 27 18 18 0 2022-01-17T10:42:29 2022-01-17T10:42:29
DL27 1 wildtype F dLight1.3b Green 29.990 73 47 21 21 0 2022-09-30T08:18:51 2022-09-30T08:18:51
DL29 1 wildtype M dLight1.3b Green 29.990 74 45 21 21 0 2022-09-30T07:43:57 2022-09-30T07:43:57
DL30 1 wildtype M dLight1.3b Green 29.990 74 38 21 21 0 2022-09-30T07:06:00 2022-09-30T07:06:00
DL31 1 wildtype M dLight1.3b Green,GreenIsosbestic 17.980 73 52 10 10 0 2023-10-20T10:33:52 2023-10-20T10:33:52
DL32 1 wildtype M dLight1.3b Green,GreenIsosbestic 17.980 73 44 10 10 0 2023-10-20T10:08:17 2023-10-20T10:08:17
Grid-9 1 DAT-Cre M Ach3.0,jRGECO1a Green,Red 11.000 97 85 0 44 45 2021-08-21T10:48:10 2021-08-21T10:48:10
Opto2 30 Dat-Cre F rDA3m Red 29.990 41 22 0 0 0 2023-03-14T16:45:16 2023-03-22T14:10:36
Opto4 13 Dat-Cre F tdTomato Red 29.990 40 24 0 0 0 2023-09-29T20:55:57 2023-10-06T20:35:41
Opto9 12 Dat-Cre F tdTomato Red 29.990 34 17 0 0 0 2023-09-29T21:46:22 2023-10-06T20:56:06
UG21 3 wildtype F tdTomato Red 30.000 30 22 0 0 0 2023-03-21T13:21:20 2023-03-21T13:42:51
The 16 subjects in the files match the archive's count. The archive lists the approach as "cell population imaging"; the files hold fiber photometry (one trace per optical fiber) and the movies of the fiber bundle it was extracted from. The subjects map onto the paper's experiments as follows, by genotype, indicator, date and the paper's code:
- DL18, DL20, DL21, DL23, DL27, DL29, DL30: C57BL/6 with dLight1.3b, one light and tone session each (Figures 4, 5, S4, S6).
- DL15: dLight1.3b, 46 reward sessions over 8 weeks (Figure 1D).
- 842: D1-cre with GCaMP7f, 470 and 405 nm interleaved at 18 Hz (Figure 1C, S2D-F).
- DL31, DL32: dLight1.3b, 470 and 415 nm interleaved at 18 Hz, light and tone (Figure S5E-H).
- Grid-9: DAT-cre with Ach3.0 and jRGECO1a at 11 Hz (Figure S2C); no claim uses it.
- Opto2: DAT-cre with rDA3m (rDAm3.0), single-fiber stimulation sessions (Figure 6).
- Opto4, Opto9: DAT-cre with tdTomato, 13 stimulation sessions each, with face videos: the tdTomato controls of Figure 6 ("26 single-site stimulations across 2 mice").
- UG21: wild type, 3 stimulation sessions, indicator given as tdTomato. Its 3 sessions and 22 included fibers match the paper's rDAm3.0 mouse without opsin ("0/22 non-stimulated sites for 3 single-site stimulations"), and the paper's code names these sessions as the photoswitching controls of Figure S7A. The indicator in the files is consistent with being a mislabel; the files cannot show which indicator was expressed.
The D1-cre mice of Figure 7 (OptoD1_6, 7, 8, 10 in the paper's code), the freely moving mice of Figure 3, the histology of Figure S1, the CT scans of Figure 2, and the 470/405 nm sessions of Figure S5A-D are not in the dandiset.
# Columns of the fiber table in every file, with their descriptions in one file: does it hold
# planned target positions as well as the localized tips (claim figs3-targeting-errors)?
print("Distinct sets of fiber table columns over the", len(files), "files:")
for cols, n in files["fiber table columns"].value_counts().items():
print(f" {n} files: {cols}")
f = V.open_nwb(V.SALIENT_SESSIONS["DL18"])
t = f["general/FiberPhotometry/FiberPhotometryTable"]
print("\nFiberPhotometryTable of", V.SALIENT_SESSIONS["DL18"], "-", V._text(t.attrs.get("description", "")))
for k in t.keys():
print(f" {k}: {V._text(t[k].attrs.get('description', ''))}")
Distinct sets of fiber table columns over the 116 files: 116 files: allen_atlas_coordinates,coordinates,dichroic_mirror,emission_filter,excitation_filter,excitation_source,id,included,indicator,location,optical_fiber,photodetector
FiberPhotometryTable of sub-DL18/sub-DL18_ses-211110_image+ophys.nwb - Contains the metadata for the fiber photometry experiment. allen_atlas_coordinates: The fiber tip coordinates (AP, ML, DV) in Allen Brain Atlas coordinates coordinates: Fiber placement in stereotactic coordinates (AP, ML, DV) mm relative to Bregma. dichroic_mirror: Link to the dichroic mirror device. emission_filter: Link to the emission filter device. excitation_filter: Link to the excitation filter device. excitation_source: Link to the excitation source device. id: included: Whether this fiber has been successfully implanted. indicator: Link to the indicator object. location: Location of fiber. optical_fiber: Link to the optical fiber device. photodetector: Link to the photodetector device.
The fiber table gives each fiber's localized tip (coordinates in mm from bregma and in the Allen
atlas), its location label and the included flag; no column holds a planned target, so the
targeting errors of Figure S3 cannot be computed from the dandiset.
def load_salient(mouse):
f = V.open_nwb(V.SALIENT_SESSIONS[mouse])
dff, times, tab = V.read_dff(f)
beh = V.read_behavior(f)
ev = V.read_events(f)
fs = 1 / np.median(np.diff(beh["t"]))
ev["on"] = V.nearest_index(beh["t"], ev.start_time.values)
ev["off"] = V.nearest_index(beh["t"], ev.stop_time.values)
return {"mouse": mouse, "sex": V._text(f["general/subject/sex"][()]), "dff_all": dff, "times": times,
"table": tab, "beh": beh, "events": ev, "fs": fs}
with ThreadPoolExecutor(7) as ex:
S = {d["mouse"]: d for d in ex.map(load_salient, V.SALIENT_SESSIONS)}
rows = []
for m, d in S.items():
ev, tab, beh = d["events"], d["table"], d["beh"]
light, tone = ev[ev.event_type == "Light"], ev[ev.event_type == "Tone"]
d["dff"] = d["dff_all"][:, tab.included.values]
d["fib"] = tab[tab.included].reset_index(drop=True)
rows.append({
"mouse": m, "sex": d["sex"], "frames": len(d["times"]), "behaviour samples": len(beh["t"]),
"fs_hz": round(d["fs"], 3), "max |fiber t - behaviour t| ms": round(1000 * np.abs(d["times"] - beh["t"]).max(), 1),
"fibers": len(tab), "included": int(tab.included.sum()), "light": len(light), "tone": len(tone),
"light s": f"{(light.stop_time - light.start_time).min():.2f}-{(light.stop_time - light.start_time).max():.2f}",
"tone s": f"{(tone.stop_time - tone.start_time).min():.2f}-{(tone.stop_time - tone.start_time).max():.2f}",
"ITI s": f"{np.diff(np.sort(ev[ev.event_type.isin(['Light', 'Tone'])].start_time)).min():.1f}-"
f"{np.diff(np.sort(ev[ev.event_type.isin(['Light', 'Tone'])].start_time)).max():.1f}",
"dff min": round(float(d["dff"].min()), 3), "dff max": round(float(d["dff"].max()), 3),
"roll m/s": f"{beh['roll'].min():.2f}..{beh['roll'].max():.2f}",
"pitch m/s": f"{beh['pitch'].min():.2f}..{beh['pitch'].max():.2f}",
"angular rad/s": f"{beh['angular'].min():.1f}..{beh['angular'].max():.1f}"})
sal = pd.DataFrame(rows)
print(sal.to_string(index=False))
print("\nTotal included fibers:", sal.included.sum(), "; sexes:", sal.sex.value_counts().to_dict())
print("Paper: 280 fibers from 7 mice (4F + 3M); 14 to 21 presentations of each modality; ITI 4 to 40 s.")
locs = pd.concat([d["fib"] for d in S.values()]).location.value_counts()
print("\nLocations of the included fibers:", locs.to_dict())
excl = pd.concat([d["table"][~d["table"].included] for d in S.values()]).location.value_counts()
print("Locations of the other fibers (top 8):", excl.head(8).to_dict())
allfib = pd.concat([d["fib"].assign(mouse=m) for m, d in S.items()])
print("\nIncluded fiber coordinates (mm from bregma): AP %.2f to %.2f, ML %.2f to %.2f, DV %.2f to %.2f"
% (allfib.AP.min(), allfib.AP.max(), allfib.ML.min(), allfib.ML.max(), allfib.DV.min(), allfib.DV.max()))
print("Fibers in the tail (AP <= -1.0 mm):", int((allfib.AP <= -1.0).sum()),
"; per mouse:", allfib[allfib.AP <= -1.0].groupby("mouse").size().to_dict())
mouse sex frames behaviour samples fs_hz max |fiber t - behaviour t| ms fibers included light tone light s tone s ITI s dff min dff max roll m/s pitch m/s angular rad/s
DL18 F 40000 40000 30.075 24.7 103 19 14 44 0.80-0.90 1.30-1.37 5.5-41.3 -0.059 0.195 -0.65..0.57 -0.38..0.73 -8.6..13.8
DL20 F 28283 28283 30.075 18.0 84 53 20 20 0.73-0.93 0.93-1.37 5.1-41.4 -0.060 0.199 -0.70..0.25 -0.24..1.23 -7.6..17.7
DL21 F 23484 23484 30.075 14.9 74 51 18 18 0.73-0.97 0.90-1.37 4.9-39.1 -0.080 0.259 -0.18..0.16 -0.26..0.19 -7.9..6.9
DL23 M 27400 27400 30.075 10.6 74 27 18 18 0.80-0.93 1.17-1.37 5.4-39.9 -0.039 0.177 -0.30..0.18 -0.28..0.20 -6.0..7.4
DL27 F 25398 25398 30.075 16.0 73 47 21 21 0.80-0.93 1.30-1.43 4.9-40.0 -0.039 0.166 -0.51..0.43 -0.20..0.72 -10.1..5.6
DL29 M 30454 30454 30.075 19.2 74 45 21 21 0.77-0.90 1.20-1.37 4.8-41.6 -0.074 0.218 -0.59..0.53 -0.30..0.81 -12.9..7.3
DL30 M 25223 25223 30.075 16.3 74 38 21 21 0.80-0.93 0.57-1.37 5.4-40.4 -0.040 0.156 -0.33..0.36 -0.36..0.43 -11.7..6.7
Total included fibers: 280 ; sexes: {'F': 4, 'M': 3}
Paper: 280 fibers from 7 mice (4F + 3M); 14 to 21 presentations of each modality; ITI 4 to 40 s.
Locations of the included fibers: {'Caudoputamen': 255, 'Nucleus accumbens': 24, 'Striatum': 1}
Locations of the other fibers (top 8): {'Caudoputamen': 87, '': 16, 'lateral ventricle': 15, 'Primary motor area Layer 6a': 12, 'Supplemental somatosensory area layer 6a': 11, 'Primary somatosensory area mouth layer 6a': 9, 'Globus pallidus external segment': 9, 'Striatum': 7}
Included fiber coordinates (mm from bregma): AP -2.12 to 1.76, ML 0.63 to 3.56, DV 1.78 to 5.49
Fibers in the tail (AP <= -1.0 mm): 35 ; per mouse: {'DL18': 2, 'DL20': 9, 'DL21': 8, 'DL23': 1, 'DL27': 5, 'DL29': 3, 'DL30': 7}
The included fibers number 280, from 4 female and 3 male mice, as in the paper; they are the fibers
the paper analysed (its code keeps fibers "in_striatum", and the per-mouse counts equal those of the
paper's results file, see fig4-link-published-results). The column describes them as "successfully
implanted"; the paper's criterion is a tip within or on the border of the striatum, and the included
fibers are labelled Caudoputamen, Nucleus accumbens or Striatum, while 87 Caudoputamen fibers are not
included. Fiber traces and behaviour have the same number of samples; the fiber series is stored with a
regular rate and the behaviour with timestamps, which differ by at most 25 ms, and every event starts
on a behaviour timestamp, so events are placed on frames by their behaviour sample, as the code does.
Velocity is in m/s and angular velocity in rad/s, with ranges plausible for a mouse on a ball. Six
mice have 18 to 21 presentations of each stimulus; DL18 has 14 light and 44 tone presentations. The
paper's code uses only 12 kHz tones and its results file has 14 for DL18; the files do not say which
of the 44 tones are 12 kHz. Light presentations last 0.7 to 1.0 s and tones about 1.33 s on the
trigger line; the code places the tone 11 frames after its trigger (speaker delay).
# Clock check: decode light against tone from the population response after each onset, with the
# true onsets and with all onsets shifted by +-5 s (labels shifted in time).
TONE_SHIFT = 11 # frames at 30 Hz (paper's code)
def decode_accuracy_fs(dff, onsets, labels, shift, w):
on = onsets + shift
ok = (on - w >= 0) & (on + w + 1 <= len(dff))
on, lab = on[ok], labels[ok]
feats = np.stack([dff[o:o + w].mean(0) - dff[o - w:o].mean(0) for o in on])
correct = 0
for i in range(len(on)): # leave-one-out nearest centroid on z-scored features
tr = np.ones(len(on), bool)
tr[i] = False
mu, sd = feats[tr].mean(0), feats[tr].std(0) + 1e-9
z = (feats - mu) / sd
c0, c1 = z[tr & (lab == 0)].mean(0), z[tr & (lab == 1)].mean(0)
pred = int(np.sum((z[i] - c1) ** 2) < np.sum((z[i] - c0) ** 2))
correct += pred == lab[i]
return correct / len(on), len(on)
def decode_accuracy(dff, onsets, labels, shift):
return decode_accuracy_fs(dff, onsets, labels, shift, 30)
rows = []
for m, d in S.items():
ev = d["events"]
st = ev[ev.event_type.isin(["Light", "Tone"])]
lab = (st.event_type == "Tone").values.astype(int)
on = st.on.values + TONE_SHIFT * lab
row = {"mouse": m}
for shift in (0, 150, -150):
acc, n = decode_accuracy(d["dff"], on, lab, shift)
row[f"shift {shift / 30:+.0f} s"] = round(acc, 2)
row["presentations"] = n
row["chance (majority)"] = round(max(lab.mean(), 1 - lab.mean()), 2)
rows.append(row)
print(pd.DataFrame(rows).to_string(index=False))
mouse shift +0 s shift +5 s shift -5 s presentations chance (majority) DL18 0.90 0.51 0.48 58 0.76 DL20 0.92 0.52 0.60 40 0.50 DL21 0.94 0.61 0.53 36 0.50 DL23 0.92 0.71 0.36 36 0.50 DL27 0.86 0.52 0.54 41 0.50 DL29 0.85 0.46 0.56 41 0.50 DL30 0.76 0.50 0.56 41 0.50
With the true onsets the stimulus type is decoded from the fibers' responses well above chance in every mouse, and with onsets shifted by 5 s it falls to about chance (for DL18, whose tones outnumber lights 44 to 14, chance is 0.76). The event table and the fiber traces share one clock, as far as this contrast can show.
def load_two_channel(path):
f = V.open_nwb(path)
out = {"path": path, "beh": V.read_behavior(f), "events": V.read_events(f)}
for color in ("Green", "GreenIsosbestic"):
dff, times, tab = V.read_dff(f, color)
out[color] = {"dff": dff[:, tab.included.values], "times": times, "fib": tab[tab.included].reset_index(drop=True),
"rate": 1 / np.median(np.diff(times))}
seg = f["processing/ophys/ImageSegmentation/PlaneSegmentation"]
out["centroids"] = seg["ROICentroids"][:]
out["included_mask"] = V.fiber_table(f, "FiberPhotometryResponseSeriesGreen").included.values
return out
with ThreadPoolExecutor(4) as ex:
TWO = dict(zip(["842-221123", "842-221118", "DL31", "DL32"], ex.map(load_two_channel, [
"sub-842/sub-842_ses-221123_ophys.nwb", "sub-842/sub-842_ses-221118_ophys.nwb",
V.ISOSBESTIC_SESSIONS["DL31"], V.ISOSBESTIC_SESSIONS["DL32"]])))
for k, d in TWO.items():
g, i = d["Green"], d["GreenIsosbestic"]
ev = d["events"]
print(f"{k}: 470 nm {g['dff'].shape} at {g['rate']:.2f} Hz from {g['times'][0]:.2f} s; isosbestic {i['dff'].shape}; "
f"behaviour samples {len(d['beh']['t'])} from {d['beh']['t'][0]:.2f} s; ROI centroids {d['centroids'].shape}; "
f"included {int(d['included_mask'].sum())}; events {ev.event_type.value_counts().to_dict()}; "
f"DF/F 470 {g['dff'].min():.3f}..{g['dff'].max():.3f}, isosbestic {i['dff'].min():.3f}..{i['dff'].max():.3f}")
for k in ("DL31", "DL32"):
d = TWO[k]
g, i = d["Green"]["times"], d["GreenIsosbestic"]["times"]
print(f"{k}: isosbestic frames lead 470 nm frames by {1000 * np.median(g[:len(i)] - i[:len(g)]):.1f} ms;"
f" behaviour timestamps differ from the isosbestic series' regular frame times by up to"
f" {1000 * np.abs(d['beh']['t'] - i[:len(d['beh']['t'])]).max():.1f} ms, a drift over the session")
ev = d["events"]
st = ev[ev.event_type.isin(["Light", "Tone"])].sort_values("start_time")
lab = (st.event_type == "Tone").values.astype(int)
on = V.nearest_index(d["beh"]["t"], st.start_time.values) + 6 * lab
for color in ("Green", "GreenIsosbestic"):
x = d[color]["dff"]
accs = [round(decode_accuracy_fs(x, on, lab, sh, 18)[0], 2) for sh in (0, 90, -90)]
print(f" {color}: light against tone decoded, true onsets / +5 s / -5 s: {accs} (chance 0.5, {len(on)} presentations)")
842-221123: 470 nm (32000, 61) at 18.00 Hz from 180.76 s; isosbestic (32000, 61); behaviour samples 32000 from 180.73 s; ROI centroids (82, 2); included 61; events {'Lick': 3138, 'Reward': 47}; DF/F 470 -0.078..0.552, isosbestic -0.026..0.051
842-221118: 470 nm (32000, 61) at 18.00 Hz from 49.75 s; isosbestic (32000, 61); behaviour samples 32000 from 49.78 s; ROI centroids (82, 2); included 61; events {'Lick': 2203, 'Reward': 48}; DF/F 470 -0.072..0.657, isosbestic -0.048..0.053
DL31: 470 nm (8176, 52) at 17.98 Hz from 3.94 s; isosbestic (8177, 52); behaviour samples 8177 from 3.91 s; ROI centroids (73, 2); included 52; events {'Light': 10, 'Tone': 10}; DF/F 470 -0.026..0.058, isosbestic -0.019..0.064
DL32: 470 nm (8575, 44) at 17.98 Hz from 7.37 s; isosbestic (8576, 44); behaviour samples 8576 from 7.34 s; ROI centroids (73, 2); included 44; events {'Tone': 10, 'Light': 10}; DF/F 470 -0.029..0.101, isosbestic -0.025..0.058
DL31: isosbestic frames lead 470 nm frames by 28.0 ms; behaviour timestamps differ from the isosbestic series' regular frame times by up to 428.8 ms, a drift over the session
Green: light against tone decoded, true onsets / +5 s / -5 s: [np.float64(0.9), np.float64(0.4), np.float64(0.65)] (chance 0.5, 20 presentations)
GreenIsosbestic: light against tone decoded, true onsets / +5 s / -5 s: [np.float64(0.65), np.float64(0.55), np.float64(0.15)] (chance 0.5, 20 presentations)
DL32: isosbestic frames lead 470 nm frames by 28.0 ms; behaviour timestamps differ from the isosbestic series' regular frame times by up to 535.3 ms, a drift over the session
Green: light against tone decoded, true onsets / +5 s / -5 s: [np.float64(0.8), np.float64(0.6), np.float64(0.65)] (chance 0.5, 20 presentations)
GreenIsosbestic: light against tone decoded, true onsets / +5 s / -5 s: [np.float64(0.65), np.float64(0.5), np.float64(0.6)] (chance 0.5, 20 presentations)
Mouse 842 has 82 fibers per channel, 61 included; the paper's Figure 1C describes a 93-fiber array, and its code for the figure reads this mouse's session of 2022-11-23. Its traces start 180 s into the behaviour clock; no claim aligns them to events. DL31 and DL32 are recorded at 18 Hz with the two wavelengths interleaved (415 nm frames 28 ms before the 470 nm frames), as the Methods describe; the behaviour has one sample per frame pair (8177 for DL31), but its timestamps drift from the regular frame times stored for the series by up to half a second over the session, so the regular rate does not describe the frame times exactly. Events are therefore placed by their behaviour sample index, the frame index the paper's code uses, for both channels. Light and tone are decoded from the 470 nm signals better with these onsets (0.90 and 0.80) than with onsets shifted by 5 s (0.40 to 0.65); with 20 presentations per mouse the contrast is weak, so the alignment of these two sessions is consistent with the data but not shown. The bundle-image ROI centroids give one row per fiber of a channel.
dl15 = files[files.subject == "DL15"].copy()
dl15["date"] = pd.to_datetime(dl15.start.str[:10])
print(dl15[["session", "start", "frames", "rate_hz", "fibers", "included", "Reward", "Lick"]].to_string(index=False))
def load_dl15(path):
f = V.open_nwb(path)
dff, times, tab = V.read_dff(f)
beh = V.read_behavior(f)
ev = V.read_events(f)
n = len(beh["t"])
lick = np.zeros(n, bool)
for _, e in ev[ev.event_type == "Lick"].iterrows():
a, b = V.nearest_index(beh["t"], np.array([e.start_time, e.stop_time]))
lick[a:b] = True
rew = np.sort(V.nearest_index(beh["t"], ev[ev.event_type == "Reward"].start_time.values))
nxt = np.r_[rew[1:] - 1, n - 1]
cons = []
for a, b in zip(rew, nxt): # first lick after each reward and before the next (getRewInfo.m)
k = np.flatnonzero(lick[a + 1:b + 1])
if len(k):
cons.append(a + 1 + k[0])
cons = np.array([c for c in cons if c + 121 <= len(dff)], int)
x = dff[:, tab.included.values]
act = np.stack([x[c:c + 121] for c in cons]) if len(cons) else np.zeros((0, 121, x.shape[1]))
return {"path": path, "start": V._text(f["session_start_time"][()])[:10], "n_rew": len(rew), "act": act,
"n_cons": len(cons), "head": np.round(dff[:200], 6)}
with ThreadPoolExecutor(10) as ex:
D15 = list(ex.map(load_dl15, dl15.path))
print("\nSessions with 25 rewards or fewer:", int((dl15.Reward <= 25).sum()), "from", dl15[dl15.Reward <= 25].date.min().date(),
"to", dl15[dl15.Reward <= 25].date.max().date(), "; with more:", int((dl15.Reward > 25).sum()),
"from", dl15[dl15.Reward > 25].date.min().date(), "to", dl15[dl15.Reward > 25].date.max().date())
print("Rewards followed by a lick before the next reward (consumption), per session:", [s_["n_cons"] for s_ in D15])
session start frames rate_hz fibers included Reward Lick 210819 2021-08-19T11:44:21 9000 30.00 64 44 16 632 210820 2021-08-20T12:23:25 9000 30.00 64 44 17 340 210822 2021-08-22T13:57:37 9000 29.99 64 44 14 973 210823 2021-08-23T11:53:29 9000 29.99 64 44 15 718 210824 2021-08-24T16:28:14 9000 29.99 64 44 17 631 210825 2021-08-25T11:34:51 9000 30.00 64 44 15 519 210826 2021-08-26T09:19:47 9000 30.00 64 44 15 318 210827 2021-08-27T09:17:55 9000 29.99 64 44 18 738 210831 2021-08-31T09:23:58 10105 29.99 64 44 15 436 210902 2021-09-02T08:44:51 12000 29.99 64 44 17 571 210903 2021-09-03T09:20:07 10424 30.00 64 44 15 673 210904 2021-09-04T09:33:17 10463 30.00 64 44 15 848 210905 2021-09-05T09:32:46 7198 30.00 64 44 12 487 210906 2021-09-06T09:26:39 8259 29.99 64 44 12 456 210908 2021-09-08T09:33:29 5618 30.00 64 44 10 336 210909 2021-09-09T09:22:37 6501 30.00 64 44 10 448 210910b 2021-09-10T09:16:19 34778 29.99 64 44 50 2724 210913 2021-09-13T09:26:41 6961 29.99 64 44 10 461 210914 2021-09-14T09:18:08 31714 30.00 64 44 60 3030 210915 2021-09-15T09:32:33 32523 30.00 64 44 60 3200 210916 2021-09-16T09:36:54 34937 30.00 64 44 61 2507 210917 2021-09-17T10:30:34 28399 29.99 64 44 60 2998 210918 2021-09-18T09:16:27 34553 29.99 64 44 62 2387 210919 2021-09-19T09:36:38 30198 30.00 64 44 60 3237 210920 2021-09-20T09:33:27 32021 30.00 64 44 60 3550 210921 2021-09-21T09:43:44 30979 29.99 64 44 0 2 210923 2021-09-23T09:44:41 32277 30.00 64 44 60 3797 210924 2021-09-24T09:50:40 31558 29.99 64 44 60 3250 210925 2021-09-25T09:45:57 33991 30.00 64 44 60 3400 210926 2021-09-26T10:07:24 31919 29.99 64 44 61 3716 210927 2021-09-27T09:43:53 31469 29.99 64 44 60 2615 210928 2021-09-28T10:06:44 32746 30.00 64 44 61 2693 210929 2021-09-29T09:55:23 32563 29.99 64 44 60 2418 210930 2021-09-30T10:13:50 20839 29.99 64 44 44 1648 211001 2021-10-01T09:47:33 32155 30.00 64 44 61 2626 211004 2021-10-04T12:41:10 37509 29.99 64 44 70 2796 211005 2021-10-05T10:11:31 42000 29.99 64 44 70 2441 211006 2021-10-06T09:37:21 33106 30.00 64 44 60 1899 211007 2021-10-07T09:29:46 34165 29.99 64 44 60 1497 211008 2021-10-08T11:39:30 31063 30.00 64 44 60 1835 211010 2021-10-10T09:39:22 30127 29.99 64 44 60 1511 211011 2021-10-11T11:42:42 33966 29.99 64 44 60 2447 211012 2021-10-12T09:39:48 33087 30.00 64 44 60 2577 211013 2021-10-13T10:54:13 31680 30.00 64 44 60 2526 211014 2021-10-14T11:06:28 33304 30.00 64 44 60 2365 211015 2021-10-15T07:56:38 32234 29.99 64 44 60 1798
Sessions with 25 rewards or fewer: 18 from 2021-08-19 to 2021-09-21 ; with more: 28 from 2021-09-10 to 2021-10-15 Rewards followed by a lick before the next reward (consumption), per session: [16, 16, 14, 15, 17, 14, 15, 18, 15, 17, 15, 15, 12, 12, 10, 10, 50, 10, 60, 60, 61, 60, 62, 60, 60, 0, 60, 60, 60, 60, 60, 60, 60, 44, 60, 70, 70, 53, 60, 60, 60, 60, 60, 60, 60, 60]
DL15 has 44 included fibers, as the 44 fibers of Figure 1D. Sessions up to 2021-09-13 have 10 to 18
rewards and nothing else; from 2021-09-10 on most sessions have 44 to 70 rewards. The paper's code
treats sessions with more than 25 rewards as a task with predicted rewards and keeps only the rewards
of one task phase (virmen_phase), a variable the files do not hold (the paper's behaviour file of
2021-09-14 on Zenodo has the ViRMEn task variables, read while preparing, not here). The session of
2021-09-21 has no reward.
opto_paths = [p for p in paths if any(s in p for s in ("Opto2", "Opto4", "Opto9", "UG21"))]
def load_opto(path):
f = V.open_nwb(path)
color = "Red"
dff, times, tab = V.read_dff(f, color)
fs = 1 / np.median(np.diff(times))
stim, n_art = V.stimulated_fiber(dff, fs)
on, off = V.laser_epochs(dff[:, stim], fs) if stim is not None else (np.array([], int), np.array([], int))
return {"path": path, "subject": path.split("/")[0][4:], "session": path.split("_ses-")[1].split("_")[0],
"dff": dff, "fs": fs, "table": tab, "stim": stim, "n_artefact": n_art, "on": on, "off": off}
with ThreadPoolExecutor(12) as ex:
OPTO = list(ex.map(load_opto, opto_paths))
op = pd.DataFrame([{
"subject": o["subject"], "session": o["session"], "frames": o["dff"].shape[0], "fs": round(o["fs"], 2),
"fibers": o["dff"].shape[1], "included": int(o["table"].included.sum()),
"stim fiber (1-based)": None if o["stim"] is None else o["stim"] + 1,
"stim included": None if o["stim"] is None else bool(o["table"].included[o["stim"]]),
"stim location": None if o["stim"] is None else o["table"].location[o["stim"]],
"fibers with artefact": o["n_artefact"], "epochs": len(o["on"]),
"epoch s": f"{(np.median(o['off'] - o['on']) / o['fs']):.2f}" if len(o["on"]) else "",
"interval s": f"{np.diff(o['on']).min() / o['fs']:.0f}-{np.diff(o['on']).max() / o['fs']:.0f}" if len(o["on"]) > 1 else "",
"max dff stim": round(float(o["dff"][:, o["stim"]].max()), 2) if o["stim"] is not None else None} for o in OPTO])
print(op.to_string(index=False))
subject session frames fs fibers included stim fiber (1-based) stim included stim location fibers with artefact epochs epoch s interval s max dff stim Opto2 230314a 9000 29.99 41 22 13 True Caudoputamen 1 8 1.00 29-50 3.90 Opto2 230314b 9000 29.99 41 22 17 True Caudoputamen 1 8 1.00 30-55 4.19 Opto2 230314c 9000 29.99 41 22 38 True Caudoputamen 1 7 1.00 30-81 5.55 Opto2 230315a 9000 30.00 41 22 11 True Caudoputamen 1 7 1.00 29-50 2.72 Opto2 230315b 9000 30.00 41 22 7 True Caudoputamen 1 8 1.00 4-51 4.64 Opto2 230315c 9000 29.99 41 22 10 True Caudoputamen 1 9 1.00 9-53 1.39 Opto2 230315d 9000 29.99 41 22 12 True Caudoputamen 1 8 1.00 26-51 3.35 Opto2 230316a 9000 30.00 41 22 14 True Caudoputamen 1 8 1.00 32-56 2.24 Opto2 230316b 9000 30.00 41 22 15 False Primary somatosensory area upper limb layer 6a 1 8 1.00 33-55 2.23 Opto2 230316c 9000 29.99 41 22 16 False Caudoputamen 1 7 1.00 35-85 1.70 Opto2 230316d 9000 29.99 41 22 18 False Primary somatosensory area mouth layer 6a 1 8 1.00 29-61 4.37 Opto2 230316e 9000 30.00 41 22 19 False supra-callosal cerebral white matter 1 8 1.00 30-60 4.61 Opto2 230316f 9000 30.00 41 22 20 True Caudoputamen 1 7 1.00 32-61 1.80 Opto2 230316g 9000 29.99 41 22 22 False Primary somatosensory area barrel field layer 6a 1 8 1.00 28-55 4.23 Opto2 230320a 9000 29.99 41 22 23 True Caudoputamen 1 8 1.00 30-46 2.44 Opto2 230320b 9000 29.99 41 22 24 True Caudoputamen 1 9 1.00 25-47 4.17 Opto2 230320c 9000 29.99 41 22 25 False corpus callosum body 1 8 1.00 32-60 3.72 Opto2 230320d 9000 30.00 41 22 26 True Caudoputamen 1 9 1.00 27-50 1.25 Opto2 230320e 9000 29.99 41 22 27 False corpus callosum body 1 10 1.00 11-57 3.51 Opto2 230320f 9000 29.99 41 22 28 False Primary somatosensory area barrel field layer 6a 13 9 1.00 29-54 3.16 Opto2 230321a 9000 29.99 41 22 39 True Caudoputamen 1 9 1.00 27-59 3.30 Opto2 230321b 9000 29.99 41 22 33 True Caudoputamen 1 9 1.00 25-60 4.64 Opto2 230321c 9000 30.00 41 22 29 False Primary somatosensory area barrel field layer 6a 1 8 1.00 32-60 1.18 Opto2 230321d 9000 30.00 41 22 31 True Caudoputamen 1 8 1.00 31-59 4.91 Opto2 230321e 9000 29.99 41 22 32 True Caudoputamen 1 8 1.00 29-55 4.28 Opto2 230321f 9000 29.99 41 22 34 False Caudoputamen 1 8 1.00 26-52 4.05 Opto2 230321g 9000 29.99 41 22 36 True Caudoputamen 1 8 1.00 32-57 3.78 Opto2 230321h 9000 30.00 41 22 37 False 1 8 1.00 32-62 5.30 Opto2 230322a 9000 29.99 41 22 40 True Caudoputamen 1 8 1.00 30-60 5.43 Opto2 230322b 9000 30.00 41 22 41 True Caudoputamen 1 8 1.00 32-55 3.39 Opto4 230929a 9000 30.00 40 24 1 True Caudoputamen 24 5 1.00 34-59 13.28 Opto4 230929b 9000 30.00 40 24 1 True Caudoputamen 30 6 1.00 35-60 10.39 Opto4 231002a 9000 30.00 40 24 6 True Caudoputamen 17 5 1.00 52-59 7.77 Opto4 231002b 9000 29.99 40 24 5 False Caudoputamen 21 6 1.00 31-52 12.65 Opto4 231002c 9000 29.99 40 24 8 False Caudoputamen 15 7 1.00 27-60 7.53 Opto4 231003a 9000 29.99 40 24 10 False lateral ventricle 11 6 1.00 32-55 8.91 Opto4 231003b 9000 30.00 40 24 12 False 1 6 1.00 39-54 33.27 Opto4 231003c 9000 30.00 40 24 14 True Caudoputamen 13 6 1.00 34-60 14.31 Opto4 231004a 9000 29.99 40 24 16 True Caudoputamen 25 5 1.00 46-60 19.24 Opto4 231004b 9000 29.99 40 24 18 True Caudoputamen 3 6 1.00 35-57 13.39 Opto4 231004c 9000 29.99 40 24 20 False lateral ventricle 11 6 1.00 37-59 18.36 Opto4 231006a 9000 29.99 40 24 22 True Nucleus accumbens 6 6 1.00 14-57 23.90 Opto4 231006b 9000 29.99 40 24 24 True Caudoputamen 18 6 1.00 37-50 12.64 Opto9 230929a 9000 29.99 34 17 1 False Field CA3 8 6 1.00 35-60 7.10 Opto9 230929b 9000 29.99 34 17 3 False fimbria 12 5 1.00 32-60 4.00 Opto9 231002a 9000 29.99 34 17 4 False Field CA1 9 7 1.00 30-55 9.00 Opto9 231002b 9000 30.00 34 17 6 False stria terminalis 4 7 1.00 34-52 13.16 Opto9 231002c 9000 29.99 34 17 8 False 1 6 1.00 37-54 5.53 Opto9 231003a 9000 30.00 34 17 10 False alveus 16 7 1.00 34-49 5.81 Opto9 231003b 9000 29.99 34 17 12 False Globus pallidus external segment 5 7 1.00 8-60 10.67 Opto9 231003c 9000 29.99 34 17 14 True Caudoputamen 14 5 1.00 38-57 3.50 Opto9 231004a 9000 29.99 34 17 18 False lateral ventricle 23 6 1.00 41-56 7.77 Opto9 231004b 9000 30.00 34 17 18 False lateral ventricle 12 6 1.00 39-54 14.80 Opto9 231004c 9000 29.99 34 17 20 False Caudoputamen 11 6 1.00 32-59 15.28 Opto9 231006a 9000 30.00 34 17 20 False Caudoputamen 19 5 1.00 48-59 2.36 UG21 230321a 9000 30.00 30 22 12 True Caudoputamen 1 8 1.00 29-76 1.19 UG21 230321b 9000 29.99 30 22 18 True Caudoputamen 1 8 1.00 27-62 2.69 UG21 230321c 9000 30.00 30 22 24 True Caudoputamen 1 8 1.00 29-68 3.00
# Recordings deposited twice: compare the first 200 DF/F frames of every session loaded here
# (all but Grid-9, the only session of its mouse).
heads = [(f"{m} salient", d["dff_all"][:200]) for m, d in S.items()]
heads += [(k, d["Green"]["dff"][:200]) for k, d in TWO.items()]
heads += [(f"DL15 {s_['start']}", s_["head"]) for s_ in D15]
heads += [(f"{o['subject']} {o['session']}", o["dff"][:200]) for o in OPTO]
same = [(a, b) for i, (a, x) in enumerate(heads) for b, y in heads[i + 1:] if x.shape == y.shape and np.allclose(x, y, equal_nan=True)]
print("Sessions compared:", len(heads), "; pairs with identical first 200 frames:", same if same else "none")
Sessions compared: 115 ; pairs with identical first 200 frames: none
The files hold no stimulation times and no record of which fiber was stimulated. As in the paper's
code (getStimTimes.m), the laser epochs are found where the stimulated fiber's DF/F exceeds 0.4
(laser light leaking into the imaged fiber), here without the code's cross-check against the trigger
line, which the files do not hold. In Opto2 and UG21 exactly one fiber carries the artefact in each
session, except Opto2 2023-03-20f, where 13 fibers do, among them the midbrain fibers 1 to 5; it is
left out as a multi-site or midbrain stimulation. UG21's stimulated fibers are 12, 18 and 24, the
numbers of the paper's code folders for this mouse (null_12, null_18, null_24), and the first
three Opto2 sessions give fibers 13, 17 and 38, matching the code's stri_13, stri_17 and stri_38
recorded on 2023-03-14. There are 7 to 11 epochs of about 1 s, 30 to 60 s apart, as the Methods
state. In the tdTomato mice (Opto4, Opto9) the artefact appears on many fibers at once; the fiber with
the largest one is taken as the stimulated fiber, which is not always included. Opto9 has 12
sessions, so the files hold 25 tdTomato stimulations where the paper reports 26.
Salient-stimulus GLM (shared by Figures 4, 5 and S4 to S6)¶
For each of the 7 mice: overall velocity and acceleration from the treadmill; for each fiber the lags of the most positive and most negative significant ITI cross-correlation with each (within 500 ms); the GLM of the Methods on the -1 to 2 s windows; and the null of 10000 refits with random onsets and shuffled labels. The Methods' cross-correlation is used; the code's version, in which negative lags are not shifted, is shown as a variant of the point estimates.
rng_glm = np.random.default_rng(10840)
N_NULL = 10000
t0 = time.time()
G = {}
for m, d in S.items():
ev = d["events"]
st = ev[ev.event_type.isin(["Light", "Tone"])].sort_values("start_time")
lab = (st.event_type == "Tone").values.astype(int)
on = st.on.values + TONE_SHIFT * lab
fs = d["fs"]
vel, acc = V.overall_velocity(d["beh"], fs)
iti = V.iti_mask(len(vel), st.on.values, st.off.values, int(round(fs)))
lags = V.loco_lags(d["dff"], vel, acc, iti, 15)
glm = V.CueGLM(d["dff"], vel, acc, lags, fs)
mu, regr = glm.fit(on, lab, keep_trials=True)
null = V.cue_glm_null(glm, lab, N_NULL, rng_glm)
lo, hi = np.percentile(null, [0.5, 99.5], axis=0)
lags_code = V.loco_lags(d["dff"], vel, acc, iti, 15, code_indexing=True)
mu_code, _ = V.CueGLM(d["dff"], vel, acc, lags_code, fs).fit(on, lab)
G[m] = {"on": on, "lab": lab, "vel": vel, "acc": acc, "lags": lags, "mu": mu, "regr": regr, "lo": lo, "hi": hi,
"mu_code": mu_code, "glm": glm}
print(f"{m}: {len(on)} presentations ({(lab == 0).sum()} light, {(lab == 1).sum()} tone), {d['dff'].shape[1]} fibers; "
f"locomotion regressors per fiber {np.bincount(glm.n_loco, minlength=5).tolist()} (0..4); "
f"windows with NaN rows {int((~np.isfinite(glm._design(on, lab)[1]).all(1)).sum())}; {time.time() - t0:.0f} s")
rows = []
for m, g in G.items():
fib = S[m]["fib"]
for c, name in ((0, "light"), (1, "tone")):
up = V.runs_of(g["mu"][c] > g["hi"][c])
down = V.runs_of(g["mu"][c] < g["lo"][c])
val, sig = V.signed_peak(g["mu"][c], up, down)
valc, _ = V.signed_peak(g["mu_code"][c], up, down)
for r in range(len(fib)):
rows.append({"mouse": m, "fiber": r, "AP": fib.AP[r], "ML": fib.ML[r], "DV": fib.DV[r], "stim": name,
"value": val[r], "sig": bool(sig[r]), "up": bool(up[r]), "down": bool(down[r]),
"value_code_lags": valc[r]})
R = pd.DataFrame(rows)
W = R.pivot_table(index=["mouse", "fiber", "AP", "ML", "DV"], columns="stim",
values=["value", "sig", "up", "down", "value_code_lags"], aggfunc="first")
W.columns = [f"{a}_{b}" if b else a for a, b in W.columns]
W = W.reset_index()
for c in ("sig_light", "up_light", "down_light", "sig_tone", "up_tone", "down_tone"):
W[c] = W[c].astype(bool)
TAIL = -1.0
W["is_tail"] = W.AP <= TAIL
print("\nfibers", len(W), "; tail", int(W.is_tail.sum()), f"; GLM and null took {time.time() - t0:.0f} s")
for name in ("light", "tone"):
print(f"{name}: value range {W['value_' + name].min():.4f} to {W['value_' + name].max():.4f} DF/F; "
f"correlation with the code-lag variant r = {np.corrcoef(W['value_' + name], W['value_code_lags_' + name])[0, 1]:.4f}")
DL18: 58 presentations (14 light, 44 tone), 19 fibers; locomotion regressors per fiber [0, 0, 2, 16, 1] (0..4); windows with NaN rows 0; 49 s
DL20: 40 presentations (20 light, 20 tone), 53 fibers; locomotion regressors per fiber [0, 0, 47, 6, 0] (0..4); windows with NaN rows 0; 96 s
DL21: 36 presentations (18 light, 18 tone), 51 fibers; locomotion regressors per fiber [0, 0, 11, 36, 4] (0..4); windows with NaN rows 0; 137 s
DL23: 36 presentations (18 light, 18 tone), 27 fibers; locomotion regressors per fiber [0, 0, 8, 17, 2] (0..4); windows with NaN rows 0; 171 s
DL27: 42 presentations (21 light, 21 tone), 47 fibers; locomotion regressors per fiber [0, 1, 28, 18, 0] (0..4); windows with NaN rows 0; 216 s
DL29: 42 presentations (21 light, 21 tone), 45 fibers; locomotion regressors per fiber [0, 0, 19, 24, 2] (0..4); windows with NaN rows 2; 260 s
DL30: 42 presentations (21 light, 21 tone), 38 fibers; locomotion regressors per fiber [1, 1, 14, 21, 1] (0..4); windows with NaN rows 0; 302 s fibers 280 ; tail 35 ; GLM and null took 302 s light: value range -0.0208 to 0.0945 DF/F; correlation with the code-lag variant r = 0.9969 tone: value range -0.0487 to 0.0617 DF/F; correlation with the code-lag variant r = 0.9888
fig4-link-published-results¶
The dandiset's DF/F traces, fibers and stimulus times are those behind the paper's salient-stimulus analysis, so that its per-fiber responses can be linked to the paper's.
Outside input: the paper's glm_spline_results.mat (Zenodo Fig04.zip), read by range requests.
pub = sio.loadmat(io.BytesIO(V.zenodo_member("Fig04.zip", "Fig04/general_functions/glm_spline_results.mat")),
squeeze_me=True, struct_as_record=False)
rows = []
for m, g in G.items():
p = pub[m]
dff = S[m]["dff"]
row = {"mouse": m, "fibers dandiset": dff.shape[1], "fibers paper": p.cue1_roi_mu.shape[1]}
for c, name in ((0, "light"), (1, "tone")):
ons = g["on"][g["lab"] == c]
eta = np.nanmean(V.windows(dff, ons, -30, 60), 0)
pm = getattr(p, f"cue{c + 1}_roi_mu")
row[f"{name} trials"] = f"{len(ons)}/{getattr(p, f'cue{c + 1}_loco_regr_act').shape[2]}"
row[f"{name} raw max|diff|"] = float(np.abs(eta - pm).max()) if eta.shape == pm.shape else np.nan
row[f"{name} spline r"] = round(np.corrcoef(g["mu"][c].ravel(), getattr(p, f"cue{c + 1}_spline_mu").ravel())[0, 1], 4)
psig = getattr(p, f"cue{c + 1}_spline_sig")
up = V.runs_of(g["mu"][c] > g["hi"][c])
down = V.runs_of(g["mu"][c] < g["lo"][c])
row[f"{name} sig agree"] = round(float(np.mean((psig[:, 0] == up) & (psig[:, 1] == down))), 3)
rows.append(row)
link = pd.DataFrame(rows)
print(link.to_string(index=False))
agree_all = sum((link[f"{n} sig agree"] * link["fibers dandiset"]).sum() for n in ("light", "tone")) / (2 * link["fibers dandiset"].sum())
print("\nSignificance agreement over all fiber and stimulus pairs:", round(float(agree_all), 3))
print("Spline estimates pooled over mice: light r = %.4f, tone r = %.4f" % (
np.corrcoef(np.concatenate([G[m]["mu"][0].ravel() for m in G]), np.concatenate([pub[m].cue1_spline_mu.ravel() for m in G]))[0, 1],
np.corrcoef(np.concatenate([G[m]["mu"][1].ravel() for m in G]), np.concatenate([pub[m].cue2_spline_mu.ravel() for m in G]))[0, 1]))
# DL18: identify the paper's 14 tones (12 kHz) among the 44 by matching its per-trial activity.
p = pub["DL18"]
dff = S["DL18"]["dff"]
tones = G["DL18"]["on"][G["DL18"]["lab"] == 1]
win = V.windows(dff, tones, -30, 60) # tones x frames x fibers
pt = p.cue2_loco_regr_act # frames x fibers x trials
match, best_r, second_r = [], [], []
for k in range(pt.shape[2]):
r = np.array([np.ma.corrcoef(np.ma.masked_invalid(pt[:, :, k].ravel()), np.ma.masked_invalid(w.ravel()))[0, 1] for w in win])
o = np.argsort(r)[::-1]
match.append(o[0]); best_r.append(r[o[0]]); second_r.append(r[o[1]])
match = np.array(match)
print("\nDL18 tones matched to the paper's 14 trials:", (match + 1).tolist(), "; distinct:", len(set(match)))
print("match r: min %.3f; best competitor r: max %.3f" % (min(best_r), max(second_r)))
eta14 = np.nanmean(win[match], 0)
print("DL18 raw tone average over the 14 matched tones against the paper's: max |diff| = %.2e" % np.abs(eta14 - p.cue2_roi_mu).max())
DL18_12K = tones[np.sort(match)]
# The GLM and null of DL18 refitted with only these 14 tones, compared with the paper's results
g = G["DL18"]
keep = (g["lab"] == 0) | np.isin(g["on"], DL18_12K)
mu18, _ = g["glm"].fit(g["on"][keep], g["lab"][keep])
null18 = V.cue_glm_null(g["glm"], g["lab"][keep], N_NULL, np.random.default_rng(18))
lo18, hi18 = np.percentile(null18, [0.5, 99.5], axis=0)
up18, down18 = V.runs_of(mu18[1] > hi18[1]), V.runs_of(mu18[1] < lo18[1])
psig = p.cue2_spline_sig
print("DL18 tone with the 14 matched tones: spline r = %.4f (44 tones: %.4f); significance agreement %.3f (44 tones: %.3f)" % (
np.corrcoef(mu18[1].ravel(), p.cue2_spline_mu.ravel())[0, 1], link.loc[link.mouse == "DL18", "tone spline r"].iloc[0],
np.mean((psig[:, 0] == up18) & (psig[:, 1] == down18)), link.loc[link.mouse == "DL18", "tone sig agree"].iloc[0]))
n_pairs = 2 * link["fibers dandiset"].sum()
agree_14 = agree_all + (np.mean((psig[:, 0] == up18) & (psig[:, 1] == down18)) - link.loc[link.mouse == "DL18", "tone sig agree"].iloc[0]) \
* link.loc[link.mouse == "DL18", "fibers dandiset"].iloc[0] / n_pairs
print("Significance agreement over all fiber and stimulus pairs with DL18's 14 tones: %.3f" % agree_14)
print("Pooled tone spline r with DL18's 14 tones: %.4f" % np.corrcoef(
np.concatenate([(mu18 if m == "DL18" else G[m]["mu"])[1].ravel() for m in G]), np.concatenate([pub[m].cue2_spline_mu.ravel() for m in G]))[0, 1])
mouse fibers dandiset fibers paper light trials light raw max|diff| light spline r light sig agree tone trials tone raw max|diff| tone spline r tone sig agree DL18 19 19 14/14 0.0 0.9984 0.895 44/14 0.014614 0.9578 0.684 DL20 53 53 20/20 0.0 0.9997 1.000 20/20 0.000000 0.9997 1.000 DL21 51 51 18/18 0.0 0.9996 1.000 18/18 0.000000 0.9982 0.961 DL23 27 27 18/18 0.0 0.9998 1.000 18/18 0.000000 0.9948 0.963 DL27 47 47 21/21 0.0 0.9994 0.979 21/21 0.000000 0.9891 0.894 DL29 45 45 21/21 0.0 0.9794 0.822 21/21 0.000000 0.9860 0.867 DL30 38 38 21/21 0.0 0.9997 0.868 21/21 0.000000 1.0000 1.000 Significance agreement over all fiber and stimulus pairs: 0.936 Spline estimates pooled over mice: light r = 0.9957, tone r = 0.9922
DL18 tones matched to the paper's 14 trials: [1, 6, 8, 9, 11, 18, 21, 25, 26, 30, 34, 38, 39, 43] ; distinct: 14 match r: min 0.952; best competitor r: max 0.693 DL18 raw tone average over the 14 matched tones against the paper's: max |diff| = 0.00e+00
DL18 tone with the 14 matched tones: spline r = 0.9977 (44 tones: 0.9578); significance agreement 0.947 (44 tones: 0.684) Significance agreement over all fiber and stimulus pairs with DL18's 14 tones: 0.945 Pooled tone spline r with DL18's 14 tones: 0.9974
Result: reproduced. For every mouse the raw triggered averages computed from the dandiset equal the paper's published ones (largest difference 0.0) for light, and for tone in the six mice other than DL18, whose file holds 44 tones where the paper's results use 14. The GLM spline estimates correlate with the published ones at r = 0.9957 for light and 0.9922 for tone over all fibers, and the significance of increases and decreases agrees for 93.6% of fiber and stimulus pairs (lowest: DL18 tone 0.684, with 44 tones in place of 14; DL29 light 0.822). DL18's 14 tones are identified one to one (match r at least 0.952, best competing tone at most 0.693), and their raw average equals the published one. Refitted with these 14 tones, DL18's tone estimates correlate with the published ones at r = 0.9977 and significance agrees for 0.947 of its fibers; with them, agreement over all pairs is 94.5% and the pooled tone r 0.9974. Per mouse, r is below 0.99 for DL29 light (0.979) and for DL27 and DL29 tone (0.989, 0.986). The dandiset's traces, fiber order and stimulus times are those of the paper's analysis; only the tone frequency of DL18 is missing.
fig4e-light-increases-widespread¶
The light evokes transient increases in dopamine release across many striatal locations.
n = len(W)
up, down = int(W.up_light.sum()), int(W.down_light.sum())
val_up = int((W.sig_light & (W.value_light > 0)).sum())
val_down = int((W.sig_light & (W.value_light < 0)).sum())
print(f"Light: fibers with a significant increase {up}/{n} ({up / n:.0%}), with a significant decrease {down}/{n};"
f" by the signed largest value: {val_up} increases, {val_down} decreases, {n - val_up - val_down} not significant")
print(W.groupby("mouse").agg(fibers=("fiber", "size"), increase=("up_light", "sum"), decrease=("down_light", "sum")).to_string())
print("Null: each fiber is tested on 31 frames against a 0.5 to 99.5 percentile band; the fraction of fibers that"
" reach 3 consecutive frames by chance is estimated below from 200 of the null refits per mouse.")
fp, fp_any, frame_out, run3 = [], [], [], []
for m, g in G.items():
sub = V.cue_glm_null(g["glm"], g["lab"], 200, np.random.default_rng(1))
fp.append(np.mean([V.runs_of(sub[i, 0] > g["hi"][0]).mean() for i in range(200)]))
fp_any.append(np.mean([(V.runs_of(sub[i, 0] > g["hi"][0]) | V.runs_of(sub[i, 0] < g["lo"][0])).mean() for i in range(200)]))
out = (sub[:, 0] > g["hi"][0]) | (sub[:, 0] < g["lo"][0]) # refits x frames x fibers
frame_out.append(out.mean())
run3.append((out[:, :-2] & out[:, 1:-1] & out[:, 2:]).mean())
print("chance rate of a significant light increase per fiber (mean over mice): %.4f" % np.mean(fp))
print("chance rate of a significant increase or decrease per fiber: %.4f" % np.mean(fp_any))
print("Paper's p < 1e-6 is 0.01 ** 3, three independent frames outside the 99%% band. In the null refits: a frame is"
" outside the band with rate %.4f, and three given consecutive frames all are with rate %.4f (0.01 ** 3 = 1e-6)"
% (np.mean(frame_out), np.mean(run3)))
Light: fibers with a significant increase 240/280 (86%), with a significant decrease 38/280; by the signed largest value: 237 increases, 11 decreases, 32 not significant
fibers increase decrease
mouse
DL18 19 17 0
DL20 53 50 0
DL21 51 34 22
DL23 27 23 0
DL27 47 44 0
DL29 45 39 12
DL30 38 33 4
Null: each fiber is tested on 31 frames against a 0.5 to 99.5 percentile band; the fraction of fibers that reach 3 consecutive frames by chance is estimated below from 200 of the null refits per mouse.
chance rate of a significant light increase per fiber (mean over mice): 0.0237 chance rate of a significant increase or decrease per fiber: 0.0426 Paper's p < 1e-6 is 0.01 ** 3, three independent frames outside the 99% band. In the null refits: a frame is outside the band with rate 0.0104, and three given consecutive frames all are with rate 0.0037 (0.01 ** 3 = 1e-6)
Result: reproduced. 240 of 280 fibers (86%) have a significant increase to light and 38 a
significant decrease; by the signed largest value, 237 are increases and 11 decreases. Every mouse has
increases in most fibers (from 34 of 51 in DL21 to 50 of 53 in DL20), and decreases are concentrated in DL21 (22 of 51) and
DL29 (12 of 45). The test's chance rate of a significant increase per fiber, from null refits, is 0.024
(about 7 of 280 fibers), and 0.043 for an increase or a decrease. This is far above the paper's stated
p < 1e-6, which is 0.01 cubed: in the null refits three given consecutive frames are all outside the band
at a rate of 0.0037, not 1e-6, because neighbouring frames of the spline estimate are correlated, and the
test looks for such a run anywhere in the 31 frames. The procedure is the paper's (its published
significance calls agree with these for 93.6% of pairs, see fig4-link-published-results), so the
difference lies in how the paper states the rate.
fig4e-light-anterior-posterior-gradient¶
Light-evoked dopamine increases follow an amplitude gradient from the anterior striatum into the posterior tail, where increases are less prevalent and smaller.
def light_gradient(w, tail):
rho, p = stats.spearmanr(w.AP, w.value_light)
t = w.AP <= tail
table = [[int((w.up_light & t).sum()), int((~w.up_light & t).sum())],
[int((w.up_light & ~t).sum()), int((~w.up_light & ~t).sum())]]
_, p_frac = stats.fisher_exact(table, alternative="less")
_, p_amp = stats.mannwhitneyu(w.value_light[t], w.value_light[~t], alternative="less")
return {"tail AP <=": tail, "n tail": int(t.sum()), "rho AP-light": round(rho, 3), "p rho": p,
"frac up tail": round(table[0][0] / max(t.sum(), 1), 2), "frac up other": round(table[1][0] / (~t).sum(), 2),
"p frac": p_frac, "median tail": round(w.value_light[t].median(), 4),
"median other": round(w.value_light[~t].median(), 4), "p amp": p_amp}
rows = [light_gradient(W, t) for t in (TAIL, -0.8, -1.3)]
print(pd.DataFrame(rows).to_string(index=False))
print("\nPer mouse (tail AP <= -1.0):")
print(pd.DataFrame([{"mouse": m, **light_gradient(W[W.mouse == m], TAIL)} for m in G]).to_string(index=False))
r_w, p_w = V.within_group_rank_corr(W.AP, W.value_light, W.mouse, 5000, np.random.default_rng(4051))
print(f"\nVariant, within mice (AP and light value ranked within each mouse, permutation within mice): r = {r_w:.3f}, p = {p_w:.2g}")
rho_c, p_c = stats.spearmanr(W.AP, W.value_code_lags_light)
print(f"\nVariant, code's lags: rho(AP, light) = {rho_c:.3f}, p = {p_c:.2g}")
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5))
for a, s in zip(ax, ("light", "tone")):
a.scatter(W.AP, W["value_" + s], c=np.where(W["sig_" + s], "k", "0.7"), s=10)
a.axvline(TAIL, ls="--", c="r")
a.axhline(0, c="0.5", lw=0.5)
a.set_xlabel("AP (mm from bregma)")
a.set_ylabel(f"{s}: signed largest DF/F")
plt.tight_layout()
plt.show()
print("Plotted ranges: AP %.2f to %.2f; light %.4f to %.4f; tone %.4f to %.4f (all points shown, black = significant)"
% (W.AP.min(), W.AP.max(), W.value_light.min(), W.value_light.max(), W.value_tone.min(), W.value_tone.max()))
tail AP <= n tail rho AP-light p rho frac up tail frac up other p frac median tail median other p amp
-1.0 35 0.168 0.004721 0.83 0.86 0.381625 0.0055 0.0117 0.001285
-0.8 46 0.168 0.004721 0.85 0.86 0.498060 0.0079 0.0118 0.001165
-1.3 21 0.168 0.004721 0.81 0.86 0.350866 0.0055 0.0113 0.011113
Per mouse (tail AP <= -1.0):
mouse tail AP <= n tail rho AP-light p rho frac up tail frac up other p frac median tail median other p amp
DL18 -1.0 2 0.691 0.001047 0.50 0.94 0.204678 0.0008 0.0080 0.052632
DL20 -1.0 9 0.445 0.000847 0.89 0.95 0.434645 0.0078 0.0167 0.005081
DL21 -1.0 8 0.069 0.632862 0.75 0.65 0.827864 0.0111 0.0130 0.175635
DL23 -1.0 1 -0.428 0.025903 1.00 0.85 1.000000 0.0333 0.0140 0.925926
DL27 -1.0 5 0.635 0.000002 0.80 0.95 0.292014 0.0033 0.0120 0.003033
DL29 -1.0 3 -0.251 0.096711 1.00 0.86 1.000000 0.0159 0.0119 0.749824
DL30 -1.0 7 -0.121 0.468734 0.86 0.87 0.661493 0.0047 0.0054 0.302828
Variant, within mice (AP and light value ranked within each mouse, permutation within mice): r = 0.152, p = 0.013 Variant, code's lags: rho(AP, light) = 0.162, p = 0.0067
Plotted ranges: AP -2.12 to 1.76; light -0.0208 to 0.0945; tone -0.0487 to 0.0617 (all points shown, black = significant)
Result: partially reproduced. The light value rises from posterior to anterior (Spearman rho = 0.168, p = 0.0047) and is smaller in the tail (median 0.0055 against 0.0117 DF/F, p = 0.0013), but increases are not less prevalent there: 83% of tail fibers against 86% elsewhere have a significant increase (p = 0.38). The same holds with the tail at -0.8 or -1.3 mm and with the code's lags (rho = 0.162). The gradient is carried by DL18, DL20 and DL27 (rho 0.45 to 0.69); three mice have a negative point estimate (DL23 -0.43, DL29 -0.25, DL30 -0.12). Ranked within each mouse and tested by permuting within mice, so that differences between mice cannot drive it, the correlation is r = 0.152 (p = 0.013).
fig4f-tone-tail-increase-anterior-decrease¶
The tone evokes the strongest dopamine increases at sites in the posterior tail, while anterior to the tail it evokes primarily decreases.
def tone_pattern(w, tail, col="value_tone"):
t = w.AP <= tail
_, p_amp = stats.mannwhitneyu(w[col][t], w[col][~t], alternative="greater")
top = w.sort_values(col, ascending=False).head(10)
ant = w[~t & w.sig_tone]
n_dec, n_inc = int((ant[col] < 0).sum()), int((ant[col] > 0).sum())
p_bin = stats.binomtest(n_dec, n_dec + n_inc, 0.5, alternative="greater").pvalue if n_dec + n_inc else np.nan
return {"tail AP <=": tail, "n tail": int(t.sum()), "median tail": round(w[col][t].median(), 4),
"median anterior": round(w[col][~t].median(), 4), "p tail>anterior": p_amp,
"top10 in tail": int((top.AP <= tail).sum()), "anterior sig dec": n_dec, "anterior sig inc": n_inc,
"p binomial": p_bin}
print(pd.DataFrame([tone_pattern(W, t) for t in (TAIL, -0.8, -1.3)]).to_string(index=False))
print("\nTail fibers with a significant tone increase / decrease:", int((W.is_tail & W.up_tone).sum()), "/", int((W.is_tail & W.down_tone).sum()),
"of", int(W.is_tail.sum()), "; anterior:", int((~W.is_tail & W.up_tone).sum()), "/", int((~W.is_tail & W.down_tone).sum()), "of", int((~W.is_tail).sum()))
print("\nPer mouse:")
print(pd.DataFrame([{"mouse": m, **tone_pattern(W[W.mouse == m], TAIL)} for m in G]).to_string(index=False))
print("\nVariant, code's lags:", tone_pattern(W.assign(value_tone=W.value_code_lags_tone), TAIL))
# Variant: DL18 with only the 12 kHz tones identified from the paper's results (fit in fig4-link-published-results)
v18, s18 = V.signed_peak(mu18[1], up18, down18)
W12 = W.copy()
W12.loc[W12.mouse == "DL18", "value_tone"] = v18
W12.loc[W12.mouse == "DL18", "sig_tone"] = s18
print("Variant, DL18 with its 14 12-kHz tones (with inputs):", tone_pattern(W12, TAIL))
tail AP <= n tail median tail median anterior p tail>anterior top10 in tail anterior sig dec anterior sig inc p binomial
-1.0 35 0.0092 -0.0026 5.529544e-13 10 103 60 0.000472
-0.8 46 0.0068 -0.0028 8.175197e-16 10 102 50 0.000015
-1.3 21 0.0102 -0.0024 4.487610e-08 5 103 72 0.011534
Tail fibers with a significant tone increase / decrease: 32 / 0 of 35 ; anterior: 62 / 108 of 245
Per mouse:
mouse tail AP <= n tail median tail median anterior p tail>anterior top10 in tail anterior sig dec anterior sig inc p binomial
DL18 -1.0 2 0.0045 0.0013 9.356725e-02 2 6 6 6.127930e-01
DL20 -1.0 9 0.0176 -0.0025 9.194148e-04 5 17 18 6.320606e-01
DL21 -1.0 8 0.0379 -0.0173 1.507625e-07 6 36 3 1.804437e-08
DL23 -1.0 1 0.0579 -0.0073 3.703704e-02 1 17 6 1.734483e-02
DL27 -1.0 5 0.0032 -0.0016 1.254026e-02 2 5 8 8.665771e-01
DL29 -1.0 3 0.0079 -0.0025 2.889359e-03 3 16 8 7.579482e-02
DL30 -1.0 7 0.0092 0.0021 7.585266e-03 4 6 11 9.282684e-01
Variant, code's lags: {'tail AP <=': -1.0, 'n tail': 35, 'median tail': np.float64(0.0095), 'median anterior': np.float64(-0.0028), 'p tail>anterior': np.float64(4.2654031389401284e-13), 'top10 in tail': 10, 'anterior sig dec': 104, 'anterior sig inc': 59, 'p binomial': np.float64(0.0002641009771453017)}
Variant, DL18 with its 14 12-kHz tones (with inputs): {'tail AP <=': -1.0, 'n tail': 35, 'median tail': np.float64(0.0092), 'median anterior': np.float64(-0.0028), 'p tail>anterior': np.float64(4.196589080170338e-13), 'top10 in tail': 10, 'anterior sig dec': 108, 'anterior sig inc': 60, 'p binomial': np.float64(0.00013182582476392674)}
Result: reproduced. Tail fibers have larger tone values than anterior fibers (median 0.0092 against -0.0026 DF/F, p = 5.5e-13), all 10 of the largest tone values are in the tail, and 32 of 35 tail fibers have a significant increase with none decreasing. Anterior to the tail, 103 fibers with a significant tone response decrease and 60 increase (binomial p = 0.00047). The amplitude comparison and the anterior excess of decreases hold at every tail boundary (-0.8, -1.0 and -1.3 mm) and with the code's lags, but at -1.3 mm only 5 of the 10 largest tone values are in the tail, below the 6 that part (a) of the criterion needs. The result is the same with DL18 restricted to its 14 12-kHz tones (108 decreases against 60). The anterior decrease is not uniform across mice: it is significant in DL21 (36 against 3) and DL23 (17 against 6), while DL18, DL20, DL27 and DL30 have as many or more anterior increases.
fig4f-decreases-lateral-ventral¶
Tone-evoked decreases are slightly more prevalent in the anterior lateral and ventral striatum, with many sites showing both increases to the light and decreases to the tone.
A = W[~W.is_tail].copy()
A["tone_dec"] = A.sig_tone & (A.value_tone < 0)
rows = []
for coord in ("ML", "DV", "AP"):
_, p = stats.mannwhitneyu(A[coord][A.tone_dec], A[coord][~A.tone_dec], alternative="greater")
rows.append({"coordinate": coord, "median, tone decrease": A[coord][A.tone_dec].median(),
"median, other": A[coord][~A.tone_dec].median(), "p (decrease > other)": p})
print("Anterior fibers:", len(A), "; with a significant tone decrease:", int(A.tone_dec.sum()))
print(pd.DataFrame(rows).to_string(index=False))
both = int((A.sig_light & (A.value_light > 0) & A.tone_dec).sum())
print(f"\nAnterior fibers with a significant light increase and a significant tone decrease: {both}/{len(A)} ({both / len(A):.0%})")
for tail in (-0.8, -1.3):
B = W[W.AP > tail].copy()
B["tone_dec"] = B.sig_tone & (B.value_tone < 0)
ps = [stats.mannwhitneyu(B[c][B.tone_dec], B[c][~B.tone_dec], alternative="greater")[1] for c in ("ML", "DV")]
print(f"Variant tail AP <= {tail}: p ML = {ps[0]:.3g}, p DV = {ps[1]:.3g}")
Anterior fibers: 245 ; with a significant tone decrease: 103
coordinate median, tone decrease median, other p (decrease > other)
ML 2.07 1.965 0.425038
DV 3.66 3.300 0.003077
AP 0.52 0.075 0.000001
Anterior fibers with a significant light increase and a significant tone decrease: 79/245 (32%)
Variant tail AP <= -0.8: p ML = 0.169, p DV = 0.0105
Variant tail AP <= -1.3: p ML = 0.834, p DV = 0.00182
Result: partially reproduced. Among 245 anterior fibers, the 103 with a significant tone decrease lie more ventrally than the others (median DV 3.66 against 3.30 mm, p = 0.0031) but not more laterally (median ML 2.07 against 1.97 mm, p = 0.43); the ventral effect holds at both other tail boundaries and the lateral one at neither. They also lie more anteriorly (p = 1e-06). 79 anterior fibers (32%) have both a significant light increase and a significant tone decrease, so the second part of the statement holds by the 20% reading.
figs4ac-patterns-consistent-across-mice¶
The modality-specific release patterns are consistent within and across mice.
rows = []
for m in G:
w = W[W.mouse == m]
rl, pl = stats.spearmanr(w.AP, w.value_light)
rt, pt_ = stats.spearmanr(w.AP, w.value_tone)
rows.append({"mouse": m, "fibers": len(w), "AP range": f"{w.AP.min():.2f}..{w.AP.max():.2f}", "tail fibers": int(w.is_tail.sum()),
"rho AP-light": round(rl, 3), "p light": round(pl, 4), "rho AP-tone": round(rt, 3), "p tone": round(pt_, 4)})
pm = pd.DataFrame(rows)
print(pm.to_string(index=False))
print(f"\nMice with a positive AP-light correlation: {(pm['rho AP-light'] > 0).sum()}/7;"
f" with a negative AP-tone correlation: {(pm['rho AP-tone'] < 0).sum()}/7")
mouse fibers AP range tail fibers rho AP-light p light rho AP-tone p tone DL18 19 -1.15..1.45 2 0.691 0.0010 -0.609 0.0057 DL20 53 -1.52..1.45 9 0.445 0.0008 -0.578 0.0000 DL21 51 -1.98..1.30 8 0.069 0.6329 -0.629 0.0000 DL23 27 -1.35..1.76 1 -0.428 0.0259 -0.441 0.0213 DL27 47 -2.12..1.41 5 0.635 0.0000 -0.648 0.0000 DL29 45 -1.47..1.35 3 -0.251 0.0967 -0.695 0.0000 DL30 38 -1.85..1.12 7 -0.121 0.4687 -0.620 0.0000 Mice with a positive AP-light correlation: 4/7; with a negative AP-tone correlation: 7/7
Result: partially reproduced. The tone pattern is consistent across mice: AP and the tone value are negatively correlated in all 7 mice (rho -0.44 to -0.70, each p < 0.03). The light gradient is not: the correlation is positive in 4 of 7 mice (DL18, DL20, DL21, DL27; significant in three) and negative in DL23 (significant), DL29 and DL30.
figs4eh-patterns-not-movement¶
Locomotor changes to the stimuli are variable and bi-directional, while the dopamine release patterns to light and tone are consistent, so the patterns are stimulus and not movement driven.
rows, split = [], []
for m, g in G.items():
dff = S[m]["dff"]
for c, name in ((0, "light"), (1, "tone")):
ons = g["on"][g["lab"] == c]
acc = V.windows(g["acc"], ons, 0, 30)[:, :, 0]
k = np.nanargmax(np.abs(acc), axis=1)
sign = np.sign(acc[np.arange(len(ons)), k])
dv = np.nanmean(V.windows(g["vel"], ons, 0, 30)[:, :, 0], 1) - np.nanmean(V.windows(g["vel"], ons, -30, -1)[:, :, 0], 1)
rows.append({"mouse": m, "stim": name, "acc trials": int((sign > 0).sum()), "dec trials": int((sign < 0).sum()),
"mean velocity change m/s": round(dv.mean(), 4), "velocity up/down": f"{(dv > 0).sum()}/{(dv < 0).sum()}"})
for part, mask in (("acc", sign > 0), ("dec", sign < 0)):
if mask.sum() > 1:
mu = np.nanmean(V.windows(dff, ons[mask], 0, 30), 0)
v = np.where(np.abs(mu.min(0)) > np.abs(mu.max(0)), mu.min(0), mu.max(0))
# variant: each split's mean minus its own mean over the 1 s before onset
mub = mu - np.nanmean(V.windows(dff, ons[mask], -30, -1), (0, 1))
vb = np.where(np.abs(mub.min(0)) > np.abs(mub.max(0)), mub.min(0), mub.max(0))
else:
v = vb = np.full(dff.shape[1], np.nan)
for r in range(dff.shape[1]):
split.append({"mouse": m, "fiber": r, "stim": name, "part": part, "value": v[r], "value_base": vb[r]})
print(pd.DataFrame(rows).to_string(index=False))
SP0 = pd.DataFrame(split)
SP0["col"] = SP0.stim + "_" + SP0.part
for value_col, title in (("value", "Raw DF/F, as the code"), ("value_base", "Variant, pre-onset baseline subtracted in each split")):
print(f"\n{title}:")
SP = SP0.pivot_table(index=["mouse", "fiber"], columns="col", values=value_col, dropna=False).reset_index()
SP.columns.name = None
SP = SP.merge(W[["mouse", "fiber", "AP", "is_tail"]], on=["mouse", "fiber"])
for name in ("light", "tone"):
a, d = SP[f"{name}_acc"], SP[f"{name}_dec"]
ok = a.notna() & d.notna()
print(f"{name}: fibers with both splits {ok.sum()}; r(acc, dec) = {np.corrcoef(a[ok], d[ok])[0, 1]:.3f}")
for part in ("acc", "dec"):
v = SP[f"{name}_{part}"]
print(f" {part}: tail median {v[SP.is_tail].median():.4f}, anterior median {v[~SP.is_tail].median():.4f}")
mouse stim acc trials dec trials mean velocity change m/s velocity up/down DL18 light 4 10 -0.2114 2/12 DL18 tone 21 23 0.0166 24/20 DL20 light 5 15 -0.0732 6/14 DL20 tone 9 11 0.0194 13/7 DL21 light 12 6 0.0099 10/8 DL21 tone 17 1 0.0365 15/3 DL23 light 9 9 0.0202 14/4 DL23 tone 16 2 0.1205 17/1 DL27 light 10 11 -0.0357 8/13 DL27 tone 8 13 -0.0947 7/14 DL29 light 0 21 -0.3337 0/21 DL29 tone 6 15 -0.0964 5/16 DL30 light 14 7 -0.0014 13/8 DL30 tone 13 8 0.0241 13/8 Raw DF/F, as the code: light: fibers with both splits 235; r(acc, dec) = 0.873 acc: tail median 0.0069, anterior median 0.0130 dec: tail median 0.0082, anterior median 0.0128 tone: fibers with both splits 229; r(acc, dec) = 0.784 acc: tail median 0.0073, anterior median -0.0031 dec: tail median 0.0098, anterior median 0.0042 Variant, pre-onset baseline subtracted in each split: light: fibers with both splits 235; r(acc, dec) = 0.755 acc: tail median 0.0054, anterior median 0.0120 dec: tail median 0.0088, anterior median 0.0132 tone: fibers with both splits 229; r(acc, dec) = 0.823 acc: tail median 0.0069, anterior median -0.0043 dec: tail median 0.0101, anterior median -0.0024
Result: partially reproduced. Locomotion after the stimuli varies across mice and presentations:
the split into presentations followed by acceleration and by deceleration ranges from 0 and 21 (DL29
light) to 21 and 23 (DL18 tone), and the mean velocity change is positive in some mice and negative in
others for both stimuli. The fiber maps of the two splits correlate (light r = 0.873 over 235 fibers, tone
r = 0.784 over 229), and the light map has positive medians in both splits. The tone pattern holds after
acceleration (tail median 0.0073, anterior median -0.0031) but not after deceleration, where the anterior
median is positive (0.0042, tail 0.0098). This test uses raw DF/F, as the paper's code
(sal_stim_bubbles_loco.m) does, without a baseline and without the locomotion regressors of the GLM.
With each split's pre-onset mean subtracted (a variant), the maps correlate at r = 0.755 (light) and
0.823 (tone), and the tone pattern holds in both splits (anterior median -0.0043 after acceleration,
-0.0024 after deceleration; tail 0.0069 and 0.0101), which points to the level before onset
in the deceleration trials as the cause of the failure under the code's measure.
figs5-isosbestic-no-pattern¶
With isosbestic (415 nm) excitation, large fluorescence changes and spatial patterns to light and tone are not observed, ruling out a hemodynamic origin of the release patterns.
rng_s5 = np.random.default_rng(5415)
rows, S5, NULL5, VAR13 = [], [], [], []
for m in ("DL31", "DL32"):
d = TWO[m]
beh, ev = d["beh"], d["events"]
st = ev[ev.event_type.isin(["Light", "Tone"])].sort_values("start_time")
for color, label in (("Green", "470"), ("GreenIsosbestic", "415")):
ch = d[color]
fs = ch["rate"]
n_fr = len(ch["dff"])
b2 = {k: beh[k][:n_fr] for k in ("roll", "pitch", "angular")} # behaviour sample k = frame pair k
vel, acc = V.overall_velocity(b2, fs)
lab = (st.event_type == "Tone").values.astype(int)
on = V.nearest_index(beh["t"], st.start_time.values) + 6 * lab # 6-frame tone delay at 18 Hz (code)
off = V.nearest_index(beh["t"], st.stop_time.values)
iti = V.iti_mask(len(vel), on - 6 * lab, off, int(round(fs)))
lags = V.loco_lags(ch["dff"], vel, acc, iti, 9)
glm = V.CueGLM(ch["dff"], vel, acc, lags, fs)
mu, _ = glm.fit(on, lab)
null = V.cue_glm_null(glm, lab, N_NULL, rng_s5)
lo, hi = np.percentile(null, [0.5, 99.5], axis=0)
for c, name in ((0, "light"), (1, "tone")):
up, down = V.runs_of(mu[c] > hi[c]), V.runs_of(mu[c] < lo[c])
val, sig = V.signed_peak(mu[c], up, down)
fib = ch["fib"]
rows.append({"mouse": m, "nm": label, "stim": name, "fibers": len(val), "trials": int((lab == c).sum()),
"median |value|": round(float(np.median(np.abs(val))), 4), "significant": int(sig.sum()),
"rho AP-value": round(stats.spearmanr(fib.AP, val)[0], 3)})
S5.append(pd.DataFrame({"mouse": m, "nm": label, "stim": name, "fiber": np.arange(len(val)), "value": val}))
# reference: |signed largest value| of the null refits (no stimulus), first 1000 refits
nv = np.where(np.abs(null[:1000, c].min(1)) > np.abs(null[:1000, c].max(1)), null[:1000, c].min(1), null[:1000, c].max(1))
NULL5.append({"mouse": m, "nm": label, "stim": name, "abs": np.abs(nv).ravel()})
# variant: 13 spline functions over the 1 s at 18 Hz (the Methods' 13 df), same lags
glm13 = V.CueGLM(ch["dff"], vel, acc, lags, fs, knot_intervals=10)
mu13, _ = glm13.fit(on, lab)
null13 = V.cue_glm_null(glm13, lab, N_NULL, rng_s5)
lo13, hi13 = np.percentile(null13, [0.5, 99.5], axis=0)
for c, name in ((0, "light"), (1, "tone")):
v13, s13 = V.signed_peak(mu13[c], V.runs_of(mu13[c] > hi13[c]), V.runs_of(mu13[c] < lo13[c]))
VAR13.append({"mouse": m, "nm": label, "stim": name, "abs": np.abs(v13), "sig": int(s13.sum())})
s5 = pd.DataFrame(rows)
print(s5.to_string(index=False))
allv = pd.concat(S5)
for name in ("light", "tone"):
a = allv[(allv.stim == name) & (allv.nm == "470")].value.abs().median()
b = allv[(allv.stim == name) & (allv.nm == "415")].value.abs().median()
n470 = s5[(s5.stim == name) & (s5.nm == "470")].significant.sum()
n415 = s5[(s5.stim == name) & (s5.nm == "415")].significant.sum()
print(f"{name}: median |value| 470 nm {a:.4f}, 415 nm {b:.4f} (ratio {b / a:.2f}); significant fibers 470 {n470}, 415 {n415}")
print("\nReference: median |value| of the null refits (random onsets), per channel")
for name in ("light", "tone"):
ref = {nm: np.median(np.concatenate([x["abs"] for x in NULL5 if x["stim"] == name and x["nm"] == nm])) for nm in ("470", "415")}
a = allv[(allv.stim == name) & (allv.nm == "470")].value.abs().median()
b = allv[(allv.stim == name) & (allv.nm == "415")].value.abs().median()
print(f"{name}: null median 470 nm {ref['470']:.4f}, 415 nm {ref['415']:.4f}; median above the null 470 nm {a - ref['470']:.4f},"
f" 415 nm {b - ref['415']:.4f}")
print("\nVariant, 13 spline functions at 18 Hz:")
for name in ("light", "tone"):
a = np.median(np.concatenate([x["abs"] for x in VAR13 if x["stim"] == name and x["nm"] == "470"]))
b = np.median(np.concatenate([x["abs"] for x in VAR13 if x["stim"] == name and x["nm"] == "415"]))
n470 = sum(x["sig"] for x in VAR13 if x["stim"] == name and x["nm"] == "470")
n415 = sum(x["sig"] for x in VAR13 if x["stim"] == name and x["nm"] == "415")
print(f"{name}: median |value| 470 nm {a:.4f}, 415 nm {b:.4f} (ratio {b / a:.2f}); significant fibers 470 {n470}, 415 {n415}")
mouse nm stim fibers trials median |value| significant rho AP-value DL31 470 light 52 10 0.0045 26 -0.250 DL31 470 tone 52 10 0.0027 9 -0.489 DL31 415 light 52 10 0.0021 1 -0.063 DL31 415 tone 52 10 0.0019 0 -0.189 DL32 470 light 44 10 0.0034 19 0.200 DL32 470 tone 44 10 0.0027 9 -0.434 DL32 415 light 44 10 0.0019 1 -0.261 DL32 415 tone 44 10 0.0018 1 0.163 light: median |value| 470 nm 0.0042, 415 nm 0.0020 (ratio 0.49); significant fibers 470 45, 415 2 tone: median |value| 470 nm 0.0027, 415 nm 0.0018 (ratio 0.68); significant fibers 470 18, 415 1 Reference: median |value| of the null refits (random onsets), per channel light: null median 470 nm 0.0025, 415 nm 0.0019; median above the null 470 nm 0.0017, 415 nm 0.0001 tone: null median 470 nm 0.0025, 415 nm 0.0019; median above the null 470 nm 0.0002, 415 nm -0.0001 Variant, 13 spline functions at 18 Hz: light: median |value| 470 nm 0.0043, 415 nm 0.0024 (ratio 0.55); significant fibers 470 36, 415 1 tone: median |value| 470 nm 0.0030, 415 nm 0.0021 (ratio 0.70); significant fibers 470 15, 415 0
Result: partially reproduced. With 415 nm excitation almost no fiber responds significantly (light 2 against 45 at 470 nm; tone 1 against 18 over 96 fibers of two mice), and the AP correlations seen at 470 nm for tone (-0.49 and -0.43) are not present at 415 nm (-0.19 and 0.16). The median absolute value at 415 nm is 0.49 of that at 470 nm for light, below the criterion's half, but 0.68 for tone. The largest-magnitude value has a floor: with random onsets its median is 0.0025 (470 nm) and 0.0019 DF/F (415 nm). The tone median at 470 nm (0.0027) is only 0.0002 above that floor, so the tone ratio compares two values near the noise and cannot fall to half; for light the 470 nm median is 0.0017 above the floor and the 415 nm median 0.0001. With 13 spline functions over the 1 s, as the Methods state for the main GLM, in place of the code's 9 at 18 Hz, the ratios are 0.55 (light) and 0.70 (tone), with 36 and 15 significant fibers at 470 nm against 1 and 0 at 415 nm. The 470 and 405 nm sessions of Figure S5A-D are not in the dandiset.
fig5ab-response-latencies¶
Latencies to half-maximum are shorter for dopamine increases to light than to tone, and shorter for increases than decreases to tone (light increase 221.6 ms, SD 42.5; tone increase 305.4 ms, SD 100.0; tone decrease 428.3 ms, SD 197.1).
FS = 30.0
lat = []
for m, g in G.items():
for c, name in ((0, "light"), (1, "tone")):
up = V.runs_of(g["mu"][c] > g["hi"][c])
down = V.runs_of(g["mu"][c] < g["lo"][c])
for r in range(g["mu"].shape[2]):
if up[r]:
lat.append({"mouse": m, "fiber": r, "stim": name, "kind": "increase", "frames": V.half_max_latency(g["mu"][c][:, r], 1)})
if down[r]:
lat.append({"mouse": m, "fiber": r, "stim": name, "kind": "decrease", "frames": V.half_max_latency(g["mu"][c][:, r], -1)})
L = pd.DataFrame(lat)
L["ms (Methods)"] = 1000 * L.frames / FS
L["ms (code, +1 frame)"] = 1000 * (L.frames + 1) / FS
summ = L.groupby(["stim", "kind"]).agg(n=("frames", "size"), no_crossing=("frames", lambda x: int(x.isna().sum())),
mean_ms=("ms (Methods)", "mean"), sd_ms=("ms (Methods)", "std"),
mean_code_ms=("ms (code, +1 frame)", "mean"), sd_code_ms=("ms (code, +1 frame)", "std"))
print(summ.round(1).to_string())
print("Paper: light increase 221.6 (SD 42.5); tone increase 305.4 (SD 100.0); tone decrease 428.3 (SD 197.1) ms")
print("\nPer mouse, mean ms (Methods):")
print(L.groupby(["mouse", "stim", "kind"])["ms (Methods)"].agg(["size", "mean"]).round(1).unstack(["stim", "kind"]).to_string())
fig, ax = plt.subplots(figsize=(6, 3))
for (s, k), col in ((("light", "increase"), "C1"), (("tone", "increase"), "C3"), (("tone", "decrease"), "C0")):
ax.hist(L[(L.stim == s) & (L.kind == k)]["ms (Methods)"].dropna(), bins=np.arange(0, 1034, 1000 / 30), alpha=0.5, label=f"{s} {k}", color=col)
ax.set_xlabel("latency to half max (ms, from onset frame)")
ax.legend()
plt.show()
print("Plotted range %.0f to %.0f ms; all values within the axes" % (L["ms (Methods)"].min(), L["ms (Methods)"].max()))
n no_crossing mean_ms sd_ms mean_code_ms sd_code_ms
stim kind
light decrease 38 0 570.2 198.0 603.5 198.0
increase 240 0 187.5 38.1 220.8 38.1
tone decrease 108 0 393.2 196.5 426.5 196.5
increase 94 1 271.7 98.7 305.0 98.7
Paper: light increase 221.6 (SD 42.5); tone increase 305.4 (SD 100.0); tone decrease 428.3 (SD 197.1) ms
Per mouse, mean ms (Methods):
size mean
stim light tone light light tone light
kind increase decrease increase decrease increase decrease increase decrease
mouse
DL18 17.0 6.0 8.0 NaN 200.0 605.6 275.0 NaN
DL20 50.0 19.0 27.0 NaN 209.3 450.9 227.2 NaN
DL21 34.0 37.0 11.0 22.0 167.6 290.1 187.9 450.0
DL23 23.0 17.0 7.0 NaN 176.8 402.0 176.2 NaN
DL27 44.0 7.0 12.0 NaN 174.2 285.7 422.2 NaN
DL29 39.0 16.0 11.0 12.0 189.7 516.7 342.4 700.0
DL30 33.0 6.0 18.0 4.0 190.9 405.6 282.4 841.7
Plotted range 33 to 967 ms; all values within the axes
Result: partially reproduced. The order of the latencies holds (light increase 187.5 ms < tone increase 271.7 ms < tone decrease 393.2 ms, counted from the onset frame as the Methods define it), and the tone decrease is within tolerance, but the light increase (expected 221.6 +- 22.2) and the tone increase (305.4 +- 30.5) are below it. Counted as the paper's code counts, one frame later, the means are 220.8, 305.0 and 426.5 ms, within 2 ms of the paper's values, with SDs of 38.1, 98.7 and 196.5 ms against 42.5, 100.0 and 197.1. The difference from the paper is the code's one-frame offset. One tone increase has no half-max crossing before its peak.
fig5g-light-latency-gradient¶
Dopamine release latencies to light follow an anterior-posterior gradient, with faster increases in the anterior striatum, consistent with wave-like dynamics.
def trial_latencies(m, c, kind):
"""Per presentation and fiber, half-max latency (frames) of the trial activity with locomotion
regressed out, where the fiber's mean response is significant and the trial exceeds the null band
for 3 consecutive frames (rise_dip_time.m)."""
g = G[m]
tau = g["glm"].tau
w = (tau >= 0) & (tau <= 30)
act = g["regr"][g["lab"] == c][:, w, :] # trials x 31 x fibers
up = V.runs_of(g["mu"][c] > g["hi"][c])
down = V.runs_of(g["mu"][c] < g["lo"][c])
out = np.full((act.shape[0], act.shape[2]), np.nan)
for t in range(act.shape[0]):
for r in range(act.shape[2]):
x = act[t, :, r]
if kind == "increase" and up[r] and V.runs_of(x > g["hi"][c][:, r])[0]:
out[t, r] = V.half_max_latency(x, 1)
if kind == "decrease" and down[r] and V.runs_of(x < g["lo"][c][:, r])[0]:
out[t, r] = V.half_max_latency(x, -1)
return out
REL = []
for m in G:
rise_l = trial_latencies(m, 0, "increase")
rise_t = trial_latencies(m, 1, "increase")
dip_t = trial_latencies(m, 1, "decrease")
either_t = np.fmin(rise_t, dip_t)
for name, mat in (("light rise", rise_l), ("tone either", either_t), ("tone rise", rise_t), ("tone dip", dip_t)):
rel = mat - np.nanmin(mat, axis=1, keepdims=True)
for r in range(mat.shape[1]):
n = int(np.isfinite(rel[:, r]).sum())
REL.append({"mouse": m, "fiber": r, "measure": name, "n_trials": n,
"rel_ms": 1000 * np.nanmean(rel[:, r]) / FS if n else np.nan,
"abs_ms": 1000 * np.nanmean(mat[:, r]) / FS if n else np.nan})
REL = pd.DataFrame(REL).merge(W[["mouse", "fiber", "AP", "ML", "DV"]], on=["mouse", "fiber"])
def gradient(measure, sign_word, rng=None):
d = REL[(REL.measure == measure) & REL.rel_ms.notna()]
rho, p = stats.spearmanr(d.AP, d.rel_ms)
print(f"{measure}: fibers with a latency {len(d)} (trials per fiber median {d.n_trials.median():.0f});"
f" Spearman rho(AP, relative latency) = {rho:.3f}, p = {p:.2g} ({sign_word})")
if rng is not None:
r_w, p_w = V.within_group_rank_corr(d.AP, d.rel_ms, d.mouse, 5000, rng)
print(f" variant, within mice (ranks within each mouse, permutation within mice): r = {r_w:.3f}, p = {p_w:.2g}")
for c in ("ML", "DV"):
r2, p2 = stats.spearmanr(d[c], d.rel_ms)
print(f" rho({c}, relative latency) = {r2:.3f}, p = {p2:.2g}")
return d, rho, p
d_light, rho_l, p_l = gradient("light rise", "negative = anterior faster", np.random.default_rng(5071))
fig, ax = plt.subplots(figsize=(5, 3))
ax.scatter(d_light.AP, d_light.rel_ms, s=10)
ax.set_xlabel("AP (mm)")
ax.set_ylabel("mean relative latency, light (ms)")
plt.show()
print("Plotted ranges: AP %.2f to %.2f; latency %.0f to %.0f ms" % (d_light.AP.min(), d_light.AP.max(), d_light.rel_ms.min(), d_light.rel_ms.max()))
light rise: fibers with a latency 240 (trials per fiber median 18); Spearman rho(AP, relative latency) = -0.285, p = 7e-06 (negative = anterior faster)
variant, within mice (ranks within each mouse, permutation within mice): r = -0.223, p = 0.0014 rho(ML, relative latency) = 0.263, p = 3.8e-05 rho(DV, relative latency) = -0.005, p = 0.93
Plotted ranges: AP -1.96 to 1.45; latency 49 to 442 ms
Result: reproduced. Over the 240 fibers with a light latency, the mean latency relative to the fastest fiber of each presentation is shorter anteriorly (Spearman rho(AP, latency) = -0.285, p = 7e-06). Latency also increases laterally (rho(ML) = 0.263, p = 3.8e-05) and not with depth (p = 0.93). Ranked within each mouse and tested by permuting within mice, which removes differences between mice and the dependence of relative latency on a mouse's number of fibers, the gradient remains (r = -0.223, p = 0.0014).
fig5h-tone-latency-gradient¶
Dopamine release latencies to tone follow a different gradient, faster in the posterior striatum than in the anterior striatum.
d_tone, rho_t, p_t = gradient("tone either", "positive = posterior faster", np.random.default_rng(5081))
for meas in ("tone rise", "tone dip"):
gradient(meas, "component")
t = d_tone.AP <= TAIL
print(f"Tone relative latency: tail median {d_tone.rel_ms[t].median():.0f} ms (n {t.sum()}), anterior median {d_tone.rel_ms[~t].median():.0f} ms (n {(~t).sum()})")
tone either: fibers with a latency 196 (trials per fiber median 18); Spearman rho(AP, relative latency) = 0.254, p = 0.00033 (positive = posterior faster)
variant, within mice (ranks within each mouse, permutation within mice): r = 0.214, p = 0.003 rho(ML, relative latency) = -0.271, p = 0.00012 rho(DV, relative latency) = 0.159, p = 0.026 tone rise: fibers with a latency 94 (trials per fiber median 18); Spearman rho(AP, relative latency) = -0.070, p = 0.5 (component) rho(ML, relative latency) = -0.189, p = 0.068 rho(DV, relative latency) = -0.066, p = 0.53 tone dip: fibers with a latency 106 (trials per fiber median 17); Spearman rho(AP, relative latency) = -0.214, p = 0.028 (component) rho(ML, relative latency) = -0.076, p = 0.44 rho(DV, relative latency) = 0.101, p = 0.3 Tone relative latency: tail median 308 ms (n 32), anterior median 368 ms (n 164)
Result: reproduced. Over the 196 fibers with a tone latency (the earlier of increase and decrease), relative latency is shorter posteriorly (rho(AP, latency) = 0.254, p = 0.00033; tail median 308 ms against 368 ms anterior). Taken separately, neither increases (rho = -0.07, p = 0.5) nor decreases (rho = -0.21, p = 0.028, the opposite sign) show this gradient, so it comes from the tail's fast increases and the anterior striatum's slower decreases, as the paper states for Figure S6A-B. Ranked within each mouse and tested by permuting within mice, the gradient remains (r = 0.214, p = 0.003).
figs6cd-gradients-individual-mice¶
The spatial latency gradients are also present in individual mice.
rows = []
for m in G:
row = {"mouse": m}
for meas, key in (("light rise", "light"), ("tone either", "tone")):
d = REL[(REL.mouse == m) & (REL.measure == meas) & REL.rel_ms.notna()]
row[f"{key} fibers"] = len(d)
if len(d) >= 5:
rho, p = stats.spearmanr(d.AP, d.rel_ms)
row[f"{key} rho"], row[f"{key} p"] = round(rho, 3), round(p, 4)
rows.append(row)
pmice = pd.DataFrame(rows)
print(pmice.to_string(index=False))
print(f"\nLight: negative in {(pmice['light rho'] < 0).sum()} of {pmice['light rho'].notna().sum()} mice with >= 5 fibers;"
f" tone: positive in {(pmice['tone rho'] > 0).sum()} of {pmice['tone rho'].notna().sum()}")
mouse light fibers light rho light p tone fibers tone rho tone p DL18 17 -0.659 0.0040 14 0.459 0.0985 DL20 50 -0.251 0.0784 43 0.453 0.0023 DL21 34 -0.000 0.9986 47 0.110 0.4625 DL23 23 0.246 0.2585 24 0.436 0.0333 DL27 44 -0.646 0.0000 17 -0.459 0.0641 DL29 39 0.074 0.6556 27 0.195 0.3288 DL30 33 -0.302 0.0873 24 0.121 0.5737 Light: negative in 4 of 7 mice with >= 5 fibers; tone: positive in 6 of 7
Result: partially reproduced. The tone gradient is positive in 6 of 7 mice (significant in DL20 and DL23); the light gradient is negative in 4 of 7 (significant in DL18 and DL27), with DL21 at 0 and DL23 and DL29 positive.
fig1c-isosbestic-control¶
In a D1-cre mouse expressing GCaMP7f in D1 spiny projection neurons, fluorescence under 470 nm excitation shows Ca2+ transients across fibers that are absent under 405 nm (isosbestic) excitation.
for key in ("842-221123", "842-221118"):
d = TWO[key]
g, i = d["Green"]["dff"], d["GreenIsosbestic"]["dff"]
n = min(len(g), len(i))
ratio = g[:n].std(0) / i[:n].std(0)
print(f"{key}: fibers {len(ratio)}; SD ratio 470/405 median {np.median(ratio):.2f}, range {ratio.min():.2f} to {ratio.max():.2f};"
f" ratio > 1 in {(ratio > 1).sum()}/{len(ratio)} ({(ratio > 1).mean():.0%}); 470 SD median {np.median(g.std(0)):.4f},"
f" 405 SD median {np.median(i.std(0)):.4f}")
big = (g[:n] > 5 * np.median(i[:n].std(0))).mean(0)
print(f" fraction of frames above 5 x the median 405 SD: 470 median {np.median(big):.4f};"
f" 405 median {np.median((i[:n] > 5 * np.median(i[:n].std(0))).mean(0)):.4f}")
big405 = (i[:n] > 5 * np.median(i[:n].std(0))).mean(0)
print(f" 405 traces with no such frame: {(big405 == 0).sum()}/{len(big405)}; 405 largest fraction {big405.max():.4f};"
f" 470 range {big.min():.4f} to {big.max():.4f}")
842-221123: fibers 61; SD ratio 470/405 median 3.01, range 1.56 to 19.41; ratio > 1 in 61/61 (100%); 470 SD median 0.0105, 405 SD median 0.0034 fraction of frames above 5 x the median 405 SD: 470 median 0.0667; 405 median 0.0000 405 traces with no such frame: 11/61; 405 largest fraction 0.0118; 470 range 0.0013 to 0.3590 842-221118: fibers 61; SD ratio 470/405 median 2.84, range 1.46 to 9.91; ratio > 1 in 61/61 (100%); 470 SD median 0.0097, 405 SD median 0.0034 fraction of frames above 5 x the median 405 SD: 470 median 0.0576; 405 median 0.0000 405 traces with no such frame: 39/61; 405 largest fraction 0.0139; 470 range 0.0005 to 0.3790
Result: reproduced. In all 61 included fibers of mouse 842 the DF/F under 470 nm varies more than under 405 nm (SD ratio median 3.01, range 1.56 to 19.41). Frames above 5 times the median 405 nm SD make up a median 6.7% of each 470 nm trace (range 0.13% to 35.9%) and a median 0.0% of each 405 nm trace (at most 1.2%; 11 of 61 traces have none). The other session gives the same (median ratio 2.84, 61 of 61 above 1).
figs2d-correlation-distance¶
Pairwise correlations of D1 SPN Ca2+ signals are generally low but variable, and drop off significantly with the distance between fiber tips.
d = TWO["842-221123"]
g = d["Green"]["dff"]
fib = d["Green"]["fib"]
cent = d["centroids"][d["included_mask"]]
C = np.corrcoef(g.T)
iu = np.triu_indices(len(fib), 1)
tip = fib[["ML", "AP", "DV"]].values
pairs = pd.DataFrame({"corr": C[iu],
"euc": np.linalg.norm(tip[:, None] - tip[None], axis=2)[iu],
"AP": np.abs(tip[:, None, 1] - tip[None, :, 1])[iu],
"ML": np.abs(tip[:, None, 0] - tip[None, :, 0])[iu],
"DV": np.abs(tip[:, None, 2] - tip[None, :, 2])[iu],
"surface": np.linalg.norm(cent[:, None] - cent[None], axis=2)[iu]})
r, p = stats.pearsonr(pairs["corr"], pairs.euc)
print(f"pairs {len(pairs)}; pairwise r median {pairs['corr'].median():.3f}, IQR {pairs['corr'].quantile(0.25):.3f} to"
f" {pairs['corr'].quantile(0.75):.3f}, range {pairs['corr'].min():.3f} to {pairs['corr'].max():.3f}")
print(f"correlation of pairwise r with tip distance: r = {r:.3f}, p = {p:.2g}; Spearman {stats.spearmanr(pairs['corr'], pairs.euc)[0]:.3f}")
q = pairs.euc <= pairs.euc.quantile(0.25)
print(f"closest quarter of pairs (<= {pairs.euc.quantile(0.25):.2f} mm): mean r {pairs['corr'][q].mean():.3f}; farthest quarter mean r"
f" {pairs['corr'][pairs.euc >= pairs.euc.quantile(0.75)].mean():.3f}")
d2 = TWO["842-221118"]
C2 = np.corrcoef(d2["Green"]["dff"].T)
print(f"Variant, session 2022-11-18: median r {np.median(C2[iu]):.3f}; r with distance {stats.pearsonr(C2[iu], pairs.euc)[0]:.3f}"
f" (p = {stats.pearsonr(C2[iu], pairs.euc)[1]:.2g})")
pairs 1830; pairwise r median 0.163, IQR -0.025 to 0.361, range -0.678 to 0.904 correlation of pairwise r with tip distance: r = -0.190, p = 2.8e-16; Spearman -0.227 closest quarter of pairs (<= 1.22 mm): mean r 0.293; farthest quarter mean r 0.133 Variant, session 2022-11-18: median r 0.192; r with distance -0.238 (p = 5.3e-25)
Result: reproduced. Pairwise correlations of the 61 fibers are low (median r = 0.163, IQR -0.025 to 0.361) and variable (range -0.678 to 0.904), and fall with distance between tips (r = -0.190 across 1830 pairs, p = 2.8e-16; mean r 0.293 in the closest quarter of pairs against 0.133 in the farthest). The other session gives the same (r = -0.238). The p-value treats pairs as independent, as the code does.
figs2ef-no-bundle-crosstalk¶
After accounting for distance in the brain, signal correlations do not depend on the distance between fibers on the bundle surface, so nearby fibers in the bundle do not cross-talk.
X = sm.add_constant(pairs[["euc", "AP", "ML", "DV", "surface"]])
fit = sm.OLS(pairs["corr"], X).fit()
print(fit.summary2().tables[1].round(5).to_string())
print(f"\nbundle-surface distance (pixels): coefficient {fit.params['surface']:.2e} per pixel, 95% CI"
f" {fit.conf_int().loc['surface', 0]:.2e} to {fit.conf_int().loc['surface', 1]:.2e}, p = {fit.pvalues['surface']:.3g}")
print(f"correlation of tip distance with surface distance: r = {stats.pearsonr(pairs.euc, pairs.surface)[0]:.3f}")
fit2 = sm.OLS(C2[iu], X).fit()
print(f"Variant, session 2022-11-18: surface coefficient {fit2.params['surface']:.2e}, p = {fit2.pvalues['surface']:.3g}")
near = pairs.surface <= pairs.surface.quantile(0.1)
print(f"Pairs adjacent on the bundle (closest 10%, <= {pairs.surface.quantile(0.1):.0f} px): mean r {pairs['corr'][near].mean():.3f},"
f" others {pairs['corr'][~near].mean():.3f}")
Coef. Std.Err. t P>|t| [0.025 0.975] const 0.31225 0.01629 19.16984 0.00000 0.28031 0.34420 euc -0.27795 0.03674 -7.56447 0.00000 -0.35001 -0.20588 AP 0.23901 0.02406 9.93508 0.00000 0.19183 0.28619 ML 0.04001 0.02385 1.67801 0.09352 -0.00675 0.08678 DV 0.10464 0.02329 4.49279 0.00001 0.05896 0.15032 surface -0.00003 0.00012 -0.29025 0.77166 -0.00027 0.00020 bundle-surface distance (pixels): coefficient -3.44e-05 per pixel, 95% CI -2.67e-04 to 1.98e-04, p = 0.772 correlation of tip distance with surface distance: r = 0.437 Variant, session 2022-11-18: surface coefficient -9.75e-05, p = 0.411 Pairs adjacent on the bundle (closest 10%, <= 42 px): mean r 0.222, others 0.156
Result: reproduced. In the paper's linear model, the distance between fibers on the bundle image has no effect once tip distances are included (coefficient -3.4e-05 per pixel, 95% CI -2.7e-04 to 2.0e-04, p = 0.77), while tip distance does (Euclidean coefficient -0.278, p < 1e-5). The other session gives p = 0.41. Fibers adjacent on the bundle correlate slightly more (0.222 against 0.156), which the model attributes to their tips also being closer (bundle and tip distances correlate at r = 0.437).
fig1d-chronic-reward-stability¶
Dopamine release to unpredicted water reward is reliable and of comparable amplitude over at least 12 weeks after implantation (median over 44 fibers, normalized to week 5).
first = min(pd.Timestamp(s["start"]) for s in D15)
def weekly(sessions):
rows = []
for s in sessions:
if len(s["act"]):
wk = 5 + (pd.Timestamp(s["start"]) - first).days // 7
rows.append((wk, s["act"]))
weeks = sorted({w for w, _ in rows})
means = {w: np.concatenate([a for ww, a in rows if ww == w]).mean(0) for w in weeks}
counts = {w: sum(len(a) for ww, a in rows if ww == w) for w in weeks}
init = means[weeks[0]][:31].max(0)
tab = pd.DataFrame([{"week": w, "rewards": counts[w], "median": np.median(means[w][:31].max(0) / init),
"q25": np.percentile(means[w][:31].max(0) / init, 25), "q75": np.percentile(means[w][:31].max(0) / init, 75),
"median peak dff": np.median(means[w][:31].max(0))} for w in weeks])
return tab
main = weekly([s for s in D15 if s["n_rew"] <= 25])
allr = weekly(D15)
print("Sessions with 25 rewards or fewer (main):")
print(main.round(3).to_string(index=False))
print("\nAll rewards of all sessions (variant):")
print(allr.round(3).to_string(index=False))
fig, ax = plt.subplots(figsize=(5, 3))
ax.errorbar(main.week, main["median"], yerr=[main["median"] - main.q25, main.q75 - main["median"]], fmt="o-k", label="<= 25 rewards")
ax.errorbar(allr.week + 0.1, allr["median"], yerr=[allr["median"] - allr.q25, allr.q75 - allr["median"]], fmt="s--", c="0.5", label="all rewards")
ax.axhspan(0.7, 1.3, color="g", alpha=0.1)
ax.set_xlabel("week (first session = week 5)")
ax.set_ylabel("peak / week-5 peak, median over fibers")
ax.legend()
plt.show()
print("Plotted range of medians and quartiles: %.2f to %.2f" % (min(main.q25.min(), allr.q25.min()), max(main.q75.max(), allr.q75.max())))
Sessions with 25 rewards or fewer (main):
week rewards median q25 q75 median peak dff
5 92 1.000 1.000 1.000 0.018
6 48 0.903 0.701 1.029 0.015
7 81 0.932 0.708 1.137 0.015
8 20 1.196 0.923 1.327 0.017
All rewards of all sessions (variant):
week rewards median q25 q75 median peak dff
5 92 1.000 1.000 1.000 0.018
6 48 0.903 0.701 1.029 0.015
7 81 0.932 0.708 1.137 0.015
8 190 0.767 0.464 0.951 0.012
9 303 0.627 0.408 0.763 0.011
10 420 0.609 0.477 0.805 0.011
11 297 0.773 0.641 0.925 0.013
12 360 0.485 0.331 0.593 0.008
13 120 0.392 0.282 0.563 0.008
Plotted range of medians and quartiles: 0.28 to 1.33
Result: partially reproduced. In the sessions whose rewards are all unpredicted, the median peak reward response of the 44 fibers stays between 0.90 and 1.20 of its week-5 value through week 8 (0.903, 0.932, 1.196), within the 0.7 to 1.3 band, but these sessions end at week 8; the 12 weeks of the claim cannot be tested because the later sessions mix predicted and unpredicted rewards and the files do not say which is which. With all rewards of all sessions the median falls to 0.77 in week 8, 0.61 in week 10 and 0.39 in week 13, which is expected if predicted rewards are included and does not test the claim.
fig6-inclusion-counts¶
In the mouse with ChR2 in dopamine axons and rDAm3.0 in the striatum, 24 recording locations lie within or on the border of the striatum; 18 of the 24 stimulated locations evoke a significant response at another location, and 22 of the 24 recording locations respond to at least one stimulation.
rng_opto = np.random.default_rng(6060)
for o in OPTO:
if o["stim"] is not None and len(o["off"]):
o["resp"] = V.stim_response(o["dff"], o["on"], o["off"], o["fs"], rng_opto, laser_trace=o["dff"][:, o["stim"]])
o2 = [o for o in OPTO if o["subject"] == "Opto2" and o["n_artefact"] == 1]
inc = o2[0]["table"].included.values
print("Opto2 sessions with one stimulated fiber:", len(o2), "; left out:", [o["session"] for o in OPTO if o["subject"] == "Opto2" and o["n_artefact"] != 1])
print("Recording locations in or bordering the striatum (included):", int(inc.sum()))
stim_inc = [o for o in o2 if inc[o["stim"]]]
print("Sessions whose stimulated fiber is included:", len(stim_inc), "; stimulated fibers (1-based):", sorted(o["stim"] + 1 for o in stim_inc))
print("Stimulated fibers not included:", sorted(o["stim"] + 1 for o in o2 if not inc[o["stim"]]))
def counts(sessions, inc_mask):
good_stim, responders = [], set()
for o in sessions:
sig = o["resp"]["sig"] & inc_mask
others = sig.copy()
others[o["stim"]] = False
if others.sum() >= 1:
good_stim.append(o)
responders |= set(np.flatnonzero(sig))
return good_stim, responders
good, responders = counts(stim_inc, inc)
print(f"\nStimulated locations with a significant response at another included location: {len(good)} of {len(stim_inc)}"
f" ({len(good) / len(stim_inc):.2f}; paper 18 of 24, 0.75)")
print(f"Included recording locations with a significant response to at least one stimulation: {len(responders)} of {int(inc.sum())}")
resp_other = set()
for o in stim_inc:
s_ = o["resp"]["sig"] & inc
s_[o["stim"]] = False
resp_other |= set(np.flatnonzero(s_))
print(f"Variant, counting only responses to another fiber's stimulation: {len(resp_other)} of {int(inc.sum())}"
f" (dropped: {sorted(r + 1 for r in responders - resp_other)})")
print("Per session: stimulated fiber, epochs, significant other fibers (increase, any):")
print(pd.DataFrame([{"session": o["session"], "stim": o["stim"] + 1, "epochs": len(o["off"]),
"sig up others": int((o["resp"]["sig_up"] & inc).sum() - o["resp"]["sig_up"][o["stim"]]),
"sig any others": int((o["resp"]["sig"] & inc).sum() - o["resp"]["sig"][o["stim"]]),
"stim fiber sig": bool(o["resp"]["sig"][o["stim"]])} for o in stim_inc]).to_string(index=False))
# Variant: the code's further exclusions (fibers 8, 9 of 35 um and broken fiber 21; 1-based)
inc_code = inc.copy()
inc_code[[7, 8, 20]] = False
good_c, resp_c = counts([o for o in stim_inc if inc_code[o["stim"]]], inc_code)
print(f"\nVariant, code's exclusions: locations {int(inc_code.sum())}; stimulated with a response elsewhere {len(good_c)} of"
f" {len([o for o in stim_inc if inc_code[o['stim']]])}; responding recording locations {len(resp_c)}")
# Variant: fibers labelled striatal but not included (candidates for "on the border")
tabO = o2[0]["table"]
striatal = tabO.location.isin(["Caudoputamen", "Striatum", "Nucleus accumbens"]).values
inc_b = inc | striatal
print(f"Variant, included or labelled striatal: locations {int(inc_b.sum())} (added fibers {(np.flatnonzero(inc_b & ~inc) + 1).tolist()},"
f" labels {tabO.location[inc_b & ~inc].tolist()})")
stim_b = [o for o in o2 if inc_b[o["stim"]]]
good_b, resp_b = counts(stim_b, inc_b)
print(f" stimulated sessions {len(stim_b)}; stimulated with a response elsewhere {len(good_b)} of {len(stim_b)};"
f" responding recording locations {len(resp_b)} of {int(inc_b.sum())}")
Opto2 sessions with one stimulated fiber: 29 ; left out: ['230320f'] Recording locations in or bordering the striatum (included): 22 Sessions whose stimulated fiber is included: 19 ; stimulated fibers (1-based): [7, 10, 11, 12, 13, 14, 17, 20, 23, 24, 26, 31, 32, 33, 36, 38, 39, 40, 41] Stimulated fibers not included: [15, 16, 18, 19, 22, 25, 27, 29, 34, 37] Stimulated locations with a significant response at another included location: 14 of 19 (0.74; paper 18 of 24, 0.75) Included recording locations with a significant response to at least one stimulation: 20 of 22 Variant, counting only responses to another fiber's stimulation: 20 of 22 (dropped: []) Per session: stimulated fiber, epochs, significant other fibers (increase, any): session stim epochs sig up others sig any others stim fiber sig 230314a 13 8 9 9 True 230314b 17 8 6 6 True 230314c 38 7 5 5 True 230315a 11 7 1 1 False 230315b 7 8 16 16 True 230315c 10 9 0 0 False 230315d 12 8 0 0 False 230316a 14 8 5 5 True 230316f 20 7 0 0 True 230320a 23 8 1 1 True 230320b 24 9 1 1 True 230320d 26 9 0 0 False 230321a 39 9 9 9 True 230321b 33 9 4 4 True 230321d 31 8 13 13 True 230321e 32 8 16 16 True 230321g 36 8 18 18 True 230322a 40 8 0 0 False 230322b 41 8 2 2 True Variant, code's exclusions: locations 19; stimulated with a response elsewhere 14 of 19; responding recording locations 18 Variant, included or labelled striatal: locations 24 (added fibers [16, 34], labels ['Caudoputamen', 'Caudoputamen']) stimulated sessions 21; stimulated with a response elsewhere 16 of 21; responding recording locations 22 of 24
Result: partially reproduced. 22 Opto2 fibers are flagged included (expected 24 +- 2.4), and 20 of
them respond to at least one stimulation (expected 22 +- 2.2), both within tolerance; 14 of the 19
stimulated included fibers evoke a significant response at another included fiber, below the expected 18
+- 1.8. As a fraction, 14 of 19 (0.74) is close to the paper's 18 of 24 (0.75); the denominators differ because
fewer included fibers were stimulated in the files. Counting only responses to another fiber's
stimulation, so that a stimulated fiber's own photoswitching does not count, still gives 20 responding
locations. Counting the two fibers labelled Caudoputamen but not flagged (16 and 34) gives exactly the
paper's 24 locations and 22 responding ones, and 16 of 21 stimulated fibers with a response elsewhere; the
files hold 29 single-fiber sessions, of which only 21 stimulate a striatal fiber, so 18 of 24 cannot be
reached. The included flag appears to mark fibers within the striatum, not those on its border.
fig6gh-release-decays-with-distance¶
Stimulating dopamine axons through one fiber evokes release at neighbouring non-stimulated sites that drops off sharply with distance in all three dimensions.
rng_fit = np.random.default_rng(6061)
tabO = o2[0]["table"]
pos = tabO[["ML", "AP", "DV"]].values
resp_fibers = np.array(sorted(responders))
pts = []
for o in good:
for r in resp_fibers:
if r != o["stim"]:
dv = pos[r] - pos[o["stim"]]
pts.append({"stim": o["stim"] + 1, "rec": r + 1, "resp": o["resp"]["resp"][r], "sig": bool(o["resp"]["sig"][r]),
"dML": dv[0], "dAP": dv[1], "dDV": dv[2], "dist": np.linalg.norm(dv)})
P6 = pd.DataFrame(pts)
print(f"points {len(P6)} ({P6.stim.nunique()} stimulation sites x {P6.rec.nunique()} recording sites); paper: 378 (18 x 21)")
a, b = V.exp_fit(P6.dist.values, P6.resp.values)
boot = []
sites = P6.stim.unique()
for _ in range(1000):
pick = rng_fit.choice(sites, len(sites), replace=True)
sub = pd.concat([P6[P6.stim == s] for s in pick])
try:
boot.append(V.exp_fit(sub.dist.values, sub.resp.values)[1])
except RuntimeError:
boot.append(np.nan)
boot = np.array(boot)
print(f"exponential fit: a = {a:.4f} DF/F, b = {b:.3f} per mm; 95% bootstrap interval of b {np.nanpercentile(boot, 2.5):.3f} to"
f" {np.nanpercentile(boot, 97.5):.3f} (failed fits {int(np.isnan(boot).sum())} of 1000)")
for ax_ in ("dML", "dAP", "dDV"):
rho, p = stats.spearmanr(P6[ax_].abs(), P6.resp)
print(f"{ax_}: Spearman rho(|distance|, response) = {rho:.3f}, p = {p:.2g}")
print(f"response range {P6.resp.min():.4f} to {P6.resp.max():.4f}; significant points {int(P6.sig.sum())}; median distance of"
f" significant {P6.dist[P6.sig].median():.2f} mm, of others {P6.dist[~P6.sig].median():.2f} mm")
bins = np.linspace(P6.dist.min(), P6.dist.max(), 12)
P6["bin"] = np.digitize(P6.dist, bins[1:-1])
print(P6.groupby("bin").agg(dist=("dist", "mean"), n=("resp", "size"), mean_resp=("resp", "mean")).round(4).to_string())
fig, ax = plt.subplots(figsize=(5, 3))
ax.scatter(P6.dist, P6.resp, s=8, c=np.where(P6.sig, "r", "0.6"))
xx = np.linspace(0, P6.dist.max(), 100)
ax.plot(xx, a * np.exp(b * xx), "k")
ax.set_xlabel("distance from stimulated fiber (mm)")
ax.set_ylabel("mean DF/F 0-0.5 s after offset")
plt.show()
print("All points shown; ranges above.")
points 266 (14 stimulation sites x 20 recording sites); paper: 378 (18 x 21)
exponential fit: a = 0.0147 DF/F, b = -1.096 per mm; 95% bootstrap interval of b -1.325 to -0.890 (failed fits 0 of 1000)
dML: Spearman rho(|distance|, response) = -0.299, p = 6.9e-07
dAP: Spearman rho(|distance|, response) = -0.268, p = 9.2e-06
dDV: Spearman rho(|distance|, response) = -0.387, p = 6.4e-11
response range -0.0008 to 0.0205; significant points 106; median distance of significant 0.95 mm, of others 1.57 mm
dist n mean_resp
bin
0 0.5400 28 0.0079
1 0.7422 24 0.0066
2 0.9196 27 0.0061
3 1.0894 36 0.0045
4 1.2751 24 0.0029
5 1.4714 31 0.0036
6 1.6445 29 0.0017
7 1.8320 31 0.0018
8 2.0290 18 0.0013
9 2.1888 11 0.0016
10 2.3774 7 0.0017
All points shown; ranges above.
Result: reproduced. Over 266 pairs of 14 stimulated and 20 recording fibers (the paper has 378, 18 by 21), the response falls with distance as a = 0.0147 DF/F, b = -1.096 per mm (95% bootstrap interval -1.325 to -0.890), and with the absolute distance along each axis (Spearman rho -0.30 ML, -0.27 AP, -0.39 DV; all p < 1e-5). The binned mean falls from 0.0079 DF/F at 0.54 mm to 0.0013 to 0.0017 beyond 2 mm, and significant responses lie closer (median 0.95 mm) than the others (1.57 mm).
fig6g-dorsal-spread¶
Evoked release drops off less with distance dorsal to the stimulated fiber than ventral to it, consistent with dorsally extending dopamine axon branches.
def two_rates(df):
out = {}
for name, m in (("dorsal", df.dDV < 0), ("ventral", df.dDV > 0)):
out[name] = V.exp_fit(df.dist[m].values, df.resp[m].values)[1]
return out
rates = two_rates(P6)
print(f"dorsal sites {int((P6.dDV < 0).sum())}, ventral sites {int((P6.dDV > 0).sum())}")
print(f"decay rate b: dorsal {rates['dorsal']:.3f}, ventral {rates['ventral']:.3f} per mm")
diffs = []
for _ in range(1000):
pick = rng_fit.choice(sites, len(sites), replace=True)
sub = pd.concat([P6[P6.stim == s] for s in pick])
try:
rr = two_rates(sub)
diffs.append(rr["dorsal"] - rr["ventral"])
except RuntimeError:
diffs.append(np.nan)
diffs = np.array(diffs)
print(f"difference dorsal - ventral {rates['dorsal'] - rates['ventral']:.3f}; 95% bootstrap interval {np.nanpercentile(diffs, 2.5):.3f}"
f" to {np.nanpercentile(diffs, 97.5):.3f} (failed {int(np.isnan(diffs).sum())})")
for name, m in (("dorsal", P6.dDV < 0), ("ventral", P6.dDV > 0)):
rho, p = stats.spearmanr(P6.dDV[m].abs(), P6.resp[m])
print(f"{name}: Spearman rho(|dDV|, response) = {rho:.3f}, p = {p:.2g}; mean response {P6.resp[m].mean():.4f}")
dorsal sites 156, ventral sites 109 decay rate b: dorsal -0.971, ventral -1.403 per mm
difference dorsal - ventral 0.432; 95% bootstrap interval -0.094 to 0.969 (failed 0) dorsal: Spearman rho(|dDV|, response) = -0.369, p = 2.1e-06; mean response 0.0043 ventral: Spearman rho(|dDV|, response) = -0.510, p = 1.5e-08; mean response 0.0035
Result: partially reproduced. Release decays more slowly dorsal to the stimulated fiber (b = -0.971 per mm, 156 pairs) than ventral to it (b = -1.403, 109 pairs), but the 95% bootstrap interval of the difference (-0.094 to 0.969) includes 0.
fig6-tdtomato-control-no-response¶
In control mice expressing tdTomato in place of rDAm3.0, single-fiber stimulation evokes no significant fluorescence increase at non-stimulated sites (0 of 48 sites, 26 stimulations in 2 mice), so the release is not a hemodynamic effect.
rows = []
for o in OPTO:
if o["subject"] not in ("Opto4", "Opto9", "Opto2") or "resp" not in o:
continue
if o["subject"] == "Opto2" and o["n_artefact"] != 1:
continue
inc_o = o["table"].included.values.copy()
oth = inc_o.copy()
oth[o["stim"]] = False
up = o["resp"]["sig_up"] & oth
rows.append({"subject": o["subject"], "session": o["session"], "stim": o["stim"] + 1, "stim included": bool(inc_o[o["stim"]]),
"artefact fibers": o["n_artefact"], "epochs": len(o["off"]), "non-stim included": int(oth.sum()),
"sig increase": int(up.sum()), "sig any": int((o["resp"]["sig"] & oth).sum()),
"max resp sig inc": round(float(o["resp"]["resp"][up].max()), 4) if up.any() else np.nan,
"median resp others": round(float(np.median(o["resp"]["resp"][oth])), 4)})
T6 = pd.DataFrame(rows)
print(T6[T6.subject != "Opto2"].to_string(index=False))
for s in ("Opto4", "Opto9", "Opto2"):
t = T6[T6.subject == s]
sites = set()
for o in OPTO:
if o["subject"] == s and "resp" in o and (s != "Opto2" or o["n_artefact"] == 1):
oth = o["table"].included.values.copy()
oth[o["stim"]] = False
sites |= set(np.flatnonzero(o["resp"]["sig_up"] & oth))
print(f"{s}: sessions {len(t)}; site-stimulation pairs tested {t['non-stim included'].sum()}; with a significant increase"
f" {t['sig increase'].sum()}; distinct sites with one {len(sites)}; median response {t['median resp others'].median():.4f}")
subject session stim stim included artefact fibers epochs non-stim included sig increase sig any max resp sig inc median resp others Opto4 230929a 1 True 24 5 23 0 0 NaN -0.0002 Opto4 230929b 1 True 30 6 23 0 0 NaN -0.0001 Opto4 231002a 6 True 17 5 23 0 0 NaN 0.0002 Opto4 231002b 5 False 21 6 24 0 0 NaN -0.0001 Opto4 231002c 8 False 15 7 24 0 0 NaN -0.0007 Opto4 231003a 10 False 11 6 24 0 0 NaN -0.0000 Opto4 231003b 12 False 1 6 24 0 0 NaN -0.0012 Opto4 231003c 14 True 13 6 23 0 0 NaN 0.0003 Opto4 231004a 16 True 25 5 23 0 0 NaN -0.0002 Opto4 231004b 18 True 3 6 23 0 0 NaN -0.0005 Opto4 231004c 20 False 11 6 24 0 0 NaN -0.0002 Opto4 231006a 22 True 6 6 23 0 0 NaN -0.0005 Opto4 231006b 24 True 18 6 23 0 0 NaN -0.0009 Opto9 230929a 1 False 8 6 17 0 0 NaN 0.0022 Opto9 230929b 3 False 12 5 17 0 0 NaN -0.0002 Opto9 231002a 4 False 9 7 17 0 0 NaN 0.0002 Opto9 231002b 6 False 4 7 17 0 0 NaN 0.0003 Opto9 231002c 8 False 1 6 17 0 0 NaN -0.0009 Opto9 231003a 10 False 16 7 17 0 1 NaN -0.0003 Opto9 231003b 12 False 5 7 17 0 0 NaN 0.0004 Opto9 231003c 14 True 14 5 16 0 0 NaN -0.0006 Opto9 231004a 18 False 23 6 17 0 0 NaN -0.0003 Opto9 231004b 18 False 12 6 17 0 0 NaN 0.0003 Opto9 231004c 20 False 11 6 17 0 0 NaN 0.0013 Opto9 231006a 20 False 19 5 17 0 0 NaN 0.0019 Opto4: sessions 13; site-stimulation pairs tested 304; with a significant increase 0; distinct sites with one 0; median response -0.0002 Opto9: sessions 12; site-stimulation pairs tested 203; with a significant increase 0; distinct sites with one 0; median response 0.0003 Opto2: sessions 29; site-stimulation pairs tested 619; with a significant increase 142; distinct sites with one 20; median response 0.0015
Result: reproduced. In the tdTomato mice no included non-stimulated fiber shows a significant increase after stimulation: 0 of 304 fiber and session pairs in Opto4 (13 sessions) and 0 of 203 in Opto9 (12 sessions), with median responses of -0.0002 and 0.0003 DF/F; one Opto9 pair has a significant decrease. The same test finds increases in 142 of 619 pairs (20 distinct fibers) in Opto2. The files hold 25 tdTomato sessions where the paper reports 26 stimulations. The stimulated fiber is taken as the one with the largest laser artefact, which can be wrong when the artefact spreads; the true stimulated fiber is then among those tested, which makes the test stricter.
figs7a-photoswitch-confined¶
In mice expressing rDAm3.0 without opsin, light through one fiber raises fluorescence at the stimulated site by photoswitching, but at no other site (0 of 22 non-stimulated sites for 3 stimulations).
for o in [o for o in OPTO if o["subject"] == "UG21"]:
inc_o = o["table"].included.values.copy()
oth = inc_o.copy()
oth[o["stim"]] = False
r = o["resp"]
print(f"UG21 {o['session']}: stimulated fiber {o['stim'] + 1} (included {bool(inc_o[o['stim']])}), epochs {len(o['off'])};"
f" stimulated fiber: response {r['resp'][o['stim']]:.4f}, significant increase {bool(r['sig_up'][o['stim']])};"
f" other included fibers {int(oth.sum())}: significant increase {int((r['sig_up'] & oth).sum())},"
f" any {int((r['sig'] & oth).sum())}, median response {np.median(r['resp'][oth]):.4f}")
UG21 230321a: stimulated fiber 12 (included True), epochs 8; stimulated fiber: response 0.0144, significant increase True; other included fibers 21: significant increase 0, any 0, median response -0.0021 UG21 230321b: stimulated fiber 18 (included True), epochs 8; stimulated fiber: response 0.0118, significant increase True; other included fibers 21: significant increase 0, any 0, median response -0.0035 UG21 230321c: stimulated fiber 24 (included True), epochs 8; stimulated fiber: response 0.0138, significant increase True; other included fibers 21: significant increase 0, any 0, median response -0.0003
Result: reproduced. In each of UG21's 3 sessions the stimulated fiber (12, 18, 24) has a significant increase after the laser (0.0118 to 0.0144 DF/F), and none of the 21 other included fibers has a significant response of either sign (median responses -0.0035 to -0.0003). The paper also required these sites to respond to reward; UG21 has no reward events in the files (data check), so that inclusion step cannot be applied.
Claims not tested here¶
The data of these claims are not in the dandiset. Except for Figure S3, the paper deposited them on Zenodo (doi:10.5281/zenodo.10272874); this attempt did not test them there.
figs1-no-tissue-damage: NeuN and DAPI images of implanted and unimplanted striatum (SuppFig01.zip).figs3-targeting-errors: the dandiset's fiber table holds the localized tips but no planned targets (data check). The paper's code (SuppFig03/B/analysis_code/planned_v_actual.m) reads the planned targets from per-mouse calibration files (calibPath), which are in none of the Zenodo archives; SuppFig03.zip holds only the CT fiber tables (MATLAB table objects) and the code. Blocked asdata_missing.fig3-freely-moving-signals: miniscope recordings in freely moving mice (Fig03.zip).fig7d-ventral-not-dorsal-jawandfig7ef-single-ventral-fiber: the D1-cre mice OptoD1_6, 7, 8 and 10, their DeepLabCut jaw tracking and stimulation files (Fig07.zip, SuppFig07.zip). None of these mice is in the dandiset.
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome |
|---|---|
| fig1c-isosbestic-control | reproduced |
| figs2d-correlation-distance | reproduced |
| figs2ef-no-bundle-crosstalk | reproduced |
| fig1d-chronic-reward-stability | partially reproduced (weeks 5 to 8 only) |
| fig4-link-published-results (with inputs) | reproduced |
| fig4e-light-increases-widespread | reproduced |
| fig4e-light-anterior-posterior-gradient | partially reproduced |
| fig4f-tone-tail-increase-anterior-decrease | reproduced |
| fig4f-decreases-lateral-ventral | partially reproduced |
| figs4ac-patterns-consistent-across-mice | partially reproduced |
| figs4eh-patterns-not-movement | partially reproduced |
| figs5-isosbestic-no-pattern | partially reproduced |
| fig5ab-response-latencies | partially reproduced |
| fig5g-light-latency-gradient | reproduced |
| fig5h-tone-latency-gradient | reproduced |
| figs6cd-gradients-individual-mice | partially reproduced |
| fig6-inclusion-counts | partially reproduced |
| fig6gh-release-decays-with-distance | reproduced |
| fig6g-dorsal-spread | partially reproduced |
| fig6-tdtomato-control-no-response | reproduced |
| figs7a-photoswitch-confined | reproduced |
| figs3 | blocked: data_missing (planned targets not public) |
| figs1, fig3, fig7d, fig7ef | blocked: data only on Zenodo, not attempted |
How the analysis differs from the paper's. The GLM, its null and the latency rule follow the
Methods and the paper's code, re-implemented in Python: the low-pass filter is a Butterworth filter in
place of MATLAB's lowpass; the locomotion lags come from a cross-correlation that shifts negative lags,
where the code does not (the code's version is shown, r = 0.997 for light and 0.989 for tone with the main values); latencies are
counted from the onset frame as in the Methods, where the code counts one frame later. DL18's tones
include all 44 presentations, as the files do not give their frequency (its 14 12-kHz tones are used as a
variant). Stimulation times and stimulated fibers in the optogenetic sessions are found from the DF/F
alone, without the trigger line. The striatal tail boundary (-1.0 mm) and several criteria (thresholds
for "many", "low", "comparable") are choices made here, before the analysis. The pooled spatial
correlations are also shown ranked within mice; the figs4eh splits with a pre-onset baseline; the figs5
GLM with 13 spline functions. The test's chance rate per fiber (0.024 for an increase) is far from the
paper's stated p < 1e-6, which assumes independent frames.
What was found about the dandiset. It holds the data of Figures 1C-D, 4, 5, 6, S2D-F, S4, S5E-H,
S6 and S7A, and for the 7 salient-stimulus mice the traces, fiber order and stimulus times are exactly those of
the paper's analysis. This match was not checked for DL15, 842, DL31, DL32 or the optogenetic mice.
Missing are: the stimulation times and the identity of the stimulated fiber in every optogenetic session;
the frequency of DL18's tones; the planned fiber targets of Figure S3 (not public); the task variables that separate unpredicted rewards in DL15's later
sessions; one tdTomato session; the data of Figures 2, 3, 7, S1 and S5A-D (on Zenodo). The included
column is described as "successfully implanted" but marks fibers within the striatum, leaving out
border fibers the optogenetic analysis counts; UG21's indicator is given as tdTomato where the paper's
counts and code identify it as the rDAm3.0 mouse without opsin; the archive's summary gives the approach
as cell population imaging; and the DL31 and DL32 series store a regular frame rate whose times drift by
up to half a second from the frame timestamps of the behaviour.
print(f"Notebook ran in {(time.time() - started) / 60:.1f} min")
Notebook ran in 15.5 min