Reproducing Chen et al. (2017) from dandiset 000015¶
Paper. Chen T-W, Li N, Daie K, Svoboda K (2017). A map of anticipatory activity in mouse motor
cortex. Neuron 94, 866-879. doi:10.1016/j.neuron.2017.05.005.
The claims file (claims/000015/10.1016-j.neuron.2017.05.005.yaml) was written from the version of
record, main text, figure legends and STAR Methods; the Supplemental Information was not in the text
read.
Dataset. DANDI:000015, version 0.220126.1914 (doi:10.48324/dandi.000015/0.220126.1914). The paper cites no DANDI version: it names a CRCNS deposit (10.6080/K0S180F9), and that DOI resolves to a placeholder record ("CRCNS.org ssc-8TBD", "This data is not yet available"). This is the only published version of the dandiset; the draft has the same 210 assets and size, so the published version is pinned. The files are streamed; nothing is downloaded.
Outside data. The paper compares its imaging with tagged PT and IT units from Li et al. (2015). Those units are on DANDI:000010, version 0.220126.1905, which is used for the link claim and for ephys-pt-it-onset-not-detected below (spike times, trial tables, behavioural events and laser power of the sessions with tagged units).
Questions. Five questions were asked of a reader (questions/000015.yaml) and answered. The
paper settles none of them; each answer wrote readings into the claims file, and every claim they
touch is run under them here: dF/F with a baseline F0 at or below zero (reading V1, the Methods'
formula, decides; V2, without the ROIs concerned, is reported), the onset rule (S1, the first bin
over threshold, and S2, the first bin from which it stays over threshold to the go cue), the
sample-epoch level and end of delay of figs4-peak-shift (L1, L2, E1, E2), the delay value of
fig6-lick-direction-spread (mean and maximum), and the method for the Li et al. (2015) units. A
sixth question, asked after the first verification and not yet answered, is over what window the
peak of figs4-peak-shift is taken; the whole extracted window is used, and the peak before the go
cue is shown.
What the dandiset holds. Two-photon imaging of 105 image planes (each deposited twice) from six mice, with the trial table, the sample, delay and go times, ROI traces and neuropil traces, an ROI table with PT, IT or unknown labels, and an area label (left ALM or left vM1, the paper's MM) per plane. It holds no plane depth, no stereotactic coordinates, no lick times, no whisker video, no wide-field imaging, no silicon probe recordings and no photoinhibition data.
| Claim | Tested here |
|---|---|
| early-lick-fraction, imaging-sample, trials-per-type, neuron-counts | Tested |
| fraction-selective-alm-mm, figs5-single-coefficient | Tested |
| fig4a-lick-direction-epochs, fig4b-object-location-memory, fig4c-outcome-after-lick, fig5-mixed-response-epoch | Tested |
| fig6-object-location-mm, fig6-lick-direction-spread (without its depth part), fig6-outcome-response-uniform, fig6h-preparatory-alm-object-mm | Tested |
| figs4-peak-shift | Tested |
| fig8c-pt-deep-only (counts only), fig8d-pt-contra-bias, fig8d-pt-earlier-than-it (all ALM IT only), fig8d-mm-posterior-preferring, fig8d-error-outnumber-correct, alm-pt-it-lick-fraction-below-20 | Tested |
| li2015-tagged-units-link (link claim added by this attempt), ephys-pt-it-onset-not-detected (on DANDI:000010) | Tested |
| fig6-lick-direction-deep-alm-first, fig8d-alm-deep-lick-direction | Blocked: no plane depth |
| fig8d-mm-object-location, fig8ef-area-explains-variance | Blocked: no plane depth (the parts that need no depth are shown, not deciding) |
| widefield-mm-to-alm, widefield-contralateral-bilateral | Blocked: no wide-field imaging |
| figs9-photoinhibition-alm-mm | Blocked: no photoinhibition data |
| fig1c-whisking-similar | Blocked: no whisker tracking |
| ephys-selective-counts, fig7d-ephys-area-distribution | Blocked: the silicon probe data are not public (placeholder DOI, not on DANDI) |
The conventions T, M, A, D, B and O cited below are those of the claims file's opening comment; the
helper module chen2017.py says where each parameter comes from.
import collections
import re
import time
import warnings
import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import chen2017 as C
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
assets = C.list_assets()
headers = C.read_headers(assets)
print(f"read {len(headers)} files of dandiset {C.DANDISET_ID} version {C.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
# Copies of one plane share subject, session, frame and ROI counts, trial table and the first and
# last 200 frames of both trace series; each plane is counted once, as its first path.
groups = collections.defaultdict(list)
for h in headers:
groups[C.content_key(h)].append(h)
planes = sorted((sorted(v, key=lambda h: h["path"])[0] for v in groups.values()), key=lambda h: h["path"])
results = C.process_planes(planes, shuffle_seed=0, start_events=("sample", "start_time"))
print(f"processed {len(planes)} planes in {time.time() - started:.0f} s")
# One row per ROI.
cells = pd.DataFrame({
"plane": np.concatenate([[i] * h["n_rois"] for i, h in enumerate(planes)]),
"animal": np.concatenate([[h["subject"]] * h["n_rois"] for h in planes]),
"area": np.concatenate([[h["area"]] * h["n_rois"] for h in planes]),
"cell_type": np.concatenate([h["cell_type"] for h in planes]),
"included": np.concatenate([h["included"] for h in planes]),
"f0_min": np.concatenate([r["T1"]["f0_min"] for r in results]),
"f0_nonpositive": np.concatenate([r["T1"]["f0_nonpositive"] for r in results]),
"dff_abs_max": np.concatenate([r["T1"]["dff_abs_max"] for r in results]),
})
area, ctype = cells.area.values, cells.cell_type.values
animal, plane = cells.animal.values, cells.plane.values
ALM, MM = area == "ALM", area == "MM"
A = {b: C.assemble(planes, results, b) for b in ("T1", "T2")}
read 210 files of dandiset 000015 version 0.220126.1914 in 316 s
processed 105 planes in 488 s
Data check¶
Each check of the metadata against the values is stated under its output.
info = pd.DataFrame([{k: h[k] for k in ("path", "identifier", "subject", "session", "area", "location",
"rate", "n_frames", "n_rois", "size")} for h in headers])
info["copies"] = [len(groups[C.content_key(h)]) for h in headers]
print("Files:", len(headers), " distinct planes:", len(planes),
" copies per plane:", dict(collections.Counter(len(v) for v in groups.values())))
print("Distinct SHA-256 digests over the 210 files:", len({h['sha256'] for h in headers}),
"; distinct identifiers:", info.identifier.nunique())
pairs = [sorted(h["identifier"] for h in v) for v in groups.values()]
same_size = all(len({h["size"] for h in v}) == 1 for v in groups.values())
print("Every pair has the same file size:", same_size)
print("First five pairs (identifiers):", pairs[:5])
print("\nPlanes by subject and area (the area is from experiment_description; ImagingPlane.location):")
pl = pd.DataFrame([{k: h[k] for k in ("subject", "session", "area", "location", "rate", "n_rois")} for h in planes])
print(pd.crosstab(pl.subject, pl.area, margins=True).to_string())
print("\nImagingPlane.location values:", dict(collections.Counter(pl.location)))
print("Imaging rates (Hz), planes:", dict(collections.Counter(pl.rate.round(4))))
print("Rate in the trace series equals the plane's imaging_rate:", all(abs(h["rate"] - h["rate_attr"]) < 1e-9 for h in planes),
"; starting_time values:", sorted({h["starting_time"] for h in planes}))
print("Units of the trace series:", sorted({str(h["rois_unit"]) for h in planes}))
print("Sessions (subject and date):", pl.groupby(["subject", "session"]).ngroups)
# Three pairs (the smallest, the median and the largest file) are read in full and every dataset
# compared.
by_size = sorted(groups.values(), key=lambda v: v[0]["size"])
for v in (by_size[0], by_size[len(by_size) // 2], by_size[-1]):
a, b = sorted(v, key=lambda h: h["path"])
n_ds, differ = C.compare_files(a, b)
print(f"Full comparison {a['path']} / {b['path']} ({a['size'] / 1e6:.0f} MB): {n_ds} datasets, differing: {differ}")
Files: 210 distinct planes: 105 copies per plane: {2: 105}
Distinct SHA-256 digests over the 210 files: 210 ; distinct identifiers: 210
Every pair has the same file size: True
First five pairs (identifiers): [['an041_2014-08-21_31', 'an041_2014-08-21_9'], ['an041_2014-08-23_20', 'an041_2014-08-23_42'], ['an041_2014-08-23_2', 'an041_2014-08-23_24'], ['an041_2014-08-23_21', 'an041_2014-08-23_43'], ['an041_2014-08-25_14', 'an041_2014-08-25_36']]
Planes by subject and area (the area is from experiment_description; ImagingPlane.location):
area ALM MM All
subject
an041 11 11 22
an043 0 9 9
an044 15 12 27
an063 16 4 20
an064 9 11 20
an066 7 0 7
All 58 47 105
ImagingPlane.location values: {'ALM': 105}
Imaging rates (Hz), planes: {13.9529: 81, 14.2038: 24}
Rate in the trace series equals the plane's imaging_rate: True ; starting_time values: [0.0]
Units of the trace series: ['a.u.']
Sessions (subject and date): 34
Full comparison sub-an044/sub-an044_ses-20140923_obj-1o335up.nwb / sub-an044/sub-an044_ses-20140923_obj-7yvntr.nwb (17 MB): 46 datasets, differing: ['file_create_date', 'identifier']
Full comparison sub-an064/sub-an064_ses-20150211_obj-128cdly.nwb / sub-an064/sub-an064_ses-20150211_obj-13oug0a.nwb (78 MB): 46 datasets, differing: ['file_create_date', 'identifier']
Full comparison sub-an041/sub-an041_ses-20140825_obj-1k6balg.nwb / sub-an041/sub-an041_ses-20140825_obj-w5c8hn.nwb (182 MB): 46 datasets, differing: ['file_create_date', 'identifier']
Every plane is deposited twice: 105 pairs of files with the same size, the same trial table and the
same first and last 200 frames of both trace series, but different identifiers (for example
an041_2014-08-23_21 and an041_2014-08-23_43) and different digests. Three pairs read in full
(above) have every dataset identical except identifier and file_create_date; the other pairs
are compared only on the parts listed. Each plane is counted once. ImagingPlane.location is "ALM" in every file,
also in the 47 planes whose experiment_description says "left_vm1"; the area used here (convention
A, reading A1) is from experiment_description, and its counts are checked against the paper's
below. The frame rates (13.95 and 14.20 Hz) agree with the stated 14 Hz. The traces are in "a.u.".
tt_all = np.concatenate([C.trial_types(h) for h in planes])
early = np.concatenate([h["early_lick"] for h in planes])
outc = np.concatenate([h["outcome"] for h in planes])
instr = np.concatenate([h["instruction"] for h in planes])
print("Trials in all planes:", len(outc))
print("outcome x early_lick:")
print(pd.crosstab(outc, early).to_string())
print("\ntrial_instruction x early_lick:")
print(pd.crosstab(instr, early).to_string())
print("\nTrial types (convention T):", dict(collections.Counter(tt_all)))
print("task values:", sorted({t for h in planes for t in h["task"]}))
rows = []
for h in planes:
ok = C.trial_types(h) != ""
rows.append({"go_minus_sample": h["go"][ok] - h["sample"][ok], "delay_minus_sample": h["delay"][ok] - h["sample"][ok],
"sample_minus_start": h["sample"] - h["trial_start"],
"frames_s_minus_last_stop": h["n_frames"] / h["rate"] - h["trial_stop"][-1]})
gs = np.concatenate([r["go_minus_sample"] for r in rows])
ds = np.concatenate([r["delay_minus_sample"] for r in rows])
ss = np.concatenate([r["sample_minus_start"] for r in rows])
fl = np.array([r["frames_s_minus_last_stop"] for r in rows])
print(f"\nAnalysed trials: go - sample {gs.min():.3f} to {gs.max():.3f} s (median {np.median(gs):.3f}); "
f"delay - sample {ds.min():.3f} to {ds.max():.3f} s")
print(f"All trials: sample - start_time {ss.min():.3f} to {ss.max():.3f} s")
print(f"Recording length (frames / rate) minus the last trial's stop_time: {fl.min():.3f} to {fl.max():.3f} s")
print("Analysed trials dropped because the window leaves the recording:", sum(r["T1"]["dropped_edge"] for r in results))
Trials in all planes: 11274
outcome x early_lick:
col_0 False True
row_0
hit 7305 0
ignore 45 0
miss 2371 0
non-performing 0 1553
trial_instruction x early_lick:
col_0 False True
row_0
left 4842 0
non-performing 0 1553
right 4879 0
Trial types (convention T): {'CL': 3729, '': 1598, 'ER': 1279, 'CR': 3576, 'EL': 1092}
task values: [np.str_('audio delay')]
Analysed trials: go - sample 4.161 to 4.202 s (median 4.200); delay - sample 1.200 to 1.200 s
All trials: sample - start_time 0.622 to 0.623 s
Recording length (frames / rate) minus the last trial's stop_time: 0.070 to 0.072 s
Analysed trials dropped because the window leaves the recording: 78
The trial table holds 'hit', 'miss', 'ignore' and 'non-performing' outcomes. Every early-lick trial
has outcome and instruction 'non-performing', so its pole position is not recorded; 'ignore' trials
(no lick after the go cue) are left out as convention T says. On every analysed trial the go cue
comes 4.2 s after the sample event and the delay starts 1.2 s after it, as the paper's sample (1.2 s)
and delay (3 s) epochs require, while the trial table's start_time is 0.623 s before sample.
So the paper's trial start, at which the claims file puts the go cue 4.2 s earlier, is the sample
event, and it is used as trial start throughout. The traces cover the trial table to within a
frame. A few analysed trials whose window runs past the end of the recording are dropped (count
above).
print("ROIs:", len(cells), " by area:", dict(collections.Counter(area)))
print("cell_type values by area:")
print(pd.crosstab(cells.area, cells.cell_type, margins=True).to_string())
print("\ncell_type by animal:")
print(pd.crosstab(cells.animal, cells.cell_type).to_string())
print("\n'included' (described as 'whether to include this roi into later analyses'):",
dict(collections.Counter(cells.included)), " by area:", cells.groupby("area").included.sum().to_dict())
with C._open(planes[0]["asset_id"]) as fh:
names = []
fh.visit(lambda n: names.append(n) if not n.startswith("specifications") else None)
ip_fields = sorted(fh["general/optophysiology/ImagingPlane"].keys())
print(f"\nEvery object in {planes[0]['path']} (except the NWB specification):")
print(names)
print("ImagingPlane fields:", ip_fields)
print("\nBehavioralEvents series, over all planes:", dict(collections.Counter(tuple(h["behavior_events"]) for h in planes)))
# Every text value (string datasets and attributes) of every plane, with the number of planes
# holding each distinct value; identifiers and dates, which differ per file, are counted only.
per_file = {"identifier", "file_create_date", "session_start_time", "timestamps_reference_time",
"general/subject/subject_id", "general/experiment_description"}
texts = collections.defaultdict(collections.Counter)
for h in planes:
for k, v in h["texts"].items():
texts[k][v] += 1
print("\nText values over the 105 planes (path: value [planes]):")
for k in sorted(texts):
if k in per_file:
print(f" {k}: {len(texts[k])} distinct values")
else:
print(f" {k}: " + "; ".join(f"{v!r} [{n}]" for v, n in texts[k].most_common()))
depth_words = re.compile(r"depth|micron|\bum\b|µm|\bz\b|coordinat|stereotax|bregma|lateral|anterior|layer|location", re.I)
hits = sorted({(k, v) for h in planes for k, v in h["texts"].items() if depth_words.search(v)})
print("\nText values naming depth, microns, z, coordinates, layer or location:", hits)
ROIs: 15431 by area: {'MM': 6594, 'ALM': 8837}
cell_type values by area:
cell_type IT PT unknown All
area
ALM 1873 385 6579 8837
MM 758 334 5502 6594
All 2631 719 12081 15431
cell_type by animal:
cell_type IT PT unknown
animal
an041 0 0 3515
an043 0 0 940
an044 0 0 2880
an063 1004 160 1813
an064 1235 527 2166
an066 392 32 767
'included' (described as 'whether to include this roi into later analyses'): {True: 8655, False: 6776} by area: {'ALM': 4694, 'MM': 3961}
Every object in sub-an041/sub-an041_ses-20140821_obj-17pzgym.nwb (except the NWB specification):
['acquisition', 'acquisition/images', 'acquisition/images/GCaMP at 940nm', 'analysis', 'file_create_date', 'general', 'general/devices', 'general/devices/two-photon microscope with Thorlabs resonant galvo scannner', 'general/experiment_description', 'general/experimenter', 'general/institution', 'general/keywords', 'general/optophysiology', 'general/optophysiology/ImagingPlane', 'general/optophysiology/ImagingPlane/description', 'general/optophysiology/ImagingPlane/excitation_lambda', 'general/optophysiology/ImagingPlane/green', 'general/optophysiology/ImagingPlane/green/description', 'general/optophysiology/ImagingPlane/green/emission_lambda', 'general/optophysiology/ImagingPlane/imaging_rate', 'general/optophysiology/ImagingPlane/indicator', 'general/optophysiology/ImagingPlane/location', 'general/related_publications', 'general/subject', 'general/subject/genotype', 'general/subject/sex', 'general/subject/species', 'general/subject/subject_id', 'general/virus', 'identifier', 'intervals', 'intervals/trials', 'intervals/trials/early_lick', 'intervals/trials/id', 'intervals/trials/outcome', 'intervals/trials/start_time', 'intervals/trials/stop_time', 'intervals/trials/task', 'intervals/trials/task_protocol', 'intervals/trials/trial_instruction', 'processing', 'processing/behavior', 'processing/behavior/BehavioralEvents', 'processing/behavior/BehavioralEvents/delay', 'processing/behavior/BehavioralEvents/delay/data', 'processing/behavior/BehavioralEvents/delay/timestamps', 'processing/behavior/BehavioralEvents/go', 'processing/behavior/BehavioralEvents/go/data', 'processing/behavior/BehavioralEvents/go/timestamps', 'processing/behavior/BehavioralEvents/sample', 'processing/behavior/BehavioralEvents/sample/data', 'processing/behavior/BehavioralEvents/sample/timestamps', 'processing/ophys', 'processing/ophys/Fluorescence', 'processing/ophys/Fluorescence/NeuropilResponseSeries', 'processing/ophys/Fluorescence/NeuropilResponseSeries/data', 'processing/ophys/Fluorescence/NeuropilResponseSeries/starting_time', 'processing/ophys/Fluorescence/RoiResponseSeries', 'processing/ophys/Fluorescence/RoiResponseSeries/data', 'processing/ophys/Fluorescence/RoiResponseSeries/rois', 'processing/ophys/Fluorescence/RoiResponseSeries/starting_time', 'processing/ophys/ImageSegmentation', 'processing/ophys/ImageSegmentation/PlaneSegmentation', 'processing/ophys/ImageSegmentation/PlaneSegmentation/cell_type', 'processing/ophys/ImageSegmentation/PlaneSegmentation/id', 'processing/ophys/ImageSegmentation/PlaneSegmentation/image_mask', 'processing/ophys/ImageSegmentation/PlaneSegmentation/included', 'processing/ophys/ImageSegmentation/PlaneSegmentation/reference_images', 'session_description', 'session_start_time', 'stimulus', 'stimulus/presentation', 'stimulus/templates', 'timestamps_reference_time']
ImagingPlane fields: ['description', 'device', 'excitation_lambda', 'green', 'imaging_rate', 'indicator', 'location']
BehavioralEvents series, over all planes: {('delay', 'go', 'sample'): 105}
Text values over the 105 planes (path: value [planes]):
acquisition/images@description: 'Structural images of a scanning' [105]
file_create_date: 105 distinct values
general/experiment_description: 2 distinct values
general/experimenter: 'Tsai-Wen Chen' [105]
general/institution: 'Janelia Research Campus' [105]
general/keywords: 'motor planning | anterior lateral cortex | medial motor cortex | ALM | MM | Two-photon imaging' [105]
general/optophysiology/ImagingPlane/description: 'Imaging session for PT and IT neurons during audio delay task' [105]
general/optophysiology/ImagingPlane/green/description: 'green channel' [105]
general/optophysiology/ImagingPlane/indicator: 'GCaMP6s' [105]
general/optophysiology/ImagingPlane/location: 'ALM' [105]
general/related_publications: 'http://dx.doi.org/10.1016/j.neuron.2017.05.005' [105]
general/subject/genotype: '' [105]
general/subject/sex: 'M' [105]
general/subject/species: 'Mus musculus' [105]
general/subject/subject_id: 6 distinct values
general/virus: '[]' [105]
identifier: 105 distinct values
intervals/trials/early_lick@description: 'varchar(32) #' [105]
intervals/trials/outcome@description: 'varchar(32) #' [105]
intervals/trials/start_time@description: 'Start time of epoch, in seconds' [105]
intervals/trials/stop_time@description: 'Stop time of epoch, in seconds' [105]
intervals/trials/task@description: 'varchar(12) # task type' [105]
intervals/trials/task_protocol@description: 'tinyint # task protocol' [105]
intervals/trials/trial_instruction@description: 'varchar(16) #' [105]
intervals/trials@colnames: 'start_time | stop_time | task | task_protocol | trial_instruction | early_lick | outcome' [105]
intervals/trials@description: 'experimental trials' [105]
processing/behavior/BehavioralEvents/delay/data@unit: 'a.u.' [105]
processing/behavior/BehavioralEvents/delay/timestamps@unit: 'seconds' [105]
processing/behavior/BehavioralEvents/delay@comments: 'no comments' [105]
processing/behavior/BehavioralEvents/delay@description: 'Time stamps of the beginning of the delay period on each trial.' [105]
processing/behavior/BehavioralEvents/go/data@unit: 'a.u.' [105]
processing/behavior/BehavioralEvents/go/timestamps@unit: 'seconds' [105]
processing/behavior/BehavioralEvents/go@comments: 'no comments' [105]
processing/behavior/BehavioralEvents/go@description: 'Time stamps of the go cue signal on each trial.' [105]
processing/behavior/BehavioralEvents/sample/data@unit: 'a.u.' [105]
processing/behavior/BehavioralEvents/sample/timestamps@unit: 'seconds' [105]
processing/behavior/BehavioralEvents/sample@comments: 'no comments' [105]
processing/behavior/BehavioralEvents/sample@description: 'Time stamps of the beginning of the sampling period on each trial.' [105]
processing/behavior@description: 'Time of behavioral events in this session' [105]
processing/ophys/Fluorescence/NeuropilResponseSeries/data@unit: 'a.u.' [105]
processing/ophys/Fluorescence/NeuropilResponseSeries/starting_time@unit: 'seconds' [105]
processing/ophys/Fluorescence/NeuropilResponseSeries@comments: 'no comments' [105]
processing/ophys/Fluorescence/NeuropilResponseSeries@description: 'average fluorescence of the neuropil surrounding the roi' [105]
processing/ophys/Fluorescence/RoiResponseSeries/data@unit: 'a.u.' [105]
processing/ophys/Fluorescence/RoiResponseSeries/rois@description: 'table region for rois in this segmentation' [105]
processing/ophys/Fluorescence/RoiResponseSeries/starting_time@unit: 'seconds' [105]
processing/ophys/Fluorescence/RoiResponseSeries@comments: 'no comments' [105]
processing/ophys/Fluorescence/RoiResponseSeries@description: 'average fluorescence of roi' [105]
processing/ophys/ImageSegmentation/PlaneSegmentation/cell_type@description: 'PT, IT, or unknown' [105]
processing/ophys/ImageSegmentation/PlaneSegmentation/image_mask@description: 'Image masks for each ROI' [105]
processing/ophys/ImageSegmentation/PlaneSegmentation/included@description: 'whether to include this roi into later analyses' [105]
processing/ophys/ImageSegmentation/PlaneSegmentation@colnames: 'cell_type | included | image_mask' [105]
processing/ophys/ImageSegmentation/PlaneSegmentation@description: 'output from segmenting the current imaging plane' [105]
processing/ophys@description: 'Processing result of imaging' [105]
session_description: 'Imaging session' [105]
session_start_time: 32 distinct values
timestamps_reference_time: 32 distinct values
Text values naming depth, microns, z, coordinates, layer or location: [('general/keywords', 'motor planning | anterior lateral cortex | medial motor cortex | ALM | MM | Two-photon imaging')]
The ROI counts are those of the paper (checked under neuron-counts). PT and IT labels are present in
three mice (an063, an064, an066), as the paper says for its retrograde labelling. The included
column marks 8,655 of 15,431 ROIs, but the paper's counts use all ROIs and nothing says what
excluded the others; all ROIs are used here, and the included ROIs are shown as a variant where the
counts matter. The files hold no plane depth: the ImagingPlane has no depth, origin or coordinates
field, and no text value of any plane names a depth, a coordinate or a layer (the search above; the
only matches are the area names). The behavioural events are sample, delay and go in every plane,
so there are no lick times. There are no whisker data and no other recordings, so the claims that
need these are blocked.
roi_trials = sum(h["n_rois"] * r["T1"]["n_trials"] for h, r in zip(planes, results))
print(f"ROI x trial baselines F0 <= 0: {int(cells.f0_nonpositive.sum())} of {roi_trials} "
f"({100 * cells.f0_nonpositive.sum() / roi_trials:.2f}%)")
print("ROIs with F0 <= 0 on at least one analysed trial:", int((cells.f0_min <= 0).sum()),
" by area:", cells[cells.f0_min <= 0].groupby("area").size().to_dict())
print("Largest |dF/F| per ROI over trials and bins, percentiles 50, 90, 99, 99.9:",
np.percentile(cells.dff_abs_max, [50, 90, 99, 99.9]).round(1))
print(f"Median over ROIs of the per-ROI median F0: {np.median(np.concatenate([r['T1']['f0_median'] for r in results])):.2f} a.u.")
ROI x trial baselines F0 <= 0: 3558 of 1369466 (0.26%)
ROIs with F0 <= 0 on at least one analysed trial: 1262 by area: {'ALM': 522, 'MM': 740}
Largest |dF/F| per ROI over trials and bins, percentiles 50, 90, 99, 99.9: [ 9.2 41.4 765.8 8828.8]
Median over ROIs of the per-ROI median F0: 4.79 a.u.
After neuropil subtraction (r = 0.7) the baseline F0 is small, and for 0.26% of ROI-trials it is zero or negative, in 1,262 ROIs. dF/F = (F - F0)/F0 then reaches thousands. The paper does not say how such trials are treated (question answered): the claims file's reading V1, dF/F as computed by the Methods' formula, decides, and V2, the ROIs with F0 <= 0 on any trial left out, is reported as a check in the claims that cite it.
Clock check. The ROI traces and the trial table share one clock if trial-type selectivity
appears only after the task events that cause it. Below, per bin, the fraction of ROIs with a
significant object-location (alpha) or lick-direction (beta) factor, with the true trial types,
with the trial types permuted within each plane (the null), and with the trials aligned to the
trial table's start_time in place of sample.
g = A["T1"]["times"]
E = A["T1"]["epochs"]
sig_true = A["T1"]["sig"]
sig_shuf = A["T1"]["shuffled"]["sig"]
p_start = np.concatenate([r["T1_start_time"]["p"] for r in results])
sig_start = p_start < C.ALPHA
fig, axes = plt.subplots(1, 2, figsize=(11, 3.2), sharey=True)
for k, name in enumerate(["object location (alpha)", "lick direction (beta)"]):
ax = axes[k]
for s, lab in ((sig_true, "trial types, aligned to sample"), (sig_shuf, "permuted trial types"),
(sig_start, "aligned to start_time")):
ax.plot(g, s[..., k].mean(0), label=lab)
for x in (0, 1.2, 4.2):
ax.axvline(x, color="k", lw=0.5)
ax.set_title(name); ax.set_xlabel("time from trial start (s)")
axes[0].set_ylabel("fraction of ROIs, p < 0.01"); axes[0].legend(fontsize=7)
plt.tight_layout(); plt.show()
for k, name in enumerate(["alpha", "beta"]):
ft, fs, fa = sig_true[..., k].mean(0), sig_shuf[..., k].mean(0), sig_start[..., k].mean(0)
print(f"{name}: true, before trial start {ft[E['pre']].mean():.4f}; max {ft.max():.4f} at {g[ft.argmax()]:.2f} s | "
f"permuted, range {fs.min():.4f} to {fs.max():.4f} | aligned to start_time, max {fa.max():.4f} at {g[fa.argmax()]:.2f} s")
first_a = g[np.flatnonzero(sig_true[..., 0].mean(0) > 0.02)[0]]
first_a_start = g[np.flatnonzero(sig_start[..., 0].mean(0) > 0.02)[0]]
print(f"First bin where alpha-significant ROIs exceed 2%: {first_a:.2f} s aligned to sample, "
f"{first_a_start:.2f} s aligned to start_time (difference {first_a_start - first_a:.2f} s)")
alpha: true, before trial start 0.0106; max 0.0426 at 1.58 s | permuted, range 0.0085 to 0.0127 | aligned to start_time, max 0.0418 at 2.22 s beta: true, before trial start 0.0105; max 0.0710 at 4.87 s | permuted, range 0.0085 to 0.0137 | aligned to start_time, max 0.0684 at 5.52 s First bin where alpha-significant ROIs exceed 2%: 0.65 s aligned to sample, 1.29 s aligned to start_time (difference 0.65 s)
Before trial start the fraction of ROIs with a significant factor is at the 1% the test allows, and
with permuted trial types it stays near 1% in every bin. With the true trial types, object-location
selectivity rises after the pole starts to move and lick-direction selectivity after the go cue;
aligned to start_time the same rise comes later on the time axis, as the offset between the
two events predicts (0.65 s for the first bin above 2%, against the 0.623 s offset, within a
frame). The traces and the trial table are consistent with one clock; the check shows
the alignment to within a few frames, not exactly.
The CRCNS deposit the paper names. The claims on the silicon probe recordings, the wide-field imaging and the photoinhibition would be tested on it if it held them.
# The deposit the paper names for its imaging and electrophysiology (Key Resources Table), read from
# its DOI record at DataCite.
rec = C.requests.get("https://api.datacite.org/dois/10.6080/K0S180F9", timeout=120).json()["data"]["attributes"]
print("DOI 10.6080/K0S180F9:", rec["publisher"], rec["publicationYear"])
print("titles:", [" ".join(t["title"].split()) for t in rec["titles"]])
print("descriptions:", [d["description"] for d in rec["descriptions"]], "; sizes:", rec["sizes"], "; formats:", rec["formats"])
print("URL:", rec.get("url"))
# Any other deposit of the paper's data on DANDI, searched by the paper's DOI and by its authors.
for q in ("10.1016/j.neuron.2017.05.005", "Chen Svoboda 2017"):
found = C.requests.get("https://api.dandiarchive.org/api/dandisets/", params={"search": q, "page_size": 100},
timeout=120).json()["results"]
print(f"DANDI search {q!r}:", [(d["identifier"], (d["most_recent_published_version"] or d["draft_version"])["name"])
for d in found])
DOI 10.6080/K0S180F9: CRCNS.org 2017 titles: ['A map of anticipatory activity in mouse motor cortex, measured with calcium imaging and silicon probes, as described in Chen et al 2017', 'CRCNS.org ssc-8TBD'] descriptions: ['(This data is not yet available but is expected to be sometime during 2017).'] ; sizes: ['TBD'] ; formats: ['HDF5', 'NWB'] URL: http://ezid.cdlib.org/id/doi:10.6080/K0S180F9
DANDI search '10.1016/j.neuron.2017.05.005': [('000015', 'A Map of Anticipatory Activity in Mouse Motor Cortex')]
DANDI search 'Chen Svoboda 2017': [('000015', 'A Map of Anticipatory Activity in Mouse Motor Cortex')]
The DOI record is a placeholder: its alternative title is "CRCNS.org ssc-8TBD", its description says the data are "not yet available", and its size is "TBD". The silicon probe recordings, wide-field imaging, photoinhibition and whisker data are therefore not public there. A search of DANDI by the paper's DOI and by its authors finds only this dandiset. Other archives were not searched here.
early-lick-fraction¶
Mice licked early, during the delay epoch, on 13.3% +- 10% of trials, and these trials were not analysed. (Expected 13.3, tolerance 10, percent.)
rows = []
for h in planes:
n = len(h["outcome"])
rows.append({"subject": h["subject"], "session": h["session"], "n": n, "early": int(h["early_lick"].sum()),
"with_lick": n - int((h["outcome"] == "ignore").sum())})
el = pd.DataFrame(rows)
out = []
for den in ("n", "with_lick"):
per_plane = 100 * el.early / el[den]
ses = el.groupby(["subject", "session"])[["early", den]].sum()
mouse = el.groupby("subject")[["early", den]].sum()
out += [
{"reading": "(a) mean over planes", "denominator": den, "mean": per_plane.mean(), "sd": per_plane.std()},
{"reading": "(b) mean over sessions", "denominator": den, "mean": (100 * ses.early / ses[den]).mean(), "sd": (100 * ses.early / ses[den]).std()},
{"reading": "(c) mean over mice", "denominator": den, "mean": (100 * mouse.early / mouse[den]).mean(), "sd": (100 * mouse.early / mouse[den]).std()},
{"reading": "(d) pooled", "denominator": den, "mean": 100 * el.early.sum() / el[den].sum(), "sd": np.nan}]
res = pd.DataFrame(out)
res["within_13.3_pm_10"] = (res["mean"] - 13.3).abs() <= 10
print(res.round(2).to_string(index=False))
print("\nPer mouse (percent of all trials):", (100 * el.groupby("subject").early.sum() / el.groupby("subject").n.sum()).round(1).to_dict())
reading denominator mean sd within_13.3_pm_10
(a) mean over planes n 13.24 9.14 True
(b) mean over sessions n 14.67 8.27 True
(c) mean over mice n 15.03 5.10 True
(d) pooled n 13.78 NaN True
(a) mean over planes with_lick 13.31 9.20 True
(b) mean over sessions with_lick 14.74 8.31 True
(c) mean over mice with_lick 15.08 5.11 True
(d) pooled with_lick 13.83 NaN True
Per mouse (percent of all trials): {'an041': 8.8, 'an043': 21.9, 'an044': 14.7, 'an063': 14.3, 'an064': 10.6, 'an066': 19.9}
Result: reproduced. Early-lick trials are 13.24% +- 9.14% of trials, averaged over planes with all trials as denominator (reading a), within 13.3 +- 10; every other reading is within tolerance too (table). They are excluded from the analysis.
imaging-sample¶
Two-photon imaging covered 12 locations in 6 mice with 8.8 +- 2.6 planes per location, 105 planes in all, each imaged for 92 +- 26 trials at 73% +- 9% correct, 3 +- 1 planes per session.
pp = pd.DataFrame([{"subject": h["subject"], "session": h["session"], "area": h["area"],
"n_all": len(h["outcome"]), "n_four": int((C.trial_types(h) != "").sum()),
"hits_four": int(np.isin(C.trial_types(h), ["CR", "CL"]).sum()),
"hits_all": int((h["outcome"] == "hit").sum())} for h in planes])
print("mice:", pp.subject.nunique(), " planes:", len(pp))
print("subject x area combinations (no location identifier is in the files):", pp.groupby(["subject", "area"]).ngroups)
print(pp.groupby(["subject", "area"]).size().to_string())
print(f"\ntrials per plane (a) four types: {pp.n_four.mean():.2f} +- {pp.n_four.std():.2f}; (b) all trials: {pp.n_all.mean():.2f} +- {pp.n_all.std():.2f}")
pc_a = 100 * pp.hits_four / pp.n_four
pc_b = 100 * pp.hits_all / pp.n_all
print(f"percent correct per plane (a) of the four types: {pc_a.mean():.2f} +- {pc_a.std():.2f}; (b) of all trials: {pc_b.mean():.2f} +- {pc_b.std():.2f}")
ps = pp.groupby(["subject", "session"]).size()
print(f"planes per session: {ps.mean():.2f} +- {ps.std():.2f}, range {ps.min()} to {ps.max()}, sessions {len(ps)}")
print("\nper mouse: planes, mean trials (a), mean percent correct (a)")
print(pp.assign(pc=pc_a).groupby("subject").agg(planes=("n_four", "size"), trials=("n_four", "mean"), pc=("pc", "mean")).round(1).to_string())
mice: 6 planes: 105
subject x area combinations (no location identifier is in the files): 10
subject area
an041 ALM 11
MM 11
an043 MM 9
an044 ALM 15
MM 12
an063 ALM 16
MM 4
an064 ALM 9
MM 11
an066 ALM 7
trials per plane (a) four types: 92.15 +- 25.98; (b) all trials: 107.37 +- 29.51
percent correct per plane (a) of the four types: 73.83 +- 9.03; (b) of all trials: 63.89 +- 10.97
planes per session: 3.09 +- 1.11, range 1 to 6, sessions 34
per mouse: planes, mean trials (a), mean percent correct (a)
planes trials pc
subject
an041 22 97.0 77.5
an043 9 95.3 77.5
an044 27 101.1 78.6
an063 20 99.6 73.7
an064 20 73.6 63.3
an066 7 70.1 69.7
Result: partially reproduced. Six mice and 105 planes, as stated. With the four analysed trial types (reading a), planes have 92.15 +- 25.98 trials at 73.83% +- 9.03% correct, both within tolerance; with all trials (reading b) there are 107.37 trials per plane, within tolerance, at 63.89% correct, just outside it. Planes per session are 3.09 +- 1.11 (range 1 to 6, 34 sessions). The 12 locations and the 8.8 planes per location cannot be counted: the files carry no location identifier, and subject and area give only 10 combinations.
trials-per-type¶
Each imaged plane had on average 34 +- 13 CR, 36 +- 14 CL, 12 +- 5 ER and 10 +- 4 EL trials.
tc = pd.DataFrame([collections.Counter(C.trial_types(h)) for h in planes]).fillna(0)[list(C.TYPES)]
exp = {"CR": (34, 13), "CL": (36, 14), "ER": (12, 5), "EL": (10, 4)}
for t in C.TYPES:
print(f"{t}: {tc[t].mean():.2f} +- {tc[t].std():.2f} (range {int(tc[t].min())} to {int(tc[t].max())}); "
f"expected {exp[t][0]} +- {exp[t][1]}; within: {abs(tc[t].mean() - exp[t][0]) <= exp[t][1]}")
print("Planes with fewer than 2 trials of some type:", int((tc.min(axis=1) < 2).sum()))
CR: 34.06 +- 13.43 (range 16 to 71); expected 34 +- 13; within: True CL: 35.51 +- 13.94 (range 16 to 81); expected 36 +- 14; within: True ER: 12.18 +- 4.71 (range 4 to 33); expected 12 +- 5; within: True EL: 10.40 +- 4.42 (range 3 to 32); expected 10 +- 4; within: True Planes with fewer than 2 trials of some type: 0
Result: reproduced. Mean trials per plane are CR 34.06 +- 13.43, CL 35.51 +- 13.94, ER 12.18 +- 4.71 and EL 10.40 +- 4.42, each within the paper's spread (34 +- 13, 36 +- 14, 12 +- 5, 10 +- 4), and the spreads match too.
neuron-counts¶
15,431 neurons were imaged across 6 mice, 8,837 in ALM and 6,594 in MM.
n_alm, n_mm = int(ALM.sum()), int(MM.sum())
print("ROIs:", len(cells), " ALM (reading A1, experiment_description left_alm):", n_alm, " MM (left_vm1):", n_mm)
print("Reading A2 (nearer area centre) cannot be run: no ROI or plane coordinates are in the files.")
print("per mouse:", cells.groupby(["animal", "area"]).size().to_dict())
ROIs: 15431 ALM (reading A1, experiment_description left_alm): 8837 MM (left_vm1): 6594
Reading A2 (nearer area centre) cannot be run: no ROI or plane coordinates are in the files.
per mouse: {('an041', 'ALM'): 1941, ('an041', 'MM'): 1574, ('an043', 'MM'): 940, ('an044', 'ALM'): 1375, ('an044', 'MM'): 1505, ('an063', 'ALM'): 2513, ('an063', 'MM'): 464, ('an064', 'ALM'): 1817, ('an064', 'MM'): 2111, ('an066', 'ALM'): 1191}
Result: reproduced. 15,431 ROIs, 8,837 in ALM and 6,594 in MM by the files' area label (reading A1), exactly the paper's numbers. Reading A2 cannot be run; A1 is used by every claim that splits ALM and MM.
fraction-selective-alm-mm¶
A large fraction of neurons distinguished between trial types, 23% in MM (1,523 of 6,594) and 31% in ALM (2,754 of 8,837).
Readings from the criterion: rule R1 (one significant bin) or R2 (two consecutive), window W1 (4.6 s before to 2.5 s after the go cue), W2 (the whole extracted window) or W3 (trial start to 2.5 s after the go cue), bins T1 or T2, area A1.
def selective_neurons(anysel, times, rule, window):
w = {"W1": (times >= -0.4 - 1e-9) & (times <= 6.7 + 1e-9), "W2": np.ones(len(times), bool),
"W3": (times >= 0) & (times <= 6.7 + 1e-9)}[window]
a = anysel[:, w]
return a.any(1) if rule == "R1" else (a[:, 1:] & a[:, :-1]).any(1)
rows = []
for b in ("T1", "T2"):
for rule in ("R1", "R2"):
for win in ("W1", "W2", "W3"):
s = selective_neurons(A[b]["any"], A[b]["times"], rule, win)
r = {"bins": b, "rule": rule, "window": win, "ALM_n": int(s[ALM].sum()), "ALM_%": 100 * s[ALM].mean(),
"MM_n": int(s[MM].sum()), "MM_%": 100 * s[MM].mean()}
r["all_within"] = (abs(r["MM_n"] - 1523) <= 152 and abs(r["ALM_n"] - 2754) <= 275 and
abs(r["MM_%"] - 23) <= 2.3 and abs(r["ALM_%"] - 31) <= 3.1)
r["ALM_higher"] = r["ALM_%"] > r["MM_%"]
v2 = cells.f0_min.values > 0
r["V2_ALM_n"], r["V2_ALM_%"] = int(s[ALM & v2].sum()), 100 * s[ALM & v2].mean()
r["V2_MM_n"], r["V2_MM_%"] = int(s[MM & v2].sum()), 100 * s[MM & v2].mean()
if b == "T1":
sn = selective_neurons(A[b]["shuffled"]["any"], A[b]["times"], rule, win)
r["null_ALM_%"], r["null_MM_%"] = 100 * sn[ALM].mean(), 100 * sn[MM].mean()
rows.append(r)
fs = pd.DataFrame(rows)
print(fs.round(1).to_string(index=False))
print("\nReadings with all four values within tolerance (V1, deciding):", int(fs.all_within.sum()))
print("V2 (ROIs with F0 <= 0 on any trial left out, a check):", int((ALM & v2).sum()), "ALM and", int((MM & v2).sum()), "MM ROIs remain")
print("\nVariants of R2, W1, T1 (not deciding):")
s = selective_neurons(A["T1"]["any"], A["T1"]["times"], "R2", "W1")
for name, m in (("all ROIs", np.ones(len(cells), bool)), ("ROIs with F0 > 0 on every trial", cells.f0_min.values > 0),
("ROIs marked included", cells.included.values)):
print(f" {name}: ALM {int(s[ALM & m].sum())}/{int((ALM & m).sum())} ({100 * s[ALM & m].mean():.1f}%), "
f"MM {int(s[MM & m].sum())}/{int((MM & m).sum())} ({100 * s[MM & m].mean():.1f}%)")
print("\nR2, W1, T1 by animal (percent selective):")
print(pd.DataFrame({"sel": s, "animal": animal, "area": area}).groupby(["animal", "area"]).sel.mean().mul(100).round(1).unstack().to_string())
print("\nPer-bin fraction of ROIs selective for any factor, T1: before trial start "
f"{A['T1']['any'][:, E['pre']].mean():.4f}, max ALM {A['T1']['any'][ALM].mean(0).max():.4f}, max MM {A['T1']['any'][MM].mean(0).max():.4f}")
bins rule window ALM_n ALM_% MM_n MM_% all_within ALM_higher V2_ALM_n V2_ALM_% V2_MM_n V2_MM_% null_ALM_% null_MM_% T1 R1 W1 7687 87.0 5408 82.0 False True 7413 89.2 5066 86.5 77.1 74.4 T1 R1 W2 8069 91.3 5741 87.1 False True 7759 93.3 5359 91.5 84.6 81.5 T1 R1 W3 7596 86.0 5297 80.3 False True 7329 88.1 4975 85.0 75.0 72.1 T1 R2 W1 2877 32.6 1755 26.6 False True 2786 33.5 1675 28.6 11.3 10.9 T1 R2 W2 3243 36.7 1957 29.7 False True 3145 37.8 1865 31.9 14.3 13.8 T1 R2 W3 2853 32.3 1736 26.3 False True 2762 33.2 1656 28.3 10.9 10.5 T2 R1 W1 7543 85.4 5290 80.2 False True 7279 87.5 4968 84.9 NaN NaN T2 R1 W2 7932 89.8 5614 85.1 False True 7641 91.9 5255 89.8 NaN NaN T2 R1 W3 7421 84.0 5169 78.4 False True 7166 86.2 4861 83.0 NaN NaN T2 R2 W1 3981 45.0 2534 38.4 False True 3859 46.4 2420 41.3 NaN NaN T2 R2 W2 4539 51.4 2852 43.3 False True 4399 52.9 2724 46.5 NaN NaN T2 R2 W3 3960 44.8 2516 38.2 False True 3838 46.2 2404 41.1 NaN NaN Readings with all four values within tolerance (V1, deciding): 0 V2 (ROIs with F0 <= 0 on any trial left out, a check): 8315 ALM and 5854 MM ROIs remain Variants of R2, W1, T1 (not deciding): all ROIs: ALM 2877/8837 (32.6%), MM 1755/6594 (26.6%) ROIs with F0 > 0 on every trial: ALM 2786/8315 (33.5%), MM 1675/5854 (28.6%) ROIs marked included: ALM 2045/4694 (43.6%), MM 1363/3961 (34.4%) R2, W1, T1 by animal (percent selective): area ALM MM animal an041 40.1 28.8 an043 NaN 29.6 an044 29.6 32.7 an063 35.2 25.0 an064 25.8 19.7 an066 28.5 NaN Per-bin fraction of ROIs selective for any factor, T1: before trial start 0.0276, max ALM 0.1560, max MM 0.0858
Result: not reproduced. No reading puts all four values within tolerance under V1, the reading that decides. Rule R1 marks 78% to 91% of ROIs as selective, and 72% to 85% with permuted trial types, so it counts chance bins. Rule R2 (permuted: 10.5% to 14.3%) with T1 bins comes closest: ALM 2,853 (32.3%) and MM 1,736 (26.3%) in window W3, and ALM 2,877 (32.6%) and MM 1,755 (26.6%) in window W1; ALM is within tolerance, MM is above it (at most 1,675 and 25.3%). With T2 bins R2 gives 45% and 38% (W1). ALM has the higher fraction under every reading. Under V2 (a check, 8,315 ALM and 5,854 MM ROIs) R2 with T1 bins gives ALM 2,762 (33.2%) and MM 1,656 (28.3%) in W3, and ALM 2,786 (33.5%) and MM 1,675 (28.6%) in W1, still above the tolerance on the MM fraction under every reading.
figs5-single-coefficient¶
The selectivity of most neurons can be attributed to one of the three coefficients, especially early in the trial, so most neurons mainly encode one behavioural variable at a time.
for b in ("T1", "T2"):
t = A[b]["times"]; Eb = A[b]["epochs"]
n_any = A[b]["any"].sum(0)
frac_pure = A[b]["pure"].any(-1).sum(0) / np.maximum(n_any, 1)
w = (t >= 0) & (t < C.GO)
print(f"{b}: bins from trial start to the go cue {w.sum()}, with pure fraction > 0.5: {(frac_pure[w] > 0.5).sum()}; "
f"range {frac_pure[w].min():.3f} to {frac_pure[w].max():.3f}; selective ROIs per bin {n_any[w].min()} to {n_any[w].max()}; "
f"after the go cue mean {frac_pure[Eb['response']].mean():.3f} (min {frac_pure[Eb['response']].min():.3f})")
if b == "T1":
plt.figure(figsize=(5, 2.6)); plt.plot(t, frac_pure); plt.axhline(0.5, color="k", lw=0.5)
plt.xlabel("time from trial start (s)"); plt.ylabel("pure / selective"); plt.title("T1"); plt.show()
sh = A["T1"]["shuffled"]; t = A["T1"]["times"]; w = (t >= 0) & (t < C.GO)
frac_null = sh["pure"].any(-1).sum(0) / np.maximum(sh["any"].sum(0), 1)
print(f"T1 with permuted trial types (null): pure / selective from trial start to the go cue {frac_null[w].min():.3f} to "
f"{frac_null[w].max():.3f} (mean {frac_null[w].mean():.3f}), bins above 0.5: {(frac_null[w] > 0.5).sum()} of {w.sum()}; "
f"selective ROIs per bin {sh['any'].sum(0)[w].min()} to {sh['any'].sum(0)[w].max()}")
T1: bins from trial start to the go cue 59, with pure fraction > 0.5: 59; range 0.845 to 0.928; selective ROIs per bin 413 to 985; after the go cue mean 0.793 (min 0.720)
T2: bins from trial start to the go cue 63, with pure fraction > 0.5: 63; range 0.826 to 0.923; selective ROIs per bin 403 to 1108; after the go cue mean 0.784 (min 0.716) T1 with permuted trial types (null): pure / selective from trial start to the go cue 0.845 to 0.911 (mean 0.876), bins above 0.5: 59 of 59; selective ROIs per bin 368 to 445
Result: reproduced. In every bin from trial start to the go cue, more than half of the selective ROIs are purely selective (lowest 0.845 with T1 bins, 0.826 with T2). After the go cue the fraction is 0.793 (T1) on average. With permuted trial types the same fraction is 0.845 to 0.911 and above one half in all 59 bins, so the criterion is met by chance-level selectivity too: the measure does not tell the result from chance, and the outcome rests on the criterion alone.
fig4a-lick-direction-epochs¶
Different neurons coded for lick direction in different epochs, some during the sample and delay epochs before licking, some during the response epoch, and some after the mouse stopped consuming water.
def fig4a_counts(pure, t):
pb = pure[..., 1]
first = np.where(pb.any(1), np.argmax(pb, 1), -1)
go_bin = np.flatnonzero(t >= C.GO)[0]
return {"delay": int(pb[:, (t >= C.SAMPLE_END) & (t < C.GO)].any(1).sum()),
"first_after_go": int((first >= go_bin).sum()),
"after_go_plus_2s": int(pb[:, t > C.GO + 2].any(1).sum())}
for b in ("T1", "T2"):
print(b, fig4a_counts(A[b]["pure"], A[b]["times"]))
print("T1 with permuted trial types (null):", fig4a_counts(A["T1"]["shuffled"]["pure"], A["T1"]["times"]))
print("Lick times are not in the files, so the median last lick cannot be computed.")
T1 {'delay': 4125, 'first_after_go': 3611, 'after_go_plus_2s': 4426}
T2 {'delay': 3825, 'first_after_go': 3514, 'after_go_plus_2s': 4155}
T1 with permuted trial types (null): {'delay': 3273, 'first_after_go': 2526, 'after_go_plus_2s': 3196}
Lick times are not in the files, so the median last lick cannot be computed.
Result: reproduced. With T1 bins, 4,125 ROIs are purely lick-direction selective in the delay epoch, 3,611 are first so after the go cue, and 4,426 more than 2 s after the go cue (T2: 3,825, 3,514 and 4,155). With permuted trial types the same counts are 3,273, 2,526 and 3,196, so each group exceeds chance, but most of each count is what chance bins alone give; the criterion asks only that such ROIs exist, which the permuted labels also satisfy. The time after the last lick cannot be tested: lick times are not in the files.
fig4b-object-location-memory¶
Many object-location-selective neurons were selective during the sample epoch, and a subset coded object location mainly in the delay epoch after the pole was removed, some up to the end of the delay.
def fig4b_counts(pure, t):
pa = pure[..., 0]
d = pa[:, (t >= C.SAMPLE_END) & (t < C.GO)].sum(1); tot = pa.sum(1)
return {"sample": int(pa[:, (t >= 0) & (t < C.SAMPLE_END)].any(1).sum()),
"most_bins_in_delay": int(((tot > 0) & (d > tot / 2)).sum()),
"last_0.5s_before_go": int(pa[:, (t >= C.GO - 0.5) & (t < C.GO)].any(1).sum())}
for b in ("T1", "T2"):
print(b, fig4b_counts(A[b]["pure"], A[b]["times"]))
print("T1 with permuted trial types (null):", fig4b_counts(A["T1"]["shuffled"]["pure"], A["T1"]["times"]))
T1 {'sample': 2148, 'most_bins_in_delay': 2487, 'last_0.5s_before_go': 1149}
T2 {'sample': 2067, 'most_bins_in_delay': 2483, 'last_0.5s_before_go': 1074}
T1 with permuted trial types (null): {'sample': 1694, 'most_bins_in_delay': 1719, 'last_0.5s_before_go': 707}
Result: reproduced. With T1 bins, 2,148 ROIs are purely object-location selective in the sample epoch, 2,487 have most of their object-location bins in the delay, and 1,149 are selective in the last 0.5 s before the go cue (T2: 2,067, 2,483, 1,074). With permuted trial types the counts are 1,694, 1,719 and 707, so each group exceeds chance; the criterion asks only that such ROIs exist, which the permuted labels also satisfy.
fig4c-outcome-after-lick¶
Outcome-selective neurons typically became selective after the answer lick; some preferred correct and others error trials, and some kept the outcome long after the trial ended.
for b in ("T1", "T2"):
t = A[b]["times"]; pg = A[b]["pure"][..., 2]; K = A[b]["coef"][..., 2]
ever = pg.any(1); first = np.argmax(pg, 1); go_bin = np.flatnonzero(t >= C.GO)[0]
after = ever & (first >= go_bin)
sign_first = np.sign(K[np.arange(len(K)), first])[ever]
print(f"{b}: purely outcome-selective ROIs {ever.sum()}, first bin after the go cue {after.sum()} "
f"({after.sum() / ever.sum():.3f}); gamma sign at the first bin: correct-preferring {(sign_first > 0).sum()}, "
f"error-preferring {(sign_first < 0).sum()}; latest outcome-selective bin {t[np.flatnonzero(pg.any(0))[-1]]:.2f} s")
if b == "T1":
s2 = (pg[:, 1:] & pg[:, :-1])
ever2 = s2.any(1); first2 = np.argmax(s2, 1)
print(f" variant, two consecutive bins (as rule R2): {ever2.sum()} ROIs, first after the go cue "
f"{(ever2 & (first2 >= go_bin)).sum() / ever2.sum():.3f}")
pgs = A["T1"]["shuffled"]["pure"][..., 2]
evs = pgs.any(1)
print(f"T1 permuted (null): {evs.sum()} ROIs, first after the go cue {(evs & (np.argmax(pgs, 1) >= go_bin)).sum() / evs.sum():.3f}")
T1: purely outcome-selective ROIs 9431, first bin after the go cue 3892 (0.413); gamma sign at the first bin: correct-preferring 4157, error-preferring 5274; latest outcome-selective bin 9.17 s variant, two consecutive bins (as rule R2): 1463 ROIs, first after the go cue 0.785 T2: purely outcome-selective ROIs 8892, first bin after the go cue 3715 (0.418); gamma sign at the first bin: correct-preferring 3974, error-preferring 4918; latest outcome-selective bin 9.14 s T1 permuted (null): 8070 ROIs, first after the go cue 0.296
Result: not reproduced. Of the ROIs purely outcome selective in any bin (9,431 with T1 bins), 41.3% have their first such bin after the go cue (T2 41.8%), below the half the criterion requires. Most of these ROIs have only chance bins: with permuted trial types 8,070 ROIs qualify, 29.6% of them first after the go cue. Requiring two consecutive bins (a variant, not deciding) leaves 1,463 ROIs, 78.5% of them first after the go cue. Both signs of gamma occur (4,157 correct-preferring and 5,274 error-preferring at the first bin), and outcome selectivity lasts to the end of the window (9.17 s after trial start).
fig5-mixed-response-epoch¶
Most neurons with mixed selectivity were seen during the response epoch; a subset responded in only one trial type and a few switched the variable they coded during the trial.
for b in ("T1", "T2"):
t = A[b]["times"]; mixed = A[b]["mixed"]; pure = A[b]["pure"]; Eb = A[b]["epochs"]
nb = mixed.sum(0)
early_obj = pure[:, t < C.GO, 0].any(1)
late_other = (pure[:, t >= C.GO, 1] | pure[:, t >= C.GO, 2]).any(1)
print(f"{b}: mixed ROI-bins {nb.sum()}, after the go cue {nb[Eb['response']].sum()} ({nb[Eb['response']].sum() / nb.sum():.3f}); "
f"ROIs purely object-location selective before the go cue and purely lick-direction or outcome selective after it: "
f"{(early_obj & late_other).sum()}")
nbs = A["T1"]["shuffled"]["mixed"].sum(0)
print(f"T1 permuted (null): mixed ROI-bins {nbs.sum()}, after the go cue {nbs[E['response']].sum() / nbs.sum():.3f} "
f"(the response window is {E['response'].mean():.3f} of the bins)")
T1: mixed ROI-bins 23363, after the go cue 17869 (0.765); ROIs purely object-location selective before the go cue and purely lick-direction or outcome selective after it: 4682 T2: mixed ROI-bins 28496, after the go cue 22115 (0.776); ROIs purely object-location selective before the go cue and purely lick-direction or outcome selective after it: 4327 T1 permuted (null): mixed ROI-bins 7298, after the go cue 0.523 (the response window is 0.500 of the bins)
Result: reproduced. 76.5% of mixed-selective ROI-bins fall after the go cue with T1 bins (77.6% with T2), against 50% of the bins in the window and 52.3% with permuted labels. 4,682 ROIs (T1) are purely object-location selective before the go cue and purely lick-direction or outcome selective after it.
fig6-object-location-mm¶
During the sample and early delay epochs, object location was represented mainly in MM.
for b in ("T1", "T2"):
t = A[b]["times"]; pa = A[b]["pure"][..., 0]
fa, fm = pa[ALM].mean(0), pa[MM].mean(0)
for end in (2.2, 2.7):
w = (t >= 0) & (t <= end + 1e-9)
print(f"{b}, trial start to {end} s: MM above ALM in {(fm[w] > fa[w]).sum()} of {w.sum()} bins; "
f"mean MM {fm[w].mean():.4f}, ALM {fa[w].mean():.4f}")
t = A["T1"]["times"]
fig, axes = plt.subplots(1, 2, figsize=(11, 3), sharey=True)
names = ["object location", "lick direction", "outcome"]
for ax, (lab, m) in zip(axes, (("ALM", ALM), ("MM", MM))):
for k in range(3):
ax.plot(t, A["T1"]["pure"][m][..., k].mean(0), label=names[k])
ax.plot(t, A["T1"]["mixed"][m].mean(0), label="mixed")
for x in (0, 1.2, 4.2):
ax.axvline(x, color="k", lw=0.5)
ax.set_title(f"{lab}, T1"); ax.set_xlabel("time from trial start (s)")
axes[0].set_ylabel("fraction of ROIs"); axes[0].legend(fontsize=7); plt.tight_layout(); plt.show()
for lab, m in (("ALM", ALM), ("MM", MM)):
print(lab, "plotted fractions range:", {n: (round(float(A['T1']['pure'][m][..., k].mean(0).min()), 4), round(float(A['T1']['pure'][m][..., k].mean(0).max()), 4)) for k, n in enumerate(names)},
"mixed", (round(float(A['T1']['mixed'][m].mean(0).min()), 4), round(float(A['T1']['mixed'][m].mean(0).max()), 4)))
pa = A["T1"]["pure"][..., 0]; w = (t >= 0) & (t <= 2.2 + 1e-9)
print("\nby animal, mean fraction purely object-location selective, trial start to 2.2 s (T1):")
print(pd.DataFrame({"f": pa[:, w].mean(1), "animal": animal, "area": area}).groupby(["animal", "area"]).f.mean().round(4).unstack().to_string())
T1, trial start to 2.2 s: MM above ALM in 27 of 31 bins; mean MM 0.0368, ALM 0.0121 T1, trial start to 2.7 s: MM above ALM in 34 of 38 bins; mean MM 0.0382, ALM 0.0125 T2, trial start to 2.2 s: MM above ALM in 29 of 33 bins; mean MM 0.0400, ALM 0.0128 T2, trial start to 2.7 s: MM above ALM in 37 of 41 bins; mean MM 0.0417, ALM 0.0133
ALM plotted fractions range: {'object location': (0.0058, 0.0199), 'lick direction': (0.0059, 0.0618), 'outcome': (0.0071, 0.0373)} mixed (0.0024, 0.0459)
MM plotted fractions range: {'object location': (0.0059, 0.0678), 'lick direction': (0.005, 0.0256), 'outcome': (0.0055, 0.0268)} mixed (0.0018, 0.0167)
by animal, mean fraction purely object-location selective, trial start to 2.2 s (T1):
area ALM MM animal an041 0.0147 0.0315 an043 NaN 0.0496 an044 0.0098 0.0474 an063 0.0126 0.0412 an064 0.0098 0.0265 an066 0.0128 NaN
Result: reproduced. The fraction of purely object-location-selective ROIs is higher in MM than in ALM in 27 of 31 bins from trial start to 2.2 s, and 34 of 38 to 2.7 s, with T1 bins (T2: 29 of 33 and 37 of 41).
fig6-lick-direction-spread¶
Lick-direction selectivity, in number of neurons and strength, kept increasing and spread across ALM layers and to a minor extent into MM, peaking after the go cue.
for b in ("T1", "T2"):
t = A[b]["times"]; Eb = A[b]["epochs"]
pb = A[b]["pure"][..., 1]; ab = np.abs(A[b]["coef_sig"][..., 1])
fA, fM, bA = pb[ALM].mean(0), pb[MM].mean(0), ab[ALM].mean(0)
print(f"{b}: (a) ALM fraction peaks at {t[fA.argmax()]:.2f} s ({fA.max():.4f}), ALM mean |beta| peaks at {t[bA.argmax()]:.2f} s "
f"({bA.max():.4f}); both after the go cue: {t[fA.argmax()] >= C.GO and t[bA.argmax()] >= C.GO}")
print(f" (c) MM maximum {fM.max():.4f} at {t[fM.argmax()]:.2f} s; MM in the delay: mean {fM[Eb['delay']].mean():.4f}, "
f"max {fM[Eb['delay']].max():.4f}; ALM maximum {fA.max():.4f}; "
f"holds with the delay mean: {fM.max() > fM[Eb['delay']].mean() and fM.max() < fA.max()}, "
f"with the delay maximum: {fM.max() > fM[Eb['delay']].max() and fM.max() < fA.max()}")
print("(b) superficial ALM cannot be separated: the files hold no plane depth.")
ab = np.abs(A["T1"]["coef_sig"][..., 1]); keep = cells.f0_min.values > 0
print(f"Variant, ROIs with F0 > 0 on every trial (T1): ALM mean |beta| peaks at {t[ab[ALM & keep].mean(0).argmax()]:.2f} s")
print("by animal, time of the ALM peak of the fraction (T1):",
{a: round(float(t[A['T1']['pure'][ALM & (animal == a)][..., 1].mean(0).argmax()]), 2) for a in np.unique(animal[ALM])})
T1: (a) ALM fraction peaks at 4.44 s (0.0618), ALM mean |beta| peaks at 4.66 s (0.1494); both after the go cue: True
(c) MM maximum 0.0256 at 5.38 s; MM in the delay: mean 0.0083, max 0.0152; ALM maximum 0.0618; holds with the delay mean: True, with the delay maximum: True
T2: (a) ALM fraction peaks at 4.67 s (0.0645), ALM mean |beta| peaks at 4.67 s (0.1531); both after the go cue: True
(c) MM maximum 0.0293 at 4.94 s; MM in the delay: mean 0.0082, max 0.0150; ALM maximum 0.0645; holds with the delay mean: True, with the delay maximum: True
(b) superficial ALM cannot be separated: the files hold no plane depth.
Variant, ROIs with F0 > 0 on every trial (T1): ALM mean |beta| peaks at 4.34 s
by animal, time of the ALM peak of the fraction (T1): {'an041': 4.67, 'an044': 4.4, 'an063': 4.14, 'an064': 4.14, 'an066': 4.07}
Result: partially reproduced. In ALM the fraction of lick-direction-selective ROIs and their mean |beta| both peak after the go cue (4.44 and 4.66 s, T1), and in MM the fraction peaks at 0.026, above its delay-epoch level (mean 0.008, maximum 0.015) and below the ALM peak (0.062), under both readings of "its value in the delay epoch" that the criterion now names (the paper does not say which) and both bin readings. The spread across ALM layers (part b) cannot be tested, as the files hold no depth.
fig6-outcome-response-uniform¶
Outcome-selective neurons were seen only during the response epoch, scattered roughly uniformly through ALM and MM; mixed neurons were also seen mainly in the response epoch in both areas.
for b in ("T1", "T2"):
t = A[b]["times"]; resp = (t >= C.GO) & (t <= C.GO + 2.5 + 1e-9)
for name, X in (("outcome", A[b]["pure"][..., 2]), ("mixed", A[b]["mixed"])):
vals = {}
for lab, m in (("ALM", ALM), ("MM", MM)):
f = X[m].mean(0)
vals[lab] = (f[t < C.GO].max(), f[t >= C.GO].max(), f[resp].mean())
a_ok = all(v[0] < v[1] / 3 for v in vals.values())
ratio = max(vals["ALM"][2], vals["MM"][2]) / min(vals["ALM"][2], vals["MM"][2])
print(f"{b} {name}: " + "; ".join(f"{k} max before go {v[0]:.4f}, max after {v[1]:.4f}, response mean {v[2]:.4f}" for k, v in vals.items())
+ f" | (a) {a_ok}; (b) ratio {ratio:.3f}, below 1.5: {ratio < 1.5}")
ps = A["T1"]["shuffled"]["pure"][..., 2]
print(f"T1 permuted (null), outcome fraction per bin: {ps.mean(0).min():.4f} to {ps.mean(0).max():.4f}")
T1 outcome: ALM max before go 0.0123, max after 0.0373, response mean 0.0282; MM max before go 0.0112, max after 0.0268, response mean 0.0193 | (a) False; (b) ratio 1.457, below 1.5: True T1 mixed: ALM max before go 0.0112, max after 0.0459, response mean 0.0313; MM max before go 0.0059, max after 0.0167, response mean 0.0110 | (a) False; (b) ratio 2.846, below 1.5: False T2 outcome: ALM max before go 0.0124, max after 0.0428, response mean 0.0302; MM max before go 0.0115, max after 0.0315, response mean 0.0209 | (a) False; (b) ratio 1.448, below 1.5: True T2 mixed: ALM max before go 0.0131, max after 0.0552, response mean 0.0358; MM max before go 0.0064, max after 0.0202, response mean 0.0127 | (a) True; (b) ratio 2.817, below 1.5: False T1 permuted (null), outcome fraction per bin: 0.0067 to 0.0105
Result: not reproduced. Part (b) holds: averaged over the response epoch the ALM and MM outcome fractions differ by a factor of 1.46 (T1) and 1.45 (T2). Part (a) fails in MM: its largest fraction before the go cue (0.0112, T1) is not below a third of the largest after it (0.0268), while in ALM it just is (0.0123 against 0.0373); T2 gives the same. The pre-cue maxima are near the range of the permuted-label fractions (0.0067 to 0.0105). Mixed ROIs, reported only, are 2.8 times as frequent in ALM as in MM over the response epoch.
fig6h-preparatory-alm-object-mm¶
Throughout the trial preparatory activity, as the fraction of lick-direction-selective neurons or the average coding coefficient, was most pronounced in ALM, while object location selectivity was mainly limited to MM.
for b in ("T1", "T2"):
t = A[b]["times"]; Eb = A[b]["epochs"]
pb, pa = A[b]["pure"][..., 1], A[b]["pure"][..., 0]
ab, aa = np.abs(A[b]["coef_sig"][..., 1]), np.abs(A[b]["coef_sig"][..., 0])
d = Eb["delay"]; w = (t >= 0) & (t < C.GO)
lick = (pb[ALM].mean(0) > pb[MM].mean(0)) & (ab[ALM].mean(0) > ab[MM].mean(0))
obj = (pa[MM].mean(0) > pa[ALM].mean(0)) & (aa[MM].mean(0) > aa[ALM].mean(0))
print(f"{b}: delay bins with ALM above MM on both lick-direction measures {lick[d].sum()} of {d.sum()} "
f"(fraction alone {(pb[ALM].mean(0) > pb[MM].mean(0))[d].sum()}, mean |beta| alone {(ab[ALM].mean(0) > ab[MM].mean(0))[d].sum()}); "
f"bins trial start to go with MM above ALM on both object measures {obj[w].sum()} of {w.sum()}")
ab = np.abs(A["T1"]["coef_sig"][..., 1]); keep = cells.f0_min.values > 0
print("Variant, ROIs with F0 > 0 on every trial (T1): delay bins with ALM mean |beta| above MM:",
int((ab[ALM & keep].mean(0) > ab[MM & keep].mean(0))[E["delay"]].sum()), "of", int(E["delay"].sum()))
T1: delay bins with ALM above MM on both lick-direction measures 24 of 42 (fraction alone 42, mean |beta| alone 24); bins trial start to go with MM above ALM on both object measures 52 of 59 T2: delay bins with ALM above MM on both lick-direction measures 30 of 45 (fraction alone 45, mean |beta| alone 30); bins trial start to go with MM above ALM on both object measures 57 of 63 Variant, ROIs with F0 > 0 on every trial (T1): delay bins with ALM mean |beta| above MM: 27 of 42
Result: reproduced. In the delay epoch ALM is above MM on both lick-direction measures in 24 of 42 bins (T1) and 30 of 45 (T2), more than half; from trial start to the go cue MM is above ALM on both object-location measures in 52 of 59 (T1) and 57 of 63 (T2) bins. The lick-direction part is carried by the fraction; on mean |beta| alone ALM leads in fewer delay bins (printed).
figs4-peak-shift¶
In two-photon imaging, the peak of activity shifted from MM to ALM over the delay.
rows = []
for b in ("T1", "T2"):
t = A[b]["times"]; Eb = A[b]["epochs"]; D = A[b]["mean_dff"]
last = (t >= C.GO - 0.5) & (t < C.GO)
for name, keep in (("V1, all ROIs", np.ones(len(cells), bool)), ("V2, F0 > 0 on every trial", cells.f0_min.values > 0)):
cA, cM = D[ALM & keep].mean(0), D[MM & keep].mean(0)
levels = {"L1 sample mean": cA[Eb["sample"]].mean(), "L2 sample max": cA[Eb["sample"]].max()}
ends = {"E1 last bin": cA[Eb["delay"]][-1], "E2 last 0.5 s": cA[last].mean()}
for ln, lv in levels.items():
for en, ev in ends.items():
rows.append({"bins": b, "dF/F": name, "MM_peak_s": t[cM.argmax()], "ALM_peak_s": t[cA.argmax()],
"MM_earlier": t[cM.argmax()] < t[cA.argmax()], "level": ln, "level_value": lv,
"end": en, "end_value": ev, "end_above_level": ev > lv})
pk = pd.DataFrame(rows)
pk["reproduced"] = pk.MM_earlier & pk.end_above_level
print(pk.round(3).to_string(index=False))
print("\nCombinations reproduced (V1):", int(pk.reproduced[pk["dF/F"].str.startswith("V1")].sum()), "of", int(pk["dF/F"].str.startswith("V1").sum()))
t = A["T1"]["times"]; D = A["T1"]["mean_dff"]; keep = cells.f0_min.values > 0
plt.figure(figsize=(5.5, 3))
for lab, m in (("ALM", ALM), ("MM", MM)):
plt.plot(t, D[m].mean(0), label=f"{lab}, all ROIs")
plt.plot(t, D[m & keep].mean(0), "--", label=f"{lab}, F0 > 0")
for x in (0, 1.2, 4.2):
plt.axvline(x, color="k", lw=0.5)
plt.xlabel("time from trial start (s)"); plt.ylabel("mean dF/F"); plt.legend(fontsize=7); plt.show()
for lab, m in (("ALM", ALM), ("MM", MM)):
print(f"{lab} plotted range: all {D[m].mean(0).min():.3f} to {D[m].mean(0).max():.3f}; F0 > 0 {D[m & keep].mean(0).min():.3f} to {D[m & keep].mean(0).max():.3f}")
print("Variant, peak searched only before the go cue or only in the delay epoch (not deciding):")
for b in ("T1", "T2"):
t = A[b]["times"]; D = A[b]["mean_dff"]; Eb = A[b]["epochs"]
for name, keep in (("V1", np.ones(len(cells), bool)), ("V2", cells.f0_min.values > 0)):
out = {}
for wname, wm in (("before go", t < C.GO), ("delay", Eb["delay"])):
for lab, m in (("ALM", ALM), ("MM", MM)):
out[f"{lab} {wname}"] = round(float(t[wm][D[m & keep].mean(0)[wm].argmax()]), 2)
print(f" {b} {name}: {out}")
bins dF/F MM_peak_s ALM_peak_s MM_earlier level level_value end end_value end_above_level reproduced T1 V1, all ROIs 6.092 4.515 False L1 sample mean 0.112 E1 last bin 0.212 True False T1 V1, all ROIs 6.092 4.515 False L1 sample mean 0.112 E2 last 0.5 s 0.205 True False T1 V1, all ROIs 6.092 4.515 False L2 sample max 0.137 E1 last bin 0.212 True False T1 V1, all ROIs 6.092 4.515 False L2 sample max 0.137 E2 last 0.5 s 0.205 True False T1 V2, F0 > 0 on every trial 4.802 4.515 False L1 sample mean 0.122 E1 last bin 0.214 True False T1 V2, F0 > 0 on every trial 4.802 4.515 False L1 sample mean 0.122 E2 last 0.5 s 0.210 True False T1 V2, F0 > 0 on every trial 4.802 4.515 False L2 sample max 0.141 E1 last bin 0.214 True False T1 V2, F0 > 0 on every trial 4.802 4.515 False L2 sample max 0.141 E2 last 0.5 s 0.210 True False T2 V1, all ROIs 6.136 4.536 False L1 sample mean 0.110 E1 last bin 0.211 True False T2 V1, all ROIs 6.136 4.536 False L1 sample mean 0.110 E2 last 0.5 s 0.204 True False T2 V1, all ROIs 6.136 4.536 False L2 sample max 0.136 E1 last bin 0.211 True False T2 V1, all ROIs 6.136 4.536 False L2 sample max 0.136 E2 last 0.5 s 0.204 True False T2 V2, F0 > 0 on every trial 4.802 4.536 False L1 sample mean 0.120 E1 last bin 0.216 True False T2 V2, F0 > 0 on every trial 4.802 4.536 False L1 sample mean 0.120 E2 last 0.5 s 0.210 True False T2 V2, F0 > 0 on every trial 4.802 4.536 False L2 sample max 0.141 E1 last bin 0.216 True False T2 V2, F0 > 0 on every trial 4.802 4.536 False L2 sample max 0.141 E2 last 0.5 s 0.210 True False Combinations reproduced (V1): 0 of 8
ALM plotted range: all -0.005 to 0.313; F0 > 0 -0.007 to 0.310
MM plotted range: all -1.150 to 0.714; F0 > 0 -0.017 to 0.268
Variant, peak searched only before the go cue or only in the delay epoch (not deciding):
T1 V1: {'ALM before go': 4.16, 'MM before go': 2.08, 'ALM delay': 4.16, 'MM delay': 2.08}
T1 V2: {'ALM before go': 4.09, 'MM before go': 1.51, 'ALM delay': 4.09, 'MM delay': 1.51}
T2 V1: {'ALM before go': 4.14, 'MM before go': 0.8, 'ALM delay': 4.14, 'MM delay': 2.6}
T2 V2: {'ALM before go': 4.14, 'MM before go': 1.53, 'ALM delay': 4.14, 'MM delay': 1.53}
Result: not reproduced. The mean dF/F of MM ROIs peaks after that of ALM (6.09 against 4.52 s with T1 bins, 6.14 against 4.54 s with T2), so the MM mean does not peak earlier, and none of the 8 combinations of bins and readings L1, L2, E1, E2 is reproduced under V1. The other part holds under all of them: the ALM mean at the end of the delay (0.212 in the last bin, 0.205 over the last 0.5 s, T1) is above its sample-epoch mean (0.112) and maximum (0.137). Under V2, without the ROIs with F0 <= 0, which dominate the MM mean, MM peaks at 4.80 s, still after ALM (4.52 s). The peak is taken over the whole extracted window, the bins of convention T. Searching only before the go cue (a variant, not deciding), MM peaks before ALM under every bin and dF/F reading: V1 2.08 against 4.16 s (T1) and 0.80 against 4.14 s (T2), V2 1.51 against 4.09 s (T1) and 1.53 against 4.14 s (T2); within the delay epoch alone V1 gives 2.08 against 4.16 s (T1) and 2.60 against 4.14 s (T2). So the outcome turns on the window of the peak, which the criterion does not name; a question has been asked, and the outcome is kept on the whole window until it is answered.
fig8c-pt-deep-only¶
IT neurons were found in both superficial and deep layers, PT neurons only in deep layers; 719 PT and 2,631 IT neurons were imaged in three mice.
n_pt, n_it = int((ctype == "PT").sum()), int((ctype == "IT").sum())
print(f"PT {n_pt} (expected 719 +- 72), IT {n_it} (expected 2631 +- 263); mice with labels: "
f"{sorted(set(animal[(ctype == 'PT') | (ctype == 'IT')]))}")
print(pd.crosstab([cells.animal, cells.area], cells.cell_type).to_string())
print("Depth of PT and IT neurons cannot be checked: the files hold no plane depth.")
PT 719 (expected 719 +- 72), IT 2631 (expected 2631 +- 263); mice with labels: ['an063', 'an064', 'an066']
cell_type IT PT unknown
animal area
an041 ALM 0 0 1941
MM 0 0 1574
an043 MM 0 0 940
an044 ALM 0 0 1375
MM 0 0 1505
an063 ALM 800 106 1607
MM 204 54 206
an064 ALM 681 247 889
MM 554 280 1277
an066 ALM 392 32 767
Depth of PT and IT neurons cannot be checked: the files hold no plane depth.
Result: partially reproduced. 719 PT and 2,631 IT ROIs in three mice, exactly the paper's counts. That PT neurons are only in deep layers and IT neurons in both cannot be tested: the files hold no plane depth.
fig8d-pt-contra-bias¶
During the delay epoch, ALM PT neurons showed a bias toward contralateral (lick right) selectivity.
rng = np.random.default_rng(81)
pt = ALM & (ctype == "PT")
for b in ("T1", "T2"):
Eb = A[b]["epochs"]
v = A[b]["coef_sig"][..., 1][:, Eb["delay"]].mean(1)
for name, keep in (("all", np.ones(len(cells), bool)), ("variant, F0 > 0", cells.f0_min.values > 0)):
g_ = pt & keep
boot = np.array([v[C.joint_resample(animal, plane, [g_], rng)[0]].mean() for _ in range(1000)])
m = v[g_].mean()
print(f"{b} {name}: ALM PT n={g_.sum()}, mean beta over the delay {m:.4f}, bootstrap SD {boot.std():.4f}, "
f"ratio {m / boot.std():.2f}, above 2 SD: {m > 2 * boot.std()}")
v = A["T1"]["coef_sig"][..., 1][:, E["delay"]].mean(1)
print("by animal (T1): " + "; ".join(f"{a} n={int((pt & (animal == a)).sum())} mean {v[pt & (animal == a)].mean():.4f}" for a in np.unique(animal[pt])))
print(f"fraction of ALM PT ROIs with positive delay-mean beta {np.mean(v[pt] > 0):.3f}, negative {np.mean(v[pt] < 0):.3f}")
T1 all: ALM PT n=385, mean beta over the delay 0.0284, bootstrap SD 0.0166, ratio 1.71, above 2 SD: False
T1 variant, F0 > 0: ALM PT n=356, mean beta over the delay 0.0339, bootstrap SD 0.0172, ratio 1.97, above 2 SD: False
T2 all: ALM PT n=385, mean beta over the delay 0.0291, bootstrap SD 0.0187, ratio 1.56, above 2 SD: False
T2 variant, F0 > 0: ALM PT n=356, mean beta over the delay 0.0346, bootstrap SD 0.0186, ratio 1.86, above 2 SD: False by animal (T1): an063 n=106 mean 0.0341; an064 n=247 mean 0.0284; an066 n=32 mean 0.0093 fraction of ALM PT ROIs with positive delay-mean beta 0.208, negative 0.187
Result: not reproduced. The mean beta of ALM PT ROIs over the delay is positive (lick right, contralateral), 0.028 with T1 bins, but only 1.71 hierarchical bootstrap SDs above zero (T2 1.56), short of the 2 SD the criterion requires; without the ROIs with F0 <= 0 it is 1.97 (T1) and 1.86 (T2). The mean is positive in each of the three mice.
fig8d-pt-earlier-than-it¶
In ALM, lick direction selectivity (|beta|) was stronger and earlier in PT than in IT neurons, with an onset difference of 1.4 +- 0.6 s (p < 0.01, bootstrap).
Reading I1 (deep ALM IT) needs depth and cannot be run; reading I2 (all ALM IT) is run with bins T1 and T2, baselines O1 and O2, onset rules S1 (first bin over the threshold) and S2 (first bin from which the curve stays over it up to the go cue), and bootstraps H1 (hierarchical) and H2 (over cells). The threshold is the baseline mean plus 6 baseline SD, and each resample finds its onsets again by the same rule.
rng = np.random.default_rng(82)
pt, it = ALM & (ctype == "PT"), ALM & (ctype == "IT")
rows = []
for b in ("T1", "T2"):
t = A[b]["times"]; ab = np.abs(A[b]["coef_sig"][..., 1])
cp, ci = ab[pt].mean(0), ab[it].mean(0)
for hier in (True, False):
# One set of resamples per bootstrap; every baseline and rule finds its onset again in each.
d = {(o, r): [] for o in ("O1", "O2") for r in ("S1", "S2")}
for _ in range(1000):
ip, ii = C.joint_resample(animal, plane, [pt, it], rng, hier)
bp, bi = ab[ip].mean(0), ab[ii].mean(0)
for (o, r), lst in d.items():
lst.append(C.onset(bi, t, o, rule=r) - C.onset(bp, t, o, rule=r))
for (o, r), lst in d.items():
dd = np.array(lst)
on_p, on_i = C.onset(cp, t, o, rule=r), C.onset(ci, t, o, rule=r)
rows.append({"bins": b, "baseline": o, "rule": r, "bootstrap": "H1" if hier else "H2",
"PT_onset": on_p, "IT_onset": on_i, "IT_minus_PT": on_i - on_p,
"boot_sd": np.nanstd(dd), "p_not_earlier": np.mean(~(dd > 0)), "no_onset": int(np.isnan(dd).sum()),
"PT_max": cp.max(), "IT_max": ci.max()})
res = pd.DataFrame(rows).sort_values(["bins", "rule", "baseline", "bootstrap"])
res["within_1.4_pm_0.6"] = (res.IT_minus_PT - 1.4).abs() <= 0.6
res["p_below_0.01"] = res.p_not_earlier < 0.01
print(res.round(3).to_string(index=False))
print("\nCombinations within 1.4 +- 0.6 s:", int(res["within_1.4_pm_0.6"].sum()), "; with p < 0.01:", int(res["p_below_0.01"].sum()),
"; both:", int((res["within_1.4_pm_0.6"] & res["p_below_0.01"]).sum()), "of", len(res))
t = A["T1"]["times"]; ab = np.abs(A["T1"]["coef_sig"][..., 1])
plt.figure(figsize=(5.5, 3))
plt.plot(t, ab[pt].mean(0), label="ALM PT"); plt.plot(t, ab[it].mean(0), label="ALM IT")
for x in (0, 1.2, 4.2):
plt.axvline(x, color="k", lw=0.5)
plt.xlabel("time from trial start (s)"); plt.ylabel("mean |beta|"); plt.legend(); plt.show()
print(f"plotted range PT {ab[pt].mean(0).min():.4f} to {ab[pt].mean(0).max():.4f}, IT {ab[it].mean(0).min():.4f} to {ab[it].mean(0).max():.4f}")
for r in ("S1", "S2"):
print(f"by animal (T1, O1, {r}): " + "; ".join(
f"{a} PT onset {C.onset(ab[pt & (animal == a)].mean(0), t, 'O1', rule=r):.2f}, IT onset {C.onset(ab[it & (animal == a)].mean(0), t, 'O1', rule=r):.2f}"
for a in np.unique(animal[pt])))
bins baseline rule bootstrap PT_onset IT_onset IT_minus_PT boot_sd p_not_earlier no_onset PT_max IT_max within_1.4_pm_0.6 p_below_0.01 T1 O1 S1 H1 0.430 0.645 0.215 1.933 0.332 9 0.229 0.084 False False T1 O1 S1 H2 0.430 0.645 0.215 1.777 0.306 0 0.229 0.084 False False T1 O2 S1 H1 1.720 0.645 -1.075 1.582 0.476 0 0.229 0.084 False False T1 O2 S1 H2 1.720 0.645 -1.075 1.544 0.517 0 0.229 0.084 False False T1 O1 S2 H1 2.150 3.942 1.792 0.498 0.319 307 0.229 0.084 True False T1 O1 S2 H2 2.150 3.942 1.792 0.382 0.170 169 0.229 0.084 True False T1 O2 S2 H1 2.222 3.512 1.290 0.515 0.149 118 0.229 0.084 True False T1 O2 S2 H2 2.222 3.512 1.290 0.439 0.021 15 0.229 0.084 True False T2 O1 S1 H1 0.467 0.667 0.200 1.728 0.385 4 0.231 0.085 False False T2 O1 S1 H2 0.467 0.667 0.200 1.469 0.359 0 0.231 0.085 False False T2 O2 S1 H1 0.467 0.667 0.200 1.598 0.467 0 0.231 0.085 False False T2 O2 S1 H2 0.467 0.667 0.200 1.333 0.487 0 0.231 0.085 False False T2 O1 S2 H1 1.734 3.535 1.801 0.541 0.129 111 0.231 0.085 True False T2 O1 S2 H2 1.734 3.535 1.801 0.415 0.003 1 0.231 0.085 True True T2 O2 S2 H1 2.134 3.535 1.401 0.533 0.065 40 0.231 0.085 True False T2 O2 S2 H2 2.134 3.535 1.401 0.434 0.006 1 0.231 0.085 True True Combinations within 1.4 +- 0.6 s: 8 ; with p < 0.01: 2 ; both: 2 of 16
plotted range PT 0.0000 to 0.2287, IT 0.0041 to 0.0842 by animal (T1, O1, S1): an063 PT onset 0.43, IT onset 3.51; an064 PT onset 2.51, IT onset 3.30; an066 PT onset 2.29, IT onset 0.65 by animal (T1, O1, S2): an063 PT onset 2.15, IT onset 3.94; an064 PT onset 2.87, IT onset 3.51; an066 PT onset nan, IT onset nan
Result: reproduced. Mean |beta| is larger in ALM PT than in ALM IT ROIs (maximum 0.229 against 0.084, T1; 0.231 against 0.085, T2). With onset rule S2, T2 bins and the bootstrap over cells (H2), the IT onset follows the PT onset by 1.80 s (baseline O1, bootstrap SD 0.42 s, p = 0.003) and 1.40 s (O2, SD 0.43 s, p = 0.006), within 1.4 +- 0.6 s at p < 0.01, as the criterion requires under one combination. The result is not robust across readings: rule S2 puts the difference within tolerance in all 8 of its combinations (1.29 to 1.80 s), but with T1 bins p is 0.021 (O2) and 0.170 (O1) under H2, and 0.065 to 0.319 with the hierarchical bootstrap (H1); rule S1, the first single bin, gives 0.20 to 0.22 s, or -1.08 s (T1, O2), with p from 0.31 to 0.52. With S2 and H1 many resamples have no onset (up to 307 of 1,000, counted as PT not earlier). Per mouse (T1, O1, S1) PT is earlier in an063 and an064 and later in an066. Reading I1 (deep IT only) could not be run. The passing combinations use T2 bins of 66.7 ms interpolated from frames of 71.7 ms (13.95 Hz) or 70.4 ms (14.20 Hz), so neighbouring T2 bins are correlated and the onset is not sharper than a frame.
fig8d-mm-posterior-preferring¶
MM neurons were more active in posterior trials, when the pole was within reach of the whiskers, than in anterior trials.
rng = np.random.default_rng(83)
for b in ("T1", "T2"):
Eb = A[b]["epochs"]
v = A[b]["coef_sig"][..., 0][:, Eb["sample"]].mean(1)
for name, keep in (("all", np.ones(len(cells), bool)), ("variant, F0 > 0", cells.f0_min.values > 0)):
g_ = MM & keep
boot = np.array([v[C.joint_resample(animal, plane, [g_], rng)[0]].mean() for _ in range(1000)])
print(f"{b} {name}: MM n={g_.sum()}, mean alpha over the sample epoch {v[g_].mean():.4f}, bootstrap SD {boot.std():.4f}, "
f"ratio {v[g_].mean() / boot.std():.2f}, above 2 SD: {v[g_].mean() > 2 * boot.std()}")
v = A["T1"]["coef_sig"][..., 0][:, E["sample"]].mean(1)
print("by animal (T1): " + "; ".join(f"{a} n={int((MM & (animal == a)).sum())} mean {v[MM & (animal == a)].mean():.4f}" for a in np.unique(animal[MM])))
sig_a = A["T1"]["sig"][..., 0][:, E["sample"]]
K = A["T1"]["coef"][..., 0][:, E["sample"]]
print(f"MM ROI-bins in the sample epoch with significant alpha: posterior-preferring {(sig_a & (K > 0))[MM].sum()}, anterior-preferring {(sig_a & (K < 0))[MM].sum()}")
T1 all: MM n=6594, mean alpha over the sample epoch 0.0232, bootstrap SD 0.0124, ratio 1.87, above 2 SD: False
T1 variant, F0 > 0: MM n=5854, mean alpha over the sample epoch 0.0266, bootstrap SD 0.0157, ratio 1.70, above 2 SD: False
T2 all: MM n=6594, mean alpha over the sample epoch 0.0242, bootstrap SD 0.0125, ratio 1.94, above 2 SD: False
T2 variant, F0 > 0: MM n=5854, mean alpha over the sample epoch 0.0270, bootstrap SD 0.0166, ratio 1.63, above 2 SD: False by animal (T1): an041 n=1574 mean 0.0269; an043 n=940 mean 0.0140; an044 n=1505 mean 0.0087; an063 n=464 mean 0.0094; an064 n=2111 mean 0.0380 MM ROI-bins in the sample epoch with significant alpha: posterior-preferring 1912, anterior-preferring 747
Result: not reproduced. The mean alpha of MM ROIs over the sample epoch is positive (posterior preferring), 0.023 with T1 bins, but 1.87 hierarchical bootstrap SDs above zero (T2 1.94), short of the 2 SD of the criterion (1.70 and 1.63 without the ROIs with F0 <= 0). It is positive in all five mice with MM planes, and significant MM ROI-bins in the sample epoch prefer the posterior pole 1,912 times against 747.
fig8d-error-outnumber-correct¶
Overall, neurons preferring error trials outnumbered those preferring correct trials.
rng = np.random.default_rng(84)
for b in ("T1", "T2"):
t = A[b]["times"]; w = (t >= C.GO) & (t <= C.GO + 2.5 + 1e-9)
pg = A[b]["pure"][:, w, 2]; P = A[b]["p"][:, w, 2]; K = A[b]["coef"][:, w, 2]
best = np.where(pg, P, 2).argmin(1)
sgn_all = np.sign(K[np.arange(len(K)), best])
v = A[b]["coef_sig"][:, w, 2].mean(1)
for name, keep in (("V1, all ROIs", np.ones(len(cells), bool)), ("V2, F0 > 0 on every trial", cells.f0_min.values > 0)):
has = pg.any(1) & keep
sgn = sgn_all[has]
n_err, n_cor = int((sgn < 0).sum()), int((sgn > 0).sum())
bt = stats.binomtest(n_err, n_err + n_cor, 0.5)
boot = np.array([v[C.joint_resample(animal, plane, [keep], rng)[0]].mean() for _ in range(1000)])
m = v[keep].mean()
print(f"{b} {name}: (a) error-preferring {n_err}, correct-preferring {n_cor}, binomial p = {bt.pvalue:.2e}; "
f"(b) mean gamma {m:.4f}, bootstrap SD {boot.std():.4f}, ratio {m / boot.std():.2f}, below -2 SD: {m < -2 * boot.std()}")
t = A["T1"]["times"]; w = (t >= C.GO) & (t <= C.GO + 2.5 + 1e-9)
v = A["T1"]["coef_sig"][:, w, 2].mean(1)
print("by animal (T1), mean gamma: " + "; ".join(f"{a} {v[animal == a].mean():.4f}" for a in np.unique(animal)))
T1 V1, all ROIs: (a) error-preferring 2825, correct-preferring 2004, binomial p = 2.78e-32; (b) mean gamma -0.0082, bootstrap SD 0.0045, ratio -1.81, below -2 SD: False
T1 V2, F0 > 0 on every trial: (a) error-preferring 2722, correct-preferring 1927, binomial p = 1.74e-31; (b) mean gamma -0.0099, bootstrap SD 0.0045, ratio -2.20, below -2 SD: True
T2 V1, all ROIs: (a) error-preferring 2730, correct-preferring 1933, binomial p = 1.52e-31; (b) mean gamma -0.0078, bootstrap SD 0.0046, ratio -1.70, below -2 SD: False
T2 V2, F0 > 0 on every trial: (a) error-preferring 2628, correct-preferring 1859, binomial p = 1.43e-30; (b) mean gamma -0.0095, bootstrap SD 0.0047, ratio -2.00, below -2 SD: False by animal (T1), mean gamma: an041 -0.0104; an043 -0.0195; an044 -0.0088; an063 -0.0089; an064 -0.0048; an066 -0.0014
Result: not reproduced. Under V1, the reading that decides, reading (a) finds error preference: among ROIs purely outcome selective after the go cue, 2,825 prefer error trials and 2,004 correct trials (T1, binomial p = 2.8e-32), and T2 agrees (2,730 against 1,933). Reading (b) does not reach its bound: the mean gamma over all ROIs is negative (-0.0082, T1), and in every mouse, but 1.81 bootstrap SDs from zero (T2 1.70). The criterion requires both readings and names no partial outcome, so the claim is not reproduced, as for the other bootstrap claims whose direction holds but whose bound does not. Under V2, a check, (a) holds (2,722 against 1,927) and (b) reaches 2.20 SD with T1 bins, past the bound, and 2.00 SD with T2, which the printed test puts just short of it.
fig8d-mm-object-location (blocked; a part without depth)¶
MM neurons showed strong object location selectivity in superficial and deep layers and in PT and IT cells, PT cells particularly in the early delay; neither PT nor IT neurons in MM showed prominent lick direction selectivity before licking.
The criterion's groups are split by depth (deep PT, deep IT, deep other, superficial IT, superficial other), and the files hold no depth, so the claim is blocked. PT neurons are deep by the paper's account (fig8c-pt-deep-only, whose depth part is itself untested), so the PT group of part (a), MM PT against ALM PT, is the one comparison that needs no depth in the files; it rests on that assumption and does not decide. The same comparisons for IT and other ROIs of the three labelled mice, without the depth split, and the quantities of parts (b) and (c) without depth, are printed for reference.
rng = np.random.default_rng(87)
labelled_mice = np.unique(animal[np.isin(ctype, ["PT", "IT"])])
lab_m = np.isin(animal, labelled_mice)
proj = np.where(np.isin(ctype, ["PT", "IT"]), ctype, "other")
v2 = cells.f0_min.values > 0
for b in ("T1", "T2"):
t = A[b]["times"]; w = (t >= 0) & (t <= 2.2 + 1e-9); ed = (t >= C.SAMPLE_END) & (t < C.SAMPLE_END + 1.0)
aa, ab = np.abs(A[b]["coef_sig"][..., 0]), np.abs(A[b]["coef_sig"][..., 1])
for vname, keep in (("V1", np.ones(len(cells), bool)), ("V2, F0 > 0 on every trial", v2)):
print(f"{b} {vname}: (a) maximum of mean |alpha| from trial start to 2.2 s, MM against ALM of the same projection "
"type (three labelled mice" + ("; hierarchical bootstrap, 1000 resamples)" if b == "T1" else ")"))
for p_ in ("PT", "IT", "other"):
gm, ga = MM & lab_m & (proj == p_) & keep, ALM & lab_m & (proj == p_) & keep
mm_, al_ = aa[gm][:, w].mean(0).max(), aa[ga][:, w].mean(0).max()
extra = ""
if b == "T1":
bm, ba = [], []
for _ in range(1000):
im, ia = C.joint_resample(animal, plane, [gm, ga], rng)
bm.append(aa[im][:, w].mean(0).max() if len(im) else np.nan)
ba.append(aa[ia][:, w].mean(0).max() if len(ia) else np.nan)
bm, ba = np.array(bm), np.array(ba)
both = np.isfinite(bm) & np.isfinite(ba)
extra = (f"; bootstrap SD {np.nanstd(bm):.4f} (MM) and {np.nanstd(ba):.4f} (ALM), resamples without MM or "
f"ALM cells {int((~both).sum())}, MM above ALM in {np.mean(bm[both] > ba[both]):.3f} of the others")
print(f" {p_}{' (part a for PT, assuming PT deep)' if p_ == 'PT' else ' (no depth split, reference)'}: "
f"MM n={gm.sum()} max {mm_:.4f}, ALM n={ga.sum()} max {al_:.4f}, MM larger: {mm_ > al_}{extra}")
first_delay = {p_: aa[MM & lab_m & (proj == p_) & keep][:, ed].mean() for p_ in ("PT", "IT", "other")}
print(" (b) reference, mean |alpha| over the first second of the delay in MM: "
+ ", ".join(f"{k} {v:.4f}" for k, v in first_delay.items()))
pre = t < C.GO
mm_pt_b = ab[MM & (proj == 'PT') & keep][:, pre].mean(0).max()
print(f" (c) reference, maximum mean |beta| before the go cue: MM PT {mm_pt_b:.4f}, "
f"MM IT {ab[MM & (proj == 'IT') & keep][:, pre].mean(0).max():.4f}, ALM PT {ab[ALM & (proj == 'PT') & keep][:, pre].mean(0).max():.4f} "
f"(the criterion's ALM deep maximum needs depth; for MM PT to pass it must exceed {2 * mm_pt_b:.4f})")
pt_mm = MM & (proj == "PT")
print(f"\nMM PT ROIs with F0 <= 0 on some trial: {int((pt_mm & ~v2).sum())} of {int(pt_mm.sum())}")
print("MM PT ROIs by animal: " + ", ".join(
f"{a} {int((pt_mm & (animal == a)).sum())} ROIs in {len(np.unique(plane[pt_mm & (animal == a)]))} planes, "
f"{int((pt_mm & ~v2 & (animal == a)).sum())} of them with F0 <= 0 on some trial" for a in labelled_mice))
T1 V1: (a) maximum of mean |alpha| from trial start to 2.2 s, MM against ALM of the same projection type (three labelled mice; hierarchical bootstrap, 1000 resamples)
PT (part a for PT, assuming PT deep): MM n=334 max 0.6233, ALM n=385 max 0.0327, MM larger: True; bootstrap SD 0.6168 (MM) and 0.0203 (ALM), resamples without MM or ALM cells 29, MM above ALM in 0.977 of the others
IT (no depth split, reference): MM n=758 max 0.0570, ALM n=1873 max 0.0274, MM larger: True; bootstrap SD 0.0250 (MM) and 0.0252 (ALM), resamples without MM or ALM cells 43, MM above ALM in 0.820 of the others
other (no depth split, reference): MM n=1483 max 0.1224, ALM n=3263 max 0.0246, MM larger: True; bootstrap SD 0.0425 (MM) and 0.0088 (ALM), resamples without MM or ALM cells 45, MM above ALM in 1.000 of the others (b) reference, mean |alpha| over the first second of the delay in MM: PT 0.2354, IT 0.0412, other 0.0945 (c) reference, maximum mean |beta| before the go cue: MM PT 0.5233, MM IT 0.0624, ALM PT 0.1119 (the criterion's ALM deep maximum needs depth; for MM PT to pass it must exceed 1.0467) T1 V2, F0 > 0 on every trial: (a) maximum of mean |alpha| from trial start to 2.2 s, MM against ALM of the same projection type (three labelled mice; hierarchical bootstrap, 1000 resamples)
PT (part a for PT, assuming PT deep): MM n=287 max 0.7037, ALM n=356 max 0.0288, MM larger: True; bootstrap SD 0.7270 (MM) and 0.0131 (ALM), resamples without MM or ALM cells 44, MM above ALM in 0.995 of the others
IT (no depth split, reference): MM n=584 max 0.0573, ALM n=1704 max 0.0124, MM larger: True; bootstrap SD 0.0258 (MM) and 0.0056 (ALM), resamples without MM or ALM cells 35, MM above ALM in 0.964 of the others
other (no depth split, reference): MM n=1165 max 0.1365, ALM n=3055 max 0.0210, MM larger: True; bootstrap SD 0.0379 (MM) and 0.0066 (ALM), resamples without MM or ALM cells 36, MM above ALM in 1.000 of the others (b) reference, mean |alpha| over the first second of the delay in MM: PT 0.2604, IT 0.0437, other 0.1044 (c) reference, maximum mean |beta| before the go cue: MM PT 0.6090, MM IT 0.0143, ALM PT 0.1184 (the criterion's ALM deep maximum needs depth; for MM PT to pass it must exceed 1.2181) T2 V1: (a) maximum of mean |alpha| from trial start to 2.2 s, MM against ALM of the same projection type (three labelled mice) PT (part a for PT, assuming PT deep): MM n=334 max 0.5859, ALM n=385 max 0.0331, MM larger: True IT (no depth split, reference): MM n=758 max 0.0536, ALM n=1873 max 0.0265, MM larger: True other (no depth split, reference): MM n=1483 max 0.1213, ALM n=3263 max 0.0244, MM larger: True (b) reference, mean |alpha| over the first second of the delay in MM: PT 0.2203, IT 0.0418, other 0.0953 (c) reference, maximum mean |beta| before the go cue: MM PT 0.4824, MM IT 0.0474, ALM PT 0.1190 (the criterion's ALM deep maximum needs depth; for MM PT to pass it must exceed 0.9648) T2 V2, F0 > 0 on every trial: (a) maximum of mean |alpha| from trial start to 2.2 s, MM against ALM of the same projection type (three labelled mice) PT (part a for PT, assuming PT deep): MM n=287 max 0.6656, ALM n=356 max 0.0292, MM larger: True IT (no depth split, reference): MM n=584 max 0.0558, ALM n=1704 max 0.0119, MM larger: True other (no depth split, reference): MM n=1165 max 0.1355, ALM n=3055 max 0.0209, MM larger: True (b) reference, mean |alpha| over the first second of the delay in MM: PT 0.2439, IT 0.0438, other 0.1056 (c) reference, maximum mean |beta| before the go cue: MM PT 0.5614, MM IT 0.0208, ALM PT 0.1287 (the criterion's ALM deep maximum needs depth; for MM PT to pass it must exceed 1.1228) MM PT ROIs with F0 <= 0 on some trial: 47 of 334 MM PT ROIs by animal: an063 54 ROIs in 4 planes, 2 of them with F0 <= 0 on some trial, an064 280 ROIs in 7 planes, 45 of them with F0 <= 0 on some trial, an066 0 ROIs in 0 planes, 0 of them with F0 <= 0 on some trial
Result: blocked (metadata insufficient), not decided here. Under the paper's account that PT neurons are deep, part (a) holds for PT: the maximum of mean |alpha| from trial start to 2.2 s is 0.623 in MM PT against 0.033 in ALM PT (T1, V1; T2 0.586 against 0.033; V2 0.704 against 0.029), larger in 0.977 of the resamples that hold both groups (29 of 1,000 lack one). The MM PT mean is uncertain: its bootstrap SD (0.617) is as large as the mean, and 280 of the 334 MM PT ROIs are in 7 planes of one mouse (an064). Without the depth split, MM IT (0.057 against 0.027) and other ROIs (0.122 against 0.025) are also above ALM. The other groups of part (a), and parts (b) and (c), need depth, so the claim stays blocked. Part (c) looks unlikely to hold for PT on the same account: the maximum of mean |beta| before the go cue is 0.523 in MM PT (T1, V1), so the ALM deep maximum would have to exceed 1.047, against 0.112 in ALM PT. The large MM PT values may come from a few ROIs: 280 of the 334 MM PT ROIs, and 45 of the 47 with F0 <= 0 on some trial, are in an064.
fig8ef-area-explains-variance (blocked; area against projection type)¶
Coding is mainly determined by cortical location; cortical area (ALM against MM) and cell depth explain most of the variance of selectivity, more than projection type.
The criterion needs the depth grouping in both of its conditions, and the files hold no depth, so the claim is blocked. Area and projection type are in the files, and their total variance explained is computed here as the criterion defines it (SSB over SST, summed over the bins from 4.6 s before to 2.5 s after the go cue and over |alpha|, |beta| and |gamma|, non-significant set to 0), under readings N1 (three labelled mice for both groupings) and N2 (all six mice for area). It does not decide: if area explains less than projection type, area is not the largest and the claim is at best partially reproduced, depending on depth; if more, depth decides whether area is the largest.
rng = np.random.default_rng(86)
names3 = ["alpha", "beta", "gamma"]
def var_explained(idx, b, w, labels):
X = A[b]["coef_sig"][idx][:, w]
ssb = sst = 0.0
per = {}
for k in range(3):
bk, tk = C.ss_between(np.abs(X[..., k]), labels)
ssb, sst = ssb + bk.sum(), sst + tk.sum()
per[f"|{names3[k]}|"] = bk.sum() / tk.sum()
bs, ts = C.ss_between(X[..., k], labels)
per[names3[k]] = bs.sum() / ts.sum()
return ssb / sst, per
for b in ("T1", "T2"):
t = A[b]["times"]; w = (t >= -0.4 - 1e-9) & (t <= 6.7 + 1e-9)
for vname, keep in (("V1", np.ones(len(cells), bool)), ("V2, F0 > 0 on every trial", v2)):
i1, iall = np.flatnonzero(lab_m & keep), np.flatnonzero(keep)
res = {"area, N1": var_explained(i1, b, w, area[i1]), "projection type, N1": var_explained(i1, b, w, proj[i1]),
"area, N2": var_explained(iall, b, w, area[iall])}
print(f"{b} {vname} (bins {w.sum()}, ROIs N1 {len(i1)}, N2 {len(iall)}):")
for k, (tot, per) in res.items():
print(f" {k}: total {tot:.4f}; per coefficient " + ", ".join(f"{c} {v:.4f}" for c, v in per.items()))
print(f" area above projection type (N1): {res['area, N1'][0] > res['projection type, N1'][0]}")
if b == "T1":
bt = []
for _ in range(1000):
ix = C.joint_resample(animal, plane, [lab_m & keep], rng)[0]
bt.append((var_explained(ix, b, w, area[ix])[0], var_explained(ix, b, w, proj[ix])[0]))
bt = np.array(bt)
print(f" hierarchical bootstrap (N1, 1000 resamples): SD area {bt[:, 0].std():.4f}, projection type "
f"{bt[:, 1].std():.4f}; resamples with area above projection type {np.mean(bt[:, 0] > bt[:, 1]):.3f}")
t = A["T1"]["times"]; w = (t >= -0.4 - 1e-9) & (t <= 6.7 + 1e-9)
print("by animal (T1, V1), projection type total: " + "; ".join(
f"{a} {var_explained(np.flatnonzero(animal == a), 'T1', w, proj[animal == a])[0]:.4f}" for a in labelled_mice))
print("by animal (T1, V1), area total: " + "; ".join(
f"{a} {var_explained(np.flatnonzero(animal == a), 'T1', w, area[animal == a])[0]:.4f}"
if len(np.unique(area[animal == a])) > 1 else f"{a} one area only" for a in labelled_mice))
T1 V1 (bins 99, ROIs N1 8096, N2 15431): area, N1: total 0.0005; per coefficient |alpha| 0.0008, alpha 0.0005, |beta| 0.0004, beta 0.0004, |gamma| 0.0003, gamma 0.0003 projection type, N1: total 0.0013; per coefficient |alpha| 0.0013, alpha 0.0012, |beta| 0.0014, beta 0.0010, |gamma| 0.0012, gamma 0.0011 area, N2: total 0.0007; per coefficient |alpha| 0.0010, alpha 0.0005, |beta| 0.0007, beta 0.0002, |gamma| 0.0002, gamma 0.0001 area above projection type (N1): False
hierarchical bootstrap (N1, 1000 resamples): SD area 0.0019, projection type 0.0011; resamples with area above projection type 0.414 T1 V2, F0 > 0 on every trial (bins 99, ROIs N1 7151, N2 14169): area, N1: total 0.0005; per coefficient |alpha| 0.0009, alpha 0.0008, |beta| 0.0004, beta 0.0007, |gamma| 0.0004, gamma 0.0004 projection type, N1: total 0.0015; per coefficient |alpha| 0.0015, alpha 0.0015, |beta| 0.0016, beta 0.0013, |gamma| 0.0014, gamma 0.0013 area, N2: total 0.0008; per coefficient |alpha| 0.0012, alpha 0.0007, |beta| 0.0008, beta 0.0003, |gamma| 0.0003, gamma 0.0001 area above projection type (N1): False
hierarchical bootstrap (N1, 1000 resamples): SD area 0.0026, projection type 0.0015; resamples with area above projection type 0.400 T2 V1 (bins 106, ROIs N1 8096, N2 15431): area, N1: total 0.0005; per coefficient |alpha| 0.0008, alpha 0.0005, |beta| 0.0004, beta 0.0005, |gamma| 0.0004, gamma 0.0002 projection type, N1: total 0.0013; per coefficient |alpha| 0.0013, alpha 0.0012, |beta| 0.0015, beta 0.0010, |gamma| 0.0013, gamma 0.0010 area, N2: total 0.0008; per coefficient |alpha| 0.0013, alpha 0.0006, |beta| 0.0009, beta 0.0002, |gamma| 0.0003, gamma 0.0001 area above projection type (N1): False
T2 V2, F0 > 0 on every trial (bins 106, ROIs N1 7151, N2 14169): area, N1: total 0.0006; per coefficient |alpha| 0.0010, alpha 0.0009, |beta| 0.0004, beta 0.0008, |gamma| 0.0004, gamma 0.0004 projection type, N1: total 0.0015; per coefficient |alpha| 0.0015, alpha 0.0015, |beta| 0.0016, beta 0.0012, |gamma| 0.0014, gamma 0.0013 area, N2: total 0.0009; per coefficient |alpha| 0.0014, alpha 0.0008, |beta| 0.0010, beta 0.0003, |gamma| 0.0003, gamma 0.0001 area above projection type (N1): False by animal (T1, V1), projection type total: an063 0.0026; an064 0.0018; an066 0.0021 by animal (T1, V1), area total: an063 0.0046; an064 0.0004; an066 one area only
Result: blocked (metadata insufficient), not decided here. Projection type explains more of the variance of selectivity than area in every reading run: total 0.0013 against 0.0005 under N1 (T1, V1; area under N2 0.0007), and the same order with T2 and with V2. Both totals are small, and the hierarchical bootstrap does not separate them (SD 0.0019 for area and 0.0011 for projection type; area above projection type in 0.414 of resamples, 0.400 under V2). If this order holds, area is not the largest of the three groupings, so the claim cannot be reproduced; whether it is partially reproduced or not depends on depth (area and depth together against projection type), which is not in the files.
alm-pt-it-lick-fraction-below-20¶
Even in ALM, less than 20% of PT or IT cells displayed lick direction selectivity at any given time, and the onset timing was variable.
pt, it = ALM & (ctype == "PT"), ALM & (ctype == "IT")
for b in ("T1", "T2"):
t = A[b]["times"]
for name, X in (("(a) purely lick-direction", A[b]["pure"][..., 1]), ("(b) beta significant", A[b]["sig"][..., 1])):
fp, fi = X[pt].mean(0), X[it].mean(0)
print(f"{b} {name}: max PT {100 * fp.max():.1f}% at {t[fp.argmax()]:.2f} s, max IT {100 * fi.max():.1f}% at {t[fi.argmax()]:.2f} s; "
f"both below 22%: {fp.max() < 0.22 and fi.max() < 0.22}")
X = A["T1"]["sig"][..., 1]
print("by animal (T1, reading b), max PT and IT: " + "; ".join(
f"{a} {100 * X[pt & (animal == a)].mean(0).max():.1f}% / {100 * X[it & (animal == a)].mean(0).max():.1f}%" for a in np.unique(animal[pt])))
T1 (a) purely lick-direction: max PT 12.7% at 4.44 s, max IT 3.3% at 4.44 s; both below 22%: True T1 (b) beta significant: max PT 15.6% at 4.44 s, max IT 6.0% at 4.87 s; both below 22%: True T2 (a) purely lick-direction: max PT 13.0% at 4.47 s, max IT 3.6% at 4.40 s; both below 22%: True T2 (b) beta significant: max PT 16.1% at 4.54 s, max IT 6.6% at 4.94 s; both below 22%: True by animal (T1, reading b), max PT and IT: an063 19.8% / 8.6%; an064 15.0% / 4.7%; an066 18.8% / 4.6%
Result: reproduced. The largest per-bin fraction of lick-direction-selective ALM PT ROIs is 12.7% (purely selective) or 15.6% (beta significant), and of ALM IT ROIs 3.3% or 6.0% (T1); T2 gives similar values (16.1% at most). All stay below 22%, also in each mouse (at most 19.8%). The variability of onset is not judged.
li2015-tagged-units-link¶
The tagged units of DANDI:000010 (Li et al. 2015) are the 45 PT and 27 IT cells the paper compares its imaging with.
li = C.li2015_tagged_units()
types = collections.Counter(t for r in li for t in r["cell_type"])
print("files:", len(li), " files with units:", sum(len(r["cell_type"]) > 0 for r in li))
print("cell_type counts over all units:", dict(types))
print("PT units by genotype:", dict(collections.Counter(r["genotype"] for r in li for t in r["cell_type"] if t == "PT")))
print("IT units by genotype:", dict(collections.Counter(r["genotype"] for r in li for t in r["cell_type"] if t == "IT")))
files: 158 files with units: 99
cell_type counts over all units: {np.str_('FS'): 124, np.str_('Pyr'): 1245, np.str_('IT'): 27, np.str_('PT'): 45}
PT units by genotype: {'rbp4-cre': 5, 'sim1-cre': 40}
IT units by genotype: {'rbp4-cre': 3, 'chr2-eyfp x tlx-cre': 8, 'tlx-cre': 16}
Result: reproduced. DANDI:000010 version 0.220126.1905 holds 45 units labelled PT and 27 labelled IT, the counts the paper gives for the electrophysiology (45/27). The units are matched by count and label only; no identifier links them to the paper's analysis.
ephys-pt-it-onset-not-detected¶
The difference in onset of lick direction selectivity between PT and IT neurons was not detected in an electrophysiology dataset recorded under similar conditions (45 PT and 27 IT cells, Li et al., 2015).
The readings are those the criterion names since the question was answered: (1) spike counts in
400 ms windows in 100 ms steps, each window placed at its centre; (2) trials of the four types
without photostimulation, early-lick trials excluded by the early_lick column, time from the
sample event (pole movement onset); (3) onset by convention O from trial start to the go cue, rules
S1 and S2, baselines O1 (windows wholly before trial start) and O2 (windows wholly within the 0.5 s
before it); (4) units of all sessions pooled, bootstrap H2 over units and H1 hierarchical over
animals, sessions and units. The 3-way ANOVA and the coefficients are those of convention M, applied
to the counts of each unit. A trial is photostimulated when a laser pulse (from the laser power
series) or a photostim_start_time event falls within it. Only the sessions with a tagged unit are
read, and of each unit only the trials within its observation intervals.
started_li = time.time()
li_rows = [r for r in li if np.isin(r["cell_type"], ["PT", "IT"]).any()]
S = C.li2015_sessions(li_rows)
print(f"read {len(S)} sessions with tagged units in {time.time() - started_li:.0f} s")
off = np.concatenate([s["sample"] - s["start"] for s in S])
print(f"trials {len(off)}, without a matched sample event {int(np.isnan(off).sum())}; "
f"sample - start_time {np.nanmin(off):.3f} to {np.nanmax(off):.3f} s")
ll = np.concatenate([s["light_laser"] for s in S]); le = np.concatenate([s["light_event"] for s in S])
print(f"photostimulated trials: by laser pulses {int(ll.sum())}, by photostim events {int(le.sum())}, both {int((ll & le).sum())}, "
f"either {int((ll | le).sum())}; sessions with a laser series {sum(s['laser_series'] > 0 for s in S)}; "
f"pulses and events outside any trial {sum(s['pulses_outside'] for s in S)}, {sum(s['events_outside'] for s in S)}")
li_tt = [C.li2015_types(s) for s in S]
early_li = np.concatenate([s["early"] for s in S])
print("early_lick values:", dict(collections.Counter(early_li)), "; outcome among early-lick trials:",
dict(collections.Counter(np.concatenate([s["outcome"][s["early"] == "early"] for s in S]))))
print("trial types used (all sessions):", dict(collections.Counter(np.concatenate(li_tt))))
rel = lambda key: np.concatenate([(s[key] - s["sample"])[t != ""] for s, t in zip(S, li_tt)])
go_rel, start_rel, stop_rel = rel("go"), rel("start"), rel("stop")
go_li = float(np.nanmedian(go_rel))
print(f"used trials: go - sample median {go_li:.3f} s, {int((np.abs(go_rel - go_li) > 0.01).sum())} trials more than 10 ms off "
f"and {int(np.isnan(go_rel).sum())} without a go event (kept, as time is taken from the sample event); "
f"stored span from {start_rel.max():.3f} to {stop_rel.min():.3f} s around the sample event in every used trial")
u_obs = [np.mean([np.any((u["obs"][:, 0] <= a + 1e-6) & (u["obs"][:, 1] >= b - 1e-6))
for a, b in zip(s["start"], s["stop"])]) for s in S for u in s["units"]]
print(f"fraction of a session's trials within a unit's observation intervals: {min(u_obs):.3f} to {max(u_obs):.3f}")
print("units:", dict(collections.Counter(u["label"] for s in S for u in s["units"])),
"; mice:", len({s["subject"] for s in S}))
read 50 sessions with tagged units in 68 s
trials 16571, without a matched sample event 1636; sample - start_time 0.524 to 0.573 s
photostimulated trials: by laser pulses 4142, by photostim events 4080, both 4070, either 4152; sessions with a laser series 44; pulses and events outside any trial 460, 76
early_lick values: {np.str_('no early'): 15251, np.str_('early'): 1320} ; outcome among early-lick trials: {np.str_('hit'): 996, np.str_('miss'): 319, np.str_('ignore'): 5}
trial types used (all sessions): {'': 6743, 'ER': 1040, 'CR': 3815, 'EL': 1252, 'CL': 3721}
used trials: go - sample median 2.601 s, 44 trials more than 10 ms off and 47 without a go event (kept, as time is taken from the sample event); stored span from -0.524 to 4.927 s around the sample event in every used trial
fraction of a session's trials within a unit's observation intervals: 0.621 to 0.956
units: {np.str_('IT'): 27, np.str_('PT'): 45} ; mice: 14
# Windows of 400 ms in 100 ms steps, edges on a grid from the sample event, wholly within the span
# stored in every used trial (spikes are stored only within the trial windows).
li_centres = np.round(np.arange(np.ceil((start_rel.max() + C.EPH_WIN / 2) / C.EPH_STEP - 1e-9),
np.floor((stop_rel.min() - C.EPH_WIN / 2) / C.EPH_STEP + 1e-9) + 1) * C.EPH_STEP, 3)
li_base = {"O1": li_centres + C.EPH_WIN / 2 <= 1e-9,
"O2": (li_centres - C.EPH_WIN / 2 >= -0.5 - 1e-9) & (li_centres + C.EPH_WIN / 2 <= 1e-9)}
print(f"window centres {li_centres[0]:.1f} to {li_centres[-1]:.1f} s; baseline windows (centres): "
+ "; ".join(f"{k} {list(li_centres[v])}" for k, v in li_base.items()))
LB, li_lab, li_animal, li_sess, li_ntr, dropped = [], [], [], [], [], []
for si, s in enumerate(S):
for u in s["units"]:
Y, types, _ = C.li2015_unit_counts(s, u, li_centres)
cnt = {t: int((types == t).sum()) for t in C.TYPES}
if min(cnt.values()) < 1:
dropped.append((s["path"], u["label"], cnt))
continue
p, coef = C.anova3(Y, types)
LB.append(np.where(p[1] < C.ALPHA, np.abs(coef[1]), 0))
li_lab.append(u["label"]); li_animal.append(s["subject"]); li_sess.append(si); li_ntr.append(len(types))
LB, li_lab, li_animal, li_sess = np.array(LB), np.array(li_lab), np.array(li_animal), np.array(li_sess)
print("units without a trial of every type (no fit, left out):", dropped)
print("units analysed:", dict(collections.Counter(li_lab)), f"; trials per unit {min(li_ntr)} to {max(li_ntr)} (median {np.median(li_ntr):.0f})")
lpt, lit = li_lab == "PT", li_lab == "IT"
cp, ci = LB[lpt].mean(0), LB[lit].mean(0)
rng = np.random.default_rng(85)
rows = []
for hier in (True, False):
d = {(o, r): [] for o in ("O1", "O2") for r in ("S1", "S2")}
for _ in range(1000):
ip, ii = C.joint_resample(li_animal, li_sess, [lpt, lit], rng, hier)
bp, bi = LB[ip].mean(0), LB[ii].mean(0)
for (o, r), lst in d.items():
on = lambda c: C.onset(c, li_centres, rule=r, go=go_li, base=li_base[o], search_end=go_li)
lst.append(on(bi) - on(bp))
for (o, r), lst in d.items():
on = lambda c: C.onset(c, li_centres, rule=r, go=go_li, base=li_base[o], search_end=go_li)
dd = np.array(lst)
rows.append({"baseline": o, "rule": r, "bootstrap": "H1" if hier else "H2", "PT_onset": on(cp), "IT_onset": on(ci),
"IT_minus_PT": on(ci) - on(cp), "boot_sd": np.nanstd(dd), "p_not_earlier": np.mean(~(dd > 0)),
"no_onset": int(np.isnan(dd).sum())})
lres = pd.DataFrame(rows).sort_values(["rule", "baseline", "bootstrap"])
lres["PT_earlier_p_below_0.01"] = lres.p_not_earlier < 0.01
print(lres.round(3).to_string(index=False))
print("\nCombinations finding PT earlier at p < 0.01:", int(lres["PT_earlier_p_below_0.01"].sum()), "of", len(lres))
for o in ("O1", "O2"):
for lab_, c in (("PT", cp), ("IT", ci)):
bm = c[li_base[o]]
print(f"{o} {lab_}: baseline mean {bm.mean():.4f}, SD {bm.std(ddof=1):.4f}, threshold {bm.mean() + 6 * bm.std(ddof=1):.4f}")
plt.figure(figsize=(5.5, 3))
plt.plot(li_centres, cp, label=f"PT (n={int(lpt.sum())})"); plt.plot(li_centres, ci, label=f"IT (n={int(lit.sum())})")
for x in (0, 1.3, go_li):
plt.axvline(x, color="k", lw=0.5)
plt.xlabel("window centre, time from the sample event (s)"); plt.ylabel("mean |beta| (spikes per window)"); plt.legend(); plt.show()
print(f"plotted range PT {cp.min():.3f} to {cp.max():.3f}, IT {ci.min():.3f} to {ci.max():.3f}")
print("by animal (O1, S1), onset PT / IT: " + "; ".join(
f"{a} {C.onset(LB[lpt & (li_animal == a)].mean(0), li_centres, base=li_base['O1'], search_end=go_li) if (lpt & (li_animal == a)).any() else np.nan:.1f} / "
f"{C.onset(LB[lit & (li_animal == a)].mean(0), li_centres, base=li_base['O1'], search_end=go_li) if (lit & (li_animal == a)).any() else np.nan:.1f}"
for a in np.unique(li_animal)))
window centres -0.3 to 4.7 s; baseline windows (centres): O1 [np.float64(-0.3), np.float64(-0.2)]; O2 [np.float64(-0.3), np.float64(-0.2)]
units without a trial of every type (no fit, left out): [('sub-219253/sub-219253_ses-20140118_obj-1m69pde_behavior+ecephys+ogen.nwb', np.str_('IT'), {'CR': 27, 'CL': 41, 'ER': 0, 'EL': 15})]
units analysed: {np.str_('IT'): 26, np.str_('PT'): 45} ; trials per unit 127 to 318 (median 193)
baseline rule bootstrap PT_onset IT_onset IT_minus_PT boot_sd p_not_earlier no_onset PT_earlier_p_below_0.01
O1 S1 H1 0.1 0.8 0.7 0.752 0.298 5 False
O1 S1 H2 0.1 0.8 0.7 0.661 0.342 0 False
O2 S1 H1 0.1 0.8 0.7 0.752 0.298 5 False
O2 S1 H2 0.1 0.8 0.7 0.661 0.342 0 False
O1 S2 H1 1.2 1.4 0.2 0.629 0.203 65 False
O1 S2 H2 1.2 1.4 0.2 0.638 0.153 8 False
O2 S2 H1 1.2 1.4 0.2 0.629 0.203 65 False
O2 S2 H2 1.2 1.4 0.2 0.638 0.153 8 False
Combinations finding PT earlier at p < 0.01: 0 of 8
O1 PT: baseline mean 0.0367, SD 0.0050, threshold 0.0664
O1 IT: baseline mean 0.0205, SD 0.0004, threshold 0.0228
O2 PT: baseline mean 0.0367, SD 0.0050, threshold 0.0664
O2 IT: baseline mean 0.0205, SD 0.0004, threshold 0.0228
plotted range PT 0.000 to 0.644, IT 0.000 to 1.210 by animal (O1, S1), onset PT / IT: 214427 2.5 / 2.3; 214429 0.9 / nan; 218453 nan / 0.8; 218457 nan / 0.8; 219030 0.1 / nan; 219031 2.4 / nan; 219033 1.4 / nan; 219036 0.0 / nan; 219037 1.6 / nan; 219038 1.4 / nan; 219247 nan / nan; 219248 nan / 1.0; 219253 nan / 1.7; 221977 nan / 2.3
Result: reproduced. No combination of readings finds the 45 PT units earlier than the 26 IT units at p < 0.01: the IT minus PT onset is 0.7 s with rule S1 (bootstrap SD 0.66 to 0.75 s, p 0.30 to 0.34) and 0.2 s with rule S2 (SD 0.63 to 0.64 s, p 0.15 to 0.20). Baselines O1 and O2 are the same two windows, as spikes are stored only from 0.52 s before the sample event, so the 8 combinations give 4 distinct results. Two overlapping windows give a weak baseline: the IT baseline SD is 0.0004 (threshold 0.0228 against a mean of 0.0205), and the PT onset under S1 (0.1 s) falls in the sample epoch, so the test carries little information. The onset difference is positive in every combination, smaller than in the imaging and not significant. One IT unit (mouse 219253) has no ER trial without photostimulation and cannot be fitted; 26 of the 27 IT units are analysed. 4,152 of 16,571 trials in these sessions had photostimulation and are left out.
Summary¶
| Claim | Outcome |
|---|---|
| early-lick-fraction | Reproduced |
| imaging-sample | Partially reproduced (locations not identifiable) |
| trials-per-type | Reproduced |
| neuron-counts | Reproduced |
| fraction-selective-alm-mm | Not reproduced |
| figs5-single-coefficient | Reproduced (the permuted null meets the criterion too) |
| fig4a-lick-direction-epochs | Reproduced |
| fig4b-object-location-memory | Reproduced |
| fig4c-outcome-after-lick | Not reproduced |
| fig5-mixed-response-epoch | Reproduced |
| fig6-object-location-mm | Reproduced |
| fig6-lick-direction-deep-alm-first | Blocked: metadata insufficient (no depth) |
| fig6-lick-direction-spread | Partially reproduced (depth part not testable) |
| fig6-outcome-response-uniform | Not reproduced |
| fig6h-preparatory-alm-object-mm | Reproduced |
| figs9-photoinhibition-alm-mm | Blocked: data missing |
| widefield-mm-to-alm | Blocked: data missing |
| widefield-contralateral-bilateral | Blocked: data missing |
| figs4-peak-shift | Not reproduced (whole window; the peak before the go cue would pass, question open) |
| fig8c-pt-deep-only | Partially reproduced (counts only) |
| fig8d-alm-deep-lick-direction | Blocked: metadata insufficient (no depth) |
| fig8d-pt-contra-bias | Not reproduced |
| fig8d-pt-earlier-than-it | Reproduced (rule S2, T2 bins, bootstrap over cells only) |
| ephys-pt-it-onset-not-detected | Reproduced |
| fig8d-mm-object-location | Blocked: metadata insufficient (no depth; if PT are deep, the PT part of (a) holds and that of (c) looks unlikely) |
| fig8d-mm-posterior-preferring | Not reproduced |
| fig8d-error-outnumber-correct | Not reproduced (reading a holds, reading b short of its bound) |
| fig8ef-area-explains-variance | Blocked: metadata insufficient (no depth; area explains less than projection type, not significantly) |
| ephys-selective-counts | Blocked: the silicon probe data are not public (placeholder DOI, not on DANDI) |
| fig7d-ephys-area-distribution | Blocked: the silicon probe data are not public (placeholder DOI, not on DANDI) |
| fig1c-whisking-similar | Blocked: data missing |
| alm-pt-it-lick-fraction-below-20 | Reproduced |
| li2015-tagged-units-link (link) | Reproduced |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. The time bins (frames on one grid, or 66.7 ms by interpolation), dF/F with baselines F0 <= 0 kept (V1), the onset rules (S1, S2) and the bootstrap details are readings from the claims file's conventions, which the paper does not fix (five questions were answered, each with "the paper does not say"); several claims use thresholds the claims file fixed itself (2 bootstrap SD, more than half of bins). Area comes from the files' label. Readings that need depth or coordinates could not be run. The PT and IT onset difference matches the paper only under one onset rule (S2), only with T2 bins and only with the bootstrap over cells. The Li et al. (2015) units are analysed with the paper's 400 ms windows on trials without photostimulation; one IT unit cannot be fitted. fig8d-pt-contra-bias and fig8d-mm-posterior-preferring, and the bootstrap part of fig8d-error-outnumber-correct, find the stated direction at 1.56 to 1.94 bootstrap SDs, short of the 2 SD bound the claims file fixed. The per-bin selectivity fractions are low (ALM purely lick-direction selective at most 6.2% with T1 bins, 6.5% with T2), and the rules that count neurons selective in any bin count many chance bins: with permuted trial types figs5-single-coefficient's criterion is met in every bin and the groups of fig4a-lick-direction-epochs and fig4b-object-location-memory exist, so those outcomes rest on their criteria, not on a difference from chance. The peak of figs4-peak-shift is taken over the whole window; before the go cue MM peaks first.
What was found about the dataset. Each of the 105 planes is deposited twice, under different
identifiers. ImagingPlane.location says ALM for every plane, also the MM ones. No text value in any file names a depth,
a coordinate or a layer. The plane depth, the cells' coordinates, the imaging location within a mouse and the lick times are not in the files,
nor are the wide-field imaging, the whisker tracking, the photoinhibition data and the silicon probe
recordings; the CRCNS DOI the paper gives for its data is an unfilled placeholder, and a search of DANDI finds no other deposit. The included
column marks 8,655 ROIs without saying why. After neuropil subtraction, 1,262 ROIs have a zero or
negative baseline on some trial. Early-lick trials carry no pole position. The counts of planes,
ROIs, areas, PT and IT cells and trials agree with the paper.
print(f"total runtime {(time.time() - started) / 60:.1f} min")
total runtime 16.8 min