Reproducing Tseng et al. (2022) from dandiset 000579¶
Paper. Tseng, Chettih, Arlt, Barroso-Luque & Harvey (2022), "Shared and specialized coding across posterior cortical areas for dynamic navigation decisions", Neuron 110, 2484-2502, doi:10.1016/j.neuron.2022.05.012.
Dataset. DANDI:000579, version 0.230728.1727 (DOI 10.48324/dandi.000579/0.230728.1727). This is the only published version; the paper cites no version (it predates the dandiset), and the draft has the same 308 assets and size. The files are streamed; nothing is downloaded.
Claims. claims/000579/10.1016-j.neuron.2022.05.012.yaml, written from the accepted author manuscript in PubMed Central (PMC9357051) with its supplement (Figures S1 to S8).
What the dataset holds. 300 behavior + two-photon sessions of 8 mice (184 in layer 2/3, 116 in layer 5) and one widefield retinotopy file per mouse. Each session file holds the trial table with the paper's behavioral LSTM outputs (association matrix, rule belief, choice bias, rule- and bias-following), the frame-aligned position and treadmill velocities, the velocity LSTM outputs (dynamic choice and cue-biased running), the area label and CCF coordinates of each neuron, and deconvolved activity and dF/F. It holds no photoinhibition sessions and no GLM fits, encoding strengths or population decoders; the decoders are fit here.
| Claim | Tested here |
|---|---|
| neuron-counts | tested |
| fig1d-switch-performance | tested |
| lstm-choice-accuracy | tested, on the LSTM predictions stored in the files |
| figs1h-lstm-beats-rl | tested, RL model refit; LSTM from the stored predictions |
| fig1i-strategy-after-switch | tested, on the stored strategy variables |
| visually-guided-rule-update | tested, on the stored association matrices |
| fig1l, fig1m, figs1m (photoinhibition) | blocked, data missing |
| fig2c-running-variability-after-switch | tested |
| fig2f-latency-strategy | tested, on the stored dynamic choice |
| fig2g-choice-commitment | tested, on the stored dynamic choice |
| cue-biased-running-uncorrelated | tested, on the stored cue-biased running |
| fig4h-activity-by-epoch | tested, all-neuron reading only |
| GLM claims (fig3i-glm-fit to fig6h-embedding-overlap, 23 claims) | blocked, too large (the timing is below) |
| fig3q-cue-decoding, figs5f-dynamic-choice-decoding, fig4d-partial-correlation, fig8cd-position-specific-decoders | tested, population decoders on a random subsample of 18 layer 2/3 sessions |
| fig8g-shattering-dimensionality | blocked, too large (the timing is below) |
| figs5l-dynamic-choice-over-cue-biased-running | blocked: it depends on well-fit-neurons, which needs the GLM |
import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
import time
import warnings
warnings.filterwarnings("ignore") # empty-slice and deprecation warnings; dropped data are counted in the cells
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import tseng
from tseng import AREAS, hboot, hboot_multi, p_two, p_unpaired, bh
T0 = time.time()
print("dandiset", tseng.DANDISET, "version", tseng.VERSION)
def pb(boot):
"""Two-tailed bootstrap p of a difference against zero, printed as a bound when no sample crosses zero."""
p = p_two(boot)
return f"< {1 / len(boot):g}" if p == 0 else f"= {p:.3g}"
def pu(b1, b2):
"""Two-tailed p of two unpaired bootstrap samples, printed as a bound when they do not overlap."""
p = p_unpaired(b1, b2)
return f"< {1 / len(b1):g}" if p == 0 else f"= {p:.3g}"
dandiset 000579 version 0.230728.1727
assets = tseng.list_assets()
sessions = assets[assets.kind == "session"].reset_index(drop=True)
# Mean epoch activity is read for the layer 2/3 sessions only (the paper analyses no layer 5 data).
res = tseng.load_many(list(sessions.asset_id),
activity_ids=list(sessions.asset_id[sessions.layer == "L23"]), workers=12)
print(f"loaded {len(res)} sessions in {time.time() - T0:.0f} s")
loaded 300 sessions in 179 s
failed = {a: r for a, r in res.items() if isinstance(r, Exception)}
print("sessions that failed to load:", len(failed), [repr(e)[:200] for e in failed.values()])
rows = []
for _, s in sessions.iterrows():
r = res[s.asset_id]
if isinstance(r, Exception):
continue
T = r["trials"]
nv = T[~T.is_vis]
rows.append(dict(asset_id=s.asset_id, mouse=s.subject, date=s.date, fov=s.fov_area, layer=s.layer,
imaging=s.imaging, n_trials=len(T), n_vis=int(T.is_vis.sum()),
fc=nv.is_correct.mean(), fc_all=T.is_correct.mean(),
n_switch=int(T.is_switch.sum()), n_neurons=len(r["neurons"]),
meta_layer=r["imaging"]["cortical_layer"], n_planes=r["imaging"]["num_of_imaging_plane"],
n_frames=r["n_frames"], stem=r["stem"], maze_len=r["maze_len"],
cue_off=1 - r["task"]["cue_delay_length"] / (r["task"]["maze_stem_length"] + r["task"]["maze_arm_length"]),
**{a: int((r["neurons"].area == a).sum()) for a in AREAS}))
D = pd.DataFrame(rows)
D["vis_frac"] = D.n_vis / D.n_trials
sessions that failed to load: 0 []
Data check¶
print("assets:", len(assets), "| session files:", len(sessions), "| widefield files:", (assets.kind == "widefield").sum())
print("mice:", sorted(D.mouse.unique(), key=int))
print("\nsessions per mouse and layer:")
print(pd.crosstab(D.mouse, D.layer, margins=True))
print("\nlayer in the file name against the file's Imaging metadata:")
print(pd.crosstab(D.layer, D.meta_layer))
print("\nimaging type:", D.imaging.value_counts().to_dict(), "| planes:", D.n_planes.value_counts().to_dict())
print("\nneurons: total", D.n_neurons.sum(), "| layer 2/3", D[D.layer == "L23"].n_neurons.sum(),
"| layer 5", D[D.layer == "L5"].n_neurons.sum())
lab = set(np.concatenate([res[a]["neurons"].area.unique() for a in D.asset_id]))
print("area labels used:", sorted(lab))
print("\ntrials per session: min", D.n_trials.min(), "median", D.n_trials.median(), "max", D.n_trials.max())
print("rule switches per session:", D.n_switch.value_counts().to_dict())
print("fraction of visually guided trials per session: median %.3f, range %.3f to %.3f, sessions above 0.2: %d"
% (D.vis_frac.median(), D.vis_frac.min(), D.vis_frac.max(), (D.vis_frac > 0.2).sum()))
print("session fraction correct, visually guided trials excluded: %.3f to %.3f; all trials: %.3f to %.3f"
% (D.fc.min(), D.fc.max(), D.fc_all.min(), D.fc_all.max()))
tp = pd.DataFrame([res[a]["task"] for a in D.asset_id])
print("\nmaze stem (virmen units), sessions:", tp.maze_stem_length.value_counts().to_dict(), "| arm:", tp.maze_arm_length.unique(),
"| m per unit:", tp.meter_per_virmen_unit.unique())
print("feedback delay s:", tp.feedback_delay_sec.unique(), "| reward delay s:", tp.reward_delay_sec.unique(),
"| ITI correct/incorrect s:", tp.iti_correct_sec.unique(), tp.iti_incorrect_sec.unique())
# Ranges of the behavioural variables used
r0 = [res[a] for a in D.asset_id]
pos = np.concatenate([r["rnn"].posF.dropna().values for r in r0])
dc = np.concatenate([r["rnn"].dyn_choice.dropna().values for r in r0])
cr = np.concatenate([r["rnn"].cue_run.dropna().values for r in r0])
print("\nforward position at the LSTM samples (m): %.3f to %.3f" % (pos.min(), pos.max()))
print("dynamic choice: %.4f to %.4f; cue-biased running: %.4f to %.4f" % (dc.min(), dc.max(), cr.min(), cr.max()))
print("velocity LSTM samples on an imaging frame timestamp: %.4f" % np.mean(np.concatenate([r["rnn"].ts_match.values for r in r0])))
vb = np.concatenate([r["velbin"][20].reshape(-1, 3, 20).transpose(0, 2, 1).reshape(-1, 3) for r in r0])
print("treadmill velocity (a.u., position-binned) 1st/50th/99th percentiles:\n", np.nanpercentile(vb, [1, 50, 99], axis=0).round(1))
allT = pd.concat([res[a]["trials"].assign(asset_id=a) for a in D.asset_id])
for c in ["rule_belief", "signed_bias", "rule_following", "bias_following", "prob_actual_cho"]:
print(f"{c}: {allT[c].min():.3f} to {allT[c].max():.3f}, NaN {allT[c].isna().sum()}")
print("trial duration (start to stop) s: median %.2f; maze epoch s: median %.2f" % (allT.trial_dur.median(), allT.maze_dur.median()))
assets: 308 | session files: 300 | widefield files: 8
mice: ['3', '4', '5', '6', '7', '8', '9', '10']
sessions per mouse and layer:
layer L23 L5 All
mouse
10 26 25 51
3 20 3 23
4 21 9 30
5 18 14 32
6 27 19 46
7 24 11 35
8 28 14 42
9 20 21 41
All 184 116 300
layer in the file name against the file's Imaging metadata:
meta_layer 2/3 5
layer
L23 184 0
L5 0 116
imaging type: {'multi': 287, 'single': 13} | planes: {4: 287, 1: 13}
neurons: total 278155 | layer 2/3 120424 | layer 5 157731
area labels used: ['A', 'AM', 'MM', 'PM', 'RSC', 'V1']
trials per session: min 320 median 408.5 max 646
rule switches per session: {2: 293, 3: 7}
fraction of visually guided trials per session: median 0.233, range 0.000 to 0.782, sessions above 0.2: 222
session fraction correct, visually guided trials excluded: 0.314 to 0.882; all trials: 0.346 to 0.894
maze stem (virmen units), sessions: {200.0: 293, 150.0: 7} | arm: [50.] | m per unit: [0.01]
feedback delay s: [1.] | reward delay s: [3.] | ITI correct/incorrect s: [2.] [4. 2.]
forward position at the LSTM samples (m): 0.050 to 2.566 dynamic choice: 0.0003 to 0.9997; cue-biased running: 0.0002 to 0.9997 velocity LSTM samples on an imaging frame timestamp: 1.0000
treadmill velocity (a.u., position-binned) 1st/50th/99th percentiles: [[ 61.5 -1261.7 -2177.6] [ 729.8 25.1 -128.8] [ 1242.1 1348.2 2099.4]] rule_belief: -0.916 to 0.935, NaN 0 signed_bias: -0.958 to 0.977, NaN 0 rule_following: -0.450 to 0.467, NaN 0 bias_following: -0.479 to 0.489, NaN 0 prob_actual_cho: 0.012 to 0.996, NaN 0 trial duration (start to stop) s: median 14.58; maze epoch s: median 8.91
Checks of the metadata against the data¶
The strategy variables are checked against the association matrix stored in the same table, with the Methods' formulas (Decision-making strategy variables): rule belief = P(rule B) - P(rule A), signed bias = P(L) - P(R), bias-following = P(actual choice) - 0.5, rule-following = P(actual choice | actual cue) - P(actual choice). The rule column is_ruleA is checked against the outcomes (rule A: black-left and white-right are correct).
chk = []
for a in D.asset_id:
T = res[a]["trials"]
nv = T[~T.is_vis]
ruleA_ok = (nv.is_cueB & nv.is_choL) | (~nv.is_cueB & ~nv.is_choL)
PL = (T.P_L_B + T.P_L_W) / 2
PR = (T.P_R_B + T.P_R_W) / 2
rbB = (T.P_L_W + T.P_R_B) / 2 - (T.P_L_B + T.P_R_W) / 2
pch = np.where(T.is_choL, PL, PR)
chk.append(dict(
outcome_matches_rule=((ruleA_ok == nv.is_correct) == True)[nv.is_ruleA].mean() * 0.5
+ ((~ruleA_ok == nv.is_correct) == True)[~nv.is_ruleA].mean() * 0.5,
rb_eq_methods=np.nanmax(np.abs(T.rule_belief - rbB)),
rb_eq_minus_methods=np.nanmax(np.abs(T.rule_belief + rbB)),
signed_bias=np.nanmax(np.abs(T.signed_bias - (PL - PR))),
bias_following=np.nanmax(np.abs(T.bias_following - (pch - 0.5))),
rule_following=np.nanmax(np.abs(T.rule_following - (T.prob_actual_cho - pch))),
r_rulebelief_ruleB=np.corrcoef(nv.rule_belief, np.where(nv.is_ruleA, -1, 1))[0, 1]))
chk = pd.DataFrame(chk)
print(chk.describe().loc[["min", "50%", "max"]].round(4).T)
min 50% max outcome_matches_rule 1.0000 1.0000 1.0000 rb_eq_methods 0.7971 1.5777 1.8692 rb_eq_minus_methods 0.0000 0.0000 0.0000 signed_bias 0.0000 0.0000 0.0000 bias_following 0.0000 0.0000 0.0000 rule_following 0.0000 0.0000 0.0000 r_rulebelief_ruleB -0.8970 -0.7252 -0.1322
The stored rule_belief equals P(rule A) - P(rule B) to the last digit in every session, and correlates negatively with the current rule being B. The Methods and the dataset's tutorial both say it is positive for rule B. The other strategy variables match the Methods exactly. Below, rule belief is recomputed from the association matrix with the Methods' sign.
Clock check. The velocity LSTM outputs are on the imaging frame clock (every sample is on a frame timestamp, printed above). That the frame-aligned behaviour and the trial table share a clock is shown by a contrast: the position-binned roll velocity separates left from right choices in the arm and not at the start of the stem, and stored dynamic choice late in the maze correlates with the trial's own reported choice far more than with the previous trial's.
rows = []
for a in D.asset_id:
r = res[a]
T = r["trials"]
v = r["velbin"][20][:, 1, :] # roll, 20 bins
L, R = T.is_choL.values, ~T.is_choL.values
late = np.nanmean(v[:, 16:19], 1)
early = np.nanmean(v[:, 0:2], 1)
q = r["rnn"]
q = q[(q.posF > 0.9 * r["maze_len"]) & q.trial.notna() & q.dyn_choice.notna()]
m = q.groupby("trial").dyn_choice.mean()
cho = T.is_choL.reindex(m.index.astype(int)).values.astype(float)
cho_shift = T.is_choL.shift(1).reindex(m.index.astype(int)).values.astype(float)
ok = ~np.isnan(cho_shift)
rows.append(dict(mouse=D.set_index("asset_id").loc[a, "mouse"],
roll_arm_LminusR=np.nanmean(late[L]) - np.nanmean(late[R]),
roll_start_LminusR=np.nanmean(early[L]) - np.nanmean(early[R]),
r_dc=stats.pearsonr(m.values, cho)[0],
r_dc_prev=stats.pearsonr(m.values[ok], cho_shift[ok])[0]))
clk = pd.DataFrame(rows)
print(clk.drop(columns="mouse").describe().loc[["mean", "min", "max"]].round(3))
print("sessions where late-maze dynamic choice correlates with the trial's own choice above 0.5:",
(clk.r_dc > 0.5).sum(), "of", len(clk), "| with the previous trial's choice:", (clk.r_dc_prev > 0.5).sum())
roll_arm_LminusR roll_start_LminusR r_dc r_dc_prev mean 929.471 156.413 0.938 0.186 min 243.004 -46.535 0.699 -0.079 max 2053.178 663.920 0.994 0.504 sessions where late-maze dynamic choice correlates with the trial's own choice above 0.5: 300 of 300 | with the previous trial's choice: 1
Further checks of what the analyses assume: that is_switch marks the first trial of the new rule (the rule column changes there and nowhere else); that a trial's start_time is maze onset (the first frame with a maze position, used for the latencies of fig2f); and how the frames of the 1 s feedback delay and the end of the maze are assigned (the epochs of fig4h).
sw_ok, sw_n, lag = 0, 0, []
for a in D.asset_id:
T = res[a]["trials"]
flip = np.r_[False, T.is_ruleA.values[1:] != T.is_ruleA.values[:-1]]
sw_n += int(T.is_switch.sum())
sw_ok += bool((flip == T.is_switch.values).all())
lag.append(T.onset_lag.values)
print(f"sessions where is_ruleA changes exactly at the is_switch trials and nowhere else: {sw_ok} of {len(D)} ({sw_n} switches)")
lag = np.concatenate(lag)
print(f"first frame with a maze position minus the trial's start_time (s): 1st/50th/99th percentiles "
f"{np.round(np.nanpercentile(lag, [1, 50, 99]), 3).tolist()}; trials more than 0.1 s apart: {(np.abs(lag) > 0.1).sum()} of {np.isfinite(lag).sum()}")
dg = pd.DataFrame([res[a]["diag"] for a in D.asset_id])
print(f"frames in the 1 s feedback delay: {dg.fd_frames.sum():,}, of which with a maze position: {dg.fd_with_pos.sum()}")
print(f"in-trial frames before the choice point: {dg.pre_frames.sum():,}; without a position: {dg.pre_no_pos.sum()}; "
f"beyond the maze end of 2.5 m: {dg.pre_beyond_maze.sum():,} (largest position {dg.max_pos.max():.3f} m)")
sessions where is_ruleA changes exactly at the is_switch trials and nowhere else: 300 of 300 (607 switches) first frame with a maze position minus the trial's start_time (s): 1st/50th/99th percentiles [-0.022, -0.004, 0.014]; trials more than 0.1 s apart: 183 of 122763 frames in the 1 s feedback delay: 3,673,732, of which with a maze position: 0 in-trial frames before the choice point: 35,156,721; without a position: 0; beyond the maze end of 2.5 m: 913,421 (largest position 2.566 m)
Result of the checks. is_ruleA changes exactly at the is_switch trials and nowhere else in all 300 sessions (607 switches), so is_switch marks the first trial of the new rule. A trial's start_time is maze onset: the first frame with a maze position falls -0.022 to 0.014 s from it (1st to 99th percentile; 183 of 122,763 trials more than 0.1 s apart). No frame of the 1 s feedback delay has a maze position, so those frames belong to the feedback/ITI epoch; 913,421 frames before the choice point lie beyond the 2.5 m maze end (up to 2.566 m) and are counted in the arm.
Session selection¶
The Methods give two thresholds on a session's overall fraction correct (Task performance analysis): above 0.5 for Figure 1D and S1B to S1F, above 0.65 for "subsequent behavioral and neural analyses", and 141 layer 2/3 sessions for the neural analyses (Classification of the sources). They also say visually guided trials are excluded "from all analyses unless mentioned otherwise" (Behavioral training), which allows two readings of fraction correct: (a) visually guided trials excluded, (b) all trials. Neither gives the paper's session counts exactly (below). Reading b is adopted, as an assumption: its layer 2/3 neuron, session and area counts are within 10% of the paper's (neuron-counts), while reading a keeps only 100 layer 2/3 sessions against 141, and the threshold that keeps 141 layer 2/3 sessions on the all-trials fraction is close to 0.65. Every analysis that depends on the selection is also run under reading a.
for thr, paper in [(0.5, "253 sessions (Figure S1B)"), (0.65, "216 sessions (Figure S1G, S1H); 141 layer 2/3 sessions")]:
for col, name in [("fc", "a, visually guided excluded"), ("fc_all", "b, all trials")]:
sel = D[D[col] > thr]
print(f"> {thr} ({name}): {len(sel)} sessions, {sel[sel.layer == 'L23'].shape[0]} layer 2/3, "
f"{sel.mouse.nunique()} mice [paper: {paper}]")
for n, col in [(253, "fc"), (216, "fc"), (253, "fc_all"), (216, "fc_all")]:
print(f"threshold that keeps {n} sessions on {col}: {np.sort(D[col].values)[::-1][n - 1]:.3f}")
l23 = D[D.layer == "L23"]
for col in ["fc", "fc_all"]:
print(f"threshold that keeps 141 layer 2/3 sessions on {col}: {np.sort(l23[col].values)[::-1][140]:.3f}")
SEL = {"a": D[D.fc > 0.65], "b": D[D.fc_all > 0.65]}
SEL50 = {"a": D[D.fc > 0.5], "b": D[D.fc_all > 0.5]}
> 0.5 (a, visually guided excluded): 291 sessions, 180 layer 2/3, 8 mice [paper: 253 sessions (Figure S1B)] > 0.5 (b, all trials): 297 sessions, 183 layer 2/3, 8 mice [paper: 253 sessions (Figure S1B)] > 0.65 (a, visually guided excluded): 154 sessions, 100 layer 2/3, 8 mice [paper: 216 sessions (Figure S1G, S1H); 141 layer 2/3 sessions] > 0.65 (b, all trials): 238 sessions, 150 layer 2/3, 8 mice [paper: 216 sessions (Figure S1G, S1H); 141 layer 2/3 sessions] threshold that keeps 253 sessions on fc: 0.576 threshold that keeps 216 sessions on fc: 0.610 threshold that keeps 253 sessions on fc_all: 0.621 threshold that keeps 216 sessions on fc_all: 0.664 threshold that keeps 141 layer 2/3 sessions on fc: 0.604 threshold that keeps 141 layer 2/3 sessions on fc_all: 0.659
neuron-counts¶
The neural analyses use 93,881 layer 2/3 neurons from 141 sessions of 8 mice, parcellated into V1 (14,373), PM (9,564), AM (14,974), MM (9,885), RSC (23,036) and A (22,049).
Counted over the neurons of each layer 2/3 session above 65% correct, with the area label stored for each neuron. The files hold one set of sources per session, which the dandiset calls neurons; their total over the 300 sessions is printed in the data check. The paper's cell-body classifier is not in the files, so whether the stored sources are the cell bodies alone is checked through that total.
exp = dict(neurons=93881, sessions=141, mice=8, V1=14373, PM=9564, AM=14974, MM=9885, RSC=23036, A=22049)
print("all 300 sessions: neurons", D.n_neurons.sum(), "(paper: 278,155 neurons in 300 sessions from 8 mice)")
out = []
for k, sel in SEL.items():
s = sel[sel.layer == "L23"]
got = dict(neurons=s.n_neurons.sum(), sessions=len(s), mice=s.mouse.nunique(), **{a: s[a].sum() for a in AREAS})
for q, v in got.items():
tol = 0 if q == "mice" else round(0.1 * exp[q])
out.append(dict(reading=k, quantity=q, measured=v, expected=exp[q], tolerance=tol,
within=abs(v - exp[q]) <= tol))
out = pd.DataFrame(out)
print(out.to_string(index=False))
s = SEL["b"][SEL["b"].layer == "L23"]
print("\nreading b, area shares:", (s[AREAS].sum() / s[AREAS].sum().sum()).round(3).to_dict())
print("paper, area shares:", {a: round(exp[a] / 93881, 3) for a in AREAS})
print("\nby mouse (reading b):")
print(s.groupby("mouse")[["n_neurons"] + AREAS].sum())
all 300 sessions: neurons 278155 (paper: 278,155 neurons in 300 sessions from 8 mice)
reading quantity measured expected tolerance within
a neurons 69707 93881 9388 False
a sessions 100 141 14 False
a mice 8 8 0 True
a V1 9039 14373 1437 False
a PM 8006 9564 956 False
a AM 10236 14974 1497 False
a MM 7289 9885 988 False
a RSC 17353 23036 2304 False
a A 17784 22049 2205 False
b neurons 97989 93881 9388 True
b sessions 150 141 14 True
b mice 8 8 0 True
b V1 15497 14373 1437 True
b PM 10313 9564 956 True
b AM 15225 14974 1497 True
b MM 10654 9885 988 True
b RSC 23728 23036 2304 True
b A 22572 22049 2205 True
reading b, area shares: {'V1': 0.158, 'PM': 0.105, 'AM': 0.155, 'MM': 0.109, 'RSC': 0.242, 'A': 0.23}
paper, area shares: {'V1': 0.153, 'PM': 0.102, 'AM': 0.159, 'MM': 0.105, 'RSC': 0.245, 'A': 0.235}
by mouse (reading b):
n_neurons V1 PM AM MM RSC A
mouse
10 5412 1634 1769 256 720 262 771
3 8912 715 1757 1071 820 3575 974
4 9435 2527 457 1022 1433 2595 1401
5 10796 2718 1709 1156 1183 2170 1860
6 11227 2664 579 1375 1703 3778 1128
7 12689 1239 278 1995 1420 4007 3750
8 24975 2080 2038 6201 2257 3557 8842
9 14543 1920 1726 2149 1118 3784 3846
Result: reproduced, resting on an assumption. The files hold 278,155 neurons from 300 sessions of 8 mice, the paper's total. With session fraction correct over all trials (reading b), 150 layer 2/3 sessions above 0.65 hold 97,989 neurons, and every count (sessions, neurons, mice and the six areas) is within its tolerance; the area shares agree with the paper's to within 0.006. With visually guided trials excluded (reading a, the Methods' general rule), 100 sessions and 69,707 neurons remain, and no count but the mice is within tolerance. Reading b is adopted as the assumption, as it is the only one consistent with the paper's counts; the paper's exact list of 141 sessions cannot be recovered from the files.
fig1d-switch-performance¶
Performance is high before a rule switch, typically falls below chance after it, and recovers gradually to high accuracy by the end of the block.
Each block after a rule switch, in sessions above 0.5 correct, is fit with the Methods' sigmoid on its non-visually guided outcomes (helper fit_block_sigmoid, t0 at the block's first trial as the Methods say, with a free midpoint as a variant). Initial and end performance are the fitted values at the block's first and last trial. Means and the end minus initial difference are taken by hierarchical bootstrap over mice, then switches (10,000 samples). Performance before the switch is the fitted end value of the preceding block.
def block_fits(sel, free_t0=False):
rows = []
for _, s in sel.iterrows():
T = res[s.asset_id]["trials"]
sw = list(np.where(T.is_switch.values)[0])
bounds = [0] + sw + [len(T)]
for bi in range(len(bounds) - 1):
blk = T.iloc[bounds[bi]:bounds[bi + 1]]
y = blk[~blk.is_vis].is_correct.values.astype(float)
if len(y) < 10:
continue
ini, end = tseng.fit_block_sigmoid(y, free_t0=free_t0)
rows.append(dict(mouse=s.mouse, asset_id=s.asset_id, block=bi, after_switch=bi > 0,
n=len(y), initial=ini, end=end, mean=y.mean()))
return pd.DataFrame(rows)
t = time.time()
BF = {k: block_fits(sel) for k, sel in SEL50.items()}
BFv = block_fits(SEL50["a"], free_t0=True)
print(f"fits took {time.time() - t:.0f} s")
fits took 177 s
rng = np.random.default_rng(101)
def summarise_switches(bf, label):
sw = bf[bf.after_switch]
pre = bf.set_index(["asset_id", "block"]).end
sw = sw.assign(before=[pre.get((a, b - 1), np.nan) for a, b in zip(sw.asset_id, sw.block)])
g = [x for _, x in sw.groupby("mouse")]
bi = hboot([x.initial.values for x in g], 10000, rng)
be = hboot([x.end.values for x in g], 10000, rng)
bb = hboot([x.before.dropna().values for x in g], 10000, rng)
bd = hboot([(x.end - x.initial).values for x in g], 10000, rng)
print(f"{label}: {len(sw)} switches from {sw.mouse.nunique()} mice")
print(f" before switch (end of previous block) {bb.mean():.3f} ± {bb.std():.3f}")
print(f" initial after switch {bi.mean():.3f} ± {bi.std():.3f}; P(boot >= 0.5) = {(bi >= 0.5).mean():.4f}")
print(f" end of block {be.mean():.3f} ± {be.std():.3f}; P(boot <= 0.7) = {(be <= 0.7).mean():.4f}")
print(f" end minus initial {bd.mean():.3f} ± {bd.std():.3f}, two-tailed p {pb(bd)}")
print(f" switches with initial below 0.5: {(sw.initial < 0.5).mean():.3f}; with end above 0.7: {(sw.end > 0.7).mean():.3f}")
return sw
SW = {}
for k in BF:
SW[k] = summarise_switches(BF[k], f"reading {k} (sessions > 0.5)")
summarise_switches(BFv, "variant, free sigmoid midpoint, reading a")
print("\nper mouse, reading a (mean initial, mean end, switches):")
print(SW["a"].groupby("mouse").agg(initial=("initial", "mean"), end=("end", "mean"), n=("initial", "size")).round(3))
print("\nblocks in reading a:", len(BF["a"]), "| sessions:", BF["a"].asset_id.nunique())
reading a (sessions > 0.5): 589 switches from 8 mice before switch (end of previous block) 0.810 ± 0.020 initial after switch 0.225 ± 0.017; P(boot >= 0.5) = 0.0000 end of block 0.730 ± 0.025; P(boot <= 0.7) = 0.1239 end minus initial 0.504 ± 0.039, two-tailed p < 0.0001 switches with initial below 0.5: 0.919; with end above 0.7: 0.565
reading b (sessions > 0.5): 601 switches from 8 mice before switch (end of previous block) 0.807 ± 0.021 initial after switch 0.226 ± 0.016; P(boot >= 0.5) = 0.0000 end of block 0.727 ± 0.026; P(boot <= 0.7) = 0.1542 end minus initial 0.501 ± 0.038, two-tailed p < 0.0001 switches with initial below 0.5: 0.918; with end above 0.7: 0.561
variant, free sigmoid midpoint, reading a: 589 switches from 8 mice
before switch (end of previous block) 0.809 ± 0.020
initial after switch 0.317 ± 0.015; P(boot >= 0.5) = 0.0000
end of block 0.730 ± 0.026; P(boot <= 0.7) = 0.1235
end minus initial 0.412 ± 0.038, two-tailed p < 0.0001
switches with initial below 0.5: 0.837; with end above 0.7: 0.589
per mouse, reading a (mean initial, mean end, switches):
initial end n
mouse
10 0.257 0.666 98
3 0.292 0.711 49
4 0.253 0.711 63
5 0.191 0.811 62
6 0.193 0.722 90
7 0.265 0.613 68
8 0.207 0.770 77
9 0.172 0.837 82
blocks in reading a: 880 | sessions: 291
Result: reproduced. Over 601 switches in 297 sessions above 0.5 (reading b), performance at the end of the block before a switch is 0.807, the fitted initial performance after the switch is 0.226 ± 0.016 (no bootstrap sample at or above 0.5), and the fitted end performance is 0.727 ± 0.026, 0.501 above the initial value (p < 0.0001). Reading a gives 0.225 and 0.730 over 589 switches, and the free-midpoint variant 0.317 and 0.730. The end value is above 0.7 as a mean but only just: 15% of bootstrap samples (12% in reading a) are at or below 0.7, and in reading a two of the 8 mice end below 0.7 (0.613 and 0.666). The switch count differs from the paper's 506 to 513.
lstm-choice-accuracy¶
The strategy LSTM predicts the mouse's choices in held-out sessions with 80.2% ± 5.1% accuracy (mean ± SD).
The files hold the association matrix of the paper's LSTM for each trial. The tutorial for the dataset and the Methods (LSTM model) describe these as the model's predicted conditional probabilities; that they are the held-out (leave-one-session-out) predictions is an assumption, as the files do not say. A trial is predicted correctly when the probability given to the actual choice for the actual cue exceeds 0.5 (the claims file's reading). Visually guided trials are excluded (Methods, LSTM model). Readings: SD over sessions or over mice; all sessions or those above 0.65.
acc = []
for _, s in D.iterrows():
T = res[s.asset_id]["trials"]
nv = T[~T.is_vis]
acc.append(dict(asset_id=s.asset_id, mouse=s.mouse, fc=s.fc, fc_all=s.fc_all,
acc=(nv.prob_actual_cho > 0.5).mean(), n=len(nv)))
acc = pd.DataFrame(acc)
rows = []
for sname, sub in [("c, all sessions", acc), ("d-a, > 0.65 (vis excluded)", acc[acc.fc > 0.65]),
("d-b, > 0.65 (all trials)", acc[acc.fc_all > 0.65])]:
pm = sub.groupby("mouse").acc.mean()
rows.append(dict(sessions=sname, n_sessions=len(sub), mean_over_sessions=sub.acc.mean(), sd_sessions=sub.acc.std(),
mean_over_mice=pm.mean(), sd_mice=pm.std(),
pooled_trials=(sub.acc * sub.n).sum() / sub.n.sum()))
rows = pd.DataFrame(rows)
print(rows.round(4).to_string(index=False))
print("expected 0.802 ± 0.051; within tolerance:", {r.sessions: (abs(r.mean_over_sessions - 0.802) <= 0.051,
abs(r.mean_over_mice - 0.802) <= 0.051) for r in rows.itertuples()})
print("\nper mouse, all sessions:", acc.groupby("mouse").acc.mean().round(3).to_dict())
sessions n_sessions mean_over_sessions sd_sessions mean_over_mice sd_mice pooled_trials
c, all sessions 300 0.7716 0.0644 0.7711 0.0432 0.7703
d-a, > 0.65 (vis excluded) 154 0.7918 0.0567 0.7797 0.0413 0.7901
d-b, > 0.65 (all trials) 238 0.7689 0.0647 0.7656 0.0474 0.7685
expected 0.802 ± 0.051; within tolerance: {'c, all sessions': (True, True), 'd-a, > 0.65 (vis excluded)': (True, True), 'd-b, > 0.65 (all trials)': (True, True)}
per mouse, all sessions: {'10': 0.74, '3': 0.762, '4': 0.705, '5': 0.826, '6': 0.753, '7': 0.752, '8': 0.813, '9': 0.817}
Result: reproduced, resting on an assumption (that the stored association matrices are the held-out predictions). Accuracy is 0.772 (mean over all 300 sessions, SD 0.064 over sessions; 0.771 over mice, SD 0.043), within 0.051 of 0.802 under every reading; the readings range from 0.766 to 0.792.
figs1h-lstm-beats-rl¶
The LSTM predicts choices better (higher log likelihood) than the best reinforcement learning model (model 4, separate learning rates for rewarded and unrewarded trials).
RL model 4 of the Methods (Reinforcement learning model; Figure S1G legend), helper rl_loglik and rl_fit: Q learning with decay toward Q0 = 0, separate learning rates for rewarded and unrewarded trials, choice perseverance, bias and lapse. Fit by maximum likelihood, leave-one-session-out within each mouse over all its sessions (as for the LSTM), and scored on the held-out session as the mean log2 likelihood of the actual choice over non-visually guided trials. The LSTM's score is the same quantity from the stored association matrix. Wilcoxon signed-rank test over sessions above 0.65 (the paper's 216).
t = time.time()
target = set(SEL["a"].asset_id) | set(SEL["b"].asset_id)
jobs, meta = [], []
for mouse, dm in D.groupby("mouse"):
dm = dm.sort_values("date")
for a in dm.asset_id:
if a not in target:
continue
train = [res[x]["trials"] for x in dm.asset_id if x != a]
jobs.append((tseng.rl_arrays(train), tseng.rl_arrays([res[a]["trials"]])))
meta.append((mouse, a))
fits = tseng.rl_loso_many(jobs, workers=os.cpu_count())
rl_rows = []
for (mouse, a), (th, runs, ll_rl) in zip(meta, fits):
nv = res[a]["trials"][~res[a]["trials"].is_vis]
ll_lstm = np.mean(np.log2(np.clip(nv.prob_actual_cho.values, 1e-12, 1)))
nlls = [r_["nll"] for r_ in runs]
rl_rows.append(dict(mouse=mouse, asset_id=a, rl=ll_rl, lstm=ll_lstm,
best_ok=runs[int(np.argmin(nlls))]["success"], n_ok=sum(r_["success"] for r_ in runs),
nll_spread=max(nlls) - min(nlls), nll_best=min(nlls),
**dict(zip(["a_rew", "a_unrew", "delta", "beta", "b0", "phi", "tau", "lam"], th))))
RL = pd.DataFrame(rl_rows)
print(f"{len(RL)} leave-one-session-out fits in {time.time() - t:.0f} s")
for k in ("a", "b"):
sub = RL[RL.asset_id.isin(SEL[k].asset_id)]
w = stats.wilcoxon(sub.lstm, sub.rl)
print(f"reading {k}: {len(sub)} sessions; LSTM {sub.lstm.mean():.3f}, RL model 4 {sub.rl.mean():.3f} (log2 likelihood per trial); "
f"LSTM higher in {(sub.lstm > sub.rl).sum()} sessions; Wilcoxon p = {w.pvalue:.3g}")
print(" per mouse, LSTM minus RL:", (sub.groupby("mouse").lstm.mean() - sub.groupby("mouse").rl.mean()).round(3).to_dict())
print("fitted parameters (median over folds):", RL[["a_rew", "a_unrew", "delta", "beta", "b0", "phi", "tau", "lam"]].median().round(3).to_dict())
print(f"optimiser: best start reported convergence in {RL.best_ok.sum()} of {len(RL)} folds; all 3 starts converged in "
f"{(RL.n_ok == 3).sum()}; spread of the training negative log likelihood across the 3 starts: median {RL.nll_spread.median():.3f}, "
f"max {RL.nll_spread.max():.3f} (best value: median {RL.nll_best.median():.0f})")
for p_, (lo_, hi_) in zip(["a_rew", "a_unrew", "delta", "beta", "b0", "phi", "tau", "lam"], tseng.RL_BOUNDS):
print(f" {p_}: median {RL[p_].median():.3f}, range {RL[p_].min():.3f} to {RL[p_].max():.3f}; folds at a bound: {((RL[p_] <= lo_ + 1e-6) | (RL[p_] >= hi_ - 1e-6)).sum()}")
238 leave-one-session-out fits in 113 s
reading a: 154 sessions; LSTM -0.667, RL model 4 -0.691 (log2 likelihood per trial); LSTM higher in 125 sessions; Wilcoxon p = 1.42e-17
per mouse, LSTM minus RL: {'10': 0.006, '3': 0.018, '4': 0.01, '5': 0.033, '6': 0.046, '7': 0.026, '8': 0.026, '9': 0.023}
reading b: 238 sessions; LSTM -0.710, RL model 4 -0.732 (log2 likelihood per trial); LSTM higher in 185 sessions; Wilcoxon p = 4.65e-22
per mouse, LSTM minus RL: {'10': 0.009, '3': 0.019, '4': 0.007, '5': 0.034, '6': 0.03, '7': 0.016, '8': 0.026, '9': 0.024}
fitted parameters (median over folds): {'a_rew': 0.158, 'a_unrew': 0.0, 'delta': 0.981, 'beta': 1.88, 'b0': 0.002, 'phi': 0.653, 'tau': 0.29, 'lam': 0.105}
optimiser: best start reported convergence in 238 of 238 folds; all 3 starts converged in 238; spread of the training negative log likelihood across the 3 starts: median 0.000, max 480.021 (best value: median 8262)
a_rew: median 0.158, range 0.004 to 0.305; folds at a bound: 0
a_unrew: median 0.000, range 0.000 to 0.297; folds at a bound: 198
delta: median 0.981, range 0.912 to 0.993; folds at a bound: 0
beta: median 1.880, range 1.476 to 50.000; folds at a bound: 50
b0: median 0.002, range -0.185 to 0.194; folds at a bound: 0
phi: median 0.653, range 0.021 to 1.115; folds at a bound: 0
tau: median 0.290, range 0.184 to 0.534; folds at a bound: 0
lam: median 0.105, range 0.000 to 0.170; folds at a bound: 5
Result: reproduced. On the 238 sessions of reading b the LSTM's held-out log2 likelihood per trial is -0.710 and that of RL model 4 is -0.732; the LSTM is higher in 185 sessions (Wilcoxon p = 4.65e-22). On the 154 sessions of reading a: -0.667 against -0.691, p = 1.42e-17. The LSTM is ahead in all 8 mice, by 0.006 to 0.046. The optimiser reported convergence for the best start in all 238 folds, and the three starts reach the same likelihood in the median fold (spread 0.000; at most 480 of a median 8,262). The learning rate for unrewarded trials sits at its lower bound of 0 in 198 folds and the inverse temperature at its upper bound of 50 in 50 folds, so in those folds the fit is limited by a bound.
Switches used for the switch-aligned strategy analyses¶
Figure 1I and 2C use 265 switches: switches whose following block ends above 70% correct (Methods, Task performance analysis: "Blocks with end performance < 70% correct were viewed as unrecovered switches and excluded"), in sessions above 0.65.
def recovered(k):
bf = block_fits(SEL[k]) if k not in ("a", "b") else BF[k][BF[k].asset_id.isin(SEL[k].asset_id)]
sw = bf[bf.after_switch & (bf.end > 0.7)]
return sw
REC = {k: recovered(k) for k in ("a", "b")}
for k, sw in REC.items():
print(f"reading {k}: {len(sw)} recovered switches from {sw.asset_id.nunique()} sessions, {sw.mouse.nunique()} mice [paper: 265]")
def switch_rows(sw):
"""Yield (mouse, trials table, index of the switch trial, index of the block end) for each switch."""
for _, r in sw.iterrows():
T = res[r.asset_id]["trials"]
sws = list(np.where(T.is_switch.values)[0])
i = sws[r.block - 1]
end = sws[r.block] if r.block < len(sws) else len(T)
yield r.mouse, T, i, end
reading a: 251 recovered switches from 147 sessions, 8 mice [paper: 265] reading b: 320 recovered switches from 205 sessions, 8 mice [paper: 265]
fig1i-strategy-after-switch¶
Before a switch mice show little choice bias and follow the rule; after a switch the probability of the actual choice falls near chance, bias-following rises moderately, and rule-following then recovers gradually.
Strategy variables of the non-visually guided trials, aligned to the switch trial (offset 0 is the first trial of the new rule). Windows of W = 10 and 20 trials before and after the switch (claims file). The end of the block is its last W non-visually guided trials. Comparisons are paired per switch, hierarchical bootstrap over mice then switches, 10,000 samples.
rng = np.random.default_rng(102)
def strat_windows(sw, W):
rows, curves = [], []
for mouse, T, i, end in switch_rows(sw):
nv = T[~T.is_vis]
pos = np.searchsorted(nv.index.values, T.index.values[i]) # first non-vis trial at or after the switch
endpos = np.searchsorted(nv.index.values, T.index.values[end - 1], side="right")
if pos < W or endpos - pos < 2 * W:
continue
b, a = nv.iloc[pos - W:pos], nv.iloc[pos:pos + W]
e = nv.iloc[endpos - W:endpos]
rows.append(dict(mouse=mouse, rf_before=b.rule_following.mean(), rf_after=a.rule_following.mean(),
rf_end=e.rule_following.mean(), bf_before=b.bias_following.mean(),
bf_after=a.bias_following.mean(), pac_before=b.prob_actual_cho.mean(),
pac_after=a.prob_actual_cho.mean(), bias_before=np.abs(b.signed_bias).mean()))
curves.append(nv.prob_actual_cho.iloc[pos - W:pos + W].values)
return pd.DataFrame(rows), np.array(curves)
for k in ("a", "b"):
for W in (10, 20):
R, C = strat_windows(REC[k], W)
g = [x for _, x in R.groupby("mouse")]
def bt(f):
return hboot([f(x).values for x in g], 10000, rng)
rf_b = bt(lambda x: x.rf_before)
d_pac = bt(lambda x: x.pac_after - x.pac_before)
d_bf = bt(lambda x: x.bf_after - x.bf_before)
d_rf = bt(lambda x: x.rf_end - x.rf_after)
cmean = C.mean(0)
pac_min = cmean[W:].min()
print(f"reading {k}, W = {W}: {len(R)} switches, {R.mouse.nunique()} mice")
print(f" (1) rule-following before {rf_b.mean():.3f} ± {rf_b.std():.3f}, p vs 0 {pb(rf_b)}; "
f"choice bias |P(L)-P(R)| before {R.bias_before.mean():.3f}")
print(f" (2) P(actual choice) before {R.pac_before.mean():.3f}, after {R.pac_after.mean():.3f}; "
f"after - before {d_pac.mean():.3f} ± {d_pac.std():.3f}, p {pb(d_pac)}; "
f"minimum of the switch-aligned mean after the switch {pac_min:.3f} (|min - 0.5| = {abs(pac_min - 0.5):.3f})")
print(f" (3) bias-following before {R.bf_before.mean():.3f}, after {R.bf_after.mean():.3f}; "
f"difference {d_bf.mean():.3f} ± {d_bf.std():.3f}, p {pb(d_bf)}")
print(f" (4) rule-following first window after {R.rf_after.mean():.3f}, end of block {R.rf_end.mean():.3f}; "
f"difference {d_rf.mean():.3f} ± {d_rf.std():.3f}, p {pb(d_rf)}; "
f"end of block vs before: {R.rf_end.mean() - R.rf_before.mean():.3f}")
if k == "a" and W == 20:
print(" per mouse:"); print(R.groupby("mouse")[["rf_before", "rf_after", "rf_end", "pac_after", "bf_before", "bf_after"]].mean().round(3))
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(np.arange(-W, W), cmean, "k")
ax.axhline(0.5, ls=":", c="gray"); ax.axvline(0, ls="--", c="r")
ax.set(xlabel="non-visually guided trials from switch", ylabel="P(actual choice)", title="reading a, W = 20")
plt.show()
print(f" plotted curve range {cmean.min():.3f} to {cmean.max():.3f}")
reading a, W = 10: 251 switches, 8 mice (1) rule-following before 0.254 ± 0.021, p vs 0 < 0.0001; choice bias |P(L)-P(R)| before 0.133 (2) P(actual choice) before 0.775, after 0.621; after - before -0.152 ± 0.015, p < 0.0001; minimum of the switch-aligned mean after the switch 0.578 (|min - 0.5| = 0.078) (3) bias-following before 0.017, after 0.054; difference 0.037 ± 0.006, p < 0.0001 (4) rule-following first window after 0.067, end of block 0.220; difference 0.151 ± 0.015, p < 0.0001; end of block vs before: -0.038
reading a, W = 20: 251 switches, 8 mice
(1) rule-following before 0.256 ± 0.020, p vs 0 < 0.0001; choice bias |P(L)-P(R)| before 0.135
(2) P(actual choice) before 0.778, after 0.609; after - before -0.167 ± 0.015, p < 0.0001; minimum of the switch-aligned mean after the switch 0.576 (|min - 0.5| = 0.076)
(3) bias-following before 0.019, after 0.064; difference 0.045 ± 0.006, p < 0.0001
(4) rule-following first window after 0.045, end of block 0.221; difference 0.174 ± 0.015, p < 0.0001; end of block vs before: -0.038
per mouse:
rf_before rf_after rf_end pac_after bf_before bf_after
mouse
10 0.186 0.023 0.142 0.575 0.028 0.052
3 0.218 0.028 0.220 0.607 0.028 0.079
4 0.192 0.030 0.160 0.566 0.005 0.036
5 0.265 0.043 0.219 0.647 0.047 0.104
6 0.227 0.040 0.202 0.610 0.019 0.070
7 0.192 0.022 0.090 0.586 0.006 0.064
8 0.292 0.060 0.257 0.618 0.007 0.058
9 0.307 0.057 0.262 0.606 0.009 0.049
plotted curve range 0.576 to 0.803
reading b, W = 10: 320 switches, 8 mice (1) rule-following before 0.227 ± 0.025, p vs 0 < 0.0001; choice bias |P(L)-P(R)| before 0.158 (2) P(actual choice) before 0.755, after 0.617; after - before -0.136 ± 0.015, p < 0.0001; minimum of the switch-aligned mean after the switch 0.575 (|min - 0.5| = 0.075) (3) bias-following before 0.026, after 0.055; difference 0.029 ± 0.008, p = 0.0004 (4) rule-following first window after 0.062, end of block 0.209; difference 0.146 ± 0.014, p < 0.0001; end of block vs before: -0.020
reading b, W = 20: 319 switches, 8 mice (1) rule-following before 0.228 ± 0.025, p vs 0 < 0.0001; choice bias |P(L)-P(R)| before 0.160 (2) P(actual choice) before 0.757, after 0.607; after - before -0.149 ± 0.017, p < 0.0001; minimum of the switch-aligned mean after the switch 0.570 (|min - 0.5| = 0.070) (3) bias-following before 0.026, after 0.064; difference 0.038 ± 0.007, p < 0.0001 (4) rule-following first window after 0.043, end of block 0.208; difference 0.164 ± 0.015, p < 0.0001; end of block vs before: -0.023
Result: reproduced. Under both selections and both windows: rule-following before the switch is 0.23 to 0.26 (p < 1e-4 against zero) with a mean choice bias of 0.13 to 0.16; the probability of the actual choice falls by 0.14 to 0.17 after the switch (p < 1e-4) and its switch-aligned minimum is 0.570 to 0.578, within 0.1 of 0.5; bias-following rises by 0.029 to 0.045 (p <= 4e-4); rule-following rises from 0.04 to 0.07 in the first window to 0.21 to 0.22 by the end of the block (p < 1e-4), close to its pre-switch level. The recovered switches number 320 (reading b) and 251 (reading a), against the paper's 265.
visually-guided-rule-update¶
Visually guided trials do not update rule belief more than normal trials in the 30 trials after a switch (0.013 ± 0.0055 against 0.0080 ± 0.0023; p = 0.34).
Rule belief is recomputed from the stored association matrix with the Methods' sign (positive for rule B) and signed toward the current (new) rule. The update caused by trial t is the change from the belief on trial t to that on trial t + 1, since the LSTM's belief on a trial is conditioned on the trials before it (Methods, LSTM model); the change from t - 1 to t is shown as a variant. Over the 30 trials from the switch trial, each switch gives the mean update over its visually guided and over its normal trials; the difference is paired per switch (switches with both kinds of trial), hierarchical bootstrap over mice then switches, 10,000 samples.
rng = np.random.default_rng(103)
def vg_update(sw, forward=True, correct_only=False):
rows = []
for mouse, T, i, end in switch_rows(sw):
rb = ((T.P_L_W + T.P_R_B) / 2 - (T.P_L_B + T.P_R_W) / 2).values
sign = np.where(T.is_ruleA.values, -1.0, 1.0)
d = (np.r_[rb[1:], np.nan] - rb) if forward else (rb - np.r_[np.nan, rb[:-1]])
d = d * sign
idx = np.arange(i, min(i + 30, end))
keep = idx if not correct_only else idx[T.is_correct.values[idx]]
vis = T.is_vis.values[keep]
dv, dn = d[keep][vis], d[keep][~vis]
dv, dn = dv[~np.isnan(dv)], dn[~np.isnan(dn)]
if len(dv) and len(dn):
rows.append(dict(mouse=mouse, vis=dv.mean(), normal=dn.mean()))
return pd.DataFrame(rows)
for k in ("a", "b"):
for fwd in (True, False):
for co in (False, True):
R = vg_update(REC[k], fwd, co)
g = [x for _, x in R.groupby("mouse")]
bv = hboot([x.vis.values for x in g], 10000, rng)
bn = hboot([x.normal.values for x in g], 10000, rng)
bd = hboot([(x.vis - x.normal).values for x in g], 10000, rng)
print(f"reading {k}, {'t to t+1' if fwd else 't-1 to t'}, {'correct trials' if co else 'all trials'}: "
f"{len(R)} switches | visually guided {bv.mean():.4f} ± {bv.std():.4f}, normal {bn.mean():.4f} ± {bn.std():.4f}, "
f"difference {bd.mean():.4f} ± {bd.std():.4f}, p {pb(bd)}")
# Split of the forward update by trial type and outcome, and the switch-aligned rule belief (reading b)
rows, curves = [], []
for mouse, T, i, end in switch_rows(REC["b"]):
rb = ((T.P_L_W + T.P_R_B) / 2 - (T.P_L_B + T.P_R_W) / 2).values
sign = np.where(T.is_ruleA.values, -1.0, 1.0)
d = (np.r_[rb[1:], np.nan] - rb) * sign
idx = np.arange(i, min(i + 30, end))
for vis in (True, False):
for cor in (True, False):
m = (T.is_vis.values[idx] == vis) & (T.is_correct.values[idx] == cor)
v = d[idx][m]
v = v[~np.isnan(v)]
if len(v):
rows.append(dict(mouse=mouse, visually_guided=vis, correct=cor, d=v.mean(), n=len(v)))
if i >= 10 and end - i >= 30:
curves.append(rb[i - 10:i + 30] * sign[i])
S5 = pd.DataFrame(rows)
print("\nforward update (t to t+1) toward the new rule over the 30 trials from the switch, reading b, by trial type and outcome:")
print(S5.groupby(["visually_guided", "correct"]).agg(mean_update=("d", "mean"), switches=("d", "size"), trials=("n", "sum")).round(4))
CV = np.array(curves).mean(0)
print(f"switch-aligned rule belief toward the new rule ({len(curves)} switches): trial -10 {CV[0]:.3f}, -1 {CV[9]:.3f}, "
f"0 {CV[10]:.3f}, +10 {CV[20]:.3f}, +20 {CV[30]:.3f}, +29 {CV[39]:.3f}; mean change per trial over trials 0 to 29 "
f"{(CV[39] - CV[10]) / 29:.4f}")
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(np.arange(-10, 30), CV, "k"); ax.axvline(0, c="r", ls="--"); ax.axhline(0, c="gray", ls=":")
ax.set(xlabel="trials from switch", ylabel="rule belief toward new rule", title="reading b")
plt.show()
print(f"plotted range {CV.min():.3f} to {CV.max():.3f}")
reading a, t to t+1, all trials: 233 switches | visually guided 0.0591 ± 0.0079, normal 0.0236 ± 0.0028, difference 0.0351 ± 0.0084, p = 0.0022
reading a, t to t+1, correct trials: 231 switches | visually guided 0.0743 ± 0.0122, normal 0.0802 ± 0.0101, difference -0.0058 ± 0.0069, p = 0.375
reading a, t-1 to t, all trials: 233 switches | visually guided 0.0309 ± 0.0060, normal 0.0287 ± 0.0032, difference 0.0022 ± 0.0053, p = 0.683
reading a, t-1 to t, correct trials: 231 switches | visually guided 0.0295 ± 0.0056, normal 0.0272 ± 0.0038, difference 0.0023 ± 0.0064, p = 0.737
reading b, t to t+1, all trials: 302 switches | visually guided 0.0545 ± 0.0075, normal 0.0215 ± 0.0029, difference 0.0329 ± 0.0074, p < 0.0001
reading b, t to t+1, correct trials: 300 switches | visually guided 0.0682 ± 0.0112, normal 0.0760 ± 0.0090, difference -0.0076 ± 0.0062, p = 0.218
reading b, t-1 to t, all trials: 302 switches | visually guided 0.0295 ± 0.0047, normal 0.0267 ± 0.0033, difference 0.0029 ± 0.0041, p = 0.481
reading b, t-1 to t, correct trials: 300 switches | visually guided 0.0281 ± 0.0047, normal 0.0258 ± 0.0035, difference 0.0022 ± 0.0056, p = 0.7
forward update (t to t+1) toward the new rule over the 30 trials from the switch, reading b, by trial type and outcome:
mean_update switches trials
visually_guided correct
False False -0.0243 320 4100
True 0.0765 320 3475
True False -0.0129 192 412
True 0.0690 300 1613
switch-aligned rule belief toward the new rule (320 switches): trial -10 -0.624, -1 -0.623, 0 -0.619, +10 -0.140, +20 0.083, +29 0.195; mean change per trial over trials 0 to 29 0.0281
plotted range -0.625 to 0.195
Result: not reproduced, resting on an assumption (the update of trial t is the change from trial t to t + 1). Under this reading visually guided trials update rule belief more than normal trials over the 30 trials after a switch: 0.0545 ± 0.0075 against 0.0215 ± 0.0029, difference 0.0329 ± 0.0074, p < 0.0001 (302 switches, reading b; reading a 0.0591 against 0.0236, p = 0.0022), so the paper's test finds the effect the claim says is absent. The change from t - 1 to t finds no difference (0.0295 against 0.0267, p = 0.48), but it is not the adopted reading and is not counted as support; under every all-trial reading both values are 2 to 5 times the paper's 0.013 and 0.0080 and outside their tolerances (the correct-only rows, which have no paper values, are larger still: 0.076 to 0.080 for normal trials). The split by outcome shows where the forward difference comes from: correct trials move belief toward the new rule (0.0765 normal, 0.0690 visually guided) and incorrect trials away from it (-0.0243 and -0.0129), and 80% of the visually guided trials are correct (1,613 of 2,025) against 46% of normal trials (3,475 of 7,575). On correct trials alone there is no difference (0.068 against 0.076, p = 0.22), as in the paper. The magnitude gap is in the stored belief itself: the switch-aligned rule belief rises from -0.619 at the switch to 0.195 at trial 29, 0.028 per trial on average, 2 to 3.5 times the paper's per-trial updates.
Photoinhibition claims (fig1l-photoinhibition-performance, fig1m-rule-following-decrease, figs1m-rule-bias-uncorrelated)¶
The dandiset's 308 files are 300 imaging sessions and 8 widefield retinotopy files (listed in the data check); none is a photoinhibition session of the VGAT-ChR2 mice, and no trial table has an inhibition column. These claims are blocked for missing data.
fig2c-running-variability-after-switch¶
Running trajectories of choice-matched trials become less correlated with one another after a rule switch.
Methods, Running trajectory correlation analysis: pitch, roll and yaw velocity of each non-visually guided trial averaged in a grid of maze positions (20 and 11 equal bins of the maze, the claims file's readings; helper load_session); Pearson correlation between two trials per velocity, averaged over the three; each trial's mean correlation with all other trials of the same choice, minus the session's mean over all choice-matched pairs (Δ trajectory correlation); mean over the 20 trials before and the 20 from the switch trial. Trials with an empty position bin are dropped. Recovered switches as above; hierarchical bootstrap, mice then switches, 10,000 samples.
rng = np.random.default_rng(104)
def delta_traj(a, nb):
T = res[a]["trials"]
v = res[a]["velbin"][nb]
keep = (~T.is_vis.values) & ~np.isnan(v).any(axis=(1, 2))
idx = np.where(keep)[0]
V = v[idx]
Z = (V - V.mean(2, keepdims=True)) / V.std(2, keepdims=True)
C = np.mean([Z[:, c] @ Z[:, c].T / nb for c in range(3)], axis=0)
cho = T.is_choL.values[idx]
same = cho[:, None] == cho[None, :]
np.fill_diagonal(same, False)
with np.errstate(invalid="ignore"):
per = np.nansum(np.where(same, C, 0), 1) / same.sum(1)
base = C[np.triu(same, 1)].mean()
out = pd.Series(np.nan, index=T.index)
out.iloc[idx] = per - base
return out, (~keep & ~T.is_vis.values).sum()
for k in ("a", "b"):
for nb in (20, 11):
cache, dropped, rows, curves = {}, 0, [], []
for mouse, T, i, end in switch_rows(REC[k]):
a = [x for x in REC[k].asset_id.unique() if res[x]["trials"] is T][0]
if a not in cache:
cache[a], n_drop = delta_traj(a, nb)
dropped += n_drop
d = cache[a]
nv = d[~T.is_vis]
pos = np.searchsorted(nv.index.values, T.index.values[i])
if pos < 20 or len(nv) - pos < 20:
continue
b, af = nv.iloc[pos - 20:pos], nv.iloc[pos:pos + 20]
if b.notna().sum() == 0 or af.notna().sum() == 0:
continue
rows.append(dict(mouse=mouse, before=b.mean(), after=af.mean()))
curves.append(nv.iloc[pos - 20:pos + 20].values)
R = pd.DataFrame(rows)
g = [x for _, x in R.groupby("mouse")]
bd = hboot([(x.after - x.before).values for x in g], 10000, rng)
print(f"reading {k}, {nb} bins: {len(R)} switches; dropped {dropped} non-visually guided trials with an empty bin; "
f"before {R.before.mean():.4f}, after {R.after.mean():.4f}, after - before {bd.mean():.4f} ± {bd.std():.4f}, "
f"p {pb(bd)}; mice with after < before: {(R.groupby('mouse').apply(lambda x: (x.after - x.before).mean()) < 0).sum()} of {R.mouse.nunique()}")
if k == "a" and nb == 20:
cm = np.nanmean(np.array(curves), 0)
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(np.arange(-20, 20), cm, "k"); ax.axvline(0, c="r", ls="--"); ax.axhline(0, c="gray", ls=":")
ax.set(xlabel="non-visually guided trials from switch", ylabel="Δ trajectory correlation", title="reading a, 20 bins")
plt.show()
print(f"plotted range {np.nanmin(cm):.4f} to {np.nanmax(cm):.4f}")
reading a, 20 bins: 248 switches; dropped 2199 non-visually guided trials with an empty bin; before 0.0382, after -0.0498, after - before -0.0869 ± 0.0120, p < 0.0001; mice with after < before: 8 of 8
plotted range -0.0863 to 0.0533
reading a, 11 bins: 251 switches; dropped 1066 non-visually guided trials with an empty bin; before 0.0397, after -0.0519, after - before -0.0903 ± 0.0127, p < 0.0001; mice with after < before: 8 of 8
reading b, 20 bins: 317 switches; dropped 2556 non-visually guided trials with an empty bin; before 0.0376, after -0.0372, after - before -0.0739 ± 0.0118, p < 0.0001; mice with after < before: 8 of 8
reading b, 11 bins: 320 switches; dropped 1412 non-visually guided trials with an empty bin; before 0.0388, after -0.0388, after - before -0.0766 ± 0.0127, p < 0.0001; mice with after < before: 8 of 8
Result: reproduced. Δ trajectory correlation is lower in the 20 trials after a switch than in the 20 before, under both bin readings and both selections: after minus before -0.077 ± 0.013 (11 bins) and -0.074 ± 0.012 (20 bins) in reading b, -0.090 and -0.087 in reading a, p < 1e-4 in each, and lower in all 8 mice.
fig2f-latency-strategy¶
Dynamic choice crosses threshold sooner on trials with higher rule-following or bias-following.
Methods, Modeling of dynamic choice: latency is the time from trial onset at which dynamic choice (the stored velocity LSTM output, P left) reaches 0.9 on left-choice trials or 0.1 on right-choice trials, divided by the session's mean trial duration. Trial onset is the trial's start time in the trial table. The trial duration of the legend (8.95 ± 2.04 s) is compared below with the maze epoch and with the whole trial (start to stop); the one that matches is used. Non-visually guided trials of sessions above 0.65. Readings: 5 and 10 quantile bins of each strategy variable (pooled over trials), and the Spearman correlation over trials. Highest against lowest bin by hierarchical bootstrap, mice then trials, 10,000 samples per bin, compared as two unpaired samples.
def latencies(sel):
rows = []
for _, s in sel.iterrows():
r = res[s.asset_id]
T = r["trials"]
q = r["rnn"].dropna(subset=["trial", "dyn_choice"])
q = q.assign(trial=q.trial.astype(int))
st = T.start_time.reindex(q.trial.values).values
q = q.assign(dt=q.t.values - st)
q = q[q.dt >= 0]
left = T.is_choL.reindex(q.trial.values).values
hit = np.where(left, q.dyn_choice.values >= 0.9, q.dyn_choice.values <= 0.1)
first = q[hit].groupby("trial").dt.min()
mdur, tdur = T.maze_dur.mean(), T.trial_dur.mean()
nv = T[~T.is_vis]
rows.append(pd.DataFrame(dict(mouse=s.mouse, asset_id=s.asset_id, lat=first.reindex(nv.index).values,
maze_dur_s=mdur, trial_dur_s=tdur, rf=nv.rule_following.values,
bf=nv.bias_following.values, has_rnn=nv.index.isin(q.trial.unique()))))
return pd.concat(rows, ignore_index=True)
LAT = {k: latencies(SEL[k]) for k in ("a", "b")}
rng = np.random.default_rng(105)
for k, L in LAT.items():
sd = L.groupby("asset_id")[["maze_dur_s", "trial_dur_s"]].first()
print(f"reading {k}: session-averaged maze duration {sd.maze_dur_s.mean():.2f} ± {sd.maze_dur_s.std():.2f} s; "
f"whole trial {sd.trial_dur_s.mean():.2f} ± {sd.trial_dur_s.std():.2f} s (expected 8.95 ± 2.04)")
L = L[L.has_rnn]
n0 = len(L)
L = L.dropna(subset=["lat"]).copy()
L["latn"] = L.lat / L.maze_dur_s
print(f" trials with LSTM output {n0}; crossing threshold {len(L)} (paper: 85,463)")
for var in ("rf", "bf"):
rho, p = stats.spearmanr(L[var], L.latn)
line = f" {var}: Spearman rho {rho:.3f} (p = {p:.2g})"
for nq in (5, 10):
L["bin"] = pd.qcut(L[var], nq, labels=False)
lo, hi = L[L.bin == 0], L[L.bin == nq - 1]
blo = hboot([x.latn.values for _, x in lo.groupby("mouse")], 10000, rng)
bhi = hboot([x.latn.values for _, x in hi.groupby("mouse")], 10000, rng)
bins = [round(float(x), 3) for x in L.groupby("bin").latn.mean().values]
line += (f" | {nq} bins: lowest {blo.mean():.3f}, highest {bhi.mean():.3f}, p {pu(bhi, blo)}, "
f"bin means {list(bins)}")
print(line)
per = L.groupby("mouse").apply(lambda x: stats.spearmanr(x.rf, x.latn)[0]).round(3)
perb = L.groupby("mouse").apply(lambda x: stats.spearmanr(x.bf, x.latn)[0]).round(3)
print(" per mouse Spearman rho, rule-following:", per.to_dict(), "\n bias-following:", perb.to_dict())
reading a: session-averaged maze duration 8.98 ± 1.93 s; whole trial 14.52 ± 1.98 s (expected 8.95 ± 2.04) trials with LSTM output 49456; crossing threshold 46708 (paper: 85,463)
rf: Spearman rho -0.139 (p = 1.3e-199) | 5 bins: lowest 0.667, highest 0.557, p = 0.0137, bin means [0.666, 0.603, 0.623, 0.608, 0.55] | 10 bins: lowest 0.700, highest 0.541, p = 0.00129, bin means [0.698, 0.635, 0.6, 0.605, 0.616, 0.63, 0.618, 0.599, 0.571, 0.528]
bf: Spearman rho -0.071 (p = 3.2e-53) | 5 bins: lowest 0.656, highest 0.587, p = 0.0576, bin means [0.655, 0.601, 0.598, 0.61, 0.586] | 10 bins: lowest 0.668, highest 0.568, p = 0.01, bin means [0.667, 0.642, 0.611, 0.592, 0.59, 0.606, 0.616, 0.605, 0.606, 0.566]
per mouse Spearman rho, rule-following: {'10': -0.053, '3': -0.105, '4': -0.057, '5': -0.067, '6': -0.103, '7': 0.146, '8': -0.128, '9': -0.115}
bias-following: {'10': -0.068, '3': -0.052, '4': -0.036, '5': -0.082, '6': -0.176, '7': -0.279, '8': -0.037, '9': -0.061}
reading b: session-averaged maze duration 9.31 ± 2.14 s; whole trial 14.92 ± 2.18 s (expected 8.95 ± 2.04)
trials with LSTM output 74425; crossing threshold 70410 (paper: 85,463)
rf: Spearman rho -0.113 (p = 4.5e-198) | 5 bins: lowest 0.673, highest 0.571, p = 0.00453, bin means [0.672, 0.589, 0.624, 0.627, 0.566] | 10 bins: lowest 0.701, highest 0.552, p = 5.85e-05, bin means [0.701, 0.643, 0.585, 0.593, 0.618, 0.629, 0.631, 0.624, 0.591, 0.541]
bf: Spearman rho -0.102 (p = 3.8e-163) | 5 bins: lowest 0.666, highest 0.571, p = 0.00113, bin means [0.666, 0.613, 0.614, 0.617, 0.569] | 10 bins: lowest 0.679, highest 0.541, p = 3.87e-05, bin means [0.679, 0.653, 0.625, 0.6, 0.604, 0.623, 0.618, 0.615, 0.601, 0.538]
per mouse Spearman rho, rule-following: {'10': -0.078, '3': -0.115, '4': -0.074, '5': -0.069, '6': -0.054, '7': 0.083, '8': -0.073, '9': -0.115}
bias-following: {'10': -0.096, '3': -0.051, '4': -0.074, '5': -0.1, '6': -0.187, '7': -0.198, '8': -0.089, '9': -0.071}
Result: reproduced, resting on an assumption (the session selection, reading b). The maze epoch matches the legend's trial duration (9.31 ± 2.14 s in reading b, 8.98 ± 1.93 s in reading a, against 8.95 ± 2.04), and trial start is maze onset (the check after the data check: median offset -0.004 s). In reading b (70,410 trials crossing threshold) normalised latency falls with rule-following (Spearman rho -0.113; highest against lowest bin 0.571 against 0.673 with 5 bins, p = 0.0045, and 0.552 against 0.701 with 10 bins, p < 0.0001) and with bias-following (rho -0.102; p = 0.0011 and p < 0.0001). In reading a every comparison holds except bias-following with 5 bins (0.587 against 0.656, p = 0.058); the hierarchical bootstrap with 8 mice is liberal near 0.05 (verifier's synthetic check), so this p value is fragile either way. One mouse of 8 has a positive rho for rule-following in both readings.
fig2g-choice-commitment¶
With high rule-following, choice commitment starts near chance and rises during the maze; with low rule-following it stays low until late; with high bias-following it is high from trial onset.
Choice commitment is the log likelihood (base 2) of the reported choice under the stored dynamic choice, averaged per trial in each maze position bin (11 or 20 equal bins of the maze, as in the claims file). Tested as in Figure S2L: per session and bin, least squares regression of choice commitment on rule-following, bias-following and local reward rate (fraction correct over the 9 adjacent non-visually guided trials, centred on the trial), with an intercept; coefficients averaged by hierarchical bootstrap over mice then sessions (10,000 samples). The last stem bin is the last bin whose centre lies in the stem (stem length from each file).
def commitment(sel, nb, exclude_current=False):
rows = []
for _, s in sel.iterrows():
r = res[s.asset_id]
T = r["trials"]
nv = T[~T.is_vis]
x = nv.is_correct.astype(float)
r9 = x.rolling(9, center=True, min_periods=5)
# centred window of 9 including the trial itself; the variant drops the trial (its 8 neighbours)
lrr = (r9.sum() - x) / (r9.count() - 1) if exclude_current else r9.mean()
q = r["rnn"].dropna(subset=["trial", "dyn_choice", "posF"])
q = q[q.trial.isin(nv.index)]
y = nv.is_choL.reindex(q.trial.astype(int).values).values.astype(float)
p = np.clip(q.dyn_choice.values, 1e-6, 1 - 1e-6)
llk = y * np.log2(p) + (1 - y) * np.log2(1 - p)
b = np.clip((q.posF.values / r["maze_len"] * nb).astype(int), 0, nb - 1)
M = pd.DataFrame(dict(trial=q.trial.astype(int).values, bin=b, ll=llk)).groupby(["trial", "bin"]).ll.mean().unstack()
centres = (np.arange(nb) + 0.5) * r["maze_len"] / nb
last_stem = int(np.where(centres < r["stem"])[0].max())
X = pd.DataFrame(dict(rf=nv.rule_following, bf=nv.bias_following, lrr=lrr)).reindex(M.index)
for bi in [0, last_stem] + list(range(nb)):
yb = M[bi] if bi in M else None
if yb is None:
continue
ok = yb.notna() & X.notna().all(1)
if ok.sum() < 20:
continue
A = np.c_[np.ones(ok.sum()), X[ok].values]
coef = np.linalg.lstsq(A, yb[ok].values, rcond=None)[0]
rows.append(dict(mouse=s.mouse, asset_id=s.asset_id, bin=bi, rel=bi / (nb - 1),
is_first=bi == 0, is_last_stem=bi == last_stem, c_rf=coef[1], c_bf=coef[2], c_lrr=coef[3]))
return pd.DataFrame(rows).drop_duplicates(["asset_id", "bin"])
rng = np.random.default_rng(106)
for k, nb, excl in [("b", 11, False), ("b", 20, False), ("a", 11, False), ("a", 20, False), ("b", 11, True), ("a", 11, True)]:
if True:
C = commitment(SEL[k], nb, excl)
f0 = C[C.is_first].set_index("asset_id")
ls = C[C.is_last_stem].set_index("asset_id")
j = f0.join(ls[["c_rf", "c_bf"]], rsuffix="_last").dropna()
g = [x for _, x in j.groupby("mouse")]
b_bf0 = hboot([x.c_bf.values for x in g], 10000, rng)
b_rf0 = hboot([x.c_rf.values for x in g], 10000, rng)
b_rfl = hboot([x.c_rf_last.values for x in g], 10000, rng)
b_d = hboot([(x.c_rf_last - x.c_rf).values for x in g], 10000, rng)
print(f"reading {k}, {nb} bins, local reward rate {'without' if excl else 'with'} the current trial, {len(j)} sessions:")
print(f" first bin: bias-following coef {b_bf0.mean():.3f} ± {b_bf0.std():.3f} (p vs 0 {pb(b_bf0)}); "
f"rule-following coef {b_rf0.mean():.3f} ± {b_rf0.std():.3f} (p vs 0 {pb(b_rf0)}); "
f"|rule| / bias = {abs(b_rf0.mean()) / b_bf0.mean():.3f}")
print(f" last stem bin: rule-following coef {b_rfl.mean():.3f} ± {b_rfl.std():.3f} (p vs 0 {pb(b_rfl)}); "
f"last minus first {b_d.mean():.3f} ± {b_d.std():.3f} (p {pb(b_d)})")
if k == "b" and nb == 11 and not excl:
mb = C.groupby("bin")[["c_rf", "c_bf", "c_lrr"]].mean()
print(" mean coefficients by bin:\n", mb.round(3).T)
fig, ax = plt.subplots(figsize=(5, 3))
for c, col in zip(["c_rf", "c_bf", "c_lrr"], ["b", "g", "crimson"]):
ax.plot(mb.index, mb[c], color=col, label=c)
ax.axhline(0, c="gray", ls=":"); ax.legend(); ax.set(xlabel="maze position bin (11), reading b", ylabel="coefficient")
plt.show()
print(f" plotted range {mb.values.min():.3f} to {mb.values.max():.3f}")
reading b, 11 bins, local reward rate with the current trial, 238 sessions: first bin: bias-following coef 0.415 ± 0.064 (p vs 0 < 0.0001); rule-following coef -0.007 ± 0.024 (p vs 0 = 0.697); |rule| / bias = 0.017 last stem bin: rule-following coef 0.850 ± 0.054 (p vs 0 < 0.0001); last minus first 0.856 ± 0.055 (p < 0.0001) mean coefficients by bin: bin 0 1 2 3 4 5 6 7 8 9 \ c_rf -0.007 0.127 0.273 0.384 0.487 0.598 0.725 0.885 0.844 0.516 c_bf 0.417 0.551 0.574 0.576 0.614 0.696 0.842 1.002 0.819 0.388 c_lrr -0.043 0.024 0.067 0.094 0.121 0.131 0.142 0.138 0.104 0.040 bin 10 c_rf 0.399 c_bf 0.267 c_lrr 0.012
plotted range -0.043 to 1.002
reading b, 20 bins, local reward rate with the current trial, 236 sessions: first bin: bias-following coef 0.348 ± 0.073 (p vs 0 < 0.0001); rule-following coef -0.013 ± 0.015 (p vs 0 = 0.366); |rule| / bias = 0.037 last stem bin: rule-following coef 0.863 ± 0.053 (p vs 0 < 0.0001); last minus first 0.876 ± 0.053 (p < 0.0001)
reading a, 11 bins, local reward rate with the current trial, 154 sessions: first bin: bias-following coef 0.339 ± 0.071 (p vs 0 < 0.0001); rule-following coef -0.000 ± 0.025 (p vs 0 = 0.95); |rule| / bias = 0.000 last stem bin: rule-following coef 0.892 ± 0.073 (p vs 0 < 0.0001); last minus first 0.891 ± 0.077 (p < 0.0001)
reading a, 20 bins, local reward rate with the current trial, 152 sessions: first bin: bias-following coef 0.280 ± 0.085 (p vs 0 = 0.001); rule-following coef -0.013 ± 0.017 (p vs 0 = 0.428); |rule| / bias = 0.048 last stem bin: rule-following coef 0.900 ± 0.074 (p vs 0 < 0.0001); last minus first 0.914 ± 0.075 (p < 0.0001)
reading b, 11 bins, local reward rate without the current trial, 238 sessions: first bin: bias-following coef 0.416 ± 0.063 (p vs 0 < 0.0001); rule-following coef -0.016 ± 0.025 (p vs 0 = 0.488); |rule| / bias = 0.038 last stem bin: rule-following coef 0.856 ± 0.050 (p vs 0 < 0.0001); last minus first 0.872 ± 0.052 (p < 0.0001)
reading a, 11 bins, local reward rate without the current trial, 154 sessions: first bin: bias-following coef 0.342 ± 0.071 (p vs 0 < 0.0001); rule-following coef -0.006 ± 0.026 (p vs 0 = 0.767); |rule| / bias = 0.017 last stem bin: rule-following coef 0.895 ± 0.069 (p vs 0 < 0.0001); last minus first 0.900 ± 0.075 (p < 0.0001)
Result: reproduced. In the first position bin the bias-following coefficient is above zero (0.28 to 0.42, p <= 0.001) and the rule-following coefficient is not different from zero (-0.016 to 0.000, p >= 0.36, and at most 0.05 of the bias coefficient); in the last stem bin the rule-following coefficient is 0.85 to 0.90 (p < 0.0001), above its first-bin value (p < 0.0001). This holds for 11 and 20 bins, both selections, and with the local reward rate taken with or without the current trial.
cue-biased-running-uncorrelated¶
Cue-biased running varies across mice and sessions and is not correlated with task performance or with the time course of dynamic choice.
Strength of cue-biased running is the session-averaged log likelihood (base 2) of the cue under the stored cue-biased running, over the maze samples of non-visually guided trials. Readings from the claims file: (a) Pearson correlation of session cue decoding against session fraction correct, pooled over sessions; (b) per-mouse regression slopes against zero across mice (bootstrap over mice, the Methods' default test, with a one-sample t-test as a variant). The logistic cue decoder of Figure S2K (mean velocities over the first 25% of the stem, 10-fold cross-validation, Methods, Logistic decoder for cue from movement) is fit per session; the sign of its roll coefficient is counted per mouse. Sessions above 0.65 under each selection reading; the fraction correct regressed on is the one that selection uses.
Result: partially reproduced. Cue-biased running varies across mice (mean session cue log likelihood from -0.92 to -0.36 by mouse) and in direction across and within mice (roll coefficient positive in 79 of 238 sessions; mouse 5 negative in all 30, mouse 10 positive in 27 and negative in 11). Per-mouse regressions of performance on cue decoding find no relation (mean slope 0.024, bootstrap p = 0.35, reading b; -0.048, p = 0.35, reading a), but the pooled correlation over sessions is positive (r = 0.224, p = 0.0005, reading b; r = 0.207, p = 0.010, reading a). The two readings the claims file names disagree, so the claim holds in part. In reading b choice decoding varies across the 7 bins of cue decoding by 0.03 and 0.04 in the first and last quadrants of the maze and by 0.09 and 0.17 in the second and third. In the second it rises from -0.817 to -0.726 over the first four bins and stays at -0.768 to -0.776 in the last three; in the third it rises from -0.675 to -0.509 over the seven bins, with a dip from the fourth to the sixth bin (-0.547 to -0.559). This part is reported, not judged.
fig4h-activity-by-epoch¶
RSC is most active in the maze stem and least active in the feedback period/ITI, when area A is more active.
Deconvolved activity of each neuron z-scored over its session (the claims file's standard reading), averaged in the maze stem (position below the stem length), the maze arm (position at or beyond the stem length, up to the choice point), and the feedback period/ITI (frames after the choice point, which have no maze position: the 1 s feedback delay falls here, as the check after the data check shows), helper _epoch_activity. Neurons of the layer 2/3 sessions above 0.65 under each selection reading. The claims file's second reading of the neurons (the 42,998 well-fit neurons) needs the GLM fits, which are not run here (see below), so only the all-neuron reading is run. Comparisons by hierarchical bootstrap over mice then neurons, 1,000 samples: between areas as two unpaired samples, within RSC as paired differences.
rng = np.random.default_rng(108)
def neuron_table(sel):
out = []
for _, s in sel[sel.layer == "L23"].iterrows():
r = res[s.asset_id]
act, neu = r["act"], r["neurons"]
assert len(act) == len(neu) and (act.plane.values == neu.plane.values).all()
out.append(pd.concat([neu.reset_index(drop=True), act.drop(columns="plane")], axis=1).assign(mouse=s.mouse))
return pd.concat(out, ignore_index=True)
for k in ("a", "b"):
N = neuron_table(SEL[k])
n0 = len(N)
N = N.dropna(subset=["stem", "arm", "fb_iti"])
print(f"reading {k}: {n0} neurons, {n0 - len(N)} dropped (no frames in an epoch or constant trace)")
B = {}
for a in AREAS:
sub = N[N.area == a]
B[a] = hboot_multi([x[["stem", "arm", "fb_iti"]].values for _, x in sub.groupby("mouse")], 1000, rng)
tab = pd.DataFrame({a: B[a].mean(0) for a in AREAS}, index=["stem", "arm", "fb_iti"]).T
fr = pd.concat([res[a]["act"][["plane", "n_stem", "n_arm", "n_fb", "ts_ratio"]].drop_duplicates("plane")
for a in SEL[k][SEL[k].layer == "L23"].asset_id])
print(f" frames per imaging plane and session: stem {fr.n_stem.min()} to {fr.n_stem.max()} (median {fr.n_stem.median():.0f}), "
f"arm {fr.n_arm.min()} to {fr.n_arm.max()} (median {fr.n_arm.median():.0f}), feedback/ITI {fr.n_fb.min()} to {fr.n_fb.max()} "
f"(median {fr.n_fb.median():.0f}); arm / stem frames {(fr.n_arm / fr.n_stem).median():.3f}")
sp = fr[fr.ts_ratio > 100]
print(f" planes whose activity timestamps run {sp.ts_ratio.median():.0f} times the behaviour clock (matched by frame index): {len(sp)}; "
f"planes with no stem frames: {(fr.n_stem == 0).sum()}")
print(tab.round(4))
print(" (a) RSC stem against each other area's stem:",
{a: f"{B['RSC'][:, 0].mean() - B[a][:, 0].mean():+.4f}, p {pu(B['RSC'][:, 0], B[a][:, 0])}" for a in AREAS if a != "RSC"})
print(f" RSC feedback/ITI {B['RSC'][:, 2].mean():.4f} against A {B['A'][:, 2].mean():.4f}, p {pu(B['RSC'][:, 2], B['A'][:, 2])}")
R = N[N.area == "RSC"]
g = [x for _, x in R.groupby("mouse")]
d1 = hboot([(x.stem - x.arm).values for x in g], 1000, rng)
d2 = hboot([(x.stem - x.fb_iti).values for x in g], 1000, rng)
d3 = hboot([(x.arm - x.fb_iti).values for x in g], 1000, rng)
print(f" (b) within RSC: stem - arm {d1.mean():+.4f} (p {pb(d1)}); stem - feedback/ITI {d2.mean():+.4f} "
f"(p {pb(d2)}); arm - feedback/ITI {d3.mean():+.4f} (p {pb(d3)})")
if k == "b":
print(" per mouse, RSC means:")
print(R.groupby("mouse")[["stem", "arm", "fb_iti"]].mean().round(3))
print(" per mouse, A feedback/ITI minus RSC feedback/ITI:",
(N[N.area == "A"].groupby("mouse").fb_iti.mean() - R.groupby("mouse").fb_iti.mean()).round(3).to_dict())
reading a: 69707 neurons, 10 dropped (no frames in an epoch or constant trace)
frames per imaging plane and session: stem 11709 to 95390 (median 16246), arm 3284 to 30425 (median 4940), feedback/ITI 11873 to 95703 (median 13259); arm / stem frames 0.297
planes whose activity timestamps run 2000 times the behaviour clock (matched by frame index): 5; planes with no stem frames: 0
stem arm fb_iti
V1 0.0041 0.0540 -0.0270
PM -0.0010 0.0656 -0.0258
AM 0.0139 0.0209 -0.0239
MM 0.0031 0.0247 -0.0135
RSC 0.0208 0.0368 -0.0369
A -0.0174 -0.0076 0.0239
(a) RSC stem against each other area's stem: {'V1': '+0.0167, p = 0.113', 'PM': '+0.0218, p = 0.0383', 'AM': '+0.0069, p = 0.432', 'MM': '+0.0177, p = 0.0762', 'A': '+0.0382, p = 4e-06'}
RSC feedback/ITI -0.0369 against A 0.0239, p < 0.001
(b) within RSC: stem - arm -0.0166 (p = 0.178); stem - feedback/ITI +0.0579 (p < 0.001); arm - feedback/ITI +0.0743 (p < 0.001) reading b: 97989 neurons, 14 dropped (no frames in an epoch or constant trace)
frames per imaging plane and session: stem 9985 to 125947 (median 16286), arm 3284 to 54295 (median 5288), feedback/ITI 11873 to 95703 (median 13530); arm / stem frames 0.312
planes whose activity timestamps run 2000 times the behaviour clock (matched by frame index): 8; planes with no stem frames: 0
stem arm fb_iti
V1 0.0011 0.0676 -0.0312
PM 0.0005 0.0724 -0.0304
AM 0.0141 0.0233 -0.0245
MM 0.0023 0.0278 -0.0136
RSC 0.0230 0.0417 -0.0415
A -0.0174 -0.0077 0.0229
(a) RSC stem against each other area's stem: {'V1': '+0.0219, p = 0.0116', 'PM': '+0.0225, p = 0.0136', 'AM': '+0.0089, p = 0.287', 'MM': '+0.0207, p = 0.016', 'A': '+0.0404, p = 4e-06'}
RSC feedback/ITI -0.0415 against A 0.0229, p < 0.001
(b) within RSC: stem - arm -0.0186 (p = 0.07); stem - feedback/ITI +0.0642 (p < 0.001); arm - feedback/ITI +0.0829 (p < 0.001)
per mouse, RSC means:
stem arm fb_iti
mouse
10 0.022 0.037 -0.056
3 0.021 -0.011 -0.024
4 -0.003 0.046 -0.024
5 -0.003 0.026 -0.008
6 0.037 0.049 -0.057
7 0.024 0.050 -0.045
8 0.023 0.055 -0.050
9 0.044 0.070 -0.062
per mouse, A feedback/ITI minus RSC feedback/ITI: {'10': 0.077, '3': 0.041, '4': 0.062, '5': 0.008, '6': 0.099, '7': 0.046, '8': 0.082, '9': 0.091}
Result: partially reproduced. Two parts of the claim hold and one does not. In the feedback period/ITI RSC is the least active area (-0.042 in reading b) and A the most active (0.023), A above RSC (p < 0.001, both readings, and in all 8 mice by 0.008 to 0.099); within RSC the feedback period/ITI is its lowest epoch (stem minus feedback/ITI 0.064, arm minus feedback/ITI 0.083, p < 0.001). But RSC is not shown to be most active in the stem: within RSC its arm mean is higher than its stem mean (0.042 against 0.023; stem minus arm -0.019, p = 0.07 in reading b, -0.017, p = 0.18 in reading a; arm higher in 7 of 8 mice), and between areas its stem mean is the highest and above V1, PM, MM and A (p <= 0.016 in reading b) but not above AM (p = 0.29; in reading a only PM and A differ). The arm here runs from the stem end to the choice point; the 1 s feedback delay has no maze position and falls in the feedback/ITI epoch (0 of 3,673,732 delay frames have a position). The well-fit reading of the neurons was not run.
Population decoders (fig3q-cue-decoding, figs5f-dynamic-choice-decoding, fig4d-partial-correlation, fig8cd-position-specific-decoders)¶
Methods, Trial phase-specific logistic decoders, helper decode_session: the simultaneously imaged neurons of a session are split into spatially adjacent, non-overlapping groups of about 100 (all layer 2/3 neurons, no fit selection); PCA to the fewest components that explain more than 90% of the variance; logistic regression with elastic net (10% L1, 90% L2) on single imaging frames (6 Hz volumes), one model per each of 11 equal position bins of the maze; 5-fold cross-validation gives a held-out prediction at every frame. Targets: cue, reported choice, and dynamic choice (the stored velocity-LSTM output, fitted as a binomial proportion). Frames are those of maze traversal on non-visually guided trials.
Readings, from the claims file: the penalty at the minimum CV deviance and the one-SE penalty; groups by k-means on neuron position or by a grid of equal-count tiles; for fig3q, a decoder's value over the stem bins or over all 11 bins. The Methods say only that lassoglm was run with 5-fold cross validation; its default partition is random observations (frames), which lets frames of one trial fall in both training and test sets, so folds of whole trials are run as a second reading (with k-means groups). The cross-position tests of fig8cd need held-out trials and use the trial folds. Areas are assigned by the area label of the neuron nearest each group's centroid. The penalty grid (10 values from the largest useful penalty down to 1e-4 of it) and the solver (FISTA, 150 iterations, checked against scikit-learn) are standard choices.
The paper has 974 decoders. Fitting them all under the three grouping and fold readings is beyond the notebook's 30-minute limit (the cell below prints the rate and the estimate for the whole pool), so the decoders are fit on a random subsample: 3 multi-plane layer 2/3 sessions of reading b per imaged field's area (seeded), 18 sessions. The single-plane layer 2/3 sessions of reading b (their number is printed below) are left out, as their frames are not volumes. The second notebook, decoders_v1_pm_offset.ipynb, fits every pool session whose field is V1 or PM (for fig3q) and the 6 pool sessions with the cue offset at 0.76 (for fig8cd); the third verification round agreed that its V1 and PM decoders decide fig3q, and its 0.76 sessions are a check of the fig8cd assumption.
rng_d = np.random.default_rng(110)
dpool = SEL["b"][(SEL["b"].layer == "L23") & (SEL["b"].imaging == "multi")]
dsess = pd.concat([x.iloc[np.sort(rng_d.permutation(len(x))[:3])] for _, x in dpool.groupby("fov")])
print(f"decoder sessions: {len(dsess)} of {len(dpool)} multi-plane layer 2/3 sessions in reading b; {dsess.mouse.nunique()} mice; "
f"{dsess.n_neurons.sum():,} neurons; sessions with the cue offset at 0.76 of the maze: {(dsess.cue_off.round(2) == 0.76).sum()} (in the pool {(dpool.cue_off.round(2) == 0.76).sum()}, by area {dpool[dpool.cue_off.round(2) == 0.76].fov.value_counts().to_dict()}); cue offsets in the decoder sessions {dsess.cue_off.round(2).value_counts().to_dict()}")
print("mice per field area:", dsess.groupby("fov").mouse.apply(lambda x: sorted(x, key=int)).to_dict())
PAIRS = [("kmeans", "frames"), ("tiles", "frames"), ("kmeans", "trials")]
t = time.time()
dec_res = tseng.decode_many(list(dsess.asset_id), PAIRS, workers=os.cpu_count())
fit_s = time.time() - t
print(f"decoders fit in {fit_s:.0f} s on {os.cpu_count()} processes: {fit_s * os.cpu_count() / len(dsess):.0f} CPU-s per session, "
f"{fit_s * os.cpu_count() / dsess.n_neurons.sum() * 1000:.1f} CPU-s per 1,000 neurons")
est = fit_s * os.cpu_count() / dsess.n_neurons.sum() * dpool.n_neurons.sum()
print(f"at this rate all {len(dpool)} pool sessions ({dpool.n_neurons.sum():,} neurons) would take {est / 3600:.1f} CPU hours, "
f"{est / os.cpu_count() / 60:.0f} min on {os.cpu_count()} processes")
print("single-plane layer 2/3 sessions of reading b, left out of the pool:",
int(((SEL["b"].layer == "L23") & (SEL["b"].imaging == "single")).sum()))
dfail = {a: r for a, r in dec_res.items() if isinstance(r, Exception)}
print("sessions that failed:", len(dfail), [repr(e)[:200] for e in dfail.values()])
mouse_of = dsess.set_index("asset_id").mouse
DEC = pd.DataFrame([r for v in dec_res.values() if not isinstance(v, Exception) for r in v if r["lam"] != "cross"])
DEC["mouse"] = DEC.asset_id.map(mouse_of)
XP = [dict(r, mouse=mouse_of[r["asset_id"]]) for v in dec_res.values() if not isinstance(v, Exception)
for r in v if r["lam"] == "cross"]
NB = tseng.DEC_NBINS
STEM_BINS = [b for b in range(NB) if (b + 0.5) / NB * D.maze_len.iloc[0] < D.stem.iloc[0]]
DEC["cue_stem"] = DEC[[f"cue_ll_{b}" for b in STEM_BINS]].mean(axis=1)
DEC["cue_all"] = DEC[[f"cue_ll_{b}" for b in range(NB)]].mean(axis=1)
for (g_, f_), x in DEC[DEC.lam == "min"].groupby(["grouping", "folds"]):
print(f"{g_} groups, {f_} folds: {len(x)} decoders; neurons per decoder {x.n_neurons.min()} to {x.n_neurons.max()} "
f"(median {x.n_neurons.median():.0f}); principal components {x.npc.min()} to {x.npc.max()} (median {x.npc.median():.0f}); "
f"by area {x.area.value_counts().reindex(AREAS).fillna(0).astype(int).to_dict()}")
k0 = DEC[(DEC.grouping == "kmeans") & (DEC.folds == "frames") & (DEC.lam == "min")]
print("paper (Figure 8D): V1 151, PM 106, AM 149, MM 108, RSC 233, A 227 decoders")
print(f"area of the neuron nearest the centroid equals the members' majority area in {(k0.area == k0.area_major).mean():.3f} of k-means decoders")
print(f"stem bins (centre in the stem): {STEM_BINS}; frames per decoder session: {DEC.groupby('asset_id').n_frames.first().min()} to "
f"{DEC.groupby('asset_id').n_frames.first().max()}; bin-target fits with too few frames or one class: "
f"{int(DEC[[f'cue_ll_{b}' for b in range(NB)]].isna().sum().sum())} cue bins")
decoder sessions: 18 of 142 multi-plane layer 2/3 sessions in reading b; 7 mice; 11,564 neurons; sessions with the cue offset at 0.76 of the maze: 0 (in the pool 6, by area {'RSC': 2, 'visA': 2, 'PM': 1, 'MM': 1}); cue offsets in the decoder sessions {1.0: 9, 0.96: 9}
mice per field area: {'AM': ['5', '7', '8'], 'MM': ['5', '7', '10'], 'PM': ['9', '10', '10'], 'RSC': ['3', '6', '8'], 'V1': ['6', '6', '10'], 'visA': ['5', '8', '9']}
decoders fit in 392 s on 4 processes: 87 CPU-s per session, 135.4 CPU-s per 1,000 neurons
at this rate all 142 pool sessions (95,988 neurons) would take 3.6 CPU hours, 54 min on 4 processes
single-plane layer 2/3 sessions of reading b, left out of the pool: 8
sessions that failed: 0 []
kmeans groups, frames folds: 116 decoders; neurons per decoder 49 to 146 (median 99); principal components 20 to 72 (median 45); by area {'V1': 19, 'PM': 12, 'AM': 16, 'MM': 18, 'RSC': 30, 'A': 21}
kmeans groups, trials folds: 116 decoders; neurons per decoder 49 to 146 (median 99); principal components 20 to 72 (median 45); by area {'V1': 19, 'PM': 12, 'AM': 16, 'MM': 18, 'RSC': 30, 'A': 21}
tiles groups, frames folds: 116 decoders; neurons per decoder 86 to 116 (median 100); principal components 25 to 59 (median 45); by area {'V1': 18, 'PM': 14, 'AM': 17, 'MM': 16, 'RSC': 30, 'A': 21}
paper (Figure 8D): V1 151, PM 106, AM 149, MM 108, RSC 233, A 227 decoders
area of the neuron nearest the centroid equals the members' majority area in 0.974 of k-means decoders
stem bins (centre in the stem): [0, 1, 2, 3, 4, 5, 6, 7, 8]; frames per decoder session: 8827 to 23283; bin-target fits with too few frames or one class: 0 cue bins
fig3q-cue-decoding¶
Cue decoded from populations of about 100 nearby neurons follows the same gradient as single-neuron cue encoding.
Cue decoding is the held-out log likelihood (base 2) of the cue per frame, averaged over a decoder's bins (stem bins or all bins), with chance at -1. Per area, hierarchical bootstrap over mice then decoders (1,000 samples); V1 against A as two unpaired bootstrap samples. Reproduced when, under every reading, V1 has the highest area mean and A the lowest, with V1 above A at p < 0.05.
rng = np.random.default_rng(111)
rows = []
for (g_, f_, lam_), sub in DEC.groupby(["grouping", "folds", "lam"]):
for bname, col in (("stem bins", "cue_stem"), ("all bins", "cue_all")):
B = {a: hboot([x[col].values for _, x in sub[sub.area == a].groupby("mouse")], 1000, rng) for a in AREAS if (sub.area == a).any()}
means = {a: B[a].mean() for a in B}
order = sorted(means, key=means.get, reverse=True)
rows.append(dict(groups=g_, folds=f_, penalty=lam_, bins=bname, **{a: round(means.get(a, np.nan), 3) for a in AREAS},
highest=order[0], lowest=order[-1], V1_minus_A=round(means["V1"] - means["A"], 3), p_V1_A=pu(B["V1"], B["A"])))
F3 = pd.DataFrame(rows)
print(F3.to_string(index=False))
print("chance (a constant prediction of 0.5): -1; readings with V1 highest and A lowest:",
int(((F3.highest == "V1") & (F3.lowest == "A")).sum()), "of", len(F3))
print("decoder values range:", round(DEC.cue_stem.min(), 3), "to", round(DEC.cue_stem.max(), 3), "(stem bins)")
print("per mouse, k-means / frames / minimum penalty, stem bins (area: mean, decoders):")
print(k0.groupby(["mouse", "area"]).cue_stem.agg(["mean", "size"]).round(3).unstack("area").to_string())
groups folds penalty bins V1 PM AM MM RSC A highest lowest V1_minus_A p_V1_A
kmeans frames 1se stem bins -0.184 -0.161 -0.295 -0.402 -0.456 -0.564 PM A 0.380 < 0.001
kmeans frames 1se all bins -0.228 -0.212 -0.337 -0.443 -0.495 -0.583 PM A 0.355 < 0.001
kmeans frames min stem bins -0.174 -0.153 -0.284 -0.390 -0.447 -0.558 PM A 0.384 < 0.001
kmeans frames min all bins -0.217 -0.206 -0.330 -0.433 -0.485 -0.575 PM A 0.358 < 0.001
kmeans trials 1se stem bins -0.211 -0.184 -0.334 -0.435 -0.493 -0.616 PM A 0.405 < 0.001
kmeans trials 1se all bins -0.265 -0.263 -0.382 -0.484 -0.536 -0.638 PM A 0.374 < 0.001
kmeans trials min stem bins -0.197 -0.179 -0.321 -0.421 -0.478 -0.606 PM A 0.409 < 0.001
kmeans trials min all bins -0.250 -0.248 -0.361 -0.469 -0.522 -0.625 PM A 0.375 < 0.001
tiles frames 1se stem bins -0.175 -0.168 -0.285 -0.397 -0.448 -0.567 PM A 0.392 < 0.001
tiles frames 1se all bins -0.206 -0.228 -0.328 -0.438 -0.488 -0.585 V1 A 0.379 < 0.001
tiles frames min stem bins -0.163 -0.162 -0.276 -0.388 -0.438 -0.554 PM A 0.390 < 0.001
tiles frames min all bins -0.200 -0.221 -0.320 -0.435 -0.477 -0.570 V1 A 0.370 < 0.001
chance (a constant prediction of 0.5): -1; readings with V1 highest and A lowest: 2 of 12
decoder values range: -0.942 to -0.035 (stem bins)
per mouse, k-means / frames / minimum penalty, stem bins (area: mean, decoders):
mean size
area A AM MM PM RSC V1 A AM MM PM RSC V1
mouse
10 NaN NaN -0.234 -0.132 NaN -0.136 NaN NaN 4.0 4.0 NaN 5.0
3 NaN NaN NaN NaN -0.474 NaN NaN NaN NaN NaN 12.0 NaN
5 -0.371 -0.187 -0.281 -0.036 -0.374 -0.044 8.0 4.0 3.0 1.0 3.0 1.0
6 NaN NaN NaN NaN -0.491 -0.233 NaN NaN NaN NaN 3.0 10.0
7 -0.459 -0.281 -0.382 -0.332 NaN -0.202 3.0 5.0 5.0 2.0 NaN 1.0
8 -0.748 -0.355 -0.552 -0.192 -0.378 NaN 6.0 7.0 4.0 2.0 10.0 NaN
9 -0.689 NaN -0.551 -0.074 -0.626 -0.128 4.0 NaN 2.0 3.0 2.0 2.0
Result: reproduced, with V1 and PM from the second notebook. On the 116 decoders of 18 sessions (a random subsample of the paper's 974), cue decoding follows a gradient from V1 and PM down through AM, MM and RSC to A under every reading: A is the lowest area in all 12 readings, and V1 is above A by 0.355 to 0.409 in log2 likelihood (p < 0.001). But V1 is the highest area in only 2 of the 12 readings (tile groups, all bins); in the other 10, PM is above V1 by 0.001 to 0.027 (PM rests on 12 to 14 decoders). On the subsample alone the criterion would hold in part. Folds of random frames and folds of whole trials give the same ranking. PM is above V1 by less than the spread its 12 to 14 decoders allow, so the subsample does not decide V1 against PM. decoders_v1_pm_offset.ipynb fits every decoder of the 41 pool sessions with a V1 or PM field (149 V1 and 71 PM decoders with k-means groups, 142 and 80 with tiles): there V1 is above PM under all 12 readings, by 0.027 to 0.062 (V1 -0.139 to -0.246, PM -0.187 to -0.282), and the highest of the other four areas here (AM, -0.276 to -0.382) is 0.118 to 0.157 below V1. With V1 and PM from the second notebook, V1 is the highest area and A the lowest under every reading, with V1 above A at p < 0.001, so the claim is reproduced; the third verification round agreed. Limits: V1 against PM is significant only for tiles over the stem bins (p = 0.039 and 0.042); the ranking takes AM, MM, RSC and A from this subsample (AM from 16 to 17 decoders in 3 mice) and V1 and PM from all their field sessions; decoders labelled PM in sessions of other fields are not fitted, so PM has 71 to 80 decoders against the paper's 106; the per-mouse and per-field checks of the second notebook are printed for one reading only (k-means groups, frame folds, minimum penalty, stem bins).
figs5f-dynamic-choice-decoding¶
From populations of about 100 neurons, dynamic choice is decoded better (Spearman 0.59 ± 0.02) than reported choice (0.26 ± 0.02), both above chance.
Spearman correlation of the held-out predictions with the true values over stem frames (position below the stem length), for the dynamic choice decoders and the reported choice decoders. Hierarchical bootstrap over mice then decoders (1,000 samples), each against zero (chance) and their paired difference. Tolerances are the paper's SEMs.
rng = np.random.default_rng(112)
for (g_, f_, lam_), sub in DEC.groupby(["grouping", "folds", "lam"]):
gg = [x for _, x in sub.groupby("mouse")]
bd = hboot([x.dyn_rho_stem.values for x in gg], 1000, rng)
br = hboot([x.choice_rho_stem.values for x in gg], 1000, rng)
bdf = hboot([(x.dyn_rho_stem - x.choice_rho_stem).values for x in gg], 1000, rng)
ok = (abs(bd.mean() - 0.59) <= 0.02, abs(br.mean() - 0.26) <= 0.02, abs(bdf.mean() - 0.33) <= 0.13)
print(f"{g_} groups, {f_} folds, {lam_} penalty, {len(sub)} decoders: dynamic {bd.mean():.3f} ± {bd.std():.3f} (p vs 0 {pb(bd)}); "
f"reported {br.mean():.3f} ± {br.std():.3f} (p vs 0 {pb(br)}); dynamic - reported {bdf.mean():.3f} ± {bdf.std():.3f} "
f"(p {pb(bdf)}); within tolerance (dynamic, reported, difference): {ok}")
print("per area, k-means / frames / minimum penalty:")
print(k0.groupby("area")[["dyn_rho_stem", "choice_rho_stem"]].agg(["mean", "size"]).reindex(AREAS).round(3).to_string())
print("per mouse:", k0.groupby("mouse")[["dyn_rho_stem", "choice_rho_stem"]].mean().round(3).to_dict("index"))
kmeans groups, frames folds, 1se penalty, 116 decoders: dynamic 0.669 ± 0.015 (p vs 0 < 0.001); reported 0.431 ± 0.020 (p vs 0 < 0.001); dynamic - reported 0.239 ± 0.011 (p < 0.001); within tolerance (dynamic, reported, difference): (np.False_, np.False_, np.True_)
kmeans groups, frames folds, min penalty, 116 decoders: dynamic 0.675 ± 0.014 (p vs 0 < 0.001); reported 0.438 ± 0.021 (p vs 0 < 0.001); dynamic - reported 0.238 ± 0.012 (p < 0.001); within tolerance (dynamic, reported, difference): (np.False_, np.False_, np.True_)
kmeans groups, trials folds, 1se penalty, 116 decoders: dynamic 0.636 ± 0.016 (p vs 0 < 0.001); reported 0.343 ± 0.025 (p vs 0 < 0.001); dynamic - reported 0.297 ± 0.016 (p < 0.001); within tolerance (dynamic, reported, difference): (np.False_, np.False_, np.True_)
kmeans groups, trials folds, min penalty, 116 decoders: dynamic 0.654 ± 0.015 (p vs 0 < 0.001); reported 0.373 ± 0.024 (p vs 0 < 0.001); dynamic - reported 0.279 ± 0.015 (p < 0.001); within tolerance (dynamic, reported, difference): (np.False_, np.False_, np.True_)
tiles groups, frames folds, 1se penalty, 116 decoders: dynamic 0.677 ± 0.014 (p vs 0 < 0.001); reported 0.435 ± 0.021 (p vs 0 < 0.001); dynamic - reported 0.241 ± 0.012 (p < 0.001); within tolerance (dynamic, reported, difference): (np.False_, np.False_, np.True_)
tiles groups, frames folds, min penalty, 116 decoders: dynamic 0.683 ± 0.014 (p vs 0 < 0.001); reported 0.444 ± 0.021 (p vs 0 < 0.001); dynamic - reported 0.240 ± 0.011 (p < 0.001); within tolerance (dynamic, reported, difference): (np.False_, np.False_, np.True_)
per area, k-means / frames / minimum penalty:
dyn_rho_stem choice_rho_stem
mean size mean size
area
V1 0.634 19 0.428 19
PM 0.693 12 0.455 12
AM 0.662 16 0.433 16
MM 0.665 18 0.419 18
RSC 0.689 30 0.427 30
A 0.706 21 0.471 21
per mouse: {'10': {'dyn_rho_stem': 0.663, 'choice_rho_stem': 0.468}, '3': {'dyn_rho_stem': 0.66, 'choice_rho_stem': 0.397}, '5': {'dyn_rho_stem': 0.709, 'choice_rho_stem': 0.484}, '6': {'dyn_rho_stem': 0.667, 'choice_rho_stem': 0.443}, '7': {'dyn_rho_stem': 0.659, 'choice_rho_stem': 0.364}, '8': {'dyn_rho_stem': 0.647, 'choice_rho_stem': 0.406}, '9': {'dyn_rho_stem': 0.747, 'choice_rho_stem': 0.529}}
Result: partially reproduced. Under every reading both decoders are above chance (p < 0.001) and dynamic choice is decoded better than reported choice (difference 0.238 to 0.297, p < 0.001, within the paper's 0.33 ± 0.13). But both values are above the paper's: dynamic choice 0.636 to 0.683 against 0.59 ± 0.02, reported choice 0.343 to 0.444 against 0.26 ± 0.02, under no reading within tolerance. Folds of whole trials give lower reported-choice values (0.343 to 0.373) than folds of frames (0.431 to 0.444). On a subsample of 116 decoders from 18 sessions.
fig4d-partial-correlation¶
Population decoders trained on reported choice are more related to dynamic choice (partial correlation 0.35 ± 0.02) than to reported choice (0.14 ± 0.02), a difference of 0.21 ± 0.02.
Per decoder, over stem frames: partial Spearman correlation (Pearson correlation of ranks, partialled) of the reported-choice decoder's held-out output with dynamic choice given reported choice, and with reported choice given dynamic choice. Hierarchical bootstrap over mice then decoders (1,000 samples). Tolerances are the paper's SEMs.
rng = np.random.default_rng(113)
for (g_, f_, lam_), sub in DEC.groupby(["grouping", "folds", "lam"]):
gg = [x for _, x in sub.groupby("mouse")]
bdy = hboot([x.pc_dyn.values for x in gg], 1000, rng)
brp = hboot([x.pc_rep.values for x in gg], 1000, rng)
bdf = hboot([(x.pc_dyn - x.pc_rep).values for x in gg], 1000, rng)
ok = (abs(bdy.mean() - 0.35) <= 0.02, abs(brp.mean() - 0.14) <= 0.02, abs(bdf.mean() - 0.21) <= 0.02)
print(f"{g_} groups, {f_} folds, {lam_} penalty, {len(sub)} decoders: with dynamic choice {bdy.mean():.3f} ± {bdy.std():.3f}; "
f"with reported choice {brp.mean():.3f} ± {brp.std():.3f}; difference {bdf.mean():.3f} ± {bdf.std():.3f} (p {pb(bdf)}); "
f"within tolerance (dynamic, reported, difference): {ok}; decoders with difference > 0: {(sub.pc_dyn > sub.pc_rep).mean():.3f}")
print("per area, k-means / frames / minimum penalty:")
print(k0.groupby("area")[["pc_dyn", "pc_rep"]].mean().reindex(AREAS).round(3).to_string())
print("per mouse, difference:", (k0.groupby("mouse").pc_dyn.mean() - k0.groupby("mouse").pc_rep.mean()).round(3).to_dict())
kmeans groups, frames folds, 1se penalty, 116 decoders: with dynamic choice 0.433 ± 0.015; with reported choice 0.204 ± 0.012; difference 0.229 ± 0.017 (p < 0.001); within tolerance (dynamic, reported, difference): (np.False_, np.False_, np.True_); decoders with difference > 0: 0.966
kmeans groups, frames folds, min penalty, 116 decoders: with dynamic choice 0.431 ± 0.014; with reported choice 0.215 ± 0.012; difference 0.217 ± 0.015 (p < 0.001); within tolerance (dynamic, reported, difference): (np.False_, np.False_, np.True_); decoders with difference > 0: 0.966
kmeans groups, trials folds, 1se penalty, 116 decoders: with dynamic choice 0.428 ± 0.017; with reported choice 0.101 ± 0.013; difference 0.326 ± 0.018 (p < 0.001); within tolerance (dynamic, reported, difference): (np.False_, np.False_, np.False_); decoders with difference > 0: 1.000
kmeans groups, trials folds, min penalty, 116 decoders: with dynamic choice 0.444 ± 0.016; with reported choice 0.130 ± 0.012; difference 0.314 ± 0.017 (p < 0.001); within tolerance (dynamic, reported, difference): (np.False_, np.True_, np.False_); decoders with difference > 0: 0.991
tiles groups, frames folds, 1se penalty, 116 decoders: with dynamic choice 0.439 ± 0.015; with reported choice 0.206 ± 0.014; difference 0.232 ± 0.018 (p < 0.001); within tolerance (dynamic, reported, difference): (np.False_, np.False_, np.False_); decoders with difference > 0: 0.966
tiles groups, frames folds, min penalty, 116 decoders: with dynamic choice 0.436 ± 0.014; with reported choice 0.217 ± 0.013; difference 0.219 ± 0.017 (p < 0.001); within tolerance (dynamic, reported, difference): (np.False_, np.False_, np.True_); decoders with difference > 0: 0.966
per area, k-means / frames / minimum penalty:
pc_dyn pc_rep
area
V1 0.388 0.221
PM 0.429 0.228
AM 0.408 0.218
MM 0.429 0.193
RSC 0.439 0.209
A 0.473 0.218
per mouse, difference: {'10': 0.158, '3': 0.209, '5': 0.184, '6': 0.255, '7': 0.222, '8': 0.226, '9': 0.267}
Result: partially reproduced. Under every reading the reported-choice decoders are more related to dynamic choice than to reported choice (difference 0.217 to 0.326, p < 0.001, positive in 97% to 100% of decoders). The partial correlation with dynamic choice is 0.428 to 0.444, above the paper's 0.35 ± 0.02 in every reading; that with reported choice is 0.101 to 0.217 (within 0.14 ± 0.02 only for trial folds at the minimum-CV penalty, 0.130); the difference is within 0.21 ± 0.02 for frame folds except tiles at the one-SE penalty (0.217 to 0.229). No reading has all three values within tolerance, so only the sign and test hold, which the criterion calls partially reproduced. On a subsample of 116 decoders from 18 sessions.
fig8cd-position-specific-decoders¶
Cue and choice decoders trained at one maze position lose accuracy when tested at other positions.
Each decoder (k-means groups, trial folds, minimum-CV penalty) trained in one position bin is applied to every bin of the held-out trials. Cue: bins with the cue present (centre before the session's cue offset, from the files' cue delay length), log likelihood (base 2). The Figure 8C legend uses only sessions with the cue offset at 0.76 of the maze (698 decoders, 100 in V1), but the files mark only 6 of the 142 sessions of the decoder pool with that offset (cue delay of 60 of 250 maze units: RSC 2, A 2, PM 1, MM 1, none with a field in V1 or AM; printed with the decoder sessions), and none of them is in the subsample, so every decoder session is used and the cue-present bins are set per session (a deviation); dynamic choice: all decoder sessions, all bins, Spearman correlation (Figure 8C, 8D legends). Per decoder, the mean over training bins of the performance at each train-test distance; the claim is tested as the mean at the largest distance minus that at zero distance, per area, hierarchical bootstrap over mice then decoders (1,000 samples), two-tailed. The 6 sessions with the cue offset at 0.76 are decoded in decoders_v1_pm_offset.ipynb.
rng = np.random.default_rng(114)
for v, label in (("cue", "cue, log2 likelihood"), ("dyn", "dynamic choice, Spearman")):
rows = []
for r in XP:
M = r[f"x_{v}"]
dist = np.abs(np.subtract.outer(np.arange(NB), np.arange(NB)))
prof = np.array([np.nanmean(M[(dist == d_) & ~np.isnan(M)]) if ((dist == d_) & ~np.isnan(M)).any() else np.nan
for d_ in range(NB)])
dmax = int(np.where(~np.isnan(prof))[0].max())
rows.append(dict(mouse=r["mouse"], area=r["area"], cue_off=r["cue_off"], d0=prof[0], dmax=prof[dmax], max_dist=dmax, **{f"dist{d_}": prof[d_] for d_ in range(NB)}))
R8 = pd.DataFrame(rows)
if v == "cue":
print("decoders by session cue offset:", R8.cue_off.round(2).value_counts().to_dict())
print(f"{label}: {len(R8)} decoders from {R8.mouse.nunique()} mice; largest distance {R8.max_dist.min()} to {R8.max_dist.max()} bins")
print(" mean performance by train-test distance (bins), all decoders:", R8[[f"dist{d_}" for d_ in range(NB)]].mean().round(3).dropna().to_dict())
for a in AREAS:
sub = R8[R8.area == a]
if len(sub) == 0:
print(f" {a}: no decoders")
continue
bdiff = hboot([(x.dmax - x.d0).values for _, x in sub.groupby("mouse")], 1000, rng)
print(f" {a}: {len(sub)} decoders, {sub.mouse.nunique()} mice; zero distance {sub.d0.mean():.3f}, largest distance {sub.dmax.mean():.3f}; "
f"largest - zero {bdiff.mean():.3f} ± {bdiff.std():.3f}, p {pb(bdiff)}; decoders with a drop {(sub.dmax < sub.d0).mean():.3f}")
decoders by session cue offset: {0.96: 72, 1.0: 44}
cue, log2 likelihood: 116 decoders from 7 mice; largest distance 10 to 10 bins
mean performance by train-test distance (bins), all decoders: {'dist0': -0.438, 'dist1': -0.536, 'dist2': -0.704, 'dist3': -0.823, 'dist4': -0.929, 'dist5': -1.047, 'dist6': -1.162, 'dist7': -1.315, 'dist8': -1.468, 'dist9': -1.578, 'dist10': -1.663}
V1: 19 decoders, 5 mice; zero distance -0.257, largest distance -2.061; largest - zero -1.884 ± 0.308, p < 0.001; decoders with a drop 1.000 PM: 12 decoders, 5 mice; zero distance -0.248, largest distance -2.105; largest - zero -1.858 ± 0.190, p < 0.001; decoders with a drop 1.000 AM: 16 decoders, 3 mice; zero distance -0.370, largest distance -1.721; largest - zero -1.373 ± 0.246, p < 0.001; decoders with a drop 1.000 MM: 18 decoders, 5 mice; zero distance -0.467, largest distance -1.532; largest - zero -1.053 ± 0.194, p < 0.001; decoders with a drop 1.000
RSC: 30 decoders, 5 mice; zero distance -0.519, largest distance -1.500; largest - zero -0.995 ± 0.157, p < 0.001; decoders with a drop 1.000
A: 21 decoders, 4 mice; zero distance -0.620, largest distance -1.352; largest - zero -0.710 ± 0.201, p < 0.001; decoders with a drop 1.000
dynamic choice, Spearman: 116 decoders from 7 mice; largest distance 10 to 10 bins
mean performance by train-test distance (bins), all decoders: {'dist0': 0.599, 'dist1': 0.553, 'dist2': 0.48, 'dist3': 0.426, 'dist4': 0.384, 'dist5': 0.342, 'dist6': 0.304, 'dist7': 0.253, 'dist8': 0.204, 'dist9': 0.161, 'dist10': 0.146}
V1: 19 decoders, 5 mice; zero distance 0.542, largest distance 0.074; largest - zero -0.471 ± 0.028, p < 0.001; decoders with a drop 1.000 PM: 12 decoders, 5 mice; zero distance 0.620, largest distance 0.138; largest - zero -0.474 ± 0.070, p < 0.001; decoders with a drop 1.000 AM: 16 decoders, 3 mice; zero distance 0.575, largest distance 0.158; largest - zero -0.415 ± 0.031, p < 0.001; decoders with a drop 1.000
MM: 18 decoders, 5 mice; zero distance 0.600, largest distance 0.143; largest - zero -0.462 ± 0.042, p < 0.001; decoders with a drop 1.000 RSC: 30 decoders, 5 mice; zero distance 0.625, largest distance 0.152; largest - zero -0.483 ± 0.050, p < 0.001; decoders with a drop 1.000
A: 21 decoders, 4 mice; zero distance 0.616, largest distance 0.199; largest - zero -0.423 ± 0.046, p < 0.001; decoders with a drop 1.000
Result: reproduced, resting on an assumption (cue decoders from all decoder sessions, with cue-present bins set per session). Decoding falls with the distance between training and test positions for both variables in every area: cue log2 likelihood from -0.438 at the training position to -1.663 ten bins away (largest minus zero distance -0.71 to -1.88 by area, p < 0.001 in each), dynamic choice Spearman correlation from 0.599 to 0.146 (-0.415 to -0.483 by area, p < 0.001 in each); every decoder drops. None of the 18 decoder sessions has the cue offset at 0.76 of the maze that the Figure 8C legend uses (their offsets are 0.96 and 1.0). On a subsample of 116 decoders from 18 sessions. The 6 pool sessions with the cue offset at 0.76 are decoded in decoders_v1_pm_offset.ipynb, under the claims file's own reading of the cue test.
Claims that need the GLM or the shattering analysis (blocked, too large)¶
The GLM claims (fig3i-glm-fit, well-fit-neurons, glm-cv-vs-test, fig3m-movement-task-gradient, fig3o-cue-gradient, figs4l-strategy-distributed, figs4p-previous-outcome, fig4b-dynamic-choice-stem, fig4c-reported-choice-feedback, figs5i-dynamic-choice-over-roll, figs5m-movement-history-control, fig4g-position-rsc, fig5e-cue-movement-distributed, fig5e-position-choice-strategy-distributed, fig6b-location-from-profiles, fig6d-three-gradients, figs6g-factor-profiles, encoding-dimensionality, fig7def-encoding-correlations, fig7g-area-decoding-strength-vs-correlation, figs8d-l2-vs-group-lasso, figs8k-robust-to-threshold, fig6h-embedding-overlap) rest on the Poisson GLM of the Methods fit to every layer 2/3 neuron. The dandiset holds no fitted coefficients, deviances or encoding strengths, so the fits would have to be repeated: per neuron, 5 cross-validation folds over 21 penalty values plus a final fit (106 fits) of a design with 972 task, 336 movement and the drift predictors on the session's frames. figs5l-dynamic-choice-over-cue-biased-running depends on well-fit-neurons and is blocked with it.
fig8g-shattering-dimensionality needs linear SVMs on every marginally balanced dichotomy of each conjunction, for 7 populations (6 areas and all areas), 10 train-test splits and 1,000 bootstrap pseudo-populations of 1,000 neurons, with two readings of choice (Methods, Quantification of shattering dimensionality; the claims file's criterion).
The timing below fits one model of each kind at the size the data give, on random data drawn for the purpose; it measures cost only, and no result is taken from it. The GLM fit uses a plain L2 penalty, faster than the group lasso of the Methods, so its estimates are lower bounds. The shattering estimate counts only the cue-by-position conjunction, so it is a lower bound too.
from sklearn.linear_model import PoissonRegressor
from sklearn.svm import SVC
from itertools import combinations
rng_t = np.random.default_rng(109)
l23 = SEL["b"][SEL["b"].layer == "L23"]
frames = int(np.median([res[a]["n_frames"] for a in l23.asset_id]) / 5) # frames per plane at 6 Hz
n_train = int(0.8 * 0.8 * frames)
Xg = rng_t.standard_normal((n_train, 1310)).astype(np.float64)
yg = rng_t.poisson(np.exp(np.log(0.2) + Xg @ rng_t.normal(0, 0.03, 1310)))
t = time.time()
PoissonRegressor(alpha=1e-3, max_iter=300).fit(Xg, yg)
t_glm = time.time() - t
n_neur = int(l23.n_neurons.sum())
ncpu = os.cpu_count()
print(f"frames per plane (median session) {frames}; one GLM fit of {n_train} frames x 1310 predictors: {t_glm:.1f} s")
for n_, what in [(n_neur, "all layer 2/3 neurons of reading b"), (2000, "a random subsample of 2,000 neurons"),
(200, "a random subsample of 200 neurons")]:
est = t_glm * 106 * n_ / 3600
print(f" {what} ({n_:,}) x 106 fits: {est:,.1f} CPU hours, {est / ncpu:,.1f} h on this machine's {ncpu} cores (lower bound)")
del Xg
# Shattering dimensionality: marginally balanced dichotomies of cue (2) x position (10): each class holds one
# condition per position and 5 of each cue, so the count is C(10, 5) / 2 (a dichotomy and its mirror are one).
n_dich = len(list(combinations(range(10), 5))) // 2
Xs = rng_t.standard_normal((20, 1000)); ys = np.r_[np.zeros(10), np.ones(10)]
t = time.time()
for _ in range(50):
SVC(kernel="linear", C=1.0).fit(Xs, ys)
t_svc = (time.time() - t) / 50
fits = n_dich * 10 * 1000 * 7 * 2
print(f"shattering, cue x position only: {n_dich} dichotomies x 10 splits x 1,000 pseudo-populations x 7 populations x 2 readings "
f"= {fits:,} SVC fits at {t_svc * 1e3:.2f} ms: {fits * t_svc / 3600:,.1f} CPU hours (lower bound: one conjunction of the "
f"several in Figure 8G, without building the trial samples)")
frames per plane (median session) 35701; one GLM fit of 22848 frames x 1310 predictors: 0.5 s all layer 2/3 neurons of reading b (97,989) x 106 fits: 1,354.7 CPU hours, 338.7 h on this machine's 4 cores (lower bound) a random subsample of 2,000 neurons (2,000) x 106 fits: 27.7 CPU hours, 6.9 h on this machine's 4 cores (lower bound) a random subsample of 200 neurons (200) x 106 fits: 2.8 CPU hours, 0.7 h on this machine's 4 cores (lower bound) shattering, cue x position only: 126 dichotomies x 10 splits x 1,000 pseudo-populations x 7 populations x 2 readings = 17,640,000 SVC fits at 0.86 ms: 4.2 CPU hours (lower bound: one conjunction of the several in Figure 8G, without building the trial samples)
Summary¶
| Claim | Outcome |
|---|---|
| neuron-counts | reproduced (assumed: session fraction correct over all trials) |
| fig1d-switch-performance | reproduced |
| lstm-choice-accuracy | reproduced (assumed: stored predictions are held out) |
| figs1h-lstm-beats-rl | reproduced |
| fig1i-strategy-after-switch | reproduced |
| visually-guided-rule-update | not reproduced (assumed: update of trial t is the change to t + 1) |
| fig1l-photoinhibition-performance | blocked, data missing |
| fig1m-rule-following-decrease | blocked, data missing |
| figs1m-rule-bias-uncorrelated | blocked, data missing |
| fig2c-running-variability-after-switch | reproduced |
| fig2f-latency-strategy | reproduced (assumed: session selection) |
| fig2g-choice-commitment | reproduced |
| cue-biased-running-uncorrelated | partially reproduced |
| fig4h-activity-by-epoch | partially reproduced |
| fig3q-cue-decoding | reproduced (V1 and PM from every V1 and PM field decoder in the second notebook; other areas from the subsample) |
| figs5f-dynamic-choice-decoding | partially reproduced (subsample of decoders) |
| fig4d-partial-correlation | partially reproduced (subsample of decoders) |
| fig8cd-position-specific-decoders | reproduced (assumed: cue decoders from all sessions; subsample of decoders; the 6 sessions with the 0.76 cue offset are in the second notebook) |
| 23 claims on the GLM fits (fig3i-glm-fit to fig6h-embedding-overlap) | blocked, too large |
| figs5l-dynamic-choice-over-cue-biased-running | blocked, too large (depends on well-fit-neurons) |
| fig8g-shattering-dimensionality | blocked, too large |
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 behavioural LSTM and the velocity LSTM were not refit; their outputs stored in the files were used. RL model 4 was refit with L-BFGS-B in Python in place of Matlab's fminunc. Sessions were selected on fraction correct over all trials, as an assumption (the Methods' rule of excluding visually guided trials gives counts far from the paper's); every analysis that depends on it also shows the other selection. Activity by epoch uses all neurons of the selected sessions, as the well-fit neurons need the GLM. The population decoders were fit on a random subsample of 18 of the 142 multi-plane layer 2/3 sessions (116 decoders against the paper's 974), with an elastic-net solver written here in place of Matlab's lassoglm, under folds of random frames and of whole trials; the cross-position cue test uses sessions whose cue offset is not the 0.76 of the Figure 8C legend. The second notebook, decoders_v1_pm_offset.ipynb, fits the decoders of all 41 pool sessions with a field in V1 or PM and of the 6 pool sessions with the 0.76 cue offset; its V1 and PM decoders decide fig3q (agreed in the third verification round), and its 0.76 sessions are a check of the fig8cd assumption. Nothing that needs the GLM or the shattering analysis was run; the timing above puts them beyond the run's limits.
What was found about the dataset.
- The stored
rule_beliefis P(rule A) - P(rule B) in all 300 sessions, the opposite sign to the Methods and to the dataset's tutorial, which both say positive means rule B. - The paper's session counts (253 sessions above 0.5, 216 and 141 layer 2/3 sessions above 0.65) cannot be recovered: with visually guided trials excluded the counts are 291, 154 and 100; over all trials 297, 238 and 150. Only the second matches the neuron counts within 10%. The files do not mark which sessions the paper used.
- The median session has 23% visually guided trials, and 222 of 300 sessions have more than 20%, where the paper says the fraction typically varied from 0 to 20%.
- In the single-plane sessions the deconvolved activity timestamps, labelled seconds, run 2000 times the behaviour clock; the planes have one frame per behaviour frame, so they are matched by index here.
- The files' cue delay marks only 6 of the 142 multi-plane layer 2/3 sessions of reading b with the cue offset at 0.76 of the maze, none with a field in V1 or AM, while the Figure 8C legend has 698 decoders (100 in V1) from such sessions.
- The photoinhibition sessions of Figures 1J to 1M and S1K to S1M are not in the dandiset.
- The dandiset holds no GLM fits, encoding strengths or population decoders.
- The total of 278,155 neurons in 300 sessions of 8 mice matches the paper exactly, and the area labels give layer 2/3 counts per area within 10% of the paper's.
print(f"notebook ran in {(time.time() - T0) / 60:.1f} min")
notebook ran in 16.2 min