Sit & Goard (2023) on dandiset 000483: single-cue rotation (Figure 1)¶
Paper. Sit KK, Goard MJ (2023). Coregistration of heading to visual cues in retrosplenial
cortex. Nature Communications 14:1992.
doi:10.1038/s41467-023-37704-5. The claims file
(claims/000483/10.1038-s41467-023-37704-5.yaml) was written from the version of record (main
text, legends, Methods; the Supplementary Information was not available).
Dataset. DANDI:000483, version 0.230421.2321. The paper cites no version of the dandiset (it cites a Dryad deposit and the dandiset by its description), so the most recent published version is pinned, and the only one: the draft has the same 128 assets and size and was last modified two months later, on metadata.
This notebook covers the single-cue controlled-rotation experiments (Figure 1). The head-rotation
and remapping claims are in rotation_remapping.ipynb, the dual-cue claims in
axons.ipynb and somata.ipynb, the decoding of Figure 6 in decoding.ipynb, and Figures S5 and S6 and the visual-cue-only files of Figure 5K in directions.ipynb.
| Claim | Here |
|---|---|
| fig1-free-heading-fraction | blocked: the free-chamber recordings are not in the dataset |
| fig1g-rotation-heading-fraction | tested |
| selection-shuffled-control | tested |
| fig1h-preferred-headings-tile | tested |
The dataset holds raw ROI fluorescence only (no neuropil trace, dF/F or inferred spikes), so dF/F
and inferred spikes are computed here, without the neuropil correction of the Methods (eq. 1).
Where each parameter comes from is in the docstring of rsc.py.
import os
for _v in ["OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"]:
os.environ[_v] = "1" # one BLAS thread, so that results do not depend on the number of cores
import time
import warnings
from multiprocessing import Pool
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import rsc
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
pd.set_option("display.max_rows", 200)
started = time.time()
cat = rsc.catalog(full=True)
print(f"dandiset {rsc.DANDISET_ID} version {rsc.DANDISET_VERSION}: {len(cat)} NWB files, "
f"read in {time.time() - started:.0f} s")
dandiset 000483 version 0.230421.2321: 128 NWB files, read in 109 s
Data check¶
The files do not say which experiment each belongs to except in their identifier; the helper
names them from it (rsc.classify). Counts of files, mice and ROIs per experiment, against the
paper's numbers:
paper_n = {"single_rotation": "16 recordings over 5 mice, 3179 neurons (Fig. 1G)",
"head_rotation": "5 recordings over 4 mice (Fig. 2)",
"chamber_rotation": "paired with the head-rotation recordings",
"remapping": "12 recordings (pairs of contexts) over 4 mice (Fig. 3)",
"adn_axons": "17 recordings over 4 mice (Fig. 4)",
"vis_axons": "7 recordings over 3 mice (Fig. 4)",
"soma_dual": "36 recordings over 6 mice (Figs. 5, 6)",
"motor_both": "Fig. S5 (rotation in both directions): 4 recordings over 2 mice",
"animal_controlled": "Fig. S6 (animal against motor control): 7 recordings over 3 mice",
"motor_cw": "Fig. S5/S6 (motor rotation, one direction)", "motor_ccw": "Fig. S5/S6",
"visual_only": "Fig. 5J, K (visual cue only rotation): number not given"}
summary = cat.groupby("experiment").agg(files=("path", "size"), mice=("subject", "nunique"),
rois=("n_rois", "sum"), frames=("n_frames", lambda x: sorted(set(x))))
summary["paper"] = summary.index.map(paper_n)
print(summary.to_string())
print("\nSubjects:", cat.subject.nunique(), "(archive summary: 20)")
print("Genotypes:", sorted(cat.genotype.unique()))
print("Sex by subject:", cat.drop_duplicates("subject").sex.value_counts().to_dict())
print("\nFluorescence rate (Hz):", sorted(cat.fluor_rate.unique()), "| heading 'rate' (Hz):",
sorted(cat.heading_rate.unique()), "| heading unit:", sorted(cat.heading_unit.astype(str).unique()))
print("Heading reference frame:", sorted(cat.heading_ref.unique()))
print("Frames of heading equal frames of fluorescence in every file:", bool((cat.n_heading == cat.n_frames).all()))
print("Imaging description:", sorted(cat.imaging_desc.unique()))
odd = cat[cat.apply(lambda r: r.subject not in r.identifier.split("_")[-1], axis=1)]
print("\nFiles whose identifier names another mouse than their subject:")
print(odd[["path", "identifier"]].to_string(index=False))
files mice rois frames paper
experiment
adn_axons 17 4 1084 [7550] 17 recordings over 4 mice (Fig. 4)
animal_controlled 5 4 1840 [7550] Fig. S6 (animal against motor control): 7 recordings over 3 mice
chamber_rotation 7 5 1436 [2550, 2600, 5100] paired with the head-rotation recordings
head_rotation 7 5 1436 [2600, 5100] 5 recordings over 4 mice (Fig. 2)
motor_both 4 3 1498 [7550, 15050] Fig. S5 (rotation in both directions): 4 recordings over 2 mice
motor_ccw 1 1 333 [7550] Fig. S5/S6
motor_cw 3 3 1084 [7550] Fig. S5/S6 (motor rotation, one direction)
remapping 24 4 5898 [7550] 12 recordings (pairs of contexts) over 4 mice (Fig. 3)
single_rotation 16 5 3179 [1950] 16 recordings over 5 mice, 3179 neurons (Fig. 1G)
soma_dual 35 6 11569 [1950, 7550] 36 recordings over 6 mice (Figs. 5, 6)
vis_axons 7 3 247 [7550] 7 recordings over 3 mice (Fig. 4)
visual_only 2 2 1137 [15050] Fig. 5J, K (visual cue only rotation): number not given
Subjects: 20 (archive summary: 20)
Genotypes: ['Slc17a7-Cre, AAV-Syn-GCaMP6s in ADN', 'Slc17a7-Cre, AAV-flex-GCaMP7b in AM/PM', 'Slc17a7-Cre, TITL-GCaMP6s', 'TITL-GCaMP6s, AAV-Syn-GCaMP6s in ADN', 'TITL-GCaMP6s, EMX1-Cre', 'TITL-GCaMP6s, Slc17a7-Cre', 'WT, AAV-Syn-GCaMP6s in AM/PM']
Sex by subject: {'M': 10, 'F': 10}
Fluorescence rate (Hz): [np.float64(10.0)] | heading 'rate' (Hz): [np.float64(0.1)] | heading unit: ['degrees']
Heading reference frame: ['cage door dead north']
Frames of heading equal frames of fluorescence in every file: True
Imaging description: ['soma recording from retrosplenial cortex']
Files whose identifier names another mouse than their subject:
path identifier
sub-SKKS089/sub-SKKS089_ses-20210429T110407_behavior+ophys.nwb 210429_KKS_SKKS089_HeadRotation_SKKS088-PlatformRotation-002
sub-SKKS089/sub-SKKS089_ses-20210429T105838_behavior+ophys.nwb 210429_KKS_SKKS089_HeadRotation_SKKS088-HeadRotation-001
sub-SKKS089/sub-SKKS089_ses-20210518T162426_behavior+ophys.nwb 210518_KKS_SKKS089_Remapping_SKKS088-Remapping-1-001
sub-SKKS089/sub-SKKS089_ses-20210518T163913_behavior+ophys.nwb 210518_KKS_SKKS089_Remapping_SKKS088-Remapping-2-002
Checks of the description against the values:
- The 16 single-cue files hold 3179 ROIs over 5 mice, exactly the 3179 neurons and 5 mice of the controlled-rotation experiment (Fig. 1G). The free-chamber experiment of Fig. 1D, E (4 recordings over 3 mice, 2665 neurons) has no files: no single-cue file is free (see below).
- The heading series states a rate of 0.1 Hz; it has one sample per imaging frame, and the fluorescence is at 10 Hz (Methods: 10 Hz), so the heading rate in the files is wrong (a metadata error). The two-photon series states about 20 Hz, but has as many frames as the fluorescence.
- Every file is described as a "soma recording", including the 24 axon recordings.
- Four files of mouse SKKS089 carry identifiers naming SKKS088.
- The files hold 20 mice, 10 male and 10 female; the paper used 26 (12 male, 14 female), so the six free-chamber and other mice not in the dataset are not described.
- There are 35 dual-cue soma recordings over 6 mice (paper: 36), 7 head-rotation pairs over 5 mice (paper: 5 over 4) and 12 remapping pairs over 4 mice (paper: 12 over 4).
single = cat[cat.experiment == "single_rotation"].reset_index(drop=True)
t0 = time.time()
recs = rsc.load_many(single.url.tolist())
with Pool(4) as pool:
for r in recs:
r["dff"] = rsc.dff(r["F"])
r["sp"] = rsc.spikes(r["dff"], pool)
del r["F"]
print(f"loaded and deconvolved {len(recs)} recordings in {time.time() - t0:.0f} s")
loaded and deconvolved 16 recordings in 24 s
rows = []
for meta, r in zip(single.itertuples(), recs):
h = rsc.clean_heading(r["heading"])
step = rsc.wrap180(np.diff(h))
segs = rsc.trial_segments(r["heading"])
rows.append(dict(subject=meta.subject, file=meta.identifier[-22:], frames=meta.n_frames, rois=meta.n_rois,
heading_nan=int(np.isnan(r["heading"]).sum()),
heading_min=round(np.nanmin(r["heading"]), 1), heading_max=round(np.nanmax(r["heading"]), 1),
median_step_deg=round(np.median(step), 2), frac_still=round(np.mean(np.abs(step) < 0.5), 3),
rpm=round(np.median(step) * 10 * 60 / 360, 2), trials=len(segs),
xy_range_mm=round(float(np.nanmax(np.abs(r["xy"]))), 1),
dff_median=round(float(np.median(r["dff"])), 3), dff_p99=round(float(np.percentile(r["dff"], 99)), 2),
spikes_per_s=round(float(r["sp"].sum(0).mean() / (meta.n_frames / 10)), 3)))
info = pd.DataFrame(rows)
print(info.to_string(index=False))
print("\nMice:", info.subject.nunique(), "| recordings:", len(info), "| ROIs:", info.rois.sum())
subject file frames rois heading_nan heading_min heading_max median_step_deg frac_still rpm trials xy_range_mm dff_median dff_p99 spikes_per_s SKKS080 _SKKS080-SingleCue-001 1950 293 0 -179.7 179.9 2.85 0.089 4.75 12 21.9 0.029 1.22 0.074 SKKS080 _SKKS080-SingleCue-002 1950 285 0 -179.9 179.9 2.86 0.084 4.76 14 21.9 0.036 1.34 0.089 SKKS080 _SKKS080-SingleCue-003 1950 313 0 -179.9 180.0 2.91 0.085 4.85 13 21.8 0.034 1.31 0.082 SKKS080 _SKKS080-SingleCue-004 1950 398 0 -179.8 179.9 3.04 0.087 5.06 13 21.9 0.035 1.31 0.083 SKKS080 _SKKS080-SingleCue-001 1950 98 0 -179.8 179.4 1.88 0.113 3.13 10 17.3 0.030 1.01 0.063 SKKS080 _SKKS080-SingleCue-002 1950 92 0 -179.1 179.6 2.32 0.081 3.87 12 17.3 0.030 0.97 0.052 SKKS081 _SKKS081-SingleCue-001 1950 264 0 -180.0 180.0 2.89 0.082 4.82 13 21.9 0.037 1.28 0.085 SKKS081 _SKKS081-SingleCue-002 1950 237 0 -179.8 179.9 3.01 0.082 5.02 15 21.9 0.034 1.19 0.079 SKKS081 _SKKS081-SingleCue-003 1950 257 0 -179.9 179.9 3.12 0.081 5.20 13 21.9 0.030 1.21 0.075 SKKS081 _SKKS081-SingleCue-004 1950 232 0 -180.0 179.8 3.18 0.094 5.31 14 21.8 0.039 1.24 0.085 SKKS089 _SKKS089-SingleCue-001 1950 116 0 -179.9 179.7 2.98 0.234 4.97 12 14.1 0.046 2.64 0.134 SKKS089 _SKKS089-SingleCue-002 1950 75 0 -180.0 179.9 3.68 0.136 6.14 16 14.2 0.049 2.37 0.109 SKKS090 _SKKS090-SingleCue-001 1950 203 0 -179.8 179.9 4.10 0.093 6.84 19 14.2 0.055 2.32 0.145 SKKS090 _SKKS090-SingleCue-002 1950 154 0 -179.8 179.9 4.38 0.088 7.30 21 14.2 0.055 2.64 0.159 SKKS091 _SKKS091-SingleCue-001 1950 78 0 -179.5 179.9 4.42 0.089 7.37 21 14.4 0.050 2.05 0.138 SKKS091 _SKKS091-SingleCue-002 1950 84 0 -180.0 179.8 4.79 0.094 7.98 23 14.2 0.047 1.94 0.123 Mice: 5 | recordings: 16 | ROIs: 3179
Every single-cue file rotates at a steady rate (median step per frame; 3.1 to 8.0 rpm), with few still frames: these are the motor-driven recordings (Methods: 3 to 7 rpm). None is a free-chamber recording, whose heading would wander. Heading is in degrees on -180 to 180, as stated.
Clock check. Heading and fluorescence have one sample per frame, but that does not show they
are aligned, and in a recording that rotates one way at constant speed a time offset only rotates
every tuning curve, which no measure within the recording can see. The files that rotate in
both directions (motor_both, Fig. S5) give a contrast: an offset of k frames between the
streams moves the tuning curves of clockwise and counterclockwise rotation in opposite
directions. The mean correlation, over ROIs, of each ROI's clockwise and counterclockwise
tuning curves is computed with the heading shifted by -10 to +10 frames against the activity;
with a shared clock it peaks at 0 or a little after (calcium and deconvolution delay).
both_dir = cat[cat.experiment == "motor_both"].reset_index(drop=True)
brecs = rsc.load_many(both_dir.url.tolist())
with Pool(4) as pool:
for r in brecs:
r["sp"] = rsc.spikes(rsc.dff(r["F"]), pool)
del r["F"]
lags = np.arange(-10, 11, 2)
rows = []
for m, r in zip(both_dir.itertuples(), brecs):
segs = rsc.trial_segments(r["heading"])
hc = rsc.clean_heading(r["heading"])
sign = np.array([np.sign(np.median(np.diff(hc[s]))) for s in segs])
cw = np.concatenate([s for s, g in zip(segs, sign) if g > 0])
ccw = np.concatenate([s for s, g in zip(segs, sign) if g < 0])
active = r["sp"].sum(0) > 0
row = {"file": m.identifier.split("_")[-1], "rois": int(active.sum()), "trials_+": int((sign > 0).sum()), "trials_-": int((sign < 0).sum())}
for lag in lags:
h = np.roll(hc, lag)
a, b = rsc.bin_tuning(r["sp"][cw][:, active], h[cw]), rsc.bin_tuning(r["sp"][ccw][:, active], h[ccw])
row[int(lag)] = round(float(np.nanmean(rsc._rowcorr(a, b))), 3)
rows.append(row)
clock = pd.DataFrame(rows)
print(clock.to_string(index=False))
print("Lag of the peak per file (frames; + = heading paired with later activity):",
[int(lags[np.argmax(clock[list(map(int, lags))].values[i])]) for i in range(len(clock))])
file rois trials_+ trials_- -10 -8 -6 -4 -2 0 2 4 6 8 10
SKKS125-DualRotation-002 399 41 42 0.000 -0.021 -0.021 0.004 0.036 0.045 0.058 0.073 0.078 0.069 0.046
SKKS126-DualRotation-002 461 39 42 0.043 0.055 0.066 0.056 0.017 -0.020 -0.037 -0.033 -0.015 0.013 0.031
SKKS127-DualDir-005 213 17 21 -0.025 -0.043 -0.057 -0.052 -0.052 -0.022 -0.014 -0.001 0.013 -0.001 -0.008
SKKS127-DualRotation-002 423 41 36 -0.042 -0.015 0.016 0.060 0.094 0.121 0.145 0.158 0.146 0.110 0.064
Lag of the peak per file (frames; + = heading paired with later activity): [6, -6, 6, 4]
The peak lag is printed per file above. A lag of a few frames after 0 is what the calcium and
deconvolution delay predict; a lag before 0 (activity leading heading) is not. Where the files
disagree in sign, the clock is consistent with the data in those that peak after 0 and is not
established in the others, and an offset that differs between two files cannot be excluded for
the other experiments either (see rotation_remapping.ipynb, fig2f).
fig1-free-heading-fraction (blocked)¶
In head-fixed mice moving a free floating chamber with a single cue, 14.0% of RSC neurons (373 of 2665) are heading selective, with no head movement.
Blocked, data_missing: the dataset holds no free-chamber recording. The 16 single-cue files are
motor-driven (table above) and hold the 3179 neurons of the controlled-rotation experiment; the
only animal-controlled files are dual-cue recordings of 2022 (Fig. S6).
fig1g-rotation-heading-fraction¶
With the chamber rotated by a motor (single cue), 14.2% of RSC neurons (452 of 3179, 16 recordings over 5 mice) are heading selective by the paper's reliability and Gaussian-fit test, a percentage similar to the free-chamber experiments.
The test of the claims file, on inferred spikes of all 3179 ROIs, with a single-term Gaussian.
The code fits the trial-mean curve as it is, with the peak location bounded to pi/2 +- pi/6
(bins 10 to 20): only a cell whose peak lies near bin 15 can be fitted well. The code also has a
center_flag that first shifts the curve's peak to bin 15, which its 'fit' preferred direction
uses. Both readings are run: (code) as the code calls it, and (centred).
trials = [rsc.trial_tuning(r["sp"], r["heading"]) for r in recs]
t0 = time.time()
with Pool(4) as pool:
sel_code = rsc.selection_many(trials, "single", 101, pool, centre=False)
sel_cent = rsc.selection_many(trials, "single", 102, pool, centre=True)
print(f"selection test, two readings, in {time.time() - t0:.0f} s")
selection test, two readings, in 143 s
def per_rec(sel):
return pd.DataFrame({"subject": single.subject, "file": single.identifier.str[-22:],
"trials": [t.shape[2] for t in trials], "rois": single.n_rois,
"reliable": [s.reliable.sum() for s in sel], "shape_ok": [s.shape_ok.sum() for s in sel],
"selective": [s.is_hd.sum() for s in sel]})
for name, sel in [("code (uncentred fit)", sel_code), ("centred fit", sel_cent)]:
t = per_rec(sel)
n, k = t.rois.sum(), t.selective.sum()
print(f"== {name}: {k} of {n} ROIs selective = {100 * k / n:.2f}% (paper 14.2% +- 1.42, 452 +- 45)")
by_mouse = t.groupby("subject")[["rois", "selective"]].sum()
by_mouse["percent"] = (100 * by_mouse.selective / by_mouse.rois).round(1)
print(by_mouse.T.to_string())
print("per recording percent: median", round(float(np.median(100 * t.selective / t.rois)), 1),
"range", round(float((100 * t.selective / t.rois).min()), 1), "to", round(float((100 * t.selective / t.rois).max()), 1))
print("\nPer recording, centred reading:")
print(per_rec(sel_cent).to_string(index=False))
print("\nROIs with fewer than 2 trials (test not run):", sum(t.shape[0] for t in trials if t.shape[2] < 2))
== code (uncentred fit): 345 of 3179 ROIs selective = 10.85% (paper 14.2% +- 1.42, 452 +- 45) subject SKKS080 SKKS081 SKKS089 SKKS090 SKKS091 rois 1479.0 990.0 191.0 357.0 162.0 selective 160.0 131.0 10.0 41.0 3.0 percent 10.8 13.2 5.2 11.5 1.9 per recording percent: median 10.2 range 1.0 to 25.3 == centred fit: 485 of 3179 ROIs selective = 15.26% (paper 14.2% +- 1.42, 452 +- 45) subject SKKS080 SKKS081 SKKS089 SKKS090 SKKS091 rois 1479.0 990.0 191.0 357.0 162.0 selective 209.0 153.0 17.0 90.0 16.0 percent 14.1 15.5 8.9 25.2 9.9 per recording percent: median 14.4 range 2.2 to 33.0 Per recording, centred reading: subject file trials rois reliable shape_ok selective SKKS080 _SKKS080-SingleCue-001 12 293 79 34 34 SKKS080 _SKKS080-SingleCue-002 14 285 132 60 60 SKKS080 _SKKS080-SingleCue-003 13 313 122 48 48 SKKS080 _SKKS080-SingleCue-004 13 398 131 58 58 SKKS080 _SKKS080-SingleCue-001 10 98 19 7 7 SKKS080 _SKKS080-SingleCue-002 12 92 16 2 2 SKKS081 _SKKS081-SingleCue-001 13 264 94 44 44 SKKS081 _SKKS081-SingleCue-002 15 237 99 36 36 SKKS081 _SKKS081-SingleCue-003 13 257 81 40 40 SKKS081 _SKKS081-SingleCue-004 14 232 85 33 33 SKKS089 _SKKS089-SingleCue-001 12 116 41 10 10 SKKS089 _SKKS089-SingleCue-002 16 75 27 7 7 SKKS090 _SKKS090-SingleCue-001 19 203 125 67 67 SKKS090 _SKKS090-SingleCue-002 21 154 68 23 23 SKKS091 _SKKS091-SingleCue-001 21 78 18 6 6 SKKS091 _SKKS091-SingleCue-002 23 84 22 10 10 ROIs with fewer than 2 trials (test not run): 0
How the count depends on the two thresholds of the reliability test (KS p and Cohen's d), centred reading, all else fixed:
allc = pd.concat(sel_cent, ignore_index=True)
for pthr, dthr in [(0.1, 0.8), (0.05, 0.8), (0.1, 0.5), (0.1, 1.0), (0.01, 1.2)]:
k = ((allc.ks_p < pthr) & (allc.cohen_d > dthr) & allc.shape_ok & (allc.a1 > 0)).sum()
print(f"KS p < {pthr}, d > {dthr}: {k} = {100 * k / len(allc):.2f}%")
KS p < 0.1, d > 0.8: 485 = 15.26% KS p < 0.05, d > 0.8: 485 = 15.26% KS p < 0.1, d > 0.5: 485 = 15.26% KS p < 0.1, d > 1.0: 444 = 13.97% KS p < 0.01, d > 1.2: 393 = 12.36%
False-positive rate of the test. The reliability step compares the split-half correlations
of 1000 random splits of the same trials with those of shifted trials by a two-sample KS test and
Cohen's d. The 1000 values are not independent (they reuse the same 13 or so trials), so the KS p
is far too small for a cell without tuning. This is the paper's own procedure: the code
(ControlledRotation.calculateHeadDirection, lines 40 to 97 at commit c6344f5) computes
bt_CC_data(iter) and bt_CC_rand(iter) over the iterations and then kstest2(bt_CC_data, bt_CC_rand) with P < 0.1 && d > 0.8, as here. The rate is measured below on synthetic untuned
cells (sparse random events, no heading dependence, 13 trials of 60 bins), and on the real cells
in selection-shuffled-control (each trial shifted once).
rng_fp = np.random.default_rng(106)
n_syn, n_tr_syn = 300, 13
frames_per_bin = 3
ev = rng_fp.random((n_syn, 60 * frames_per_bin, n_tr_syn)) < 0.02 # 2% of frames carry an event
amp = rng_fp.exponential(0.3, ev.shape) * ev
syn = amp.reshape(n_syn, 60, frames_per_bin, n_tr_syn).mean(2)
syn = (np.roll(syn, 1, 1) + syn + np.roll(syn, -1, 1)) / 3 # same 3-bin smoothing as bin_tuning
with Pool(4) as pool:
fp = {c: pd.concat(rsc.selection_many([syn], "single", 107, pool, centre=c), ignore_index=True) for c in [True, False]}
for c, df in fp.items():
print(f"synthetic untuned cells, {'centred' if c else 'code (uncentred)'} fit: reliable {100 * df.reliable.mean():.1f}%, "
f"pass the whole test {100 * df.is_hd.mean():.1f}% of {len(df)}")
synthetic untuned cells, centred fit: reliable 15.7%, pass the whole test 4.7% of 300 synthetic untuned cells, code (uncentred) fit: reliable 15.7%, pass the whole test 1.7% of 300
Result: reproduced, under the centred fit (an assumption). With the centred fit 485 of 3179 ROIs (15.26%) are heading selective, within 14.2 +- 1.42% and 452 +- 45; by mouse 8.9 to 25.2%. With the fit as the code calls it (uncentred), 345 (10.85%), outside the tolerance. The two readings give different outcomes; the centred one is adopted because only it lets cells of every preferred heading pass (see fig1h below: the uncentred fit puts 93.6% of selected cells in one quadrant). Stricter reliability thresholds lower the count (KS p < 0.01 and d > 1.2: 393, 12.36%). The test is liberal: synthetic untuned cells pass at 4.7% (centred) and 1.7% (uncentred), and 15.7% pass its reliability step; real cells with each trial shifted once pass at 5.19% (selection-shuffled-control). This is the paper's own test, so the count is comparable with the paper's, but about a third of the selected cells could pass by chance.
selection-shuffled-control¶
Run on pre-shuffled data, the heading-selection test passes no cells, so the count of heading-selective cells is not due to chance.
Two readings of "pre-shuffled". (code) As in the commented calculateHeadDirection_shuff of
Cue.m that the claims file cites: in every iteration the "real" trials are themselves circularly
shifted afresh (activity_rand2), and the real fit uses the shifted trials of the last
iteration. The "real" and the null split-half correlations then come from the same distribution,
so no cell can pass, whatever the data: this reading cannot show whether the count is due to
chance. (once) Each trial is shifted once by its own random amount before the whole test is run,
as the Methods' words ("repeating these procedures using pre-shuffled data") and the criterion's
("circularly shifting each trial's tuning curve ... before the test") say. Only this reading
measures how many untuned cells the procedure passes. Both fit readings are run under (code), the
adopted centred fit under (once).
t0 = time.time()
with Pool(4) as pool:
ctrl = {"each, centred": rsc.selection_many(trials, "single", 103, pool, centre=True, pre_shuffle="each"),
"each, code": rsc.selection_many(trials, "single", 104, pool, centre=False, pre_shuffle="each"),
"once, centred (variant)": rsc.selection_many(trials, "single", 105, pool, centre=True, pre_shuffle="once")}
print(f"control in {time.time() - t0:.0f} s")
control in 99 s
n_all = sum(len(s) for s in sel_cent)
for name, sel in ctrl.items():
k = sum(s.is_hd.sum() for s in sel)
rel_k = sum(s.reliable.sum() for s in sel)
print(f"{name}: {k} of {n_all} pass = {100 * k / n_all:.2f}% (reliable {rel_k}); "
f"per mouse {dict(zip(*np.unique(np.repeat(single.subject, [len(s) for s in sel])[np.concatenate([s.is_hd.values for s in sel])], return_counts=True)))}")
print("Thresholds of the claims file: none passing = reproduced; up to 1% = partially; above 1% = not reproduced")
each, centred: 0 of 3179 pass = 0.00% (reliable 0); per mouse {}
each, code: 0 of 3179 pass = 0.00% (reliable 0); per mouse {}
once, centred (variant): 165 of 3179 pass = 5.19% (reliable 478); per mouse {'SKKS080': np.int64(78), 'SKKS081': np.int64(60), 'SKKS089': np.int64(8), 'SKKS090': np.int64(13), 'SKKS091': np.int64(6)}
Thresholds of the claims file: none passing = reproduced; up to 1% = partially; above 1% = not reproduced
Result: not reproduced, under the (once) reading (an assumption). Shifting each trial once before the test, as the Methods describe, 165 of 3179 ROIs (5.19%) still pass, above the 1% limit of the criterion; 478 pass the reliability step. Synthetic untuned cells pass at a similar rate (4.7%, under fig1g). Run as the commented code does it, with the trials re-shifted in every iteration, 0 of 3179 pass under both fit readings, but under that reading no cell can pass whatever the data, as its "real" and null values come from the same distribution. Which reading the claims file intends is asked of a reader.
fig1h-preferred-headings-tile¶
The preferred headings of heading-selective cells tile the whole circumference of the chamber, including headings with the cue outside the visual field, with a moderate overrepresentation near 0 degrees, where the cue is centered in front of the mouse.
Preferred heading: the cross-validated peak (for each trial, the peak bin of the mean of the other trials; their circular mean). The dataset's heading is relative to "cage door dead north"; where the cue sits in that frame is not stated, so heading 0 of the dataset is not known to be the cue in front of the mouse. The share near dataset heading 0 is given, and where the densest 90 degree window lies.
The code's uncentred fit bounds the Gaussian's peak to pi/2 +- pi/6 on an axis whose 0 is the first bin (-180 degrees), that is, to dataset headings -120 to -60. That is why the cells it passes cluster near -90 (below): it is a property of the fit's bounds, not evidence of where the cue is. The share near -90 is printed as well, as the heading the code's bounds favour.
def tiling(sel, label):
pref = np.concatenate([rsc.cv_preferred(t)[s.is_hd.values] for t, s in zip(trials, sel)])
subj = np.concatenate([np.repeat(m, s.is_hd.sum()) for m, s in zip(single.subject, sel)])
bins = np.clip(((pref + 180) // 6).astype(int), 0, 59)
occ = np.bincount(bins, minlength=60)
quad = [np.mean(np.abs(rsc.wrap180(pref - c)) <= 45) for c in (0, 90, 180, -90)]
near0 = np.mean(np.abs(pref) <= 45)
cs = np.arange(-180, 180, 6)
shares = [np.mean(np.abs(rsc.wrap180(pref - c)) <= 45) for c in cs]
print(f"== {label}: {len(pref)} cells")
print(" (a) bins of 6 degrees holding at least one preferred heading:", int((occ > 0).sum()), "of 60; empty bins at",
(-180 + 6 * np.flatnonzero(occ == 0)).tolist())
print(" (b) share per 90-degree quadrant centred on 0, 90, 180, -90:", np.round(quad, 3).tolist())
print(f" share within 45 degrees of dataset heading 0: {near0:.3f} (uniform 0.25); of -90: "
f"{np.mean(np.abs(rsc.wrap180(pref + 90)) <= 45):.3f}; of +90: {np.mean(np.abs(rsc.wrap180(pref - 90)) <= 45):.3f}")
print(f" densest 90-degree window centred at {cs[int(np.argmax(shares))]} degrees with {max(shares):.3f}; "
f"sparsest at {cs[int(np.argmin(shares))]} with {min(shares):.3f}")
print(" Rayleigh test of uniformity: R =", round(float(np.abs(np.mean(np.exp(1j * np.deg2rad(pref))))), 3),
"p =", f"{np.exp(-len(pref) * np.abs(np.mean(np.exp(1j * np.deg2rad(pref)))) ** 2):.2g}",
"| mean direction", round(float(np.rad2deg(np.angle(np.mean(np.exp(1j * np.deg2rad(pref)))))), 1))
print(" counts in 30-degree bins from -180:", np.histogram(pref, bins=12, range=(-180, 180))[0].tolist())
for m in np.unique(subj):
pm = pref[subj == m]
print(f" {m}: n={len(pm)}, 30-degree bins:", np.histogram(pm, bins=12, range=(-180, 180))[0].tolist())
return pref
pref_cent = tiling(sel_cent, "centred reading")
pref_code = tiling(sel_code, "code reading")
== centred reading: 485 cells (a) bins of 6 degrees holding at least one preferred heading: 60 of 60; empty bins at [] (b) share per 90-degree quadrant centred on 0, 90, 180, -90: [0.247, 0.181, 0.184, 0.4] share within 45 degrees of dataset heading 0: 0.247 (uniform 0.25); of -90: 0.400; of +90: 0.181 densest 90-degree window centred at -102 degrees with 0.404; sparsest at 162 with 0.151 Rayleigh test of uniformity: R = 0.149 p = 2e-05 | mean direction -80.5 counts in 30-degree bins from -180: [31, 49, 76, 65, 34, 15, 63, 35, 35, 39, 19, 24] SKKS080: n=209, 30-degree bins: [23, 23, 34, 29, 17, 11, 31, 16, 5, 5, 4, 11] SKKS081: n=153, 30-degree bins: [7, 21, 32, 28, 15, 4, 21, 14, 0, 4, 5, 2] SKKS089: n=17, 30-degree bins: [1, 4, 0, 0, 1, 0, 9, 1, 0, 0, 0, 1] SKKS090: n=90, 30-degree bins: [0, 1, 10, 8, 1, 0, 2, 1, 23, 26, 9, 9] SKKS091: n=16, 30-degree bins: [0, 0, 0, 0, 0, 0, 0, 3, 7, 4, 1, 1] == code reading: 345 cells (a) bins of 6 degrees holding at least one preferred heading: 25 of 60; empty bins at [-180, -174, -30, -24, -18, -12, -6, 6, 12, 18, 24, 30, 36, 42, 48, 54, 60, 66, 72, 78, 84, 90, 96, 102, 108, 114, 120, 126, 132, 138, 144, 150, 156, 162, 168] (b) share per 90-degree quadrant centred on 0, 90, 180, -90: [0.026, 0.0, 0.043, 0.936] share within 45 degrees of dataset heading 0: 0.026 (uniform 0.25); of -90: 0.936; of +90: 0.000 densest 90-degree window centred at -96 degrees with 0.939; sparsest at 54 with 0.000 Rayleigh test of uniformity: R = 0.906 p = 1.1e-123 | mean direction -94.3 counts in 30-degree bins from -180: [4, 45, 148, 113, 33, 0, 1, 0, 0, 0, 0, 1] SKKS080: n=160, 30-degree bins: [4, 21, 70, 50, 15, 0, 0, 0, 0, 0, 0, 0] SKKS081: n=131, 30-degree bins: [0, 18, 55, 44, 14, 0, 0, 0, 0, 0, 0, 0] SKKS089: n=10, 30-degree bins: [0, 3, 2, 2, 2, 0, 1, 0, 0, 0, 0, 0] SKKS090: n=41, 30-degree bins: [0, 3, 18, 17, 2, 0, 0, 0, 0, 0, 0, 1] SKKS091: n=3, 30-degree bins: [0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0]
fig, ax = plt.subplots(1, 2, figsize=(9, 3))
for a, p, t in zip(ax, [pref_cent, pref_code], ["centred", "code"]):
a.hist(p, bins=60, range=(-180, 180), color="0.4")
a.set_xlabel("preferred heading (deg, dataset frame)"); a.set_ylabel("cells"); a.set_title(t)
plt.tight_layout(); plt.show()
print("All preferred headings lie in -180 to 180; none outside the axes.")
All preferred headings lie in -180 to 180; none outside the axes.
Result: partially reproduced, under the centred fit (an assumption). The 485 preferred headings tile the circle: all 60 bins of 6 degrees hold at least one, and the quadrants hold 18.1 to 40.0%. The overrepresentation cannot be placed relative to the cue: the files give heading relative to the cage door, not the cue. Within 45 degrees of the dataset's heading 0 lie 24.7% of the cells (uniform 25%). The densest 90-degree window, centred at -102 degrees, holds 40.4% (Rayleigh p = 2e-05), and it differs between mice (SKKS090 peaks at +60 to +120 and SKKS091 at +60 to +90, in 30-degree bins). Were the cue at dataset heading -90, the share near it would be 40.0%; nothing in the files says it is. With the code's uncentred fit, 25 of 60 bins hold a cell and 93.6% lie in one quadrant, near -90, where the fit's bounds put the peak.
Summary¶
| Claim | Outcome |
|---|---|
| fig1-free-heading-fraction | blocked: data missing (no free-chamber recordings) |
| fig1g-rotation-heading-fraction | reproduced, assumed (centred fit): 15.26% (485 of 3179); the test passes 4.7% of synthetic untuned cells |
| selection-shuffled-control | not reproduced, assumed (each trial shifted once): 5.19% pass; 0% under the code's reading, which cannot fail |
| fig1h-preferred-headings-tile | partially reproduced, assumed: tiling holds; the cue's position in the heading frame is not in the files |
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: dF/F and spikes are computed from raw fluorescence without neuropil correction (the files hold no neuropil trace), with the Python port of OASIS. The selection test and its Gaussian fits are reimplemented in numpy. The trial-mean curve is centred before fitting, an assumption; the code's own call does not centre it, and both readings are shown. The shuffled control shifts each trial once before the test, an assumption; the code's reading is shown too.
What was found about the dataset: no free-chamber recordings; the experiment of each file is named only in its identifier; the heading series states a rate of 0.1 Hz for 10 Hz data; every file, axon recordings included, is described as a soma recording; four SKKS089 files carry identifiers naming SKKS088; the position of the cue in the heading frame is not given; the files hold 20 mice where the paper used 26. The clock check peaks after 0 in three of the four files that rotate both ways (4 to 6 frames) and before 0 in one (-6 frames).