Reproducing Hulsey et al. (2024) from dandiset 000678¶
Paper. Hulsey D, Zumwalt K, Mazzucato L, McCormick DA, Jaramillo S (2024). Decision-making dynamics are predicted by arousal and uninstructed movements. Cell Reports 43(1), 113709. doi:10.1016/j.celrep.2024.113709. The claims file was written from the accepted author manuscript in PubMed Central (PMC11016285) with its supplement; this attempt read the same text and the paper's analysis code (Zenodo 10.5281/zenodo.10306018, and github.com/sjara/uobrainflex at commit ef48985192e4af332c64f9d2e70ddcfc78659040).
Dandiset. DANDI:000678, version 0.231004.2146, streamed from the archive. This is the version the paper cites (Data and code availability, and the Key resources table). A later version, 0.261020.0029, has the same 391 assets with the same asset identifiers and sizes; it differs only in the archive metadata (keywords, affiliations, ORCIDs), so the choice does not change the data.
Claims. claims/000678/10.1016-j.celrep.2024.113709.yaml (29 claims). The GLM-HMM and the analyses after it follow the paper's code, reimplemented in hulsey.py; where the claims file leaves a detail open, the passage used is cited in the section or in the helper's docstring.
| Claim | Tested here |
|---|---|
| cohort-inclusion | tested |
| training-duration | blocked: data missing (no sessions before stage 5) |
| distractor-response-rate | tested |
| figs2-n-states | tested |
| fig1g-example-mouse-prediction | tested |
| fig2a-state-types | tested |
| fig2b-optimal-occupancy | tested |
| fig2c-dwell-times | tested |
| fig2d-prediction-accuracy | tested |
| fig3c-inverted-u-pupil | tested |
| fig3d-optimal-pupil-modality | tested |
| fig3e-u-disengaged-pupil | tested |
| figs4-r2-individual-movement | tested |
| fig3f-movement-index-inverted-u | tested |
| fig3f-movement-index-u-disengaged | tested |
| figs5-movement-index-correlations | tested |
| fig3-distance-from-optimal-pupil | tested |
| fig3-raw-face-energy-trend | tested |
| fig3-raw-locomotion-no-difference | tested |
| fig4a-lower-variability | tested |
| fig4b-variability-at-transitions | tested |
| fig4d-pupil-movement-decoupling | tested |
| figs6-individual-measure-decoupling | tested |
| fig5c-state-classification | tested (fewer shuffle runs) |
| figs7-pupil-best-raw-measure | tested (fewer shuffle runs) |
| fig5c-movement-index-best | tested (fewer shuffle runs) |
| figs8a-regression-weights | tested |
| figs8b-pupil-movement-interaction | tested |
| figs1-iti-independence | tested |
The notebook prints its own run time at the end of the decoding section. On 4 cores it took 15.7 minutes in the run this text was checked against: 2.0 to stream the 391 files, 9.1 for the cross-validated GLM-HMM fits and 4.1 for the state decoders. A run on the workflow's runner took 2.3, 12.9 and 5.7 minutes for these stages. Random seeds are fixed, and each claim's section has its own generator.
import os
for _v in ["OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"]:
os.environ[_v] = "1" # many small fits run in parallel processes; one linear-algebra thread each
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import hulsey
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
assets = hulsey.list_assets()
print(f"dandiset {hulsey.DANDISET} version {hulsey.VERSION}: {len(assets)} assets, "
f"{assets['size'].sum() / 1e9:.2f} GB")
trials_all, alltrials, sessions = hulsey.load_all(assets)
print(f"streamed {len(sessions)} files in {(time.time() - started) / 60:.1f} min")
dandiset 000678 version 0.231004.2146: 391 assets, 8.11 GB
streamed 391 files in 2.0 min
Data check¶
Each file holds one session: a trials table, the measures of the face camera (online and post hoc pupil, face and whisker motion energy, sampled at about 29.5 Hz), wheel speed (about 100 Hz), and lick and reward times. The checks of the metadata against the values are stated under each output.
meta = hulsey.version_metadata()
print("Archive metadata: name:", meta["name"])
print(" numberOfSubjects:", meta.get("assetsSummary", {}).get("numberOfSubjects"),
"| numberOfFiles:", meta.get("assetsSummary", {}).get("numberOfFiles"))
print("\nSessions by subject, task (session_type), training stage:")
by = sessions.assign(task=sessions.session_type).groupby(["subject", "task"]).agg(
sessions=("path", "size"), S5=("stage", lambda x: int((x == "S5").sum())),
S6=("stage", lambda x: int((x == "S6").sum())),
no_post_hoc_pupil=("pupil_raw_median", lambda x: int(x.isna().sum())))
print(by.to_string())
print("\nTask per subject from session_type agrees with target_modality of the trials (0 auditory, 1 visual):",
bool((sessions.groupby("subject").modality.nunique() == 1).all()
and ((sessions.session_type == "visual") == (sessions.modality == 1)).all()))
print("Training stages present:", sorted(sessions.stage.unique()))
print("Rigs:", sessions.rig.value_counts().to_dict())
print("\nTrials in the files:", len(alltrials), "| target trials:", len(trials_all),
"| distractor trials:", int(sessions.n_distractor.sum()))
print("Outcome codes by stimulus type (0 target, 1 distractor):")
print(pd.crosstab(alltrials.stimulus_type, alltrials.outcome).to_string())
print("\nStimulus values (code's inpt) per task:")
print(trials_all.assign(task=np.where(trials_all.target_modality == 0, "auditory", "visual"))
.groupby("task").inpt.apply(lambda x: sorted(np.round(x.unique(), 3))).to_string())
print("Visual difficulty codes present:", sorted(alltrials.visual_stim_difficulty.dropna().unique()),
"| auditory difficulty values:", sorted(np.round(alltrials.auditory_stim_difficulty.dropna().unique(), 2)))
print("Trials table columns:", list(alltrials.columns))
has_ang = alltrials[alltrials.target_modality == 1].groupby("path").visual_gabor_angle.apply(lambda x: x.notna().any())
print(f"Visual-task sessions: {len(has_ang)}, of which {int(has_ang.sum())} have a visual_gabor_angle column")
vsets = trials_all[trials_all.target_modality == 1].groupby("path").inpt.apply(lambda x: tuple(np.round(np.unique(np.abs(x)), 2)))
print("Visual stimulus values |inpt| per session (get_stim_vals), sessions per set:", vsets.value_counts().to_dict())
print("Stimulus duration (ms):", alltrials.stimulus_duration.unique())
iti = alltrials.groupby("path").start_time.diff().dropna()
print(f"Start-to-start interval (s): median {iti.median():.2f}, 5-95% {iti.quantile(.05):.2f} to {iti.quantile(.95):.2f}")
print("\nSampling: face camera", sessions.face_rate.round(2).value_counts().to_dict(),
"Hz; wheel", sessions.run_rate.round(1).value_counts().to_dict(), "Hz")
print("Raw running_speed per session (unit stated 'cm/s'): median of session medians "
f"{sessions.run_raw_median.median():.3f}, median of session 99th percentiles {sessions.run_raw_p99.median():.3f}, "
f"largest 99th percentile {sessions.run_raw_p99.max():.3f}")
print("Raw post hoc pupil diameter (pixels): session medians from "
f"{sessions.pupil_raw_median.min():.1f} to {sessions.pupil_raw_median.max():.1f}")
cols = ["post_hoc_pupil_diameter", "face_energy", "running_speed", "post_hoc_pupil_std10", "face_std_10", "running_std10"]
print("\nMeasures at stimulus onset after the code's filtering and normalisation (target trials):")
print(trials_all[cols].describe().T[["count", "mean", "min", "50%", "max"]].round(3).to_string())
print("\nDuplicate recordings: identical face motion energy arrays", int(sessions.fp_face.duplicated(keep=False).sum()),
"| identical trial tables (stimulus times and outcomes)", int(sessions.fp_trials.duplicated(keep=False).sum()))
Archive metadata: name: Pupil and movement measurements during mouse auditory and visual discrimination
numberOfSubjects: 13 | numberOfFiles: 391
Sessions by subject, task (session_type), training stage:
sessions S5 S6 no_post_hoc_pupil
subject task
BW041 visual 25 15 10 0
BW045 auditory 84 10 74 1
BW046 visual 34 18 16 14
BW048 visual 35 20 15 0
BW051 visual 35 25 10 4
BW053 auditory 32 23 9 3
BW058 visual 12 12 0 2
BW059 auditory 25 18 7 5
BW061 auditory 23 23 0 1
BW062 auditory 17 17 0 4
BW064 auditory 21 12 9 2
BW067 auditory 23 23 0 0
BW071 auditory 25 25 0 2
Task per subject from session_type agrees with target_modality of the trials (0 auditory, 1 visual): True
Training stages present: ['S5', 'S6']
Rigs: {'Box 2': 157, 'Box1': 63, 'Box 4': 61, 'WF2': 59, 'Box 3': 35, 'Box 1': 13, 'Boox1': 3}
Trials in the files: 214985 | target trials: 167948 | distractor trials: 47037
Outcome codes by stimulus type (0 target, 1 distractor):
outcome 1.0 2.0 4.0 8.0 16.0
stimulus_type
0.0 109907 26879 0 0 31162
1.0 0 0 2446 44591 0
Stimulus values (code's inpt) per task:
task
auditory [-0.97, -0.625, -0.43, 0.43, 0.625, 0.97]
visual [-1.0, -0.8, -0.6, -0.2, 0.2, 0.6, 0.8, 1.0]
Visual difficulty codes present: [np.float64(0.0), np.float64(1.0), np.float64(2.0)] | auditory difficulty values: [np.float64(0.62), np.float64(0.75), np.float64(0.98)]
Trials table columns: ['start_time', 'stop_time', 'target_modality', 'cue_time', 'cue_ID', 'cue_duration', 'stimulus_type', 'stimulus_time', 'stimulus_duration', 'target_port', 'auditory_stim_id', 'auditory_stim_band', 'auditory_stim_difficulty', 'visual_stim_id', 'visual_stim_oreintation', 'visual_stim_difficulty', 'outcome', 'choice', 'aud_stim', 'vis_stim', 'stage', 'path', 'visual_gabor_angle']
Visual-task sessions: 141, of which 123 have a visual_gabor_angle column
Visual stimulus values |inpt| per session (get_stim_vals), sessions per set: {(np.float64(0.2), np.float64(0.6), np.float64(1.0)): 118, (np.float64(0.6), np.float64(0.8), np.float64(1.0)): 23}
Stimulus duration (ms): [1200.]
Start-to-start interval (s): median 7.66, 5-95% 6.11 to 14.69
Sampling: face camera {29.49: 331, 29.51: 52, 29.52: 7, 30.36: 1} Hz; wheel {100.3: 337, 100.1: 35, 100.2: 16, 100.0: 3} Hz
Raw running_speed per session (unit stated 'cm/s'): median of session medians 0.072, median of session 99th percentiles 0.342, largest 99th percentile 0.667
Raw post hoc pupil diameter (pixels): session medians from 21.5 to 84.7
Measures at stimulus onset after the code's filtering and normalisation (target trials):
count mean min 50% max
post_hoc_pupil_diameter 151452.0 0.741 0.030 0.756 1.000
face_energy 166907.0 0.460 0.003 0.437 0.976
running_speed 167568.0 0.136 -0.169 0.119 0.658
post_hoc_pupil_std10 145891.0 0.033 0.003 0.024 0.359
face_std_10 161926.0 0.074 0.004 0.062 0.328
running_std10 161333.0 0.059 0.000 0.052 0.282
Duplicate recordings: identical face motion energy arrays 0 | identical trial tables (stimulus times and outcomes) 0
The archive gives 13 subjects and 391 files; the files hold 13 mice, 8 on the auditory task and 5 on the visual task, all in training stages 5 and 6. No recording is deposited twice. running_speed is labelled cm/s, but its values (session medians around 0.07, 99th percentiles up to about 0.67) fit m/s, the unit of the paper (locomotion speed around 0.14 m/s) and of its code, which uses the values as they are; the values are taken as m/s and the label is recorded as a metadata error. Some sessions have a post hoc pupil series that is all missing; the paper's analyses that need pupil drop those trials, as here. The visual difficulty column is described as five angles, but only codes 0 to 2 occur. Most visual-task files also have a visual_gabor_angle column; for these, loadbehavior.get_stim_vals maps the codes 0, 1 and 2 to the session's distinct values of |angle - 45| / 45, largest first, and otherwise to 1, 0.8 and 0.6. The helper does the same, so the visual inputs are 1, 0.6 and 0.2 in most sessions (the sets per session are printed above).
The next cell checks that the measures share the clock of the trial table, with a contrast between conditions: if they do, the measures change after stimulus onset differently on trials with a response (a lick) than on trials without, and the difference shrinks when the same window is taken 4 s later, between trials.
resp = trials_all.hmm_choice != 2
clock = []
for subj, g in trials_all.groupby("subject"):
r = g.hmm_choice != 2
row = {"subject": subj}
for k in ["face_energy", "running_speed", "post_hoc_pupil_diameter"]:
for w in ["d0", "d4"]:
row[f"{k[:4]}_{w}"] = g.loc[r, f"{k}_{w}"].mean() - g.loc[~r, f"{k}_{w}"].mean()
clock.append(row)
clock = pd.DataFrame(clock).set_index("subject")
print("Response minus no-response trials: change of each measure from stimulus onset to 1.5 s after (d0),")
print("and over the same 1.5 s starting 4 s after onset (d4). Face and pupil normalised units, running m/s.")
print(clock.round(4).to_string())
for k in ["face", "runn", "post"]:
print(f"{k}: mice with a larger response contrast at onset than 4 s later (absolute values): "
f"{int((clock[k + '_d0'].abs() > clock[k + '_d4'].abs()).sum())} of {len(clock)}; "
f"median d0 {clock[k + '_d0'].median():.4f}, d4 {clock[k + '_d4'].median():.4f}")
Response minus no-response trials: change of each measure from stimulus onset to 1.5 s after (d0),
and over the same 1.5 s starting 4 s after onset (d4). Face and pupil normalised units, running m/s.
face_d0 face_d4 runn_d0 runn_d4 post_d0 post_d4
subject
BW041 0.0056 -0.0227 -0.0233 0.0090 -0.0033 0.0074
BW045 0.0588 0.0122 -0.0296 0.0040 0.0075 -0.0032
BW046 0.0252 -0.0095 -0.1593 0.0242 -0.0520 0.0148
BW048 0.0352 -0.0082 -0.1329 0.0807 -0.0395 0.0129
BW051 -0.0714 0.0108 -0.1641 0.0509 -0.0287 0.0195
BW053 0.0815 -0.0180 -0.0026 0.0025 0.0461 0.0007
BW058 0.0153 0.0063 -0.1134 0.0416 -0.0157 0.0061
BW059 0.0515 -0.0086 -0.3026 0.0868 -0.0261 -0.0050
BW061 -0.0791 -0.0013 -0.1167 -0.0223 -0.0200 -0.0003
BW062 -0.1024 -0.0052 -0.1024 0.0050 0.0118 -0.0194
BW064 0.0584 -0.0285 -0.1068 0.0082 -0.0040 -0.0062
BW067 -0.0423 0.0233 -0.1305 0.0829 -0.0084 0.0062
BW071 -0.0043 0.0010 -0.1426 0.0118 -0.0183 0.0014
face: mice with a larger response contrast at onset than 4 s later (absolute values): 12 of 13; median d0 0.0153, d4 -0.0052
runn: mice with a larger response contrast at onset than 4 s later (absolute values): 13 of 13; median d0 -0.1167, d4 0.0118
post: mice with a larger response contrast at onset than 4 s later (absolute values): 10 of 13; median d0 -0.0157, d4 0.0014
Sessions and trials (claim cohort-inclusion)¶
Of 31 mice, 13 met the inclusion criteria (8 auditory, 5 visual task), with 381 of their 500 stage 5 and 6 sessions included.
The criterion is the STAR Methods' (Training and inclusion criteria): stage 5 or 6, at least 100 rewards with at least 20% at each port, and at least 10 such sessions per mouse. Rewards are counted from the reward_left and reward_right series; the hits of the trial table are shown as a second count. The sessions that pass are the sessions of every later analysis.
s = sessions.copy()
s["rewards"] = s.reward_left + s.reward_right
s["min_port_frac"] = np.minimum(s.reward_left, s.reward_right) / s.rewards
hits = alltrials[alltrials.outcome == 1].groupby("path").target_port.agg(hits="size", hits_left=lambda x: int((x == 0).sum()))
s = s.merge(hits, left_on="path", right_index=True, how="left")
s["hit_min_frac"] = np.minimum(s.hits_left, s.hits - s.hits_left) / s.hits
s["pass_rewards"] = s.stage.isin(["S5", "S6"]) & (s.rewards >= 100) & (s.min_port_frac >= 0.2)
s["pass_hits"] = s.stage.isin(["S5", "S6"]) & (s.hits >= 100) & (s.hit_min_frac >= 0.2)
per_mouse = s.groupby("subject").pass_rewards.sum()
s["included"] = s.pass_rewards & s.subject.map(per_mouse >= 10)
print("Sessions in the dandiset:", len(s))
print("Pass by rewards:", int(s.pass_rewards.sum()), "| pass by hits:", int(s.pass_hits.sum()),
"| either:", int((s.pass_rewards | s.pass_hits).sum()))
print("Mice with at least 10 passing sessions:", int((per_mouse >= 10).sum()), "of", len(per_mouse))
mods = s.groupby("subject").session_type.first()
print("Mice by task:", mods.value_counts().to_dict())
print("Included sessions (rewards):", int(s.included.sum()))
print("\nSessions that fail either count:")
print(s.loc[~s.pass_rewards | ~s.pass_hits, ["path", "stage", "rewards", "min_port_frac", "hits", "hit_min_frac"]]
.round(3).to_string(index=False))
print("\nIncluded sessions per mouse:", s[s.included].groupby("subject").size().to_dict())
included = s[s.included]
trials = trials_all[trials_all.path.isin(included.path)].copy()
print("Target trials in included sessions:", len(trials))
Sessions in the dandiset: 391
Pass by rewards: 376 | pass by hits: 375 | either: 377
Mice with at least 10 passing sessions: 13 of 13
Mice by task: {'auditory': 8, 'visual': 5}
Included sessions (rewards): 376
Sessions that fail either count:
path stage rewards min_port_frac hits hit_min_frac
sub-BW041/sub-BW041_ses-20210905T142058.nwb S5 341 0.109 340 0.106
sub-BW045/sub-BW045_ses-20211108T091227.nwb S6 92 0.402 93 0.398
sub-BW045/sub-BW045_ses-20211116T093954.nwb S6 99 0.384 99 0.384
sub-BW045/sub-BW045_ses-20211117T085905.nwb S6 93 0.419 93 0.419
sub-BW045/sub-BW045_ses-20211118T084349.nwb S6 98 0.459 98 0.459
sub-BW045/sub-BW045_ses-20211211T141311.nwb S6 99 0.404 99 0.404
sub-BW046/sub-BW046_ses-20211116T093954.nwb S6 87 0.402 87 0.402
sub-BW046/sub-BW046_ses-20211122T113133.nwb S6 99 0.394 99 0.394
sub-BW053/sub-BW053_ses-20211228T140451.nwb S5 220 0.195 222 0.203
sub-BW053/sub-BW053_ses-20211231T155717.nwb S5 97 0.454 97 0.454
sub-BW053/sub-BW053_ses-20220121T144423.nwb S6 97 0.392 97 0.392
sub-BW059/sub-BW059_ses-20220220T101110.nwb S5 102 0.431 97 0.402
sub-BW062/sub-BW062_ses-20220802T094351.nwb S5 66 0.409 66 0.409
sub-BW064/sub-BW064_ses-20220714T114217.nwb S5 78 0.218 78 0.218
sub-BW064/sub-BW064_ses-20220815T094736.nwb S6 86 0.302 86 0.302
sub-BW064/sub-BW064_ses-20220816T102856.nwb S6 76 0.276 76 0.276
sub-BW071/sub-BW071_ses-20220710T124035.nwb S5 150 0.487 79 0.481
Included sessions per mouse: {'BW041': 24, 'BW045': 79, 'BW046': 32, 'BW048': 35, 'BW051': 35, 'BW053': 29, 'BW058': 12, 'BW059': 25, 'BW061': 23, 'BW062': 16, 'BW064': 18, 'BW067': 23, 'BW071': 25}
Target trials in included sessions: 164503
Result: Partially reproduced. The dandiset holds 13 mice, 8 on the auditory and 5 on the visual task, all of which have at least 10 passing sessions, as the paper reports. Of its 391 sessions, 376 pass the criteria counted from the reward series (375 from the hits), not 381; the paper's 500 stage 5 and 6 sessions and its 18 excluded mice are not in the dandiset, so which 381 sessions were analysed cannot be told. The 376 sessions are used below.
Training before stage 5 (claim training-duration)¶
Mice that reached proficiency took 21 +/- 8 training sessions to reach stage 5.
The dandiset holds only stage 5 and 6 sessions (data check), so the sessions before stage 5 cannot be counted. Blocked: data missing.
Distractor responses (claim distractor-response-rate)¶
In stage 6, responses to distractors of the untrained modality were below 10%, with no discernible change in task performance.
Distractor trials have outcome 4 (false_alarm, a lick) or 8 (correct_reject, no lick); the rate is the fraction with a lick. All stage 6 sessions in the dandiset are used, and the included ones are shown too.
d = alltrials[(alltrials.stimulus_type == 1)].merge(s[["path", "subject", "included"]], on="path")
d = d[d.stage == "S6"]
d["lick"] = (d.outcome == 4).astype(float)
for label, dd in [("all stage 6 sessions", d), ("included stage 6 sessions", d[d.included])]:
pooled = dd.groupby("subject").lick.mean()
per_sess = dd.groupby(["subject", "path"]).lick.mean().groupby("subject").mean()
tab = pd.DataFrame({"sessions": dd.groupby("subject").path.nunique(), "distractor trials": dd.groupby("subject").size(),
"(a) pooled": pooled.round(4), "(b) mean over sessions": per_sess.round(4),
"max session": dd.groupby(["subject", "path"]).lick.mean().groupby("subject").max().round(4)})
print(label)
print(tab.to_string())
print(f"(c) pooled over all mice: {dd.lick.mean():.4f} over {len(dd)} trials; mice below 0.10 under (a): "
f"{int((pooled < 0.1).sum())} of {len(pooled)}, under (b): {int((per_sess < 0.1).sum())} of {len(per_sess)}\n")
st6 = trials_all.merge(s[["path", "included"]], on="path")
perf = st6[st6.hmm_choice != 2].assign(correct=lambda x: x.outcome == 1).groupby(["subject", "stage"]).correct.mean().unstack()
print("Not tested, for reference: fraction correct of responded target trials by stage")
print(perf.round(3).to_string())
all stage 6 sessions
sessions distractor trials (a) pooled (b) mean over sessions max session
subject
BW041 10 3719 0.0347 0.0334 0.0644
BW045 74 20216 0.0561 0.0554 0.1522
BW046 16 4977 0.0203 0.0178 0.0597
BW048 15 4738 0.0441 0.0441 0.0793
BW051 10 3961 0.0134 0.0135 0.0403
BW053 9 2498 0.1621 0.1525 0.3233
BW059 7 4158 0.0808 0.0827 0.1402
BW064 9 2770 0.0285 0.0306 0.0742
(c) pooled over all mice: 0.0520 over 47037 trials; mice below 0.10 under (a): 7 of 8, under (b): 7 of 8
included stage 6 sessions
sessions distractor trials (a) pooled (b) mean over sessions max session
subject
BW041 10 3719 0.0347 0.0334 0.0644
BW045 69 19460 0.0558 0.0547 0.1522
BW046 14 4648 0.0211 0.0192 0.0597
BW048 15 4738 0.0441 0.0441 0.0793
BW051 10 3961 0.0134 0.0135 0.0403
BW053 8 2356 0.1702 0.1680 0.3233
BW059 7 4158 0.0808 0.0827 0.1402
BW064 7 2226 0.0279 0.0305 0.0742
(c) pooled over all mice: 0.0524 over 45266 trials; mice below 0.10 under (a): 7 of 8, under (b): 7 of 8
Not tested, for reference: fraction correct of responded target trials by stage stage S5 S6 subject BW041 0.936 0.943 BW045 0.798 0.786 BW046 0.792 0.801 BW048 0.836 0.837 BW051 0.888 0.886 BW053 0.760 0.709 BW058 0.793 NaN BW059 0.825 0.830 BW061 0.707 NaN BW062 0.798 NaN BW064 0.816 0.762 BW067 0.709 NaN BW071 0.838 NaN
Result: Reproduced. Under reading (c), 5.2% of the 47,037 stage 6 distractor trials had a lick (5.2% on the included sessions). Per mouse, 7 of the 8 mice with stage 6 sessions are below 10% under readings (a) and (b); BW053 is at 16.2% pooled (15.3% over sessions). The fraction correct of responded target trials changes between stages by at most 5.4 points (BW064, 0.816 to 0.762), shown but not tested.
GLM-HMM: number of states (claim figs2-n-states)¶
Cross-validated GLM-HMMs fit separately to each of the 13 mice gave between 3 and 5 states per mouse.
Per mouse, the target trials of the included sessions (distractor trials dropped, no-response trials kept) with inputs [stimulus value, 1] and choices left, right or none (claims file, M1 and M2). The GLM-HMM is the paper's ssm model, reimplemented in hulsey.py (its docstring gives the source of each part). 5-fold cross-validation over sessions, 1 to 7 states, MLE and MAP fits (2_n_states.py). 10 random initialisations per fold and number of states, as the paper. The number of states is the plateau rule of 3_generate_hmms.py on the test log likelihood per trial (mean over folds, best initialisation): reading (b) on the MLE curve alone, the code's; reading (a) the larger of the MLE and MAP plateaus.
N_INIT_CV = 10
t0 = time.time()
subjects = sorted(trials.subject.unique())
seqs = {subj: hulsey.sequences(trials[trials.subject == subj]) for subj in subjects}
tasks, keys = [], []
for i, subj in enumerate(subjects):
Xs, ys = seqs[subj]
tt, kk = hulsey.cv_tasks(Xs, ys, N_INIT_CV, seed=1000 + i)
tasks += tt
keys += [(subj,) + k for k in kk]
order = np.argsort([-sum(len(x) for x in t[1]) * t[0] for t in tasks]) # largest fits first
res = hulsey.run_parallel(hulsey.cv_test_ll, [tasks[j] for j in order])
cv = {keys[j]: r for j, r in zip(order, res)}
print(f"{len(tasks)} cross-validation fits in {(time.time() - t0) / 60:.1f} min; "
f"{sum(r[2] == 1000 for r in res)} stopped at 1000 iterations")
n_states, curves = {}, {}
rows = []
for subj in subjects:
row = {"subject": subj, "task": mods[subj], "trials": sum(len(x) for x in seqs[subj][0]),
"sessions": len(seqs[subj][0])}
for use_map, name in [(False, "MLE"), (True, "MAP")]:
a = np.full((N_INIT_CV, 7, 5), np.nan)
for (sb, f, K, z, m), v in cv.items():
if sb == subj and m == use_map:
a[z, K - 1, f] = v[1]
n, curve = hulsey.select_states(a)
curves[(subj, name)] = curve
row[name] = n
row["(a) max"] = max(row["MLE"], row["MAP"])
row["(b) MLE"] = row["MLE"]
rows.append(row)
nst = pd.DataFrame(rows).set_index("subject")
print(nst.to_string())
for r in ["(a) max", "(b) MLE"]:
print(f"reading {r}: range {nst[r].min()} to {nst[r].max()}; counts {nst[r].value_counts().sort_index().to_dict()}")
print("\nTest log likelihood per trial (MLE, mean over folds, best initialisation), 1 to 7 states:")
print(pd.DataFrame({subj: curves[(subj, 'MLE')] for subj in subjects}, index=range(1, 8)).T.round(4).to_string())
print("\nThe same for MAP:")
print(pd.DataFrame({subj: curves[(subj, 'MAP')] for subj in subjects}, index=range(1, 8)).T.round(4).to_string())
n_states = nst["(b) MLE"].to_dict()
9100 cross-validation fits in 9.4 min; 38 stopped at 1000 iterations
task trials sessions MLE MAP (a) max (b) MLE
subject
BW041 visual 11542 24 3 3 3 3
BW045 auditory 23179 79 4 4 4 4
BW046 visual 13795 32 4 4 4 4
BW048 visual 14109 35 4 4 4 4
BW051 visual 13935 35 3 3 3 3
BW053 auditory 11858 29 4 3 4 4
BW058 visual 5623 12 3 3 3 3
BW059 auditory 13184 25 4 3 4 4
BW061 auditory 13793 23 4 3 4 4
BW062 auditory 11142 16 3 3 3 3
BW064 auditory 7034 18 3 3 3 3
BW067 auditory 13603 23 4 3 4 4
BW071 auditory 11706 25 4 4 4 4
reading (a) max: range 3 to 4; counts {3: 5, 4: 8}
reading (b) MLE: range 3 to 4; counts {3: 5, 4: 8}
Test log likelihood per trial (MLE, mean over folds, best initialisation), 1 to 7 states:
1 2 3 4 5 6 7
BW041 -0.7114 -0.4855 -0.4228 -0.4306 -0.4413 -0.4530 -0.4664
BW045 -0.8578 -0.6413 -0.6110 -0.5928 -0.5928 -0.5957 -0.6003
BW046 -0.7954 -0.6040 -0.5666 -0.5413 -0.5462 -0.5531 -0.5620
BW048 -0.7124 -0.6079 -0.5533 -0.5367 -0.5358 -0.5359 -0.5428
BW051 -0.6546 -0.5534 -0.5362 -0.5375 -0.5412 -0.5477 -0.5568
BW053 -0.9080 -0.6980 -0.6860 -0.6850 -0.6859 -0.6976 -0.7091
BW058 -0.7975 -0.6477 -0.6312 -0.6331 -0.6497 -0.6757 -0.7040
BW059 -0.7503 -0.6729 -0.6505 -0.6436 -0.6504 -0.6512 -0.6575
BW061 -1.0153 -0.7655 -0.6686 -0.6646 -0.6745 -0.6815 -0.6861
BW062 -0.8491 -0.7061 -0.6694 -0.6747 -0.6799 -0.6899 -0.7050
BW064 -0.9800 -0.9033 -0.8807 -0.8920 -0.9128 -0.9337 -0.9581
BW067 -0.9500 -0.7672 -0.6886 -0.6834 -0.6859 -0.6917 -0.6991
BW071 -0.8312 -0.5530 -0.5342 -0.5282 -0.5291 -0.5382 -0.5516
The same for MAP:
1 2 3 4 5 6 7
BW041 -0.7009 -0.4701 -0.4086 -0.4198 -0.4386 -0.4618 -0.4861
BW045 -0.8525 -0.6335 -0.6018 -0.5872 -0.5899 -0.5989 -0.6146
BW046 -0.7866 -0.5916 -0.5558 -0.5318 -0.5433 -0.5629 -0.5870
BW048 -0.7038 -0.5945 -0.5384 -0.5250 -0.5295 -0.5424 -0.5613
BW051 -0.6463 -0.5405 -0.5228 -0.5263 -0.5362 -0.5518 -0.5746
BW053 -0.8973 -0.6823 -0.6676 -0.6728 -0.6815 -0.7049 -0.7345
BW058 -0.7742 -0.6150 -0.5957 -0.6041 -0.6396 -0.6889 -0.7474
BW059 -0.7409 -0.6579 -0.6295 -0.6310 -0.6446 -0.6631 -0.6871
BW061 -1.0063 -0.7524 -0.6554 -0.6587 -0.6800 -0.6985 -0.7239
BW062 -0.8381 -0.6910 -0.6553 -0.6641 -0.6803 -0.7040 -0.7319
BW064 -0.9612 -0.8711 -0.8422 -0.8553 -0.8860 -0.9266 -0.9666
BW067 -0.9408 -0.7521 -0.6720 -0.6712 -0.6788 -0.6961 -0.7205
BW071 -0.8207 -0.5381 -0.5205 -0.5187 -0.5315 -0.5551 -0.5828
Result: Reproduced. Under both readings every mouse has 3 or 4 states (five mice with 3, eight with 4), inside the claimed 3 to 5; no mouse has 5 states. The MAP plateau is never larger than the MLE plateau, so reading (a) gives the same counts as reading (b).
Final models and state labels¶
As 3_generate_hmms.py: per mouse, 10 MLE fits on all sessions with the selected number of states (reading (b), which is also reading (a) when they agree), the one with the highest log likelihood kept; trials in a state when its posterior is at least 0.8 (flex_hmm.get_posterior_probs); states ordered optimal, disengaged, others (flex_hmm.permute_hmm). Sub-optimal states are then matched to the four templates of Figure2.py, on the sessions in which the optimal state occurs.
t0 = time.time()
tasks, owner = [], []
for i, subj in enumerate(subjects):
Xs, ys = seqs[subj]
tt = hulsey.final_tasks(Xs, ys, n_states[subj], 10, seed=2000 + i)
tasks += tt
owner += [subj] * len(tt)
res = hulsey.run_parallel(hulsey.fit_one, tasks)
models = {}
for subj in subjects:
rr = [res[j] for j in range(len(res)) if owner[j] == subj]
models[subj] = rr[int(np.argmax([r[1] for r in rr]))][0]
print(f"final fits in {(time.time() - t0) / 60:.1f} min")
post_cols = []
parts = []
for subj in subjects:
g = trials[trials.subject == subj]
m = models[subj]
p0 = hulsey.posterior_states(m, g)
ps, _ = hulsey.psychos(g, p0.hmm_state.values)
perm = hulsey.permutation_optimal_disengaged(ps)
m.permute(perm)
parts.append(hulsey.posterior_states(m, g))
post = pd.concat(parts)
trials = trials.drop(columns=[c for c in trials.columns if c.startswith("state_") or c == "hmm_state"], errors="ignore")
trials = trials.join(post)
# sessions in which the optimal state occurs (Figure2.py) and with more than 10 optimal trials (Figure3.py etc.)
opt_n = trials.groupby("path").hmm_state.apply(lambda x: int((x == 0).sum()))
trials["sess_has_opt"] = trials.path.map(opt_n > 0)
trials["sess_opt10"] = trials.path.map(opt_n > 10) & trials.path.map(included.set_index("path").pupil_raw_median.notna())
types = {}
rows = []
for subj in subjects:
g = trials[(trials.subject == subj) & trials.sess_has_opt]
ps, inp = hulsey.psychos(g, g.hmm_state.values)
lab, sim = hulsey.state_types(ps)
lab_all, _ = hulsey.state_types(ps, fixed_first_two=False)
types[subj] = lab
occ = g.hmm_state.value_counts(normalize=True).sort_index()
rows.append({"subject": subj, "states": n_states[subj], "inputs": len(inp),
"labels": [hulsey.STATE_NAMES[k] for k in lab],
"nearest of all six templates": [hulsey.STATE_NAMES[k] for k in lab_all],
"occupancy": [round(occ.get(k, 0), 3) for k in range(n_states[subj])],
"indeterminate": round(g.hmm_state.isna().mean(), 3)})
print(pd.DataFrame(rows).set_index("subject").to_string())
trials["state_type"] = np.nan
for subj in subjects:
sel = trials.subject == subj
st = trials.loc[sel, "hmm_state"]
trials.loc[sel, "state_type"] = st.map(lambda k: types[subj][int(k)] if k == k else np.nan)
print("\nSessions with the optimal state:", int(trials.groupby("path").sess_has_opt.first().sum()),
"| with more than 10 optimal trials and a post hoc pupil (M5):", int(trials.groupby("path").sess_opt10.first().sum()),
"of", trials.path.nunique())
m5 = trials[trials.sess_opt10]
print("Trials in the M5 sessions:", len(m5), "| with all six measures:", len(m5.dropna(subset=["post_hoc_pupil_diameter",
"post_hoc_pupil_std10", "face_energy", "face_std_10", "running_speed", "running_std10"])),
"| of these in a discrete state:", int(m5.dropna(subset=["post_hoc_pupil_diameter", "post_hoc_pupil_std10", "face_energy",
"face_std_10", "running_speed", "running_std10"]).hmm_state.notna().sum()))
print("\nTrials by state (0 optimal, 1 disengaged, 2 or more other; NaN indeterminate):",
trials.hmm_state.value_counts(dropna=False).sort_index().to_dict())
final fits in 0.1 min
states inputs labels nearest of all six templates occupancy indeterminate subject BW041 3 8 [optimal, disengaged, avoid left] [optimal, disengaged, avoid left] [0.678, 0.128, 0.194] 0.070 BW045 4 6 [optimal, disengaged, bias right, avoid right] [optimal, disengaged, bias right, avoid right] [0.422, 0.173, 0.319, 0.086] 0.238 BW046 4 6 [optimal, disengaged, bias right, bias left] [optimal, disengaged, bias right, bias left] [0.592, 0.082, 0.118, 0.207] 0.124 BW048 4 6 [optimal, disengaged, avoid left, bias left] [optimal, disengaged, avoid left, bias left] [0.719, 0.048, 0.095, 0.139] 0.108 BW051 3 6 [optimal, disengaged, avoid right] [optimal, disengaged, avoid right] [0.616, 0.05, 0.334] 0.143 BW053 4 6 [optimal, disengaged, bias left, bias right] [optimal, disengaged, bias left, bias right] [0.546, 0.224, 0.114, 0.116] 0.129 BW058 3 6 [optimal, disengaged, avoid right] [optimal, disengaged, avoid right] [0.745, 0.095, 0.161] 0.052 BW059 4 6 [optimal, disengaged, no trials, avoid left] [optimal, disengaged, no trials, avoid left] [0.703, 0.076, 0, 0.22] 0.171 BW061 4 6 [optimal, disengaged, bias left, bias right] [optimal, disengaged, bias left, bias right] [0.423, 0.273, 0.026, 0.278] 0.071 BW062 3 6 [optimal, disengaged, avoid left] [optimal, disengaged, avoid left] [0.539, 0.112, 0.348] 0.172 BW064 3 6 [optimal, disengaged, avoid right] [optimal, disengaged, avoid right] [0.43, 0.331, 0.239] 0.223 BW067 4 6 [optimal, disengaged, bias left, bias right] [optimal, disengaged, bias left, bias right] [0.491, 0.142, 0.33, 0.037] 0.104 BW071 4 6 [optimal, disengaged, bias left, bias right] [optimal, disengaged, bias left, bias right] [0.611, 0.174, 0.091, 0.125] 0.102 Sessions with the optimal state: 352 | with more than 10 optimal trials and a post hoc pupil (M5): 313 of 376
Trials in the M5 sessions: 136640 | with all six measures: 129568 | of these in a discrete state: 111714
Trials by state (0 optimal, 1 disengaged, 2 or more other; NaN indeterminate): {0.0: 75695, 1.0: 20620, 2.0: 30246, 3.0: 15656, nan: 22286}
# a state with no discrete trial in the sessions where the optimal state occurs is typed on all sessions
for subj in subjects:
if (types[subj] == -1).any():
g = trials[trials.subject == subj]
ps, _ = hulsey.psychos(g, g.hmm_state.values)
lab_all_sess, _ = hulsey.state_types(ps)
print(f"{subj}: states with no discrete trial in sessions with the optimal state:",
list(np.where(types[subj] == -1)[0]), "| typed on all sessions:",
[hulsey.STATE_NAMES[k] for k in lab_all_sess],
"| trials per state over all sessions:", g.hmm_state.value_counts().sort_index().to_dict())
types[subj] = np.where(types[subj] == -1, lab_all_sess, types[subj])
sel = trials.subject == subj
trials.loc[sel, "state_type"] = trials.loc[sel, "hmm_state"].map(lambda k: types[subj][int(k)] if k == k else np.nan)
BW059: states with no discrete trial in sessions with the optimal state: [np.int64(2)] | typed on all sessions: ['optimal', 'disengaged', 'avoid right', 'avoid left'] | trials per state over all sessions: {0.0: 7420, 1.0: 820, 2.0: 432, 3.0: 2327}
State types (claim fig2a-state-types)¶
Six stereotyped performance states emerged across mice. Every mouse had an optimal and a disengaged state and at least one sub-optimal state (bias left or right, avoid left or right).
rows = []
for subj in subjects:
lab = types[subj]
rows.append({"subject": subj, "states": len(lab), "optimal": int((lab == 0).sum()), "disengaged": int((lab == 1).sum()),
"other states": int((lab >= 2).sum()),
"other states matched to a sub-optimal template": int(((lab >= 2) & (lab <= 5)).sum())})
ft = pd.DataFrame(rows).set_index("subject")
print(ft.to_string())
ok = (ft.optimal == 1) & (ft.disengaged == 1) & (ft["other states"] >= 1) & (ft["other states"] == ft["other states matched to a sub-optimal template"])
print("\nMice with an optimal, a disengaged and at least one other state, each other state typed:", int(ok.sum()), "of", len(ft))
counts = {hulsey.STATE_NAMES[k]: int(sum((types[s] == k).any() for s in subjects)) for k in range(2, 6)}
print("Mice with each sub-optimal type:", counts, "(Figure 2A: bias left 6, avoid right 4, bias right 5, avoid left 5)")
print("Sub-optimal states in all:", int(sum((types[s] >= 2).sum() for s in subjects)))
print("\nDistance of each state's psychometric (sessions with the optimal state) to the six templates of Figure2.py;"
" the smallest per row is the nearest template:")
for subj in ["BW061", "BW048"]:
g = trials[(trials.subject == subj) & trials.sess_has_opt]
ps, inp = hulsey.psychos(g, g.hmm_state.values)
_, sim = hulsey.state_types(ps, fixed_first_two=False)
print(f"{subj}, labels {[hulsey.STATE_NAMES[k] for k in types[subj]]}:")
print(pd.DataFrame(sim, columns=hulsey.STATE_NAMES[:6]).round(2).to_string())
for k in range(len(ps)):
print(f" state {k} ({hulsey.STATE_NAMES[types[subj][k]]}), P(left), P(right), P(no response) by input {np.round(inp, 2).tolist()}:")
print(" ", np.round(ps[k], 2).tolist())
states optimal disengaged other states other states matched to a sub-optimal template
subject
BW041 3 1 1 1 1
BW045 4 1 1 2 2
BW046 4 1 1 2 2
BW048 4 1 1 2 2
BW051 3 1 1 1 1
BW053 4 1 1 2 2
BW058 3 1 1 1 1
BW059 4 1 1 2 2
BW061 4 1 1 2 2
BW062 3 1 1 1 1
BW064 3 1 1 1 1
BW067 4 1 1 2 2
BW071 4 1 1 2 2
Mice with an optimal, a disengaged and at least one other state, each other state typed: 13 of 13
Mice with each sub-optimal type: {'bias left': 6, 'avoid right': 5, 'bias right': 6, 'avoid left': 4} (Figure 2A: bias left 6, avoid right 4, bias right 5, avoid left 5)
Sub-optimal states in all: 21
Distance of each state's psychometric (sessions with the optimal state) to the six templates of Figure2.py; the smallest per row is the nearest template:
BW061, labels ['optimal', 'disengaged', 'bias left', 'bias right']:
optimal disengaged bias left avoid right bias right avoid left
0 2.52 11.95 6.57 7.48 5.48 6.98
1 10.59 2.92 10.71 6.86 10.36 6.66
2 4.14 11.66 2.53 5.83 9.81 9.97
3 4.13 11.86 9.85 9.95 2.29 6.04
state 0 (optimal), P(left), P(right), P(no response) by input [-0.97, -0.62, -0.43, 0.43, 0.62, 0.97]:
[[0.87, 0.73, 0.64, 0.26, 0.15, 0.07], [0.13, 0.27, 0.36, 0.74, 0.84, 0.92], [0.0, 0.0, 0.01, 0.0, 0.0, 0.02]]
state 1 (disengaged), P(left), P(right), P(no response) by input [-0.97, -0.62, -0.43, 0.43, 0.62, 0.97]:
[[0.11, 0.09, 0.1, 0.11, 0.11, 0.12], [0.14, 0.16, 0.12, 0.13, 0.12, 0.15], [0.74, 0.75, 0.78, 0.76, 0.77, 0.74]]
state 2 (bias left), P(left), P(right), P(no response) by input [-0.97, -0.62, -0.43, 0.43, 0.62, 0.97]:
[[0.96, 0.97, 0.98, 0.74, 0.68, 0.4], [0.04, 0.03, 0.02, 0.22, 0.29, 0.5], [0.0, 0.0, 0.0, 0.04, 0.03, 0.1]]
state 3 (bias right), P(left), P(right), P(no response) by input [-0.97, -0.62, -0.43, 0.43, 0.62, 0.97]:
[[0.47, 0.3, 0.24, 0.03, 0.02, 0.01], [0.49, 0.69, 0.76, 0.97, 0.97, 0.99], [0.04, 0.01, 0.01, 0.01, 0.01, 0.0]]
BW048, labels ['optimal', 'disengaged', 'avoid left', 'bias left']:
optimal disengaged bias left avoid right bias right avoid left
0 1.37 11.97 5.67 6.49 6.36 6.85
1 11.94 0.21 11.81 6.07 11.98 6.08
2 5.68 8.18 11.14 10.90 4.68 2.96
3 5.82 10.56 1.69 4.68 11.75 11.69
state 0 (optimal), P(left), P(right), P(no response) by input [-1.0, -0.6, -0.2, 0.2, 0.6, 1.0]:
[[0.97, 0.95, 0.83, 0.33, 0.07, 0.03], [0.03, 0.05, 0.17, 0.67, 0.93, 0.97], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]
state 1 (disengaged), P(left), P(right), P(no response) by input [-1.0, -0.6, -0.2, 0.2, 0.6, 1.0]:
[[0.03, 0.0, 0.0, 0.04, 0.02, 0.01], [0.0, 0.0, 0.01, 0.0, 0.0, 0.0], [0.97, 1.0, 0.99, 0.96, 0.98, 0.99]]
state 2 (avoid left), P(left), P(right), P(no response) by input [-1.0, -0.6, -0.2, 0.2, 0.6, 1.0]:
[[0.18, 0.14, 0.06, 0.03, 0.01, 0.0], [0.25, 0.26, 0.37, 0.8, 0.99, 0.99], [0.57, 0.6, 0.58, 0.16, 0.0, 0.0]]
state 3 (bias left), P(left), P(right), P(no response) by input [-1.0, -0.6, -0.2, 0.2, 0.6, 1.0]:
[[1.0, 0.99, 0.98, 0.89, 0.69, 0.61], [0.0, 0.0, 0.0, 0.01, 0.04, 0.07], [0.0, 0.01, 0.02, 0.1, 0.27, 0.32]]
Result: Reproduced. All 13 mice have an optimal state, a disengaged state and at least one other state, and every other state matches a sub-optimal template (21 sub-optimal states). The mice per sub-optimal type are bias left 6, avoid right 5, bias right 6 and avoid left 4, against 6, 4, 5 and 5 in Figure 2A; these counts do not decide the claim. In BW059 one state has no trial with a posterior of at least 0.8 in the sessions where the optimal state occurs; it is typed on all sessions (avoid right), an assumption. Under M3 each other state is matched only to the four sub-optimal templates, so the matching itself cannot fail. The distances to all six templates printed for BW061 and BW048 show each of their states nearest the template of its label (for example BW048 state 2, avoid left, 2.96 against 4.68 for the next nearest).
Occupancy (claim fig2b-optimal-occupancy)¶
Mice spent the largest proportion of trials in the optimal state (around half), more than in the disengaged, sub-optimal or indeterminate states; 13.1% +/- 1.6% of trials were indeterminate.
Per mouse, trials of the sessions in which the optimal state occurs, pooled (Figure2.py).
occ = []
for subj in subjects:
g = trials[(trials.subject == subj) & trials.sess_has_opt]
h = g.hmm_state
occ.append({"subject": subj, "optimal": (h == 0).mean(), "disengaged": (h == 1).mean(),
"sub-optimal": (h > 1).mean(), "indeterminate": h.isna().mean(), "trials": len(g)})
occ = pd.DataFrame(occ).set_index("subject")
print((occ.drop(columns="trials") * 100).round(1).assign(trials=occ.trials).to_string())
m_ = occ.drop(columns="trials").mean() * 100
se_ = occ.drop(columns="trials").sem() * 100
print("\nmean +/- SEM over mice (%):", {k: f"{m_[k]:.1f} +/- {se_[k]:.1f}" for k in m_.index})
for other in ["disengaged", "sub-optimal", "indeterminate"]:
w = stats.wilcoxon(occ.optimal, occ[other])
print(f"optimal vs {other}: optimal larger in {int((occ.optimal > occ[other]).sum())} of 13 mice, "
f"Wilcoxon signed-rank p = {w.pvalue:.2g}")
optimal disengaged sub-optimal indeterminate trials
subject
BW041 63.0 11.9 18.1 7.0 11542
BW045 32.1 13.2 30.9 23.8 22472
BW046 51.9 7.2 28.5 12.4 13795
BW048 64.1 4.3 20.8 10.8 14109
BW051 52.8 4.3 28.6 14.3 13935
BW053 47.6 19.5 20.0 12.9 10600
BW058 70.6 9.0 15.2 5.2 5623
BW059 58.3 6.3 18.3 17.1 12737
BW061 39.3 25.4 28.3 7.1 7807
BW062 44.7 9.3 28.8 17.2 11142
BW064 33.4 25.7 18.6 22.3 7034
BW067 44.0 12.7 32.9 10.4 12197
BW071 54.9 15.6 19.3 10.2 9952
mean +/- SEM over mice (%): {'optimal': '50.5 +/- 3.3', 'disengaged': '12.7 +/- 2.0', 'sub-optimal': '23.7 +/- 1.7', 'indeterminate': '13.1 +/- 1.6'}
optimal vs disengaged: optimal larger in 13 of 13 mice, Wilcoxon signed-rank p = 0.00024
optimal vs sub-optimal: optimal larger in 13 of 13 mice, Wilcoxon signed-rank p = 0.00024
optimal vs indeterminate: optimal larger in 13 of 13 mice, Wilcoxon signed-rank p = 0.00024
Result: Reproduced. The optimal state holds the largest share of trials (50.5% +/- 3.3%), more than the disengaged (12.7%), sub-optimal (23.7%) and indeterminate (13.1%) states, larger in 13 of 13 mice for each (Wilcoxon signed-rank p = 0.00024 for all three). The indeterminate share, 13.1% +/- 1.6%, matches the expected 13.1% +/- 1.6%.
Dwell times (claim fig2c-dwell-times)¶
Epochs of the optimal state lasted longer than disengaged and sub-optimal epochs (about 37 against about 15 trials), and the median transition between discrete states took 3.4 +/- 0.3 trials.
Epochs from flex_hmm.get_dwell_times on the sessions in which the optimal state occurs (Figure2.py), with (y) the code's dwell and (x) the run length. A disengaged epoch longer than 40 trials (code's dwell) that ends a session is left out (Figure2.py). The per-mouse median is the first length at which the cumulative distribution exceeds 0.499. get_dwell_times also adds a last epoch of dwell 0 when a session's last trial differs from the one before it, which adds a zero to the indeterminate dwells; the code's output is used, and a variant without these epochs is shown.
def cdf_median(v):
v = np.asarray(v)
for i in range(500):
if np.mean(v <= i) > 0.499:
return i
return np.nan
ep_rows, ind_rows = [], []
for subj in subjects:
g = trials[(trials.subject == subj) & trials.sess_has_opt].sort_values(["path", "trial_in_session"])
for path, sg in g.groupby("path"):
dt = hulsey.dwell_table(sg.hmm_state.values)
if len(dt) == 0:
continue
gap = dt.first_trial.values - np.r_[0, dt.last_trial.values[:-1]]
ind_rows += [{"subject": subj, "gap": x} for x in gap]
drop_last = dt.state.iloc[-1] == 1 and dt.dwell.iloc[-1] > 40
dt = dt.iloc[:-1] if drop_last else dt
ep_rows.append(dt.assign(subject=subj, path=path))
ep = pd.concat(ep_rows, ignore_index=True)
ep["group"] = np.where(ep.state == 0, "optimal", np.where(ep.state == 1, "disengaged", "sub-optimal"))
ind = pd.DataFrame(ind_rows)
res_dw = {}
for col, cname in [("dwell", "(y) code dwell"), ("run", "(x) run length")]:
med = ep.groupby(["subject", "group"])[col].apply(cdf_median).unstack()
mean = ep.groupby(["subject", "group"])[col].mean().unstack()
pooled = ep.groupby("group")[col].mean()
res_dw[cname] = (med, mean, pooled)
print(cname)
print(" (a) mean over mice of per-mouse median:", med.mean().round(1).to_dict())
print(" (b) mean over mice of per-mouse mean: ", mean.mean().round(1).to_dict())
print(" (c) mean over all epochs pooled: ", pooled.round(1).to_dict())
for other in ["disengaged", "sub-optimal"]:
for nm, tab in [("medians", med), ("means", mean)]:
w = stats.wilcoxon(tab.optimal, tab[other])
print(f" optimal vs {other} ({nm}): optimal longer in {int((tab.optimal > tab[other]).sum())} of 13, p = {w.pvalue:.2g}")
ind_med = ind.groupby("subject").gap.apply(cdf_median)
print("\nMedian indeterminate dwell per mouse:", ind_med.to_dict())
print(f"mean +/- SEM over mice: {ind_med.mean():.2f} +/- {ind_med.sem():.2f} trials")
print("Per-mouse medians (code dwell):")
print(res_dw["(y) code dwell"][0].assign(indeterminate=ind_med).to_string())
print("Epochs per group:", ep.group.value_counts().to_dict())
# get_dwell_times adds a last epoch of dwell 0 (run 1) when a session's last trial differs from the one before
# it; it also adds a zero gap to the indeterminate dwells. Variant without these epochs and gaps:
spur = ep.dwell == 0
print(f"\nEpochs of dwell 0 (the extra last epoch): {int(spur.sum())} of {len(ep)}")
ind2 = []
for subj in subjects:
g = trials[(trials.subject == subj) & trials.sess_has_opt].sort_values(["path", "trial_in_session"])
for path, sg in g.groupby("path"):
dt = hulsey.dwell_table(sg.hmm_state.values)
dt = dt[dt.dwell > 0]
if len(dt) == 0:
continue
gap = dt.first_trial.values - np.r_[0, dt.last_trial.values[:-1]]
ind2 += [{"subject": subj, "gap": x} for x in gap]
ind_med2 = pd.DataFrame(ind2).groupby("subject").gap.apply(cdf_median)
print(f"Without them: median indeterminate dwell {ind_med2.mean():.2f} +/- {ind_med2.sem():.2f} trials")
ep2 = ep[~spur]
for col, cname in [("dwell", "(y) code dwell"), ("run", "(x) run length")]:
med2 = ep2.groupby(["subject", "group"])[col].apply(cdf_median).unstack()
print(f" {cname}, (a) mean over mice of per-mouse median:", med2.mean().round(1).to_dict())
(y) code dwell
(a) mean over mice of per-mouse median: {'disengaged': 14.0, 'optimal': 43.8, 'sub-optimal': 19.1}
(b) mean over mice of per-mouse mean: {'disengaged': 21.1, 'optimal': 74.1, 'sub-optimal': 30.1}
(c) mean over all epochs pooled: {'disengaged': 14.6, 'optimal': 50.3, 'sub-optimal': 20.7}
optimal vs disengaged (medians): optimal longer in 13 of 13, p = 0.00024
optimal vs disengaged (means): optimal longer in 13 of 13, p = 0.00024
optimal vs sub-optimal (medians): optimal longer in 13 of 13, p = 0.00024
optimal vs sub-optimal (means): optimal longer in 13 of 13, p = 0.00024
(x) run length
(a) mean over mice of per-mouse median: {'disengaged': 13.8, 'optimal': 42.8, 'sub-optimal': 18.1}
(b) mean over mice of per-mouse mean: {'disengaged': 20.8, 'optimal': 73.2, 'sub-optimal': 29.2}
(c) mean over all epochs pooled: {'disengaged': 14.2, 'optimal': 49.3, 'sub-optimal': 19.7}
optimal vs disengaged (medians): optimal longer in 13 of 13, p = 0.00024
optimal vs disengaged (means): optimal longer in 13 of 13, p = 0.00024
optimal vs sub-optimal (medians): optimal longer in 13 of 13, p = 0.00024
optimal vs sub-optimal (means): optimal longer in 13 of 13, p = 0.00024
Median indeterminate dwell per mouse: {'BW041': 2, 'BW045': 4, 'BW046': 3, 'BW048': 3, 'BW051': 4, 'BW053': 4, 'BW058': 2, 'BW059': 4, 'BW061': 6, 'BW062': 3, 'BW064': 3, 'BW067': 1, 'BW071': 4}
mean +/- SEM over mice: 3.31 +/- 0.35 trials
Per-mouse medians (code dwell):
group disengaged optimal sub-optimal indeterminate
subject
BW041 13 25 11 2
BW045 13 17 10 4
BW046 11 28 13 3
BW048 9 27 11 3
BW051 13 24 17 4
BW053 16 72 17 4
BW058 11 55 16 2
BW059 8 16 11 4
BW061 32 106 91 6
BW062 14 16 11 3
BW064 8 12 11 3
BW067 6 40 8 1
BW071 28 131 21 4
Epochs per group: {'sub-optimal': 1920, 'optimal': 1534, 'disengaged': 876}
Epochs of dwell 0 (the extra last epoch): 3 of 4330
Without them: median indeterminate dwell 3.31 +/- 0.35 trials
(y) code dwell, (a) mean over mice of per-mouse median: {'disengaged': 14.0, 'optimal': 43.8, 'sub-optimal': 19.1}
(x) run length, (a) mean over mice of per-mouse median: {'disengaged': 13.8, 'optimal': 42.8, 'sub-optimal': 18.1}
Result: Partially reproduced. Optimal epochs last longer than disengaged and sub-optimal epochs in 13 of 13 mice under every reading (p = 0.00024), and the median indeterminate dwell, 3.31 +/- 0.35 trials, is within 3.4 +/- 0.3 (3.31 also without the 3 extra epochs of dwell 0). The optimal and sub-optimal dwells are longer than the paper's 37 and about 15 under every reading, while the disengaged dwell is within tolerance under readings (a) and (c): reading (a) gives 43.8 trials optimal, 14.0 disengaged and 19.1 sub-optimal with the code dwell, and 42.8, 13.8 and 18.1 as run length; the optimal dwell is above 40.7 and the sub-optimal above 16.5 in both, and readings (b) and (c) give 49 to 74 trials for the optimal state.
Choice prediction (claims fig2d-prediction-accuracy and fig1g-example-mouse-prediction)¶
The best-fit GLM-HMM predicted 75.7% +/- 1.7% of choices, a 10.8% +/- 1.3% improvement over a one-state (static performance) model.
As 6_predictive_accuracy.py: per mouse, 5 folds, each holding out a contiguous block of 20% of every session's trials, the order of the blocks permuted per session; on the training trials, concatenated as one sequence, a one-state model and a model with the mouse's number of states are fit (flex_hmm.choice_hmm_fit, MAP), 5 initialisations each, the best by test log likelihood kept; each test trial's choice is predicted from the state the transition matrix makes most likely after the filtered state of the preceding test trials.
rng_pred = np.random.default_rng(6)
t0 = time.time()
tasks, meta = [], []
folds = {}
for subj in subjects:
g = trials[trials.subject == subj].sort_values(["path", "trial_in_session"])
X = np.c_[g.inpt.values, np.ones(len(g))]
y = g.hmm_choice.values.astype(int)
sess_idx = pd.factorize(g.path)[0]
tr_i, te_i = hulsey.stratified_kfold_random_order(sess_idx, 5, rng_pred)
folds[subj] = (X, y, tr_i, te_i)
for f in range(5):
for K in (1, n_states[subj]):
for z in range(5):
tasks.append((K, X[tr_i[f]], y[tr_i[f]], X[te_i[f]], y[te_i[f]], int(rng_pred.integers(2**31))))
meta.append((subj, f, K))
res = hulsey.run_parallel(hulsey.fit_map_one, tasks)
print(f"{len(tasks)} fits in {(time.time() - t0) / 60:.1f} min")
acc = {}
for subj in subjects:
X, y, tr_i, te_i = folds[subj]
for K in (1, n_states[subj]):
a = []
for f in range(5):
cand = [res[j] for j in range(len(res)) if meta[j] == (subj, f, K)]
best = cand[int(np.argmax([c[1] for c in cand]))][0]
pred = hulsey.predict_choices(best, X[te_i[f]], y[te_i[f]])
a.append(np.mean(pred == y[te_i[f]]))
acc[(subj, "one state" if K == 1 else "best")] = np.mean(a) * 100
acc = pd.Series(acc).unstack()
acc["states"] = pd.Series(n_states)
acc["gain (points)"] = acc.best - acc["one state"]
acc["gain (relative %)"] = (acc.best / acc["one state"] - 1) * 100
print(acc.round(2).to_string())
print(f"\nbest model: {acc.best.mean():.1f} +/- {acc.best.sem():.1f}%; one state: {acc['one state'].mean():.1f} +/- {acc['one state'].sem():.1f}%")
print(f"(a) gain in points: {acc['gain (points)'].mean():.1f} +/- {acc['gain (points)'].sem():.1f}; "
f"(b) relative gain: {acc['gain (relative %)'].mean():.1f} +/- {acc['gain (relative %)'].sem():.1f}%")
print(f"best better in {int((acc['gain (points)'] > 0).sum())} of 13; Wilcoxon signed-rank p = {stats.wilcoxon(acc.best, acc['one state']).pvalue:.2g}; "
f"rank-sum (variant, O4) p = {stats.ranksums(acc.best, acc['one state']).pvalue:.2g}")
650 fits in 0.1 min
best one state states gain (points) gain (relative %) BW041 87.90 71.36 3 16.54 23.18 BW045 73.42 66.31 4 7.11 10.72 BW046 78.08 70.18 4 7.90 11.26 BW048 79.74 75.38 4 4.35 5.77 BW051 79.27 70.28 3 8.99 12.79 BW053 73.91 62.21 4 11.69 18.80 BW058 78.11 70.17 3 7.93 11.31 BW059 76.05 70.12 4 5.92 8.45 BW061 74.06 54.52 4 19.53 35.83 BW062 73.16 62.28 3 10.88 17.47 BW064 62.34 48.57 3 13.78 28.37 BW067 73.26 56.33 4 16.94 30.07 BW071 80.92 69.02 4 11.91 17.26 best model: 76.2 +/- 1.6%; one state: 65.1 +/- 2.2% (a) gain in points: 11.0 +/- 1.3; (b) relative gain: 17.8 +/- 2.5% best better in 13 of 13; Wilcoxon signed-rank p = 0.00024; rank-sum (variant, O4) p = 0.00025
Result: Reproduced. The selected GLM-HMM predicts 76.2% +/- 1.6% of choices against 65.1% for one state, a gain of 11.0 +/- 1.3 points (17.8% relative), better in 13 of 13 mice (signed-rank p = 0.00024; rank-sum p = 0.00025).
For the example mouse, a 4-state GLM-HMM was selected; it predicted 78% of choices against 70% for a one-state model.
The example mouse is BW046 (Figure1.py, 7_example_mouse_decoding.py).
print("BW046: states selected", n_states["BW046"], "(reading (a):", nst.loc["BW046", "(a) max"], ")")
print(f"BW046: one state {acc.loc['BW046', 'one state']:.1f}%, best model {acc.loc['BW046', 'best']:.1f}%")
BW046: states selected 4 (reading (a): 4 ) BW046: one state 70.2%, best model 78.1%
Result: Reproduced. BW046 has 4 states under both readings, and its choices are predicted at 78.1% by the 4-state model against 70.2% by the one-state model.
Pupil diameter and the optimal state (claims fig3c-inverted-u-pupil, fig3d-optimal-pupil-modality, fig3e-u-disengaged-pupil)¶
The probability of being in the optimal state follows an inverted-U function of pupil diameter, highest at intermediate diameters, in nearly all mice.
Sessions with more than 10 optimal-state trials and a post hoc pupil (claims file, M5). As 4_optimal_pupils.py: discrete-state trials with a pupil value are binned by pupil diameter in 0.02 bins (lower, upper]; per fold of a shuffled 5-fold split, P(optimal) per bin with more than 2 training or 3 test trials, polynomials of degree 1 to 7 fit to the training bins and scored by r2 on the test bins; the degree is the elbow of the mean r2 (largest height above the line from degree 1 to degree 7). As Figure3.py: P(state) per 0.02 bin [lower, upper) over trials with all six measures, for bins with more than 5 trials (a state absent from a bin is left missing, as the code does), the bins below the first one with a value merged into it; the polynomial of the selected degree is fit to P(optimal) and the optimal pupil is its maximum on a 0.001 grid within the range of the bins.
MEAS6 = ["post_hoc_pupil_diameter", "post_hoc_pupil_std10", "face_energy", "face_std_10", "running_speed", "running_std10"]
def fig3_set(subj, sess_col="sess_opt10"):
"""Figure3.py / FigureS3_S4.py / FigureS8.py trials: the mouse's sessions in the set, all six
measures present, the movement index z-scored over these trials."""
g = trials[(trials.subject == subj) & trials[sess_col]].dropna(subset=MEAS6).copy()
g["movement"] = stats.zscore(g.face_energy) + stats.zscore(g.running_speed)
g["movement_std"] = stats.zscore(g.face_std_10) + stats.zscore(g.running_std10)
return g
def p_state_bins(values, state, lo, hi, step, min_n=5, merge_low=True):
"""Figure3.py: P(optimal), P(disengaged), P(other) per bin [p, p + step) with more than min_n
trials; a state absent from a bin stays NaN; bins below the first with an optimal value are merged into it."""
st = np.where(state > 2, 2, state)
edges = np.arange(lo, hi, step)
out = np.full((len(edges), 3), np.nan)
for k, p in enumerate(edges):
sel = (values >= p) & (values < p + step)
n = sel.sum()
if n > min_n:
for s_ in np.unique(st[sel]):
out[k, int(s_)] = np.mean(st[sel] == s_)
if merge_low and np.isfinite(out[:, 0]).any():
k = np.where(np.isfinite(out[:, 0]))[0][0]
sel = (values >= lo) & (values < lo + (k + 1) * step)
if sel.sum() > min_n:
for s_ in np.unique(st[sel]):
out[k, int(s_)] = np.mean(st[sel] == s_)
return edges, out
def cv_r2(values, opt, edges_lower, step, rng, degrees=range(1, 8)):
"""4_optimal_pupils.py: 5-fold shuffled split of trials; binned P(optimal) in (lower, upper]
with more than 2 training and 3 test trials per bin; test r2 of polynomial fits per degree."""
from sklearn.model_selection import KFold
from sklearn.metrics import r2_score
kf = KFold(n_splits=5, shuffle=True, random_state=int(rng.integers(2**31)))
r2 = np.full((len(degrees), 5), np.nan)
def binp(v, o, mn):
y = np.full(len(edges_lower), np.nan)
for i, lo in enumerate(edges_lower):
sel = (v > lo) & (v <= lo + step)
if sel.sum() > mn:
y[i] = o[sel].mean()
return y
for f, (a, b) in enumerate(kf.split(values)):
ytr, yte = binp(values[a], opt[a], 2), binp(values[b], opt[b], 3)
x = np.asarray(edges_lower, dtype=np.float32)
xtr, ytr2 = x[np.isfinite(ytr)], ytr[np.isfinite(ytr)]
xte, yte2 = x[np.isfinite(yte)], yte[np.isfinite(yte)]
for j, d in enumerate(degrees):
c = np.polyfit(xtr, ytr2, d)
r2[j, f] = r2_score(yte2, np.poly1d(c)(xte))
return r2.mean(axis=1)
def elbow_degree(mean_r2):
"""4_optimal_pupils.py: the degree (2 to 6) whose mean r2 lies furthest above the line from degree 1 to 7."""
x1, y1, x2, y2 = 0, mean_r2[0], 6, mean_r2[6]
b = (x1 * y2 - x2 * y1) / (x1 - x2)
m = (y1 - y2) / (x1 - x2)
signs = [mean_r2[ii] - (m * ii + b) for ii in range(1, 6)]
return int(np.argmax(signs) + 2)
rng3 = np.random.default_rng(3)
pupil_rows, curves_pupil = [], {}
for subj in subjects:
g4 = trials[(trials.subject == subj) & trials.sess_opt10].dropna(subset=["post_hoc_pupil_diameter", "hmm_state"])
r2 = cv_r2(g4.post_hoc_pupil_diameter.values, (g4.hmm_state.values == 0).astype(float), np.arange(0, 1, .02), .02, rng3)
deg = elbow_degree(r2)
g = fig3_set(subj).dropna(subset=["hmm_state"])
edges, ps3 = p_state_bins(g.post_hoc_pupil_diameter.values, g.hmm_state.values, 0, 1, .02)
curves_pupil[subj] = (edges, ps3)
y = ps3[:, 0]
x2, y2 = edges[np.isfinite(y)], y[np.isfinite(y)]
xs = np.arange(x2[0], x2[-1], .001)
fit = np.poly1d(np.polyfit(x2.astype(np.float32), y2.astype(np.float32), deg))(xs)
opt_pupil = xs[np.argmax(fit)]
c2 = np.polyfit(x2, y2, 2)
pupil_rows.append({"subject": subj, "task": mods[subj], "trials": len(g), "bins": len(x2),
"range": f"{x2[0]:.2f}-{x2[-1]:.2f}", "degree": deg, "r2 deg2": r2[1],
"optimal pupil": opt_pupil, "fit at max": fit.max(), "fit at low end": fit[0], "fit at high end": fit[-1],
"(i) interior max, lower ends": bool(0 < np.argmax(fit) < len(fit) - 1 and fit[0] < fit.max() and fit[-1] < fit.max()),
"(ii) quad coef": c2[0]})
pupil_tab = pd.DataFrame(pupil_rows).set_index("subject")
print(pupil_tab.round(3).to_string())
n_i = int(pupil_tab["(i) interior max, lower ends"].sum())
n_ii = int((pupil_tab["(ii) quad coef"] < 0).sum())
print(f"\n(i) maximum inside the range with lower values at both ends: {n_i} of 13; (ii) negative quadratic term: {n_ii} of 13")
print("Readings of 'nearly all': (a) at least 12 of 13, (b) at least 11 of 13:",
{"(i)": (n_i >= 12, n_i >= 11), "(ii)": (n_ii >= 12, n_ii >= 11)})
optimal_pupil = pupil_tab["optimal pupil"]
task trials bins range degree r2 deg2 optimal pupil fit at max fit at low end fit at high end (i) interior max, lower ends (ii) quad coef
subject
BW041 visual 10470 24 0.50-0.96 2 0.700 0.699 0.751 0.379 0.111 True -9.430
BW045 auditory 14768 30 0.38-0.96 3 0.712 0.562 0.569 0.161 0.120 True -4.043
BW046 visual 7345 19 0.60-0.98 2 0.715 0.807 0.661 -0.004 0.197 True -15.597
BW048 visual 11888 23 0.54-0.98 2 0.749 0.680 0.824 0.660 0.074 True -8.384
BW051 visual 9419 22 0.56-0.98 2 0.403 0.801 0.748 0.275 0.489 True -8.162
BW053 auditory 8148 33 0.02-0.94 2 0.751 0.385 0.681 0.340 -0.103 True -2.558
BW058 visual 4542 21 0.52-0.94 2 0.532 0.739 0.785 0.141 0.245 True -13.388
BW059 auditory 8520 22 0.54-0.96 2 0.788 0.689 0.752 0.475 -0.156 True -12.472
BW061 auditory 7021 21 0.56-0.96 3 0.570 0.663 0.671 0.430 0.152 True -6.092
BW062 auditory 6224 24 0.50-0.96 2 0.871 0.710 0.753 0.192 -0.042 True -12.787
BW064 auditory 4718 20 0.60-0.98 6 0.797 0.635 0.754 0.680 0.196 True -0.507
BW067 auditory 10604 22 0.52-0.94 3 0.853 0.618 0.775 0.466 0.046 True -6.260
BW071 auditory 8047 21 0.56-0.96 3 0.534 0.662 0.872 0.318 0.274 True -7.018
(i) maximum inside the range with lower values at both ends: 13 of 13; (ii) negative quadratic term: 13 of 13
Readings of 'nearly all': (a) at least 12 of 13, (b) at least 11 of 13: {'(i)': (True, True), '(ii)': (True, True)}
fig, axes = plt.subplots(2, 7, figsize=(18, 5), sharey=True)
for ax, subj in zip(axes.ravel(), subjects):
edges, ps3 = curves_pupil[subj]
for k, c in zip(range(3), ["teal", "firebrick", "orange"]):
ax.plot(edges, ps3[:, k], ".-", color=c, lw=1)
ax.axvline(optimal_pupil[subj], color="k", ls=":")
ax.set_title(subj, fontsize=9)
ax.set_xlim(0, 1)
axes.ravel()[-1].axis("off")
fig.supxlabel("pupil diameter (fraction of session maximum)")
fig.supylabel("P(state): optimal teal, disengaged red, other orange")
plt.tight_layout()
plt.show()
print("All plotted values lie in [0, 1] on both axes; dotted line: optimal pupil.")
All plotted values lie in [0, 1] on both axes; dotted line: optimal pupil.
Result: Reproduced. In 13 of 13 mice the fit of P(optimal) on pupil diameter has its maximum inside the data range with lower values at both ends, and in 13 of 13 a second-order fit has a negative quadratic term, so "nearly all" holds under every reading.
The optimal pupil diameter is smaller in mice performing the auditory task than in those performing the visual task.
aud = optimal_pupil[pupil_tab.task == "auditory"]
vis = optimal_pupil[pupil_tab.task == "visual"]
rs = stats.ranksums(aud, vis)
print("auditory:", aud.round(3).to_dict())
print("visual:", vis.round(3).to_dict())
print(f"median auditory {aud.median():.3f}, visual {vis.median():.3f}; Wilcoxon rank-sum p = {rs.pvalue:.3g}; "
f"Mann-Whitney U exact p = {stats.mannwhitneyu(aud, vis).pvalue:.3g}")
auditory: {'BW045': 0.562, 'BW053': 0.385, 'BW059': 0.689, 'BW061': 0.663, 'BW062': 0.71, 'BW064': 0.635, 'BW067': 0.618, 'BW071': 0.662}
visual: {'BW041': 0.699, 'BW046': 0.807, 'BW048': 0.68, 'BW051': 0.801, 'BW058': 0.739}
median auditory 0.649, visual 0.739; Wilcoxon rank-sum p = 0.0128; Mann-Whitney U exact p = 0.0109
Result: Reproduced. The optimal pupil is smaller in auditory-task mice (median 0.649 of the session maximum) than in visual-task mice (0.739); Wilcoxon rank-sum p = 0.013 (exact Mann-Whitney p = 0.011).
Across mice, aligned to each mouse's optimal pupil, the probability of the optimal state is an inverted U and the probability of the disengaged state a U of pupil diameter; this holds also on trials when the mice were stationary.
As Figure3.py: each mouse's P(state) per 0.02 bin shifted so that the bin nearest its optimal pupil is at 0, and averaged over mice. The shape test is the claims file's: on bins with data from at least 7 of 13 mice, a second-order fit with a negative (optimal) or positive (disengaged) quadratic term and its vertex inside the range. Stationary trials, readings (a) locomotion speed exactly 0 and (b) absolute speed below 0.02 m/s, with the same optimal pupils; (b') speed below 0.02 m/s with backward movement counted as stationary, the rule of the first version of this notebook, is shown for comparison.
def shifted_mean(sel_fn):
shifted = np.full((len(subjects), 3, 50), np.nan)
x = np.arange(0, 1, .02)
for i, subj in enumerate(subjects):
g = fig3_set(subj).dropna(subset=["hmm_state"])
g = g[sel_fn(g)]
if len(g) == 0:
continue
_, ps3 = p_state_bins(g.post_hoc_pupil_diameter.values, g.hmm_state.values, 0, 1, .02)
ob = np.argmin(np.abs(x - optimal_pupil[subj]))
for j, b in enumerate(range(ob - 25, ob + 25)):
if 0 <= b < 50:
shifted[i, :, j] = ps3[b]
return shifted
def shape(curve_by_mouse, k, sign):
n = np.sum(np.isfinite(curve_by_mouse[:, k, :]), axis=0)
mean = np.nanmean(curve_by_mouse[:, k, :], axis=0)
keep = n >= 7
xb = np.arange(-50, 50, 2)[keep]
if keep.sum() < 3:
return dict(bins=int(keep.sum()), quad=np.nan, vertex=np.nan, holds=False)
c = np.polyfit(xb, mean[keep], 2)
vx = -c[1] / (2 * c[0])
return dict(bins=int(keep.sum()), range=f"{xb[0]} to {xb[-1]}", quad=c[0], vertex=round(vx, 1),
holds=bool(np.sign(c[0]) == sign and xb[0] < vx < xb[-1]))
shift_res = {}
for name, fn in [("all trials", lambda g: np.ones(len(g), bool)),
("(a) stationary, speed == 0", lambda g: g.running_speed.values == 0),
("(b) stationary, |speed| < 0.02 m/s", lambda g: np.abs(g.running_speed.values) < 0.02),
("(b') speed < 0.02 m/s, backward included", lambda g: g.running_speed.values < 0.02)]:
sh = shifted_mean(fn)
shift_res[name] = sh
o, d = shape(sh, 0, -1), shape(sh, 1, +1)
nt = sum(int(fn(fig3_set(s_).dropna(subset=["hmm_state"])).sum()) for s_ in subjects)
print(f"{name} ({nt} trials): optimal {o}\n disengaged {d}")
nback = sum(int((fig3_set(s_).dropna(subset=["hmm_state"]).running_speed.values <= -0.02).sum()) for s_ in subjects)
print(f"Trials with speed at or below -0.02 m/s (backward), counted as stationary by (b') only: {nback}")
nnear = sum(int(((v := fig3_set(s_).dropna(subset=["hmm_state"]).running_speed.abs().values) > 0).sum() - (v >= 1e-9).sum()) for s_ in subjects)
print(f"Trials with speed not exactly 0 but below 1e-9 m/s in absolute value (zero up to rounding in the filter): {nnear}")
fig, ax = plt.subplots(1, 4, figsize=(18, 3.5), sharey=True)
for a, (name, sh) in zip(ax, shift_res.items()):
for k, c in zip(range(3), ["teal", "firebrick", "orange"]):
a.plot(np.arange(-50, 50, 2), np.nanmean(sh[:, k, :], axis=0), color=c)
a.set_title(name, fontsize=9)
a.set_xlabel("pupil minus optimal pupil (% of max)")
ax[0].set_ylabel("mean P(state) over mice")
plt.tight_layout()
plt.show()
print("x from -50 to 48 % of max; y in [0, 1].")
all trials (111714 trials): optimal {'bins': 23, 'range': '-14 to 30', 'quad': np.float64(-0.0006669451053081474), 'vertex': np.float64(-0.9), 'holds': True}
disengaged {'bins': 23, 'range': '-14 to 30', 'quad': np.float64(0.000523755085708531), 'vertex': np.float64(-1.6), 'holds': True}
(a) stationary, speed == 0 (1821 trials): optimal {'bins': 0, 'quad': nan, 'vertex': nan, 'holds': False}
disengaged {'bins': 0, 'quad': nan, 'vertex': nan, 'holds': False}
(b) stationary, |speed| < 0.02 m/s (32495 trials): optimal {'bins': 21, 'range': '-14 to 26', 'quad': np.float64(-0.0004066500841704846), 'vertex': np.float64(-3.3), 'holds': True}
disengaged {'bins': 21, 'range': '-14 to 26', 'quad': np.float64(0.00030294345210026875), 'vertex': np.float64(-6.1), 'holds': True}
(b') speed < 0.02 m/s, backward included (32709 trials): optimal {'bins': 21, 'range': '-14 to 26', 'quad': np.float64(-0.0004047777523347985), 'vertex': np.float64(-3.3), 'holds': True}
disengaged {'bins': 21, 'range': '-14 to 26', 'quad': np.float64(0.0003036137994277513), 'vertex': np.float64(-6.0), 'holds': True}
Trials with speed at or below -0.02 m/s (backward), counted as stationary by (b') only: 214
Trials with speed not exactly 0 but below 1e-9 m/s in absolute value (zero up to rounding in the filter): 9869
x from -50 to 48 % of max; y in [0, 1].
Result: Reproduced. Aligned to each mouse's optimal pupil and averaged, P(optimal) has a negative and P(disengaged) a positive quadratic term with the vertex inside the range (23 bins with at least 7 mice), and both shapes hold on stationary trials with absolute speed below 0.02 m/s (reading (b), 32,495 trials, 21 bins; vertices -3.3 and -6.1). Counting the 214 backward trials at or below -0.02 m/s as stationary (b') changes the vertices by at most 0.1. Reading (a), speed exactly 0, leaves 1,821 trials and no bin with 7 mice, so the shapes cannot be tested under it. A speed of exactly 0 depends on rounding in the filter: 9,869 more trials have a speed below 1e-9 m/s that is not exactly 0, and a run on another machine counted 1,879 trials under (a). Readings (b) and (b') do not depend on this.
Movement measures and the optimal state (claims figs4-r2-individual-movement, fig3f-movement-index-inverted-u, fig3f-movement-index-u-disengaged, figs5-movement-index-correlations)¶
Second-order cross-validated regression fits of optimal-state probability are poorer for face motion energy (r2 0.49 +/- 0.08) and locomotion speed (0.44 +/- 0.07) than for pupil diameter (0.72 +/- 0.03).
The cross-validated r2 of a degree-2 fit as 4_optimal_pupils.py computes it (above), for pupil on its own trial set, and for face motion energy, locomotion speed and the movement index on the Figure3.py trial set (sessions of M5, six measures present). Bins, O5: (a) the plotting bins of FigureS3_S4.py and Figure3.py (face 0.05 from 0 to 1, locomotion 0.02 m/s from -0.02 to 0.30, movement index 0.25 from -5 to 5); (b) 50 equal-width bins over each mouse's range.
rng4 = np.random.default_rng(4)
BINS_A = {"face_energy": (0, 1, .05), "running_speed": (-.02, .301, .02), "movement": (-5, 5, .25)}
r2rows = []
for subj in subjects:
g = fig3_set(subj).dropna(subset=["hmm_state"])
opt = (g.hmm_state.values == 0).astype(float)
row = {"subject": subj, "pupil": pupil_tab.loc[subj, "r2 deg2"]}
for meas, (lo, hi, step) in BINS_A.items():
v = g[meas].values
row[f"{meas} (a)"] = cv_r2(v, opt, np.arange(lo, hi, step), step, rng4, degrees=[1, 2])[1]
stepb = (v.max() - v.min()) / 50
row[f"{meas} (b)"] = cv_r2(v, opt, v.min() - 1e-9 + np.arange(50) * stepb, stepb, rng4, degrees=[1, 2])[1]
r2rows.append(row)
r2tab = pd.DataFrame(r2rows).set_index("subject")
print(r2tab.round(3).to_string())
print("\nmean +/- SEM over mice:", {c: f"{r2tab[c].mean():.4f} +/- {r2tab[c].sem():.4f}" for c in r2tab.columns})
for c in r2tab.columns[1:]:
w = stats.wilcoxon(r2tab[c], r2tab.pupil)
print(f"{c} vs pupil: difference {np.mean(r2tab[c] - r2tab.pupil):+.2f}, lower in {int((r2tab[c] < r2tab.pupil).sum())} of 13, p = {w.pvalue:.3g}")
pupil face_energy (a) face_energy (b) running_speed (a) running_speed (b) movement (a) movement (b)
subject
BW041 0.700 0.535 0.726 -0.169 0.250 0.701 0.575
BW045 0.712 0.628 0.355 0.794 0.681 0.480 0.375
BW046 0.715 -0.108 0.186 0.529 0.214 0.323 0.147
BW048 0.749 0.481 0.124 0.781 0.351 0.386 0.113
BW051 0.403 0.701 0.650 0.789 0.477 0.836 0.756
BW053 0.751 0.234 0.351 0.411 -0.090 0.393 0.458
BW058 0.532 0.847 0.819 0.623 0.252 0.747 0.760
BW059 0.788 0.070 0.027 0.553 0.469 0.260 0.627
BW061 0.570 0.717 0.738 0.753 0.770 0.798 0.782
BW062 0.871 0.119 0.078 0.843 0.719 0.767 0.693
BW064 0.797 -0.053 -0.063 0.139 0.461 0.386 0.284
BW067 0.853 0.730 0.733 0.431 0.197 0.669 0.670
BW071 0.534 0.669 0.598 0.899 0.616 0.713 0.694
mean +/- SEM over mice: {'pupil': '0.6904 +/- 0.0388', 'face_energy (a)': '0.4284 +/- 0.0919', 'face_energy (b)': '0.4093 +/- 0.0874', 'running_speed (a)': '0.5673 +/- 0.0855', 'running_speed (b)': '0.4127 +/- 0.0688', 'movement (a)': '0.5737 +/- 0.0569', 'movement (b)': '0.5334 +/- 0.0652'}
face_energy (a) vs pupil: difference -0.26, lower in 9 of 13, p = 0.11
face_energy (b) vs pupil: difference -0.28, lower in 8 of 13, p = 0.0574
running_speed (a) vs pupil: difference -0.12, lower in 7 of 13, p = 0.414
running_speed (b) vs pupil: difference -0.28, lower in 10 of 13, p = 0.0105
movement (a) vs pupil: difference -0.12, lower in 8 of 13, p = 0.216
movement (b) vs pupil: difference -0.16, lower in 9 of 13, p = 0.168
# Variant: the r2 above comes from one shuffled 5-fold split per mouse, as in 4_optimal_pupils.py.
# Mean over 10 splits per mouse, to show how much the split moves it.
rng4b = np.random.default_rng(41)
N_SPLITS = 10
rep_rows = []
for subj in subjects:
g4 = trials[(trials.subject == subj) & trials.sess_opt10].dropna(subset=["post_hoc_pupil_diameter", "hmm_state"])
o4 = (g4.hmm_state.values == 0).astype(float)
pr = [cv_r2(g4.post_hoc_pupil_diameter.values, o4, np.arange(0, 1, .02), .02, rng4b, degrees=[1, 2])[1] for _ in range(N_SPLITS)]
row = {"subject": subj, "pupil": np.mean(pr), "pupil min": np.min(pr), "pupil max": np.max(pr)}
g = fig3_set(subj).dropna(subset=["hmm_state"])
opt = (g.hmm_state.values == 0).astype(float)
for meas, (lo, hi, step) in BINS_A.items():
v = g[meas].values
row[f"{meas} (a)"] = np.mean([cv_r2(v, opt, np.arange(lo, hi, step), step, rng4b, degrees=[1, 2])[1] for _ in range(N_SPLITS)])
stepb = (v.max() - v.min()) / 50
row[f"{meas} (b)"] = np.mean([cv_r2(v, opt, v.min() - 1e-9 + np.arange(50) * stepb, stepb, rng4b, degrees=[1, 2])[1] for _ in range(N_SPLITS)])
rep_rows.append(row)
r2rep = pd.DataFrame(rep_rows).set_index("subject")
print(f"Mean r2 over {N_SPLITS} shuffled splits per mouse (pupil min and max over the splits):")
print(r2rep.round(3).to_string())
print("\nmean +/- SEM over mice:", {c: f"{r2rep[c].mean():.4f} +/- {r2rep[c].sem():.4f}" for c in r2rep.columns if c not in ("pupil min", "pupil max")})
for c in [c for c in r2rep.columns if "(" in c]:
w = stats.wilcoxon(r2rep[c], r2rep.pupil)
print(f"{c} vs pupil: difference {np.mean(r2rep[c] - r2rep.pupil):+.2f}, lower in {int((r2rep[c] < r2rep.pupil).sum())} of 13, p = {w.pvalue:.3g}")
Mean r2 over 10 shuffled splits per mouse (pupil min and max over the splits):
pupil pupil min pupil max face_energy (a) face_energy (b) running_speed (a) running_speed (b) movement (a) movement (b)
subject
BW041 0.696 0.652 0.744 0.556 0.664 -0.078 0.238 0.684 0.541
BW045 0.679 0.611 0.706 0.521 0.383 0.788 0.692 0.463 0.446
BW046 0.710 0.643 0.790 0.438 0.189 0.496 0.216 0.208 0.068
BW048 0.774 0.720 0.819 0.463 0.269 0.772 0.357 0.361 0.125
BW051 0.342 0.001 0.562 0.735 0.670 0.794 0.530 0.837 0.773
BW053 0.737 0.694 0.755 0.226 0.355 0.322 -0.097 0.419 0.478
BW058 0.478 0.382 0.548 0.875 0.838 0.628 0.247 0.765 0.759
BW059 0.783 0.754 0.801 -0.228 0.017 0.598 0.487 0.382 0.516
BW061 0.636 0.604 0.651 0.757 0.738 0.771 0.760 0.798 0.782
BW062 0.877 0.867 0.890 0.165 -0.116 0.837 0.716 0.764 0.676
BW064 0.803 0.721 0.848 -0.032 -0.047 0.226 0.477 0.390 0.271
BW067 0.839 0.817 0.853 0.740 0.748 0.404 0.154 0.663 0.645
BW071 0.507 0.470 0.549 0.696 0.598 0.867 0.657 0.718 0.700
mean +/- SEM over mice: {'pupil': '0.6816 +/- 0.0431', 'face_energy (a)': '0.4547 +/- 0.0927', 'face_energy (b)': '0.4083 +/- 0.0907', 'running_speed (a)': '0.5711 +/- 0.0791', 'running_speed (b)': '0.4180 +/- 0.0716', 'movement (a)': '0.5733 +/- 0.0577', 'movement (b)': '0.5216 +/- 0.0666'}
face_energy (a) vs pupil: difference -0.23, lower in 9 of 13, p = 0.146
face_energy (b) vs pupil: difference -0.27, lower in 9 of 13, p = 0.0681
running_speed (a) vs pupil: difference -0.11, lower in 8 of 13, p = 0.34
running_speed (b) vs pupil: difference -0.26, lower in 9 of 13, p = 0.0134
movement (a) vs pupil: difference -0.11, lower in 9 of 13, p = 0.216
movement (b) vs pupil: difference -0.16, lower in 9 of 13, p = 0.127
Result: Partially reproduced. Pupil diameter gives 0.690 +/- 0.039, within 0.72 +/- 0.03. Face motion energy gives 0.428 +/- 0.092 with the plotting bins (a), within tolerance, and 0.409 with 50 bins (b), just outside; locomotion speed gives 0.567 under (a), outside (0.570 in a run on another machine, where the speeds at zero round differently, see fig3e), and 0.413 under (b), within. So no bin reading puts all three values in tolerance. Both measures are below pupil on average, but face is not lower at p < 0.05 under either reading (p = 0.11 and 0.057), and locomotion only under (b) (p = 0.011; (a) p = 0.41). Over 10 shuffled splits per mouse, pupil gives 0.682 +/- 0.043 (BW051 from 0.001 to 0.562), and the tests come out the same way (face p = 0.15 and 0.068, locomotion p = 0.34 and 0.013); the single-split r2 of a mouse can move by more than 0.1 with the split.
Optimal-state probability peaks at intermediate values of the movement index, an inverted U similar to pupil diameter (cross-validated r2 0.63 +/- 0.06, not different from pupil).
The movement index is the sum of the z-scored face motion energy and locomotion speed over the trials of the Figure3.py set (M4). Its r2 and the test against pupil are in the table above (rows movement). The shape is tested on the mean P(optimal) curve over the 0.25 bins of Figure3.py, as for pupil.
mv = np.full((len(subjects), 3, 40), np.nan)
for i, subj in enumerate(subjects):
g = fig3_set(subj).dropna(subset=["hmm_state"])
_, ps3 = p_state_bins(g.movement.values, g.hmm_state.values, -5, 5, .25, merge_low=False)
mv[i] = ps3.T
def shape_x(curves, k, sign, xb_all):
n = np.sum(np.isfinite(curves[:, k, :]), axis=0)
mean = np.nanmean(curves[:, k, :], axis=0)
keep = n >= 7
xb = xb_all[keep]
c = np.polyfit(xb, mean[keep], 2)
vx = -c[1] / (2 * c[0])
return dict(bins=int(keep.sum()), range=f"{xb[0]} to {xb[-1]}", quad=round(c[0], 4), vertex=round(vx, 2),
peak_bin_of_mean=float(xb[np.argmax(mean[keep])] if sign < 0 else xb[np.argmin(mean[keep])]),
holds=bool(np.sign(c[0]) == sign and xb[0] < vx < xb[-1]))
xmv = np.arange(-5, 5, .25)
print("optimal:", shape_x(mv, 0, -1, xmv))
for c in ["movement (a)", "movement (b)"]:
print(f"{c}: r2 {r2tab[c].mean():.4f} +/- {r2tab[c].sem():.3f}; vs pupil difference {np.mean(r2tab[c] - r2tab.pupil):+.3f}, "
f"p = {stats.wilcoxon(r2tab[c], r2tab.pupil).pvalue:.3g}")
fig, ax = plt.subplots(figsize=(5, 3.5))
for k, c in zip(range(3), ["teal", "firebrick", "orange"]):
ax.plot(xmv, np.nanmean(mv[:, k, :], axis=0), color=c)
ax.set_xlabel("movement index (bins of 0.25)")
ax.set_ylabel("mean P(state) over mice")
plt.tight_layout()
plt.show()
print("x from -5 to 4.75; y in [0, 1]; bins with fewer than 7 mice are included in the plot but not in the fit.")
optimal: {'bins': 33, 'range': '-4.25 to 3.75', 'quad': np.float64(-0.0239), 'vertex': np.float64(0.28), 'peak_bin_of_mean': 0.0, 'holds': True}
movement (a): r2 0.5737 +/- 0.057; vs pupil difference -0.117, p = 0.216
movement (b): r2 0.5334 +/- 0.065; vs pupil difference -0.157, p = 0.168
x from -5 to 4.75; y in [0, 1]; bins with fewer than 7 mice are included in the plot but not in the fit.
Result: Reproduced. The mean P(optimal) over movement-index bins is an inverted U (negative quadratic term, vertex at 0.28, inside -4.25 to 3.75). Its r2 is 0.574 +/- 0.057 with the plotting bins (a), within 0.63 +/- 0.06, and 0.533 with 50 bins (b), outside; it does not differ from pupil's under either reading (p = 0.22 and 0.17). The (a) value is close to the lower bound of 0.57: over 10 shuffled splits per mouse it is 0.573 (b: 0.522).
The probability of the disengaged state is a U-shaped function of the movement index.
print("disengaged:", shape_x(mv, 1, +1, xmv))
disengaged: {'bins': 34, 'range': '-4.25 to 4.0', 'quad': np.float64(0.0225), 'vertex': np.float64(0.08), 'peak_bin_of_mean': -0.25, 'holds': True}
Result: Reproduced. The mean P(disengaged) over movement-index bins has a positive quadratic term with its vertex (0.08) inside the range of bins with at least 7 mice (-4.25 to 4.0).
The movement index is equally correlated with face motion energy and with locomotion speed, and highly correlated with both in each mouse.
cr = []
for subj in subjects:
g = fig3_set(subj)
cr.append({"subject": subj, "r face": stats.pearsonr(g.movement, g.face_energy)[0],
"r locomotion": stats.pearsonr(g.movement, g.running_speed)[0]})
cr = pd.DataFrame(cr).set_index("subject")
cr["difference"] = cr["r face"] - cr["r locomotion"]
print(cr.round(4).to_string())
print(f"\nlargest |difference| {cr.difference.abs().max():.4f}; smallest r {cr[['r face', 'r locomotion']].min().min():.3f}; "
f"(a) both above 0.7 in every mouse: {bool((cr[['r face', 'r locomotion']] > .7).all().all())}; "
f"(b) above 0.5: {bool((cr[['r face', 'r locomotion']] > .5).all().all())}")
r face r locomotion difference subject BW041 0.7441 0.7441 -0.0 BW045 0.7862 0.7862 0.0 BW046 0.7310 0.7310 -0.0 BW048 0.7307 0.7307 0.0 BW051 0.8319 0.8319 -0.0 BW053 0.8047 0.8047 0.0 BW058 0.8547 0.8547 0.0 BW059 0.7353 0.7353 -0.0 BW061 0.8221 0.8221 -0.0 BW062 0.7942 0.7942 0.0 BW064 0.7373 0.7373 0.0 BW067 0.7605 0.7605 0.0 BW071 0.6737 0.6737 -0.0 largest |difference| 0.0000; smallest r 0.674; (a) both above 0.7 in every mouse: False; (b) above 0.5: True
Result: Reproduced. In every mouse the movement index correlates equally with face motion energy and with locomotion speed (difference below 0.0001), with r from 0.674 to 0.855, so reading (b), both above 0.5, holds; reading (a), above 0.7, fails in one mouse (BW071, 0.674).
Raw measures by state (claims fig3-distance-from-optimal-pupil, fig3-raw-face-energy-trend, fig3-raw-locomotion-no-difference)¶
Per mouse, the mean at stimulus onset over the trials of each state, sub-optimal states pooled; sessions read both ways of O3: (a) sessions with more than 10 optimal-state trials and a post hoc pupil (Figure3.py), (b) all included sessions. Each measure's trials are those where it is present.
def by_state(col, sess_col, fn=lambda g, c: g[c]):
rows = []
for subj in subjects:
g = trials[(trials.subject == subj)]
if sess_col:
g = g[g[sess_col]]
v = fn(g, col)
st = g.hmm_state
rows.append({"subject": subj, "optimal": v[st == 0].mean(), "disengaged": v[st == 1].mean(), "sub-optimal": v[st > 1].mean()})
return pd.DataFrame(rows).set_index("subject")
def report(tab, scale, larger=None):
m_, se_ = tab.mean() * scale, tab.sem() * scale
print(" mean +/- SEM:", {k: f"{m_[k]:.3f} +/- {se_[k]:.3f}" for k in tab.columns})
for other in ["disengaged", "sub-optimal"]:
w = stats.wilcoxon(tab.optimal, tab[other])
print(f" optimal vs {other}: difference {np.mean(tab.optimal - tab[other]) * scale:+.3f}, "
f"optimal lower in {int((tab.optimal < tab[other]).sum())} of 13, p = {w.pvalue:.2g}")
dist = lambda g, c: (g[c] - g.subject.map(optimal_pupil)).abs()
for label, sc in [("(a) sessions with >10 optimal trials and pupil", "sess_opt10"), ("(b) all included sessions", None)]:
print("Absolute difference from the optimal pupil (% of session maximum),", label)
report(by_state("post_hoc_pupil_diameter", sc, dist), 100)
Absolute difference from the optimal pupil (% of session maximum), (a) sessions with >10 optimal trials and pupil
mean +/- SEM: {'optimal': '8.804 +/- 0.964', 'disengaged': '15.955 +/- 1.506', 'sub-optimal': '11.795 +/- 1.298'}
optimal vs disengaged: difference -7.152, optimal lower in 13 of 13, p = 0.00024
optimal vs sub-optimal: difference -2.991, optimal lower in 13 of 13, p = 0.00024
Absolute difference from the optimal pupil (% of session maximum), (b) all included sessions
mean +/- SEM: {'optimal': '8.806 +/- 0.964', 'disengaged': '15.830 +/- 1.484', 'sub-optimal': '11.450 +/- 1.182'}
optimal vs disengaged: difference -7.024, optimal lower in 13 of 13, p = 0.00024
optimal vs sub-optimal: difference -2.644, optimal lower in 12 of 13, p = 0.0012
Result: Partially reproduced. The distance from the optimal pupil is smallest in the optimal state and lower than in the disengaged (13 of 13 mice, p = 0.00024 under both session readings) and sub-optimal states (p = 0.00024 under (a), 0.0012 under (b)), but the values are larger than the paper's: 8.8%, 16.0% and 11.8% of the session maximum under (a) against 7%, 12.6% and 9%.
for label, sc in [("(a)", "sess_opt10"), ("(b)", None)]:
print("Face motion energy (a.u.)", label)
report(by_state("face_energy", sc), 1)
Face motion energy (a.u.) (a)
mean +/- SEM: {'optimal': '0.482 +/- 0.024', 'disengaged': '0.441 +/- 0.030', 'sub-optimal': '0.448 +/- 0.025'}
optimal vs disengaged: difference +0.041, optimal lower in 5 of 13, p = 0.19
optimal vs sub-optimal: difference +0.034, optimal lower in 2 of 13, p = 0.017
Face motion energy (a.u.) (b)
mean +/- SEM: {'optimal': '0.485 +/- 0.024', 'disengaged': '0.435 +/- 0.029', 'sub-optimal': '0.453 +/- 0.027'}
optimal vs disengaged: difference +0.050, optimal lower in 4 of 13, p = 0.08
optimal vs sub-optimal: difference +0.032, optimal lower in 3 of 13, p = 0.048
Result: Reproduced. Under reading (a), face motion energy is 0.482 +/- 0.024 a.u. in the optimal state, 0.441 in the disengaged and 0.448 in the sub-optimal states, all within tolerance; the optimal mean is the largest, above sub-optimal with p = 0.017, and above disengaged with p = 0.19 (paper 0.06). Under reading (b) the values are 0.485, 0.435 and 0.453 (p = 0.048 and 0.08).
for label, sc in [("(a)", "sess_opt10"), ("(b)", None)]:
print("Locomotion speed (m/s)", label)
report(by_state("running_speed", sc), 1)
Locomotion speed (m/s) (a)
mean +/- SEM: {'optimal': '0.136 +/- 0.024', 'disengaged': '0.137 +/- 0.026', 'sub-optimal': '0.136 +/- 0.024'}
optimal vs disengaged: difference -0.002, optimal lower in 5 of 13, p = 0.84
optimal vs sub-optimal: difference -0.001, optimal lower in 5 of 13, p = 0.89
Locomotion speed (m/s) (b)
mean +/- SEM: {'optimal': '0.137 +/- 0.023', 'disengaged': '0.132 +/- 0.024', 'sub-optimal': '0.139 +/- 0.024'}
optimal vs disengaged: difference +0.005, optimal lower in 5 of 13, p = 0.68
optimal vs sub-optimal: difference -0.002, optimal lower in 6 of 13, p = 0.95
Result: Reproduced. Locomotion speed is 0.136, 0.137 and 0.136 m/s in the optimal, disengaged and sub-optimal states (reading (a)), and neither test finds a difference (p = 0.84 and 0.89; reading (b) p = 0.68 and 0.95).
Variability of the measures (claims fig4a-lower-variability, fig4b-variability-at-transitions)¶
The 10-trial variability of pupil diameter, face motion energy and locomotion speed is lower in the optimal state than in the disengaged and sub-optimal states.
As Figure4.py: all included sessions, trials with the three 10-trial standard deviations and the pupil present; per mouse the mean of each standard deviation over the trials of each state, sub-optimal states pooled.
STD3 = ["post_hoc_pupil_std10", "face_std_10", "running_std10"]
fig4 = {}
for subj in subjects:
g = trials[trials.subject == subj].dropna(subset=STD3 + ["post_hoc_pupil_diameter"]).copy()
g["movement"] = stats.zscore(g.face_energy) + stats.zscore(g.running_speed)
fig4[subj] = g
n_tests = 0
for meas in STD3:
tab = pd.DataFrame({subj: {"optimal": g.loc[g.hmm_state == 0, meas].mean(), "disengaged": g.loc[g.hmm_state == 1, meas].mean(),
"sub-optimal": g.loc[g.hmm_state > 1, meas].mean()} for subj, g in fig4.items()}).T
print(meas, "per mouse:")
print(tab.round(4).T.to_string())
print(meas, "mean over mice:", tab.mean().round(4).to_dict())
for other in ["disengaged", "sub-optimal"]:
w = stats.wilcoxon(tab.optimal, tab[other])
n_tests += (w.pvalue < 0.05) and (tab.optimal.mean() < tab[other].mean())
print(f" optimal vs {other}: difference {np.mean(tab.optimal - tab[other]):+.4f}, optimal lower in "
f"{int((tab.optimal < tab[other]).sum())} of 13, p = {w.pvalue:.2g}")
print("tests with the optimal state lower at p < 0.05:", int(n_tests), "of 6")
post_hoc_pupil_std10 per mouse:
BW041 BW045 BW046 BW048 BW051 BW053 BW058 BW059 BW061 BW062 BW064 BW067 BW071
optimal 0.0302 0.0462 0.0217 0.0242 0.0286 0.0482 0.0241 0.0177 0.0188 0.0281 0.0275 0.0161 0.0186
disengaged 0.0558 0.0501 0.0301 0.0499 0.0463 0.0953 0.0702 0.0308 0.0288 0.0288 0.0317 0.0256 0.0457
sub-optimal 0.0422 0.0468 0.0236 0.0257 0.0358 0.0622 0.0622 0.0205 0.0240 0.0267 0.0283 0.0189 0.0259
post_hoc_pupil_std10 mean over mice: {'optimal': 0.0269, 'disengaged': 0.0453, 'sub-optimal': 0.0341}
optimal vs disengaged: difference -0.0184, optimal lower in 13 of 13, p = 0.00024
optimal vs sub-optimal: difference -0.0071, optimal lower in 12 of 13, p = 0.0012
face_std_10 per mouse:
BW041 BW045 BW046 BW048 BW051 BW053 BW058 BW059 BW061 BW062 BW064 BW067 BW071
optimal 0.0612 0.1007 0.0435 0.0425 0.0700 0.0991 0.0669 0.0763 0.0571 0.0733 0.0689 0.0500 0.0497
disengaged 0.0902 0.0856 0.0551 0.0842 0.0957 0.0939 0.0956 0.0937 0.0853 0.1042 0.0839 0.0903 0.1203
sub-optimal 0.0839 0.0973 0.0460 0.0500 0.0953 0.1131 0.0965 0.0647 0.0637 0.0844 0.0721 0.0532 0.0684
face_std_10 mean over mice: {'optimal': 0.0661, 'disengaged': 0.0906, 'sub-optimal': 0.076}
optimal vs disengaged: difference -0.0245, optimal lower in 11 of 13, p = 0.0024
optimal vs sub-optimal: difference -0.0099, optimal lower in 11 of 13, p = 0.017
running_std10 per mouse:
BW041 BW045 BW046 BW048 BW051 BW053 BW058 BW059 BW061 BW062 BW064 BW067 BW071
optimal 0.0220 0.0363 0.0572 0.0519 0.0358 0.0449 0.0555 0.0680 0.0749 0.0516 0.0675 0.0805 0.0419
disengaged 0.0168 0.0684 0.0844 0.0862 0.0548 0.0290 0.0868 0.1351 0.1139 0.0882 0.0944 0.1115 0.0590
sub-optimal 0.0207 0.0448 0.0704 0.0814 0.0588 0.0301 0.0767 0.0835 0.0657 0.0579 0.0808 0.0787 0.0470
running_std10 mean over mice: {'optimal': 0.0529, 'disengaged': 0.0791, 'sub-optimal': 0.0613}
optimal vs disengaged: difference -0.0262, optimal lower in 11 of 13, p = 0.0012
optimal vs sub-optimal: difference -0.0083, optimal lower in 9 of 13, p = 0.057
tests with the optimal state lower at p < 0.05: 5 of 6
Result: Partially reproduced. All six differences have the stated sign, and five of the six tests give p < 0.05; locomotion speed variability in the optimal state against the sub-optimal states does not (p = 0.057).
Variability of the measures decreases on entering the optimal state and increases on leaving it.
As Figure4.py: optimal epochs from flex_hmm.get_dwell_times on all included sessions with a code dwell above 9; entries from a discrete non-optimal epoch and exits to one; the 10-trial standard deviation over the 10 trials before and the 10 after the boundary, averaged per mouse; the test compares position 10 (the last trial before the entry, or the last optimal trial) with position 20.
enter_means, exit_means = {m: [] for m in STD3}, {m: [] for m in STD3}
n_trans = {}
for subj in subjects:
g = trials[trials.subject == subj].sort_values(["path", "trial_in_session"])
E = {m: [] for m in STD3}
X_ = {m: [] for m in STD3}
ne = nx = 0
for path, sg in g.groupby("path"):
dt = hulsey.dwell_table(sg.hmm_state.values)
for SI in np.where(dt.state.values == 0)[0]:
if dt.dwell.iloc[SI] <= 9:
continue
start, end = dt.first_trial.iloc[SI], dt.last_trial.iloc[SI]
prior = dt.state.iloc[SI - 1] if SI > 0 else np.nan
nxt = dt.state.iloc[SI + 1] if SI < len(dt) - 1 else np.nan
for m in STD3:
v = sg[m].values
a, b = v[start - 10:start + 10], v[end - 10:end + 10]
if prior > 0 and len(a) == 20:
E[m].append(a)
if nxt > 0 and len(b) == 20:
X_[m].append(b)
ne += prior > 0
nx += nxt > 0
for m in STD3:
enter_means[m].append(np.nanmean(E[m], axis=0))
exit_means[m].append(np.nanmean(X_[m], axis=0))
n_trans[subj] = (ne, nx)
print("transitions per mouse (entries, exits):", n_trans)
fig, ax = plt.subplots(1, 3, figsize=(13, 3))
for i, m in enumerate(STD3):
en, ex = np.array(enter_means[m]), np.array(exit_means[m])
we, wx = stats.wilcoxon(en[:, 9], en[:, 19]), stats.wilcoxon(ex[:, 9], ex[:, 19])
print(f"{m}: entry {np.mean(en[:, 9]):.4f} -> {np.mean(en[:, 19]):.4f} (lower after in {int((en[:, 19] < en[:, 9]).sum())} of 13, p = {we.pvalue:.2g}); "
f"exit {np.mean(ex[:, 9]):.4f} -> {np.mean(ex[:, 19]):.4f} (higher after in {int((ex[:, 19] > ex[:, 9]).sum())} of 13, p = {wx.pvalue:.2g})")
ax[i].plot(np.arange(-10, 10), (en - en.mean(1, keepdims=True)).mean(0), label="entry")
ax[i].plot(np.arange(-10, 10), (ex - ex.mean(1, keepdims=True)).mean(0), label="exit")
ax[i].axvline(0, color="k", lw=.5)
ax[i].set_title(m, fontsize=9)
ax[0].legend()
ax[0].set_ylabel("mean-subtracted 10-trial SD")
ax[1].set_xlabel("trials from the boundary")
plt.tight_layout()
plt.show()
transitions per mouse (entries, exits): {'BW041': (np.int64(52), np.int64(50)), 'BW045': (np.int64(85), np.int64(125)), 'BW046': (np.int64(55), np.int64(69)), 'BW048': (np.int64(73), np.int64(79)), 'BW051': (np.int64(63), np.int64(67)), 'BW053': (np.int64(16), np.int64(32)), 'BW058': (np.int64(9), np.int64(19)), 'BW059': (np.int64(74), np.int64(79)), 'BW061': (np.int64(10), np.int64(18)), 'BW062': (np.int64(65), np.int64(67)), 'BW064': (np.int64(39), np.int64(43)), 'BW067': (np.int64(29), np.int64(46)), 'BW071': (np.int64(15), np.int64(31))}
post_hoc_pupil_std10: entry 0.0319 -> 0.0283 (lower after in 7 of 13, p = 0.45); exit 0.0301 -> 0.0346 (higher after in 11 of 13, p = 0.0012)
face_std_10: entry 0.0750 -> 0.0680 (lower after in 9 of 13, p = 0.11); exit 0.0713 -> 0.0780 (higher after in 8 of 13, p = 0.13)
running_std10: entry 0.0620 -> 0.0545 (lower after in 9 of 13, p = 0.04); exit 0.0553 -> 0.0640 (higher after in 11 of 13, p = 0.0012)
Result: Partially reproduced. All three measures' variability falls on entering the optimal state and rises on leaving it on average. Of the four comparisons the figure marks, pupil exit (p = 0.0012), locomotion entry (p = 0.04) and locomotion exit (p = 0.0012) give p < 0.05, and face entry does not (p = 0.11).
Coupling of pupil and movement (claims fig4d-pupil-movement-decoupling, figs6-individual-measure-decoupling)¶
The correlation between pupil diameter and the movement index is lower during the optimal state than during the disengaged and sub-optimal states.
As Figure4.py: the trial set of fig4a, per mouse and state the Pearson r between movement index and pupil diameter.
rr = pd.DataFrame({subj: {"optimal": np.corrcoef(g.loc[g.hmm_state == 0, "movement"], g.loc[g.hmm_state == 0, "post_hoc_pupil_diameter"])[0, 1],
"disengaged": np.corrcoef(g.loc[g.hmm_state == 1, "movement"], g.loc[g.hmm_state == 1, "post_hoc_pupil_diameter"])[0, 1],
"sub-optimal": np.corrcoef(g.loc[g.hmm_state > 1, "movement"], g.loc[g.hmm_state > 1, "post_hoc_pupil_diameter"])[0, 1]}
for subj, g in fig4.items()}).T
print(rr.round(3).to_string())
print("mean:", rr.mean().round(3).to_dict())
for a_, b_ in [("optimal", "disengaged"), ("optimal", "sub-optimal"), ("disengaged", "sub-optimal")]:
print(f"{a_} vs {b_}: difference {np.mean(rr[a_] - rr[b_]):+.3f}, first lower in {int((rr[a_] < rr[b_]).sum())} of 13, "
f"p = {stats.wilcoxon(rr[a_], rr[b_]).pvalue:.2g}")
optimal disengaged sub-optimal
BW041 0.207 0.532 0.542
BW045 0.516 0.604 0.610
BW046 -0.195 0.029 0.034
BW048 0.057 0.550 0.082
BW051 -0.134 0.537 0.168
BW053 0.103 0.538 0.300
BW058 0.059 0.658 0.484
BW059 0.275 0.506 0.306
BW061 -0.114 0.243 0.229
BW062 0.265 0.639 0.548
BW064 0.011 0.116 0.112
BW067 -0.173 0.255 -0.051
BW071 0.121 0.498 0.105
mean: {'optimal': 0.077, 'disengaged': 0.439, 'sub-optimal': 0.267}
optimal vs disengaged: difference -0.362, first lower in 13 of 13, p = 0.00024
optimal vs sub-optimal: difference -0.190, first lower in 12 of 13, p = 0.00049
disengaged vs sub-optimal: difference +0.172, first lower in 3 of 13, p = 0.0081
Result: Reproduced. The pupil-movement correlation is 0.077 in the optimal state against 0.439 in the disengaged and 0.267 in the sub-optimal states, lower in 13 and 12 of 13 mice (p = 0.00024 and 0.00049).
The correlations of pupil diameter with face motion energy and with locomotion speed are also reduced during the optimal state.
As FigureS6.py: sessions of M5 with the six measures present; per mouse and state the Pearson r.
for a_m, b_m in [("face_energy", "post_hoc_pupil_diameter"), ("running_speed", "post_hoc_pupil_diameter"), ("running_speed", "face_energy")]:
tab = {}
for subj in subjects:
g = fig3_set(subj)
tab[subj] = {nm: np.corrcoef(g.loc[sel, a_m], g.loc[sel, b_m])[0, 1]
for nm, sel in [("optimal", g.hmm_state == 0), ("disengaged", g.hmm_state == 1), ("sub-optimal", g.hmm_state > 1)]}
tab = pd.DataFrame(tab).T
print(f"r({a_m}, {b_m}) per mouse:")
print(tab.round(3).T.to_string())
print(f"r({a_m}, {b_m}) mean over mice:", tab.mean().round(3).to_dict())
for other in ["disengaged", "sub-optimal"]:
print(f" optimal vs {other}: difference {np.mean(tab.optimal - tab[other]):+.3f}, optimal lower in "
f"{int((tab.optimal < tab[other]).sum())} of 13, p = {stats.wilcoxon(tab.optimal, tab[other]).pvalue:.2g}")
r(face_energy, post_hoc_pupil_diameter) per mouse:
BW041 BW045 BW046 BW048 BW051 BW053 BW058 BW059 BW061 BW062 BW064 BW067 BW071
optimal 0.271 0.200 -0.077 0.063 -0.231 0.119 0.068 0.300 -0.167 0.001 0.299 -0.265 0.450
disengaged 0.339 0.479 0.018 0.396 0.556 0.408 0.603 0.381 0.006 0.620 0.211 0.391 0.441
sub-optimal 0.539 0.332 0.068 0.143 0.035 0.260 0.377 0.550 0.272 0.401 0.262 -0.108 0.523
r(face_energy, post_hoc_pupil_diameter) mean over mice: {'optimal': 0.079, 'disengaged': 0.373, 'sub-optimal': 0.281}
optimal vs disengaged: difference -0.294, optimal lower in 11 of 13, p = 0.0024
optimal vs sub-optimal: difference -0.202, optimal lower in 12 of 13, p = 0.00049
r(running_speed, post_hoc_pupil_diameter) per mouse:
BW041 BW045 BW046 BW048 BW051 BW053 BW058 BW059 BW061 BW062 BW064 BW067 BW071
optimal 0.047 0.674 -0.212 0.017 0.050 0.049 0.018 0.072 0.019 0.436 -0.333 0.005 -0.388
disengaged 0.519 0.506 0.024 0.448 0.457 0.473 0.590 0.507 0.071 0.551 -0.013 0.082 0.317
sub-optimal 0.271 0.651 -0.020 -0.012 0.204 0.220 0.502 -0.094 -0.033 0.433 -0.072 0.048 -0.324
r(running_speed, post_hoc_pupil_diameter) mean over mice: {'optimal': 0.035, 'disengaged': 0.349, 'sub-optimal': 0.136}
optimal vs disengaged: difference -0.314, optimal lower in 12 of 13, p = 0.0017
optimal vs sub-optimal: difference -0.101, optimal lower in 8 of 13, p = 0.068
r(running_speed, face_energy) per mouse:
BW041 BW045 BW046 BW048 BW051 BW053 BW058 BW059 BW061 BW062 BW064 BW067 BW071
optimal 0.018 0.223 0.087 0.042 0.246 0.239 0.120 -0.019 0.181 0.162 -0.119 0.037 -0.439
disengaged 0.321 0.318 0.066 0.171 0.745 0.306 0.638 0.572 0.426 0.675 0.312 0.361 0.451
sub-optimal 0.186 0.300 0.054 0.149 0.442 0.305 0.675 0.092 0.067 0.162 0.025 0.132 -0.201
r(running_speed, face_energy) mean over mice: {'optimal': 0.06, 'disengaged': 0.412, 'sub-optimal': 0.184}
optimal vs disengaged: difference -0.352, optimal lower in 12 of 13, p = 0.00049
optimal vs sub-optimal: difference -0.124, optimal lower in 11 of 13, p = 0.01
Result: Reproduced. The correlation of pupil with face motion energy is lower in the optimal state than in the disengaged (p = 0.0024) and sub-optimal states (p = 0.00049); with locomotion speed it is lower than in the disengaged state (p = 0.0017) and not significantly lower than in the sub-optimal states (p = 0.068).
Trial history and evoked pupil (claim figs1-iti-independence)¶
The inter-trial interval is long enough that trial-evoked responses do not carry over; correlations of the measures at trial start with the prior ITI are weak and centered on zero, and evoked pupil responses do not differ across performance states.
As FigureS1.py on all included sessions: the ITI is the difference between consecutive trial start times within a session; per mouse, Pearson r of each measure at stimulus onset with the preceding ITI (the pupil value at onset is used for the code's trace sample 1 ms later); the evoked response is the signed larger of the largest rise and fall of pupil within 2.5 s of onset; per mouse the median per group, and a Kruskal-Wallis test over mice.
itr, ev = [], []
for subj in subjects:
g = trials[trials.subject == subj].sort_values(["path", "trial_in_session"]).copy()
g["iti"] = g.groupby("path").start_time.diff()
g = g.dropna(subset=["evoked_pupil", "post_hoc_pupil_diameter", "iti", "running_speed", "face_energy"])
itr.append({"subject": subj, **{m: stats.pearsonr(g.iti, g[m])[0] for m in ["post_hoc_pupil_diameter", "running_speed", "face_energy"]}})
ev.append({"subject": subj, "optimal": g.loc[g.hmm_state == 0, "evoked_pupil"].median() * 100,
"disengaged": g.loc[g.hmm_state == 1, "evoked_pupil"].median() * 100,
"sub-optimal": g.loc[g.hmm_state > 1, "evoked_pupil"].median() * 100,
"indeterminate": g.loc[g.hmm_state.isna(), "evoked_pupil"].median() * 100})
itr = pd.DataFrame(itr).set_index("subject")
ev = pd.DataFrame(ev).set_index("subject")
print("r with the preceding ITI:")
print(itr.round(3).to_string())
for m in itr.columns:
print(f"{m}: median r {itr[m].median():+.3f}, Wilcoxon against 0 p = {stats.wilcoxon(itr[m]).pvalue:.2g}; "
f"(a) |median| < 0.1: {abs(itr[m].median()) < .1}, (b) < 0.2: {abs(itr[m].median()) < .2}")
print("\nmedian evoked pupil response per mouse (% of max):")
print(ev.round(2).to_string())
kw4 = stats.kruskal(*[ev[c] for c in ev.columns])
kw3 = stats.kruskal(*[ev[c] for c in ["optimal", "disengaged", "sub-optimal"]])
print(f"Kruskal-Wallis, four groups: H = {kw4.statistic:.2f}, p = {kw4.pvalue:.2g}; three discrete states (variant): p = {kw3.pvalue:.2g}")
print("group means:", ev.mean().round(2).to_dict())
r with the preceding ITI:
post_hoc_pupil_diameter running_speed face_energy
subject
BW041 -0.150 0.058 -0.070
BW045 -0.031 -0.011 -0.123
BW046 0.079 0.058 -0.081
BW048 0.081 0.134 -0.026
BW051 -0.044 -0.007 0.138
BW053 -0.276 -0.085 -0.070
BW058 -0.063 -0.056 0.123
BW059 0.036 -0.009 0.020
BW061 0.012 0.006 0.006
BW062 -0.002 -0.036 -0.065
BW064 -0.229 0.142 -0.119
BW067 -0.089 -0.153 -0.003
BW071 -0.058 0.074 -0.044
post_hoc_pupil_diameter: median r -0.044, Wilcoxon against 0 p = 0.13; (a) |median| < 0.1: True, (b) < 0.2: True
running_speed: median r -0.007, Wilcoxon against 0 p = 0.89; (a) |median| < 0.1: True, (b) < 0.2: True
face_energy: median r -0.044, Wilcoxon against 0 p = 0.27; (a) |median| < 0.1: True, (b) < 0.2: True
median evoked pupil response per mouse (% of max):
optimal disengaged sub-optimal indeterminate
subject
BW041 1.33 1.26 1.52 -0.76
BW045 3.88 2.43 2.81 3.43
BW046 -9.39 -2.83 -6.40 -7.49
BW048 -5.47 2.33 -3.89 -4.42
BW051 -6.27 1.47 -3.67 -4.40
BW053 6.84 5.25 8.56 7.98
BW058 -4.47 -1.03 2.72 -3.05
BW059 -1.83 2.93 1.67 1.03
BW061 -3.25 1.30 -2.27 -2.52
BW062 2.82 1.00 1.16 1.89
BW064 0.85 1.41 1.32 1.49
BW067 -1.76 1.54 -1.03 -1.34
BW071 1.23 3.35 0.04 1.53
Kruskal-Wallis, four groups: H = 4.38, p = 0.22; three discrete states (variant): p = 0.15
group means: {'optimal': -1.19, 'disengaged': 1.57, 'sub-optimal': 0.2, 'indeterminate': -0.51}
Result: Reproduced. The median correlations of pupil, locomotion and face motion energy with the preceding ITI are -0.044, -0.007 and -0.044, none different from zero (p = 0.13, 0.89, 0.27), within reading (a). The median evoked pupil responses do not differ across the four groups (Kruskal-Wallis p = 0.22; three discrete states p = 0.15).
Regression of state probabilities on the measures (claims figs8a-regression-weights, figs8b-pupil-movement-interaction)¶
In linear regressions of state probability on the measures, the quadratic weights of pupil diameter and of locomotion speed are negative for the optimal state and positive for the disengaged state, as are the weights of their 10-trial SDs, while face motion energy weights center on zero.
As FigureS8.py: sessions of M5 with the six measures; per mouse, OLS of each state's posterior on the min-max scaled design (constant, pupil P, face F, locomotion L, their squares and products, their 10-trial SDs, and PFL; the design scaled again); second model on pupil, movement index M, their squares and product, the pupil SD and the movement SD index.
import statsmodels.api as sm
from sklearn.preprocessing import MinMaxScaler, PolynomialFeatures
names1 = ["1", "P", "F", "L", "P*P", "P*F", "P*L", "F*F", "F*L", "L*L", "Pstd", "Fstd", "Lstd", "P*F*L"]
names2 = ["1", "P", "M", "P*P", "P*M", "M*M", "Pstd", "Mstd"]
W1, W2 = [], []
for subj in subjects:
g = fig3_set(subj)
mm = MinMaxScaler()
X = mm.fit_transform(g[["post_hoc_pupil_diameter", "face_energy", "running_speed"]].values)
S = mm.fit_transform(g[["post_hoc_pupil_std10", "face_std_10", "running_std10"]].values)
X2 = mm.fit_transform(np.c_[g.post_hoc_pupil_diameter.values, g.movement.values])
pf = PolynomialFeatures(degree=2)
A1 = mm.fit_transform(np.c_[pf.fit_transform(X), S, X[:, 0] * X[:, 1] * X[:, 2]])
A1[:, 0] = 1
A2 = mm.fit_transform(np.c_[pf.fit_transform(X2), g.post_hoc_pupil_std10.values, g.movement_std.values])
A2[:, 0] = 1
for k in range(n_states[subj]):
yk = g[f"state_{k + 1}_prob"].values
grp = "optimal" if k == 0 else "disengaged" if k == 1 else "sub-optimal"
W1.append({"subject": subj, "state": grp, **dict(zip(names1, sm.OLS(yk, A1).fit().params))})
W2.append({"subject": subj, "state": grp, **dict(zip(names2, sm.OLS(yk, A2).fit().params))})
W1, W2 = pd.DataFrame(W1), pd.DataFrame(W2)
rows = []
for grp in ["optimal", "disengaged"]:
w = W1[W1.state == grp]
for term in ["P*P", "L*L", "F*F", "Pstd", "Lstd", "Fstd"]:
rows.append({"state": grp, "term": term, "median": w[term].median(), "mean": w[term].mean(),
"positive": int((w[term] > 0).sum()), "p (vs 0)": stats.wilcoxon(w[term]).pvalue})
s8 = pd.DataFrame(rows)
print(s8.round(3).to_string(index=False))
print("p (vs 0) with two significant digits:", [f"{r.state} {r.term} {r['p (vs 0)']:.2g}" for _, r in s8.iterrows()])
state term median mean positive p (vs 0) optimal P*P -1.114 -1.248 0 0.000 optimal L*L -1.098 -1.135 1 0.001 optimal F*F -0.023 -0.117 6 0.839 optimal Pstd -0.260 -0.404 2 0.008 optimal Lstd -0.332 -0.273 4 0.027 optimal Fstd -0.035 -0.073 6 0.455 disengaged P*P 0.560 0.895 12 0.000 disengaged L*L 0.644 0.772 13 0.000 disengaged F*F 0.219 0.160 9 0.455 disengaged Pstd 0.193 0.297 9 0.033 disengaged Lstd 0.197 0.211 11 0.008 disengaged Fstd 0.129 0.110 9 0.191 p (vs 0) with two significant digits: ['optimal P*P 0.00024', 'optimal L*L 0.0012', 'optimal F*F 0.84', 'optimal Pstd 0.0081', 'optimal Lstd 0.027', 'optimal Fstd 0.45', 'disengaged P*P 0.00049', 'disengaged L*L 0.00024', 'disengaged F*F 0.45', 'disengaged Pstd 0.033', 'disengaged Lstd 0.0081', 'disengaged Fstd 0.19']
Result: Reproduced. For the optimal state, the quadratic and SD weights of pupil and locomotion are negative (Wilcoxon against zero p at most 0.027) and positive for the disengaged state (p at most 0.033), while the face weights do not differ from zero (p from 0.191 to 0.839), under reading (a) and so under (b).
With the movement index, the quadratic movement weight is negative for the optimal state, and the pupil-by-movement interaction is the only consistent weight for sub-optimal states (0.42 +/- 0.11), differing from its weight for the optimal state (-0.95 +/- 0.21).
for grp in ["optimal", "disengaged", "sub-optimal"]:
w = W2[W2.state == grp]
print(grp, f"({len(w)} fits): mean +/- SEM", {t: f"{w[t].mean():.2f} +/- {w[t].sem():.2f}" for t in names2[1:]})
print(" Wilcoxon vs 0 p:", {t: round(stats.wilcoxon(w[t]).pvalue, 3) for t in names2[1:]})
pm_s, pm_o = W2[W2.state == "sub-optimal"]["P*M"], W2[W2.state == "optimal"]["P*M"]
print(f"P*M sub-optimal {pm_s.mean():.2f} +/- {pm_s.sem():.2f} vs optimal {pm_o.mean():.2f} +/- {pm_o.sem():.2f}; "
f"rank-sum p = {stats.ranksums(pm_s, pm_o).pvalue:.2g}")
print(f"M*M optimal mean {W2[W2.state == 'optimal']['M*M'].mean():.2f}")
optimal (13 fits): mean +/- SEM {'P': '1.36 +/- 0.41', 'M': '1.99 +/- 0.38', 'P*P': '-1.26 +/- 0.35', 'P*M': '-0.95 +/- 0.21', 'M*M': '-1.30 +/- 0.31', 'Pstd': '-0.38 +/- 0.15', 'Mstd': '-0.32 +/- 0.13'}
Wilcoxon vs 0 p: {'P': np.float64(0.002), 'M': np.float64(0.0), 'P*P': np.float64(0.0), 'P*M': np.float64(0.0), 'M*M': np.float64(0.0), 'Pstd': np.float64(0.021), 'Mstd': np.float64(0.033)}
disengaged (13 fits): mean +/- SEM {'P': '-0.81 +/- 0.25', 'M': '-1.46 +/- 0.25', 'P*P': '0.92 +/- 0.20', 'P*M': '0.29 +/- 0.21', 'M*M': '1.29 +/- 0.23', 'Pstd': '0.29 +/- 0.15', 'Mstd': '0.28 +/- 0.10'}
Wilcoxon vs 0 p: {'P': np.float64(0.001), 'M': np.float64(0.0), 'P*P': np.float64(0.0), 'P*M': np.float64(0.216), 'M*M': np.float64(0.0), 'Pstd': np.float64(0.057), 'Mstd': np.float64(0.013)}
sub-optimal (21 fits): mean +/- SEM {'P': '-0.34 +/- 0.29', 'M': '-0.33 +/- 0.22', 'P*P': '0.21 +/- 0.24', 'P*M': '0.41 +/- 0.10', 'M*M': '0.00 +/- 0.22', 'Pstd': '0.06 +/- 0.05', 'Mstd': '0.02 +/- 0.07'}
Wilcoxon vs 0 p: {'P': np.float64(0.179), 'M': np.float64(0.229), 'P*P': np.float64(0.355), 'P*M': np.float64(0.0), 'M*M': np.float64(0.892), 'Pstd': np.float64(0.562), 'Mstd': np.float64(0.585)}
P*M sub-optimal 0.41 +/- 0.10 vs optimal -0.95 +/- 0.21; rank-sum p = 3.8e-06
M*M optimal mean -1.30
Result: Reproduced. The pupil-by-movement weight is 0.41 +/- 0.10 over the 21 sub-optimal fits (expected 0.42 +/- 0.11) and -0.95 +/- 0.21 over the optimal fits (expected -0.95 +/- 0.21), different by a rank-sum test (p = 3.8e-06), and the mean quadratic movement weight of the optimal fits is negative (-1.30). It is the only sub-optimal weight that differs from zero (p < 0.001).
Decoding the state from the measures (claims fig5c-state-classification, figs7-pupil-best-raw-measure, fig5c-movement-index-best)¶
For every mouse, support vector machines on the value and 10-trial variability of the measures predict the performance state above chance, with 88.2% +/- 1.8% accuracy for optimal against disengaged and 77.7% +/- 2% for optimal against sub-optimal states.
As 5_decoding.py: sessions of M5; discrete-state trials with the six measures; the movement index and its SD z-scored over these trials; for each mouse, optimal against each other state when the two classes give more than 30 trials; at most 1500 trials per class, balanced, sampled at random; features min-max scaled; RBF SVC; 5-fold shuffled KFold; accuracy is the mean over folds and runs. In a run, the condition's columns are permuted across the training trials. The paper runs 1000 runs per condition; here the unshuffled condition needs one run (it does not vary), the shuffled conditions 3 runs, and the null (all columns shuffled) 20 runs, to fit the time limit (see the deviations). z = (accuracy - null mean) / null SD.
N_RUN_COND, N_RUN_NULL = 3, 20
rng5 = np.random.default_rng(5)
t0 = time.time()
clf = []
tasks, tmeta = [], []
for subj in subjects:
g = trials[(trials.subject == subj) & trials.sess_opt10].dropna(subset=["hmm_state"] + MEAS6).copy()
g["movement"] = stats.zscore(g.face_energy) + stats.zscore(g.running_speed)
g["movement_std"] = stats.zscore(g.face_std_10) + stats.zscore(g.running_std10)
for st in sorted(g.hmm_state.unique())[1:]:
o = g[g.hmm_state == 0].sample(frac=1, random_state=int(rng5.integers(2**31)))
t_ = g[g.hmm_state == st].sample(frac=1, random_state=int(rng5.integers(2**31)))
n = min(len(o), len(t_), 1500)
o, t_ = o.iloc[:n], t_.iloc[:n]
y = np.r_[np.zeros(n), np.ones(n)]
if len(y) <= 30:
continue
X = MinMaxScaler().fit_transform(np.r_[o[hulsey.DECODE_MEASURES].values, t_[hulsey.DECODE_MEASURES].values])
cid = len(clf)
clf.append({"subject": subj, "state": int(st), "type": "disengaged" if st == 1 else "sub-optimal",
"template": hulsey.STATE_NAMES[types[subj][int(st)]], "trials per class": n})
from sklearn.model_selection import KFold
for f, (a, b) in enumerate(KFold(5, shuffle=True, random_state=int(rng5.integers(2**31))).split(y)):
for cond, cols in hulsey.SHUFFLE.items():
nr = 1 if cond == 13 else N_RUN_NULL if cond == 14 else N_RUN_COND
for r in range(nr):
tasks.append((X[a], y[a], X[b], y[b], cols, int(rng5.integers(2**31))))
tmeta.append((cid, cond, r, f))
res = hulsey.run_parallel(hulsey.svm_task, tasks)
print(f"{len(clf)} classifications, {len(tasks)} SVM fits in {(time.time() - t0) / 60:.1f} min")
acc5 = pd.DataFrame(tmeta, columns=["cid", "cond", "run", "fold"]).assign(acc=res)
per_run = acc5.groupby(["cid", "cond", "run"]).acc.mean() # mean over folds
null = per_run.xs(14, level="cond").groupby("cid").agg(["mean", "std", lambda x: np.percentile(x, 95)])
null.columns = ["null mean", "null sd", "null p95"]
cond_acc = per_run.groupby(["cid", "cond"]).mean().unstack()
clf = pd.DataFrame(clf).join(null)
for cond in cond_acc.columns:
clf[f"acc {cond}"] = cond_acc[cond] * 100
clf[f"z {cond}"] = (cond_acc[cond] - clf["null mean"]) / clf["null sd"]
show = clf[["subject", "state", "type", "template", "trials per class", "acc 13", "acc 0", "null mean", "null sd", "z 13", "z 0"]].copy()
show["null mean"] *= 100
show["null sd"] *= 100
print(show.round(2).to_string(index=False))
for cond, label in [(13, "eight features unshuffled (condition 13)"), (0, "'All measures' of Figure5.py (condition 0)")]:
for tp in ["disengaged", "sub-optimal"]:
c = clf[clf.type == tp][f"acc {cond}"]
print(f"{label}, optimal vs {tp}: {c.mean():.1f} +/- {c.sem():.1f}% over {len(c)} classifications")
z = clf[f"z {cond}"]
above95 = (clf[f"acc {cond}"] / 100 > clf["null p95"])
print(f" significant: (a) z > 1.98 {int((z > 1.98).sum())} of {len(z)}, (b) z > 1.645 {int((z > 1.645).sum())}, "
f"(c) above the null's 95th percentile {int(above95.sum())}")
33 classifications, 8415 SVM fits in 4.1 min subject state type template trials per class acc 13 acc 0 null mean null sd z 13 z 0 BW041 1 disengaged disengaged 1357 88.17 87.60 50.60 4.09 9.18 9.04 BW041 2 sub-optimal avoid left 1500 74.90 74.74 49.79 3.28 7.66 7.61 BW045 1 disengaged disengaged 1500 82.73 82.14 50.26 3.83 8.49 8.34 BW045 2 sub-optimal bias right 1500 59.53 58.39 49.57 1.41 7.07 6.26 BW045 3 sub-optimal avoid right 1296 86.65 86.25 49.39 7.02 5.31 5.25 BW046 1 disengaged disengaged 625 77.12 76.32 51.06 2.77 9.40 9.11 BW046 2 sub-optimal bias right 899 78.92 78.57 50.06 2.22 13.00 12.84 BW046 3 sub-optimal bias left 1336 70.13 68.82 49.76 1.76 11.59 10.84 BW048 1 disengaged disengaged 585 88.21 87.66 50.38 4.18 9.05 8.92 BW048 2 sub-optimal avoid left 1134 81.39 81.01 49.71 4.12 7.69 7.60 BW048 3 sub-optimal bias left 1500 78.07 77.49 49.01 2.71 10.71 10.49 BW051 1 disengaged disengaged 396 94.06 93.90 49.97 9.04 4.88 4.86 BW051 2 sub-optimal avoid right 1500 73.70 73.00 49.19 4.61 5.32 5.17 BW053 1 disengaged disengaged 1500 88.90 88.58 51.05 3.13 12.10 12.00 BW053 2 sub-optimal bias left 925 83.62 82.13 50.14 4.08 8.20 7.83 BW053 3 sub-optimal bias right 972 71.55 70.47 50.38 2.37 8.92 8.46 BW058 1 disengaged disengaged 456 91.45 91.05 50.38 5.75 7.14 7.07 BW058 2 sub-optimal avoid right 835 89.70 89.28 50.92 5.21 7.44 7.36 BW059 1 disengaged disengaged 642 95.56 95.48 51.98 6.96 6.26 6.25 BW059 3 sub-optimal avoid left 1500 69.47 68.90 49.34 2.44 8.25 8.02 BW061 1 disengaged disengaged 1500 92.63 92.48 49.42 5.74 7.52 7.50 BW061 2 sub-optimal bias left 173 92.19 91.52 49.02 6.41 6.74 6.64 BW061 3 sub-optimal bias right 1500 75.57 74.83 49.41 3.09 8.46 8.22 BW062 1 disengaged disengaged 777 93.69 93.37 50.49 5.54 7.80 7.74 BW062 2 sub-optimal avoid left 1500 77.20 76.97 47.51 3.48 8.54 8.48 BW064 1 disengaged disengaged 1500 77.33 76.16 49.22 2.79 10.08 9.66 BW064 2 sub-optimal avoid right 1164 73.07 71.89 49.55 2.42 9.72 9.24 BW067 1 disengaged disengaged 1500 93.07 92.71 50.14 4.60 9.33 9.25 BW067 2 sub-optimal bias left 1500 77.80 77.57 51.24 3.71 7.16 7.09 BW067 3 sub-optimal bias right 376 91.63 91.31 46.86 6.36 7.04 6.99 BW071 1 disengaged disengaged 1365 93.19 92.92 52.10 4.37 9.41 9.35 BW071 2 sub-optimal bias left 764 78.60 77.07 49.31 2.96 9.90 9.38 BW071 3 sub-optimal bias right 1033 84.08 83.85 49.52 3.90 8.87 8.81 eight features unshuffled (condition 13), optimal vs disengaged: 88.9 +/- 1.7% over 13 classifications eight features unshuffled (condition 13), optimal vs sub-optimal: 78.4 +/- 1.8% over 20 classifications significant: (a) z > 1.98 33 of 33, (b) z > 1.645 33, (c) above the null's 95th percentile 33 'All measures' of Figure5.py (condition 0), optimal vs disengaged: 88.5 +/- 1.8% over 13 classifications 'All measures' of Figure5.py (condition 0), optimal vs sub-optimal: 77.7 +/- 1.9% over 20 classifications significant: (a) z > 1.98 33 of 33, (b) z > 1.645 33, (c) above the null's 95th percentile 33
Result: Reproduced. With the eight features unshuffled, accuracy is 88.9% +/- 1.7% for optimal against disengaged (13 classifications) and 78.4% +/- 1.8% against sub-optimal states (20), within tolerance; the 'All measures' condition gives 88.5% and 77.7%. All 33 classifications are significant under every reading. Here 20 sub-optimal classifications are possible against 19 in the paper.
Pupil diameter (value and variability) is the best single raw measure for state classification, better than face motion energy and only slightly better than locomotion speed.
Readings (a), Figure S7A: keep one measure's value and SD, shuffle all others (conditions 2 pupil, 3 face, 4 locomotion); (b), Figure S7B: shuffle one measure's value and SD with the movement index (conditions 6, 7, 8), comparing the drop from condition 13. Quantity (i) accuracy, (ii) z.
for reading, (cp, cf, cl), drop in [("(a) keep one", (2, 3, 4), False), ("(b) shuffle one", (6, 7, 8), True)]:
for q in ["acc", "z"]:
P, F, L = (clf[f"{q} {c}"] for c in (cp, cf, cl))
if drop:
base = clf[f"{q} 13"]
P, F, L = base - P, base - F, base - L # a larger drop means more information
wf, wl = stats.wilcoxon(P, F), stats.wilcoxon(P, L)
print(f"{reading}, {q}{' drop' if drop else ''}: pupil {P.mean():.2f}, face {F.mean():.2f}, locomotion {L.mean():.2f}; "
f"pupil > face in {int((P > F).sum())} of {len(P)}, p = {wf.pvalue:.2g}; pupil > locomotion in {int((P > L).sum())}, p = {wl.pvalue:.2g}")
(a) keep one, acc: pupil 71.51, face 66.09, locomotion 68.10; pupil > face in 24 of 33, p = 0.0011; pupil > locomotion in 23, p = 0.053 (a) keep one, z: pupil 5.21, face 4.02, locomotion 4.37; pupil > face in 24 of 33, p = 0.0013; pupil > locomotion in 23, p = 0.039 (b) shuffle one, acc drop: pupil 7.71, face 4.45, locomotion 5.59; pupil > face in 25 of 33, p = 0.00029; pupil > locomotion in 23, p = 0.029 (b) shuffle one, z drop: pupil 2.15, face 1.40, locomotion 1.63; pupil > face in 25 of 33, p = 0.00056; pupil > locomotion in 23, p = 0.022
Result: Reproduced under reading (a). Keeping one measure, pupil gives 71.5% against 66.1% for face motion energy (pupil higher in 24 of 33, p = 0.0011) and 68.1% for locomotion speed (higher on average, p = 0.053; 0.051 in a run on another machine); reading (b) agrees (face p = 0.00029, locomotion p = 0.029).
The movement index is the best predictor of performance state, encoding more information about states than pupil diameter.
Main reading (Figure 5C legend): keep the movement index value and SD, shuffle all others ("M only"), against keeping only the pupil value and SD (condition 2). Variants, as the criterion asks: condition 5 of 5_decoding.py (movement index kept with the pupil value) and condition 1 (face and locomotion shuffled), the label "Movement index" of the final Figure5.py. Both keep pupil information in the movement-index classifier, so they compare pupil and movement together against pupil alone and cannot show that the movement index carries more information than pupil.
for label, cm in [("main: M only", "M only"), ("variant: condition 5", 5), ("variant: condition 1", 1)]:
for q in ["acc", "z"]:
M, P = clf[f"{q} {cm}"], clf[f"{q} 2"]
w = stats.wilcoxon(M, P)
print(f"{label}, {q}: movement {M.mean():.2f} vs pupil {P.mean():.2f}; movement higher in {int((M > P).sum())} of {len(M)}, p = {w.pvalue:.2g}")
print(f"\nnotebook run time from the first cell: {(time.time() - started) / 60:.1f} min")
main: M only, acc: movement 68.63 vs pupil 71.51; movement higher in 14 of 33, p = 0.086 main: M only, z: movement 4.51 vs pupil 5.21; movement higher in 14 of 33, p = 0.068 variant: condition 5, acc: movement 75.67 vs pupil 71.51; movement higher in 29 of 33, p = 1.3e-05 variant: condition 5, z: movement 6.33 vs pupil 5.21; movement higher in 29 of 33, p = 3.1e-05 variant: condition 1, acc: movement 77.86 vs pupil 71.51; movement higher in 33 of 33, p = 2.3e-10 variant: condition 1, z: movement 6.97 vs pupil 5.21; movement higher in 33 of 33, p = 2.3e-10 notebook run time from the first cell: 16.1 min
Result: Not reproduced. Under the claim's main reading, keeping only the movement index and its SD gives 68.6% against 71.5% for the pupil alone, higher in 14 of 33 classifications (p = 0.086; z p = 0.068): the movement index is not better than pupil. The code's conditions shown as variants give the movement index higher accuracy (condition 5: 75.7%, p = 1.3e-05; condition 1: 77.9%, p = 2.3e-10), but they keep pupil information as well and do not test the claim.
Summary¶
| Claim | Outcome |
|---|---|
| cohort-inclusion | partially reproduced (13 mice, 8 and 5 by task; 376 sessions pass, not 381) |
| training-duration | blocked: data missing |
| distractor-response-rate | reproduced (5.2% overall; one mouse at 16%) |
| figs2-n-states | reproduced (3 or 4 states; no mouse with 5) |
| fig1g-example-mouse-prediction | reproduced (4 states; 70.2% and 78.1%) |
| fig2a-state-types | reproduced, resting on an assumption for one state of BW059 |
| fig2b-optimal-occupancy | reproduced (indeterminate 13.1%) |
| fig2c-dwell-times | partially reproduced (tests and indeterminate dwell hold; optimal and sub-optimal dwells longer than 37 and 15) |
| fig2d-prediction-accuracy | reproduced (76.2%, gain 11.0 points) |
| fig3c-inverted-u-pupil | reproduced |
| fig3d-optimal-pupil-modality | reproduced (p = 0.013) |
| fig3e-u-disengaged-pupil | reproduced |
| figs4-r2-individual-movement | partially reproduced (lower on average; values and tests do not all hold under one reading) |
| fig3f-movement-index-inverted-u | reproduced (r2 0.574, close to the lower bound 0.57) |
| fig3f-movement-index-u-disengaged | reproduced |
| figs5-movement-index-correlations | reproduced |
| fig3-distance-from-optimal-pupil | partially reproduced (order and tests hold; values larger) |
| fig3-raw-face-energy-trend | reproduced |
| fig3-raw-locomotion-no-difference | reproduced |
| fig4a-lower-variability | partially reproduced (5 of 6 tests) |
| fig4b-variability-at-transitions | partially reproduced (signs hold; 3 of 4 marked tests) |
| fig4d-pupil-movement-decoupling | reproduced |
| figs6-individual-measure-decoupling | reproduced |
| fig5c-state-classification | reproduced (88.9% and 78.4%) |
| figs7-pupil-best-raw-measure | reproduced |
| fig5c-movement-index-best | not reproduced (movement index 68.6% against pupil 71.5%) |
| figs8a-regression-weights | reproduced |
| figs8b-pupil-movement-interaction | reproduced |
| figs1-iti-independence | reproduced |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. The GLM-HMM of the paper's ssm fork is reimplemented (hulsey.py), fit by the same EM with the same priors, initialisation and stopping rule; its random initial values differ from ssm's. The decoders use 3 shuffle runs per shuffled condition and 20 null runs, not 1000, to keep the notebook within its time limit. The sessions are the 376 that pass the paper's inclusion criteria, where the paper reports 381. Where the paper's state-typing code would score a state with no discrete trials as equal to every template, that state is typed on all sessions (one state of BW059).
What was found about the dataset. The dandiset holds the 13 included mice and 391 sessions of stages 5 and 6, more than the 381 the paper analysed, with no way to tell which were left out; the excluded mice, the other stage 5 and 6 sessions and the training sessions before stage 5 are not in it. running_speed is labelled cm/s but holds m/s. 38 of the 391 sessions have a post hoc pupil series that is entirely missing. The rig names are spelled several ways ('Box 1', 'Box1', 'Boox1'). No recording is deposited twice, and the measures share the clock of the trial table (data check).