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. The paper is not open access and has no copy in PubMed Central, so the claims and criteria were written from its abstract only. Its public Data availability statement was read to find the linked data and code. After the first round of verification, the public legends of Extended Data Figs 5 to 10 were read for the number of units the paper used (61 PTupper, 69 PTlower) and for the single-unit decoding measure of Extended Data Fig. 9; both are used in variants only, and the criteria were not changed.
Dandiset. DANDI:000006, version 0.220126.1855. The paper cites no DANDI version (it points to CRCNS and Figshare). This is the only published version, and it is pinned for that reason. The draft holds the same 53 files with the same sizes; it differs only in metadata added later (ORCID, affiliations, funders, anatomy terms and a link to the DataLad source), none of which the analysis uses.
What the dandiset holds. Spike-sorted units from silicon probe recordings in the anterior
lateral motor cortex (ALM) of 12 mice in a delayed response task (lick left or right after a go
cue, depending on the position of a pole presented in the sample epoch), with trial tables and
lick times. Tagged pyramidal tract (PT) units carry the label PTupper (8 mice) or PTlower
(4 mice). The single-cell and bulk RNA-seq, the anatomy and the axonal reconstructions of the
paper are not in it.
Code read. The paper's Data availability statement points to the Figshare record
doi:10.25378/janelia.7007846 (version 1). Its
Ephys.zip holds the MATLAB analysis code (defineDataset.m, loadAndProcessData.m, simplePlots.m,
Projections.m, CDCorrelation.m, decodingAccuracy.m, VarianceExplained.m and helpers) and the
processed data it reads. The code names the PTlower animals "MP" (medulla-projecting) and the
PTupper animals "TP" (thalamus-projecting), and fixes the epochs, the trial and unit selection and
the selectivity test used here (see alm.py). The processed data are read in the data check, to
show how the files' spike times were built and which of the paper's units the dandiset holds, and
in variants that rerun the tests on the paper's 61 and 69 units.
| Claim | Tested |
|---|---|
| pt-types-transcriptomic | no |
| pt-types-shared-targets | no |
| pt-upper-thalamus-projection | no |
| pt-lower-medulla-projection | no |
| pt-preparatory-activity | yes |
| pt-upper-early-preparatory | yes |
| pt-upper-persistent | yes |
| pt-lower-late-and-motor | yes |
| pt-types-differ-in-timing | yes |
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", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
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 58 s
Data check¶
Each check of the metadata against the values is stated under its output.
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", [])])
rows = []
for s in sessions:
u = s["units"]
rows.append({"subject": s["subject"], "session": s["session"], "hemisphere": s["hemisphere"],
"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())
only_tagged = files[files.unidentified == 0]
print("\nFiles without unidentified units:", len(only_tagged), "; units in them", only_tagged.units.min(), "to", only_tagged.units.max())
print("Electrode x, y, z all zero in every file:", all(s["meta"]["electrodes_xyz_zero"] for s in sessions))
print("Subjects with both PTupper and PTlower units:",
int(((by_subject.PTupper > 0) & (by_subject.PTlower > 0)).sum()))
units = pd.concat([s["units"].assign(subject=s["subject"], session=s["session"]) for s in sessions])
print("\ncell_type counts:", units.cell_type.value_counts().to_dict())
print("Tagged label exactly when id > 1000:", bool(((units.id > 1000) == (units.cell_type != "unidentified")).all()))
print("quality values (all units):", units.quality.value_counts().to_dict())
print("quality values (tagged units):", units[units.cell_type != "unidentified"].quality.value_counts().to_dict())
print("\nDepth (um) by cell_type:")
print(units.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')]
Files 53 ; subjects 12
subject session hemisphere trials units PTupper PTlower unidentified nwb sex genotype session_start
anm369962 20170309 left 274 3 0 3 0 2.0.2 U '' 2017-03-09T00:00:00
anm369962 20170310 right 275 79 0 5 74 2.0.2 U '' 2017-03-10T00:00:00
anm369962 20170313 right 354 71 0 5 66 2.0.2 U '' 2017-03-13T00:00:00
anm369962 20170314 right 306 51 0 4 47 2.0.2 U '' 2017-03-14T00:00:00
anm369962 20170316 right 252 2 0 2 0 2.0.2 U '' 2017-03-16T00:00:00
anm369962 20170317 right 274 59 0 2 57 2.0.2 U '' 2017-03-17T00:00:00
anm369963 20170226 left 295 49 0 2 47 2.0.2 U '' 2017-02-26T00:00:00
anm369963 20170227 left 358 42 0 2 40 2.0.2 U '' 2017-02-27T00:00:00
anm369963 20170228 left 178 1 0 1 0 2.0.2 U '' 2017-02-28T00:00:00
anm369963 20170301 left 233 3 0 3 0 2.0.2 U '' 2017-03-01T00:00:00
anm369963 20170302 left 193 2 0 2 0 2.0.2 U '' 2017-03-02T00:00:00
anm369963 20170306 right 257 61 0 1 60 2.0.2 U '' 2017-03-06T00:00:00
anm369963 20170309 right 245 2 0 2 0 2.0.2 U '' 2017-03-09T00:00:00
anm369963 20170310 right 273 4 0 4 0 2.0.2 U '' 2017-03-10T00:00:00
anm369964 20170320 left 319 55 0 5 50 2.0.2 U '' 2017-03-20T00:00:00
anm369964 20170321 left 386 3 0 3 0 2.0.2 U '' 2017-03-21T00:00:00
anm369964 20170322 left 288 40 0 2 38 2.0.2 U '' 2017-03-22T00:00:00
anm369964 20170323 left 459 2 0 2 0 2.0.2 U '' 2017-03-23T00:00:00
anm372793 20170504 left 304 1 1 0 0 2.0.2 U '' 2017-05-04T00:00:00
anm372793 20170508 left 340 1 1 0 0 2.0.2 U '' 2017-05-08T00:00:00
anm372793 20170512 left 282 2 2 0 0 2.0.2 U '' 2017-05-12T00:00:00
anm372793 20170513 left 368 47 1 0 46 2.0.2 U '' 2017-05-13T00:00:00
anm372793 20170514 left 270 67 3 0 64 2.0.2 U '' 2017-05-14T00:00:00
anm372794 20170621 left 248 4 0 4 0 2.0.2 U '' 2017-06-21T00:00:00
anm372794 20170622 left 340 4 0 4 0 2.0.2 U '' 2017-06-22T00:00:00
anm372794 20170624 left 248 50 0 3 47 2.0.2 U '' 2017-06-24T00:00:00
anm372794 20170625 left 323 56 0 5 51 2.0.2 U '' 2017-06-25T00:00:00
anm372794 20170626 left 328 34 0 1 33 2.0.2 U '' 2017-06-26T00:00:00
anm372794 20170627 left 385 4 0 4 0 2.0.2 U '' 2017-06-27T00:00:00
anm372795 20170714 left 219 4 4 0 0 2.0.2 U '' 2017-07-14T00:00:00
anm372795 20170715 left 203 1 1 0 0 2.0.2 U '' 2017-07-15T00:00:00
anm372795 20170716 left 273 4 4 0 0 2.0.2 U '' 2017-07-16T00:00:00
anm372795 20170718 left 263 2 2 0 0 2.0.2 U '' 2017-07-18T00:00:00
anm372797 20170615 left 224 76 3 0 73 2.0.2 U '' 2017-06-15T00:00:00
anm372797 20170617 left 234 1 1 0 0 2.0.2 U '' 2017-06-17T00:00:00
anm372904 20170615 left 345 1 1 0 0 2.0.2 U '' 2017-06-15T00:00:00
anm372904 20170616 left 314 1 1 0 0 2.0.2 U '' 2017-06-16T00:00:00
anm372904 20170617 left 304 1 1 0 0 2.0.2 U '' 2017-06-17T00:00:00
anm372904 20170618 left 334 70 1 0 69 2.0.2 U '' 2017-06-18T00:00:00
anm372904 20170619 left 254 66 1 0 65 2.0.2 U '' 2017-06-19T00:00:00
anm372905 20170715 left 224 1 1 0 0 2.0.2 U '' 2017-07-15T00:00:00
anm372905 20170716 left 277 2 2 0 0 2.0.2 U '' 2017-07-16T00:00:00
anm372905 20170717 left 291 1 1 0 0 2.0.2 U '' 2017-07-17T00:00:00
anm372906 20170608 left 306 57 4 0 53 2.0.2 U '' 2017-06-08T00:00:00
anm372906 20170610 left 315 1 1 0 0 2.0.2 U '' 2017-06-10T00:00:00
anm372906 20170611 left 227 1 1 0 0 2.0.2 U '' 2017-06-11T00:00:00
anm372906 20170612 left 301 1 1 0 0 2.0.2 U '' 2017-06-12T00:00:00
anm372907 20170608 left 319 2 2 0 0 2.0.2 U '' 2017-06-08T00:00:00
anm372907 20170610 left 310 1 1 0 0 2.0.2 U '' 2017-06-10T00:00:00
anm372907 20170612 left 312 1 1 0 0 2.0.2 U '' 2017-06-12T00:00:00
anm372907 20170613 left 139 1 1 0 0 2.0.2 U '' 2017-06-13T00:00:00
anm372909 20170520 left 241 1 1 0 0 2.0.2 U '' 2017-05-20T00:00:00
anm372909 20170522 left 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 ; units in them 1 to 4
Electrode x, y, z all zero in every file: True
Subjects with both PTupper and PTlower units: 0
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 only one tagged type:
PTupper in 8 mice and PTlower in 4, so the type comparison is also a comparison between animals.
The code's grouping of animals into "MP" and "TP" (defineDataset.m) puts the four PTlower mice in
"MP" (medulla-projecting) and the eight PTupper mice in "TP" (thalamus-projecting). That is how
the types are linked to projections here; the files themselves do not name the projection
target. Tagged units are exactly the units with ids above 1000. Unit quality has a mix of
vocabularies (single, Excellent, Good, Fair, Poor) and some empty values ([], 0,
blank), which are not described. Sex is U and genotype empty in every file, and the session
start times are at midnight. 34 files hold only the tagged units of the session (1 to 4 units),
with no unidentified units, and electrode positions are all zero.
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())
odd = (trials.cue - trials.pole_out - 1.5).abs() > 0.1
print("\nTrials 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 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 dandiset's layout implies but does not
state. The pole comes in 2.5 s before the go cue and starts to leave 1.5 s before it in nearly all
trials, so the code's epochs (sample from -2.5 s, delay from -1.3 s) fit the events. The trials
where the delay is not 1.5 s are nearly all early-lick trials, which are left out. stop_time is
NaN throughout, so trials have no stated end. Consecutive trial starts are a median 6.4 s apart
(5.4 to 7.9 s from the 1st to the 99th percentile), 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 often falls after the next trial's start. Photostimulation was flagged on 2
trials only, and is_good is false for 651 trials; the code does not use is_good (a variant
below drops those trials).
# 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()):
fl = left[np.searchsorted(left, c, "right")] if np.searchsorted(left, c, "right") < len(left) else np.inf
fr = right[np.searchsorted(right, c, "right")] if np.searchsorted(right, c, "right") < 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("Licks are stated in", {s["meta"]["lick_units"]["left"] for s in sessions}, "(timestamps)")
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("\nLeave-one-out decoding accuracy of direction (correct trials):")
print(decode.round(3).to_string(index=False))
print(f"\nMedian accuracy: response {decode.response.median():.3f}; response with shifted labels "
f"{decode['response, labels shifted'].median():.3f}; 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))
Licks are stated in {'Seconds'} (timestamps)
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 accuracy of direction (correct trials): subject session units trials response response, labels shifted baseline p (response > 0.5) anm369962 20170309 3 191 0.880 0.602 0.576 0.000 anm369962 20170310 79 186 0.935 0.568 0.559 0.000 anm369962 20170313 71 262 0.973 0.588 0.473 0.000 anm369962 20170314 51 239 0.979 0.588 0.556 0.000 anm369962 20170316 2 200 0.700 0.538 0.510 0.000 anm369962 20170317 59 226 0.965 0.533 0.562 0.000 anm369963 20170226 49 236 0.987 0.568 0.530 0.000 anm369963 20170227 42 286 0.979 0.570 0.510 0.000 anm369963 20170228 1 127 0.472 0.449 0.559 0.761 anm369963 20170301 3 162 0.988 0.617 0.519 0.000 anm369963 20170302 2 125 0.664 0.480 0.464 0.000 anm369963 20170306 61 202 0.985 0.552 0.500 0.000 anm369963 20170309 2 160 0.938 0.569 0.538 0.000 anm369963 20170310 4 219 0.945 0.557 0.484 0.000 anm369964 20170320 55 261 0.992 0.608 0.521 0.000 anm369964 20170321 3 334 0.740 0.538 0.518 0.000 anm369964 20170322 40 249 0.996 0.593 0.482 0.000 anm369964 20170323 2 421 0.506 0.489 0.549 0.423 anm372793 20170504 1 270 0.596 0.550 0.537 0.001 anm372793 20170508 1 284 0.482 0.518 0.567 0.743 anm372793 20170512 2 224 0.844 0.612 0.509 0.000 anm372793 20170513 47 300 0.987 0.532 0.557 0.000 anm372793 20170514 67 213 0.981 0.552 0.573 0.000 anm372794 20170621 4 158 0.854 0.525 0.544 0.000 anm372794 20170622 4 267 0.764 0.539 0.569 0.000 anm372794 20170624 50 208 0.981 0.585 0.529 0.000 anm372794 20170625 56 280 0.982 0.518 0.579 0.000 anm372794 20170626 34 257 0.981 0.580 0.553 0.000 anm372794 20170627 4 313 0.671 0.508 0.537 0.000 anm372795 20170714 4 193 0.953 0.508 0.451 0.000 anm372795 20170715 1 164 0.665 0.524 0.476 0.000 anm372795 20170716 4 235 0.966 0.570 0.583 0.000 anm372795 20170718 2 157 0.834 0.471 0.541 0.000 anm372797 20170615 76 150 0.980 0.553 0.587 0.000 anm372797 20170617 1 115 0.496 0.079 0.426 0.574 anm372904 20170615 1 308 0.893 0.573 0.516 0.000 anm372904 20170616 1 284 0.187 0.572 0.525 1.000 anm372904 20170617 1 250 0.764 0.548 0.496 0.000 anm372904 20170618 70 286 0.986 0.582 0.542 0.000 anm372904 20170619 66 203 0.970 0.547 0.532 0.000 anm372905 20170715 1 161 0.596 0.522 0.540 0.009 anm372905 20170716 2 211 0.919 0.600 0.431 0.000 anm372905 20170717 1 230 0.848 0.572 0.491 0.000 anm372906 20170608 57 223 0.991 0.556 0.556 0.000 anm372906 20170610 1 225 0.693 0.524 0.538 0.000 anm372906 20170611 1 159 0.604 0.528 0.497 0.005 anm372906 20170612 1 257 0.556 0.494 0.502 0.040 anm372907 20170608 2 228 0.917 0.581 0.548 0.000 anm372907 20170610 1 218 0.844 0.587 0.509 0.000 anm372907 20170612 1 236 0.746 0.500 0.492 0.000 anm372907 20170613 1 88 0.489 0.563 0.500 0.625 anm372909 20170520 1 167 0.605 0.509 0.479 0.004 anm372909 20170522 53 213 0.991 0.465 0.521 0.000 Median accuracy: 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 (11689 of 12191) and on the other side in nearly all incorrect trials (1832 of 1885), with a median latency of 0.12 s; with the labels of the next trial the match falls to 5381 of 12191 correct trials. Spikes share it too: direction is decoded from the response epoch with a median accuracy of 0.89, against 0.53 in the baseline epoch and 0.55 with the labels of the next trial. The seven sessions listed above, each with one or two units, could not be checked this way; for them the clock is assumed. Two single-unit sessions give accuracies far below 0.5 (0.19 with true labels in anm372904 2017-06-16, 0.08 with shifted labels in anm372797 2017-06-17), which a leave-one-out decoder can produce from a feature without information; they are not read as evidence either way.
# How the spike times are stored.
rows = []
for s in sessions:
t = s["trials"]
cue = t.cue_abs.to_numpy()
for _, u in s["units"].iterrows():
raw = np.asarray(u.spikes)
flags = alm.overlap_flags(raw)
sorted_spikes = 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)), "flagged": int(flags.sum())}
for name in ("baseline", "sample", "delay", "response"):
a, b = alm.EPOCHS[name]
inside = np.zeros(len(raw), bool)
lo = np.searchsorted(sorted_spikes, cue + a)
hi = np.searchsorted(sorted_spikes, cue + b)
order = np.argsort(raw, kind="stable")
for l, h in zip(lo, hi):
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("Backward steps per unit, median", storage["backward steps"].median(), "; max", storage["backward steps"].max())
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%})")
# The paper's processed data for one session (anm372906, 2017-06-08), read from Figshare.
mat = alm.read_figshare_processed()
s = next(x for x in sessions if x["subject"] == "anm372906" and x["session"] == "20170608")
per_unit, spans = alm.compare_with_processed(s, mat, 0)
print(f"\nProcessed file {alm.FIGSHARE_CHECK_FILE.split('/')[-1]}, session 1: {len(per_unit)} units; "
f"units whose stored spike times equal trial start + go cue + processed time: "
f"{int(per_unit.equal_in_stored_order.sum())}")
spans["length"] = spans["last"] - spans["first"]
spans["first - pole_in"] = spans["first"] - spans["pole_in"]
print("Span of each trial's processed spike times, seconds from the go cue (percentiles 0, 5, 50, 95, 100):")
print(spans[["first", "last", "length", "first - pole_in", "start_spacing"]].quantile([0, 0.05, 0.5, 0.95, 1]).round(3).T.to_string())
print("First 12 trials:")
print(spans.head(12).round(3).to_string())
Units with spike times out of order: 1120 of 1149 Backward steps per unit, median 38.0 ; max 156 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%)
Processed file data_structure_anm372906_WR40_psth.mat, session 1: 57 units; units whose stored spike times equal trial start + go cue + processed time: 57
Span of each trial's processed spike times, seconds from the go cue (percentiles 0, 5, 50, 95, 100):
0.00 0.05 0.50 0.95 1.00
first -6.710 -3.473 -3.471 -0.772 -0.761
last -0.322 2.912 2.925 5.622 5.628
length 6.360 6.385 6.395 6.399 6.400
first - pole_in -0.973 -0.972 -0.970 1.729 8.566
start_spacing 5.918 6.399 6.400 6.400 6.882
First 12 trials:
first last pole_in start_spacing length first - pole_in
0 -0.766 5.625 -2.978 6.401 6.391 2.211
1 -0.768 5.626 -2.500 6.399 6.394 1.733
2 -0.772 5.626 -2.500 6.401 6.398 1.729
3 -0.761 5.624 -2.501 6.399 6.385 1.739
4 -0.772 5.627 -2.501 6.400 6.399 1.729
5 -0.769 5.626 -2.501 6.400 6.395 1.731
6 -0.770 5.626 -2.501 6.400 6.397 1.730
7 -0.771 5.620 -2.501 6.400 6.391 1.730
8 -0.769 5.622 -2.501 6.400 6.391 1.732
9 -0.764 5.626 -2.501 6.400 6.390 1.736
10 -3.472 2.927 -2.501 6.400 6.399 -0.971
11 -3.715 2.669 -2.743 6.400 6.383 -0.971
The spike times are not sorted in 1120 of 1149 units, which breaks the NWB convention for
spike_times. In the session compared here with the paper's processed data (all sessions are
compared further below), every unit's stored spike times are exactly trial start + go cue time + the processed spike time relative to the go
cue, in the processed order (trial by trial). In the processed data each trial's spikes cover a
window of 6.4 s, from about 0.97 s before the pole comes in; the first ten trials of that session
are shifted by about 2.7 s, and the code leaves out the first ten trials of every unit. As trial
starts in the files are 6.4 s apart but the pole comes in at different times, consecutive trials'
windows overlap on the file's clock, and spikes of one trial land in the next trial's early
epochs. The stored order shows part of this: the flagged spikes above come before, in storage, a
spike of a later trial at an earlier time. Among the spikes counted in each epoch on the file's
clock, 11.7% are flagged in the baseline epoch and under 2% in the sample, delay and response
epochs. Overlaps that leave
no backward step are not detected, so the flagged share is a lower bound.
The analysis below uses the file's clock as published (spikes counted in each epoch around each trial's go cue). A variant leaves out the flagged spikes.
records = alm.tagged_units(sessions)
units_df = pd.DataFrame([{k: v for k, v in r.items() if k not in ("rates", "is_right", "cue", "spike_times")} for r in records])
print("Tagged units:", len(units_df), units_df.cell_type.value_counts().to_dict())
print("\nInclusion (the code's rule):")
print(pd.crosstab(units_df.cell_type, units_df.included).to_string())
print("\nReasons for leaving units out (a unit can have several):")
for t in alm.TYPES:
reasons = units_df[(units_df.cell_type == t) & ~units_df.included].reason.str.split(", ").explode()
print(" ", t, reasons.value_counts().to_dict())
inc = units_df[units_df.included]
print("\nIncluded units by animal:")
print(inc.groupby(["cell_type", "subject"]).size().to_string())
print("\nIncluded units by session:")
print(inc.groupby(["cell_type", "subject", "session"]).size().unstack(0).fillna(0).astype(int).to_string())
print("\nCorrect trials per included unit: left", int(inc.n_correct_L.min()), "to", int(inc.n_correct_L.max()),
"; right", int(inc.n_correct_R.min()), "to", int(inc.n_correct_R.max()))
print("Spikes per trial of included units:", round(inc.spikes_per_trial.min(), 2), "to", round(inc.spikes_per_trial.max(), 2))
print("Depth of included units (um):", inc.groupby("cell_type").depth.describe()[["count", "min", "50%", "max"]].to_string())
Tagged units: 117 {'PTlower': 71, 'PTupper': 46}
Inclusion (the code's rule):
included False True
cell_type
PTlower 15 56
PTupper 2 44
Reasons for leaving units out (a unit can have several):
PTupper {'performance': 1, 'too few correct trials': 1}
PTlower {'low rate': 7, 'quality': 6, 'performance': 2, 'too few correct trials': 1}
Included units by animal:
cell_type subject
PTlower anm369962 20
anm369963 13
anm369964 9
anm372794 14
PTupper anm372793 8
anm372795 11
anm372797 3
anm372904 5
anm372905 4
anm372906 7
anm372907 4
anm372909 2
Included units by session:
cell_type PTlower PTupper
subject session
anm369962 20170309 3 0
20170310 4 0
20170313 5 0
20170314 4 0
20170316 2 0
20170317 2 0
anm369963 20170226 2 0
20170227 2 0
20170228 1 0
20170301 3 0
20170309 2 0
20170310 3 0
anm369964 20170320 4 0
20170321 2 0
20170322 2 0
20170323 1 0
anm372793 20170504 0 1
20170508 0 1
20170512 0 2
20170513 0 1
20170514 0 3
anm372794 20170621 4 0
20170622 1 0
20170624 3 0
20170625 3 0
20170627 3 0
anm372795 20170714 0 4
20170715 0 1
20170716 0 4
20170718 0 2
anm372797 20170615 0 3
anm372904 20170615 0 1
20170616 0 1
20170617 0 1
20170618 0 1
20170619 0 1
anm372905 20170715 0 1
20170716 0 2
20170717 0 1
anm372906 20170608 0 4
20170610 0 1
20170611 0 1
20170612 0 1
anm372907 20170608 0 2
20170610 0 1
20170612 0 1
anm372909 20170520 0 1
20170522 0 1
Correct trials per included unit: left 62 to 203 ; right 56 to 199
Spikes per trial of included units: 3.98 to 126.71
Depth of included units (um): count min 50% max
cell_type
PTlower 56.0 630.0 977.5 1600.0
PTupper 44.0 430.0 617.5 930.0
Each session's tagged units are counted as separate units, as in the code. Units recorded in the same animal on different days may be the same neuron; the files have no identifier that would show it, so results are also given by animal. The depth is the nominal depth of the unit along the probe; the files give no layer per unit (every electrode group says "cortical_layer: 5").
# The paper's processed data (Ephys.zip on Figshare), every session: how the files' spike times
# were built, the code's inclusion rule on the paper's own units, and which of them the dandiset holds.
started = time.time()
print("Processed data files (_psth.mat) in Ephys.zip, all read below:", len(alm.processed_members()))
proc_records, storage_all, shifts = alm.read_processed(sessions)
proc = pd.DataFrame([{k: v for k, v in r.items() if k not in ("rates", "is_right")} for r in proc_records])
print(f"Processed data read in {time.time() - started:.0f} s: {len(proc)} tagged units in "
f"{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; units of the processed data 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 (median per session):")
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))
print("\nThe code's inclusion rule on the processed data (the Extended Data legends give n = 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())
both = proc.merge(units_df, on=["subject", "session", "id"], suffixes=(" processed", " dandiset"))
print(f"\nTagged units in both: {len(both)} of the dandiset's {len(units_df)}; same cell type "
f"{int((both['cell_type processed'] == both['cell_type dandiset']).sum())}; same inclusion decision "
f"{int((both['included processed'] == both['included dandiset']).sum())}; same reasons "
f"{int((both['reason processed'] == both['reason dandiset']).sum())}")
print("Reasons for leaving out, dandiset units:")
for t in alm.TYPES:
reasons = units_df[(units_df.cell_type == t) & ~units_df.included].reason.str.split(", ").explode()
print(" ", t, reasons.value_counts().to_dict())
missing = proc[~proc.in_dandiset]
print("\nProcessed sessions with tagged units that are not in the dandiset:")
print(missing.groupby(["cell_type", "subject", "session"]).agg(tagged=("id", "size"), included=("included", "sum")).to_string())
print("Included units in them:", missing[missing.included].cell_type.value_counts().to_dict())
print("\nDandiset description:", meta["description"])
Processed data files (_psth.mat) in Ephys.zip, all read below: 14
Processed data read in 55 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; units of the processed data not in the file: 0
Start of each trial's processed spike window, seconds from the pole coming in (median per session):
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
The code's inclusion rule on the processed data (the Extended Data legends give n = 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
Tagged units in both: 117 of the dandiset's 117; same cell type 117; same inclusion decision 117; same reasons 117
Reasons for leaving out, dandiset units:
PTupper {'performance': 1, 'too few correct trials': 1}
PTlower {'low rate': 7, 'quality': 6, 'performance': 2, 'too few correct trials': 1}
Processed sessions with tagged units that are not in the dandiset:
tagged included
cell_type subject session
PTlower anm369962 20170315 6 6
anm369963 20170308 4 3
anm372794 20170623 4 4
PTupper anm372793 20170505 3 3
20170511 2 2
anm372795 20170713 1 1
20170717 2 2
anm372905 20170713 3 3
anm372907 20170609 3 3
20170611 3 3
anm372909 20170519 2 0
Included units in them: {'PTupper': 17, 'PTlower': 13}
Dandiset 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.
The paper's processed data hold 150 tagged units in 64 sessions. In all 53 sessions of the dandiset, every unit's stored spike times (1149 of 1149 units) equal trial start + go cue time + the processed spike time, in the processed order. The files' spike times were therefore built the same way in every session, including the seven sessions whose clock the decoding check could not establish. In all 53 sessions the spike windows of the first ten trials start later, relative to the pole coming in, than those of the other trials (session medians 1.73 to 5.84 s after the pole comes in, against -0.97 to -0.16 s); the code leaves out the first ten trials of each unit.
The code's inclusion rule applied to the processed data keeps 61 PTupper and 69 PTlower units, the numbers in the Extended Data legends, so the rule used here gives the paper's selection. All 117 tagged units of the dandiset are in the processed data, with the same cell type, the same inclusion decision and the same reasons, so the 15 PTlower units left out here were left out by the paper's code too. The difference in counts comes from the sessions deposited: 11 processed sessions with tagged units are not in the dandiset, and they hold 17 included PTupper and 13 included PTlower units (44 + 17 = 61, 56 + 13 = 69). The dandiset's description says that some files were found to be broken and were left out. Variants below rerun each test on the paper's 130 units, from the processed data.
# Per unit: selectivity p-value and d' in each epoch (correct trials, left against right).
def unit_table(recs):
rows = []
for r in recs:
if not r["included"]:
continue
row = {"subject": r["subject"], "session": r["session"], "id": r["id"], "cell_type": r["cell_type"],
"depth": r["depth"]}
for name, rates in r["rates"].items():
row[f"p {name}"] = alm.ranksum_p(rates, r["is_right"])
row[f"d {name}"] = alm.dprime(rates, r["is_right"])
rows.append(row)
return pd.DataFrame(rows)
def shuffle_fraction(recs, cell_types, epoch, rng, n=1000, threshold=alm.SELECTIVITY_P):
# Fraction of included units of `cell_types` selective in `epoch`, and its null from
# shuffling each unit's trial labels.
obs, null = [], []
for r in recs:
if not r["included"] or r["cell_type"] not in cell_types:
continue
rates, right = r["rates"][epoch], r["is_right"]
obs.append(alm.ranksum_p(rates, right) < threshold)
null.append(alm.shuffled_ranksum_p(rates, right, alm.label_shuffles(len(right), n, rng)) < threshold)
obs, null = np.array(obs), np.array(null).mean(axis=0)
return obs.mean(), obs.sum(), len(obs), null, alm.null_p(obs.mean(), null)
main = unit_table(records)
records_overlap = alm.tagged_units(sessions, drop_overlap=True)
records_good = alm.tagged_units(sessions, exclude_bad=True)
variant_overlap = unit_table(records_overlap)
variant_good = unit_table(records_good)
# Check the vectorised shuffle test against scipy on one unit.
r0 = next(r for r in records if r["included"])
perm = alm.label_shuffles(len(r0["is_right"]), 3, np.random.default_rng(0))
fast = alm.shuffled_ranksum_p(r0["rates"]["delay"], r0["is_right"], perm)
slow = [alm.ranksum_p(r0["rates"]["delay"], r0["is_right"][p]) for p in perm]
print("Shuffle test against scipy on one unit:", np.round(fast, 6), np.round(slow, 6))
fast_d = alm.shuffled_dprime(r0["rates"]["delay"], r0["is_right"], perm)
slow_d = [alm.dprime(r0["rates"]["delay"], r0["is_right"][p]) for p in perm]
print("Shuffled d' against direct:", np.round(fast_d, 6), np.round(slow_d, 6))
print("\nIncluded units: main", len(main), "; without flagged spikes", len(variant_overlap),
"; without is_good false trials", len(variant_good))
print("\nFraction of included units selective (p < 0.01) per epoch:")
print(main.groupby("cell_type")[[f"p {e}" for e in alm.EPOCHS]].agg(lambda p: (p < 0.01).mean()).round(3).to_string())
print("\nSame, without flagged spikes:")
print(variant_overlap.groupby("cell_type")[[f"p {e}" for e in alm.EPOCHS]].agg(lambda p: (p < 0.01).mean()).round(3).to_string())
print("\nSame, without is_good false trials:")
print(variant_good.groupby("cell_type")[[f"p {e}" for e in alm.EPOCHS]].agg(lambda p: (p < 0.01).mean()).round(3).to_string())
Shuffle test against scipy on one unit: [0.373262 0.647188 0.79109 ] [0.373262 0.647188 0.79109 ]
Shuffled d' against direct: [ 0.08357 -0.070126 -0.001129] [ 0.08357 -0.070126 -0.001129]
Included units: main 100 ; without flagged spikes 100 ; without is_good false trials 100
Fraction of included units selective (p < 0.01) per epoch:
p baseline p sample p delay p late delay p response
cell_type
PTlower 0.018 0.482 0.768 0.714 0.643
PTupper 0.045 0.614 0.659 0.614 0.841
Same, without flagged spikes:
p baseline p sample p delay p late delay p response
cell_type
PTlower 0.036 0.482 0.786 0.714 0.643
PTupper 0.045 0.591 0.659 0.614 0.841
Same, without is_good false trials:
p baseline p sample p delay p late delay p response
cell_type
PTlower 0.018 0.482 0.768 0.714 0.643
PTupper 0.045 0.614 0.659 0.614 0.841
# Units in the sessions whose clock the decoding check could not establish (the data check).
unchecked = main.merge(not_shown[["subject", "session"]], on=["subject", "session"])
print("Included units in those sessions:", len(unchecked))
print(unchecked[["subject", "session", "id", "cell_type"] + [f"p {e}" for e in alm.EPOCHS]].round(4).to_string(index=False))
# Units recorded in one animal on different days may be the same neuron. Depths of the included
# units by session, and a subset with one session per animal (the session with the most included
# tagged units, the earliest if tied), used as a variant below.
print("\nDepths (um) of included units by animal and session:")
print(main.groupby(["cell_type", "subject", "session"]).depth.agg(lambda d: ", ".join(f"{x:.0f}" for x in sorted(d))).to_string())
near = [bool(((main.subject == u.subject) & (main.session != u.session) & ((main.depth - u.depth).abs() <= 20)).any())
for u in main.itertuples()]
print("Included units with a unit of the same animal within 20 um on another day:", sum(near), "of", len(main))
counts = main.groupby(["subject", "session"]).size().rename("n").reset_index()
chosen = counts.sort_values(["subject", "n", "session"], ascending=[True, False, True]).groupby("subject").head(1)
chosen_keys = set(zip(chosen.subject, chosen.session))
one = main[[k in chosen_keys for k in zip(main.subject, main.session)]]
records_one = [r for r in records if (r["subject"], r["session"]) in chosen_keys]
print("One session per animal:", one.groupby("cell_type").size().to_dict(), "units in", len(chosen_keys), "sessions")
# The same per-unit table on the paper's units from the processed data.
proc_main = unit_table(proc_records)
print("Paper's units (processed data), included:", proc_main.groupby("cell_type").size().to_dict())
print("\nFraction selective (p < 0.01) per epoch, paper's units:")
print(proc_main.groupby("cell_type")[[f"p {e}" for e in alm.EPOCHS]].agg(lambda p: (p < 0.01).mean()).round(3).to_string())
print("\nSame, one session per animal:")
print(one.groupby("cell_type")[[f"p {e}" for e in alm.EPOCHS]].agg(lambda p: (p < 0.01).mean()).round(3).to_string())
Included units in those sessions: 5
subject session id cell_type p baseline p sample p delay p late delay p response
anm369963 20170228 1001 PTlower 0.0844 0.0007 0.0000 0.0000 0.6282
anm369964 20170323 1002 PTlower 0.6529 0.0032 0.0692 0.0349 1.0000
anm372793 20170508 1001 PTupper 0.0366 0.0000 0.0000 0.0000 0.7023
anm372904 20170616 1001 PTupper 0.4581 0.0000 0.1946 0.4307 0.9734
anm372906 20170612 1001 PTupper 0.7065 0.0011 0.0000 0.0000 0.0045
Depths (um) of included units by animal and session:
cell_type subject session
PTlower anm369962 20170309 730, 980, 1005
20170310 840, 865, 890, 940
20170313 1000, 1150, 1250, 1425, 1600
20170314 1105, 1130, 1230, 1255
20170316 980, 1105
20170317 1235, 1460
anm369963 20170226 705, 905
20170227 930, 1005
20170228 780
20170301 925, 975, 975
20170309 995, 1520
20170310 1205, 1255, 1480
anm369964 20170320 880, 930, 1030, 1030
20170321 805, 855
20170322 930, 1030
20170323 1135
anm372794 20170621 630, 655, 830, 905
20170622 780
20170624 780, 880, 955
20170625 765, 765, 1015
20170627 830, 980, 1005
PTupper anm372793 20170504 430
20170508 755
20170512 630, 705
20170513 830
20170514 640, 740, 740
anm372795 20170714 480, 505, 555, 580
20170715 760
20170716 730, 780, 855, 905
20170718 800, 850
anm372797 20170615 430, 455, 480
anm372904 20170615 480
20170616 535
20170617 580
20170618 580
20170619 740
anm372905 20170715 580
20170716 655, 680
20170717 530
anm372906 20170608 480, 505, 555, 555
20170610 605
20170611 580
20170612 630
anm372907 20170608 580, 630
20170610 655
20170612 930
anm372909 20170520 755
20170522 530
Included units with a unit of the same animal within 20 um on another day: 44 of 100
One session per animal: {'PTlower': 16, 'PTupper': 20} units in 12 sessions
Paper's units (processed data), included: {'PTlower': 69, 'PTupper': 61}
Fraction selective (p < 0.01) per epoch, paper's units:
p baseline p sample p delay p late delay p response
cell_type
PTlower 0.029 0.507 0.710 0.652 0.551
PTupper 0.033 0.656 0.656 0.623 0.770
Same, one session per animal:
p baseline p sample p delay p late delay p response
cell_type
PTlower 0.00 0.25 0.688 0.625 0.438
PTupper 0.05 0.40 0.750 0.600 0.800
The five included units in the sessions whose clock the decoding check could not establish are listed with their p-values. Their spike times equal those of the processed data (data check), so they are on the same clock as the other units. Units recorded in one animal on different days may be the same neuron. The files hold no waveforms, and 44 of the 100 included units lie within 20 um of a unit of the same animal recorded on another day, so repeats cannot be ruled out. Each claim below therefore also gives the result on one session per animal (20 PTupper and 16 PTlower units in 12 sessions), and on the paper's 61 PTupper and 69 PTlower units from the processed data. These variants do not change the outcomes.
The baseline epoch, before the pole comes in, is a reference: a selective fraction there near 1% would mean the trial labels carry no information before the sample epoch.
pt-preparatory-activity¶
Descending PT neurons in layer 5 of motor cortex carry preparatory activity, which predicts the upcoming movement before it is initiated.
Criterion: over all included tagged units, the fraction selective in the delay epoch (-1.3 to 0 s) is above chance by the label-shuffle null (p < 0.01).
rng = np.random.default_rng(1)
frac, k, n, null, p = shuffle_fraction(records, alm.TYPES, "delay", rng)
print(f"All tagged units, delay epoch: {k} of {n} selective ({frac:.3f}); shuffle null median {np.median(null):.3f}, "
f"99th percentile {np.quantile(null, 0.99):.3f}; p = {p:.4f}")
for t in alm.TYPES:
f, k, n, null, p = shuffle_fraction(records, (t,), "delay", rng)
print(f" {t}: {k} of {n} ({f:.3f}); null 99th percentile {np.quantile(null, 0.99):.3f}; p = {p:.4f}")
print("\nBy animal (delay epoch, p < 0.01):")
print(main.assign(sel=main["p delay"] < 0.01).groupby(["cell_type", "subject"]).sel.agg(["sum", "size", "mean"]).round(3).to_string())
print("\nThreshold of the selectivity test, delay epoch, all tagged units:")
for thr in (0.05, 0.01, 0.001):
print(f" p < {thr}: {(main['p delay'] < thr).mean():.3f} ({(main['p delay'] < thr).sum()} of {len(main)})")
for name, table in (("without flagged spikes", variant_overlap), ("without is_good false trials", variant_good)):
print(f"Variant {name}: {(table['p delay'] < 0.01).sum()} of {len(table)} ({(table['p delay'] < 0.01).mean():.3f})")
print("\nDepth of delay-selective units (um):", main[main["p delay"] < 0.01].depth.describe()[["min", "50%", "max"]].to_dict())
All tagged units, delay epoch: 72 of 100 selective (0.720); shuffle null median 0.010, 99th percentile 0.040; p = 0.0010
PTupper: 29 of 44 (0.659); null 99th percentile 0.045; p = 0.0010
PTlower: 43 of 56 (0.768); null 99th percentile 0.054; p = 0.0010
By animal (delay epoch, p < 0.01):
sum size mean
cell_type subject
PTlower anm369962 16 20 0.800
anm369963 11 13 0.846
anm369964 6 9 0.667
anm372794 10 14 0.714
PTupper anm372793 5 8 0.625
anm372795 9 11 0.818
anm372797 2 3 0.667
anm372904 4 5 0.800
anm372905 1 4 0.250
anm372906 4 7 0.571
anm372907 3 4 0.750
anm372909 1 2 0.500
Threshold of the selectivity test, delay epoch, all tagged units:
p < 0.05: 0.760 (76 of 100)
p < 0.01: 0.720 (72 of 100)
p < 0.001: 0.670 (67 of 100)
Variant without flagged spikes: 73 of 100 (0.730)
Variant without is_good false trials: 72 of 100 (0.720)
Depth of delay-selective units (um): {'min': 430.0, '50%': 845.0, 'max': 1600.0}
# Variants: one session per animal, and the paper's units from the processed data.
rng = np.random.default_rng(11)
for name, recs, table in (("one session per animal", records_one, one), ("paper's units (processed data)", proc_records, proc_main)):
f, k, n, null, p = shuffle_fraction(recs, alm.TYPES, "delay", rng)
per_type = "; ".join(f"{t} {(table[table.cell_type == t]['p delay'] < 0.01).sum()} of {(table.cell_type == t).sum()}" for t in alm.TYPES)
print(f"Variant {name}: delay epoch {k} of {n} ({f:.3f}); null 99th percentile {np.quantile(null, 0.99):.3f}; p = {p:.4f} ({per_type})")
Variant one session per animal: delay epoch 26 of 36 (0.722); null 99th percentile 0.056; p = 0.0010 (PTupper 15 of 20; PTlower 11 of 16)
Variant paper's units (processed data): delay epoch 89 of 130 (0.685); null 99th percentile 0.031; p = 0.0010 (PTupper 40 of 61; PTlower 49 of 69)
Result: reproduced. 72 of 100 included tagged units are selective in the delay epoch (p < 0.01), against a shuffle null whose 99th percentile is 0.040 (p = 0.001, the smallest 1000 shuffles allow). It holds in both types (PTupper 29 of 44, PTlower 43 of 56) and in every animal (25% to 85% of its units), at test thresholds from 0.05 to 0.001 (76 to 67 units), without flagged spikes (73) and without the trials marked not good (72). With one session per animal 26 of 36 units are selective (p = 0.001), and among the paper's 130 units from the processed data 89 are (p = 0.001). The selective units lie at nominal depths of 430 to 1600 um; the files give no layer per unit, so "layer 5" is not checked.
pt-upper-early-preparatory¶
The thalamus-projecting PT type produces early preparatory activity.
Criterion: among included PTupper units, the fraction selective in the sample epoch (-2.5 to -2.1 s) is above chance by the label-shuffle null (p < 0.01).
rng = np.random.default_rng(2)
frac, k, n, null, p = shuffle_fraction(records, ("PTupper",), "sample", rng)
print(f"PTupper, sample epoch: {k} of {n} selective ({frac:.3f}); shuffle null median {np.median(null):.3f}, "
f"99th percentile {np.quantile(null, 0.99):.3f}; p = {p:.4f}")
f0, k0, n0, null0, p0 = shuffle_fraction(records, ("PTupper",), "baseline", rng)
print(f"PTupper, baseline epoch (reference): {k0} of {n0} ({f0:.3f}); p = {p0:.4f}")
print("\nBy animal (sample epoch, p < 0.01):")
up = main[main.cell_type == "PTupper"]
print(up.assign(sel=up["p sample"] < 0.01).groupby("subject").sel.agg(["sum", "size", "mean"]).round(3).to_string())
for thr in (0.05, 0.01, 0.001):
print(f"Threshold p < {thr}: {(up['p sample'] < thr).sum()} of {len(up)} ({(up['p sample'] < thr).mean():.3f})")
for name, recs in (("without flagged spikes", records_overlap), ("without is_good false trials", records_good)):
f, k, n, null, p = shuffle_fraction(recs, ("PTupper",), "sample", rng)
print(f"Variant {name}: {k} of {n} ({f:.3f}); p = {p:.4f}")
PTupper, sample epoch: 27 of 44 selective (0.614); shuffle null median 0.000, 99th percentile 0.045; p = 0.0010
PTupper, baseline epoch (reference): 2 of 44 (0.045); p = 0.0789
By animal (sample epoch, p < 0.01):
sum size mean
subject
anm372793 4 8 0.500
anm372795 8 11 0.727
anm372797 2 3 0.667
anm372904 4 5 0.800
anm372905 1 4 0.250
anm372906 3 7 0.429
anm372907 4 4 1.000
anm372909 1 2 0.500
Threshold p < 0.05: 35 of 44 (0.795)
Threshold p < 0.01: 27 of 44 (0.614)
Threshold p < 0.001: 19 of 44 (0.432)
Variant without flagged spikes: 26 of 44 (0.591); p = 0.0010
Variant without is_good false trials: 27 of 44 (0.614); p = 0.0010
# Variants: one session per animal, and the paper's units from the processed data.
rng = np.random.default_rng(12)
for name, recs in (("one session per animal", records_one), ("paper's units (processed data)", proc_records)):
f, k, n, null, p = shuffle_fraction(recs, ("PTupper",), "sample", rng)
f0, k0, n0, _, p0 = shuffle_fraction(recs, ("PTupper",), "baseline", rng)
print(f"Variant {name}: PTupper sample epoch {k} of {n} ({f:.3f}); null 99th percentile {np.quantile(null, 0.99):.3f}; "
f"p = {p:.4f}; baseline {k0} of {n0} (p = {p0:.4f})")
Variant one session per animal: PTupper sample epoch 8 of 20 (0.400); null 99th percentile 0.100; p = 0.0010; baseline 1 of 20 (p = 0.1648)
Variant paper's units (processed data): PTupper sample epoch 40 of 61 (0.656); null 99th percentile 0.049; p = 0.0010; baseline 2 of 61 (p = 0.1259)
Result: reproduced. 27 of 44 included PTupper units (0.614) are selective in the sample epoch, against a null 99th percentile of 0.045 (p = 0.001); in the baseline epoch, before the pole comes in, 2 of 44 are (p = 0.079). Every PTupper animal has at least one sample-selective unit (25% to 100% of its units). The fraction is 0.795 at p < 0.05 and 0.432 at p < 0.001, and 26 of 44 units are selective without flagged spikes. With one session per animal 8 of 20 units are selective (null 99th percentile 0.100, p = 0.001), and among the paper's 61 PTupper units 40 are (p = 0.001).
pt-upper-persistent¶
The early preparatory activity of the thalamus-projecting PT type persists until the movement is initiated.
Criterion: (a) the fraction of included PTupper units selective in the late delay epoch (-0.4 to 0 s) is above chance (p < 0.01); (b) the cosine similarity of the sample and late delay d' vectors of PTupper units is above that of 95% of label shuffles (p < 0.05). Both: reproduced; one: partially reproduced.
rng = np.random.default_rng(3)
def similarity_test(recs, cell_type, rng, n=1000):
# Cosine similarity of sample and late delay d' over included units of a type, and its null
# from shuffling each unit's trial labels (the same shuffle for both epochs).
d_s, d_l, null_s, null_l = [], [], [], []
for r in recs:
if not r["included"] or r["cell_type"] != cell_type:
continue
right = r["is_right"]
perms = alm.label_shuffles(len(right), n, rng)
d_s.append(alm.dprime(r["rates"]["sample"], right))
d_l.append(alm.dprime(r["rates"]["late delay"], right))
null_s.append(alm.shuffled_dprime(r["rates"]["sample"], right, perms))
null_l.append(alm.shuffled_dprime(r["rates"]["late delay"], right, perms))
d_s, d_l = np.array(d_s), np.array(d_l)
null_s, null_l = np.array(null_s), np.array(null_l)
null = (null_s * null_l).sum(axis=0) / (np.linalg.norm(null_s, axis=0) * np.linalg.norm(null_l, axis=0))
obs = alm.cosine(d_s, d_l)
return obs, null, alm.null_p(obs, null), d_s, d_l
frac, k, n, null, p = shuffle_fraction(records, ("PTupper",), "late delay", rng)
print(f"(a) PTupper, late delay: {k} of {n} selective ({frac:.3f}); null 99th percentile {np.quantile(null, 0.99):.3f}; p = {p:.4f}")
obs, null, p_sim, d_s, d_l = similarity_test(records, "PTupper", rng)
print(f"(b) PTupper, cosine similarity of sample and late delay d': {obs:.3f}; shuffle null median {np.median(null):.3f}, "
f"95th percentile {np.quantile(null, 0.95):.3f}; p = {p_sim:.4f}")
r, p_r = stats.pearsonr(d_s, d_l)
print(f"Variant, the code's Pearson correlation of sample and late delay d' (Projections.m): r = {r:.3f}, p = {p_r:.4g}, n = {len(d_s)}")
sel_s = (np.array([alm.ranksum_p(x["rates"]["sample"], x["is_right"]) for x in records if x["included"] and x["cell_type"] == "PTupper"]) < 0.01)
same_sign = np.sign(d_s) == np.sign(d_l)
print(f"Units selective in the sample epoch whose late delay d' has the same sign: {int((sel_s & same_sign).sum())} of {int(sel_s.sum())}")
print("\nBy animal (late delay, p < 0.01):")
up = main[main.cell_type == "PTupper"]
print(up.assign(sel=up["p late delay"] < 0.01).groupby("subject").sel.agg(["sum", "size", "mean"]).round(3).to_string())
for name, recs in (("without flagged spikes", records_overlap), ("without is_good false trials", records_good)):
f, k, n, _, p = shuffle_fraction(recs, ("PTupper",), "late delay", rng)
o, nl, ps, _, _ = similarity_test(recs, "PTupper", rng)
print(f"Variant {name}: (a) {k} of {n} ({f:.3f}), p = {p:.4f}; (b) cosine {o:.3f}, p = {ps:.4f}")
print(f"Range of plotted d': sample {d_s.min():.2f} to {d_s.max():.2f}; late delay {d_l.min():.2f} to {d_l.max():.2f}")
(a) PTupper, late delay: 27 of 44 selective (0.614); null 99th percentile 0.046; p = 0.0010
(b) PTupper, cosine similarity of sample and late delay d': 0.495; shuffle null median 0.182, 95th percentile 0.433; p = 0.0160
Variant, the code's Pearson correlation of sample and late delay d' (Projections.m): r = 0.496, p = 0.0006251, n = 44
Units selective in the sample epoch whose late delay d' has the same sign: 17 of 27
By animal (late delay, p < 0.01):
sum size mean
subject
anm372793 5 8 0.625
anm372795 7 11 0.636
anm372797 2 3 0.667
anm372904 4 5 0.800
anm372905 1 4 0.250
anm372906 5 7 0.714
anm372907 2 4 0.500
anm372909 1 2 0.500
Variant without flagged spikes: (a) 27 of 44 (0.614), p = 0.0010; (b) cosine 0.497, p = 0.0110
Variant without is_good false trials: (a) 27 of 44 (0.614), p = 0.0010; (b) cosine 0.495, p = 0.0170 Range of plotted d': sample -1.21 to 1.61; late delay -2.32 to 2.96
# Variants: one session per animal, and the paper's units from the processed data.
rng = np.random.default_rng(13)
for name, recs in (("one session per animal", records_one), ("paper's units (processed data)", proc_records)):
f, k, n, _, p = shuffle_fraction(recs, ("PTupper",), "late delay", rng)
o, nl, ps, ds, dl = similarity_test(recs, "PTupper", rng)
r, p_r = stats.pearsonr(ds, dl)
print(f"Variant {name}: (a) {k} of {n} ({f:.3f}), p = {p:.4f}; (b) cosine {o:.3f}, null 95th percentile "
f"{np.quantile(nl, 0.95):.3f}, p = {ps:.4f}; Pearson r = {r:.3f} (p = {p_r:.4g})")
Variant one session per animal: (a) 12 of 20 (0.600), p = 0.0010; (b) cosine 0.670, null 95th percentile 0.547, p = 0.0060; Pearson r = 0.709 (p = 0.0004691)
Variant paper's units (processed data): (a) 38 of 61 (0.623), p = 0.0010; (b) cosine 0.514, null 95th percentile 0.379, p = 0.0020; Pearson r = 0.518 (p = 1.929e-05)
Result: reproduced. (a) 27 of 44 PTupper units are selective in the late delay (null 99th percentile 0.046, p = 0.001). (b) The cosine similarity of the sample and late delay d' vectors is 0.495, against a shuffle null median of 0.182 and 95th percentile of 0.433 (p = 0.016); the code's Pearson correlation gives r = 0.496 (p = 0.0006). Of the 27 sample-selective units, 17 keep the sign of their d' in the late delay. Both parts hold without flagged spikes (p = 0.001 and 0.011), without the trials marked not good (p = 0.001 and 0.017), with one session per animal ((a) 12 of 20, p = 0.001; (b) cosine 0.670, p = 0.006) and on the paper's 61 PTupper units ((a) 38 of 61, p = 0.001; (b) cosine 0.514, p = 0.002).
pt-lower-late-and-motor¶
The medulla-projecting PT type mainly produces late preparatory activity and motor commands.
Criterion: among included PTlower units, (a) the fraction selective in the late delay and (b) the fraction selective in the response epoch are each larger than the fraction selective in the sample epoch, by exact one-sided McNemar tests at p < 0.05. Both: reproduced; one: partially reproduced.
def mcnemar_table(table, cell_type, first, second, threshold=0.01):
t = table[table.cell_type == cell_type]
a, b = t[f"p {first}"] < threshold, t[f"p {second}"] < threshold
p = alm.mcnemar_one_sided(int((a & ~b).sum()), int((~a & b).sum()))
return a.mean(), b.mean(), int((a & ~b).sum()), int((~a & b).sum()), len(t), p
for later in ("late delay", "response"):
fa, fb, only_a, only_b, n, p = mcnemar_table(main, "PTlower", "sample", later)
print(f"PTlower: sample {fa:.3f} against {later} {fb:.3f} (n = {n}); selective only in sample {only_a}, "
f"only in {later} {only_b}; one-sided McNemar p = {p:.4g}")
low = main[main.cell_type == "PTlower"]
print("\nBy animal (fraction selective, p < 0.01):")
print(low.groupby("subject")[["p sample", "p late delay", "p response"]].agg(lambda p: (p < 0.01).mean()).round(3)
.join(low.groupby("subject").size().rename("units")).to_string())
for name, table in (("without flagged spikes", variant_overlap), ("without is_good false trials", variant_good)):
for later in ("late delay", "response"):
fa, fb, only_a, only_b, n, p = mcnemar_table(table, "PTlower", "sample", later)
print(f"Variant {name}: sample {fa:.3f} against {later} {fb:.3f}; p = {p:.4g}")
for thr in (0.05, 0.001):
for later in ("late delay", "response"):
fa, fb, only_a, only_b, n, p = mcnemar_table(main, "PTlower", "sample", later, thr)
print(f"Threshold p < {thr}: sample {fa:.3f} against {later} {fb:.3f}; McNemar p = {p:.4g}")
# Variant reading of "motor commands": a change of rate in the response epoch from the baseline
# epoch, either direction, over all correct trials (paired Wilcoxon signed-rank, p < 0.01).
rows = []
for r in records:
if r["included"]:
diff = r["rates"]["response"] - r["rates"]["baseline"]
p = stats.wilcoxon(diff).pvalue if np.any(diff != 0) else 1.0
rows.append({"cell_type": r["cell_type"], "changed": p < 0.01, "up": np.mean(diff) > 0})
change = pd.DataFrame(rows)
print("\nVariant: response epoch rate differs from baseline (p < 0.01):")
print(change.groupby("cell_type").agg(fraction=("changed", "mean"), n=("changed", "size"),
increased=("up", lambda u: (u & change.loc[u.index, "changed"]).sum())).round(3).to_string())
k_low, n_low = int(change[change.cell_type == "PTlower"].changed.sum()), int((change.cell_type == "PTlower").sum())
k_up, n_up = int(change[change.cell_type == "PTupper"].changed.sum()), int((change.cell_type == "PTupper").sum())
print(f"PTlower {k_low}/{n_low} against PTupper {k_up}/{n_up}; one-sided Fisher p (PTlower > PTupper) = "
f"{stats.fisher_exact([[k_low, n_low - k_low], [k_up, n_up - k_up]], alternative='greater').pvalue:.4g}")
PTlower: sample 0.482 against late delay 0.714 (n = 56); selective only in sample 7, only in late delay 20; one-sided McNemar p = 0.009579
PTlower: sample 0.482 against response 0.643 (n = 56); selective only in sample 13, only in response 22; one-sided McNemar p = 0.08773
By animal (fraction selective, p < 0.01):
p sample p late delay p response units
subject
anm369962 0.300 0.900 0.850 20
anm369963 0.615 0.769 0.615 13
anm369964 0.667 0.556 0.444 9
anm372794 0.500 0.500 0.500 14
Variant without flagged spikes: sample 0.482 against late delay 0.714; p = 0.009579
Variant without flagged spikes: sample 0.482 against response 0.643; p = 0.08773
Variant without is_good false trials: sample 0.482 against late delay 0.714; p = 0.009579
Variant without is_good false trials: sample 0.482 against response 0.643; p = 0.08773
Threshold p < 0.05: sample 0.589 against late delay 0.750; McNemar p = 0.05388
Threshold p < 0.05: sample 0.589 against response 0.714; McNemar p = 0.1405
Threshold p < 0.001: sample 0.357 against late delay 0.607; McNemar p = 0.003305
Threshold p < 0.001: sample 0.357 against response 0.589; McNemar p = 0.01754
Variant: response epoch rate differs from baseline (p < 0.01):
fraction n increased
cell_type
PTlower 0.946 56 13
PTupper 0.841 44 10
PTlower 53/56 against PTupper 37/44; one-sided Fisher p (PTlower > PTupper) = 0.07957
# Variants: one session per animal, and the paper's units from the processed data.
for name, table in (("one session per animal", one), ("paper's units (processed data)", proc_main)):
for thr in (0.01, 0.05, 0.001):
for later in ("late delay", "response"):
fa, fb, only_a, only_b, n, p = mcnemar_table(table, "PTlower", "sample", later, thr)
print(f"Variant {name}, p < {thr}: sample {fa:.3f} against {later} {fb:.3f} (n = {n}); "
f"only sample {only_a}, only {later} {only_b}; McNemar p = {p:.4g}")
low_p = proc_main[proc_main.cell_type == "PTlower"]
print("\nPaper's units, PTlower by animal (fraction selective, p < 0.01):")
print(low_p.groupby("subject")[["p sample", "p late delay", "p response"]].agg(lambda p: (p < 0.01).mean()).round(3)
.join(low_p.groupby("subject").size().rename("units")).to_string())
Variant one session per animal, p < 0.01: sample 0.250 against late delay 0.625 (n = 16); only sample 0, only late delay 6; McNemar p = 0.01562
Variant one session per animal, p < 0.01: sample 0.250 against response 0.438 (n = 16); only sample 2, only response 5; McNemar p = 0.2266
Variant one session per animal, p < 0.05: sample 0.312 against late delay 0.688 (n = 16); only sample 1, only late delay 7; McNemar p = 0.03516
Variant one session per animal, p < 0.05: sample 0.312 against response 0.625 (n = 16); only sample 2, only response 7; McNemar p = 0.08984
Variant one session per animal, p < 0.001: sample 0.188 against late delay 0.562 (n = 16); only sample 0, only late delay 6; McNemar p = 0.01562
Variant one session per animal, p < 0.001: sample 0.188 against response 0.438 (n = 16); only sample 2, only response 6; McNemar p = 0.1445
Variant paper's units (processed data), p < 0.01: sample 0.507 against late delay 0.652 (n = 69); only sample 10, only late delay 20; McNemar p = 0.04937
Variant paper's units (processed data), p < 0.01: sample 0.507 against response 0.551 (n = 69); only sample 18, only response 21; McNemar p = 0.3746
Variant paper's units (processed data), p < 0.05: sample 0.580 against late delay 0.710 (n = 69); only sample 9, only late delay 18; McNemar p = 0.06104
Variant paper's units (processed data), p < 0.05: sample 0.580 against response 0.667 (n = 69); only sample 13, only response 19; McNemar p = 0.1885
Variant paper's units (processed data), p < 0.001: sample 0.362 against late delay 0.580 (n = 69); only sample 8, only late delay 23; McNemar p = 0.005337
Variant paper's units (processed data), p < 0.001: sample 0.362 against response 0.536 (n = 69); only sample 13, only response 25; McNemar p = 0.03648
Paper's units, PTlower by animal (fraction selective, p < 0.01):
p sample p late delay p response units
subject
anm369962 0.462 0.846 0.692 26
anm369963 0.562 0.625 0.625 16
anm369964 0.667 0.556 0.333 9
anm372794 0.444 0.444 0.389 18
Result: partially reproduced. Among 56 PTlower units, more are selective in the late delay than in the sample epoch (0.714 against 0.482; one-sided McNemar p = 0.0096), but more being selective in the response epoch (0.643) is not significant (p = 0.088). Nearly half of the PTlower units are already selective in the sample epoch. The pattern rests mostly on one animal: anm369962 goes from 0.30 in the sample epoch to 0.90 in the late delay, while the other three stay between 0.44 and 0.77 in all three epochs without a common rise. With the test at p < 0.001 both parts hold (p = 0.0033 and 0.018); at p < 0.05 neither does (p = 0.054 and 0.14). The rate in the response epoch differs from baseline in 53 of 56 PTlower units, and in 37 of 44 PTupper units (one-sided Fisher p = 0.080 for PTlower above PTupper). The variants give the same split. With one session per animal (16 units), part (a) holds (0.250 against 0.625, p = 0.016) and part (b) does not (0.438, p = 0.23). On the paper's 69 PTlower units, part (a) holds (0.507 against 0.652, p = 0.049) and part (b) does not (0.551, p = 0.37). On those units both parts hold at p < 0.001 (p = 0.0053 and 0.036), and neither does at p < 0.05 (p = 0.061 and 0.19).
pt-types-differ-in-timing¶
The two PT types differ in when they carry information about the movement; the thalamus-projecting type carries early preparatory activity, the medulla-projecting type mainly late preparatory activity and motor commands.
Criterion: (a) the fraction selective in the sample epoch is larger for PTupper than PTlower; (b) the fraction selective in the response epoch is larger for PTlower than PTupper; one-sided Fisher exact tests at p < 0.05. Both: reproduced; one: partially reproduced.
def fisher(table, epoch, larger, smaller, threshold=0.01):
a = table[table.cell_type == larger][f"p {epoch}"] < threshold
b = table[table.cell_type == smaller][f"p {epoch}"] < threshold
p = stats.fisher_exact([[a.sum(), (~a).sum()], [b.sum(), (~b).sum()]], alternative="greater").pvalue
return a.mean(), a.sum(), len(a), b.mean(), b.sum(), len(b), p
for label, table in (("main", main), ("without flagged spikes", variant_overlap), ("without is_good false trials", variant_good)):
fa, ka, na, fb, kb, nb, pa = fisher(table, "sample", "PTupper", "PTlower")
ga, la, ma, gb, lb, mb, pb = fisher(table, "response", "PTlower", "PTupper")
print(f"{label}: (a) sample PTupper {ka}/{na} ({fa:.3f}) against PTlower {kb}/{nb} ({fb:.3f}), p = {pa:.4g}; "
f"(b) response PTlower {la}/{ma} ({ga:.3f}) against PTupper {lb}/{mb} ({gb:.3f}), p = {pb:.4g}")
print("\nAll epochs, one-sided Fisher p for PTupper > PTlower and for PTlower > PTupper:")
for e in alm.EPOCHS:
_, _, _, _, _, _, p_up = fisher(main, e, "PTupper", "PTlower")
_, _, _, _, _, _, p_low = fisher(main, e, "PTlower", "PTupper")
print(f" {e:10s} PTupper {(main[main.cell_type == 'PTupper'][f'p {e}'] < 0.01).mean():.3f} "
f"PTlower {(main[main.cell_type == 'PTlower'][f'p {e}'] < 0.01).mean():.3f} p(up>low) {p_up:.4g} p(low>up) {p_low:.4g}")
by_animal = main.groupby(["cell_type", "subject"])[["p sample", "p response"]].agg(lambda p: (p < 0.01).mean())
by_animal["units"] = main.groupby(["cell_type", "subject"]).size()
print("\nBy animal:")
print(by_animal.round(3).to_string())
for e, larger, smaller in (("sample", "PTupper", "PTlower"), ("response", "PTlower", "PTupper")):
x, y = by_animal.loc[larger][f"p {e}"], by_animal.loc[smaller][f"p {e}"]
print(f"Animals, {e}: {larger} median {x.median():.3f} (n {len(x)}) against {smaller} {y.median():.3f} (n {len(y)}); "
f"one-sided Mann-Whitney p = {stats.mannwhitneyu(x, y, alternative='greater').pvalue:.4g}")
main: (a) sample PTupper 27/44 (0.614) against PTlower 27/56 (0.482), p = 0.134; (b) response PTlower 36/56 (0.643) against PTupper 37/44 (0.841), p = 0.9936
without flagged spikes: (a) sample PTupper 26/44 (0.591) against PTlower 27/56 (0.482), p = 0.1895; (b) response PTlower 36/56 (0.643) against PTupper 37/44 (0.841), p = 0.9936
without is_good false trials: (a) sample PTupper 27/44 (0.614) against PTlower 27/56 (0.482), p = 0.134; (b) response PTlower 36/56 (0.643) against PTupper 37/44 (0.841), p = 0.9936
All epochs, one-sided Fisher p for PTupper > PTlower and for PTlower > PTupper:
baseline PTupper 0.045 PTlower 0.018 p(up>low) 0.4095 p(low>up) 0.9181
sample PTupper 0.614 PTlower 0.482 p(up>low) 0.134 p(low>up) 0.935
delay PTupper 0.659 PTlower 0.768 p(up>low) 0.9229 p(low>up) 0.164
late delay PTupper 0.614 PTlower 0.714 p(up>low) 0.899 p(low>up) 0.198
response PTupper 0.841 PTlower 0.643 p(up>low) 0.02213 p(low>up) 0.9936
By animal:
p sample p response units
cell_type subject
PTlower anm369962 0.300 0.850 20
anm369963 0.615 0.615 13
anm369964 0.667 0.444 9
anm372794 0.500 0.500 14
PTupper anm372793 0.500 0.750 8
anm372795 0.727 1.000 11
anm372797 0.667 0.333 3
anm372904 0.800 0.600 5
anm372905 0.250 1.000 4
anm372906 0.429 1.000 7
anm372907 1.000 0.750 4
anm372909 0.500 1.000 2
Animals, sample: PTupper median 0.583 (n 8) against PTlower 0.558 (n 4); one-sided Mann-Whitney p = 0.3036
Animals, response: PTlower median 0.558 (n 4) against PTupper 0.875 (n 8); one-sided Mann-Whitney p = 0.9295
# Variants: one session per animal, and the paper's units from the processed data.
for label, table in (("one session per animal", one), ("paper's units (processed data)", proc_main)):
fa, ka, na, fb, kb, nb, pa = fisher(table, "sample", "PTupper", "PTlower")
ga, la, ma, gb, lb, mb, pb = fisher(table, "response", "PTlower", "PTupper")
_, _, _, _, _, _, pr = fisher(table, "response", "PTupper", "PTlower")
print(f"Variant {label}: (a) sample PTupper {ka}/{na} ({fa:.3f}) against PTlower {kb}/{nb} ({fb:.3f}), p = {pa:.4g}; "
f"(b) response PTlower {la}/{ma} ({ga:.3f}) against PTupper {lb}/{mb} ({gb:.3f}), p = {pb:.4g} (reverse p = {pr:.4g})")
# The measure of Extended Data Fig. 9 (decodingAccuracy.m): each unit's best single-threshold
# accuracy in classifying correct lick-left and lick-right trials from its sample epoch rate;
# also with rates rectified at the baseline (the mean over both trial types of the baseline epoch
# rate). The legend gives 15 of 61 PTupper against 3 of 69 PTlower at 70% or more, and a mean
# accuracy of 64.4 +/- 1.0% against 58.9 +/- 0.6% (two-sided Mann-Whitney p = 1e-4); rectified,
# 15 of 61 against 0 of 69 (62.7 +/- 1.1% against 56.7 +/- 0.4%, p = 9e-7).
def accuracies(recs, rectified):
rows = []
for r in recs:
if not r["included"]:
continue
x = r["rates"]["sample"]
if rectified:
base = r["rates"]["baseline"]
x = np.maximum(x - np.mean([base[~r["is_right"]].mean(), base[r["is_right"]].mean()]), 0)
rows.append({"cell_type": r["cell_type"], "subject": r["subject"], "acc": alm.roc_accuracy(x, r["is_right"])})
return pd.DataFrame(rows)
for label, recs in (("dandiset", records), ("paper's units (processed data)", proc_records)):
for rectified in (False, True):
acc = accuracies(recs, rectified)
up, low = acc[acc.cell_type == "PTupper"].acc, acc[acc.cell_type == "PTlower"].acc
k_up, k_low = int((up >= 0.7).sum()), int((low >= 0.7).sum())
p_f = stats.fisher_exact([[k_up, len(up) - k_up], [k_low, len(low) - k_low]], alternative="greater").pvalue
p_mw = stats.mannwhitneyu(up, low, alternative="two-sided").pvalue
print(f"{label}{', rectified' if rectified else ''}: accuracy >= 0.70 PTupper {k_up}/{len(up)} ({k_up / len(up):.3f}) "
f"against PTlower {k_low}/{len(low)} ({k_low / len(low):.3f}), one-sided Fisher p = {p_f:.4g}; "
f"mean {up.mean():.3f} +/- {up.sem():.3f} against {low.mean():.3f} +/- {low.sem():.3f}, Mann-Whitney p = {p_mw:.3g}; "
f"range {acc.acc.min():.3f} to {acc.acc.max():.3f}")
Variant one session per animal: (a) sample PTupper 8/20 (0.400) against PTlower 4/16 (0.250), p = 0.2783; (b) response PTlower 7/16 (0.438) against PTupper 16/20 (0.800), p = 0.9957 (reverse p = 0.02831) Variant paper's units (processed data): (a) sample PTupper 40/61 (0.656) against PTlower 35/69 (0.507), p = 0.06245; (b) response PTlower 38/69 (0.551) against PTupper 47/61 (0.770), p = 0.9977 (reverse p = 0.00692)
dandiset: accuracy >= 0.70 PTupper 7/44 (0.159) against PTlower 1/56 (0.018), one-sided Fisher p = 0.01248; mean 0.636 +/- 0.013 against 0.585 +/- 0.006, Mann-Whitney p = 0.000537; range 0.527 to 0.902 dandiset, rectified: accuracy >= 0.70 PTupper 6/44 (0.136) against PTlower 0/56 (0.000), one-sided Fisher p = 0.005922; mean 0.617 +/- 0.013 against 0.560 +/- 0.004, Mann-Whitney p = 2.7e-05; range 0.515 to 0.902
paper's units (processed data): accuracy >= 0.70 PTupper 12/61 (0.197) against PTlower 1/69 (0.014), one-sided Fisher p = 0.0004846; mean 0.638 +/- 0.011 against 0.586 +/- 0.006, Mann-Whitney p = 3.61e-05; range 0.510 to 0.902 paper's units (processed data), rectified: accuracy >= 0.70 PTupper 12/61 (0.197) against PTlower 0/69 (0.000), one-sided Fisher p = 6.045e-05; mean 0.626 +/- 0.011 against 0.567 +/- 0.004, Mann-Whitney p = 2.27e-06; range 0.510 to 0.902
# The code's comparison (CDCorrelation.m): dot product of the unit-normalised sample and late
# delay coding directions per type, over 100 bootstraps that resample animals, then units within
# them, then trials within units.
rng = np.random.default_rng(5)
def bootstrap_similarity(recs, cell_type, rng, n_boot=100):
pool = [r for r in recs if r["included"] and r["cell_type"] == cell_type]
animals = sorted({r["subject"] for r in pool})
out = []
for _ in range(n_boot):
chosen = rng.choice(animals, len(animals), replace=True)
cand = [i for a in chosen for i, r in enumerate(pool) if r["subject"] == a]
picks = rng.choice(cand, len(pool), replace=True)
d_s, d_l = [], []
for i in picks:
r = pool[i]
right = r["is_right"]
idx_l = rng.choice(np.flatnonzero(~right), (~right).sum())
idx_r = rng.choice(np.flatnonzero(right), right.sum())
idx = np.r_[idx_l, idx_r]
lab = np.r_[np.zeros(len(idx_l), bool), np.ones(len(idx_r), bool)]
d_s.append(alm.dprime(r["rates"]["sample"][idx], lab))
d_l.append(alm.dprime(r["rates"]["late delay"][idx], lab))
out.append(alm.cosine(np.array(d_s), np.array(d_l)))
return np.array(out)
boots = {t: bootstrap_similarity(records, t, rng) for t in alm.TYPES}
for t, b in boots.items():
print(f"{t}: sample to late delay similarity {np.nanmean(b):.3f} +/- {np.nanstd(b):.3f} (mean +/- sd over {len(b)} bootstraps)")
print(f"Fraction of bootstraps with PTupper above PTlower: {np.mean(boots['PTupper'] > boots['PTlower']):.2f}")
for t in alm.TYPES:
o, _, p, _, _ = similarity_test(records, t, rng)
print(f"{t}: cosine without bootstrap {o:.3f}, label-shuffle p = {p:.4f}")
PTupper: sample to late delay similarity 0.473 +/- 0.214 (mean +/- sd over 100 bootstraps) PTlower: sample to late delay similarity 0.216 +/- 0.204 (mean +/- sd over 100 bootstraps) Fraction of bootstraps with PTupper above PTlower: 0.83
PTupper: cosine without bootstrap 0.495, label-shuffle p = 0.0070
PTlower: cosine without bootstrap 0.177, label-shuffle p = 0.4276
# Fraction of included units selective (p < 0.01) in 200 ms windows, every 100 ms, by type.
starts = np.round(np.arange(-3.3, 1.01, 0.1), 2)
curves = {}
for t in alm.TYPES:
pool = [r for r in records if r["included"] and r["cell_type"] == t]
curves[t] = np.array([np.mean([alm.ranksum_p(alm.epoch_rates(r["spike_times"], r["cue"], (a, a + 0.2)), r["is_right"]) < 0.01
for r in pool]) for a in starts])
centres = starts + 0.1
fig, ax = plt.subplots(figsize=(7, 3.5))
colours = {"PTupper": "#2a78d6", "PTlower": "#eb6834"}
for t in alm.TYPES:
n = sum(r["included"] and r["cell_type"] == t for r in records)
ax.plot(centres, curves[t], color=colours[t], lw=2, label=f"{t} ({alm.PROJECTION[t]}, n = {n})")
ax.axhline(0.01, color="#52514e", lw=1, ls=":", label="test level (0.01)")
for x in (-2.5, -1.5, -1.3, 0):
ax.axvline(x, color="#c3c2b7", lw=1)
ax.set_xlabel("Time from go cue (s); lines at pole in, pole out, delay start, go cue")
ax.set_ylabel("Fraction selective")
ax.legend(frameon=False, fontsize=8)
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
plt.show()
print("Window centres", centres.min(), "to", centres.max(), "s; all points within the axes.")
for t in alm.TYPES:
print(f"{t}: fraction range {curves[t].min():.3f} to {curves[t].max():.3f}; peak at {centres[np.argmax(curves[t])]:.1f} s")
print(pd.DataFrame(curves, index=centres).iloc[::5].round(3).to_string())
Window centres -3.1999999999999997 to 1.1 s; all points within the axes.
PTupper: fraction range 0.000 to 0.841; peak at 0.2 s
PTlower: fraction range 0.000 to 0.714; peak at 0.0 s
PTupper PTlower
-3.2 0.023 0.000
-2.7 0.000 0.000
-2.2 0.500 0.554
-1.7 0.455 0.589
-1.2 0.614 0.482
-0.7 0.568 0.679
-0.2 0.636 0.661
0.3 0.727 0.661
0.8 0.659 0.661
Result: not reproduced. (a) PTupper units are selective in the sample epoch more often than PTlower units (0.614 against 0.482), but not significantly (Fisher p = 0.134). (b) In the response epoch PTlower units are selective less often than PTupper units (0.643 against 0.841), the reverse of the claim (p = 0.99 for the stated direction, 0.022 for the reverse). Over animals neither difference holds (Mann-Whitney p = 0.30 and 0.93). The code's coding-direction comparison points the way the abstract describes, but the difference is not significant: the sample to late delay similarity is 0.473 +/- 0.214 for PTupper and 0.216 +/- 0.204 for PTlower, with PTupper higher in only 83 of 100 bootstraps. Without the bootstrap, PTupper's similarity is above its shuffle null (p = 0.007) and PTlower's is not (0.177, p = 0.43). With one session per animal, (a) is not significant (0.400 against 0.250, p = 0.28) and (b) is reversed (0.438 against 0.800). On the paper's 130 units from the processed data, (a) is not significant (0.656 against 0.507, p = 0.062) and (b) is reversed (0.551 against 0.770; p = 0.007 for the reverse), so the reversal does not come from the units missing from the dandiset. The paper's own measure of early information, single-unit decoding accuracy in the sample epoch (Extended Data Fig. 9), goes the abstract's way. At 70% accuracy or more, it finds 7 of 44 PTupper against 1 of 56 PTlower units here (Fisher p = 0.012), and 12 of 61 against 1 of 69 on the paper's units (p = 0.0005), where the legend gives 15 of 61 against 3 of 69. The mean accuracy on the paper's units is 0.638 against 0.586 (legend 0.644 and 0.589). With rates rectified at baseline, it finds 12 of 61 against 0 of 69 (legend 15 of 61 against 0 of 69). The means match the legend, but this attempt does not explain the smaller counts. These variants support part (a) under the paper's measure but not under the criterion's. For part (b), the variant reading of "motor commands" in the result above, a response epoch rate that differs from baseline, goes the claim's way without reaching significance (53 of 56 PTlower against 37 of 44 PTupper, Fisher p = 0.080). The reversal of part (b) thus shows that PTlower units are less often selective for lick direction in the response epoch; it does not show that they lack activity related to the movement. The variants do not change the outcome, which the criterion decides. Each animal has one type only, so a difference between types is also a difference between 8 and 4 animals.
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 |
|---|---|---|
| pt-types-transcriptomic | dandiset | blocked: data missing |
| pt-types-shared-targets | dandiset | blocked: data missing |
| pt-upper-thalamus-projection | dandiset | blocked: data missing |
| pt-lower-medulla-projection | dandiset | blocked: data missing |
| pt-preparatory-activity | dandiset | reproduced |
| pt-upper-early-preparatory | dandiset | reproduced |
| pt-upper-persistent | dandiset | reproduced |
| pt-lower-late-and-motor | dandiset | partially reproduced |
| pt-types-differ-in-timing | dandiset | not reproduced |
How the analysis differs from the paper's. The claims and their criteria were written from the abstract only, so they and the tests of significance are this attempt's; the paper's own statistics are unknown. The public Extended Data legends 5 to 10, read later, are used in variants only. The epochs, trial selection, unit inclusion and selectivity test follow the paper's code on Figshare. Epoch rates are spike counts over the epoch, where the code averages a causally smoothed PSTH. Spikes are counted on the file's clock, where some spikes of one trial fall into the next trial's early epochs; the code worked on spikes already assigned to trials. Leaving out the spikes that the stored order shows to overlap changes no outcome. The tests of chance (label shuffles), of change between epochs (McNemar) and between types (Fisher) are choices made here. The dandiset holds 44 of the paper's 61 PTupper and 56 of its 69 PTlower units; variants rerun every test on the paper's 130 units from its processed data and give the same outcomes.
What was found about the dandiset.
- It holds the electrophysiology only: the transcriptomic, anatomical and axonal data behind four of the abstract's findings are elsewhere (GEO GSE119182 and the Figshare record).
- Spike times are not sorted in 1120 of 1149 units. They were written as trial start + go cue time
- the spike's time from the go cue, trial by trial; as each trial's spikes cover 6.4 s and trial starts are a median 6.4 s apart (5.4 to 7.9 s from the 1st to the 99th percentile) while the pole comes in at different times, consecutive trials overlap on the file's clock.
- 11 sessions of the paper's processed data that hold tagged units are not in the dandiset. They hold 17 of the paper's 61 PTupper and 13 of its 69 PTlower units. The description says some files were broken and left out, without naming them.
- In every session the spike windows of the first ten trials start later, relative to the pole coming in, than those of the other trials, in the files as in the processed data.
stop_timeof the trial table is NaN in every file; the times of the trial events are relative to the trial start, which the files do not say.- The projection class of the tagged types is not stated in the files;
PTupperandPTlowerare linked to thalamus- and medulla-projecting neurons only through the paper's code. - Unit quality mixes vocabularies and has undescribed empty values; sex is unknown and genotype empty for every mouse; electrode positions are all zero, so depth comes only from the unit table.