Shin et al. (2025) from dandiset 000248: the Neuropixels claims¶
Paper. Shin H, Ogando MB, Abdeladim L, ... Lecoq J, Adesnik H (2025). Recurrent pattern completion drives the neocortical representation of sensory inference. Nature Neuroscience 28, 2319-2329. doi:10.1038/s41593-025-02055-5. The claims file was written from the version of record (main text, figure legends, Methods and Extended Data legends); the Supplementary Information was not available to it.
Dandiset. DANDI:000248, version
0.240502.2344. The paper cites the dandiset without a version. This is its only published version,
so it is pinned. The draft has the same 78 assets and size; it was last modified at the time of
publication (2024-05-02). The paper's analysis code reads a copy named 00248_v240130 (January
2024), before this version; no earlier version was published, so a difference between the two
cannot be checked.
What the dandiset holds. The 12 Neuropixels sessions (12 mice), each a _ogen.nwb file with
the units of all six probes, the stimulus tables, the receptive field and size blocks, and in 10
sessions eye tracking; and 66 _probe-N_ecephys.nwb files with LFP, which are not used. The 2p
imaging, 2p mesoscope and 2p holography data of the paper are not in it.
Method sources. The claims file gives each method; where it leaves a detail open, the detail
is taken from the paper or from the analysis code the claims file cites
(github.com/hs13/IC_paper_code at commit 50e04ea), and the helper module illusion.py says which.
| Claims | Status |
|---|---|
| 13 Neuropixels claims (unit counts, functional subsets, receptive fields, fixed gaze, eye, cross-correlograms, population responses) | Tested here |
| 5 decoding claims (fig3a, fig3c, fig3d, ed7b, ed7a) | Tested in decoding_shin_2025.ipynb, with further seeds in decoding_checks_shin_2025.ipynb (reading a) and decoding_seed1_shin_2025.ipynb, decoding_seed2_shin_2025.ipynb (reading b) |
| supp3-decoder-types | decoding_checks_shin_2025.ipynb: the SVMs are run; blocked (too_large) for the network |
| 21 claims on the 2p imaging, mesoscope and holography data | Blocked (data_missing): not in the dandiset |
import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import time
import warnings
from itertools import combinations
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from joblib import Parallel, delayed
from scipy import stats
import illusion as il
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
N_JOBS = 4
started = time.time()
sessions = il.load_all()
print(f"loaded {len(sessions)} sessions of dandiset {il.DANDISET_ID} version {il.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
loaded 12 sessions of dandiset 000248 version 0.240502.2344 in 133 s
Data check¶
Each check of the description against the values is stated under its output.
assets = pd.DataFrame(il.all_assets())
assets["subject"] = assets.path.str.split("/").str[0]
assets["kind"] = np.where(assets.path.str.endswith("_ogen.nwb"), "session (_ogen)", "probe LFP (_ecephys)")
print("Assets:", len(assets), "; total size (GB):", round(assets["size"].sum() / 1e9, 1))
print(pd.crosstab(assets.subject, assets.kind).to_string())
info = pd.DataFrame([{"subject": s["subject"], "session": s["session_id"], "genotype": s["genotype"].split("/")[0],
"age": s["age"], "units_in_file": s["n_units_file"],
"visual_area_units": len(s["units"]), "eye_tracking": s["eye"] is not None}
for s in sessions])
print()
print(info.to_string(index=False))
print("\nMice by genotype:", info.genotype.value_counts().to_dict())
print("Paper, Methods (Animals): Sst-IRES-Cre;Ai32 n = 10 and PV-IRES-Cre;Ai32 n = 4; Results and legends: 12 sessions from 12 mice.")
Assets: 78 ; total size (GB): 186.0
kind probe LFP (_ecephys) session (_ogen)
subject
sub-619293 0 1
sub-619296 6 1
sub-620333 6 1
sub-620334 6 1
sub-625545 6 1
sub-625554 6 1
sub-625555 6 1
sub-630506 6 1
sub-631510 6 1
sub-631570 6 1
sub-633229 6 1
sub-637484 6 1
subject session genotype age units_in_file visual_area_units eye_tracking
619293 1184980079 Sst-IRES-Cre P141D 2136 709 False
619296 1187930705 Sst-IRES-Cre P154D 1918 676 True
620333 1188137866 Sst-IRES-Cre P148D 2593 807 True
620334 1189887297 Sst-IRES-Cre P154D 2092 669 True
625545 1182865981 Sst-IRES-Cre P89D 2793 1123 True
625554 1181330601 Pvalb-IRES-Cre P82D 2930 953 True
625555 1183070926 Pvalb-IRES-Cre P90D 2621 1269 True
630506 1192952695 Sst-IRES-Cre P92D 2517 799 True
631510 1196157974 Sst-IRES-Cre P99D 2386 740 True
631570 1194857009 Pvalb-IRES-Cre P92D 1789 688 True
633229 1199247593 Sst-IRES-Cre P101D 3026 1211 True
637484 1208667752 Sst-IRES-Cre P92D 2373 934 False
Mice by genotype: {'Sst-IRES-Cre': 9, 'Pvalb-IRES-Cre': 3}
Paper, Methods (Animals): Sst-IRES-Cre;Ai32 n = 10 and PV-IRES-Cre;Ai32 n = 4; Results and legends: 12 sessions from 12 mice.
The dandiset holds 12 sessions from 12 mice, as the Results and legends state. One file per session
(_ogen.nwb) holds the units of all probes; sub-619293 has no probe LFP files. The genotypes are 9
Sst-IRES-Cre;Ai32 and 3 Pvalb-IRES-Cre;Ai32 mice, against the 10 and 4 (14 mice) of the Methods'
Animals list: the Methods list two more mice than the analyses use, and the split differs.
rows = []
for s in sessions:
wd, q, ar = s["waveform_duration_all"], s["quality_all"], s["area_all"]
rows.append({"subject": s["subject"], "units": len(wd), "noise_label": int(np.sum(q == "noise")),
"in_6_visual_areas": int(np.sum(ar != None)),
"RS_all_areas": int(np.sum(wd > 0.4)),
**{f"{a}_RS": int(np.sum((s["units"].area == a) & s["units"].rs)) for a in il.AREAS}})
units_tab = pd.DataFrame(rows)
print(units_tab.to_string(index=False))
allwd = np.concatenate([s["waveform_duration_all"] for s in sessions])
print("\nwaveform_duration (ms): min %.3f, median %.3f, max %.3f; fraction > 0.4 ms: %.3f" %
(allwd.min(), np.median(allwd), allwd.max(), np.mean(allwd > 0.4)))
print("quality values:", sorted(set(np.concatenate([s["quality_all"] for s in sessions]))))
labs = pd.concat([s["unit_labels_all"] for s in sessions]).groupby(level=0).sum()
vis = labs[[il.area_of(x) is not None for x in labs.index]]
print("\nVisual cortex labels used (units):", vis.to_dict())
print("Other labels with 'VIS':", [x for x in labs.index if "VIS" in x and il.area_of(x) is None])
print("Spike times sorted within units (first 50 units per session):", all(s["sorted_ok"] for s in sessions))
subject units noise_label in_6_visual_areas RS_all_areas V1_RS LM_RS RL_RS AL_RS PM_RS AM_RS
619293 2136 488 709 1490 89 0 156 106 151 107
619296 1918 360 676 1398 168 4 87 94 95 93
620333 2593 476 807 1834 114 119 128 91 102 129
620334 2092 455 669 1334 68 183 100 0 90 82
625545 2793 523 1123 2126 190 181 122 160 173 134
625554 2930 461 953 2116 179 0 136 139 193 183
625555 2621 632 1269 2083 314 0 208 159 214 198
630506 2517 515 799 1889 58 131 140 121 208 0
631510 2386 546 740 1603 129 105 75 87 70 96
631570 1789 436 688 1278 142 74 0 105 158 109
633229 3026 539 1211 2161 181 172 130 179 210 142
637484 2373 531 934 1590 172 122 279 0 67 96
waveform_duration (ms): min 0.014, median 0.577, max 2.692; fraction > 0.4 ms: 0.716
quality values: [np.str_('good'), np.str_('noise')]
Visual cortex labels used (units): {'VISal1': 9, 'VISal2/3': 390, 'VISal4': 279, 'VISal5': 555, 'VISal6a': 233, 'VISal6b': 27, 'VISam1': 25, 'VISam2/3': 595, 'VISam4': 194, 'VISam5': 409, 'VISam6a': 251, 'VISam6b': 173, 'VISl1': 45, 'VISl2/3': 280, 'VISl4': 245, 'VISl5': 496, 'VISl6': 4, 'VISl6a': 235, 'VISl6b': 18, 'VISp1': 19, 'VISp2/3': 629, 'VISp4': 415, 'VISp5': 641, 'VISp6a': 422, 'VISp6b': 34, 'VISpm1': 23, 'VISpm2/3': 675, 'VISpm4': 308, 'VISpm5': 653, 'VISpm6a': 334, 'VISpm6b': 18, 'VISrl1': 36, 'VISrl2/3': 503, 'VISrl4': 577, 'VISrl5': 522, 'VISrl6a': 256, 'VISrl6b': 50}
Other labels with 'VIS': ['VIS2/3', 'VIS4', 'VIS5', 'VIS6a', 'VIS6b']
Spike times sorted within units (first 50 units per session): True
Units carry a waveform duration in ms (0.014 to 2.7 ms, a bimodal range as expected for
trough-to-peak widths), a peak channel whose CCF label gives the area and layer, and a quality
label of good or noise. The Methods say that "all units not classified as noise were packaged
into Neurodata Without Borders files"; the files nevertheless hold units labelled noise (counts
above). The analyses below keep them, as the paper's code does (it applies no quality filter), and
the counts in the claims below show the paper's numbers include them. LM has fewer than 50 RS units
in four sessions (none in three), and RL none in one, as the Fig. 3d legend's session counts imply.
rows = []
for s in sessions:
for blk, b in s["blocks"].items():
codes, n = np.unique(b["codes"], return_counts=True)
rows.append({"subject": s["subject"], "block": blk, "table_rows": b["n_rows"], "trials": len(b["codes"]),
"ITI_rows_with_nonzero_frame": b["iti_rows_nonzero_frame"],
"n_IC1": int(np.sum(b["codes"] == il.IC1)), "n_other_min": int(n.min()), "n_other_max": int(n.max()),
"dur_s_min": b["duration"].min().round(3), "dur_s_max": b["duration"].max().round(3),
"iti_s_median": np.median(b["iti"]).round(3)})
blk_tab = pd.DataFrame(rows)
print(blk_tab.groupby("block").agg({"trials": ["min", "max"], "table_rows": ["min", "max"], "ITI_rows_with_nonzero_frame": "max",
"n_IC1": ["min", "max"], "n_other_min": "min", "n_other_max": "max",
"dur_s_min": "min", "dur_s_max": "max", "iti_s_median": ["min", "max"]}).to_string())
b = sessions[0]["blocks"]["ICwcfg1"]
codes, n = np.unique(b["codes"], return_counts=True)
print("\nICwcfg1 trials per image code, first session:", dict(zip(codes.tolist(), n.tolist())))
rf = pd.DataFrame([{"subject": s["subject"], "rf_rows": s["rf"]["n_rows"], "rf_trials": len(s["rf"]["pos"]),
"classic_per_position": np.bincount(s["rf"]["pos"][~s["rf"]["inverse"]], minlength=9).tolist(),
"inverse_per_position": np.bincount(s["rf"]["pos"][s["rf"]["inverse"]], minlength=9).tolist(),
"row_dur_s": (round(s["rf"]["row_dur"].min(), 3), round(s["rf"]["row_dur"].max(), 3)),
"directions_per_trial": s["rf"]["directions_per_trial"][0],
"longest_spontaneous_s": round(s["spont_interval"][1] - s["spont_interval"][0], 1)}
for s in sessions])
print(rf.to_string(index=False))
trials table_rows ITI_rows_with_nonzero_frame n_IC1 n_other_min n_other_max dur_s_min dur_s_max iti_s_median
min max min max max min max min max min max min max
block
ICkcfg0 660 660 1380 1380 0 30 30 30 30 0.4 0.401 0.4 0.4
ICkcfg1 660 660 1380 1380 0 30 30 30 30 0.4 0.401 0.4 0.4
ICwcfg0 660 660 1380 1380 0 30 30 30 30 0.4 0.401 0.4 0.4
ICwcfg1 5300 5300 10720 10720 0 400 400 50 400 0.4 0.401 0.4 0.4
ICwcfg1 trials per image code, first session: {0: 400, 101: 400, 105: 400, 106: 400, 107: 400, 109: 400, 110: 400, 111: 400, 506: 50, 511: 50, 1105: 400, 1109: 400, 1201: 400, 1299: 400, 1301: 50, 1302: 50, 1303: 50, 1304: 50, 1305: 50, 1306: 50, 1307: 50, 1308: 50}
subject rf_rows rf_trials classic_per_position inverse_per_position row_dur_s directions_per_trial longest_spontaneous_s
619293 720 180 [10, 10, 10, 10, 10, 10, 10, 10, 10] [10, 10, 10, 10, 10, 10, 10, 10, 10] (0.25, 0.284) (0.0, 45.0, 90.0, 135.0) 312.3
619296 720 180 [10, 10, 10, 10, 10, 10, 10, 10, 10] [10, 10, 10, 10, 10, 10, 10, 10, 10] (0.25, 0.284) (0.0, 45.0, 90.0, 135.0) 312.3
620333 720 180 [10, 10, 10, 10, 10, 10, 10, 10, 10] [10, 10, 10, 10, 10, 10, 10, 10, 10] (0.25, 0.284) (0.0, 45.0, 90.0, 135.0) 312.3
620334 720 180 [10, 10, 10, 10, 10, 10, 10, 10, 10] [10, 10, 10, 10, 10, 10, 10, 10, 10] (0.25, 0.284) (0.0, 45.0, 90.0, 135.0) 312.3
625545 720 180 [10, 10, 10, 10, 10, 10, 10, 10, 10] [10, 10, 10, 10, 10, 10, 10, 10, 10] (0.25, 0.284) (0.0, 45.0, 90.0, 135.0) 312.3
625554 720 180 [10, 10, 10, 10, 10, 10, 10, 10, 10] [10, 10, 10, 10, 10, 10, 10, 10, 10] (0.25, 0.284) (0.0, 45.0, 90.0, 135.0) 312.3
625555 720 180 [10, 10, 10, 10, 10, 10, 10, 10, 10] [10, 10, 10, 10, 10, 10, 10, 10, 10] (0.25, 0.284) (0.0, 45.0, 90.0, 135.0) 312.3
630506 720 180 [10, 10, 10, 10, 10, 10, 10, 10, 10] [10, 10, 10, 10, 10, 10, 10, 10, 10] (0.25, 0.284) (0.0, 45.0, 90.0, 135.0) 312.3
631510 720 180 [10, 10, 10, 10, 10, 10, 10, 10, 10] [10, 10, 10, 10, 10, 10, 10, 10, 10] (0.25, 0.3) (0.0, 45.0, 90.0, 135.0) 312.3
631570 720 180 [10, 10, 10, 10, 10, 10, 10, 10, 10] [10, 10, 10, 10, 10, 10, 10, 10, 10] (0.25, 0.284) (0.0, 45.0, 90.0, 135.0) 312.3
633229 720 180 [10, 10, 10, 10, 10, 10, 10, 10, 10] [10, 10, 10, 10, 10, 10, 10, 10, 10] (0.25, 0.284) (0.0, 45.0, 90.0, 135.0) 312.3
637484 720 180 [10, 10, 10, 10, 10, 10, 10, 10, 10] [10, 10, 10, 10, 10, 10, 10, 10, 10] (0.25, 0.284) (0.0, 45.0, 90.0, 135.0) 312.3
Each IC block table alternates image rows with ITI rows whose frame is 0, as the code assumes; no
ITI row carries an image. ICwcfg1 has 5,300 trials per session: 400 of blank, X, TC1, IC1, LC1,
TC2, LC2, IC2, TRE1, TRE2, XRE1 and XRE2, and 50 of IRE1, IRE2 and the eight segment images. The
Methods say IC and LC images were shown 400 times and each other image at least 50 times; the
data agree, with TRE, XRE, X and TC images at 400. The other blocks have 30 trials of each of the
22 images. Images last 0.400 s with a 0.4 s ITI, as stated. The receptive field block has 180
trials (four 0.25 s directions each, 0 s ITI) at 9 positions for standard and inverse patches,
10 each, as the Methods state. The stimulus tables carry the image index (frame) without image
names; the names are the code's list, and the image templates stored in the files are consistent
with it (shown at the end of this section).
# Clock check: decode the four training images from V1 in the 400 ms after onset, in the 400 ms
# before onset, and after shuffling the labels. If spike times and the trial table share a clock,
# only the first is above chance.
def clock_check(i, subject, b, m):
rng = np.random.default_rng(100 + i)
shuffled = b["codes"].copy()
sel = np.isin(shuffled, [il.IC1, il.LC1, il.LC2, il.IC2])
shuffled[sel] = rng.permutation(shuffled[sel])
return {"subject": subject,
"post_0_400ms": il.decode(b["R"][:, m], b["codes"], [106, 107, 110, 111], [1105], seed=i)["acc"],
"pre_-400_0ms": il.decode(b["R_pre"][:, m], b["codes"], [106, 107, 110, 111], [1105], seed=i)["acc"],
"post_shuffled_labels": il.decode(b["R"][:, m], shuffled, [106, 107, 110, 111], [1105], seed=i)["acc"],
"rate_post_Hz": b["R"][:, m].mean(), "rate_pre_Hz": b["R_pre"][:, m].mean()}
def v1_block(s):
m = ((s["units"].area == "V1") & s["units"].rs).values
b = s["blocks"]["ICwcfg1"]
return {"R": b["R"][:, m], "R_pre": b["R_pre"][:, m], "codes": b["codes"]}, np.ones(m.sum(), bool)
clock = pd.DataFrame(Parallel(n_jobs=N_JOBS)(delayed(clock_check)(i, s["subject"], *v1_block(s)) for i, s in enumerate(sessions)))
print(clock.round(3).to_string(index=False))
print("chance = 0.25")
subject post_0_400ms pre_-400_0ms post_shuffled_labels rate_post_Hz rate_pre_Hz 619293 0.268 0.274 0.260 5.966 5.990 619296 0.893 0.252 0.261 3.571 3.425 620333 0.684 0.226 0.265 3.557 3.106 620334 0.708 0.258 0.244 6.326 5.116 625545 0.788 0.240 0.249 3.265 2.662 625554 0.791 0.263 0.236 4.017 3.594 625555 0.524 0.232 0.235 2.435 2.291 630506 0.746 0.236 0.243 7.033 5.566 631510 0.825 0.236 0.259 3.351 3.194 631570 0.808 0.257 0.243 3.145 2.799 633229 0.857 0.258 0.234 3.480 3.228 637484 0.554 0.277 0.238 3.081 2.797 chance = 0.25
In every session the four training images are decoded from V1 spikes after onset and not from the 400 ms before onset or with shuffled labels (values above), so the spike times and the stimulus table share a clock. The eye tracking timestamps are checked below.
rows = []
for s in sessions:
e = s["eye"]
if e is None:
rows.append({"subject": s["subject"], "eye": False})
continue
rows.append({"subject": s["subject"], "eye": True, "frame_rate_Hz": round(e["hz"], 2), "frames": e["n_frames"],
"nan_timestamps": e["n_nan_timestamps"], "nan_is_last": e["nan_timestamp_last"],
"stated_unit_pupil_xy": e["units"]["data"], "pupil_x_1_99": np.round(e["x_range"], 0).tolist(),
"pupil_y_1_99": np.round(e["y_range"], 0).tolist(), "pupil_width_median": round(e["width_median"], 1),
"eye_width_median": round(e["eye_width_median"], 1), "frac_nan_pupil": round(e["frac_nan"], 3),
"frac_likely_blink": round(e["frac_blink"], 3),
"first_frame_after_onset_s": np.round(np.percentile(e["ICwcfg1"]["lag_first_frame"], [0, 100]), 4).tolist()})
eye_tab = pd.DataFrame(rows)
print(eye_tab.to_string(index=False))
subject eye frame_rate_Hz frames nan_timestamps nan_is_last stated_unit_pupil_xy pupil_x_1_99 pupil_y_1_99 pupil_width_median eye_width_median frac_nan_pupil frac_likely_blink first_frame_after_onset_s 619293 False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 619296 True 60.02 538128.0 1.0 True meters [325.0, 359.0] [222.0, 247.0] 27.2 160.0 0.001 0.009 [0.0, 0.0167] 620333 True 60.02 539586.0 1.0 True meters [348.0, 399.0] [237.0, 275.0] 23.0 166.2 0.001 0.021 [0.0, 0.0167] 620334 True 60.02 542540.0 1.0 True meters [272.0, 347.0] [243.0, 265.0] 29.5 163.6 0.010 0.027 [0.0, 0.0167] 625545 True 60.02 538553.0 1.0 True meters [321.0, 382.0] [227.0, 253.0] 25.7 146.1 0.007 0.010 [0.0, 0.0167] 625554 True 60.02 540449.0 1.0 True meters [322.0, 387.0] [240.0, 260.0] 33.6 150.8 0.001 0.018 [0.0, 0.0167] 625555 True 60.02 538050.0 1.0 True meters [376.0, 418.0] [261.0, 283.0] 18.9 128.6 0.006 0.012 [0.0, 0.0167] 630506 True 60.02 538426.0 1.0 True meters [305.0, 394.0] [244.0, 294.0] 21.7 158.0 0.010 0.031 [0.0, 0.0167] 631510 True 60.02 543400.0 1.0 True meters [270.0, 327.0] [223.0, 249.0] 18.4 154.2 0.042 0.070 [0.0, 0.0167] 631570 True 60.02 545249.0 1.0 True meters [310.0, 363.0] [227.0, 246.0] 19.5 142.5 0.019 0.026 [0.0, 0.0167] 633229 True 60.02 542132.0 1.0 True meters [347.0, 401.0] [258.0, 281.0] 19.2 151.4 0.041 0.067 [0.0, 0.0167] 637484 False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
Ten sessions have eye tracking, as the Extended Data Fig. 2 legend states. The eye video is sampled at 60 Hz (timestamps 16.7 ms apart), where the Methods say "Videos of the eye and body were acquired at 30 Hz"; the paper's code assumes 60 Hz. The pupil position is stated in meters, but its values (pupil center near 300 to 400, pupil half-width near 20, eye half-width near 150) are camera pixels, as the code's comments say. No conversion to visual degrees is in the files. One timestamp per session is NaN, the last one, and is not used. Every IC trial has an eye frame within 17 ms of its onset, so eye frames cover the trial table; that the eye clock matches the stimulus clock is consistent with this, not shown by it.
# The image templates stored in the first session, with the code's names for each frame index.
import h5py, remfile
f = h5py.File(remfile.File(il.session_assets()[0]["url"]), "r")
g = f["stimulus/templates/ICwcfg1_presentations"]
names = ["blank", "X", "TC1", "IC1", "LC1", "TC2", "LC2", "IC2", "IRE1", "IRE2", "TRE1", "TRE2", "XRE1", "XRE2",
"In BR", "In BL", "In TL", "In TR", "Out BR", "Out BL", "Out TL", "Out TR"]
fig, axs = plt.subplots(4, 6, figsize=(12, 5.5))
for k, ax in enumerate(axs.flat):
ax.axis("off")
if k < 22:
img = g[f"ICwcfg1_presentations{k}"][::10, ::10]
ax.imshow(img, cmap="gray", vmin=0, vmax=255)
ax.set_title(f"{k}: {names[k]} ({il.IC_TRIAL_TYPES[k]})", fontsize=7)
plt.tight_layout()
plt.show()
print("template size:", g["ICwcfg1_presentations0"].shape, "; values", g["ICwcfg1_presentations0"][:].min(), "to",
g["ICwcfg1_presentations0"][:].max())
# Orientation of the real bar of the XRE images: agreement of each XRE template with the IRE templates
# (the IC images with their real bar outlined) in the central 200 x 200 pixels, where the bar crosses the gap.
T = {k: g[f"ICwcfg1_presentations{k}"][:] for k in [8, 9, 12, 13]}
h, w = T[8].shape[:2]
ctr = (slice(h // 2 - 100, h // 2 + 100), slice(w // 2 - 100, w // 2 + 100))
for k in [12, 13]:
print(f"{names[k]}: fraction of central pixels equal to IRE1 {np.mean(T[k][ctr] == T[8][ctr]):.3f}, to IRE2 {np.mean(T[k][ctr] == T[9][ctr]):.3f}")
template size: (1200, 1920) ; values 0 to 255 XRE1: fraction of central pixels equal to IRE1 0.924, to IRE2 0.656 XRE2: fraction of central pixels equal to IRE1 0.656, to IRE2 0.924
The templates match the code's names: index 3 and 7 show the two illusory bars (IC1 at 135 degrees, IC2 at 45 degrees), 4 and 6 the recombined LC images, 8 and 9 the IC images with the real bar outlined (IRE), 10 and 11 the TRE images, 12 and 13 the XRE images, and 14 to 21 single inducer segments. The real bar of XRE1 has the orientation of IC1's illusory bar and that of XRE2 the orientation of IC2's: in the central 200 x 200 pixels, XRE1 agrees with IRE1 in 0.924 of pixels and with IRE2 in 0.656, and XRE2 the reverse.
# Recordings deposited twice: compare spike counts, unit counts and trial order between sessions.
sig = []
for s in sessions:
sig.append((s["subject"], s["n_units_file"], int(s["units"].n_spikes_bins.sum()),
tuple(s["blocks"]["ICwcfg1"]["codes"][:200]), round(s["blocks"]["ICwcfg1"]["start"][0], 3)))
pairs = [(a[0], b[0]) for a, b in combinations(sig, 2) if a[1] == b[1] or a[2] == b[2] or a[3] == b[3] or a[4] == b[4]]
print("Pairs sharing unit count, spike count, trial order or block start:", pairs if pairs else "none")
print("Trial order identical across sessions:", len({x[3] for x in sig}) == 1)
Pairs sharing unit count, spike count, trial order or block start: none Trial order identical across sessions: False
No two sessions share their unit count, spike count or first block start; every session is a distinct recording. The trial order of ICwcfg1 is the same in every session (one stimulus script), which the Methods' "Trial order within each block was randomized" does not contradict.
results-npx-units-per-session¶
The Neuropixels sessions hold 150 +/- 20 V1 RS units and 597 +/- 46 higher visual area RS units per session (mean +/- s.e.m., 12 sessions).
v1 = units_tab["V1_RS"]
hva = units_tab[[f"{a}_RS" for a in il.AREAS[1:]]].sum(axis=1)
print("V1 RS per session:", v1.tolist())
print("HVA RS per session:", hva.tolist())
print("V1 RS: total %d, mean %.1f +/- %.1f (s.e.m.)" % (v1.sum(), v1.mean(), v1.sem()))
print("HVA RS: total %d, mean %.1f +/- %.1f (s.e.m.)" % (hva.sum(), hva.mean(), hva.sem()))
print("Paper: 150 +/- 20 and 597 +/- 46; Methods: 1,804 V1 RS units")
good = [(int(((s["units"].area == "V1") & s["units"].rs & (s["units"].quality == "good")).sum()),
int(((s["units"].area != "V1") & s["units"].rs & (s["units"].quality == "good")).sum())) for s in sessions]
good = np.array(good)
print("Variant without units labelled noise: V1 %.1f +/- %.1f, HVA %.1f +/- %.1f" %
(good[:, 0].mean(), stats.sem(good[:, 0]), good[:, 1].mean(), stats.sem(good[:, 1])))
V1 RS per session: [89, 168, 114, 68, 190, 179, 314, 58, 129, 142, 181, 172] HVA RS per session: [520, 373, 569, 455, 770, 651, 779, 600, 433, 446, 833, 564] V1 RS: total 1804, mean 150.3 +/- 19.8 (s.e.m.) HVA RS: total 6993, mean 582.8 +/- 43.3 (s.e.m.) Paper: 150 +/- 20 and 597 +/- 46; Methods: 1,804 V1 RS units Variant without units labelled noise: V1 131.4 +/- 18.3, HVA 509.9 +/- 40.7
Result: reproduced. V1 holds 150.3 +/- 19.8 RS units per session (1,804 in all) and the higher visual areas 582.8 +/- 43.3, both within the paper's 150 +/- 20 and 597 +/- 46.
methods-npx-ic-encoder-count¶
Across the ICwcfg0 and ICwcfg1 blocks, 159 of 1,804 V1 RS units are IC-encoders.
IC-encoders as in the claims file (Kruskal-Wallis across blank, IC1, LC1, LC2, IC2 at p < 0.05, Tukey-Kramer against blank, AUROC sign; code analyzeStaticICtxi.m), on the IC response (0 to 400 ms), in each block with that block's IC1 and IC2 images.
groups = {}
for s in sessions:
u = s["units"]
for blk in ["ICwcfg1", "ICwcfg0"]:
b = s["blocks"][blk]
g = il.ic_groups(b["R"], b["codes"])
g.index = u.index
groups[(s["subject"], blk)] = pd.concat([u[["area", "layer", "rs", "quality"]], g], axis=1)
def v1rs(subject, blk):
g = groups[(subject, blk)]
return g[(g.area == "V1") & g.rs]
per = pd.DataFrame([{"subject": s["subject"], "V1_RS": len(v1rs(s["subject"], "ICwcfg1")),
"IC_enc_ICwcfg1": int(v1rs(s["subject"], "ICwcfg1").ic_enc.sum()),
"IC_enc_ICwcfg0": int(v1rs(s["subject"], "ICwcfg0").ic_enc.sum()),
"IC_enc_either": int((v1rs(s["subject"], "ICwcfg1").ic_enc | v1rs(s["subject"], "ICwcfg0").ic_enc).sum())}
for s in sessions])
print(per.to_string(index=False))
tot = per.sum(numeric_only=True)
print("\nV1 RS units:", tot.V1_RS)
print("(a) IC-encoder in at least one block:", tot.IC_enc_either)
print("(b) sum of the two blocks' counts:", tot.IC_enc_ICwcfg1 + tot.IC_enc_ICwcfg0)
print("Paper: 159 of 1,804 (tolerance 15.9)")
subject V1_RS IC_enc_ICwcfg1 IC_enc_ICwcfg0 IC_enc_either 619293 89 0 0 0 619296 168 7 6 12 620333 114 6 5 11 620334 68 2 9 11 625545 190 15 15 29 625554 179 9 8 16 625555 314 9 8 17 630506 58 3 1 4 631510 129 8 7 15 631570 142 4 4 8 633229 181 10 13 23 637484 172 6 8 13 V1 RS units: 1804 (a) IC-encoder in at least one block: 159 (b) sum of the two blocks' counts: 163 Paper: 159 of 1,804 (tolerance 15.9)
Result: reproduced. 159 V1 RS units are IC-encoders in at least one of ICwcfg1 and ICwcfg0 (reading a), the paper's number; the sum of the two blocks' counts (reading b) is 163, also within 159 +/- 15.9.
ed6b-npx-v1-group-counts¶
The Neuropixels data hold 87 V1 IC-encoders and 636 V1 segment responders (the sets whose layer distribution Extended Data Fig. 6b shows).
Readings of the criterion: (a) RS units, ICwcfg1; (b) all units not classified as noise, ICwcfg1;
(c) RS units, union of ICwcfg1 and ICwcfg0. For (b), the Methods say "All units not classified as
noise were packaged into Neurodata Without Borders files and analyzed in this paper", so on the
paper's own description (b) is every unit in the files; the files also carry a quality label
noise on some units, so (b) is run both ways: (b1) all units in the files, (b2) without the units
labelled noise.
def v1all(subject, blk):
g = groups[(subject, blk)]
return g[g.area == "V1"]
rows = []
for name, sel in [("(a) RS, ICwcfg1", lambda sub: (v1rs(sub, "ICwcfg1"), None)),
("(b1) all units in the files, ICwcfg1", lambda sub: (v1all(sub, "ICwcfg1"), None)),
("(b2) units not labelled noise, ICwcfg1", lambda sub: (v1all(sub, "ICwcfg1").query("quality != 'noise'"), None)),
("(c) RS, ICwcfg1 or ICwcfg0", lambda sub: (v1rs(sub, "ICwcfg1"), v1rs(sub, "ICwcfg0")))]:
ic = seg = 0
for s in sessions:
g1, g0 = sel(s["subject"])
if g0 is None:
ic += int(g1.ic_enc.sum()); seg += int(g1.seg_resp.sum())
else:
ic += int((g1.ic_enc | g0.ic_enc).sum()); seg += int((g1.seg_resp | g0.seg_resp).sum())
rows.append({"reading": name, "IC_encoders": ic, "segment_responders": seg,
"IC_within_87+/-8.7": abs(ic - 87) <= 8.7, "seg_within_636+/-63.6": abs(seg - 636) <= 63.6})
print(pd.DataFrame(rows).to_string(index=False))
g = pd.concat([v1all(s["subject"], "ICwcfg1") for s in sessions])
lay = pd.crosstab(g.layer, [g.ic_enc, g.seg_resp])
for name, m in [("IC-encoders", g.ic_enc), ("segment responders", g.seg_resp), ("all V1 units", g.ic_enc | ~g.ic_enc)]:
print(f"\nLayer proportions, {name} (reading b1), n = {int(m.sum())}:",
(g[m].layer.value_counts(normalize=True).sort_index().round(3)).to_dict())
reading IC_encoders segment_responders IC_within_87+/-8.7 seg_within_636+/-63.6
(a) RS, ICwcfg1 79 526 True False
(b1) all units in the files, ICwcfg1 87 636 True True
(b2) units not labelled noise, ICwcfg1 79 566 True False
(c) RS, ICwcfg1 or ICwcfg0 159 684 False True
Layer proportions, IC-encoders (reading b1), n = 87: {'2/3': 0.356, '4': 0.149, '5': 0.322, '6a': 0.172}
Layer proportions, segment responders (reading b1), n = 636: {'2/3': 0.281, '4': 0.225, '5': 0.321, '6a': 0.162, '6b': 0.011}
Layer proportions, all V1 units (reading b1), n = 2160: {'1': 0.009, '2/3': 0.291, '4': 0.192, '5': 0.297, '6a': 0.195, '6b': 0.016}
Result: reproduced, resting on an assumption. Reading (b1), all V1 units in the files in ICwcfg1, gives 87 IC-encoders and 636 segment responders, the paper's two numbers exactly. The assumption is that the criterion's "all units not classified as noise" means the units of the files, which the Methods say hold all units not classified as noise; dropping the units the files label noise (b2) gives 79 and 566, and 566 is outside 636 +/- 63.6. Reading (a) gives 79 and 526, the counts of Fig. 2d, and reading (c) 159 and 684; neither puts both counts within tolerance.
ed6cd-npx-ic-encoders-all-areas¶
IC-encoders are found in every visual area, in the Neuropixels and in the calcium imaging data.
Neuropixels part only: ICwcfg1, RS units of each area, sessions with at least 50 RS units in the area. Visually responsive: Kruskal-Wallis across blank, IC1, LC1, LC2, IC2 at p < 0.05.
rows = []
for s in sessions:
g = groups[(s["subject"], "ICwcfg1")]
for a in il.AREAS:
x = g[(g.area == a) & g.rs]
if len(x) < 50:
continue
vr = x.p_kw_BK < 0.05
rows.append({"subject": s["subject"], "area": a, "RS": len(x), "IC_enc": int(x.ic_enc.sum()),
"vis_resp": int(vr.sum()), "IC_enc_frac_all": x.ic_enc.mean(),
"IC_enc_frac_visresp": x.ic_enc[vr].mean() if vr.any() else np.nan})
prop = pd.DataFrame(rows)
summary = prop.groupby("area").agg(sessions=("subject", "size"), IC_encoders=("IC_enc", "sum"), RS_units=("RS", "sum"),
frac_all_mean=("IC_enc_frac_all", "mean"), frac_all_sem=("IC_enc_frac_all", "sem"),
frac_visresp_mean=("IC_enc_frac_visresp", "mean"),
sessions_with_any=("IC_enc", lambda v: int((v > 0).sum()))).reindex(il.AREAS)
print(summary.round(4).to_string())
print("\nAreas with at least one IC-encoder:", [a for a in il.AREAS if summary.loc[a, "IC_encoders"] > 0])
sessions IC_encoders RS_units frac_all_mean frac_all_sem frac_visresp_mean sessions_with_any area V1 12 79 1804 0.0428 0.0059 0.0550 11 LM 8 48 1087 0.0417 0.0062 0.0568 8 RL 11 73 1561 0.0481 0.0080 0.0862 11 AL 10 64 1241 0.0502 0.0067 0.0707 10 PM 12 59 1731 0.0341 0.0044 0.0787 12 AM 11 49 1369 0.0350 0.0083 0.0582 9 Areas with at least one IC-encoder: ['V1', 'LM', 'RL', 'AL', 'PM', 'AM']
Result: reproduced (Neuropixels part). All six areas hold IC-encoders in ICwcfg1: from 48 (LM) to 79 (V1) per area, in 8 to 12 sessions per area, 3.4% to 5.0% of RS units on average. The calcium imaging part of the claim is not in the dandiset.
fig1e-npx-exclusive-center¶
In the Neuropixels data, 18 of 1,804 V1 RS units respond exclusively to the grating patch at the center position (the illusory gap region) of the receptive field mapping block.
Per RF trial (four directions, 1 s), the rate 0 to 1000 ms after onset, standard patches only. Reading (a), the code: signed-rank of each position's 10 trial rates against the unit's mean rate in the longest spontaneous block; the center is the position of largest mean response and the only one surviving the code's Holm step. Reading (b), the Methods: rank-sum of each position's trial rates against the rates in 1-s windows tiled over the longest spontaneous (gray screen) block (312 windows), Holm over the 9 positions, the center the only significant position (with a mean above the spontaneous mean, as "responded" implies).
rfres = {}
for s in sessions:
u = s["units"]
m = ((u.area == "V1") & u.rs).values
d = il.rf_positions(s["rf"]["R"][:, m], s["rf"]["pos"], s["rf"]["inverse"], s["spont_rate"][m], s["spont_windows"][:, m])
d.index = u.index[m]
rfres[s["subject"]] = d
rf_tab = pd.DataFrame([{"subject": k, "V1_RS": len(d), "best_center": int((d.best == 0).sum()),
"excl_center_a": int(d.excl_a.sum()), "excl_center_b": int(d.excl_b.sum()),
"kw_sig": int((d.p_kw_rf < 0.05).sum())} for k, d in rfres.items()])
print(rf_tab.to_string(index=False))
t = rf_tab.sum(numeric_only=True)
print(f"\nV1 RS units {t.V1_RS}; exclusively center-responsive: (a) {t.excl_center_a}, (b) {t.excl_center_b}; paper 18 (tolerance 1.8)")
allrf = pd.concat(rfres.values())
for name, m in [("center best and p_center (a) < 0.05/9", (allrf.best == 0) & (allrf.p_center_a * 9 < 0.05)),
("center best and p_center (b) < 0.05/9", (allrf.best == 0) & (allrf.p_center_b * 9 < 0.05)),
("(a) and (b) both", allrf.excl_a & allrf.excl_b),
("(a), units not labelled noise", allrf.excl_a & np.concatenate([s["units"].quality.values[((s["units"].area == "V1") & s["units"].rs).values] != "noise" for s in sessions]))]:
print(f"variant, {name}: {int(m.sum())}")
print("\nSteps before the exclusivity rule (all V1 RS units):")
steps = pd.DataFrame([
{"step": "center is the position of largest mean response", "units": int((allrf.best == 0).sum())},
{"step": "(a) center best and its signed-rank p < 0.05 (uncorrected)", "units": int(((allrf.best == 0) & (allrf.p_center_a < 0.05)).sum())},
{"step": "(a) center best and at least one position passes the Holm step", "units": int(((allrf.best == 0) & (allrf.n_holm_a >= 1)).sum())},
{"step": "(a) ... exactly one position passes", "units": int(((allrf.best == 0) & (allrf.n_holm_a == 1)).sum())},
{"step": "(a) ... and it is the center (= exclusive, reading a)", "units": int(allrf.excl_a.sum())},
{"step": "(b) center rank-sum p < 0.05 (uncorrected)", "units": int((allrf.p_center_b < 0.05).sum())},
{"step": "(b) center passes Holm over 9 positions", "units": int(allrf.center_holm_b.sum())},
{"step": "(b) ... and mean above the spontaneous windows", "units": int((allrf.center_holm_b & allrf.center_above_b).sum())},
{"step": "(b) ... and no other position passes (= exclusive, reading b)", "units": int(allrf.excl_b.sum())}])
print(steps.to_string(index=False))
for rd in ["a", "b"]:
x = allrf[(allrf.best == 0) & (allrf[f"n_holm_{rd}"] >= 1)] if rd == "a" else allrf[allrf.center_holm_b & allrf.center_above_b]
among = "units with the center best and at least one position passing" if rd == "a" else "units whose center passes, above the spontaneous mean"
print(f"Reading ({rd}), number of positions passing the Holm step among {among}:",
x[f"n_holm_{rd}"].value_counts().sort_index().to_dict())
subject V1_RS best_center excl_center_a excl_center_b kw_sig
619293 89 16 0 1 3
619296 168 89 7 8 59
620333 114 9 0 0 70
620334 68 19 0 3 35
625545 190 39 4 4 106
625554 179 85 11 12 152
625555 314 52 2 4 62
630506 58 9 0 0 38
631510 129 65 6 6 70
631570 142 35 5 3 85
633229 181 42 7 10 115
637484 172 14 2 1 109
V1 RS units 1804; exclusively center-responsive: (a) 44, (b) 52; paper 18 (tolerance 1.8)
variant, center best and p_center (a) < 0.05/9: 231
variant, center best and p_center (b) < 0.05/9: 273
variant, (a) and (b) both: 5
variant, (a), units not labelled noise: 37
Steps before the exclusivity rule (all V1 RS units):
step units
center is the position of largest mean response 474
(a) center best and its signed-rank p < 0.05 (uncorrected) 310
(a) center best and at least one position passes the Holm step 383
(a) ... exactly one position passes 68
(a) ... and it is the center (= exclusive, reading a) 44
(b) center rank-sum p < 0.05 (uncorrected) 951
(b) center passes Holm over 9 positions 765
(b) ... and mean above the spontaneous windows 602
(b) ... and no other position passes (= exclusive, reading b) 52
Reading (a), number of positions passing the Holm step among units with the center best and at least one position passing: {1: 68, 2: 57, 3: 46, 4: 33, 5: 30, 6: 21, 7: 19, 8: 10, 9: 99}
Reading (b), number of positions passing the Holm step among units whose center passes, above the spontaneous mean: {1: 52, 2: 78, 3: 86, 4: 68, 5: 49, 6: 39, 7: 38, 8: 44, 9: 148}
Result: not reproduced. The V1 RS count is 1,804, as in the paper, but 44 units are exclusively center-responsive under reading (a) and 52 under reading (b), against 18 +/- 1.8. The steps above show where the counts form: under (a), 474 units have the center as best position, 383 of them have at least one position passing the Holm step, 68 have exactly one position passing, and for 44 that position is the center; under (b), 602 have the center passing the Holm step with a mean above the spontaneous windows, and 52 of them have no other position passing. The selection is far wider than the paper's under both readings.
fig1f-npx-ic-ire-responsive¶
Of the exclusively center-responsive V1 RS units, 66.7% respond substantially to the IC image and 77.8% to the real-edge (IRE) image of their preferred IRE orientation, although the IC image has no contrast in their receptive field.
Units of fig1e under each of its readings. Preferred IRE orientation: the IRE image with the strongest mean IC response; the IC image is the one with the illusory bar of that orientation. Test (a): the code's Kruskal-Wallis across blank, IC and IRE with Tukey-Kramer against blank and AUROC > 0.5; test (b): rank-sum against blank with AUROC > 0.5. Block (i): ICwcfg1, two orientations; block (ii): ICwcfg1 and ICwcfg0, four orientations, each image tested in its block.
def ic_ire_table(sel_by_subject, ic_trials=None):
"""Per selected unit and (test, block reading): IC- and IRE-responsiveness, preferred orientations.
ic_trials: optional dict subject -> dict block -> boolean mask of trials to keep (fixed gaze)."""
rows = []
for s in sessions:
idx = sel_by_subject.get(s["subject"], [])
if len(idx) == 0:
continue
pos = s["units"].index.get_indexer(idx)
for blocks in [("ICwcfg1",), ("ICwcfg1", "ICwcfg0")]:
for test in ["a", "b"]:
res = {}
for blk in blocks:
b = s["blocks"][blk]
keep = np.ones(len(b["codes"]), bool) if ic_trials is None else ic_trials[s["subject"]][blk]
R, c = b["R"][keep][:, pos], b["codes"][keep]
for k, (icc, irec) in enumerate([(il.IC1, il.IRE1), (il.IC2, il.IRE2)], start=1):
icr, irer = il.ic_ire_responsive(R, c, icc, irec, test)
res[(blk, k)] = {"ire_mean": R[c == irec].mean(axis=0), "ic_mean": R[c == icc].mean(axis=0),
"ic_resp": icr, "ire_resp": irer}
keys = list(res)
ire_means = np.stack([res[k]["ire_mean"] for k in keys], axis=1)
ic_means = np.stack([res[k]["ic_mean"] for k in keys], axis=1)
pi = ire_means.argmax(axis=1)
pc = ic_means.argmax(axis=1)
for j, unit in enumerate(idx):
kk = keys[pi[j]]
rows.append({"subject": s["subject"], "unit": unit, "blocks": "+".join(blocks), "test": test,
"pref_ire_ori": il.BAR_ORIENTATION[kk], "pref_ic_ori": il.BAR_ORIENTATION[keys[pc[j]]],
"ic_resp": bool(res[kk]["ic_resp"][j]), "ire_resp": bool(res[kk]["ire_resp"][j])})
return pd.DataFrame(rows)
def frac_table(t, n_units):
out = []
for (blocks, test), x in t.groupby(["blocks", "test"]):
ire = x[x.ire_resp]
out.append({"blocks": blocks, "test": test, "units": n_units, "IC_resp_%": 100 * x.ic_resp.mean(),
"IRE_resp_%": 100 * x.ire_resp.mean(), "IRE_resp_units": len(ire),
"ori_match_%": 100 * (ire.pref_ic_ori == ire.pref_ire_ori).mean() if len(ire) else np.nan})
return pd.DataFrame(out)
fig1f = {}
for rd in ["a", "b"]:
sel = {k: d.index[d[f"excl_{rd}"]] for k, d in rfres.items()}
n = sum(len(v) for v in sel.values())
t = ic_ire_table(sel)
fig1f[rd] = t
ft = frac_table(t, n)
ft.insert(0, "fig1e_reading", rd)
print(ft.round(1).to_string(index=False))
print("\nPaper: IC-responsive 66.7% (tolerance 6.67), IRE-responsive 77.8% (7.78); orientation match 78.6% (7.86)")
fig1e_reading blocks test units IC_resp_% IRE_resp_% IRE_resp_units ori_match_%
a ICwcfg1 a 44 50.0 63.6 28 75.0
a ICwcfg1 b 44 56.8 63.6 28 75.0
a ICwcfg1+ICwcfg0 a 44 38.6 77.3 34 61.8
a ICwcfg1+ICwcfg0 b 44 54.5 79.5 35 62.9
fig1e_reading blocks test units IC_resp_% IRE_resp_% IRE_resp_units ori_match_%
b ICwcfg1 a 52 53.8 57.7 30 93.3
b ICwcfg1 b 52 55.8 63.5 33 90.9
b ICwcfg1+ICwcfg0 a 52 38.5 65.4 34 61.8
b ICwcfg1+ICwcfg0 b 52 44.2 67.3 35 60.0
Paper: IC-responsive 66.7% (tolerance 6.67), IRE-responsive 77.8% (7.78); orientation match 78.6% (7.86)
Result: not reproduced. No combination of fig1e reading, test and block puts both fractions within tolerance. The IC-responsive fraction ranges from 38.5% to 56.8%, below 66.7 +/- 6.67 in every combination; the IRE-responsive fraction is within 77.8 +/- 7.78 only with both blocks under fig1e reading (a) (77.3% and 79.5%). These units are those of fig1e, which did not reproduce.
fig1g-npx-orientation-match¶
For IRE-responsive exclusively center-responsive V1 RS units, the preferred orientation of the illusory bar matches that of the real bar in 78.6% of units.
The ori_match_% column above gives the match for every reading; the table below adds the chance
level of each reading (1/number of orientations) and a shuffle of the preferred IC orientations
among the same units (1,000 shuffles, 95th percentile).
rng = np.random.default_rng(7)
rows = []
for rd, t in fig1f.items():
for (blocks, test), x in t.groupby(["blocks", "test"]):
ire = x[x.ire_resp]
if len(ire) == 0:
continue
match = (ire.pref_ic_ori.values == ire.pref_ire_ori.values).mean()
null = [np.mean(rng.permutation(ire.pref_ic_ori.values) == ire.pref_ire_ori.values) for _ in range(1000)]
rows.append({"fig1e_reading": rd, "blocks": blocks, "test": test, "IRE_resp_units": len(ire),
"match_%": 100 * match, "shuffle_95th_%": 100 * np.percentile(null, 95),
"chance_%": 100 / (2 if blocks == "ICwcfg1" else 4)})
print(pd.DataFrame(rows).round(1).to_string(index=False))
fig1e_reading blocks test IRE_resp_units match_% shuffle_95th_% chance_%
a ICwcfg1 a 28 75.0 67.9 50.0
a ICwcfg1 b 28 75.0 67.9 50.0
a ICwcfg1+ICwcfg0 a 34 61.8 41.2 25.0
a ICwcfg1+ICwcfg0 b 35 62.9 40.0 25.0
b ICwcfg1 a 30 93.3 66.7 50.0
b ICwcfg1 b 33 90.9 66.7 50.0
b ICwcfg1+ICwcfg0 a 34 61.8 38.2 25.0
b ICwcfg1+ICwcfg0 b 35 60.0 37.1 25.0
Result: not reproduced. The criterion judges the match under the reading that reproduces fig1f, and none does. The match ranges from 60.0% to 93.3% across readings; only fig1e (a) with block ICwcfg1 (75.0%) is within 78.6 +/- 7.86. Every reading is above its shuffle 95th percentile (37.1% to 67.9%).
ed2b-npx-fixed-gaze¶
Limiting the Neuropixels analysis to fixed-gaze trials (pupil within 8 degrees of its modal position) leaves 11 exclusively center-responsive V1 RS units, and their responses to ICs and orientation match do not change the conclusion.
The files give pupil position in camera pixels only, so reading (a) (8 visual degrees) cannot be run; reading (b), the code's 20 pixels, is used. Fixed-gaze trial: the largest distance of the pupil center from its modal position, from the first eye frame after onset to the end of the trial (1 s for RF trials, 0.4 s for IC trials), is below 20 pixels and no frame is a likely blink (code, analyzeMouseEye.m). A session is valid when its fixed-gaze RF trials cover all 9 positions (criterion; Extended Data Fig. 2 legend, "at least one trial in each position"). Which RF trials count is not stated, so two rules are run: (s) standard trials only, the trials the exclusivity test uses; (r) any RF trial, standard or inverse. fig1e, fig1f and fig1g are then repeated on fixed-gaze trials of the valid sessions under each rule; a valid session with no fixed-gaze standard trial at a position gets p = 1 there.
fix_rf, fix_ic, valid_by = {}, {}, {"s": [], "r": []}
rows = []
for s in sessions:
e = s["eye"]
if e is None:
continue
fr = (e["RF"]["maxdist"] < 20) & ~e["RF"]["blink"]
fi = (e["ICwcfg1"]["maxdist"] < 20) & ~e["ICwcfg1"]["blink"]
fix_rf[s["subject"]] = fr
fix_ic[s["subject"]] = {"ICwcfg1": fi, "ICwcfg0": (e["ICwcfg0"]["maxdist"] < 20) & ~e["ICwcfg0"]["blink"]}
cl = ~s["rf"]["inverse"]
cov_s = np.unique(s["rf"]["pos"][fr & cl])
cov_r = np.unique(s["rf"]["pos"][fr])
if len(cov_s) == 9:
valid_by["s"].append(s["subject"])
if len(cov_r) == 9:
valid_by["r"].append(s["subject"])
rows.append({"subject": s["subject"], "fixed_RF_trials_%": 100 * fr.mean(), "fixed_standard_RF_trials": int((fr & cl).sum()),
"fixed_inverse_RF_trials": int((fr & ~cl).sum()), "positions_standard": len(cov_s), "positions_any": len(cov_r),
"valid_s": len(cov_s) == 9, "valid_r": len(cov_r) == 9, "fixed_ICwcfg1_trials_%": 100 * fi.mean(),
"standard_trials_per_position": np.bincount(s["rf"]["pos"][fr & cl], minlength=9).tolist()})
print(pd.DataFrame(rows).round(1).to_string(index=False))
print(f"\nSessions with eye tracking: {len(fix_rf)} (paper 10); valid fixed-gaze RF sessions: (s) {len(valid_by['s'])}, "
f"(r) {len(valid_by['r'])} (paper 9)")
rf_fix = {}
for s in sessions:
if s["subject"] not in valid_by["r"]:
continue
u = s["units"]
m = ((u.area == "V1") & u.rs).values
d = il.rf_positions(s["rf"]["R"][:, m], s["rf"]["pos"], s["rf"]["inverse"], s["spont_rate"][m],
s["spont_windows"][:, m], trials=fix_rf[s["subject"]])
d.index = u.index[m]
rf_fix[s["subject"]] = d
rng = np.random.default_rng(11)
rows = []
for rule in ["s", "r"]:
for rd in ["a", "b"]:
sel = {k: d.index[d[f"excl_{rd}"]] for k, d in rf_fix.items() if k in valid_by[rule]}
n = sum(len(v) for v in sel.values())
print(f"session rule ({rule}), fixed-gaze exclusively center-responsive V1 RS units, fig1e reading ({rd}): {n} "
f"(paper 11, tolerance 1.1); per session {({k: len(v) for k, v in sel.items()})}")
t = ic_ire_table(sel, ic_trials=fix_ic)
for (blocks, test), x in t.groupby(["blocks", "test"]):
ire = x[x.ire_resp]
match = (ire.pref_ic_ori.values == ire.pref_ire_ori.values).mean() if len(ire) else np.nan
null = [np.mean(rng.permutation(ire.pref_ic_ori.values) == ire.pref_ire_ori.values) for _ in range(1000)] if len(ire) else [np.nan]
rows.append({"session_rule": rule, "fig1e_reading": rd, "blocks": blocks, "test": test, "units": n,
"IC_resp_%": 100 * x.ic_resp.mean(), "IRE_resp_%": 100 * x.ire_resp.mean(), "IRE_resp_units": len(ire),
"match_%": 100 * match, "shuffle_95th_%": 100 * np.percentile(null, 95)})
print(pd.DataFrame(rows).round(1).to_string(index=False))
subject fixed_RF_trials_% fixed_standard_RF_trials fixed_inverse_RF_trials positions_standard positions_any valid_s valid_r fixed_ICwcfg1_trials_% standard_trials_per_position 619296 91.7 82 83 9 9 True True 98.3 [9, 9, 9, 9, 10, 10, 8, 9, 9] 620333 28.3 29 22 9 9 True True 84.2 [1, 3, 4, 4, 4, 4, 2, 3, 4] 620334 53.3 47 49 9 9 True True 72.6 [6, 6, 4, 5, 6, 5, 5, 5, 5] 625545 85.6 77 77 9 9 True True 84.2 [9, 9, 9, 8, 6, 8, 9, 10, 9] 625554 5.6 3 7 3 7 False False 83.3 [1, 0, 0, 1, 0, 0, 0, 0, 1] 625555 58.9 53 53 9 9 True True 94.8 [5, 6, 4, 8, 9, 5, 4, 5, 7] 630506 36.7 30 36 8 9 False True 88.8 [3, 4, 5, 3, 4, 2, 0, 6, 3] 631510 58.9 53 53 9 9 True True 86.2 [6, 5, 6, 7, 8, 3, 5, 6, 7] 631570 78.3 68 73 9 9 True True 92.6 [8, 9, 7, 7, 6, 8, 9, 7, 7] 633229 36.1 29 36 9 9 True True 83.2 [4, 4, 3, 3, 1, 3, 3, 5, 3] Sessions with eye tracking: 10 (paper 10); valid fixed-gaze RF sessions: (s) 8, (r) 9 (paper 9)
session rule (s), fixed-gaze exclusively center-responsive V1 RS units, fig1e reading (a): 10 (paper 11, tolerance 1.1); per session {'619296': 6, '620333': 0, '620334': 0, '625545': 4, '625555': 0, '631510': 0, '631570': 0, '633229': 0}
session rule (s), fixed-gaze exclusively center-responsive V1 RS units, fig1e reading (b): 41 (paper 11, tolerance 1.1); per session {'619296': 7, '620333': 0, '620334': 5, '625545': 3, '625555': 2, '631510': 9, '631570': 4, '633229': 11}
session rule (r), fixed-gaze exclusively center-responsive V1 RS units, fig1e reading (a): 10 (paper 11, tolerance 1.1); per session {'619296': 6, '620333': 0, '620334': 0, '625545': 4, '625555': 0, '630506': 0, '631510': 0, '631570': 0, '633229': 0}
session rule (r), fixed-gaze exclusively center-responsive V1 RS units, fig1e reading (b): 41 (paper 11, tolerance 1.1); per session {'619296': 7, '620333': 0, '620334': 5, '625545': 3, '625555': 2, '630506': 0, '631510': 9, '631570': 4, '633229': 11}
session_rule fig1e_reading blocks test units IC_resp_% IRE_resp_% IRE_resp_units match_% shuffle_95th_%
s a ICwcfg1 a 10 50.0 80.0 8 100.0 75.0
s a ICwcfg1 b 10 60.0 90.0 9 100.0 77.8
s a ICwcfg1+ICwcfg0 a 10 50.0 70.0 7 71.4 57.1
s a ICwcfg1+ICwcfg0 b 10 60.0 80.0 8 75.0 50.0
s b ICwcfg1 a 41 46.3 65.9 27 88.9 66.7
s b ICwcfg1 b 41 48.8 73.2 30 86.7 66.7
s b ICwcfg1+ICwcfg0 a 41 41.5 78.0 32 46.9 37.5
s b ICwcfg1+ICwcfg0 b 41 51.2 80.5 33 45.5 39.4
r a ICwcfg1 a 10 50.0 80.0 8 100.0 75.0
r a ICwcfg1 b 10 60.0 90.0 9 100.0 77.8
r a ICwcfg1+ICwcfg0 a 10 50.0 70.0 7 71.4 57.1
r a ICwcfg1+ICwcfg0 b 10 60.0 80.0 8 75.0 50.0
r b ICwcfg1 a 41 46.3 65.9 27 88.9 66.7
r b ICwcfg1 b 41 48.8 73.2 30 86.7 66.7
r b ICwcfg1+ICwcfg0 a 41 41.5 78.0 32 46.9 40.6
r b ICwcfg1+ICwcfg0 b 41 51.2 80.5 33 45.5 39.4
Result: reproduced, resting on an assumption. Ten sessions have eye tracking, as in the paper. Which fixed-gaze RF trials must cover the 9 positions is not stated: counting any RF trial (rule r) gives 9 valid sessions, the paper's number, while counting standard trials only (rule s) gives 8, because session 630506 has no fixed-gaze standard trial at one position. Rule (r) is adopted, as the criterion and the legend speak of RF trials without restricting them to standard patches. Under rule (r) and fig1e reading (a), 10 exclusively center-responsive units remain (within 11 +/- 1.1), and their orientation match (71.4% to 100% across tests and blocks) is above the shuffle 95th percentile (50.0% to 77.8%) in every combination. Under fig1e reading (b), 41 units remain, outside tolerance. The units are the same 10 under rule (s), so the count rests on the session rule only through the session count. The fixed-gaze units rest on the fig1e selection, which did not reproduce at 18 units (44 with all trials), so the match of the fixed-gaze count with the paper's 11 may be partly coincidental.
ed2f-npx-rf-match¶
Receptive field positions mapped with fixed-gaze trials match those mapped with all trials for 81% of the 398 V1 RS units with significant receptive fields.
RF position: the position of largest mean response over standard RF trials. Significant RF: Kruskal-Wallis across the 9 positions at p < 0.05. Readings: significant with (a) all trials, (b) fixed-gaze trials, (c) both. Sessions: those with eye tracking whose fixed-gaze trials cover the 9 positions (the valid sessions of ed2b, under rule s and rule r); the 10 sessions with eye tracking, with the RF position taken over the positions that have fixed-gaze trials, are shown as a variant.
def rf_match(subjects):
rows = []
for s in sessions:
if s["subject"] not in subjects:
continue
u = s["units"]
m = ((u.area == "V1") & u.rs).values
allt = rfres[s["subject"]]
fx = rf_fix[s["subject"]] if s["subject"] in rf_fix else il.rf_positions(
s["rf"]["R"][:, m], s["rf"]["pos"], s["rf"]["inverse"], s["spont_rate"][m], s["spont_windows"][:, m],
trials=fix_rf[s["subject"]])
rows.append(pd.DataFrame({"subject": s["subject"], "best_all": allt.best.values, "best_fix": fx.best.values,
"sig_all": allt.p_kw_rf.values < 0.05, "sig_fix": fx.p_kw_rf.values < 0.05}))
return pd.concat(rows)
for name, subs in [("valid sessions, rule (s)", valid_by["s"]), ("valid sessions, rule (r)", valid_by["r"]),
("all sessions with eye tracking (variant)", list(fix_rf))]:
t = rf_match(subs)
print(f"{name}: {t.subject.nunique()} sessions, {len(t)} V1 RS units")
for rd, m in [("(a) significant with all trials", t.sig_all), ("(b) significant with fixed-gaze trials", t.sig_fix),
("(c) both", t.sig_all & t.sig_fix)]:
x = t[m]
print(f" {rd}: n = {len(x)}, match = {np.mean(x.best_all == x.best_fix):.3f}")
print("Paper: 0.81 (tolerance 0.081) among 398 (39.8) units")
valid sessions, rule (s): 8 sessions, 1306 V1 RS units (a) significant with all trials: n = 602, match = 0.733 (b) significant with fixed-gaze trials: n = 414, match = 0.812 (c) both: n = 374, match = 0.832 valid sessions, rule (r): 9 sessions, 1364 V1 RS units (a) significant with all trials: n = 640, match = 0.723 (b) significant with fixed-gaze trials: n = 434, match = 0.804 (c) both: n = 392, match = 0.821
all sessions with eye tracking (variant): 10 sessions, 1543 V1 RS units (a) significant with all trials: n = 792, match = 0.676 (b) significant with fixed-gaze trials: n = 434, match = 0.804 (c) both: n = 392, match = 0.821 Paper: 0.81 (tolerance 0.081) among 398 (39.8) units
Result: reproduced. On the 8 sessions of rule (s), reading (b) gives n = 414 with a match of 0.812 and reading (c) n = 374 with 0.832; on the 9 sessions of rule (r), (b) gives n = 434 with 0.804 and (c) n = 392 with 0.821. All are within 398 +/- 39.8 and 0.81 +/- 0.081. Reading (a) gives n = 602 or 640, outside tolerance.
ed2g-npx-pupil-area¶
Pupil area during the presentation does not differ across the images of the IC block.
Pupil area z-scored over the session; per trial its mean over the 0 to 0.4 s presentation (end and maximum as variants); per session and image the mean over trials; Friedman test with sessions as blocks. Readings: (a) IC1, LC1, LC2, IC2 of ICwcfg1; (b) all 22 images of ICwcfg1.
eye_sessions = [s for s in sessions if s["eye"] is not None]
names22 = il.IC_TRIAL_TYPES
rows = []
for measure in ["area_mean", "area_end", "area_max"]:
for rd, imgs in [("(a) IC1, LC1, LC2, IC2", [il.IC1, il.LC1, il.LC2, il.IC2]), ("(b) all 22 images", list(names22))]:
M = np.array([[np.nanmean(s["eye"]["ICwcfg1"][measure][s["blocks"]["ICwcfg1"]["codes"] == c]) for c in imgs]
for s in eye_sessions])
p = stats.friedmanchisquare(*M.T).pvalue
rows.append({"measure": measure, "reading": rd, "sessions": len(M), "friedman_p": p,
"range_of_image_means_z": np.ptp(M.mean(axis=0))})
print(pd.DataFrame(rows).round(4).to_string(index=False))
print("Paper: p = 0.7220, Friedman test, 10 sessions")
measure reading sessions friedman_p range_of_image_means_z area_mean (a) IC1, LC1, LC2, IC2 10 0.3561 0.0110 area_mean (b) all 22 images 10 0.8979 0.0451 area_end (a) IC1, LC1, LC2, IC2 10 0.5399 0.0075 area_end (b) all 22 images 10 0.3473 0.0511 area_max (a) IC1, LC1, LC2, IC2 10 0.9484 0.0670 area_max (b) all 22 images 10 0.0630 0.4638 Paper: p = 0.7220, Friedman test, 10 sessions
Result: reproduced. The Friedman test across sessions gives p = 0.356 for IC1, LC1, LC2 and IC2 (reading a) and p = 0.898 for all 22 images (reading b), both at or above 0.05 (paper p = 0.722). The end and maximum variants give p from 0.063 to 0.948.
ed2a-npx-gaze-confined¶
Eye gaze is confined to a small region that fits within the illusory gap, with no systematic difference in pupil position among IC1, LC1, LC2 and IC2 presentations.
Pupil center during the 0 to 0.4 s presentation of IC1, LC1, LC2 and IC2 in ICwcfg1; per session the image means of x and of y; a Friedman test with sessions as blocks (the criterion). Per-session Kruskal-Wallis tests across the four images are a variant. Confinement: the extent, in camera pixels, of the region where each session's 2D density of pupil positions (2-pixel bins, over the whole session) is at least 0.01 of its peak, a reading of the criterion's "0.01 isoline".
imgs = [il.IC1, il.LC1, il.LC2, il.IC2]
rows, kw = [], []
for coord in ["x_mean", "y_mean"]:
M = np.array([[np.nanmean(s["eye"]["ICwcfg1"][coord][s["blocks"]["ICwcfg1"]["codes"] == c]) for c in imgs]
for s in eye_sessions])
rows.append({"coordinate": coord[0], "friedman_p": stats.friedmanchisquare(*M.T).pvalue,
"max_diff_between_image_means_px": np.ptp(M - M.mean(axis=1, keepdims=True), axis=1).max()})
for s in eye_sessions:
c = s["blocks"]["ICwcfg1"]["codes"]
v = s["eye"]["ICwcfg1"][coord]
kw.append({"subject": s["subject"], "coordinate": coord[0],
"kw_p": stats.kruskal(*[v[(c == k) & np.isfinite(v)] for k in imgs]).pvalue})
print(pd.DataFrame(rows).round(4).to_string(index=False))
kw = pd.DataFrame(kw)
print("\nVariant, per-session Kruskal-Wallis across the 4 images: p < 0.05 in", int((kw.kw_p < 0.05).sum()), "of", len(kw),
"session x coordinate tests; smallest p", kw.kw_p.min().round(4))
ext = []
for s in eye_sessions:
xy = s["eye"]["density_xy"]
H, xe, ye = np.histogram2d(xy[:, 0], xy[:, 1], bins=[np.arange(-150, 151, 2)] * 2)
inside = H >= 0.01 * H.max()
xs, ys = xe[:-1][inside.any(axis=1)], ye[:-1][inside.any(axis=0)]
ext.append({"subject": s["subject"], "x_extent_px": xs.max() - xs.min() + 2, "y_extent_px": ys.max() - ys.min() + 2,
"frac_within_20px_of_mode": np.mean(np.hypot(xy[:, 0], xy[:, 1]) < 20)})
print(pd.DataFrame(ext).round(3).to_string(index=False))
coordinate friedman_p max_diff_between_image_means_px
x 0.4717 2.2437
y 0.3735 0.6590
Variant, per-session Kruskal-Wallis across the 4 images: p < 0.05 in 1 of 20 session x coordinate tests; smallest p 0.0258
subject x_extent_px y_extent_px frac_within_20px_of_mode
619296 50 32 0.974
620333 64 50 0.824
620334 92 28 0.799
625545 88 34 0.829
625554 88 26 0.860
625555 54 26 0.966
630506 124 24 0.832
631510 76 42 0.904
631570 72 26 0.934
633229 74 30 0.855
Result: reproduced (the part that can be judged). The Friedman test of the per-session image means gives p = 0.472 for x and 0.374 for y, so no systematic difference in pupil position across IC1, LC1, LC2 and IC2; the largest difference between image means in a session is 2.2 pixels. Whether the gaze region fits within the illusory gap cannot be judged: its extent (24 to 124 pixels by session) is in camera pixels and the files give no conversion to degrees.
fig2d-npx-ccg-hva-input¶
V1 IC-encoders (n = 79) are more likely than V1 segment responders (n = 526) to receive putative excitatory connections from higher visual areas, by cross-correlogram analysis.
Spike trains of RS units over the ICwcfg1 block in 1-ms bins; the raw CCG for every pair of a V1 IC-encoder or segment responder (RS, ICwcfg1) and an HVA RS unit of the same session; the CCG of the smoothed trains is subtracted; a connection from the HVA unit when the corrected CCG peak at lags 1 to 10 ms with the HVA unit leading exceeds 7 times the flank s.d. (50 to 100 ms). Readings of the criterion: smoothing (a) 25-ms boxcar, (b) 25-ms Gaussian with s.d. 5 ms; flank s.d. from (i) the raw, (ii) the corrected CCG; measure per V1 unit (x) fraction, (y) number of HVA units connected to it.
t0 = time.time()
ccg_rows = []
for s in sessions:
u = s["units"]
g = groups[(s["subject"], "ICwcfg1")]
tg = list(g.index[(g.area == "V1") & g.rs & (g.ic_enc | g.seg_resp)])
src = list(u.index[(u.area != "V1") & u.rs])
counts = il.ccg_counts([s["ccg_bins"][i] for i in tg], [s["ccg_bins"][i] for i in src], 124)
for sm in ["box", "gauss"]:
for fl in ["raw", "corrected"]:
con = il.connections(counts, 124, sm, fl)
for k, i in enumerate(tg):
ccg_rows.append({"subject": s["subject"], "unit": i, "group": "IC-encoder" if g.loc[i, "ic_enc"] else "segment responder",
"smoothing": sm, "flank": fl, "n_conn": int(con[k].sum()), "frac_conn": con[k].mean(),
"n_hva": len(src)})
ccg = pd.DataFrame(ccg_rows)
print(f"CCGs computed in {time.time() - t0:.0f} s")
x = ccg[(ccg.smoothing == "box") & (ccg.flank == "raw")]
print("Units in the analysis:", x.group.value_counts().to_dict(), "(paper 79 and 526)")
rows = []
for (sm, fl), x in ccg.groupby(["smoothing", "flank"]):
for meas in ["frac_conn", "n_conn"]:
a, b = x[x.group == "IC-encoder"][meas], x[x.group == "segment responder"][meas]
rows.append({"smoothing": sm, "flank_sd_from": fl, "measure": meas, "IC_enc_median": a.median(), "seg_median": b.median(),
"IC_enc_mean": a.mean(), "seg_mean": b.mean(), "IC_enc_frac_with_any": (a > 0).mean(),
"seg_frac_with_any": (b > 0).mean(),
"ranksum_p_tie_corrected": stats.mannwhitneyu(a, b, alternative="two-sided").pvalue,
"ranksum_p_no_tie_corr": stats.ranksums(a, b).pvalue,
"P(IC-enc > seg)": stats.mannwhitneyu(a, b).statistic / (len(a) * len(b))})
ccg_tab = pd.DataFrame(rows)
print(ccg_tab.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
CCGs computed in 35 s
Units in the analysis: {'segment responder': 526, 'IC-encoder': 79} (paper 79 and 526)
smoothing flank_sd_from measure IC_enc_median seg_median IC_enc_mean seg_mean IC_enc_frac_with_any seg_frac_with_any ranksum_p_tie_corrected ranksum_p_no_tie_corr P(IC-enc > seg)
box corrected frac_conn 0.002681 0.001765 0.007114 0.003452 0.7848 0.6464 0.0003913 0.0005049 0.6213
box corrected n_conn 2 1 5.152 2.16 0.7848 0.6464 0.0003357 0.0005023 0.6213
box raw frac_conn 0.002198 0.0012 0.00429 0.001985 0.6835 0.5038 0.0001667 0.0003704 0.6241
box raw n_conn 1 1 3.114 1.228 0.6835 0.5038 9.825e-05 0.0002636 0.6272
gauss corrected frac_conn 0.002681 0.001536 0.006278 0.002638 0.7595 0.5913 5.236e-05 8.643e-05 0.6369
gauss corrected n_conn 2 1 4.557 1.622 0.7595 0.5913 4.351e-05 8.68e-05 0.6368
gauss raw frac_conn 0.002198 0 0.003889 0.001639 0.6835 0.4753 1.646e-05 5.487e-05 0.6406
gauss raw n_conn 1 0 2.81 1.008 0.6835 0.4753 1.009e-05 4.241e-05 0.6427
x = ccg[(ccg.smoothing == "box") & (ccg.flank == "raw")]
print("Per session, reading (a)(i), mean fraction of HVA units connected:")
print(x.pivot_table(index="subject", columns="group", values="frac_conn", aggfunc=["mean", "count"]).round(4).to_string())
fig, ax = plt.subplots(figsize=(4, 3))
ax.boxplot([x[x.group == "IC-encoder"].frac_conn, x[x.group == "segment responder"].frac_conn], showfliers=False)
ax.set_xticks([1, 2], ["IC-encoders", "segment responders"])
ax.set_ylabel("fraction of HVA RS units connected")
plt.tight_layout(); plt.show()
print("plotted range (no fliers); all values from", x.frac_conn.min(), "to", round(x.frac_conn.max(), 4))
Per session, reading (a)(i), mean fraction of HVA units connected:
mean count
group IC-encoder segment responder IC-encoder segment responder
subject
619293 NaN 0.0038 NaN 2.0
619296 0.0004 0.0014 7.0 56.0
620333 0.0006 0.0012 6.0 39.0
620334 0.0011 0.0012 2.0 26.0
625545 0.0037 0.0028 15.0 77.0
625554 0.0036 0.0010 9.0 89.0
625555 0.0151 0.0048 9.0 24.0
630506 0.0028 0.0009 3.0 29.0
631510 0.0017 0.0015 8.0 36.0
631570 0.0028 0.0035 4.0 59.0
633229 0.0055 0.0021 10.0 55.0
637484 0.0030 0.0018 6.0 34.0
plotted range (no fliers); all values from 0.0 to 0.0501
# The criterion as written: under every combination, a higher median for IC-encoders and rank-sum p < 0.05.
lit = ccg_tab[["smoothing", "flank_sd_from", "measure", "IC_enc_median", "seg_median", "ranksum_p_tie_corrected"]].copy()
lit["higher_median"] = lit.IC_enc_median > lit.seg_median
lit["p_below_0.05"] = lit.ranksum_p_tie_corrected < 0.05
lit["meets_rule"] = lit.higher_median & lit["p_below_0.05"]
print(lit.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
n = int(lit.meets_rule.sum())
print(f"Combinations meeting the rule: {n} of {len(lit)}; literal verdict:",
"reproduced" if n == len(lit) else "not reproduced" if n == 0 else "readings differ (blocked, methods_underspecified)")
smoothing flank_sd_from measure IC_enc_median seg_median ranksum_p_tie_corrected higher_median p_below_0.05 meets_rule
box corrected frac_conn 0.002681 0.001765 0.0003913 True True True
box corrected n_conn 2 1 0.0003357 True True True
box raw frac_conn 0.002198 0.0012 0.0001667 True True True
box raw n_conn 1 1 9.825e-05 False True False
gauss corrected frac_conn 0.002681 0.001536 5.236e-05 True True True
gauss corrected n_conn 2 1 4.351e-05 True True True
gauss raw frac_conn 0.002198 0 1.646e-05 True True True
gauss raw n_conn 1 0 1.009e-05 True True True
Combinations meeting the rule: 7 of 8; literal verdict: readings differ (blocked, methods_underspecified)
Result: blocked (methods_underspecified). The groups hold 79 IC-encoders and 526 segment responders, the paper's numbers. In all 8 combinations of smoothing, flank and measure, IC-encoders have the larger mean and P(IC-encoder > segment responder) of 0.62 to 0.64, with a tie-corrected rank-sum p from 1.0e-05 to 3.9e-04. The medians are higher for IC-encoders in 7 combinations and equal (1 and 1) for the number of connections with boxcar smoothing and the raw flank. The criterion requires a higher median under every combination, and the cell above shows that 7 of 8 combinations meet it. Equal medians are not a higher median, so the readings differ, which the criterion sends to blocked (methods_underspecified). The effect the claim describes is present under every combination.
Decoding¶
The decoding claims (fig3a, fig3c, fig3d, ed7b, ed7a) are in decoding_shin_2025.ipynb, and the
checks of the decoder over cross-validation seeds, its C bound and the other decoder types
(supp3-decoder-types) in decoding_checks_shin_2025.ipynb. Each notebook loads the data again, so
that each fits within the notebook time limit.
ed6ef-npx-ic-not-greater¶
In the Neuropixels data, IC images do not evoke larger population average responses than LC images in any area, averaged over all units or over visually responsive units.
Sessions with at least 50 RS units in the area; ICwcfg1; per session the IC response averaged over (e) all RS units, (f) visually responsive RS units (Kruskal-Wallis across blank, IC1, LC1, LC2, IC2 at p < 0.05), on IC trials (IC1 and IC2) and on LC trials (LC1 and LC2); one-tailed signed-rank test of IC > LC across sessions.
paper_e = {"V1": 0.9614, "LM": 0.9805, "RL": 0.7676, "AL": 0.7217, "PM": 0.9895, "AM": 0.8398}
paper_f = {"V1": 0.9788, "LM": 0.9805, "RL": 0.7114, "AL": 0.7842, "PM": 0.9866, "AM": 0.6499}
rows = []
for s in sessions:
u = s["units"]
b = s["blocks"]["ICwcfg1"]
g = groups[(s["subject"], "ICwcfg1")]
for a in il.AREAS:
m = ((u.area == a) & u.rs).values
if m.sum() < 50:
continue
R = b["R"][:, m]
icr = R[np.isin(b["codes"], [il.IC1, il.IC2])].mean(axis=0)
lcr = R[np.isin(b["codes"], [il.LC1, il.LC2])].mean(axis=0)
vr = (g.p_kw_BK.values[m] < 0.05)
rows.append({"subject": s["subject"], "area": a, "IC_all": icr.mean(), "LC_all": lcr.mean(),
"IC_vr": icr[vr].mean(), "LC_vr": lcr[vr].mean(), "n_vr": int(vr.sum())})
pop = pd.DataFrame(rows)
out = []
for a in il.AREAS:
x = pop[pop.area == a]
out.append({"area": a, "sessions": len(x),
"e_IC_minus_LC_Hz": (x.IC_all - x.LC_all).mean(), "e_p_IC>LC": il.signrank_p(x.IC_all, x.LC_all, "greater"), "e_paper": paper_e[a],
"f_IC_minus_LC_Hz": (x.IC_vr - x.LC_vr).mean(), "f_p_IC>LC": il.signrank_p(x.IC_vr, x.LC_vr, "greater"), "f_paper": paper_f[a],
"f_units_mean": x.n_vr.mean()})
pop_tab = pd.DataFrame(out)
print(pop_tab.round(4).to_string(index=False))
print("Areas with p < 0.05 in (e):", pop_tab[pop_tab["e_p_IC>LC"] < 0.05].area.tolist(), "; in (f):", pop_tab[pop_tab["f_p_IC>LC"] < 0.05].area.tolist())
area sessions e_IC_minus_LC_Hz e_p_IC>LC e_paper f_IC_minus_LC_Hz f_p_IC>LC f_paper f_units_mean V1 12 -0.0544 0.9614 0.9614 -0.0782 0.9788 0.9788 110.4167 LM 8 -0.0613 0.9805 0.9805 -0.0688 0.9805 0.9805 104.3750 RL 11 -0.0020 0.7676 0.7676 -0.0047 0.7114 0.7114 77.5455 AL 10 -0.0055 0.7217 0.7217 -0.0130 0.7842 0.7842 88.0000 PM 12 -0.0325 0.9895 0.9895 -0.0640 0.9866 0.9866 73.6667 AM 11 -0.0201 0.8608 0.8398 0.0048 0.6499 0.6499 63.0909 Areas with p < 0.05 in (e): [] ; in (f): []
Result: reproduced. No area gives p < 0.05 in either panel: the one-tailed p for IC > LC ranges from 0.72 to 0.99 over all RS units and from 0.65 to 0.99 over visually responsive units. Eleven of the 12 p-values equal the paper's to four decimals; AM in panel (e) gives 0.8608 against 0.8398. The mean IC minus LC difference is between -0.078 and 0.005 spikes/s.
Claims on the 2p imaging, 2p mesoscope and 2p holography data¶
These 21 claims (fig1ik, fig2b, fig2c, ed4c, ed5cd, fig3e (two), ed7b-2p, ed5e, ed4hi, ed3, fig4a,
fig4b, fig5c, fig5, fig5d, fig5e, ed9b, ed10, fig6d, fig6e) rest on recordings that are not in this
dandiset; the paper names no public archive for them (its Data availability lists only this
dandiset for the Neuropixels data and "additional postprocessed data" in a Code Ocean capsule,
10.24433/CO.6659254.v2, which does not say what it holds and could not be reached from this run).
They are blocked with data_missing.
print("Files in the pinned version:", len(assets), "; kinds:", assets.kind.value_counts().to_dict())
print("Any file name mentioning ophys, imaging or 2p:", any(k in p for p in assets.path for k in ["ophys", "imaging", "2p"]))
print(f"Notebook ran in {(time.time() - started) / 60:.1f} min")
Files in the pinned version: 78 ; kinds: {'probe LFP (_ecephys)': 66, 'session (_ogen)': 12}
Any file name mentioning ophys, imaging or 2p: False
Notebook ran in 11.9 min
Summary¶
| Claim | Outcome |
|---|---|
| results-npx-units-per-session | reproduced |
| methods-npx-ic-encoder-count | reproduced |
| ed6b-npx-v1-group-counts | reproduced (assumption: the files' units are the units not classified as noise) |
| ed6cd-npx-ic-encoders-all-areas | reproduced (Neuropixels part) |
| fig1e-npx-exclusive-center | not reproduced (44 or 52 units against 18) |
| fig1f-npx-ic-ire-responsive | not reproduced |
| fig1g-npx-orientation-match | not reproduced |
| ed2b-npx-fixed-gaze | reproduced (assumption: any RF trial counts for position coverage) |
| ed2f-npx-rf-match | reproduced |
| ed2g-npx-pupil-area | reproduced |
| ed2a-npx-gaze-confined | reproduced (no systematic difference; the extent in degrees cannot be judged) |
| fig2d-npx-ccg-hva-input | blocked: methods_underspecified (7 of 8 combinations meet the criterion; in the eighth the medians are equal, 1 and 1; the effect is present in all 8) |
| ed6ef-npx-ic-not-greater | reproduced |
| fig3a, fig3c, fig3d, ed7b, ed7a, supp3 | see decoding_shin_2025.ipynb, decoding_checks_shin_2025.ipynb and the two seed notebooks |
| 21 claims on 2p imaging, mesoscope and holography | 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 exclusive center selection (fig1e) uses the two readings of the criterion, as the code's analyzeRFCI.m was not compared line by line; the fixed-gaze criterion uses the code's 20 camera pixels, as the files give no conversion to visual degrees, and counts any RF trial towards the coverage of the 9 positions; the cross-correlogram correction is computed from the raw counts and the per-pair normalisation is left out (it does not change the comparison of a peak with the flank of the same pair).
What was found about the dataset (attempt record, findings): the files hold units labelled noise although the Methods say none were packaged; the Methods list 14 mice where the files and the Results have 12; the pupil position is labelled in meters but is in camera pixels; the eye video is 60 Hz where the Methods say 30 Hz; the IC block tables carry image indices without names; and the 2p data are not in the dandiset.