Reproducing Pagan et al. (2025) from dandiset 001550¶
Paper. Pagan M, Tang VD, Aoi MC, Pillow JW, Mante V, Sussillo D, Brody CD (2025). Individual variability of neural computations underlying flexible decisions. Nature 639, 421-429. doi:10.1038/s41586-024-08433-6. This attempt used the open access full text (PubMed Central PMC11903320, CC BY), with the Methods and the Extended Data legends.
Dandiset. DANDI:001550, version 0.260916.0811. The paper cites no DANDI version: its data statement points to a GitHub repository, and the dandiset was made later from the same behavioural files. This is the only published version, so it is pinned; the draft has the same 16,115 assets and size.
What the dandiset holds. The behaviour of 26 rats: 16,115 sessions in 11 protocols, each with a trials table (context, outcome, generative evidence strengths, and the time of every pulse) and, in some TaskSwitch6 sessions, optogenetic epochs. The electrophysiology (3,495 units) is not in it. Files are streamed from the archive: each file is read whole into memory and opened there.
Outside inputs. The paper's data repository
Brody-Lab/flexible_decision_making_rats
at commit 0ea8fdb7c56d755d99155171cc314dbbc24393ac: its per-session behavioural tables (to test
the link to the dandiset), its trial counts and saved regression weights (as reference values),
and its neural data and fitted per-neuron pulse kernels (for the neural claims, which the dandiset
cannot test). Everything is fetched from GitHub at that commit, into memory.
| Claim | Data | Tested |
|---|---|---|
| session-link | with inputs | yes |
| subjects | dandiset | yes |
| trials-per-session | dandiset | yes |
| blocks-per-session | dandiset | yes |
| session-duration | dandiset | yes |
| block-min-30 | dandiset | yes |
| fig1c-20-rats | dandiset | yes |
| fig1c-context-selection | dandiset | yes |
| ed1e-switch-four-trials | dandiset | yes |
| ed2g-learning-2-5-months | dandiset | yes |
| fig4-converging | dandiset | yes |
| fig4-loc-frq-uncorrelated | dandiset | yes |
| ed3g-halfsplit | dandiset | yes |
| ed4-bilateral-fof-impairs | dandiset | yes |
| neural-link | with inputs | yes |
| units-per-rat | elsewhere (GitHub) | yes |
| fig3h-neural-uncorrelated | elsewhere (GitHub) | yes |
| choice-axes-parallel | elsewhere (GitHub) | yes |
| fig5d-behaviour-neural | with inputs | yes |
| fig1e-choice-axis | elsewhere (GitHub) | blocked: methods underspecified |
| fig1f-irrelevant-encoded | elsewhere (GitHub) | blocked: not runnable |
| fig3a-trained-rnn-svm | elsewhere (GitHub) | blocked: not runnable |
| fig5c-engineered-rnn | elsewhere (GitHub) | blocked: not runnable |
The claims and their criteria are in claims/001550/10.1038-s41586-024-08433-6.yaml; the method
choices common to all claims are in its opening comment and in the docstring of taskswitch.py.
import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import time
import warnings
from concurrent.futures import ThreadPoolExecutor
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import taskswitch as ts
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
assets = ts.asset_table()
task_assets = assets[assets.protocol.isin(ts.TASK_PROTOCOLS)].reset_index(drop=True)
results = ts.load_many(task_assets.url)
other_assets = assets[assets.protocol.str.startswith("TaskSwitch") & ~assets.protocol.isin(ts.TASK_PROTOCOLS)].reset_index(drop=True)
other_results = ts.load_many(other_assets.url, keep_bins=False)
print(f"dandiset {ts.DANDISET_ID} version {ts.DANDISET_VERSION}: {len(assets)} files; "
f"read {len(results)} task sessions and {len(other_results)} other TaskSwitch sessions "
f"in {time.time() - started:.0f} s")
def session_table(table, res):
s = pd.DataFrame([r["session"] for r in res])
return pd.concat([table[["subject", "protocol", "session", "date", "size"]].reset_index(drop=True), s], axis=1)
S = session_table(task_assets, results)
S_other = session_table(other_assets, other_results)
S["readable"] = S.missing.fillna("") == ""
S["selected"] = S.readable & (S.valid >= ts.MIN_VALID) & (S.hit >= ts.MIN_HIT) & (S.fsi >= ts.MIN_FSI)
dandiset 001550 version 0.260916.0811: 16115 files; read 15413 task sessions and 440 other TaskSwitch sessions in 66 s
Data check¶
Each check of the metadata against the values is stated under its output.
print("Files per protocol:")
print(assets.protocol.value_counts().to_string())
print("\nSubjects:", assets.subject.nunique(), "| files per subject, min to max:",
assets.subject.value_counts().min(), "to", assets.subject.value_counts().max())
print("Total size of the task-protocol files:", round(task_assets["size"].sum() / 1e9, 1), "GB")
print("\nTask sessions without one of the required trial columns:", int((~S.readable).sum()))
print(S.loc[~S.readable, ["subject", "session", "missing", "n_trials"]].to_string(index=False))
ok = S[S.readable]
print("\nSessions with a hit_history column:", int(ok.has_hit_history.sum()), "of", len(ok),
"| sessions where it disagrees with the outcome derived from result_history:",
int((ok.hit_history_mismatch > 0).sum()))
print("dtype of HistorySection_gammafreq_history:", ok.hist_gfreq_dtype.value_counts().to_dict(),
"| sessions where it differs from gamma_freq:", int((ok.hist_gfreq_mismatch > 0).sum()))
print("Sex:", ok.sex.value_counts().to_dict(), "| strain:", ok.strain.value_counts().to_dict())
Files per protocol:
protocol
TaskSwitch4 13777
TaskSwitch6 1550
TaskSwitch2 201
ProAnti3Marino 173
TaskSwitch3 111
TaskSwitch4double 96
TaskSwitch4new 86
PBups 73
TaskSwitch 19
ProAnti3 16
TaskSwitch4repeat 13
Subjects: 26 | files per subject, min to max: 433 to 827
Total size of the task-protocol files: 20.1 GB
Task sessions without one of the required trial columns: 13
subject session missing n_trials
H113 TaskSwitch4-190110a right_lo 2
H113 TaskSwitch4-190122b HistorySection_result_history,left_lo 1
P049 TaskSwitch4-180223a HistorySection_result_history,right_hi 1
P055 TaskSwitch4-160507a right_hi,left_hi 2
P055 TaskSwitch4-160820a HistorySection_result_history,left_hi 1
P055 TaskSwitch4-170205a HistorySection_result_history,right_hi 1
P055 TaskSwitch4-171026a HistorySection_result_history 1
P055 TaskSwitch4-180516a HistorySection_result_history,right_hi,right_lo 1
P059 TaskSwitch4-160827a HistorySection_result_history,right_lo,left_lo 1
P072 TaskSwitch4-180227c HistorySection_result_history,left_hi 1
P095 TaskSwitch4-171130a left_lo 2
P124 TaskSwitch6-190726a right_hi 3
P135 TaskSwitch4-190514a HistorySection_result_history,left_hi,right_lo,left_lo 1
Sessions with a hit_history column: 15305 of 15400 | sessions where it disagrees with the outcome derived from result_history: 0
dtype of HistorySection_gammafreq_history: {'float64': 14830, 'int16': 351, 'uint8': 219} | sessions where it differs from gamma_freq: 0
Sex: {'M': 15400} | strain: {'Long Evans': 15400}
Thirteen task sessions lack the result codes or a pulse column, and are left out of every count
below. Where both exist, hit_history agrees with the outcome derived from result_history on
every trial, so the outcome is taken from result_history, which every other session has. Some
files store HistorySection_gammafreq_history, a signed strength, as an unsigned or 16-bit
integer; in every one of them its values equal the float gamma_freq, so no sign was lost.
trials_all = pd.concat([r["trials"].assign(sid=i) for i, r in enumerate(results) if r["trials"] is not None],
ignore_index=True)
valid_all = trials_all[trials_all.hit.notna()]
print("Trials:", len(trials_all), "| valid:", len(valid_all),
"| result codes:", trials_all.result.value_counts().sort_index().to_dict())
print("Context labels:", trials_all.task.value_counts().to_dict(), "| side labels:", trials_all.side.value_counts().to_dict())
print("gamma_dir values (valid trials):", sorted(valid_all.gdir.unique().round(3).tolist()))
print("gamma_freq values (valid trials):", sorted(valid_all.gfreq.unique().round(3).tolist()))
print("duration (s) of the pulse train, valid trials:", valid_all.duration.value_counts().head(5).to_dict())
loose = S[S.readable & S.loose]
print("\nSessions whose pulse times were binned:", len(loose))
print("Pulse times over valid trials (s): min", round(loose.pulse_min.min(), 5), "max", round(loose.pulse_max.max(), 5),
"| pulses outside 0 to 1.3 s:", int(loose.pulses_outside.sum()))
print("Pulses per valid trial: median", round(loose.pulses_per_trial.median(), 2),
"range", round(loose.pulses_per_trial.min(), 2), "to", round(loose.pulses_per_trial.max(), 2),
f"(the Methods' 40 Hz over 1.3 s gives {40 * 1.3:.0f})")
Trials: 11545151 | valid: 7611092 | result codes: {1: 5757103, 2: 1853193, 3: 3367532, 4: 566527, 5: 790, 6: 6}
Context labels: {'d': 6195573, 'f': 5349578} | side labels: {'r': 5819368, 'l': 5725783}
gamma_dir values (valid trials): [-4.0, -3.5, -3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0]
gamma_freq values (valid trials): [-4.0, -3.5, -3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0]
duration (s) of the pulse train, valid trials: {1.3: 7609091, 10.0: 2001}
Sessions whose pulse times were binned: 12530
Pulse times over valid trials (s): min 0.0 max 1.3 | pulses outside 0 to 1.3 s: 0
Pulses per valid trial: median 52.19 range 48.9 to 55.15 (the Methods' 40 Hz over 1.3 s gives 52)
Pulse times lie within 0 to 1.3 s of the start of the pulse train, and there are about 52 per
valid trial, the 40 Hz over 1.3 s of the Methods. The generative strengths take six values,
plus 0 for the training stages that lack one feature. The pulse columns' description says
seconds from the centre poke; the optogenetic epochs below start at cpoke_start_time or 0.65 s
after it, which is consistent with the pulse train starting at the centre poke.
# Do the pulse columns, the correct side and the choice describe the same trials? Contrast the
# true pairing with each trial's pulses paired with the previous trial's side and choice.
rng = np.random.default_rng(1)
rows = []
for i in rng.choice(np.flatnonzero(S.selected), 200, replace=False):
r = results[i]
tr = r["trials"][r["trials"].hit.notna()].reset_index(drop=True)
p = r["pulses"].sum(2).astype(int) # right_hi, left_hi, right_lo, left_lo
loc = (p[:, 0] + p[:, 2]) - (p[:, 1] + p[:, 3])
frq = (p[:, 0] + p[:, 1]) - (p[:, 2] + p[:, 3])
rel = np.where(tr.task == "d", loc, frq)
right = (tr.side == "r").values
keep = rel != 0
rows.append(dict(side_true=np.mean((rel[keep] > 0) == right[keep]),
side_shift=np.mean((rel[1:][keep[1:]] > 0) == right[:-1][keep[1:]]),
choice_true=np.mean((rel[keep] > 0) == (tr.choice.values[keep] == 1)),
choice_shift=np.mean((rel[1:][keep[1:]] > 0) == (tr.choice.values[:-1][keep[1:]] == 1)),
gdir_sign=np.mean(np.sign(loc[tr.gdir != 0]) == np.sign(tr.gdir[tr.gdir != 0]))))
chk = pd.DataFrame(rows)
print("200 selected sessions; fraction of valid trials (relevant pulse difference not 0) where:")
print(chk.describe().loc[["mean", "min", "max"]].round(3).to_string())
# Optogenetic epochs against the trial clock.
eps = []
for i, r in enumerate(results):
if r["opto"] is not None and r["trials"] is not None:
eps.append(ts.opto_trials(r["trials"], r["opto"]).assign(subject=S.subject[i], sid=i))
eps = pd.concat(eps, ignore_index=True)
print("\nOptogenetic epochs:", len(eps), "in", eps.sid.nunique(), "sessions of", sorted(eps.subject.unique()))
print("Epochs whose start falls inside a trial:", int(eps.in_trial.sum()))
print("Start minus cpoke_start_time (s), rounded:", eps.offset.round(3).value_counts().head(6).to_dict())
print("Epoch length (s):", eps.duration.value_counts().to_dict())
print("Sites:", eps.sites.value_counts().to_dict(), "| power (mW):", eps.power.value_counts().to_dict(),
"| wavelength (nm):", eps.wavelength.value_counts().to_dict())
200 selected sessions; fraction of valid trials (relevant pulse difference not 0) where:
side_true side_shift choice_true choice_shift gdir_sign
mean 1.000 0.504 0.786 0.493 1.000
min 0.961 0.441 0.701 0.427 0.977
max 1.000 0.570 0.905 0.549 1.000
Optogenetic epochs: 42596 in 237 sessions of ['P116', 'P124', 'P127', 'P131']
Epochs whose start falls inside a trial: 42596
Start minus cpoke_start_time (s), rounded: {0.0: 31540, 0.65: 11056}
Epoch length (s): {0.65: 22076, 1.3: 20520}
Sites: {'bilateral': 14298, 'left': 14229, 'right': 14069} | power (mW): {25.0: 42596} | wavelength (nm): {473.0: 42596}
In selected sessions, the side with more relevant pulses is the stored correct side on every
trial where the counts differ, and the choice follows the relevant pulses on most trials; paired
with the previous trial instead, both fall to about one half. So the pulse columns, the correct
side and the derived choice refer to the same trials, and the sign of the location pulse
difference agrees with gamma_dir. Every optogenetic epoch starts inside a trial, at the centre
poke or 0.65 s after it, and lasts 1.3 or 0.65 s, as the column description of opto_type says
(full trial, first or second half). The files give the wavelength as 473 nm; Extended Data Figure
4 says 450 nm. Epochs are in four rats, where the paper reports three.
sel = S[S.selected]
print("Sessions that pass the Methods selection (>= 100 valid trials, >= 0.7 correct, feature selection index >= 0.7):",
len(sel), "of", int(S.readable.sum()))
print(sel.protocol.value_counts().to_dict())
per_rat = S[S.readable].groupby("subject").agg(sessions=("session", "size"), selected=("selected", "sum"),
valid=("valid", "sum"))
per_rat["valid_selected"] = sel.groupby("subject").valid.sum()
print(per_rat.astype(int).to_string())
Sessions that pass the Methods selection (>= 100 valid trials, >= 0.7 correct, feature selection index >= 0.7): 8849 of 15400
{'TaskSwitch4': 7917, 'TaskSwitch6': 858, 'TaskSwitch4new': 74}
sessions selected valid valid_selected
subject
H113 624 466 249273 203339
P007 363 214 218781 148259
P013 461 242 387436 233734
P049 747 62 259022 23014
P055 814 706 512241 480165
P059 639 69 268802 25847
P072 580 233 323287 177844
P074 668 234 340403 149761
P078 500 264 261584 147043
P095 607 366 316247 197061
P100 719 553 355776 283348
P101 579 260 218400 100290
P102 629 465 372177 290572
P116 621 213 272067 107940
P120 547 465 373323 319079
P123 588 401 283478 202690
P124 657 365 197444 115123
P127 605 430 268024 195545
P130 618 521 335692 287944
P131 677 227 259382 104187
P135 591 367 365305 251475
P150 616 412 402184 292255
P159 427 288 186698 135369
P167 488 368 171723 128935
P187 488 378 194224 153441
P189 547 280 218119 128309
session-link¶
The dandiset's sessions can be linked to the paper's published per-session behavioural tables, and the linked records agree on trial counts, outcomes, evidence and pulse timing.
rng = np.random.default_rng(101)
rats = sorted(assets.subject.unique())
blobs = {}
with ThreadPoolExecutor(8) as ex:
for rat, b in zip(rats, ex.map(lambda r: ts.paper_file(f"data/Sessions_trials/{r}_sessions.mat"), rats)):
blobs[rat] = b
from scipy.io import loadmat
import io as _io
paper_rows = []
for rat in rats:
m = loadmat(_io.BytesIO(blobs[rat]), squeeze_me=True, struct_as_record=False, variable_names=["sessions"])
for k, s in enumerate(np.atleast_1d(m["sessions"])):
paper_rows.append(dict(subject=rat, k=k, file=s.file, session=ts.paper_session_name(s.file),
valid_p=s.valid, hit_p=s.hit, good=s.good, bdir=s.bdir, bfreq=s.bfreq))
Pp = pd.DataFrame(paper_rows)
Pp["has_valid"] = Pp.valid_p.notna() & (Pp.valid_p > 0)
Pp["selected_p"] = (Pp.hit_p >= 0.7) & ((Pp.bdir + Pp.bfreq) / 2 >= 0.7) & (Pp.valid_p >= 100) & (Pp.good == 1)
M = Pp.merge(S[S.readable], on=["subject", "session"], how="outer", indicator=True)
print("Paper sessions:", len(Pp), "| with valid trials:", int(Pp.has_valid.sum()),
"| found in the dandiset (readable):", int((M._merge == "both").sum()))
found = M[(M._merge != "right_only") & M.has_valid.fillna(False)]
frac_found = (found._merge == "both").mean()
print(f"Paper sessions with valid trials found in the dandiset: {frac_found:.4f}")
print("Paper sessions not found:", found.loc[found._merge == "left_only", ["subject", "session", "valid_p"]].to_string(index=False))
both = M[(M._merge == "both") & M.has_valid.fillna(False)]
eq_valid = (both.valid_p == both.valid).mean()
eq_hit = np.isclose(both.hit_p, both.hit).mean()
rel_ok = np.isclose(both.bdir, both.rel_loc, atol=1e-6) & np.isclose(both.bfreq, both.rel_frq, atol=1e-6)
print(f"Matched sessions: {len(both)}; valid trials equal in {eq_valid:.4f}, fraction correct equal in {eq_hit:.4f}, "
f"both relative weights within 1e-6 in {rel_ok.mean():.4f}")
mis = both[~rel_ok]
capped_loc, capped_frq = np.minimum(mis.rel_loc, 1), np.minimum(mis.rel_frq, 1)
ok_capped = np.isclose(mis.bdir, capped_loc, atol=1e-6) & np.isclose(mis.bfreq, capped_frq, atol=1e-6)
print(f"\nSessions whose relative weights differ: {len(mis)}; of them equal once this attempt's values are capped at 1: "
f"{int(ok_capped.sum())}; differing sessions that pass the paper's selection: {int(mis.selected_p.sum())}")
rest = mis[~ok_capped]
print("Differing after capping:", len(rest), "| with a missing value on either side:",
int((rest[["bdir", "bfreq", "rel_loc", "rel_frq"]].isna().any(axis=1)).sum()),
"| valid trials in them: median", rest.valid.median(), "max", rest.valid.max(),
"| of them in the paper's selection:", int(rest.selected_p.sum()))
na_p = rest[["bdir", "bfreq"]].isna().any(axis=1)
na_d = rest[["rel_loc", "rel_frq"]].isna().any(axis=1)
print("Of them, missing in the paper's table only:", int((na_p & ~na_d).sum()), "| in this attempt's values only:",
int((~na_p & na_d).sum()), "| on both sides:", int((na_p & na_d).sum()), "| on neither side:", int((~na_p & ~na_d).sum()))
print(rest[["subject", "session", "valid", "hit", "bdir", "rel_loc", "bfreq", "rel_frq"]].head(12).round(4).to_string(index=False))
both_present = rest[~na_p & ~na_d]
d = np.maximum(np.abs(both_present.bdir - np.minimum(both_present.rel_loc, 1)), np.abs(both_present.bfreq - np.minimum(both_present.rel_frq, 1)))
print(f"\nDiffering with both values present: {len(both_present)} | largest difference after capping: median {d.median():.4f}, "
f"max {d.max():.4f} | valid trials: median {both_present.valid.median():.0f}, range {both_present.valid.min()} to {both_present.valid.max()} "
f"| in the paper's selection: {int(both_present.selected_p.sum())}")
print("Their rats:", both_present.subject.value_counts().to_dict())
print(both_present.assign(diff=d)[["subject", "session", "valid", "hit", "bdir", "rel_loc", "bfreq", "rel_frq", "diff"]]
.sort_values("diff", ascending=False).head(12).round(4).to_string(index=False))
only_d = M[M._merge == "right_only"]
print("\nTask sessions in the dandiset and not in the paper's tables:", len(only_d), "| of them pass the selection:",
int(only_d.selected.sum()))
print(only_d[["subject", "session", "valid", "hit", "fsi", "selected"]].round(3).to_string(index=False))
extra = M[(M._merge == "both") & M.selected.astype(bool) & ~M.selected_p.fillna(False).astype(bool)]
print("\nSelected here, in the paper's tables, not selected there:")
print(extra[["subject", "session", "valid", "hit", "fsi", "valid_p", "hit_p", "good", "bdir", "bfreq", "rel_loc", "rel_frq"]].round(4).to_string(index=False))
xs = pd.crosstab(M.selected.fillna(False).rename("dandiset selection"), M.selected_p.fillna(False).rename("paper selection with good flag"))
print("\n", xs.to_string())
Paper sessions: 15487 | with valid trials: 15278 | found in the dandiset (readable): 15394
Paper sessions with valid trials found in the dandiset: 0.9999
Paper sessions not found: subject session valid_p
P189 TaskSwitch6-210516a 27.0
Matched sessions: 15277; valid trials equal in 1.0000, fraction correct equal in 1.0000, both relative weights within 1e-6 in 0.8983
Sessions whose relative weights differ: 1553; of them equal once this attempt's values are capped at 1: 861; differing sessions that pass the paper's selection: 650
Differing after capping: 692 | with a missing value on either side: 632 | valid trials in them: median 201.0 max 981.0 | of them in the paper's selection: 3
Of them, missing in the paper's table only: 5 | in this attempt's values only: 0 | on both sides: 627 | on neither side: 60
subject session valid hit bdir rel_loc bfreq rel_frq
H113 TaskSwitch4-170518a 473.0 0.7928 NaN NaN NaN NaN
H113 TaskSwitch4-170519a 609.0 0.5501 NaN NaN NaN NaN
H113 TaskSwitch4-170520a 475.0 0.7221 NaN NaN NaN NaN
H113 TaskSwitch4-170521a 645.0 0.8295 NaN NaN NaN NaN
H113 TaskSwitch4-170522a 644.0 0.8556 NaN NaN NaN NaN
H113 TaskSwitch4-170523a 662.0 0.8625 NaN NaN NaN NaN
H113 TaskSwitch4-170525a 914.0 0.8775 NaN NaN NaN NaN
H113 TaskSwitch4-170526a 613.0 0.8874 NaN NaN NaN NaN
H113 TaskSwitch4-170703a 30.0 0.7667 NaN 0.6875 NaN NaN
H113 TaskSwitch4-170716a 16.0 0.7500 NaN 568.2172 NaN NaN
H113 TaskSwitch4-171120a 7.0 0.1429 NaN NaN NaN NaN
H113 TaskSwitch4-180313a 244.0 0.8074 0.6786 0.6786 0.8283 0.8301
Differing with both values present: 60 | largest difference after capping: median 0.8732, max 13.7986 | valid trials: median 198, range 39.0 to 916.0 | in the paper's selection: 3
Their rats: {'P189': 8, 'P127': 6, 'P131': 5, 'P013': 4, 'P123': 4, 'P135': 4, 'P150': 4, 'P167': 4, 'P101': 3, 'P159': 3, 'H113': 2, 'P100': 2, 'P187': 2, 'P049': 1, 'P055': 1, 'P059': 1, 'P072': 1, 'P078': 1, 'P095': 1, 'P102': 1, 'P120': 1, 'P130': 1}
subject session valid hit bdir rel_loc bfreq rel_frq diff
P131 TaskSwitch4-180106a 180.0 0.5111 1.0000 -12.7986 0.2313 0.2313 13.7986
P100 TaskSwitch4-170610a 594.0 0.5101 0.2715 0.2715 1.0000 -11.9991 12.9991
P095 TaskSwitch4-180713a 149.0 0.3826 1.0000 3.6079 0.0000 -7.7570 7.7570
P167 TaskSwitch4-181212a 361.0 0.5291 1.0000 -3.6637 1.0000 2.1512 4.6637
P078 TaskSwitch4-180313a 94.0 0.4894 1.0000 -3.4661 0.0000 21.0193 4.4661
P189 TaskSwitch6-200819a 176.0 0.5568 1.0000 -2.6819 0.3961 0.3961 3.6819
P123 TaskSwitch4-190105a 425.0 0.3412 0.0000 -3.2861 1.0000 1.0185 3.2861
P127 TaskSwitch4-180701a 47.0 0.5745 1.0000 1.1313 1.0000 -1.1439 2.1439
P189 TaskSwitch6-210105a 39.0 0.6923 1.0000 1.2727 1.0000 -1.0416 2.0416
P131 TaskSwitch4-190411a 112.0 0.5536 1.0000 1.3711 0.0000 -1.7742 1.7742
P189 TaskSwitch6-191003a 196.0 0.5459 1.0000 -0.5823 0.6196 0.6196 1.5823
P135 TaskSwitch4-171028a 146.0 0.4247 1.0000 2.7815 0.0000 -1.5536 1.5536
Task sessions in the dandiset and not in the paper's tables: 6 | of them pass the selection: 5
subject session valid hit fsi selected
P187 TaskSwitch6-210117a 247.0 0.826 0.828 True
P187 TaskSwitch6-210121a 333.0 0.829 0.830 True
P187 TaskSwitch6-210122a 404.0 0.809 0.755 True
P187 TaskSwitch6-210224a 458.0 0.707 0.767 True
P189 TaskSwitch6-210204a 25.0 0.800 NaN False
P189 TaskSwitch6-210226a 609.0 0.818 0.771 True
Selected here, in the paper's tables, not selected there:
subject session valid hit fsi valid_p hit_p good bdir bfreq rel_loc rel_frq
P059 TaskSwitch4-170331a 213.0 0.7089 0.7527 213.0 0.7089 1.0 1.0 0.3935 1.1119 0.3935
paper selection with good flag False True
dandiset selection
False 6644 0
True 6 8843
# Per-trial comparison on a seeded sample of 20 selected sessions found in both.
cand = M[(M._merge == "both") & M.selected & M.selected_p.fillna(False)]
sample = cand.iloc[rng.choice(len(cand), 20, replace=False)]
url_of = dict(zip(task_assets.subject + "/" + task_assets.session, task_assets.url))
rows = []
for _, row in sample.iterrows():
tr_p = loadmat(_io.BytesIO(blobs[row.subject]), squeeze_me=True, struct_as_record=False,
variable_names=["trials"])["trials"][int(row.k)]
d = ts.load_session(url_of[f"{row.subject}/{row.session}"], bin_s=0.01, n_bins=130, all_bins=True)
tr = d["trials"][d["trials"].hit.notna()].reset_index(drop=True)
p = d["pulses"].astype(int)
mine = dict(vecr=p[:, 0] + p[:, 2], vecl=p[:, 1] + p[:, 3], vechi=p[:, 0] + p[:, 1], veclo=p[:, 2] + p[:, 3])
n = len(tr)
same_n = n == len(np.atleast_1d(tr_p.choice_all))
agree = (same_n
and np.array_equal(list(str(tr_p.task_all)), tr.task.tolist())
and np.array_equal(list(str(tr_p.side_all)), tr.side.tolist()))
per_trial = np.ones(n, bool) if same_n else np.zeros(n, bool)
if same_n:
per_trial &= np.array(list(str(tr_p.task_all))) == tr.task.values
per_trial &= np.array(list(str(tr_p.side_all))) == tr.side.values
per_trial &= np.isclose(tr_p.gdir_all, tr.gdir) & np.isclose(tr_p.gfreq_all, tr.gfreq)
per_trial &= np.atleast_1d(tr_p.choice_all) == tr.choice.values
per_trial &= np.atleast_1d(tr_p.hits_all) == tr.hit.values
nta_eq = np.mean(np.atleast_1d(tr_p.nta_all) == tr.nta.values)
bins_eq = np.mean([np.mean(getattr(tr_p, f"{k}_all") == mine[k]) for k in mine])
else:
nta_eq = bins_eq = np.nan
rows.append(dict(subject=row.subject, session=row.session, trials=n, same_count=same_n,
trials_agree=per_trial.mean(), nta_agree=nta_eq, bins_10ms_agree=bins_eq))
sample_cmp = pd.DataFrame(rows)
print(sample_cmp.round(5).to_string(index=False))
print("\nOver the sample: trials agreeing on context, side, strengths, choice and outcome:",
round(np.average(sample_cmp.trials_agree, weights=sample_cmp.trials), 5),
"| position in block:", round(np.average(sample_cmp.nta_agree, weights=sample_cmp.trials), 5),
"| 10 ms pulse bins:", round(np.average(sample_cmp.bins_10ms_agree, weights=sample_cmp.trials), 5))
# Every matched session, not only the sample: the per-trial fields, and the pulse counts in 50 ms
# bins (the paper's 10 ms bins summed in fives) for the sessions whose pulses were binned here.
idx_of = dict(zip(S.subject + "/" + S.session, S.index))
tot = dict(sessions=0, empty=0, same_count=0, trials=0, trials_agree=0, nta_agree=0, binned_trials=0, bins=0, bins_agree=0)
for rat in rats:
trs = np.atleast_1d(loadmat(_io.BytesIO(blobs[rat]), squeeze_me=True, struct_as_record=False,
variable_names=["trials"])["trials"])
for _, row in both[both.subject == rat].iterrows():
i = idx_of[f"{rat}/{row.session}"]
tr_p, tr = trs[int(row.k)], results[i]["trials"]
tr = tr[tr.hit.notna()]
n = len(tr)
tot["sessions"] += 1
if np.atleast_1d(tr_p.choice_all).size == 0:
tot["empty"] += 1
continue
if n != np.atleast_1d(tr_p.choice_all).size:
continue
tot["same_count"] += 1
ok = ((np.array(list(str(tr_p.task_all))) == tr.task.values) & (np.array(list(str(tr_p.side_all))) == tr.side.values)
& np.isclose(np.atleast_1d(tr_p.gdir_all), tr.gdir) & np.isclose(np.atleast_1d(tr_p.gfreq_all), tr.gfreq)
& (np.atleast_1d(tr_p.choice_all) == tr.choice.values) & (np.atleast_1d(tr_p.hits_all) == tr.hit.values))
tot["trials"] += n; tot["trials_agree"] += int(ok.sum())
tot["nta_agree"] += int((np.atleast_1d(tr_p.nta_all) == tr.nta.values).sum())
if results[i]["pulses"] is not None:
p = results[i]["pulses"].astype(int)
mine = dict(vecr=p[:, 0] + p[:, 2], vecl=p[:, 1] + p[:, 3], vechi=p[:, 0] + p[:, 1], veclo=p[:, 2] + p[:, 3])
for key, v in mine.items():
theirs = np.asarray(getattr(tr_p, f"{key}_all")).reshape(n, -1).reshape(n, 26, 5).sum(2)
tot["bins"] += v.size; tot["bins_agree"] += int((theirs == v).sum())
tot["binned_trials"] += n
print(f"\nAll {tot['sessions']} matched sessions: {tot['empty']} have an empty per-trial record in the paper's tables; of the "
f"other {tot['sessions'] - tot['empty']}, {tot['same_count']} have the same number of valid trials; over their "
f"{tot['trials']} valid trials, context, side, strengths, choice and outcome agree on {tot['trials_agree'] / tot['trials']:.5f} "
f"and position in block on {tot['nta_agree'] / tot['trials']:.5f}; over the {tot['binned_trials']} trials of the "
f"sessions with binned pulses, 50 ms pulse bins agree on {tot['bins_agree'] / tot['bins']:.5f}")
num = ts.paper_mat("data/numsessions_all_rats.mat")
numhits = pd.Series(num["numhits"], index=num["rats"])
cmp_hits = pd.DataFrame(dict(paper_numhits=numhits, dandiset=per_rat.valid_selected)).fillna(0).astype(int)
cmp_hits["difference"] = cmp_hits.dandiset - cmp_hits.paper_numhits
print("\nValid trials in selected sessions, per rat, against the paper's numsessions_all_rats.mat:")
print(cmp_hits.to_string())
print("Rats with equal counts:", int((cmp_hits.difference == 0).sum()), "of", len(cmp_hits))
del blobs
subject session trials same_count trials_agree nta_agree bins_10ms_agree P124 TaskSwitch4-181130a 271 True 1.0 1.0 1.00000 P135 TaskSwitch4-180724a 823 True 1.0 1.0 0.99997 P100 TaskSwitch4-180816a 429 True 1.0 1.0 0.99995 P055 TaskSwitch4-170122a 549 True 1.0 1.0 1.00000 P127 TaskSwitch4-180922a 482 True 1.0 1.0 0.99994 P101 TaskSwitch4-180418a 395 True 1.0 1.0 0.99996 P130 TaskSwitch4-180504a 285 True 1.0 1.0 0.99997 P100 TaskSwitch4-181024a 493 True 1.0 1.0 0.99997 P189 TaskSwitch6-191122a 588 True 1.0 1.0 0.99995 P159 TaskSwitch4-190315a 780 True 1.0 1.0 0.99996 P072 TaskSwitch4-180413a 774 True 1.0 1.0 0.99999 P095 TaskSwitch4-180206a 631 True 1.0 1.0 0.99994 P095 TaskSwitch4-180911a 381 True 1.0 1.0 0.99998 P189 TaskSwitch6-210518a 540 True 1.0 1.0 0.99990 P124 TaskSwitch4-171119a 279 True 1.0 1.0 0.99997 P187 TaskSwitch6-191223a 192 True 1.0 1.0 1.00000 P167 TaskSwitch4-191218a 113 True 1.0 1.0 1.00000 P135 TaskSwitch4-181013a 980 True 1.0 1.0 0.99994 P100 TaskSwitch4-181013a 445 True 1.0 1.0 0.99998 P007 TaskSwitch4-161130a 880 True 1.0 1.0 0.99998 Over the sample: trials agreeing on context, side, strengths, choice and outcome: 1.0 | position in block: 1.0 | 10 ms pulse bins: 0.99996
All 15277 matched sessions: 627 have an empty per-trial record in the paper's tables; of the other 14650, 14650 have the same number of valid trials; over their 7447692 valid trials, context, side, strengths, choice and outcome agree on 1.00000 and position in block on 1.00000; over the 6532005 trials of the sessions with binned pulses, 50 ms pulse bins agree on 0.99996
Valid trials in selected sessions, per rat, against the paper's numsessions_all_rats.mat:
paper_numhits dandiset difference
H113 203339 203339 0
P007 148259 148259 0
P013 233734 233734 0
P049 23014 23014 0
P055 480165 480165 0
P059 25634 25847 213
P072 177844 177844 0
P074 149761 149761 0
P078 147043 147043 0
P095 197061 197061 0
P100 283348 283348 0
P101 100290 100290 0
P102 290572 290572 0
P116 107940 107940 0
P120 319079 319079 0
P123 202690 202690 0
P124 115123 115123 0
P127 195545 195545 0
P130 287944 287944 0
P131 104187 104187 0
P135 251475 251475 0
P150 292255 292255 0
P159 135369 135369 0
P167 128935 128935 0
P187 151999 153441 1442
P189 127700 128309 609
Rats with equal counts: 23 of 26
Result: reproduced. Of the 15,278 paper sessions with valid trials, 15,277 (0.9999) are in the dandiset; the one missing (P189 TaskSwitch6-210516a) has 27 valid trials. In all 15,277 matched sessions the number of valid trials and the fraction correct are equal. In the 20 sampled sessions every valid trial agrees on context, correct side, strengths, choice and outcome, the position in block agrees on every trial, and 0.99996 of the 10 ms pulse bins agree. Over every matched session the same holds: 627 have an empty per-trial record in the paper's tables, and in the other 14,650 the valid trial counts are equal, the 7,447,692 valid trials agree on every per-trial field and on position in block (1.00000), and 0.99996 of the 50 ms pulse bins agree in the sessions whose pulses were binned (6,532,005 trials).
Both relative weights agree to 1e-6 in 0.898 of matched sessions. Of the 1,553 that differ, 861
agree once this attempt's values are capped at 1, which the paper's table does. Of the 692 left,
627 have a missing value on both sides and 5 in the paper's table only; the other 60 differ with
both values present, by a median of 0.87 after capping (39 to 916 valid trials), and were not
explained. The 12 largest of these differences are in sessions with a fraction correct of 0.34 to
0.69. Only 3 of the 692 are in the paper's selection. The
Methods selection made from the dandiset alone keeps all 8,843 sessions the paper's code keeps,
with its good flag, plus 6 more: 5 P187 and P189 sessions absent from the paper's tables, and
P059 TaskSwitch4-170331a, where the location weight in the location context is 1.112 here and
1.0 in the paper's table, so its index falls to 0.697 there. The per-rat trial counts equal the
paper's for 23 of 26 rats and exceed them by 213, 1,442 and 609 trials for P059, P187 and P189.
subjects¶
The study used 26 male Long-Evans rats, aged between 6 and 24 months.
dob = S[S.readable].groupby("subject").date_of_birth.agg(lambda x: sorted(set(x)))
print("Subjects:", assets.subject.nunique(), "| subjects with more than one date of birth:", int((dob.str.len() > 1).sum()))
age = (pd.to_datetime(sel.session_start.str[:10]) - pd.to_datetime(sel.date_of_birth.str[:10])).dt.days / 30.44
print(f"Age at the {len(sel)} selected sessions (months): min {age.min():.1f}, median {age.median():.1f}, max {age.max():.1f}")
inside = ((age >= 5.4) & (age <= 26.4)).mean()
print(f"Fraction within 5.4 to 26.4 months: {inside:.4f}; within 6 to 24: {((age >= 6) & (age <= 24)).mean():.4f}")
by_rat = pd.DataFrame(dict(min=age.groupby(sel.subject).min(), max=age.groupby(sel.subject).max())).round(1)
print(by_rat.T.to_string())
print("Rats with a selected session below 5.4 months:", int((by_rat["min"] < 5.4).sum()), "| above 26.4 months:", int((by_rat["max"] > 26.4).sum()))
# Plausibility of the dates of birth: each rat's age at its first session in the dandiset (any
# protocol), at its first TaskSwitch session and at its first selected session.
day = pd.to_datetime(assets.date, format="%y%m%d")
first_any = day.groupby(assets.subject).min()
first_ts = day[assets.protocol.str.startswith("TaskSwitch")].groupby(assets.subject[assets.protocol.str.startswith("TaskSwitch")]).min()
first_sel = pd.to_datetime(sel.date, format="%y%m%d").groupby(sel.subject).min()
birth = pd.to_datetime(dob.str[0].str[:10])
ages = pd.DataFrame(dict(date_of_birth=birth.dt.date, first_session=first_any.dt.date,
age_first_session=(first_any - birth).dt.days / 30.44,
age_first_taskswitch=(first_ts - birth).dt.days / 30.44,
age_first_selected=(first_sel - birth).dt.days / 30.44)).round(1)
print(ages.to_string())
print("Age at the first session in the dandiset (months): min", ages.age_first_session.min(), "median",
ages.age_first_session.median(), "max", ages.age_first_session.max(),
"| rats younger than 2 months at their first session:", int((ages.age_first_session < 2).sum()))
Subjects: 26 | subjects with more than one date of birth: 0
Age at the 8849 selected sessions (months): min 3.0, median 15.4, max 39.1
Fraction within 5.4 to 26.4 months: 0.9130; within 6 to 24: 0.8441
subject H113 P007 P013 P049 P055 P059 P072 P074 P078 P095 P100 P101 P102 P116 P120 P123 P124 P127 P130 P131 P135 P150 P159 P167 P187 P189
min 4.4 11.4 8.7 7.5 3.6 4.6 4.8 10.9 3.3 4.2 4.8 5.0 4.2 6.4 3.0 3.5 4.7 3.1 4.1 5.7 6.7 4.1 3.7 3.4 4.2 18.0
max 26.0 24.2 27.5 28.7 29.3 25.2 23.5 25.2 19.7 24.3 29.4 23.4 25.0 27.6 22.0 23.0 27.3 23.9 25.4 28.3 23.8 25.6 18.5 20.7 24.3 39.1
Rats with a selected session below 5.4 months: 18 | above 26.4 months: 8
date_of_birth first_session age_first_session age_first_taskswitch age_first_selected
subject
H113 2017-02-21 2017-05-14 2.7 2.8 4.4
P007 2015-02-03 2015-03-26 1.7 4.4 11.4
P013 2015-04-07 2015-05-28 1.7 1.7 8.7
P049 2015-12-15 2016-02-21 2.2 2.4 7.5
P055 2016-02-02 2016-03-27 1.8 2.5 3.6
P059 2016-03-29 2016-05-25 1.9 2.1 4.6
P072 2016-05-03 2016-06-27 1.8 2.0 4.8
P074 2016-05-03 2016-06-27 1.8 2.0 10.9
P078 2016-07-12 2016-09-01 1.7 1.9 3.3
P095 2017-02-21 2017-04-10 1.6 1.9 4.2
P100 2017-02-21 2017-04-10 1.6 1.8 4.8
P101 2017-02-21 2017-04-10 1.6 1.9 5.0
P102 2017-02-21 2017-04-10 1.6 1.7 4.2
P116 2017-04-25 2017-06-25 2.0 2.4 6.4
P120 2017-06-06 2017-07-24 1.6 1.9 3.0
P123 2017-06-06 2017-07-26 1.6 1.9 3.5
P124 2017-06-06 2017-07-26 1.6 1.9 4.7
P127 2017-06-06 2017-07-25 1.6 1.9 3.1
P130 2017-06-27 2017-08-21 1.8 2.2 4.1
P131 2017-06-27 2017-08-21 1.8 2.2 5.7
P135 2017-06-27 2017-08-21 1.8 2.2 6.7
P150 2018-01-09 2018-03-01 1.7 2.0 4.1
P159 2018-03-27 2018-05-17 1.7 1.9 3.7
P167 2018-05-29 2018-07-18 1.6 1.9 3.4
P187 2019-02-21 2019-04-25 2.1 2.4 4.2
P189 2018-02-21 2019-04-25 14.1 14.4 18.0
Age at the first session in the dandiset (months): min 1.6 median 1.7 max 14.1 | rats younger than 2 months at their first session: 21
Result: partially reproduced. The dandiset holds 26 subjects, all male Long Evans. Ages at the 8,849 selected sessions run from 3.0 to 39.1 months (median 15.4); 0.913 of sessions lie within the widened 5.4 to 26.4 months and 0.844 within 6 to 24. Eighteen rats have selected sessions before 5.4 months (the youngest at 3.0) and eight after 26.4 (P189 up to 39.1).
The dates of birth put 25 rats at 1.6 to 2.7 months at their first session in the dandiset (21 below 2 months), and 3.0 to 11.4 months at their first selected session. P189 is the exception: 14.1 months at its first session, which falls on the same day as P187's, whose date of birth is 2019-02-21 against P189's 2018-02-21. P189's date of birth may be a year early; with 2019-02-21 its first session would be at about 2.1 months. For the other rats the ages at the first session are those of rats that start training young, so the selected sessions below 6 months are likely real, and the Methods' 6 to 24 months does not describe them.
trials-per-session, blocks-per-session, session-duration¶
Rats performed on average 542 trials per session. On average, rats switched across 14.6 context blocks per session. Behavioural sessions lasted 2 to 4 hours.
per_session = trials_all.groupby("sid").agg(trials=("hit", "size"), valid=("hit", "count"),
blocks=("block", "max"), first=("start", "min"), last=("stop", "max"))
per_session["blocks"] += 1
per_session["switches"] = per_session.blocks - 1
per_session["hours"] = (per_session["last"] - per_session["first"]) / 3600
per_session["selected"] = S.selected.reindex(per_session.index).values
groups = {"selected sessions": per_session[per_session.selected], "all task sessions": per_session}
rows = []
for name, g in groups.items():
rows.append(dict(sessions=name, n=len(g), mean_trials=g.trials.mean(), mean_valid=g.valid.mean(),
mean_blocks=g.blocks.mean(), mean_switches=g.switches.mean(),
median_hours=g.hours.median(), frac_2_to_4h=((g.hours >= 2) & (g.hours <= 4)).mean(),
frac_1p8_to_4p4h=((g.hours >= 1.8) & (g.hours <= 4.4)).mean(),
hours_5_95=f"{g.hours.quantile(.05):.2f} to {g.hours.quantile(.95):.2f}"))
summary_sessions = pd.DataFrame(rows).set_index("sessions")
print(summary_sessions.round(3).T.to_string())
unread = S[~S.readable]
n_all = len(per_session) + len(unread)
print(f"\nWith the {len(unread)} unreadable sessions ({unread.n_trials.min()} to {unread.n_trials.max()} trials) added: reading (c) "
f"{(per_session.trials.sum() + unread.n_trials.sum()) / n_all:.1f} trials; reading (d), counting them as 0 valid trials, "
f"{per_session.valid.sum() / n_all:.1f}")
print("\nPaper: 542 trials (tolerance 487.8 to 596.2); 14.6 blocks (13.14 to 16.06); 2 to 4 h (median within 1.8 to 4.4 h)")
sessions selected sessions all task sessions n 8849 15400 mean_trials 817.851 749.685 mean_valid 551.765 494.227 mean_blocks 16.975 13.898 mean_switches 15.975 12.898 median_hours 2.499 2.569 frac_2_to_4h 0.559 0.556 frac_1p8_to_4p4h 0.65 0.65 hours_5_95 1.26 to 4.41 1.25 to 4.48 With the 13 unreadable sessions (1 to 3 trials) added: reading (c) 749.1 trials; reading (d), counting them as 0 valid trials, 493.8 Paper: 542 trials (tolerance 487.8 to 596.2); 14.6 blocks (13.14 to 16.06); 2 to 4 h (median within 1.8 to 4.4 h)
fig, ax = plt.subplots(1, 3, figsize=(12, 3))
g = per_session[per_session.selected]
ax[0].hist(g.trials, 60); ax[0].set_xlabel("trials per selected session")
ax[1].hist(g.blocks, np.arange(0, 60)); ax[1].set_xlabel("context blocks per selected session")
ax[2].hist(g.hours.clip(0, 8), 60); ax[2].set_xlabel("session duration (h, clipped at 8)")
plt.tight_layout(); plt.show()
print("Ranges: trials", g.trials.min(), "to", g.trials.max(), "| blocks", g.blocks.min(), "to", g.blocks.max(),
"| hours", round(g.hours.min(), 2), "to", round(g.hours.max(), 2), "; sessions above 8 h:", int((g.hours > 8).sum()))
Ranges: trials 130 to 2543 | blocks 3 to 59 | hours 0.26 to 7.74 ; sessions above 8 h: 0
Result: trials-per-session reproduced. Valid trials per session average 551.8 over the selected sessions (reading b) and 494.2 over all task sessions (reading d), both within 487.8 to 596.2; all trials average 817.9 and 749.7 (readings a and c), outside. The 13 unreadable sessions (1 to 3 trials) are left out of readings (c) and (d); added, they give 749.1 and 493.8.
Result: blocks-per-session reproduced. The number of switches averages 15.98 over the selected sessions and the number of blocks 13.90 over all task sessions, both within 13.14 to 16.06; blocks over selected sessions (16.98) and switches over all sessions (12.90) are outside.
Result: session-duration reproduced. The median duration is 2.50 h over the selected sessions and 2.57 h over all task sessions, within 1.8 to 4.4 h. Only 0.559 and 0.556 of sessions last 2 to 4 h, and 5% to 95% of selected sessions span 1.26 to 4.41 h.
block-min-30¶
Block switches occurred after a minimum of 30 trials per block, and when a local estimate of performance reached 80% correct.
rows = []
for i in np.flatnonzero(S.selected):
tr = results[i]["trials"]
last = tr.block.max()
for b, g in tr[tr.block < last].groupby("block"):
v = g.hit.dropna()
rows.append(dict(sid=i, all=len(g), valid=len(v), perf_last30=v.tail(30).mean()))
blocks = pd.DataFrame(rows)
print("Non-final blocks in selected sessions:", len(blocks))
print(f"Blocks with >= 30 trials (all trials): {(blocks['all'] >= 30).mean():.4f}; (valid trials): {(blocks.valid >= 30).mean():.4f}")
print("Block length, all trials: percentiles 1, 5, 50:", blocks["all"].quantile([.01, .05, .5]).tolist(),
"| valid trials:", blocks.valid.quantile([.01, .05, .5]).tolist())
print("Fraction correct over the last 30 valid trials before a switch: median", round(blocks.perf_last30.median(), 3),
"| fraction at or above 0.8:", round((blocks.perf_last30 >= 0.8).mean(), 3))
Non-final blocks in selected sessions: 141359 Blocks with >= 30 trials (all trials): 1.0000; (valid trials): 0.8170 Block length, all trials: percentiles 1, 5, 50: [30.0, 32.0, 44.0] | valid trials: [29.0, 29.0, 30.0] Fraction correct over the last 30 valid trials before a switch: median 0.828 | fraction at or above 0.8: 0.601
Result: reproduced. Counted over all trials, every one of the 141,359 non-final blocks in the selected sessions has at least 30 trials (fraction 1.0000; 1st percentile 30, median 44); counted over valid trials, 0.817 do. Over the last 30 valid trials before a switch the median fraction correct is 0.828, and 0.601 of blocks are at or above 0.8, so the 80% rule is not reached on that window in 40% of blocks; the task's own estimate is not described.
fig1c-20-rats and fig1c-context-selection¶
Twenty rats each contributed more than 120,000 trials after training. After training, choices are mostly driven by location in the location context and by frequency in the frequency context.
counts = per_rat.valid_selected.sort_values(ascending=False)
rats120 = counts[counts > 120000].index.tolist()
print("Rats with more than 120,000 valid trials in selected sessions:", len(rats120),
"| with at least 120,000 (Figure1c.m uses >=):", int((counts >= 120000).sum()))
print(counts.to_string())
Rats with more than 120,000 valid trials in selected sessions: 20 | with at least 120,000 (Figure1c.m uses >=): 20 subject P055 480165.0 P120 319079.0 P150 292255.0 P102 290572.0 P130 287944.0 P100 283348.0 P135 251475.0 P013 233734.0 H113 203339.0 P123 202690.0 P095 197061.0 P127 195545.0 P072 177844.0 P187 153441.0 P074 149761.0 P007 148259.0 P078 147043.0 P159 135369.0 P167 128935.0 P189 128309.0 P124 115123.0 P116 107940.0 P131 104187.0 P101 100290.0 P059 25847.0 P049 23014.0
Result: reproduced. Twenty rats have more than 120,000 valid trials in the selected sessions (P189 is the twentieth, with 128,309; P124 is next, with 115,123).
rows = []
for rat in rats120:
v = pd.concat([results[i]["trials"] for i in np.flatnonzero(S.selected & (S.subject == rat))])
v = v[v.hit.notna()]
rl, rf, fsi = ts.relative_weights(v.choice.values.astype(int), v.task.values, v.gdir.values, v.gfreq.values)
rows.append(dict(rat=rat, trials=len(v), rel_location_in_LOC=rl, rel_frequency_in_FRQ=rf, index=fsi))
ctx_sel = pd.DataFrame(rows).set_index("rat")
print(ctx_sel.round(3).to_string())
both_ok = (ctx_sel.rel_location_in_LOC > 0.5) & (ctx_sel.rel_frequency_in_FRQ > 0.5)
print("\nRats with both relative weights above 0.5:", int(both_ok.sum()), "of", len(ctx_sel),
"| lowest:", round(ctx_sel[["rel_location_in_LOC", "rel_frequency_in_FRQ"]].min().min(), 3))
trials rel_location_in_LOC rel_frequency_in_FRQ index rat P055 480165 0.836 0.825 0.831 P120 319079 0.860 0.836 0.848 P150 292255 0.797 0.789 0.793 P102 290572 0.856 0.779 0.817 P130 287944 0.857 0.870 0.863 P100 283348 0.841 0.767 0.804 P135 251475 0.830 0.788 0.809 P013 233734 0.795 0.750 0.772 H113 203339 0.834 0.765 0.800 P123 202690 0.844 0.786 0.815 P095 197061 0.804 0.724 0.764 P127 195545 0.829 0.802 0.816 P072 177844 0.810 0.775 0.792 P187 153441 0.842 0.804 0.823 P074 149761 0.770 0.731 0.750 P007 148259 0.804 0.724 0.764 P078 147043 0.808 0.788 0.798 P159 135369 0.805 0.783 0.794 P167 128935 0.844 0.760 0.802 P189 128309 0.825 0.765 0.795 Rats with both relative weights above 0.5: 20 of 20 | lowest: 0.724
Result: reproduced. In all 20 rats the relevant feature carries more than half of the weight in both contexts: the relative weight of location in the location context is 0.770 to 0.860 and of frequency in the frequency context 0.724 to 0.870.
ed1e-switch-four-trials¶
Rats switch the feature they attend to within about four trials of a new context block.
def switch_curves(mask, min_n=30):
rows = []
for rat in sorted(S.subject.unique()):
ids = np.flatnonzero(mask & (S.subject == rat))
if not len(ids):
continue
v = pd.concat([results[i]["trials"] for i in ids])
v = v[v.hit.notna() & (v.block > 0)]
for new in ("f", "d"):
for k in range(1, 31):
g = v[(v.task == new) & (v.nta == k)]
if len(g) < min_n:
continue
X = np.column_stack([np.ones(len(g)), g.gdir, g.gfreq])
w = ts._irls(X, g.choice.values.astype(float))
rel = w[2] / (w[1] + w[2]) if new == "f" else w[1] / (w[1] + w[2])
rows.append(dict(rat=rat, new=new, k=k, rel=rel, n=len(g)))
return pd.DataFrame(rows)
def switch_summary(c):
m = c.groupby(["new", "k"]).rel.mean().unstack(0)
plateau = m.loc[10:30].mean()
return m, plateau
main_mask = S.selected
code_mask = S.readable & (S.valid >= 100) & (S.hit >= 0.65) & (S.fsi >= 0.65)
for name, mask in (("Methods selection (main)", main_mask), ("code's 0.65 selection (variant)", code_mask)):
c = switch_curves(mask)
m, plateau = switch_summary(c)
print(f"{name}: rats {c.rat.nunique()}; relevant relative weight by position in the new block "
f"(new block f = switch from location to frequency):")
print(m.loc[[1, 2, 3, 4, 5, 6, 8, 10, 20, 30]].round(3).T.to_string())
print("plateau (positions 10 to 30):", plateau.round(3).to_dict())
for k in (4, 5):
print(f" position {k}: fraction of plateau", (m.loc[k] / plateau).round(3).to_dict())
if name.startswith("Methods"):
switch_main = (c, m, plateau)
c, m, plateau = switch_main
# Per rat at positions 1 to 5 (Methods selection): relevant relative weight and trials fitted.
per_rat_k = c[c.k <= 5].pivot_table(index="rat", columns=["new", "k"], values=["rel", "n"])
print("\nPer rat, positions 1 to 5 (Methods selection):")
print(per_rat_k["rel"].round(2).to_string())
print("trials per fit, min and max over rats:", per_rat_k["n"].min().min(), "to", per_rat_k["n"].max().max(),
"| rats with a fit at each position:", c[c.k <= 5].groupby(["new", "k"]).rat.nunique().unstack(0).to_dict())
# Bootstrap over rats of the fraction of plateau at positions 4 and 5.
R = c.pivot_table(index=["new", "k"], columns="rat", values="rel")
rng_sw = np.random.default_rng(301)
boot = []
for _ in range(2000):
pick = rng_sw.choice(R.columns, len(R.columns), replace=True)
mm = R[pick].mean(axis=1).unstack(0)
pl = mm.loc[10:30].mean()
boot.append({(new, k): mm.loc[k, new] / pl[new] for new in ("d", "f") for k in (4, 5)})
boot = pd.DataFrame(boot)
print("\nBootstrap over rats (2000 resamples), fraction of plateau: 2.5th, 50th and 97.5th percentiles; share of resamples >= 0.9")
print(pd.DataFrame(dict(p2_5=boot.quantile(0.025), p50=boot.median(), p97_5=boot.quantile(0.975),
share_ge_0_9=(boot >= 0.9).mean())).round(3).to_string())
for min_n in (10, 100):
cc_ = switch_curves(main_mask, min_n)
mv_, pv_ = switch_summary(cc_)
print(f"Minimum trials per fit {min_n} (main uses 30): position 4 fraction of plateau",
(mv_.loc[4] / pv_).round(3).to_dict(), "| position 5", (mv_.loc[5] / pv_).round(3).to_dict(),
"| rat x position fits dropped:", 26 * 2 * 30 - len(cc_))
fig, ax = plt.subplots(figsize=(5, 3))
for new, lab in (("f", "location to frequency"), ("d", "frequency to location")):
ax.plot(m.index, m[new], label=lab)
ax.set_xlabel("position in new block (all trials)"); ax.set_ylabel("relevant relative weight"); ax.legend(); plt.show()
print("Plotted range:", round(m.min().min(), 3), "to", round(m.max().max(), 3))
Methods selection (main): rats 26; relevant relative weight by position in the new block (new block f = switch from location to frequency):
k 1 2 3 4 5 6 8 10 20 30
new
d 0.502 0.642 0.735 0.776 0.816 0.832 0.866 0.871 0.851 0.828
f 0.560 0.683 0.729 0.760 0.771 0.775 0.781 0.787 0.787 0.787
plateau (positions 10 to 30): {'d': 0.854, 'f': 0.789}
position 4: fraction of plateau {'d': 0.91, 'f': 0.963}
position 5: fraction of plateau {'d': 0.956, 'f': 0.977}
code's 0.65 selection (variant): rats 26; relevant relative weight by position in the new block (new block f = switch from location to frequency):
k 1 2 3 4 5 6 8 10 20 30
new
d 0.486 0.636 0.704 0.759 0.793 0.806 0.838 0.847 0.829 0.809
f 0.552 0.670 0.716 0.743 0.750 0.757 0.771 0.772 0.771 0.771
plateau (positions 10 to 30): {'d': 0.831, 'f': 0.773}
position 4: fraction of plateau {'d': 0.914, 'f': 0.961}
position 5: fraction of plateau {'d': 0.955, 'f': 0.97}
Per rat, positions 1 to 5 (Methods selection):
new d f
k 1 2 3 4 5 1 2 3 4 5
rat
H113 0.46 0.68 0.78 0.87 0.92 0.57 0.69 0.75 0.75 0.76
P007 0.50 0.62 0.76 0.70 0.74 0.59 0.70 0.72 0.72 0.76
P013 0.44 0.61 0.66 0.73 0.74 0.60 0.68 0.71 0.71 0.72
P049 0.75 0.78 0.96 0.75 0.92 0.63 0.64 0.75 0.79 0.71
P055 0.51 0.66 0.76 0.81 0.85 0.57 0.72 0.78 0.81 0.84
P059 0.76 0.63 0.90 0.74 0.80 0.59 0.68 0.69 0.71 0.74
P072 0.80 0.94 0.96 0.88 0.92 0.60 0.71 0.74 0.76 0.79
P074 0.58 0.65 0.72 0.77 0.76 0.53 0.66 0.67 0.73 0.78
P078 0.34 0.49 0.61 0.68 0.70 0.54 0.71 0.72 0.77 0.81
P095 0.49 0.59 0.66 0.76 0.83 0.54 0.68 0.72 0.69 0.73
P100 0.38 0.57 0.67 0.72 0.77 0.52 0.68 0.71 0.75 0.75
P101 0.49 0.68 0.78 0.87 0.90 0.61 0.69 0.75 0.72 0.74
P102 0.33 0.55 0.65 0.68 0.79 0.56 0.71 0.71 0.76 0.79
P116 0.63 0.86 0.83 0.90 0.93 0.57 0.65 0.73 0.76 0.74
P120 0.46 0.65 0.72 0.79 0.83 0.62 0.72 0.78 0.81 0.84
P123 0.37 0.55 0.64 0.76 0.78 0.49 0.69 0.74 0.79 0.77
P124 0.43 0.58 0.63 0.73 0.69 0.51 0.65 0.68 0.73 0.76
P127 0.33 0.48 0.59 0.72 0.70 0.47 0.66 0.73 0.78 0.80
P130 0.28 0.54 0.71 0.75 0.84 0.58 0.74 0.81 0.85 0.86
P131 0.80 0.79 0.93 0.91 1.01 0.55 0.65 0.69 0.73 0.74
P135 0.43 0.61 0.68 0.74 0.79 0.49 0.67 0.74 0.80 0.78
P150 0.31 0.49 0.58 0.70 0.78 0.63 0.72 0.75 0.76 0.76
P159 0.56 0.69 0.72 0.83 0.80 0.62 0.71 0.72 0.81 0.78
P167 0.49 0.63 0.70 0.75 0.78 0.48 0.60 0.69 0.74 0.72
P187 0.56 0.71 0.78 0.82 0.85 0.52 0.67 0.76 0.77 0.81
P189 0.55 0.67 0.73 0.82 0.82 0.59 0.68 0.72 0.77 0.78
trials per fit, min and max over rats: 103.0 to 6096.0 | rats with a fit at each position: {'d': {1: 26, 2: 26, 3: 26, 4: 26, 5: 26}, 'f': {1: 26, 2: 26, 3: 26, 4: 26, 5: 26}}
Bootstrap over rats (2000 resamples), fraction of plateau: 2.5th, 50th and 97.5th percentiles; share of resamples >= 0.9
p2_5 p50 p97_5 share_ge_0_9
(d, 4) 0.874 0.909 0.946 0.678
(d, 5) 0.916 0.956 1.000 0.998
(f, 4) 0.949 0.962 0.981 1.000
(f, 5) 0.967 0.977 0.989 1.000
Minimum trials per fit 10 (main uses 30): position 4 fraction of plateau {'d': 0.91, 'f': 0.963} | position 5 {'d': 0.956, 'f': 0.977} | rat x position fits dropped: 0
Minimum trials per fit 100 (main uses 30): position 4 fraction of plateau {'d': 0.91, 'f': 0.963} | position 5 {'d': 0.956, 'f': 0.977} | rat x position fits dropped: 0
Plotted range: 0.502 to 0.877
Result: reproduced. At position 4 of a new block the across-rat relative weight of the newly relevant feature is 0.91 (frequency to location) and 0.963 (location to frequency) of its plateau under the Methods selection, and 0.914 and 0.961 under the code's looser selection. At position 1 it is 0.502 and 0.560, against plateaus of 0.854 and 0.789. All 26 rats have a fit at every position from 1 to 5, with 103 to 6,096 trials each, so the minimum of 30 trials drops nothing (10 and 100 give the same values). The position 4 reading is thin for the switch to location: bootstrapping the rats gives 0.874 to 0.946 (95%), and 0.678 of resamples reach 0.9. The position 5 reading, which the criterion also allows, is firm: 0.956 and 0.977 of plateau, with 0.998 and 1.000 of resamples at 0.9 or above.
ed2g-learning-2-5-months¶
Most rats learned the task in two to five months.
all_sessions = pd.concat([S[S.readable], S_other[S_other.missing.fillna("") == ""]], ignore_index=True)
all_sessions["day"] = pd.to_datetime(all_sessions.date, format="%y%m%d")
assets["day"] = pd.to_datetime(assets.date, format="%y%m%d")
start_any = assets.groupby("subject").day.min()
start_ts = assets[assets.protocol.str.startswith("TaskSwitch")].groupby("subject").day.min()
crit_a = all_sessions[(all_sessions.valid >= 100) & (all_sessions.fsi >= 0.7)].groupby("subject").day.min()
crit_b = all_sessions[(all_sessions.valid >= 100) & (all_sessions.fsi >= 0.7) & (all_sessions.hit >= 0.7)].groupby("subject").day.min()
learn = pd.DataFrame({f"{s}_{c}": (cv - sv).dt.days / 30.44
for s, sv in (("i_first_session", start_any), ("ii_first_taskswitch", start_ts))
for c, cv in (("a_index", crit_a), ("b_selection", crit_b))})
print(learn.round(2).to_string())
print("\nFraction of rats within 1.8 to 5.5 months, per reading:")
print(((learn >= 1.8) & (learn <= 5.5)).sum().astype(str) + " of " + learn.notna().sum().astype(str))
print("Within 2 to 5 months:", (((learn >= 2) & (learn <= 5)).sum()).to_dict())
print("Protocols of each rat's first session:", assets.sort_values("day").groupby("subject").protocol.first().value_counts().to_dict())
i_first_session_a_index i_first_session_b_selection ii_first_taskswitch_a_index ii_first_taskswitch_b_selection
subject
H113 1.74 1.74 1.61 1.61
P007 8.18 9.69 5.49 7.00
P013 4.93 4.93 4.93 4.93
P049 2.60 5.29 2.43 5.12
P055 1.81 1.81 1.12 1.12
P059 2.69 2.69 2.50 2.50
P072 1.68 2.99 1.51 2.83
P074 3.42 9.07 3.25 8.90
P078 1.64 1.64 1.41 1.41
P095 2.60 2.63 2.27 2.30
P100 1.84 3.19 1.58 2.92
P101 3.32 3.42 3.02 3.12
P102 2.37 2.66 2.23 2.53
P116 3.81 4.37 3.42 3.98
P120 1.45 1.45 1.08 1.08
P123 1.84 1.84 1.54 1.54
P124 3.06 3.06 2.76 2.76
P127 1.45 1.45 1.12 1.12
P130 1.41 2.33 1.05 1.97
P131 2.23 3.88 1.87 3.52
P135 0.85 4.93 0.49 4.57
P150 1.87 2.43 1.51 2.07
P159 1.87 2.00 1.68 1.81
P167 1.77 1.77 1.54 1.54
P187 2.17 2.17 1.87 1.87
P189 3.12 3.94 2.83 3.65
Fraction of rats within 1.8 to 5.5 months, per reading:
i_first_session_a_index 17 of 26
i_first_session_b_selection 19 of 26
ii_first_taskswitch_a_index 13 of 26
ii_first_taskswitch_b_selection 17 of 26
dtype: str
Within 2 to 5 months: {'i_first_session_a_index': 12, 'i_first_session_b_selection': 16, 'ii_first_taskswitch_a_index': 10, 'ii_first_taskswitch_b_selection': 13}
Protocols of each rat's first session: {'ProAnti3Marino': 24, 'ProAnti3': 1, 'TaskSwitch': 1}
Result: reproduced. Counting from each rat's first session and taking the first session with an index of at least 0.7 (readings i and a), 17 of 26 rats took 1.8 to 5.5 months; with the full selection (i, b) 19 of 26, and from the first TaskSwitch session 13 (ii, a, exactly half) and 17 (ii, b). Within the unwidened 2 to 5 months the counts are 12, 16, 10 and 13. Times under reading (i, a) run from 0.85 to 8.18 months.
Behavioural kernels¶
The kernels of fig4-converging, fig4-loc-frq-uncorrelated, ed3g-halfsplit and fig5d-behaviour-neural are fitted here. The main fit uses the code's smoothness penalty with its weight chosen per rat and context by 3-fold cross-validation over 20, 200 and 2000; variants use the code's fixed 200, and for the recorded rats the code's looser selection (0.65 correct, 0.6 index). The cross-validation folds are drawn at random over trials, as the claims file fixes, so trials of one session can fall in different folds; this only affects which penalty weight is chosen, and the fixed-200 variant shows how much that matters for each claim.
def rat_data(rat, mask):
X, y, task = [], [], []
for i in np.flatnonzero(mask & (S.subject == rat)):
r = results[i]
if r["pulses"] is None:
raise RuntimeError(f"pulses not kept for {S.session[i]}")
tr = r["trials"][r["trials"].hit.notna()]
X.append(ts.evidence_features(r["pulses"]))
y.append(tr.choice.values)
task.append(tr.task.values)
return np.vstack(X), np.concatenate(y), np.concatenate(task)
rats100 = counts[counts >= 100000].index.tolist()
rats_fit = sorted(set(rats100) | set(ts.PAPER_EPHYS_RATS))
code_ephys_mask = S.readable & (S.valid >= 100) & (S.hit >= 0.65) & (S.fsi >= 0.6)
t0 = time.time()
def fit_rat(rat):
seed = int(rat[1:]) + 1000 # each rat's folds have their own seed
X, y, task = rat_data(rat, S.selected)
out = dict(rat=rat, n=len(y), main=ts.kernels(X, y, task, seed=seed))
out["fixed200"] = ts.kernels(X, y, task, lam=200.0)
if rat in rats100:
halves = []
for h in (0, 1):
Xh, yh, th = [], [], []
for ctx in ("d", "f"):
idx = np.flatnonzero(task == ctx)
cut = int(round(len(idx) / 2))
part = idx[:cut] if h == 0 else idx[cut:]
Xh.append(X[part]); yh.append(y[part]); th.append(task[part])
Xh, yh, th = np.vstack(Xh), np.concatenate(yh), np.concatenate(th)
halves.append((ts.kernels(Xh, yh, th, lam={c: out["main"][f"lam_{c}"] for c in "df"}),
ts.kernels(Xh, yh, th, lam=200.0)))
out["halves"] = halves
if rat in ts.PAPER_EPHYS_RATS:
Xc, yc, tc = rat_data(rat, code_ephys_mask)
out["code_sel_n"] = len(yc)
out["code_sel_200"] = ts.kernels(Xc, yc, tc, lam=200.0)
out["code_sel_cv"] = ts.kernels(Xc, yc, tc, seed=seed)
return out
with ThreadPoolExecutor(4) as ex:
fits = {f["rat"]: f for f in ex.map(fit_rat, rats_fit)}
print(f"fitted {len(fits)} rats in {time.time() - t0:.0f} s; every fit converged (a failed fit raises)")
print(pd.DataFrame({r: dict(trials=f["n"], lam_LOC=f["main"]["lam_d"], lam_FRQ=f["main"]["lam_f"])
for r, f in fits.items()}).T.to_string())
fitted 26 rats in 188 s; every fit converged (a failed fit raises)
trials lam_LOC lam_FRQ
H113 203339.0 200.0 20.0
P007 148259.0 200.0 200.0
P013 233734.0 200.0 200.0
P049 23014.0 2000.0 2000.0
P055 480165.0 20.0 20.0
P059 25847.0 200.0 200.0
P072 177844.0 20.0 200.0
P074 149761.0 20.0 200.0
P078 147043.0 200.0 200.0
P095 197061.0 20.0 200.0
P100 283348.0 20.0 20.0
P101 100290.0 200.0 200.0
P102 290572.0 200.0 200.0
P116 107940.0 200.0 200.0
P120 319079.0 200.0 200.0
P123 202690.0 200.0 200.0
P124 115123.0 200.0 200.0
P127 195545.0 200.0 200.0
P130 287944.0 20.0 20.0
P131 104187.0 200.0 200.0
P135 251475.0 20.0 20.0
P150 292255.0 200.0 200.0
P159 135369.0 20.0 200.0
P167 128935.0 20.0 200.0
P187 153441.0 20.0 200.0
P189 128309.0 2000.0 200.0
fig4-converging¶
Differential behavioural kernels predominantly converge towards the end of the pulse train, that is the difference between contexts in the weight of a feature shrinks for later pulses.
rats18 = [r for r in rats120 if int(r[1:]) >= 15]
print("Figure 4 rats (> 120,000 valid trials, not numbered below 15):", len(rats18), rats18)
saved = {}
for rat in rats18:
p = ts.paper_mat(f"data/logistic_regression_results/logreg_best_{rat}_5.mat")
saved[rat] = ts.slopes(dict(dd=p["thetad"][1:27], fd=p["thetad"][27:], df=p["thetaf"][1:27], ff=p["thetaf"][27:]))
sl = pd.DataFrame({r: dict(loc=ts.slopes(fits[r]["main"])[0], frq=ts.slopes(fits[r]["main"])[1],
loc_200=ts.slopes(fits[r]["fixed200"])[0], frq_200=ts.slopes(fits[r]["fixed200"])[1],
loc_paper_weights=saved[r][0], frq_paper_weights=saved[r][1]) for r in rats18}).T
print(sl.round(3).to_string())
for a, b in (("loc", "frq"), ("loc_200", "frq_200")):
v = np.r_[sl[a], sl[b]]
print(f"{a[:3]}+{b[:3]} ({'main' if a == 'loc' else 'fixed 200'}): negative {int((v < 0).sum())} of {len(v)}, "
f"median {np.median(v):.3f}, range {v.min():.3f} to {v.max():.3f}, within -0.8 to -0.15: {int(((v >= -0.8) & (v <= -0.15)).sum())}")
mine = np.r_[sl.loc_200, sl.frq_200]; theirs = np.r_[sl.loc_paper_weights, sl.frq_paper_weights]
print(f"Fixed-200 slope indices against those from the paper's saved weights: largest absolute difference "
f"{np.max(np.abs(mine - theirs)):.4f}, r = {np.corrcoef(mine, theirs)[0, 1]:.5f}")
Figure 4 rats (> 120,000 valid trials, not numbered below 15): 18 ['P055', 'P120', 'P150', 'P102', 'P130', 'P100', 'P135', 'H113', 'P123', 'P095', 'P127', 'P072', 'P187', 'P074', 'P078', 'P159', 'P167', 'P189']
loc frq loc_200 frq_200 loc_paper_weights frq_paper_weights P055 -0.569 -0.510 -0.626 -0.523 -0.626 -0.522 P120 -0.441 -0.502 -0.441 -0.502 -0.441 -0.502 P150 -0.302 -0.640 -0.302 -0.640 -0.302 -0.640 P102 -0.725 -0.493 -0.725 -0.493 -0.726 -0.493 P130 -0.257 -0.283 -0.288 -0.279 -0.288 -0.279 P100 -0.512 -0.329 -0.639 -0.333 -0.638 -0.333 P135 -0.355 -0.473 -0.390 -0.482 -0.389 -0.482 H113 -0.539 -0.532 -0.529 -0.484 -0.527 -0.485 P123 -0.774 -0.621 -0.774 -0.621 -0.775 -0.621 P095 -0.426 -0.312 -0.614 -0.373 -0.613 -0.373 P127 -0.233 -0.631 -0.233 -0.631 -0.234 -0.631 P072 -0.603 -0.330 -0.664 -0.405 -0.664 -0.405 P187 -0.428 -0.287 -0.494 -0.301 -0.512 -0.294 P074 -0.543 -0.189 -0.586 -0.258 -0.585 -0.258 P078 -0.437 -0.240 -0.437 -0.240 -0.436 -0.239 P159 -0.390 -0.397 -0.506 -0.489 -0.506 -0.489 P167 -0.424 -0.582 -0.539 -0.635 -0.539 -0.635 P189 -0.097 -0.238 -0.270 -0.224 -0.263 -0.216 loc+frq (main): negative 36 of 36, median -0.432, range -0.774 to -0.097, within -0.8 to -0.15: 35 loc+frq (fixed 200): negative 36 of 36, median -0.491, range -0.774 to -0.224, within -0.8 to -0.15: 36 Fixed-200 slope indices against those from the paper's saved weights: largest absolute difference 0.0177, r = 0.99973
fig, ax = plt.subplots(1, 2, figsize=(10, 3))
lo, hi = np.inf, -np.inf
for rat in rats18:
k = fits[rat]["main"]
for a, (rel, irr) in zip(ax, (("dd", "df"), ("ff", "fd"))):
d = (k[rel] - k[irr]) / np.max(np.abs(k[rel] - k[irr]))
a.plot(ts.BIN_CENTRES, d, lw=0.8)
lo, hi = min(lo, d.min()), max(hi, d.max())
for a, t in zip(ax, ("location", "frequency")):
a.set_title(f"differential kernel, {t} (normalised)"); a.set_xlabel("time in pulse train (s)")
plt.tight_layout(); plt.show()
print("Plotted range:", round(lo, 3), "to", round(hi, 3))
Plotted range: -0.186 to 1.0
Result: reproduced. All 36 slope indices of the 18 rats are negative (median -0.432, range -0.774 to -0.097; 35 within the code's axis range of -0.8 to -0.15). With the code's fixed penalty all 36 are negative as well (median -0.491), and these agree with the indices from the paper's saved weights to at most 0.0177 (r = 0.99973).
fig4-loc-frq-uncorrelated¶
Across rats, the shape of the differential behavioural kernel for one feature does not predict the shape for the other feature.
Result: reproduced. Across the 18 rats the location and frequency slope indices are not significantly correlated: r = 0.180 (95% CI -0.313 to 0.597), p = 0.475; with the fixed penalty r = 0.092, p = 0.717.
ed3g-halfsplit¶
The behavioural slope index computed from the first half of a rat's trials is correlated with the index from the second half (r = 0.58, p = 0.000013).
rows = []
for rat in rats100:
(m1, f1), (m2, f2) = fits[rat]["halves"]
rows.append(dict(rat=rat, loc1=ts.slopes(m1)[0], frq1=ts.slopes(m1)[1], loc2=ts.slopes(m2)[0], frq2=ts.slopes(m2)[1],
loc1_200=ts.slopes(f1)[0], frq1_200=ts.slopes(f1)[1], loc2_200=ts.slopes(f2)[0], frq2_200=ts.slopes(f2)[1]))
hs = pd.DataFrame(rows).set_index("rat")
print(hs.round(3).to_string())
for suf, name in (("", "main"), ("_200", "fixed 200")):
a = np.r_[hs[f"loc1{suf}"], hs[f"frq1{suf}"]]; b = np.r_[hs[f"loc2{suf}"], hs[f"frq2{suf}"]]
r, p, lo, hi = pearson_ci(a, b)
print(f"{name}: {len(rats100)} rats, {len(a)} points: r = {r:.3f} (95% CI {lo:.3f} to {hi:.3f}), p = {p:.3g}")
for f in ("loc", "frq"):
r1, p1 = stats.pearsonr(hs[f"{f}1{suf}"], hs[f"{f}2{suf}"])
print(f" {f} alone: r = {r1:.3f}, p = {p1:.3g}")
print("Paper: r = 0.58 (tolerance 0.522 to 0.638), p = 0.000013")
loc1 frq1 loc2 frq2 loc1_200 frq1_200 loc2_200 frq2_200 rat P055 -0.620 -0.486 -0.500 -0.533 -0.681 -0.518 -0.555 -0.530 P120 -0.472 -0.533 -0.424 -0.479 -0.472 -0.533 -0.424 -0.479 P150 0.206 -0.574 -0.536 -0.700 0.206 -0.574 -0.536 -0.700 P102 -0.850 -0.486 -0.624 -0.462 -0.850 -0.486 -0.624 -0.462 P130 -0.093 -0.163 -0.379 -0.381 -0.116 -0.161 -0.420 -0.379 P100 -0.585 -0.201 -0.464 -0.436 -0.715 -0.196 -0.582 -0.448 P135 0.032 -0.488 -0.487 -0.460 0.043 -0.489 -0.600 -0.469 P013 0.120 -0.953 0.305 -0.609 0.120 -0.953 0.305 -0.609 H113 -0.817 -0.483 -0.200 -0.536 -0.688 -0.464 -0.334 -0.511 P123 -0.663 -0.573 -0.851 -0.672 -0.663 -0.573 -0.851 -0.672 P095 -0.471 -0.236 -0.403 -0.370 -0.723 -0.252 -0.560 -0.474 P127 -0.359 -0.519 -0.154 -0.751 -0.359 -0.519 -0.154 -0.751 P072 -0.690 -0.267 -0.508 -0.394 -0.753 -0.334 -0.561 -0.474 P187 -0.403 -0.212 -0.454 -0.319 -0.473 -0.245 -0.512 -0.348 P074 -0.468 -0.136 -0.623 -0.230 -0.488 -0.168 -0.698 -0.320 P007 0.158 -0.407 0.255 -0.368 0.158 -0.407 0.255 -0.368 P078 -0.434 -0.030 -0.433 -0.428 -0.434 -0.030 -0.433 -0.428 P159 -0.355 -0.427 -0.435 -0.386 -0.447 -0.522 -0.577 -0.476 P167 -0.356 -0.627 -0.508 -0.492 -0.479 -0.710 -0.612 -0.536 P189 0.226 -0.234 -0.307 -0.245 0.049 -0.198 -0.450 -0.251 P124 -0.719 -0.083 -0.479 -0.014 -0.719 -0.083 -0.479 -0.014 P116 -0.447 0.044 -0.724 -0.198 -0.447 0.044 -0.724 -0.198 P131 0.202 -0.151 -0.382 -0.386 0.202 -0.151 -0.382 -0.386 P101 -0.078 -0.599 -0.422 -0.485 -0.078 -0.599 -0.422 -0.485 main: 24 rats, 48 points: r = 0.552 (95% CI 0.318 to 0.723), p = 4.71e-05 loc alone: r = 0.470, p = 0.0206 frq alone: r = 0.736, p = 4.19e-05 fixed 200: 24 rats, 48 points: r = 0.577 (95% CI 0.350 to 0.740), p = 1.77e-05 loc alone: r = 0.519, p = 0.00933 frq alone: r = 0.722, p = 6.71e-05 Paper: r = 0.58 (tolerance 0.522 to 0.638), p = 0.000013
Result: reproduced. Over 24 rats and 48 values the first-half and second-half slope indices correlate with r = 0.552 (95% CI 0.318 to 0.723), p = 4.7e-05, within 0.522 to 0.638; with the code's fixed penalty r = 0.577, p = 1.8e-05.
ed4-bilateral-fof-impairs¶
Bilateral optogenetic silencing of the frontal orienting fields during the stimulus impairs performance of the task (three rats).
rows = []
opto_sessions = S[S.readable & (S.valid >= 100)].index.intersection(eps.sid.unique())
for i in opto_sessions:
tr = results[i]["trials"].reset_index(drop=True)
e = eps[(eps.sid == i) & eps.in_trial]
lab = pd.Series("control", index=tr.index, dtype=object)
for _, ep in e.iterrows():
lab[ep.trial] = f"{ep.sites}_{ep.window}"
tr = tr.assign(cond=lab.values)
tr = tr[tr.hit.notna()]
rows.append(tr.assign(subject=S.subject[i]))
ot = pd.concat(rows, ignore_index=True)
n_sessions = S.loc[opto_sessions].groupby("subject").size()
print("Opto sessions with >= 100 valid trials per rat:", n_sessions.to_dict())
tab = ot.groupby(["subject", "cond"]).hit.agg(["size", "mean"]).unstack(1)
print(tab["mean"].round(3).to_string()); print(tab["size"].fillna(0).astype(int).to_string())
top3 = n_sessions.sort_values(ascending=False).index[:3].tolist()
res_rows = []
for rat in n_sessions.index:
g = ot[ot.subject == rat]
ctl = g[g.cond == "control"].hit
for name, cond in (("full", g.cond == "bilateral_full"), ("any", g.cond.str.startswith("bilateral"))):
b = g[cond].hit
if len(b) == 0:
res_rows.append(dict(rat=rat, epochs=name, n_bilateral=0)); continue
table = [[int(b.sum()), int(len(b) - b.sum())], [int(ctl.sum()), int(len(ctl) - ctl.sum())]]
p = stats.fisher_exact(table)[1]
res_rows.append(dict(rat=rat, epochs=name, n_bilateral=len(b), acc_bilateral=b.mean(), n_control=len(ctl),
acc_control=ctl.mean(), difference=b.mean() - ctl.mean(), p=p))
bil = pd.DataFrame(res_rows)
# Was the laser connected on the trials of each label? (valid trials, per rat)
print("\nValid trials by label, opto_connected and opto_type:")
print(pd.crosstab([ot.subject, ot.cond], [ot.opto_connected, ot.opto_type]).to_string())
conn_rows = []
for rat in n_sessions.index:
g = ot[(ot.subject == rat) & (ot.opto_connected == 1)]
b, ctl = g[g.cond == "bilateral_full"].hit, g[g.cond == "control"].hit
p = stats.fisher_exact([[int(b.sum()), int(len(b) - b.sum())], [int(ctl.sum()), int(len(ctl) - ctl.sum())]])[1]
conn_rows.append(dict(rat=rat, n_bilateral=len(b), acc_bilateral=b.mean(), n_control=len(ctl), acc_control=ctl.mean(),
difference=b.mean() - ctl.mean(), p=p))
print("\nBilateral full-trial epochs against control, laser-connected trials only:")
print(pd.DataFrame(conn_rows).round(4).to_string(index=False))
print("\n", bil.round(4).to_string(index=False))
print("\nThree rats with the most opto sessions:", top3)
for name in ("full", "any"):
for reading, rr in (("three rats", top3), ("all rats", list(n_sessions.index))):
g = bil[(bil.epochs == name) & bil.rat.isin(rr)]
sig = ((g.difference < 0) & (g.p < 0.05)).sum()
print(f" {name} bilateral epochs, {reading}: lower with p < 0.05 in {sig} of {len(rr)}")
Opto sessions with >= 100 valid trials per rat: {'P116': 47, 'P124': 74, 'P127': 4, 'P131': 88}
cond bilateral_first bilateral_full bilateral_second control left_first left_full left_second right_first right_full right_second
subject
P116 0.693 0.639 0.683 0.733 0.699 0.636 0.687 0.678 0.611 0.655
P124 0.667 0.712 0.627 0.687 0.667 0.678 0.664 0.723 0.708 0.650
P127 0.500 0.603 0.250 0.634 0.600 0.594 0.500 0.333 0.629 0.375
P131 0.707 0.686 0.699 0.727 0.709 0.699 0.736 0.720 0.694 0.719
cond bilateral_first bilateral_full bilateral_second control left_first left_full left_second right_first right_full right_second
subject
P116 264 873 265 8462 246 943 278 233 895 252
P124 448 712 442 9072 445 674 432 423 664 452
P127 4 126 8 741 5 96 6 6 105 8
P131 1179 1988 1234 25308 1172 1889 1192 1142 2014 1239
Valid trials by label, opto_connected and opto_type:
opto_connected 0 1
opto_type First Half Full Trial Second Half First Half Full Trial Second Half
subject cond
P116 bilateral_first 0 0 0 264 0 0
bilateral_full 0 0 0 0 873 0
bilateral_second 0 0 0 0 0 265
control 2 4 0 1492 5440 1524
left_first 0 0 0 246 0 0
left_full 0 0 0 0 943 0
left_second 0 0 0 0 0 278
right_first 0 0 0 233 0 0
right_full 0 0 0 0 895 0
right_second 0 0 0 0 0 252
P124 bilateral_first 0 0 0 448 0 0
bilateral_full 0 0 0 0 712 0
bilateral_second 0 0 0 0 0 442
control 3 5 4 2568 3959 2533
left_first 0 0 0 445 0 0
left_full 0 0 0 0 674 0
left_second 0 0 0 0 0 432
right_first 0 0 0 423 0 0
right_full 0 0 0 0 664 0
right_second 0 0 0 0 0 452
P127 bilateral_first 0 0 0 4 0 0
bilateral_full 0 0 0 0 126 0
bilateral_second 0 0 0 0 0 8
control 0 0 0 39 656 46
left_first 0 0 0 5 0 0
left_full 0 0 0 0 96 0
left_second 0 0 0 0 0 6
right_first 0 0 0 6 0 0
right_full 0 0 0 0 105 0
right_second 0 0 0 0 0 8
P131 bilateral_first 0 0 0 1179 0 0
bilateral_full 0 0 0 0 1988 0
bilateral_second 0 0 0 0 0 1234
control 0 7 32 6907 11607 6755
left_first 0 0 0 1172 0 0
left_full 0 0 0 0 1889 0
left_second 0 0 0 0 0 1192
right_first 0 0 0 1142 0 0
right_full 0 0 0 0 2014 0
right_second 0 0 0 0 0 1239
Bilateral full-trial epochs against control, laser-connected trials only:
rat n_bilateral acc_bilateral n_control acc_control difference p
P116 873 0.6392 8456 0.7329 -0.0937 0.0000
P124 712 0.7121 9060 0.6865 0.0255 0.1655
P127 126 0.6032 741 0.6343 -0.0311 0.5495
P131 1988 0.6861 25269 0.7269 -0.0407 0.0001
rat epochs n_bilateral acc_bilateral n_control acc_control difference p
P116 full 873 0.6392 8462 0.7330 -0.0939 0.0000
P116 any 1402 0.6576 8462 0.7330 -0.0754 0.0000
P124 full 712 0.7121 9072 0.6868 0.0252 0.1656
P124 any 1602 0.6760 9072 0.6868 -0.0108 0.3975
P127 full 126 0.6032 741 0.6343 -0.0311 0.5495
P127 any 138 0.5797 741 0.6343 -0.0546 0.2503
P131 full 1988 0.6861 25308 0.7269 -0.0408 0.0001
P131 any 4401 0.6953 25308 0.7269 -0.0316 0.0000
Three rats with the most opto sessions: ['P131', 'P124', 'P116']
full bilateral epochs, three rats: lower with p < 0.05 in 2 of 3
full bilateral epochs, all rats: lower with p < 0.05 in 2 of 4
any bilateral epochs, three rats: lower with p < 0.05 in 2 of 3
any bilateral epochs, all rats: lower with p < 0.05 in 2 of 4
Result: partially reproduced. Four rats have optogenetic sessions, where the paper reports three. With bilateral full-trial epochs, accuracy is lower than on control trials of the same sessions with p < 0.05 in 2 of the 3 rats with the most opto sessions: P116 (0.639 against 0.733, p below 1e-4) and P131 (0.686 against 0.727, p = 0.0001). P124 is higher (0.712 against 0.687, p = 0.17), and P127 is lower but not significantly (4 sessions, 0.603 against 0.634, p = 0.55). Any bilateral epoch gives the same counts. Every trial with an epoch has the laser connected, and the few control trials without it change nothing: with laser-connected trials only, P116, P124 and P131 give the same accuracies and p values.
neural-link, units-per-rat¶
The neurons in the paper's published neural data can be linked to the dandiset's rats and sessions; they come from 199 sessions of 7 rats. The recordings yielded 3,495 putative single units, 2,047, 832, 258, 94 and 54 from five rats in FOF and 112 and 98 from two rats in mPFC.
t0 = time.time()
rat_numbers = ts.paper_mat("data/rat_for_each_neuron.mat")["rat_numbers"].astype(int)
n_cells = len(rat_numbers)
def read_cell(n):
m = ts.paper_mat(f"data/neural_data/cell_{n}.mat", variable_names=["behav", "rat_name", "brain_region"])
b = m["behav"]
return dict(cell=n, rat=str(m["rat_name"]), region=str(m["brain_region"]), n_trials=int(b.nTrials),
sym=ts.trial_symbols(list(str(b.task)), np.atleast_1d(b.gdir), np.atleast_1d(b.gfreq), np.atleast_1d(b.choice)))
with ThreadPoolExecutor(16) as ex:
cells = pd.DataFrame(list(ex.map(read_cell, range(1, n_cells + 1))))
print(f"read {len(cells)} cell files in {time.time() - t0:.0f} s")
print("Cells per rat and region:")
print(pd.crosstab(cells.rat, cells.region).to_string())
print("Units in total:", len(cells), "| FOF:", int((cells.region == "FOF").sum()), "| mPFC:", int((cells.region == "mPFC").sum()))
print("rat_for_each_neuron agrees with the rat named in each cell file:",
bool((("P" + pd.Series(rat_numbers).astype(str).str.zfill(3)).values == cells.rat.values).all()))
read 3495 cell files in 25 s Cells per rat and region: region FOF mPFC rat P049 54 0 P055 258 0 P059 832 0 P095 94 0 P100 2047 0 P101 0 112 P102 0 98 Units in total: 3495 | FOF: 3285 | mPFC: 210 rat_for_each_neuron agrees with the rat named in each cell file: True
# Each rat's valid trials, session by session, as one string; a cell matches where its trials appear.
sym_sessions = {}
for i, r in enumerate(results):
if r["trials"] is None or S.subject[i] not in ts.PAPER_EPHYS_RATS:
continue
v = r["trials"][r["trials"].hit.notna()]
sym_sessions.setdefault(S.subject[i], []).append((S.session[i], ts.trial_symbols(v.task, v.gdir, v.gfreq, v.choice)))
matches = []
for _, c in cells.iterrows():
found = [sid for sid, s in sym_sessions.get(c.rat, []) if c.sym and c.sym in s]
matches.append(found)
cells["n_match"] = [len(f) for f in matches]
cells["session"] = [f[0] if len(f) == 1 else None for f in matches]
print("Rats of the neural data in the dandiset:", sorted(set(cells.rat) & set(assets.subject)))
print("Cells matching exactly one session:", int((cells.n_match == 1).sum()), "of", len(cells),
f"({(cells.n_match == 1).mean():.4f}); none: {int((cells.n_match == 0).sum())}; more than one: {int((cells.n_match > 1).sum())}")
print("Distinct matched sessions:", cells.session.nunique(), "| per rat:", cells.groupby("rat").session.nunique().to_dict())
print("Cells without a match, per rat:", cells[cells.n_match == 0].groupby("rat").size().to_dict(),
"| their trial counts:", cells[cells.n_match == 0].n_trials.describe()[["min", "50%", "max"]].to_dict())
print("Matched sessions by protocol:", pd.Series([s.split("-")[0] for s in cells.session.dropna().unique()]).value_counts().to_dict())
print("Paper: 199 sessions (tolerance 179.1 to 218.9)")
# Variant: the cell's trials as an ordered subsequence of one session's valid trials (gaps allowed).
gap = [m if m else ts.gap_match(c.sym, sym_sessions.get(c.rat, [])) for m, (_, c) in zip(matches, cells.iterrows())]
cells["n_gap"] = [len(g) for g in gap]
cells["session_gap"] = [g[0] if len(g) == 1 else None for g in gap]
print("\nVariant, gaps allowed: cells matching exactly one session", int((cells.n_gap == 1).sum()), "of", len(cells),
f"({(cells.n_gap == 1).mean():.4f}); none: {int((cells.n_gap == 0).sum())}; more than one: {int((cells.n_gap > 1).sum())}")
print("Distinct sessions:", cells.session_gap.nunique(), "| per rat:", cells.groupby("rat").session_gap.nunique().to_dict())
unm = cells[cells.n_gap == 0]
print("Cells still unmatched, per rat:", unm.groupby("rat").size().to_dict())
# Variant: every TaskSwitch protocol of the recorded rats, adding TaskSwitch4repeat and others.
sym_all = {k: list(v) for k, v in sym_sessions.items()}
for i, r in enumerate(other_results):
if r["trials"] is None or other_assets.subject[i] not in ts.PAPER_EPHYS_RATS:
continue
v = r["trials"][r["trials"].hit.notna()]
sym_all[other_assets.subject[i]].append((other_assets.session[i], ts.trial_symbols(v.task, v.gdir, v.gfreq, v.choice)))
print("\nSessions searched per rat, three task protocols:", {k: len(v) for k, v in sym_sessions.items()},
"| every TaskSwitch protocol:", {k: len(v) for k, v in sym_all.items()})
m_all = [[sid for sid, s in sym_all.get(c.rat, []) if c.sym and c.sym in s] for _, c in cells.iterrows()]
g_all = [m if m else ts.gap_match(c.sym, sym_all.get(c.rat, [])) for m, (_, c) in zip(m_all, cells.iterrows())]
cells["n_match_all"] = [len(f) for f in m_all]
cells["n_gap_all"] = [len(f) for f in g_all]
cells["session_all"] = [f[0] if len(f) == 1 else None for f in m_all]
cells["session_gap_all"] = [f[0] if len(f) == 1 else None for f in g_all]
for col, name in (("n_match_all", "contiguous"), ("n_gap_all", "gaps allowed")):
sess_col = "session_all" if col == "n_match_all" else "session_gap_all"
print(f"Every TaskSwitch protocol, {name}: cells matching exactly one session {int((cells[col] == 1).sum())} "
f"({(cells[col] == 1).mean():.4f}); none {int((cells[col] == 0).sum())}; more than one {int((cells[col] > 1).sum())}; "
f"distinct sessions {cells[sess_col].nunique()}")
print("Sessions matched (gaps allowed) by protocol:", pd.Series([s.split("-")[0] for s in cells.session_gap_all.dropna().unique()]).value_counts().to_dict())
# Per rat: cells, matches, sessions and the dates of the matched sessions (every TaskSwitch protocol).
cells["date_all"] = cells.session_gap_all.str[-7:-1]
per_rat_cells = cells.groupby("rat").agg(
cells=("cell", "size"), contiguous=("n_match_all", lambda x: int((x == 1).sum())),
gaps_allowed=("n_gap_all", lambda x: int((x == 1).sum())), unmatched=("n_gap_all", lambda x: int((x == 0).sum())),
sessions_contiguous=("session_all", "nunique"), sessions_gaps=("session_gap_all", "nunique"),
first_date=("date_all", "min"), last_date=("date_all", "max"), trials_median=("n_trials", "median"))
per_rat_cells["unmatched_trials_median"] = cells[cells.n_gap_all == 0].groupby("rat").n_trials.median()
print(per_rat_cells.to_string())
# The cells that match no session: how much of each falls inside one dandiset session. Their
# trials are cut into windows of 8, and the share of windows found in the best session is given.
rows = []
for _, c in cells[cells.n_gap_all == 0].iterrows():
w = [c.sym[k:k + 8] for k in range(0, len(c.sym) - 7, 8)]
best = max((np.mean([x in s for x in w]), sid) for sid, s in sym_all[c.rat])
rows.append(dict(rat=c.rat, trials=c.n_trials, share=best[0], session=best[1]))
unm_all = pd.DataFrame(rows)
print(f"\nUnmatched cells: {len(unm_all)}; share of their 8-trial windows found in their best session: median "
f"{unm_all.share.median():.2f}, range {unm_all.share.min():.2f} to {unm_all.share.max():.2f}; "
f"best sessions not among those matched: {len(set(unm_all.session) - set(cells.session_gap_all.dropna()))} "
f"of {unm_all.session.nunique()}")
Rats of the neural data in the dandiset: ['P049', 'P055', 'P059', 'P095', 'P100', 'P101', 'P102']
Cells matching exactly one session: 2371 of 3495 (0.6784); none: 1124; more than one: 0
Distinct matched sessions: 110 | per rat: {'P049': 6, 'P055': 10, 'P059': 10, 'P095': 36, 'P100': 71, 'P101': 23, 'P102': 21}
Cells without a match, per rat: {'P049': 2, 'P059': 168, 'P095': 34, 'P100': 835, 'P101': 36, 'P102': 49} | their trial counts: {'min': 81.0, '50%': 409.0, 'max': 798.0}
Matched sessions by protocol: {'TaskSwitch4': 98, 'TaskSwitch6': 12}
Paper: 199 sessions (tolerance 179.1 to 218.9)
Variant, gaps allowed: cells matching exactly one session 3269 of 3495 (0.9353); none: 226; more than one: 0
Distinct sessions: 112 | per rat: {'P049': 6, 'P055': 10, 'P059': 10, 'P095': 39, 'P100': 72, 'P101': 28, 'P102': 26}
Cells still unmatched, per rat: {'P059': 11, 'P095': 8, 'P100': 189, 'P101': 2, 'P102': 16}
Sessions searched per rat, three task protocols: {'P049': 747, 'P055': 814, 'P059': 639, 'P095': 607, 'P100': 719, 'P101': 579, 'P102': 629} | every TaskSwitch protocol: {'P049': 747, 'P055': 814, 'P059': 639, 'P095': 607, 'P100': 732, 'P101': 579, 'P102': 629}
Every TaskSwitch protocol, contiguous: cells matching exactly one session 2425 (0.6938); none 1070; more than one 0; distinct sessions 113
Every TaskSwitch protocol, gaps allowed: cells matching exactly one session 3349 (0.9582); none 146; more than one 0; distinct sessions 115
Sessions matched (gaps allowed) by protocol: {'TaskSwitch4': 100, 'TaskSwitch6': 12, 'TaskSwitch4repeat': 3}
cells contiguous gaps_allowed unmatched sessions_contiguous sessions_gaps first_date last_date trials_median unmatched_trials_median
rat
P049 54 52 54 0 6 6 180402 180508 244.0 NaN
P055 258 258 258 0 10 10 180205 180516 427.0 NaN
P059 832 664 821 11 10 10 180402 180508 319.0 303.0
P095 94 60 86 8 36 39 180810 190204 480.5 394.0
P100 2047 1266 1938 109 74 75 180713 190716 396.0 421.0
P101 112 76 110 2 23 28 180823 190123 286.5 424.5
P102 98 49 82 16 21 26 180823 190204 587.0 540.5
Unmatched cells: 146; share of their 8-trial windows found in their best session: median 0.43, range 0.02 to 0.94; best sessions not among those matched: 3 of 64
Result: neural-link not reproduced. The seven rats are in the dandiset and every cell names one of them, but only 2,371 of 3,495 cells (0.678) have trials equal to a contiguous run of one session's valid trials, from 110 distinct sessions against the paper's 199 (tolerance 179 to 219). With gaps allowed (the variant), 3,269 cells (0.935) match exactly one session and none matches more than one, so the neural files leave out some of a session's valid trials and the rule chosen here is stricter than the data need; but these cells still come from only 112 dandiset sessions, and 226 cells, 189 of them P100's, match none.
Searching every TaskSwitch protocol of the recorded rats (adding P100's 13 TaskSwitch4repeat sessions) raises the contiguous matches to 2,425 (0.694) and the matches with gaps to 3,349 (0.958), from 113 and 115 sessions; 3 of the new sessions are TaskSwitch4repeat. The 146 cells still unmatched share a median of 0.43 of their 8-trial windows with their best session, and all but 3 of those best sessions are among the 115 already matched, so they do not point to missing sessions. Per rat the matched sessions are 6, 10, 10, 39, 75, 28 and 26, recorded from February 2018 to July 2019. Under no search does the number of sessions approach 199.
Result: units-per-rat reproduced. The cell files hold 3,495 units, 3,285 in FOF (2,047, 832, 258, 94 and 54 in P100, P059, P055, P095 and P049) and 210 in mPFC (112 in P101 and 98 in P102), as the Methods' per-rat counts say. The Methods' "3,495 ... recorded from 5 rats in FOF" counts the mPFC units as well.
fig3h-neural-uncorrelated, choice-axes-parallel¶
Neural slope indices vary across rats and features, with no apparent correlation between the location and frequency indices across the seven recorded rats. Choice axes estimated separately in the two contexts are essentially parallel (average angle 1.6 degrees).
Result: reproduced. The neural location and frequency slope indices from the paper's fitted kernels are not significantly correlated across the seven rats (r = -0.177, 95% CI -0.821 to 0.665, p = 0.704), and they vary: location -0.221 to 1.728, frequency 0.465 to 1.500. The Methods say the differential pulse responses "only displayed positive (or zero) slope indices"; three location indices are negative here, P049 -0.221, P101 -0.045 and P095 -0.007. The claim's criterion does not test this.
ang = pd.DataFrame({r: ts.choice_axis_angle(ker_ctx["data_all_trials_loc"], ker_ctx["data_all_trials_freq"], rat_numbers, int(r[1:]))
for r in ts.PAPER_EPHYS_RATS}).T.astype(float)
ang["full_folded"] = np.minimum(ang.full, 180 - ang.full)
print(ang[["projected", "full", "full_folded", "var3", "units"]].round(3).to_string())
print(f"\nMean over 7 rats: projected on the first 3 PCs and folded (the code's plotted value, main) {ang.projected.mean():.2f} degrees "
f"(range {ang.projected.min():.2f} to {ang.projected.max():.2f}); full-space angle {ang.full.mean():.2f}, folded {ang.full_folded.mean():.2f}; "
f"paper 1.6 (tolerance 1.44 to 1.76)")
print("Rats with a projected angle within 1.44 to 1.76:", int(ang.projected.between(1.44, 1.76).sum()),
"| below 1.76:", int((ang.projected < 1.76).sum()))
projected full full_folded var3 units P049 5.532 24.750 24.750 0.904 54.0 P055 7.914 36.254 36.254 0.846 258.0 P059 4.488 27.080 27.080 0.899 832.0 P095 1.819 13.388 13.388 0.939 94.0 P100 1.121 19.815 19.815 0.884 2047.0 P101 4.245 32.892 32.892 0.891 112.0 P102 8.004 128.192 51.808 0.845 98.0 Mean over 7 rats: projected on the first 3 PCs and folded (the code's plotted value, main) 4.73 degrees (range 1.12 to 8.00); full-space angle 40.34, folded 29.43; paper 1.6 (tolerance 1.44 to 1.76) Rats with a projected angle within 1.44 to 1.76: 0 | below 1.76: 1
Result: not reproduced. ExtDataFigure5ag.m computes two angles and plots the second: the two context choice axes projected on the first 3 principal components of both contexts' choice kernels (which explain 0.85 to 0.94 of their variance), folded to at most 90 degrees. On the paper's own fitted per-context kernels these angles are 1.12 to 8.00 degrees (P100 1.12, P095 1.82, P101 4.25, P059 4.49, P049 5.53, P055 7.91, P102 8.00), mean 4.73, against the paper's 1.6 (tolerance 1.44 to 1.76); no rat lies in the tolerance and one lies below it. The angle in the full space, which the code computes first and then overwrites, is 13.4 to 128.2 degrees, mean 40.3 (29.4 folded). The axes are close to parallel in the 3-dimensional space but not as close as reported. The paper's test of each angle against 0 ("not different from 0 in 6 of 7 rats") is the permutation test the Methods describe under "Estimating the choice axis": the context labels of the trials are shuffled and the axes recomputed. It is not in the repository's code and needs the per-neuron kernels refitted for every shuffle, so it was not run.
fig5d-behaviour-neural¶
The behavioural slope index and the neural slope index are anti-correlated across rats and features (r = -0.73, P < 0.01), and for location (r = -0.71) and frequency (r = -0.71) alone (P < 0.1).
half2 = slice(13, 26) # bins 14 to 26, the second half of the stimulus
beh = pd.DataFrame({r: dict(loc=ts.slopes(fits[r]["main"], half2)[0], frq=ts.slopes(fits[r]["main"], half2)[1],
loc_code=ts.slopes(fits[r]["code_sel_200"], half2)[0], frq_code=ts.slopes(fits[r]["code_sel_200"], half2)[1],
trials=fits[r]["n"], trials_code=fits[r]["code_sel_n"]) for r in ts.PAPER_EPHYS_RATS}).T
saved_e = ts.paper_mat("data/logistic_regression_results_ephys_rats.mat")
beh["loc_paper_saved"], beh["frq_paper_saved"] = saved_e["location"], saved_e["frequency"]
print(pd.concat([beh, neural.add_prefix("neural_")], axis=1).round(3).to_string())
def corr_report(bl, bf, name):
r, p, lo, hi = pearson_ci(np.r_[bl, bf], np.r_[neural["loc"], neural["frq"]])
rl, pl = stats.pearsonr(bl, neural["loc"]); rf, pf = stats.pearsonr(bf, neural["frq"])
print(f"{name}: both r = {r:.3f} (95% CI {lo:.3f} to {hi:.3f}), p = {p:.3g}; location r = {rl:.3f} (p = {pl:.3g}); "
f"frequency r = {rf:.3f} (p = {pf:.3g})")
corr_report(beh["loc"], beh["frq"], "main (Methods selection, cross-validated penalty)")
corr_report(beh["loc_code"], beh["frq_code"], "variant (code's 0.65/0.6 selection, penalty 200)")
corr_report(beh["loc_paper_saved"], beh["frq_paper_saved"], "reference (paper's saved behavioural indices)")
beh["loc_code_cv"] = [ts.slopes(fits[r]["code_sel_cv"], half2)[0] for r in beh.index]
beh["frq_code_cv"] = [ts.slopes(fits[r]["code_sel_cv"], half2)[1] for r in beh.index]
corr_report(beh["loc_code_cv"], beh["frq_code_cv"], "variant (code's 0.65/0.6 selection, cross-validated penalty)")
mixed = beh.index.isin(["P049", "P059"])
corr_report(beh["loc"].where(~mixed, beh["loc_code_cv"]), beh["frq"].where(~mixed, beh["frq_code_cv"]),
"variant (Methods selection, P049 and P059 from the code's selection)")
loo = []
for r in beh.index:
k = beh.index != r
loo.append(dict(left_out=r, r=stats.pearsonr(np.r_[beh["loc"][k], beh["frq"][k]], np.r_[neural["loc"][k], neural["frq"][k]])[0]))
loo = pd.DataFrame(loo).set_index("left_out")
print("Main, leaving out one rat (12 points):", loo.r.round(3).to_dict(), "| within -0.803 to -0.657:",
int(loo.r.between(-0.803, -0.657).sum()), "of 7")
print("Paper: both -0.73 (tolerance -0.803 to -0.657), P < 0.01; location -0.71, frequency -0.71, P < 0.1")
fig, ax = plt.subplots(figsize=(4, 3.5))
ax.scatter(beh["loc"], neural["loc"], c="k", label="location"); ax.scatter(beh["frq"], neural["frq"], c="b", label="frequency")
ax.set_xlabel("behavioural slope index (dandiset)"); ax.set_ylabel("neural slope index (paper kernels)"); ax.legend(); plt.show()
print("Plotted ranges: behaviour", round(min(beh['loc'].min(), beh['frq'].min()), 3), "to", round(max(beh['loc'].max(), beh['frq'].max()), 3),
"| neural", round(min(neural['loc'].min(), neural['frq'].min()), 3), "to", round(max(neural['loc'].max(), neural['frq'].max()), 3))
loc frq loc_code frq_code trials trials_code loc_paper_saved frq_paper_saved neural_loc neural_frq neural_units
P049 0.131 -0.407 -0.432 -0.147 23014.0 135001.0 -0.415 -0.138 -0.221 0.610 54
P055 -0.966 -0.803 -0.831 -0.731 480165.0 498391.0 -0.833 -0.727 0.649 0.596 258
P059 1.292 -1.136 -0.613 -1.013 25847.0 162317.0 -0.608 -1.017 0.373 1.500 832
P095 -0.454 -0.788 0.540 -0.641 197061.0 269734.0 0.535 -0.641 -0.007 0.465 94
P100 -1.398 -1.143 -1.425 -0.960 283348.0 318846.0 -1.418 -0.961 1.312 1.030 2047
P101 0.085 -1.054 -0.129 -1.040 100290.0 172819.0 -0.128 -1.045 -0.045 1.482 112
P102 -1.328 -0.655 -0.994 -0.666 290572.0 328818.0 -0.989 -0.665 1.728 0.531 98
main (Methods selection, cross-validated penalty): both r = -0.697 (95% CI -0.896 to -0.264), p = 0.0056; location r = -0.675 (p = 0.0963); frequency r = -0.772 (p = 0.042)
variant (code's 0.65/0.6 selection, penalty 200): both r = -0.765 (95% CI -0.921 to -0.394), p = 0.00144; location r = -0.766 (p = 0.0448); frequency r = -0.729 (p = 0.0629)
reference (paper's saved behavioural indices): both r = -0.768 (95% CI -0.923 to -0.402), p = 0.00132; location r = -0.769 (p = 0.0432); frequency r = -0.730 (p = 0.0625)
variant (code's 0.65/0.6 selection, cross-validated penalty): both r = -0.830 (95% CI -0.945 to -0.535), p = 0.000239; location r = -0.885 (p = 0.00804); frequency r = -0.719 (p = 0.0684)
variant (Methods selection, P049 and P059 from the code's selection): both r = -0.834 (95% CI -0.946 to -0.543), p = 0.000213; location r = -0.924 (p = 0.00296); frequency r = -0.625 (p = 0.134)
Main, leaving out one rat (12 points): {'P049': -0.657, 'P055': -0.717, 'P059': -0.895, 'P095': -0.739, 'P100': -0.655, 'P101': -0.66, 'P102': -0.67} | within -0.803 to -0.657: 4 of 7
Paper: both -0.73 (tolerance -0.803 to -0.657), P < 0.01; location -0.71, frequency -0.71, P < 0.1
Plotted ranges: behaviour -1.398 to 1.292 | neural -0.221 to 1.728
Result: reproduced. Between the behavioural slope index from the dandiset (Methods selection, second half of the stimulus) and the neural slope index from the paper's kernels, r = -0.697 over 14 values (95% CI -0.896 to -0.264), p = 0.0056, within -0.803 to -0.657. For location alone r = -0.675 (p = 0.096) and for frequency alone r = -0.772 (p = 0.042). With the code's selection and penalty r = -0.765 (p = 0.0014), and the paper's own saved behavioural indices give -0.768.
The result is close to the edge of the tolerance and rests on two rats with few trials. P049 and P059 have 23,014 and 25,847 trials under the Methods selection against 135,001 and 162,317 under the code's, and their location indices change sign between the two (0.131 against -0.432; 1.292 against -0.613). Taking those two rats from the code's selection gives r = -0.834, and the code's selection with the cross-validated penalty r = -0.830, both more negative than the tolerance. Leaving out one rat at a time gives -0.655 to -0.895, inside the tolerance for 4 of 7 rats. The five 14-point correlations printed above are all negative with p < 0.01.
print(f"total run time {time.time() - started:.0f} s")
total run time 413 s
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Data | Outcome |
|---|---|---|
| session-link | with inputs | reproduced |
| subjects | dandiset | partially reproduced |
| trials-per-session | dandiset | reproduced |
| blocks-per-session | dandiset | reproduced |
| session-duration | dandiset | reproduced |
| block-min-30 | dandiset | reproduced |
| fig1c-20-rats | dandiset | reproduced |
| fig1c-context-selection | dandiset | reproduced |
| ed1e-switch-four-trials | dandiset | reproduced |
| ed2g-learning-2-5-months | dandiset | reproduced |
| fig4-converging | dandiset | reproduced |
| fig4-loc-frq-uncorrelated | dandiset | reproduced |
| ed3g-halfsplit | dandiset | reproduced |
| ed4-bilateral-fof-impairs | dandiset | partially reproduced |
| neural-link | with inputs | not reproduced |
| units-per-rat | elsewhere | reproduced |
| fig3h-neural-uncorrelated | elsewhere | reproduced |
| choice-axes-parallel | elsewhere | not reproduced |
| fig5d-behaviour-neural | with inputs | reproduced |
| fig1e-choice-axis | elsewhere | blocked: methods underspecified |
| fig1f-irrelevant-encoded | elsewhere | blocked: not runnable |
| fig3a-trained-rnn-svm | elsewhere | blocked: not runnable |
| fig5c-engineered-rnn | elsewhere | blocked: not runnable |
How the analysis differs from the paper's. The session selection is the Methods' (the code's
extra good flag is not in the dandiset, and never removes a session that passes the other
criteria). The feature selection index uses the stored generative log rate ratios; the relative weights
equal the paper's per-session values to 1e-6 in 0.898 of matched sessions, 861 more agree once
capped at 1 as the paper's table does, 627 are missing on both sides, and 60 differ with both
values present. The behavioural kernels use the code's smoothness penalty, read
as the Methods' "ridge regression", with its weight chosen by 3-fold cross-validation over 20,
200 and 2000 (grid and folds chosen here); the code's fixed 200 is shown as a variant and gives
slope indices within 0.018 of those from the paper's saved weights. For the recorded rats the
Methods selection keeps far fewer trials of P049 and P059 than the code's looser one, which is
shown as a variant. Several numbers in the Methods (trials and blocks per session, session
length) do not say which sessions or trials they count; the readings are listed in the claims
file. The neural claims use the paper's fitted per-neuron kernels from its GitHub repository, not
a refit from spikes. The choice-axis angle is the one ExtDataFigure5ag.m plots, in the space of
the first 3 principal components of both contexts' choice kernels. The mTDR and RNN analyses were not run.
What was found about the dandiset. It holds the behaviour of all 26 rats and links to the paper's published per-session tables: 15,277 of 15,278 paper sessions with valid trials are in it, with equal valid trial counts and fraction correct, and 6 dandiset task sessions are not in the tables. In the 14,650 matched sessions with a per-trial record in the paper's tables every valid trial agrees on context, side, strengths, choice and outcome, and 0.99996 of 50 ms pulse bins agree; 627 have no per-trial record there to compare. The paper's trial selection can be rebuilt from the dandiset alone, adding 6 sessions. It does not hold the electrophysiology, which is only in the paper's GitHub repository, as MATLAB files per cell with no session identifier; most cells can be placed in a dandiset session only by matching their trial sequences. The optogenetic epochs state 473 nm where the paper says 450 nm, and are present in four rats where the paper reports three, without saying which three were analysed. Dates of birth give ages of 3 to 39 months at the analysed sessions, where the Methods say 6 to 24; they put 25 rats at 1.6 to 2.7 months at their first session, and P189's date of birth looks a year early. Thirteen task sessions lack the result codes or a pulse column, and one paper session with 27 valid trials is not in the dandiset. The cells of the paper's neural files match 110 dandiset sessions as contiguous runs of trials and 112 with gaps allowed (113 and 115 when every TaskSwitch protocol is searched), not the 199 the paper reports, and 146 to 226 match none. The choice-axis angles computed from the paper's own kernels by its own code average 4.7 degrees, not the 1.6 it reports.