Reproducing Li, Namburi, Olson et al. (2022) from dandiset 000689¶
Paper. Li H, Namburi P, Olson JM, et al. (2022). Neurotensin orchestrates valence assignment in the amygdala. Nature 608, 586-592. doi:10.1038/s41586-022-04964-y. The claims file was written from the accepted author manuscript in PubMed Central (PMC9583860): main text, Methods with the per-panel statistics, and the legends of Figures 1 to 5 and Extended Data Figures 1 to 10.
Dandiset. DANDI:000689, version 0.240530.1923. The paper cites no DANDI version (it predates the dandiset, created in October 2023). There are two published versions, 0.240524.2012 and 0.240530.1923; the data check below shows that both list the same 85 assets with the same asset identifiers, so they do not differ in content. The most recent, 0.240530.1923, is pinned, as the protocol asks when the paper cites none.
This notebook tests the electrophysiology, behavior and NT sensor photometry claims that the files can support, and shows why the others are blocked.
What the dandiset holds. 31 sessions of BLA recordings during the three-cue discrimination task (one per mouse, spike times of 1902 units, trial events, DeepLabCut pose in 29 of them), two pose-only files for other sessions of two of those mice, 32 videos, and 20 photometry sessions of the GRAB-NTS1.0 NT sensor in the BLA (5 mice, early and late acquisition of tone-sucrose and tone-airpuff learning). No file says which mice are CRISPR and which are controls, nor which projection (NAc or CeA) was tagged in each mouse, the units are not the paper's sorted single units, and the pose has no timestamps, so its clock is fitted to the trial table (data check).
| Claim | Tested here | Outcome basis |
|---|---|---|
| fig3k-nt-sensor-sucrose-increase | yes | NT sensor files |
| fig3l-nt-sensor-airpuff-decrease | yes | NT sensor files |
| fig3-nt-sensor-early-dip | yes | NT sensor files |
| fig4-recording-count | yes | units table |
| fig4c-photoidentification | run, blocked, metadata_insufficient | light events and units; no fake-block or projection labels |
| fig4d-seven-clusters | run, blocked, methods_underspecified | units table |
| fig5dh-behavior-clusters | yes | pose, on a clock fitted to the port entries |
| fig2c, ed8d, ed8fg, fig4e, ed8i, fig4fg (2), fig5a, fig5b, fig5c, ed8e, fig5fi, ed8os | blocked, metadata_insufficient | no group or projection labels |
| fig1c, ed1c, ed2, ed3, ed5f, fig1g, fig1h, fig2d, fig2g, fig2h, fig3bc, fig3e, fig3f, ed6, fig3h, ed9, ed9kl, ed10b, ed10 | blocked, data_missing | experiments not in the dandiset |
import time, json, warnings
from concurrent.futures import ThreadPoolExecutor
import numpy as np
import pandas as pd
import requests
import matplotlib.pyplot as plt
from scipy import stats
from scipy.cluster.hierarchy import linkage, fcluster
from scipy.spatial.distance import pdist
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
import nt_bla as nb
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 220); pd.set_option("display.max_columns", 40); pd.set_option("display.max_rows", 200)
started = time.time()
A = nb.assets()
print(f"dandiset {nb.DANDISET} version {nb.VERSION}: {len(A)} assets")
dandiset 000689 version 0.240530.1923: 85 assets
Data check¶
Each check of the metadata against the values is stated under its output.
# The two published versions: same assets?
other = []
url = f"{nb.API}/dandisets/{nb.DANDISET}/versions/0.240524.2012/assets/?page_size=200"
while url:
r = requests.get(url, timeout=60).json(); other += r["results"]; url = r.get("next")
a1 = {a["path"]: (a["asset_id"], a["size"]) for a in A}
a2 = {a["path"]: (a["asset_id"], a["size"]) for a in other}
print("assets in 0.240530.1923:", len(a1), " in 0.240524.2012:", len(a2),
" same paths:", set(a1) == set(a2), " same asset ids and sizes:", a1 == a2)
kind = lambda p: ("photometry nwb" if "ophys" in p else "ecephys nwb" if "ecephys" in p and p.endswith(".nwb")
else "pose-only nwb" if p.endswith(".nwb") else "video")
inv = pd.DataFrame({"path": list(a1), "GB": [a1[p][1] / 1e9 for p in a1]})
inv["kind"] = inv.path.map(kind)
inv["subject"] = inv.path.str.split("/").str[0]
print(inv.groupby("kind").agg(files=("path", "size"), subjects=("subject", "nunique"), GB=("GB", "sum")).round(1).to_string())
assets in 0.240530.1923: 85 in 0.240524.2012: 85 same paths: True same asset ids and sizes: True
files subjects GB
kind
ecephys nwb 31 31 534.0
photometry nwb 20 5 0.1
pose-only nwb 2 2 0.6
video 32 30 40.2
The two published versions hold the same assets, so the choice between them cannot change a result. The archive's summary counts 36 subjects: 31 electrophysiology mice and 5 photometry mice.
# Electrophysiology: spike times and trial tables of the 31 sessions (raw voltage is not read).
EP = sorted(a["path"] for a in A if a["path"].endswith(".nwb") and "ecephys" in a["path"])
t = time.time()
with ThreadPoolExecutor(4) as ex:
E = list(ex.map(nb.load_ephys, EP))
print(f"loaded {len(E)} sessions in {time.time() - t:.0f} s")
rows = []
for e in E:
tr = e["trials"]
r = {"subject": e["subject"], "session": e["session_id"][-22:], "genotype": e["genotype"],
"units": len(e["spikes"]), "channels": e["n_channels"], "dur_s": round(e["duration"])}
for k, v in nb.CS_TYPES.items():
r[k] = len(tr[v][0]) if v in tr else 0
r["sucrose_deliv"] = len(tr.get("reward_delivery", ([],))[0])
r["shock_relay"] = len(tr.get("shock_relay", ([],))[0])
r["light"] = len(tr.get("phototagging", ([],))[0])
cs = np.concatenate([tr[v][0] for v in nb.CS_TYPES.values() if v in tr])
r["task_end_s"] = round(cs.max()); r["light_start_s"] = round(tr["phototagging"][0].min())
rows.append(r)
ses = pd.DataFrame(rows)
print(ses.to_string(index=False))
print("\nmice:", ses.subject.nunique(), " sessions:", len(ses), " units:", ses.units.sum(),
" genotype values:", sorted(ses.genotype.unique()))
print("sessions without neutral CS trials:", list(ses.subject[ses.neutral == 0]))
print("CS counts (sucrose, neutral, shock) in the other sessions:",
ses[ses.neutral > 0][["sucrose", "neutral", "shock"]].apply(tuple, axis=1).value_counts().to_dict())
loaded 31 sessions in 111 s
subject session genotype units channels dur_s sucrose neutral shock sucrose_deliv shock_relay light task_end_s light_start_s
H01 9-07-13-11-59-42-DiscS Wild type 54 32 11888 60 30 30 59 30 854 10916 11130
H02 9-07-14-11-14-47-DiscS Wild type 64 32 12091 60 30 30 60 30 854 10934 11188
H07 9-07-13-18-44-41-DiscS Wild type 64 32 12789 60 30 30 60 30 856 10934 12008
H08 9-07-14-18-11-14-DiscS Wild type 60 32 12048 60 30 30 60 30 854 10951 11221
H09 9-08-30-16-32-13-Disc4 Wild type 62 32 12224 75 0 45 72 45 854 10986 11089
H10 9-08-30-16-32-20-Disc4 Wild type 64 32 12221 75 0 45 75 45 854 10980 11089
H11 9-08-30-12-58-04-Disc4 Wild type 64 32 12271 75 0 45 75 45 854 10974 11498
H12 9-08-30-20-05-12-Disc4 Wild type 64 32 11847 75 0 45 73 45 854 10982 11088
H13 9-08-30-09-21-02-Disc4 Wild type 64 32 12261 75 0 45 63 45 854 10969 11391
H14 9-08-30-09-21-05-Disc4 Wild type 64 32 12261 73 0 45 70 45 854 10966 11392
H17 1-14-10-52-Disc5-20khz Wild type 64 32 12118 60 30 30 60 30 854 10989 11015
H19 -17-44-38-Disc5-3,5khz Wild type 64 32 12434 60 30 30 60 30 854 11001 11154
H22 9-12-27-13-58-44-Disc3 Wild type 64 32 12026 60 30 30 60 30 854 10978 11252
H23 9-12-29-09-55-50-Disc4 Wild type 64 32 12287 60 30 30 27 31 854 10990 11251
H24 9-12-27-10-25-11-Disc3 Wild type 46 32 12451 60 30 30 60 30 859 10986 11538
H25 20-10-06-39-Disc4-3.5k Wild type 58 32 12307 60 30 30 56 31 854 11043 11358
H27 -20-10-07-01-Disc4-20k Wild type 38 32 12300 60 30 30 60 31 854 11034 11353
H28 -19-14-27-39-Disc3-20k Wild type 64 32 12114 60 30 30 60 30 854 11025 11238
H29 -20-13-42-59-Disc4-20k Wild type 64 32 12831 60 30 30 60 30 854 11009 11435
H31 0-04-30-09-56-31-Disc7 Wild type 64 32 12186 60 30 30 60 30 852 11032 3
H32 0-04-23-13-59-47-Disc5 Wild type 64 32 11643 60 30 30 60 30 852 11035 3
H33 0-04-15-12-44-39-Disc5 Wild type 64 32 11638 60 30 30 60 30 854 10988 3
H34 0-04-15-12-44-32-Disc5 Wild type 64 32 11690 60 30 28 60 30 854 10995 13
H35 09-57-13-Disc7-8,22,23 Wild type 58 32 12152 60 30 30 60 30 854 10991 2
H36 0-04-24-12-47-33-Disc6 Wild type 64 32 11617 60 30 30 60 30 854 11001 2
H37 29-10-04-51-Disc5-9-11 Wild type 58 32 11590 60 30 30 58 30 854 10994 2
H38 0-04-29-10-04-41-Disc5 Wild type 64 32 11610 60 30 30 60 30 854 11006 2
H39 0-04-29-13-33-40-Disc5 Wild type 64 32 11535 60 30 30 42 30 854 11012 2
H40 0-04-29-13-33-54-Disc5 Wild type 64 32 11525 60 30 30 60 30 854 11000 2
H41 04-29-16-55-26-Disc5-9 Wild type 62 32 11532 59 30 30 60 30 854 10992 2
H42 04-29-16-55-04-Disc5-5 Wild type 62 32 11558 60 30 30 60 30 854 11016 29
mice: 31 sessions: 31 units: 1902 genotype values: ['Wild type']
sessions without neutral CS trials: ['H09', 'H10', 'H11', 'H12', 'H13', 'H14']
CS counts (sucrose, neutral, shock) in the other sessions: {(60, 30, 30): 23, (60, 30, 28): 1, (59, 30, 30): 1}
Checked against the paper. The session description says sucrose, neutral and shock CSs at 2:1:1, and 25 sessions hold 60, 30 and 30 CSs (some with one or two fewer), as the Methods give (120 trials). Six sessions (H09 to H14, recorded on 2019-08-30, "Disc4") hold 75 sucrose and 45 shock CSs and no neutral CS: the description does not fit them. Every file gives the genotype as "Wild type"; none says whether the mouse is in the CRISPR or the control group, nor whether Cav2-Cre was put in the NAc or the CeA, which the group comparisons and the phototagging need. There are 31 mice where the paper reports 30. H28 and H32 also have a pose-only file for a later discrimination day (Disc4, Disc6) than their recording (Disc3, Disc5), so their recorded session may not be the last day, which is the only day the paper analyses. Each session holds about 850 light events (the phototagging sessions); no "fake" session is labelled.
# Are the units sorted single units? Refractory violations and spikes shared between units.
def coinc(a, b, w=1e-4):
if len(a) == 0 or len(b) < 2: return 0.0
j = np.clip(np.searchsorted(b, a), 1, len(b) - 1)
return float(np.mean(np.minimum(np.abs(b[j] - a), np.abs(b[j - 1] - a)) < w))
def unit_quality(e):
sp = e["spikes"]
isi = [np.mean(np.diff(s) < 0.001) if len(s) > 1 else np.nan for s in sp]
allsp = np.concatenate(sp)
shared = 1 - len(np.unique(np.round(allsp * 30000).astype(np.int64))) / len(allsp)
nn = min(len(sp), 16)
pair = [coinc(sp[i], sp[j]) for i in range(nn) for j in range(nn) if i != j]
return {"subject": e["subject"], "units": len(sp), "median_frac_isi_lt_1ms": np.nanmedian(isi),
"units_with_isi_viol_gt_5pct": int(np.sum(np.array(isi) > 0.05)),
"frac_spike_samples_shared": shared, "median_pair_coincidence_0.1ms": np.median(pair)}
with ThreadPoolExecutor(4) as ex:
q = pd.DataFrame(list(ex.map(unit_quality, E)))
print(q.round(3).to_string(index=False))
print("\nover all units: median fraction of ISIs below 1 ms = %.3f; units with more than 5%% = %d of %d"
% (q.median_frac_isi_lt_1ms.median(), q.units_with_isi_viol_gt_5pct.sum(), q.units.sum()))
print("units per session: %s" % q.units.value_counts().sort_index().to_dict())
subject units median_frac_isi_lt_1ms units_with_isi_viol_gt_5pct frac_spike_samples_shared median_pair_coincidence_0.1ms
H01 54 0.265 54 0.201 0.078
H02 64 0.190 61 0.101 0.028
H07 64 0.168 63 0.136 0.045
H08 60 0.217 59 0.168 0.057
H09 62 0.241 61 0.202 0.074
H10 64 0.351 64 0.267 0.158
H11 64 0.255 62 0.178 0.056
H12 64 0.289 62 0.194 0.093
H13 64 0.315 62 0.216 0.090
H14 64 0.243 64 0.408 0.281
H17 64 0.048 12 0.652 0.647
H19 64 0.046 23 0.166 0.049
H22 64 0.048 31 0.102 0.021
H23 64 0.045 26 0.102 0.031
H24 46 0.031 10 0.191 0.048
H25 58 0.030 8 0.257 0.061
H27 38 0.046 14 0.449 0.056
H28 64 0.061 39 0.076 0.109
H29 64 0.039 25 0.141 0.084
H31 64 0.046 27 0.130 0.110
H32 64 0.032 17 0.114 0.043
H33 64 0.037 21 0.238 0.089
H34 64 0.074 49 0.133 0.105
H35 58 0.058 32 0.137 0.098
H36 64 0.026 22 0.105 0.064
H37 58 0.025 20 0.057 0.016
H38 64 0.033 18 0.116 0.032
H39 64 0.036 15 0.267 0.075
H40 64 0.023 12 0.081 0.031
H41 62 0.032 9 0.479 0.070
H42 62 0.028 19 0.081 0.028
over all units: median fraction of ISIs below 1 ms = 0.046; units with more than 5% = 1061 of 1902
units per session: {38: 1, 46: 1, 54: 1, 58: 3, 60: 1, 62: 3, 64: 21}
Checked against the paper. The Methods say single units were sorted in Offline Sorter, and the paper analyses 685 neurons from 30 mice (about 23 per mouse). The files hold 1902 units, 64 in 21 of 31 sessions (two per channel of the 32-channel optrode). In 1061 units more than 5% of the inter-spike intervals are below 1 ms (a median of 17 to 35% per session in H01 to H14), which a well-isolated single neuron does not show, and in each session 6 to 65% of spike samples recur in another unit. These are not the paper's sorted single units: they look like unsorted or multi-unit threshold crossings. The claims that use units are tested on them as they are, and this is recorded as a finding.
# Every text field of each recording file, the units and electrodes columns, and the archive's
# per-asset metadata: is there any group (CRISPR, control) or projection (NAc, CeA) label?
# Also every object under acquisition/ and processing/, with its type, and for image containers
# the description and the shape of each image.
import re
def contents(h, path):
objs, images = [], []
for top in ("acquisition", "processing"):
for k in h[top]:
o = h[top][k]
nt = nb._s(o.attrs.get("neurodata_type", ""))
objs.append(f"{top}/{k.split(':')[0]} ({nt})")
if nt == "ProcessingModule":
objs += [f"{top}/{k}/{j} ({nb._s(o[j].attrs.get('neurodata_type', ''))})" for j in o]
if nt == "Images":
images.append({"subject": path.split("/")[0][4:], "container": k,
"description": nb._s(o.attrs.get("description", ""))[:110],
"n_images": sum(hasattr(o[j], "shape") for j in o),
"first images": "; ".join(f"{j} {tuple(o[j].shape)}" for j in list(o)[:2] if hasattr(o[j], "shape"))})
return objs, images
def text_fields(path):
h = nb.open_h5(path)
out = {}
def visit(name, obj):
if isinstance(obj, nb.h5py.Dataset) and obj.dtype.kind in "OSU" and obj.size <= 64 and not name.startswith(("units/", "intervals/", "acquisition/", "processing/")):
v = obj[()]
out[name] = " | ".join(str(nb._s(x)) for x in np.atleast_1d(v) if isinstance(nb._s(x), str))
h["general"].visititems(lambda n, o: visit("general/" + n, o))
out["session_description"] = nb._s(h["session_description"][()])
out["units columns"] = ", ".join(sorted(h["units"].keys()))
out["electrodes columns"] = ", ".join(sorted(h["general/extracellular_ephys/electrodes"].keys()))
out["electrode locations"] = ", ".join(sorted(set(nb._s(x) for x in h["general/extracellular_ephys/electrodes/location"][:])))
return out, contents(h, path)
with ThreadPoolExecutor(4) as ex:
TF, CT = zip(*ex.map(text_fields, EP))
fields = sorted(set().union(*TF))
for f in fields:
vals = sorted(set(t.get(f, "(absent)") for t in TF))
shown = vals if f != "general/experiment_description" else [v[:80] + "..." for v in vals]
print(f"{f}: {len(vals)} distinct value(s): {shown[:6]}{' ...' if len(vals) > 6 else ''}")
pat = re.compile(r"crispr|control|scrambl|sgrna|guide|nac|accumbens|cea|central amyg|cav2|projection|group", re.I)
hits = sorted({(f, m.group(0)) for t in TF for f, v in t.items() if f != "general/experiment_description" for m in pat.finditer(v)})
print("\nlabel words in the files' fields (experiment_description, the paper's abstract, left out):", hits)
def asset_subject(a):
r = requests.get(f"{nb.API}/dandisets/{nb.DANDISET}/versions/{nb.VERSION}/assets/{a['asset_id']}/", timeout=60).json()
return {k: str(v) for k, v in r.get("wasAttributedTo", [{}])[0].items() if k not in ("identifier",)}
with ThreadPoolExecutor(8) as ex:
am = pd.DataFrame(list(ex.map(asset_subject, [a for a in A if a["path"] in EP])))
print("\narchive asset metadata, subject fields:", list(am.columns))
print({c: sorted(am[c].unique().tolist())[:4] for c in am.columns if c not in ("age",)})
counts = pd.Series([x for objs, _ in CT for x in objs]).value_counts()
print("\nobjects under acquisition/ and processing/ (number of the 31 files holding each):")
print(counts.to_string())
imgs = pd.DataFrame([r for _, im in CT for r in im])
print(f"\nimage containers: {len(imgs)} in {imgs.subject.nunique()} files; files without: {[e['subject'] for e, (_, im) in zip(E, CT) if not im]}")
print(imgs.to_string(index=False, max_colwidth=110))
print("image names in all containers:", sorted({re.sub(r"\d+", "#", n.split(" ")[0]) for s in imgs["first images"] for n in s.split("; ")}))
electrode locations: 1 distinct value(s): ['unknown']
electrodes columns: 1 distinct value(s): ['channel_name, gain_to_uV, group, group_name, id, location, offset_to_uV']
general/experiment_description: 1 distinct value(s): ['The ability to associate temporally segregated information and assign positive o...']
general/experimenter: 1 distinct value(s): ['Li, Hao']
general/extracellular_ephys/electrodes/channel_name: 1 distinct value(s): ['CH1 | CH2 | CH3 | CH4 | CH5 | CH6 | CH7 | CH8 | CH9 | CH10 | CH11 | CH12 | CH13 | CH14 | CH15 | CH16 | CH17 | CH18 | CH19 | CH20 | CH21 | CH22 | CH23 | CH24 | CH25 | CH26 | CH27 | CH28 | CH29 | CH30 | CH31 | CH32']
general/extracellular_ephys/electrodes/group: 1 distinct value(s): ['']
general/extracellular_ephys/electrodes/group_name: 1 distinct value(s): ['ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup | ElectrodeGroup']
general/extracellular_ephys/electrodes/location: 1 distinct value(s): ['unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown | unknown']
general/institution: 1 distinct value(s): ['Salk Institute for Biological Studies']
general/lab: 1 distinct value(s): ['Tye']
general/related_publications: 2 distinct value(s): ['doi: 10.1038/s41586-022-04964-y', 'https://doi.org/10.1038/s41586-022-04964-y']
general/session_id: 31 distinct value(s): ['1-2019-07-13-11-59-42-DiscS', '10-2019-08-30-16-32-20-Disc4', '11-2019-08-30-12-58-04-Disc4', '12-2019-08-30-20-05-12-Disc4', '13-2019-08-30-09-21-02-Disc4', '14-2019-08-30-09-21-05-Disc4'] ...
general/subject/age: 2 distinct value(s): ['P109D', 'P60D']
general/subject/genotype: 1 distinct value(s): ['Wild type']
general/subject/sex: 1 distinct value(s): ['M']
general/subject/species: 1 distinct value(s): ['Mus musculus']
general/subject/strain: 1 distinct value(s): ['C57BL/6J']
general/subject/subject_id: 31 distinct value(s): ['H01', 'H02', 'H07', 'H08', 'H09', 'H10'] ...
session_description: 2 distinct value(s): ['Auto-generated by neuroconv', 'Discrimination task in which three types of trials with distinct auditory tones predicting sucrose, shock or no outcome (sucrose, shock and neutral conditioned stimulus (CS), respectively) were presented randomly at a 2:1:1 ratio']
units columns: 1 distinct value(s): ['id, spike_times, spike_times_index, unit_name']
label words in the files' fields (experiment_description, the paper's abstract, left out): [('electrodes columns', 'group'), ('general/extracellular_ephys/electrodes/group_name', 'Group')]
archive asset metadata, subject fields: ['age', 'sex', 'species', 'genotype', 'schemaKey']
{'sex': ["{'name': 'Male', 'schemaKey': 'SexType', 'identifier': 'http://purl.obolibrary.org/obo/PATO_0000384'}"], 'species': ["{'name': 'Mus musculus - House mouse', 'schemaKey': 'SpeciesType', 'identifier': 'http://purl.obolibrary.org/obo/NCBITaxon_10090'}"], 'genotype': ['Wild type'], 'schemaKey': ['Participant']}
objects under acquisition/ and processing/ (number of the 31 files holding each):
acquisition/ElectricalSeries (ElectricalSeries) 31
acquisition/Video (ImageSeries) 30
processing/behavior (ProcessingModule) 29
processing/behavior/PoseEstimation (PoseEstimation) 29
acquisition/COFAImages (Images) 27
image containers: 27 in 27 files; files without: ['H01', 'H02', 'H07', 'H08']
subject container description n_images first images
H09 COFAImages The confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H9_20x_PVT.oif. 8 GrayScaleImage1Depth1 (2420, 1992); GrayScaleImage1Depth2 (2420, 1992)
H10 COFAImages The 0 confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H10_20x_PVT.oif. 8 GrayScaleImage1Depth1 (2418, 2440); GrayScaleImage1Depth2 (2418, 2440)
H11 COFAImages The 1 confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H11_20x_PVT.oif. 8 GrayScaleImage1Depth1 (2866, 1988); GrayScaleImage1Depth2 (2866, 1988)
H12 COFAImages The 2 confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H12_slidePVT3_slice3_zstack_40x_PVT. 12 GrayScaleImage1Composite (3764, 2448); GrayScaleImage1Depth1 (3764, 2448)
H13 COFAImages The 3 confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H13_slidePVT4_slice6_zstack_40x_PVT. 18 GrayScaleImage1Composite (2418, 2900); GrayScaleImage1Depth1 (2418, 2900)
H14 COFAImages The 4 confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H14_slidePVT4_slice6_zstack_40x_PVT. 12 GrayScaleImage1Depth1 (3764, 2896); GrayScaleImage1Depth2 (3764, 2896)
H17 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H17PVT_40x.oif. 42 GrayScaleImage1Composite (3328, 3786); GrayScaleImage1Depth1 (3328, 3786)
H19 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H19PVT_40x.oif. 27 GrayScaleImage1Composite (3772, 3348); GrayScaleImage1Depth1 (3772, 3348)
H22 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H22PVT_40x.oif. 15 GrayScaleImage1Composite (3772, 3338); GrayScaleImage1Depth1 (3772, 3338)
H23 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H23PVT_40x.oif. 18 GrayScaleImage1Composite (3764, 3338); GrayScaleImage1Depth1 (3764, 3338)
H24 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H24PVT_40x.oif. 15 GrayScaleImage1Composite (3766, 3798); GrayScaleImage1Depth1 (3766, 3798)
H25 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H25PVT_40x.oif. 27 GrayScaleImage1Composite (3772, 3784); GrayScaleImage1Depth1 (3772, 3784)
H27 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H27PVT_40x.oif. 32 GrayScaleImage1Composite (4660, 4254); GrayScaleImage1Depth1 (4660, 4254)
H28 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H28PVT_40x.oif. 28 GrayScaleImage1Composite (3774, 3782); GrayScaleImage1Depth1 (3774, 3782)
H29 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H29PVT_40x.oif. 32 GrayScaleImage1Composite (3766, 3338); GrayScaleImage1Depth1 (3766, 3338)
H31 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H31PVT_40x.oif. 28 GrayScaleImage1Composite (3316, 3340); GrayScaleImage1Depth1 (3316, 3340)
H32 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H32PVT_40x.oif. 32 GrayScaleImage1Composite (3770, 2888); GrayScaleImage1Depth1 (3770, 2888)
H33 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H33PVT_40x.oif. 40 GrayScaleImage1Composite (3778, 2894); GrayScaleImage1Depth1 (3778, 2894)
H34 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H34PVT_40x.oif. 24 GrayScaleImage1Composite (3320, 3340); GrayScaleImage1Depth1 (3320, 3340)
H35 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H35PVT_40x.oif. 32 GrayScaleImage1Composite (3318, 2890); GrayScaleImage1Depth1 (3318, 2890)
H36 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H36PVT_40x.oif. 24 GrayScaleImage1Composite (3318, 2444); GrayScaleImage1Depth1 (3318, 2444)
H37 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H37PVT_40x.oif. 32 GrayScaleImage1Composite (4224, 4224); GrayScaleImage1Depth1 (4224, 4224)
H38 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H38PVT_40x.oif. 24 GrayScaleImage1Composite (3766, 3786); GrayScaleImage1Depth1 (3766, 3786)
H39 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H39PVT_40x.oif. 24 GrayScaleImage1Composite (3766, 3790); GrayScaleImage1Depth1 (3766, 3790)
H40 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H40PVT_40x.oif. 24 GrayScaleImage1Composite (3766, 3782); GrayScaleImage1Depth1 (3766, 3782)
H41 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H41PVT_40x.oif. 28 GrayScaleImage1Composite (2880, 3330); GrayScaleImage1Depth1 (2880, 3330)
H42 COFAImages The PVT confocal images from FLUOVIEW FV1000 version 1.2.6.0 extracted from H42PVT_40x.oif. 24 GrayScaleImage1Composite (3320, 3788); GrayScaleImage1Depth1 (3320, 3788)
image names in all containers: ['GrayScaleImage#Composite', 'GrayScaleImage#Depth#']
No field of the files, and no subject field of the archive's metadata, names a CRISPR or control group, a guide RNA, or the NAc or CeA (the words found, if any, are listed above with their field). The units table holds only spike times and a unit name; the electrodes table gives no location.
27 of the 31 files (all but H01, H02, H07 and H08) also hold a stack of confocal images (COFAImages, 8 to 42 images each), whose description names an Olympus FLUOVIEW file of the PVT (for example H17PVT_40x.oif). The images are named only by depth (GrayScaleImage1Depth#, and a composite), with no channel, probe or group label.
# Clock check for the spikes and the trial table: CS-evoked change in population rate at the
# true CS onsets against the same onsets shifted 10 s earlier (inside the inter-trial interval).
def pop_resp(e, shift):
on = np.concatenate([e["trials"][v][0] for v in nb.CS_TYPES.values() if v in e["trials"]]) + shift
out = []
for sp in e["spikes"]:
p = nb.psth(sp, on, 1.0, 1.0, 0.1)
b = p[:10].mean()
out.append((p[10:13].mean() - b) / (p[:10].std(ddof=1) + 1e-9))
return np.median(np.abs(out))
ck = pd.DataFrame([{"subject": e["subject"], "true_onsets": pop_resp(e, 0.0), "shifted_-10s": pop_resp(e, -10.0)} for e in E])
print("median |z| of the 0 to 0.3 s response across units, per session:")
print(ck.round(2).to_string(index=False))
print("sessions where the true onsets give the larger response: %d of %d" % ((ck.true_onsets > ck["shifted_-10s"]).sum(), len(ck)))
median |z| of the 0 to 0.3 s response across units, per session:
subject true_onsets shifted_-10s
H01 8.25 0.41
H02 2.92 0.50
H07 1.94 0.50
H08 1.16 0.49
H09 3.75 1.39
H10 7.76 0.58
H11 1.99 0.93
H12 5.99 0.44
H13 2.12 0.55
H14 1.01 0.50
H17 0.00 0.00
H19 4.70 0.42
H22 0.50 0.65
H23 1.66 0.67
H24 1.05 0.48
H25 0.71 0.60
H27 0.98 1.02
H28 0.63 0.58
H29 0.96 0.60
H31 2.17 0.34
H32 3.53 0.47
H33 1.89 0.51
H34 0.41 0.46
H35 4.36 0.63
H36 2.28 0.51
H37 1.18 0.52
H38 1.86 0.42
H39 0.48 0.51
H40 2.08 0.42
H41 1.16 0.70
H42 2.16 0.99
sessions where the true onsets give the larger response: 26 of 31
The units respond more at the true CS onsets than at onsets shifted into the inter-trial interval in the sessions counted above, which is consistent with the spikes and the trial table sharing a clock (a check that supports, without proving, the alignment). It fails in H17, whose spike times end at 249 s of a 12118 s recording (next cell): its units cannot be aligned with its trials. In H22, H27, H34 and H39 the two values are close, so the alignment is not shown there.
# Span of the spike times against the length of the recording, per session.
span = pd.DataFrame([{"subject": e["subject"], "first_spike_s": round(min(s[0] for s in e["spikes"] if len(s)), 1),
"last_spike_s": round(max(s[-1] for s in e["spikes"] if len(s)), 1), "recording_s": round(e["duration"], 1),
"task_start_s": round(min(e["trials"][v][0].min() for v in nb.CS_TYPES.values() if v in e["trials"]), 1),
"task_end_s": round(max(e["trials"][v][0].max() for v in nb.CS_TYPES.values() if v in e["trials"]), 1),
"light_before_task": int(np.sum(e["trials"]["phototagging"][0] < min(e["trials"][v][0].min() for v in nb.CS_TYPES.values() if v in e["trials"]))),
"light_after_task": int(np.sum(e["trials"]["phototagging"][0] > max(e["trials"][v][0].max() for v in nb.CS_TYPES.values() if v in e["trials"])))} for e in E])
print(span.to_string(index=False))
print("sessions whose spikes stop before 90% of the recording:", list(span.subject[span.last_spike_s < 0.9 * span.recording_s]))
print("sessions with light events before the task:", list(span.subject[span.light_before_task > 0]))
# Light events form two blocks per session: split at the largest gap between events.
for e in E[:3] + E[-3:]:
x = e["trials"]["phototagging"][0]; g = np.argmax(np.diff(x))
print(f'{e["subject"]}: blocks of {g + 1} and {len(x) - g - 1} light events, gap {np.diff(x)[g]:.0f} s')
subject first_spike_s last_spike_s recording_s task_start_s task_end_s light_before_task light_after_task
H01 0.0 11887.6 11887.8 216.4 10916.4 0 854
H02 0.0 12090.3 12090.6 233.8 10933.8 0 854
H07 0.0 12788.0 12788.5 233.8 10933.8 0 856
H08 0.0 12047.7 12048.2 251.1 10951.1 0 854
H09 0.0 12220.5 12223.7 285.6 10985.7 0 854
H10 0.0 12213.5 12220.9 279.5 10979.5 0 854
H11 0.0 12270.6 12270.7 273.5 10973.5 0 854
H12 0.0 11846.3 11846.5 282.4 10982.4 0 854
H13 0.0 12249.0 12261.5 269.0 10969.0 0 854
H14 0.0 12259.5 12261.4 266.1 10966.1 0 854
H17 0.0 249.2 12117.8 289.2 10989.3 0 854
H19 0.0 11969.4 12434.4 301.0 11001.1 0 854
H22 0.0 12026.1 12026.2 278.2 10978.3 0 854
H23 0.0 12285.7 12286.8 290.1 10990.2 0 854
H24 0.0 12451.0 12451.5 285.8 10985.8 0 859
H25 0.0 12255.5 12306.8 342.6 11042.6 0 854
H27 0.0 12249.9 12300.0 333.5 11033.5 0 854
H28 0.0 11993.3 12113.8 325.1 11025.1 0 854
H29 0.1 12468.7 12831.4 308.9 11009.0 0 854
H31 0.1 12165.5 12186.5 332.0 11032.0 386 426
H32 0.0 11594.4 11643.5 335.3 11035.3 390 426
H33 0.0 11624.3 11637.6 288.1 10988.1 353 426
H34 0.0 11682.6 11690.5 294.9 10994.9 349 426
H35 0.0 12126.4 12152.4 290.5 10990.5 356 426
H36 0.0 11552.1 11616.6 300.5 11000.6 366 426
H37 0.0 11585.3 11590.0 294.1 10994.2 360 426
H38 0.0 11599.5 11609.6 306.0 11006.0 368 426
H39 0.0 11431.7 11534.6 311.9 11012.0 368 426
H40 0.0 11421.8 11524.5 299.7 10999.8 366 426
H41 0.0 11510.6 11531.7 291.7 10991.7 357 426
H42 0.0 11536.6 11558.0 315.5 11015.5 354 426
sessions whose spikes stop before 90% of the recording: ['H17']
sessions with light events before the task: ['H31', 'H32', 'H33', 'H34', 'H35', 'H36', 'H37', 'H38', 'H39', 'H40', 'H41', 'H42']
H01: blocks of 428 and 426 light events, gap 48 s
H02: blocks of 428 and 426 light events, gap 69 s
H07: blocks of 430 and 426 light events, gap 60 s
H40: blocks of 428 and 426 light events, gap 10714 s
H41: blocks of 428 and 426 light events, gap 10803 s
H42: blocks of 428 and 426 light events, gap 10801 s
H17's spike times all fall within the first 249 s of a recording of 12118 s, while its trial events cover the whole session: its spike times are not on the trial clock (a file error). Every session holds about 854 light events in two blocks of about 427, which fits the Methods' two light sessions (a "fake" one with the patch cord disconnected, then the real one), but no label says which block is the fake session. In 12 sessions (H31 to H42) one block comes before the task and one after it; the Methods put both after the task.
# Does the pose share the trial clock? The pose series give only a start time (0 s) and a nominal
# 15 frames/s, no timestamps. At each port entry (beam break, trial table) the head should be at the
# port. (1) With the nominal clock, the fixed lag (pose time minus trial time) that minimises the
# median head-to-port distance at port entries, in each third of the session. (2) A linear map,
# pose time = T + offset + factor * (t - T), fitted to all port entries (nb.fit_pose_clock), and
# the best extra lag left in each third under it (0 when the map holds throughout).
def pose_clock(p):
h = nb.open_h5(p)
subj = p.split("/")[0][4:]
if "processing/behavior/PoseEstimation" not in h:
return subj, None
pose = nb.load_pose(h); r, t0 = pose["rate"], pose["t0"]
m = nb.pose_measures(pose)
dist = m["distance_to_port"]
tt = h["intervals/trials/trial_type"][:].astype(str); st = h["intervals/trials/start_time"][:]
pe = np.sort(st[tt == "port_entry"])
out = {"subject": subj, "pose_s": round(len(dist) / r), "port_entries": len(pe)}
lags = np.arange(-200, 200.01, 0.2)
for k, seg in enumerate(np.array_split(pe, 3)):
med = nb._median_at(dist, r, t0, seg[None, :] + lags[:, None])
out[f"lag{k + 1}_s"] = round(float(lags[np.nanargmin(med)]), 1)
c = nb.fit_pose_clock(dist, r, t0, pe)
out.update({"fit_offset_s": round(c["offset"], 1), "fit_factor": round(c["factor"], 4),
"resid1_s": c["resid_lag_thirds"][0], "resid2_s": c["resid_lag_thirds"][1], "resid3_s": c["resid_lag_thirds"][2],
"dist_fit": round(c["med_dist_fit"], 1), "dist_lag0": round(c["med_dist_lag0"], 1), "dist_all_frames": round(c["med_dist_all"], 1)})
return subj, (out, {"measures": m, "rate": r, "t0": t0, "clock": c})
t = time.time()
with ThreadPoolExecutor(4) as ex:
PC = dict(ex.map(pose_clock, EP))
pl = pd.DataFrame([v[0] for v in PC.values() if v is not None])
print(f"pose read and clocks fitted in {time.time() - t:.0f} s; sessions without pose: {[k for k, v in PC.items() if v is None]}")
print(pl.round(1).assign(fit_factor=pl.fit_factor.round(4)).to_string(index=False))
spread = pl[["lag1_s", "lag2_s", "lag3_s"]].max(axis=1) - pl[["lag1_s", "lag2_s", "lag3_s"]].min(axis=1)
resid = pl[["resid1_s", "resid2_s", "resid3_s"]].abs().max(axis=1)
pl["aligned"] = resid <= 1.0
print("\nnominal clock: median |lag| %.1f s; sessions whose best lag changes by more than 5 s within the session: %d of %d"
% (np.median(np.abs(pl[["lag1_s", "lag2_s", "lag3_s"]].values)), int((spread > 5).sum()), len(pl)))
print("linear map: sessions with every third within 1 s of the map: %d of %d; not: %s"
% (pl.aligned.sum(), len(pl), list(pl.subject[~pl.aligned])))
print("fitted factor in the aligned sessions: median %.4f, range %.4f to %.4f; median distance at port entries %.1f px under the map, %.1f px at lag 0, %.1f px over all frames"
% (pl.fit_factor[pl.aligned].median(), pl.fit_factor[pl.aligned].min(), pl.fit_factor[pl.aligned].max(),
pl.dist_fit[pl.aligned].median(), pl.dist_lag0[pl.aligned].median(), pl.dist_all_frames[pl.aligned].median()))
print("aligned sessions with a factor between 0.987 and 0.989: %d of %d" % (pl.fit_factor[pl.aligned].between(0.987, 0.989).sum(), pl.aligned.sum()))
pose read and clocks fitted in 157 s; sessions without pose: ['H28', 'H32']
subject pose_s port_entries lag1_s lag2_s lag3_s fit_offset_s fit_factor resid1_s resid2_s resid3_s dist_fit dist_lag0 dist_all_frames
H01 11024 619 -6.0 -31.4 -77.0 -28.2 0.9876 0.2 -0.0 -0.2 28.3 51.6 63.9
H02 11079 520 -27.8 -41.8 -122.8 -58.7 0.9875 -0.0 -0.0 0.4 39.3 100.3 105.1
H07 11140 554 -8.8 33.6 -106.8 -55.8 0.9875 -0.2 -0.0 0.2 36.8 93.1 94.1
H08 11092 1640 -37.6 -73.4 -114.4 -73.7 0.9877 -0.0 -0.0 -0.0 29.0 34.7 35.1
H09 10982 885 -40.2 -83.8 -119.2 -73.7 0.9874 -0.0 -0.0 0.2 36.4 87.3 86.5
H10 11042 1640 -30.2 -82.8 -105.4 -82.9 0.9875 -0.0 -0.0 -0.0 27.7 31.8 34.8
H11 10947 2171 177.0 -192.2 148.6 243.3 1.0030 3.4 19.2 -0.0 37.6 39.1 39.2
H12 11151 1597 47.0 4.2 -35.4 2.7 0.9877 0.2 -0.0 -0.2 24.6 36.2 39.0
H13 11274 576 2.4 -65.4 -62.4 -56.4 0.9876 -0.0 -0.0 -0.0 31.1 53.7 62.4
H14 11272 854 175.2 168.2 -80.2 292.7 0.9590 0.6 -2.2 -11.6 24.5 29.7 34.9
H17 10891 1914 -32.8 -75.2 -104.6 -75.3 0.9876 -0.0 3.6 -0.0 32.2 38.7 40.0
H19 11007 1456 -40.8 -79.6 -127.4 -83.9 0.9876 -0.4 -0.0 -0.0 36.1 67.2 89.6
H22 10963 1328 -30.0 -61.6 -114.0 -62.8 0.9883 -0.0 -0.0 -0.2 43.7 105.7 114.2
H23 11117 125 79.2 136.2 -80.8 169.8 0.9703 -1.6 -2.0 -11.6 37.8 41.9 49.3
H24 11418 1558 186.8 -51.4 -102.0 -69.3 0.9881 -0.0 -0.0 -0.0 38.4 40.4 46.8
H25 11345 708 33.4 -52.0 10.2 -49.7 1.0109 -14.8 -0.0 5.0 46.4 50.3 97.2
H27 11345 1614 -15.4 -15.4 -15.4 -15.5 1.0000 -0.0 -0.0 -0.0 45.2 80.8 96.0
H29 11327 714 84.8 127.4 -83.0 151.5 0.9655 -9.8 -0.0 -0.6 151.4 206.2 212.0
H31 11586 1082 -88.2 -137.8 -158.4 -141.2 0.9884 -0.2 -0.0 0.2 40.7 136.7 187.4
H33 11028 580 -37.2 -54.4 -192.4 -64.8 0.9884 -0.2 -0.0 -0.0 44.5 137.2 137.2
H34 11027 1582 -47.0 -103.0 -111.2 -83.4 0.9885 0.4 -0.0 0.4 36.9 158.8 149.0
H35 11586 646 -106.8 -58.2 -105.6 -71.3 0.9884 -0.0 -0.2 -0.0 41.9 232.6 181.3
H36 10987 1394 -37.4 -71.8 -132.0 -81.3 0.9885 -0.0 -0.0 0.2 39.8 173.1 179.6
H37 11062 1876 -41.4 -68.0 -112.6 -74.3 0.9884 -0.2 -0.0 -0.0 39.6 113.3 137.0
H38 11060 899 -55.8 -87.6 -136.6 -93.6 0.9885 0.4 -0.0 -0.0 38.3 209.4 223.0
H39 10909 606 -171.0 -89.2 -27.4 -129.7 0.9721 -0.6 -2.6 -14.0 154.9 226.4 220.9
H40 10913 1226 -102.6 -79.4 186.0 207.1 0.9939 -0.2 -0.8 2.8 166.3 203.4 221.3
H41 10946 686 -50.6 -106.2 -126.4 -92.8 0.9884 -0.2 -0.0 -0.2 39.7 204.9 205.9
H42 10939 1065 -81.8 -104.4 -131.0 -110.4 0.9884 -0.0 -0.2 -0.0 35.8 198.5 180.9
nominal clock: median |lag| 80.2 s; sessions whose best lag changes by more than 5 s within the session: 28 of 29
linear map: sessions with every third within 1 s of the map: 21 of 29; not: ['H11', 'H14', 'H17', 'H23', 'H25', 'H29', 'H39', 'H40']
fitted factor in the aligned sessions: median 0.9883, range 0.9874 to 1.0000; median distance at port entries 38.3 px under the map, 100.3 px at lag 0, 105.1 px over all frames
aligned sessions with a factor between 0.987 and 0.989: 20 of 21
# The video assets: frame count, frame rate and duration from each file's container header
# (streamed, no frame is decoded), against the pose series and the electrophysiology recording.
import av
def video_info(a):
out = {"subject": a["path"].split("/")[0][4:], "file": a["path"].split("/")[1][:40], "container": a["path"][-3:]}
try:
rf = nb.RangeFile(f"{nb.API}/assets/{a['asset_id']}/download/")
f = av.open(rf)
v = f.streams.video[0]
out.update({"frames": v.frames, "rate": float(v.average_rate), "duration_s": round(f.duration / 1e6, 1) if f.duration else np.nan,
"MB_read": round(rf.fetched / 1e6, 1)})
f.close()
except Exception as ex:
out["error"] = str(ex)[:60]
return out
VID = [a for a in A if not a["path"].endswith(".nwb")]
with ThreadPoolExecutor(8) as ex:
vi = pd.DataFrame(list(ex.map(video_info, VID)))
rec = {e["subject"]: e["duration"] for e in E}
vi["recording_s"] = vi.subject.map(rec).round(1)
vi["pose_s"] = vi.subject.map(dict(zip(pl.subject, pl.pose_s)))
print(vi.sort_values("subject").to_string(index=False))
print("mkv files report no frame count in their header (0); their rate is the container's nominal one.")
subject file container frames rate duration_s MB_read recording_s pose_s
H01 sub-H01_ses-1-2019-07-13-11-59-42-DiscS_ mp4 165366 15.0 11024.4 2.0 11887.8 11024.0
H02 sub-H02_ses-2-2019-07-14-11-14-47-DiscS_ mp4 166185 15.0 11079.0 1.9 12090.6 11079.0
H07 sub-H07_ses-7-2019-07-13-18-44-41-DiscS_ mp4 179366 15.0 11957.7 2.1 12788.5 11140.0
H08 sub-H08_ses-8-2019-07-14-18-11-14-DiscS_ mp4 166376 15.0 11091.7 2.0 12048.2 11092.0
H09 sub-H09_ses-9-2019-08-30-16-32-13-Disc4_ mp4 187658 15.0 12510.5 2.0 12223.7 10982.0
H10 sub-H10_ses-10-2019-08-30-16-32-20-Disc4 mp4 187638 15.0 12509.2 2.2 12220.9 11042.0
H11 sub-H11_ses-11-2019-08-30-12-58-04-Disc4 mp4 170599 15.0 11373.3 1.8 12270.7 10947.0
H13 sub-H13_ses-13-2019-08-30-09-21-02-Disc4 mp4 169115 15.0 11274.3 1.8 12261.5 11274.0
H14 sub-H14_ses-14-2019-08-30-09-21-05-Disc4 mp4 169078 15.0 11271.9 2.0 12261.4 11272.0
H17 sub-H17_ses-H17-2019-10-31-14-10-52-Disc mp4 163360 15.0 10890.7 1.8 12117.8 10891.0
H19 sub-H19_ses-H19-2019-10-31-17-44-38-Disc mp4 165108 15.0 11007.2 1.9 12434.4 11007.0
H22 sub-H22_ses-H22-2019-12-27-13-58-44-Disc mp4 164449 15.0 10963.3 1.2 12026.2 10963.0
H23 sub-H23_ses-H23-2019-12-29-09-55-50-Disc mp4 166763 15.0 11117.5 1.4 12286.8 11117.0
H24 sub-H24_ses-H24-2019-12-27-10-25-11-Disc mp4 171278 15.0 11418.5 1.4 12451.5 11418.0
H25 sub-H25_ses-H25-2020-02-20-10-06-39-Disc mkv 0 30.0 11344.5 0.1 12306.8 11345.0
H27 sub-H27_ses-H27-2020-02-20-10-07-01-Disc mkv 0 30.0 11344.6 0.1 12300.0 11345.0
H28 sub-H28_ses-H28-2020-02-19-14-27-39-Disc mkv 0 30.0 11240.3 0.1 12113.8 NaN
H28 sub-H28_ses-28-Disc4DLC-resnet50-Hao-Med mkv 0 30.0 11327.0 0.1 12113.8 NaN
H29 sub-H29_ses-H29-2020-02-20-13-42-59-Disc mkv 0 30.0 11326.6 0.1 12831.4 11327.0
H31 sub-H31_ses-H31-2020-04-30-09-56-31-Disc mp4 173857 15.0 11590.5 1.3 12186.5 11586.0
H32 sub-H32_ses-32-Disc6DLC-resnet50-Hao-Med mp4 165802 15.0 11053.5 1.1 11643.5 NaN
H32 sub-H32_ses-H32-2020-04-23-13-59-47-Disc mp4 165880 15.0 11058.7 1.3 11643.5 NaN
H33 sub-H33_ses-H33-2020-04-15-12-44-39-Disc mp4 165486 15.0 11032.4 1.2 11637.6 11028.0
H34 sub-H34_ses-H34-2020-04-15-12-44-32-Disc mp4 165492 15.0 11032.8 1.2 11690.5 11027.0
H35 sub-H35_ses-H35-2020-04-30-09-57-13-Disc mp4 173882 15.0 11592.1 1.2 12152.4 11586.0
H36 sub-H36_ses-H36-2020-04-24-12-47-33-Disc mp4 164916 15.0 10994.4 1.1 11616.6 10987.0
H37 sub-H37_ses-H37-2020-04-29-10-04-51-Disc mp4 165971 15.0 11064.7 1.2 11590.0 11062.0
H38 sub-H38_ses-H38-2020-04-29-10-04-41-Disc mp4 165979 15.0 11065.3 1.2 11609.6 11060.0
H39 sub-H39_ses-H39-2020-04-29-13-33-40-Disc mp4 163737 15.0 10915.8 1.1 11534.6 10909.0
H40 sub-H40_ses-H40-2020-04-29-13-33-54-Disc mp4 163736 15.0 10915.7 1.4 11524.5 10913.0
H41 sub-H41_ses-H41-2020-04-29-16-55-26-Disc mp4 164227 15.0 10948.5 1.4 11531.7 10946.0
H42 sub-H42_ses-H42-2020-04-29-16-55-04-Disc mp4 164228 15.0 10948.5 1.1 11558.0 10939.0
mkv files report no frame count in their header (0); their rate is the container's nominal one.
With the nominal clock (15 frames per second from 0 s) the head is closest to the port up to 192 s away from the port entries (median 80 s), and the best lag drifts within the session in 28 of 29 sessions (columns lag1 to lag3). A linear map with an offset and a rate factor removes the drift in 21 of the 29 sessions: under it, the best extra lag in each third is within 1 s of zero, 20 of those 21 fitted factors lie between 0.9874 and 0.9885 (H27 is at 1.0000), and the median head-to-port distance at port entries is 38 pixels against 105 over all frames. A common factor of about 0.988 is what a video whose true frame rate is about 1.2% below the nominal 15 per second would give. In H11, H14, H17, H23, H25, H29, H39 and H40 the map leaves a lag of 2.8 to 19.2 s in some third, or the distance stays near that over all frames (H29, H39, H40).
The video headers (next output) give a constant 15 frames per second for the mp4 files and durations within 10 s of the pose series, except H07, H09 and H10, whose videos run 818 to 1528 s longer than their pose, and H11, 426 s longer; the mkv files give no frame count. So the drift is not visible in the files' own metadata; the original videos were .asf files (PoseEstimation original_videos), whose timing is lost. The pose can be put on the trial clock in the 21 aligned sessions, by a map estimated from the data, not given by the files. This is an assumption for the behavior claim, which uses only those sessions. The map is consistent with the data, not shown by them.
# Photometry: the 20 NT sensor sessions.
PP = sorted(a["path"] for a in A if "ophys" in a["path"])
t = time.time()
with ThreadPoolExecutor(4) as ex:
S = list(ex.map(nb.load_photometry, PP))
print(f"loaded {len(S)} sessions in {time.time() - t:.0f} s")
ph = []
for s in S:
sid = s["session_id"]
lick_rel = np.concatenate([s["licks"][(s["licks"] > c - 2) & (s["licks"] < c + 6)] - c for c in s["cs_on"]])
ph.append({"subject": s["subject"], "task": "sucrose" if "Reward" in sid else "airpuff",
"phase": "early" if "early" in sid else "late", "fluorophore": s["fluorophore"][:14],
"fiber": s["fiber_location"], "nm": s["wavelengths"], "frames": s["n_frames"],
"dt_ms": tuple(round(x * 1000, 1) for x in s["frame_dt"]),
"phase_means": tuple(round(x, 3) for x in s["phase_means"]),
"ph415": s["phase415"], "ph470": s["phase470"],
"agree": round(min(s["phase_agree415"], s["phase_agree470"]), 3),
"labels_moved": s["labels_moved415"] + s["labels_moved470"],
"phase_means_after_irregular": tuple(round(x, 3) for x in s["phase_means_after_irregular"]),
"cues": len(s["cs_on"]), "cue_s": round(float(np.median(s["cs_dur"])), 2),
"min_gap_s": round(float(np.diff(s["cs_on"]).min()), 1),
"lick_frames_-2_0": int(((lick_rel >= -2) & (lick_rel < 0)).sum()),
"lick_frames_2_4": int(((lick_rel >= 2) & (lick_rel < 4)).sum())})
ph = pd.DataFrame(ph)
print(ph.to_string(index=False))
loaded 20 sessions in 45 s
subject task phase fluorophore fiber nm frames dt_ms phase_means ph415 ph470 agree labels_moved phase_means_after_irregular cues cue_s min_gap_s lick_frames_-2_0 lick_frames_2_4
M-0 airpuff early neurotensin-fl BLA [415.0, 470.0] 138159 (16.5, 16.8) (0.201, 0.012, 0.039) 1 2 1.000 0 (nan, nan, nan) 29 2.0 4.7 8 13
M-0 airpuff late neurotensin-fl BLA [415.0, 470.0] 251595 (16.5, 16.8) (0.202, 0.013, 0.039) 1 2 1.000 0 (nan, nan, nan) 60 2.0 29.7 2367 2338
M-0 sucrose early neurotensin-fl BLA [415.0, 470.0] 281228 (16.5, 16.8) (0.133, 0.015, 0.033) 1 2 1.000 0 (nan, nan, nan) 59 2.0 29.6 103 289
M-0 sucrose late neurotensin-fl BLA [415.0, 470.0] 247903 (16.5, 16.8) (0.198, 0.013, 0.04) 1 2 1.000 0 (nan, nan, nan) 59 2.0 29.7 139 489
M-1 airpuff early neurotensin-fl BLA [415.0, 470.0] 268252 (16.5, 16.8) (0.2, 0.015, 0.042) 1 2 1.000 0 (nan, nan, nan) 60 2.0 29.7 2400 2373
M-1 airpuff late neurotensin-fl BLA [415.0, 470.0] 257196 (16.5, 16.8) (0.2, 0.015, 0.043) 1 2 1.000 0 (nan, nan, nan) 60 2.0 29.6 0 0
M-1 sucrose early neurotensin-fl BLA [415.0, 470.0] 258206 (16.6, 16.8) (0.199, 0.015, 0.045) 1 2 1.000 0 (nan, nan, nan) 59 2.0 29.7 49 146
M-1 sucrose late neurotensin-fl BLA [415.0, 470.0] 251391 (16.5, 16.8) (0.201, 0.017, 0.045) 1 2 1.000 0 (nan, nan, nan) 60 2.0 29.7 654 1405
M-3 airpuff early neurotensin-fl BLA [415.0, 470.0] 313891 (11.3, 16.8) (0.201, 0.019, 0.048) 1 2 1.000 0 (nan, nan, nan) 60 2.0 29.7 7 0
M-3 airpuff late neurotensin-fl BLA [415.0, 470.0] 251673 (16.5, 16.8) (0.2, 0.015, 0.05) 1 2 1.000 0 (nan, nan, nan) 60 2.0 29.7 70 129
M-3 sucrose early neurotensin-fl BLA [415.0, 470.0] 245497 (16.5, 16.8) (0.196, 0.024, 0.082) 1 2 1.000 0 (nan, nan, nan) 59 2.0 29.7 0 130
M-3 sucrose late neurotensin-fl BLA [415.0, 470.0] 266708 (16.5, 16.8) (0.198, 0.026, 0.059) 1 2 1.000 0 (nan, nan, nan) 59 2.0 29.7 109 1417
M-4 airpuff early neurotensin-fl BLA [415.0, 470.0] 263618 (16.5, 16.8) (0.201, 0.014, 0.033) 1 2 1.000 0 (nan, nan, nan) 60 2.0 29.7 4 86
M-4 airpuff late neurotensin-fl BLA [415.0, 470.0] 273117 (11.4, 16.8) (0.2, 0.01, 0.03) 1 2 1.000 0 (nan, nan, nan) 60 2.0 29.7 103 113
M-4 sucrose early neurotensin-fl BLA [415.0, 470.0] 249451 (16.5, 16.8) (0.133, 0.012, 0.027) 1 2 1.000 0 (nan, nan, nan) 59 2.0 29.5 76 564
M-4 sucrose late neurotensin-fl BLA [415.0, 470.0] 256056 (16.5, 16.8) (0.2, 0.013, 0.033) 1 2 1.000 0 (nan, nan, nan) 60 2.0 29.6 72 538
M-5 airpuff early neurotensin-fl BLA [415.0, 470.0] 263757 (16.5, 16.8) (0.201, 0.022, 0.053) 1 2 1.000 0 (nan, nan, nan) 60 2.0 29.7 1150 1243
M-5 airpuff late neurotensin-fl BLA [415.0, 470.0] 248811 (16.5, 16.8) (0.202, 0.023, 0.054) 1 2 1.000 0 (nan, nan, nan) 60 2.0 29.7 0 0
M-5 sucrose early neurotensin-fl BLA [415.0, 470.0] 291288 (16.5, 16.8) (0.131, 0.021, 0.048) 1 2 1.000 0 (nan, nan, nan) 61 2.0 29.5 52 560
M-5 sucrose late neurotensin-fl BLA [415.0, 470.0] 250581 (2.3, 16.8) (0.2, 0.019, 0.057) 1 2 0.989 165289 (0.2, 0.019, 0.057) 59 2.0 29.6 50 328
Checked against the paper. The files hold the GRAB-NTS1.0 sensor ("neurotensin-fluorescent sensor") in the BLA of five wild-type mice, with 415 and 470 nm excitation, which is the experiment of Figure 3j-l (N = 5; the Methods exclude one of six mice for fiber placement, and the subjects are Mouse1, 2, 4, 5 and 6). Frames come every 16.7 ms (60 Hz), not at the Methods' 40 Hz: one series interleaves three frame types. The LED labels mark every 415 nm frame and every 470 nm frame (a third of the frames each); these sit on a fixed phase of three ("agree" is the fraction of labelled frames on the most common phase). In Mouse6's late sucrose session (M-5) the labels move by one frame after a 2.3 ms frame interval ("labels_moved"), while the data keep their order: the mean of each phase after that interval ("phase_means_after_irregular") is the same as over the session, with the near-constant third frame type still on phase 0. Following the labels there would read that constant frame type as the 470 nm channel, so the channels are taken from the frame order, at the phase the labels give before the first irregular interval. The third frame type is a near-constant value around 0.13 to 0.20 that no metadata describes. Each session holds about 60 cues lasting 1.9 to 2 s, about 30 s apart. The exception is Mouse1's early airpuff session (M-0): it holds 29 cues, and its last three come 4.7 and 7.4 s apart, one of them a 4.7 s bout, which looks like a disrupted end of the session. These cues are among the "all trials" readings; the "first 16 trials" readings do not reach them. The Methods give a 4 s CS with the US 2 s after onset and about half of the CSs unpaired; the files mark one kind of cue, 2 s long, and no US. In sucrose sessions licking rises after the cue ends (frames with a lick 2 to 4 s after onset against the 2 s before), which places the sucrose at the end of the 2 s cue, as the Methods' US time of 2 s gives; no file marks whether a cue was paired or unpaired, nor the airpuff. All cues are therefore analysed as CS-US trials, an assumption recorded with the claims.
# Per-trial sucrose check: in the late sucrose sessions, how many cues are followed by licking
# within 2 to 5 s of onset (a paired trial should be, an unpaired one need not)?
for s in S:
if "Reward" in s["session_id"]:
n = [np.sum((s["licks"] >= c + 2) & (s["licks"] < c + 5)) > 0 for c in s["cs_on"]]
print(f'{s["session_id"]:34s} cues followed by a lick within 2-5 s: {sum(n)} of {len(n)}')
RewardLearning_earlyAcq_Mouse1 cues followed by a lick within 2-5 s: 50 of 59 RewardLearning_lateAcq_Mouse1 cues followed by a lick within 2-5 s: 57 of 59 RewardLearning_earlyAcq_Mouse2 cues followed by a lick within 2-5 s: 34 of 59 RewardLearning_lateAcq_Mouse2 cues followed by a lick within 2-5 s: 55 of 60 RewardLearning_earlyAcq_Mouse4 cues followed by a lick within 2-5 s: 17 of 59 RewardLearning_lateAcq_Mouse4 cues followed by a lick within 2-5 s: 58 of 59 RewardLearning_earlyAcq_Mouse5 cues followed by a lick within 2-5 s: 54 of 59 RewardLearning_lateAcq_Mouse5 cues followed by a lick within 2-5 s: 53 of 60 RewardLearning_earlyAcq_Mouse6 cues followed by a lick within 2-5 s: 24 of 61 RewardLearning_lateAcq_Mouse6 cues followed by a lick within 2-5 s: 49 of 59
In the late sucrose sessions most cues are followed by licking, so the files' cues are mostly, and perhaps all, paired with sucrose; the early sessions, before learning, cannot tell.
Photometry claims (Figure 3j-l)¶
Processing as the claims' criterion gives it (Methods, Fiber photometry, Data analysis): 415 and 470 nm channels median-filtered over 200 ms and averaged into 200 ms (5 Hz) bins from -5 to 10 s around each CS onset; the trial-averaged 470 nm trace regressed on the trial-averaged 415 nm trace over -5 to 0 s (the 415 nm fit subtracted from the 470 nm trace), the residuals z-scored against -5 to 0 s. Windows: CS 0 to 2 s, US 2 to 7 s (Figure 3 legend). Early and late are the files' early and late acquisition sessions. Readings run, as the criterion lists them: (a) first 16 trials or (b) all trials; (i) area under the z-scored curve or (ii) mean z. One more detail the Methods leave open is run both ways: the z-score uses each trial's own baseline ("trial") or the session's pooled baseline ("session"). The regression is fitted on the trial-averaged baseline, as the Methods say; the fitted slope, intercept and correlation of each session are printed, and a fit on the baseline bins of every trial pooled is run as a variant, since the averaged baseline may hold too little variance to fit.
TRIALS, FITS = {}, []
for s in S:
key = (s["subject"], "sucrose" if "Reward" in s["session_id"] else "airpuff", "early" if "early" in s["session_id"] else "late")
for zm in ("trial", "session"):
tb, Z, keep, fit = nb.photometry_trials(s, zmode=zm, return_fit=True)
TRIALS[key + (zm,)] = (tb, Z)
tb, Z, keep, fitp = nb.photometry_trials(s, zmode="trial", fit="pooled", return_fit=True)
TRIALS[key + ("trial, pooled fit",)] = (tb, Z)
FITS.append({"mouse": key[0], "task": key[1], "phase": key[2], "slope": fit[0], "intercept": fit[1], "r": fit[2],
"pooled_slope": fitp[0], "pooled_intercept": fitp[1], "pooled_r": fitp[2]})
print("470 nm on 415 nm regression over -5 to 0 s, per session (trial-averaged fit, and pooled-trials variant):")
print(pd.DataFrame(FITS).round(3).to_string(index=False))
print("trials kept per session (all cues have a full -5 to 10 s window):",
{f"{k[0]} {k[1]} {k[2]}": v[1].shape[0] for k, v in TRIALS.items() if k[3] == "trial"})
mice = sorted({k[0] for k in TRIALS})
def measure(task, phase, zm, sub, meas, lo, hi):
out = []
for m in mice:
tb, Z = TRIALS[(m, task, phase, zm)]
Z = Z[:16] if sub == "first16" else Z
out.append(nb.window_measure(tb, Z, lo, hi, meas).mean())
return np.array(out)
def compare(task, zmodes=("trial", "session")):
rows = []
for zm in zmodes:
for sub in ("first16", "all"):
for meas in ("auc", "mean"):
r = {"zscore": zm, "trials": sub, "measure": meas}
for win, (lo, hi) in (("US", (2, 7)), ("CS", (0, 2))):
e, l = measure(task, "early", zm, sub, meas, lo, hi), measure(task, "late", zm, sub, meas, lo, hi)
tt = stats.ttest_rel(l, e)
d = l - e
ci = stats.t.interval(0.95, len(d) - 1, d.mean(), stats.sem(d))
r.update({f"{win}_late-early": d.mean(), f"{win}_ci_low": ci[0], f"{win}_ci_high": ci[1],
f"{win}_t": tt.statistic, f"{win}_p": tt.pvalue, f"{win}_mice_in_direction": None})
r[f"{win}_mice_up"] = int((d > 0).sum())
rows.append(r)
return pd.DataFrame(rows).drop(columns=["US_mice_in_direction", "CS_mice_in_direction"])
470 nm on 415 nm regression over -5 to 0 s, per session (trial-averaged fit, and pooled-trials variant):
mouse task phase slope intercept r pooled_slope pooled_intercept pooled_r
M-0 airpuff early 0.375 0.034 0.248 1.434 0.021 0.995
M-0 airpuff late 0.805 0.028 0.467 1.487 0.019 0.994
M-0 sucrose early 0.928 0.019 0.563 1.249 0.014 0.990
M-0 sucrose late 1.054 0.026 0.745 1.433 0.021 0.996
M-1 airpuff early 1.613 0.018 0.877 1.150 0.025 0.994
M-1 airpuff late 1.203 0.025 0.868 1.218 0.025 0.996
M-1 sucrose early 0.939 0.030 0.936 1.670 0.020 0.991
M-1 sucrose late 0.737 0.032 0.480 1.432 0.020 0.997
M-3 airpuff early 2.331 0.002 0.972 1.166 0.025 0.992
M-3 airpuff late 2.056 0.019 0.972 1.444 0.029 0.989
M-3 sucrose early 1.458 0.047 0.614 2.058 0.033 0.990
M-3 sucrose late 1.857 0.012 0.966 1.689 0.016 0.993
M-4 airpuff early 0.865 0.021 0.819 1.249 0.016 0.997
M-4 airpuff late 1.604 0.015 0.943 1.148 0.019 0.985
M-4 sucrose early 1.559 0.008 0.821 1.099 0.014 0.998
M-4 sucrose late 2.018 0.008 0.935 1.229 0.018 0.958
M-5 airpuff early 2.189 0.005 0.978 1.432 0.022 0.989
M-5 airpuff late 2.185 0.004 0.958 1.382 0.022 0.980
M-5 sucrose early 1.712 0.012 0.903 1.615 0.014 0.986
M-5 sucrose late 2.168 0.017 0.988 1.246 0.034 0.938
trials kept per session (all cues have a full -5 to 10 s window): {'M-0 airpuff early': 29, 'M-0 airpuff late': 60, 'M-0 sucrose early': 59, 'M-0 sucrose late': 59, 'M-1 airpuff early': 60, 'M-1 airpuff late': 60, 'M-1 sucrose early': 59, 'M-1 sucrose late': 60, 'M-3 airpuff early': 60, 'M-3 airpuff late': 60, 'M-3 sucrose early': 59, 'M-3 sucrose late': 59, 'M-4 airpuff early': 60, 'M-4 airpuff late': 60, 'M-4 sucrose early': 59, 'M-4 sucrose late': 60, 'M-5 airpuff early': 60, 'M-5 airpuff late': 60, 'M-5 sucrose early': 61, 'M-5 sucrose late': 59}
fig3k-nt-sensor-sucrose-increase¶
The BLA NT sensor response to sucrose delivery increases after reward learning. (Paper: US window late against early, paired t-test t(4) = 8.709, P = 0.001, effect 4.14, CI 2.744 to 5.538; CS window t(4) = 1.788, P = 0.1483.) Reproduced when the US comparison has p < 0.05 with late higher under every reading.
res_suc = compare("sucrose")
print(res_suc.round(4).to_string(index=False))
print("\nPer mouse, US window AUC (trial z-score, first 16 trials):")
print(pd.DataFrame({"mouse": mice, "early": measure("sucrose", "early", "trial", "first16", "auc", 2, 7),
"late": measure("sucrose", "late", "trial", "first16", "auc", 2, 7)}).round(2).to_string(index=False))
ok = (res_suc["US_p"] < 0.05) & (res_suc["US_late-early"] > 0)
print("\nReadings with US p < 0.05 and late higher: %d of %d; per-trial z-score %d of 4, session z-score %d of 4"
% (ok.sum(), len(res_suc), ok[res_suc.zscore == "trial"].sum(), ok[res_suc.zscore == "session"].sum()))
print("\nVariant, regression fitted on the pooled baseline bins of every trial (per-trial z-score):")
print(compare("sucrose", ("trial, pooled fit",))[["zscore", "trials", "measure", "US_late-early", "US_ci_low", "US_ci_high", "US_t", "US_p", "US_mice_up"]].round(4).to_string(index=False))
zscore trials measure US_late-early US_ci_low US_ci_high US_t US_p US_mice_up CS_late-early CS_ci_low CS_ci_high CS_t CS_p CS_mice_up
trial first16 auc 2.9428 -0.7733 6.6589 2.1987 0.0928 4 0.6614 -2.2730 3.5958 0.6258 0.5654 3
trial first16 mean 0.5886 -0.1547 1.3318 2.1987 0.0928 4 0.3307 -1.1365 1.7979 0.6258 0.5654 3
trial all auc 3.5294 1.8236 5.2352 5.7445 0.0046 5 0.5455 -3.0270 4.1180 0.4240 0.6934 3
trial all mean 0.7059 0.3647 1.0470 5.7445 0.0046 5 0.2728 -1.5135 2.0590 0.4240 0.6934 3
session first16 auc -2.9182 -9.5474 3.7111 -1.2222 0.2887 2 -1.2509 -3.8352 1.3334 -1.3439 0.2501 2
session first16 mean -0.5836 -1.9095 0.7422 -1.2222 0.2887 2 -0.6255 -1.9176 0.6667 -1.3439 0.2501 2
session all auc 0.1574 -0.0985 0.4133 1.7074 0.1629 4 -0.0076 -0.2250 0.2098 -0.0967 0.9276 2
session all mean 0.0315 -0.0197 0.0827 1.7074 0.1629 4 -0.0038 -0.1125 0.1049 -0.0967 0.9276 2
Per mouse, US window AUC (trial z-score, first 16 trials):
mouse early late
M-0 -7.85 -3.59
M-1 -8.94 -6.50
M-3 0.52 2.71
M-4 -0.67 6.33
M-5 1.44 0.27
Readings with US p < 0.05 and late higher: 2 of 8; per-trial z-score 2 of 4, session z-score 0 of 4
Variant, regression fitted on the pooled baseline bins of every trial (per-trial z-score):
zscore trials measure US_late-early US_ci_low US_ci_high US_t US_p US_mice_up
trial, pooled fit first16 auc 2.5313 -5.7652 10.8277 0.8471 0.4446 4
trial, pooled fit first16 mean 0.5063 -1.1530 2.1655 0.8471 0.4446 4
trial, pooled fit all auc 1.4543 -6.4857 9.3942 0.5085 0.6378 4
trial, pooled fit all mean 0.2909 -1.2971 1.8788 0.5085 0.6378 4
fig, axes = plt.subplots(1, 2, figsize=(10, 3.2), sharey=True)
for ax, task in zip(axes, ("sucrose", "airpuff")):
for phase, col in (("early", "tab:blue"), ("late", "tab:red")):
tr = np.array([TRIALS[(m, task, phase, "trial")][1][:16].mean(0) for m in mice])
tb = TRIALS[(mice[0], task, phase, "trial")][0] + 0.1
ax.plot(tb, tr.mean(0), color=col, label=phase)
ax.fill_between(tb, tr.mean(0) - stats.sem(tr), tr.mean(0) + stats.sem(tr), color=col, alpha=0.2)
print(f"{task} {phase}: plotted mean z from {tr.mean(0).min():.2f} to {tr.mean(0).max():.2f}")
ax.axvspan(2, 7, color="0.9"); ax.axvline(0, color="k", lw=0.5)
ax.set_title(f"NT sensor, {task} (first 16 trials, mean of 5 mice)"); ax.set_xlabel("s from CS onset")
axes[0].set_ylabel("z"); axes[0].legend(); plt.tight_layout(); plt.show()
sucrose early: plotted mean z from -1.26 to 0.27 sucrose late: plotted mean z from -0.76 to 0.45 airpuff early: plotted mean z from -0.63 to 0.59 airpuff late: plotted mean z from -1.74 to 0.23
Result: not reproduced. Under the per-trial z-score the US response to sucrose is larger in the late session under all four readings of the criterion, but p < 0.05 holds only with all trials: AUC +3.53 (95% CI 1.82 to 5.24), t(4) = 5.74, p = 0.0046, 5 of 5 mice higher, inside the paper's CI of 2.744 to 5.538. With the first 16 trials it is AUC +2.94, p = 0.093, 4 of 5 mice higher. AUC and mean give the same t and p, so the two passing readings are one result, and the criterion's "p < 0.05 with late higher under every reading" holds in 2 of its 4 readings: the criterion is not met. The significant result also rests on the regression fitted on the trial-averaged baseline, whose correlation is 0.48 to 0.75 in four sucrose sessions (M-0 early and late, M-1 late, M-3 early); fitted on the pooled baseline bins of every trial (correlation at least 0.938) no reading reaches p < 0.05 (p >= 0.44), and under the session-pooled z-score none does either (p >= 0.16). The CS window does not change under any reading (p >= 0.25), as in the paper. The per-trial z-score is the reading adopted (an assumption: the Methods z-score "using data a baseline window of -5 s to 0 s", most simply read as each trial's own window). The criterion has no partial condition, so with p < 0.05 failing in half its readings, and in every variant of the processing, the claim is not reproduced; only the direction holds throughout.
fig3l-nt-sensor-airpuff-decrease¶
The BLA NT sensor response to airpuff decreases after punishment learning. (Paper: US window t(4) = 3.996, P = 0.0162, effect -3.657, CI -6.083 to -0.9583; CS t(4) = 2.137, P = 0.0994.) Reproduced when the US comparison has p < 0.05 with late lower under every reading.
res_air = compare("airpuff")
print(res_air.round(4).to_string(index=False))
print("\nPer mouse, US window AUC (trial z-score, first 16 trials):")
print(pd.DataFrame({"mouse": mice, "early": measure("airpuff", "early", "trial", "first16", "auc", 2, 7),
"late": measure("airpuff", "late", "trial", "first16", "auc", 2, 7)}).round(2).to_string(index=False))
print("\nReadings with US p < 0.05 and late lower: %d of %d" % (((res_air["US_p"] < 0.05) & (res_air["US_late-early"] < 0)).sum(), len(res_air)))
print("\nVariant, regression fitted on the pooled baseline bins of every trial (per-trial z-score):")
print(compare("airpuff", ("trial, pooled fit",))[["zscore", "trials", "measure", "US_late-early", "US_ci_low", "US_ci_high", "US_t", "US_p", "US_mice_up"]].round(4).to_string(index=False))
zscore trials measure US_late-early US_ci_low US_ci_high US_t US_p US_mice_up CS_late-early CS_ci_low CS_ci_high CS_t CS_p CS_mice_up
trial first16 auc -5.4694 -12.3523 1.4135 -2.2062 0.0920 0 -0.8698 -2.3015 0.5619 -1.6868 0.1669 1
trial first16 mean -1.0939 -2.4705 0.2827 -2.2062 0.0920 0 -0.4349 -1.1508 0.2809 -1.6868 0.1669 1
trial all auc -4.6570 -10.0465 0.7325 -2.3991 0.0744 1 -0.7447 -1.6415 0.1521 -2.3055 0.0824 0
trial all mean -0.9314 -2.0093 0.1465 -2.3991 0.0744 1 -0.3724 -0.8208 0.0761 -2.3055 0.0824 0
session first16 auc -1.1811 -11.1338 8.7715 -0.3295 0.7583 2 -0.3223 -4.4390 3.7945 -0.2173 0.8386 2
session first16 mean -0.2362 -2.2268 1.7543 -0.3295 0.7583 2 -0.1611 -2.2195 1.8972 -0.2173 0.8386 2
session all auc -0.6074 -1.3853 0.1705 -2.1679 0.0960 0 -0.1427 -0.3445 0.0591 -1.9637 0.1210 0
session all mean -0.1215 -0.2771 0.0341 -2.1679 0.0960 0 -0.0714 -0.1722 0.0295 -1.9637 0.1210 0
Per mouse, US window AUC (trial z-score, first 16 trials):
mouse early late
M-0 -1.36 -2.51
M-1 -6.08 -18.11
M-3 -0.33 -3.33
M-4 10.77 -0.07
M-5 -0.27 -0.61
Readings with US p < 0.05 and late lower: 0 of 8
Variant, regression fitted on the pooled baseline bins of every trial (per-trial z-score):
zscore trials measure US_late-early US_ci_low US_ci_high US_t US_p US_mice_up
trial, pooled fit first16 auc 0.2685 -3.6100 4.1470 0.1922 0.8569 3
trial, pooled fit first16 mean 0.0537 -0.7220 0.8294 0.1922 0.8569 3
trial, pooled fit all auc -1.5799 -4.6941 1.5343 -1.4085 0.2318 2
trial, pooled fit all mean -0.3160 -0.9388 0.3069 -1.4085 0.2318 2
Result: not reproduced. Under the per-trial z-score the US response to airpuff is lower in the late session under all four readings (AUC -5.47 with the first 16 trials, 5 of 5 mice lower, and -4.66 with all trials; the paper's effect is -3.657), but no reading reaches p < 0.05 (smallest p = 0.074, t(4) = -2.40). The session-pooled z-score gives p >= 0.096, and the pooled-baseline fit p >= 0.23. The CS window does not change significantly under any reading (p >= 0.08).
fig3-nt-sensor-early-dip¶
Early in acquisition, the BLA NT sensor signal transiently decreases after both sucrose and airpuff. Reproduced when the mean z-scored signal over mice in the US window (2 to 7 s from CS onset) is below zero in early trials for both sucrose and airpuff.
rows = []
for zm in ("trial", "session"):
for sub in ("first16", "all"):
for meas in ("mean", "auc"):
r = {"zscore": zm, "trials": sub, "measure": meas}
for task in ("sucrose", "airpuff"):
v = measure(task, "early", zm, sub, meas, 2, 7)
r[f"{task}_mean"] = v.mean(); r[f"{task}_mice_below_0"] = int((v < 0).sum())
r[f"{task}_p_vs_0"] = stats.ttest_1samp(v, 0).pvalue
rows.append(r)
dip = pd.DataFrame(rows)
print(dip.round(3).to_string(index=False))
print("\nReadings with both early means below zero: %d of %d" % (((dip.sucrose_mean < 0) & (dip.airpuff_mean < 0)).sum(), len(dip)))
print("\nEarly sessions, per mouse, mean z 2-7 s (trial z-score, first 16 trials):")
print(pd.DataFrame({"mouse": mice, "sucrose": measure("sucrose", "early", "trial", "first16", "mean", 2, 7),
"airpuff": measure("airpuff", "early", "trial", "first16", "mean", 2, 7)}).round(2).to_string(index=False))
# Variant for the word "transient" (asked of a reader): the lowest point, 2 to 7 s, of the trace
# averaged over trials and then over mice, and of each mouse's trial-averaged trace.
print("\nVariant: minimum of the early trace in 2 to 7 s (z), mean of 5 mice and per mouse:")
for zm in ("trial", "session"):
for sub in ("first16", "all"):
r = {}
for task in ("sucrose", "airpuff"):
tb = TRIALS[(mice[0], task, "early", zm)][0] + 0.1
w = (tb > 2) & (tb < 7)
tr = np.array([(TRIALS[(m, task, "early", zm)][1][:16] if sub == "first16" else TRIALS[(m, task, "early", zm)][1]).mean(0) for m in mice])
r[task] = f"mean-trace min {tr.mean(0)[w].min():+.2f} at {tb[w][np.argmin(tr.mean(0)[w])]:.1f} s; per mouse {np.round(tr[:, w].min(1), 2).tolist()}"
print(f"{zm:8s} {sub:8s} sucrose: {r['sucrose']}\n airpuff: {r['airpuff']}")
zscore trials measure sucrose_mean sucrose_mice_below_0 sucrose_p_vs_0 airpuff_mean airpuff_mice_below_0 airpuff_p_vs_0
trial first16 mean -0.620 3 0.231 0.109 4 0.854
trial first16 auc -3.099 3 0.231 0.544 4 0.854
trial all mean -0.399 4 0.351 0.192 4 0.726
trial all auc -1.994 4 0.351 0.960 4 0.726
session first16 mean 0.344 2 0.557 -0.397 3 0.573
session first16 auc 1.719 2 0.557 -1.987 3 0.573
session all mean -0.018 4 0.655 0.014 3 0.641
session all auc -0.088 4 0.655 0.068 3 0.641
Readings with both early means below zero: 0 of 8
Early sessions, per mouse, mean z 2-7 s (trial z-score, first 16 trials):
mouse sucrose airpuff
M-0 -1.57 -0.27
M-1 -1.79 -1.22
M-3 0.10 -0.07
M-4 -0.13 2.15
M-5 0.29 -0.05
Variant: minimum of the early trace in 2 to 7 s (z), mean of 5 mice and per mouse:
trial first16 sucrose: mean-trace min -1.09 at 4.3 s; per mouse [-3.0, -2.44, -1.07, -0.82, -0.24]
airpuff: mean-trace min -0.63 at 3.7 s; per mouse [-3.81, -2.47, -0.77, -0.11, -0.27]
trial all sucrose: mean-trace min -0.82 at 4.3 s; per mouse [-2.64, -1.67, -1.59, -0.39, 0.65]
airpuff: mean-trace min -0.36 at 3.7 s; per mouse [-2.94, -1.48, -1.04, 0.15, -0.27]
session first16 sucrose: mean-trace min +0.31 at 4.3 s; per mouse [0.7, 1.33, 1.23, -1.45, -0.43]
airpuff: mean-trace min -0.46 at 3.7 s; per mouse [0.47, -1.58, -1.46, 1.4, -1.43]
session all sucrose: mean-trace min -0.05 at 4.3 s; per mouse [-0.27, -0.07, -0.05, -0.03, 0.08]
airpuff: mean-trace min -0.03 at 3.7 s; per mouse [-0.23, -0.09, -0.03, -0.02, -0.03]
Result: not reproduced. In the early sessions the mean over mice of the US-window signal is below zero for sucrose under the per-trial z-score (-0.62 with the first 16 trials, -0.40 with all trials) but above zero for airpuff (+0.11 and +0.19; M-4 at +2.15), so the criterion, which needs both below zero, is not met under any per-trial reading. Under the session-pooled z-score the signs change with the trial subset (0 of 8 readings in all). No mean differs from zero (p >= 0.23). Read as a brief dip, the paper's word "transient", the mouse-averaged trace does go below zero within 2 to 7 s for both outcomes under the per-trial z-score (lowest -1.09 for sucrose and -0.63 for airpuff, first 16 trials); whether the criterion's window mean is the right measure for a transient decrease is asked of a reader.
Electrophysiology claims (Figure 4)¶
fig4-recording-count¶
685 BLA neurons were recorded from 30 mice during the three-cue discrimination task. Units from the last discrimination-task session of each mouse, with baseline firing of at least 0.3 Hz (Methods). Whether 685 is counted before or after that exclusion is not said; both counts are shown and either may match (expected 685 +/- 68.5 neurons, 30 +/- 3 mice).
Baseline rate: the mean rate in the 1 s before each CS onset, over all CSs of the session (the PSTH baseline of the Methods); the whole-session rate is shown as a variant.
base_rate, sess_rate, unit_rows = {}, {}, []
for e in E:
on = np.sort(np.concatenate([e["trials"][v][0] for v in nb.CS_TYPES.values() if v in e["trials"]]))
for u, sp in enumerate(e["spikes"]):
b = np.mean(np.searchsorted(sp, on) - np.searchsorted(sp, on - 1.0))
unit_rows.append({"subject": e["subject"], "unit": u, "baseline_hz": b, "session_hz": len(sp) / e["duration"]})
U = pd.DataFrame(unit_rows)
n_all, n_base, n_sess = len(U), int((U.baseline_hz >= 0.3).sum()), int((U.session_hz >= 0.3).sum())
print(f"units before exclusion: {n_all}; with baseline >= 0.3 Hz: {n_base}; with session rate >= 0.3 Hz: {n_sess}")
print(f"mice: {U.subject.nunique()}")
print("per mouse (all units / kept by baseline):")
print(U.groupby("subject").agg(units=("unit", "size"), kept=("baseline_hz", lambda x: int((x >= 0.3).sum()))).T.to_string())
print(f"\nexpected 685 +/- 68.5 neurons: before exclusion {abs(n_all - 685) <= 68.5}, after {abs(n_base - 685) <= 68.5}")
print(f"expected 30 +/- 3 mice: {abs(U.subject.nunique() - 30) <= 3}")
print("units per mouse: median %.0f (paper 685 / 30 = 22.8)" % U.groupby("subject").size().median())
ratio = ses[["sucrose", "neutral", "shock"]]
near = (ratio.neutral > 0) & ((ratio.sucrose - 2 * ratio.neutral).abs() <= 2) & ((ratio.shock - ratio.neutral).abs() <= 2)
print(f"sessions with sucrose, neutral and shock CSs near 2:1:1 (within 2 trials): {int(near.sum())} of {len(ses)}; not: {list(ses.subject[~near])}")
units before exclusion: 1902; with baseline >= 0.3 Hz: 1739; with session rate >= 0.3 Hz: 1884 mice: 31 per mouse (all units / kept by baseline): subject H01 H02 H07 H08 H09 H10 H11 H12 H13 H14 H17 H19 H22 H23 H24 H25 H27 H28 H29 H31 H32 H33 H34 H35 H36 H37 H38 H39 H40 H41 H42 units 54 64 64 60 62 64 64 64 64 64 64 64 64 64 46 58 38 64 64 64 64 64 64 58 64 58 64 64 64 62 62 kept 47 61 64 60 62 56 55 60 43 64 0 63 52 63 46 58 38 63 61 53 64 64 64 51 64 48 64 63 64 62 62 expected 685 +/- 68.5 neurons: before exclusion False, after False expected 30 +/- 3 mice: True units per mouse: median 64 (paper 685 / 30 = 22.8) sessions with sucrose, neutral and shock CSs near 2:1:1 (within 2 trials): 25 of 31; not: ['H09', 'H10', 'H11', 'H12', 'H13', 'H14']
Result: partially reproduced. The number of mice (31) is within the tolerance of 30. The unit count is not: 1902 units, 1739 with a baseline rate of at least 0.3 Hz (1884 by whole-session rate), against the paper's 685; neither is within 68.5 of 685. The units are not sorted single units (data check), which fits a median of 64 per mouse instead of 23. The 2:1:1 ratio of the statement holds in 25 of 31 sessions; H09 to H14 have no neutral CS.
fig4d-seven-clusters¶
Unsupervised hierarchical clustering of BLA responses to the sucrose, neutral and shock CSs gives seven functional clusters. Per neuron, 20 ms PSTHs from -1 to 1 s around CS onset, z-scored with the 1 s baseline, smoothed with a Gaussian-weighted moving average (smoothing factor 0.85), the three CS traces concatenated; Ward linkage with correlation distance cut at 35% of the maximum linkage distance; readings (a) Ward on the correlation distance matrix, (b) Ward with Euclidean distance on centred, unit-norm rows. Reproduced when one reading gives seven clusters.
Units with a baseline rate below 0.3 Hz are excluded (Methods). The six sessions without neutral CSs cannot give the concatenated trace and are left out; so are units whose baseline SD is zero. The "smoothing factor of 0.85" sets the window through a MATLAB heuristic the paper does not give, so the smoothing window is run over a range (in 20 ms bins, smoothing each CS trace separately).
feats, owner = [], []
dropped = {"low_rate": 0, "no_neutral_session": 0, "zero_sd": 0}
for e in E:
on = {k: e["trials"][v][0] for k, v in nb.CS_TYPES.items() if v in e["trials"]}
allon = np.sort(np.concatenate(list(on.values())))
for u, sp in enumerate(e["spikes"]):
b = np.mean(np.searchsorted(sp, allon) - np.searchsorted(sp, allon - 1.0))
if b < 0.3: dropped["low_rate"] += 1; continue
if "neutral" not in on: dropped["no_neutral_session"] += 1; continue
z = np.array([nb.zscore_psth(nb.psth(sp, on[k])) for k in ("sucrose", "neutral", "shock")])
if np.isnan(z).any(): dropped["zero_sd"] += 1; continue
feats.append(z); owner.append(e["subject"])
F = np.array(feats); owner = np.array(owner)
print("units clustered:", len(F), "from", len(set(owner)), "mice; dropped:", dropped)
def cluster(F, window):
X = nb.gauss_smooth(F, window).reshape(len(F), -1)
Za = linkage(pdist(X, "correlation"), "ward")
ca = fcluster(Za, 0.35 * Za[:, 2].max(), "distance")
Xc = X - X.mean(1, keepdims=True); Xc /= np.linalg.norm(Xc, axis=1, keepdims=True)
Zb = linkage(Xc, "ward")
cb = fcluster(Zb, 0.35 * Zb[:, 2].max(), "distance")
return ca, cb
rows = []
CL = {}
for w in (1, 3, 5, 7, 9, 11, 15, 21, 31, 41, 51, 75, 101):
ca, cb = cluster(F, w)
CL[w] = (ca, cb)
rows.append({"window_bins": w, "window_ms": w * 20, "clusters_a_corr_ward": ca.max(), "clusters_b_unitnorm_ward": cb.max(),
"sizes_a": sorted(np.bincount(ca)[1:].tolist(), reverse=True)})
cl = pd.DataFrame(rows)
print(cl.to_string(index=False))
print("\nwindows giving 7 clusters: reading (a)", list(cl.window_bins[cl.clusters_a_corr_ward == 7]),
" reading (b)", list(cl.window_bins[cl.clusters_b_unitnorm_ward == 7]))
units clustered: 1399 from 24 mice; dropped: {'low_rate': 163, 'no_neutral_session': 340, 'zero_sd': 0}
window_bins window_ms clusters_a_corr_ward clusters_b_unitnorm_ward sizes_a
1 20 9 11 [338, 243, 192, 183, 162, 83, 78, 61, 59]
3 60 12 14 [383, 292, 146, 105, 94, 84, 63, 57, 52, 44, 41, 38]
5 100 11 13 [364, 303, 119, 108, 105, 100, 97, 69, 58, 40, 36]
7 140 12 13 [381, 283, 125, 89, 88, 76, 73, 71, 67, 58, 48, 40]
9 180 11 13 [320, 291, 130, 129, 113, 91, 86, 81, 60, 58, 40]
11 220 10 12 [575, 191, 125, 107, 86, 83, 76, 64, 51, 41]
15 300 8 12 [601, 227, 137, 118, 101, 88, 68, 59]
21 420 9 11 [340, 238, 211, 153, 109, 106, 105, 77, 60]
31 620 6 8 [393, 249, 237, 208, 192, 120]
41 820 7 9 [385, 229, 188, 187, 178, 118, 114]
51 1020 6 9 [315, 268, 263, 263, 210, 80]
75 1500 6 8 [342, 285, 276, 201, 183, 112]
101 2020 8 9 [250, 245, 175, 158, 153, 148, 138, 132]
windows giving 7 clusters: reading (a) [41] reading (b) []
# Cut-off sensitivity for one window (11 bins = 220 ms): clusters as a function of the cut fraction.
Xw = nb.gauss_smooth(F, 11).reshape(len(F), -1)
Za = linkage(pdist(Xw, "correlation"), "ward")
print("reading (a), window 11 bins: clusters by cut fraction of the maximum linkage distance")
print({f: int(fcluster(Za, f * Za[:, 2].max(), "distance").max()) for f in (0.25, 0.3, 0.35, 0.4, 0.45, 0.5)})
# Per mouse composition of the clusters, reading (a), window 11.
ca = CL[11][0]
print("\nunits per mouse in each cluster (reading (a), window 11):")
print(pd.crosstab(owner, ca).to_string())
fig, ax = plt.subplots(figsize=(7, 4))
order = np.argsort(ca, kind="stable")
im = ax.imshow(Xw[order], aspect="auto", cmap="RdBu_r", vmin=-5, vmax=5, interpolation="nearest")
ax.set_xticks([0, 50, 100, 150, 200, 250]); ax.set_xticklabels(["S -1", "S 0", "N -1", "N 0", "Sh -1", "Sh 0"])
ax.set_ylabel("unit (sorted by cluster)"); plt.colorbar(im, label="z (clipped at +/-5)"); plt.show()
print("plotted z from %.1f to %.1f; %.1f%% of values outside +/-5" % (Xw.min(), Xw.max(), 100 * np.mean(np.abs(Xw) > 5)))
# Clusters dominated by single mice, and the count without the sessions whose CS alignment is not
# shown (H22, H27, H34, H39; data check), for every window, reading (a).
for w in (11, 41):
ct = pd.crosstab(owner, CL[w][0])
share = (ct.max(axis=0) / ct.sum(axis=0)).round(2)
print(f"window {w} bins, reading (a): largest share of one mouse in each cluster: {dict(zip(ct.columns, share))}; mouse: {dict(zip(ct.columns, ct.idxmax(axis=0)))}")
keep = ~np.isin(owner, ["H22", "H27", "H34", "H39"])
print("\nwithout H22, H27, H34 and H39 (%d units): clusters (a), (b) by window:" % keep.sum(),
{w: tuple(int(c.max()) for c in cluster(F[keep], w)) for w in (1, 5, 11, 21, 31, 41, 51, 75, 101)})
reading (a), window 11 bins: clusters by cut fraction of the maximum linkage distance
{0.25: 14, 0.3: 13, 0.35: 10, 0.4: 7, 0.45: 7, 0.5: 5}
units per mouse in each cluster (reading (a), window 11):
col_0 1 2 3 4 5 6 7 8 9 10
row_0
H01 0 0 0 3 0 0 38 6 0 0
H02 2 0 1 35 7 0 8 1 0 7
H07 0 0 0 14 0 5 9 14 0 22
H08 0 0 0 1 0 17 41 1 0 0
H19 0 0 0 0 1 33 4 2 0 23
H22 8 3 35 0 2 0 0 1 0 3
H23 2 0 7 0 0 7 7 4 29 7
H24 0 0 0 0 6 0 22 0 18 0
H25 1 0 0 5 47 0 0 5 0 0
H27 0 38 0 0 0 0 0 0 0 0
H28 2 0 2 0 2 0 54 0 0 3
H29 40 0 5 0 2 0 8 0 0 6
H31 2 0 0 1 3 0 42 0 0 5
H32 0 0 0 2 7 2 39 8 0 6
H33 58 0 0 1 1 0 0 4 0 0
H34 1 0 0 2 6 0 52 0 0 3
H35 0 0 0 4 1 0 46 0 0 0
H36 0 0 0 4 10 11 28 8 0 3
H37 6 0 7 5 6 0 8 2 0 14
H38 0 0 2 4 23 0 27 2 4 2
H39 0 0 1 0 29 0 33 0 0 0
H40 1 0 4 1 13 0 30 13 0 2
H41 0 0 0 0 18 0 44 0 0 0
H42 2 0 0 1 7 1 35 15 0 1
plotted z from -8.6 to 23.4; 0.7% of values outside +/-5
window 11 bins, reading (a): largest share of one mouse in each cluster: {1: 0.46, 2: 0.93, 3: 0.55, 4: 0.42, 5: 0.25, 6: 0.43, 7: 0.09, 8: 0.17, 9: 0.57, 10: 0.21}; mouse: {1: 'H33', 2: 'H27', 3: 'H22', 4: 'H02', 5: 'H25', 6: 'H19', 7: 'H28', 8: 'H42', 9: 'H23', 10: 'H19'}
window 41 bins, reading (a): largest share of one mouse in each cluster: {1: 0.23, 2: 0.18, 3: 0.11, 4: 0.25, 5: 0.29, 6: 0.2, 7: 0.26}; mouse: {1: 'H07', 2: 'H02', 3: 'H28', 4: 'H41', 5: 'H25', 6: 'H23', 7: 'H33'}
without H22, H27, H34 and H39 (1182 units): clusters (a), (b) by window: {1: (9, 11), 5: (11, 12), 11: (8, 12), 21: (7, 9), 31: (6, 9), 41: (6, 9), 51: (6, 7), 75: (6, 6), 101: (6, 7)}
Result: blocked (methods_underspecified). With the paper's rule (Ward linkage, correlation distance, cut at 35% of the maximum linkage distance) the 1399 units of 24 mice give 8 to 12 clusters under reading (a) and 11 to 14 under reading (b) for smoothing windows of 20 to 420 ms, and 6 to 8 under (a) for windows of 620 ms to 2 s; seven clusters come only under reading (a) with an 820 ms window. The count depends on the smoothing window, which the paper sets by MATLAB's "smoothing factor of 0.85" without the window it produced, so no reading can be preferred; the count also moves with the cut fraction (14 at 25%, 10 at 35%, 7 at 40%, 5 at 50% for a 220 ms window). Without H22, H27, H34 and H39, whose CS alignment is not shown, seven clusters come under (a) at 420 ms and under (b) at 1020 ms and 2 s. Some clusters are single sessions (at 220 ms one cluster is 93% H27), so a count here would not show functional types. The units are not the paper's sorted neurons.
fig4c-photoidentification¶
BLA-NAc and BLA-CeA neurons are identified by short-latency responses to light at the end of each session. The Methods' rule (1 ms PSTH 50 ms around light onset, z-scored against the 50 ms before; excitation above z = 1.98 within 6 ms, none 7 to 50 ms after with a 6 ms moving average, none in the fake session); readings (a) all pulse onsets of the trains and the onset of the long light, (b) the 1 Hz pulses only. Reproduced when the counts are within 10% of 39 BLA-NAc and 32 BLA-CeA neurons under one reading.
No file says which projection (NAc or CeA) carries the Cav2-Cre in each mouse, so the two counts cannot be split. Their sum can still be tested: both counts are within 10% only if the total is within 63.9 to 78.1 (0.9 and 1.1 times 71). The fake-session test needs to know which of the two light blocks is the fake one. The Methods put the fake session first ("a phototagging session was conducted after a 'fake' phototagging session in which the patch cord was disconnected", Methods, In vivo electrophysiology), so block 1 is taken as the fake one; the reverse is run as a reading, and for H31 to H42, where one block comes before the task, either order per session is allowed (the lowest total over those choices is printed). Light events of 5 and 10 ms are the pulses; the files' 100 ms events are the only longer light, and their onsets are added in reading (a).
tag_rows = []
for e in E:
on, off = e["trials"]["phototagging"]
dur = off - on
g = np.argmax(np.diff(on))
block = [np.arange(len(on)) <= g, np.arange(len(on)) > g]
def pulses(b, reading):
p = on[block[b] & (dur <= 0.011)]
if reading == "a":
return np.sort(np.r_[p, on[block[b] & (np.abs(dur - 0.1) < 0.01)]])
return p[(np.r_[np.inf, np.diff(p)] > 0.9) & (np.r_[np.diff(p), np.inf] > 0.9)]
base = U[U.subject == e["subject"]].baseline_hz.values
r = {"subject": e["subject"], "units": len(e["spikes"])}
for rd in ("a", "b"):
P = [pulses(0, rd), pulses(1, rd)]
r[f"{rd}_pulses"] = len(P[1])
# nb.photoresponsive(sp, real, fake) is "responsive to real and not to fake"; computed once per block.
resp = [np.array([nb.photoresponsive(sp, P[b])[0] for sp in e["spikes"]]) for b in (0, 1)]
for name, (fk, rl) in {"fake1": (0, 1), "fake2": (1, 0)}.items():
ok = resp[rl] & ~resp[fk]
r[f"{rd}_{name}"] = int(ok.sum()); r[f"{rd}_{name}_0.3Hz"] = int((ok & (base >= 0.3)).sum())
r[f"{rd}_resp_block1"] = int(resp[0].sum()); r[f"{rd}_resp_block2"] = int(resp[1].sum())
tag_rows.append(r)
tg = pd.DataFrame(tag_rows)
print(tg.to_string(index=False))
late = tg.subject.isin(list(span.subject[span.light_before_task > 0]))
print("\nlight event durations (s) and counts, first session:", dict(zip(*np.unique(np.round(E[0]["trials"]["phototagging"][1] - E[0]["trials"]["phototagging"][0], 3), return_counts=True))))
print("\ntagged units in total (paper 39 + 32 = 71; both counts within 10% needs 63.9 to 78.1):")
for rd in ("a", "b"):
for sfx in ("", "_0.3Hz"):
f1, f2 = tg[f"{rd}_fake1{sfx}"], tg[f"{rd}_fake2{sfx}"]
lowest = int(f1[~late].sum() + np.minimum(f1[late], f2[late]).sum())
print(f"reading ({rd}){' baseline >= 0.3 Hz' if sfx else ''}: block 1 fake {f1.sum()}, block 2 fake {f2.sum()}, "
f"block 1 fake with either order in H31-H42 (lowest) {lowest}; "
f"responsive without the fake test: block 1 {tg[f'{rd}_resp_block1'].sum()}, block 2 {tg[f'{rd}_resp_block2'].sum()}")
for rd in ("a", "b"):
print(f"reading ({rd}): units responsive to block 1 {tg[f'{rd}_resp_block1'].sum()} of {tg.units.sum()}, to block 2 {tg[f'{rd}_resp_block2'].sum()} of {tg.units.sum()}")
print("units responsive to both blocks (reading a), share of those responsive to block 2: %.2f"
% ((tg.a_resp_block2 - tg.a_fake1).sum() / tg.a_resp_block2.sum()))
print("mice with no tagged unit (reading a, block 1 fake):", list(tg.subject[tg.a_fake1 == 0]))
subject units a_pulses a_fake1 a_fake1_0.3Hz a_fake2 a_fake2_0.3Hz a_resp_block1 a_resp_block2 b_pulses b_fake1 b_fake1_0.3Hz b_fake2 b_fake2_0.3Hz b_resp_block1 b_resp_block2
H01 54 426 0 0 1 1 53 52 150 4 4 5 5 48 47
H02 64 426 4 4 2 2 45 47 150 8 7 4 4 47 51
H07 64 426 0 0 10 10 63 53 150 1 1 5 5 61 57
H08 60 426 0 0 4 4 60 56 150 0 0 16 16 60 44
H09 62 426 2 2 1 1 59 60 150 0 0 5 5 61 56
H10 64 426 0 0 11 11 64 53 150 0 0 20 19 64 44
H11 64 426 13 12 16 16 47 44 150 17 14 13 13 30 34
H12 64 426 1 1 2 2 63 62 150 3 3 6 6 60 57
H13 64 426 4 4 6 6 38 36 150 11 9 11 6 33 33
H14 64 426 0 0 3 3 33 30 150 0 0 0 0 0 0
H17 64 426 0 0 0 0 0 0 150 0 0 0 0 0 0
H19 64 426 0 0 41 41 64 23 150 0 0 38 38 62 24
H22 64 426 17 13 23 18 29 23 150 20 16 19 18 28 29
H23 64 426 3 3 32 31 34 5 150 9 9 30 30 45 24
H24 46 428 8 8 0 0 32 40 151 0 0 0 0 46 46
H25 58 426 0 0 10 10 51 41 150 1 1 1 1 56 56
H27 38 426 4 4 1 1 27 30 150 3 3 0 0 35 38
H28 64 426 2 2 6 6 38 34 150 0 0 8 7 56 48
H29 64 426 2 1 3 3 44 43 150 0 0 6 6 61 55
H31 64 426 4 2 15 14 46 35 150 0 0 0 0 64 64
H32 64 426 23 23 1 1 38 60 150 1 1 0 0 63 64
H33 64 426 10 10 2 2 42 50 150 4 4 0 0 60 64
H34 64 426 9 9 9 9 42 42 150 12 12 2 2 46 56
H35 58 426 0 0 6 5 52 46 150 5 4 10 7 45 40
H36 64 426 10 10 4 4 41 47 150 2 2 5 5 53 50
H37 58 426 8 4 4 4 31 35 150 14 9 1 1 36 49
H38 64 426 2 2 9 9 48 41 150 6 6 3 3 52 55
H39 64 426 4 3 17 17 60 47 150 0 0 0 0 64 64
H40 64 426 4 4 6 6 34 32 150 5 5 6 6 40 39
H41 62 426 7 7 8 8 45 44 150 0 0 10 10 62 52
H42 62 426 5 5 3 3 30 32 150 7 7 9 9 36 34
light event durations (s) and counts, first session: {np.float64(0.005): np.int64(480), np.float64(0.01): np.int64(300), np.float64(0.1): np.int64(72), np.float64(0.228): np.int64(1), np.float64(0.279): np.int64(1)}
tagged units in total (paper 39 + 32 = 71; both counts within 10% needs 63.9 to 78.1):
reading (a): block 1 fake 146, block 2 fake 256, block 1 fake with either order in H31-H42 (lowest) 104; responsive without the fake test: block 1 1353, block 2 1243
reading (a) baseline >= 0.3 Hz: block 1 fake 133, block 2 fake 248, block 1 fake with either order in H31-H42 (lowest) 95; responsive without the fake test: block 1 1353, block 2 1243
reading (b): block 1 fake 133, block 2 fake 233, block 1 fake with either order in H31-H42 (lowest) 102; responsive without the fake test: block 1 1474, block 2 1374
reading (b) baseline >= 0.3 Hz: block 1 fake 117, block 2 fake 222, block 1 fake with either order in H31-H42 (lowest) 91; responsive without the fake test: block 1 1474, block 2 1374
reading (a): units responsive to block 1 1353 of 1902, to block 2 1243 of 1902
reading (b): units responsive to block 1 1474 of 1902, to block 2 1374 of 1902
units responsive to both blocks (reading a), share of those responsive to block 2: 0.88
mice with no tagged unit (reading a, block 1 fake): ['H01', 'H07', 'H08', 'H10', 'H14', 'H17', 'H19', 'H25', 'H35']
Result: blocked, metadata insufficient. With the fake-session test the rule tags 91 to 256 units in total, depending on the reading of the pulses, which block is the fake one, and the 0.3 Hz exclusion: the fewest are 91 (1 Hz pulses, either block order allowed in H31 to H42, baseline at least 0.3 Hz), and with the Methods' order (block 1 fake) 146 and 133 for readings (a) and (b). All are above 78.1. But the test cannot be applied as the Methods mean it. Block 1, taken as the fake session, evokes a short-latency response in 1353 of 1902 units (reading a), and block 2 in 1243; 88% of the units responsive to block 2 also respond to block 1. A session with the patch cord disconnected should evoke no light response, so either neither block is the fake session as the Methods describe it, or the responses are light artefacts on units that are not sorted. The 91 to 256 "tagged" units are then units that cross the threshold in one block and not in the other, and their count does not test the claim. Without labels for the fake and real blocks (and for the projection of each mouse), the claim cannot be tested.
Behavior claim (Figure 5d-h)¶
fig5dh-behavior-clusters¶
Unsupervised trial-by-trial clustering of pose-derived behavior gives seven clusters for sucrose trials and five for shock trials. Eight measures for the 300 frames (20 s) after CS onset; tsfresh features (about 6000); PCA, t-SNE on the top 30 PCs, a kernel density over the t-SNE map (seaborn kdeplot), local maxima as watershed markers; sucrose and shock trials separately. t-SNE perplexity, kernel bandwidth and seed are not given: the defaults are used with several seeds. Reproduced when the default settings give 7 and 5 clusters for most seeds.
The trials come from the sessions where the fitted pose clock holds (data check), with the CS onsets of the trial table mapped onto the pose clock. Choices the paper does not make, each shown: tsfresh 0.16.0 does not install with this Python (3.13), so the current tsfresh with its default, comprehensive feature set is used (6264 features per trial, against the paper's "~6000"); features that are not finite or constant are dropped. Whether the features are scaled before PCA is not said; both are run. The z-scored features are the reading adopted, since on the raw features the PCA sees only the few features with the largest scales and the cluster counts move with changes of 1e-9 in the features (shown below). "The defaults of scikit-learn" are read two ways: the current defaults (PCA initialisation, automatic learning rate), as installed, and those of the scikit-learn releases contemporary with tsfresh 0.16.0 (2020: random initialisation, learning rate 200). The kernel density is seaborn's default (scipy gaussian_kde, Scott's bandwidth, 200 x 200 grid, 3 bandwidths beyond the data), and each local maximum of the grid is one watershed region.
# Pose measures for the 300 frames after each sucrose and each shock CS, aligned sessions only.
BT = {"sucrose": [], "shock": []}; BW = {"sucrose": [], "shock": []}; lost = {"sucrose": 0, "shock": 0}
Ed = {e["subject"]: e for e in E}
for subj in pl.subject[pl.aligned]:
d = PC[subj][1]
for k in BT:
on = Ed[subj]["trials"][nb.CS_TYPES[k]][0]
X, kept = nb.behavior_trials(d["measures"], d["clock"], d["rate"], d["t0"], on)
BT[k].append(X); BW[k] += [subj] * len(X); lost[k] += len(on) - len(kept)
BT = {k: np.concatenate(v) for k, v in BT.items()}; BW = {k: np.array(v) for k, v in BW.items()}
print({k: v.shape for k, v in BT.items()}, "trials dropped (frames beyond the pose):", lost)
print("trials per mouse:", pd.DataFrame({k: pd.Series(BW[k]).value_counts() for k in BW}).fillna(0).astype(int).T.to_string())
for k in BT:
print(k, "measure ranges:", {m: (round(float(np.nanmin(BT[k][:, j])), 1), round(float(np.nanmax(BT[k][:, j])), 1)) for j, m in enumerate(nb.POSE_MEASURES)})
{'sucrose': (1319, 8, 300), 'shock': (688, 8, 300)} trials dropped (frames beyond the pose): {'sucrose': 0, 'shock': 0}
trials per mouse: H01 H02 H07 H08 H09 H10 H12 H13 H19 H22 H24 H27 H31 H33 H34 H35 H36 H37 H38 H41 H42
sucrose 60 60 60 60 75 75 75 75 60 60 60 60 60 60 60 60 60 60 60 59 60
shock 30 30 30 30 45 45 45 45 30 30 30 30 30 30 28 30 30 30 30 30 30
sucrose measure ranges: {'velocity': (0.0, 10223.0), 'distance_to_port': (1.9, 596.8), 'direction': (-3.1, 3.1), 'rel_direction': (-180.0, 180.0), 'acceleration': (-151274.8, 153204.8), 'd_distance_to_port': (-8044.7, 8053.2), 'd_direction': (-47.1, 47.1), 'd_rel_direction': (-2699.8, 2699.2)}
shock measure ranges: {'velocity': (0.0, 9844.5), 'distance_to_port': (0.7, 603.2), 'direction': (-3.1, 3.1), 'rel_direction': (-180.0, 180.0), 'acceleration': (-146097.2, 147608.1), 'd_distance_to_port': (-7873.5, 8410.2), 'd_direction': (-47.1, 47.1), 'd_rel_direction': (-2699.2, 2699.9)}
import logging
from tsfresh import extract_features
from tsfresh.feature_extraction import ComprehensiveFCParameters
logging.getLogger("tsfresh").setLevel(logging.ERROR)
FEAT = {}
for k, X in BT.items():
t = time.time()
df = pd.DataFrame({"id": np.repeat(np.arange(len(X)), X.shape[2]), "time": np.tile(np.arange(X.shape[2]), len(X))})
for j, m in enumerate(nb.POSE_MEASURES):
df[m] = X[:, j, :].ravel()
with warnings.catch_warnings():
warnings.simplefilter("ignore")
Fk = extract_features(df, column_id="id", column_sort="time", default_fc_parameters=ComprehensiveFCParameters(),
n_jobs=4, disable_progressbar=True)
V = Fk.replace([np.inf, -np.inf], np.nan)
V = V.loc[:, V.notna().all() & (V.std() > 0)]
FEAT[k] = V
print(f"{k}: {len(X)} trials, {Fk.shape[1]} features, {V.shape[1]} finite and not constant ({time.time() - t:.0f} s)")
sucrose: 1319 trials, 6264 features, 5931 finite and not constant (443 s)
shock: 688 trials, 6264 features, 5927 finite and not constant (232 s)
from scipy import ndimage
from threadpoolctl import threadpool_info
def tsne_counts(PCs, seeds=range(5), store=None, **kw):
"""Density peaks of the t-SNE map per seed: all maxima, and maxima above 5% of the highest."""
counts, counts5 = [], []
for seed in seeds:
Y = TSNE(2, random_state=seed, **kw).fit_transform(PCs)
n, D = nb.density_peaks(Y); counts.append(n)
counts5.append(int(((D == ndimage.maximum_filter(D, size=3, mode="constant", cval=-np.inf)) & (D > 0.05 * D.max())).sum()))
if store is not None:
store[seed] = Y
return counts, counts5
target = {"sucrose": 7, "shock": 5}
rows = []
EMB = {}
for k, V in FEAT.items():
for scaling in ("z-scored", "raw"):
M = ((V - V.mean()) / V.std()).values if scaling == "z-scored" else V.values
pca = PCA(30, random_state=0).fit(M)
PCs = pca.transform(M)
for defaults, kw in (("current", {}), ("2020", {"init": "random", "learning_rate": 200.0})):
store = {}
counts, counts5 = tsne_counts(PCs, store=store, **kw)
EMB.update({(k, scaling, defaults, s): Y for s, Y in store.items()})
rows.append({"trials": k, "features": scaling, "tsne_defaults": defaults, "pca_solver": pca._fit_svd_solver,
"var_30pcs": round(pca.explained_variance_ratio_.sum(), 3), "clusters_by_seed": counts,
"clusters_peaks_above_5pct": counts5, "target": target[k], "seeds_at_target": sum(c == target[k] for c in counts)})
bc = pd.DataFrame(rows)
print(bc.to_string(index=False))
both = {(f, d): (bc[(bc.trials == "sucrose") & (bc.features == f) & (bc.tsne_defaults == d)].seeds_at_target.iloc[0] >= 3,
bc[(bc.trials == "shock") & (bc.features == f) & (bc.tsne_defaults == d)].seeds_at_target.iloc[0] >= 3)
for f in ("z-scored", "raw") for d in ("current", "2020")}
print("\n(sucrose 7 in most seeds, shock 5 in most seeds) per reading:", both)
print("readings with both:", [r for r, v in both.items() if all(v)])
for k in FEAT:
raw_var = FEAT[k].var().sort_values(ascending=False)
print(f"{k}, raw features: share of the total variance in the largest 1, 3 and 10 features: "
f"{raw_var.iloc[0] / raw_var.sum():.4f}, {raw_var.iloc[:3].sum() / raw_var.sum():.4f}, {raw_var.iloc[:10].sum() / raw_var.sum():.4f}; largest: {raw_var.index[0]}")
fig, axes = plt.subplots(1, 2, figsize=(9, 4))
for ax, k in zip(axes, ("sucrose", "shock")):
Y = EMB[(k, "z-scored", "current", 0)]
ax.scatter(Y[:, 0], Y[:, 1], s=3, c=pd.factorize(BW[k])[0], cmap="tab20")
ax.set_title(f"{k} trials, t-SNE (z-scored features, current defaults, seed 0)", fontsize=9); ax.set_xticks([]); ax.set_yticks([])
plt.tight_layout(); plt.show()
print("colours are mice; all points are plotted")
# Sensitivity to numerical detail that the seeds do not fix, with the current defaults: an exact PCA
# (svd_solver="full") in place of the randomized one PCA(30) chooses here, and each feature
# multiplied by (1 + 1e-9 * standard normal noise) before scaling. Seeds 0 to 2.
print("\nBLAS and OpenMP libraries (the PCA's arithmetic depends on them):",
[(i.get("internal_api"), i.get("version"), i.get("architecture"), i.get("num_threads")) for i in threadpool_info()])
rng_b = np.random.default_rng(5)
for k, V in FEAT.items():
noisy = V.values * (1 + 1e-9 * rng_b.standard_normal(V.shape))
for scaling in ("z-scored", "raw"):
for name, (Mv, solver) in {"exact PCA": (V.values, "full"), "features x (1 + 1e-9 noise)": (noisy, "auto")}.items():
Mv = (Mv - Mv.mean(0)) / Mv.std(0, ddof=1) if scaling == "z-scored" else Mv
counts, _ = tsne_counts(PCA(30, svd_solver=solver, random_state=0).fit_transform(Mv), seeds=range(3))
print(f"{k}, {scaling} features, {name}: clusters by seed {counts}")
# Variant: the adopted reading (z-scored features, current defaults) without H08, H10 and H24,
# whose head-to-port distance under the fitted clock is least below that over all frames.
weak = ["H08", "H10", "H24"]
print("\nhead-to-port distance (px) under the map / over all frames:", {s: (float(pl.set_index("subject").dist_fit[s]), float(pl.set_index("subject").dist_all_frames[s])) for s in weak})
for k, V in FEAT.items():
Vk = V[~np.isin(BW[k], weak)]
Vk = Vk.loc[:, Vk.std() > 0]
counts, _ = tsne_counts(PCA(30, random_state=0).fit_transform(((Vk - Vk.mean()) / Vk.std()).values))
print(f"{k} without {weak}: {len(Vk)} trials, clusters by seed {counts}, at target {sum(c == target[k] for c in counts)} of 5")
print(f"\nnotebook ran in {(time.time() - started) / 60:.1f} min")
trials features tsne_defaults pca_solver var_30pcs clusters_by_seed clusters_peaks_above_5pct target seeds_at_target
sucrose z-scored current randomized 0.421 [9, 9, 9, 9, 9] [9, 9, 9, 9, 9] 7 0
sucrose z-scored 2020 randomized 0.421 [7, 7, 8, 8, 7] [7, 7, 8, 8, 7] 7 3
sucrose raw current randomized 1.000 [7, 7, 7, 7, 7] [7, 7, 7, 7, 7] 7 5
sucrose raw 2020 randomized 1.000 [6, 4, 6, 6, 4] [6, 4, 6, 6, 4] 7 0
shock z-scored current randomized 0.409 [4, 4, 4, 4, 4] [4, 4, 4, 4, 4] 5 0
shock z-scored 2020 randomized 0.409 [4, 4, 4, 5, 3] [4, 4, 4, 5, 3] 5 1
shock raw current randomized 1.000 [5, 5, 5, 5, 5] [5, 5, 5, 5, 5] 5 5
shock raw 2020 randomized 1.000 [3, 4, 6, 4, 3] [3, 4, 6, 4, 3] 5 0
(sucrose 7 in most seeds, shock 5 in most seeds) per reading: {('z-scored', 'current'): (np.False_, np.False_), ('z-scored', '2020'): (np.True_, np.False_), ('raw', 'current'): (np.True_, np.True_), ('raw', '2020'): (np.False_, np.False_)}
readings with both: [('raw', 'current')]
sucrose, raw features: share of the total variance in the largest 1, 3 and 10 features: 0.6309, 0.9530, 1.0000; largest: acceleration__time_reversal_asymmetry_statistic__lag_1
shock, raw features: share of the total variance in the largest 1, 3 and 10 features: 0.6414, 0.9753, 1.0000; largest: acceleration__time_reversal_asymmetry_statistic__lag_1
colours are mice; all points are plotted
BLAS and OpenMP libraries (the PCA's arithmetic depends on them): [('openmp', None, None, 4), ('openblas', '0.3.34.106.0', 'SkylakeX', 4), ('openblas', '0.3.31.dev', 'SkylakeX', 4)]
sucrose, z-scored features, exact PCA: clusters by seed [7, 7, 7]
sucrose, z-scored features, features x (1 + 1e-9 noise): clusters by seed [7, 7, 7]
sucrose, raw features, exact PCA: clusters by seed [9, 9, 9]
sucrose, raw features, features x (1 + 1e-9 noise): clusters by seed [8, 8, 8]
shock, z-scored features, exact PCA: clusters by seed [3, 3, 3]
shock, z-scored features, features x (1 + 1e-9 noise): clusters by seed [3, 3, 3]
shock, raw features, exact PCA: clusters by seed [5, 5, 5]
shock, raw features, features x (1 + 1e-9 noise): clusters by seed [5, 5, 5]
head-to-port distance (px) under the map / over all frames: {'H08': (29.0, 35.1), 'H10': (27.7, 34.8), 'H24': (38.4, 46.8)}
sucrose without ['H08', 'H10', 'H24']: 1124 trials, clusters by seed [5, 5, 5, 5, 5], at target 0 of 5
shock without ['H08', 'H10', 'H24']: 583 trials, clusters by seed [3, 3, 3, 3, 3], at target 0 of 5 notebook ran in 24.4 min
Result: partially reproduced. On 1319 sucrose and 688 shock trials from the 21 sessions with an aligned pose clock, no reading gives both counts: each of the four readings meets one count and not the other. With z-scored features, sucrose trials give 7 clusters for 5 of 5 seeds under the installed defaults (7, 7, 8, 7, 8 under the 2020 defaults, 3 of 5), and shock trials give 3 for 5 of 5 (4, 3, 4, 3, 3 under the 2020 defaults), never 5. With raw features, sucrose trials give 6 for 5 of 5 (7, 5, 7, 5, 6 under the 2020 defaults) and shock trials 5 for 5 of 5 (5, 5, 5, 6, 7 under the 2020 defaults, 3 of 5). The raw-feature counts are not stable: on the raw features the 3 largest features carry 95% (sucrose) and 98% (shock) of the variance and the 10 largest all of it, and with the installed defaults an exact PCA gives 8 and 5 clusters, and multiplying each feature by 1 + 1e-9 noise gives 7 and 7. With z-scored features the same changes give 7 and 4 clusters: the sucrose count stays at 7 and the shock count moves from 3 to 4, below 5 either way. The reading adopted is z-scored features with the installed defaults (7 and 3): the sucrose count holds and the shock count does not, as with the 2020 defaults; with raw features it is the reverse. Without H08, H10 and H24, whose pose clock is checked least well (head-to-port distance 29.0, 27.7 and 38.4 px under the map against 35.1, 34.8 and 46.8 px over all frames), the adopted reading gives 5 and 4. The trials rest on a pose clock fitted to the data, and the count of the sucrose trials also moves with the sessions included.
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 |
|---|---|
| fig3k-nt-sensor-sucrose-increase | not reproduced (higher late in every reading, but p = 0.093 with the first 16 trials; p = 0.0046 with all trials only under the per-trial z-score and trial-averaged fit) |
| fig3l-nt-sensor-airpuff-decrease | not reproduced (right direction, p = 0.074 at best) |
| fig3-nt-sensor-early-dip | not reproduced (sucrose below zero, airpuff above) |
| fig4-recording-count | partially reproduced (31 mice against 30; 1902 or 1739 units against 685) |
| fig4c-photoidentification | blocked, metadata_insufficient (the "fake" block evokes responses in most units, and no block or projection labels) |
| fig4d-seven-clusters | blocked, methods_underspecified (7 clusters only with an 820 ms smoothing window) |
| fig5dh-behavior-clusters | partially reproduced (7 sucrose clusters, 3 shock clusters, adopted reading; no reading gives both) |
| 13 claims needing group or projection labels | blocked, metadata_insufficient |
| 19 claims on experiments not in the dandiset | blocked, data_missing |
How the analysis differs from the paper's. The photometry cues are taken as CS-US trials with the US 2 s after onset, since the files mark neither the US nor unpaired cues; the photometry z-score is run with per-trial and session baselines; the clustering smoothing window is run over a range since MATLAB's smoothing factor cannot be reproduced; the units are the files' units, not sorted single units; the fake phototagging block is taken as the first one, as the Methods order them; the pose is put on the trial clock by a linear map fitted to the port entries, in 21 of 29 sessions; tsfresh is the current version, not 0.16.0, even under the 2020 t-SNE defaults, and the t-SNE defaults and feature scaling are run both ways.
What was found about the dataset. The files carry no CRISPR or control label and no tagged projection, which most of the paper's electrophysiology findings need; the units are not sorted single units (frequent sub-millisecond intervals, spikes shared between units, 64 units in most sessions); H17's spike times stop at 249 s of a 12118 s recording; the pose has no timestamps and on its nominal clock of 15 frames per second the head reaches the port up to 192 s (median 80 s) away from the port entries, a drift that a rate factor of about 0.988 removes in 21 sessions; there are 31 mice where the paper has 30, and two of them may not be recorded on their last day; six sessions have no neutral CS though the description says 2:1:1; the two light blocks are not labelled fake and real, and both evoke short-latency responses in most units (1353 and 1243 of 1902); 27 files hold confocal images of the PVT with no channel or group label; the photometry series interleaves three frame types at 60 Hz in one series, in one session the LED labels move by one frame against the data, the cue events last 2 s where the Methods give 4 s, and no event marks the US or unpaired cues. Only the NT sensor experiment of the photometry is deposited; the behavioral cohorts, the axon-terminal and somatic photometry, the anatomy and the slice physiology are not.