Reproducing Economo et al. (2018) from dandiset 000006¶
Paper. Economo MN, Viswanathan S, Tasic B, et al. (2018). Distinct descending motor cortex pathways and their roles in movement. Nature 563, 79-84. doi:10.1038/s41586-018-0642-9.
Claims. The claims file claims/000006/10.1038-s41586-018-0642-9.yaml was written by a reader
from the version of record (main text, Methods, Extended Data legends and the Reporting Summary)
and checked by a second reader. It decides what is tested and how; this notebook does not change
a claim or add a method detail it lacks. Where the file leaves a detail open, it names readings,
and the notebook runs every one. Where it gives no reading for something the test needs, the
question is in questions/000006.yaml, and the notebook runs each reading it can see. A reader
answered this attempt's questions on 2026-10-20 and the claims file was updated with the answers;
this run takes them up, and each claim's section says what its answer changed.
This is the second attempt on the dandiset (000006-economo-2018-02). It replaces
000006-economo-2018-01, whose claims were written from the abstract only. The claims file was
rewritten from the full text; where a claim kept its id, the attempt record compares the outcomes.
Dandiset. DANDI:000006, version 0.220126.1855, streamed from the DANDI API. Why this version: the paper cites no DANDI version (it points to CRCNS and Figshare), and 0.220126.1855 is the only published version, so it is the most recent one and it is pinned. The draft has the same 53 files; the earlier attempt found it differs only in metadata the analysis does not use. The earlier attempt used the same version, so every difference from its results comes from the claims file, not the data.
Outside data. The Figshare record of the paper, doi:10.25378/janelia.7007846 (version 1), holds the paper's code and processed electrophysiology (Ephys.zip), the sequencing data (Seq.zip), the axon reconstructions (Axons.zip) and anatomy images (Anatomy.zip). The claims file names it as the home of the claims whose data are elsewhere, which are tested on it. The processed electrophysiology is read in the data check and in variants on the paper's own units; no outcome on the dandiset's data rests on it.
| Claim | Data | Tested |
|---|---|---|
| pt-types-transcriptomic | elsewhere (Seq.zip) | yes |
| pt-upper-thalamus-projection | elsewhere (Axons.zip) | yes |
| pt-lower-medulla-projection | elsewhere (Axons.zip) | yes |
| pt-types-shared-targets | elsewhere (Axons.zip, Anatomy.zip) | first part; second part no: blocked |
| pt-retrograde-sublaminae | elsewhere (Anatomy.zip) | no: blocked |
| pt-projection-cluster-mapping | elsewhere (Seq.zip) | yes |
| pt-types-account-for-all-pt | elsewhere (Axons.zip) | yes |
| pt-marker-genes | dataset (not public) | no: blocked |
| pt-bulk-sc-fold-change-correlation | elsewhere (Seq.zip) | yes |
| pt-tagged-units | dandiset | yes |
| pt-depth | dandiset | yes |
| pt-spike-rates-similar | dandiset | yes |
| pt-lower-heterogeneous-rates | dandiset | yes |
| pt-lower-bursting | dandiset | yes |
| pt-selective-any-epoch | dandiset | yes |
| pt-variance-explained | dandiset | yes |
| pt-preparatory-activity | dandiset | yes |
| pt-upper-early-preparatory | dandiset | yes |
| pt-upper-persistent | dandiset | yes |
| pt-lower-cdearly-decay | dandiset | yes |
| pt-cd-stability | dandiset | yes |
| pt-single-neuron-early-selectivity | dandiset | yes |
| pt-upper-decoding | dandiset | yes |
| pt-upper-decoding-rectified | dandiset | yes |
| pt-cdgo-onset | dandiset | yes |
| pt-lower-late-and-motor | dandiset | yes |
| reaction-time | dandiset | no: blocked (no video) |
| pt-cd-change-at-go-cue | dandiset | yes |
| pt-lower-movement-termination | dandiset | yes |
| pt-cohort-control | dandiset | yes |
| pt-types-differ-in-timing | dandiset (summary of six claims) | yes |
Readings shared by the electrophysiology claims (claims file, opening comment). Unit selection: (a) the code's per-unit rule, which the claims file now settles as the selection the claims use, and (b) the Methods' session rule, a variant. Trial range of the go-cue aligned rates: (a) the unit's range in the code, settled the same way, and (b) all trials, a variant. Still open, in the bootstrap: what a mouse is, (a) the 4 and 8 mice of the Methods, (b) the code's animal numbers (6 for PTlower); and whether the resampled trials enter the projected rates, (a) no, as the code, (b) yes. Outcomes are decided with selection (a) and range (a) under every open reading; the variants are printed next to them.
import io
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import alm
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
sessions = alm.load_all()
meta = alm.dandiset_metadata()
print(f"loaded {len(sessions)} files from dandiset {alm.DANDISET_ID} version {alm.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
loaded 53 files from dandiset 000006 version 0.220126.1855 in 49 s
Data check¶
Each check of the metadata against the values is stated under its output. The last cell of this section gives the status of each finding of the earlier attempt.
summary = meta["assetsSummary"]
print("Archive summary: subjects", summary["numberOfSubjects"], "; files", summary["numberOfFiles"],
"; species", [s["name"] for s in summary["species"]], "; variables", summary["variableMeasured"])
print("Related resources:", [(r.get("identifier"), r["relation"]) for r in meta.get("relatedResource", [])])
print("Description:", meta.get("description", "")[:600])
rows = []
for s in sessions:
u = s["units"]
rows.append({"subject": s["subject"], "session": s["session"], "hemisphere": s["hemisphere"],
"animal (code)": alm.animal_numbers(s), "trials": len(s["trials"]), "units": len(u),
"PTupper": int((u.cell_type == "PTupper").sum()), "PTlower": int((u.cell_type == "PTlower").sum()),
"unidentified": int((u.cell_type == "unidentified").sum()),
"nwb": s["meta"]["nwb_version"], "sex": s["meta"]["sex"], "genotype": repr(s["meta"]["genotype"]),
"session_start": s["meta"]["session_start"][:19], "paper": s["meta"]["related_publications"]})
files = pd.DataFrame(rows)
print("\nFiles", len(files), "; subjects", files.subject.nunique())
print(files.drop(columns=["paper"]).to_string(index=False))
print("\nrelated_publications in every file:", files.paper.unique().tolist())
by_subject = files.groupby("subject")[["trials", "PTupper", "PTlower", "unidentified"]].sum()
by_subject["sessions"] = files.groupby("subject").size()
by_subject["hemispheres"] = files.groupby("subject").hemisphere.agg(lambda h: "/".join(sorted(set(h))))
print("\n", by_subject.to_string())
print("\nFiles without unidentified units:", int((files.unidentified == 0).sum()))
print("Electrode x, y, z all zero in every file:", all(s["meta"]["electrodes_xyz_zero"] for s in sessions))
print("Electrode group locations:", sorted({s["location"] for s in sessions}))
print("Subjects with both PTupper and PTlower units:", int(((by_subject.PTupper > 0) & (by_subject.PTlower > 0)).sum()))
print("Units table columns in every file:", sorted({tuple(s["meta"]["units_columns"]) for s in sessions}))
print("Any acquisition object named like a video or image, in any file:", any(s["meta"]["has_video"] for s in sessions))
units_all = pd.concat([s["units"].assign(subject=s["subject"], session=s["session"]) for s in sessions])
print("\ncell_type counts:", units_all.cell_type.value_counts().to_dict())
print("Tagged label exactly when id > 1000:", bool(((units_all.id > 1000) == (units_all.cell_type != "unidentified")).all()))
print("quality values (all units):", units_all.quality.value_counts().to_dict())
print("quality values (tagged units):", units_all[units_all.cell_type != "unidentified"].quality.value_counts().to_dict())
print("\nDepth (um) by cell_type:")
print(units_all.groupby("cell_type").depth.describe()[["count", "min", "25%", "50%", "75%", "max"]].to_string())
Archive summary: subjects 12 ; files 53 ; species ['House mouse'] ; variables ['Units', 'ElectrodeGroup', 'BehavioralEvents']
Related resources: [('doi:10.1038/s41586-018-0642-9', 'dcite:IsDescribedBy')]
Description: Extracellular electrophysiology recordings performed on mouse anterior lateral motor cortex (ALM) in delay response task. Neural activity from two neuron populations, pyramidal track upper and lower, were characterized, in relation to movement execution. Some files, as originally (re)distributed from e.g. http://datasets.datalad.org/?dir=/labs/svoboda/Economo_2018 were found to be broken and would not be able available among reorganized files under sub-* directories.
Files 53 ; subjects 12
subject session hemisphere animal (code) trials units PTupper PTlower unidentified nwb sex genotype session_start
anm369962 20170309 left anm369962 274 3 0 3 0 2.0.2 U '' 2017-03-09T00:00:00
anm369962 20170310 right anm369962 right 275 79 0 5 74 2.0.2 U '' 2017-03-10T00:00:00
anm369962 20170313 right anm369962 right 354 71 0 5 66 2.0.2 U '' 2017-03-13T00:00:00
anm369962 20170314 right anm369962 right 306 51 0 4 47 2.0.2 U '' 2017-03-14T00:00:00
anm369962 20170316 right anm369962 right 252 2 0 2 0 2.0.2 U '' 2017-03-16T00:00:00
anm369962 20170317 right anm369962 right 274 59 0 2 57 2.0.2 U '' 2017-03-17T00:00:00
anm369963 20170226 left anm369963 295 49 0 2 47 2.0.2 U '' 2017-02-26T00:00:00
anm369963 20170227 left anm369963 358 42 0 2 40 2.0.2 U '' 2017-02-27T00:00:00
anm369963 20170228 left anm369963 178 1 0 1 0 2.0.2 U '' 2017-02-28T00:00:00
anm369963 20170301 left anm369963 233 3 0 3 0 2.0.2 U '' 2017-03-01T00:00:00
anm369963 20170302 left anm369963 193 2 0 2 0 2.0.2 U '' 2017-03-02T00:00:00
anm369963 20170306 right anm369963 right 257 61 0 1 60 2.0.2 U '' 2017-03-06T00:00:00
anm369963 20170309 right anm369963 right 245 2 0 2 0 2.0.2 U '' 2017-03-09T00:00:00
anm369963 20170310 right anm369963 right 273 4 0 4 0 2.0.2 U '' 2017-03-10T00:00:00
anm369964 20170320 left anm369964 319 55 0 5 50 2.0.2 U '' 2017-03-20T00:00:00
anm369964 20170321 left anm369964 386 3 0 3 0 2.0.2 U '' 2017-03-21T00:00:00
anm369964 20170322 left anm369964 288 40 0 2 38 2.0.2 U '' 2017-03-22T00:00:00
anm369964 20170323 left anm369964 459 2 0 2 0 2.0.2 U '' 2017-03-23T00:00:00
anm372793 20170504 left anm372793 304 1 1 0 0 2.0.2 U '' 2017-05-04T00:00:00
anm372793 20170508 left anm372793 340 1 1 0 0 2.0.2 U '' 2017-05-08T00:00:00
anm372793 20170512 left anm372793 282 2 2 0 0 2.0.2 U '' 2017-05-12T00:00:00
anm372793 20170513 left anm372793 368 47 1 0 46 2.0.2 U '' 2017-05-13T00:00:00
anm372793 20170514 left anm372793 270 67 3 0 64 2.0.2 U '' 2017-05-14T00:00:00
anm372794 20170621 left anm372794 248 4 0 4 0 2.0.2 U '' 2017-06-21T00:00:00
anm372794 20170622 left anm372794 340 4 0 4 0 2.0.2 U '' 2017-06-22T00:00:00
anm372794 20170624 left anm372794 248 50 0 3 47 2.0.2 U '' 2017-06-24T00:00:00
anm372794 20170625 left anm372794 323 56 0 5 51 2.0.2 U '' 2017-06-25T00:00:00
anm372794 20170626 left anm372794 328 34 0 1 33 2.0.2 U '' 2017-06-26T00:00:00
anm372794 20170627 left anm372794 385 4 0 4 0 2.0.2 U '' 2017-06-27T00:00:00
anm372795 20170714 left anm372795 219 4 4 0 0 2.0.2 U '' 2017-07-14T00:00:00
anm372795 20170715 left anm372795 203 1 1 0 0 2.0.2 U '' 2017-07-15T00:00:00
anm372795 20170716 left anm372795 273 4 4 0 0 2.0.2 U '' 2017-07-16T00:00:00
anm372795 20170718 left anm372795 263 2 2 0 0 2.0.2 U '' 2017-07-18T00:00:00
anm372797 20170615 left anm372797 224 76 3 0 73 2.0.2 U '' 2017-06-15T00:00:00
anm372797 20170617 left anm372797 234 1 1 0 0 2.0.2 U '' 2017-06-17T00:00:00
anm372904 20170615 left anm372904 345 1 1 0 0 2.0.2 U '' 2017-06-15T00:00:00
anm372904 20170616 left anm372904 314 1 1 0 0 2.0.2 U '' 2017-06-16T00:00:00
anm372904 20170617 left anm372904 304 1 1 0 0 2.0.2 U '' 2017-06-17T00:00:00
anm372904 20170618 left anm372904 334 70 1 0 69 2.0.2 U '' 2017-06-18T00:00:00
anm372904 20170619 left anm372904 254 66 1 0 65 2.0.2 U '' 2017-06-19T00:00:00
anm372905 20170715 left anm372905 224 1 1 0 0 2.0.2 U '' 2017-07-15T00:00:00
anm372905 20170716 left anm372905 277 2 2 0 0 2.0.2 U '' 2017-07-16T00:00:00
anm372905 20170717 left anm372905 291 1 1 0 0 2.0.2 U '' 2017-07-17T00:00:00
anm372906 20170608 left anm372906 306 57 4 0 53 2.0.2 U '' 2017-06-08T00:00:00
anm372906 20170610 left anm372906 315 1 1 0 0 2.0.2 U '' 2017-06-10T00:00:00
anm372906 20170611 left anm372906 227 1 1 0 0 2.0.2 U '' 2017-06-11T00:00:00
anm372906 20170612 left anm372906 301 1 1 0 0 2.0.2 U '' 2017-06-12T00:00:00
anm372907 20170608 left anm372907 319 2 2 0 0 2.0.2 U '' 2017-06-08T00:00:00
anm372907 20170610 left anm372907 310 1 1 0 0 2.0.2 U '' 2017-06-10T00:00:00
anm372907 20170612 left anm372907 312 1 1 0 0 2.0.2 U '' 2017-06-12T00:00:00
anm372907 20170613 left anm372907 139 1 1 0 0 2.0.2 U '' 2017-06-13T00:00:00
anm372909 20170520 left anm372909 241 1 1 0 0 2.0.2 U '' 2017-05-20T00:00:00
anm372909 20170522 left anm372909 265 53 1 0 52 2.0.2 U '' 2017-05-22T00:00:00
related_publications in every file: ['doi:10.1038/s41586-018-0642-9']
trials PTupper PTlower unidentified sessions hemispheres
subject
anm369962 1735 0 21 244 6 left/right
anm369963 2032 0 17 147 8 left/right
anm369964 1452 0 12 88 4 left
anm372793 1564 8 0 110 5 left
anm372794 1872 0 21 131 6 left
anm372795 958 11 0 0 4 left
anm372797 458 4 0 73 2 left
anm372904 1551 5 0 134 5 left
anm372905 792 4 0 0 3 left
anm372906 1149 7 0 53 4 left
anm372907 1080 5 0 0 4 left
anm372909 506 2 0 52 2 left
Files without unidentified units: 34
Electrode x, y, z all zero in every file: True
Electrode group locations: ['brain_region: ALM; brain_subregion: N/A; cortical_layer: 5; hemisphere: left; brain_location_full_name: N/A', 'brain_region: ALM; brain_subregion: N/A; cortical_layer: 5; hemisphere: right; brain_location_full_name: N/A']
Subjects with both PTupper and PTlower units: 0
Units table columns in every file: [('cell_type', 'depth', 'electrodes', 'electrodes_index', 'id', 'quality', 'spike_times', 'spike_times_index')]
Any acquisition object named like a video or image, in any file: False
cell_type counts: {'unidentified': 1032, 'PTlower': 71, 'PTupper': 46}
Tagged label exactly when id > 1000: True
quality values (all units): {'single': 534, '[]': 194, 'Good': 159, 'Fair': 87, 'Excellent': 80, '0': 64, 'Poor': 29, ' ': 2}
quality values (tagged units): {'Excellent': 57, 'Good': 54, 'Fair': 6}
Depth (um) by cell_type:
count min 25% 50% 75% max
cell_type
PTlower 71.0 555.0 872.5 975.0 1052.5 1600.0
PTupper 46.0 430.0 540.0 630.0 740.0 930.0
unidentified 1032.0 280.0 655.0 855.0 1015.0 1630.0
The files match the archive's summary: 53 files and 12 mice. Each mouse has one tagged type only:
PTupper in 8 mice and PTlower in 4, so the type comparison is also a comparison between animals.
The PTupper mice were recorded in the left hemisphere only, so for them the two readings of what a
mouse is coincide; two PTlower mice (anm369962 and anm369963) were also recorded in the right
hemisphere, which the code counts as separate animals. The code's grouping of animals into "MP"
(medulla-projecting) and "TP" (thalamus-projecting) in defineDataset.m puts the four PTlower mice
in "MP" and the eight PTupper mice in "TP"; the files themselves do not name the projection target.
Tagged units are exactly the units with ids above 1000. The units table holds a unit depth but no
session manipulator depth, and no file holds video. Unit quality mixes vocabularies and has empty
values, sex is U and genotype empty in every file, the electrode positions are all zero, and
every electrode group says layer 5.
trials = pd.concat([s["trials"].assign(subject=s["subject"], session=s["session"]) for s in sessions], ignore_index=True)
print("Trials", len(trials))
print(pd.crosstab(trials.response, trials.type, margins=True).to_string())
print("\nis_good:", trials.is_good.value_counts().to_dict(), "; stim_present:", trials.stim.value_counts().to_dict())
print("stop_time: all NaN in every file:", bool(trials.stop.isna().all()))
rel = pd.DataFrame({"pole_in - cue": trials.pole_in - trials.cue, "pole_out - cue": trials.pole_out - trials.cue,
"cue from start": trials.cue, "start spacing": trials.groupby(["subject", "session"]).start.diff()})
print("\nTimes in seconds (percentiles 0, 1, 50, 99, 100):")
print(rel.quantile([0, 0.01, 0.5, 0.99, 1]).round(3).T.to_string())
nxt = trials.groupby(["subject", "session"]).start.shift(-1)
print("\nTrials whose go cue (start + cue_start_time) falls after the next trial's start_time:",
int((trials.cue_abs > nxt).sum()), "of", int(nxt.notna().sum()))
odd = (trials.cue - trials.pole_out - 1.5).abs() > 0.1
print("Trials whose go cue is not 1.5 s (within 0.1 s) after pole_out:", int(odd.sum()),
"; by response:", trials[odd].response.value_counts().to_dict())
Trials 15149
type lick left lick right All
response
correct 6233 5958 12191
early lick 353 314 667
incorrect 1012 873 1885
no response 221 185 406
All 7819 7330 15149
is_good: {True: 14498, False: 651} ; stim_present: {False: 15147, True: 2}
stop_time: all NaN in every file: True
Times in seconds (percentiles 0, 1, 50, 99, 100):
0.00 0.01 0.50 0.99 1.00
pole_in - cue -88.329 -3.750 -2.501 -2.500 -2.395
pole_out - cue -87.329 -2.558 -1.500 -1.500 -1.395
cue from start 6.166 6.285 6.822 8.535 114.100
start spacing 4.918 5.400 6.400 7.900 17.870
Trials whose go cue (start + cue_start_time) falls after the next trial's start_time: 13306 of 15096
Trials whose go cue is not 1.5 s (within 0.1 s) after pole_out: 396 ; by response: {'early lick': 394, 'no response': 2}
The trial table holds start_time on the session clock and pole_in_time, pole_out_time and
cue_start_time in seconds from the trial start, as the values show though the files do not say
so. The pole comes in 2.5 s before the go cue and starts to leave 1.5 s before it in nearly all
trials, which fits the code's epochs. stop_time is NaN throughout. Consecutive trial starts are a
median 6.4 s apart, while the go cue comes 6.3 to 8.5 s after the start, so on the file's clock a
trial's go cue mostly falls after the next trial's start.
# Clock of the licks: first lick in the 2 s after the go cue, against the trial's direction,
# with the true labels and with each trial given the label of the next trial.
rows = []
for s in sessions:
t = s["trials"]
left, right = np.sort(s["licks"]["left"]), np.sort(s["licks"]["right"])
for shift in (0, 1):
labels = np.roll(t.type.to_numpy(), -shift)
for i, c in enumerate(t.cue_abs.to_numpy()):
jl, jr = np.searchsorted(left, c, "right"), np.searchsorted(right, c, "right")
fl = left[jl] if jl < len(left) else np.inf
fr = right[jr] if jr < len(right) else np.inf
first = min(fl, fr)
side = "none" if first - c > 2 else ("left" if fl < fr else "right")
rows.append({"shift": shift, "response": t.response.iloc[i],
"first lick": "none" if side == "none" else ("trial side" if side in labels[i] else "other side"),
"latency": first - c if side != "none" else np.nan})
licks = pd.DataFrame(rows)
print("First lick after the go cue, true labels (shift 0) and labels of the next trial (shift 1):")
print(pd.crosstab([licks["shift"], licks.response], licks["first lick"]).to_string())
lat = licks[(licks["shift"] == 0) & (licks.response == "correct")].latency
print(f"\nCorrect trials: first lick latency median {lat.median():.3f} s, 5th to 95th percentile "
f"{lat.quantile(0.05):.3f} to {lat.quantile(0.95):.3f} s")
# Clock of the spikes: decode the direction of correct trials from all units of a session in
# the response epoch, with true labels and with labels taken from the next trial.
rows = []
for s in sessions:
acc, n = alm.decode_direction(s, alm.EPOCHS["response"])
acc_shift, _ = alm.decode_direction(s, alm.EPOCHS["response"], shift=1)
acc_base, _ = alm.decode_direction(s, alm.EPOCHS["baseline"])
rows.append({"subject": s["subject"], "session": s["session"], "units": len(s["units"]), "trials": n,
"response": acc, "response, labels shifted": acc_shift, "baseline": acc_base,
"p (response > 0.5)": stats.binomtest(round(acc * n), n, 0.5, alternative="greater").pvalue})
decode = pd.DataFrame(rows)
print(f"\nLeave-one-out decoding of direction, median accuracy over {len(decode)} sessions: response "
f"{decode.response.median():.3f}; response with shifted labels {decode['response, labels shifted'].median():.3f}; "
f"baseline {decode.baseline.median():.3f}")
not_shown = decode[decode["p (response > 0.5)"] >= 0.01]
print("Sessions where the response-epoch decoding is not above chance (p >= 0.01):", len(not_shown))
print(not_shown[["subject", "session", "units", "response"]].round(3).to_string(index=False))
First lick after the go cue, true labels (shift 0) and labels of the next trial (shift 1):
first lick none other side trial side
shift response
0 correct 21 481 11689
early lick 18 139 510
incorrect 45 1832 8
no response 398 3 5
1 correct 21 6789 5381
early lick 18 343 306
incorrect 45 1074 766
no response 398 5 3
Correct trials: first lick latency median 0.122 s, 5th to 95th percentile 0.034 to 0.372 s
Leave-one-out decoding of direction, median accuracy over 53 sessions: response 0.893; response with shifted labels 0.552; baseline 0.529 Sessions where the response-epoch decoding is not above chance (p >= 0.01): 7 subject session units response anm369963 20170228 1 0.472 anm369964 20170323 2 0.506 anm372793 20170508 1 0.482 anm372797 20170617 1 0.496 anm372904 20170616 1 0.187 anm372906 20170612 1 0.556 anm372907 20170613 1 0.489
Licks share the trial clock: after the go cue the first lick is on the trial's side in nearly all correct trials and on the other side in nearly all incorrect trials, and the match falls to chance with the labels of the next trial. Spikes share it too: direction is decoded from the response epoch well above chance and above the baseline epoch and the shifted labels. The sessions listed, each with one or two units, could not be checked this way; their spike times equal those of the processed data (below), which puts them on the same clock.
# How the spike times and the lick times are stored.
rows = []
for s in sessions:
cue = s["trials"].cue_abs.to_numpy()
for _, u in s["units"].iterrows():
raw = np.asarray(u.spikes)
flags = alm.overlap_flags(raw)
srt = np.sort(raw)
row = {"subject": s["subject"], "session": s["session"], "id": u.id, "cell_type": u.cell_type,
"spikes": len(raw), "unsorted": bool(np.any(np.diff(raw) < 0)),
"backward steps": int(np.sum(np.diff(raw) < 0)), "trials": len(cue), "flagged": int(flags.sum())}
order = np.argsort(raw, kind="stable")
for name in ("baseline", "sample", "delay", "response"):
a, b = alm.EPOCHS[name]
inside = np.zeros(len(raw), bool)
for l, h in zip(np.searchsorted(srt, cue + a), np.searchsorted(srt, cue + b)):
inside[order[l:h]] = True
row[f"in {name}"] = int(inside.sum())
row[f"flagged in {name}"] = int((inside & flags).sum())
rows.append(row)
storage = pd.DataFrame(rows)
print("Units with spike times out of order:", int(storage.unsorted.sum()), "of", len(storage))
print("Units with more backward steps than trials:", int((storage["backward steps"] > storage.trials).sum()))
print("Spikes stored before an earlier-timed spike (flagged):", int(storage.flagged.sum()), "of", int(storage.spikes.sum()),
f"({storage.flagged.sum() / storage.spikes.sum():.2%})")
for name in ("baseline", "sample", "delay", "response"):
n, f = storage[f"in {name}"].sum(), storage[f"flagged in {name}"].sum()
print(f" spikes counted in the {name} epoch: {n}; flagged {f} ({f / n:.2%})")
lick_rows = [{"side": side, "unsorted": bool(np.any(np.diff(s["licks"][side]) < 0)),
"backward steps": int(np.sum(np.diff(s["licks"][side]) < 0)), "licks": len(s["licks"][side])}
for s in sessions for side in ("left", "right")]
lick_storage = pd.DataFrame(lick_rows)
print("\nLick time series out of order:", int(lick_storage.unsorted.sum()), "of", len(lick_storage),
"; backward steps", int(lick_storage["backward steps"].sum()), "in", int(lick_storage.licks.sum()), "licks")
Units with spike times out of order: 1120 of 1149 Units with more backward steps than trials: 0 Spikes stored before an earlier-timed spike (flagged): 370517 of 15418565 (2.40%) spikes counted in the baseline epoch: 1026782; flagged 119584 (11.65%) spikes counted in the sample epoch: 1012913; flagged 16359 (1.62%) spikes counted in the delay epoch: 3232765; flagged 37324 (1.15%) spikes counted in the response epoch: 1397243; flagged 10231 (0.73%) Lick time series out of order: 105 of 106 ; backward steps 14360 in 180679 licks
The spike times are not sorted in most units, which breaks the NWB convention for spike_times,
and no unit has more backward steps than trials, as expected if they are stored trial by trial.
The processed data below show how: trial start + go cue time + the spike's time from the go cue,
trial by trial. As consecutive trials' spike windows overlap on the file's clock, spikes of one trial
land in the next trial's early epochs; the flagged share above is a lower bound. The lick times are
stored the same way and are also out of order. The analysis uses the file's clock as published: a
trial's spikes are those within the grid around its go cue (3.4 s before to 3.0 s after), and its
last lick is rebuilt from the stored order (alm.last_licks), checked against the processed data
below.
# The paper's processed data (Ephys.zip on Figshare), every session: how the files were built, the
# code's inclusion rule on the paper's units, the last licks, and which of the paper's units the
# dandiset holds.
t0 = time.time()
members = alm.processed_members()
print("Processed data files (_psth.mat) in Ephys.zip, all read below:", len(members))
proc_units, storage_all, shifts, lick_check, proc_unidentified, lick_agree = alm.read_processed(sessions)
proc = pd.DataFrame([{k: v for k, v in r.items() if k not in ("go", "lick", "is_right", "in_range", "isi_times", "isi_trials")}
for r in proc_units])
print(f"read in {time.time() - t0:.0f} s: {len(proc)} tagged units in {proc.groupby(['subject', 'session']).ngroups} sessions")
print("\nStored spike times equal to trial start + go cue + processed time, in stored order:",
int(storage_all.equal_in_stored_order.sum()), "of", len(storage_all), "units in",
storage_all.groupby(["subject", "session"]).ngroups, "sessions; processed units not in the file:",
int((~storage_all.in_file).sum()))
print("Start of each trial's processed spike window, seconds from the pole coming in (session medians):")
print(shifts[["first 10: first - pole_in", "others: first - pole_in"]].describe().round(3).to_string())
print("Sessions where the first ten trials' windows start later than the others':",
int((shifts["first 10: first - pole_in"] > shifts["others: first - pole_in"]).sum()), "of", len(shifts))
tot = lick_check[["correct trials", "equal (within 5 ms)", "missing in files"]].sum()
print(f"\nLast lick after the go cue, rebuilt from the files against the processed trial properties: "
f"{tot['equal (within 5 ms)']} of {tot['correct trials']} correct trials equal within 5 ms "
f"({tot['equal (within 5 ms)'] / tot['correct trials']:.3f}); none found in {tot['missing in files']}")
print("Sessions with the lowest agreement:")
print(lick_check.assign(frac=lick_check["equal (within 5 ms)"] / lick_check["correct trials"])
.sort_values("frac").head(5).round(3).to_string(index=False))
print("\nThe code's inclusion rule on the processed data (the paper: 61 PTupper, 69 PTlower):")
print(pd.crosstab(proc.cell_type, proc.included, margins=True).to_string())
print("\nBy whether the session is in the dandiset:")
print(pd.crosstab([proc.cell_type, proc.in_dandiset], proc.included).to_string())
missing = proc[~proc.in_dandiset].groupby(["subject", "session"]).size()
print("Processed sessions with tagged units that are not in the dandiset:", len(missing))
print("\nUnidentified good-quality units in the processed data (the paper: 495 PTupper cohort, 511 PTlower cohort):")
print(proc_unidentified.groupby(["cohort", "in_dandiset"])["unidentified good"].agg(["sum", "count"])
.rename(columns={"sum": "units", "count": "sessions"}).to_string())
Processed data files (_psth.mat) in Ephys.zip, all read below: 14
read in 42 s: 150 tagged units in 64 sessions
Stored spike times equal to trial start + go cue + processed time, in stored order: 1149 of 1149 units in 53 sessions; processed units not in the file: 0
Start of each trial's processed spike window, seconds from the pole coming in (session medians):
first 10: first - pole_in others: first - pole_in
count 53.000 53.000
mean 2.190 -0.839
std 0.906 0.204
min 1.730 -0.971
25% 1.732 -0.970
50% 1.773 -0.928
75% 2.260 -0.811
max 5.837 -0.155
Sessions where the first ten trials' windows start later than the others': 53 of 53
Last lick after the go cue, rebuilt from the files against the processed trial properties: 11955 of 12190 correct trials equal within 5 ms (0.981); none found in 21
Sessions with the lowest agreement:
subject session correct trials equal (within 5 ms) missing in files max |difference| (s) frac
anm369963 20170309 170 158 3 3.017 0.929
anm372793 20170512 231 219 0 2.408 0.948
anm372795 20170718 160 152 4 1.282 0.950
anm372906 20170608 229 218 0 1.652 0.952
anm372793 20170508 291 278 0 1.698 0.955
The code's inclusion rule on the processed data (the paper: 61 PTupper, 69 PTlower):
included False True All
cell_type
PTlower 16 69 85
PTupper 4 61 65
All 20 130 150
By whether the session is in the dandiset:
included False True
cell_type in_dandiset
PTlower False 1 13
True 15 56
PTupper False 2 17
True 2 44
Processed sessions with tagged units that are not in the dandiset: 11
Unidentified good-quality units in the processed data (the paper: 495 PTupper cohort, 511 PTlower cohort):
units sessions
cohort in_dandiset
PTlower cohort False 117 3
True 430 24
PTupper cohort False 285 9
True 232 29
In all 53 sessions of the dandiset every unit's stored spike times equal trial start + go cue time
- the processed spike time, in the processed order. In every session the spike windows of the first ten trials start later, relative to the pole, than the others'; the code leaves out the first ten trials of each unit. The code's inclusion rule on the processed data keeps the paper's 61 PTupper and 69 PTlower units; the dandiset lacks 11 processed sessions with tagged units, and with them 17 included PTupper and 13 included PTlower units, and it also lacks the unidentified units of those sessions. The last lick rebuilt from the files matches the processed one in nearly every correct trial.
# Units of the dandiset under both readings of the unit inclusion, and the dandiset's units against
# the processed data.
t0 = time.time()
units = alm.build_units(sessions)
udf = pd.DataFrame([{k: u[k] for k in ("subject", "session", "hemisphere", "animal", "cohort", "id", "cell_type",
"quality", "depth", "included", "session_ok", "reason", "spikes_per_trial",
"n_correct_L", "n_correct_R")} for u in units])
print(f"built {len(units)} units in {time.time() - t0:.0f} s (all tagged units; unidentified units of good quality)")
tagged = udf[udf.cell_type.isin(alm.TYPES)]
print("Tagged units in the dandiset:", len(tagged), tagged.cell_type.value_counts().to_dict(), "(the paper sorted 143)")
print("\nReading (a), the code's per-unit rule:")
print(pd.crosstab(tagged.cell_type, tagged.included).to_string())
for t in alm.TYPES:
reasons = tagged[(tagged.cell_type == t) & ~tagged.included].reason.str.split(", ").explode()
print(" ", t, "left out for:", reasons.value_counts().to_dict())
print("\nReading (b), the Methods' session rule:")
print(pd.crosstab(tagged.cell_type, tagged.session_ok).to_string())
both = proc.merge(tagged, on=["subject", "session", "id"], suffixes=(" processed", " dandiset"))
print(f"\nTagged units in both the dandiset and the processed data: {len(both)} of {len(tagged)}; same cell type "
f"{int((both['cell_type processed'] == both['cell_type dandiset']).sum())}; same inclusion "
f"{int((both['included processed'] == both['included dandiset']).sum())}; same depth "
f"{int((both['depth processed'] == both['depth dandiset']).sum())}")
for coh in ("PTupper cohort", "PTlower cohort"):
g = udf[(udf.cell_type == "unidentified") & (udf.cohort == coh)]
print(f"Unidentified good-quality units, {coh}: {len(g)} in {g.groupby(['subject', 'session']).ngroups} sessions; "
f"reading (a) {int(g.included.sum())}, reading (b) {int(g.session_ok.sum())}; left hemisphere, reading (a) "
f"{int((g.included & (g.hemisphere == 'left')).sum())}")
print("\nIncluded tagged units by animal (reading (a)):")
print(tagged[tagged.included].groupby(["cell_type", "animal"]).size().to_string())
built 779 units in 3 s (all tagged units; unidentified units of good quality)
Tagged units in the dandiset: 117 {'PTlower': 71, 'PTupper': 46} (the paper sorted 143)
Reading (a), the code's per-unit rule:
included False True
cell_type
PTlower 15 56
PTupper 2 44
PTupper left out for: {'performance': 1, 'too few correct trials': 1}
PTlower left out for: {'low rate': 7, 'quality': 6, 'performance': 2, 'too few correct trials': 1}
Reading (b), the Methods' session rule:
session_ok False True
cell_type
PTlower 8 63
PTupper 4 42
Tagged units in both the dandiset and the processed data: 117 of 117; same cell type 117; same inclusion 117; same depth 117
Unidentified good-quality units, PTupper cohort: 232 in 6 sessions; reading (a) 226, reading (b) 232; left hemisphere, reading (a) 226
Unidentified good-quality units, PTlower cohort: 430 in 12 sessions; reading (a) 403, reading (b) 430; left hemisphere, reading (a) 212
Included tagged units by animal (reading (a)):
cell_type animal
PTlower anm369962 3
anm369962 right 17
anm369963 8
anm369963 right 5
anm369964 9
anm372794 14
PTupper anm372793 8
anm372795 11
anm372797 3
anm372904 5
anm372905 4
anm372906 7
anm372907 4
anm372909 2
Neither reading of the unit inclusion gives the paper's counts on the dandiset: the code's rule keeps 44 PTupper and 56 PTlower units, the session rule 42 and 63. The unidentified units are fewer than the paper's too: 232 of good quality in 6 sessions of the PTupper cohort (the paper: 495 in 15) and 430 in 12 sessions of the PTlower cohort (511 in 14). Every tagged unit the dandiset holds is in the processed data, with the same cell type, inclusion decision and depth. Neither reading matches, and the claims file settles which selection the other claims use: the code's per-unit rule, which gives the paper's 61 and 69 on the processed data. The session rule is a variant, shown next to it. Units recorded in one animal on different days may be the same neuron; the files hold no waveforms or identifier that would show it, so units are counted per session, as in the code.
# Status of each finding of the earlier attempt (000006-economo-2018-01), checked again here.
status = pd.DataFrame([
{"finding": "1. data_missing: processed sessions with tagged units not in the dandiset",
"check": f"{len(missing)} sessions; included units missing PTupper "
f"{int((~proc.in_dandiset & proc.included & (proc.cell_type == 'PTupper')).sum())}, PTlower "
f"{int((~proc.in_dandiset & proc.included & (proc.cell_type == 'PTlower')).sum())}",
"status": "still present" if len(missing) else "fixed"},
{"finding": "2. data_missing: no transcriptomic or anatomical data in the dandiset",
"check": f"variables measured {summary['variableMeasured']}",
"status": "still present" if all("RNA" not in v and "Image" not in v for v in summary["variableMeasured"]) else "check"},
{"finding": "3. format_error: spike times not sorted, trials overlap",
"check": f"{int(storage.unsorted.sum())} of {len(storage)} units unsorted; {storage.flagged.sum() / storage.spikes.sum():.2%} of spikes flagged",
"status": "still present" if storage.unsorted.any() else "fixed"},
{"finding": "4. format_error: stop_time NaN; relative trial times undocumented; go cue after next start",
"check": f"stop_time all NaN {bool(trials.stop.isna().all())}; go cue after next start in {int((trials.cue_abs > nxt).sum())} trials",
"status": "still present" if trials.stop.isna().all() else "fixed"},
{"finding": "5. metadata_insufficient: PTupper and PTlower not described",
"check": f"cell_type descriptions in the files: {sorted({s['meta']['cell_type_description'] for s in sessions})}",
"status": "still present"},
{"finding": "6. metadata_insufficient: quality vocabularies, sex U, genotype empty, positions zero, layer 5",
"check": f"quality values {sorted(units_all.quality.unique().tolist())}; sex {sorted(files.sex.unique().tolist())}; "
f"positions zero {all(s['meta']['electrodes_xyz_zero'] for s in sessions)}",
"status": "still present"},
{"finding": "7. other: first ten trials' spike windows shifted in the processed data",
"check": f"{int((shifts['first 10: first - pole_in'] > shifts['others: first - pole_in']).sum())} of {len(shifts)} sessions",
"status": "still present"},
])
print(status.to_string(index=False))
finding check status
1. data_missing: processed sessions with tagged units not in the dandiset 11 sessions; included units missing PTupper 17, PTlower 13 still present
2. data_missing: no transcriptomic or anatomical data in the dandiset variables measured ['Units', 'ElectrodeGroup', 'BehavioralEvents'] still present
3. format_error: spike times not sorted, trials overlap 1120 of 1149 units unsorted; 2.40% of spikes flagged still present
4. format_error: stop_time NaN; relative trial times undocumented; go cue after next start stop_time all NaN True; go cue after next start in 13306 trials still present
5. metadata_insufficient: PTupper and PTlower not described cell_type descriptions in the files: ['cell type (e.g. wide width, narrow width spiking)'] still present
6. metadata_insufficient: quality vocabularies, sex U, genotype empty, positions zero, layer 5 quality values [' ', '0', 'Excellent', 'Fair', 'Good', 'Poor', '[]', 'single']; sex ['U']; positions zero True still present
7. other: first ten trials' spike windows shifted in the processed data 53 of 53 sessions still present
The version is the one the earlier attempt read, and nothing in it has changed: all seven findings
of 000006-economo-2018-01 are still present, as the table shows. Finding 5 is shown by the
cell_type description, which is generic. This attempt adds findings: the unidentified units of the
missing sessions are missing too; the lick times are stored trial by trial and out of order, as the
spike times are; the files hold no session manipulator depth and no video.
Per-unit measures and the shared bootstrap¶
Every per-unit measure is computed once, under both readings of the trial range, from the smoothed rates (5 ms bins, causal 50 ms boxcar). The hierarchical bootstrap (1,000 iterations, 50 trials of each type per unit) is run for each population under each reading of the unit selection, the trial range and what a mouse is; each run gives both readings of the projection. Outcomes use selection (a) and range (a); the other runs are variants. Each run has its own seeded generator. The same is done on the paper's 61 and 69 units from the processed data, as a variant.
t0 = time.time()
rows = []
for i, u in enumerate(units):
for tr in ("unit", "all"):
rows.append({"index": i, "range": tr, **{k: u[k] for k in ("subject", "session", "animal", "cohort", "hemisphere",
"id", "cell_type", "quality", "depth", "included",
"session_ok")},
"burst fraction": alm.isi_burst_fraction(u["raw_spikes"], np.diff(u["raw_spikes"]) < 0),
"burst fraction, sorted": alm.isi_burst_fraction(np.sort(u["raw_spikes"]), np.zeros(max(len(u["raw_spikes"]) - 1, 0), bool)),
**alm.unit_measures(u, tr)})
per_unit = pd.DataFrame(rows)
rows = []
for i, u in enumerate(proc_units):
rows.append({"index": i, "range": "unit", **{k: u[k] for k in ("subject", "session", "animal", "id", "cell_type", "depth",
"manipulator_depth", "included", "in_dandiset")},
"burst fraction": alm.isi_burst_fraction(u["isi_times"], np.diff(u["isi_trials"]) != 0),
**alm.unit_measures(u, "unit")})
per_unit_proc = pd.DataFrame(rows)
print(f"per-unit measures in {time.time() - t0:.0f} s")
SELECT = {"code": "included", "methods": "session_ok"}
READINGS = [(s, r) for s in ("code", "methods") for r in ("unit", "all")]
def pop_table(cell_types, selection, trial_range, frame=None):
# Units of the given types under a reading of the selection and of the trial range.
t = per_unit if frame is None else frame
t = t[(t.range == trial_range) & t.cell_type.isin(cell_types)]
return t[t[SELECT[selection]]] if frame is None else t[t.included]
for s, r in READINGS:
t = pop_table(alm.TYPES, s, r)
print(f"selection {s}, trial range {r}: {t.cell_type.value_counts().to_dict()}; trials per unit (left) "
f"{int(t.n_left.min())} to {int(t.n_left.max())}, (right) {int(t.n_right.min())} to {int(t.n_right.max())}")
per-unit measures in 10 s
selection code, trial range unit: {'PTlower': 56, 'PTupper': 44}; trials per unit (left) 62 to 203, (right) 56 to 199
selection code, trial range all: {'PTlower': 56, 'PTupper': 44}; trials per unit (left) 69 to 218, (right) 63 to 211
selection methods, trial range unit: {'PTlower': 63, 'PTupper': 42}; trials per unit (left) 35 to 209, (right) 37 to 202
selection methods, trial range all: {'PTlower': 63, 'PTupper': 42}; trials per unit (left) 69 to 218, (right) 63 to 211
t0 = time.time()
boots = {}
for ct in alm.TYPES:
for sel in ("code", "methods"):
pop = alm.selected(units, ct, sel)
for tr in ("unit", "all"):
for mk in ("subject", "animal"):
if ct == "PTupper" and mk == "animal":
boots[(ct, sel, tr, mk)] = boots[(ct, sel, tr, "subject")] # the same animals
continue
key = (ct, sel, tr, mk)
boots[key] = alm.summarise(alm.bootstrap(pop, mk, tr, alm.seeded("bootstrap " + " ".join(key)), lick=True))
boots_proc = {}
for ct in alm.TYPES:
pop = [u for u in proc_units if u["included"] and u["cell_type"] == ct]
for mk in ("subject", "animal"):
boots_proc[(ct, mk)] = alm.summarise(alm.bootstrap(pop, mk, "unit", alm.seeded(f"bootstrap processed {ct} {mk}"), lick=True))
print(f"bootstraps in {time.time() - t0:.0f} s: {len(boots)} on the dandiset, {len(boots_proc)} on the processed data")
BOOT_KEYS = [(sel, tr, mk, pr) for sel in ("code", "methods") for tr in ("unit", "all")
for mk in ("subject", "animal") for pr in ("proj a", "proj b")]
def sel_summary(ct, sel, tr, mk, pr, cd, align="go", processed=False):
return (boots_proc[(ct, mk)] if processed else boots[(ct, sel, tr, mk)])[(align, pr, cd)]
def reading_label(sel, tr, mk, pr):
return f"sel {sel}, range {tr}, mouse {mk}, {pr}"
bootstraps in 71 s: 16 on the dandiset, 4 on the processed data
# Variant: the processed data restricted to the tagged units the dandiset also holds, under the code's
# rule; it has the dandiset's units and the paper's per-trial rates (section "Missing units against the
# file format").
t0 = time.time()
per_unit_same = per_unit_proc[per_unit_proc.in_dandiset]
boots_same = {}
for ct in alm.TYPES:
pop = [u for u in proc_units if u["included"] and u["in_dandiset"] and u["cell_type"] == ct]
boots_same[ct] = alm.summarise(alm.bootstrap(pop, "subject", "unit", alm.seeded(f"bootstrap processed same units {ct}")))
print(f"processed data, units the dandiset holds: {per_unit_same[per_unit_same.included].cell_type.value_counts().to_dict()}; "
f"bootstraps in {time.time() - t0:.0f} s")
processed data, units the dandiset holds: {'PTlower': 56, 'PTupper': 44}; bootstraps in 5 s
# Selectivity along each coding direction, go-cue aligned (main reading), with the 5-95% band.
fig, axes = plt.subplots(1, 3, figsize=(13, 3.4), sharex=True)
colours = {"PTupper": "#1b9e77", "PTlower": "#d95f02"}
ranges = []
for ax, cd in zip(axes, ("early", "late", "go")):
for ct in alm.TYPES:
s = sel_summary(ct, "code", "unit", "subject", "proj a", cd)
ax.fill_between(alm.TIME, s["q5"], s["q95"], color=colours[ct], alpha=0.25, lw=0)
ax.plot(alm.TIME, s["q50"], color=colours[ct], lw=1.5, label=f"{ct} ({alm.PROJECTION[ct]})")
ranges.append((cd, ct, float(s["q5"].min()), float(s["q95"].max())))
for x in (-2.5, -1.3, 0):
ax.axvline(x, color="#999999", lw=0.8)
ax.axhline(0, color="#444444", lw=0.8, ls=":")
ax.set_title(f"CD{cd}")
ax.set_xlabel("Time from go cue (s)")
axes[0].set_ylabel("Selectivity (a.u.)")
axes[0].legend(frameon=False, fontsize=8)
for ax in axes:
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
plt.show()
print("Plotted ranges (5th to 95th percentile), all within the axes:")
for cd, ct, lo, hi in ranges:
print(f" CD{cd} {ct}: {lo:.2f} to {hi:.2f}")
Plotted ranges (5th to 95th percentile), all within the axes: CDearly PTupper: -2.20 to 9.36 CDearly PTlower: -5.55 to 13.82 CDlate PTupper: -3.75 to 13.38 CDlate PTlower: -2.16 to 22.55 CDgo PTupper: -1.68 to 10.27 CDgo PTlower: -2.96 to 27.70
Claims with data elsewhere: sequencing and axon reconstructions¶
These claims rest on data the dandiset does not hold. The claims file names the paper's Figshare record (doi:10.25378/janelia.7007846, version 1) as their home; Seq.zip and Axons.zip are read from it in memory. They are counted apart from the claims on the dandiset.
t0 = time.time()
seq = alm.read_seq()
anno = seq["anno"]
print(f"Seq.zip read in {time.time() - t0:.0f} s: annotation {anno.shape}, single-cell expression {seq['sc'].shape} "
f"(genes x cells), bulk {seq['bulk'].shape} (genes x samples)")
print("Regions:", anno.region_label.value_counts().to_dict(), "; classes:", anno.class_label.value_counts().to_dict())
print("cluster_metadata.csv: clusters", len(seq["cluster_metadata"]), "; labels of its PT-like clusters:",
seq["cluster_metadata"][seq["cluster_metadata"].cluster_label.str.contains("L5b ALM")].cluster_label.tolist())
neurons = alm.read_axons()
print("Axons.zip: objects", len(neurons), "; axons", sum(n.Name.endswith("-A") for n in neurons),
"; dendrites", sum(n.Name.endswith("-D") for n in neurons))
anatomy = alm.figshare_zip(alm.FIGSHARE_ANATOMY_ZIP).namelist()
folders = pd.Series([n.split("/")[1] for n in anatomy if n.count("/") >= 2 and not n.startswith("__MACOSX")]).value_counts()
print("Anatomy.zip folders (files):", folders.to_dict())
print("File types in Anatomy.zip:", pd.Series([n.rsplit(".", 1)[-1] for n in anatomy if "." in n.split("/")[-1]]).value_counts().to_dict())
Seq.zip read in 5 s: annotation (24411, 152), single-cell expression (45768, 368) (genes x cells), bulk (48526, 22) (genes x samples)
Regions: {'VISp': 14739, 'ALM': 9672} ; classes: {'Glutamatergic': 11905, 'GABAergic': 10534, 'Non-Neuronal': 1383, 'Low Quality': 589}
cluster_metadata.csv: clusters 142 ; labels of its PT-like clusters: ['L5b ALM Kcnip1', 'L5b ALM Pappa2_1', 'L5b ALM Pappa2_2']
Axons.zip: objects 41 ; axons 22 ; dendrites 19
Anatomy.zip folders (files): {'BL777_Medulla_ChR2-YFP_Ch01': 122, 'BL772_ThalMed_anterograde_Ch01': 111, 'BL772_ThalMed_anterograde_Ch02': 111, 'BL648_Thal_Med_H2b_Ch01': 109, 'BL648_Thal_Med_H2b_Ch02': 109, 'BL780_Thalamus_ChR2-YFP_Ch01': 101}
File types in Anatomy.zip: {'jpg': 657}
Seq.zip holds the full single-cell annotation of the Tasic et al. (2018) taxonomy (cluster, class,
region and retrograde injection site per cell), the single-cell expression of the PT cells and the
bulk RNA-seq counts. Its cluster_metadata.csv is from an earlier version of the taxonomy (other
cluster names), so the annotation is used. Axons.zip holds 22 reconstructed neurons, as axon and
dendrite objects, with the structure of every node. Anatomy.zip holds only image stacks of six
brains: no cell counts, no co-labelling annotations and no smFISH data.
pt-types-transcriptomic¶
In ALM scRNA-seq, glutamatergic neurons form 27 clusters and PT neurons map to three of them (Slco2a1, and the closely related Npsr1 and Hpgd), and with axonal reconstructions this suggests that PT neurons can be divided into two cell types.
Criterion: the cluster count is exact (tolerance 0). Counting rule (claims file, from a reader's answer): Fig. 1a names 26 glutamatergic clusters and leaves out the Cajal-Retzius cluster, and the Methods give 27, so the count is of the ALM glutamatergic clusters not named for VISp, the Cajal-Retzius cluster included. The other counts are variants.
alm_cells = anno[anno.region_label == "ALM"]
glu = alm_cells[alm_cells.class_label == "Glutamatergic"]
clusters = glu.cluster_label.value_counts()
share = anno[anno.class_label == "Glutamatergic"].groupby("cluster_label").region_label.agg(lambda r: (r == "ALM").mean())
readings = {
"criterion: clusters with an ALM cell not named for VISp, Cajal-Retzius included": [c for c in clusters.index if "VISp" not in c],
"variant: as the criterion, without the Cajal-Retzius cluster (the clusters of Fig. 1a)": [c for c in clusters.index if "VISp" not in c and not c.startswith("CR")],
"variant: every glutamatergic cluster with an ALM cell": clusters.index,
"variant: glutamatergic clusters with most of their cells in ALM": share[share > 0.5].index,
}
for name, cl in readings.items():
print(f"{name}: {len(cl)}")
print("\nGlutamatergic ALM cells:", len(glu), "; cells per cluster under the criterion's rule:")
print(clusters[[c for c in clusters.index if "VISp" not in c]].to_string())
pt_cells = alm_cells[alm_cells.subclass_label == "L5 PT"]
print("\nALM cells of the L5 PT subclass by cluster:", pt_cells.cluster_label.value_counts().to_dict())
criterion: clusters with an ALM cell not named for VISp, Cajal-Retzius included: 27
variant: as the criterion, without the Cajal-Retzius cluster (the clusters of Fig. 1a): 26
variant: every glutamatergic cluster with an ALM cell: 45
variant: glutamatergic clusters with most of their cells in ALM: 25
Glutamatergic ALM cells: 4539 ; cells per cluster under the criterion's rule:
cluster_label
L5 IT ALM Tnc 704
L5 IT ALM Pld5 412
L5 IT ALM Cbln4 Fezf2 387
L5 NP ALM Trhr Nefl 294
L6 CT ALM Cpa6 258
L6 NP ALM Trh 241
L5 IT ALM Gkn1 Pcdh19 230
L5 IT ALM Tmem163 Dmrtb1 224
L6 IT ALM Tgfb1 215
L5 PT ALM Slco2a1 203
L5 IT ALM Lypd1 Gpr88 184
L6 IT ALM Oprk1 172
L2/3 IT ALM Sla 124
L5 PT ALM Hpgd 122
L5 IT ALM Cpa6 Gpr88 121
L2/3 IT ALM Macc1 Lrg1 119
L6 CT Nxph2 Sla 85
L5 IT ALM Tmem163 Arhgap25 72
L2/3 IT ALM Ptrf 70
L5 IT ALM Npw 53
L6b ALM Olfr111 Nxph1 51
L6b ALM Olfr111 Spon1 46
L5 PT ALM Npsr1 43
L6b Col8a1 Rprm 30
CR Lhx5 11
L6b Hsd17b2 9
L6b P2ry12 8
ALM cells of the L5 PT subclass by cluster: {'L5 PT ALM Slco2a1': 203, 'L5 PT ALM Hpgd': 122, 'L5 PT ALM Npsr1': 43}
Result: reproduced. Under the criterion's rule the ALM glutamatergic cells of the annotation fall in 27 clusters, the paper's number. The variants give 26 without the Cajal-Retzius cluster (the clusters of Fig. 1a), 45 counting every glutamatergic cluster with an ALM cell, and 25 counting the clusters with most of their cells in ALM. The 368 ALM cells of the L5 PT subclass map to three clusters: Slco2a1 (203), Hpgd (122) and Npsr1 (43). That the reconstructions suggest two types is an interpretation and is not tested. The answer fixed the counting rule, which the earlier run had taken as a reading of the reported number; the outcome is unchanged.
pt-upper-thalamus-projection, pt-lower-medulla-projection, pt-types-shared-targets and pt-types-account-for-all-pt¶
pt-upper-thalamus-projection: Of 12 reconstructed PT neurons, 8 innervate the thalamus.
pt-lower-medulla-projection: The other 4 of 12 reconstructed PT neurons bypass the thalamus and branch in the medulla, including the reticular nuclei with premotor neurons for orofacial movements.
pt-types-shared-targets: Both PT types project to several shared targets; all reconstructed cells innervate the superior colliculus, and both types send collaterals to the subthalamic nucleus, zona incerta and midbrain reticular nucleus.
pt-types-account-for-all-pt: All 18 partially and fully reconstructed PT neurons project to the superior colliculus and none lacks projections to both thalamus and medulla, suggesting that most if not all PT neurons belong to the two types.
Criteria (claims file, from a reader's answers). The 12 neurons and their groups are those of Supplementary data table 1, the same as the hand grouping of AxonReconstructions.m. Innervating the thalamus is read as (a) any axon terminus there, (b) any axon length there, (c) more termini in the thalamus than in the medulla; bypassing it as (a) no terminus, (b) no axon, (c) more termini in the medulla than in the thalamus. The 8 and the 4 are reported counts, so each holds if it matches under one reading. pt-types-shared-targets: the first part (all cells reach the superior colliculus) is stated of the reconstructions and tested on them; the second part is stated in the Extended Data Fig. 2e legend of the bulk labelling, its intended test is those images, and the paper gives no method to analyse them, so the test on the reconstructions is a variant. pt-types-account-for-all-pt: the 18 neurons listed in Supplementary data table 1, with projecting read as (a) any terminus, (b) any axon length. Termini and axon length per structure are counted as AxonReconstructions.m counts them.
structures = ["Thalamus", "Medulla", "Intermediate reticular nucleus", "Superior colliculus motor related",
"Superior colliculus sensory related", "Subthalamic nucleus", "Zona incerta", "Midbrain reticular nucleus"]
axons = alm.axon_targets(neurons, structures)
cols = ["neuron", "soma", "group"] + [c for c in axons.columns if c.endswith("termini")]
print("Axon termini per structure (and its subdivisions):")
print(axons[cols].to_string(index=False))
print("\nAxon length (mm):")
print(axons[["neuron", "group"] + [c for c in axons.columns if c.endswith("length")]].round(1).to_string(index=False))
axons["SC termini"] = axons["Superior colliculus motor related termini"] + axons["Superior colliculus sensory related termini"]
axons["SC length"] = axons["Superior colliculus motor related length"] + axons["Superior colliculus sensory related length"]
twelve = axons[axons.group != ""]
print("\nThe code's 12 neurons:", len(twelve), twelve.group.value_counts().to_dict())
Axon termini per structure (and its subdivisions):
neuron soma group Thalamus termini Medulla termini Intermediate reticular nucleus termini Superior colliculus motor related termini Superior colliculus sensory related termini Subthalamic nucleus termini Zona incerta termini Midbrain reticular nucleus termini
AA0066 ACAd5 7 0 0 0 0 0 66 5
AA0121 ORBl5 11 0 0 17 0 0 0 9
AA0114 MOp5 thalamus 15 10 0 80 0 8 1 37
AA0131 MOp5 1 3 0 8 0 0 1 7
AA0132 MOp5 1 0 0 18 0 0 0 24
AA0133 MOp6a medulla 0 83 15 3 0 0 5 4
AA0135 MOp5 4 0 0 3 0 0 0 2
AA0169 MOp6a 0 16 13 10 0 0 0 4
AA0261 MOp5 thalamus 89 0 0 28 0 4 2 59
AA0119 AId5 14 0 0 32 0 0 0 31
AA0134 MOp6a 0 13 6 0 0 1 0 3
AA0001 SSp-m5 13 0 0 0 0 0 0 4
AA0011 MOs5 thalamus 51 7 0 51 0 9 2 53
AA0115 MOs5 thalamus 36 3 0 17 0 0 0 32
AA0122 MOs6a thalamus 39 0 0 26 0 4 0 27
AA0181 MOs5 thalamus 61 5 0 6 0 15 1 9
AA0182 MOs5 thalamus 33 8 0 22 0 0 0 14
AA0245 MOs5 thalamus 75 1 0 17 0 13 1 25
AA0250 MOs5 33 15 0 15 0 1 0 85
AA0012 MOs6a medulla 1 117 48 19 0 1 2 17
AA0179 MOs6a medulla 7 52 1 17 0 2 13 9
AA0180 MOs6a medulla 1 130 13 3 0 2 0 9
Axon length (mm):
neuron group Thalamus length Medulla length Intermediate reticular nucleus length Superior colliculus motor related length Superior colliculus sensory related length Subthalamic nucleus length Zona incerta length Midbrain reticular nucleus length
AA0066 3.0 0.0 0.0 0.0 0.0 0.0 13.9 3.1
AA0121 5.8 0.0 0.0 4.8 0.0 0.3 0.4 10.0
AA0114 thalamus 10.6 2.3 0.0 23.8 0.0 1.5 0.3 23.4
AA0131 0.2 5.6 0.4 3.3 0.0 0.3 0.5 4.1
AA0132 2.2 0.0 0.0 5.3 0.0 0.0 0.0 12.9
AA0133 medulla 0.0 30.2 6.7 1.3 0.0 0.2 0.5 5.7
AA0135 1.6 0.0 0.0 1.4 0.0 0.0 0.0 3.6
AA0169 0.0 10.0 5.3 4.0 0.0 0.0 0.0 10.9
AA0261 thalamus 25.3 0.1 0.0 5.5 0.0 2.3 1.4 15.9
AA0119 10.2 0.0 0.0 14.0 0.0 0.1 0.9 20.9
AA0134 0.5 16.6 4.2 0.0 0.0 0.3 0.0 3.5
AA0001 6.9 0.0 0.0 0.0 0.0 0.0 0.0 2.4
AA0011 thalamus 21.0 1.7 0.0 18.7 0.0 3.0 1.5 21.8
AA0115 thalamus 17.8 1.5 0.0 4.5 0.0 0.0 1.2 13.7
AA0122 thalamus 29.9 0.2 0.0 8.0 0.0 1.2 0.7 10.8
AA0181 thalamus 21.7 1.0 0.0 3.6 0.0 2.5 0.6 10.3
AA0182 thalamus 19.6 5.0 0.2 11.4 0.0 0.0 0.0 13.7
AA0245 thalamus 37.5 0.2 0.0 5.3 0.0 3.9 1.5 7.5
AA0250 14.4 5.6 0.1 9.4 0.0 0.1 0.3 32.2
AA0012 medulla 0.3 37.9 12.7 6.8 0.0 1.5 0.6 6.4
AA0179 medulla 4.7 29.2 2.1 12.0 0.0 1.1 5.8 8.6
AA0180 medulla 0.3 51.0 6.6 0.8 0.0 0.6 0.0 13.0
The code's 12 neurons: 12 {'thalamus': 8, 'medulla': 4}
def innervates(frame, s, reading):
if reading == "termini":
return frame[f"{s} termini"] > 0
return frame[f"{s} length"] > 0
print("pt-upper-thalamus-projection and pt-lower-medulla-projection (the 12 neurons of Supplementary data table 1):")
print("Same neurons and groups as AxonReconstructions.m:", sorted(twelve.neuron) == sorted(alm.AXON_GROUPS))
dominant = twelve["Thalamus termini"] > twelve["Medulla termini"]
for name, th, me in (("(a) any terminus", innervates(twelve, "Thalamus", "termini"), innervates(twelve, "Medulla", "termini")),
("(b) any length", innervates(twelve, "Thalamus", "length"), innervates(twelve, "Medulla", "length")),
("(c) more termini than in the other", dominant, twelve["Medulla termini"] > twelve["Thalamus termini"])):
print(f" {name}: innervate the thalamus {int(th.sum())} of 12 (the paper 8; thalamus group {int(th[twelve.group == 'thalamus'].sum())} of 8); "
f"bypass the thalamus {int((~th).sum())} of 12 (the paper 4; medulla group {int((~th)[twelve.group == 'medulla'].sum())} of 4); "
f"medulla group in the medulla {int(me[twelve.group == 'medulla'].sum())} of 4")
print(" (c) the 8 are the thalamus group:", bool((dominant == (twelve.group == "thalamus")).all()))
for reading in ("termini", "length"):
print(f" medulla group in the intermediate reticular nucleus ({reading}): "
f"{int(innervates(twelve[twelve.group == 'medulla'], 'Intermediate reticular nucleus', reading).sum())} of 4")
print(" Medulla group, thalamic termini:", twelve[twelve.group == "medulla"]["Thalamus termini"].tolist(),
"; thalamus group, medullary termini:", twelve[twelve.group == "thalamus"]["Medulla termini"].tolist())
print("\npt-types-shared-targets (the 12 neurons):")
for reading in ("termini", "length"):
sc = innervates(twelve, "SC", reading)
shared = {s: {g: int(innervates(twelve[twelve.group == g], s, reading).sum()) for g in ("thalamus", "medulla")}
for s in ("Subthalamic nucleus", "Zona incerta", "Midbrain reticular nucleus")}
print(f" {reading}: superior colliculus {int(sc.sum())} of 12; variant on the reconstructions, neurons of each group reaching:", shared)
listed = {**alm.AXON_GROUPS, **alm.AXON_PARTIAL}
eighteen = axons[axons.neuron.isin(listed)]
print("\npt-types-account-for-all-pt: the 18 neurons of Supplementary data table 1; found in Axons.zip:", len(eighteen),
"; the same as the 18 with somata in motor cortex:", sorted(eighteen.neuron) == sorted(axons[axons.soma.str.startswith("MO")].neuron))
print(" Neurons of Axons.zip not among the 18:", axons[~axons.neuron.isin(listed)][["neuron", "soma"]].to_dict("records"))
for name, frame in (("the 18", eighteen), ("variant, all 22", axons)):
for reading in ("termini", "length"):
sc = innervates(frame, "SC", reading)
neither = ~innervates(frame, "Thalamus", reading) & ~innervates(frame, "Medulla", reading)
print(f" {name}, n = {len(frame)}, ({'a' if reading == 'termini' else 'b'}) {reading}: superior colliculus {int(sc.sum())} of {len(frame)} "
f"(not: {frame[~sc].neuron.tolist()}); lacking both thalamus and medulla {int(neither.sum())} "
f"({frame[neither].neuron.tolist()})")
a134 = axons[axons.neuron == "AA0134"].iloc[0]
print(" AA0134 (Supplementary data table 1, medulla-projecting, partial): soma", a134.soma, "; superior colliculus termini",
int(a134["SC termini"]), "length", round(a134["SC length"], 2), "mm; thalamus termini", int(a134["Thalamus termini"]),
"; medulla termini", int(a134["Medulla termini"]), "; medulla length", round(a134["Medulla length"], 1), "mm")
pt-upper-thalamus-projection and pt-lower-medulla-projection (the 12 neurons of Supplementary data table 1):
Same neurons and groups as AxonReconstructions.m: True
(a) any terminus: innervate the thalamus 11 of 12 (the paper 8; thalamus group 8 of 8); bypass the thalamus 1 of 12 (the paper 4; medulla group 1 of 4); medulla group in the medulla 4 of 4
(b) any length: innervate the thalamus 11 of 12 (the paper 8; thalamus group 8 of 8); bypass the thalamus 1 of 12 (the paper 4; medulla group 1 of 4); medulla group in the medulla 4 of 4
(c) more termini than in the other: innervate the thalamus 8 of 12 (the paper 8; thalamus group 8 of 8); bypass the thalamus 4 of 12 (the paper 4; medulla group 4 of 4); medulla group in the medulla 4 of 4
(c) the 8 are the thalamus group: True
medulla group in the intermediate reticular nucleus (termini): 4 of 4
medulla group in the intermediate reticular nucleus (length): 4 of 4
Medulla group, thalamic termini: [0, 1, 7, 1] ; thalamus group, medullary termini: [10, 0, 7, 3, 0, 5, 8, 1]
pt-types-shared-targets (the 12 neurons):
termini: superior colliculus 12 of 12; variant on the reconstructions, neurons of each group reaching: {'Subthalamic nucleus': {'thalamus': 6, 'medulla': 3}, 'Zona incerta': {'thalamus': 5, 'medulla': 3}, 'Midbrain reticular nucleus': {'thalamus': 8, 'medulla': 4}}
length: superior colliculus 12 of 12; variant on the reconstructions, neurons of each group reaching: {'Subthalamic nucleus': {'thalamus': 6, 'medulla': 4}, 'Zona incerta': {'thalamus': 7, 'medulla': 3}, 'Midbrain reticular nucleus': {'thalamus': 8, 'medulla': 4}}
pt-types-account-for-all-pt: the 18 neurons of Supplementary data table 1; found in Axons.zip: 18 ; the same as the 18 with somata in motor cortex: True
Neurons of Axons.zip not among the 18: [{'neuron': 'AA0066', 'soma': 'ACAd5'}, {'neuron': 'AA0121', 'soma': 'ORBl5'}, {'neuron': 'AA0119', 'soma': 'AId5'}, {'neuron': 'AA0001', 'soma': 'SSp-m5'}]
the 18, n = 18, (a) termini: superior colliculus 17 of 18 (not: ['AA0134']); lacking both thalamus and medulla 0 ([])
the 18, n = 18, (b) length: superior colliculus 17 of 18 (not: ['AA0134']); lacking both thalamus and medulla 0 ([])
variant, all 22, n = 22, (a) termini: superior colliculus 19 of 22 (not: ['AA0066', 'AA0134', 'AA0001']); lacking both thalamus and medulla 0 ([])
variant, all 22, n = 22, (b) length: superior colliculus 19 of 22 (not: ['AA0066', 'AA0134', 'AA0001']); lacking both thalamus and medulla 0 ([])
AA0134 (Supplementary data table 1, medulla-projecting, partial): soma MOp6a ; superior colliculus termini 0 length 0.0 mm; thalamus termini 0 ; medulla termini 13 ; medulla length 16.6 mm
# Where the axon of AA0134 runs: the structures of its axon nodes (the 25 most frequent) and of its termini.
a134 = next(n for n in neurons if n.Name == "AA0134-A")
node_names = pd.Series([str(a.safeName) for a in np.atleast_1d(a134.Anatomy)])
ends134 = ~np.isin(a134.Nodes[:, 4].astype(int), a134.Nodes[:, 3].astype(int))
print("AA0134 axon nodes:", len(node_names), "; termini:", int(ends134.sum()))
print("Nodes by structure:", node_names.value_counts().head(25).to_dict())
print("Termini by structure:", node_names[ends134].value_counts().to_dict())
print("Any node in a structure named for the superior colliculus:", bool(node_names.str.contains("colliculus", case=False).any()))
AA0134 axon nodes: 2368 ; termini: 36
Nodes by structure: {'Intermediate reticular nucleus': 235, 'medial longitudinal fascicle': 209, 'pyramid': 188, 'Gigantocellular reticular nucleus': 149, 'Midbrain reticular nucleus': 141, 'medial lemniscus': 138, 'Parvicellular reticular nucleus': 121, 'Pontine reticular nucleus caudal part': 113, 'Pontine gray': 100, 'Medullary reticular nucleus': 95, 'Medulla': 72, 'Magnocellular reticular nucleus': 66, 'Red nucleus': 66, 'tectospinal pathway': 55, 'cerebal peduncle': 50, 'corticospinal tract': 48, 'Inferior olivary complex': 47, 'Nucleus raphe magnus': 44, 'Caudoputamen': 43, 'Superior central nucleus raphe': 38, 'Nucleus raphe pallidus': 31, 'internal capsule': 29, 'fiber tracts': 28, 'Substantia nigra reticular part': 26, 'facial nerve': 24}
Termini by structure: {'Intermediate reticular nucleus': 6, 'medial longitudinal fascicle': 4, 'Parvicellular reticular nucleus': 4, 'Pontine gray': 4, 'Midbrain reticular nucleus': 3, 'medial lemniscus': 3, 'Tegmental reticular nucleus': 2, 'Medullary reticular nucleus': 1, 'Medulla': 1, 'Gigantocellular reticular nucleus': 1, 'Pontine reticular nucleus caudal part': 1, 'tectospinal pathway': 1, 'facial nerve': 1, 'Pons': 1, 'Superior central nucleus raphe': 1, 'Red nucleus': 1, 'Subthalamic nucleus': 1}
Any node in a structure named for the superior colliculus: False
Result: reproduced. Under reading (c), 8 of the 12 neurons have more termini in the thalamus than in the medulla, the paper's 8, and they are exactly the thalamus group of Supplementary data table 1. Readings (a) and (b) give 11 of 12, as three of the four medulla-group neurons have a few thalamic termini (1, 7 and 1) and thalamic axon. The 8 is a reported count, so it holds under the reading that matches. The answer named the readings and the neurons; the earlier run had the same counts and was blocked, as the claims file then gave no rule.
Result: reproduced. Under reading (c), 4 of the 12 neurons have more termini in the medulla than in the thalamus, the paper's 4, all of them the medulla group, and all 4 have termini and axon in the medulla and in the intermediate reticular nucleus. Under readings (a) and (b) only 1 of the 4 lacks thalamic termini and axon. Six of the eight thalamus-group neurons also have medullary termini (1 to 10). The 4 is a reported count, so it holds under the reading that matches. The answer named the readings; the earlier run was blocked, as the claims file then gave no rule.
Result: blocked (methods underspecified). The first part holds: all 12 neurons have axon and termini in the superior colliculus, all of it in the motor-related part. The second part's intended test is the bulk labelling images of Extended Data Fig. 2e, for which the paper gives no method of analysis, so it is not tested. On the reconstructions, a variant, neurons of both groups reach the subthalamic nucleus (6 thalamus-group and 3 medulla-group neurons by termini, 6 and 4 by axon length), the zona incerta (5 and 3; 7 and 3) and the midbrain reticular nucleus (8 and 4 under both). The answer made the test on the reconstructions a variant; the earlier run had the claim reproduced on it.
Result: partially reproduced. The 18 neurons of Supplementary data table 1 are all in Axons.zip, and they are the 18 with somata in motor cortex that the earlier run used. None lacks both thalamus and medulla (0 of 18), but 17 of 18 reach the superior colliculus: AA0134, a partial medulla-projecting neuron of the table, has neither termini nor axon there (it has 13 medullary termini). Both readings of projecting (termini, axon length) agree. Over all 22 neurons of the file, a variant, 19 reach it. The answer fixed which 18 neurons are meant; the outcome is unchanged. The axon of AA0134 (second cell) has 2368 nodes and 36 termini, and no node lies in a structure named for the superior colliculus; it runs through the midbrain reticular nucleus (141 nodes) and the tectospinal pathway (55) to the intermediate reticular nucleus (235 nodes, 6 termini) and other reticular nuclei of the pons and medulla.
pt-projection-cluster-mapping¶
All 62 medulla-projecting PT neurons map to the Slco2a1 cluster, and 91 of 94 thalamus-projecting PT neurons map to the Npsr1 and Hpgd clusters.
Criterion (claims file, from a reader's answer): the paper does not say which annotation field
defines the injection groups of Fig. 2f. Readings, (a) the grouping of the paper's code for the same
cells (SC_correlations_across_clusters.ipynb), which matches inj_roi_label to "IRN" (medulla),
"VAL" (thalamus), "PG" (pons) and "SC|SCm" (superior colliculus); (b) inj_primary_label grouped
with the Allen ontology (medulla IRN, GRN, PARN; thalamus VAL, PF, MD; pons PG, PRNc; superior
colliculus SC, SCm). The counts (62 exact; 91 +/- 9.1) are reported numbers, so the claim holds
when both match under one reading; the group sizes are compared with the Methods' n (62, 94, 100,
97). PT neurons are the ALM cells of the L5 PT subclass.
site = pt_cells.inj_roi_label.astype(str)
print("Reading (a), the code's grouping (SC_correlations_across_clusters.ipynb, grep on inj_roi_label):")
for g, pattern in (("medulla", "IRN"), ("thalamus", "VAL"), ("pons", "PG"), ("superior colliculus", "SC|SCm")):
cells = pt_cells[site.str.contains(pattern)]
counts = cells.cluster_label.str.replace("L5 PT ALM ", "", regex=False).value_counts().to_dict()
print(f" {g} ({pattern}): {len(cells)} cells; {counts}")
groups = {"medulla": {"IRN", "GRN", "PARN"}, "thalamus": {"VAL", "PF", "MD"}, "superior colliculus": {"SC", "SCm"},
"pons": {"PG", "PRNc"}}
site = pt_cells.inj_primary_label.astype(str).str.replace("ALM : ", "", regex=False)
print("Reading (b), inj_primary_label with the Allen ontology:")
for g, members in groups.items():
cells = pt_cells[site.isin(members)]
counts = cells.cluster_label.str.replace("L5 PT ALM ", "", regex=False).value_counts().to_dict()
print(f" {g}: {len(cells)} cells; {counts}")
print("\nThe paper: medulla 62, all Slco2a1; thalamus 94, Npsr1 19 and Hpgd 72 (91); superior colliculus 97 "
"(Slco2a1 43, Npsr1 15, Hpgd 39); pons 100 (Slco2a1 86, Npsr1 5, Hpgd 9)")
Reading (a), the code's grouping (SC_correlations_across_clusters.ipynb, grep on inj_roi_label):
medulla (IRN): 58 cells; {'Slco2a1': 58}
thalamus (VAL): 92 cells; {'Hpgd': 71, 'Npsr1': 18, 'Slco2a1': 3}
pons (PG): 99 cells; {'Slco2a1': 85, 'Hpgd': 9, 'Npsr1': 5}
superior colliculus (SC|SCm): 97 cells; {'Slco2a1': 43, 'Hpgd': 39, 'Npsr1': 15}
Reading (b), inj_primary_label with the Allen ontology:
medulla: 62 cells; {'Slco2a1': 62}
thalamus: 94 cells; {'Hpgd': 72, 'Npsr1': 19, 'Slco2a1': 3}
superior colliculus: 92 cells; {'Slco2a1': 43, 'Hpgd': 36, 'Npsr1': 13}
pons: 100 cells; {'Slco2a1': 86, 'Hpgd': 9, 'Npsr1': 5}
The paper: medulla 62, all Slco2a1; thalamus 94, Npsr1 19 and Hpgd 72 (91); superior colliculus 97 (Slco2a1 43, Npsr1 15, Hpgd 39); pons 100 (Slco2a1 86, Npsr1 5, Hpgd 9)
Result: reproduced. Under reading (b), all 62 medulla-labelled PT neurons are in Slco2a1 and 91 of the 94 thalamus-labelled ones are in Npsr1 (19) or Hpgd (72), the paper's numbers. Its pons group matches the Methods too (100: 86, 5 and 9), and its superior colliculus group has 92 cells against 97. Under reading (a), the code's grouping, the medulla group has 58 cells, all Slco2a1, and 89 of the 92 thalamus-labelled cells are in Npsr1 (18) or Hpgd (71), so the medulla count does not match; its superior colliculus group matches (97: 43, 15 and 39) and its pons group has 99 cells. The answer named the code's grouping as a reading; the outcome is unchanged and rests on reading (b).
pt-bulk-sc-fold-change-correlation¶
Fold changes between PTupper and PTlower of the 200 most discriminative genes correlate between scRNA-seq and bulk RNA-seq (Pearson's r = 0.87).
Criterion (claims file, from a reader's answer): fold change is the mean log2(CPM + 1) of the Npsr1 and Hpgd cells minus that of the Slco2a1 cells, in single cells and in bulk, as the Extended Data Fig. 4i legend says; the 200 genes (100 at each end of the binary detection accuracy in single cells) and the cell and replicate groups (six thalamus against six medulla samples) are those of bulkVsSC.m on Seq.zip. The code applies log2(x + 1) to the bulk counts as stored, a variant. Reproduced when r is within 0.87 +/- 0.087.
print("Sample key:", seq["key"].label.tolist())
print("Single-cell values stored as CPM (every cell sums to 1e6):", bool(np.allclose(seq["sc"].sum(axis=0), 1e6)))
for scale, name in (("cpm", "criterion, log2(CPM + 1) as the legend"), ("counts", "variant, bulk values as stored, as bulkVsSC.m")):
r = alm.bulk_vs_sc(seq, scale)
print(f"{name}: r = {r['r']:.3f} (p = {r['p']:.2g}) over {r['kept']} genes of {r['genes']} in both; "
f"r over all those genes {r['r all genes']:.3f}; single cells {r['cells thal']} Npsr1/Hpgd and {r['cells med']} Slco2a1")
Sample key: ['SL1-GFP_Thal', 'SL1-tdtomato_IRN', 'SL1-GFP_Thal_Outlier', 'SL1-tdtomato_IRN ', 'SL1-tdtomato_IRN', 'SL1-GFP_Thal', 'SL1-tdtomato_IRN', 'SL1-GFP_Thal', 'SL1-GFP_Thal', 'SL1-tdtomato_IRN', 'SL1-GFP_Thal', 'SL1-tdtomato_IRN ', 'SL1-GFP_Thal', 'SL1-tdtomato_IRN', 'Mouse_qRNA', 'Mouse_qRNA', 'SC', 'SC', 'SC', 'Pons', 'Pons', 'Pons'] Single-cell values stored as CPM (every cell sums to 1e6): True
criterion, log2(CPM + 1) as the legend: r = 0.874 (p = 4.1e-64) over 200 genes of 31002 in both; r over all those genes 0.516; single cells 165 Npsr1/Hpgd and 203 Slco2a1
variant, bulk values as stored, as bulkVsSC.m: r = 0.867 (p = 8.1e-62) over 200 genes of 31002 in both; r over all those genes 0.458; single cells 165 Npsr1/Hpgd and 203 Slco2a1
Result: reproduced. With the bulk samples as counts per million, as the legend says (the single-cell values are stored as CPM), r = 0.874 over the 200 genes, within 0.87 +/- 0.087. The code's variant, with the bulk counts as stored, gives 0.867. Over all 31,002 genes in both data sets r is 0.516 (0.458 in the variant), so the agreement lies in the discriminative genes. The answer made the legend's scale the main reading, where the earlier run led with the code's; the outcome is unchanged.
pt-retrograde-sublaminae and pt-marker-genes (blocked)¶
pt-retrograde-sublaminae: Thalamus-projecting PT neurons lie in upper L5b and medulla-projecting PT neurons in lower L5b; few neurons are co-labelled from thalamus and medulla (2.2%, 22 of 984) while many are co-labelled from thalamus and superior colliculus (77.1%, 687 of 890).
pt-marker-genes: smFISH confirms Slco2a1 as a PTlower marker and Npnt as a PTupper marker.
Anatomy.zip (listed in the data check above) holds image stacks of six brains (thalamus and
medulla H2B labelling, anterograde labelling, ChR2 expression) without cell counts or
co-labelling annotations. Neither the paper nor its code says how labelled and co-labelled cells
were counted (claims file, from a reader's answer), so the claim stays methods underspecified
unless the counts can be shown not to depend on it. pt-marker-genes is now data: dataset
(claims file, from a reader's answer): the paper names archives only for the bulk RNA-seq and the
electrophysiology, and no file of the dandiset or of the Figshare record holds smFISH images or
puncta counts.
print("Anatomy.zip folders:", sorted(folders.index.tolist()))
print("Any folder named for the superior colliculus:", any("sc" in f.lower().replace("_", " ").split() or "colliculus" in f.lower() for f in folders.index))
print("Any folder or file named for smFISH, RNAscope, Npnt or Slco2a1:",
any(k in n.lower() for n in anatomy for k in ("fish", "rnascope", "npnt", "slco2a1")))
print("Any table (csv, xls, txt, mat) in Anatomy.zip:", [n for n in anatomy if n.lower().endswith((".csv", ".xls", ".xlsx", ".txt", ".mat"))])
Anatomy.zip folders: ['BL648_Thal_Med_H2b_Ch01', 'BL648_Thal_Med_H2b_Ch02', 'BL772_ThalMed_anterograde_Ch01', 'BL772_ThalMed_anterograde_Ch02', 'BL777_Medulla_ChR2-YFP_Ch01', 'BL780_Thalamus_ChR2-YFP_Ch01'] Any folder named for the superior colliculus: False Any folder or file named for smFISH, RNAscope, Npnt or Slco2a1: False Any table (csv, xls, txt, mat) in Anatomy.zip: []
# The image stacks of Anatomy.zip, from the archive's directory (no image is opened): images per folder
# and their compressed sizes. The files hold images only, with no scale or cell annotation beside them.
z_anatomy = alm.figshare_zip(alm.FIGSHARE_ANATOMY_ZIP)
for folder in sorted(folders.index):
jpgs = sorted(n for n in anatomy if n.split("/")[1:2] == [folder] and n.lower().endswith(".jpg"))
mb = np.array([z_anatomy.getinfo(n).file_size for n in jpgs]) / 1e6
print(f"{folder}: {len(jpgs)} images, {mb.min():.2f} to {mb.max():.2f} MB each; first {jpgs[0].split('/')[-1]}")
BL648_Thal_Med_H2b_Ch01: 108 images, 0.78 to 1.84 MB each; first 648_jSV101_H2B_SL1_3xPT_jpgs0000.jpg BL648_Thal_Med_H2b_Ch02: 108 images, 1.00 to 2.07 MB each; first BL648_Thal_Med_H2b_Ch020000.jpg BL772_ThalMed_anterograde_Ch01: 110 images, 0.81 to 1.85 MB each; first 772_jSV108_172124_1002_PSjpgs0001.jpg BL772_ThalMed_anterograde_Ch02: 110 images, 0.63 to 2.12 MB each; first 772_jSV108_172124_1002_PSjpgs0001.jpg BL777_Medulla_ChR2-YFP_Ch01: 121 images, 0.75 to 1.81 MB each; first 777_ME15_ALM_med_WR43_172x124_1013_PSjpgs0000.jpg BL780_Thalamus_ChR2-YFP_Ch01: 100 images, 0.95 to 1.93 MB each; first 780_ME16_ALM_med_WR42_172x124_1025_PSjpgs0000.jpg
Result: blocked. pt-retrograde-sublaminae (methods underspecified): the only data are image stacks (657 jpg files, no tables), no folder is named for superior colliculus labelling, which the 77.1% needs, and neither the paper nor its code gives a method for counting labelled and co-labelled cells. The answer confirmed this; the outcome is unchanged. pt-marker-genes (data missing): no smFISH images or puncta counts are in the dandiset or in any file of the Figshare record, and the paper names no other archive for them. The answer changed the claim to data dataset, so it is blocked as data missing, where the earlier run used other. The two H2B folders of brain BL648, its thalamus and medulla retrograde labelling, hold 108 images each (0.78 to 2.07 MB); nothing beside them gives a scale, cell positions or counts, so the co-labelling count would need a cell detection method, which the paper does not give. Opened outside the notebook while preparing (not run here, and deciding nothing), one image of each folder is 8600 by 6200 pixels, greyscale.
pt-tagged-units¶
The recordings yield 61 identified PTupper cells from 8 mice and 69 identified PTlower neurons from 4 mice after exclusions (130 of 143 tagged units kept).
Criterion: reproduced when one reading of the unit selection gives both unit counts within tolerance (61 +/- 6.1, 69 +/- 6.9) and the mouse counts exactly (8, 4). When neither reading matches, the other claims use the code's per-unit rule (claims file, from a reader's answer).
rows = []
for sel, col in SELECT.items():
for ct in alm.TYPES:
t = tagged[(tagged.cell_type == ct) & tagged[col]]
rows.append({"reading": sel, "type": ct, "units": len(t), "mice": t.subject.nunique()})
counts = pd.DataFrame(rows)
print(counts.to_string(index=False))
print("Tagged units before exclusion:", len(tagged), "(the paper: 143)")
ok = {sel: bool((abs(counts[(counts.reading == sel) & (counts.type == "PTupper")].units.iloc[0] - 61) <= 6.1)
and (abs(counts[(counts.reading == sel) & (counts.type == "PTlower")].units.iloc[0] - 69) <= 6.9)
and counts[(counts.reading == sel) & (counts.type == "PTupper")].mice.iloc[0] == 8
and counts[(counts.reading == sel) & (counts.type == "PTlower")].mice.iloc[0] == 4) for sel in SELECT}
print("Within tolerance under the reading:", ok)
print("Processed data, the code's rule:", proc[proc.included].cell_type.value_counts().to_dict(),
"; mice", proc[proc.included].groupby("cell_type").subject.nunique().to_dict())
reading type units mice
code PTupper 44 8
code PTlower 56 4
methods PTupper 42 8
methods PTlower 63 4
Tagged units before exclusion: 117 (the paper: 143)
Within tolerance under the reading: {'code': False, 'methods': False}
Processed data, the code's rule: {'PTlower': 69, 'PTupper': 61} ; mice {'PTlower': 4, 'PTupper': 8}
# The processed data hold more tagged units than the 143 the paper sorted: by quality and the code's rule.
print("Tagged units in the processed data:", len(proc), "(the paper: 143 tagged units, 130 kept)")
print(pd.crosstab([proc.cell_type, proc.quality], proc.included, margins=True).to_string())
print("Of good quality (Excellent, Good, single):", int(proc.quality.isin(alm.GOOD_QUALITY).sum()),
"; reasons for leaving units out:", proc[~proc.included].reason.str.split(", ").explode().value_counts().to_dict())
Tagged units in the processed data: 150 (the paper: 143 tagged units, 130 kept)
included False True All
cell_type quality
PTlower Excellent 3 47 50
Fair 7 0 7
Good 6 22 28
PTupper Excellent 2 22 24
Good 2 39 41
All 20 130 150
Of good quality (Excellent, Good, single): 143 ; reasons for leaving units out: {'quality': 7, 'low rate': 7, 'too few correct trials': 4, 'performance': 3}
Result: not reproduced. The dandiset holds 117 tagged units (the paper sorted 143). The code's rule keeps 44 PTupper and 56 PTlower units, and the session rule keeps 42 and 63, from 8 and 4 mice under both. Only the session rule's 63 PTlower units are within tolerance (69 +/- 6.9); no reading gives both counts. On the processed data the code's rule gives 61 and 69 units from 8 and 4 mice. The difference is the 11 sessions missing from the dandiset (data check). Following the answer, the other claims use the code's rule, with 44 PTupper and 56 PTlower units, and the session rule is a variant. The outcome is unchanged. The processed data hold 150 tagged units (second cell): 143 of good quality (Excellent or Good), the paper's 143, and 7 PTlower units of Fair quality, which the code's quality filter leaves out; of the 143 the code's rule keeps 130, the paper's number.
pt-depth¶
Identified PTupper cells are at more superficial recording depths than PTlower neurons.
Criterion: reproduced when the mean depth of included PTupper units is smaller than that of PTlower units, judged on each unit's depth field, which Fig. 4e plots (claims file, from a reader's answer); the session manipulator depth is a variant, shown where available. Shown with all units and without the right-hemisphere PTlower recordings.
for sel in SELECT:
t = pop_table(alm.TYPES, sel, "unit")
for label, frame in (("all units", t), ("without right-hemisphere PTlower", t[~((t.cell_type == "PTlower") & (t.hemisphere == "right"))])):
d = frame.groupby("cell_type").depth.agg(["mean", "std", "count"]).round(1)
print(f"selection {sel}, {label}: PTupper {d.loc['PTupper', 'mean']} +/- {d.loc['PTupper', 'std']} um (n = {int(d.loc['PTupper', 'count'])}); "
f"PTlower {d.loc['PTlower', 'mean']} +/- {d.loc['PTlower', 'std']} um (n = {int(d.loc['PTlower', 'count'])}); "
f"PTupper shallower: {d.loc['PTupper', 'mean'] < d.loc['PTlower', 'mean']}")
print("\nThe files hold no session manipulator depth (units columns in the data check). Variant on the processed data:")
pp = per_unit_proc[per_unit_proc.included]
for label, col in (("unit depth field", "depth"), ("manipulator depth", "manipulator_depth")):
for sub, frame in (("all units", pp), ("without right-hemisphere PTlower", pp[~pp.animal.isin(["29.5", "30.5"])])):
d = frame.groupby("cell_type")[col].mean().round(1)
print(f" {label}, {sub}: PTupper {d['PTupper']} um, PTlower {d['PTlower']} um")
selection code, all units: PTupper 636.0 +/- 132.2 um (n = 44); PTlower 1004.1 +/- 217.3 um (n = 56); PTupper shallower: True selection code, without right-hemisphere PTlower: PTupper 636.0 +/- 132.2 um (n = 44); PTlower 891.6 +/- 120.8 um (n = 34); PTupper shallower: True selection methods, all units: PTupper 627.0 +/- 128.4 um (n = 42); PTlower 994.0 +/- 202.7 um (n = 63); PTupper shallower: True selection methods, without right-hemisphere PTlower: PTupper 627.0 +/- 128.4 um (n = 42); PTlower 901.1 +/- 122.2 um (n = 40); PTupper shallower: True The files hold no session manipulator depth (units columns in the data check). Variant on the processed data: unit depth field, all units: PTupper 631.6 um, PTlower 1037.9 um unit depth field, without right-hemisphere PTlower: PTupper 631.6 um, PTlower 895.1 um manipulator depth, all units: PTupper 1117.1 um, PTlower 1373.8 um manipulator depth, without right-hemisphere PTlower: PTupper 1117.1 um, PTlower 1188.9 um
Result: reproduced. With the code's selection, included PTupper units lie at a mean unit depth of 636 +/- 132 um and PTlower units at 1004 +/- 217 um, and at 892 um without the right-hemisphere recordings. The session rule, a variant, gives 627 and 994 um (901 um). The files hold no session manipulator depth; on the processed data, a variant, it orders the types the same way (1117 against 1374 um, and 1189 um without the right hemisphere). The answer made the unit depth decide; the outcome is unchanged.
pt-spike-rates-similar¶
Baseline and trial-averaged peak spike rates do not differ between PTupper (baseline 4.3 +/- 3.5 Hz, peak 17.7 +/- 13.8 Hz, median +/- s.d.) and PTlower cells (5.2 +/- 4.8 Hz, 17.2 +/- 22 Hz).
Criterion: a claim of absence, judged on two-sided Student's t-tests between the populations, for baseline and peak. Reproduced when neither test has p < 0.05 and the four medians are within 10%; partially reproduced when the tests hold and a median is out of tolerance; not reproduced when either test has p < 0.05. Selection (a) and range (a) decide; the rest are variants.
expected = {("PTupper", "baseline rate"): 4.3, ("PTupper", "peak rate"): 17.7, ("PTlower", "baseline rate"): 5.2, ("PTlower", "peak rate"): 17.2}
def rate_rows(frame, label):
out = []
for m in ("baseline rate", "peak rate"):
up, low = frame[frame.cell_type == "PTupper"][m], frame[frame.cell_type == "PTlower"][m]
p = stats.ttest_ind(up, low).pvalue
out.append({"reading": label, "measure": m, "PTupper median": up.median(), "PTupper sd": up.std(),
"PTlower median": low.median(), "PTlower sd": low.std(), "t-test p": p,
"medians in tolerance": all(abs(x.median() - expected[(ct, m)]) <= 0.1 * expected[(ct, m)]
for ct, x in (("PTupper", up), ("PTlower", low)))})
return out
rows = []
for s, r in READINGS:
rows += rate_rows(pop_table(alm.TYPES, s, r), f"sel {s}, range {r}")
rows += rate_rows(pop_table(alm.TYPES, None, "unit", per_unit_proc), "processed data (variant)")
rates_table = pd.DataFrame(rows)
print(rates_table.round(3).to_string(index=False))
reading measure PTupper median PTupper sd PTlower median PTlower sd t-test p medians in tolerance
sel code, range unit baseline rate 4.933 3.672 6.689 5.849 0.074 False
sel code, range unit peak rate 18.369 11.177 20.958 22.916 0.030 False
sel code, range all baseline rate 4.877 3.560 6.593 5.750 0.069 False
sel code, range all peak rate 17.737 10.372 19.629 21.894 0.028 False
sel methods, range unit baseline rate 4.862 3.755 5.513 5.975 0.435 False
sel methods, range unit peak rate 18.369 11.169 15.172 22.564 0.356 False
sel methods, range all baseline rate 4.678 3.641 5.363 5.891 0.440 True
sel methods, range all peak rate 17.737 10.358 14.375 21.552 0.336 False
processed data (variant) baseline rate 4.265 3.528 5.182 4.827 0.178 True
processed data (variant) peak rate 17.647 13.818 17.167 22.044 0.142 True
Result: not reproduced. With the code's selection and the unit's range, peak rates differ between the types (two-sided Student's t-test p = 0.030), so the claim of absence fails; baseline rates do not differ (p = 0.074), and neither p exceeds the paper's bound of 0.1. The medians are 4.93 Hz (baseline) and 18.37 Hz (peak) for PTupper and 6.69 and 20.96 Hz for PTlower; only the PTupper peak is within tolerance. The s.d. are 3.67, 11.18, 5.85 and 22.92 Hz. With all trials, a variant, p = 0.069 and 0.028. The session rule, a variant, gives p = 0.34 to 0.44, with medians out of tolerance. On the paper's units from the processed data all four medians are within tolerance (4.27, 17.65, 5.18 and 17.17 Hz) and neither test finds a difference (p = 0.18 and 0.14). The answer fixed the selection; the earlier run was blocked because the two selections gave different outcomes. On the processed data for the same 44 and 56 units the dandiset holds (section "Missing units against the file format"), peak rates also differ (p = 0.03) and baseline rates do not (p = 0.085), so the difference from the paper's units comes from the missing sessions, not from the files. The files raise every unit's baseline rate (median 0.99 Hz for PTupper and 0.97 Hz for PTlower), as spikes of neighbouring trials fall in the baseline epoch.
pt-lower-heterogeneous-rates¶
Spike rates are more heterogeneous across the PTlower population than the PTupper population, for baseline (P = 0.02) and peak rates (P = 3 x 10^-4).
Criterion: rates as defined in pt-spike-rates-similar; the claim depends on pt-tagged-units, not on that claim's outcome (claims file, from a reader's answer); readings (a) Bartlett's test, (b) the code's two-sided F test (vartest2); the outcome is the one they give. Reproduced when, for both baseline and peak, the PTlower variance is larger and the test has p < 0.05; partially reproduced when that holds for one.
def ftest(a, b):
f = np.var(a, ddof=1) / np.var(b, ddof=1)
cdf = stats.f.cdf(f, len(a) - 1, len(b) - 1)
return 2 * min(cdf, 1 - cdf)
rows = []
frames = [(f"sel {s}, range {r}", pop_table(alm.TYPES, s, r)) for s, r in READINGS]
frames.append(("processed data (variant)", pop_table(alm.TYPES, None, "unit", per_unit_proc)))
for label, frame in frames:
for m in ("baseline rate", "peak rate"):
up, low = frame[frame.cell_type == "PTupper"][m], frame[frame.cell_type == "PTlower"][m]
rows.append({"reading": label, "measure": m, "PTupper var": up.var(), "PTlower var": low.var(),
"PTlower larger": low.var() > up.var(), "Bartlett p": stats.bartlett(low, up).pvalue,
"F test p": ftest(low, up)})
print(pd.DataFrame(rows).round(4).to_string(index=False))
reading measure PTupper var PTlower var PTlower larger Bartlett p F test p
sel code, range unit baseline rate 13.4825 34.2147 True 0.0019 0.0020
sel code, range unit peak rate 124.9289 525.1627 True 0.0000 0.0000
sel code, range all baseline rate 12.6723 33.0650 True 0.0014 0.0015
sel code, range all peak rate 107.5823 479.3607 True 0.0000 0.0000
sel methods, range unit baseline rate 14.1015 35.6970 True 0.0019 0.0021
sel methods, range unit peak rate 124.7425 509.1443 True 0.0000 0.0000
sel methods, range all baseline rate 13.2597 34.7022 True 0.0014 0.0014
sel methods, range all peak rate 107.2792 464.4954 True 0.0000 0.0000
processed data (variant) baseline rate 12.4481 23.3021 True 0.0140 0.0142
processed data (variant) peak rate 190.9359 485.9449 True 0.0003 0.0003
Result: reproduced. With the code's selection and the unit's range, PTlower rates are more variable than PTupper rates for baseline (variance 34.2 against 13.5 Hz^2; Bartlett p = 0.0019, F test p = 0.0020) and for peak (525 against 125 Hz^2; both p < 0.0001), under both tests. The variants (all trials, the session rule) give the same (baseline p 0.0014 to 0.0021, peak p < 0.0001). On the processed data the tests give p = 0.014 and 0.0003, near the paper's 0.02 and 3 x 10^-4. The answer moved the dependency from pt-spike-rates-similar, which blocked the claim in the earlier run, to pt-tagged-units. pt-tagged-units is not reproduced on the dandiset, as 11 sessions are missing, so the test runs on 44 and 56 units, not the paper's 61 and 69.
pt-lower-bursting¶
A substantial proportion of PTlower neurons fire bursts (18.8%), which are rare among PTupper cells (3.3%).
Criterion: a bursting cell has more than 10% of its inter-spike intervals shorter than 5 ms, intervals across a trial boundary counted as long, over all its spikes. Reproduced when both percentages are within 10% and the two-sided Fisher exact test has p < 0.05; partially reproduced when the test holds and a percentage is out of tolerance. The files do not mark trial boundaries; in the stored order (trial by trial) a boundary is taken where the time steps backwards, and the sorted spike times are a variant.
rows = []
frames = [(f"sel {s}", pop_table(alm.TYPES, s, "unit"), "burst fraction") for s in SELECT]
frames += [(f"sel {s}, sorted spike times", pop_table(alm.TYPES, s, "unit"), "burst fraction, sorted") for s in SELECT]
frames.append(("processed data (variant)", pop_table(alm.TYPES, None, "unit", per_unit_proc), "burst fraction"))
for label, frame, col in frames:
b = frame[col] > alm.BURST_FRACTION
up, low = b[frame.cell_type == "PTupper"], b[frame.cell_type == "PTlower"]
p = stats.fisher_exact([[low.sum(), (~low).sum()], [up.sum(), (~up).sum()]]).pvalue
rows.append({"reading": label, "PTupper bursting": f"{up.sum()} of {len(up)}", "PTupper %": 100 * up.mean(),
"PTlower bursting": f"{low.sum()} of {len(low)}", "PTlower %": 100 * low.mean(), "Fisher p": p,
"in tolerance": abs(100 * up.mean() - 3.3) <= 0.33 and abs(100 * low.mean() - 18.8) <= 1.88})
print(pd.DataFrame(rows).round(4).to_string(index=False))
f = pop_table(alm.TYPES, "code", "unit")
print("\nBurst fraction (% of intervals under 5 ms) by type, reading (a):")
print((100 * f.groupby("cell_type")["burst fraction"].describe()[["min", "50%", "max"]]).round(2).to_string())
reading PTupper bursting PTupper % PTlower bursting PTlower % Fisher p in tolerance
sel code 2 of 44 4.5455 11 of 56 19.6429 0.0350 False
sel methods 2 of 42 4.7619 12 of 63 19.0476 0.0419 False
sel code, sorted spike times 2 of 44 4.5455 11 of 56 19.6429 0.0350 False
sel methods, sorted spike times 2 of 42 4.7619 12 of 63 19.0476 0.0419 False
processed data (variant) 2 of 61 3.2787 13 of 69 18.8406 0.0058 True
Burst fraction (% of intervals under 5 ms) by type, reading (a):
min 50% max
cell_type
PTlower 0.0 0.75 39.76
PTupper 0.0 0.43 17.35
Result: partially reproduced. With the code's selection, 11 of 56 PTlower units (19.6%) and 2 of 44 PTupper units (4.5%) burst (two-sided Fisher p = 0.035). The PTlower share is within 18.8 +/- 1.88% and the test holds, but the PTupper share is above 3.3 +/- 0.33%, as 2 units of 44 are 4.5%. Sorting the spike times gives the same counts; the session rule, a variant, gives 12 of 63 (19.0%) and 2 of 42 (4.8%), p = 0.042. On the paper's units the shares are 13 of 69 (18.8%) and 2 of 61 (3.3%), p = 0.0058: the paper's numbers. The answer fixed the selection; the outcome is unchanged. The processed data for the same units give the same counts, 2 of 44 and 11 of 56, so the PTupper share is out of tolerance because units are missing: with 44 units one unit is 2.3%, and 3.3 +/- 0.33% cannot be reached.
pt-selective-any-epoch¶
Most recorded PT neurons (122 of 130, 94%) have significant trial-type selectivity in at least one task epoch.
Criterion: per unit and epoch, the mean rate of each correct trial, two-sided Mann-Whitney test at p < 0.01; a unit counts when any of the sample (-2.5 to -2.1 s), delay and response epochs is significant; the full sample epoch (-2.5 to -1.5 s) is a reading. The count holds if within 122 +/- 12.2 under one reading.
rows = []
frames = [(f"sel {s}, range {r}", pop_table(alm.TYPES, s, r)) for s, r in READINGS]
frames.append(("processed data (variant)", pop_table(alm.TYPES, None, "unit", per_unit_proc)))
for label, frame in frames:
for sample in ("sample", "full sample"):
sig = (frame[[f"p {sample}", "p delay", "p response"]] < alm.SELECTIVITY_P).any(axis=1)
rows.append({"reading": label, "sample epoch": sample, "selective": int(sig.sum()), "units": len(frame),
"fraction": sig.mean(), "count in tolerance": abs(sig.sum() - 122) <= 12.2})
print(pd.DataFrame(rows).round(3).to_string(index=False))
reading sample epoch selective units fraction count in tolerance
sel code, range unit sample 99 100 0.990 False
sel code, range unit full sample 99 100 0.990 False
sel code, range all sample 96 100 0.960 False
sel code, range all full sample 97 100 0.970 False
sel methods, range unit sample 97 105 0.924 False
sel methods, range unit full sample 100 105 0.952 False
sel methods, range all sample 95 105 0.905 False
sel methods, range all full sample 99 105 0.943 False
processed data (variant) sample 121 130 0.931 True
processed data (variant) full sample 125 130 0.962 True
Result: not reproduced. With the code's selection and the unit's range, 99 of 100 units are selective in at least one epoch, under both readings of the sample epoch. The dandiset holds only 100 included units, so the count cannot reach 122 +/- 12.2. The fraction is 0.99, against the paper's 94%. The variants give 95 to 100 selective units of 100 or 105 (0.90 to 0.97). Of the paper's 130 units, 121 (first 400 ms of the sample epoch) or 125 (full sample epoch) are selective, within tolerance. The answer fixed the selection; the outcome is unchanged. On the processed data for the 100 units the dandiset holds, 97 are selective; the count is bounded by the number of units the dandiset holds.
pt-variance-explained¶
Trial-averaged activity patterns are nearly as diverse within each PT population as across all neurons recorded in ALM; cell type explains relatively little of the variance.
Criterion (Methods): trial-averaged lick right and lick left rates of each neuron, concatenated and z-scored; PCA per PT population; variance explained by 1 to 6 components, against 1,000 size-matched samples drawn without replacement from the unidentified units of the same cohort. Reproduced when, for both populations, the median ratio of PT to unidentified variance explained is at most 1.3 at every number of components. The code's own procedure is a variant.
rows = []
for s, r in READINGS:
for ct, cohort in (("PTupper", "PTupper cohort"), ("PTlower", "PTlower cohort")):
pt = alm.selected(units, ct, s)
pool = alm.selected(units, "unidentified", s, cohort)
ve_pt = alm.variance_explained(alm.concatenated(pt, r))
x_pool = alm.concatenated(pool, r)
rng = alm.seeded(f"variance explained {ct} {s} {r}")
ve_u = np.array([alm.variance_explained(x_pool[:, rng.choice(x_pool.shape[1], len(pt), replace=False)])
for _ in range(alm.N_BOOT)])
ratio = np.median(ve_pt[None, :] / ve_u, axis=0)
rows.append({"reading": f"sel {s}, range {r}", "type": ct, "n": len(pt), "pool": len(pool),
**{f"k={k + 1}": ratio[k] for k in range(6)}, "max ratio": ratio.max(), "at most 1.3": bool(ratio.max() <= 1.3)})
ve_table = pd.DataFrame(rows)
print("Median ratio of PT to unidentified variance explained, by number of components:")
print(ve_table.round(3).to_string(index=False))
Median ratio of PT to unidentified variance explained, by number of components:
reading type n pool k=1 k=2 k=3 k=4 k=5 k=6 max ratio at most 1.3
sel code, range unit PTupper 44 226 1.055 1.081 1.062 1.072 1.072 1.060 1.081 True
sel code, range unit PTlower 56 403 1.200 1.228 1.204 1.162 1.150 1.115 1.228 True
sel code, range all PTupper 44 226 1.033 1.073 1.054 1.066 1.067 1.054 1.073 True
sel code, range all PTlower 56 403 1.214 1.235 1.211 1.167 1.151 1.115 1.235 True
sel methods, range unit PTupper 42 232 1.060 1.097 1.066 1.078 1.083 1.075 1.097 True
sel methods, range unit PTlower 63 430 1.362 1.254 1.186 1.141 1.115 1.084 1.362 False
sel methods, range all PTupper 42 232 1.043 1.096 1.064 1.076 1.078 1.066 1.096 True
sel methods, range all PTlower 63 430 1.354 1.250 1.187 1.138 1.114 1.084 1.354 False
# The code's procedure (VarianceExplained.m), main reading only: 100 ms smoothing, baseline
# subtraction and division by s.d., the first 1000 bins, trials resampled for every unit,
# variance explained per unit and averaged; 1,000 iterations.
rows = []
for ct, cohort in (("PTupper", "PTupper cohort"), ("PTlower", "PTlower cohort")):
pt = alm.selected(units, ct, "code")
pool = alm.selected(units, "unidentified", "code", cohort)
v_pt = alm.code_variance_explained(pt, "unit", alm.seeded(f"code variance explained {ct}"))
v_u = alm.code_variance_explained(pt, "unit", alm.seeded(f"code variance explained pool {ct}"), pool=pool, size=len(pt))
ratio = np.median(v_pt.mean(axis=0)[None, :] / v_u, axis=0)
z = (v_pt.mean(axis=0) - v_u.mean(axis=0)) / v_pt.std(axis=0, ddof=1)
rows.append({"type": ct, **{f"k={k + 1}": ratio[k] for k in range(6)},
"components with z > 1.645": [k + 1 for k in range(6) if z[k] > 1.645]})
print("Code's variant, median ratio by number of components:")
print(pd.DataFrame(rows).round(3).to_string(index=False))
Code's variant, median ratio by number of components: type k=1 k=2 k=3 k=4 k=5 k=6 components with z > 1.645 PTupper 1.146 1.146 1.126 1.119 1.104 1.090 [1, 2, 3, 4, 5, 6] PTlower 1.208 1.220 1.207 1.158 1.130 1.094 [1, 2, 3, 4, 5, 6]
Result: reproduced. With the code's selection and the unit's range, the median ratio of PT to unidentified variance explained is at most 1.08 for PTupper and 1.23 for PTlower over 1 to 6 components, so it stays at or below 1.3 throughout. All trials, a variant, gives at most 1.07 and 1.24. The session rule, a variant, gives a PTlower ratio for one component of 1.36 and 1.35, above the bound. On the figure's asterisks: the ratio is above 1 at every number of components in both populations, and the code's own procedure (ratios 1.09 to 1.22) puts the PT populations above the unidentified at every number of components (z > 1.645). The unidentified pools are smaller than the paper's (226 and 403 units against 495 and 511). The answer fixed the selection; the earlier run was blocked because the session rule exceeded the bound.
pt-preparatory-activity¶
Along the late delay coding direction (CDlate) selectivity is substantial in both PT types, with 40 of 61 PTupper and 44 of 69 PTlower neurons selective.
Criterion: (a) units whose late delay epoch rates (-0.4 to 0 s) differ between trial types by a two-sided Mann-Whitney test at p < 0.01, both counts within tolerance (40 +/- 4, 44 +/- 4.4); (b) selectivity along each population's CDlate is significant (5th percentile above zero) over the late delay epoch in both populations, read as (i) at every bin from -0.4 to 0 s, or (ii) for the mean over that window. Reproduced when both parts hold, partially reproduced when one does.
rows = []
frames = [(f"sel {s}, range {r}", pop_table(alm.TYPES, s, r)) for s, r in READINGS]
frames.append(("processed data (variant)", pop_table(alm.TYPES, None, "unit", per_unit_proc)))
for label, frame in frames:
sig = frame["p late delay"] < alm.SELECTIVITY_P
up, low = int(sig[frame.cell_type == "PTupper"].sum()), int(sig[frame.cell_type == "PTlower"].sum())
rows.append({"reading": label, "PTupper selective": f"{up} of {(frame.cell_type == 'PTupper').sum()}",
"PTlower selective": f"{low} of {(frame.cell_type == 'PTlower').sum()}",
"both in tolerance": abs(up - 40) <= 4 and abs(low - 44) <= 4.4})
print("(a) late delay selective units:")
print(pd.DataFrame(rows).to_string(index=False))
late_bins = alm.bins(alm.CD_WINDOWS["late"], "cd")
rows = []
for key in BOOT_KEYS:
row = {"reading": reading_label(*key)}
for ct in alm.TYPES:
s = sel_summary(ct, *key, "late")
row[f"{ct} (i) every bin"] = bool((s["q5"][late_bins] > 0).all())
row[f"{ct} (ii) mean, p"] = float((s["means"]["late delay"] <= 0).mean())
rows.append(row)
for ct in alm.TYPES:
for mk in ("subject", "animal"):
for pr in ("proj a", "proj b"):
s = sel_summary(ct, None, None, mk, pr, "late", processed=True)
rows.append({"reading": f"processed data, {ct}, mouse {mk}, {pr}", f"{ct} (i) every bin": bool((s["q5"][late_bins] > 0).all()),
f"{ct} (ii) mean, p": float((s["means"]["late delay"] <= 0).mean())})
print("\n(b) CDlate selectivity in the late delay epoch (p: fraction of iterations with the window mean at or below zero):")
print(pd.DataFrame(rows).to_string(index=False))
(a) late delay selective units:
reading PTupper selective PTlower selective both in tolerance
sel code, range unit 29 of 44 38 of 56 False
sel code, range all 28 of 44 37 of 56 False
sel methods, range unit 28 of 42 41 of 63 False
sel methods, range all 27 of 42 39 of 63 False
processed data (variant) 40 of 61 44 of 69 True
(b) CDlate selectivity in the late delay epoch (p: fraction of iterations with the window mean at or below zero):
reading PTupper (i) every bin PTupper (ii) mean, p PTlower (i) every bin PTlower (ii) mean, p
sel code, range unit, mouse subject, proj a True 0.0 True 0.0
sel code, range unit, mouse subject, proj b True 0.0 True 0.0
sel code, range unit, mouse animal, proj a True 0.0 True 0.0
sel code, range unit, mouse animal, proj b True 0.0 True 0.0
sel code, range all, mouse subject, proj a True 0.0 True 0.0
sel code, range all, mouse subject, proj b True 0.0 True 0.0
sel code, range all, mouse animal, proj a True 0.0 True 0.0
sel code, range all, mouse animal, proj b True 0.0 True 0.0
sel methods, range unit, mouse subject, proj a True 0.0 True 0.0
sel methods, range unit, mouse subject, proj b True 0.0 True 0.0
sel methods, range unit, mouse animal, proj a True 0.0 True 0.0
sel methods, range unit, mouse animal, proj b True 0.0 True 0.0
sel methods, range all, mouse subject, proj a True 0.0 True 0.0
sel methods, range all, mouse subject, proj b True 0.0 True 0.0
sel methods, range all, mouse animal, proj a True 0.0 True 0.0
sel methods, range all, mouse animal, proj b True 0.0 True 0.0
processed data, PTupper, mouse subject, proj a True 0.0 NaN NaN
processed data, PTupper, mouse subject, proj b True 0.0 NaN NaN
processed data, PTupper, mouse animal, proj a True 0.0 NaN NaN
processed data, PTupper, mouse animal, proj b True 0.0 NaN NaN
processed data, PTlower, mouse subject, proj a NaN NaN True 0.0
processed data, PTlower, mouse subject, proj b NaN NaN True 0.0
processed data, PTlower, mouse animal, proj a NaN NaN True 0.0
processed data, PTlower, mouse animal, proj b NaN NaN True 0.0
Result: partially reproduced. (a) With the code's selection and the unit's range, 29 of 44 PTupper and 38 of 56 PTlower units are selective in the late delay, below 40 +/- 4 and 44 +/- 4.4. The variants give 27 or 28 and 37 to 41. (b) CDlate selectivity is above zero at every bin of the late delay, and for its mean (0 of 1,000 iterations at or below zero), in both populations under the four open readings of the bootstrap and under the variants. On the paper's 130 units, (a) gives 40 of 61 and 44 of 69, the paper's counts, and (b) holds. The answer fixed the selection; the outcome is unchanged. The processed data for the same units give the dandiset's counts, 29 of 44 and 38 of 56; as fractions (65.9% and 67.9%) they are close to those of the paper's units (65.6% and 63.8%).
pt-upper-early-preparatory¶
Along the early coding direction (CDearly, first 400 ms of the sample epoch), selectivity is larger and more consistent across trials in the PTupper population than in the PTlower population.
Criterion: CDearly selectivity with the shared bootstrap, averaged over (a) -2.5 to -2.1 s or (b) -2.5 to -1.5 s. Larger: (i) the bootstrap median of the PTupper mean exceeds that of PTlower; (ii) in addition the PTupper 5th percentile exceeds the PTlower median. More consistent: a smaller s.d. over median of the bootstrap. Reproduced when both hold, partially reproduced when larger holds alone, not reproduced when larger does not hold. The paper does not say which rule for larger decides (claims file, from a reader's answer): both stay, and the claim is methods underspecified when they give different outcomes.
def early_rows(get, label):
out = []
for win, wname in (("CDearly window", "(a) -2.5 to -2.1 s"), ("sample epoch", "(b) -2.5 to -1.5 s")):
up, low = get("PTupper")["means"][win], get("PTlower")["means"][win]
mu, ml = np.median(up), np.median(low)
cu, cl = up.std() / mu, low.std() / ml
out.append({"reading": label, "window": wname, "PTupper median": mu, "PTupper 5th pct": np.percentile(up, 5),
"PTlower median": ml, "larger (i)": mu > ml, "larger (ii)": mu > ml and np.percentile(up, 5) > ml,
"PTupper sd/median": cu, "PTlower sd/median": cl, "more consistent": cu < cl})
return out
rows = []
for key in BOOT_KEYS:
rows += early_rows(lambda ct: sel_summary(ct, *key, "early"), reading_label(*key))
for mk in ("subject", "animal"):
for pr in ("proj a", "proj b"):
rows += early_rows(lambda ct: sel_summary(ct, None, None, mk, pr, "early", processed=True), f"processed data, mouse {mk}, {pr}")
early_table = pd.DataFrame(rows)
print(early_table.round(3).to_string(index=False))
reading window PTupper median PTupper 5th pct PTlower median larger (i) larger (ii) PTupper sd/median PTlower sd/median more consistent
sel code, range unit, mouse subject, proj a (a) -2.5 to -2.1 s 3.970 1.946 2.091 True False 0.314 0.249 False
sel code, range unit, mouse subject, proj a (b) -2.5 to -1.5 s 3.708 1.767 2.211 True False 0.343 0.214 False
sel code, range unit, mouse subject, proj b (a) -2.5 to -2.1 s 4.119 2.052 2.461 True False 0.320 0.296 False
sel code, range unit, mouse subject, proj b (b) -2.5 to -1.5 s 3.758 1.747 2.367 True False 0.351 0.283 False
sel code, range unit, mouse animal, proj a (a) -2.5 to -2.1 s 3.970 1.946 2.080 True False 0.314 0.254 False
sel code, range unit, mouse animal, proj a (b) -2.5 to -1.5 s 3.708 1.767 2.182 True False 0.343 0.221 False
sel code, range unit, mouse animal, proj b (a) -2.5 to -2.1 s 4.119 2.052 2.441 True False 0.320 0.299 False
sel code, range unit, mouse animal, proj b (b) -2.5 to -1.5 s 3.758 1.747 2.355 True False 0.351 0.276 False
sel code, range all, mouse subject, proj a (a) -2.5 to -2.1 s 3.806 1.888 1.983 True False 0.295 0.259 False
sel code, range all, mouse subject, proj a (b) -2.5 to -1.5 s 3.481 1.753 2.010 True False 0.310 0.222 False
sel code, range all, mouse subject, proj b (a) -2.5 to -2.1 s 3.938 1.952 2.365 True False 0.309 0.310 True
sel code, range all, mouse subject, proj b (b) -2.5 to -1.5 s 3.584 1.749 2.174 True False 0.323 0.301 False
sel code, range all, mouse animal, proj a (a) -2.5 to -2.1 s 3.806 1.888 1.984 True False 0.295 0.261 False
sel code, range all, mouse animal, proj a (b) -2.5 to -1.5 s 3.481 1.753 1.997 True False 0.310 0.226 False
sel code, range all, mouse animal, proj b (a) -2.5 to -2.1 s 3.938 1.952 2.327 True False 0.309 0.305 False
sel code, range all, mouse animal, proj b (b) -2.5 to -1.5 s 3.584 1.749 2.146 True False 0.323 0.285 False
sel methods, range unit, mouse subject, proj a (a) -2.5 to -2.1 s 3.928 1.873 1.858 True True 0.333 0.249 False
sel methods, range unit, mouse subject, proj a (b) -2.5 to -1.5 s 3.633 1.710 1.942 True False 0.356 0.222 False
sel methods, range unit, mouse subject, proj b (a) -2.5 to -2.1 s 4.018 1.978 2.173 True False 0.342 0.266 False
sel methods, range unit, mouse subject, proj b (b) -2.5 to -1.5 s 3.624 1.717 2.082 True False 0.369 0.274 False
sel methods, range unit, mouse animal, proj a (a) -2.5 to -2.1 s 3.928 1.873 1.868 True True 0.333 0.257 False
sel methods, range unit, mouse animal, proj a (b) -2.5 to -1.5 s 3.633 1.710 1.933 True False 0.356 0.232 False
sel methods, range unit, mouse animal, proj b (a) -2.5 to -2.1 s 4.018 1.978 2.193 True False 0.342 0.268 False
sel methods, range unit, mouse animal, proj b (b) -2.5 to -1.5 s 3.624 1.717 2.094 True False 0.369 0.270 False
sel methods, range all, mouse subject, proj a (a) -2.5 to -2.1 s 3.860 1.804 1.667 True True 0.310 0.288 False
sel methods, range all, mouse subject, proj a (b) -2.5 to -1.5 s 3.544 1.713 1.746 True False 0.332 0.246 False
sel methods, range all, mouse subject, proj b (a) -2.5 to -2.1 s 4.053 2.032 2.013 True True 0.311 0.303 False
sel methods, range all, mouse subject, proj b (b) -2.5 to -1.5 s 3.631 1.657 1.890 True False 0.343 0.304 False
sel methods, range all, mouse animal, proj a (a) -2.5 to -2.1 s 3.860 1.804 1.667 True True 0.310 0.294 False
sel methods, range all, mouse animal, proj a (b) -2.5 to -1.5 s 3.544 1.713 1.752 True False 0.332 0.252 False
sel methods, range all, mouse animal, proj b (a) -2.5 to -2.1 s 4.053 2.032 2.037 True False 0.311 0.303 False
sel methods, range all, mouse animal, proj b (b) -2.5 to -1.5 s 3.631 1.657 1.898 True False 0.343 0.305 False
processed data, mouse subject, proj a (a) -2.5 to -2.1 s 3.937 2.596 1.694 True True 0.221 0.234 True
processed data, mouse subject, proj a (b) -2.5 to -1.5 s 3.673 2.477 1.819 True True 0.234 0.192 False
processed data, mouse subject, proj b (a) -2.5 to -2.1 s 4.125 2.744 2.042 True True 0.221 0.261 True
processed data, mouse subject, proj b (b) -2.5 to -1.5 s 3.770 2.538 1.994 True True 0.236 0.243 True
processed data, mouse animal, proj a (a) -2.5 to -2.1 s 3.947 2.656 1.729 True True 0.223 0.247 True
processed data, mouse animal, proj a (b) -2.5 to -1.5 s 3.665 2.465 1.846 True True 0.239 0.211 False
processed data, mouse animal, proj b (a) -2.5 to -2.1 s 4.120 2.776 2.059 True True 0.230 0.287 True
processed data, mouse animal, proj b (b) -2.5 to -1.5 s 3.772 2.538 2.009 True True 0.241 0.262 True
Result: blocked (methods underspecified). With the code's selection and the unit's range, PTupper CDearly selectivity is larger by rule (i) in all 8 combinations of the open readings and windows (bootstrap medians 3.71 to 4.12 against 2.08 to 2.46), but by rule (ii) in none, as the PTupper 5th percentile (1.75 to 2.05) is below the PTlower median. It is more consistent in none (s.d. over median 0.31 to 0.35 for PTupper against 0.21 to 0.30). So rule (i) gives partially reproduced and rule (ii) not reproduced, in all 8. The answer kept both rules, so the claim stays blocked. The selection and trial range it fixed leave out the one combination where rule (i) gave reproduced in the earlier run (all trials, now a variant). Among the variants, rule (ii) holds in 5 of 24 combinations, all with the session rule. On the paper's units larger holds under both rules in all 8 combinations, and more consistent in 6 of 8.
pt-upper-persistent¶
In the PTupper population, selectivity along CDearly remains constant through the sample and delay epochs until the go cue.
Criterion: PTupper CDearly selectivity, shared bootstrap. Readings: (a) the 5th percentile is above zero at every bin from -2.1 s to the go cue; (b) the mean over -0.4 to 0 s is significant (fewer than 5% of iterations at or below zero). Reproduced when the readings agree that it holds; the ratio of the mean in -0.4 to 0 s to the mean in -2.5 to -2.1 s is reported.
span = (alm.TIME >= -2.1) & (alm.TIME < 0)
def persistent_row(s, label):
ratio = np.median(s["means"]["late delay"]) / np.median(s["means"]["CDearly window"])
return {"reading": label, "(a) every bin -2.1 to 0 s": bool((s["q5"][span] > 0).all()),
"bins not significant": int((s["q5"][span] <= 0).sum()), "of": int(span.sum()),
"(b) p, -0.4 to 0 s": float((s["means"]["late delay"] <= 0).mean()),
"(b) holds": bool((s["means"]["late delay"] <= 0).mean() < 0.05), "late / early": ratio}
rows = [persistent_row(sel_summary("PTupper", *key, "early"), reading_label(*key)) for key in BOOT_KEYS if key[2] == "subject"]
rows += [persistent_row(sel_summary("PTupper", None, None, "subject", pr, "early", processed=True), f"processed data, {pr}")
for pr in ("proj a", "proj b")]
print("(PTupper: the two readings of what a mouse is coincide.)")
print(pd.DataFrame(rows).round(3).to_string(index=False))
(PTupper: the two readings of what a mouse is coincide.)
reading (a) every bin -2.1 to 0 s bins not significant of (b) p, -0.4 to 0 s (b) holds late / early
sel code, range unit, mouse subject, proj a False 68 420 0.053 False 0.777
sel code, range unit, mouse subject, proj b False 160 420 0.058 False 0.739
sel code, range all, mouse subject, proj a False 118 420 0.067 False 0.739
sel code, range all, mouse subject, proj b False 165 420 0.062 False 0.702
sel methods, range unit, mouse subject, proj a False 44 420 0.049 True 0.797
sel methods, range unit, mouse subject, proj b False 107 420 0.051 False 0.775
sel methods, range all, mouse subject, proj a False 81 420 0.058 False 0.776
sel methods, range all, mouse subject, proj b False 200 420 0.062 False 0.714
processed data, proj a True 0 420 0.007 True 0.797
processed data, proj b True 0 420 0.011 True 0.758
Result: blocked (methods underspecified), as it depends on pt-upper-early-preparatory. Its own test does not hold with the code's selection and the unit's range. The PTupper 5th percentile is at or below zero in 68 (projection a) and 160 (projection b) of the 420 bins from -2.1 s to the go cue, and the mean over -0.4 to 0 s is at or below zero in 5.3% and 5.8% of iterations, so neither reading holds. The selectivity in the last 400 ms is 0.78 and 0.74 of that in the CDearly window. Among the variants, reading (b) holds in one (4.9%, session rule, unit's range, projection a). On the paper's units both readings hold: no bin is at or below zero, and p = 0.007 and 0.011. The answers left its dependency blocked, so the outcome is unchanged. On the processed data for the same 44 PTupper units neither reading holds either (85 and 153 bins not significant, p = 0.060 and 0.066), so its own test fails here because units are missing, not because of the files.
pt-lower-cdearly-decay¶
In the PTlower population, selectivity along CDearly decays rapidly, with no significant selectivity in the 400 ms before the go cue (P = 0.34, bootstrap).
Criterion: a claim of absence, judged on the paper's test. PTlower CDearly selectivity, shared bootstrap; p is the fraction of iterations with the mean over -0.4 to 0 s at or below zero. Reproduced when p >= 0.05; not reproduced when p < 0.05. Compared with 0.34; the selectivity in that window is given.
rows = []
for key in BOOT_KEYS:
s = sel_summary("PTlower", *key, "early")
m = s["means"]["late delay"]
rows.append({"reading": reading_label(*key), "p": float((m <= 0).mean()), "median": np.median(m),
"5th pct": np.percentile(m, 5), "95th pct": np.percentile(m, 95),
"early window median": np.median(s["means"]["CDearly window"])})
for mk in ("subject", "animal"):
for pr in ("proj a", "proj b"):
s = sel_summary("PTlower", None, None, mk, pr, "early", processed=True)
m = s["means"]["late delay"]
rows.append({"reading": f"processed data, mouse {mk}, {pr}", "p": float((m <= 0).mean()), "median": np.median(m),
"5th pct": np.percentile(m, 5), "95th pct": np.percentile(m, 95),
"early window median": np.median(s["means"]["CDearly window"])})
print(pd.DataFrame(rows).round(3).to_string(index=False))
reading p median 5th pct 95th pct early window median
sel code, range unit, mouse subject, proj a 0.348 0.818 -3.006 4.035 2.091
sel code, range unit, mouse subject, proj b 0.377 0.673 -3.089 4.139 2.461
sel code, range unit, mouse animal, proj a 0.379 0.736 -3.219 3.932 2.080
sel code, range unit, mouse animal, proj b 0.395 0.555 -3.310 3.940 2.441
sel code, range all, mouse subject, proj a 0.304 0.912 -2.725 3.901 1.983
sel code, range all, mouse subject, proj b 0.346 0.821 -2.938 4.059 2.365
sel code, range all, mouse animal, proj a 0.356 0.833 -2.954 3.593 1.984
sel code, range all, mouse animal, proj b 0.378 0.635 -3.221 3.562 2.327
sel methods, range unit, mouse subject, proj a 0.355 0.737 -2.963 3.606 1.858
sel methods, range unit, mouse subject, proj b 0.390 0.615 -3.138 3.489 2.173
sel methods, range unit, mouse animal, proj a 0.356 0.710 -3.139 3.257 1.868
sel methods, range unit, mouse animal, proj b 0.385 0.571 -3.267 3.257 2.193
sel methods, range all, mouse subject, proj a 0.336 0.874 -2.338 3.056 1.667
sel methods, range all, mouse subject, proj b 0.372 0.679 -2.637 2.950 2.013
sel methods, range all, mouse animal, proj a 0.319 0.878 -2.578 3.344 1.667
sel methods, range all, mouse animal, proj b 0.339 0.702 -2.749 3.252 2.037
processed data, mouse subject, proj a 0.458 0.149 -3.207 2.732 1.694
processed data, mouse subject, proj b 0.480 0.113 -3.325 2.836 2.042
processed data, mouse animal, proj a 0.455 0.230 -3.401 2.947 1.729
processed data, mouse animal, proj b 0.474 0.116 -3.450 2.812 2.059
Result: reproduced. With the code's selection and the unit's range, the PTlower CDearly selectivity in the 400 ms before the go cue is at or below zero in 35% to 40% of iterations (p = 0.348 to 0.395) under the four open readings, against the paper's P = 0.34. Its median there is 0.56 to 0.82, against 2.08 to 2.46 in the CDearly window. The variants give p = 0.30 to 0.39. On the paper's units p = 0.46 to 0.48. The answer fixed the selection; the outcome is unchanged.
pt-cd-stability¶
The coding direction stays more similar across the sample and delay epochs in the PTupper population than in the PTlower population, whose delay-epoch coding directions are largely orthogonal to those of the sample epoch.
Criterion: (a) Pearson's r across the units of a population between the CDearly and CDlate weights (Projections.m, plotEarlyLateCorr); (b) the inner product of the unit-normalised CDs of the sample and late delay epochs, mean over the bootstrap (CDCorrelation.m, compareEpochs: 100 iterations, all trials up to 500). Reproduced when PTupper's value exceeds PTlower's under both measures.
rows = []
for s, r in READINGS:
row = {"reading": f"sel {s}, range {r}"}
for ct in alm.TYPES:
pop = alm.selected(units, ct, s)
e, l = alm.coding_direction(pop, r, alm.CD_WINDOWS["early"]), alm.coding_direction(pop, r, alm.CD_WINDOWS["late"])
rr = stats.pearsonr(e, l)
row[f"{ct} r"], row[f"{ct} p"] = rr[0], rr[1]
row["PTupper higher"] = row["PTupper r"] > row["PTlower r"]
rows.append(row)
row = {"reading": "processed data (variant)"}
for ct in alm.TYPES:
pop = [u for u in proc_units if u["included"] and u["cell_type"] == ct]
rr = stats.pearsonr(alm.coding_direction(pop, "unit", alm.CD_WINDOWS["early"]), alm.coding_direction(pop, "unit", alm.CD_WINDOWS["late"]))
row[f"{ct} r"], row[f"{ct} p"] = rr[0], rr[1]
row["PTupper higher"] = row["PTupper r"] > row["PTlower r"]
rows.append(row)
print("(a) Pearson's r between CDearly and CDlate weights:")
print(pd.DataFrame(rows).round(4).to_string(index=False))
rows = []
for s, r in READINGS:
for mk in ("subject", "animal"):
row = {"reading": f"sel {s}, range {r}, mouse {mk}"}
for ct in alm.TYPES:
sim = alm.epoch_similarity(alm.selected(units, ct, s), mk, r, alm.seeded(f"epoch similarity {ct} {s} {r} {mk}"))
row[f"{ct} mean"], row[f"{ct} 1.96 s.e."] = sim.mean(), 1.96 * sim.std(ddof=1) / np.sqrt(len(sim))
row["PTupper higher"] = row["PTupper mean"] > row["PTlower mean"]
rows.append(row)
for mk in ("subject", "animal"):
row = {"reading": f"processed data, mouse {mk}"}
for ct in alm.TYPES:
pop = [u for u in proc_units if u["included"] and u["cell_type"] == ct]
sim = alm.epoch_similarity(pop, mk, "unit", alm.seeded(f"epoch similarity processed {ct} {mk}"))
row[f"{ct} mean"], row[f"{ct} 1.96 s.e."] = sim.mean(), 1.96 * sim.std(ddof=1) / np.sqrt(len(sim))
row["PTupper higher"] = row["PTupper mean"] > row["PTlower mean"]
rows.append(row)
print("\n(b) sample to late delay inner product of unit-normalised CDs (100 iterations):")
print(pd.DataFrame(rows).round(3).to_string(index=False))
(a) Pearson's r between CDearly and CDlate weights:
reading PTupper r PTupper p PTlower r PTlower p PTupper higher
sel code, range unit 0.4898 0.0007 0.1838 0.1750 True
sel code, range all 0.4863 0.0008 0.2290 0.0895 True
sel methods, range unit 0.5088 0.0006 0.1749 0.1703 True
sel methods, range all 0.5032 0.0007 0.2092 0.0998 True
processed data (variant) 0.5141 0.0000 0.1552 0.2029 True
(b) sample to late delay inner product of unit-normalised CDs (100 iterations):
reading PTupper mean PTupper 1.96 s.e. PTlower mean PTlower 1.96 s.e. PTupper higher
sel code, range unit, mouse subject 0.458 0.048 0.223 0.050 True
sel code, range unit, mouse animal 0.407 0.059 0.192 0.059 True
sel code, range all, mouse subject 0.427 0.054 0.219 0.055 True
sel code, range all, mouse animal 0.466 0.054 0.219 0.057 True
sel methods, range unit, mouse subject 0.406 0.057 0.187 0.055 True
sel methods, range unit, mouse animal 0.423 0.058 0.207 0.057 True
sel methods, range all, mouse subject 0.455 0.053 0.143 0.055 True
sel methods, range all, mouse animal 0.415 0.060 0.176 0.052 True
processed data, mouse subject 0.451 0.050 0.172 0.048 True
processed data, mouse animal 0.509 0.033 0.188 0.050 True
Result: reproduced. With the code's selection and the unit's range: (a) CDearly and CDlate weights correlate in PTupper (r = 0.49, p = 0.0007) and less in PTlower (r = 0.18, p = 0.18). (b) The sample to late delay inner product is 0.46 and 0.41 for PTupper against 0.22 and 0.19 for PTlower (the two readings of what a mouse is), with intervals (mean +/- 1.96 s.e. over 100 iterations) that do not overlap. Both measures hold under every variant (r 0.49 to 0.51 against 0.17 to 0.23; inner products 0.41 to 0.47 against 0.14 to 0.22) and on the paper's units. The PTlower values are low; the paper does not quantify "largely orthogonal". The answer fixed the selection; the outcome is unchanged.
pt-single-neuron-early-selectivity¶
Individual PTupper neurons have higher early selectivity than PTlower neurons (0.41 +/- 0.04 against 0.25 +/- 0.02, mean +/- s.e.m., P < 0.01).
Criterion: measures from the first 400 ms of the sample epoch, (a) the absolute d' (the unit's CD weight before normalisation), (b) the absolute d' at the accuracy-maximising threshold of the single-unit decoder; tests (i) two-sided Student's t-test, (ii) two-sided Mann-Whitney test, both run. Trials are those of the unit's range (claims file, from a reader's answer); all trials is a variant. Reproduced when, under one measure, both means are within the s.e.m. and PTupper is higher with p < 0.01 under both tests; partially reproduced when PTupper is higher with p < 0.01 under both tests but a mean is out of tolerance under every measure.
rows = []
frames = [(f"sel {s}, range {r}", pop_table(alm.TYPES, s, r)) for s, r in READINGS]
frames.append(("processed data (variant)", pop_table(alm.TYPES, None, "unit", per_unit_proc)))
for label, frame in frames:
for col, name in (("early d'", "(a) |d'|"), ("roc d'", "(b) |d'| at best threshold")):
x = frame.assign(v=frame[col].abs())
finite = np.isfinite(x.v)
x = x[finite]
up, low = x[x.cell_type == "PTupper"].v, x[x.cell_type == "PTlower"].v
rows.append({"reading": label, "measure": name, "infinite values dropped": int((~finite).sum()),
"PTupper mean": up.mean(), "PTupper sem": up.sem(), "PTlower mean": low.mean(), "PTlower sem": low.sem(),
"t p": stats.ttest_ind(up, low).pvalue, "MW p": stats.mannwhitneyu(up, low, alternative="two-sided").pvalue,
"means in tolerance": abs(up.mean() - 0.41) <= 0.04 and abs(low.mean() - 0.25) <= 0.02})
print(pd.DataFrame(rows).round(4).to_string(index=False))
reading measure infinite values dropped PTupper mean PTupper sem PTlower mean PTlower sem t p MW p means in tolerance
sel code, range unit (a) |d'| 0 0.4100 0.0507 0.2365 0.0236 0.0012 0.0041 True
sel code, range unit (b) |d'| at best threshold 10 0.8190 0.0808 0.6197 0.0512 0.0343 0.0569 False
sel code, range all (a) |d'| 0 0.3621 0.0470 0.2231 0.0220 0.0051 0.0260 False
sel code, range all (b) |d'| at best threshold 11 0.7612 0.0748 0.5962 0.0477 0.0589 0.0892 False
sel methods, range unit (a) |d'| 0 0.4162 0.0529 0.2288 0.0223 0.0004 0.0019 False
sel methods, range unit (b) |d'| at best threshold 15 0.8240 0.0827 0.6074 0.0483 0.0202 0.0580 False
sel methods, range all (a) |d'| 0 0.3703 0.0488 0.2135 0.0209 0.0013 0.0085 False
sel methods, range all (b) |d'| at best threshold 15 0.7652 0.0765 0.5787 0.0451 0.0309 0.0518 False
processed data (variant) (a) |d'| 0 0.4136 0.0415 0.2475 0.0218 0.0004 0.0019 True
processed data (variant) (b) |d'| at best threshold 14 0.8358 0.0663 0.6436 0.0484 0.0198 0.0265 False
Result: reproduced. With the code's selection and the unit's range, measure (a) gives 0.410 +/- 0.051 for PTupper and 0.237 +/- 0.024 for PTlower (mean +/- s.e.m.), both within tolerance (0.41 +/- 0.04, 0.25 +/- 0.02), and PTupper is higher with p = 0.0012 (Student's t-test) and 0.0041 (Mann-Whitney), both below 0.01. Measure (b) gives 0.82 against 0.62, with 10 infinite values dropped, t-test p = 0.034 and Mann-Whitney p = 0.057; a match under one measure counts. Among the variants, all trials gives 0.362 and 0.223 with Mann-Whitney p = 0.026, and the session rule gives 0.416 and 0.229 (PTlower out of tolerance) with p = 0.0004 and 0.0019. On the paper's units measure (a) gives 0.414 and 0.248 (p = 0.0004 and 0.0019), the paper's numbers. The answers fixed the selection and the trial range; the earlier run was blocked because these readings gave three outcomes.
pt-upper-decoding¶
Single PTupper neurons decode trial type better than PTlower neurons in the first 400 ms of the sample epoch; 24.6% (15 of 61) of PTupper and 4.4% (3 of 69) of PTlower neurons reach 70% accuracy, mean accuracy is 64.4 +/- 1.0% against 58.9 +/- 0.6%, and the ten most discriminative neurons are all PTupper.
Criterion: per unit, the best single-threshold accuracy from the mean rate of each correct trial in -2.5 to -2.1 s. "At least 70%" is read as (a) accuracy >= 0.70 or (b) the code's cut, above the mean + 1 s.d. of all PT accuracies; the counts hold if they match under one. Reproduced when all four numbers are within tolerance (15 +/- 1.5, 3 +/- 0.3, 64.4 +/- 1.0%, 58.9 +/- 0.6%), the two-sided Mann-Whitney test has p < 0.05 with PTupper higher, and the ten most accurate units are all PTupper; partially reproduced when the test holds and some of the rest do not.
def decoding_rows(frame, label, col, n_top):
up, low = frame[frame.cell_type == "PTupper"][col], frame[frame.cell_type == "PTlower"][col]
allv = frame[col]
cut = allv.mean() + allv.std()
top = frame.sort_values(col, ascending=False).head(n_top)
return {"reading": label, "PTupper >= 0.70": int((up >= 0.70).sum()), "PTlower >= 0.70": int((low >= 0.70).sum()),
"cut (b)": cut, "PTupper > cut": int((up > cut).sum()), "PTlower > cut": int((low > cut).sum()),
"n PTupper": len(up), "n PTlower": len(low),
"PTupper mean %": 100 * up.mean(), "PTupper sem %": 100 * up.sem(), "PTlower mean %": 100 * low.mean(),
"PTlower sem %": 100 * low.sem(), "MW p": stats.mannwhitneyu(up, low, alternative="two-sided").pvalue,
f"PTupper in top {n_top}": int((top.cell_type == "PTupper").sum())}
rows = [decoding_rows(pop_table(alm.TYPES, s, r), f"sel {s}, range {r}", "accuracy", 10) for s, r in READINGS]
rows.append(decoding_rows(pop_table(alm.TYPES, None, "unit", per_unit_proc), "processed data (variant)", "accuracy", 10))
print(pd.DataFrame(rows).round(4).T.to_string())
0 1 2 3 4 reading sel code, range unit sel code, range all sel methods, range unit sel methods, range all processed data (variant) PTupper >= 0.70 9 7 9 7 15 PTlower >= 0.70 3 1 3 1 3 cut (b) 0.6833 0.674 0.681 0.6718 0.686 PTupper > cut 9 9 9 9 15 PTlower > cut 3 3 3 3 3 n PTupper 44 44 42 42 61 n PTlower 56 56 63 63 69 PTupper mean % 64.282 63.5136 64.4295 63.7107 64.3781 PTupper sem % 1.2546 1.189 1.3076 1.2333 1.036 PTlower mean % 58.7142 58.3463 58.3372 57.9912 58.8171 PTlower sem % 0.6591 0.6306 0.6611 0.6228 0.6028 MW p 0.0001 0.0001 0.0 0.0 0.0 PTupper in top 10 8 8 8 8 10
Result: partially reproduced. With the code's selection and the unit's range, mean accuracy is higher in PTupper units (Mann-Whitney p = 0.0001 to four decimals): 64.3 +/- 1.3% against 58.7 +/- 0.7%, both within tolerance. But only 9 of 44 PTupper units reach 70%, and 9 pass the code's cut (0.683), against 15 +/- 1.5; 3 PTlower units reach 70% and 3 pass the cut, the paper's 3; and 8 of the 10 most accurate units are PTupper, not all 10. The variants give 7 to 9 PTupper units, means of 63.5% to 64.4% and 58.0% to 58.3%, and 8 of 10. On the paper's units: 15 of 61 and 3 of 69, means of 64.4% and 58.8%, and all ten best units PTupper, the paper's numbers. The answer fixed the selection; the outcome is unchanged. The processed data for the same units give 9 PTupper and 3 PTlower units at 70% and 8 PTupper units in the top 10, as the dandiset, so the shortfall comes from the missing units.
pt-upper-decoding-rectified¶
Decoding from spike rates rectified at baseline, no PTlower neuron reaches 70% accuracy while the same 24.6% of PTupper neurons do, PTupper neurons make up 20 of the 21 most predictive, and mean accuracy is 62.7 +/- 1.1% against 56.7 +/- 0.4%.
Criterion: as pt-upper-decoding, after subtracting the unit's baseline (the mean over the two trial types of their mean rate in -3.1 to -2.7 s) and setting negative values to zero; under reading (b) the cut is recomputed from the rectified accuracies. Reproduced when all numbers are within tolerance (0, 15 +/- 1.5, 20 +/- 2 of 21, 62.7 +/- 1.1%, 56.7 +/- 0.4%) and the Mann-Whitney test has p < 0.05 with PTupper higher; partially reproduced when the test holds and a number is out of tolerance.
rows = [decoding_rows(pop_table(alm.TYPES, s, r), f"sel {s}, range {r}", "rectified accuracy", 21) for s, r in READINGS]
rows.append(decoding_rows(pop_table(alm.TYPES, None, "unit", per_unit_proc), "processed data (variant)", "rectified accuracy", 21))
print(pd.DataFrame(rows).round(4).T.to_string())
0 1 2 3 4 reading sel code, range unit sel code, range all sel methods, range unit sel methods, range all processed data (variant) PTupper >= 0.70 8 6 8 6 15 PTlower >= 0.70 0 0 0 0 0 cut (b) 0.6544 0.6483 0.6509 0.6449 0.6639 PTupper > cut 11 10 11 11 16 PTlower > cut 1 1 1 1 1 n PTupper 44 44 42 42 61 n PTlower 56 56 63 63 69 PTupper mean % 61.8081 61.3283 61.9631 61.4958 62.7314 PTupper sem % 1.2698 1.1948 1.3261 1.2453 1.0574 PTlower mean % 56.1981 56.1348 55.9024 55.8492 56.6024 PTlower sem % 0.4808 0.4781 0.4542 0.4424 0.4278 MW p 0.0 0.0001 0.0 0.0 0.0 PTupper in top 21 16 16 16 17 20
Result: partially reproduced. With the code's selection and the unit's range, no PTlower unit reaches 70% with rectified rates, as in the paper, and the Mann-Whitney test holds (p = 0.0000 to four decimals). But only 8 of 44 PTupper units reach 70% (11 pass the recomputed cut of 0.654, as does 1 PTlower unit), against 15 +/- 1.5; 16 of the 21 most predictive units are PTupper, against 20 +/- 2; and while the PTupper mean, 61.8 +/- 1.3%, is within 62.7 +/- 1.1%, the PTlower mean, 56.2 +/- 0.5%, is below 56.7 +/- 0.4%. The variants give 6 to 8 PTupper units at 70%, 16 or 17 of 21, and PTlower means of 55.8% to 56.1%. On the paper's units: 15 and 0 units at 70%, 20 of 21, and means of 62.7% and 56.6%, the paper's numbers. The answer fixed the selection; the outcome is unchanged. The processed data for the same units give 9 PTupper and 0 PTlower units at 70%, 15 PTupper units in the top 21, and means of 62.33% and 56.84%: the counts fall short because units are missing, and the PTlower mean is out of tolerance on the files and within it on the same units processed. The files' baseline rates, the rectification threshold, are higher, which may explain the difference, but the notebook does not isolate the baseline as its cause.
pt-cdgo-onset¶
Along CDgo, significant selectivity emerges 24 ms after the go cue in PTlower neurons, faster than in PTupper neurons (48 ms), consistent with a role in movement initiation.
Criterion: onset is the first time after the go cue at which the 5th percentile of the bootstrap selectivity along CDgo is above zero; tolerance one 5 ms bin. Readings: the bin's time stamp (a) centre, (b) end, (c) start; (i) the first significant bin, (ii) the first of three consecutive significant bins. Reproduced when, under one of these readings, both onsets are within tolerance and PTlower is earlier; partially reproduced when PTlower is earlier but an onset is out of tolerance.
STAMPS = [(st, k) for st in ("centre", "end", "start") for k in (1, 3)]
def onset_rows(get, label):
out = []
for st, k in STAMPS:
lo, up = alm.onset(get("PTlower")["q5"], st, k), alm.onset(get("PTupper")["q5"], st, k)
out.append({"reading": label, "stamp": st, "consecutive": k, "PTlower ms": 1000 * lo, "PTupper ms": 1000 * up,
"PTlower earlier": lo < up, "in tolerance": abs(1000 * lo - 24) <= 5 and abs(1000 * up - 48) <= 5})
return out
rows = []
for key in BOOT_KEYS:
rows += onset_rows(lambda ct: sel_summary(ct, *key, "go"), reading_label(*key))
onset_table = pd.DataFrame(rows)
print(onset_table.round(1).to_string(index=False))
rows = []
for mk in ("subject", "animal"):
for pr in ("proj a", "proj b"):
rows += onset_rows(lambda ct: sel_summary(ct, None, None, mk, pr, "go", processed=True), f"processed data, mouse {mk}, {pr}")
print("\nVariant, the paper's units from the processed data:")
print(pd.DataFrame(rows).round(1).to_string(index=False))
reading stamp consecutive PTlower ms PTupper ms PTlower earlier in tolerance
sel code, range unit, mouse subject, proj a centre 1 37.5 77.5 True False
sel code, range unit, mouse subject, proj a centre 3 37.5 77.5 True False
sel code, range unit, mouse subject, proj a end 1 40.0 80.0 True False
sel code, range unit, mouse subject, proj a end 3 40.0 80.0 True False
sel code, range unit, mouse subject, proj a start 1 35.0 75.0 True False
sel code, range unit, mouse subject, proj a start 3 35.0 75.0 True False
sel code, range unit, mouse subject, proj b centre 1 37.5 82.5 True False
sel code, range unit, mouse subject, proj b centre 3 37.5 82.5 True False
sel code, range unit, mouse subject, proj b end 1 40.0 85.0 True False
sel code, range unit, mouse subject, proj b end 3 40.0 85.0 True False
sel code, range unit, mouse subject, proj b start 1 35.0 80.0 True False
sel code, range unit, mouse subject, proj b start 3 35.0 80.0 True False
sel code, range unit, mouse animal, proj a centre 1 32.5 77.5 True False
sel code, range unit, mouse animal, proj a centre 3 32.5 77.5 True False
sel code, range unit, mouse animal, proj a end 1 35.0 80.0 True False
sel code, range unit, mouse animal, proj a end 3 35.0 80.0 True False
sel code, range unit, mouse animal, proj a start 1 30.0 75.0 True False
sel code, range unit, mouse animal, proj a start 3 30.0 75.0 True False
sel code, range unit, mouse animal, proj b centre 1 32.5 82.5 True False
sel code, range unit, mouse animal, proj b centre 3 32.5 82.5 True False
sel code, range unit, mouse animal, proj b end 1 35.0 85.0 True False
sel code, range unit, mouse animal, proj b end 3 35.0 85.0 True False
sel code, range unit, mouse animal, proj b start 1 30.0 80.0 True False
sel code, range unit, mouse animal, proj b start 3 30.0 80.0 True False
sel code, range all, mouse subject, proj a centre 1 37.5 82.5 True False
sel code, range all, mouse subject, proj a centre 3 37.5 82.5 True False
sel code, range all, mouse subject, proj a end 1 40.0 85.0 True False
sel code, range all, mouse subject, proj a end 3 40.0 85.0 True False
sel code, range all, mouse subject, proj a start 1 35.0 80.0 True False
sel code, range all, mouse subject, proj a start 3 35.0 80.0 True False
sel code, range all, mouse subject, proj b centre 1 37.5 82.5 True False
sel code, range all, mouse subject, proj b centre 3 37.5 82.5 True False
sel code, range all, mouse subject, proj b end 1 40.0 85.0 True False
sel code, range all, mouse subject, proj b end 3 40.0 85.0 True False
sel code, range all, mouse subject, proj b start 1 35.0 80.0 True False
sel code, range all, mouse subject, proj b start 3 35.0 80.0 True False
sel code, range all, mouse animal, proj a centre 1 37.5 82.5 True False
sel code, range all, mouse animal, proj a centre 3 37.5 82.5 True False
sel code, range all, mouse animal, proj a end 1 40.0 85.0 True False
sel code, range all, mouse animal, proj a end 3 40.0 85.0 True False
sel code, range all, mouse animal, proj a start 1 35.0 80.0 True False
sel code, range all, mouse animal, proj a start 3 35.0 80.0 True False
sel code, range all, mouse animal, proj b centre 1 37.5 82.5 True False
sel code, range all, mouse animal, proj b centre 3 37.5 82.5 True False
sel code, range all, mouse animal, proj b end 1 40.0 85.0 True False
sel code, range all, mouse animal, proj b end 3 40.0 85.0 True False
sel code, range all, mouse animal, proj b start 1 35.0 80.0 True False
sel code, range all, mouse animal, proj b start 3 35.0 80.0 True False
sel methods, range unit, mouse subject, proj a centre 1 37.5 77.5 True False
sel methods, range unit, mouse subject, proj a centre 3 37.5 77.5 True False
sel methods, range unit, mouse subject, proj a end 1 40.0 80.0 True False
sel methods, range unit, mouse subject, proj a end 3 40.0 80.0 True False
sel methods, range unit, mouse subject, proj a start 1 35.0 75.0 True False
sel methods, range unit, mouse subject, proj a start 3 35.0 75.0 True False
sel methods, range unit, mouse subject, proj b centre 1 42.5 77.5 True False
sel methods, range unit, mouse subject, proj b centre 3 42.5 77.5 True False
sel methods, range unit, mouse subject, proj b end 1 45.0 80.0 True False
sel methods, range unit, mouse subject, proj b end 3 45.0 80.0 True False
sel methods, range unit, mouse subject, proj b start 1 40.0 75.0 True False
sel methods, range unit, mouse subject, proj b start 3 40.0 75.0 True False
sel methods, range unit, mouse animal, proj a centre 1 37.5 77.5 True False
sel methods, range unit, mouse animal, proj a centre 3 37.5 77.5 True False
sel methods, range unit, mouse animal, proj a end 1 40.0 80.0 True False
sel methods, range unit, mouse animal, proj a end 3 40.0 80.0 True False
sel methods, range unit, mouse animal, proj a start 1 35.0 75.0 True False
sel methods, range unit, mouse animal, proj a start 3 35.0 75.0 True False
sel methods, range unit, mouse animal, proj b centre 1 37.5 77.5 True False
sel methods, range unit, mouse animal, proj b centre 3 37.5 77.5 True False
sel methods, range unit, mouse animal, proj b end 1 40.0 80.0 True False
sel methods, range unit, mouse animal, proj b end 3 40.0 80.0 True False
sel methods, range unit, mouse animal, proj b start 1 35.0 75.0 True False
sel methods, range unit, mouse animal, proj b start 3 35.0 75.0 True False
sel methods, range all, mouse subject, proj a centre 1 37.5 77.5 True False
sel methods, range all, mouse subject, proj a centre 3 37.5 77.5 True False
sel methods, range all, mouse subject, proj a end 1 40.0 80.0 True False
sel methods, range all, mouse subject, proj a end 3 40.0 80.0 True False
sel methods, range all, mouse subject, proj a start 1 35.0 75.0 True False
sel methods, range all, mouse subject, proj a start 3 35.0 75.0 True False
sel methods, range all, mouse subject, proj b centre 1 37.5 77.5 True False
sel methods, range all, mouse subject, proj b centre 3 37.5 77.5 True False
sel methods, range all, mouse subject, proj b end 1 40.0 80.0 True False
sel methods, range all, mouse subject, proj b end 3 40.0 80.0 True False
sel methods, range all, mouse subject, proj b start 1 35.0 75.0 True False
sel methods, range all, mouse subject, proj b start 3 35.0 75.0 True False
sel methods, range all, mouse animal, proj a centre 1 37.5 77.5 True False
sel methods, range all, mouse animal, proj a centre 3 37.5 77.5 True False
sel methods, range all, mouse animal, proj a end 1 40.0 80.0 True False
sel methods, range all, mouse animal, proj a end 3 40.0 80.0 True False
sel methods, range all, mouse animal, proj a start 1 35.0 75.0 True False
sel methods, range all, mouse animal, proj a start 3 35.0 75.0 True False
sel methods, range all, mouse animal, proj b centre 1 37.5 77.5 True False
sel methods, range all, mouse animal, proj b centre 3 37.5 77.5 True False
sel methods, range all, mouse animal, proj b end 1 40.0 80.0 True False
sel methods, range all, mouse animal, proj b end 3 40.0 80.0 True False
sel methods, range all, mouse animal, proj b start 1 35.0 75.0 True False
sel methods, range all, mouse animal, proj b start 3 35.0 75.0 True False
Variant, the paper's units from the processed data:
reading stamp consecutive PTlower ms PTupper ms PTlower earlier in tolerance
processed data, mouse subject, proj a centre 1 37.5 57.5 True False
processed data, mouse subject, proj a centre 3 37.5 57.5 True False
processed data, mouse subject, proj a end 1 40.0 60.0 True False
processed data, mouse subject, proj a end 3 40.0 60.0 True False
processed data, mouse subject, proj a start 1 35.0 55.0 True False
processed data, mouse subject, proj a start 3 35.0 55.0 True False
processed data, mouse subject, proj b centre 1 37.5 62.5 True False
processed data, mouse subject, proj b centre 3 37.5 62.5 True False
processed data, mouse subject, proj b end 1 40.0 65.0 True False
processed data, mouse subject, proj b end 3 40.0 65.0 True False
processed data, mouse subject, proj b start 1 35.0 60.0 True False
processed data, mouse subject, proj b start 3 35.0 60.0 True False
processed data, mouse animal, proj a centre 1 32.5 62.5 True False
processed data, mouse animal, proj a centre 3 32.5 62.5 True False
processed data, mouse animal, proj a end 1 35.0 65.0 True False
processed data, mouse animal, proj a end 3 35.0 65.0 True False
processed data, mouse animal, proj a start 1 30.0 60.0 True False
processed data, mouse animal, proj a start 3 30.0 60.0 True False
processed data, mouse animal, proj b centre 1 32.5 67.5 True False
processed data, mouse animal, proj b centre 3 32.5 77.5 True False
processed data, mouse animal, proj b end 1 35.0 70.0 True False
processed data, mouse animal, proj b end 3 35.0 80.0 True False
processed data, mouse animal, proj b start 1 30.0 65.0 True False
processed data, mouse animal, proj b start 3 30.0 75.0 True False
# The 5th percentile of CDgo selectivity in the first 100 ms after the go cue (bin centres), selection
# (a), range (a), mouse (a), projection (a): on the dandiset, on the processed data for the units the
# dandiset holds, and on all the paper's units.
first = (alm.TIME > 0) & (alm.TIME <= 0.1)
tab = {"time ms": 1000 * alm.TIME[first]}
for ct in alm.TYPES:
tab[f"{ct} files"] = sel_summary(ct, "code", "unit", "subject", "proj a", "go")["q5"][first]
tab[f"{ct} processed, same units"] = boots_same[ct][("go", "proj a", "go")]["q5"][first]
tab[f"{ct} processed, all"] = sel_summary(ct, None, None, "subject", "proj a", "go", processed=True)["q5"][first]
print(pd.DataFrame(tab).round(2).to_string(index=False))
time ms PTupper files PTupper processed, same units PTupper processed, all PTlower files PTlower processed, same units PTlower processed, all
2.5 -0.48 -0.58 -0.52 -1.97 -1.42 -1.26
7.5 -0.50 -0.63 -0.52 -1.98 -1.40 -1.25
12.5 -0.47 -0.64 -0.53 -2.20 -1.50 -1.29
17.5 -0.35 -0.70 -0.46 -1.92 -1.49 -1.38
22.5 -0.32 -0.46 -0.31 -1.56 -1.23 -1.10
27.5 -0.32 -0.58 -0.40 -0.91 -0.89 -0.64
32.5 -0.48 -0.62 -0.43 -0.29 -0.09 -0.08
37.5 -0.45 -0.70 -0.55 0.29 0.29 0.29
42.5 -0.31 -0.61 -0.41 0.70 0.70 0.75
47.5 -0.33 -0.45 -0.16 1.06 1.07 1.11
52.5 -0.32 -0.41 -0.12 1.40 1.54 1.50
57.5 -0.17 -0.36 0.03 2.04 2.17 1.96
62.5 -0.19 -0.38 0.07 2.71 2.59 2.35
67.5 -0.15 -0.26 0.24 3.01 3.01 2.63
72.5 -0.21 -0.30 0.18 3.16 3.25 2.86
77.5 0.08 -0.04 0.52 3.28 3.55 3.32
82.5 0.19 0.07 0.71 3.61 3.82 3.47
87.5 0.38 0.17 0.92 4.16 4.11 3.82
92.5 0.42 0.25 1.00 4.56 4.46 4.30
97.5 0.52 0.32 1.03 4.76 4.76 4.50
Result: partially reproduced. With the code's selection and the unit's range, PTlower selectivity along CDgo becomes significant before PTupper's under all 24 combinations of the open readings: 30 to 40 ms after the go cue, against 75 to 85 ms. Under none is either onset within one bin of the paper's 24 and 48 ms. The variants give the same order. On the paper's units the onsets are 30 to 40 ms for PTlower and 55 to 80 ms for PTupper, also outside tolerance. The answer fixed the selection; the outcome is unchanged. In the first 100 ms after the go cue (second cell), the PTlower 5th percentile first rises above zero in the 37.5 ms bin in the files, in the processed data for the same units and in all the paper's units alike, so its distance from 24 ms does not come from the data but from the method shared by all three; what in it is not tested here. The PTupper 5th percentile first rises above zero at 77.5 ms in the files, 82.5 ms in the same units processed and 57.5 ms in all the paper's units, so most of its delay comes from the missing units.
pt-lower-late-and-motor¶
PTlower neurons carry late preparatory activity (selectivity along CDlate) and a movement command; along the go cue coding direction (CDgo) selectivity is larger, emerges faster and persists longer in PTlower than in PTupper.
Criterion: late preparatory activity is part (b) of pt-preparatory-activity for PTlower. Along CDgo: larger, (a) the bootstrap median of the mean over 0 to 0.4 s, (b) the peak of the median selectivity after the go cue, PTlower above PTupper; faster, from pt-cdgo-onset; persisted longer, the last time of the first run of significant bins after the go cue is later for PTlower. Reproduced when PTlower CDlate selectivity holds and larger, faster and longer all hold; partially reproduced when some do.
after = alm.TIME > 0
rows = []
for key in BOOT_KEYS:
lo, up = sel_summary("PTlower", *key, "go"), sel_summary("PTupper", *key, "go")
late = sel_summary("PTlower", *key, "late")
row = {"reading": reading_label(*key),
"CDlate (i)": bool((late["q5"][late_bins] > 0).all()), "CDlate (ii)": bool((late["means"]["late delay"] <= 0).mean() < 0.05),
"mean 0-0.4 s PTlower": np.median(lo["means"]["response"]), "PTupper": np.median(up["means"]["response"]),
"peak PTlower": lo["q50"][after].max(), "peak PTupper": up["q50"][after].max()}
row["larger (a)"] = row["mean 0-0.4 s PTlower"] > row["PTupper"]
row["larger (b)"] = row["peak PTlower"] > row["peak PTupper"]
faster = [alm.onset(lo["q5"], st, k) < alm.onset(up["q5"], st, k) for st, k in STAMPS]
row["faster (all stamps)"] = all(faster)
row["run end PTlower s"] = alm.first_run_end(lo["q5"])
row["run end PTupper s"] = alm.first_run_end(up["q5"])
row["longer"] = row["run end PTlower s"] > row["run end PTupper s"]
rows.append(row)
motor_table = pd.DataFrame(rows)
print(motor_table.round(3).to_string(index=False))
print("\nRun ends under the other time stamps differ by half a bin and give the same order:",
all((alm.first_run_end(sel_summary("PTlower", *k, "go")["q5"], st) > alm.first_run_end(sel_summary("PTupper", *k, "go")["q5"], st))
== motor_table.longer.iloc[i] for i, k in enumerate(BOOT_KEYS) for st in ("end", "start")))
rows = []
for mk in ("subject", "animal"):
for pr in ("proj a", "proj b"):
lo, up = (sel_summary(ct, None, None, mk, pr, "go", processed=True) for ct in ("PTlower", "PTupper"))
rows.append({"reading": f"processed data, mouse {mk}, {pr}", "mean 0-0.4 s PTlower": np.median(lo["means"]["response"]),
"PTupper": np.median(up["means"]["response"]), "peak PTlower": lo["q50"][after].max(),
"peak PTupper": up["q50"][after].max(), "run end PTlower s": alm.first_run_end(lo["q5"]),
"run end PTupper s": alm.first_run_end(up["q5"])})
print("\nVariant, the paper's units:")
print(pd.DataFrame(rows).round(3).to_string(index=False))
reading CDlate (i) CDlate (ii) mean 0-0.4 s PTlower PTupper peak PTlower peak PTupper larger (a) larger (b) faster (all stamps) run end PTlower s run end PTupper s longer
sel code, range unit, mouse subject, proj a True True 13.251 3.938 16.946 6.183 True True True 1.442 1.308 True
sel code, range unit, mouse subject, proj b True True 13.405 4.038 17.006 6.242 True True True 1.438 1.232 True
sel code, range unit, mouse animal, proj a True True 13.730 3.938 17.295 6.183 True True True 1.608 1.308 True
sel code, range unit, mouse animal, proj b True True 13.810 4.038 17.392 6.242 True True True 1.568 1.232 True
sel code, range all, mouse subject, proj a True True 11.796 3.520 15.113 5.658 True True True 1.418 1.308 True
sel code, range all, mouse subject, proj b True True 12.055 3.641 15.259 5.836 True True True 1.412 1.288 True
sel code, range all, mouse animal, proj a True True 12.111 3.520 15.446 5.658 True True True 1.552 1.308 True
sel code, range all, mouse animal, proj b True True 12.219 3.641 15.630 5.836 True True True 1.438 1.288 True
sel methods, range unit, mouse subject, proj a True True 13.980 4.071 17.599 6.292 True True True 1.442 1.308 True
sel methods, range unit, mouse subject, proj b True True 13.983 4.190 17.797 6.350 True True True 1.428 1.298 True
sel methods, range unit, mouse animal, proj a True True 14.166 4.071 17.720 6.292 True True True 1.552 1.308 True
sel methods, range unit, mouse animal, proj b True True 14.352 4.190 17.780 6.350 True True True 1.478 1.298 True
sel methods, range all, mouse subject, proj a True True 12.555 3.593 15.954 5.611 True True True 1.422 1.302 True
sel methods, range all, mouse subject, proj b True True 12.765 3.682 15.967 5.800 True True True 1.412 1.288 True
sel methods, range all, mouse animal, proj a True True 12.264 3.593 15.213 5.611 True True True 1.428 1.302 True
sel methods, range all, mouse animal, proj b True True 12.346 3.682 15.401 5.800 True True True 1.422 1.288 True
Run ends under the other time stamps differ by half a bin and give the same order: True
Variant, the paper's units:
reading mean 0-0.4 s PTlower PTupper peak PTlower peak PTupper run end PTlower s run end PTupper s
processed data, mouse subject, proj a 13.134 5.231 16.355 7.796 1.458 0.728
processed data, mouse subject, proj b 13.151 5.320 16.429 7.844 1.442 0.602
processed data, mouse animal, proj a 13.607 5.263 16.930 7.845 1.622 0.602
processed data, mouse animal, proj b 13.649 5.303 16.910 7.921 1.608 0.068
Result: reproduced. With the code's selection and the unit's range, under the four open readings, PTlower CDlate selectivity holds (every bin of the late delay and the window mean are significant). Along CDgo, PTlower selectivity is larger (median over 0 to 0.4 s 13.3 to 13.8, against 3.9 to 4.0; peak 16.9 to 17.4, against 6.2), emerges faster (pt-cdgo-onset, under every time stamp), and persists longer: its first significant run ends 1.44 to 1.61 s after the go cue, against 1.23 to 1.31 s. On the paper's units the same holds, with the PTupper run ending 0.07 to 0.73 s after the go cue. The answer fixed the selection; the outcome is unchanged.
reaction-time (blocked)¶
The first detectable movement occurs about 50 ms after the go cue (99% confidence interval 38 to 56 ms).
The test needs the 400 Hz side-view video. The data check shows that no file holds video; the acquisition group of every file holds the lick times only.
print("Acquisition objects in the files:", sorted({k for s in sessions for k in s["meta"]["acquisition"]}),
"; any video or image object:", any(s["meta"]["has_video"] for s in sessions))
Acquisition objects in the files: ['lick_times'] ; any video or image object: False
Result: blocked (data missing). The files hold lick times only. There is no video from which to measure the reaction time, and the paper deposits none elsewhere.
pt-cd-change-at-go-cue¶
The coding direction changes rapidly at the go cue in the PTlower population, and more slowly, over several hundred milliseconds, in the PTupper population.
Criterion: the inner product of the unit-normalised CD at each bin with CDlate, normalised to its mean in the 400 ms before the go cue, 1,000 bootstrap iterations; the time after the go cue at which the bootstrap median first falls below (a) 0.5 or (b) 0.75. Reproduced when PTlower's time is shorter than PTupper's and PTupper's is at least 0.2 s; partially reproduced when PTlower's is shorter only. The code's version (the mean over the delay epoch's bins, 100 iterations) is a variant. Neither the paper nor its code fixes the level (claims file, from a reader's answer): both stay, and the claim is methods underspecified when they give different outcomes.
t0 = time.time()
window = (-0.5, 2.0)
tsub = alm.TIME[(alm.TIME >= window[0]) & (alm.TIME <= window[1])]
pre = (tsub > -0.4) & (tsub <= 0)
def fall_times(corr):
norm = corr / corr[:, pre].mean(axis=1, keepdims=True)
med = np.median(norm, axis=0)
out = {}
for level in (0.5, 0.75):
below = np.flatnonzero((tsub > 0) & (med < level))
out[level] = float(tsub[below[0]]) if len(below) else np.nan
return out, med
rows = []
curves = {}
for s, r in READINGS:
for mk in ("subject", "animal"):
row = {"reading": f"sel {s}, range {r}, mouse {mk}"}
for ct in alm.TYPES:
if ct == "PTupper" and mk == "animal":
ft = rows[-1][f"{ct} fall"]
else:
corr = alm.binned_cd_correlation(alm.selected(units, ct, s), mk, r, alm.seeded(f"cd change {ct} {s} {r} {mk}"), window)
ft, med = fall_times(corr)
curves[(ct, s, r, mk)] = med
row[f"{ct} fall"] = ft
for level in (0.5, 0.75):
lo, up = row["PTlower fall"][level], row["PTupper fall"][level]
row[f"level {level}: PTlower s"], row[f"level {level}: PTupper s"] = lo, up
row[f"level {level}: holds"] = bool(lo < up and up >= 0.2)
row[f"level {level}: PTlower shorter"] = bool(lo < up)
rows.append(row)
change_table = pd.DataFrame(rows).drop(columns=["PTupper fall", "PTlower fall"])
print(change_table.round(4).to_string(index=False))
rows = []
for ct in alm.TYPES:
corr = alm.binned_cd_correlation(alm.selected(units, ct, "code"), "subject", "unit", alm.seeded(f"cd change code {ct}"),
window, n_iter=100, n_trials=500, reference="delay")
ft, _ = fall_times(corr)
rows.append({"variant": "code (delay epoch reference, 100 iterations)", "type": ct, "below 0.5 s": ft[0.5], "below 0.75 s": ft[0.75]})
for ct in alm.TYPES:
pop = [u for u in proc_units if u["included"] and u["cell_type"] == ct]
ft, _ = fall_times(alm.binned_cd_correlation(pop, "subject", "unit", alm.seeded(f"cd change processed {ct}"), window))
rows.append({"variant": "processed data", "type": ct, "below 0.5 s": ft[0.5], "below 0.75 s": ft[0.75]})
print("\nVariants:")
print(pd.DataFrame(rows).round(4).to_string(index=False))
print(f"({time.time() - t0:.0f} s)")
reading level 0.5: PTlower s level 0.5: PTupper s level 0.5: holds level 0.5: PTlower shorter level 0.75: PTlower s level 0.75: PTupper s level 0.75: holds level 0.75: PTlower shorter
sel code, range unit, mouse subject 0.3675 0.3725 True True 0.0625 0.1925 False True
sel code, range unit, mouse animal 0.2725 0.3725 True True 0.0625 0.1925 False True
sel code, range all, mouse subject 0.2475 0.3725 True True 0.0575 0.1775 False True
sel code, range all, mouse animal 0.2325 0.3725 True True 0.0575 0.1775 False True
sel methods, range unit, mouse subject 0.2225 0.3975 True True 0.0625 0.1925 False True
sel methods, range unit, mouse animal 0.2225 0.3975 True True 0.0625 0.1925 False True
sel methods, range all, mouse subject 0.2225 0.3825 True True 0.0625 0.1675 False True
sel methods, range all, mouse animal 0.2225 0.3825 True True 0.0625 0.1675 False True
Variants:
variant type below 0.5 s below 0.75 s
code (delay epoch reference, 100 iterations) PTupper 0.3825 0.1975
code (delay epoch reference, 100 iterations) PTlower 0.0725 0.0575
processed data PTupper 0.3825 0.1875
processed data PTlower 0.2375 0.0625
(38 s)
fig, ax = plt.subplots(figsize=(6, 3.2))
for ct in alm.TYPES:
med = curves[(ct, "code", "unit", "subject")]
ax.plot(tsub, med, color=colours[ct], lw=1.5, label=ct)
for level in (0.5, 0.75):
ax.axhline(level, color="#999999", lw=0.8, ls=":")
ax.axvline(0, color="#444444", lw=0.8)
ax.set_xlabel("Time from go cue (s)")
ax.set_ylabel("Correlation with CDlate (normalised)")
ax.legend(frameon=False, fontsize=8)
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
plt.show()
for ct in alm.TYPES:
med = curves[(ct, "code", "unit", "subject")]
print(f"{ct}: median normalised correlation from {med.min():.2f} to {med.max():.2f}, all within the axes")
PTupper: median normalised correlation from -0.08 to 1.02, all within the axes PTlower: median normalised correlation from 0.05 to 1.02, all within the axes
Result: blocked (methods underspecified). With the code's selection and the unit's range, the correlation with CDlate falls sooner in PTlower than in PTupper under both levels and both readings of what a mouse is. At level 0.5, PTlower falls below it at 0.37 s (4 and 8 mice) or 0.27 s (the code's animals) and PTupper at 0.37 s, so the claim holds (reproduced), though with the 4 and 8 mice the two differ by one 5 ms bin. At level 0.75, PTlower falls below it at 0.06 s and PTupper at 0.19 s, under the 0.2 s that "several hundred milliseconds" asks (partially reproduced). The answer kept both levels, as the paper and its code fix none, so the claim stays blocked. The code's variant (PTlower 0.07 and 0.06 s, PTupper 0.38 and 0.20 s) and the paper's units (0.24 and 0.06 s against 0.38 and 0.19 s) give the same order. The crossing times are taken from the bootstrap median with no interval, so the one-bin difference at level 0.5 has no measure of uncertainty.
pt-lower-movement-termination¶
In PTlower neurons, selectivity along both CDlate and CDgo ceases with the offset of movement (the last lick), together with an abrupt change in the coding direction; this transition is not seen in the PTupper population.
Criterion: rates aligned to the last lick-port contact after the go cue on correct trials, within each unit's trial range, projected on the go-cue CDlate and CDgo without baseline subtraction, shared bootstrap. Readings: (a) the 5th percentile is above zero for most bins of the 200 ms before the last lick and for none of the bins from 0.2 to 0.4 s after it; (b) the same with 400 ms before and 0.1 to 0.5 s after. Reproduced when this holds for PTlower along both CDs; partially reproduced for one CD. The criterion fixes the unit's trial range, so the runs with that range are used. The paper gives no windows, so both stay (claims file, from a reader's answer). The last lick is the code's (getLickInfo), the latest left or right lick after the go cue in the trial. Here it is rebuilt from the files (data check), a deviation that the criterion accepts when the outcome is the same with the trials whose rebuilt last lick disagrees with the processed data left out; the second cell shows this.
WINDOWS = {"(a)": ((-0.2, 0.0), (0.2, 0.4)), "(b)": ((-0.4, 0.0), (0.1, 0.5))}
def termination(s, w):
(b0, b1), (a0, a1) = WINDOWS[w]
before = (alm.TIME > b0) & (alm.TIME <= b1)
aft = (alm.TIME >= a0) & (alm.TIME <= a1)
frac_before = float((s["q5"][before] > 0).mean())
n_after = int((s["q5"][aft] > 0).sum())
return frac_before, n_after, frac_before > 0.5 and n_after == 0
rows = []
for key in [k for k in BOOT_KEYS if k[1] == "unit"]:
for ct in alm.TYPES:
for cd in ("late", "go"):
s = sel_summary(ct, *key, cd, align="lick")
row = {"reading": reading_label(*key), "type": ct, "CD": cd}
for w in WINDOWS:
fb, na, ok = termination(s, w)
row[f"{w} before: share significant"], row[f"{w} after: significant bins"], row[f"{w} holds"] = fb, na, ok
rows.append(row)
term_table = pd.DataFrame(rows)
print(term_table.round(3).to_string(index=False))
rows = []
for mk in ("subject", "animal"):
for pr in ("proj a", "proj b"):
for ct in alm.TYPES:
for cd in ("late", "go"):
s = sel_summary(ct, None, None, mk, pr, cd, align="lick", processed=True)
row = {"reading": f"processed data, mouse {mk}, {pr}", "type": ct, "CD": cd}
for w in WINDOWS:
fb, na, ok = termination(s, w)
row[f"{w} before"], row[f"{w} after bins"], row[f"{w} holds"] = fb, na, ok
rows.append(row)
print("\nVariant, the paper's units with their own last licks (processed data):")
print(pd.DataFrame(rows).round(3).to_string(index=False))
reading type CD (a) before: share significant (a) after: significant bins (a) holds (b) before: share significant (b) after: significant bins (b) holds
sel code, range unit, mouse subject, proj a PTupper late 0.000 0 False 0.000 0 False
sel code, range unit, mouse subject, proj a PTupper go 1.000 9 False 1.000 13 False
sel code, range unit, mouse subject, proj a PTlower late 1.000 9 False 1.000 29 False
sel code, range unit, mouse subject, proj a PTlower go 1.000 0 True 1.000 15 False
sel code, range unit, mouse subject, proj b PTupper late 0.000 0 False 0.000 0 False
sel code, range unit, mouse subject, proj b PTupper go 0.825 4 False 0.912 4 False
sel code, range unit, mouse subject, proj b PTlower late 1.000 8 False 1.000 26 False
sel code, range unit, mouse subject, proj b PTlower go 1.000 0 True 1.000 13 False
sel code, range unit, mouse animal, proj a PTupper late 0.000 0 False 0.000 0 False
sel code, range unit, mouse animal, proj a PTupper go 1.000 9 False 1.000 13 False
sel code, range unit, mouse animal, proj a PTlower late 1.000 9 False 0.962 25 False
sel code, range unit, mouse animal, proj a PTlower go 1.000 16 False 1.000 37 False
sel code, range unit, mouse animal, proj b PTupper late 0.000 0 False 0.000 0 False
sel code, range unit, mouse animal, proj b PTupper go 0.825 4 False 0.912 4 False
sel code, range unit, mouse animal, proj b PTlower late 0.925 8 False 0.925 23 False
sel code, range unit, mouse animal, proj b PTlower go 1.000 11 False 1.000 31 False
sel methods, range unit, mouse subject, proj a PTupper late 0.000 0 False 0.000 0 False
sel methods, range unit, mouse subject, proj a PTupper go 1.000 10 False 1.000 15 False
sel methods, range unit, mouse subject, proj a PTlower late 0.750 7 False 0.750 16 False
sel methods, range unit, mouse subject, proj a PTlower go 1.000 2 False 1.000 18 False
sel methods, range unit, mouse subject, proj b PTupper late 0.000 0 False 0.000 0 False
sel methods, range unit, mouse subject, proj b PTupper go 0.925 8 False 0.962 9 False
sel methods, range unit, mouse subject, proj b PTlower late 0.600 5 False 0.662 8 False
sel methods, range unit, mouse subject, proj b PTlower go 1.000 1 False 1.000 13 False
sel methods, range unit, mouse animal, proj a PTupper late 0.000 0 False 0.000 0 False
sel methods, range unit, mouse animal, proj a PTupper go 1.000 10 False 1.000 15 False
sel methods, range unit, mouse animal, proj a PTlower late 0.675 8 False 0.588 21 False
sel methods, range unit, mouse animal, proj a PTlower go 1.000 13 False 1.000 33 False
sel methods, range unit, mouse animal, proj b PTupper late 0.000 0 False 0.000 0 False
sel methods, range unit, mouse animal, proj b PTupper go 0.925 8 False 0.962 9 False
sel methods, range unit, mouse animal, proj b PTlower late 0.500 7 False 0.488 15 False
sel methods, range unit, mouse animal, proj b PTlower go 1.000 8 False 1.000 26 False
Variant, the paper's units with their own last licks (processed data):
reading type CD (a) before (a) after bins (a) holds (b) before (b) after bins (b) holds
processed data, mouse subject, proj a PTupper late 0.0 0 False 0.000 0 False
processed data, mouse subject, proj a PTupper go 0.0 0 False 0.000 0 False
processed data, mouse subject, proj a PTlower late 1.0 8 False 1.000 24 False
processed data, mouse subject, proj a PTlower go 1.0 0 True 1.000 13 False
processed data, mouse subject, proj b PTupper late 0.0 0 False 0.000 0 False
processed data, mouse subject, proj b PTupper go 0.0 0 False 0.000 0 False
processed data, mouse subject, proj b PTlower late 1.0 7 False 1.000 23 False
processed data, mouse subject, proj b PTlower go 1.0 0 True 1.000 12 False
processed data, mouse animal, proj a PTupper late 0.0 0 False 0.000 0 False
processed data, mouse animal, proj a PTupper go 0.0 0 False 0.000 0 False
processed data, mouse animal, proj a PTlower late 1.0 8 False 1.000 23 False
processed data, mouse animal, proj a PTlower go 1.0 15 False 1.000 35 False
processed data, mouse animal, proj b PTupper late 0.0 0 False 0.000 0 False
processed data, mouse animal, proj b PTupper go 0.0 0 False 0.000 0 False
processed data, mouse animal, proj b PTlower late 1.0 7 False 0.975 19 False
processed data, mouse animal, proj b PTlower go 1.0 11 False 1.000 31 False
# The rebuilt last lick (criterion): the outcome with the trials whose rebuilt last lick disagrees with
# the processed data's (more than 5 ms apart, or missing) left out of every unit, coding directions and
# projections alike; selection (a), range (a).
t0 = time.time()
def lick_checked(u):
agree = lick_agree[(u["subject"], u["session"])][u["trial_index"]]
return {**u, "in_range": u["in_range"] & agree}
rows = []
for ct in alm.TYPES:
pop = alm.selected(units, ct, "code")
kept = [lick_checked(u) for u in pop]
n_all = sum(int(u["in_range"].sum()) for u in pop)
n_kept = sum(int(u["in_range"].sum()) for u in kept)
print(f"{ct}: trials in the units' ranges {n_all}, left out {n_all - n_kept} ({(n_all - n_kept) / n_all:.1%}); "
f"fewest trials of one type in a unit {min(min(len(r) for r in alm.trial_rows(u, 'unit')) for u in kept)}")
for mk in ("subject", "animal"):
if ct == "PTupper" and mk == "animal":
continue
b = alm.summarise(alm.bootstrap(kept, mk, "unit", alm.seeded(f"last lick checked {ct} {mk}"), lick=True))
for pr in ("proj a", "proj b"):
for cd in ("late", "go"):
row = {"reading": f"mouse {mk}, {pr}", "type": ct, "CD": cd}
for w in WINDOWS:
fb, na, ok = termination(b[("lick", pr, cd)], w)
full = term_table[(term_table.reading == reading_label("code", "unit", mk, pr)) & (term_table.type == ct) & (term_table.CD == cd)].iloc[0]
row[f"{w} before"], row[f"{w} after bins"], row[f"{w} holds"] = fb, na, ok
row[f"{w} holds, all trials"] = bool(full[f"{w} holds"])
rows.append(row)
checked = pd.DataFrame(rows)
print(checked.round(3).to_string(index=False))
print("Same holds/fails as with all trials, every row and window:",
bool(all((checked[f"{w} holds"] == checked[f"{w} holds, all trials"]).all() for w in WINDOWS)))
print(f"({time.time() - t0:.0f} s)")
PTupper: trials in the units' ranges 9252, left out 158 (1.7%); fewest trials of one type in a unit 61
PTlower: trials in the units' ranges 12566, left out 257 (2.0%); fewest trials of one type in a unit 56
reading type CD (a) before (a) after bins (a) holds (a) holds, all trials (b) before (b) after bins (b) holds (b) holds, all trials mouse subject, proj a PTupper late 0.000 0 False False 0.000 0 False False mouse subject, proj a PTupper go 1.000 10 False False 1.000 18 False False mouse subject, proj b PTupper late 0.000 0 False False 0.000 0 False False mouse subject, proj b PTupper go 0.900 4 False False 0.950 4 False False mouse subject, proj a PTlower late 1.000 8 False False 0.975 26 False False mouse subject, proj a PTlower go 1.000 1 False True 1.000 16 False False mouse subject, proj b PTlower late 1.000 7 False False 0.962 20 False False mouse subject, proj b PTlower go 1.000 0 True True 1.000 12 False False mouse animal, proj a PTlower late 0.850 8 False False 0.812 20 False False mouse animal, proj a PTlower go 1.000 18 False False 1.000 39 False False mouse animal, proj b PTlower late 0.825 7 False False 0.800 15 False False mouse animal, proj b PTlower go 1.000 14 False False 1.000 34 False False Same holds/fails as with all trials, every row and window: False (12 s)
# How often the lick-aligned rates are cut: a trial's spikes reach 3.0 s after its go cue, so bins
# later than that after the last lick count no spikes. Share of correct trials in the included units'
# ranges whose rebuilt last lick is more than 2.5, 2.6 or 2.8 s after the go cue (bins up to 0.5, 0.4
# and 0.2 s after it are then cut).
for ct in alm.TYPES:
rel = np.concatenate([u["last_lick_rel"][u["in_range"]] for u in alm.selected(units, ct, "code")])
rel = rel[np.isfinite(rel)]
print(f"{ct}: {len(rel)} unit-trials; last lick median {np.median(rel):.3f} s, 95th percentile "
f"{np.percentile(rel, 95):.3f} s after the go cue; later than 2.5 / 2.6 / 2.8 s: "
+ " / ".join(f"{np.mean(rel > x):.2%}" for x in (2.5, 2.6, 2.8)))
PTupper: 9241 unit-trials; last lick median 1.414 s, 95th percentile 1.800 s after the go cue; later than 2.5 / 2.6 / 2.8 s: 0.16% / 0.14% / 0.13% PTlower: 12545 unit-trials; last lick median 1.365 s, 95th percentile 1.734 s after the go cue; later than 2.5 / 2.6 / 2.8 s: 0.17% / 0.13% / 0.10%
# The coding direction around the last lick (the paper gives no measure of the abrupt change):
# per bin, d' over all trials of each unit's last-lick aligned rates, unit-normalised; inner products
# between the mean direction 0.3 to 0.1 s before the last lick and each later bin.
def lick_cd(pop):
t_sel = (alm.TIME >= -1.0) & (alm.TIME <= 1.0)
vecs = []
for u in pop:
rl, rr = alm.trial_rows(u, "unit")
sl = alm.smooth(u["lick"][rl].toarray())[:, t_sel] / alm.DT
sr = alm.smooth(u["lick"][rr].toarray())[:, t_sel] / alm.DT
den = np.sqrt(sl.var(axis=0, ddof=1) + sr.var(axis=0, ddof=1))
vecs.append(np.where(den > 0, (sr.mean(axis=0) - sl.mean(axis=0)) / np.where(den > 0, den, 1), 0))
v = np.array(vecs)
v = v / (np.linalg.norm(v, axis=0, keepdims=True) + 1e-30)
return alm.TIME[t_sel], v
for ct in alm.TYPES:
t_l, v = lick_cd(alm.selected(units, ct, "code"))
ref = v[:, (t_l >= -0.3) & (t_l <= -0.1)].mean(axis=1)
ref /= np.linalg.norm(ref)
corr = ref @ v
print(f"{ct}: inner product with the direction 0.3 to 0.1 s before the last lick, at "
+ ", ".join(f"{x:+.1f} s {corr[np.argmin(abs(t_l - x))]:.2f}" for x in (-0.2, 0.0, 0.1, 0.2, 0.3, 0.5)))
PTupper: inner product with the direction 0.3 to 0.1 s before the last lick, at -0.2 s 0.95, +0.0 s 0.89, +0.1 s 0.88, +0.2 s 0.82, +0.3 s 0.79, +0.5 s 0.68
PTlower: inner product with the direction 0.3 to 0.1 s before the last lick, at -0.2 s 0.98, +0.0 s 0.94, +0.1 s 0.84, +0.2 s 0.70, +0.3 s 0.52, +0.5 s 0.40
Result: blocked (methods underspecified). With the code's selection and the unit's range, PTlower CDlate selectivity is significant in 92% to 100% of the bins before the last lick, but it stays significant after it: 8 or 9 bins from 0.2 to 0.4 s and 23 to 29 from 0.1 to 0.5 s, so it does not cease under either window. Along CDgo it ceases under window (a) with the 4 and 8 mice (no significant bin after, under both projections), which gives partially reproduced; with the code's animals 11 or 16 bins stay significant under window (a), and under window (b) 13 to 37 stay, which gives not reproduced. The window and the reading of what a mouse is decide, and the answer kept both windows. With the trials whose rebuilt last lick disagrees with the processed data left out (1.7% of PTupper and 2.0% of PTlower trials), the outcome is the same in every combination but one: with the 4 and 8 mice and projection (a), 1 bin of PTlower CDgo stays significant from 0.2 to 0.4 s after the last lick, where none did with all trials. So in that combination the rebuilt last lick is not shown to leave the outcome unchanged, as the criterion asks. PTupper CDlate selectivity is not significant before the last lick under any reading. The PTlower coding direction changes more after the last lick than PTupper's: its inner product with the direction before the last lick is 0.52 at +0.3 s, against 0.79. For PTlower, the paper's units with their own last licks give the same pattern as the dandiset with all trials. For PTupper they do not: on the dandiset, with the code's selection, PTupper CDgo selectivity is significant in 82.5% to 100% of the bins before the last lick and in 4 to 13 bins after it, while on the paper's units no bin before or after it is significant. A trial's spikes reach 3.0 s after its go cue, so lick-aligned bins later than that count no spikes; the last lick is more than 2.5 s after the go cue in 0.16% of PTupper and 0.17% of PTlower unit-trials (third cell), so this cut does not bear on the result.
pt-cohort-control¶
Differences in coding between PTupper and PTlower neurons are not explained by the animal cohorts; unidentified neurons recorded in the two cohorts give similar coding direction projections, nearly indistinguishable once depth-matched.
Criterion: unidentified units of each cohort (counts 495, 511 and 276 +/- 10%, before or after the inclusion rule), depth-matched by keeping the PTlower cohort's left ALM recordings. Reproduced when, for the depth-matched populations, the PTlower-cohort population's mean CDearly selectivity over -2.5 to -2.1 s lies inside the PTupper-cohort population's 5-95% bootstrap interval, and the two CDgo onsets (pt-cdgo-onset rule) differ by no more than the difference between the PT types; partially reproduced when one holds. The z-scored variant of the PT populations is shown.
t0 = time.time()
print("Unidentified units (good quality / reading (a) / reading (b)):")
for coh, left, paper in (("PTupper cohort", False, 495), ("PTlower cohort", False, 511), ("PTlower cohort", True, 276)):
g = udf[(udf.cell_type == "unidentified") & (udf.cohort == coh) & ((udf.hemisphere == "left") | (not left))]
print(f" {coh}{', left ALM only' if left else ''}: {len(g)} / {int(g.included.sum())} / {int(g.session_ok.sum())} "
f"(the paper {paper}; within 10%: {any(abs(x - paper) <= 0.1 * paper for x in (len(g), g.included.sum(), g.session_ok.sum()))})")
cohort_boots = {}
for s in ("code", "methods"):
for r in ("unit", "all"):
for name, coh, left, mks in (("PTupper cohort", "PTupper cohort", False, ("subject",)),
("PTlower cohort", "PTlower cohort", False, ("subject", "animal")),
("PTlower cohort, left", "PTlower cohort", True, ("subject",))):
pop = alm.selected(units, "unidentified", s, coh, left_only=left)
for mk in mks:
key = (name, s, r, mk)
cohort_boots[key] = alm.summarise(alm.bootstrap(pop, mk, r, alm.seeded("cohort " + " ".join(key))))
print(f"cohort bootstraps in {time.time() - t0:.0f} s")
Unidentified units (good quality / reading (a) / reading (b)): PTupper cohort: 232 / 226 / 232 (the paper 495; within 10%: False) PTlower cohort: 430 / 403 / 430 (the paper 511; within 10%: False) PTlower cohort, left ALM only: 226 / 212 / 226 (the paper 276; within 10%: False)
cohort bootstraps in 196 s
rows = []
for s in ("code", "methods"):
for r in ("unit", "all"):
for pr in ("proj a", "proj b"):
up = cohort_boots[("PTupper cohort", s, r, "subject")][("go", pr, "early")]
low = cohort_boots[("PTlower cohort, left", s, r, "subject")][("go", pr, "early")]
lo_all = cohort_boots[("PTlower cohort", s, r, "subject")][("go", pr, "early")]
q5, q95 = np.percentile(up["means"]["CDearly window"], [5, 95])
m_low = np.median(low["means"]["CDearly window"])
row = {"reading": f"sel {s}, range {r}, {pr}", "PTupper cohort 5%": q5, "95%": q95,
"PTlower cohort left median": m_low, "PTlower cohort all median": np.median(lo_all["means"]["CDearly window"]),
"early inside": bool(q5 <= m_low <= q95)}
gu = cohort_boots[("PTupper cohort", s, r, "subject")][("go", pr, "go")]
gl = cohort_boots[("PTlower cohort, left", s, r, "subject")][("go", pr, "go")]
ok = []
for st, k in STAMPS:
d_unid = abs(alm.onset(gu["q5"], st, k) - alm.onset(gl["q5"], st, k))
d_pt = abs(alm.onset(sel_summary("PTupper", s, r, "subject", pr, "go")["q5"], st, k)
- alm.onset(sel_summary("PTlower", s, r, "subject", pr, "go")["q5"], st, k))
ok.append(d_unid <= d_pt)
if (st, k) == ("centre", 1):
row["unidentified onsets ms (PTupper cohort, PTlower left)"] = (1000 * alm.onset(gu["q5"]), 1000 * alm.onset(gl["q5"]))
row["PT onset difference ms"] = 1000 * d_pt
row["go onsets within PT difference (all stamps)"] = all(ok)
row["go onsets within PT difference (any stamp)"] = any(ok)
rows.append(row)
cohort_table = pd.DataFrame(rows)
print(cohort_table.round(3).to_string(index=False))
rows = []
for s in ("code", "methods"):
for r in ("unit", "all"):
for pr in ("proj a", "proj b"):
a = cohort_boots[("PTlower cohort", s, r, "subject")][("go", pr, "early")]["means"]["CDearly window"]
b = cohort_boots[("PTlower cohort", s, r, "animal")][("go", pr, "early")]["means"]["CDearly window"]
rows.append({"reading": f"sel {s}, range {r}, {pr}", "PTlower cohort (all) median, mouse subject": np.median(a),
"mouse animal": np.median(b)})
print("\nPTlower cohort, both hemispheres, both readings of what a mouse is:")
print(pd.DataFrame(rows).round(3).to_string(index=False))
reading PTupper cohort 5% 95% PTlower cohort left median PTlower cohort all median early inside unidentified onsets ms (PTupper cohort, PTlower left) PT onset difference ms go onsets within PT difference (all stamps) go onsets within PT difference (any stamp)
sel code, range unit, proj a 3.965 8.857 9.822 7.269 False (42.5, 37.5) 40.0 True True
sel code, range unit, proj b 4.202 9.291 10.072 7.615 False (42.5, 37.5) 45.0 True True
sel code, range all, proj a 4.014 8.496 8.817 6.589 False (42.5, 37.5) 45.0 True True
sel code, range all, proj b 4.282 9.027 9.056 7.017 False (42.5, 37.5) 45.0 True True
sel methods, range unit, proj a 4.044 8.833 9.413 6.924 False (42.5, 37.5) 40.0 True True
sel methods, range unit, proj b 4.322 9.294 9.639 7.249 False (42.5, 37.5) 35.0 True True
sel methods, range all, proj a 3.839 8.336 8.407 6.294 False (42.5, 37.5) 40.0 True True
sel methods, range all, proj b 4.104 8.819 8.665 6.674 True (42.5, 37.5) 40.0 True True
PTlower cohort, both hemispheres, both readings of what a mouse is:
reading PTlower cohort (all) median, mouse subject mouse animal
sel code, range unit, proj a 7.269 7.503
sel code, range unit, proj b 7.615 7.808
sel code, range all, proj a 6.589 6.535
sel code, range all, proj b 7.017 6.932
sel methods, range unit, proj a 6.924 7.236
sel methods, range unit, proj b 7.249 7.577
sel methods, range all, proj a 6.294 6.343
sel methods, range all, proj b 6.674 6.738
# Variant: PT populations with rates z-scored on the epoch before the stimulus (Extended Data Fig. 8,
# bottom row). Each unit's rates are divided by the s.d. over its trials of the baseline epoch rate
# (a standard choice; the legend does not say which s.d.); the mean cancels in the selectivity.
rows = []
for ct in alm.TYPES:
pop = alm.selected(units, ct, "code")
sd = []
for u in pop:
rl, rr = alm.trial_rows(u, "unit")
x = alm.trial_means(u["go"], np.r_[rl, rr], alm.EPOCHS["baseline"])
sd.append(1 / x.std(ddof=1) if x.std(ddof=1) > 0 else 0.0)
z = alm.summarise(alm.bootstrap(pop, "subject", "unit", alm.seeded(f"z-scored {ct}"), scales=np.array(sd)))
rows.append({"type": ct, "CDearly window median": np.median(z[("go", "proj a", "early")]["means"]["CDearly window"]),
"5%": np.percentile(z[("go", "proj a", "early")]["means"]["CDearly window"], 5),
"CDgo onset ms": 1000 * alm.onset(z[("go", "proj a", "go")]["q5"])})
print("z-scored PT populations (selection (a), range (a), mouse (a), projection (a)):")
print(pd.DataFrame(rows).round(3).to_string(index=False))
z-scored PT populations (selection (a), range (a), mouse (a), projection (a)): type CDearly window median 5% CDgo onset ms PTupper 0.831 0.439 92.5 PTlower 0.290 0.189 37.5
Result: blocked (methods underspecified), as it depends on pt-upper-early-preparatory. Its own test, with the code's selection and the unit's range, gives partially reproduced under both projections. The CDgo onsets of the depth-matched unidentified populations (42.5 and 37.5 ms) differ by 5 ms, less than the 40 and 45 ms between the PT types, under every time stamp. But the PTlower cohort's left ALM population has a larger CDearly selectivity (median 9.8 and 10.1) than the PTupper cohort, outside its 5-95% interval (4.0 to 8.9, and 4.2 to 9.3). The dandiset holds 232 unidentified good-quality units of the PTupper cohort, 430 of the PTlower cohort and 226 of its left ALM recordings (226, 403 and 212 after the code's rule), against the paper's 495, 511 and 276, all out of tolerance, as the units of the missing sessions are missing. The z-scored PT populations keep the PT difference (CDearly 0.83 against 0.29; CDgo onset 92.5 against 37.5 ms). The answer fixed the selection, which leaves out the one combination of the earlier run where the CDearly part held (session rule, all trials, projection b); the outcome is unchanged. The CDearly part fails because the PTlower cohort is more selective than the PTupper cohort, while among the PT types PTupper is the more selective (3.97 against 2.09).
pt-types-differ-in-timing¶
The two PT types have different roles in time; PTupper populations carry early preparatory activity that persists until movement, while PTlower populations mainly carry late preparatory activity and movement commands.
Criterion: the summary of the six claims it depends on (pt-upper-early-preparatory, pt-upper-persistent, pt-lower-cdearly-decay, pt-single-neuron-early-selectivity, pt-lower-late-and-motor, pt-cdgo-onset). Reproduced when all six are reproduced; partially reproduced when at least one claim on each side (PTupper early and persistent; PTlower late and motor) is reproduced; not reproduced otherwise.
print("The six claims, main reading (selection (a), range (a), mouse (a), projection (a)), as tested above:")
main_key = ("code", "unit", "subject", "proj a")
e = early_table[early_table.reading == reading_label(*main_key)]
print(" pt-upper-early-preparatory: larger (i)", e["larger (i)"].tolist(), "; more consistent", e["more consistent"].tolist())
pr_ = persistent_row(sel_summary("PTupper", *main_key, "early"), "main")
print(" pt-upper-persistent: (a)", pr_["(a) every bin -2.1 to 0 s"], "; (b) p", round(pr_["(b) p, -0.4 to 0 s"], 3))
f = pop_table(alm.TYPES, "code", "unit")
up_, low_ = f[f.cell_type == "PTupper"]["early d'"].abs(), f[f.cell_type == "PTlower"]["early d'"].abs()
print(f" pt-single-neuron-early-selectivity, measure (a): {up_.mean():.3f} against {low_.mean():.3f}; t p "
f"{stats.ttest_ind(up_, low_).pvalue:.4f}; Mann-Whitney p {stats.mannwhitneyu(up_, low_, alternative='two-sided').pvalue:.4f}")
print(" pt-lower-cdearly-decay: p", round(float((sel_summary("PTlower", *main_key, "early")["means"]["late delay"] <= 0).mean()), 3))
print(" pt-lower-late-and-motor:", motor_table[motor_table.reading == reading_label(*main_key)][["CDlate (i)", "larger (a)", "larger (b)", "faster (all stamps)", "longer"]].to_dict("records"))
print(" pt-cdgo-onset:", onset_table[(onset_table.reading == reading_label(*main_key)) & (onset_table.stamp == "centre") & (onset_table.consecutive == 1)][["PTlower ms", "PTupper ms"]].to_dict("records"))
The six claims, main reading (selection (a), range (a), mouse (a), projection (a)), as tested above:
pt-upper-early-preparatory: larger (i) [True, True] ; more consistent [False, False]
pt-upper-persistent: (a) False ; (b) p 0.053
pt-single-neuron-early-selectivity, measure (a): 0.410 against 0.237; t p 0.0012; Mann-Whitney p 0.0041
pt-lower-cdearly-decay: p 0.348
pt-lower-late-and-motor: [{'CDlate (i)': True, 'larger (a)': True, 'larger (b)': True, 'faster (all stamps)': True, 'longer': True}]
pt-cdgo-onset: [{'PTlower ms': 37.5, 'PTupper ms': 77.5}]
Result: blocked (methods underspecified). Two of the six claims it summarises are blocked: pt-upper-early-preparatory, and pt-upper-persistent through its dependency on the first, so this claim is blocked for the same reason. Of the others, pt-single-neuron-early-selectivity, pt-lower-cdearly-decay and pt-lower-late-and-motor are reproduced and pt-cdgo-onset is partially reproduced. pt-upper-persistent's own test does not hold with the code's selection and the unit's range (p = 0.053 for the mean over the last 400 ms). The answers changed pt-single-neuron-early-selectivity from blocked to reproduced; this claim's outcome is unchanged. Without the dependency rule the criterion would give partially reproduced, with pt-single-neuron-early-selectivity reproduced on the PTupper side and pt-lower-late-and-motor on the PTlower side; the claims file makes the two blocked claims part of this one, so it is blocked.
Missing units against the file format (variant)¶
Two things separate the dandiset's results from the paper's: the 11 missing sessions, and rates
rebuilt from the files' spike times on the file's clock (where spikes of a neighbouring trial are
counted). The paper's processed data restricted to the tagged units the dandiset also holds
(per_unit_same, boots_same) have the dandiset's units and the paper's per-trial rates. A
difference between the dandiset and these comes from the files; a difference between these and
all the paper's units comes from the missing units. Selection (a), range (a), mouse (a). No outcome
rests on this section.
frames3 = [("dandiset files", pop_table(alm.TYPES, "code", "unit")),
("processed, same units", pop_table(alm.TYPES, None, "unit", per_unit_same)),
("processed, all units", pop_table(alm.TYPES, None, "unit", per_unit_proc))]
rows = []
for label, f in frames3:
up, low = f[f.cell_type == "PTupper"], f[f.cell_type == "PTlower"]
r = {"population": label, "units (PTupper, PTlower)": (len(up), len(low))}
for m in ("baseline rate", "peak rate"):
r[f"{m}: medians Hz"] = (round(up[m].median(), 2), round(low[m].median(), 2))
r[f"{m}: t-test p"] = round(stats.ttest_ind(up[m], low[m]).pvalue, 4)
b = f["burst fraction"] > alm.BURST_FRACTION
r["bursting (PTupper, PTlower)"] = (int(b[up.index].sum()), int(b[low.index].sum()))
sig = (f[["p sample", "p delay", "p response"]] < alm.SELECTIVITY_P).any(axis=1)
r["selective in any epoch"] = f"{int(sig.sum())} of {len(f)} ({sig.mean():.1%})"
late = f["p late delay"] < alm.SELECTIVITY_P
r["late delay selective"] = (f"{int(late[up.index].sum())} of {len(up)} ({late[up.index].mean():.1%})",
f"{int(late[low.index].sum())} of {len(low)} ({late[low.index].mean():.1%})")
a_up, a_low = up["early d'"].abs(), low["early d'"].abs()
r["early |d'| means"] = (round(a_up.mean(), 3), round(a_low.mean(), 3))
r["early |d'|: t-test p"] = round(stats.ttest_ind(a_up, a_low).pvalue, 4)
for col, n_top in (("accuracy", 10), ("rectified accuracy", 21)):
d = decoding_rows(f, label, col, n_top)
r[f"{col}: >= 70% (PTupper, PTlower)"] = (d["PTupper >= 0.70"], d["PTlower >= 0.70"])
r[f"{col}: means %"] = (round(d["PTupper mean %"], 2), round(d["PTlower mean %"], 2))
r[f"{col}: PTupper in top {n_top}"] = d[f"PTupper in top {n_top}"]
rows.append(r)
print("Per-unit measures:")
print(pd.DataFrame(rows).set_index("population").T.to_string())
def boot_get(population, ct, pr, cd):
if population == "dandiset files":
return sel_summary(ct, "code", "unit", "subject", pr, cd)
if population == "processed, all units":
return sel_summary(ct, None, None, "subject", pr, cd, processed=True)
return boots_same[ct][("go", pr, cd)]
rows = []
for label, _ in frames3:
for pr in ("proj a", "proj b"):
e = early_rows(lambda ct: boot_get(label, ct, pr, "early"), label)[0]
per = persistent_row(boot_get(label, "PTupper", pr, "early"), label)
lo, upg = boot_get(label, "PTlower", pr, "go"), boot_get(label, "PTupper", pr, "go")
rows.append({"population": label, "projection": pr,
"CDearly medians (PTupper, PTlower)": (round(e["PTupper median"], 2), round(e["PTlower median"], 2)),
"larger (i)": e["larger (i)"], "larger (ii)": e["larger (ii)"],
"sd/median (PTupper, PTlower)": (round(e["PTupper sd/median"], 3), round(e["PTlower sd/median"], 3)),
"persistent: bins not significant": per["bins not significant"],
"persistent (b) p": round(per["(b) p, -0.4 to 0 s"], 3),
"PTlower decay p": round(float((boot_get(label, "PTlower", pr, "early")["means"]["late delay"] <= 0).mean()), 3),
"CDgo onsets ms (PTlower, PTupper)": (1000 * alm.onset(lo["q5"]), 1000 * alm.onset(upg["q5"])),
"CDgo run end s (PTlower, PTupper)": (alm.first_run_end(lo["q5"]), alm.first_run_end(upg["q5"]))})
print("\nShared bootstrap (mouse (a); bin centres, first significant bin):")
print(pd.DataFrame(rows).to_string(index=False))
# The same units in the files and in the processed data: how much the file's clock moves each unit's rates.
cols = ["baseline rate", "peak rate", "accuracy", "rectified accuracy", "early d'"]
same = (pop_table(alm.TYPES, "code", "unit")[["subject", "session", "id", "cell_type"] + cols]
.merge(pop_table(alm.TYPES, None, "unit", per_unit_same)[["subject", "session", "id"] + cols],
on=["subject", "session", "id"], suffixes=(" files", " processed")))
print(f"\nThe same units, files against processed data ({len(same)} units matched):")
rows = []
for ct in alm.TYPES:
g = same[same.cell_type == ct]
for c in cols:
d = g[f"{c} files"] - g[f"{c} processed"]
rows.append({"type": ct, "measure": c, "median files": g[f"{c} files"].median(),
"median processed": g[f"{c} processed"].median(), "median difference": d.median(),
"min difference": d.min(), "max difference": d.max(), "units higher in files": int((d > 0).sum()),
"of": len(g)})
print(pd.DataFrame(rows).round(3).to_string(index=False))
Per-unit measures:
population dandiset files processed, same units processed, all units
units (PTupper, PTlower) (44, 56) (44, 56) (61, 69)
baseline rate: medians Hz (4.93, 6.69) (4.26, 5.66) (4.27, 5.18)
baseline rate: t-test p 0.0742 0.0847 0.1782
peak rate: medians Hz (18.37, 20.96) (17.94, 19.41) (17.65, 17.17)
peak rate: t-test p 0.0302 0.0325 0.1418
bursting (PTupper, PTlower) (2, 11) (2, 11) (2, 13)
selective in any epoch 99 of 100 (99.0%) 97 of 100 (97.0%) 121 of 130 (93.1%)
late delay selective (29 of 44 (65.9%), 38 of 56 (67.9%)) (29 of 44 (65.9%), 38 of 56 (67.9%)) (40 of 61 (65.6%), 44 of 69 (63.8%))
early |d'| means (0.41, 0.237) (0.409, 0.256) (0.414, 0.247)
early |d'|: t-test p 0.0012 0.0044 0.0004
accuracy: >= 70% (PTupper, PTlower) (9, 3) (9, 3) (15, 3)
accuracy: means % (64.28, 58.71) (64.2, 58.93) (64.38, 58.82)
accuracy: PTupper in top 10 8 8 10
rectified accuracy: >= 70% (PTupper, PTlower) (8, 0) (9, 0) (15, 0)
rectified accuracy: means % (61.81, 56.2) (62.33, 56.84) (62.73, 56.6)
rectified accuracy: PTupper in top 21 16 15 20
Shared bootstrap (mouse (a); bin centres, first significant bin):
population projection CDearly medians (PTupper, PTlower) larger (i) larger (ii) sd/median (PTupper, PTlower) persistent: bins not significant persistent (b) p PTlower decay p CDgo onsets ms (PTlower, PTupper) CDgo run end s (PTlower, PTupper)
dandiset files proj a (3.97, 2.09) True False (0.314, 0.249) 68 0.053 0.348 (37.5, 77.5) (1.4425, 1.3075)
dandiset files proj b (4.12, 2.46) True False (0.32, 0.296) 160 0.058 0.377 (37.5, 82.5) (1.4375, 1.2325)
processed, same units proj a (3.89, 1.84) True True (0.323, 0.259) 85 0.060 0.474 (37.5, 82.5) (1.5625, 1.3025)
processed, same units proj b (4.01, 2.19) True False (0.333, 0.28) 153 0.066 0.493 (37.5, 87.5) (1.4425, 1.2175)
processed, all units proj a (3.94, 1.69) True True (0.221, 0.234) 0 0.007 0.458 (37.5, 57.5) (1.4575, 0.7275)
processed, all units proj b (4.13, 2.04) True True (0.221, 0.261) 0 0.011 0.480 (37.5, 62.5) (1.4425, 0.6025)
The same units, files against processed data (100 units matched):
type measure median files median processed median difference min difference max difference units higher in files of
PTupper baseline rate 4.933 4.256 0.986 0.166 3.638 44 44
PTupper peak rate 18.369 17.936 0.168 -0.667 3.636 28 44
PTupper accuracy 0.621 0.618 0.000 -0.011 0.020 17 44
PTupper rectified accuracy 0.586 0.594 0.000 -0.060 0.020 13 44
PTupper early d' -0.158 -0.147 0.002 -0.115 0.041 25 44
PTlower baseline rate 6.689 5.663 0.966 0.140 6.473 56 56
PTlower peak rate 20.958 19.414 0.198 -1.649 4.004 35 56
PTlower accuracy 0.568 0.568 0.000 -0.038 0.017 22 56
PTlower rectified accuracy 0.554 0.557 -0.002 -0.057 0.017 19 56
PTlower early d' -0.105 -0.139 0.012 -0.040 0.188 40 56
What the split shows. On the same 100 units, the files and the processed data agree on the per-unit measures except the baseline rate: every unit's baseline rate is higher in the files (median by 0.99 Hz for PTupper and 0.97 Hz for PTlower, 0.14 to 6.47 Hz), as spikes of neighbouring trials fall in the baseline epoch, while decoding accuracy differs by a median of 0.000 and early d' by 0.002 and 0.012. The results that fail on the dandiset fail on the same units from the processed data too: peak rates differ (t-test p = 0.03), 9 PTupper units reach 70% accuracy and 8 of the 10 best are PTupper, the late delay counts are 29 and 38, PTupper CDearly selectivity is not significant in 85 and 153 bins before the go cue (p = 0.060 and 0.066 over the last 400 ms), and the PTupper CDgo onset is 82.5 and 87.5 ms. On all the paper's units they hold, so these failures come from the missing units. The files move the rectified PTlower mean (56.20% against 56.84% on the same units), the rectified PTupper count at 70% (8 against 9), rule (ii) of larger CDearly selectivity under projection (a) (fails on the files, holds on the same units processed) and the PTlower CDearly decay p (0.348 against 0.474). The PTlower CDgo onset is 37.5 ms in all three.
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Data | Outcome | Before the answers | Earlier attempt (-01) |
|---|---|---|---|---|
| pt-types-transcriptomic | elsewhere | reproduced | reproduced | blocked (data missing) |
| pt-upper-thalamus-projection | elsewhere | reproduced | blocked (methods underspecified) | blocked (data missing) |
| pt-lower-medulla-projection | elsewhere | reproduced | blocked (methods underspecified) | blocked (data missing) |
| pt-types-shared-targets | elsewhere | blocked (methods underspecified) | reproduced | blocked (data missing) |
| pt-retrograde-sublaminae | elsewhere | blocked (methods underspecified) | blocked (methods underspecified) | new |
| pt-projection-cluster-mapping | elsewhere | reproduced | reproduced | new |
| pt-types-account-for-all-pt | elsewhere | partially reproduced | partially reproduced | new |
| pt-marker-genes | dataset | blocked (data missing) | blocked (other) | new |
| pt-bulk-sc-fold-change-correlation | elsewhere | reproduced | reproduced | new |
| pt-tagged-units | dandiset | not reproduced | not reproduced | new |
| pt-depth | dandiset | reproduced | reproduced | new |
| pt-spike-rates-similar | dandiset | not reproduced | blocked (methods underspecified) | new |
| pt-lower-heterogeneous-rates | dandiset | reproduced | blocked (by dependency) | new |
| pt-lower-bursting | dandiset | partially reproduced | partially reproduced | new |
| pt-selective-any-epoch | dandiset | not reproduced | not reproduced | new |
| pt-variance-explained | dandiset | reproduced | blocked (methods underspecified) | new |
| pt-preparatory-activity | dandiset | partially reproduced | partially reproduced | reproduced |
| pt-upper-early-preparatory | dandiset | blocked (methods underspecified) | blocked (methods underspecified) | reproduced |
| pt-upper-persistent | dandiset | blocked (methods underspecified, by dependency) | blocked (by dependency) | reproduced |
| pt-lower-cdearly-decay | dandiset | reproduced | reproduced | new |
| pt-cd-stability | dandiset | reproduced | reproduced | new |
| pt-single-neuron-early-selectivity | dandiset | reproduced | blocked (methods underspecified) | new |
| pt-upper-decoding | dandiset | partially reproduced | partially reproduced | new |
| pt-upper-decoding-rectified | dandiset | partially reproduced | partially reproduced | new |
| pt-cdgo-onset | dandiset | partially reproduced | partially reproduced | new |
| pt-lower-late-and-motor | dandiset | reproduced | reproduced | partially reproduced |
| reaction-time | dandiset | blocked (data missing) | blocked (data missing) | new |
| pt-cd-change-at-go-cue | dandiset | blocked (methods underspecified) | blocked (methods underspecified) | new |
| pt-lower-movement-termination | dandiset | blocked (methods underspecified) | blocked (methods underspecified) | new |
| pt-cohort-control | dandiset | blocked (methods underspecified, by dependency) | blocked (methods underspecified) | new |
| pt-types-differ-in-timing | dandiset | blocked (methods underspecified, by dependency) | blocked (by dependency) | not reproduced |
"Before the answers" is this attempt's run before a reader answered its questions. The answers settled the unit selection (the code's per-unit rule) and the trial range (the unit's range), so the electrophysiology outcomes are now decided under those, with the other readings as variants; this unblocks pt-spike-rates-similar, pt-variance-explained and pt-single-neuron-early-selectivity. The answers moved the dependency of pt-lower-heterogeneous-rates to pt-tagged-units, named the readings and neurons of the reconstruction claims, made the test of pt-types-shared-targets on the reconstructions a variant, made pt-marker-genes a claim on data the paper deposited nowhere, and kept both rules or levels or windows for pt-upper-early-preparatory, pt-cd-change-at-go-cue and pt-lower-movement-termination, which stay blocked. The earlier attempt's outcomes were reached with criteria written from the abstract; the data are the same version.
How the analysis differs from the paper's. The files hold spike times, not the paper's per-trial rates. The rates are rebuilt on the file's clock: a trial's spikes are those within 3.4 s before and 3.0 s after its go cue, so where trials overlap on that clock, spikes of a neighbouring trial are counted (11.7% of the spikes in the baseline epoch are known to be out of order, under 2% in the later epochs). The last lick of each trial is rebuilt from the stored order of the lick times; it matches the processed data's in 98.1% of correct trials, and leaving out the trials where it does not changes the movement termination outcome in one combination of readings. The dandiset lacks 11 of the paper's sessions, so the code's unit rule keeps 44 PTupper and 56 PTlower units, not 61 and 69, and the unidentified pools are smaller than the paper's: 232 units in the PTupper cohort against the paper's 495, and 430 in the PTlower cohort against 511. The same code run on the paper's 130 units from its processed data gives the paper's numbers for the unit counts, bursting, late delay selective counts, single-unit decoding, early selectivity and spike rates, and later CDgo onsets than the paper's (30 to 40 and 55 to 80 ms against 24 and 48 ms). The same code on the processed data for the 100 units the dandiset holds (section "Missing units against the file format") separates the two causes: the failures of counts and tests (unit counts, the spike rate test, decoding counts, late delay counts, PTupper persistence, the PTupper CDgo onset) come from the missing units; the files raise every unit's baseline rate by a median of about 1 Hz, and the rectified PTlower decoding mean is out of tolerance on the files (56.20%) and within it on the same units processed (56.84%), which the higher baseline, the rectification threshold, may explain but the notebook does not isolate; and the PTlower CDgo onset is 37.5 ms on all three, so its distance from 24 ms comes from the method here, not the data.
What was found about the dandiset. All seven findings of the earlier attempt are still present, as the version has not changed: 11 sessions of the paper's data are missing (17 PTupper and 13 PTlower included units), the data hold no sequencing or anatomy, spike times are not sorted and trials overlap on the file's clock, the trial table has no stop times and does not say that its event times are relative to the trial start, the PT labels are not described, the quality vocabulary, sex, genotype, electrode positions and layer are uninformative, and the first ten trials of each session are shifted in the paper's processed data. New in this attempt: the missing sessions also hold 285 PTupper-cohort and 117 PTlower-cohort unidentified units (the dandiset has 232 and 430, against the paper's 495 and 511); the lick times are stored trial by trial and out of order in 105 of 106 series, so licks cannot be assigned to trials directly; the files hold no session manipulator depth; and there is no video.
print(f"Notebook ran in {time.time() - started:.0f} s")
Notebook ran in 580 s