Reproducing Davatolhagh et al. (2025) from dandiset 001712¶
Paper. Davatolhagh MF, Couto J, Melin M, Oesch LT, Findling C, Morales P, Hong I, International Brain Laboratory, Churchland AK (2025). Altered use of prior expectations and modified neural dynamics in a mouse model of autism. bioRxiv, version 1. doi:10.1101/2025.09.12.675910. This attempt used the full text of the preprint (CC-BY-NC-ND).
Dandiset. DANDI:001712, version draft. The
paper cites no version, and the dandiset has no published version: the archive lists only the
draft. The draft is used, with its modification time (2026-05-28T10:04:28Z) and the ids and sizes
of the three assets read, which the next cell checks against the archive. A draft can change; if
the check reports a difference, the numbers below may not match the attempt record.
What the dandiset holds. One session (2023-11-14) of one mouse, FD-28, in the biased choice-world task with widefield imaging: a processed NWB file, a raw NWB file and three camera videos. The paper analyses 11 wild-type and 8 Nrxn1a knockout mice with 10 sessions each. The genotype of FD-28 is not stated in the files or the dandiset's metadata, nor in the public IBL database (queried by hand during the attempt, not in this notebook). The findings that compare genotypes cannot be tested; those that the paper states for mice of both genotypes are tested on this session.
Code read. The paper links churchlandlab/chiCa for its encoding analyses (read at commit
d98b3c6) and jcouto/wfield for preprocessing (commit 0befe16). The dandiset links the conversion
code catalystneuro/IBL-widefield-to-nwb (commit 4a6fb06). Choices that follow them are given in
the helper module wfprior.py.
| Claim | Data | Status |
|---|---|---|
| fig1c-psychometric: rightward choices rise with right contrast | dandiset | tested |
| fig1d-psychometric-params-genotype | dandiset | blocked: data missing |
| fig1e-rt-contrast: longer reaction times at lower contrast | dandiset | tested |
| fig1f-rt-genotype | dandiset | blocked: data missing |
| fig1g-learning-genotype | dandiset | blocked: data missing |
| fig2bc-prior-shift: more rightward choices in right-likely blocks | dandiset | tested |
| fig2bd-ko-smaller-shift | dandiset | blocked: data missing |
| fig2ef-block-adaptation: choices adapt after a block switch | dandiset | tested |
| fig2g-ko-slower-adaptation | dandiset | blocked: data missing |
| fig3b-contralateral-visual: lateralized visual responses | dandiset | tested |
| fig3c-ko-elevated-activity | dandiset | blocked: data missing |
| fig3d-ko-higher-correlation | dandiset | blocked: data missing |
| fig3e-ko-higher-variance | dandiset | blocked: data missing |
| fig4c-full-cvr2: the encoding model explains over half the variance | dandiset | tested |
| fig4cd-movement-dominates: movement predictors dominate | dandiset | tested |
| fig4cd-ko-more-movement | dandiset | blocked: data missing |
| figs2-ko-elevated-motion | dandiset | blocked: data missing |
| figs3-ko-higher-unique | dandiset | blocked: data missing |
In [1]:
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import optimize, stats
import wfprior as wf
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
info = wf.draft_info()
print(f"dandiset {wf.DANDISET_ID} version {wf.DANDISET_VERSION}, modified {info['modified']}")
for path, (aid, size) in sorted(info["assets"].items()):
print(f" {size / 1e9:7.2f} GB {aid} {path}")
diff = wf.check_assets(info)
print("Assets differ from those this entry was written against:" if diff else
"The assets read match those this entry was written against.", *diff, sep="\n")
f = wf.open_nwb(wf.PROCESSED)
trials = wf.load_trials(f)
svd = wf.load_svd(f)
frame_t = svd["t"]
movement = wf.load_movement(f, frame_t)
print(f"loaded trials, imaging components and movement in {time.time() - started:.0f} s")
dandiset 001712 version draft, modified 2026-05-28 10:04:28.686333+00:00
1.05 GB 3102badd-5260-498a-aaaf-acec7d7960e7 sub-FD-28/sub-FD-28_ses-81f90b18-e61c-4d32-bbce-3e0c5f33f06c_desc-processed_behavior+ophys.nwb
33.35 GB d99e80d1-1b9e-4a83-836c-dec4a46e4cce sub-FD-28/sub-FD-28_ses-81f90b18-e61c-4d32-bbce-3e0c5f33f06c_desc-raw_behavior+ophys.nwb
0.61 GB 07e84de1-50e0-443d-9ca0-caa5d7a2636d sub-FD-28/sub-FD-28_ses-81f90b18-e61c-4d32-bbce-3e0c5f33f06c_ecephys+image/sub-FD-28_ses-81f90b18-e61c-4d32-bbce-3e0c5f33f06c_VideoBodyCamera.mp4
3.50 GB 8ee6d630-f8d1-4400-aa82-782341c16309 sub-FD-28/sub-FD-28_ses-81f90b18-e61c-4d32-bbce-3e0c5f33f06c_ecephys+image/sub-FD-28_ses-81f90b18-e61c-4d32-bbce-3e0c5f33f06c_VideoLeftCamera.mp4
1.10 GB 18493935-8c8c-4f8b-b04f-802e8a7cdef5 sub-FD-28/sub-FD-28_ses-81f90b18-e61c-4d32-bbce-3e0c5f33f06c_ecephys+image/sub-FD-28_ses-81f90b18-e61c-4d32-bbce-3e0c5f33f06c_VideoRightCamera.mp4
The assets read match those this entry was written against.
loaded trials, imaging components and movement in 23 s
Data check¶
Each check of the metadata against the values is stated under its output.
In [2]:
import re
meta = {k: wf.text(f, "general/" + k) for k in ("lab", "institution", "protocol", "session_id", "keywords")}
subj = {k: wf.text(f, "general/subject/" + k) for k in ("subject_id", "sex", "strain", "date_of_birth", "species")}
print("Session:", wf.text(f, "session_start_time"), "| lab:", meta["lab"], "| protocol:", meta["protocol"])
print("Subject:", subj)
print("Subject fields:", sorted(f["general/subject"].keys()), "(no genotype field)")
print("Keywords in the file:", meta["keywords"])
print("Keywords in the dandiset's metadata (DANDI API):", info["keywords"])
print("Gene names in the dandiset's description:", sorted(set(re.findall(r"CNTN\w*|Cntn\w*|Nrxn\w*|NRXN\w*", info["description"]))))
for plane in ("ImagingPlaneCalcium", "ImagingPlaneIsosbestic"):
p = f["general/optophysiology/" + plane]
print(plane, "| indicator:", wf.text(f, f"general/optophysiology/{plane}/indicator"),
"| excitation:", p["excitation_lambda"][()], "nm | location:", wf.text(f, f"general/optophysiology/{plane}/location"))
print("Assets: files per subject", pd.Series([p.split("/")[0] for p in info["assets"]]).value_counts().to_dict())
Session: 2023-11-14T12:10:33.757227-08:00 | lab: Churchland_ucla | protocol: _iblrig_tasks_biasedChoiceWorld
Subject: {'subject_id': 'FD-28', 'sex': 'F', 'strain': 'C57BL/6', 'date_of_birth': '2022-12-08T00:00:00-08:00', 'species': 'Mus musculus'}
Subject fields: ['date_of_birth', 'description', 'expected_water_ml', 'remaining_water_ml', 'sex', 'species', 'strain', 'subject_id', 'uuid', 'weight'] (no genotype field)
Keywords in the file: ['widefield calcium imaging', 'decision-making', 'prior information', 'autism spectrum disorder', 'Neurexin1alpha', 'Nrxn1alpha knockout', 'cortical dynamics', 'GCaMP6s', 'dorsal cortex', 'synaptic signaling']
Keywords in the dandiset's metadata (DANDI API): ['CNTNALP2', 'widefield imaging', 'decision-making', 'prior expectations', 'sensory perception', 'calcium imaging']
Gene names in the dandiset's description: ['CNTNALP2']
ImagingPlaneCalcium | indicator: GCaMP6f | excitation: 470.0 nm | location: Whole Brain
ImagingPlaneIsosbestic | indicator: GCaMP6f | excitation: 405.0 nm | location: Whole Brain
Assets: files per subject {'sub-FD-28': 5}
One subject and one session. The subject record has no genotype, and its description is the paper's general description of the cross (KO and wild-type littermates). The paper (Methods) and the file's keywords give the indicator as GCaMP6s, while both imaging planes say GCaMP6f. The dandiset's own keywords and description name "CNTNALP2" (printed above), a different gene from the Nrxn1a of the paper and the file.
In [3]:
print("Trials:", len(trials))
for c in ["gabor_stimulus_contrast", "gabor_stimulus_side", "mouse_wheel_choice", "block_type", "probability_left",
"is_mouse_rewarded", "reward_volume_uL"]:
print(f" {c}: {trials[c].value_counts(dropna=False).sort_index().to_dict()}")
print("Blocks:", trials.block_index.nunique(), "| first block length:", int((trials.block_index == 0).sum()),
"| biased block lengths:", trials[trials.block_index > 0].groupby("block_index").size().describe()[["min", "50%", "max"]].to_dict())
side_given_block = pd.crosstab(trials.block_type, trials.gabor_stimulus_side, normalize="index").round(3)
print("Stimulus side by block (fraction):", side_given_block.to_dict("index"))
ok = trials.mouse_wheel_choice != "none"
correct = ((trials.choice_right == 1) & (trials.gabor_stimulus_side == "right")) | ((trials.choice_right == 0) & (trials.gabor_stimulus_side == "left"))
zero = trials.gabor_stimulus_contrast == 0
print("Rewarded equals correct choice (counter-clockwise = right) on non-zero contrast trials:",
bool((trials.is_mouse_rewarded[~zero] == correct[~zero]).all()),
"| on 0% trials:", f"{(trials.is_mouse_rewarded[zero] == correct[zero]).mean():.3f}")
times = ["start_time", "quiescence_period", "gabor_stimulus_onset_time", "auditory_cue_time", "wheel_movement_onset_time",
"choice_registration_time", "feedback_time", "gabor_stimulus_offset_time", "stop_time"]
print(trials[times].describe().T[["count", "min", "max"]].round(3).to_string())
d = trials.gabor_stimulus_onset_time - trials.auditory_cue_time
print(f"Stimulus onset minus go cue (s): median {d.median():.4f}, range {d.min():.4f} to {d.max():.4f}")
q = trials.gabor_stimulus_onset_time - trials.start_time
print(f"Stimulus onset minus trial start (s): median {q.median():.3f}; quiescence period range {trials.quiescence_period.min():.3f} to {trials.quiescence_period.max():.3f}")
print(f"Reaction time (s): {trials.rt.notna().sum()} finite, {(trials.rt <= 0).sum()} at or before stimulus onset, "
f"median {trials.rt.median():.3f}, range {trials.rt.min():.3f} to {trials.rt.max():.3f}")
g = trials.wheel_movement_onset_time - trials.auditory_cue_time
print(f"Wheel movement onset minus go cue (s): {int((g < 0).sum())} of {int(g.notna().sum())} negative, minimum {g.min():.3f}")
print(f"Feedback minus choice registration (s): median {(trials.feedback_time - trials.choice_registration_time).median():.4f}")
Trials: 913
gabor_stimulus_contrast: {0.0: 178, 6.25: 195, 12.5: 184, 25.0: 167, 100.0: 189}
gabor_stimulus_side: {'left': 428, 'right': 485}
mouse_wheel_choice: {'clockwise': 308, 'counter_clockwise': 603, 'none': 2}
block_type: {'left_block': 367, 'right_block': 456, 'unbiased': 90}
probability_left: {0.2: 456, 0.5: 90, 0.8: 367}
is_mouse_rewarded: {False: 226, True: 687}
reward_volume_uL: {0.0: 226, 1.5: 687}
Blocks: 16 | first block length: 90 | biased block lengths: {'min': 9.0, '50%': 54.0, 'max': 96.0}
Stimulus side by block (fraction): {'left_block': {'left': 0.793, 'right': 0.207}, 'right_block': {'left': 0.215, 'right': 0.785}, 'unbiased': {'left': 0.433, 'right': 0.567}}
Rewarded equals correct choice (counter-clockwise = right) on non-zero contrast trials: True | on 0% trials: 1.000
count min max
start_time 913.0 9.005 3993.066
quiescence_period 913.0 0.400 0.700
gabor_stimulus_onset_time 913.0 9.612 3993.582
auditory_cue_time 913.0 9.527 3993.583
wheel_movement_onset_time 899.0 9.547 3993.647
choice_registration_time 913.0 9.824 3993.833
feedback_time 913.0 9.824 3993.834
gabor_stimulus_offset_time 913.0 10.878 3995.899
stop_time 913.0 11.378 3996.399
Stimulus onset minus go cue (s): median -0.0010, range -0.0015 to 0.0844
Stimulus onset minus trial start (s): median 0.675; quiescence period range 0.400 to 0.700
Reaction time (s): 899 finite, 149 at or before stimulus onset, median 0.069, range -0.182 to 16.771
Wheel movement onset minus go cue (s): 151 of 899 negative, minimum -0.183
Feedback minus choice registration (s): median 0.0001
The contrasts are in percent (0 to 100) and the five values of the description occur with similar counts. The first block is the 90 unbiased trials the description gives, and the stimulus side follows the block probabilities (about 0.8 on the likely side). A counter-clockwise turn is the correct response to a right stimulus on every trial with a visible stimulus, which confirms the coding of rightward choices. The go cue and the photodiode-detected stimulus onset are close to simultaneous, as the description says (offsets printed above). The wheel movement onset is described as the first movement after the go cue, but some onsets fall before the go cue and the stimulus onset (counts above); those trials are left out of reaction-time measures.
In [4]:
dt = np.diff(frame_t)
print(f"Imaging: {len(frame_t)} frames per channel, {frame_t[0]:.2f} to {frame_t[-1]:.2f} s, "
f"interval median {np.median(dt) * 1000:.3f} ms ({1 / np.median(dt):.2f} Hz), range {dt.min() * 1000:.3f} to {dt.max() * 1000:.3f} ms")
print("Trials span", round(trials.start_time.min(), 2), "to", round(trials.stop_time.max(), 2), "s, inside the imaging")
raw = wf.open_nwb(wf.RAW)
raw_t = raw["acquisition/OnePhotonSeriesCalcium/timestamps"][:]
print("Raw 470 nm frame timestamps equal the processed component timestamps:", np.array_equal(raw_t, frame_t),
f"(raw frames {raw['acquisition/OnePhotonSeriesCalcium/data'].shape})")
print("Raw imaging rate stated:", raw["general/optophysiology/ImagingPlaneCalcium/imaging_rate"][()], "Hz")
u = svd["u"].reshape(len(svd["u"]), -1).astype(np.float64)
gram = u @ u.T
print(f"Spatial components: {svd['u'].shape}, Gram diagonal {np.diag(gram).min():.5f} to {np.diag(gram).max():.5f}, "
f"largest off-diagonal {np.abs(gram - np.diag(np.diag(gram))).max():.2e}, every pixel nonzero: {bool((np.abs(u).sum(0) > 0).all())}")
svt = svd["svt"]
recon = svt[::5000] @ u
print(f"Temporal components: {svt.shape}; reconstructed dF/F over {len(recon)} sampled frames: 1st to 99th percentile "
f"{np.percentile(recon, 1):.4f} to {np.percentile(recon, 99):.4f}")
print(f"Mean image: {svd['mean_image'].shape}, range {svd['mean_image'].min():.0f} to {svd['mean_image'].max():.0f}")
print("Lab metadata groups in the processed file:", sorted(f["general"].keys()))
Imaging: 124886 frames per channel, 1.42 to 4000.69 s, interval median 32.033 ms (31.22 Hz), range 32.000 to 32.033 ms Trials span 9.01 to 3996.4 s, inside the imaging
Raw 470 nm frame timestamps equal the processed component timestamps: True (raw frames (124886, 540, 640)) Raw imaging rate stated: 15.0 Hz
Spatial components: (200, 540, 640), Gram diagonal 1.00000 to 1.00026, largest off-diagonal 1.57e-05, every pixel nonzero: True Temporal components: (124886, 200); reconstructed dF/F over 25 sampled frames: 1st to 99th percentile -0.0276 to 0.0457 Mean image: (540, 640), range 534 to 50601 Lab metadata groups in the processed file: ['devices', 'experiment_description', 'experimenter', 'ibl_metadata', 'institution', 'keywords', 'lab', 'optophysiology', 'protocol', 'related_publications', 'session_id', 'source_script', 'subject']
The Methods give 30 frames per second per channel; the timestamps give 31.2 Hz, with regular
intervals, and the raw file's imaging plane states 15 Hz. The raw and the processed file hold the same recording (same timestamps), so it is
counted once. The spatial components are orthonormal, so variance in pixel space is the sum of the
variances of the temporal components. They cover every pixel of the frame: there is no brain mask.
The conversion code's documentation says the processed file holds landmarks and a registration to
the Allen atlas (lab_meta_data localization, atlas_registration); neither file holds them, so
no pixel can be assigned to a cortical area. The reconstructed values are of the order of a few
percent, consistent with dF/F.
In [5]:
print("Movement variables per imaging frame:")
print(movement.describe().T[["mean", "std", "min", "max"]].round(3).to_string())
for cam in ("Left", "Right", "Body"):
t = f[f"processing/motion_energy/{cam}CameraMotionEnergy/timestamps"][:]
print(f"{cam} camera: {len(t)} frames, {1 / np.median(np.diff(t)):.2f} Hz")
print("Units: wheel", f["processing/wheel/WheelVelocitySmoothed/data"].attrs["unit"],
"| paws", f["processing/pose_estimation/LeftCamera/PoseEstimationSeriesLeftPaw/data"].attrs["unit"],
"| motion energy", f["processing/motion_energy/BodyCameraMotionEnergy/data"].attrs["unit"])
for cam in ("Left", "Right"):
for paw in ("LeftPaw", "RightPaw"):
conf = f[f"processing/pose_estimation/{cam}Camera/PoseEstimationSeries{paw}/confidence"][:]
print(f" {cam} camera {paw}: fraction of frames with likelihood below 0.9: {(conf < 0.9).mean():.3f}")
# Where the ROI of each motion energy series lies, against the median position of the nose tip and the pupil
# centre (mean of the four pupil points) in the same camera, from frames with likelihood 0.9 or more.
import re
def median_xy(cam, part):
g = f[f"processing/pose_estimation/{cam}Camera/PoseEstimationSeries{part}"]
xy = g["data"][:]
return np.nanmedian(xy[g["confidence"][:] >= 0.9], axis=0)
for cam in ("Left", "Right", "Body"):
desc = f[f"processing/motion_energy/{cam}CameraMotionEnergy"].attrs["description"]
w_, h_, x0, y0 = map(int, re.search(r"(\d+) pixels wide, (\d+) pixels tall, top-left corner at \((\d+), (\d+)\)", desc).groups())
print(f"{cam} camera ROI: x {x0} to {x0 + w_}, y {y0} to {y0 + h_}")
if cam == "Body":
print(" tail start", median_xy(cam, "TailStart").round())
continue
nose = median_xy(cam, "NoseTip")
pupil = np.mean([median_xy(cam, f"RightPupil{s}") for s in ("Top", "Bottom", "Left", "Right")], axis=0)
print(f" nose tip {nose.round()}, pupil centre {pupil.round()}; ROI centre ({x0 + w_ / 2:.0f}, {y0 + h_ / 2:.0f}); "
f"ROI centre between the two in x: {min(nose[0], pupil[0]) < x0 + w_ / 2 < max(nose[0], pupil[0])}, "
f"below the pupil: {y0 + h_ / 2 > pupil[1]}")
Movement variables per imaging frame:
mean std min max
wheel_velocity -0.182 0.848 -7.026 6.478
paw_left_leftcam 343.419 546.761 0.183 8266.333
paw_right_leftcam 472.724 788.476 0.156 10671.467
paw_left_rightcam 298.151 420.719 0.684 4750.300
paw_right_rightcam 224.384 432.715 0.005 5833.166
motion_energy_left 7.642 3.415 1.505 117.863
motion_energy_right 3.500 2.049 0.764 47.627
motion_energy_body 2.420 2.366 0.348 67.980
Left camera: 239237 frames, 59.76 Hz
Right camera: 600397 frames, 150.00 Hz
Body camera: 119964 frames, 30.00 Hz
Units: wheel rad/s | paws px | motion energy a.u.
Left camera LeftPaw: fraction of frames with likelihood below 0.9: 0.115
Left camera RightPaw: fraction of frames with likelihood below 0.9: 0.004
Right camera LeftPaw: fraction of frames with likelihood below 0.9: 0.182
Right camera RightPaw: fraction of frames with likelihood below 0.9: 0.017
Left camera ROI: x 179 to 333, y 165 to 267
nose tip [110. 223.], pupil centre [398. 111.]; ROI centre (256, 216); ROI centre between the two in x: True, below the pupil: True Right camera ROI: x 478 to 564, y 104 to 161
nose tip [609. 106.], pupil centre [485. 62.]; ROI centre (521, 132); ROI centre between the two in x: True, below the pupil: True Body camera ROI: x 138 to 345, y 44 to 254
tail start [308. 221.]
The camera rates match the Methods (left 60 Hz, right 150 Hz, body 30 Hz). The ROI motion energy descriptions give the ROI's position in pixels but not the body part. In each side camera the ROI centre lies between the nose tip and the pupil and below the pupil (printed above), where the whisker pad is. The body camera's ROI contains the median position of the tail start, so it covers the body. These are taken to be the whisker pads and the body, as the paper's model uses; the positions agree with that, but the files do not name the body parts.
In [6]:
# Clock checks. (1) The wheel turns the way the choice says, after the movement onset of the trial table.
wt = movement.wheel_velocity.to_numpy()
on = wf.frame_index(frame_t, trials.wheel_movement_onset_time.fillna(trials.choice_registration_time).to_numpy())
lags = np.arange(-15, 16)
vel = np.stack([wt[on + k] for k in lags])
sign = np.where(trials.choice_right == 1, 1, -1)
signed = (vel * sign[None, :])[:, ok.to_numpy()]
print("Wheel velocity (rad/s) times choice sign (+1 right), median over trials, by frame from movement onset:")
print(pd.Series(np.median(signed, axis=1), index=lags).round(3).to_string())
print("Wheel position reference frame:", wf.text(f, "processing/wheel/WheelPosition/reference_frame"))
after = vel[15 + 3:15 + 10][:, ok.to_numpy()].mean(0)
ccw = (trials.mouse_wheel_choice == "counter_clockwise").to_numpy()[ok.to_numpy()]
print(f"Mean wheel velocity 0.1 to 0.3 s after movement onset: negative on {(after[ccw] < 0).mean():.3f} of counter-clockwise choices, "
f"positive on {(after[~ccw] > 0).mean():.3f} of clockwise choices")
# (2) Licks follow rewarded feedback.
licks = f["processing/lick_times/EventsLickTimes/timestamps"][:]
def lick_rate(t0, a, b):
return np.array([((licks >= t + a) & (licks < t + b)).sum() / (b - a) for t in t0])
fb = trials.feedback_time.to_numpy()
rew = trials.is_mouse_rewarded.to_numpy()
for name, (a, b) in {"-1 to 0 s": (-1, 0), "0 to 1 s": (0, 1)}.items():
print(f"Lick rate {name} from feedback: rewarded {lick_rate(fb[rew], a, b).mean():.2f}/s, "
f"unrewarded {lick_rate(fb[~rew], a, b).mean():.2f}/s")
Wheel velocity (rad/s) times choice sign (+1 right), median over trials, by frame from movement onset: -15 -0.001 -14 -0.001 -13 -0.002 -12 -0.002 -11 -0.001 -10 -0.002 -9 -0.002 -8 -0.002 -7 -0.002 -6 -0.003 -5 -0.003 -4 -0.002 -3 -0.003 -2 -0.005 -1 -0.014 0 -0.134 1 -0.369 2 -0.636 3 -1.001 4 -1.431 5 -1.570 6 -1.591 7 -1.589 8 -1.527 9 -1.407 10 -1.242 11 -1.003 12 -0.757 13 -0.552 14 -0.417 15 -0.270
Wheel position reference frame: Initial angle at session start is defined as zero. Counter-clockwise is positive. Mean wheel velocity 0.1 to 0.3 s after movement onset: negative on 0.955 of counter-clockwise choices, positive on 0.919 of clockwise choices Lick rate -1 to 0 s from feedback: rewarded 1.58/s, unrewarded 0.96/s Lick rate 0 to 1 s from feedback: rewarded 20.37/s, unrewarded 0.19/s
In [7]:
# (3) Imaging and the trial table. Two contrasts, each as the norm in pixel space of the difference between the mean
# responses of two trial groups, by frame from the event, baseline -0.6 to -0.3 s before it. Null: the labels permuted
# within strata (200 permutations), which keeps any difference that comes from the strata.
# feedback: rewarded against unrewarded trials, strata = absolute contrast;
# stimulus: right against left stimulus at 25 and 100% contrast, strata = block type.
strong = (trials.gabor_stimulus_contrast >= 25).to_numpy()
right = (trials.gabor_stimulus_side == "right").to_numpy()
onset = wf.frame_index(frame_t, trials.gabor_stimulus_onset_time.to_numpy())
feedback = wf.frame_index(frame_t, trials.feedback_time.to_numpy())
fps = 1 / np.median(np.diff(frame_t))
clk_lags = np.arange(-19, 17)
def clock_check(event_frames, label, strata, rng):
seg = np.stack([svt[o + clk_lags] - svt[o - 19:o - 9].mean(0) for o in event_frames])
def course(lab):
return np.linalg.norm(seg[lab].mean(0) - seg[~lab].mean(0), axis=1)
null = []
for _ in range(200):
perm = label.copy()
for s_ in np.unique(strata):
m = strata == s_
perm[m] = rng.permutation(label[m])
null.append(course(perm))
out = pd.DataFrame({"s": np.round(clk_lags / fps, 3), "true": course(label), "null_95th": np.percentile(null, 95, axis=0)}).set_index("s")
above = out.true > out.null_95th
run = above & above.shift(-1, fill_value=False) & above.shift(-2, fill_value=False)
first = out.index[(out.index > -0.3) & run.to_numpy()]
return out, int(above[out.index < 0].sum()), (first[0] if len(first) else None)
rng_clk = np.random.default_rng(7)
fb_course, fb_pre, fb_first = clock_check(feedback, trials.is_mouse_rewarded.to_numpy(), trials.gabor_stimulus_contrast.to_numpy(), rng_clk)
st_course, st_pre, st_first = clock_check(onset[strong], right[strong], trials.block_type.to_numpy()[strong], rng_clk)
both = pd.concat({"feedback: rewarded vs not": fb_course, "stimulus: right vs left": st_course}, axis=1)
print(f"Feedback: {int(trials.is_mouse_rewarded.sum())} rewarded, {int((~trials.is_mouse_rewarded).sum())} unrewarded trials; "
f"stimulus: {int((strong & right).sum())} right, {int((strong & ~right).sum())} left trials")
print(both.round(3).to_string())
for name, pre, first in (("feedback", fb_pre, fb_first), ("stimulus", st_pre, st_first)):
print(f"{name}: frames before the event above the null 95th percentile {pre} of {int((clk_lags < 0).sum())}; "
f"first of three consecutive frames above it at {first} s")
lick_lat = [licks[(licks > t) & (licks < t + 1)].min() - t for t in fb[rew] if ((licks > t) & (licks < t + 1)).any()]
print(f"First lick after rewarded feedback (video clock): median {np.median(lick_lat):.3f} s over {len(lick_lat)} trials")
rts = trials.rt[strong & ok.to_numpy() & (trials.rt > 0).to_numpy()]
print(f"Movement onset after the stimulus on the stimulus-contrast trials: median {rts.median():.3f} s")
Feedback: 687 rewarded, 226 unrewarded trials; stimulus: 196 right, 160 left trials
feedback: rewarded vs not stimulus: right vs left
true null_95th true null_95th
s
-0.609 0.274 0.630 0.461 0.760
-0.577 0.211 0.545 0.427 0.625
-0.545 0.242 0.382 0.404 0.499
-0.513 0.119 0.322 0.328 0.417
-0.481 0.076 0.273 0.158 0.427
-0.448 0.151 0.287 0.268 0.389
-0.416 0.234 0.348 0.333 0.434
-0.384 0.230 0.442 0.387 0.484
-0.352 0.201 0.475 0.269 0.572
-0.320 0.259 0.569 0.470 0.754
-0.288 0.337 0.603 0.526 0.890
-0.256 0.506 0.709 0.445 1.086
-0.224 0.619 0.777 0.557 1.261
-0.192 0.684 0.818 0.725 1.354
-0.160 0.830 0.929 0.717 1.351
-0.128 0.791 0.925 0.807 1.328
-0.096 0.708 0.990 0.892 1.383
-0.064 0.720 1.112 0.731 1.410
-0.032 0.691 1.241 0.884 1.645
0.000 0.678 1.289 0.897 1.772
0.032 0.773 1.399 1.063 1.994
0.064 0.667 1.493 0.975 2.008
0.096 1.003 1.472 0.923 1.891
0.128 1.174 1.525 0.739 1.656
0.160 1.336 1.567 0.920 1.487
0.192 1.509 1.564 1.117 1.490
0.224 1.845 1.594 1.285 1.460
0.256 2.235 1.535 1.435 1.435
0.288 2.532 1.598 1.621 1.572
0.320 3.145 1.650 1.823 1.639
0.352 3.904 1.726 2.030 1.686
0.384 4.441 1.810 2.153 1.712
0.416 5.063 1.869 2.248 1.839
0.448 5.587 1.945 2.337 2.034
0.481 5.907 1.921 2.397 2.100
0.513 6.089 1.922 2.507 2.125
feedback: frames before the event above the null 95th percentile 0 of 19; first of three consecutive frames above it at 0.224 s
stimulus: frames before the event above the null 95th percentile 0 of 19; first of three consecutive frames above it at 0.256 s
First lick after rewarded feedback (video clock): median 0.216 s over 663 trials
Movement onset after the stimulus on the stimulus-contrast trials: median 0.086 s
The wheel starts to turn at the trial table's movement onset, and licking rises after rewarded
feedback and not after unrewarded feedback, so the wheel, the lick times and the trial table share
a clock. The direction disagrees with the descriptions: the wheel series say counter-clockwise is
positive, but trials labelled counter_clockwise turn the wheel with negative velocity. The choice
labels agree with the rewards (above), so one of the two sign conventions is misdescribed; the
analysis uses the choice labels.
For both contrasts the difference stays below the null's 95th percentile on every frame before the event (0 of 19). After feedback the rewarded-unrewarded difference exceeds the null from 0.224 s, close to the median first lick after a reward on the video clock (0.216 s), as expected if licking drives the cortical signal and the imaging, the cameras and the trial table share a clock. After the stimulus the right-left difference exceeds the null only from 0.256 s, later than a visual latency; on these trials the median movement onset is 0.086 s, so the side difference may come from the choice movement. The checks are consistent with a shared clock to within a few frames; they do not show a visual response at its expected latency.
fig1c-psychometric¶
Mice of both genotypes show sigmoidal psychometric curves, with the proportion of rightward choices increasing as the stimulus contrast increases on the right.
In [8]:
ch = trials[ok]
psy = ch.groupby("signed_contrast").choice_right.agg(["mean", "size"]).rename(columns={"mean": "p_right", "size": "trials"})
print(psy.round(3).to_string())
rho, p_rho = stats.spearmanr(psy.index, psy.p_right)
span = psy.p_right.loc[100] - psy.p_right.loc[-100]
print(f"Spearman rho of P(right) with signed contrast over {len(psy)} levels: {rho:.3f} (p = {p_rho:.2g})")
print(f"P(right) at +100% minus at -100%: {span:.3f}")
fit_all = wf.fit_psychometric(ch.signed_contrast.to_numpy(), ch.choice_right.to_numpy())
print("Lapse-and-erf fit, all choice trials:", {k: round(float(v), 3) for k, v in fit_all.items()})
rng = np.random.default_rng(101)
boot = []
for _ in range(200):
s = ch.sample(len(ch), replace=True, random_state=rng.integers(1 << 31))
boot.append(wf.fit_psychometric(s.signed_contrast.to_numpy(), s.choice_right.to_numpy())["sigma"])
print(f"sigma (contrast %, inverse slope): {fit_all['sigma']:.2f}, bootstrap 95% interval {np.percentile(boot, 2.5):.2f} to {np.percentile(boot, 97.5):.2f}")
print("By block (variant):")
for b, g in ch.groupby("block_type"):
r = wf.fit_psychometric(g.signed_contrast.to_numpy(), g.choice_right.to_numpy())
print(f" {b}: n = {len(g)}, P(right) at -100 / +100: {g[g.signed_contrast == -100].choice_right.mean():.3f} / "
f"{g[g.signed_contrast == 100].choice_right.mean():.3f}, bias {r['bias']:.2f}, sigma {r['sigma']:.2f}")
p_right trials signed_contrast -100.00 0.112 80 -25.00 0.253 79 -12.50 0.461 89 -6.25 0.527 93 0.00 0.695 177 6.25 0.824 102 12.50 0.874 95 25.00 0.977 87 100.00 1.000 109 Spearman rho of P(right) with signed contrast over 9 levels: 1.000 (p = 0) P(right) at +100% minus at -100%: 0.887
Lapse-and-erf fit, all choice trials: {'bias': -6.769, 'sigma': 17.043, 'lapse_left': 0.116, 'lapse_right': 0.0, 'nll': 401.207, 'n': 911.0}
sigma (contrast %, inverse slope): 17.04, bootstrap 95% interval 13.98 to 20.47 By block (variant): left_block: n = 365, P(right) at -100 / +100: 0.077 / 1.000, bias 2.86, sigma 25.03
right_block: n = 456, P(right) at -100 / +100: 0.278 / 1.000, bias -16.54, sigma 18.39 unbiased: n = 90, P(right) at -100 / +100: 0.000 / 1.000, bias -10.37, sigma 24.24
In [9]:
fig, ax = plt.subplots(figsize=(4.5, 3.5))
cc = np.linspace(-100, 100, 201)
ax.plot(psy.index, psy.p_right, "ko")
ax.plot(cc, wf.psychometric(cc, *[fit_all[k] for k in ("bias", "sigma", "lapse_left", "lapse_right")]), "k-")
ax.set_xlabel("signed contrast (%)")
ax.set_ylabel("P(rightward choice)")
ax.set_ylim(0, 1)
ax.set_title("FD-28, all blocks")
plt.show()
print("Plotted P(right) range:", round(psy.p_right.min(), 3), "to", round(psy.p_right.max(), 3), "; none outside the axes")
Plotted P(right) range: 0.112 to 1.0 ; none outside the axes
Result: reproduced. The proportion of rightward choices rises at every step of signed contrast, from 0.112 at -100% to 1.000 at +100% (Spearman rho 1.000 over nine levels; span 0.887). The lapse-and-erf fit gives sigma 17.04% contrast (bootstrap 95% interval 13.98 to 20.47), a bias of -6.77% (toward rightward choices) and a left lapse of 0.116. This is one session of one mouse of unknown genotype.
fig1e-rt-contrast¶
Reaction time (time to first movement) is longer at lower contrasts in both genotypes.
In [10]:
rt_ok = ok & (trials.rt > 0)
r = trials[rt_ok]
print(f"Trials used: {len(r)} of {len(trials)} ({int((~ok).sum())} without a choice, {int((ok & ~(trials.rt > 0)).sum())} with no "
"movement onset after the stimulus)")
tab = r.groupby("gabor_stimulus_contrast").rt.describe()[["count", "25%", "50%", "75%"]]
print(tab.round(3).to_string())
rho, p = stats.spearmanr(r.gabor_stimulus_contrast, r.rt)
print(f"Spearman rho of |contrast| with reaction time: {rho:.3f}, p = {p:.2g}")
print(f"Median at 0%: {tab.loc[0.0, '50%']:.3f} s; at 100%: {tab.loc[100.0, '50%']:.3f} s")
for name, g in (("correct", r[correct[rt_ok]]), ("incorrect", r[~correct[rt_ok]])):
rr, pp = stats.spearmanr(g.gabor_stimulus_contrast, g.rt)
print(f"Variant, {name} trials only: n = {len(g)}, rho {rr:.3f}, p = {pp:.2g}, median RT {g.rt.median():.3f} s")
for lim in (0.08, 0.15):
g = r[r.rt > lim]
rr, pp = stats.spearmanr(g.gabor_stimulus_contrast, g.rt)
print(f"Variant, reaction times above {lim} s: n = {len(g)}, rho {rr:.3f}, p = {pp:.2g}")
bins = [-np.inf, 0, 0.08, 0.15, 0.4, np.inf]
print("Movement onset minus stimulus onset on choice trials, counts by interval (s):",
{str(k): int(v) for k, v in pd.cut(trials.rt[ok], bins).value_counts(sort=False).items()}, f"| no onset: {int(trials.rt[ok].isna().sum())}")
# Variant: response time, choice registration minus stimulus onset, on every choice trial.
resp = (trials.choice_registration_time - trials.gabor_stimulus_onset_time)[ok]
tab_r = resp.groupby(trials.gabor_stimulus_contrast[ok]).describe()[["count", "25%", "50%", "75%"]]
print("Response time (choice registration minus stimulus onset) by contrast:")
print(tab_r.round(3).to_string())
rr, pp = stats.spearmanr(trials.gabor_stimulus_contrast[ok], resp)
print(f"Variant, response time on all {int(ok.sum())} choice trials: rho {rr:.3f}, p = {pp:.2g}; median at 0% {tab_r.loc[0.0, '50%']:.3f} s, "
f"at 100% {tab_r.loc[100.0, '50%']:.3f} s; range {resp.min():.3f} to {resp.max():.3f} s")
rr, pp = stats.spearmanr(r.gabor_stimulus_contrast, resp[rt_ok[ok]])
print(f"Variant, response time on the {len(r)} trials of the criterion: rho {rr:.3f}, p = {pp:.2g}")
Trials used: 750 of 913 (2 without a choice, 161 with no movement onset after the stimulus)
count 25% 50% 75%
gabor_stimulus_contrast
0.00 147.0 0.034 0.089 0.172
6.25 157.0 0.035 0.091 0.173
12.50 152.0 0.034 0.095 0.174
25.00 136.0 0.031 0.089 0.195
100.00 158.0 0.032 0.081 0.146
Spearman rho of |contrast| with reaction time: -0.053, p = 0.15
Median at 0%: 0.089 s; at 100%: 0.081 s
Variant, correct trials only: n = 562, rho -0.050, p = 0.23, median RT 0.093 s
Variant, incorrect trials only: n = 188, rho -0.163, p = 0.026, median RT 0.077 s
Variant, reaction times above 0.08 s: n = 404, rho -0.091, p = 0.068
Variant, reaction times above 0.15 s: n = 226, rho -0.226, p = 0.00062
Movement onset minus stimulus onset on choice trials, counts by interval (s): {'(-inf, 0.0]': 149, '(0.0, 0.08]': 346, '(0.08, 0.15]': 178, '(0.15, 0.4]': 127, '(0.4, inf]': 99} | no onset: 12
Response time (choice registration minus stimulus onset) by contrast:
count 25% 50% 75%
gabor_stimulus_contrast
0.00 177.0 0.245 0.326 0.939
6.25 195.0 0.237 0.333 0.555
12.50 184.0 0.235 0.331 0.642
25.00 166.0 0.239 0.320 0.536
100.00 189.0 0.223 0.281 0.380
Variant, response time on all 911 choice trials: rho -0.134, p = 5.2e-05; median at 0% 0.326 s, at 100% 0.281 s; range 0.069 to 34.794 s
Variant, response time on the 750 trials of the criterion: rho -0.148, p = 4.5e-05
In [11]:
fig, axs = plt.subplots(1, 2, figsize=(9, 3.2))
edges = np.arange(-0.2, 1.0001, 0.02)
axs[0].hist(trials.rt[ok].dropna(), bins=edges, color="0.4")
axs[0].axvline(0, color="r", lw=0.8); axs[0].set_xlabel("movement onset minus stimulus onset (s)"); axs[0].set_ylabel("trials")
axs[1].hist(resp, bins=edges, color="0.4")
axs[1].axvline(0, color="r", lw=0.8); axs[1].set_xlabel("choice registration minus stimulus onset (s)")
plt.show()
for name, v in (("movement onset", trials.rt[ok].dropna()), ("choice registration", resp)):
print(f"{name}: {int(((v < edges[0]) | (v > edges[-1])).sum())} of {len(v)} trials outside -0.2 to 1.0 s")
movement onset: 48 of 899 trials outside -0.2 to 1.0 s choice registration: 144 of 911 trials outside -0.2 to 1.0 s
Result: not reproduced. Over 750 trials the Spearman correlation between absolute contrast and reaction time is -0.053 (p = 0.15), so the test of the criterion fails; the median at 0% (0.089 s) does exceed that at 100% (0.081 s). Reaction times are short at every contrast (medians 0.081 to 0.095 s). Of the 899 choice trials with a movement onset, 149 have it at or before the stimulus onset and 346 within 0.08 s after it, so many onsets in the trial table may mark movements begun before the stimulus could guide them. The result depends on this. With reaction times above 0.15 s only (226 trials) the correlation is -0.226 (p = 0.00062). The response time, from stimulus onset to choice registration, is shorter at higher contrast over all 911 choice trials (rho -0.134, p = 5.2e-05; medians 0.326 s at 0% and 0.281 s at 100%). These variants do not change the outcome, as the criterion and the paper's "time to first movement" use the movement onset. The paper's text does not say how the first movement is detected or how movements before the stimulus are treated.
fig2bc-prior-shift¶
Mice use the block prior; rightward choices are more frequent in right-likely blocks than in left-likely blocks, most clearly when contrast is low, shifting the psychometric curve toward the likely side.
In [12]:
biased = trials[ok & (trials.block_type != "unbiased")]
tab = biased.pivot_table(index="signed_contrast", columns="block_type", values="choice_right", aggfunc=["mean", "size"])
shift = pd.DataFrame({"p_right_left_block": tab[("mean", "left_block")], "p_right_right_block": tab[("mean", "right_block")],
"n_left_block": tab[("size", "left_block")], "n_right_block": tab[("size", "right_block")]})
shift["right_minus_left"] = shift.p_right_right_block - shift.p_right_left_block
print(shift.round(3).to_string())
z = biased[biased.signed_contrast == 0]
a = int(z[(z.block_type == "right_block")].choice_right.sum()); n_a = int((z.block_type == "right_block").sum())
b = int(z[(z.block_type == "left_block")].choice_right.sum()); n_b = int((z.block_type == "left_block").sum())
odds, p0 = stats.fisher_exact([[a, n_a - a], [b, n_b - b]], alternative="greater")
print(f"0% contrast: P(right) {a}/{n_a} = {a / n_a:.3f} in right-likely blocks, {b}/{n_b} = {b / n_b:.3f} in left-likely blocks; "
f"difference {a / n_a - b / n_b:.3f}, one-sided Fisher p = {p0:.2g}")
fits = {bt: wf.fit_psychometric(g.signed_contrast.to_numpy(), g.choice_right.to_numpy()) for bt, g in trials[ok].groupby("block_type")}
print(pd.DataFrame(fits).T[["bias", "sigma", "lapse_left", "lapse_right", "n"]].astype(float).round(3).to_string())
print(f"Bias, left-likely minus right-likely block: {fits['left_block']['bias'] - fits['right_block']['bias']:.2f} % contrast")
low = biased[biased.gabor_stimulus_contrast <= 12.5]
for bt in ("right_block", "left_block"):
g = low[low.block_type == bt]
print(f"Contrast 12.5% or less, {bt}: P(right) {g.choice_right.mean():.3f} (n = {len(g)})")
p_right_left_block p_right_right_block n_left_block n_right_block right_minus_left signed_contrast -100.00 0.077 0.278 52 18 0.201 -25.00 0.164 0.500 55 18 0.336 -12.50 0.333 0.792 57 24 0.458 -6.25 0.431 0.842 65 19 0.411 -0.00 0.507 0.831 73 89 0.325 6.25 0.562 0.926 16 68 0.364 12.50 0.611 0.957 18 69 0.345 25.00 0.818 1.000 11 68 0.182 100.00 1.000 1.000 18 83 0.000 0% contrast: P(right) 74/89 = 0.831 in right-likely blocks, 37/73 = 0.507 in left-likely blocks; difference 0.325, one-sided Fisher p = 9.2e-06
bias sigma lapse_left lapse_right n left_block 2.863 25.034 0.072 0.0 365.0 right_block -16.543 18.392 0.285 0.0 456.0 unbiased -10.372 24.241 0.000 0.0 90.0 Bias, left-likely minus right-likely block: 19.41 % contrast Contrast 12.5% or less, right_block: P(right) 0.885 (n = 269) Contrast 12.5% or less, left_block: P(right) 0.454 (n = 229)
In [13]:
fig, ax = plt.subplots(figsize=(4.5, 3.5))
for bt, colour in (("left_block", "tab:blue"), ("right_block", "tab:orange"), ("unbiased", "0.5")):
g = trials[ok & (trials.block_type == bt)].groupby("signed_contrast").choice_right.mean()
ax.plot(g.index, g.values, "o", color=colour, label=bt)
ax.plot(cc, wf.psychometric(cc, *[fits[bt][k] for k in ("bias", "sigma", "lapse_left", "lapse_right")]), color=colour)
ax.set_xlabel("signed contrast (%)")
ax.set_ylabel("P(rightward choice)")
ax.set_ylim(0, 1)
ax.legend(fontsize=8)
plt.show()
print("All plotted proportions lie in 0 to 1; none outside the axes")
All plotted proportions lie in 0 to 1; none outside the axes
Result: reproduced. At 0% contrast the mouse chose right on 0.831 of trials (74 of 89) in right-likely blocks and 0.507 (37 of 73) in left-likely blocks, a difference of 0.325 (one-sided Fisher p = 9.2e-06). The difference is 0.182 to 0.458 at every signed contrast except +100%, where both blocks give 1.000, and the psychometric bias moves by 19.41% contrast between left-likely (2.86) and right-likely (-16.54) blocks. Whether the shift is smaller than in wild-type mice cannot be tested on one mouse of unknown genotype.
fig2ef-block-adaptation¶
After a block switch, mice adapt their choices toward the newly likely side, and the proportion of choices toward the block side rises from an initial value A to an asymptote B with a time constant tau.
In [14]:
def switch_curve(sel):
rows = []
for bi in sorted(trials.block_index.unique()):
if bi < 2: # switches between biased blocks: from block 1 to 2 onward
continue
g = trials[(trials.block_index == bi)]
toward_right = g.block_type.iloc[0] == "right_block"
for k, (i, t) in enumerate(g.iterrows()):
if k >= 40:
break
if not np.isfinite(t.choice_right) or not sel(t):
continue
rows.append({"block": bi, "t": k, "toward": float(t.choice_right == (1.0 if toward_right else 0.0))})
return pd.DataFrame(rows)
def model(t, A, B, tau):
return B - (B - A) * np.exp(-t / tau)
def adaptation(df, label):
curve = df.groupby("t").toward.mean()
(A, B, tau), cov = optimize.curve_fit(model, curve.index.to_numpy(float), curve.to_numpy(), p0=[0.5, 0.8, 5],
bounds=([0, 0, 0.01], [1, 1, 200]))
early = df[df.t <= 4].toward; late = df[(df.t >= 20) & (df.t <= 39)].toward
_, p = stats.fisher_exact([[late.sum(), len(late) - late.sum()], [early.sum(), len(early) - early.sum()]], alternative="greater")
print(f"{label}: {df.block.nunique()} switches, {len(df)} trials; fit A = {A:.3f}, B = {B:.3f}, tau = {tau:.2f} trials "
f"(SE {np.sqrt(np.diag(cov))[2]:.2f}); toward block side trials 0-4 {early.mean():.3f} (n = {len(early)}), "
f"trials 20-39 {late.mean():.3f} (n = {len(late)}), one-sided Fisher p = {p:.2g}")
return curve, (A, B, tau)
sw_all = switch_curve(lambda t: True)
print("Block lengths of the switches used:", trials[trials.block_index >= 2].groupby("block_index").size().tolist())
curve_all, par_all = adaptation(sw_all, "All contrasts")
curve_low, par_low = adaptation(switch_curve(lambda t: t.gabor_stimulus_contrast <= 12.5), "Variant, contrast 12.5% or less")
print("Per switch, fraction toward the block side in trials 0-4 and 20-39:")
print(sw_all.assign(phase=np.where(sw_all.t <= 4, "0-4", np.where(sw_all.t >= 20, "20-39", "5-19")))
.pivot_table(index="block", columns="phase", values="toward", aggfunc="mean").round(2).to_string())
Block lengths of the switches used: [57, 80, 53, 54, 86, 38, 33, 59, 26, 77, 47, 43, 65, 9] All contrasts: 14 switches, 506 trials; fit A = 0.480, B = 1.000, tau = 25.12 trials (SE 16.14); toward block side trials 0-4 0.543 (n = 70), trials 20-39 0.827 (n = 237), one-sided Fisher p = 2.7e-06 Variant, contrast 12.5% or less: 14 switches, 316 trials; fit A = 0.354, B = 1.000, tau = 21.57 trials (SE 11.64); toward block side trials 0-4 0.452 (n = 42), trials 20-39 0.820 (n = 150), one-sided Fisher p = 5.4e-06 Per switch, fraction toward the block side in trials 0-4 and 20-39: phase 0-4 20-39 5-19 block 2 0.4 0.70 0.73 3 0.8 1.00 0.87 4 0.2 0.80 0.60 5 1.0 0.95 1.00 6 0.2 0.80 0.33 7 0.6 0.67 0.93 8 0.4 0.54 0.40 9 0.6 0.95 1.00 10 0.0 0.50 0.53 11 0.8 1.00 0.73 12 0.2 0.95 0.53 13 1.0 0.95 0.87 14 0.6 0.60 0.20 15 0.8 NaN 0.50
In [15]:
# Bootstrap over switches: switches drawn with replacement and the curve refitted (500 draws).
rng_sw = np.random.default_rng(202)
blocks_ = sw_all.block.unique()
boot_par = []
for _ in range(500):
d = pd.concat([sw_all[sw_all.block == b] for b in rng_sw.choice(blocks_, len(blocks_))])
c = d.groupby("t").toward.mean()
try:
boot_par.append(optimize.curve_fit(model, c.index.to_numpy(float), c.to_numpy(), p0=[0.5, 0.8, 5], bounds=([0, 0, 0.01], [1, 1, 200]))[0])
except RuntimeError:
boot_par.append([np.nan] * 3)
boot_par = pd.DataFrame(boot_par, columns=["A", "B", "tau"])
print(f"Bootstrap over switches: {int(boot_par.tau.isna().sum())} of {len(boot_par)} fits failed")
print(boot_par.quantile([0.025, 0.5, 0.975]).round(3).to_string())
print(f"Fraction of draws with B > A: {(boot_par.B > boot_par.A).mean():.3f}; with B at its bound (0.999 or more): {(boot_par.B >= 0.999).mean():.3f}; "
f"with tau of 40 trials or more: {(boot_par.tau >= 40).mean():.3f}")
# Variant: maximum-likelihood fit on the single trials (Bernoulli), which weights each trial number by its count.
t_, y_ = sw_all.t.to_numpy(float), sw_all.toward.to_numpy()
def nll_ad(p):
q = np.clip(model(t_, *p), 1e-6, 1 - 1e-6)
return -np.sum(y_ * np.log(q) + (1 - y_) * np.log(1 - q))
ml = min((optimize.minimize(nll_ad, [a0, 0.9, s0], method="L-BFGS-B", bounds=[(0, 1), (0, 1), (0.01, 200)])
for a0 in (0.3, 0.6) for s0 in (3, 10, 30)), key=lambda r_: r_.fun)
print(f"Variant, Bernoulli maximum-likelihood fit on {len(sw_all)} trials: A = {ml.x[0]:.3f}, B = {ml.x[1]:.3f}, tau = {ml.x[2]:.2f} trials")
print("Trials per trial number after the switch:", sw_all.groupby("t").size().agg(["min", "max"]).to_dict())
Bootstrap over switches: 0 of 500 fits failed
A B tau
0.025 0.279 0.836 6.410
0.500 0.460 1.000 23.870
0.975 0.643 1.000 38.369
Fraction of draws with B > A: 1.000; with B at its bound (0.999 or more): 0.632; with tau of 40 trials or more: 0.018
Variant, Bernoulli maximum-likelihood fit on 506 trials: A = 0.479, B = 1.000, tau = 25.03 trials
Trials per trial number after the switch: {'min': 10, 'max': 14}
In [16]:
fig, ax = plt.subplots(figsize=(5, 3.5))
tt = np.arange(40)
ax.plot(curve_all.index, curve_all.values, "ko", ms=3, label="all contrasts")
ax.plot(tt, model(tt, *par_all), "k-")
ax.plot(curve_low.index, curve_low.values, "o", color="tab:green", ms=3, label="contrast <= 12.5%")
ax.plot(tt, model(tt, *par_low), color="tab:green")
ax.set_xlabel("trials after block switch")
ax.set_ylabel("P(choice toward block side)")
ax.set_ylim(0, 1)
ax.legend(fontsize=8)
plt.show()
print("Plotted proportions range", round(min(curve_all.min(), curve_low.min()), 3), "to", round(max(curve_all.max(), curve_low.max()), 3))
Plotted proportions range 0.167 to 1.0
Result: reproduced. Over 14 switches between biased blocks, the proportion of choices toward the new block side is 0.543 in trials 0 to 4 and 0.827 in trials 20 to 39 (one-sided Fisher p = 2.7e-06). The fit gives A = 0.480, B = 1.000 (at its upper bound) and tau = 25.12 trials with a standard error of 16.14, so the time constant is poorly determined in one session. Resampling the switches (500 draws) gives 95% intervals of 0.279 to 0.643 for A, 0.836 to 1.000 for B and 6.41 to 38.37 trials for tau; B exceeds A in every draw, B is at its bound in 0.632 of them and tau is 40 trials or more in 0.018. A maximum-likelihood fit on the single trials, which weights each trial number by its count (10 to 14 trials), gives A = 0.479, B = 1.000 and tau = 25.03. At contrasts of 12.5% or less the pattern is the same (0.452 against 0.820; tau 21.57).
fig3b-contralateral-visual¶
Visual cortex (VISp) responds robustly to contralateral visual stimuli.
The files carry no atlas registration, so the test is in image coordinates (see the claims file).
In [17]:
mean_img = svd["mean_image"]
H, W = mean_img.shape
cols = np.arange(int(W * 0.3), int(W * 0.7))
midline = int(cols[np.argmin(mean_img[:, cols].mean(0))])
print(f"Midline estimate: column {midline} of {W} (lowest mean fluorescence in columns {cols[0]} to {cols[-1]})")
ev = wf.evoked(svt, onset) # trials x components
uim = svd["u"].reshape(len(svd["u"]), -1).astype(np.float64)
def peak_test(sel_trials):
rmap = (ev[sel_trials & strong & right].mean(0) - ev[sel_trials & strong & ~right].mean(0)) @ uim
res = {}
for name, m in (("right_minus_left", rmap), ("left_minus_right", -rmap)):
p = int(np.argmax(m)); y, x = divmod(p, W)
res[name] = {"row": y, "col": x, "side": "left of midline" if x < midline else "right of midline", "value": m[p], "pixel": p}
return rmap, res
all_trials = np.ones(len(trials), bool)
rmap, peaks = peak_test(all_trials)
for name, r in peaks.items():
vals = ev @ uim[:, r["pixel"]]
a, b = (vals[strong & right], vals[strong & ~right]) if name == "right_minus_left" else (vals[strong & ~right], vals[strong & right])
_, p = stats.mannwhitneyu(a, b, alternative="greater")
print(f"{name}: peak at row {r['row']}, column {r['col']} ({r['side']}), difference {r['value']:.4f} dF/F; "
f"defining side {np.median(a):.4f} against other side {np.median(b):.4f} (median dF/F), one-sided Mann-Whitney p = {p:.2g}")
print("Peaks on opposite sides of the midline:", peaks["right_minus_left"]["side"] != peaks["left_minus_right"]["side"])
print(f"Mean image at the two peaks: {mean_img.ravel()[peaks['right_minus_left']['pixel']]:.0f} and "
f"{mean_img.ravel()[peaks['left_minus_right']['pixel']]:.0f}; median over the frame {np.median(mean_img):.0f}, 5th percentile {np.percentile(mean_img, 5):.0f}")
for side, sel in (("right", strong & right), ("left", strong & ~right)):
m = ev[sel].mean(0) @ uim
p = int(np.argmax(m)); y, x = divmod(p, W)
print(f"Peak of the {side}-stimulus evoked map itself: row {y}, column {x}, {m[p]:.4f} dF/F")
Midline estimate: column 316 of 640 (lowest mean fluorescence in columns 192 to 447)
right_minus_left: peak at row 0, column 0 (left of midline), difference 0.0502 dF/F; defining side -0.0424 against other side -0.1129 (median dF/F), one-sided Mann-Whitney p = 0.053 left_minus_right: peak at row 539, column 2 (left of midline), difference 0.0222 dF/F; defining side 0.1086 against other side 0.0511 (median dF/F), one-sided Mann-Whitney p = 0.14 Peaks on opposite sides of the midline: False Mean image at the two peaks: 534 and 926; median over the frame 4223, 5th percentile 1538 Peak of the right-stimulus evoked map itself: row 539, column 0, 0.0697 dF/F Peak of the left-stimulus evoked map itself: row 539, column 0, 0.0898 dF/F
In [18]:
# Variant: peaks chosen on odd trials, tested on even trials (the main test chooses and tests on the same trials).
idx = np.arange(len(trials))
_, pk_odd = peak_test(idx % 2 == 1)
for name, r in pk_odd.items():
vals = ev @ uim[:, r["pixel"]]
sel_even = idx % 2 == 0
a, b = (vals[sel_even & strong & right], vals[sel_even & strong & ~right]) if name == "right_minus_left" else (vals[sel_even & strong & ~right], vals[sel_even & strong & right])
_, p = stats.mannwhitneyu(a, b, alternative="greater")
print(f"Odd-trial {name} peak: row {r['row']}, column {r['col']} ({r['side']}); even trials: {np.median(a):.4f} against {np.median(b):.4f}, p = {p:.2g}")
# Variant: correct trials only and 100% contrast only.
for label, sel in (("correct trials", correct.to_numpy()), ("100% contrast", (trials.gabor_stimulus_contrast == 100).to_numpy())):
m = (ev[sel & strong & right].mean(0) - ev[sel & strong & ~right].mean(0)) @ uim
pr, pl = int(np.argmax(m)), int(np.argmax(-m))
print(f"Variant, {label}: right-minus-left peak column {pr % W}, left-minus-right peak column {pl % W} (midline {midline})")
# Variant: peaks only among pixels brighter than the median of the mean image and 20 pixels or more from the frame edge.
yy, xx = np.mgrid[:H, :W]
mask = ((mean_img > np.median(mean_img)) & (yy >= 20) & (yy < H - 20) & (xx >= 20) & (xx < W - 20)).ravel()
print(f"Variant mask: {mask.sum()} of {mask.size} pixels")
masked_peaks = {}
for name, m in (("right_minus_left", rmap), ("left_minus_right", -rmap)):
p = int(np.flatnonzero(mask)[np.argmax(m[mask])]); y, x = divmod(p, W)
masked_peaks[name] = p
vals = ev @ uim[:, p]
a, b = (vals[strong & right], vals[strong & ~right]) if name == "right_minus_left" else (vals[strong & ~right], vals[strong & right])
_, pv = stats.mannwhitneyu(a, b, alternative="greater")
print(f" masked {name}: peak row {y}, column {x} ({'left' if x < midline else 'right'} of midline), {m[p]:.4f} dF/F; "
f"{np.median(a):.4f} against {np.median(b):.4f}, p = {pv:.2g}")
fig, axs = plt.subplots(1, 3, figsize=(12, 3.2))
axs[0].imshow(mean_img, cmap="gray"); axs[0].axvline(midline, color="y"); axs[0].set_title("mean image, midline")
lim = np.percentile(np.abs(rmap), 99.5)
im = axs[1].imshow(rmap.reshape(H, W), cmap="RdBu_r", vmin=-lim, vmax=lim); axs[1].axvline(midline, color="k", lw=0.5)
for name, r in peaks.items():
axs[1].plot(r["col"], r["row"], "k+" if name == "right_minus_left" else "kx", ms=10)
axs[1].set_title("right minus left stimulus, 0.1-0.4 s"); plt.colorbar(im, ax=axs[1], fraction=0.04)
pix = peaks["right_minus_left"]["pixel"]; pix2 = peaks["left_minus_right"]["pixel"]
for p_, lab in ((pix, "+ peak"), (pix2, "x peak")):
for sel, ls, stim in ((strong & right, "-", "right"), (strong & ~right, "--", "left")):
tr_ = np.stack([np.stack([svt[o + k] - svt[o - 9:o].mean(0) for k in lags]) @ uim[:, p_] for o in onset[sel]]).mean(0)
axs[2].plot(lags / fps, tr_, ls, label=f"{lab}, {stim} stimulus")
axs[2].axvline(0, color="0.5"); axs[2].set_xlabel("s from stimulus onset"); axs[2].set_ylabel("dF/F"); axs[2].legend(fontsize=7)
plt.show()
print(f"Map colour range clipped at +/-{lim:.4f} dF/F (99.5th percentile of |value|); {(np.abs(rmap) > lim).sum()} pixels beyond it")
Odd-trial right_minus_left peak: row 0, column 0 (left of midline); even trials: -0.0039 against -0.0834, p = 0.027 Odd-trial left_minus_right peak: row 392, column 479 (right of midline); even trials: 0.0171 against 0.0048, p = 2.7e-09 Variant, correct trials: right-minus-left peak column 0, left-minus-right peak column 2 (midline 316) Variant, 100% contrast: right-minus-left peak column 0, left-minus-right peak column 2 (midline 316) Variant mask: 172437 of 345600 pixels masked right_minus_left: peak row 459, column 199 (left of midline), 0.0119 dF/F; 0.0177 against 0.0037, p = 8.1e-15 masked left_minus_right: peak row 386, column 481 (right of midline), 0.0113 dF/F; 0.0172 against 0.0057, p = 1.6e-17
Map colour range clipped at +/-0.0098 dF/F (99.5th percentile of |value|); 1728 pixels beyond it
In [19]:
# Variant: the masked search on trials whose wheel movement starts more than 0.4 s after the stimulus, so that no
# movement falls in the 0.1 to 0.4 s window of the evoked map. The mask is the one above, set from the mean image only.
slow = (trials.rt > 0.4).to_numpy()
print(f"Slow trials at 25 or 100% contrast: {int((slow & strong & right).sum())} right, {int((slow & strong & ~right).sum())} left stimulus")
rmap_slow = (ev[slow & strong & right].mean(0) - ev[slow & strong & ~right].mean(0)) @ uim
for name, m in (("right_minus_left", rmap_slow), ("left_minus_right", -rmap_slow)):
p = int(np.flatnonzero(mask)[np.argmax(m[mask])]); y, x = divmod(p, W)
vals = ev @ uim[:, p]
a, b = (vals[slow & strong & right], vals[slow & strong & ~right]) if name == "right_minus_left" else (vals[slow & strong & ~right], vals[slow & strong & right])
_, pv = stats.mannwhitneyu(a, b, alternative="greater")
print(f" slow trials, masked {name}: peak row {y}, column {x} ({'left' if x < midline else 'right'} of midline), {m[p]:.4f} dF/F; "
f"{np.median(a):.4f} against {np.median(b):.4f}, p = {pv:.2g}")
for name, p in masked_peaks.items():
vals = ev @ uim[:, p]
a, b = (vals[slow & strong & right], vals[slow & strong & ~right]) if name == "right_minus_left" else (vals[slow & strong & ~right], vals[slow & strong & right])
_, pv = stats.mannwhitneyu(a, b, alternative="greater")
print(f" slow trials at the all-trial masked {name} peak (column {p % W}): {np.median(a):.4f} against {np.median(b):.4f}, p = {pv:.2g}")
# Time courses at the two all-trial masked peaks: defining-side stimulus minus other-side stimulus, baseline -0.3 to 0 s before
# the stimulus, aligned to the stimulus and to the movement onset, on trials at 25 or 100% with a movement onset after the stimulus.
moved = strong & (trials.rt > 0).to_numpy()
mo = wf.frame_index(frame_t, trials.wheel_movement_onset_time.fillna(trials.gabor_stimulus_onset_time).to_numpy())
tc_lags = np.arange(-10, 14)
rows_tc = {}
for name, p in masked_peaks.items():
base = np.array([svt[o - 9:o].mean(0) for o in onset]) @ uim[:, p]
for align, ev_frames in (("stimulus", onset), ("movement", mo)):
tr_ = np.stack([svt[o + tc_lags] @ uim[:, p] for o in ev_frames]) - base[:, None]
hi_side, lo_side = (moved & right, moved & ~right) if name == "right_minus_left" else (moved & ~right, moved & right)
diff = tr_[hi_side].mean(0) - tr_[lo_side].mean(0)
se = np.sqrt(tr_[hi_side].var(0, ddof=1) / hi_side.sum() + tr_[lo_side].var(0, ddof=1) / lo_side.sum())
rows_tc[(f"{name} col {p % W}", f"{align}: diff")] = diff
rows_tc[(f"{name} col {p % W}", f"{align}: z")] = diff / se
tc = pd.DataFrame(rows_tc, index=np.round(tc_lags / fps, 3))
tc.index.name = "s from event"
print(f"Trials: {int((moved & right).sum())} right, {int((moved & ~right).sum())} left stimulus; median movement onset {trials.rt[moved].median():.3f} s after the stimulus")
print(tc.round(4).to_string())
Slow trials at 25 or 100% contrast: 13 right, 14 left stimulus
slow trials, masked right_minus_left: peak row 270, column 550 (right of midline), 0.0343 dF/F; 0.0205 against -0.0136, p = 0.0062
slow trials, masked left_minus_right: peak row 364, column 72 (left of midline), 0.0277 dF/F; 0.0179 against -0.0244, p = 0.012
slow trials at the all-trial masked right_minus_left peak (column 199): 0.0147 against 0.0038, p = 0.02
slow trials at the all-trial masked left_minus_right peak (column 481): 0.0138 against 0.0025, p = 0.091
Trials: 162 right, 132 left stimulus; median movement onset 0.086 s after the stimulus
right_minus_left col 199 left_minus_right col 481
stimulus: diff stimulus: z movement: diff movement: z stimulus: diff stimulus: z movement: diff movement: z
s from event
-0.320 0.0006 0.5062 0.0036 2.7799 -0.0007 -0.7779 -0.0001 -0.1238
-0.288 0.0007 0.7513 0.0040 3.1975 -0.0003 -0.4031 -0.0002 -0.2556
-0.256 0.0007 1.0060 0.0036 2.8888 -0.0001 -0.2156 0.0002 0.2481
-0.224 0.0009 1.5254 0.0036 2.9301 -0.0003 -0.7030 0.0001 0.0717
-0.192 0.0004 0.9269 0.0034 2.6487 -0.0001 -0.3499 0.0002 0.1527
-0.160 -0.0000 -0.0888 0.0027 2.0503 0.0002 0.4798 0.0010 0.9009
-0.128 -0.0000 -0.1239 0.0035 2.4503 -0.0001 -0.1927 0.0008 0.6893
-0.096 -0.0006 -1.0166 0.0038 2.4881 0.0001 0.2932 0.0016 1.2511
-0.064 -0.0011 -1.5577 0.0039 2.3947 0.0005 0.8838 0.0028 1.9751
-0.032 -0.0008 -0.8969 0.0050 2.9163 0.0001 0.1772 0.0036 2.4093
0.000 -0.0006 -0.5036 0.0061 3.3999 -0.0001 -0.0550 0.0049 3.1375
0.032 -0.0008 -0.5848 0.0082 4.3485 -0.0002 -0.2199 0.0060 3.6513
0.064 -0.0002 -0.1581 0.0104 5.3440 0.0004 0.3095 0.0082 4.8305
0.096 0.0021 1.5150 0.0131 6.6520 0.0022 1.8512 0.0102 5.9770
0.128 0.0049 3.2980 0.0155 7.8488 0.0047 3.7236 0.0124 7.0913
0.160 0.0073 4.5457 0.0169 8.4819 0.0079 5.8278 0.0139 7.9738
0.192 0.0098 5.7695 0.0178 8.9154 0.0103 7.3039 0.0148 8.6172
0.224 0.0121 6.9217 0.0186 9.1566 0.0122 8.2610 0.0153 8.7202
0.256 0.0137 7.5639 0.0195 9.4459 0.0137 8.8074 0.0154 8.6561
0.288 0.0148 8.0398 0.0204 9.7053 0.0146 9.2550 0.0153 8.4459
0.320 0.0163 8.7111 0.0209 9.6209 0.0149 9.1463 0.0152 8.1114
0.352 0.0173 8.9188 0.0210 9.6420 0.0151 8.9292 0.0151 8.0278
0.384 0.0180 9.1728 0.0204 9.1794 0.0153 8.8502 0.0155 7.9247
0.416 0.0185 9.0764 0.0200 8.8582 0.0156 8.5046 0.0153 7.6472
In [20]:
# Variant: movements other than the wheel in the window of the evoked map. For each movement variable, its mean over 0.1 to
# 0.4 s after the stimulus minus its mean over -0.3 to 0 s (the windows of the imaging), compared between right and left
# stimulus trials at 25 or 100% contrast, on all of them and on the slow trials above (two-sided Mann-Whitney).
mv_ev = pd.DataFrame(wf.evoked(movement.to_numpy(), onset), columns=movement.columns)
rows_mv = []
for col in movement.columns:
for label, sel in (("all", strong), ("slow", slow & strong)):
a_, b_ = mv_ev[col][sel & right], mv_ev[col][sel & ~right]
rows_mv.append({"variable": col, "trials": label, "n right/left": f"{len(a_)}/{len(b_)}",
"median change, right stimulus": a_.median(), "median change, left stimulus": b_.median(),
"fraction rising": (mv_ev[col][sel] > 0).mean(), "p, right vs left": stats.mannwhitneyu(a_, b_).pvalue})
print(pd.DataFrame(rows_mv).set_index(["variable", "trials"]).round(4).to_string())
n right/left median change, right stimulus median change, left stimulus fraction rising p, right vs left
variable trials
wheel_velocity all 196/160 -1.7350 0.7222 0.3820 0.0000
slow 13/14 0.0000 -0.0000 0.5185 0.3958
paw_left_leftcam all 196/160 317.8015 289.6250 0.9382 0.5046
slow 13/14 21.5077 111.1769 0.7778 0.1389
paw_right_leftcam all 196/160 653.4738 865.7304 0.9551 0.0000
slow 13/14 25.7006 192.5895 0.6667 0.1040
paw_left_rightcam all 196/160 605.9570 646.6317 0.9045 0.9831
slow 13/14 -71.7499 2.8570 0.4444 0.5441
paw_right_rightcam all 196/160 592.3854 487.4757 0.9607 0.0005
slow 13/14 -11.2796 82.4070 0.7037 0.0143
motion_energy_left all 196/160 2.8773 4.3556 0.9129 0.0003
slow 13/14 -0.2376 0.6129 0.4815 0.1389
motion_energy_right all 196/160 2.0919 2.0238 0.9073 0.1091
slow 13/14 -0.0840 0.4230 0.4815 0.3693
motion_energy_body all 196/160 2.2043 2.1182 0.9045 0.0610
slow 13/14 0.0982 0.1020 0.5556 0.9806
Result: not reproduced. Under the criterion fixed beforehand, the peaks of the right-minus-left and left-minus-right maps both fall at the left edge of the frame (columns 0 and 2, rows 0 and 539), on the same side of the midline (column 316), with trial tests at p = 0.053 and 0.14. These corner pixels are darker in the mean image (534 and 926) than 95% of the frame (5th percentile 1538), so they lie outside the imaged brain; the spatial components cover the whole frame with no mask. In a variant restricted to pixels brighter than the median of the mean image and 20 pixels from the edge, the right-minus-left peak lies left of the midline (column 199) and the left-minus-right peak right of it (column 481), with p = 8.1e-15 and 1.6e-17.
The variants do not show that this lateralization is visual. On the trials at 25 or 100% contrast with a movement onset after the stimulus (162 right, 132 left), the side difference at the two masked peaks exceeds 3 standard errors from 0.128 s after the stimulus (z 3.3 and 3.7), while the median movement onset is 0.086 s. Aligned to the movement onset, the difference at the column 199 peak is already above zero 0.32 s before the onset (z 2.8), and at the column 481 peak from 0.03 s before it (z 2.4). Only 27 trials (13 right, 14 left) have no movement before 0.4 s. On them the difference at the two masked peaks points the same way (p = 0.02 and 0.091), but a masked peak search on these trials alone puts the right-minus-left peak right of the midline (column 550) and the left-minus-right peak left of it (column 72), the reverse of the search on all trials (p = 0.0062 and 0.012). This reversal is not explained here.
The slow trials exclude only wheel movement. Over the same 0.1 to 0.4 s window, paw speed and camera motion energy rise on most trials at 25 or 100% contrast, and several of them differ between right and left stimulus trials: wheel velocity, the right paw in both side cameras and the left camera's motion energy (p = 0.0005 or less). On the 27 slow trials, paw speed rises on 0.44 to 0.78 of trials and motion energy on 0.48 to 0.56, and the right paw in the right camera still differs by side (p = 0.014). So paw and face movements are not excluded on the slow trials, and over all trials they differ with the stimulus side in this window. Without an atlas registration neither the hemisphere nor VISp can be named. The variants do not change the outcome.
Encoding model¶
The design follows the claims file and wfprior.py. The body camera's video and motion-energy
components are computed here by streaming the video from the archive.
In [21]:
t0 = time.time()
vid = wf.video_components(f, frame_t, "Body")
print(f"Body camera: {vid['n_video_frames']} decoded frames, {vid['n_timestamps']} timestamps, {vid['pixels']} pixels after binning, "
f"in {time.time() - t0:.0f} s; top {wf.N_VIDEO_COMPONENTS} components explain {vid['video_explained']:.3f} of video variance "
f"and {vid['motion_explained']:.3f} of motion-energy variance")
ak = wf.fit_action_kernel(trials)
print(f"Action-kernel learning rate {ak['rate']:.2f} (negative log-likelihood {ak['nll']:.1f}; at rate 0.01 {ak['nll_grid'][0]:.1f}, at 0.99 {ak['nll_grid'][-1]:.1f}); "
f"weights (bias, contrast, sqrt contrast, prior) {np.round(ak['weights'], 2).tolist()}")
print(f"Correlation of subjective prior with block (0.5 - probability_left): {np.corrcoef(ak['prior'], 0.5 - trials.probability_left)[0, 1]:.3f}")
valid = ok.to_numpy()
rows_frame, rows_tp, rows_trial = wf.segment_frames(trials, frame_t, valid)
points = wf.STIM_WINDOW[1] - wf.STIM_WINDOW[0] + wf.FEEDBACK_WINDOW[1] - wf.FEEDBACK_WINDOW[0]
mov_basic = movement.to_numpy()
mov_full = np.hstack([mov_basic, vid["video"], vid["motion"]])
X, index, analog = wf.build_design(trials, rows_trial, rows_tp, points, ak["prior"], mov_full, rows_frame)
Y = svt[rows_frame]
task_cols = np.concatenate([index[k] for k in ("choice", "stim_left", "stim_right", "outcome", "true_prior", "subjective_prior")])
mov_cols = index["movement_columns"]
print(f"Valid trials {valid.sum()} of {len(trials)}; rows {X.shape[0]} ({points} frames per trial); columns {X.shape[1]}: "
f"intercept {len(index['intercept'])}, task {len(task_cols)}, movement {len(mov_cols)} "
f"({mov_basic.shape[1]} basic + {vid['video'].shape[1]} video + {vid['motion'].shape[1]} motion energy)")
print(f"Rows in which a frame appears in both windows of a trial: {len(rows_frame) - len(np.unique(np.stack([rows_frame, rows_trial]), axis=1)[0])}")
u_flat = svd["u"][:, ::2, ::2].reshape(len(svd["u"]), -1).T.astype(np.float64)
print(f"Per-pixel cvR2 over every second pixel in each direction: {u_flat.shape[0]} pixels")
Body camera: 119964 decoded frames, 119964 timestamps, 1280 pixels after binning, in 83 s; top 200 components explain 0.996 of video variance and 0.969 of motion-energy variance
Action-kernel learning rate 0.40 (negative log-likelihood 321.0; at rate 0.01 402.0, at 0.99 343.4); weights (bias, contrast, sqrt contrast, prior) [0.54, 0.84, 2.8, 3.59] Correlation of subjective prior with block (0.5 - probability_left): 0.604 Valid trials 911 of 913; rows 63770 (70 frames per trial); columns 898: intercept 70, task 420, movement 408 (8 basic + 200 video + 200 motion energy)
Rows in which a frame appears in both windows of a trial: 8971 Per-pixel cvR2 over every second pixel in each direction: 86400 pixels
In [22]:
t0 = time.time()
rng_fit = np.random.default_rng(404)
test_folds = wf.folds(rows_tp, rng_fit)
all_cols = np.concatenate([task_cols, mov_cols])
models = {
"full": [], "task shuffled": task_cols, "movement shuffled": mov_cols, "all shuffled": all_cols,
}
fits = {name: wf.fit_model(X, Y, analog, test_folds, cols, np.random.default_rng(404 + i), u_flat) for i, (name, cols) in enumerate(models.items())}
# The intercept-only model: the time-point columns alone, task and movement columns dropped.
icp = index["intercept"]
fits["intercept only"] = wf.fit_model(X[:, icp], Y, [], test_folds, [], np.random.default_rng(403), u_flat)
# Variant: without the body video components (columns dropped).
keep = np.setdiff1d(np.arange(X.shape[1]), mov_cols[mov_basic.shape[1]:])
Xb = X[:, keep]
pos = {c: i for i, c in enumerate(keep)}
mov_b = np.array([pos[c] for c in mov_cols[:mov_basic.shape[1]]]); task_b = np.array([pos[c] for c in task_cols])
an_b = [pos[c] for c in analog if c in pos]
fits["no video: full"] = wf.fit_model(Xb, Y, an_b, test_folds, [], np.random.default_rng(420), u_flat)
fits["no video: movement shuffled"] = wf.fit_model(Xb, Y, an_b, test_folds, mov_b, np.random.default_rng(421), u_flat)
fits["no video: task shuffled"] = wf.fit_model(Xb, Y, an_b, test_folds, task_b, np.random.default_rng(422), u_flat)
# Variant: each frame once per trial (the frames in both windows of a trial kept in the stimulus window only), and folds
# that keep each trial's rows together, so no frame of a test trial is in the training rows.
_, first = np.unique(np.stack([rows_trial, rows_frame], 1), axis=0, return_index=True)
keep_rows = np.sort(first)
trial_test = wf.trial_folds(rows_trial[keep_rows], np.random.default_rng(430))
Xg, Yg = X[keep_rows], Y[keep_rows]
for i, (name, cols) in enumerate((("full", []), ("task shuffled", task_cols), ("movement shuffled", mov_cols))):
fits[f"trial folds: {name}"] = wf.fit_model(Xg, Yg, analog, trial_test, cols, np.random.default_rng(431 + i), u_flat)
fits["trial folds: intercept only"] = wf.fit_model(Xg[:, icp], Yg, [], trial_test, [], np.random.default_rng(434), u_flat)
print(f"Trial-fold variant: {len(keep_rows)} of {len(rows_frame)} rows, {len(np.unique(rows_trial))} trials in {wf.N_FOLDS} folds")
readings = {"a pooled pixel space": "pooled", "b mean over pixels": "pixel_mean", "c mean over components": "component_mean"}
table = pd.DataFrame({name: {r: fit[k] for r, k in readings.items()} for name, fit in fits.items()}).T
print(f"cvR2 by model and reading, mean over {wf.N_FOLDS} folds ({time.time() - t0:.0f} s):")
print(table.round(4).to_string())
print("Fold range of the full model:", {r: (round(fits['full']['folds'][k].min(), 4), round(fits['full']['folds'][k].max(), 4)) for r, k in readings.items()})
print(f"Ridge penalties chosen per component and fold (grid {wf.ALPHAS[0]:g} to {wf.ALPHAS[-1]:g}):")
print(pd.DataFrame({k: {"min": v["alphas"].min(), "median": np.median(v["alphas"]), "max": v["alphas"].max(),
"at grid minimum": (v["alphas"] == wf.ALPHAS[0]).mean(), "at grid maximum": (v["alphas"] == wf.ALPHAS[-1]).mean()}
for k, v in fits.items()}).T.round(3).to_string())
# Variant: reading (b) over the pixels of the mask of the fig3b variant (brighter than the median of the mean image and 20
# pixels or more from the edge), at the resolution of the per-pixel maps.
mask_half = mask.reshape(H, W)[::2, ::2].ravel()
print(f"Reading (b) over the {int(mask_half.sum())} masked pixels of {mask_half.size}:",
{k: round(float(v["pixel_map"][mask_half].mean()), 4) for k, v in fits.items()})
Trial-fold variant: 54799 of 63770 rows, 911 trials in 10 folds
cvR2 by model and reading, mean over 10 folds (139 s):
a pooled pixel space b mean over pixels c mean over components
full 0.4709 0.4218 0.1264
task shuffled 0.4540 0.4067 0.1131
movement shuffled 0.3671 0.3028 0.0590
all shuffled 0.3405 0.2796 0.0394
intercept only 0.3481 0.2886 0.0486
no video: full 0.4140 0.3591 0.0757
no video: movement shuffled 0.3715 0.3078 0.0656
no video: task shuffled 0.3934 0.3417 0.0604
trial folds: full 0.4224 0.3629 0.0472
trial folds: task shuffled 0.4051 0.3472 0.0351
trial folds: movement shuffled 0.3739 0.3032 0.0524
trial folds: intercept only 0.3541 0.2888 0.0447
Fold range of the full model: {'a pooled pixel space': (np.float64(0.466), np.float64(0.4777)), 'b mean over pixels': (np.float64(0.4148), np.float64(0.4313)), 'c mean over components': (np.float64(0.1236), np.float64(0.1282))}
Ridge penalties chosen per component and fold (grid 0.1 to 1e+06):
min median max at grid minimum at grid maximum
full 3.162 56.234 5623.413 0.0 0.000
task shuffled 17.783 177.828 5623.413 0.0 0.000
movement shuffled 3.162 56.234 1000000.000 0.0 0.005
all shuffled 31.623 177.828 1000000.000 0.0 0.125
intercept only 1.000 31.623 562.341 0.0 0.000
no video: full 3.162 56.234 1778.279 0.0 0.000
no video: movement shuffled 3.162 56.234 316.228 0.0 0.000
no video: task shuffled 31.623 177.828 5623.413 0.0 0.000
trial folds: full 1.778 56.234 5623.413 0.0 0.000
trial folds: task shuffled 17.783 177.828 10000.000 0.0 0.000
trial folds: movement shuffled 1.778 56.234 1000000.000 0.0 0.013
trial folds: intercept only 1.778 31.623 1000.000 0.0 0.000
Reading (b) over the 43076 masked pixels of 86400: {'full': 0.4335, 'task shuffled': 0.418, 'movement shuffled': 0.3055, 'all shuffled': 0.2811, 'intercept only': 0.2899, 'no video: full': 0.3678, 'no video: movement shuffled': 0.3104, 'no video: task shuffled': 0.3496, 'trial folds: full': 0.3739, 'trial folds: task shuffled': 0.3574, 'trial folds: movement shuffled': 0.3067, 'trial folds: intercept only': 0.2906}
fig4c-full-cvr2¶
The full linear encoding model explains more than half of single-trial cortical variance (cvR2 0.60 +/- 0.02 in KO and 0.54 +/- 0.02 in WT).
In [23]:
lo, hi = 0.573 - 0.087, 0.573 + 0.087
rows = []
for r in readings:
full, icpt = table.loc["full", r], table.loc["intercept only", r]
rows.append({"reading": r, "full": full, "intercept only": icpt, "full minus intercept": full - icpt,
"full in range": lo <= full <= hi, "minus intercept in range": lo <= full - icpt <= hi})
cv = pd.DataFrame(rows).set_index("reading")
print(f"Tolerance range {lo:.3f} to {hi:.3f}")
print(cv.round(4).to_string())
print(f"Full minus the model with all task and movement columns shuffled: " + ", ".join(f"{r} {table.loc['full', r] - table.loc['all shuffled', r]:.4f}" for r in readings))
print("Variant without body video components:")
print(pd.DataFrame({"full": table.loc["no video: full"], "minus intercept": table.loc["no video: full"] - table.loc["intercept only"]}).round(4).to_string())
print("Variant, trial folds without repeated frames:")
print(pd.DataFrame({"full": table.loc["trial folds: full"], "minus intercept": table.loc["trial folds: full"] - table.loc["trial folds: intercept only"]}).round(4).to_string())
bm = {k: float(fits[k]["pixel_map"][mask_half].mean()) for k in ("full", "intercept only", "trial folds: full", "trial folds: intercept only")}
print(f"Variant, reading (b) over masked pixels: full {bm['full']:.4f}, minus intercept {bm['full'] - bm['intercept only']:.4f}; "
f"trial folds full {bm['trial folds: full']:.4f}, minus intercept {bm['trial folds: full'] - bm['trial folds: intercept only']:.4f}")
pm = fits["full"]["pixel_map"]
print(f"Per-pixel cvR2 of the full model: median {np.median(pm):.3f}, 5th to 95th percentile {np.percentile(pm, 5):.3f} to {np.percentile(pm, 95):.3f}")
Tolerance range 0.486 to 0.660
full intercept only full minus intercept full in range minus intercept in range
reading
a pooled pixel space 0.4709 0.3481 0.1227 False False
b mean over pixels 0.4218 0.2886 0.1332 False False
c mean over components 0.1264 0.0486 0.0778 False False
Full minus the model with all task and movement columns shuffled: a pooled pixel space 0.1304, b mean over pixels 0.1422, c mean over components 0.0870
Variant without body video components:
full minus intercept
a pooled pixel space 0.4140 0.0659
b mean over pixels 0.3591 0.0705
c mean over components 0.0757 0.0271
Variant, trial folds without repeated frames:
full minus intercept
a pooled pixel space 0.4224 0.0683
b mean over pixels 0.3629 0.0741
c mean over components 0.0472 0.0025
Variant, reading (b) over masked pixels: full 0.4335, minus intercept 0.1436; trial folds full 0.3739, minus intercept 0.0832
Per-pixel cvR2 of the full model: median 0.424, 5th to 95th percentile 0.237 to 0.587
In [24]:
fig, axs = plt.subplots(1, 2, figsize=(9, 3.2))
im = axs[0].imshow(pm.reshape(H // 2, W // 2), vmin=0, vmax=1, cmap="viridis"); plt.colorbar(im, ax=axs[0], fraction=0.04)
axs[0].set_title("full model cvR2 per pixel")
im = axs[1].imshow((pm - fits["intercept only"]["pixel_map"]).reshape(H // 2, W // 2), vmin=0, vmax=1, cmap="viridis")
plt.colorbar(im, ax=axs[1], fraction=0.04); axs[1].set_title("full minus intercept only")
plt.show()
d = pm - fits["intercept only"]["pixel_map"]
print(f"Pixels outside 0 to 1: full {int(((pm < 0) | (pm > 1)).sum())}, difference {int(((d < 0) | (d > 1)).sum())} of {pm.size}")
Pixels outside 0 to 1: full 0, difference 0 of 86400
Result: partially reproduced. The full model's cvR2 is 0.471 pooled in pixel space, 0.422
as the mean over pixels and 0.126 as the mean over components; with the intercept-only model
(the time-point columns alone) subtracted, 0.123, 0.133 and 0.078. None is within 0.486 to 0.660;
reading (a) is above 0.3, which the criterion counts as partial. No penalty is at either end of
the grid for the full model. Four variants do not change the outcome. Over the masked pixels of
the fig3b variant, reading (b) is 0.434. With folds that keep each trial's frames together and
each frame once per trial, the full model's cvR2 falls to 0.422, 0.363 and 0.047, and to 0.068,
0.074 and 0.003 above the intercept-only model, so much of its gain over the intercept under the
chiCa folds, which share frames of a trial between training and test rows, does not carry over to
new trials. With the video components of the left and right cameras added, as in the paper, the
notebook encoding_all_cameras.ipynb gives a full cvR2 of 0.574 (a) and 0.525 (b), within the range,
and 0.208 (c); with folds by trial it is 0.472, 0.405 and 0.046, below the range again. The claims
file fixed the body-camera design, so this does not decide the outcome. Without the body camera's components the
pooled cvR2 falls from 0.471 to 0.414.
The paper is ambiguous on the reading. Its Results say the model accounts "for more than half of the total observed single trial variance" and give 0.60 and 0.54 as the full model's cvR2, while the Figure 4C legend says the full model is shown "with trial time regressor subtracted out". Under the subtracted reading the values here are 0.123 and 0.133 with the body camera and 0.226 and 0.236 with all three cameras (readings a and b), far below the range. The criterion's reading, without the subtraction, decides the outcome.
fig4cd-movement-dominates¶
Movement predictors account for most of the explained variance; shuffling them reduces cvR2 by much more than shuffling task predictors (delta R2 for movement 0.40 +/- 0.02 in KO and 0.33 +/- 0.02 in WT; task shuffle beta -0.02).
In [25]:
lo, hi = 0.3685 - 0.0715, 0.3685 + 0.0715
rows = []
for r in readings:
dm = table.loc["full", r] - table.loc["movement shuffled", r]
dt_ = table.loc["full", r] - table.loc["task shuffled", r]
k = readings[r]
fm = fits["full"]["folds"][k] - fits["movement shuffled"]["folds"][k]
ft = fits["full"]["folds"][k] - fits["task shuffled"]["folds"][k]
rows.append({"reading": r, "movement dR2": dm, "task dR2": dt_, "movement in range": lo <= dm <= hi,
"task below movement/5": dt_ < dm / 5, "movement folds min": fm.min(), "task folds max": ft.max()})
dd = pd.DataFrame(rows).set_index("reading")
print(f"Tolerance range for movement dR2: {lo:.4f} to {hi:.4f}")
print(dd.round(4).to_string())
print("Variant without body video components:")
for r in readings:
print(f" {r}: movement dR2 {table.loc['no video: full', r] - table.loc['no video: movement shuffled', r]:.4f}, "
f"task dR2 {table.loc['no video: full', r] - table.loc['no video: task shuffled', r]:.4f}")
print("Variant, trial folds without repeated frames:")
for r in readings:
print(f" {r}: movement dR2 {table.loc['trial folds: full', r] - table.loc['trial folds: movement shuffled', r]:.4f}, "
f"task dR2 {table.loc['trial folds: full', r] - table.loc['trial folds: task shuffled', r]:.4f}")
bmm = {k: float(fits[k]["pixel_map"][mask_half].mean()) for k in ("full", "movement shuffled", "task shuffled")}
print(f"Variant, reading (b) over masked pixels: movement dR2 {bmm['full'] - bmm['movement shuffled']:.4f}, task dR2 {bmm['full'] - bmm['task shuffled']:.4f}")
print("Explained variance above intercept only, and the share left after each shuffle:")
for r in readings:
above = table.loc["full", r] - table.loc["intercept only", r]
print(f" {r}: full minus intercept {above:.4f}; movement-shuffled minus intercept {table.loc['movement shuffled', r] - table.loc['intercept only', r]:.4f}; "
f"task-shuffled minus intercept {table.loc['task shuffled', r] - table.loc['intercept only', r]:.4f}")
Tolerance range for movement dR2: 0.2970 to 0.4400
movement dR2 task dR2 movement in range task below movement/5 movement folds min task folds max
reading
a pooled pixel space 0.1038 0.0168 False True 0.0958 0.0193
b mean over pixels 0.1190 0.0151 False True 0.1103 0.0171
c mean over components 0.0673 0.0133 False True 0.0651 0.0142
Variant without body video components:
a pooled pixel space: movement dR2 0.0424, task dR2 0.0206
b mean over pixels: movement dR2 0.0513, task dR2 0.0174
c mean over components: movement dR2 0.0101, task dR2 0.0153
Variant, trial folds without repeated frames:
a pooled pixel space: movement dR2 0.0485, task dR2 0.0173
b mean over pixels: movement dR2 0.0597, task dR2 0.0156
c mean over components: movement dR2 -0.0052, task dR2 0.0121
Variant, reading (b) over masked pixels: movement dR2 0.1280, task dR2 0.0155
Explained variance above intercept only, and the share left after each shuffle:
a pooled pixel space: full minus intercept 0.1227; movement-shuffled minus intercept 0.0189; task-shuffled minus intercept 0.1059
b mean over pixels: full minus intercept 0.1332; movement-shuffled minus intercept 0.0142; task-shuffled minus intercept 0.1181
c mean over components: full minus intercept 0.0778; movement-shuffled minus intercept 0.0105; task-shuffled minus intercept 0.0645
Result: partially reproduced. Shuffling the movement predictors lowers cvR2 by 0.104
(pooled), 0.119 (mean over pixels) and 0.067 (mean over components); shuffling the task predictors
lowers it by 0.017, 0.015 and 0.013. Movement dominates under every reading and in every fold
(smallest movement drop 0.065, largest task drop 0.019), but the movement drop is below the range
0.297 to 0.440 under every reading. The body camera's video components carry most of it: without
them the movement drop is 0.042 pooled against 0.021 for task, and 0.010 against 0.015 over
components. With folds by trial and each frame once per trial, the movement drop is 0.049, 0.060
and -0.005 against 0.017, 0.016 and 0.012 for task, so over components task exceeds movement,
which under the criterion would count as not reproduced; the criterion's design uses the chiCa
folds, so the outcome stays. With all three cameras (encoding_all_cameras.ipynb) the movement
drop is 0.216, 0.232 and 0.161, still below the range, against 0.008, 0.008 and 0.007 for task;
with folds by trial, 0.109, 0.115 and 0.008 against 0.008, 0.008 and 0.000.
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Data | Outcome |
|---|---|---|
| fig1c-psychometric | dandiset | reproduced |
| fig1d-psychometric-params-genotype | dandiset | blocked: data missing |
| fig1e-rt-contrast | dandiset | not reproduced |
| fig1f-rt-genotype | dandiset | blocked: data missing |
| fig1g-learning-genotype | dandiset | blocked: data missing |
| fig2bc-prior-shift | dandiset | reproduced |
| fig2bd-ko-smaller-shift | dandiset | blocked: data missing |
| fig2ef-block-adaptation | dandiset | reproduced |
| fig2g-ko-slower-adaptation | dandiset | blocked: data missing |
| fig3b-contralateral-visual | dandiset | not reproduced |
| fig3c-ko-elevated-activity | dandiset | blocked: data missing |
| fig3d-ko-higher-correlation | dandiset | blocked: data missing |
| fig3e-ko-higher-variance | dandiset | blocked: data missing |
| fig4c-full-cvr2 | dandiset | partially reproduced |
| fig4cd-movement-dominates | dandiset | partially reproduced |
| fig4cd-ko-more-movement | dandiset | blocked: data missing |
| figs2-ko-elevated-motion | dandiset | blocked: data missing |
| figs3-ko-higher-unique | dandiset | blocked: data missing |
How the analysis differs from the paper's. One session of one mouse, where the paper uses 10
sessions from each of 19 mice; the genotype is unknown, so values reported per genotype are
compared with a range spanning both groups. The encoding model has no video components from the left
and right cameras (left out for runtime; the notebook encoding_all_cameras.ipynb adds them as a
variant), uses Lightning Pose paw speeds where the paper uses
DeepLabCut, takes the ROI motion energies to be the whisker pads and body, and uses an action-kernel
prior fitted here. The trial windows, the ridge penalty selection (generalized cross-validation) and
the three readings of cvR2 are choices made here, as the paper's widefield code is not public; the
lab's chiCa code sets the rest, including folds that split the frames of a trial between training
and test rows. Reaction time is the trial table's wheel movement onset, which falls within 0.08 s
of the stimulus or before it on more than half of the trials. VISp cannot be located without an atlas registration, so the visual
response is tested in image coordinates.
What was found about the dandiset. It holds one session of one of the paper's 19 mice, and no genotype for its mouse. It has no published version. The atlas landmarks and registration that the conversion documentation describes are not in either file. The imaging planes give GCaMP6f where the paper and the file's keywords give GCaMP6s, the raw file states an imaging rate of 15 Hz where the timestamps give 31.2 Hz per channel, the dandiset's keywords name "CNTNALP2" where the paper and files name Nrxn1a, and the wheel's stated sign convention (counter-clockwise positive) is the opposite of what the choice labels imply. The trial table, wheel, licks and imaging are consistent with a shared clock.
In [26]:
print(f"Notebook ran in {(time.time() - started) / 60:.1f} min")
Notebook ran in 5.1 min