Siegle, Jia et al. (2021): the passive viewing survey, from the paper's DANDI deposits¶
Paper. Siegle JH, Jia X, Durand S, et al. (2021). Survey of spiking in the mouse visual system reveals functional hierarchy. Nature 592, 86-92. doi:10.1038/s41586-020-03171-x. The claims file was written from the NIH author manuscript in PubMed Central (PMC10399640): main text, Methods, the legends of Figures 1 to 4 and Extended Data Figures 1 to 10, and the availability statements.
Target dandiset. DANDI:000253,
version 0.240923.1441, the most recent published version (the paper cites no version of it; it
predates the dandiset). Its title and description are those of the OpenScope Global/Local Oddball
project; it names this paper as IsSupplementTo. The first section shows that it holds a
different experiment, so every claim the claims file marks elsewhere is tested on the archives
the paper names for its data.
Where the paper's data are. The Data availability statement deposits the 58 passive viewing experiments on DANDI as dandisets 000021 and 000022. Pinned here: 000021 version 0.251116.2246 (Brain Observatory 1.1 stimulus set) and 000022 version 0.251116.2247 (Functional Connectivity stimulus set). The paper cites no version; these are the only published versions of each (the other is the draft), so they are pinned. The per-unit tables and the behaviour metrics tables of the paper's GitHub repository (AllenInstitute/neuropixels_platform_paper at commit 4f8bcc26) are read for the claims that the claims file says rest on them, and for the link claim.
This notebook tests the claims that do not need cross-correlograms. The CCG claims (Figure 2 and
Extended Data Figures 7, 8 and 10c to e) are in ccg_siegle_2021.ipynb.
| Claim | Data | Here |
|---|---|---|
| methods-sessions-passing-qc | 000021, 000022 | tested |
| methods-probes-ccf | 000021, 000022 | tested |
| ed1jk-units-per-experiment | 000021, 000022 | tested |
| fig1g-receptive-fields | 000021, 000022 | tested |
| link-unit-tables (added link) | 000021, 000022 and the paper's per-unit tables | tested |
| fig3c, ed9a, fig3f, ed9d, fig3i, ed9e, ed9fg, fig3l, ed9h, fig3-unit-counts, ed9bc, ed9ij | 000021, 000022 | tested |
| ed10b-layers | 000021, 000022 with the paper's layer table | tested |
| anatomical-hierarchy-order | the paper's code | tested (order only) |
| fig4e, fig4f, ed9k, ed9lmn, fig4h, ed9p | the paper's behaviour metrics tables | tested |
| methods-gaze-stability | 000021, 000022 | tested |
| fig4b-behaviour-performance | not deposited | blocked, data missing |
| CCG claims (nine) | 000021, 000022 | in ccg_siegle_2021.ipynb |
import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import io
import re
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
from scipy import stats
from statsmodels.stats.multitest import multipletests
import siegle as S
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()
def repo_csv(name):
return pd.read_csv(io.StringIO(requests.get(f"{S.CODE_RAW}/data/{name}", timeout=60).text))
The target dandiset, 000253¶
Its assets and metadata from the DANDI API, and the contents of one file, read without its data.
tgt = S.list_assets(*S.TARGET)
meta253 = requests.get(f"{S.API}/dandisets/{S.TARGET[0]}/versions/{S.TARGET[1]}/").json()
print("name:", meta253["name"])
print("related:", [(r["relation"], r.get("name")) for r in meta253.get("relatedResource", [])])
print("assets:", len(tgt), " subjects:", tgt.path.str.split("/").str[0].nunique(),
" sessions:", tgt.path.str.extract(r"ses-(\d+)")[0].nunique(),
" size (TB):", round(tgt.size_gb.sum() / 1000, 2))
print("files by kind:", tgt.path.str.extract(r"_(ogen|probe-\d)")[0].str.replace(r"probe-\d", "probe", regex=True).value_counts().to_dict())
f253 = S.open_asset(tgt[tgt.path.str.endswith("_ogen.nwb")].iloc[0].asset_id)
print("\nstimulus:", f253["general/stimulus"][()].decode())
print("interval tables:", sorted(f253["intervals"].keys()))
print("processing:", sorted(f253["processing"].keys()))
print("units:", f253["units/id"].shape[0], "columns:", sorted(k for k in f253["units"].keys() if not k.endswith("_index")))
f253.close()
name: Allen Institute Openscope - Global/Local Oddball project
related: [('dcite:IsSupplementedBy', 'OpenScope DataBook'), ('dcite:IsSupplementedBy', 'OpenScope DataPortal'), ('dcite:IsSupplementTo', 'Survey of spiking in the mouse visual system reveals functional hierarchy ')]
assets: 98 subjects: 14 sessions: 14 size (TB): 0.68
files by kind: {'probe': 84, 'ogen': 14}
stimulus: OpenScopeGlobalLocalOddball interval tables: ['create_receptive_field_mapping_presentations', 'init_grating_presentations', 'init_intermission_presentations', 'invalid_times', 'spontaneous_presentations'] processing: ['optotagging', 'running', 'stimulus']
units: 2446 columns: ['PT_ratio', 'amplitude', 'amplitude_cutoff', 'cluster_id', 'cumulative_drift', 'd_prime', 'firing_rate', 'id', 'isi_violations', 'isolation_distance', 'l_ratio', 'max_drift', 'nn_hit_rate', 'nn_miss_rate', 'peak_channel_id', 'presence_ratio', 'quality', 'recovery_slope', 'repolarization_slope', 'silhouette_score', 'snr', 'spike_amplitudes', 'spike_times', 'spread', 'velocity_above', 'velocity_below', 'waveform_duration', 'waveform_halfwidth', 'waveform_mean']
The files of 000253 hold an oddball grating experiment (init_grating, receptive field mapping
with 1,944 presentations, spontaneous blocks) in 14 mice, with no full-field flashes, no 81-location
Gabor block of 3,645 presentations and no drifting gratings of 75 repeats, which the paper's
analyses use. The paper's recordings are not in it; the claims are tested on 000021 and 000022.
Loading the paper's 58 sessions¶
For each session file of 000021 and 000022: the subject, the units table with its peak channel's structure and CCF coordinates, the stimulus tables, the gaze in degrees, and, streamed for the non-noise units of the eight visual areas and the hippocampus only, the spike times, from which the helper computes each unit's receptive field test and area, times to first spike, drifting grating rate and modulation index, and the two timescales. Four sessions are read at a time; one session takes 30 to 60 s.
assets = pd.concat([S.list_assets(d, v) for d, v in S.SOURCES.items()])
sessions_assets = assets[~assets.path.str.contains("probe")].reset_index(drop=True)
print("assets per dandiset:", assets.groupby("dandiset").size().to_dict(),
" session files:", sessions_assets.groupby("dandiset").size().to_dict())
res = S.load_all(sessions_assets)
units = pd.concat([r[0] for r in res], ignore_index=True)
meta = pd.DataFrame([r[1] for r in res])
gaze = {r[1]["subject"]: r[2] for r in res}
print(f"loaded {len(meta)} sessions, {len(units)} units, in {(time.time() - started) / 60:.1f} min")
assets per dandiset: {'000021': 214, '000022': 169} session files: {'000021': 32, '000022': 26}
loaded 58 sessions, 136556 units, in 8.9 min
Data check¶
print("Sessions by dandiset and stimulus set:")
print(pd.crosstab(meta.dandiset, meta.session_type).to_string())
print("\nSubjects:", meta.subject.nunique(), " sessions:", len(meta))
print("\nGenotypes:")
print(meta.genotype.value_counts().to_string())
age_days = meta.age.str.extract(r"P(\d+)D")[0].astype(int)
print("\nSex:", meta.sex.value_counts().to_dict(), " age range (days):", age_days.min(), "to", age_days.max())
print("\nProbes per session:", meta.n_probes.value_counts().sort_index().to_dict(), " total:", meta.n_probes.sum())
print("\nUnits by quality:", units.quality.value_counts().to_dict())
print("Non-noise units by structure (eight areas and hippocampus):")
good = units[units.quality == "good"]
print(good.location.value_counts().reindex(S.AREAS + S.HPC).to_string())
print("\nUnits whose peak channel has no structure label:", int((units.location == "").sum()),
" no CCF coordinate:", int(units.x.isna().sum()))
Sessions by dandiset and stimulus set:
session_type brain_observatory_1.1 functional_connectivity
dandiset
000021 32 0
000022 0 26
Subjects: 58 sessions: 58
Genotypes:
genotype
wt/wt 30
Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt 12
Pvalb-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt 8
Vip-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt 8
Sex: {'M': 45, 'F': 13} age range (days): 91 to 142
Probes per session: {4: 2, 5: 12, 6: 44} total: 332
Units by quality: {'good': 99844, 'noise': 36712}
Non-noise units by structure (eight areas and hippocampus):
location
LGd 2591
VISp 8605
VISl 4941
VISrl 6031
LP 4873
VISal 6468
VISpm 4270
VISam 6211
CA1 17196
CA3 3211
DG 5882
Units whose peak channel has no structure label: 1333 no CCF coordinate: 21130
print("Stimulus tables per session (counts, median durations in s):")
cols = ["n_flashes", "flash_duration", "flash_interval", "n_gabors", "gabor_duration", "gabor_positions", "dg_stimulus", "n_dg"]
print(meta.groupby("session_type")[cols].agg(lambda s: sorted(set(np.round(s, 3) if s.dtype.kind == "f" else s))).T.to_string())
print("\nFlash colours (dark -1, light +1) per session:", meta.flash_colors.astype(str).value_counts().to_dict())
print("\nSessions whose electrode y and z columns are equal on every channel:", int(meta.y_equals_z.sum()), "of", len(meta))
print("CCF ranges over units (um): x", np.nanmin(units.x), "to", np.nanmax(units.x),
"; y", np.nanmin(units.y), "to", np.nanmax(units.y), "; z", np.nanmin(units.z), "to", np.nanmax(units.z))
print("\nSessions with gaze mapping:", int(meta.has_gaze.sum()))
m = units[units.p_value_rf.notna()]
print("\nUnits measured (non-noise, eight areas and hippocampus):", len(m))
desc = m[["p_value_rf", "area_rf", "tfs_all_250", "tfs_spont", "rate_pre", "firing_rate_dg", "mod_idx_dg",
"timescale_it", "timescale_ac", "firing_rate", "snr"]].describe().T
print(desc[["count", "min", "50%", "max"]].round(4).to_string())
measured = units[(units.quality == "good") & units.location.isin(S.AREAS + S.HPC)]
print("\nunits with fewer than 10 spikes over the Gabor presentations (no p_value_rf, as the AllenSDK):",
int((measured.n_spikes_gabor < 10).sum()), "of", len(measured))
print("\nSpikes read:", int(meta.spikes_read.sum()), " seconds per session:", round(meta.seconds.min()), "to", round(meta.seconds.max()))
Stimulus tables per session (counts, median durations in s):
session_type brain_observatory_1.1 functional_connectivity
n_flashes [150] [150]
flash_duration [0.25] [0.25]
flash_interval [2.002] [2.002]
n_gabors [3645] [3645]
gabor_duration [0.25] [0.25]
gabor_positions [81] [81]
dg_stimulus [drifting_gratings_presentations] [drifting_gratings_75_repeats_presentations]
n_dg [628, 630] [600]
Flash colours (dark -1, light +1) per session: {"{'-1.0': np.int64(75), '1.0': np.int64(75)}": 58}
Sessions whose electrode y and z columns are equal on every channel: 58 of 58
CCF ranges over units (um): x 6487.0 to 9727.0 ; y -1250.0 to 4890.0 ; z -1250.0 to 4890.0
Sessions with gaze mapping: 50
Units measured (non-noise, eight areas and hippocampus): 65945
count min 50% max
p_value_rf 65945.0 0.0000 0.0550 1.0000
area_rf 65945.0 100.0000 500.0000 6400.0000
tfs_all_250 62294.0 0.0300 0.1040 0.2490
tfs_spont 61625.0 0.0300 0.1110 0.2490
rate_pre 65945.0 0.0000 1.7067 185.2267
firing_rate_dg 65945.0 0.0000 3.0831 115.6250
mod_idx_dg 65785.0 0.0000 0.8220 9.7124
timescale_it 23442.0 1.0000 1.0002 1000.0000
timescale_ac 26128.0 1.0000 61.9964 1000.0000
firing_rate 65945.0 0.0019 3.2154 98.8355
snr 65945.0 0.0356 2.1246 10.3526
units with fewer than 10 spikes over the Gabor presentations (no p_value_rf, as the AllenSDK): 7619 of 70279
Spikes read: 3341589748 seconds per session: 23 to 85
Checks of the description against the values. The 58 sessions split 32 and 26 between the Brain Observatory 1.1 and Functional Connectivity sets, one session per mouse, as the Methods give. Every session has 150 full-field flashes (75 dark, 75 light) of 250 ms every 2 s, and a Gabor block of 81 locations with 250 ms presentations, as the Methods describe; the Functional Connectivity sessions have the drifting gratings of 4 orientations at 2 Hz with 75 repeats per contrast. Times are in seconds and spike times share the clock of the stimulus tables: the receptive field test finds significant fields in most visual area units and almost none in the hippocampus, and flash latencies cluster in tens of ms (below), which a clock offset would not give. The electrodes table gives the CCF position in x, y and z, but y and z are equal on every channel of every session, so one of the three coordinates (the left-right one, z) is not the CCF value; x (anterior-posterior) is used below only to say whether a probe was mapped.
methods-sessions-passing-qc¶
Of 87 mice entering the recording step, 61 passed session-level quality control and 58 resulted in NWB files; these 58 passive viewing experiments are the ones analysed.
Counts of sessions and of mice by genotype, from the session metadata (tolerance zero).
geno = meta.genotype.str.split("/").str[0].str.replace("-IRES-Cre", "", regex=False)
counts = {"sessions": len(meta), "wt": int((geno == "wt").sum()), "Pvalb": int((geno == "Pvalb").sum()),
"Sst": int((geno == "Sst").sum()), "Vip": int((geno == "Vip").sum())}
print(counts, " expected: sessions 58, wt 30, Pvalb 8, Sst 12, Vip 8")
print("Ai32 in every Cre line:", bool(meta[geno != "wt"].genotype.str.contains("Ai32").all()))
{'sessions': 58, 'wt': 30, 'Pvalb': 8, 'Sst': 12, 'Vip': 8} expected: sessions 58, wt 30, Pvalb 8, Sst 12, Vip 8
Ai32 in every Cre line: True
Result: reproduced. The two dandisets hold 58 sessions of 58 mice: 30 wild type, 8 Pvalb, 12 Sst and 8 Vip Cre lines crossed with Ai32, the counts of the Methods. The 87 and 61 are pipeline counts before release and are not judged.
methods-probes-ccf¶
284 of the 332 probes were mapped to the CCFv3.
As ExtDataFigure1jk.py counts it: per probe, the units passing ISI violations < 0.5, amplitude cutoff < 0.1 and presence ratio > 0.95 (no quality filter in that loop); a probe is mapped when the sum of their anterior-posterior CCF coordinates is positive (missing values skipped, as pandas does). Reading (a), the 58 released sessions; reading (b) cannot be counted on released data.
filt = (units.isi_violations < 0.5) & (units.amplitude_cutoff < 0.1) & (units.presence_ratio > 0.95)
probes = units.groupby(["session", "probe_id"]).size().rename("units").to_frame()
probes["units_filtered"] = units[filt].groupby(["session", "probe_id"]).size()
probes["ap_sum"] = units[filt].groupby(["session", "probe_id"]).x.sum()
probes = probes.fillna(0)
n_probes = int(meta.n_probes.sum())
print("probes (devices) in the 58 sessions:", n_probes, " probes with units:", len(probes))
print("probes mapped (sum of AP coordinates > 0):", int((probes.ap_sum > 0).sum()), " expected 284 of 332")
print("probes whose units all lack a CCF coordinate:", int((units.groupby(["session", "probe_id"]).x.count() == 0).sum()))
print("unmapped probes per session:", (probes.ap_sum <= 0).groupby(level=0).sum().value_counts().sort_index().to_dict())
probes (devices) in the 58 sessions: 332 probes with units: 332
probes mapped (sum of AP coordinates > 0): 284 expected 284 of 332
probes whose units all lack a CCF coordinate: 48
unmapped probes per session: {0: 46, 1: 4, 5: 4, 6: 4}
Result: reproduced. The 58 sessions hold 332 probes, and 284 of them have a positive sum of anterior-posterior CCF coordinates over their filtered units; the 48 others have no CCF coordinate on any unit. Both counts equal the paper's under reading (a).
ed1jk-units-per-experiment¶
After quality filtering, 682 +/- 144 units were recorded per experiment, 119 +/- 48 per probe and 56 +/- 30 per visual area, with 6.1 +/- 1.1 visual areas sampled per experiment.
As ExtDataFigure1jk.py: the three filters and quality good per experiment and per area, the three filters only per probe; areas with at least one unit; mean and numpy SD.
f3 = (units.amplitude_cutoff < 0.1) & (units.presence_ratio > 0.95) & (units.isi_violations < 0.5)
sel1 = f3 & (units.quality == "good")
per_exp = units[sel1].groupby("session").size().reindex(meta.path.str.extract(r"ses-(\d+)")[0]).fillna(0)
per_probe = units[f3].groupby(["session", "probe_id"]).size()
per_probe = per_probe.reindex(units.groupby(["session", "probe_id"]).size().index).fillna(0)
pa = units[sel1 & units.location.isin(S.AREAS)].groupby(["session", "location"]).size()
areas_per_exp = pa.groupby(level=0).size().reindex(per_exp.index).fillna(0)
out = pd.DataFrame({"mean": [per_exp.mean(), per_probe.mean(), pa.mean(), areas_per_exp.mean()],
"sd": [per_exp.std(ddof=0), per_probe.std(ddof=0), pa.std(ddof=0), areas_per_exp.std(ddof=0)],
"n": [len(per_exp), len(per_probe), len(pa), len(areas_per_exp)],
"paper": ["682 +/- 144", "119 +/- 48", "56 +/- 30", "6.1 +/- 1.1"]},
index=["units per experiment", "units per probe", "units per area per experiment", "areas per experiment"])
print(out.round(2).to_string())
print("\nunits per experiment by stimulus set:")
st_of = meta.set_index(meta.path.str.extract(r"ses-(\d+)")[0]).session_type
print(per_exp.groupby(st_of).agg(["mean", "std", "min", "max"]).round(1).to_string())
print("probes counted here:", len(per_probe), " (probes with any unit)")
mean sd n paper
units per experiment 686.17 144.58 58 682 +/- 144
units per probe 128.78 51.06 332 119 +/- 48
units per area per experiment 56.17 30.34 352 56 +/- 30
areas per experiment 6.07 1.13 58 6.1 +/- 1.1
units per experiment by stimulus set:
mean std min max
session_type
brain_observatory_1.1 678.9 144.3 439 954
functional_connectivity 695.1 150.1 411 999
probes counted here: 332 (probes with any unit)
Result: reproduced. Units per experiment 686.2 +/- 144.6 (paper 682 +/- 144), per probe 128.8 +/- 51.1 (119 +/- 48), per area per experiment 56.2 +/- 30.3 (56 +/- 30), and 6.07 +/- 1.13 areas per experiment (6.1 +/- 1.1); all five values are within their tolerances.
fig1g-receptive-fields¶
Units in the eight visual areas were highly visually responsive, 60% having a significant spatial receptive field within the monitor, against only 1.4% of units in the hippocampus (CA1, CA3 and dentate gyrus).
p < 0.01 by the chi-square shuffle test (1,000 shuffles, stopped early once p >= 0.01 is settled); all non-noise units, and, as the AllenSDK's p_value_rf (ReceptiveFieldMapping, minimum_spike_count 10), no p-value for a unit with fewer than 10 spikes over the Gabor presentations, which then counts as not significant, as Figure1g.py counts it (its denominator is every unit). The 60% pooled as Figure1g.py; the 1.4% under (a) pooled and (b) the mean across mice of each mouse's fraction. Variants: without the 10-spike rule (the round-1 result), and with the units without a p-value left out of the denominator.
g = units[(units.quality == "good")]
vis = g[g.location.isin(S.AREAS)]
hpc = g[g.location.isin(S.HPC)]
sig = lambda d, col="p_value_rf": (d[col] < 0.01)
print("eight areas, pooled: %d of %d = %.4f" % (sig(vis).sum(), len(vis), sig(vis).mean()))
print("hippocampus (a) pooled: %d of %d = %.4f" % (sig(hpc).sum(), len(hpc), sig(hpc).mean()))
per_mouse = hpc.groupby("subject").apply(lambda d: sig(d).mean())
print("hippocampus (b) mean over %d mice of each mouse's fraction = %.4f (median %.4f, max %.4f)"
% (len(per_mouse), per_mouse.mean(), per_mouse.median(), per_mouse.max()))
print("expected 0.60 +/- 0.06 and 0.014 +/- 0.0014")
print("\nunits with fewer than 10 Gabor spikes (no p-value):", int((vis.n_spikes_gabor < 10).sum()), "visual,",
int((hpc.n_spikes_gabor < 10).sum()), "hippocampal; of these significant without the rule:",
int(sig(vis[vis.n_spikes_gabor < 10], "p_value_rf_all").sum()), "visual,", int(sig(hpc[hpc.n_spikes_gabor < 10], "p_value_rf_all").sum()), "hippocampal")
pm_all = hpc.groupby("subject").apply(lambda d: sig(d, "p_value_rf_all").mean())
print("variant without the 10-spike rule: eight areas %.4f; hippocampus (a) %.4f, (b) %.4f"
% (sig(vis, "p_value_rf_all").mean(), sig(hpc, "p_value_rf_all").mean(), pm_all.mean()))
vv, hh = vis[vis.p_value_rf.notna()], hpc[hpc.p_value_rf.notna()]
print("variant, units without a p-value left out: eight areas %.4f (%d units); hippocampus (a) %.4f (%d units), (b) %.4f"
% (sig(vv).mean(), len(vv), sig(hh).mean(), len(hh), hh.groupby("subject").apply(lambda d: sig(d).mean()).mean()))
print("\nper area, pooled fraction and number of units:")
print(g[g.location.isin(S.AREAS + S.HPC)].groupby("location").apply(lambda d: pd.Series({"fraction": sig(d).mean(), "units": len(d)})).round(3).to_string())
print("threshold sensitivity, eight areas: p<0.001 %.3f, p<0.005 %.3f, p<0.02 %.3f" %
tuple((vis.p_value_rf < t).mean() for t in (0.001, 0.005, 0.02)))
print("threshold sensitivity, hippocampus (a): p<0.005 %.4f, p<0.02 %.4f" % tuple((hpc.p_value_rf < t).mean() for t in (0.005, 0.02)))
print("by stimulus set, eight areas:", vis.groupby("session_type").apply(lambda d: round(sig(d).mean(), 3)).to_dict())
eight areas, pooled: 26694 of 43990 = 0.6068 hippocampus (a) pooled: 415 of 26289 = 0.0158 hippocampus (b) mean over 52 mice of each mouse's fraction = 0.0171 (median 0.0138, max 0.0658) expected 0.60 +/- 0.06 and 0.014 +/- 0.0014 units with fewer than 10 Gabor spikes (no p-value): 4267 visual, 3352 hippocampal; of these significant without the rule: 45 visual, 9 hippocampal variant without the 10-spike rule: eight areas 0.6078; hippocampus (a) 0.0161, (b) 0.0174 variant, units without a p-value left out: eight areas 0.6720 (39723 units); hippocampus (a) 0.0181 (22937 units), (b) 0.0194 per area, pooled fraction and number of units:
fraction units
location
CA1 0.014 17196.0
CA3 0.024 3211.0
DG 0.016 5882.0
LGd 0.522 2591.0
LP 0.451 4873.0
VISal 0.634 6468.0
VISam 0.607 6211.0
VISl 0.706 4941.0
VISp 0.688 8605.0
VISpm 0.562 4270.0
VISrl 0.575 6031.0
threshold sensitivity, eight areas: p<0.001 0.554, p<0.005 0.590, p<0.02 0.625
threshold sensitivity, hippocampus (a): p<0.005 0.0103, p<0.02 0.0258
by stimulus set, eight areas: {'brain_observatory_1.1': 0.621, 'functional_connectivity': 0.59}
Result: partially reproduced. 60.7% of the 43,990 non-noise units of the eight areas have p < 0.01 (paper 60%, within tolerance). In the hippocampus 1.58% do pooled (reading a) and 1.71% on the mean over 52 mice (reading b); both are outside 1.4% +/- 0.14%, though the contrast with the visual areas is as large as the paper's. The AllenSDK's rule of at least 10 Gabor spikes removes 45 visual and 9 hippocampal units that were significant without it (60.8%, 1.61% and 1.74% without the rule). Leaving the units without a p-value out of the denominator gives 67.2%, 1.81% and 1.94%, so the judgement does not change either way.
link-unit-tables¶
The paper's per-unit tables (data/time_to_first_spike.csv, data/timescale_metrics.csv,
data/layer_info.csv) are keyed by ecephys_unit_id. They are matched to the units of 000021 and
000022 by the units table id. The check: the share of the tables' ids found, and agreement of the
values recomputed here with the tables' values, and the layer codes against the structures. The
timescale table was written by Figure3/timescale_calculation.py (response decay timescale over 40
to 290 ms after onset), whose procedure is recomputed. No script of the repository writes
time_to_first_spike.csv: Figure3/get_flash_PSTH.py reads it (as ttf0), computes first spikes from
30 ms within 200 ms, merges the two and writes a different file, time_to_first_spike_20200312.csv.
So the latency table's procedure is not known, and it is compared with every window computed here.
tfs_csv = repo_csv("time_to_first_spike.csv")
ts_csv = repo_csv("timescale_metrics.csv")
layer_csv = repo_csv("layer_info.csv")
ids = set(units.unit_id)
good_ids = set(units[units.quality == "good"].unit_id)
for name, t in (("time_to_first_spike", tfs_csv), ("timescale_metrics", ts_csv), ("layer_info", layer_csv)):
print(f"{name}: {len(t)} rows, {t.ecephys_unit_id.isin(ids).mean():.4f} of ids in the dandisets, "
f"{t.ecephys_unit_id.isin(good_ids).mean():.4f} among non-noise units")
print("non-noise units of the dandisets absent from layer_info:", len(good_ids - set(layer_csv.ecephys_unit_id)))
a = units.merge(tfs_csv, left_on="unit_id", right_on="ecephys_unit_id")
a = a[a.tfs_all_200.notna() & a.time_to_first_spike_fl.notna()]
d = (a.tfs_all_200 - a.time_to_first_spike_fl).abs() * 1000
print(f"\ntime to first spike, 200 ms window, {len(a)} units with both: within 1 ms {np.mean(d <= 1):.4f}, "
f"median |difference| {d.median():.2f} ms, 95th percentile {d.quantile(0.95):.2f} ms, "
f"median signed difference {((a.tfs_all_200 - a.time_to_first_spike_fl) * 1000).median():.2f} ms")
for c in ("tfs_all_250", "tfs_dark_250", "tfs_light_250", "tfs_all_500"):
dd = (a[c] - a.time_to_first_spike_fl).abs() * 1000
print(f" {c}: within 1 ms {np.mean(dd <= 1):.4f}, median |difference| {dd.median():.2f} ms")
b = units.merge(ts_csv, left_on="unit_id", right_on="ecephys_unit_id", suffixes=("", "_csv"))
b = b[b.timescale_ac_code.notna() & b.timescale_ac_csv.notna()]
rel = (b.timescale_ac_code - b.timescale_ac_csv).abs() / b.timescale_ac_csv.abs()
print(f"\nresponse decay timescale (40-290 ms), {len(b)} units with both: within 1% {np.mean(rel <= 0.01):.4f}")
print(f"units in timescale_metrics: {len(ts_csv)}; with a value here: {ts_csv.ecephys_unit_id.isin(units[units.timescale_ac_code.notna()].unit_id).sum()}")
bi = b[b.timescale_it_code.notna() & b.timescale_it_csv.notna()]
reli = (bi.timescale_it_code - bi.timescale_it_csv).abs() / bi.timescale_it_csv.abs()
print(f"intrinsic timescale (code window), {len(bi)} units: within 1% {np.mean(reli <= 0.01):.4f}")
c = units.merge(layer_csv, left_on="unit_id", right_on="ecephys_unit_id")
c["cortex"] = c.location.str.startswith("VIS")
print("\nlayer codes by structure class (rows: layer code; cortex = structure name starting with VIS):")
print(pd.crosstab(c.cortical_layer, c.cortex).to_string())
print("units with layer > 0 outside VIS* structures, by structure:",
c[(c.cortical_layer > 0) & ~c.cortex].location.value_counts().head(8).to_dict())
time_to_first_spike: 99180 rows, 1.0000 of ids in the dandisets, 1.0000 among non-noise units
timescale_metrics: 23524 rows, 1.0000 of ids in the dandisets, 1.0000 among non-noise units
layer_info: 99180 rows, 1.0000 of ids in the dandisets, 1.0000 among non-noise units
non-noise units of the dandisets absent from layer_info: 664
time to first spike, 200 ms window, 62483 units with both: within 1 ms 0.1613, median |difference| 4.00 ms, 95th percentile 37.00 ms, median signed difference -0.50 ms
tfs_all_250: within 1 ms 0.1125, median |difference| 8.00 ms
tfs_dark_250: within 1 ms 0.1362, median |difference| 12.00 ms
tfs_light_250: within 1 ms 0.1297, median |difference| 12.00 ms
tfs_all_500: within 1 ms 0.0299, median |difference| 64.00 ms
response decay timescale (40-290 ms), 23352 units with both: within 1% 0.9930
units in timescale_metrics: 23524; with a value here: 23352
intrinsic timescale (code window), 21122 units: within 1% 0.9955
layer codes by structure class (rows: layer code; cortex = structure name starting with VIS):
cortex False True
cortical_layer
0 44390 5746
1 11185 141
2 156 9247
3 3155 0
4 0 6470
5 0 13154
6 0 5536
units with layer > 0 outside VIS* structures, by structure: {'CA1': 11180, 'CA3': 3148, 'CA2': 156, 'DG': 10, 'HPF': 1, 'ProS': 1}
# One session, unit by unit: each flash's first spike from 30 ms, and the latency under other procedures
one = sessions_assets[sessions_assets.dandiset == "000021"].iloc[0]
f1 = S.open_asset(one.asset_id)
u1 = f1["units"]
ids1, ix1 = u1["id"][:], u1["spike_times_index"][:]
ix0 = np.r_[0, ix1[:-1]]
fl1 = S._intervals(f1, "flashes_presentations", ["color"])
fon1 = fl1.start_time.values
tab = tfs_csv.set_index("ecephys_unit_id").time_to_first_spike_fl
rows1 = np.where(pd.Series(ids1).isin(tab.dropna().index).values)[0]
st1 = u1["spike_times"][ix0[rows1.min()]:ix1[rows1.max()]]
T1 = []
for k in rows1:
st = st1[ix0[k] - ix0[rows1.min()]:ix1[k] - ix0[rows1.min()]]
j = np.searchsorted(st, fon1 + 0.030)
t = np.full(len(fon1), np.nan)
okk = j < len(st)
t[okk] = st[j[okk]] - fon1[okk]
T1.append(t)
f1.close()
T1 = np.array(T1)
tb1 = tab[ids1[rows1]].values
print(f"session {one.path}: {len(rows1)} units with a value in the paper's table, {len(fon1)} flashes")
grid = []
for end in (0.10, 0.12, 0.15, 0.17, 0.19, 0.20, 0.22, 0.25, 0.30, 0.50):
for trials, m in (("all", slice(None)), ("first 75", slice(0, 75)), ("last 75", slice(75, None))):
x = np.where(T1 < end, T1, np.nan)[:, m]
for stat, fn in (("median", np.nanmedian), ("mean", np.nanmean)):
with warnings.catch_warnings():
warnings.simplefilter("ignore")
v = fn(np.floor(x * 1000) / 1000, axis=1)
dd = np.abs(v - tb1) * 1000
grid.append({"end_ms": int(end * 1000), "flashes": trials, "statistic": stat,
"within_1ms": np.nanmean(dd <= 1), "median_abs_diff_ms": np.nanmedian(dd)})
grid = pd.DataFrame(grid).sort_values("within_1ms", ascending=False)
print("procedures tried, best five by share of units within 1 ms of the table:")
print(grid.head(5).round(3).to_string(index=False))
print("share within 1 ms, over all", len(grid), "procedures: %.3f to %.3f" % (grid.within_1ms.min(), grid.within_1ms.max()))
for k in range(3):
t = np.round(T1[k] * 1000, 1)
print(f"\nunit {ids1[rows1[k]]}: table {tb1[k] * 1000:.1f} ms; median within 200 ms {np.nanmedian(np.floor(t[t < 200])):.1f}, "
f"within 250 ms {np.nanmedian(np.floor(t[t < 250])):.1f}; trials with a spike before 200 ms {int((t < 200).sum())} of {len(t)}")
print(" first spike per flash from 30 ms (ms, in presentation order; None when there is none before 250 ms):",
[v if v < 250 else None for v in t.tolist()])
session sub-699733573/sub-699733573_ses-715093703.nwb: 1591 units with a value in the paper's table, 150 flashes
procedures tried, best five by share of units within 1 ms of the table:
end_ms flashes statistic within_1ms median_abs_diff_ms
200 all median 0.145 4.5
190 all median 0.137 5.0
220 all median 0.136 5.0
170 all median 0.111 7.0
200 last 75 median 0.101 7.0
share within 1 ms, over all 60 procedures: 0.011 to 0.145
unit 950922146: table 105.0 ms; median within 200 ms 102.0, within 250 ms 107.5; trials with a spike before 200 ms 34 of 150
first spike per flash from 30 ms (ms, in presentation order; None when there is none before 250 ms): [None, None, None, 51.2, None, None, None, None, None, None, None, None, 246.7, None, 220.9, None, None, None, None, None, None, 183.1, 144.4, None, 94.6, None, None, 110.8, 59.7, None, None, 57.9, 240.4, None, 191.2, None, None, None, None, 117.1, None, 96.7, 243.1, 99.7, 47.0, 39.3, None, 105.0, None, None, None, None, 114.7, None, 38.5, None, None, None, 105.4, None, None, None, None, None, None, None, 57.3, None, None, None, 199.6, None, None, 131.1, None, None, None, 31.0, None, None, None, None, None, None, None, None, None, 79.0, None, None, None, None, 174.4, None, None, None, 158.6, 105.5, None, 124.2, 119.0, None, None, None, None, None, None, None, 62.7, None, None, None, None, 226.3, 61.6, None, None, None, None, None, None, None, 156.7, None, None, None, None, None, None, None, 39.7, None, None, None, None, None, 43.8, 41.7, None, None, None, 123.0, None, None, None, 249.1, None, None, None, None]
unit 950922020: table 102.0 ms; median within 200 ms 104.0, within 250 ms 114.5; trials with a spike before 200 ms 62 of 150
first spike per flash from 30 ms (ms, in presentation order; None when there is none before 250 ms): [None, 221.2, None, None, 184.4, 54.8, 243.4, None, None, None, 104.1, 180.0, 216.3, None, 179.1, None, 148.8, 88.7, None, None, None, 129.9, None, 227.1, 88.8, 78.2, 99.3, 233.6, 210.0, 55.8, 77.3, None, 40.3, 118.1, 151.2, 116.4, None, None, 45.7, None, None, 57.7, None, 183.3, None, None, None, None, None, None, None, None, 129.9, None, 90.8, None, 241.9, 65.3, None, 53.6, None, None, None, None, None, None, 50.8, None, 43.1, None, 81.9, 127.7, None, None, None, 53.0, 34.9, 124.1, None, 137.0, 194.7, 139.8, None, 59.4, 42.0, 81.1, 45.4, None, None, 161.0, 134.8, 204.9, 49.5, 104.6, 139.9, None, None, None, 194.2, None, None, None, 83.4, None, None, None, 197.4, 40.8, None, 113.3, 34.8, None, 106.8, None, None, None, None, None, None, 138.7, None, 170.7, 132.8, None, 31.6, 142.9, None, None, 170.9, None, None, None, None, 56.1, 33.9, None, 108.9, None, 135.6, None, None, None, 153.4, None, None, 89.4, None, None, None, None]
unit 950921996: table 105.0 ms; median within 200 ms 109.0, within 250 ms 124.0; trials with a spike before 200 ms 54 of 150
first spike per flash from 30 ms (ms, in presentation order; None when there is none before 250 ms): [None, None, None, 88.5, 186.8, 209.3, None, 104.3, 76.0, None, None, None, 225.4, None, 122.1, 188.1, 48.1, 43.5, None, None, None, None, None, None, None, 239.5, None, None, 111.1, None, None, None, 69.9, None, None, None, None, None, None, None, 77.1, None, None, None, 149.5, 80.4, 46.1, 30.3, None, None, None, None, 30.4, 143.5, None, None, 126.4, None, None, 96.7, None, None, 179.8, None, None, None, None, 224.5, 175.7, None, None, 78.8, 51.6, None, 127.1, 172.9, None, 79.6, None, 186.4, None, None, 52.2, None, 128.0, 241.1, 160.9, 98.7, 185.9, None, None, None, 104.8, 147.9, 191.6, None, None, 126.8, 84.1, 83.5, None, None, None, None, None, None, None, 118.5, 241.8, 184.5, 46.6, 112.2, None, None, None, None, None, None, None, 107.7, None, 190.9, None, None, None, 60.1, None, None, 164.8, None, 240.1, 198.9, None, 170.3, 196.9, 176.1, None, None, None, 228.6, None, 85.8, 59.7, None, None, 97.1, None, None, 81.0, None]
Result: partially reproduced. Every id of the three tables is a non-noise unit of the dandisets (664 non-noise units are absent from the layer table). The timescales agree: 99.3% of 23,352 units within 1% for the response decay timescale and 99.6% of 21,122 for the intrinsic timescale. The time to first spike does not: with a 200 ms window 16.1% of 62,483 units are within 1 ms (median absolute difference 4.0 ms, median signed difference -0.5 ms), and the 250 ms, dark, light and 500 ms readings agree less (3.0% to 13.6%). On one session, unit by unit, none of 60 procedures (window ends from 100 to 500 ms, all flashes or one half, median or mean) brings more than 14.5% of 1,591 units within 1 ms of the table, so the table's procedure is not found; the script that would show it is not in the repository. The layer codes 1, 2 and 3 are also given to CA1, CA2 and CA3 units, so code 2 (L2/3) is not cortical only; with cortical_depth > 0, which the layer analysis requires, only VIS* units remain (listed under ed10b).
Figure 3: unit selection and the four metrics¶
The selection of Figure3.py: per area, p_value_rf < 0.01, area_rf < 2500, snr > 1, firing_rate_dg > 0.1 Hz; then flash time to first spike < 100 ms for all four metrics; then, for the response decay timescale, 1 < timescale < 300 ms, more than 50 spikes and SD of the estimate < 20. Area means are over units (log10 for the modulation index), and Pearson's r is against the eight hierarchy scores. Pairs of areas: two-sided rank-sum with 1e-5 added to each p, Benjamini-Hochberg over the 28 pairs at 0.05, as Figure3.py.
HS = np.array([S.HIERARCHY[a] for a in S.AREAS])
U3 = units[units.location.isin(S.AREAS) & (units.quality == "good")].copy()
base = (U3.p_value_rf < 0.01) & (U3.area_rf < 2500) & (U3.snr > 1) & (U3.firing_rate_dg > 0.1)
def area_corr(values, sel, f=lambda v: v, label=""):
means = pd.Series({a: np.nanmean(f(values[sel & (U3.location == a)].values)) for a in S.AREAS})
n = pd.Series({a: int((sel & (U3.location == a) & values.notna()).sum()) for a in S.AREAS})
r, p = stats.pearsonr(HS, means.values)
return means, n, r, p
def pairs_test(values, sel, f=lambda v: v):
vals = {a: f(values[sel & (U3.location == a)].dropna().values) for a in S.AREAS}
rows = []
for i, a1 in enumerate(S.AREAS):
for a2 in S.AREAS[i + 1:]:
z, p = stats.ranksums(vals[a1], vals[a2])
rows.append({"pair": f"{S.NAMES[a1]}-{S.NAMES[a2]}", "p": p + 1e-5})
t = pd.DataFrame(rows)
t["reject"], t["p_bh"] = multipletests(t.p.values, alpha=0.05, method="fdr_bh")[:2]
return t
def judge_pairs(t, expected_ns):
ns = set(t.pair[~t.reject])
exp = {p if p in set(t.pair) else "-".join(p.split("-")[::-1]) for p in expected_ns}
disagree = (ns - exp) | (exp - ns)
return sorted(ns), sorted(disagree)
print("units per area: all non-noise, with RF, after the filters")
cnt = pd.DataFrame({"non_noise": U3.groupby("location").size(),
"p_rf<0.01": U3[U3.p_value_rf < 0.01].groupby("location").size(),
"after_filter": U3[base].groupby("location").size()}).reindex(S.AREAS)
cnt.index = [S.NAMES[a] for a in cnt.index]
print(cnt.to_string())
other3 = (U3.p_value_rf < 0.01) & (U3.area_rf < 2500) & (U3.snr > 1)
print("units passing the other three filters:", int(other3.sum()), "; of these removed by firing_rate_dg > 0.1 Hz:",
int((other3 & ~(U3.firing_rate_dg > 0.1)).sum()))
print("mice per area after the filter:", {S.NAMES[a]: U3[base & (U3.location == a)].subject.nunique() for a in S.AREAS})
units per area: all non-noise, with RF, after the filters
non_noise p_rf<0.01 after_filter
LGN 2591 1352 1233
V1 8605 5922 5581
LM 4941 3486 2911
RL 6031 3465 3122
LP 4873 2197 2073
AL 6468 4103 3755
PM 4270 2400 2229
AM 6211 3769 3555
units passing the other three filters: 25019 ; of these removed by firing_rate_dg > 0.1 Hz: 560
mice per area after the filter: {'LGN': 32, 'V1': 55, 'LM': 42, 'RL': 50, 'LP': 40, 'AL': 44, 'PM': 35, 'AM': 50}
fig3c-latency-hierarchy¶
The mean time to first spike after a full-field flash in each area is correlated with its anatomical hierarchy score.
Readings of the open detail (which flashes): (a) all, (b) dark, (c) light; window end 250 ms (Methods), with 500 ms (functions/time_to_first_spike.py default) as the claims file's variant, and 200 ms (the window of Figure3/get_flash_PSTH.py, which does not write the paper's table) and the paper's table as further variants. Each reading's own latency is used for the < 100 ms cut.
rows = []
for col in ["tfs_all_250", "tfs_dark_250", "tfs_light_250", "tfs_all_500", "tfs_dark_500", "tfs_light_500", "tfs_all_200"]:
sel = base & (U3[col] < 0.1)
means, n, r, p = area_corr(U3[col] * 1000, sel)
rows.append({"reading": col, "r": r, "p": p, "units": int(n.sum()), **{S.NAMES[a]: means[a] for a in S.AREAS}})
U3 = U3.merge(tfs_csv, left_on="unit_id", right_on="ecephys_unit_id", how="left").set_index(U3.index)
sel = base & (U3.time_to_first_spike_fl < 0.1)
means, n, r, p = area_corr(U3.time_to_first_spike_fl * 1000, sel)
rows.append({"reading": "paper table time_to_first_spike.csv", "r": r, "p": p, "units": int(n.sum()), **{S.NAMES[a]: means[a] for a in S.AREAS}})
lat = pd.DataFrame(rows).set_index("reading")
print("mean time to first spike per area (ms), Pearson r and p against hierarchy score (expected r 0.95 +/- 0.095):")
print(lat.round(4).to_string())
lat_sel = base & (U3.tfs_all_250 < 0.1)
print("\nper area range of unit latencies (ms), reading (a):",
{S.NAMES[a]: (round(U3.tfs_all_250[lat_sel & (U3.location == a)].min() * 1000), round(U3.tfs_all_250[lat_sel & (U3.location == a)].max() * 1000)) for a in S.AREAS})
by_set = []
for stp in ("brain_observatory_1.1", "functional_connectivity"):
m_, n_, r_, p_ = area_corr(U3.tfs_all_250 * 1000, lat_sel & (U3.session_type == stp))
by_set.append({"set": stp, "r": r_, "p": p_, "units": int(n_.sum())})
print("reading (a) by stimulus set:", pd.DataFrame(by_set).round(4).to_dict("records"))
mean time to first spike per area (ms), Pearson r and p against hierarchy score (expected r 0.95 +/- 0.095):
r p units LGN V1 LM RL LP AL PM AM
reading
tfs_all_250 0.9591 0.0002 13344 66.4089 66.0935 69.2344 71.4664 71.5003 71.1346 75.1614 75.0562
tfs_dark_250 0.9044 0.0020 11708 62.9946 61.9246 63.9122 67.9350 67.6542 66.5561 72.6798 71.9158
tfs_light_250 0.9329 0.0007 14635 66.4857 65.6572 71.0752 72.0624 72.3894 72.5477 73.5040 74.0889
tfs_all_500 0.8021 0.0166 6527 70.2665 67.2058 70.7875 70.4127 72.4994 71.5684 75.0506 72.9228
tfs_dark_500 0.9548 0.0002 6504 62.9072 63.6804 66.1338 67.8960 68.8208 67.5137 72.6911 71.6830
tfs_light_500 0.9465 0.0004 7547 65.7989 67.7294 71.7007 70.5670 73.7807 72.5449 73.8722 73.9128
tfs_all_200 0.9364 0.0006 17163 64.4900 64.4537 66.9964 70.0583 68.6017 68.9797 73.2659 73.4453
paper table time_to_first_spike.csv 0.9531 0.0002 15849 62.0256 63.3343 66.1658 69.0097 67.9994 67.9524 72.8890 72.3628
per area range of unit latencies (ms), reading (a): {'LGN': (34, 100), 'V1': (30, 100), 'LM': (30, 100), 'RL': (31, 100), 'LP': (37, 100), 'AL': (30, 100), 'PM': (35, 100), 'AM': (32, 100)}
reading (a) by stimulus set: [{'set': 'brain_observatory_1.1', 'r': 0.9614, 'p': 0.0001, 'units': 7082}, {'set': 'functional_connectivity', 'r': 0.891, 'p': 0.003, 'units': 6262}]
fig, ax = plt.subplots(1, 1, figsize=(4, 3))
ax.plot(HS, lat.loc["tfs_all_250", [S.NAMES[a] for a in S.AREAS]].values.astype(float), "o")
for a in S.AREAS:
ax.annotate(S.NAMES[a], (S.HIERARCHY[a], lat.loc["tfs_all_250", S.NAMES[a]]), fontsize=8)
ax.set_xlabel("anatomical hierarchy score")
ax.set_ylabel("mean time to first spike (ms)")
plt.tight_layout()
print("plotted range (ms):", lat.loc["tfs_all_250", [S.NAMES[a] for a in S.AREAS]].astype(float).agg(["min", "max"]).round(2).to_dict())
plotted range (ms): {'min': 66.09, 'max': 75.16}
Result: reproduced. With all flashes and the Methods' 250 ms window (13,344 units), the area means range from 66.1 ms (V1) and 66.4 ms (LGN) to 75.2 ms (PM) and 75.1 ms (AM), and r = 0.959, p = 0.0002 (paper 0.95). Dark flashes give r = 0.904 and light flashes 0.933, both within tolerance with p < 0.01, so the readings agree. The paper's own latency table gives r = 0.953.
ed9a-latency-pairs¶
Latencies differ significantly between all pairs of areas except LGN-V1, RL-LP, LP-AL and AM-PM.
t = pairs_test(U3.tfs_all_250 * 1000, lat_sel)
ns, dis = judge_pairs(t, ["LGN-V1", "RL-LP", "LP-AL", "PM-AM"])
print(t.round(6).to_string(index=False))
print("\nnot significant:", ns)
print("pairs that disagree with the claim:", len(dis), dis)
for col in ("tfs_dark_250", "tfs_light_250", "tfs_all_500"):
tt = pairs_test(U3[col] * 1000, base & (U3[col] < 0.1))
print(f"variant {col}: not significant {judge_pairs(tt, [])[0]}")
pair p reject p_bh LGN-V1 0.534526 False 0.575643 LGN-LM 0.000036 True 0.000051 LGN-RL 0.000010 True 0.000015 LGN-LP 0.000010 True 0.000015 LGN-AL 0.000010 True 0.000015 LGN-PM 0.000010 True 0.000015 LGN-AM 0.000010 True 0.000015 V1-LM 0.000010 True 0.000015 V1-RL 0.000010 True 0.000015 V1-LP 0.000010 True 0.000015 V1-AL 0.000010 True 0.000015 V1-PM 0.000010 True 0.000015 V1-AM 0.000010 True 0.000015 LM-RL 0.000132 True 0.000168 LM-LP 0.000104 True 0.000139 LM-AL 0.000712 True 0.000867 LM-PM 0.000010 True 0.000015 LM-AM 0.000010 True 0.000015 RL-LP 0.866404 False 0.898493 RL-AL 0.467074 False 0.523123 RL-PM 0.000010 True 0.000015 RL-AM 0.000010 True 0.000015 LP-AL 0.380376 False 0.443772 LP-PM 0.000010 True 0.000015 LP-AM 0.000010 True 0.000015 AL-PM 0.000010 True 0.000015 AL-AM 0.000010 True 0.000015 PM-AM 0.901965 False 0.901965 not significant: ['LGN-V1', 'LP-AL', 'PM-AM', 'RL-AL', 'RL-LP'] pairs that disagree with the claim: 1 ['RL-AL'] variant tfs_dark_250: not significant ['LGN-LM', 'LGN-V1', 'LP-AL', 'PM-AM', 'RL-LP'] variant tfs_light_250: not significant ['AL-PM', 'LGN-V1', 'LM-LP', 'LM-RL', 'LP-AL', 'PM-AM', 'RL-AL', 'RL-LP'] variant tfs_all_500: not significant ['AL-AM', 'LGN-AL', 'LGN-LM', 'LGN-RL', 'LM-AL', 'LM-LP', 'LM-RL', 'LP-AL', 'LP-AM', 'RL-AL']
Result: partially reproduced. Of the 28 pairs, 27 agree with the claim: LGN-V1, RL-LP, LP-AL and PM-AM are not significant, but RL-AL is not significant either (corrected p = 0.52). With dark flashes LGN-LM is the one pair that disagrees.
fig3f-rf-size-hierarchy¶
Receptive field size increases with anatomical hierarchy score.
means, n, r, p = area_corr(U3.area_rf, lat_sel)
print("mean RF area (deg2) per area:", {S.NAMES[a]: round(means[a], 1) for a in S.AREAS})
print("units:", {S.NAMES[a]: int(n[a]) for a in S.AREAS})
print(f"Pearson r = {r:.4f}, p = {p:.3g} (expected 0.97 +/- 0.097)")
m2, n2, r2, p2 = area_corr(U3.area_rf, base)
print(f"variant without the latency cut: r = {r2:.4f}, p = {p2:.3g}")
for stp in ("brain_observatory_1.1", "functional_connectivity"):
m_, n_, r_, p_ = area_corr(U3.area_rf, lat_sel & (U3.session_type == stp))
print(f"{stp}: r = {r_:.4f}, p = {p_:.3g}")
print("range of unit RF areas (deg2):", U3.area_rf[lat_sel].min(), "to", U3.area_rf[lat_sel].max())
mean RF area (deg2) per area: {'LGN': np.float64(543.9), 'V1': np.float64(579.0), 'LM': np.float64(732.8), 'RL': np.float64(710.6), 'LP': np.float64(688.2), 'AL': np.float64(754.0), 'PM': np.float64(846.3), 'AM': np.float64(905.1)}
units: {'LGN': 851, 'V1': 2385, 'LM': 1630, 'RL': 1741, 'LP': 1588, 'AL': 2254, 'PM': 1106, 'AM': 1789}
Pearson r = 0.9577, p = 0.000183 (expected 0.97 +/- 0.097)
variant without the latency cut: r = 0.9737, p = 4.47e-05
brain_observatory_1.1: r = 0.9570, p = 0.000193
functional_connectivity: r = 0.9125, p = 0.00157
range of unit RF areas (deg2): 100.0 to 2400.0
Result: reproduced. Mean RF area rises from 544 deg2 (LGN) to 905 deg2 (AM); r = 0.958, p = 0.0002 (paper 0.97).
ed9d-rf-size-pairs¶
Receptive field sizes differ significantly between all pairs of areas except LM-RL.
t = pairs_test(U3.area_rf, lat_sel)
ns, dis = judge_pairs(t, ["LM-RL"])
print(t.round(6).to_string(index=False))
print("\nnot significant:", ns)
print("pairs that disagree with the claim:", len(dis), dis)
pair p reject p_bh LGN-V1 0.000010 True 0.000013 LGN-LM 0.000010 True 0.000013 LGN-RL 0.000010 True 0.000013 LGN-LP 0.000010 True 0.000013 LGN-AL 0.000010 True 0.000013 LGN-PM 0.000010 True 0.000013 LGN-AM 0.000010 True 0.000013 V1-LM 0.000010 True 0.000013 V1-RL 0.000010 True 0.000013 V1-LP 0.000309 True 0.000360 V1-AL 0.000010 True 0.000013 V1-PM 0.000010 True 0.000013 V1-AM 0.000010 True 0.000013 LM-RL 0.408960 False 0.408960 LM-LP 0.000011 True 0.000013 LM-AL 0.363684 False 0.377154 LM-PM 0.000010 True 0.000013 LM-AM 0.000010 True 0.000013 RL-LP 0.000040 True 0.000049 RL-AL 0.062632 False 0.067450 RL-PM 0.000010 True 0.000013 RL-AM 0.000010 True 0.000013 LP-AL 0.000010 True 0.000013 LP-PM 0.000010 True 0.000013 LP-AM 0.000010 True 0.000013 AL-PM 0.000010 True 0.000013 AL-AM 0.000010 True 0.000013 PM-AM 0.001088 True 0.001219 not significant: ['LM-AL', 'LM-RL', 'RL-AL'] pairs that disagree with the claim: 2 ['LM-AL', 'RL-AL']
Result: partially reproduced. LM-RL is not significant as stated, and two further pairs, LM-AL (corrected p = 0.38) and RL-AL (0.067), are not significant either; 26 of 28 pairs agree.
fig3i-modulation-index-hierarchy¶
The modulation index of responses to drifting gratings is highest in LGN and decreases along the hierarchy.
log10 = lambda v: np.log10(v)
means, n, r, p = area_corr(U3.mod_idx_dg, lat_sel, log10)
print("mean log10 MI per area:", {S.NAMES[a]: round(means[a], 3) for a in S.AREAS})
print("units:", {S.NAMES[a]: int(n[a]) for a in S.AREAS})
print(f"Pearson r = {r:.4f}, p = {p:.3g} (expected -0.89 +/- 0.089); highest area: {S.NAMES[means.idxmax()]}")
m2, n2, r2, p2 = area_corr(U3.mod_idx_dg, lat_sel)
print(f"variant, mean of MI: r = {r2:.4f}, p = {p2:.3g}; highest {S.NAMES[m2.idxmax()]}")
for stp in ("brain_observatory_1.1", "functional_connectivity"):
m_, n_, r_, p_ = area_corr(U3.mod_idx_dg, lat_sel & (U3.session_type == stp), log10)
print(f"{stp}: r = {r_:.4f}, p = {p_:.3g}, highest {S.NAMES[m_.idxmax()]}")
print("units with MI not computed or zero:", int((U3.mod_idx_dg[lat_sel].isna() | (U3.mod_idx_dg[lat_sel] <= 0)).sum()))
print("range of MI:", U3.mod_idx_dg[lat_sel].min(), "to", U3.mod_idx_dg[lat_sel].max())
mean log10 MI per area: {'LGN': np.float64(0.628), 'V1': np.float64(0.41), 'LM': np.float64(0.09), 'RL': np.float64(0.083), 'LP': np.float64(0.091), 'AL': np.float64(0.009), 'PM': np.float64(-0.079), 'AM': np.float64(-0.124)}
units: {'LGN': 851, 'V1': 2385, 'LM': 1630, 'RL': 1741, 'LP': 1588, 'AL': 2254, 'PM': 1106, 'AM': 1789}
Pearson r = -0.9535, p = 0.000243 (expected -0.89 +/- 0.089); highest area: LGN
variant, mean of MI: r = -0.9343, p = 0.000674; highest LGN
brain_observatory_1.1: r = -0.9106, p = 0.00167, highest LGN
functional_connectivity: r = -0.9631, p = 0.000123, highest LGN units with MI not computed or zero: 0 range of MI: 7.032151919312814e-05 to 9.71244661568323
Result: reproduced. Mean log10 MI is highest in LGN (0.628) and lowest in AM (-0.124); r = -0.954, p = 0.0002 (paper -0.89).
ed9e-modulation-index-pairs¶
Modulation indices differ significantly between all pairs of areas except RL-AL and AM-PM.
t = pairs_test(U3.mod_idx_dg.where(U3.mod_idx_dg > 0), lat_sel, log10)
ns, dis = judge_pairs(t, ["RL-AL", "PM-AM"])
print(t.round(6).to_string(index=False))
print("\nnot significant:", ns)
print("pairs that disagree with the claim:", len(dis), dis)
pair p reject p_bh LGN-V1 0.000010 True 0.000014 LGN-LM 0.000010 True 0.000014 LGN-RL 0.000010 True 0.000014 LGN-LP 0.000010 True 0.000014 LGN-AL 0.000010 True 0.000014 LGN-PM 0.000010 True 0.000014 LGN-AM 0.000010 True 0.000014 V1-LM 0.000010 True 0.000014 V1-RL 0.000010 True 0.000014 V1-LP 0.000010 True 0.000014 V1-AL 0.000010 True 0.000014 V1-PM 0.000010 True 0.000014 V1-AM 0.000010 True 0.000014 LM-RL 0.607602 False 0.630106 LM-LP 0.878786 False 0.878786 LM-AL 0.000102 True 0.000124 LM-PM 0.000010 True 0.000014 LM-AM 0.000010 True 0.000014 RL-LP 0.498849 False 0.537221 RL-AL 0.000497 True 0.000580 RL-PM 0.000010 True 0.000014 RL-AM 0.000010 True 0.000014 LP-AL 0.000073 True 0.000093 LP-PM 0.000010 True 0.000014 LP-AM 0.000010 True 0.000014 AL-PM 0.000011 True 0.000014 AL-AM 0.000010 True 0.000014 PM-AM 0.107435 False 0.120327 not significant: ['LM-LP', 'LM-RL', 'PM-AM', 'RL-LP'] pairs that disagree with the claim: 4 ['LM-LP', 'LM-RL', 'RL-AL', 'RL-LP']
Result: not reproduced. PM-AM is not significant as stated, but RL-AL is significant (corrected p = 0.0006), and LM-RL, LM-LP and RL-LP are not; 4 of 28 pairs disagree, more than the two the criterion allows for a partial result.
ed9fg-intrinsic-timescale-not-correlated¶
The mean intrinsic timescale of spontaneous activity in each area is not correlated with the visual hierarchy.
Units of timescale_calculation.py (the Figure 3 filter, no latency cut); kept when the SD of the estimate < 100 and at least 100 spikes. Variant from the code: 1.0 to 2.0 s after onset and 5 < timescale < 400 ms.
Result: reproduced. The mean intrinsic timescale is not correlated with the hierarchy: r = -0.158, p = 0.71 (paper -0.24, 0.57). The code's window gives r = -0.290, p = 0.49. The test is weak: 8,834 of the 15,958 kept units sit at the 1 ms lower bound of the fit (median 1.00 ms), so most fits under the Methods' reading are degenerate; without them r = 0.082, p = 0.85.
fig3l-response-decay-timescale-hierarchy¶
Higher-order areas have a longer response decay timescale during the flash response, and the area mean correlates with hierarchy score.
def ac_sel(col, err, cnt):
return lat_sel & (U3[col] < 300) & (U3[col] > 1) & (U3[cnt] > 50) & (U3[err] < 20)
ac_main = ac_sel("timescale_ac", "err_ac", "spike_count_ac")
means, n, r, p = area_corr(U3.timescale_ac, ac_main)
print("mean response decay timescale (ms), 0-250 ms:", {S.NAMES[a]: round(means[a], 1) for a in S.AREAS})
print("units:", {S.NAMES[a]: int(n[a]) for a in S.AREAS}, " total", int(n.sum()))
print(f"Pearson r = {r:.4f}, p = {p:.3g} (expected 0.86 +/- 0.086)")
ac_code = ac_sel("timescale_ac_code", "err_ac_code", "spike_count_ac_code")
m2, n2, r2, p2 = area_corr(U3.timescale_ac_code, ac_code)
print("variant 40-290 ms (code):", {S.NAMES[a]: round(m2[a], 1) for a in S.AREAS})
print(f" r = {r2:.4f}, p = {p2:.3g}, units {int(n2.sum())}")
U3 = U3.merge(ts_csv[["ecephys_unit_id", "timescale_ac", "err_ac", "spike_count_ac"]].rename(columns=lambda c: c + "_csv" if c != "ecephys_unit_id" else "tsid"),
left_on="unit_id", right_on="tsid", how="left").set_index(U3.index)
csv_sel = (base & (U3.time_to_first_spike_fl < 0.1) & (U3.timescale_ac_csv < 300) & (U3.timescale_ac_csv > 1)
& (U3.spike_count_ac_csv > 50) & (U3.err_ac_csv < 20))
m3, n3, r3, p3 = area_corr(U3.timescale_ac_csv, csv_sel)
print(f"variant, the paper's tables for timescale and latency: r = {r3:.4f}, p = {p3:.3g}, units {int(n3.sum())}")
mean response decay timescale (ms), 0-250 ms: {'LGN': np.float64(44.6), 'V1': np.float64(41.4), 'LM': np.float64(43.0), 'RL': np.float64(49.5), 'LP': np.float64(47.1), 'AL': np.float64(45.6), 'PM': np.float64(53.3), 'AM': np.float64(54.6)}
units: {'LGN': 513, 'V1': 1061, 'LM': 780, 'RL': 867, 'LP': 730, 'AL': 1119, 'PM': 301, 'AM': 471} total 5842
Pearson r = 0.8038, p = 0.0162 (expected 0.86 +/- 0.086)
variant 40-290 ms (code): {'LGN': np.float64(35.8), 'V1': np.float64(38.6), 'LM': np.float64(40.0), 'RL': np.float64(44.9), 'LP': np.float64(42.1), 'AL': np.float64(40.6), 'PM': np.float64(49.9), 'AM': np.float64(49.3)}
r = 0.8763, p = 0.00431, units 7467
variant, the paper's tables for timescale and latency: r = 0.8548, p = 0.00685, units 7798
Result: reproduced. With the Methods' 0 to 250 ms window, r = 0.804, p = 0.016 (paper 0.86; within tolerance), with PM and AM highest (53.3 and 54.6 ms). The code's 40 to 290 ms window gives r = 0.876, and the paper's table 0.855.
ed9h-response-decay-pairs¶
Response decay timescales differ significantly between all pairs of areas except LM-AL and AM-PM.
t = pairs_test(U3.timescale_ac, ac_main)
ns, dis = judge_pairs(t, ["LM-AL", "PM-AM"])
print(t.round(6).to_string(index=False))
print("\nnot significant:", ns)
print("pairs that disagree with the claim:", len(dis), dis)
tc = pairs_test(U3.timescale_ac_code, ac_code)
print("variant 40-290 ms, not significant:", judge_pairs(tc, [])[0])
pair p reject p_bh LGN-V1 0.000364 True 0.000537 LGN-LM 0.062979 False 0.067823 LGN-RL 0.000054 True 0.000088 LGN-LP 0.015998 True 0.019476 LGN-AL 0.814266 False 0.814266 LGN-PM 0.000010 True 0.000024 LGN-AM 0.000010 True 0.000024 V1-LM 0.042538 True 0.049627 V1-RL 0.000010 True 0.000024 V1-LP 0.000010 True 0.000024 V1-AL 0.000012 True 0.000025 V1-PM 0.000010 True 0.000024 V1-AM 0.000010 True 0.000024 LM-RL 0.000010 True 0.000024 LM-LP 0.000011 True 0.000025 LM-AL 0.013973 True 0.018180 LM-PM 0.000010 True 0.000024 LM-AM 0.000010 True 0.000024 RL-LP 0.057244 False 0.064113 RL-AL 0.000014 True 0.000026 RL-PM 0.014285 True 0.018180 RL-AM 0.000037 True 0.000064 LP-AL 0.008236 True 0.011530 LP-PM 0.000133 True 0.000207 LP-AM 0.000010 True 0.000024 AL-PM 0.000010 True 0.000024 AL-AM 0.000010 True 0.000024 PM-AM 0.267709 False 0.277625 not significant: ['LGN-AL', 'LGN-LM', 'PM-AM', 'RL-LP'] pairs that disagree with the claim: 4 ['LGN-AL', 'LGN-LM', 'LM-AL', 'RL-LP']
variant 40-290 ms, not significant: ['LGN-V1', 'LM-AL', 'PM-AM']
Result: not reproduced. With the Methods' window PM-AM is not significant as stated, but LM-AL is significant (corrected p = 0.018), and LGN-LM, LGN-AL and RL-LP are not; 4 of 28 pairs disagree. With the code's 40 to 290 ms window only LGN-V1 disagrees.
fig3-unit-counts¶
The Figure 3 analyses use 15,713 units from 58 mice, and the response decay timescale 7,837 units.
print("judged, the claims' procedure (latency within 250 ms; decay timescale over 0 to 250 ms):")
print(" readings of the 15,713: (a) after the latency cut", int(lat_sel.sum()), "; (b) before it", int(base.sum()),
" (expected 15,713 +/- 1,571)")
print(" mice in (a):", U3[lat_sel].subject.nunique())
print(" response decay timescale units:", int(ac_main.sum()), " (expected 7,837 +/- 784)")
print("variants, not judged:")
print(" without the firing_rate_dg > 0.1 Hz cut (its definition is assumed): after the latency cut",
int((other3 & (U3.tfs_all_250 < 0.1)).sum()), "; before it", int(other3.sum()))
print(" latency within 500 ms (the claims file's variant): after the latency cut", int((base & (U3.tfs_all_500 < 0.1)).sum()))
print(" decay timescale over 40 to 290 ms with the 250 ms latency cut:", int(ac_code.sum()))
code_lat = base & (U3.tfs_all_200 < 0.1)
code_ac = (code_lat & (U3.timescale_ac_code < 300) & (U3.timescale_ac_code > 1) & (U3.spike_count_ac_code > 50)
& (U3.err_ac_code < 20))
print(" latency within 200 ms (get_flash_PSTH.py) and decay timescale over 40 to 290 ms (timescale_calculation.py):",
"after the latency cut", int(code_lat.sum()), "; decay timescale", int(code_ac.sum()))
print(" the paper's tables for latency and timescale:", int(csv_sel.sum()),
"; after the latency cut from the paper's table:", int((base & (U3.time_to_first_spike_fl < 0.1)).sum()))
judged, the claims' procedure (latency within 250 ms; decay timescale over 0 to 250 ms): readings of the 15,713: (a) after the latency cut 13344 ; (b) before it 24459 (expected 15,713 +/- 1,571) mice in (a): 58 response decay timescale units: 5842 (expected 7,837 +/- 784) variants, not judged: without the firing_rate_dg > 0.1 Hz cut (its definition is assumed): after the latency cut 13476 ; before it 25019 latency within 500 ms (the claims file's variant): after the latency cut 6527 decay timescale over 40 to 290 ms with the 250 ms latency cut: 7467 latency within 200 ms (get_flash_PSTH.py) and decay timescale over 40 to 290 ms (timescale_calculation.py): after the latency cut 17163 ; decay timescale 8381 the paper's tables for latency and timescale: 7798 ; after the latency cut from the paper's table: 15849
Result: not reproduced. With the claims' procedure (latency within the Methods' 250 ms; response decay timescale over the Methods' 0 to 250 ms) the latency cut keeps 13,344 units (reading a) and the unit table before it 24,459 (reading b), against 15,713 +/- 1,571, and the decay timescale 5,842 against 7,837 +/- 784: no count is within tolerance. Variants, not judged: with a 200 ms latency window and the 40 to 290 ms decay window of the paper's scripts, 17,163 and 8,381, and with the paper's own tables 15,849 and 7,798, all within tolerance. Without the firing_rate_dg > 0.1 Hz cut, whose definition is assumed, (a) is 13,476 and (b) 25,019: the cut removes 560 units and does not change the outcome. The first revision adopted the scripts' windows on the reading that get_flash_PSTH.py wrote the paper's latency table; it reads that table and does not write it (link-unit-tables), so that basis is gone and the claims file's reading is kept.
ed9bc-spontaneous-latency¶
Time to first spike measured in the grey pre-stimulus interval is not correlated with hierarchy score but is strongly negatively correlated with baseline firing rate.
The same units as fig3c (reading (a)). Baseline rate readings: (a) the mean rate in the same 1 s grey interval, (b) the unit table's overall rate.
means, n, r, p = area_corr(U3.tfs_spont * 1000, lat_sel)
print("mean spontaneous time to first spike (ms):", {S.NAMES[a]: round(means[a], 1) for a in S.AREAS})
print(f"against hierarchy: r = {r:.4f}, p = {p:.3g} (paper r = 0.57, p = 0.14)")
rate_a = pd.Series({a: U3.rate_pre[lat_sel & (U3.location == a)].mean() for a in S.AREAS})
rate_b = pd.Series({a: U3.firing_rate[lat_sel & (U3.location == a)].mean() for a in S.AREAS})
for lab, rate in (("(a) pre-flash rate", rate_a), ("(b) overall rate", rate_b)):
rr, pp = stats.pearsonr(rate.values, means.values)
print(f"against baseline rate {lab}: r = {rr:.4f}, p = {pp:.3g} (expected -0.98 +/- 0.098); rates {rate.round(2).to_dict()}")
print("units without a spontaneous latency:", int(U3.tfs_spont[lat_sel].isna().sum()), "of", int(lat_sel.sum()))
mean spontaneous time to first spike (ms): {'LGN': np.float64(94.9), 'V1': np.float64(109.7), 'LM': np.float64(111.8), 'RL': np.float64(105.0), 'LP': np.float64(104.1), 'AL': np.float64(110.2), 'PM': np.float64(107.0), 'AM': np.float64(103.4)}
against hierarchy: r = 0.2903, p = 0.485 (paper r = 0.57, p = 0.14)
against baseline rate (a) pre-flash rate: r = -0.9904, p = 2.22e-06 (expected -0.98 +/- 0.098); rates {'LGd': 9.32, 'VISp': 5.16, 'VISl': 4.37, 'VISrl': 6.24, 'LP': 6.59, 'VISal': 4.4, 'VISpm': 5.93, 'VISam': 6.44}
against baseline rate (b) overall rate: r = -0.9862, p = 6.43e-06 (expected -0.98 +/- 0.098); rates {'LGd': 10.7, 'VISp': 6.97, 'VISl': 6.02, 'VISrl': 7.8, 'LP': 7.7, 'VISal': 6.27, 'VISpm': 7.52, 'VISam': 8.18}
units without a spontaneous latency: 156 of 13344
Result: reproduced. The spontaneous time to first spike is not correlated with the hierarchy (r = 0.290, p = 0.49; paper 0.57, 0.14) and is correlated with the pre-flash rate (r = -0.991, p = 2e-06; paper -0.98). The overall rate gives r = -0.987.
ed9ij-firing-rate-not-correlated¶
Overall firing rate does not correlate with hierarchy score.
Result: reproduced. Mean log10 firing rate is not correlated with the hierarchy: r = -0.327, p = 0.43; the untransformed mean gives r = -0.363, p = 0.38.
ed10b-layers¶
Compared with deep layers (5 and 6), superficial layers (2/3 and 4) have an earlier time to first spike, smaller receptive fields, a higher modulation index and a faster response decay timescale.
As the significance block of ExtDataFigure10.py, with the paper's layer table (linked above): units with cortical_depth > 0 and the Figure 3 filters, cumulative cuts; superficial layers 2 and 4 (the codes the figure names L2/3 and L4), deep 5 and 6; difference of means; rank-sum.
L = U3.merge(layer_csv, left_on="unit_id", right_on="ecephys_unit_id", how="left", suffixes=("", "_l")).set_index(U3.index)
Lall = units[(units.quality == "good")].merge(layer_csv, left_on="unit_id", right_on="ecephys_unit_id")
print("layer codes among non-noise units with depth > 0:", Lall[Lall.cortical_depth > 0].cortical_layer.value_counts().sort_index().to_dict())
# all cortical units with depth: the pool of the code is all structures, so the Figure 3 frame (eight areas) is widened
Ucx = units[(units.quality == "good") & units.p_value_rf.notna()].merge(layer_csv, left_on="unit_id", right_on="ecephys_unit_id")
print("units measured here with depth > 0:", int((Ucx.cortical_depth > 0).sum()),
"; their structures:", Ucx[Ucx.cortical_depth > 0].location.value_counts().to_dict())
def layer_test(D, sup, deep, tfs="tfs_all_250", ts=("timescale_ac", "err_ac", "spike_count_ac"), mi=np.log10):
sel = (D.p_value_rf < 0.01) & (D.cortical_depth > 0) & (D.area_rf < 2500) & (D.snr > 1) & (D.firing_rate_dg > 0.1)
out = []
specs = [("time to first spike (ms)", tfs, lambda v: v * 1000, "deep-sup"),
("RF area (deg2)", "area_rf", lambda v: v, "deep-sup"),
("log10 MI" if mi is np.log10 else "MI", "mod_idx_dg", mi, "sup-deep"),
("response decay timescale (ms)", ts[0], lambda v: v, "deep-sup")]
for k, (name, col, f, sign) in enumerate(specs):
if k == 0:
sel = sel & (D[tfs] < 0.1)
if k == 3:
sel = sel & (D[ts[0]] < 300) & (D[ts[0]] > 1) & (D[ts[2]] > 50) & (D[ts[1]] < 20)
a = f(D[col][sel & D.cortical_layer.isin(sup)].values)
b = f(D[col][sel & D.cortical_layer.isin(deep)].values)
a, b = a[np.isfinite(a)], b[np.isfinite(b)]
diff = (b.mean() - a.mean()) if sign == "deep-sup" else (a.mean() - b.mean())
out.append({"metric": name, "difference": diff, "as": sign, "p": stats.ranksums(a, b)[1],
"n_sup": len(a), "n_deep": len(b), "median_diff": (np.median(b) - np.median(a)) if sign == "deep-sup" else (np.median(a) - np.median(b))})
return pd.DataFrame(out)
print("\nmain (layers 2/3 and 4 against 5 and 6, all cortical units with a depth):")
print(layer_test(Ucx, [2, 4], [5, 6]).to_string(index=False))
print("expected: 2.59 ms, 109 deg2, 0.09, 6.6 ms (each +/- 10%)")
print("\nvariant, raw MI in place of log10 MI (the other three rows are the main result again):")
print(layer_test(Ucx, [2, 4], [5, 6], mi=lambda v: v).to_string(index=False))
print("\nvariant, the code's split (layer < 5 against >= 5):")
print(layer_test(Ucx, [1, 2, 3, 4], [5, 6]).to_string(index=False))
print("\nvariant, the six cortical areas only:")
print(layer_test(Ucx[Ucx.location.isin(S.AREAS[1:4] + S.AREAS[5:])], [2, 4], [5, 6]).to_string(index=False))
print("\nvariant, first spike within 200 ms (get_flash_PSTH.py) and decay timescale over 40-290 ms (timescale_calculation.py):")
print(layer_test(Ucx, [2, 4], [5, 6], "tfs_all_200", ("timescale_ac_code", "err_ac_code", "spike_count_ac_code")).to_string(index=False))
Ucsv = Ucx.merge(tfs_csv[["ecephys_unit_id", "time_to_first_spike_fl"]], on="ecephys_unit_id", how="left").merge(
ts_csv[["ecephys_unit_id", "timescale_ac", "err_ac", "spike_count_ac"]].rename(
columns={"timescale_ac": "ts_csv", "err_ac": "err_csv", "spike_count_ac": "cnt_csv"}), on="ecephys_unit_id", how="left")
print("\nvariant, the paper's tables for latency and timescale (as ExtDataFigure10.py reads them):")
print(layer_test(Ucsv, [2, 4], [5, 6], "time_to_first_spike_fl", ("ts_csv", "err_csv", "cnt_csv")).to_string(index=False))
layer codes among non-noise units with depth > 0:
{1: 141, 2: 9247, 4: 6470, 5: 13154, 6: 5536}
units measured here with depth > 0:
30638 ; their structures: {'VISp': 6490, 'VISam': 5208, 'VISal': 5078, 'VISrl': 4435, 'VISl': 3675, 'VISpm': 3009, 'VIS': 1348, 'VISmma': 598, 'VISli': 590, 'VISmmp': 207}
main (layers 2/3 and 4 against 5 and 6, all cortical units with a depth):
metric difference as p n_sup n_deep median_diff
time to first spike (ms) 1.337767 deep-sup 3.906916e-05 4430 5656 1.500000
RF area (deg2) 105.131760 deep-sup 7.780026e-27 4430 5656 100.000000
log10 MI 0.117054 sup-deep 7.446111e-21 4430 5656 0.213403
response decay timescale (ms) 7.711249 deep-sup 2.559392e-32 1838 2372 8.499093
expected: 2.59 ms, 109 deg2, 0.09, 6.6 ms (each +/- 10%)
variant, raw MI in place of log10 MI (the other three rows are the main result again):
metric difference as p n_sup n_deep median_diff
time to first spike (ms) 1.337767 deep-sup 3.906916e-05 4430 5656 1.500000
RF area (deg2) 105.131760 deep-sup 7.780026e-27 4430 5656 100.000000
MI 0.587319 sup-deep 7.446111e-21 4430 5656 0.623844
response decay timescale (ms) 7.711249 deep-sup 2.559392e-32 1838 2372 8.499093
variant, the code's split (layer < 5 against >= 5):
metric difference as p n_sup n_deep median_diff
time to first spike (ms) 1.310884 deep-sup 5.552925e-05 4445 5656 1.500000
RF area (deg2) 104.135612 deep-sup 2.295162e-26 4445 5656 100.000000
log10 MI 0.116381 sup-deep 1.184137e-20 4445 5656 0.213129
response decay timescale (ms) 7.625358 deep-sup 1.010511e-31 1843 2372 8.475208
variant, the six cortical areas only:
metric difference as p n_sup n_deep median_diff
time to first spike (ms) 1.247134 deep-sup 2.272293e-04 4125 5171 1.500000
RF area (deg2) 109.424668 deep-sup 4.584805e-27 4125 5171 100.000000
log10 MI 0.106689 sup-deep 9.584783e-16 4125 5171 0.199359
response decay timescale (ms) 7.577990 deep-sup 8.295942e-30 1722 2201 8.295328
variant, first spike within 200 ms (get_flash_PSTH.py) and decay timescale over 40-290 ms (timescale_calculation.py):
metric difference as p n_sup n_deep median_diff
time to first spike (ms) 2.797684 deep-sup 8.122645e-24 6165 7115 4.000000
RF area (deg2) 113.232316 deep-sup 1.089164e-42 6165 7115 100.000000
log10 MI 0.126691 sup-deep 7.276410e-33 6165 7115 0.223976
response decay timescale (ms) 6.386411 deep-sup 1.939687e-35 2758 3491 6.917262
variant, the paper's tables for latency and timescale (as ExtDataFigure10.py reads them):
metric difference as p n_sup n_deep median_diff
time to first spike (ms) 2.550243 deep-sup 6.445638e-19 5615 6615 4.000000
RF area (deg2) 110.299873 deep-sup 7.120542e-37 5615 6615 100.000000
log10 MI 0.124270 sup-deep 3.847460e-28 5615 6615 0.212794
response decay timescale (ms) 6.692416 deep-sup 1.721902e-33 2395 2970 7.454013
Result: not reproduced. With the claims' procedure (the Methods' windows of the Figure 3 metrics) all four differences have the stated sign with p < 1e-4, but only the RF area is within tolerance: deep minus superficial 1.34 ms in latency (2.59 +/- 0.26), 105 deg2 in RF area (109 +/- 10.9) and 7.71 ms in decay timescale (6.6 +/- 0.66), and superficial minus deep 0.117 in log10 MI (0.09 +/- 0.009); 1 of 4 parts, below the two the criterion needs for a partial result. Raw MI gives a difference of 0.587, further from 0.09. The code's layer split and the six areas alone change the differences by at most 0.09 ms, 4.3 deg2, 0.010 and 0.13 ms. Variants, not judged: with the scripts' windows (200 ms; 40 to 290 ms) 2.80 ms, 113 deg2, 0.127 and 6.39 ms, and with the paper's tables 2.55 ms, 110 deg2, 0.124 and 6.69 ms, 3 of 4 parts each.
anatomical-hierarchy-order¶
In the anatomical hierarchy recomputed from ref. 3, LGN is lowest, then V1, with LM, RL, LP and AL at intermediate levels and PM and AM at the top.
The scores are read from Figure3/Figure3.py at the pinned commit.
src = requests.get(f"{S.CODE_RAW}/Figure3/Figure3.py", timeout=60).text
block = re.search(r"hierarchy_score = \{(.*?)\}", src, re.S).group(1)
scores = {k: float(v) for k, v in re.findall(r"'(\w+)'\s*:\s*(-?[\d.]+)", block)}
order = sorted(scores, key=scores.get)
print("scores from the code:", scores)
print("order, low to high:", [S.NAMES[a] for a in order])
ok = (order[0] == "LGd" and order[1] == "VISp" and set(order[2:6]) == {"VISl", "VISrl", "LP", "VISal"}
and set(order[6:]) == {"VISpm", "VISam"})
print("LGN lowest, V1 next, LM/RL/LP/AL intermediate, PM/AM top:", ok)
print("same scores as the helper uses:", scores == S.HIERARCHY)
scores from the code: {'LGd': -0.515, 'VISp': -0.357, 'VISl': -0.093, 'VISrl': -0.059, 'LP': 0.105, 'VISal': 0.152, 'VISpm': 0.327, 'VISam': 0.441}
order, low to high: ['LGN', 'V1', 'LM', 'RL', 'LP', 'AL', 'PM', 'AM']
LGN lowest, V1 next, LM/RL/LP/AL intermediate, PM/AM top: True
same scores as the helper uses: True
Result: reproduced. The scores in the paper's code order the areas LGN, V1, LM, RL, LP, AL, PM, AM, as stated. The bootstrap test of RL against LM needs the connectivity data of ref. 3 and was not run.
Figure 4 and Extended Data Figure 9k to p: the paper's behaviour metrics tables¶
data/change_modulation_data.csv (one row per unit, 4,057 units) and data/decoding_data.csv
(one row per session and area) at the pinned commit; per-region means and Pearson against the
hierarchy scores, as Figure4.py.
cm = repo_csv("change_modulation_data.csv")
dec = repo_csv("decoding_data.csv")
REG = ["LGd", "V1", "LM", "RL", "LP", "AL", "PM", "AM"]
print("change modulation table:", cm.shape, " units per region:", cm.Region.value_counts().reindex(REG).to_dict())
print("decoding table:", dec.shape, " sessions per region:", dec.Region.value_counts().reindex(REG).to_dict())
print(cm.describe().T[["count", "min", "50%", "max"]].round(3).to_string())
def region_corr(df, col, f=np.nanmean):
m = pd.Series({r: f(df[col][df.Region == r].values) for r in REG})
r, p = stats.pearsonr(HS, m.values)
return m, r, p
change modulation table: (4057, 13) units per region: {'LGd': 297, 'V1': 915, 'LM': 724, 'RL': 483, 'LP': 311, 'AL': 482, 'PM': 314, 'AM': 531}
decoding table: (77, 4) sessions per region: {'LGd': 3, 'V1': 16, 'LM': 14, 'RL': 10, 'LP': 5, 'AL': 10, 'PM': 7, 'AM': 12}
count min 50% max
Unnamed: 0 4057.0 0.000 2028.000 4056.000
Time To First Spike 4056.0 30.853 66.129 146.422
Change Modulation Active 4057.0 -1.000 0.229 1.000
Change Modulation Passive 4057.0 -0.743 0.147 0.971
Change Modulation Hit 4056.0 -0.805 0.226 1.000
Change Modulation Miss 4050.0 -1.000 0.137 1.000
Pre-change Response Active 4057.0 -46.557 0.998 75.907
Pre-change Response Passive 4057.0 -29.930 0.933 60.151
Change Response Active 4057.0 -30.270 2.563 112.292
Change Response Passive 4057.0 -25.558 1.878 83.969
Baseline Rate Active 4057.0 0.000 2.221 134.077
Baseline Rate Passive 4057.0 0.000 3.699 117.509
fig4e-behaviour-latency-hierarchy¶
During behaviour the first spike latency of image responses is correlated with the anatomical hierarchy score.
m, r, p = region_corr(cm, "Time To First Spike")
print("mean time to first spike (ms):", m.round(2).to_dict())
print(f"Pearson r = {r:.4f}, p = {p:.3g}")
mean time to first spike (ms): {'LGd': 60.09, 'V1': 65.63, 'LM': 68.4, 'RL': 68.59, 'LP': 70.57, 'AL': 67.98, 'PM': 77.16, 'AM': 74.11}
Pearson r = 0.9171, p = 0.00134
Result: reproduced. Mean latency ranges from 60.1 ms (LGd) to 77.2 ms (PM); r = 0.917, p = 0.0013.
fig4f-change-modulation-hierarchy¶
During active behaviour the change modulation index is positive in every area and increases along the hierarchy from LGN to AM, and it is larger during active behaviour than during passive replay.
ma, ra, pa_ = region_corr(cm, "Change Modulation Active")
mp, rp, pp = region_corr(cm, "Change Modulation Passive")
print("active means:", ma.round(4).to_dict())
print("passive means:", mp.round(4).to_dict())
print(f"active: Pearson r = {ra:.4f}, p = {pa_:.3g} (expected 0.83 +/- 0.083); passive r = {rp:.4f}, p = {pp:.3g}")
print("active mean > 0 in every area:", bool((ma > 0).all()))
print("(a) active > passive in every area:", bool((ma > mp).all()), " areas:", int((ma > mp).sum()), "of 8")
print("(b) active > passive in the mean over areas:", ma.mean() > mp.mean(), round(ma.mean(), 4), round(mp.mean(), 4))
active means: {'LGd': 0.1112, 'V1': 0.2078, 'LM': 0.2277, 'RL': 0.1971, 'LP': 0.239, 'AL': 0.2181, 'PM': 0.312, 'AM': 0.2511}
passive means: {'LGd': 0.0985, 'V1': 0.1388, 'LM': 0.1413, 'RL': 0.1081, 'LP': 0.1268, 'AL': 0.1551, 'PM': 0.1941, 'AM': 0.1582}
active: Pearson r = 0.8279, p = 0.0112 (expected 0.83 +/- 0.083); passive r = 0.7297, p = 0.0399
active mean > 0 in every area: True
(a) active > passive in every area: True areas: 8 of 8
(b) active > passive in the mean over areas: True 0.2205 0.1401
Result: reproduced. Active CMI: r = 0.828, p = 0.011 (paper 0.83, 0.011); positive in all eight areas (0.111 to 0.312); above the passive mean in all eight areas (mean 0.221 against 0.140).
ed9k-cmi-hit-not-miss¶
The change modulation index correlates with hierarchy on hit trials but not on miss trials.
mh, rh, ph = region_corr(cm, "Change Modulation Hit")
mm, rm, pm = region_corr(cm, "Change Modulation Miss")
print(f"hit: r = {rh:.4f}, p = {ph:.3g} (expected 0.85 +/- 0.085); miss: r = {rm:.4f}, p = {pm:.3g} (paper 0.51, 0.2)")
print("hit means:", mh.round(4).to_dict())
print("miss means:", mm.round(4).to_dict())
hit: r = 0.8531, p = 0.00707 (expected 0.85 +/- 0.085); miss: r = 0.5098, p = 0.197 (paper 0.51, 0.2)
hit means: {'LGd': 0.1005, 'V1': 0.2026, 'LM': 0.2263, 'RL': 0.2122, 'LP': 0.2399, 'AL': 0.2141, 'PM': 0.3154, 'AM': 0.261}
miss means: {'LGd': 0.0735, 'V1': 0.1475, 'LM': 0.1378, 'RL': 0.1072, 'LP': 0.101, 'AL': 0.125, 'PM': 0.1852, 'AM': 0.1328}
Result: reproduced. Hit trials r = 0.853, p = 0.007; miss trials r = 0.510, p = 0.197, as the paper reports.
ed9lmn-other-responses-not-correlated¶
Baseline firing rate, the response to the pre-change image and the response to the change image during the task are not correlated with hierarchy score.
Result: reproduced. None of the three active measures is correlated with the hierarchy (baseline r = -0.584, p = 0.13; pre-change r = -0.239, p = 0.57; change r = 0.107, p = 0.80). Nor are the passive ones (p = 0.079 to 0.66).
fig4h-decoder-behaviour-hierarchy¶
The trial-wise covariation of change decoder output with the mouse's choices increases along the hierarchy, starting with no correlation in LGN.
col = "Correleation of decoder prediction and mouse behavior"
m = dec.groupby("Region")[col].mean().reindex(REG)
sem = dec.groupby("Region")[col].agg(lambda v: np.std(v) / np.sqrt(len(v))).reindex(REG)
r, p = stats.pearsonr(HS, m.values)
print("mean per region:", m.round(4).to_dict())
print(f"Pearson r = {r:.4f}, p = {p:.3g} (expected 0.88 +/- 0.088)")
print(f"LGN mean +/- s.e.m.: {m['LGd']:.4f} +/- {sem['LGd']:.4f}; interval includes zero: {m['LGd'] - sem['LGd'] <= 0 <= m['LGd'] + sem['LGd']}")
mean per region: {'LGd': -0.0148, 'V1': 0.1073, 'LM': 0.1335, 'RL': 0.2514, 'LP': 0.2264, 'AL': 0.1889, 'PM': 0.2326, 'AM': 0.2791}
Pearson r = 0.8812, p = 0.00383 (expected 0.88 +/- 0.088)
LGN mean +/- s.e.m.: -0.0148 +/- 0.0319; interval includes zero: True
Result: reproduced. r = 0.881, p = 0.0038 (paper 0.88, 0.004); LGN -0.015 +/- 0.032 includes zero.
ed9p-decoders-above-chance¶
Change decoders trained on 20 units performed better than chance in all eight areas, and decoder accuracy was not correlated with hierarchy score.
acc = dec.groupby("Region")["Decoder accuracy"]
t = pd.DataFrame({"mean": acc.mean(), "sem": acc.agg(lambda v: np.std(v) / np.sqrt(len(v))), "n": acc.size(),
"t_p_one_sided": acc.agg(lambda v: stats.ttest_1samp(v, 0.5, alternative="greater").pvalue)}).reindex(REG)
t["(a) mean-sem > 0.5"] = t["mean"] - t["sem"] > 0.5
t["(b) p < 0.05"] = t.t_p_one_sided < 0.05
print(t.round(5).to_string())
r, p = stats.pearsonr(HS, t["mean"].values)
print(f"accuracy against hierarchy: r = {r:.4f}, p = {p:.3g}")
mean sem n t_p_one_sided (a) mean-sem > 0.5 (b) p < 0.05 Region LGd 0.68849 0.00543 3 0.00062 True True V1 0.80053 0.00996 16 0.00000 True True LM 0.80585 0.01332 14 0.00000 True True RL 0.76132 0.01706 10 0.00000 True True LP 0.74952 0.04250 5 0.00315 True True AL 0.81956 0.01437 10 0.00000 True True PM 0.83545 0.01971 7 0.00000 True True AM 0.78223 0.01749 12 0.00000 True True accuracy against hierarchy: r = 0.5662, p = 0.143
Result: reproduced. Every area is above 0.5 under both readings (means 0.688 to 0.835; the smallest, LGN with 3 sessions, has one-sided p = 0.0006), and accuracy is not correlated with the hierarchy (r = 0.566, p = 0.14).
methods-gaze-stability¶
Across 50 mice with eye tracking, 95% of gaze locations fell within 6.4 +/- 2.1 degrees of the mean, with a maximum of 13.6 degrees.
Filtered gaze in degrees (screen_coordinates_spherical), frames with a finite value; per mouse the
95th percentile of the distance from the mean gaze: (a) radial, (b) the larger of the two axes'
95th percentiles of the absolute deviation. Raw gaze is a variant.
rows = []
for subj, g in gaze.items():
for kind, xy in g.items():
rad, ax_, n = S.gaze_radius(xy)
rows.append({"subject": subj, "kind": kind, "radial": rad, "per_axis": ax_, "frames": n, "of": len(xy)})
gz = pd.DataFrame(rows)
for kind in ("filtered", "raw"):
d = gz[gz.kind == kind]
print(f"{kind}: mice {len(d)}; (a) radial {d.radial.mean():.2f} +/- {d.radial.std(ddof=0):.2f} (sd ddof 0), "
f"sd ddof 1 {d.radial.std():.2f}, max {d.radial.max():.2f}; (b) per axis {d.per_axis.mean():.2f} +/- {d.per_axis.std(ddof=0):.2f}, max {d.per_axis.max():.2f}")
d = gz[gz.kind == "filtered"]
print("frames kept per mouse (fraction):", round((d.frames / d["of"]).min(), 3), "to", round((d.frames / d["of"]).max(), 3))
print("radial values, sorted:", np.round(np.sort(d.radial.values), 2).tolist())
filtered: mice 50; (a) radial 6.59 +/- 2.18 (sd ddof 0), sd ddof 1 2.20, max 14.29; (b) per axis 6.21 +/- 2.13, max 13.37 raw: mice 50; (a) radial 6.99 +/- 2.31 (sd ddof 0), sd ddof 1 2.33, max 15.96; (b) per axis 6.51 +/- 2.27, max 15.06 frames kept per mouse (fraction): 0.536 to 0.972 radial values, sorted: [2.34, 3.55, 3.92, 4.03, 4.22, 4.39, 4.88, 4.99, 5.06, 5.1, 5.19, 5.23, 5.27, 5.32, 5.32, 5.38, 5.57, 5.61, 5.64, 5.7, 5.73, 5.79, 5.93, 6.0, 6.06, 6.21, 6.21, 6.29, 6.32, 6.76, 7.04, 7.08, 7.19, 7.37, 7.38, 7.4, 7.76, 7.83, 7.86, 7.87, 7.89, 8.03, 8.11, 8.45, 8.65, 8.8, 8.9, 9.42, 14.27, 14.29]
Result: reproduced. Over the 50 mice with gaze mapping, the radial 95% distance is 6.59 +/- 2.18 degrees with a maximum of 14.29 (paper 6.4 +/- 2.1 and 13.6, both within tolerance); the per-axis reading gives 6.21 +/- 2.13 and 13.37.
Summary¶
Claims tested in this notebook; the CCG claims are summarized in ccg_siegle_2021.ipynb.
| Claim | Data | Outcome |
|---|---|---|
| methods-sessions-passing-qc | 000021, 000022 | reproduced |
| methods-probes-ccf | 000021, 000022 | reproduced (284 of 332) |
| ed1jk-units-per-experiment | 000021, 000022 | reproduced |
| fig1g-receptive-fields | 000021, 000022 | partially reproduced (60.7%; hippocampus 1.58%) |
| link-unit-tables | 000021, 000022 and the paper's tables | partially reproduced (ids and timescales match; latencies do not) |
| fig3c-latency-hierarchy | 000021, 000022 | reproduced, r = 0.959 |
| ed9a-latency-pairs | 000021, 000022 | partially reproduced (27 of 28 pairs) |
| fig3f-rf-size-hierarchy | 000021, 000022 | reproduced, r = 0.958 |
| ed9d-rf-size-pairs | 000021, 000022 | partially reproduced (26 of 28) |
| fig3i-modulation-index-hierarchy | 000021, 000022 | reproduced, r = -0.954 |
| ed9e-modulation-index-pairs | 000021, 000022 | not reproduced (24 of 28) |
| ed9fg-intrinsic-timescale-not-correlated | 000021, 000022 | reproduced (p = 0.71; most fits at the 1 ms bound) |
| fig3l-response-decay-timescale-hierarchy | 000021, 000022 | reproduced, r = 0.804 |
| ed9h-response-decay-pairs | 000021, 000022 | not reproduced (24 of 28) |
| fig3-unit-counts | 000021, 000022 | not reproduced (13,344 or 24,459, and 5,842); the scripts' windows give counts within tolerance |
| ed9bc-spontaneous-latency | 000021, 000022 | reproduced |
| ed9ij-firing-rate-not-correlated | 000021, 000022 | reproduced |
| ed10b-layers | 000021, 000022 with the layer table | not reproduced (1 of 4 parts, the RF area); the scripts' windows give 3 of 4 |
| anatomical-hierarchy-order | the paper's code | reproduced (order only) |
| fig4e, fig4f, ed9k, ed9lmn, fig4h, ed9p | the paper's behaviour tables | reproduced |
| methods-gaze-stability | 000021, 000022 | reproduced |
| fig4b-behaviour-performance | not deposited | blocked, data missing |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. The paper took the receptive field, modulation index and drifting grating rate of each unit from the AllenSDK's precomputed metrics, which are not in the NWB files; they are recomputed here from the spike times as the Methods describe them. The preferred grating condition (the stimulus condition with the highest mean count) and firing_rate_dg (spikes over all grating presentations divided by their duration) are assumptions; firing_rate_dg enters only as the > 0.1 Hz cut, which removes 560 of 25,019 units (fig3-unit-counts). Latencies and response decay timescales follow the Methods' windows (250 ms; 0 to 250 ms); the scripts' windows (200 ms in get_flash_PSTH.py; 40 to 290 ms in timescale_calculation.py) are shown as variants, and they change the unit counts, the layer differences and some pairwise tests. The paper's latency table matches none of the procedures tried (at most 16% of units within 1 ms), and no script of the repository writes it. A unit needs at least 10 spikes over the Gabor presentations for a receptive field p-value, as in the AllenSDK.
What was found about the data. Dandiset 000253 is an OpenScope oddball experiment and holds none of
the paper's recordings, though its metadata names the paper as IsSupplementTo. The paper's data
are in 000021 and 000022, which match the paper's counts of sessions, genotypes and probes. In
every session the electrodes' y and z columns are identical, so the left-right CCF coordinate
is lost. The units tables carry no stimulus metrics (p_value_rf, area_rf, mod_idx_dg), which the
paper's figures use. The behaviour experiments of Figure 4 are not deposited as recordings; only
the paper's per-unit and per-session metric tables are public.