Reproducing Chang et al. (2023) from dandiset 000232¶
Paper. Chang Y-T, Finkel EA, Xu D, O'Connor DH. Rule-based modulation of a sensorimotor transformation across cortical areas. bioRxiv, doi:10.1101/2023.08.21.554194. This attempt used the full text of preprint version 2 (2024-03-09, CC-BY), as served by eLife for its reviewed preprint 92620v2, whose source is bioRxiv version 2; bioRxiv refused the download. Table 1 and the table of coding-dimension trial types are images in that copy and were not read.
Dandiset. DANDI:000232, version 0.240510.2038. The preprint cites no DANDI version (its data statement says "available upon request"); the paper's code repository cites this version's DOI, and it is the only published version. The draft holds the same 179 files (same assets and blobs) and differs only in contributor metadata.
What the dandiset holds. 49 silicon-probe sessions of the cross-modal selection task (S1, S2, MM, ALM and AMM, one area per session; trials and spike times of the units kept after the paper's quality filter) and 130 optogenetic behaviour sessions (cross-modal task and tactile detection task).
Code read. The paper's analysis code, https://github.com/YitingChang/cross_modal_task at commit 227631ee133709df979b9914c9031cfa992c3736 (MATLAB; preprocessing, data arrays, one folder per figure). Where the Methods leave a choice open, the analysis follows the code; where they disagree, the Methods' choice is the main result and the code's a variant. Nothing from outside the dandiset is used as data.
| Claim | Tested |
|---|---|
| fig1e-recording-counts | tested |
| fig6-session-counts | tested |
| fig1d-performance | tested |
| fig1d-outcome-fractions | tested |
| fig1h-thit-tcr-selective | tested |
| fig2c-rule-prestim-units | tested |
| fig2d-stimulus-prestim-units | tested |
| fig2e-stimulus-poststim-units | tested |
| fig2f-prestim-predicts-touch-response | tested |
| fig3b-rule-decoding | tested |
| fig3c-rule-decoding-motor-stronger | tested |
| fig3de-stimulus-decoding-controls | tested |
| figs3-svm-rf-decoding | tested |
| fig4c-prestim-distance-predicts-divergence | tested |
| fig4e-subspace-overlap | tested |
| fig4f-prestim-overlap-predicts-poststim | tested |
| figs4ef-coding-dimensions | tested |
| figs4h-choice-right-left | tested |
| fig5b-first-hit-switches | tested |
| fig5de-motor-states-track-switch | tested |
| figs5-sensory-states-no-tracking | tested |
| fig6c-pre-inhibition-tactile | tested |
| figs6b-pre-inhibition-hit-rate | tested |
| fig6d-post-inhibition-tactile | tested |
| figs7gh-visual-sensitivity | tested |
| fig6-ppc-inhibition | tested |
| fig6fg-detection-task | tested |
| fig6ij-laser-only-licks | tested |
| fig1e-probe-locations | blocked (data missing: no histology or electrode coordinates) |
In [1]:
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from sklearn.ensemble import RandomForestClassifier
from sklearn.svm import SVC
import xmodal as xm
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)
BLUE, GRAY, ORANGE = "#2a78d6", "#8a8984", "#eb6834"
started = time.time()
assets = xm.asset_urls()
files = xm.load_all(assets, spikes=True)
ephys = [f for f in files if f["kind"] == "ephys"]
opto = [f for f in files if f["kind"] == "opto"]
print(f"dandiset {xm.DANDISET_ID} version {xm.DANDISET_VERSION}: {len(assets)} NWB files, "
f"{len(ephys)} with units, {len(opto)} without; loaded in {time.time() - started:.0f} s")
dandiset 000232 version 0.240510.2038: 179 NWB files, 49 with units, 130 without; loaded in 249 s
Data check¶
Every check of the description against the values is stated under its output.
In [2]:
info = pd.DataFrame([{"path": f["path"], "subject": f["subject"], "genotype": f["genotype"], "sex": f["sex"],
"kind": f["kind"], "description": f["description"], "trials": len(f["trials"]),
"units": f.get("n_units", 0)} for f in files])
print("Files by experiment description:")
print(info.groupby("description").agg(files=("path", "size"), mice=("subject", "nunique")).to_string())
print("\nSubjects:")
print(info.groupby("subject").agg(genotype=("genotype", "first"), sex=("sex", "first"),
files=("path", "size"), ephys=("kind", lambda k: (k == "ephys").sum())).to_string())
print("\nArchive summary: 18 subjects, 179 files; files hold", info.subject.nunique(), "subjects,", len(info), "files")
Files by experiment description:
files mice
description
electrophysiological recording of ALM during the cross-modal sensory selection task 13 8
electrophysiological recording of AMM during the cross-modal sensory selection task 9 7
electrophysiological recording of MM during the cross-modal sensory selection task 9 7
electrophysiological recording of S1 during the cross-modal sensory selection task 10 6
electrophysiological recording of S2 during the cross-modal sensory selection task 8 5
optogenetic inhibition of ALM during the cross-modal sensory selection task 10 4
optogenetic inhibition of ALM during the tactile detection task 4 3
optogenetic inhibition of AMM during the cross-modal sensory selection task 21 7
optogenetic inhibition of AMM during the tactile detection task 4 3
optogenetic inhibition of MM during the cross-modal sensory selection task 11 4
optogenetic inhibition of MM during the tactile detection task 5 3
optogenetic inhibition of PPC during the cross-modal sensory selection task 21 7
optogenetic inhibition of PPC during the tactile detection task 4 3
optogenetic inhibition of S1/S2 during the cross-modal sensory selection task 10 4
optogenetic inhibition of S1/S2 during the tactile detection task 3 3
sham experiment for optogentic inhibition during the cross-modal sensory selection task 28 7
sham experiment for optogentic inhibition during the tactile detection task 9 5
Subjects:
genotype sex files ephys
subject
EF0147 Emx-Cre-Ai32 M 2 2
EF0148 Emx-Cre-Ai32 M 5 5
EF0150 Emx-Cre-Ai32 M 7 7
EF0151 Emx-Cre-Ai32 M 2 2
JL005 VGAT-ChR2 M 3 3
YT053 PV-Cre-Ai32 M 16 0
YT071 VGAT-ChR2 M 3 3
YT080 PV-Cre-Ai32 M 14 3
YT081 PV-Cre-Ai32 M 15 4
YT083 PV-Cre-Ai32 F 18 1
YT084 PV-Cre-Ai32 F 26 8
YT085 PV-Cre-Ai32 M 22 9
YT086 PV-Cre-Ai32 M 17 2
YT091 PV-Cre-Ai32 M 10 0
YT092 PV-Cre-Ai32 F 6 0
YT093 PV-Cre-Ai32 M 4 0
YT094 PV-Cre-Ai32 M 4 0
YT095 PV-Cre-Ai32 F 5 0
Archive summary: 18 subjects, 179 files; files hold 18 subjects, 179 files
The descriptions name one recorded area per ephys file and one inhibited area (or sham) and one task per optogenetic file. The 18 subjects and 179 files agree with the archive's summary.
In [3]:
rows = []
for f in ephys:
rows.append({"session": f["path"].split("/")[1][:31], "subject": f["subject"], "area": f["area"],
"group_location": f["group_location"], "electrode_locations": ",".join(f["electrode_locations"]),
"units": f["n_units"], "unit_columns": ",".join(f["unit_columns"]), "trials": len(f["trials"])})
eph_info = pd.DataFrame(rows)
print("Area in the description equals the electrode group's location in every file:",
bool((eph_info.area == eph_info.group_location).all()))
print("Electrode-level locations:", eph_info.electrode_locations.unique().tolist())
print("Units table columns:", eph_info.unit_columns.unique().tolist())
print("\nSessions, mice and units per area:")
print(eph_info.groupby("area").agg(sessions=("session", "size"), mice=("subject", "nunique"),
units=("units", "sum"), units_min=("units", "min"), units_max=("units", "max"),
trials_min=("trials", "min"), trials_max=("trials", "max")).to_string())
f0 = ephys[0]
print("\nTrials table columns (ephys):", f0["trial_columns"])
print("Files whose electrode x, y and z are all NaN:", sum(f["electrode_xyz_all_nan"] for f in ephys), "of", len(ephys))
print("First unit ids of each file start at 0:", all(f["unit_ids_head"][0] == 0 for f in ephys), "; e.g.", f0["unit_ids_head"], "...")
Area in the description equals the electrode group's location in every file: True
Electrode-level locations: ['unknown']
Units table columns: ['spike_times,obs_intervals']
Sessions, mice and units per area:
sessions mice units units_min units_max trials_min trials_max
area
ALM 13 8 256 9 35 163 368
AMM 9 7 110 2 28 160 358
MM 9 7 140 4 39 185 263
S1 10 6 177 4 28 210 315
S2 8 5 162 5 30 212 330
Trials table columns (ephys): ['start_time', 'stop_time', 'blockType', 'trialType', 'response', 'visStimType', 'somStimType', 'somOnset_times', 'somOffset_times', 'visOnset_times', 'visOffset_times', 'stimOnset_times', 'rLickOnset_times', 'rLickOffset_times', 'lLickOnset_times', 'lLickOffset_times', 'firstLick_times']
Files whose electrode x, y and z are all NaN: 49 of 49
First unit ids of each file start at 0: True ; e.g. [0, 1, 2, 3, 4] ...
In [4]:
# Trial coding: block, stimulus and response against the lick times.
rows = []
for f in ephys:
t = xm.code_trials(f["trials"])
first_side = []
for r, l, s0 in zip(t.r_licks, t.l_licks, t.stim_on):
rr, ll = r[(r - s0 > 0.1) & (r - s0 < 2)], l[(l - s0 > 0.1) & (l - s0 < 2)]
first_side.append(3 if len(rr) + len(ll) == 0 else 1 if len(ll) == 0 or (len(rr) and rr[0] < ll[0]) else 2)
first_side = np.array(first_side)
lick_pre = xm.licks_in(t, -1, 0.1)
gaps = t.start.to_numpy()[1:] - t.stop.to_numpy()[:-1]
rows.append({"session": f["path"].split("/")[1][:31], "trials": len(t), "cue_trials": int((t.b == 0).sum()),
"no_stimulus": int((t.s == 0).sum()), "response_values": ",".join(map(str, sorted(t.response.unique()))),
"response_eq_first_lick_side": round(float(np.mean(first_side[t.d > 0] == t.d[t.d > 0])), 3),
"trials_lick_-1_to_0.1s": int(lick_pre.sum()),
"stim_minus_start_min_s": round(float((t.stim_on - t.start).min()), 3),
"stop_minus_stim_min_s": round(float((t.stop - t.stim_on).min()), 3),
"gaps_over_10ms": int((gaps > 0.01).sum()), "gap_total_s": round(float(gaps[gaps > 0.01].sum()), 0)})
coding = pd.DataFrame(rows)
print(coding.describe().T[["min", "50%", "max"]].to_string())
allt = pd.concat([xm.code_trials(f["trials"]).assign(area=f["area"]) for f in ephys], ignore_index=True)
print("\nBlock type x trial type x response, all ephys sessions:")
print(pd.crosstab([allt.block_type, allt.trial_type], allt.response).to_string())
print("\nStimulus types:", allt.vis_stim_type.unique().tolist(), allt.som_stim_type.unique().tolist())
min 50% max
trials 160.000 236.000 368.000
cue_trials 2.000 4.000 6.000
no_stimulus 0.000 0.000 0.000
response_eq_first_lick_side 0.948 1.000 1.000
trials_lick_-1_to_0.1s 0.000 0.000 0.000
stim_minus_start_min_s 1.049 1.077 1.140
stop_minus_stim_min_s 3.477 3.852 4.116
gaps_over_10ms 8.000 29.000 64.000
gap_total_s 77.000 303.000 835.000
Block type x trial type x response, all ephys sessions:
response 0.0 1.0 2.0
block_type trial_type
Change Stim_Som_NoCue 0 94 0
Stim_Vis_NoCue 0 0 91
Visual Stim_Som_NoCue 2338 279 390
Stim_Vis_NoCue 722 28 2235
Whisker Stim_Som_NoCue 721 2093 56
Stim_Vis_NoCue 2151 476 306
Stimulus types: ['Amp0p5Dur0p150'] ['SineAmp0p875Freq20Cyc3', 'SineAmp0p8375Freq20Cyc3', 'SineAmp1p000Freq20Cyc3', 'SineAmp0p750Freq20Cyc3']
response agrees with the side of the first lick between 0.1 and 2 s after the stimulus in at least the
minimum fraction printed (the remainder are licks the behavioural system timed differently); 1 is a right
lick, 2 a left lick and 0 no lick. "Change" rows are the cue trials (water given on the 9th trial of a
block), left out of all analyses as in the paper's code. No trial in any ephys file has a lick between -1 and 0.1 s of the stimulus: the trials the
Methods remove for pre-stimulus licking (about 25%) are already absent, and the gaps between trials mark
where they were. The dandiset does not say this. The stimulus comes at least 1 s after the trial start and
the trial ends at least 3.4 s after it, so every analysis window lies inside a trial.
In [5]:
rows = []
for f in ephys:
t = f["trials"]
dur = float((t.stop - t.start).sum())
rates = np.array([len(s) / dur for s in f["spikes"]])
obs_is_trials = all(np.allclose(o, np.c_[t.start, t.stop], atol=2e-3) for o in f["obs_intervals"])
starts, stops = t.start.to_numpy() - 2e-3, t.stop.to_numpy() + 2e-3
in_trial = np.mean([np.mean(np.searchsorted(starts, sp, side="right") - 1 == np.searchsorted(stops, sp, side="left"))
for sp in f["spikes"]])
rows.append({"session": f["path"].split("/")[1][:31], "area": f["area"], "units": f["n_units"],
"sorted": all(np.all(np.diff(s) >= 0) for s in f["spikes"]),
"obs_intervals_are_trials": obs_is_trials, "spikes_in_trials": round(float(in_trial), 4),
"rate_min_hz": rates.min().round(2), "rate_max_hz": rates.max().round(1),
"spike_hash": hash(np.concatenate(f["spikes"][:3]).tobytes())})
units_check = pd.DataFrame(rows)
print(units_check.agg({"sorted": "all", "obs_intervals_are_trials": "all", "spikes_in_trials": "min",
"rate_min_hz": "min", "rate_max_hz": "max", "units": "sum"}).to_string())
print("Sessions whose first three spike trains are identical to another session's:",
int(units_check.spike_hash.duplicated().sum()))
sorted True obs_intervals_are_trials True spikes_in_trials 0.9995 rate_min_hz 0.16 rate_max_hz 49.1 units 845 Sessions whose first three spike trains are identical to another session's: 0
Spike times are sorted, each unit's observation intervals are the trial intervals, and spikes are stored only inside trials. No unit has a mean rate within trials above 50 Hz, consistent with the code's removal of such units, and the unit counts per area equal the paper's after its quality filter (see fig1e-recording-counts), so the Units tables hold the selected units only. They carry no quality metrics, waveforms, channel or depth, and unit ids run from 0, so the units the paper names by cluster number (for example unit 127 of YT084, 2020-02-01) cannot be found. Electrode x, y and z are NaN in all 49 files and electrode locations "unknown". No two sessions share spike trains, so no recording is deposited twice.
In [6]:
# Clock check: decoding the stimulus (tactile against visual, correct trials) in 100 ms windows around its onset.
rng_clock = np.random.default_rng(101)
windows = [(-0.3, -0.2), (-0.2, -0.1), (-0.1, 0.0), (0.0, 0.1), (0.1, 0.2)]
rows = []
for f in ephys:
t = xm.code_trials(f["trials"])
sel, block, stim = xm.correct_trials(t)
on = t.stim_on.to_numpy()[sel]
acc = {}
for lo, hi in windows:
x = xm.binned_rates(f["spikes"], on, lo, hi, hi - lo)[:, :, 0]
acc[f"{lo:+.1f}"] = xm.cv_accuracy(x, stim, rng_clock)
x = xm.binned_rates(f["spikes"], on, 0.0, 0.1, 0.1)[:, :, 0]
acc["shuffled_0.0"] = xm.cv_accuracy(x, rng_clock.permutation(stim), rng_clock)
rows.append({"session": f["path"].split("/")[1][:31], "area": f["area"], "units": f["n_units"], **acc})
clock = pd.DataFrame(rows)
print(clock.groupby("area")[[f"{lo:+.1f}" for lo, _ in windows] + ["shuffled_0.0"]].median().round(3).to_string())
clock["jump"] = clock["+0.0"] - clock["-0.1"]
print("\nSessions where accuracy 0 to 100 ms exceeds -100 to 0 ms by at least 0.05:",
int((clock.jump >= 0.05).sum()), "of", len(clock))
print(clock.loc[clock.jump < 0.05, ["session", "area", "units", "-0.1", "+0.0", "+0.1"]].round(3).to_string())
-0.3 -0.2 -0.1 +0.0 +0.1 shuffled_0.0
area
ALM 0.528 0.529 0.527 0.720 0.736 0.546
AMM 0.531 0.555 0.536 0.642 0.760 0.531
MM 0.497 0.515 0.517 0.653 0.745 0.510
S1 0.504 0.539 0.482 0.830 0.744 0.503
S2 0.506 0.541 0.510 0.887 0.774 0.562
Sessions where accuracy 0 to 100 ms exceeds -100 to 0 ms by at least 0.05: 43 of 49
session area units -0.1 +0.0 +0.1
1 sub-EF0147_ses-20190206T150628_ ALM 10 0.527 0.539 0.648
33 sub-YT084_ses-20200202T162019_b AMM 2 0.539 0.535 0.544
34 sub-YT084_ses-20200203T174908_b AMM 5 0.545 0.581 0.545
39 sub-YT085_ses-20200310T113706_b MM 4 0.528 0.565 0.745
46 sub-YT085_ses-20200322T133303_b AMM 2 0.552 0.541 0.639
48 sub-YT086_ses-20200320T171513_b MM 5 0.630 0.649 0.649
In 43 of 49 sessions the stimulus is decoded from 0 to 100 ms after the onset in the trial table, at a median accuracy of 0.64 to 0.89 by area, and not from the 100 ms windows before it (medians 0.48 to 0.56, shuffled 0.50 to 0.56), which shows that spike times and trial times share a clock. In the six sessions listed the rise in the first 100 ms is below 0.05. Three of them (EF0147 2019-02-06, YT085 2020-03-10 and 2020-03-22) rise in the next window, 100 to 200 ms (0.64 to 0.75); for the other three, YT084 2020-02-02 and 2020-02-03 (AMM, 2 and 5 units) and YT086 2020-03-20 (MM, 5 units, 0.63 before and 0.65 after), the shared clock is not established by this check.
In [7]:
# Optogenetic files: conditions, laser timing, trial labels and the code's performance check.
rows = []
for f in opto:
t = f["trials"]
rel = np.round((t.opto_on - t.stim_on).to_numpy(), 1)
dur = np.round((t.opto_off - t.opto_on).to_numpy(), 1)
licks = xm.opto_licks(t, f["task"])
rows.append({"session": f["path"].split("/")[1][:31], "mouse": f["subject"], "genotype": f["genotype"],
"task": f["task"], "condition": f["condition"], "trials": len(t),
"laser_onsets": ",".join(sorted({str(v) for v in rel[~np.isnan(rel)]})),
"laser_durations": ",".join(sorted({str(v) for v in dur[~np.isnan(dur)]})),
"laser_trials": int((~np.isnan(rel)).sum()),
"start_eq_stim": bool(np.allclose(t.start, t.stim_on)),
"overlapping_trials": int((t.start.to_numpy()[1:] < t.stop.to_numpy()[:-1]).sum()),
"passes_code_check": xm.cross_modal_passed(licks) if f["task"] == "cross-modal" else None,
"licks": licks})
optos = pd.DataFrame(rows)
print(optos.groupby(["task", "condition"]).agg(sessions=("session", "size"), mice=("mouse", "nunique"),
genotypes=("genotype", lambda g: ",".join(sorted(set(g)))),
laser_trials_median=("laser_trials", "median"),
passes_code_check=("passes_code_check", lambda p: int(p.sum()) if p.notna().all() else "not run")).to_string())
print("\nLaser onset relative to the reference time (s):", sorted(set(",".join(optos.laser_onsets).split(","))))
print("Laser durations (s):", sorted(set(",".join(optos.laser_durations).split(","))))
onsets = np.concatenate([np.round((f["trials"].opto_on - f["trials"].stim_on).to_numpy(), 1) for f in opto])
durs = np.concatenate([np.round((f["trials"].opto_off - f["trials"].opto_on).to_numpy(), 1) for f in opto])
print("Laser trials by onset:", pd.Series(onsets[~np.isnan(onsets)]).value_counts().to_dict(),
"; by duration:", pd.Series(durs[~np.isnan(durs)]).value_counts().to_dict())
print("start_time equals stimulus onset in every file:", bool(optos.start_eq_stim.all()),
"; files with overlapping trial intervals:", int((optos.overlapping_trials > 0).sum()))
ot = pd.concat([f["trials"].assign(task=f["task"]) for f in opto], ignore_index=True)
print("\nTrial types by task (rows) with somatosensory / visual onsets present:")
print(ot.groupby(["task", "trial_type"]).agg(trials=("stim_on", "size"), som=("som_on", lambda x: x.notna().sum()),
vis=("vis_on", lambda x: x.notna().sum()),
laser=("opto_on", lambda x: x.notna().sum())).to_string())
print("\nDuplicate trial tables among optogenetic files:",
int(pd.Series([hash(f["trials"].stim_on.to_numpy().tobytes()) for f in opto]).duplicated().sum()))
sessions mice genotypes laser_trials_median passes_code_check
task condition
cross-modal ALM 10 4 PV-Cre-Ai32 125.0 10
AMM 21 7 PV-Cre-Ai32 132.0 21
MM 11 4 PV-Cre-Ai32 99.0 11
PPC 21 7 PV-Cre-Ai32 127.0 21
S1/S2 10 4 PV-Cre-Ai32 124.0 10
sham 28 7 PV-Cre-Ai32 130.5 27
detection ALM 4 3 PV-Cre-Ai32 109.0 not run
AMM 4 3 PV-Cre-Ai32 111.0 not run
MM 5 3 PV-Cre-Ai32 111.0 not run
PPC 4 3 PV-Cre-Ai32 102.0 not run
S1/S2 3 3 PV-Cre-Ai32 89.0 not run
sham 9 5 PV-Cre-Ai32 93.0 not run
Laser onset relative to the reference time (s): ['-0.0', '-0.8', '-1.6', '-2.7', '0.0']
Laser durations (s): ['0.8', '2.0', '9.5']
Laser trials by onset: {0.0: 8146, -0.8: 7980, -2.7: 1, -1.6: 1} ; by duration: {2.0: 8114, 0.8: 8013, 9.5: 1}
start_time equals stimulus onset in every file: True ; files with overlapping trial intervals: 130
Trial types by task (rows) with somatosensory / visual onsets present:
trials som vis laser
task trial_type
cross-modal Stim_Som 3977 3977 0 3767
Stim_Som_NoCue 9480 9480 0 2
Stim_Vis 3306 0 3306 3147
Stim_Vis_NoCue 8163 0 8163 0
catch 6076 0 0 6076
detection Stim_Som 1398 1398 0 1398
Stim_Som_NoCue 3261 3261 0 1
Stim_Vis 816 0 0 816
Stim_Vis_NoCue 1943 0 0 0
catch 921 0 0 921
Duplicate trial tables among optogenetic files: 0
Laser onsets fall at -0.8 s (pre-stimulus inhibition, 0.8 s long) or 0 s (from the stimulus, 2 s long)
from the reference time, as the Methods state; on laser-only ("catch") trials the reference time is the
laser offset for short trains and the onset for long ones. Two laser trials start at -1.6 and -2.7 s and
one lasts 9.5 s; as in the code, any onset before 0 counts as pre-stimulus. In these files start_time equals the
stimulus onset and trial intervals overlap, so start_time and stop_time do not delimit trials. In the
detection task, trials without a whisker stimulus are labelled "Stim_Vis" though no visual stimulus is
given (no visual onset). The code's performance check passes all cross-modal sessions but one sham
session, so the files hold the sessions the paper kept (the counts match, fig6-session-counts), and the
one sham session that fails here is kept as deposited.
In [8]:
# Per-session arrays for the ephys analyses, and the code's session checks.
sessions = [xm.ephys_session(f) for f in ephys]
main = [s for s in sessions if s["area"] in xm.AREAS]
rows = []
for s in main:
t = s["trials"]
n = {cat: int(xm.is_category(t, key).sum()) for key, cat in xm.CATEGORIES.items()}
rows.append({"session": s["path"].split("/")[1][:31], "area": s["area"], **n})
counts = pd.DataFrame(rows)
needed = [c for c in xm.CATEGORIES.values() if c not in ("TTV", "VVT")]
print("Sessions with fewer than 2 trials of a category the code requires:", int((counts[needed] < 2).any(axis=1).sum()))
print("Trials per category over the 40 sessions (min, median, max):")
print(counts[list(xm.CATEGORIES.values())].agg(["min", "median", "max"]).T.to_string())
print("\nRates used: 100 ms bins", main[0]["r100"].shape[2], "; smoothed 10 ms bins", main[0]["r10"].shape[2],
"; max smoothed rate over all sessions", round(max(float(s["r10"].max()) for s in main), 1), "Hz")
Sessions with fewer than 2 trials of a category the code requires: 0
Trials per category over the 40 sessions (min, median, max):
min median max
TTT 26.0 40.0 76.0
TTV 0.0 1.0 7.0
TTN 3.0 13.5 41.0
TVT 3.0 7.0 27.0
TVV 2.0 5.5 12.0
TVN 24.0 44.5 61.0
VTT 2.0 6.0 13.0
VTV 2.0 7.0 16.0
VTN 19.0 44.0 76.0
VVT 0.0 0.0 3.0
VVV 20.0 40.0 86.0
VVN 2.0 14.5 35.0
Rates used: 100 ms bins 2 ; smoothed 10 ms bins 25 ; max smoothed rate over all sessions 318.5 Hz
fig1e-recording-counts¶
Single units were recorded in S1 (6 mice, 10 sessions, 177 neurons), S2 (5, 8, 162), MM (7, 9, 140) and ALM (8, 13, 256), one area per session, and in AMM (110 neurons).
In [9]:
paper = pd.DataFrame({"mice": [6, 5, 7, 8, None], "sessions": [10, 8, 9, 13, None], "units": [177, 162, 140, 256, 110]},
index=["S1", "S2", "MM", "ALM", "AMM"])
found = eph_info.groupby("area").agg(mice=("subject", "nunique"), sessions=("session", "size"), units=("units", "sum"))
table = paper.join(found, lsuffix="_paper", rsuffix="_dandiset")
print(table.to_string())
print("All stated counts match:", bool(all((table[c + "_paper"].isna() | (table[c + "_paper"] == table[c + "_dandiset"])).all()
for c in ("mice", "sessions", "units"))))
mice_paper sessions_paper units_paper mice_dandiset sessions_dandiset units_dandiset S1 6.0 10.0 177 6 10 177 S2 5.0 8.0 162 5 8 162 MM 7.0 9.0 140 7 9 140 ALM 8.0 13.0 256 8 13 256 AMM NaN NaN 110 7 9 110 All stated counts match: True
Result: reproduced. Sessions, mice and units per area match the paper exactly: S1 10 sessions, 6 mice, 177 units; S2 8, 5, 162; MM 9, 7, 140; ALM 13, 8, 256; AMM 110 units (9 sessions, 7 mice).
fig6-session-counts¶
Optogenetic sessions and mice per condition in the cross-modal and detection tasks (Figure 6 legend).
In [10]:
paper6 = {("cross-modal", "S1/S2"): (4, 10), ("cross-modal", "MM"): (4, 11), ("cross-modal", "ALM"): (4, 10),
("cross-modal", "sham"): (7, 28), ("cross-modal", "AMM"): (7, 21), ("cross-modal", "PPC"): (7, 21),
("detection", "sham"): (5, 9), ("detection", "S1/S2"): (3, 3), ("detection", "MM"): (3, 5),
("detection", "ALM"): (3, 4), ("detection", "AMM"): (3, 4), ("detection", "PPC"): (3, 4)}
g = optos.groupby(["task", "condition"]).agg(mice=("mouse", "nunique"), sessions=("session", "size"))
rows = [{"task": k[0], "condition": k[1], "mice_paper": v[0], "mice": g.loc[k, "mice"],
"sessions_paper": v[1], "sessions": g.loc[k, "sessions"]} for k, v in paper6.items()]
t6 = pd.DataFrame(rows)
print(t6.to_string(index=False))
print("Counts matching:", int((t6.mice == t6.mice_paper).sum() + (t6.sessions == t6.sessions_paper).sum()), "of", 2 * len(t6))
task condition mice_paper mice sessions_paper sessions cross-modal S1/S2 4 4 10 10 cross-modal MM 4 4 11 11 cross-modal ALM 4 4 10 10 cross-modal sham 7 7 28 28 cross-modal AMM 7 7 21 21 cross-modal PPC 7 7 21 21 detection sham 5 5 9 9 detection S1/S2 3 3 3 3 detection MM 3 3 5 5 detection ALM 3 3 4 4 detection AMM 3 3 4 4 detection PPC 3 3 4 4 Counts matching: 24 of 24
Result: reproduced. All 24 counts of optogenetic sessions and mice per condition match the Figure 6 legend.
fig1d-performance and fig1d-outcome-fractions¶
Performance was about 75% correct under both rules (touch 74 ± 1%, light 75 ± 1%, 12 mice), with no difference (paired t-test, p = 0.56). Outcome fractions, touch against light: hit 36/35%, correct rejection 38/40%, miss 13/14%, false alarm 13/12%.
Reading (a) uses the trials as deposited, reading (b) drops the last 20 trials of each session (Methods).
In [11]:
def block_outcomes(t):
c = t.category[t.category != ""]
out = {}
for blk, hit, cr, miss, fa in (("touch", "TTT", "TVN", "TTN", ["TTV", "TVT", "TVV"]),
("light", "VVV", "VTN", "VVN", ["VTT", "VTV", "VVT"])):
n = c.str[0].eq("T" if blk == "touch" else "V").sum()
out.update({f"{blk}_hit": (c == hit).sum() / n, f"{blk}_cr": (c == cr).sum() / n,
f"{blk}_miss": (c == miss).sum() / n, f"{blk}_fa": c.isin(fa).sum() / n})
out[f"{blk}_perf"] = out[f"{blk}_hit"] + out[f"{blk}_cr"]
return out
perf = pd.DataFrame([{"subject": s["subject"], "area": s["area"], "reading": r,
**block_outcomes(s["trials"] if r == "a" else s["trials"].iloc[:-20])}
for s in main for r in ("a", "b")])
gap_rows = []
for s in main:
t = s["trials"]
gap = t.start.to_numpy()[1:] - t.stop.to_numpy()[:-1] > 0.01
gap_rows.append({"last_20": gap[-19:].mean(), "earlier": gap[:-19].mean(), "last_row_gap_before": bool(gap[-1])})
gap_rows = pd.DataFrame(gap_rows)
print("Fraction of consecutive rows separated by removed trials: within the last 20 rows", gap_rows.last_20.mean().round(3),
"; earlier rows", gap_rows.earlier.mean().round(3), "(mean over the 40 sessions)\n")
fig1d = {}
for r, g in perf.groupby("reading"):
m = g.groupby("subject").mean(numeric_only=True)
test = stats.ttest_rel(m.touch_perf, m.light_perf)
fig1d[r] = (m, test)
print(f"Reading ({r}): {len(m)} mice; percent, mean and SEM over mice")
print(pd.DataFrame({"mean": 100 * m.mean(), "sem": 100 * m.sem(), "min_mouse": 100 * m.min(),
"max_mouse": 100 * m.max()}).round(1).T.to_string())
print(f" paired t-test touch against light performance: t = {test.statistic:.3f}, p = {test.pvalue:.3f}\n")
Fraction of consecutive rows separated by removed trials: within the last 20 rows 0.045 ; earlier rows 0.127 (mean over the 40 sessions)
Reading (a): 12 mice; percent, mean and SEM over mice
touch_hit touch_cr touch_miss touch_fa touch_perf light_hit light_cr light_miss light_fa light_perf
mean 35.6 38.1 13.0 13.3 73.6 35.2 39.6 13.6 11.6 74.8
sem 1.1 1.0 1.1 1.0 0.8 1.8 1.2 1.4 0.6 1.5
min_mouse 28.4 31.7 5.9 6.9 68.6 25.3 30.9 5.0 8.6 69.1
max_mouse 43.2 42.6 19.4 19.3 77.1 47.7 46.2 19.9 15.5 86.1
paired t-test touch against light performance: t = -0.599, p = 0.561
Reading (b): 12 mice; percent, mean and SEM over mice
touch_hit touch_cr touch_miss touch_fa touch_perf light_hit light_cr light_miss light_fa light_perf
mean 35.7 37.1 13.2 14.0 72.8 35.9 38.3 13.1 12.7 74.2
sem 1.1 1.0 1.2 1.1 0.9 1.8 1.2 1.5 0.7 1.6
min_mouse 28.3 30.9 5.9 7.4 66.2 26.6 28.9 5.6 8.1 65.0
max_mouse 43.2 41.9 20.4 19.3 76.4 47.8 45.3 20.4 17.0 84.9
paired t-test touch against light performance: t = -0.604, p = 0.558
In [12]:
expected = {"touch_perf": (74, 1), "light_perf": (75, 1), "touch_hit": (36, 1), "light_hit": (35, 2),
"touch_cr": (38, 1), "light_cr": (40, 1), "touch_miss": (13, 1), "light_miss": (14, 1),
"touch_fa": (13, 1), "light_fa": (12, 1)}
for r, (m, test) in fig1d.items():
within = {k: abs(100 * m[k].mean() - v) <= tol for k, (v, tol) in expected.items()}
print(f"Reading ({r}) within tolerance:", within)
print(f" fractions within tolerance: {sum(within[k] for k in expected if 'perf' not in k)} of 8")
print("\nPer area (reading a, mean over sessions):")
print((100 * perf[perf.reading == "a"].groupby("area")[["touch_perf", "light_perf"]].mean()).round(1).to_string())
Reading (a) within tolerance: {'touch_perf': np.True_, 'light_perf': np.True_, 'touch_hit': np.True_, 'light_hit': np.True_, 'touch_cr': np.True_, 'light_cr': np.True_, 'touch_miss': np.True_, 'light_miss': np.True_, 'touch_fa': np.True_, 'light_fa': np.True_}
fractions within tolerance: 8 of 8
Reading (b) within tolerance: {'touch_perf': np.False_, 'light_perf': np.True_, 'touch_hit': np.True_, 'light_hit': np.True_, 'touch_cr': np.True_, 'light_cr': np.False_, 'touch_miss': np.True_, 'light_miss': np.True_, 'touch_fa': np.False_, 'light_fa': np.True_}
fractions within tolerance: 6 of 8
Per area (reading a, mean over sessions):
touch_perf light_perf
area
ALM 72.3 73.7
MM 71.9 78.0
S1 73.1 73.7
S2 74.4 74.6
Result: reproduced. Under reading (a), the trials as deposited, performance is 73.6 ± 0.8% in respond-to-touch and 74.8 ± 1.5% in respond-to-light blocks (mean and SEM over 12 mice), both within 1 point of the paper, and the paired t-test gives p = 0.561 (paper 0.56). Reading (b), without the last 20 trials, gives 72.8 and 74.2% (p = 0.558); its touch value is outside tolerance. Removed trials (gaps between rows) separate 4.5% of consecutive rows within the last 20 rows of a session against 12.7% earlier, this does not tell whether 20 trials were removed after the last row.
Result: reproduced. Under reading (a) all eight outcome fractions are within the paper's SEM (touch against light: hit 35.6 against 35.2%, correct rejection 38.1 against 39.6%, miss 13.0 against 13.6%, false alarm 13.3 against 11.6%); reading (b) puts 6 of 8 within.
Single-neuron discrimination (fig1h, fig2c, fig2d, fig2e, fig2f)¶
18 to 33% of neurons discriminate tactile hits (tHit) from tactile correct rejections (tCR) in the first 150 ms (S1 28.8, S2 17.9, MM 32.9, ALM 23.4%), mostly with AUC above 0.5 (fig1h). Before the stimulus, S1 4.5, S2 2.5, MM 21.4, ALM 10.2% discriminate the rule (fig2c); almost none discriminate the stimulus (fig2d; 0, 0, 0, 0.8%); after it, 37.3, 35.2, 15, 14.1% do (fig2e). Among rule-discriminating neurons, the pre-stimulus rule preference predicts the early touch-response modulation (fig2f; MM r = 0.93, 30 neurons; ALM r = 0.83, 26 neurons).
The bootstrap resamples are shared by all units and bins of a session (the code draws them separately); each unit's interval depends only on its own resamples.
In [13]:
rng_roc = np.random.default_rng(202)
N_BOOT = 1000
roc = []
for s in main:
t = s["trials"]
hit, cr, vh, vcr = (xm.is_category(t, k) for k in (xm.T_HIT, xm.T_CR, xm.V_HIT, xm.V_CR))
n = s["n_units"]
b1h = xm.bootstrap_auc(s["r10"][hit][:, :, 10:25].reshape(hit.sum(), -1),
s["r10"][cr][:, :, 10:25].reshape(cr.sum(), -1), N_BOOT, rng_roc).reshape(N_BOOT, n, 15)
rule = xm.bootstrap_auc(s["r100"][hit | vcr][:, :, 0], s["r100"][cr | vh][:, :, 0], N_BOOT, rng_roc)
stim_pre = xm.bootstrap_auc(s["r100"][hit | cr][:, :, 0], s["r100"][vh | vcr][:, :, 0], N_BOOT, rng_roc)
stim_post = xm.bootstrap_auc(s["r100"][hit | cr][:, :, 1], s["r100"][vh | vcr][:, :, 1], N_BOOT, rng_roc)
roc.append({"area": s["area"], "subject": s["subject"], "n": n, "trials": (hit.sum(), cr.sum(), vh.sum(), vcr.sum()),
"b1h": b1h, "rule": rule, "stim_pre": stim_pre, "stim_post": stim_post})
units = []
for area in xm.AREAS:
R = [r for r in roc if r["area"] == area]
n_area = sum(r["n"] for r in R)
b1h = np.concatenate([r["b1h"] for r in R], axis=1)
sig1h = xm.bonferroni_significant(b1h.reshape(N_BOOT, -1), n_area).reshape(n_area, 15)
df = pd.DataFrame({"area": area, "subject": np.concatenate([[r["subject"]] * r["n"] for r in R]),
"selective": xm.consecutive(sig1h, 3), "selective_2bins": xm.consecutive(sig1h, 2),
"selective_4bins": xm.consecutive(sig1h, 4),
"auc_0_150": b1h.mean(axis=(0, 2)), "auc_0_50": b1h[:, :, :5].mean(axis=(0, 2))})
for name in ("rule", "stim_pre", "stim_post"):
b = np.concatenate([r[name] for r in R], axis=1)
df[name + "_sig"] = xm.bonferroni_significant(b, n_area)
df[name + "_sig_uncorrected"] = xm.bonferroni_significant(b, 1)
df[name + "_auc"] = b.mean(axis=0)
units.append(df)
units = pd.concat(units, ignore_index=True)
print("Trials per class (tHit, tCR, vHit, vCR), min over sessions:", np.min([r["trials"] for r in roc], axis=0))
Trials per class (tHit, tCR, vHit, vCR), min over sessions: [26 19 20 24]
In [14]:
summary = units.groupby("area", sort=False).agg(
neurons=("selective", "size"), fig1h_pct=("selective", "mean"), fig2c_pct=("rule_sig", "mean"),
fig2d_pct=("stim_pre_sig", "mean"), fig2e_pct=("stim_post_sig", "mean"))
for c in summary.columns[1:]:
summary[c] = (100 * summary[c]).round(1)
summary["fig1h_paper"] = [28.8, 17.9, 32.9, 23.4]
summary["fig2c_paper"] = [4.5, 2.5, 21.4, 10.2]
summary["fig2d_paper"] = [0, 0, 0, 0.8]
summary["fig2e_paper"] = [37.3, 35.2, 15.0, 14.1]
summary["fig2c_neurons"] = units.groupby("area", sort=False).rule_sig.sum()
print(summary.to_string())
sel = units[units.selective]
print(f"\nfig1h: selective neurons with mean AUC above 0.5: {(sel.auc_0_150 > 0.5).sum()} of {len(sel)} "
f"({100 * (sel.auc_0_150 > 0.5).mean():.1f}%); by area:",
sel.groupby("area", sort=False).auc_0_150.apply(lambda a: round(100 * (a > 0.5).mean(), 1)).to_dict())
print("\nVariants: percent selective with 2 or 4 consecutive bins, and fig2 percentages without Bonferroni correction")
print((100 * units.groupby("area", sort=False)[["selective_2bins", "selective_4bins", "rule_sig_uncorrected",
"stim_pre_sig_uncorrected", "stim_post_sig_uncorrected"]].mean()).round(1).to_string())
print("\nBy mouse (percent; neurons):")
print(units.groupby(["area", "subject"], sort=False).agg(n=("selective", "size"), fig1h=("selective", "mean"),
fig2c=("rule_sig", "mean"), fig2e=("stim_post_sig", "mean"))
.assign(fig1h=lambda d: (100 * d.fig1h).round(0), fig2c=lambda d: (100 * d.fig2c).round(0),
fig2e=lambda d: (100 * d.fig2e).round(0)).to_string())
neurons fig1h_pct fig2c_pct fig2d_pct fig2e_pct fig1h_paper fig2c_paper fig2d_paper fig2e_paper fig2c_neurons
area
S1 177 30.5 5.6 0.0 36.2 28.8 4.5 0.0 37.3 10
S2 162 16.7 2.5 0.0 35.8 17.9 2.5 0.0 35.2 4
MM 140 36.4 20.0 0.0 13.6 32.9 21.4 0.0 15.0 28
ALM 256 23.0 9.8 0.8 14.5 23.4 10.2 0.8 14.1 25
fig1h: selective neurons with mean AUC above 0.5: 165 of 191 (86.4%); by area: {'S1': 92.6, 'S2': 100.0, 'MM': 74.5, 'ALM': 84.7}
Variants: percent selective with 2 or 4 consecutive bins, and fig2 percentages without Bonferroni correction
selective_2bins selective_4bins rule_sig_uncorrected stim_pre_sig_uncorrected stim_post_sig_uncorrected
area
S1 36.2 22.6 16.9 4.0 53.7
S2 26.5 15.4 17.3 4.9 44.4
MM 39.3 30.7 37.9 3.6 32.9
ALM 30.9 18.0 26.2 5.9 28.9
By mouse (percent; neurons):
n fig1h fig2c fig2e
area subject
S1 EF0150 8 62.0 25.0 62.0
YT080 22 32.0 5.0 23.0
YT081 55 27.0 4.0 20.0
YT084 49 29.0 4.0 51.0
YT085 30 23.0 3.0 47.0
YT086 13 46.0 15.0 31.0
S2 EF0148 41 2.0 0.0 59.0
EF0150 34 29.0 6.0 47.0
YT080 5 40.0 20.0 60.0
YT084 25 28.0 0.0 40.0
YT085 57 12.0 2.0 9.0
MM JL005 13 38.0 0.0 31.0
YT071 44 45.0 25.0 16.0
YT081 21 33.0 24.0 10.0
YT083 5 20.0 0.0 20.0
YT084 39 36.0 23.0 10.0
YT085 13 8.0 15.0 0.0
YT086 5 60.0 20.0 20.0
ALM EF0147 22 41.0 18.0 27.0
EF0148 32 19.0 12.0 9.0
EF0150 26 35.0 31.0 8.0
EF0151 29 17.0 10.0 10.0
JL005 67 19.0 7.0 16.0
YT081 13 8.0 0.0 8.0
YT084 50 22.0 2.0 16.0
YT085 17 29.0 0.0 18.0
In [15]:
# How much the percentages move with the bootstrap draws: the same analysis with four other seeds.
reps = []
for seed in (11, 12, 13, 14):
rng_rep = np.random.default_rng(seed)
per_area = {a: {"b1h": [], "rule": [], "stim_pre": [], "stim_post": []} for a in xm.AREAS}
for s in main:
t = s["trials"]
hit, cr, vh, vcr = (xm.is_category(t, k) for k in (xm.T_HIT, xm.T_CR, xm.V_HIT, xm.V_CR))
n = s["n_units"]
d = per_area[s["area"]]
d["b1h"].append(xm.bootstrap_auc(s["r10"][hit][:, :, 10:25].reshape(hit.sum(), -1),
s["r10"][cr][:, :, 10:25].reshape(cr.sum(), -1), N_BOOT, rng_rep).reshape(N_BOOT, n, 15))
d["rule"].append(xm.bootstrap_auc(s["r100"][hit | vcr][:, :, 0], s["r100"][cr | vh][:, :, 0], N_BOOT, rng_rep))
d["stim_pre"].append(xm.bootstrap_auc(s["r100"][hit | cr][:, :, 0], s["r100"][vh | vcr][:, :, 0], N_BOOT, rng_rep))
d["stim_post"].append(xm.bootstrap_auc(s["r100"][hit | cr][:, :, 1], s["r100"][vh | vcr][:, :, 1], N_BOOT, rng_rep))
for a, d in per_area.items():
b1h = np.concatenate(d["b1h"], axis=1)
na = b1h.shape[1]
row = {"seed": seed, "area": a,
"fig1h": 100 * xm.consecutive(xm.bonferroni_significant(b1h.reshape(N_BOOT, -1), na).reshape(na, 15)).mean()}
for k, name in (("rule", "fig2c"), ("stim_pre", "fig2d"), ("stim_post", "fig2e")):
row[name] = 100 * xm.bonferroni_significant(np.concatenate(d[k], axis=1), na).mean()
reps.append(row)
main_seed = pd.DataFrame({"seed": 202, "area": units.area, "fig1h": 100 * units.selective,
"fig2c": 100 * units.rule_sig, "fig2d": 100 * units.stim_pre_sig,
"fig2e": 100 * units.stim_post_sig}).groupby("area", sort=False).mean().reset_index()
reps = pd.concat([main_seed, pd.DataFrame(reps)], ignore_index=True)
print("Percent of neurons over the main seed and four other seeds:")
print(reps.pivot(index="seed", columns="area", values=["fig1h", "fig2c"]).round(1).to_string())
avg = reps.groupby("area", sort=False)[["fig1h", "fig2c", "fig2d", "fig2e"]].agg(["mean", "min", "max"])
print("\nMean, min and max over the five seeds:")
print(avg.round(1).to_string())
tol1h = {"S1": (28.8, 2.88), "S2": (17.9, 1.79), "MM": (32.9, 3.29), "ALM": (23.4, 2.34)}
tol2c = {"S1": (4.5, 0.56), "S2": (2.5, 0.62), "MM": (21.4, 2.14), "ALM": (10.2, 1.02)}
for name, tol in (("fig1h", tol1h), ("fig2c", tol2c)):
print(f"{name}, five-seed mean within tolerance:", {a: bool(abs(avg.loc[a, (name, "mean")] - v) <= d) for a, (v, d) in tol.items()},
"; seeds within tolerance per area:", {a: int((abs(reps[reps.area == a][name] - v) <= d).sum()) for a, (v, d) in tol.items()})
Percent of neurons over the main seed and four other seeds:
fig1h fig2c
area ALM MM S1 S2 ALM MM S1 S2
seed
11.0 23.0 35.0 29.4 18.5 9.0 20.0 5.1 3.1
12.0 22.7 36.4 28.2 19.1 10.9 21.4 5.1 3.1
13.0 21.9 35.0 30.5 17.9 11.3 20.0 3.4 3.1
14.0 25.0 35.0 29.9 18.5 8.6 22.1 5.6 2.5
202.0 23.0 36.4 30.5 16.7 9.8 20.0 5.6 2.5
Mean, min and max over the five seeds:
fig1h fig2c fig2d fig2e
mean min max mean min max mean min max mean min max
area
S1 29.7 28.2 30.5 5.0 3.4 5.6 0.0 0.0 0.0 35.9 35.0 36.2
S2 18.1 16.7 19.1 2.8 2.5 3.1 0.2 0.0 0.6 35.7 35.2 36.4
MM 35.6 35.0 36.4 20.7 20.0 22.1 0.6 0.0 0.7 16.1 13.6 17.1
ALM 23.1 21.9 25.0 9.9 8.6 11.3 0.9 0.8 1.2 14.1 12.9 15.2
fig1h, five-seed mean within tolerance: {'S1': True, 'S2': True, 'MM': True, 'ALM': True} ; seeds within tolerance per area: {'S1': 5, 'S2': 5, 'MM': 3, 'ALM': 5}
fig2c, five-seed mean within tolerance: {'S1': True, 'S2': True, 'MM': True, 'ALM': True} ; seeds within tolerance per area: {'S1': 0, 'S2': 5, 'MM': 5, 'ALM': 2}
In [16]:
fig, axes = plt.subplots(1, 4, figsize=(12, 2.6), sharey=True)
for ax, area in zip(axes, xm.AREAS):
u = units[units.area == area]
ax.hist(u.auc_0_150, bins=np.arange(0, 1.01, 0.05), color=GRAY, label="all")
ax.hist(u.auc_0_150[u.selective], bins=np.arange(0, 1.01, 0.05), color=BLUE, label="selective")
ax.set_title(f"{area} ({100 * u.selective.mean():.1f}% selective)")
ax.set_xlabel("AUC tHit vs tCR, 0-150 ms")
axes[0].set_ylabel("neurons")
axes[0].legend(frameon=False)
plt.tight_layout()
plt.show()
print("Mean AUC range over neurons:", units.auc_0_150.min().round(3), "to", units.auc_0_150.max().round(3), "(all inside the axes)")
Mean AUC range over neurons: 0.241 to 0.933 (all inside the axes)
Result: partially reproduced. S1 30.5%, S2 16.7% and ALM 23.0% are within tolerance; MM is 36.4%, above 32.9 ± 3.29 and above the stated range widened by 10% (up to 36.3%). 165 of 191 selective neurons (86.4%) have mean AUC above 0.5. The outcome is judged on the main seed. Over the main seed and four others, MM gives 35.0 to 36.4% and is within tolerance in 3 of 5; the five-seed means (S1 29.7, S2 18.1, MM 35.6, ALM 23.1%) are all within tolerance. With 1000 resamples the Bonferroni interval is the minimum and maximum of the resamples, so the MM excess over the tolerance depends on the draw.
Result: partially reproduced. S2 2.5%, MM 20.0% and ALM 9.8% are within tolerance, S1 is 5.6% (10 neurons) against 4.5 ± 0.56%. MM and ALM each exceed S1 and S2. The outcome is judged on the main seed. Over five seeds S1 gives 3.4 to 5.6% and is within tolerance in none (9 neurons, 5.1%, is just above it), and ALM is within tolerance in 2 of 5; the five-seed means (S1 5.0, S2 2.8, MM 20.7, ALM 9.9%) are all within tolerance.
Result: reproduced. Before the stimulus 0% of S1, S2 and MM neurons and 0.8% of ALM neurons discriminate tactile from visual trials, as in the paper.
Result: reproduced. After the stimulus S1 36.2%, S2 35.8%, MM 13.6% and ALM 14.5% of neurons discriminate the stimulus, all within tolerance of 37.3, 35.2, 15 and 14.1%.
In [17]:
rows = []
for area in xm.AREAS:
u = units[(units.area == area) & units.rule_sig]
r, p = stats.pearsonr(u.rule_auc, u.auc_0_50) if len(u) > 2 else (np.nan, np.nan)
rho, p_s = stats.spearmanr(u.rule_auc, u.auc_0_50) if len(u) > 2 else (np.nan, np.nan)
ua = units[units.area == area]
r_all, p_all = stats.pearsonr(ua.rule_auc, ua.auc_0_50)
rows.append({"area": area, "neurons": len(u), "pearson_r": r, "p": p, "spearman_rho": rho, "p_spearman": p_s,
"r_all_neurons": r_all, "p_all": p_all,
"x_range": f"{u.rule_auc.min():.2f}-{u.rule_auc.max():.2f}" if len(u) else "",
"y_range": f"{u.auc_0_50.min():.2f}-{u.auc_0_50.max():.2f}" if len(u) else ""})
f2f = pd.DataFrame(rows)
f2f["paper_r"] = [0.69, 0.91, 0.93, 0.83]
f2f["paper_neurons"] = [8, 4, 30, 26]
print(f2f.round(4).to_string(index=False))
area neurons pearson_r p spearman_rho p_spearman r_all_neurons p_all x_range y_range paper_r paper_neurons S1 10 0.6807 0.0303 0.8182 0.0038 0.4592 0.0 0.38-0.67 0.39-0.84 0.69 8 S2 4 0.9015 0.0985 1.0000 0.0000 0.4968 0.0 0.32-0.70 0.42-0.92 0.91 4 MM 28 0.9578 0.0000 0.9431 0.0000 0.8131 0.0 0.11-0.81 0.20-0.73 0.93 30 ALM 25 0.8405 0.0000 0.8585 0.0000 0.6118 0.0 0.31-0.76 0.41-0.73 0.83 26
Result: reproduced. Among rule-discriminating neurons, r = 0.958 in MM (28 neurons, p < 0.001) and 0.841 in ALM (25 neurons, p < 0.001), within tolerance of 0.93 and 0.83. S1 gives r = 0.681 (10 neurons, p = 0.030) and S2 0.902 (4 neurons, p = 0.0985).
fig3b-rule-decoding¶
The rule can be decoded from pre-stimulus population activity; median LDA accuracy, true against shuffled, S1 0.61/0.5, S2 0.62/0.53, MM 0.7/0.52, ALM 0.68/0.55.
In [18]:
rng_dec = np.random.default_rng(303)
rows = []
for s in main:
sel, block, stim = xm.correct_trials(s["trials"])
for name, y, b in (("rule", block, 0), ("stimulus_pre", stim, 0), ("stimulus_post", stim, 1)):
x = s["r100"][sel][:, :, b]
rows.append({"area": s["area"], "subject": s["subject"], "units": s["n_units"], "what": name,
"true": xm.cv_accuracy(x, y, rng_dec), "shuffled": xm.cv_accuracy(x, rng_dec.permutation(y), rng_dec)})
dec = pd.DataFrame(rows)
med = dec.groupby(["what", "area"], sort=False)[["true", "shuffled"]].median().round(3)
print("Median accuracy over sessions:")
print(med.to_string())
r = dec[dec.what == "rule"]
print("\nRule decoding per session: range", r.true.min().round(3), "to", r.true.max().round(3),
"; sessions above their shuffled value:", int((r.true > r.shuffled).sum()), "of", len(r))
print(r.groupby("area", sort=False).apply(lambda g: pd.Series({"sessions": len(g), "min": g.true.min(), "max": g.true.max(),
"wilcoxon_p_vs_shuffled": stats.wilcoxon(g.true, g.shuffled).pvalue})).round(3).to_string())
bal = pd.DataFrame([{"area": s["area"], "larger_class": max(np.mean(xm.correct_trials(s["trials"])[1] == v)
for v in np.unique(xm.correct_trials(s["trials"])[1]))} for s in main])
print("\nShare of the larger block among correct trials, median (min to max) per area, next to the shuffled median:")
print(bal.groupby("area", sort=False).larger_class.agg(["median", "min", "max"]).assign(
shuffled_median=r.groupby("area", sort=False).shuffled.median()).round(3).to_string())
print("\nBy mouse (median rule accuracy):")
print(r.groupby(["area", "subject"], sort=False).true.median().round(3).to_string())
Median accuracy over sessions:
true shuffled
what area
rule ALM 0.682 0.535
stimulus_pre ALM 0.535 0.514
stimulus_post ALM 0.706 0.504
rule S2 0.617 0.532
stimulus_pre S2 0.520 0.511
stimulus_post S2 0.882 0.534
rule S1 0.606 0.564
stimulus_pre S1 0.491 0.500
stimulus_post S1 0.839 0.529
rule MM 0.708 0.526
stimulus_pre MM 0.500 0.503
stimulus_post MM 0.663 0.513
Rule decoding per session: range 0.49 to 0.853 ; sessions above their shuffled value: 37 of 40
sessions min max wilcoxon_p_vs_shuffled
area
ALM 13.0 0.533 0.801 0.000
S2 8.0 0.524 0.721 0.008
S1 10.0 0.490 0.740 0.055
MM 9.0 0.572 0.853 0.004
Share of the larger block among correct trials, median (min to max) per area, next to the shuffled median:
median min max shuffled_median
area
ALM 0.552 0.503 0.656 0.535
S2 0.556 0.507 0.600 0.532
S1 0.556 0.500 0.627 0.564
MM 0.529 0.509 0.578 0.526
By mouse (median rule accuracy):
area subject
ALM EF0147 0.695
S2 EF0148 0.607
ALM EF0148 0.683
S2 EF0150 0.617
S1 EF0150 0.667
ALM EF0150 0.771
EF0151 0.741
JL005 0.680
MM JL005 0.572
YT071 0.838
S2 YT080 0.721
S1 YT080 0.570
YT081 0.606
MM YT081 0.818
ALM YT081 0.560
MM YT083 0.665
YT084 0.853
S1 YT084 0.542
ALM YT084 0.593
S2 YT084 0.584
MM YT085 0.662
S1 YT085 0.649
S2 YT085 0.612
ALM YT085 0.533
S1 YT086 0.740
MM YT086 0.708
Result: reproduced. Median rule-decoding accuracy is S1 0.606, S2 0.617, MM 0.708 and ALM 0.682, all within 10% of the paper, against shuffled medians of 0.564, 0.532, 0.526 and 0.535. 37 of 40 sessions exceed their shuffled value. The shuffled medians (0.526 to 0.564) are about equal to the larger block's share of correct trials (medians 0.529 to 0.556 by area).
fig3c-rule-decoding-motor-stronger and fig3de-stimulus-decoding-controls¶
Rule decoding CIs do not overlap shuffled ones in MM and ALM and overlap in S1 and S2; the stimulus is decoded after, not before, its onset.
In [19]:
rng_boot = np.random.default_rng(304)
rows = []
for name, target, b in (("rule", "block", 0), ("stimulus_pre", "stimulus", 0), ("stimulus_post", "stimulus", 1)):
for area in xm.AREAS:
bb = xm.bootstrap_decoding([s for s in main if s["area"] == area], target, b, 200, rng_boot)
lo, hi = np.quantile(bb, [0.025, 0.975], axis=0)
rows.append({"what": name, "area": area, "true_mean": bb[:, 0].mean(), "true_lo": lo[0], "true_hi": hi[0],
"shuffled_mean": bb[:, 1].mean(), "shuffled_lo": lo[1], "shuffled_hi": hi[1],
"overlap": bool(lo[0] <= hi[1] and lo[1] <= hi[0])})
boot_dec = pd.DataFrame(rows)
print(boot_dec.round(3).to_string(index=False))
print("\nPaper, rule: S1 [0.52,0.61] vs [0.49,0.55]; S2 [0.53,0.63] vs [0.49,0.54]; MM [0.60,0.73] vs [0.48,0.53]; ALM [0.56,0.66] vs [0.49,0.54]")
what area true_mean true_lo true_hi shuffled_mean shuffled_lo shuffled_hi overlap
rule S1 0.570 0.530 0.614 0.518 0.493 0.545 True
rule S2 0.574 0.528 0.619 0.514 0.488 0.539 True
rule MM 0.667 0.601 0.731 0.499 0.472 0.523 False
rule ALM 0.610 0.560 0.655 0.517 0.495 0.544 False
stimulus_pre S1 0.495 0.464 0.527 0.503 0.479 0.529 True
stimulus_pre S2 0.509 0.473 0.550 0.513 0.487 0.542 True
stimulus_pre MM 0.499 0.464 0.536 0.502 0.476 0.527 True
stimulus_pre ALM 0.506 0.472 0.534 0.510 0.486 0.532 True
stimulus_post S1 0.783 0.728 0.831 0.502 0.477 0.527 False
stimulus_post S2 0.805 0.701 0.900 0.509 0.483 0.537 False
stimulus_post MM 0.642 0.578 0.728 0.502 0.477 0.529 False
stimulus_post ALM 0.660 0.603 0.720 0.511 0.488 0.535 False
Paper, rule: S1 [0.52,0.61] vs [0.49,0.55]; S2 [0.53,0.63] vs [0.49,0.54]; MM [0.60,0.73] vs [0.48,0.53]; ALM [0.56,0.66] vs [0.49,0.54]
In [20]:
fig, axes = plt.subplots(1, 3, figsize=(11, 2.8), sharey=True)
for ax, name in zip(axes, ("rule", "stimulus_pre", "stimulus_post")):
d = boot_dec[boot_dec.what == name].reset_index(drop=True)
x = np.arange(len(d))
ax.errorbar(x - 0.1, d.true_mean, yerr=[d.true_mean - d.true_lo, d.true_hi - d.true_mean], fmt="o", color=BLUE, label="true")
ax.errorbar(x + 0.1, d.shuffled_mean, yerr=[d.shuffled_mean - d.shuffled_lo, d.shuffled_hi - d.shuffled_mean],
fmt="o", color=GRAY, label="shuffled")
ax.axhline(0.5, color="#cccccc", lw=1, zorder=0)
ax.set_xticks(x, d.area)
ax.set_title(name.replace("_", " "))
axes[0].set_ylabel("accuracy, 95% CI")
axes[0].legend(frameon=False)
plt.tight_layout()
plt.show()
print("Plotted range:", boot_dec[["true_lo", "shuffled_lo"]].min().min().round(3), "to",
boot_dec[["true_hi", "shuffled_hi"]].max().max().round(3))
Plotted range: 0.464 to 0.9
Result: reproduced. The bootstrap CIs of true and shuffled accuracy do not overlap in MM (0.601 to 0.731 against 0.472 to 0.523) and ALM (0.560 to 0.655 against 0.495 to 0.544) and overlap in S1 (0.530 to 0.614 against 0.493 to 0.545) and S2 (0.528 to 0.619 against 0.488 to 0.539), close to the paper's intervals.
Result: reproduced. Before the stimulus the true and shuffled CIs overlap in all four areas (true means 0.495 to 0.509); after it they do not overlap in any (true means 0.642 to 0.805).
figs3-svm-rf-decoding¶
SVM and random forests decode the rule about as well as LDA (medians in the claims file).
In [21]:
rng_svm = np.random.default_rng(305)
class LinearSVM:
def fit(self, x, y):
self.model = SVC(kernel="linear", C=1.0).fit(x, y)
return self
def predict(self, x):
return self.model.predict(x)
rows = []
for s in main:
t = s["trials"]
sel, block, _ = xm.correct_trials(t)
x = s["r100"][sel][:, :, 0]
ref = s["r100"][(t.category != "").to_numpy()][:, :, 0]
xs = xm.soft_normalize(x.T, ref.T).T
xz = (x - ref.mean(0)) / (ref.std(0, ddof=1) + 0.01)
yp = rng_svm.permutation(block)
row = {"area": s["area"], "units": s["n_units"],
"svm": xm.cv_accuracy(xs, block, rng_svm, make_model=LinearSVM),
"svm_shuffled": xm.cv_accuracy(xs, yp, rng_svm, make_model=LinearSVM),
"svm_zscore": xm.cv_accuracy(xz, block, rng_svm, make_model=LinearSVM),
"svm_zscore_shuffled": xm.cv_accuracy(xz, yp, rng_svm, make_model=LinearSVM)}
for name, y in (("rf", block), ("rf_shuffled", yp)):
forest = RandomForestClassifier(500, max_features="sqrt", oob_score=True, n_jobs=4,
random_state=int(rng_svm.integers(2**31))).fit(xs, y)
row[name] = forest.oob_score_
rows.append(row)
svm = pd.DataFrame(rows)
med = svm.groupby("area").median().loc[xm.AREAS].drop(columns="units")
med["svm_paper"] = [0.60, 0.61, 0.71, 0.65]
med["rf_paper"] = [0.59, 0.60, 0.65, 0.70]
print(med.round(3).to_string())
within = sum(abs(med.loc[a, m] - med.loc[a, m + "_paper"]) <= 0.1 * med.loc[a, m + "_paper"] for a in xm.AREAS for m in ("svm", "rf"))
print("\nTrue medians within 10% of the paper:", within, "of 8; all above shuffled:",
bool(((med.svm > med.svm_shuffled) & (med.rf > med.rf_shuffled)).all()))
print("Per-session range: SVM", svm.svm.min().round(3), "-", svm.svm.max().round(3), "; RF", svm.rf.min().round(3), "-", svm.rf.max().round(3))
svm svm_shuffled svm_zscore svm_zscore_shuffled rf rf_shuffled svm_paper rf_paper area S1 0.606 0.556 0.575 0.511 0.586 0.481 0.60 0.59 S2 0.589 0.562 0.603 0.515 0.606 0.533 0.61 0.60 MM 0.658 0.545 0.701 0.539 0.676 0.545 0.71 0.65 ALM 0.648 0.535 0.642 0.537 0.655 0.503 0.65 0.70 True medians within 10% of the paper: 8 of 8; all above shuffled: True Per-session range: SVM 0.506 - 0.866 ; RF 0.455 - 0.914
Result: reproduced. All eight true medians are within 10% of the paper (SVM 0.606, 0.589, 0.658, 0.648; random forest 0.586, 0.606, 0.676, 0.655 for S1, S2, MM, ALM) and all exceed their shuffled medians.
fig4c-prestim-distance-predicts-divergence¶
Sessions with a larger pre-stimulus tHit-tCR distance show a larger distance after the stimulus (permutation test, p < 0.001, 40 sessions).
In [22]:
rng_dist = np.random.default_rng(404)
def split_test(dist, n_perm=1000):
pre = dist[:, :10].mean(1)
top = pre >= np.median(pre)
post = dist[:, 10:25]
stat = np.linalg.norm(post[top].mean(0) - post[~top].mean(0))
null = np.array([np.linalg.norm(post[p].mean(0) - post[~p].mean(0))
for p in (rng_dist.permutation(len(pre)) < top.sum() for _ in range(n_perm))])
return top, post[top].mean(), post[~top].mean(), stat, float(np.mean(null >= stat))
dist = np.array([xm.trajectory_distance(s) for s in main])
dist_code = np.array([xm.trajectory_distance(s, offset=0.0) for s in main])
for name, d in (("Methods normalization (range + 5)", dist), ("code normalization (range)", dist_code)):
top, a, b, stat, p = split_test(d)
print(f"{name}: top {top.sum()} sessions mean post distance {a:.3f}, bottom {(~top).sum()} {b:.3f}; "
f"statistic {stat:.3f}, permutation p = {p:.3f}")
areas = pd.Series([s["area"] for s in main])
print(" top-group sessions by area:", areas[top].value_counts().to_dict())
print("Distance range over sessions and bins:", dist.min().round(3), "to", dist.max().round(3))
pre, post = dist[:, :10].mean(1), dist[:, 10:25].mean(1)
print("Spearman correlation of pre and post mean distance over 40 sessions:", stats.spearmanr(pre, post))
Methods normalization (range + 5): top 20 sessions mean post distance 1.171, bottom 20 0.966; statistic 0.830, permutation p = 0.004
top-group sessions by area: {'ALM': 8, 'MM': 6, 'S1': 5, 'S2': 1}
code normalization (range): top 20 sessions mean post distance 1.428, bottom 20 1.112; statistic 1.251, permutation p = 0.001
top-group sessions by area: {'ALM': 9, 'MM': 5, 'S1': 4, 'S2': 2}
Distance range over sessions and bins: 0.279 to 2.16
Spearman correlation of pre and post mean distance over 40 sessions: SignificanceResult(statistic=np.float64(0.5369606003752346), pvalue=np.float64(0.00035391093709066643))
Result: partially reproduced. The top 20 sessions by pre-stimulus distance have a larger mean post-stimulus distance than the bottom 20 (1.171 against 0.966), with permutation p = 0.004, not below 0.001. The code's normalization gives p = 0.001.
fig4e-subspace-overlap and fig4f-prestim-overlap-predicts-poststim¶
tCR overlap minus control tHit overlap: S1 -0.23, S2 -0.23, MM -0.36, ALM -0.35 (paired t-tests); pre- and post-stimulus tCR overlaps correlate, r = 0.68, 31 sessions.
In [23]:
rng_ov = np.random.default_rng(505)
rows = []
for s in main:
o = xm.subspace_overlaps(s, rng_ov, n_splits=100)
single = [xm.subspace_overlaps(s, rng_ov, n_splits=1, ref_offset=5.0, cmp_offset=0.01, per_neuron_split=True)
for _ in range(200)]
rows.append({"area": s["area"], "units": s["n_units"], **o,
"code_control": [v["control_post"] for v in single], "code_cr": [v["cr_post"] for v in single],
"code_cr_pre": [v["cr_pre"] for v in single]})
ov = pd.DataFrame(rows)
def overlap_tests(df, ctrl="control_post", cr="cr_post"):
out = []
for area in xm.AREAS:
g = df[df.area == area]
d = g[cr] - g[ctrl]
out.append({"area": area, "sessions": len(g), "difference": d.mean(), "p_paired_t": stats.ttest_rel(g[cr], g[ctrl]).pvalue,
"sessions_negative": int((d < 0).sum())})
return pd.DataFrame(out)
main_ov = overlap_tests(ov[ov.units >= 10])
main_ov["paper"] = [-0.23, -0.23, -0.36, -0.35]
main_ov["within_10pct"] = (main_ov.difference - main_ov.paper).abs() <= 0.1 * main_ov.paper.abs()
print("Main (at least 10 units, 100 equal splits, offset 5):")
print(main_ov.round(4).to_string(index=False))
print("\nVariant, code's session rule (more than 10 units):")
print(overlap_tests(ov[ov.units > 10]).round(4).to_string(index=False))
# The code's single per-neuron split with offsets 5 and 0.01, repeated 200 times: distribution of each area's mean difference.
code_runs = []
g10 = ov[ov.units > 10]
for k in range(200):
for area in xm.AREAS:
g = g10[g10.area == area]
code_runs.append({"run": k, "area": area,
"difference": np.mean([c[k] - h[k] for c, h in zip(g.code_cr, g.code_control)]),
"p": stats.ttest_rel([c[k] for c in g.code_cr], [h[k] for h in g.code_control]).pvalue})
code_runs = pd.DataFrame(code_runs)
print("\nVariant, code's procedure (one split per run, 200 runs): difference 2.5, 50, 97.5 percentiles, and runs with p < 0.05")
print(code_runs.groupby("area", sort=False).agg(p2_5=("difference", lambda d: d.quantile(0.025)), median=("difference", "median"),
p97_5=("difference", lambda d: d.quantile(0.975)),
runs_p_below_0_05=("p", lambda p: int((p < 0.05).sum()))).round(3).to_string())
print("Overlap range:", ov[["control_post", "cr_post", "cr_pre"]].min().min().round(3), "to", ov[["control_post", "cr_post", "cr_pre"]].max().max().round(3))
Main (at least 10 units, 100 equal splits, offset 5): area sessions difference p_paired_t sessions_negative paper within_10pct S1 8 -0.2908 0.0007 8 -0.23 False S2 7 -0.2181 0.0254 6 -0.23 True MM 5 -0.4658 0.0029 5 -0.36 False ALM 12 -0.4224 0.0000 12 -0.35 False Variant, code's session rule (more than 10 units): area sessions difference p_paired_t sessions_negative S1 8 -0.2908 0.0007 8 S2 7 -0.2181 0.0254 6 MM 5 -0.4658 0.0029 5 ALM 11 -0.4198 0.0000 11
Variant, code's procedure (one split per run, 200 runs): difference 2.5, 50, 97.5 percentiles, and runs with p < 0.05
p2_5 median p97_5 runs_p_below_0_05
area
S1 -0.327 -0.286 -0.239 200
S2 -0.257 -0.222 -0.168 181
MM -0.497 -0.445 -0.376 200
ALM -0.473 -0.420 -0.359 200
Overlap range: 0.171 to 0.988
Result: partially reproduced. The tCR overlap is lower than the control tHit overlap in all four areas (S1 -0.291, p = 0.0007; S2 -0.218, p = 0.025; MM -0.466, p = 0.003; ALM -0.422, p < 0.0001), but only S2 is within 10% of the paper's difference; MM and ALM are larger in magnitude. The code's single-split procedure gives medians of -0.286, -0.222, -0.445 and -0.420 over 200 runs, and its 95% ranges exclude the paper's MM (-0.36) and ALM (-0.35) values.
In [24]:
for label, g in (("at least 10 units (main)", ov[ov.units >= 10]), ("more than 10 units (code)", ov[ov.units > 10])):
r, p = stats.pearsonr(g.cr_pre, g.cr_post)
print(f"{label}: {len(g)} sessions, Pearson r = {r:.3f}, p = {p:.2g}")
for area in xm.AREAS:
a = g[g.area == area]
print(f" {area}: {len(a)} sessions, r = {stats.pearsonr(a.cr_pre, a.cr_post)[0]:.2f}")
g = ov[ov.units > 10]
rs = [stats.pearsonr([c[k] for c in g.code_cr_pre], [c[k] for c in g.code_cr])[0] for k in range(200)]
print("Code's single-split procedure, 200 runs: r 2.5, 50, 97.5 percentiles", np.quantile(rs, [0.025, 0.5, 0.975]).round(3))
at least 10 units (main): 32 sessions, Pearson r = 0.688, p = 1.3e-05 S1: 8 sessions, r = 0.36 S2: 7 sessions, r = 0.70 MM: 5 sessions, r = 0.83 ALM: 12 sessions, r = 0.59 more than 10 units (code): 31 sessions, Pearson r = 0.686, p = 2e-05 S1: 8 sessions, r = 0.36 S2: 7 sessions, r = 0.70 MM: 5 sessions, r = 0.83 ALM: 11 sessions, r = 0.59 Code's single-split procedure, 200 runs: r 2.5, 50, 97.5 percentiles [0.573 0.653 0.734]
Result: reproduced. Over 32 sessions with at least 10 units, r = 0.688 (p = 1.3e-05); with the code's 31 sessions r = 0.686. The paper gives 0.68.
figs4ef-coding-dimensions and figs4h-choice-right-left¶
Stimulus CDs are aligned between rules in S1, S2 and MM but not ALM; choice CDs are not aligned in any area. Right-against-left choice CDs are still not well aligned.
For figs4h the criterion's "needed condition" was first read as each of the four rare categories (left
licks in touch blocks, right licks in light blocks, by stimulus); that leaves 1 of 10 S1 sessions, none in S2
and MM and 1 of 13 in ALM, so the criterion cannot decide. As the code averages over the stimulus types that exist (nanmean), the right-
and left-lick trials are pooled over stimuli and need at least 2 each; the first reading is shown as a
variant.
In [25]:
rng_cd = np.random.default_rng(606)
cd_rows, cd_boot = [], {}
for area in xm.AREAS:
ss = [s for s in main if s["area"] == area]
variants = {"main": (ss, "pooled"), "each_category": (ss, "each")}
if any(s["n_units"] <= 10 for s in ss):
variants["more_than_10_units"] = ([s for s in ss if s["n_units"] > 10], "pooled")
for vname, (sub, rule) in variants.items():
df = xm.bootstrap_coding_dimensions(sub, 1000 if vname == "main" else 300, rng_cd, rl_rule=rule)
cd_boot[(area, vname)] = df
for k in ("stim", "choice", "rl"):
if vname == "each_category" and k != "rl":
continue
tq, sq = df[k + "_true"].quantile([0.025, 0.975]), df[k + "_shuffled"].quantile([0.025, 0.975])
cd_rows.append({"area": area, "variant": vname, "sessions": len(sub), "cd": k,
"true": df[k + "_true"].mean(), "true_lo": tq.iloc[0], "true_hi": tq.iloc[1],
"shuffled": df[k + "_shuffled"].mean(), "shuffled_lo": sq.iloc[0], "shuffled_hi": sq.iloc[1],
"aligned": bool(tq.iloc[0] > sq.iloc[1] or sq.iloc[0] > tq.iloc[1]),
"resamples_with_value": int(df[k + "_true"].notna().sum())})
cd = pd.DataFrame(cd_rows)
print(cd.round(3).to_string(index=False))
rl_sessions = {area: sum(min(sum(xm.is_category(s["trials"], (b, st, d)).sum() for st in (1, 2)) for b in (1, 2) for d in (1, 2)) >= 2
for s in main if s["area"] == area) for area in xm.AREAS}
print("\nSessions with at least 2 right- and 2 left-lick trials in each block (pooled over stimuli):", rl_sessions)
rl_each = {area: sum(min(xm.is_category(s["trials"], k).sum() for k in ((1, 1, 2), (1, 2, 2), (2, 1, 1), (2, 2, 1))) >= 2
for s in main if s["area"] == area) for area in xm.AREAS}
print("Sessions with at least 2 trials in each of the four rare categories (first reading):", rl_each,
"of", {area: sum(s["area"] == area for s in main) for area in xm.AREAS})
area variant sessions cd true true_lo true_hi shuffled shuffled_lo shuffled_hi aligned resamples_with_value
S1 main 10 stim 0.557 0.397 0.714 0.252 0.137 0.384 True 1000
S1 main 10 choice 0.302 0.169 0.450 0.258 0.142 0.399 False 1000
S1 main 10 rl 0.438 0.286 0.585 0.243 0.128 0.376 False 1000
S1 each_category 10 rl 0.376 0.044 0.734 0.226 0.027 0.600 False 193
S1 more_than_10_units 8 stim 0.506 0.342 0.649 0.210 0.109 0.329 True 300
S1 more_than_10_units 8 choice 0.283 0.148 0.440 0.215 0.110 0.322 False 300
S1 more_than_10_units 8 rl 0.413 0.249 0.576 0.209 0.111 0.340 False 300
S2 main 8 stim 0.620 0.431 0.781 0.230 0.114 0.361 True 1000
S2 main 8 choice 0.249 0.132 0.399 0.230 0.121 0.360 False 1000
S2 main 8 rl 0.422 0.226 0.616 0.226 0.117 0.352 False 1000
S2 each_category 8 rl NaN NaN NaN NaN NaN NaN False 0
S2 more_than_10_units 7 stim 0.617 0.403 0.777 0.209 0.107 0.335 True 300
S2 more_than_10_units 7 choice 0.230 0.112 0.362 0.212 0.116 0.331 False 300
S2 more_than_10_units 7 rl 0.356 0.202 0.515 0.206 0.114 0.330 False 300
MM main 9 stim 0.366 0.218 0.521 0.282 0.150 0.441 False 1000
MM main 9 choice 0.304 0.165 0.465 0.286 0.154 0.448 False 1000
MM main 9 rl 0.445 0.268 0.621 0.275 0.153 0.424 False 1000
MM each_category 9 rl NaN NaN NaN NaN NaN NaN False 0
MM more_than_10_units 5 stim 0.295 0.139 0.502 0.201 0.090 0.349 False 300
MM more_than_10_units 5 choice 0.236 0.111 0.402 0.207 0.089 0.353 False 300
MM more_than_10_units 5 rl 0.516 0.287 0.704 0.205 0.091 0.353 False 300
ALM main 13 stim 0.348 0.230 0.474 0.228 0.140 0.323 False 1000
ALM main 13 choice 0.289 0.180 0.410 0.226 0.143 0.326 False 1000
ALM main 13 rl 0.307 0.196 0.424 0.223 0.138 0.321 False 1000
ALM each_category 13 rl 0.284 0.024 0.631 0.238 0.020 0.580 False 195
ALM more_than_10_units 11 stim 0.332 0.199 0.461 0.212 0.121 0.320 False 300
ALM more_than_10_units 11 choice 0.269 0.161 0.380 0.209 0.119 0.313 False 300
ALM more_than_10_units 11 rl 0.289 0.183 0.410 0.213 0.123 0.321 False 300
Sessions with at least 2 right- and 2 left-lick trials in each block (pooled over stimuli): {'S1': np.int64(10), 'S2': np.int64(8), 'MM': np.int64(9), 'ALM': np.int64(13)}
Sessions with at least 2 trials in each of the four rare categories (first reading): {'S1': np.int64(1), 'S2': np.int64(0), 'MM': np.int64(0), 'ALM': np.int64(1)} of {'S1': 10, 'S2': 8, 'MM': 9, 'ALM': 13}
In [26]:
m = cd[cd.variant == "main"].set_index(["cd", "area"])
expected_aligned = {("stim", "S1"): True, ("stim", "S2"): True, ("stim", "MM"): True, ("stim", "ALM"): False,
**{("choice", a): False for a in xm.AREAS}}
match = {k: m.loc[k, "aligned"] == v for k, v in expected_aligned.items()}
print("figs4ef judgements matching the paper:", sum(match.values()), "of 8;", {f"{k[0]} {k[1]}": bool(v) for k, v in match.items()})
paper_rl = {"S1": 0.39, "S2": 0.32, "MM": 0.35, "ALM": 0.28}
for a in xm.AREAS:
v = m.loc[("rl", a), "true"]
print(f"figs4h {a}: true mean {v:.3f} (paper {paper_rl[a]}), below 0.5: {v < 0.5}, within 0.1: {abs(v - paper_rl[a]) <= 0.1}")
figs4ef judgements matching the paper: 7 of 8; {'stim S1': True, 'stim S2': True, 'stim MM': False, 'stim ALM': True, 'choice S1': True, 'choice S2': True, 'choice MM': True, 'choice ALM': True}
figs4h S1: true mean 0.438 (paper 0.39), below 0.5: True, within 0.1: True
figs4h S2: true mean 0.422 (paper 0.32), below 0.5: True, within 0.1: False
figs4h MM: true mean 0.445 (paper 0.35), below 0.5: True, within 0.1: True
figs4h ALM: true mean 0.307 (paper 0.28), below 0.5: True, within 0.1: True
In [27]:
# figs4h S2: the bootstrap mean of the right-left dot product with four other seeds (1000 resamples each).
ss = [s for s in main if s["area"] == "S2"]
s2 = [cd_boot[("S2", "main")]["rl_true"].mean()]
for seed in (61, 62, 63, 64):
s2.append(xm.bootstrap_coding_dimensions(ss, 1000, np.random.default_rng(seed))["rl_true"].mean())
print("S2 right-left true mean, main seed then four others:", np.round(s2, 3),
"; within 0.1 of 0.32:", [bool(abs(v - 0.32) <= 0.1) for v in s2])
S2 right-left true mean, main seed then four others: [0.422 0.42 0.424 0.424 0.425] ; within 0.1 of 0.32: [False, True, False, False, False]
Result: partially reproduced. Seven of eight judgements match: stimulus CDs are aligned in S1 (0.557 against shuffled up to 0.384) and S2 (0.620 against up to 0.361) and not in ALM, and choice CDs are not aligned in any area, but the MM stimulus CD is not aligned (0.366, CI 0.218 to 0.521, against shuffled 0.150 to 0.441), where the paper reports 0.48 against 0.24.
Result: partially reproduced. Right-against-left choice CDs have true means below 0.5 in all four areas (S1 0.438, S2 0.422, MM 0.445, ALM 0.307), and S1, MM and ALM are within 0.1 of the paper; S2 is 0.422 against 0.32, more than 0.1 above it. Over the main seed and four others the S2 mean is 0.420 to 0.425 and within 0.1 of the paper in 1 of 5 (0.420); the outcome is judged on the main seed. The CIs overlap the shuffled ones in all four areas.
fig5b-first-hit-switches¶
88 switches into touch blocks and 91 into light blocks; mice switched before or immediately after the cue.
In [28]:
blocks = []
for s in main:
_, b = xm.transition_info(s["trials"])
blocks.append(b.assign(session=s["path"], area=s["area"]))
blocks = pd.concat(blocks, ignore_index=True)
nonfirst = blocks[~blocks.first_block]
with_hit = nonfirst[nonfirst.first_hit.notna()]
print("Blocks of each type, all:", blocks.groupby("block").size().to_dict(),
"; first blocks of a session:", blocks[blocks.first_block].groupby("block").size().to_dict())
print("Switches (non-first blocks) with a hit:", with_hit.groupby("block").size().to_dict(), "(paper: touch 88, light 91)")
c = with_hit[with_hit.cue.notna()].copy()
c["rel"] = c.first_hit.astype(float) - c.cue.astype(float)
c["at_or_before_cue_plus_1"] = c.rel <= 1
print("First hit relative to the cue row (rows of the file, removed trials not counted):")
print(c.groupby("block").agg(blocks_with_cue=("rel", "size"), fraction_by_cue_plus_1=("at_or_before_cue_plus_1", "mean"),
median_rel=("rel", "median"), min_rel=("rel", "min"), max_rel=("rel", "max")).round(3).to_string())
print("Blocks with a hit but no cue row in the file:", int(with_hit.cue.isna().sum()))
print("\nBlocks with a hit, first blocks included:", blocks[blocks.first_hit.notna()].groupby("block").size().to_dict(),
"; first blocks with a hit:", blocks[blocks.first_block & blocks.first_hit.notna()].groupby("block").size().to_dict())
est = []
for s in main:
_, b = xm.transition_info(s["trials"])
num = xm.estimated_trial_numbers(s["trials"])
b = b[~b.first_block & b.first_hit.notna() & b.cue.notna()]
est.append(pd.DataFrame({"block": b.block, "rel_rows": b.first_hit.astype(int) - b.cue.astype(int),
"rel_trials": num[b.first_hit.astype(int)] - num[b.cue.astype(int)]}))
est = pd.concat(est, ignore_index=True)
print("First hit relative to the cue in estimated trial numbers (removed trials counted from the gaps, xm.estimated_trial_numbers):")
print(est.groupby("block").agg(blocks=("rel_trials", "size"), fraction_by_cue_plus_1=("rel_trials", lambda x: (x <= 1).mean()),
median=("rel_trials", "median"), min=("rel_trials", "min"), max=("rel_trials", "max"),
differs_from_rows=("rel_rows", lambda x: (x != est.loc[x.index, "rel_trials"]).mean())).round(3).to_string())
print("Trials from block start to first hit (rows): median", (with_hit.first_hit.astype(float) - with_hit.start).median() + 1)
Blocks of each type, all: {'light': 91, 'touch': 88} ; first blocks of a session: {'light': 16, 'touch': 24}
Switches (non-first blocks) with a hit: {'light': 75, 'touch': 63} (paper: touch 88, light 91)
First hit relative to the cue row (rows of the file, removed trials not counted):
blocks_with_cue fraction_by_cue_plus_1 median_rel min_rel max_rel
block
light 66 0.652 1.0 -8.0 15.0
touch 56 0.571 1.0 -8.0 12.0
Blocks with a hit but no cue row in the file: 16
Blocks with a hit, first blocks included: {'light': 91, 'touch': 87} ; first blocks with a hit: {'light': 16, 'touch': 24}
First hit relative to the cue in estimated trial numbers (removed trials counted from the gaps, xm.estimated_trial_numbers):
blocks fraction_by_cue_plus_1 median min max differs_from_rows
block
light 66 0.621 1.0 -8 21 0.227
touch 56 0.554 1.0 -9 13 0.250
Trials from block start to first hit (rows): median 9.0
Result: not reproduced. Blocks after the first of a session that contain a hit, as the criterion and the code's GetFirstHitInfo define them, number 63 (touch) and 75 (light), outside tolerance of 88 and 91, so two of the three conditions fail. The third holds: more than half of the first hits come before the cue trial or on the trial after it (touch 0.571, light 0.652 in rows of the file; 0.554 and 0.621 in estimated trial numbers, which count the removed trials). The paper's numbers equal other counts: all blocks, first blocks included, are 88 and 91, and blocks with a hit, first blocks included, are 87 and 91.
fig5de-motor-states-track-switch and figs5-sensory-states-no-tracking¶
Kendall's tau of the touch-block probability over transitions: touch to light MM -0.39 [-0.67, -0.11], ALM -0.29 [-0.54, -0.04]; light to touch MM 0.37 [0.07, 0.63], ALM 0.27 [0.03, 0.49]; S1 and S2 CIs include 0.
In [29]:
rng_tr = np.random.default_rng(707)
true_tau = pd.DataFrame([{"area": s["area"], **{k: v for k, v in xm.transition_taus(s, rng_tr).items() if not k.endswith("values")}}
for s in main])
print("Mean tau over sessions, 10-fold estimate (no resampling):")
print(true_tau.groupby("area", sort=False)[["to_light", "to_touch", "pre_touch", "pre_light"]].mean().round(3).to_string())
rows = []
for area in xm.AREAS:
ss = [s for s in main if s["area"] == area]
boot = []
for b in range(500):
vals = [xm.transition_taus(ss[j], rng_tr, bootstrap=True) for j in rng_tr.integers(len(ss), size=len(ss))]
boot.append((np.nanmean([v["to_light"] for v in vals]), np.nanmean([v["to_touch"] for v in vals])))
boot = np.array(boot)
for i, name in enumerate(("to_light", "to_touch")):
lo, hi = np.nanquantile(boot[:, i], [0.025, 0.975])
rows.append({"area": area, "transition": name, "sessions": len(ss), "mean": np.nanmean(boot[:, i]), "lo": lo, "hi": hi,
"excludes_0": bool(lo > 0 or hi < 0)})
taus = pd.DataFrame(rows)
taus["paper"] = ["-0.16 [-0.43, 0.1]", "0.21 [-0.07, 0.5]", "-0.15 [-0.54, 0.21]", "0.25 [-0.17, 0.63]",
"-0.39 [-0.67, -0.11]", "0.37 [0.07, 0.63]", "-0.29 [-0.54, -0.04]", "0.27 [0.03, 0.49]"]
print("\nBootstrap (500 resamples of sessions, neurons and trials):")
print(taus.round(3).to_string(index=False))
Mean tau over sessions, 10-fold estimate (no resampling):
to_light to_touch pre_touch pre_light
area
ALM -0.487 0.436 0.597 0.391
S2 -0.333 0.583 0.537 0.425
S1 -0.267 0.333 0.551 0.439
MM -0.370 0.519 0.672 0.351
Bootstrap (500 resamples of sessions, neurons and trials): area transition sessions mean lo hi excludes_0 paper S1 to_light 10 -0.163 -0.433 0.117 False -0.16 [-0.43, 0.1] S1 to_touch 10 0.219 -0.067 0.500 False 0.21 [-0.07, 0.5] S2 to_light 8 -0.140 -0.500 0.250 False -0.15 [-0.54, 0.21] S2 to_touch 8 0.259 -0.105 0.620 False 0.25 [-0.17, 0.63] MM to_light 9 -0.401 -0.667 -0.111 True -0.39 [-0.67, -0.11] MM to_touch 9 0.382 0.074 0.630 True 0.37 [0.07, 0.63] ALM to_light 13 -0.302 -0.538 -0.051 True -0.29 [-0.54, -0.04] ALM to_touch 13 0.268 0.026 0.513 True 0.27 [0.03, 0.49]
In [30]:
n_trans = []
for s in main:
f, _ = xm.transition_info(s["trials"])
n_trans.append({"area": s["area"], "early": int(f.early.sum()), "late": int(f.late.sum()), "train": int(f.not_trans.sum())})
print("Transition trials per session (median, min):")
print(pd.DataFrame(n_trans).groupby("area", sort=False).agg(["median", "min"]).to_string())
Transition trials per session (median, min):
early late train
median min median min median min
area
ALM 6.0 2 16.0 5 197.0 146
S2 5.0 1 11.0 3 222.0 183
S1 6.5 3 20.5 11 207.0 172
MM 5.0 2 21.0 10 198.0 161
Result: reproduced. All four CIs exclude 0 with the paper's sign: touch to light MM -0.401 [-0.667, -0.111] and ALM -0.302 [-0.538, -0.051]; light to touch MM 0.382 [0.074, 0.630] and ALM 0.268 [0.026, 0.513], close to the paper's values.
Result: reproduced. All four S1 and S2 CIs include 0 (S1 -0.163 [-0.433, 0.117] and 0.219 [-0.067, 0.500]; S2 -0.140 [-0.500, 0.250] and 0.259 [-0.105, 0.620]). The 10-fold estimates without resampling are of the same size in S1 and S2 as in MM and ALM (S2 to touch 0.583, MM 0.519), so the absence rests on the wide CIs from 8 and 10 sessions.
Optogenetic inhibition (fig6c, figs6b, fig6d, figs7gh, fig6-ppc, fig6fg, fig6ij)¶
Changes are with laser minus without, per session, averaged within mouse and then over mice; 95% CIs from 10000 hierarchical resamples. "methods" resamples hit and false-alarm trials; "code" keeps the observed false-alarm rate, as the code does.
In [31]:
rng_opto = np.random.default_rng(808)
def counts_for(sub, key):
return (np.array([len(L[key]) for L in sub.licks]), np.array([int(np.sum(L[key])) for L in sub.licks]))
def change(sub, task, measure, o):
if task == "cross-modal":
hit, fa = {"tactile": ((1, 1, 1), (1, 2, 1)), "visual": ((2, 2, 2), (2, 1, 2)), "tactile_hit": ((1, 1, 1), None)}[measure]
c = {"h0": counts_for(sub, hit + (1,)), "ho": counts_for(sub, hit + (o,))}
if fa:
c.update({"f0": counts_for(sub, fa + (1,)), "fo": counts_for(sub, fa + (o,))})
else:
c = {"h0": counts_for(sub, (1, 1)), "ho": counts_for(sub, (1, o)), "f0": counts_for(sub, (2, 1)), "fo": counts_for(sub, (2, o))}
formula = (lambda p: (p["ho"] - p["fo"]) - (p["h0"] - p["f0"])) if "f0" in c else (lambda p: p["ho"] - p["h0"])
mouse = sub.mouse.to_numpy()
b = xm.hierarchical_bootstrap(mouse, c, formula, 10000, rng_opto)
bc = xm.hierarchical_bootstrap(mouse, c, formula, 10000, rng_opto, resample={"h0", "ho"})
return {"point": xm.point_estimate(mouse, c, formula), "boot_mean": b.mean(),
"lo": np.quantile(b, 0.025), "hi": np.quantile(b, 0.975),
"code_lo": np.quantile(bc, 0.025), "code_hi": np.quantile(bc, 0.975),
"sessions": len(sub), "mice": sub.mouse.nunique(), "trials_with_laser": int(c["ho"][0].sum())}
rows = []
for task, measures in (("cross-modal", ("tactile", "tactile_hit", "visual")), ("detection", ("sensitivity",))):
for measure in measures:
for o, window in ((2, "pre"), (3, "post")):
for cond in ("sham", "S1/S2", "MM", "ALM", "AMM", "PPC"):
sub = optos[(optos.task == task) & (optos.condition == cond)]
rows.append({"task": task, "measure": measure, "window": window, "condition": cond, **change(sub, task, measure, o)})
opto_res = pd.DataFrame(rows)
opto_res["ci_excludes_0"] = np.where(opto_res.hi < 0, "below", np.where(opto_res.lo > 0, "above", "includes 0"))
opto_res["code_ci_excludes_0"] = np.where(opto_res.code_hi < 0, "below", np.where(opto_res.code_lo > 0, "above", "includes 0"))
print(opto_res.round(3).to_string(index=False))
task measure window condition point boot_mean lo hi code_lo code_hi sessions mice trials_with_laser ci_excludes_0 code_ci_excludes_0 cross-modal tactile pre sham 0.011 0.010 -0.126 0.144 -0.118 0.134 28 7 297 includes 0 includes 0 cross-modal tactile pre S1/S2 -0.178 -0.179 -0.468 0.094 -0.457 0.080 10 4 98 includes 0 includes 0 cross-modal tactile pre MM -0.390 -0.388 -0.658 -0.123 -0.642 -0.135 11 4 88 below below cross-modal tactile pre ALM -0.373 -0.372 -0.707 0.005 -0.682 -0.018 10 4 73 includes 0 below cross-modal tactile pre AMM -0.062 -0.062 -0.216 0.077 -0.193 0.061 21 7 201 includes 0 includes 0 cross-modal tactile pre PPC -0.232 -0.232 -0.408 -0.052 -0.399 -0.060 21 7 174 below below cross-modal tactile post sham 0.028 0.029 -0.111 0.168 -0.099 0.159 28 7 300 includes 0 includes 0 cross-modal tactile post S1/S2 -0.389 -0.393 -0.699 -0.080 -0.677 -0.092 10 4 87 below below cross-modal tactile post MM -0.487 -0.487 -0.612 -0.327 -0.598 -0.334 11 4 100 below below cross-modal tactile post ALM -0.402 -0.402 -0.568 -0.253 -0.545 -0.274 10 4 98 below below cross-modal tactile post AMM -0.306 -0.307 -0.516 -0.111 -0.495 -0.118 21 7 196 below below cross-modal tactile post PPC -0.218 -0.218 -0.466 -0.009 -0.435 -0.034 21 7 195 below below cross-modal tactile_hit pre sham 0.022 0.021 -0.076 0.108 -0.074 0.110 28 7 297 includes 0 includes 0 cross-modal tactile_hit pre S1/S2 -0.147 -0.147 -0.303 0.010 -0.306 0.010 10 4 98 includes 0 includes 0 cross-modal tactile_hit pre MM -0.318 -0.318 -0.492 -0.152 -0.495 -0.151 11 4 88 below below cross-modal tactile_hit pre ALM -0.227 -0.227 -0.463 -0.007 -0.466 0.002 10 4 73 below includes 0 cross-modal tactile_hit pre AMM -0.048 -0.049 -0.158 0.053 -0.158 0.055 21 7 201 includes 0 includes 0 cross-modal tactile_hit pre PPC -0.205 -0.204 -0.365 -0.045 -0.359 -0.049 21 7 174 below below cross-modal tactile_hit post sham 0.050 0.050 -0.041 0.141 -0.042 0.139 28 7 300 includes 0 includes 0 cross-modal tactile_hit post S1/S2 -0.341 -0.340 -0.566 -0.098 -0.564 -0.105 10 4 87 below below cross-modal tactile_hit post MM -0.544 -0.546 -0.748 -0.296 -0.748 -0.294 11 4 100 below below cross-modal tactile_hit post ALM -0.471 -0.472 -0.714 -0.215 -0.712 -0.220 10 4 98 below below cross-modal tactile_hit post AMM -0.304 -0.304 -0.468 -0.147 -0.471 -0.153 21 7 196 below below cross-modal tactile_hit post PPC 0.038 0.037 -0.149 0.193 -0.149 0.195 21 7 195 includes 0 includes 0 cross-modal visual pre sham -0.061 -0.061 -0.186 0.063 -0.178 0.054 28 7 200 includes 0 includes 0 cross-modal visual pre S1/S2 -0.100 -0.099 -0.282 0.070 -0.271 0.050 10 4 60 includes 0 includes 0 cross-modal visual pre MM -0.383 -0.382 -0.603 -0.173 -0.591 -0.188 11 4 76 below below cross-modal visual pre ALM -0.207 -0.207 -0.494 0.050 -0.498 0.043 10 4 60 includes 0 includes 0 cross-modal visual pre AMM -0.221 -0.221 -0.356 -0.082 -0.346 -0.085 21 7 173 below below cross-modal visual pre PPC -0.173 -0.174 -0.333 -0.031 -0.309 -0.052 21 7 142 below below cross-modal visual post sham -0.081 -0.082 -0.237 0.066 -0.228 0.059 28 7 190 includes 0 includes 0 cross-modal visual post S1/S2 -0.197 -0.198 -0.460 0.045 -0.448 0.046 10 4 60 includes 0 includes 0 cross-modal visual post MM -0.673 -0.672 -0.839 -0.518 -0.825 -0.529 11 4 70 below below cross-modal visual post ALM -0.631 -0.631 -0.870 -0.385 -0.839 -0.391 10 4 91 below below cross-modal visual post AMM -0.542 -0.543 -0.725 -0.353 -0.717 -0.353 21 7 144 below below cross-modal visual post PPC -0.301 -0.302 -0.555 -0.070 -0.538 -0.074 21 7 165 below below detection sensitivity pre sham 0.066 0.066 -0.100 0.240 -0.071 0.211 9 5 207 includes 0 includes 0 detection sensitivity pre S1/S2 -0.038 -0.037 -0.373 0.243 -0.337 0.178 3 3 63 includes 0 includes 0 detection sensitivity pre MM -0.094 -0.094 -0.270 0.087 -0.230 0.055 5 3 130 includes 0 includes 0 detection sensitivity pre ALM -0.113 -0.114 -0.320 0.112 -0.281 0.067 4 3 83 includes 0 includes 0 detection sensitivity pre AMM -0.005 -0.007 -0.264 0.198 -0.239 0.158 4 3 94 includes 0 includes 0 detection sensitivity pre PPC -0.004 -0.003 -0.175 0.159 -0.119 0.103 4 3 90 includes 0 includes 0 detection sensitivity post sham 0.034 0.035 -0.152 0.225 -0.130 0.200 9 5 206 includes 0 includes 0 detection sensitivity post S1/S2 -0.448 -0.448 -0.643 -0.251 -0.593 -0.304 3 3 69 below below detection sensitivity post MM -0.434 -0.435 -0.628 -0.233 -0.593 -0.259 5 3 121 below below detection sensitivity post ALM -0.266 -0.265 -0.471 -0.065 -0.446 -0.098 4 3 80 below below detection sensitivity post AMM -0.106 -0.107 -0.331 0.127 -0.301 0.103 4 3 97 includes 0 includes 0 detection sensitivity post PPC -0.079 -0.079 -0.282 0.106 -0.238 0.059 4 3 98 includes 0 includes 0
In [32]:
# Monte Carlo spread of the CI bounds: the same hierarchical bootstraps (10000 resamples) with ten other seeds.
judged = [("cross-modal", "tactile", "pre", c) for c in ("S1/S2", "MM", "ALM", "sham", "AMM", "PPC")] + \
[("cross-modal", "tactile_hit", "pre", c) for c in ("S1/S2", "MM", "ALM")] + \
[("cross-modal", "tactile", "post", c) for c in ("S1/S2", "MM", "ALM", "PPC")] + \
[("cross-modal", "visual", w, c) for w in ("pre", "post") for c in ("S1/S2", "MM", "ALM", "AMM", "PPC")] + \
[("detection", "sensitivity", w, c) for w in ("pre", "post") for c in ("S1/S2", "MM", "ALM", "AMM", "PPC")]
main_rng = rng_opto
mc_rows = []
for seed in range(21, 31):
rng_opto = np.random.default_rng(seed)
for task, measure, window, cond in judged:
sub = optos[(optos.task == task) & (optos.condition == cond)]
r = change(sub, task, measure, 2 if window == "pre" else 3)
mc_rows.append({"task": task, "measure": measure, "window": window, "condition": cond, "seed": seed,
"lo": r["lo"], "hi": r["hi"], "code_lo": r["code_lo"], "code_hi": r["code_hi"]})
rng_opto = main_rng
mc = pd.DataFrame(mc_rows)
side = lambda lo, hi: np.where(hi < 0, "below", np.where(lo > 0, "above", "includes 0"))
mc["call"], mc["code_call"] = side(mc.lo, mc.hi), side(mc.code_lo, mc.code_hi)
g = mc.groupby(["task", "measure", "window", "condition"], sort=False)
mc_sum = pd.DataFrame({"lo_mean": g.lo.mean(), "lo_sd": g.lo.std(), "hi_mean": g.hi.mean(), "hi_sd": g.hi.std(),
"hi_min": g.hi.min(), "hi_max": g.hi.max(), "seeds_below_0": g.call.apply(lambda x: int((x == "below").sum())),
"code_hi_mean": g.code_hi.mean(), "code_hi_sd": g.code_hi.std(),
"code_seeds_below_0": g.code_call.apply(lambda x: int((x == "below").sum()))})
main_call = opto_res.set_index(["task", "measure", "window", "condition"])
mc_sum["main_seed_call"] = main_call.loc[mc_sum.index, "ci_excludes_0"].to_numpy()
mc_sum["calls_that_differ_from_main"] = g.call.apply(lambda x: int((x != main_call.loc[x.name, "ci_excludes_0"]).sum())).to_numpy()
print("CI bounds over ten other seeds (mean and SD of each bound; the SD is the Monte Carlo error of one run):")
print(mc_sum.round(3).to_string())
CI bounds over ten other seeds (mean and SD of each bound; the SD is the Monte Carlo error of one run):
lo_mean lo_sd hi_mean hi_sd hi_min hi_max seeds_below_0 code_hi_mean code_hi_sd code_seeds_below_0 main_seed_call calls_that_differ_from_main
task measure window condition
cross-modal tactile pre S1/S2 -0.473 0.005 0.099 0.004 0.094 0.105 0 0.075 0.002 0 includes 0 0
MM -0.661 0.002 -0.121 0.003 -0.128 -0.118 10 -0.133 0.003 10 below 0
ALM -0.710 0.004 -0.001 0.003 -0.006 0.001 7 -0.018 0.004 10 includes 0 7
sham -0.128 0.003 0.145 0.002 0.141 0.147 0 0.134 0.002 0 includes 0 0
AMM -0.212 0.002 0.076 0.001 0.073 0.078 0 0.062 0.001 0 includes 0 0
PPC -0.410 0.003 -0.051 0.002 -0.054 -0.048 10 -0.059 0.003 10 below 0
tactile_hit pre S1/S2 -0.304 0.002 0.013 0.002 0.010 0.016 0 0.011 0.002 0 includes 0 0
MM -0.496 0.003 -0.152 0.002 -0.155 -0.150 10 -0.152 0.003 10 below 0
ALM -0.465 0.002 -0.004 0.002 -0.007 -0.002 10 -0.003 0.002 9 below 0
tactile post S1/S2 -0.692 0.004 -0.075 0.004 -0.080 -0.067 10 -0.091 0.004 10 below 0
MM -0.611 0.001 -0.328 0.004 -0.336 -0.324 10 -0.341 0.003 10 below 0
ALM -0.568 0.002 -0.253 0.002 -0.256 -0.249 10 -0.272 0.002 10 below 0
PPC -0.467 0.004 -0.011 0.003 -0.017 -0.008 10 -0.032 0.001 10 below 0
visual pre S1/S2 -0.286 0.003 0.069 0.001 0.067 0.071 0 0.051 0.001 0 includes 0 0
MM -0.606 0.003 -0.176 0.003 -0.179 -0.172 10 -0.183 0.003 10 below 0
ALM -0.501 0.005 0.048 0.004 0.043 0.059 0 0.041 0.002 0 includes 0 0
AMM -0.356 0.001 -0.078 0.002 -0.082 -0.075 10 -0.086 0.002 10 below 0
PPC -0.332 0.002 -0.030 0.002 -0.034 -0.026 10 -0.053 0.001 10 below 0
post S1/S2 -0.461 0.003 0.049 0.003 0.046 0.053 0 0.043 0.002 0 includes 0 0
MM -0.840 0.002 -0.520 0.001 -0.523 -0.518 10 -0.526 0.003 10 below 0
ALM -0.873 0.004 -0.382 0.003 -0.386 -0.377 10 -0.394 0.003 10 below 0
AMM -0.725 0.002 -0.350 0.001 -0.352 -0.348 10 -0.358 0.003 10 below 0
PPC -0.553 0.003 -0.066 0.003 -0.072 -0.060 10 -0.071 0.004 10 below 0
detection sensitivity pre S1/S2 -0.372 0.003 0.246 0.002 0.242 0.250 0 0.179 0.001 0 includes 0 0
MM -0.272 0.002 0.087 0.002 0.085 0.091 0 0.056 0.002 0 includes 0 0
ALM -0.321 0.003 0.117 0.002 0.114 0.119 0 0.068 0.002 0 includes 0 0
AMM -0.261 0.005 0.197 0.002 0.192 0.199 0 0.158 0.001 0 includes 0 0
PPC -0.179 0.002 0.159 0.002 0.155 0.161 0 0.101 0.001 0 includes 0 0
post S1/S2 -0.650 0.003 -0.250 0.002 -0.254 -0.247 10 -0.306 0.002 10 below 0
MM -0.625 0.003 -0.231 0.003 -0.236 -0.226 10 -0.257 0.003 10 below 0
ALM -0.470 0.003 -0.067 0.002 -0.071 -0.063 10 -0.101 0.001 10 below 0
AMM -0.332 0.002 0.130 0.003 0.123 0.134 0 0.104 0.003 0 includes 0 0
PPC -0.278 0.002 0.106 0.003 0.100 0.110 0 0.060 0.001 0 includes 0 0
In [33]:
def show(task, measure, window, conds, paper):
d = opto_res[(opto_res.task == task) & (opto_res.measure == measure) & (opto_res.window == window)].set_index("condition").loc[conds]
d = d[["sessions", "mice", "point", "lo", "hi", "ci_excludes_0", "code_lo", "code_hi", "code_ci_excludes_0"]].copy()
d["paper_ci"] = paper
print(f"{task}, {measure}, laser {window}-stimulus")
print(d.round(3).to_string())
print()
show("cross-modal", "tactile", "pre", ["S1/S2", "MM", "ALM", "sham", "AMM"],
["[-0.47, 0.08]", "[-0.65, -0.12]", "[-0.68, -0.02]", "[-0.12, 0.13]", "[-0.19, 0.06]"])
cross-modal, tactile, laser pre-stimulus
sessions mice point lo hi ci_excludes_0 code_lo code_hi code_ci_excludes_0 paper_ci
condition
S1/S2 10 4 -0.178 -0.468 0.094 includes 0 -0.457 0.080 includes 0 [-0.47, 0.08]
MM 11 4 -0.390 -0.658 -0.123 below -0.642 -0.135 below [-0.65, -0.12]
ALM 10 4 -0.373 -0.707 0.005 includes 0 -0.682 -0.018 below [-0.68, -0.02]
sham 28 7 0.011 -0.126 0.144 includes 0 -0.118 0.134 includes 0 [-0.12, 0.13]
AMM 21 7 -0.062 -0.216 0.077 includes 0 -0.193 0.061 includes 0 [-0.19, 0.06]
Result: partially reproduced. Pre-stimulus inhibition of MM lowers tactile sensitivity (-0.390, CI -0.658 to -0.123), and S1/S2, sham and AMM CIs include 0; for ALM the change is -0.373 but its CI reaches 0.005 when false alarms are also resampled. The outcome is judged on the main seed. The ALM upper bound is at the edge: over ten other seeds it averages -0.001 (SD 0.003, -0.006 to 0.001) and is below 0 in 7 of 10, so the ALM interval touches 0 under the Methods' resampling. With the code's resampling the ALM CI is -0.682 to -0.018, as in the paper (-0.68 to -0.02), and below 0 in all ten other seeds.
In [34]:
show("cross-modal", "tactile_hit", "pre", ["S1/S2", "MM", "ALM"], ["[-0.31, 0.01]", "[-0.49, -0.15]", "[-0.46, -0.07]"])
cross-modal, tactile_hit, laser pre-stimulus
sessions mice point lo hi ci_excludes_0 code_lo code_hi code_ci_excludes_0 paper_ci
condition
S1/S2 10 4 -0.147 -0.303 0.010 includes 0 -0.306 0.010 includes 0 [-0.31, 0.01]
MM 11 4 -0.318 -0.492 -0.152 below -0.495 -0.151 below [-0.49, -0.15]
ALM 10 4 -0.227 -0.463 -0.007 below -0.466 0.002 includes 0 [-0.46, -0.07]
Result: reproduced. The tactile hit rate falls with MM (-0.318, CI -0.492 to -0.152) and ALM (-0.227, -0.463 to -0.007) inhibition, and the S1/S2 CI includes 0 (-0.303 to 0.010). The ALM upper bound is close to 0 but below it in all ten other seeds (mean -0.004, SD 0.002); the "code" column, the same procedure with another draw, gives 0.002 here and is below 0 in 9 of 10 other seeds.
In [35]:
show("cross-modal", "tactile", "post", ["S1/S2", "MM", "ALM"], ["[-0.67, -0.08]", "[-0.6, -0.34]", "[-0.55, -0.27]"])
cross-modal, tactile, laser post-stimulus
sessions mice point lo hi ci_excludes_0 code_lo code_hi code_ci_excludes_0 paper_ci
condition
S1/S2 10 4 -0.389 -0.699 -0.080 below -0.677 -0.092 below [-0.67, -0.08]
MM 11 4 -0.487 -0.612 -0.327 below -0.598 -0.334 below [-0.6, -0.34]
ALM 10 4 -0.402 -0.568 -0.253 below -0.545 -0.274 below [-0.55, -0.27]
Result: reproduced. Inhibition from the stimulus onset lowers tactile sensitivity for S1/S2 (-0.389, CI -0.699 to -0.080), MM (-0.487, -0.612 to -0.327) and ALM (-0.402, -0.568 to -0.253).
In [36]:
show("cross-modal", "visual", "pre", ["MM", "AMM", "ALM", "S1/S2"], ["[-0.6, -0.18]", "[-0.35, -0.08]", "[-0.5, 0.04]", "[-0.27, 0.05]"])
show("cross-modal", "visual", "post", ["MM", "ALM", "AMM", "S1/S2"], ["[-0.83, -0.52]", "[-0.84, -0.39]", "[-0.72, -0.36]", "[-0.44, 0.04]"])
cross-modal, visual, laser pre-stimulus
sessions mice point lo hi ci_excludes_0 code_lo code_hi code_ci_excludes_0 paper_ci
condition
MM 11 4 -0.383 -0.603 -0.173 below -0.591 -0.188 below [-0.6, -0.18]
AMM 21 7 -0.221 -0.356 -0.082 below -0.346 -0.085 below [-0.35, -0.08]
ALM 10 4 -0.207 -0.494 0.050 includes 0 -0.498 0.043 includes 0 [-0.5, 0.04]
S1/S2 10 4 -0.100 -0.282 0.070 includes 0 -0.271 0.050 includes 0 [-0.27, 0.05]
cross-modal, visual, laser post-stimulus
sessions mice point lo hi ci_excludes_0 code_lo code_hi code_ci_excludes_0 paper_ci
condition
MM 11 4 -0.673 -0.839 -0.518 below -0.825 -0.529 below [-0.83, -0.52]
ALM 10 4 -0.631 -0.870 -0.385 below -0.839 -0.391 below [-0.84, -0.39]
AMM 21 7 -0.542 -0.725 -0.353 below -0.717 -0.353 below [-0.72, -0.36]
S1/S2 10 4 -0.197 -0.460 0.045 includes 0 -0.448 0.046 includes 0 [-0.44, 0.04]
Result: reproduced. All eight judgements match: pre-stimulus MM (-0.603 to -0.173) and AMM (-0.356 to -0.082) CIs are below 0, ALM and S1/S2 include 0; post-stimulus MM, ALM and AMM are below 0 and S1/S2 includes 0 (-0.460 to 0.045).
In [37]:
for measure, window, paper in (("tactile", "pre", "[-0.4, -0.06]"), ("tactile", "post", "[-0.43, -0.04]"),
("visual", "pre", "[-0.31, -0.05]"), ("visual", "post", "[-0.54, -0.07]")):
show("cross-modal", measure, window, ["PPC"], [paper])
cross-modal, tactile, laser pre-stimulus
sessions mice point lo hi ci_excludes_0 code_lo code_hi code_ci_excludes_0 paper_ci
condition
PPC 21 7 -0.232 -0.408 -0.052 below -0.399 -0.06 below [-0.4, -0.06]
cross-modal, tactile, laser post-stimulus
sessions mice point lo hi ci_excludes_0 code_lo code_hi code_ci_excludes_0 paper_ci
condition
PPC 21 7 -0.218 -0.466 -0.009 below -0.435 -0.034 below [-0.43, -0.04]
cross-modal, visual, laser pre-stimulus
sessions mice point lo hi ci_excludes_0 code_lo code_hi code_ci_excludes_0 paper_ci
condition
PPC 21 7 -0.173 -0.333 -0.031 below -0.309 -0.052 below [-0.31, -0.05]
cross-modal, visual, laser post-stimulus
sessions mice point lo hi ci_excludes_0 code_lo code_hi code_ci_excludes_0 paper_ci
condition
PPC 21 7 -0.301 -0.555 -0.07 below -0.538 -0.074 below [-0.54, -0.07]
Result: reproduced. PPC inhibition lowers tactile sensitivity before (-0.408 to -0.052) and after (-0.466 to -0.009) the stimulus onset and visual sensitivity before (-0.333 to -0.031) and after (-0.555 to -0.070). The closest bound, tactile post-stimulus, is below 0 in all ten other seeds (mean -0.011, SD 0.003).
In [38]:
show("detection", "sensitivity", "pre", ["S1/S2", "MM", "ALM", "AMM", "PPC", "sham"],
["[-0.34, 0.18]", "[-0.23, 0.06]", "[-0.28, 0.07]", "[-0.25, 0.16]", "[-0.12, 0.1]", ""])
show("detection", "sensitivity", "post", ["S1/S2", "MM", "ALM", "AMM", "PPC", "sham"],
["[-0.59, -0.31]", "[-0.59, -0.25]", "[-0.45, -0.1]", "[-0.17, 0.18]", "[-0.14, 0.3]", ""])
detection, sensitivity, laser pre-stimulus
sessions mice point lo hi ci_excludes_0 code_lo code_hi code_ci_excludes_0 paper_ci
condition
S1/S2 3 3 -0.038 -0.373 0.243 includes 0 -0.337 0.178 includes 0 [-0.34, 0.18]
MM 5 3 -0.094 -0.270 0.087 includes 0 -0.230 0.055 includes 0 [-0.23, 0.06]
ALM 4 3 -0.113 -0.320 0.112 includes 0 -0.281 0.067 includes 0 [-0.28, 0.07]
AMM 4 3 -0.005 -0.264 0.198 includes 0 -0.239 0.158 includes 0 [-0.25, 0.16]
PPC 4 3 -0.004 -0.175 0.159 includes 0 -0.119 0.103 includes 0 [-0.12, 0.1]
sham 9 5 0.066 -0.100 0.240 includes 0 -0.071 0.211 includes 0
detection, sensitivity, laser post-stimulus
sessions mice point lo hi ci_excludes_0 code_lo code_hi code_ci_excludes_0 paper_ci
condition
S1/S2 3 3 -0.448 -0.643 -0.251 below -0.593 -0.304 below [-0.59, -0.31]
MM 5 3 -0.434 -0.628 -0.233 below -0.593 -0.259 below [-0.59, -0.25]
ALM 4 3 -0.266 -0.471 -0.065 below -0.446 -0.098 below [-0.45, -0.1]
AMM 4 3 -0.106 -0.331 0.127 includes 0 -0.301 0.103 includes 0 [-0.17, 0.18]
PPC 4 3 -0.079 -0.282 0.106 includes 0 -0.238 0.059 includes 0 [-0.14, 0.3]
sham 9 5 0.034 -0.152 0.225 includes 0 -0.130 0.200 includes 0
Result: reproduced. In the detection task, pre-stimulus inhibition CIs include 0 for S1/S2, MM and ALM, and post-stimulus CIs are below 0 (S1/S2 -0.643 to -0.251, MM -0.628 to -0.233, ALM -0.471 to -0.065); AMM and PPC include 0 in both windows.
In [39]:
rows = []
for task in ("cross-modal", "detection"):
for duration in ("short", "long"):
for cond in ("sham", "S1/S2", "MM", "ALM", "AMM", "PPC"):
sub = optos[(optos.task == task) & (optos.condition == cond)]
if task == "cross-modal":
iti, rw = ((1, 3, 1, 1), (1, 3, 1, 2)) if duration == "short" else ((1, 4, 1, 1), (1, 4, 1, 3))
else:
iti, rw = ((4, 1), (4, 3)) if duration == "short" else ((5, 1), (5, 3))
c = {"iti": counts_for(sub, iti), "rw": counts_for(sub, rw)}
f = lambda p: p["rw"] - p["iti"]
b = xm.hierarchical_bootstrap(sub.mouse.to_numpy(), c, f, 10000, rng_opto)
lo, hi = np.quantile(b, [0.025, 0.975])
rows.append({"task": task, "laser": duration, "condition": cond, "laser_only_trials": int(c["rw"][0].sum()),
"p_iti": c["iti"][1].sum() / max(c["iti"][0].sum(), 1), "p_response": c["rw"][1].sum() / max(c["rw"][0].sum(), 1),
"point": xm.point_estimate(sub.mouse.to_numpy(), c, f), "lo": lo, "hi": hi,
"includes_0": bool(lo <= 0 <= hi)})
for cond in ("sham", "S1/S2", "MM", "ALM", "AMM", "PPC"):
sub = optos[(optos.task == "detection") & (optos.condition == cond)]
c = {"iti": tuple(a + b for a, b in zip(counts_for(sub, (4, 1)), counts_for(sub, (5, 1)))),
"rw": tuple(a + b for a, b in zip(counts_for(sub, (4, 3)), counts_for(sub, (5, 3))))}
f = lambda p: p["rw"] - p["iti"]
b = xm.hierarchical_bootstrap(sub.mouse.to_numpy(), c, f, 10000, rng_opto)
lo, hi = np.quantile(b, [0.025, 0.975])
rows.append({"task": "detection", "laser": "pooled", "condition": cond, "laser_only_trials": int(c["rw"][0].sum()),
"p_iti": c["iti"][1].sum() / max(c["iti"][0].sum(), 1), "p_response": c["rw"][1].sum() / max(c["rw"][0].sum(), 1),
"point": xm.point_estimate(sub.mouse.to_numpy(), c, f), "lo": lo, "hi": hi, "includes_0": bool(lo <= 0 <= hi)})
laser_only = pd.DataFrame(rows)
paper_lo = {("cross-modal", "S1/S2"): "[-0.18, 0.24]", ("cross-modal", "MM"): "[-0.31, 0.03]", ("cross-modal", "ALM"): "[-0.24, 0.16]",
("detection", "S1/S2"): "[-0.13, 0.31]", ("detection", "MM"): "[-0.06, 0.45]", ("detection", "ALM"): "[-0.18, 0.34]"}
laser_only["paper_ci"] = [paper_lo.get((t, c), "") for t, c in zip(laser_only.task, laser_only.condition)]
print(laser_only.round(3).to_string(index=False))
key = laser_only[laser_only.condition.isin(["S1/S2", "MM", "ALM"]) & (laser_only.laser != "pooled")]
print("\nS1/S2, MM and ALM CIs including 0:", int(key.includes_0.sum()), "of", len(key))
task laser condition laser_only_trials p_iti p_response point lo hi includes_0 paper_ci cross-modal short sham 403 0.218 0.213 -0.000 -0.084 0.090 True cross-modal short S1/S2 163 0.196 0.202 0.005 -0.117 0.128 True [-0.18, 0.24] cross-modal short MM 124 0.242 0.226 0.005 -0.144 0.156 True [-0.31, 0.03] cross-modal short ALM 131 0.168 0.305 0.114 -0.058 0.284 True [-0.24, 0.16] cross-modal short AMM 308 0.198 0.214 -0.021 -0.192 0.108 True cross-modal short PPC 298 0.141 0.168 -0.001 -0.113 0.101 True cross-modal long sham 429 0.184 0.200 0.014 -0.093 0.114 True cross-modal long S1/S2 173 0.179 0.225 0.033 -0.180 0.245 True [-0.18, 0.24] cross-modal long MM 138 0.210 0.109 -0.143 -0.312 0.030 True [-0.31, 0.03] cross-modal long ALM 149 0.221 0.242 -0.040 -0.243 0.152 True [-0.24, 0.16] cross-modal long AMM 332 0.154 0.145 0.012 -0.095 0.146 True cross-modal long PPC 331 0.160 0.411 0.096 -0.102 0.290 True detection short sham 236 0.347 0.394 0.053 -0.092 0.200 True detection short S1/S2 62 0.290 0.435 0.140 -0.042 0.309 True [-0.13, 0.31] detection short MM 141 0.333 0.553 0.258 0.045 0.483 False [-0.06, 0.45] detection short ALM 99 0.384 0.586 0.184 -0.037 0.406 True [-0.18, 0.34] detection short AMM 120 0.325 0.467 0.135 -0.090 0.356 True detection short PPC 97 0.320 0.505 0.200 0.026 0.383 False detection long sham 220 0.386 0.423 0.062 -0.131 0.277 True detection long S1/S2 61 0.377 0.361 -0.014 -0.197 0.167 True [-0.13, 0.31] detection long MM 146 0.377 0.411 0.024 -0.240 0.257 True [-0.06, 0.45] detection long ALM 111 0.360 0.324 -0.042 -0.309 0.231 True [-0.18, 0.34] detection long AMM 123 0.374 0.504 0.108 -0.077 0.282 True detection long PPC 99 0.394 0.596 0.200 0.017 0.383 False detection pooled sham 456 0.366 0.408 0.051 -0.057 0.164 True detection pooled S1/S2 123 0.333 0.398 0.060 -0.074 0.194 True [-0.13, 0.31] detection pooled MM 287 0.355 0.481 0.149 0.020 0.274 False [-0.06, 0.45] detection pooled ALM 210 0.371 0.448 0.068 -0.114 0.291 True [-0.18, 0.34] detection pooled AMM 243 0.350 0.486 0.125 -0.005 0.257 True detection pooled PPC 196 0.357 0.551 0.201 0.061 0.353 False S1/S2, MM and ALM CIs including 0: 11 of 12
In [40]:
# Monte Carlo spread of the laser-only CIs: the same hierarchical bootstraps (10000 resamples) with ten other seeds.
lo_rows = []
for seed in range(21, 31):
rng_lo = np.random.default_rng(seed)
for task in ("cross-modal", "detection"):
for duration in ("short", "long"):
for cond in ("S1/S2", "MM", "ALM"):
sub = optos[(optos.task == task) & (optos.condition == cond)]
if task == "cross-modal":
iti, rw = ((1, 3, 1, 1), (1, 3, 1, 2)) if duration == "short" else ((1, 4, 1, 1), (1, 4, 1, 3))
else:
iti, rw = ((4, 1), (4, 3)) if duration == "short" else ((5, 1), (5, 3))
c = {"iti": counts_for(sub, iti), "rw": counts_for(sub, rw)}
b = xm.hierarchical_bootstrap(sub.mouse.to_numpy(), c, lambda p: p["rw"] - p["iti"], 10000, rng_lo)
lo, hi = np.quantile(b, [0.025, 0.975])
lo_rows.append({"task": task, "laser": duration, "condition": cond, "seed": seed, "lo": lo, "hi": hi,
"includes_0": bool(lo <= 0 <= hi)})
lo_mc = pd.DataFrame(lo_rows)
g = lo_mc.groupby(["task", "laser", "condition"], sort=False)
lo_sum = pd.DataFrame({"lo_mean": g.lo.mean(), "lo_sd": g.lo.std(), "lo_min": g.lo.min(), "lo_max": g.lo.max(),
"hi_mean": g.hi.mean(), "hi_sd": g.hi.std(), "seeds_including_0": g.includes_0.sum()})
main_lo = key.set_index(["task", "laser", "condition"]).includes_0
lo_sum["main_seed_includes_0"] = main_lo.loc[lo_sum.index].to_numpy()
print("Laser-only CI bounds over ten other seeds:")
print(lo_sum.round(3).to_string())
Laser-only CI bounds over ten other seeds:
lo_mean lo_sd lo_min lo_max hi_mean hi_sd seeds_including_0 main_seed_includes_0
task laser condition
cross-modal short S1/S2 -0.119 0.002 -0.122 -0.117 0.130 0.002 10 True
MM -0.140 0.002 -0.142 -0.136 0.153 0.001 10 True
ALM -0.060 0.002 -0.063 -0.057 0.289 0.002 10 True
long S1/S2 -0.180 0.003 -0.187 -0.175 0.246 0.002 10 True
MM -0.312 0.003 -0.317 -0.308 0.033 0.002 10 True
ALM -0.244 0.003 -0.248 -0.240 0.156 0.001 10 True
detection short S1/S2 -0.038 0.004 -0.044 -0.034 0.311 0.002 10 True
MM 0.045 0.004 0.040 0.053 0.483 0.004 0 False
ALM -0.036 0.002 -0.037 -0.031 0.404 0.003 10 True
long S1/S2 -0.200 0.004 -0.206 -0.192 0.171 0.003 10 True
MM -0.240 0.002 -0.243 -0.236 0.252 0.002 10 True
ALM -0.314 0.001 -0.314 -0.312 0.235 0.004 10 True
Result: not reproduced. 11 of the 12 CIs include 0, but on short laser-only trials of the detection task MM inhibition raises the right-lick probability (0.258, CI 0.045 to 0.483); over ten other seeds its lower bound is 0.040 to 0.053, and the other eleven CIs include 0 in every seed. All six cross-modal CIs and all long-laser CIs include 0. The windows are the code's (GetBehavInfo and GetBehavInfo_tDetection): in the cross-modal task licks on short laser-only trials are scored 0 to 2 s after the laser offset, in the detection task -0.8 to 1.2 s from it, which covers the laser. The paper's cross-modal CIs match the long-laser rows (MM -0.312 to 0.030 against -0.31 to 0.03); its detection CIs match neither the short, the long nor the pooled rows (the code's "all laser trials"; pooled MM 0.020 to 0.274 also excludes 0). A rise in licking does not bear on the paper's point that inhibition does not prevent licking.
Summary¶
| Claim | Outcome |
|---|---|
| fig1e-recording-counts | reproduced |
| fig6-session-counts | reproduced |
| fig1d-performance | reproduced |
| fig1d-outcome-fractions | reproduced |
| fig1h-thit-tcr-selective | partially reproduced |
| fig2c-rule-prestim-units | partially reproduced |
| fig2d-stimulus-prestim-units | reproduced |
| fig2e-stimulus-poststim-units | reproduced |
| fig2f-prestim-predicts-touch-response | reproduced |
| fig3b-rule-decoding | reproduced |
| fig3c-rule-decoding-motor-stronger | reproduced |
| fig3de-stimulus-decoding-controls | reproduced |
| figs3-svm-rf-decoding | reproduced |
| fig4c-prestim-distance-predicts-divergence | partially reproduced |
| fig4e-subspace-overlap | partially reproduced |
| fig4f-prestim-overlap-predicts-poststim | reproduced |
| figs4ef-coding-dimensions | partially reproduced |
| figs4h-choice-right-left | partially reproduced |
| fig5b-first-hit-switches | not reproduced |
| fig5de-motor-states-track-switch | reproduced |
| figs5-sensory-states-no-tracking | reproduced |
| fig6c-pre-inhibition-tactile | partially reproduced |
| figs6b-pre-inhibition-hit-rate | reproduced |
| fig6d-post-inhibition-tactile | reproduced |
| figs7gh-visual-sensitivity | reproduced |
| fig6-ppc-inhibition | reproduced |
| fig6fg-detection-task | reproduced |
| fig6ij-laser-only-licks | not reproduced |
| fig1e-probe-locations | 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 analysis is in Python, not MATLAB. Bootstrap resamples for AUC are shared by all units and bins of a session. The decoding CIs use 200 resamples per area, the transition CIs 500 and the coding-dimension CIs 1000. The transition classifier is LDA with Ledoit-Wolf shrinkage in place of the code's optimized LDA. Subspace overlaps average 100 equal splits of the tactile hits; the code's single split per neuron is a variant. Where the Methods and code disagree (normalization offsets, the session rule of the subspace analysis, resampling of false alarms in the optogenetic CIs), the Methods' choice is the main result and the code's is shown. The first-hit trial numbers are counted in rows of the files, which lack the removed trials, and in trial numbers estimated from the gaps between rows. Outcomes that rest on a bootstrap are judged on the main seed. Other seeds are shown for the single-neuron percentages and the figs4h S2 mean (five seeds each) and the optogenetic CIs (ten seeds). Over the main seed and four others, fig1h MM is within tolerance in 3 of 5, fig2c S1 in 0 of 5 and ALM in 2 of 5, and figs4h S2 in 1 of 5. On the five-seed means all four areas of fig1h and fig2c are within tolerance. Over ten other seeds the fig6c ALM interval under the Methods' resampling is below 0 in 7 of 10, and no other Methods-resampling call of the optogenetic CIs checked over seeds changes side; the figs6b ALM interval with the code's resampling includes 0 on the main draw (upper bound 0.002) and is below 0 in 9 of 10 other seeds. The fig6ij laser-only calls keep their side in all ten seeds.
What was found about the dandiset. Its counts of sessions, mice, units and optogenetic sessions
match the paper exactly, and most results reproduce closely, so the dandiset holds the data the paper
analysed. The ephys files hold only the trials and units the paper kept: trials with licks between -1 and
0.1 s of the stimulus and units failing the quality filter are absent, which the dandiset does not say.
Units carry no quality metrics, channel, depth or original cluster id; electrode coordinates are NaN and
locations "unknown", and no histology is included. In the behaviour files start_time equals the stimulus
onset and trials overlap; in the detection task, trials without a whisker stimulus are labelled
"Stim_Vis". Trial table columns have no descriptions, and codes such as response (1 right, 2 left,
0 none) and blockType ("Change" for the cue trial) are explained only by the paper's code.
In [41]:
print(f"Notebook ran in {(time.time() - started) / 60:.1f} minutes")
Notebook ran in 12.8 minutes