Sit & Goard (2023) on dandiset 000483: decoding heading from RSC somata (Figure 6)¶
Paper. Sit KK, Goard MJ (2023), Nature Communications 14:1992, doi:10.1038/s41467-023-37704-5; claims file written from the version of record without its supplement.
Dataset. DANDI:000483, version
0.230421.2321, the only published version (the paper cites none). See single_cue.ipynb for the
data check of all files.
| Claim | Here |
|---|---|
| fig6a-population-light-on | tested |
| fig6a-population-light-off | tested |
| fig6b-heading-cells | tested |
| fig6b-landmark-cells | tested |
| fig6b-alignment-cells | tested |
| fig6b-180-errors | tested |
| fig6-heading-below-population | tested |
The tuned somata and their classes are rebuilt here exactly as in somata.ipynb (same files,
seeds and code), so this notebook can run on its own; the class counts are printed so that they
can be compared with those of somata.ipynb. The classes found there differ from the paper's
(fig5-class-proportions, not reproduced), so the class decoders here decode other groups of cells
than the paper's.
Decoding follows the claims file's opening comment (Methods, "Decoding RSC somatic activity";
Statistical information), with its readings: the input, (a) dF/F or (b) inferred spikes; the
decoder, (a) a naive Bayes decoder on binarized activity (Etter et al. 2020), trained on even and
tested on odd time points, or (b) the lab's HeadingDecoder (github.com/ucsb-goard-lab/HeadingDecoder
at commit ca9842e: activity times row-rescaled leave-one-trial-out tuning curves, normalized,
filtered over 10 samples with its temporalFilter, maximum taken); the t-test, (a) paired over
trials or (b) unpaired; chance, (a) 0.18 (Methods) or (b) 0.10 (Figure 6A legend). All readings
are run.
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
from sklearn.mixture import GaussianMixture
import rsc
warnings.filterwarnings("ignore")
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
started = time.time()
cat = rsc.catalog(full=False)
sub = cat[cat.experiment == "soma_dual"].reset_index(drop=True)
recs = rsc.load_many(sub.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"dandiset {rsc.DANDISET_ID} version {rsc.DANDISET_VERSION}: {len(sub)} files loaded and deconvolved "
f"in {time.time() - started:.0f} s")
dandiset 000483 version 0.230421.2321: 35 files loaded and deconvolved in 205 s
rows = []
for m, r in zip(sub.itertuples(), recs):
L, D = rsc.light_dark_frames(m.n_frames)
r["seg_l"], r["seg_d"] = rsc.segments_within(r["heading"], L), rsc.segments_within(r["heading"], D)
r["seg_all"] = rsc.trial_segments(r["heading"])
r["blocks"] = (L, D)
per_l = [sum(1 for s in r["seg_l"] if s[0] >= b[0] and s[-1] <= b[-1]) for b in L]
per_d = [sum(1 for s in r["seg_d"] if s[0] >= b[0] and s[-1] <= b[-1]) for b in D]
rows.append(dict(subject=m.subject, identifier=m.identifier[:6] + " " + m.identifier.split("_")[-1], frames=m.n_frames,
rois=m.n_rois, light_trials_per_block=per_l, dark_trials_per_block=per_d,
median_trial_frames=int(np.median([len(s) for s in r["seg_l"]])) if r["seg_l"] else 0,
dff_range=f"{r['dff'].min():.2f}..{r['dff'].max():.1f}", spikes_max=round(float(r["sp"].max()), 2)))
info = pd.DataFrame(rows)
print(info.to_string(index=False))
print(f"{len(sub)} recordings, {sub.subject.nunique()} mice, {sub.n_rois.sum()} ROIs; dF/F is unitless, spikes are OASIS "
"event amplitudes in dF/F units")
subject identifier frames rois light_trials_per_block dark_trials_per_block median_trial_frames dff_range spikes_max SKKS080 201123 SKKS080-DualCue-003 1950 112 [5] [6] 156 -0.57..4.1 1.78 SKKS080 201123 SKKS080-DualCue-004 1950 110 [6] [6] 142 -0.49..3.4 1.58 SKKS088 210118 SKKS088-DualCue-001 1950 329 [6] [2] 132 -0.55..5.0 1.53 SKKS088 210118 SKKS088-DualCue-002 1950 240 [0] [3] 0 -0.68..5.5 1.63 SKKS088 210118 SKKS088-DualCue-003 1950 286 [1] [2] 145 -0.64..5.2 1.77 SKKS088 210216 SKKS088-DualCue-001 7550 411 [7, 6, 5] [6, 6, 7] 183 -0.55..4.9 1.76 SKKS088 210216 SKKS088-DualCue-002 7550 289 [6, 6, 7] [7, 5, 6] 182 -0.63..5.2 2.09 SKKS089 201205 SKKS089-DualCue-003 1950 115 [10] [10] 90 -0.59..6.6 2.68 SKKS089 201205 SKKS089-DualCue-004 1950 82 [9] [10] 85 -0.58..5.3 2.51 SKKS089 210118 SKKS089-DualCue-002 1950 307 [3] [2] 128 -0.58..5.4 4.87 SKKS089 210216 SKKS089-DualCue-001 7550 435 [5, 6, 6] [7, 6, 6] 182 -0.69..5.4 2.13 SKKS089 210216 SKKS089-DualCue-002 7550 459 [7, 5, 6] [6, 7, 6] 181 -0.67..5.8 2.32 SKKS089 210312 SKKS089-RSCSoma-3rpm-001 7550 498 [6, 5, 6] [6, 6, 6] 187 -0.78..7.8 3.87 SKKS089 210312 SKKS089-RSCSoma-3rpm-002 7550 430 [7, 6, 6] [5, 5, 6] 186 -0.77..6.8 3.21 SKKS089 210312 SKKS089-RSCSoma-3rpm-003 7550 518 [6, 5, 7] [6, 6, 5] 187 -0.86..7.1 3.62 SKKS090 201205 SKKS090-DualCue-003 1950 172 [11] [11] 80 -0.73..6.6 3.45 SKKS090 201205 SKKS090-DualCue-004 1950 145 [10] [11] 79 -0.63..6.3 2.52 SKKS090 210203 SKKS090-DualCue-001 7550 164 [6, 7, 6] [6, 6, 6] 180 -0.70..7.5 3.42 SKKS090 210216 SKKS090-DualCue-001 7550 368 [6, 5, 6] [6, 5, 6] 192 -0.56..5.8 2.38 SKKS090 210216 SKKS090-DualCue-002 7550 549 [5, 6, 5] [6, 6, 6] 191 -0.64..6.6 4.20 SKKS090 210312 SKKS090-RSCSoma-3rpm-001 7550 479 [5, 6, 7] [6, 6, 6] 185 -0.76..8.0 4.17 SKKS090 210312 SKKS090-RSCSoma-3rpm-002 7550 549 [7, 6, 6] [6, 6, 6] 184 -0.83..8.7 3.38 SKKS090 210312 SKKS090-RSCSoma-3rpm-003 7550 586 [6, 5, 5] [6, 6, 6] 184 -0.85..8.0 3.44 SKKS091 210203 SKKS091-DualCue-001 7550 162 [7, 6, 6] [6, 7, 6] 180 -0.58..5.8 1.93 SKKS091 210203 SKKS091-DualCue-002 7550 158 [7, 6, 6] [6, 6, 6] 180 -0.62..6.8 3.49 SKKS091 210216 SKKS091-DualCue-001 7550 326 [6, 7, 6] [7, 6, 6] 177 -0.66..6.9 2.98 SKKS091 210216 SKKS091-DualCue-002 7550 317 [7, 6, 6] [6, 6, 6] 177 -0.72..6.7 3.13 SKKS091 210311 SKKS091-RSCSoma-3rpm-001 7550 346 [6, 6, 6] [5, 5, 5] 185 -0.74..9.4 3.37 SKKS091 210311 SKKS091-RSCSoma-3rpm-002 7550 290 [6, 6, 6] [5, 5, 6] 185 -0.76..6.9 2.59 SKKS091 210311 SKKS091-RSCSoma-3rpm-003 7550 307 [6, 6, 6] [7, 5, 6] 183 -0.73..6.4 3.14 SKKS100 210216 SKKS100-DualCue-001 7550 496 [6, 6, 6] [6, 7, 6] 178 -0.81..6.6 2.66 SKKS100 210216 SKKS100-DualCue-002 7550 440 [6, 6, 7] [7, 5, 6] 178 -0.71..6.2 2.21 SKKS100 210311 SKKS100-RSCSoma-3rpm-001 7550 395 [6, 6, 6] [6, 5, 5] 183 -0.68..7.4 3.15 SKKS100 210311 SKKS100-RSCSoma-3rpm-002 7550 373 [6, 6, 5] [6, 7, 6] 186 -0.77..8.0 2.76 SKKS100 210311 SKKS100-RSCSoma-3rpm-003 7550 326 [5, 5, 6] [6, 6, 6] 185 -0.67..7.8 2.50 35 recordings, 6 mice, 11569 ROIs; dF/F is unitless, spikes are OASIS event amplitudes in dF/F units
Tuned somata and their classes (as in somata.ipynb)¶
Heading-selection test with the two-term Gaussian on all trials (centred fit, seed 306), light-on
and light-off curves aligned with the light-on shifts, Gaussian-fit features, three-component GMM
(random state 304) on the selective somata, archetypes and class names by the rule of
somata.ipynb (no outlier step).
t0 = time.time()
trials_all = [rsc.trial_tuning(r["sp"], r["heading"], r["seg_all"]) for r in recs]
with Pool(4) as pool:
sel = rsc.selection_many(trials_all, "double", 306, pool, centre=True)
def light_dark_aligned(i):
r = recs[i]
lt = rsc.trial_tuning(r["sp"], r["heading"], r["seg_l"])
dt = rsc.trial_tuning(r["sp"], r["heading"], r["seg_d"])
if lt.shape[2] == 0:
return None, None
sh = rsc.cv_shifts(lt) if lt.shape[2] > 1 else (rsc.SHIFT_POINT - lt[:, :, 0].argmax(1))[:, None]
la = rsc.apply_shifts(lt, sh).mean(2)
sh_all = (rsc.SHIFT_POINT - lt.mean(2).argmax(1))[:, None]
da = rsc.apply_shifts(dt, np.repeat(sh_all, max(dt.shape[2], 1), 1)).mean(2) if dt.shape[2] else np.full_like(la, np.nan)
return la, da
aligned = [light_dark_aligned(i) for i in range(len(recs))]
soma_idx = [i for i in range(len(recs)) if aligned[i][0] is not None]
L, D, cell_rec, cell_roi = [], [], [], []
for i in soma_idx:
m = sel[i].is_hd.values & ~np.isnan(aligned[i][1]).any(1)
L.append(aligned[i][0][m]); D.append(aligned[i][1][m])
cell_rec += [i] * int(m.sum()); cell_roi += np.flatnonzero(m).tolist()
cell_rec, cell_roi = np.array(cell_rec), np.array(cell_roi)
both = rsc.rescale(np.hstack([np.vstack(L), np.vstack(D)]))
pl, _ = rsc.fit_gauss(both[:, :60], "double", do_rescale=False)
pd_, _ = rsc.fit_gauss(both[:, 60:], "double", do_rescale=False)
X = np.column_stack([pl[:, [1, 2, 5, 6]], pd_[:, [1, 2, 5, 6]]])
g = GaussianMixture(3, covariance_type="full", n_init=20, max_iter=1000, random_state=304).fit(X)
def archetype(coef):
a1, a2, c1, c2 = coef
return a1 * np.exp(-((rsc.X - np.pi / 2) / c1) ** 2) + a2 * np.exp(-((rsc.X - 3 * np.pi / 2) / c2) ** 2)
def cosdist(A, b):
with np.errstate(invalid="ignore", divide="ignore"):
return 1 - (A @ b) / (np.linalg.norm(A, axis=1) * np.linalg.norm(b))
dist = np.column_stack([cosdist(both[:, :60], archetype(mu[:4])) + cosdist(both[:, 60:], archetype(mu[4:])) for mu in g.means_])
cls = np.nanargmin(np.nan_to_num(dist, nan=9), 1)
bimod = [min(mu[0], mu[1]) / max(mu[0], mu[1]) for mu in g.means_]
dark_amp = [max(mu[4], mu[5]) for mu in g.means_]
heading_c = int(np.argmin(bimod))
rest = [c for c in range(3) if c != heading_c]
landmark_c = rest[int(np.argmin([dark_amp[c] for c in rest]))]
names = {heading_c: "heading", landmark_c: "landmark", [c for c in rest if c != landmark_c][0]: "alignment"}
print(f"selection, fits and clustering in {time.time() - t0:.0f} s")
print(f"tuned somata with light-on and light-off trials: {len(cls)} from {len(soma_idx)} recordings")
print("class counts:", {names[c]: int((cls == c).sum()) for c in range(3)}, "(compare somata.ipynb, fig5-class-proportions)")
selection, fits and clustering in 451 s
tuned somata with light-on and light-off trials: 1419 from 34 recordings
class counts: {'heading': 202, 'landmark': 67, 'alignment': 1150} (compare somata.ipynb, fig5-class-proportions)
Pseudopopulation trials¶
Methods: the trials are the full rotations before and after each rest that separates light-on from light-off; each trial is resampled to 200 bins and trials are aligned across recordings into pseudopopulations. Here, in each recording, the last 5 light-on trials before and the first 5 light-off trials after each of the three light-to-dark transitions of a 7550-frame recording (15 trials each; HeadingDecoder takes up to 15 a side). The 1950-frame recordings have one transition and some have fewer than 5 trials a side, so they are left out (their cells are counted below). Each trial is resampled to 200 bins on its heading (bins every 1.8 degrees from -180 to 180) in place of time, so that bin k is the same heading in every recording; at the motor's constant speed the two differ only by speed fluctuations (a standard choice). Trial k of each recording is put with trial k of the others.
K = 5
N_RS = 200
grid = -180 + (np.arange(N_RS) + 0.5) * 360 / N_RS
def resample(act, h, seg):
hu = np.rad2deg(np.unwrap(np.deg2rad(h[seg])))
sgn = 1 if hu[-1] >= hu[0] else -1
hu = sgn * hu
target = hu[0] + np.mod(sgn * grid - hu[0], 360)
return np.stack([np.interp(target, hu, act[seg][:, n]) for n in range(act.shape[1])], 1)
def trials_around(i):
r = recs[i]
L, D = r["blocks"]
on, off = [], []
for bl, bd in zip(L, D):
sl = [s for s in r["seg_l"] if s[0] >= bl[0] and s[-1] <= bl[-1]]
sd = [s for s in r["seg_d"] if s[0] >= bd[0] and s[-1] <= bd[-1]]
if len(sl) < K or len(sd) < K:
return None
on += sl[-K:]; off += sd[:K]
return on, off
use = [i for i in soma_idx if sub.n_frames[i] == 7550 and trials_around(i) is not None]
inuse = np.isin(cell_rec, use)
print(f"recordings used: {len(use)} of {len(soma_idx)} (7550-frame with {K} trials on each side of all 3 transitions); "
f"tuned cells {inuse.sum()} of {len(cls)}; left out {int((~inuse).sum())}")
print("cells used by class:", {names[c]: int(((cls == c) & inuse).sum()) for c in range(3)},
"by mouse:", pd.Series(sub.subject[cell_rec[inuse]].values).value_counts().to_dict())
pop = {}
for inp in ["dff", "sp"]:
on_all, off_all = [], []
for i in use:
on, off = trials_around(i)
act = recs[i][inp][:, cell_roi[cell_rec == i]]
hc = rsc.clean_heading(recs[i]["heading"])
on_all.append(np.stack([resample(act, hc, s) for s in on]))
off_all.append(np.stack([resample(act, hc, s) for s in off]))
pop[inp] = {"on": np.concatenate(on_all, 2), "off": np.concatenate(off_all, 2)} # trials x 200 bins x cells
ucls = cls[inuse]
print("pseudopopulation shape (trials, bins, cells):", pop["sp"]["on"].shape)
recordings used: 25 of 34 (7550-frame with 5 trials on each side of all 3 transitions); tuned cells 1331 of 1419; left out 88
cells used by class: {'heading': 179, 'landmark': 50, 'alignment': 1102} by mouse: {'SKKS089': 528, 'SKKS100': 266, 'SKKS091': 242, 'SKKS090': 182, 'SKKS088': 113}
pseudopopulation shape (trials, bins, cells): (15, 200, 1331)
Decoders¶
- (a) Naive Bayes (Etter et al. 2020): activity binarized (spikes: above 0; dF/F: above its mean plus 2 SD over the recording's trials, a standard choice), P(active | heading state) per cell for 60 states of 6 degrees from the even bins of every trial, with a +1/+2 prior; the odd bins are decoded as the state of maximum summed log-likelihood.
- (b) HeadingDecoder: for each trial, each cell's tuning curve (60 bins) from the other trials of
the same condition, rescaled 0 to 1; the heading distribution is activity times tuning
weights, normalized, filtered over 10 samples (
temporalFilter); the maximum is taken. Activity below 0 (dF/F) is set to 0, as the product needs non-negative weights.
Accuracy per trial: the share of decoded bins within 18 degrees of the true heading. Each bootstrap iteration samples 332 cells with replacement from the group (all tuned cells: all three classes, so in their proportions); 1000 iterations; mean and 2.5 and 97.5 percentiles over iterations. p of light-on against light-off is one minus the share of iterations in which a t-test over trials gives p < 0.05 ((a) paired, (b) unpaired). Errors at 180 degrees: the share of light-on bins with an error of 162 degrees or more. Chance is computed with each decoder by decoding with each cell's test activity circularly shifted along the trial by its own random amount (the tuning to heading is broken, the activity kept).
NS = 60
state = np.clip(((grid + 180) // 6).astype(int), 0, NS - 1)
centres = -177 + 6 * np.arange(NS)
def err_deg(dec, true):
return np.abs((dec - true + 180) % 360 - 180)
def prep_bayes(A, inp):
"""Log-likelihood weights from the even bins; binarized odd bins to decode."""
n = A.shape[2]
flat = A.reshape(-1, n)
th = np.zeros(n) if inp == "sp" else flat.mean(0) + 2 * flat.std(0)
B = (A > th).astype(np.float32)
tr, te = B[:, 0::2], B[:, 1::2]
s_tr = state[0::2]
p = np.stack([(tr[:, s_tr == s].sum((0, 1)) + 1) / (tr[:, s_tr == s].shape[0] * tr[:, s_tr == s].shape[1] + 2)
for s in range(NS)], 1)
return dict(kind="bayes", W=(np.log(p) - np.log(1 - p)).astype(np.float32), B=np.log(1 - p).astype(np.float32),
X=te, true=grid[1::2])
def prep_hd(A, inp):
"""Leave-one-trial-out tuning weights per trial (trials x cells x 60), rescaled 0 to 1 per cell."""
A = np.clip(A, 0, None)
S = np.stack([A[:, state == s].mean(1) for s in range(NS)], 2)
loo = (S.sum(0)[None] - S) / (A.shape[0] - 1)
lo = loo.min(2, keepdims=True); rg = loo.max(2, keepdims=True) - lo
W = np.where(rg > 0, (loo - lo) / np.where(rg > 0, rg, 1), 0)
return dict(kind="hd", W=W.astype(np.float32), X=A.astype(np.float32), true=grid)
def temporal_filter(P, w=10):
"""HeadingDecoder.temporalFilter: product of (1 + p) over a window of w samples, renormalized."""
Lg = np.log1p(P)
c = np.concatenate([np.zeros_like(Lg[:, :1]), np.cumsum(Lg, 1)], 1)
F = P.shape[1]; i = np.arange(F)
a, b = np.clip(i - w // 2, 0, F), np.clip(i + w // 2, 0, F)
out = np.expm1(c[:, b] - c[:, a])
return out / out.sum(2, keepdims=True)
def decode(Pp, idx, X=None):
X = Pp["X"] if X is None else X
if Pp["kind"] == "bayes":
sc = X[:, :, idx] @ Pp["W"][idx] + Pp["B"][idx].sum(0)
else:
sc = np.matmul(X[:, :, idx], Pp["W"][:, idx])
tot = sc.sum(2, keepdims=True)
sc = temporal_filter(np.where(tot > 0, sc / np.where(tot > 0, tot, 1), 1 / NS))
return err_deg(centres[sc.argmax(2)], Pp["true"][None])
groups = {"all tuned": np.ones(len(ucls), bool), **{names[c]: ucls == c for c in range(3)}}
print("cells per group:", {k: int(v.sum()) for k, v in groups.items()})
cells per group: {'all tuned': 1331, 'heading': 179, 'landmark': 50, 'alignment': 1102}
t0 = time.time()
N_BOOT, N_SAMPLE = 1000, 332
results, per_iter = [], {}
chance = []
for inp in ["dff", "sp"]:
for dec_name, prep in [("(a) Bayes", prep_bayes), ("(b) HeadingDecoder", prep_hd)]:
P = {cond: prep(pop[inp][cond], inp) for cond in ["on", "off"]}
rng = np.random.default_rng(601)
# chance: test activity of each cell circularly shifted along the trial by its own random amount
ch = []
for _ in range(20):
idx = rng.choice(len(ucls), N_SAMPLE, replace=True)
Xs = P["on"]["X"]
sh = rng.integers(Xs.shape[1], size=(Xs.shape[0], Xs.shape[2]))
ar = (np.arange(Xs.shape[1])[None, :, None] - sh[:, None, :]) % Xs.shape[1]
ch.append(float((decode(P["on"], idx, np.take_along_axis(Xs, ar, 1)) <= 18).mean()))
chance.append(dict(input=inp, decoder=dec_name, chance_mean=round(np.mean(ch), 3), chance_max=round(np.max(ch), 3)))
for gname, m in groups.items():
pool_ = np.flatnonzero(m)
acc = {"on": np.zeros(N_BOOT), "off": np.zeros(N_BOOT)}
e180 = np.zeros(N_BOOT); err_on = np.zeros(N_BOOT)
pt_paired = np.zeros(N_BOOT, bool); pt_unpaired = np.zeros(N_BOOT, bool)
for it in range(N_BOOT):
idx = rng.choice(pool_, N_SAMPLE, replace=True)
a = {}
for cond in ["on", "off"]:
e = decode(P[cond], idx)
a[cond] = (e <= 18).mean(1)
acc[cond][it] = a[cond].mean()
if cond == "on":
e180[it] = (e >= 162).mean(); err_on[it] = e.mean()
pt_paired[it] = stats.ttest_rel(a["on"], a["off"]).pvalue < 0.05
pt_unpaired[it] = stats.ttest_ind(a["on"], a["off"]).pvalue < 0.05
per_iter[(inp, dec_name, gname)] = acc
results.append(dict(input=inp, decoder=dec_name, group=gname,
on_mean=acc["on"].mean(), on_lo=np.percentile(acc["on"], 2.5), on_hi=np.percentile(acc["on"], 97.5),
off_mean=acc["off"].mean(), off_lo=np.percentile(acc["off"], 2.5), off_hi=np.percentile(acc["off"], 97.5),
p_paired=1 - pt_paired.mean(), p_unpaired=1 - pt_unpaired.mean(),
off_lower=bool(acc["off"].mean() < acc["on"].mean()),
err180_on=e180.mean(), mean_err_on=err_on.mean()))
res = pd.DataFrame(results).round(3)
print(f"decoding in {time.time() - t0:.0f} s")
print(res.to_string(index=False))
print("\nChance accuracy with each decoder (20 draws of shifted test activity, light-on):")
print(pd.DataFrame(chance).to_string(index=False))
print("Chance readings of the claims file: (a) 0.18 (Methods), (b) 0.10 (Figure 6A legend; 36/360)")
decoding in 171 s input decoder group on_mean on_lo on_hi off_mean off_lo off_hi p_paired p_unpaired off_lower err180_on mean_err_on dff (a) Bayes all tuned 0.708 0.665 0.747 0.569 0.525 0.612 0.016 0.061 True 0.034 26.925 dff (a) Bayes heading 0.429 0.399 0.457 0.377 0.350 0.403 0.697 0.964 True 0.070 54.409 dff (a) Bayes landmark 0.396 0.377 0.414 0.207 0.192 0.220 0.000 0.000 True 0.120 60.465 dff (a) Bayes alignment 0.720 0.680 0.761 0.580 0.538 0.619 0.011 0.053 True 0.033 25.769 dff (b) HeadingDecoder all tuned 0.202 0.158 0.251 0.137 0.101 0.178 0.346 0.361 True 0.138 81.206 dff (b) HeadingDecoder heading 0.130 0.113 0.153 0.109 0.091 0.131 0.939 0.915 True 0.104 88.125 dff (b) HeadingDecoder landmark 0.201 0.171 0.234 0.103 0.086 0.120 0.003 0.001 True 0.115 79.021 dff (b) HeadingDecoder alignment 0.203 0.160 0.251 0.151 0.110 0.201 0.530 0.527 True 0.148 81.794 sp (a) Bayes all tuned 0.712 0.678 0.746 0.572 0.535 0.608 0.000 0.002 True 0.055 31.609 sp (a) Bayes heading 0.415 0.390 0.441 0.341 0.317 0.364 0.170 0.534 True 0.091 61.531 sp (a) Bayes landmark 0.383 0.361 0.401 0.187 0.175 0.199 0.000 0.000 True 0.136 66.258 sp (a) Bayes alignment 0.711 0.675 0.745 0.581 0.545 0.613 0.001 0.001 True 0.057 31.741 sp (b) HeadingDecoder all tuned 0.215 0.178 0.257 0.144 0.112 0.176 0.231 0.261 True 0.149 82.373 sp (b) HeadingDecoder heading 0.146 0.130 0.164 0.124 0.103 0.146 0.907 0.934 True 0.106 87.329 sp (b) HeadingDecoder landmark 0.247 0.225 0.267 0.123 0.115 0.131 0.000 0.000 True 0.151 77.356 sp (b) HeadingDecoder alignment 0.219 0.181 0.259 0.149 0.116 0.183 0.237 0.233 True 0.159 82.748 Chance accuracy with each decoder (20 draws of shifted test activity, light-on): input decoder chance_mean chance_max dff (a) Bayes 0.102 0.128 dff (b) HeadingDecoder 0.101 0.120 sp (a) Bayes 0.101 0.115 sp (b) HeadingDecoder 0.101 0.135 Chance readings of the claims file: (a) 0.18 (Methods), (b) 0.10 (Figure 6A legend; 36/360)
fig, ax = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
for a, (inp, dec_name) in zip(ax, [("sp", "(a) Bayes"), ("sp", "(b) HeadingDecoder")]):
for k, gname in enumerate(groups):
acc = per_iter[(inp, dec_name, gname)]
a.errorbar([k - 0.1], [acc["on"].mean()], yerr=[[acc["on"].mean() - np.percentile(acc["on"], 2.5)], [np.percentile(acc["on"], 97.5) - acc["on"].mean()]], fmt="o", color="0.2")
a.errorbar([k + 0.1], [acc["off"].mean()], yerr=[[acc["off"].mean() - np.percentile(acc["off"], 2.5)], [np.percentile(acc["off"], 97.5) - acc["off"].mean()]], fmt="s", color="0.6")
a.axhline(0.10, ls=":", color="k"); a.axhline(0.18, ls="--", color="k")
a.set_xticks(range(len(groups))); a.set_xticklabels(list(groups), rotation=20); a.set_title(f"spikes, {dec_name}")
ax[0].set_ylabel("decoder accuracy (light-on dark, light-off grey)"); ax[0].set_ylim(0, 1)
plt.tight_layout(); plt.show()
print("All accuracies lie in", round(float(res[["on_lo", "off_lo"]].min().min()), 3), "to", round(float(res[["on_hi", "off_hi"]].max().max()), 3), "; none outside the axes.")
All accuracies lie in 0.086 to 0.761 ; none outside the axes.
Checks of the decoders (asked for in verification)¶
Three checks bear on the choice of decoder reading.
- Both decoders as written here are run on synthetic cells whose tuning is known, with as many trials and bins as the real pseudopopulation, and the HeadingDecoder is also run on the real data without its temporal filter. If the reimplementation works, it decodes the synthetic cells well, and its low accuracy on the real cells comes from the data or from the reading, not from a coding error.
- The resampled heading bins are 1.8 degrees apart and each is interpolated from the two nearest frames, so neighbouring even (training) and odd (test) bins share frames. The paper does the same: it resamples each trial to 200 bins and then splits "even and odd timepoints for training and testing" (Methods, Decoding RSC somatic activity, "Preprocessing and alignment" and "Decoding strategy"), so reading (a) follows the paper. Variant (a') trains the naive Bayes decoder on the even trials and tests it on the odd trials, so that no frame and no trial is in both.
- Variant (a'') keeps the paper's split of even and odd time points but removes the interpolation: each trial is cut into 100 heading bins of 3.6 degrees, each holding the one frame nearest to it, so the even and odd bins hold different frames (the share they still share is printed). Neighbouring frames remain correlated through the calcium dynamics, as in any split of time points within a trial.
Both variants use the same bootstrap and enter every claim's table below as further readings of (a).
# (1) Both decoders on synthetic cells of known tuning, shaped like the pseudopopulation (spikes)
rng = np.random.default_rng(603)
n_tr_syn = pop["sp"]["on"].shape[0]
syn_rows = []
for shape in ["one peak", "two peaks 180 deg apart, the second half as high"]:
pref = rng.uniform(-180, 180, N_SAMPLE)
dd = np.deg2rad(grid[:, None] - pref[None])
rate = np.exp(3 * (np.cos(dd) - 1))
if shape != "one peak":
rate = rate + 0.5 * np.exp(3 * (np.cos(dd - np.pi) - 1))
for gain in [0.2, 1.0]:
A = rng.poisson(gain * rate[None], size=(n_tr_syn, len(grid), N_SAMPLE)).astype(float)
idx = np.arange(N_SAMPLE)
row_ = dict(tuning=shape, mean_events_per_bin=round(float(A.mean()), 3))
for dec_name, prep in [("Bayes", prep_bayes), ("HeadingDecoder", prep_hd)]:
e = decode(prep(A, "sp"), idx)
row_[dec_name] = round(float((e <= 18).mean()), 3)
row_[dec_name + " errors >= 162"] = round(float((e >= 162).mean()), 3)
syn_rows.append(row_)
print("Synthetic cells (332 cells, von Mises tuning, Poisson events; accuracy within 18 degrees):")
print(pd.DataFrame(syn_rows).to_string(index=False))
print("Real cells: share of light-on spike bins above 0 (the Bayes binarization):", round(float((pop["sp"]["on"] > 0).mean()), 3))
# HeadingDecoder without its temporal filter, on the real light-on spikes of all tuned cells (100 draws)
def decode_hd_nofilter(Pp, idx):
sc = np.matmul(Pp["X"][:, :, idx], Pp["W"][:, idx])
return err_deg(centres[sc.argmax(2)], Pp["true"][None])
rng = np.random.default_rng(604)
Ph = prep_hd(pop["sp"]["on"], "sp")
acc_f, acc_nf = [], []
for _ in range(100):
idx = rng.choice(len(ucls), N_SAMPLE, replace=True)
acc_f.append((decode(Ph, idx) <= 18).mean()); acc_nf.append((decode_hd_nofilter(Ph, idx) <= 18).mean())
print(f"HeadingDecoder on real spikes, all tuned, light-on (100 draws): with temporal filter {np.mean(acc_f):.3f}, "
f"without {np.mean(acc_nf):.3f}")
Synthetic cells (332 cells, von Mises tuning, Poisson events; accuracy within 18 degrees):
tuning mean_events_per_bin Bayes Bayes errors >= 162 HeadingDecoder HeadingDecoder errors >= 162
one peak 0.049 0.743 0.000 0.938 0.000
one peak 0.243 0.982 0.000 0.999 0.000
two peaks 180 deg apart, the second half as high 0.073 0.663 0.085 0.908 0.008
two peaks 180 deg apart, the second half as high 0.365 0.979 0.006 0.973 0.000
Real cells: share of light-on spike bins above 0 (the Bayes binarization): 0.16
HeadingDecoder on real spikes, all tuned, light-on (100 draws): with temporal filter 0.216, without 0.197
# (2) Naive Bayes trained on the even trials and tested on the odd trials (no frame shared), and
# (3) on even and odd time points without interpolation: each of 100 heading bins (3.6 degrees) takes the frame
# nearest to it, so the even and odd bins of a trial hold different frames, as in the paper's split of time points
def prep_bayes_general(A, inp, st, true, split):
n = A.shape[2]
flat = A.reshape(-1, n)
th = np.zeros(n) if inp == "sp" else flat.mean(0) + 2 * flat.std(0)
B = (A > th).astype(np.float32)
if split == "trials":
tr, te, s_tr, true_te = B[0::2], B[1::2], st, true
else:
tr, te, s_tr, true_te = B[:, 0::2], B[:, 1::2], st[0::2], true[1::2]
p = np.stack([(tr[:, s_tr == s].sum((0, 1)) + 1) / (tr[:, s_tr == s].shape[0] * tr[:, s_tr == s].shape[1] + 2)
for s in range(NS)], 1)
return dict(kind="bayes", W=(np.log(p) - np.log(1 - p)).astype(np.float32), B=np.log(1 - p).astype(np.float32),
X=te, true=true_te)
N_NF = 100
grid_nf = -180 + (np.arange(N_NF) + 0.5) * 360 / N_NF
state_nf = np.clip(((grid_nf + 180) // 6).astype(int), 0, NS - 1)
def nearest_frames(h, seg):
hu = np.rad2deg(np.unwrap(np.deg2rad(h[seg])))
sgn = 1 if hu[-1] >= hu[0] else -1
hu = sgn * hu
target = hu[0] + np.mod(sgn * grid_nf - hu[0], 360)
return seg[np.abs(hu[None, :] - target[:, None]).argmin(1)]
pop_nf, shared = {}, []
for inp in ["dff", "sp"]:
on_all, off_all = [], []
for i in use:
on, off = trials_around(i)
act = recs[i][inp][:, cell_roi[cell_rec == i]]
hc = rsc.clean_heading(recs[i]["heading"])
fr_on = [nearest_frames(hc, s) for s in on]
fr_off = [nearest_frames(hc, s) for s in off]
if inp == "sp":
shared += [len(np.intersect1d(f[0::2], f[1::2])) / (N_NF // 2) for f in fr_on + fr_off]
on_all.append(np.stack([act[f] for f in fr_on]))
off_all.append(np.stack([act[f] for f in fr_off]))
pop_nf[inp] = {"on": np.concatenate(on_all, 2), "off": np.concatenate(off_all, 2)}
print(f"nearest-frame bins: share of odd (test) bins whose frame is also an even (training) bin's frame: "
f"mean {np.mean(shared):.3f}, max {np.max(shared):.3f} over {len(shared)} trials")
def bayes_boot(make, label, seed):
rows_ = []
for inp in ["dff", "sp"]:
P = {cond: make(inp, cond) for cond in ["on", "off"]}
rng = np.random.default_rng(seed)
for gname, m in groups.items():
pool_ = np.flatnonzero(m)
acc = {"on": np.zeros(N_BOOT), "off": np.zeros(N_BOOT)}
e180 = np.zeros(N_BOOT); pt_paired = np.zeros(N_BOOT, bool); pt_unpaired = np.zeros(N_BOOT, bool)
for it in range(N_BOOT):
idx = rng.choice(pool_, N_SAMPLE, replace=True)
a = {}
for cond in ["on", "off"]:
e = decode(P[cond], idx)
a[cond] = (e <= 18).mean(1)
acc[cond][it] = a[cond].mean()
if cond == "on":
e180[it] = (e >= 162).mean()
pt_paired[it] = stats.ttest_rel(a["on"], a["off"]).pvalue < 0.05
pt_unpaired[it] = stats.ttest_ind(a["on"], a["off"]).pvalue < 0.05
per_iter[(inp, label, gname)] = acc
rows_.append(dict(input=inp, decoder=label, group=gname,
on_mean=acc["on"].mean(), on_lo=np.percentile(acc["on"], 2.5), on_hi=np.percentile(acc["on"], 97.5),
off_mean=acc["off"].mean(), off_lo=np.percentile(acc["off"], 2.5), off_hi=np.percentile(acc["off"], 97.5),
p_paired=1 - pt_paired.mean(), p_unpaired=1 - pt_unpaired.mean(),
off_lower=bool(acc["off"].mean() < acc["on"].mean()), err180_on=e180.mean(), mean_err_on=np.nan))
return rows_
t0 = time.time()
rows_tr = bayes_boot(lambda inp, cond: prep_bayes_general(pop[inp][cond], inp, state, grid, "trials"), "(a') Bayes, trial split", 602)
rows_tr += bayes_boot(lambda inp, cond: prep_bayes_general(pop_nf[inp][cond], inp, state_nf, grid_nf, "points"),
"(a'') Bayes, nearest frames", 605)
res_tr = pd.DataFrame(rows_tr).round(3)
print(f"variant decoding in {time.time() - t0:.0f} s; trial split: {pop['sp']['on'][0::2].shape[0]} training and "
f"{pop['sp']['on'][1::2].shape[0]} test trials per condition; nearest frames: {N_NF} bins per trial, even to train, odd to test")
print(res_tr.drop(columns="mean_err_on").to_string(index=False))
print("\nDifference from the even/odd-bin split of (a) (variant minus (a)), light-on and light-off means:")
cmp_ = res_tr.merge(res[res.decoder == "(a) Bayes"], on=["input", "group"], suffixes=("", "_bins"))
print(pd.DataFrame(dict(variant=cmp_.decoder, input=cmp_.input, group=cmp_.group, on=(cmp_.on_mean - cmp_.on_mean_bins).round(3),
off=(cmp_.off_mean - cmp_.off_mean_bins).round(3))).to_string(index=False))
res = pd.concat([res, res_tr], ignore_index=True)
nearest-frame bins: share of odd (test) bins whose frame is also an even (training) bin's frame: mean 0.002, max 0.040 over 750 trials
variant decoding in 63 s; trial split: 8 training and 7 test trials per condition; nearest frames: 100 bins per trial, even to train, odd to test
input decoder group on_mean on_lo on_hi off_mean off_lo off_hi p_paired p_unpaired off_lower err180_on
dff (a') Bayes, trial split all tuned 0.255 0.210 0.301 0.167 0.137 0.202 0.348 0.284 True 0.170
dff (a') Bayes, trial split heading 0.158 0.131 0.187 0.162 0.139 0.184 0.991 0.988 False 0.157
dff (a') Bayes, trial split landmark 0.253 0.228 0.278 0.138 0.126 0.150 0.014 0.000 True 0.179
dff (a') Bayes, trial split alignment 0.259 0.215 0.308 0.160 0.129 0.192 0.250 0.158 True 0.168
sp (a') Bayes, trial split all tuned 0.239 0.207 0.274 0.140 0.116 0.165 0.047 0.055 True 0.201
sp (a') Bayes, trial split heading 0.161 0.139 0.182 0.133 0.116 0.151 0.857 0.955 True 0.124
sp (a') Bayes, trial split landmark 0.230 0.209 0.251 0.104 0.093 0.115 0.000 0.000 True 0.194
sp (a') Bayes, trial split alignment 0.239 0.206 0.272 0.136 0.114 0.159 0.023 0.023 True 0.202
dff (a'') Bayes, nearest frames all tuned 0.532 0.487 0.579 0.406 0.361 0.452 0.105 0.175 True 0.050
dff (a'') Bayes, nearest frames heading 0.328 0.297 0.357 0.280 0.253 0.305 0.806 0.917 True 0.078
dff (a'') Bayes, nearest frames landmark 0.352 0.333 0.371 0.165 0.148 0.184 0.000 0.000 True 0.135
dff (a'') Bayes, nearest frames alignment 0.548 0.503 0.596 0.414 0.367 0.464 0.065 0.140 True 0.048
sp (a'') Bayes, nearest frames all tuned 0.296 0.263 0.332 0.206 0.177 0.233 0.062 0.049 True 0.136
sp (a'') Bayes, nearest frames heading 0.202 0.177 0.225 0.173 0.151 0.195 0.837 0.817 True 0.137
sp (a'') Bayes, nearest frames landmark 0.250 0.224 0.275 0.132 0.117 0.147 0.000 0.000 True 0.171
sp (a'') Bayes, nearest frames alignment 0.295 0.261 0.331 0.197 0.171 0.229 0.022 0.021 True 0.133
Difference from the even/odd-bin split of (a) (variant minus (a)), light-on and light-off means:
variant input group on off
(a') Bayes, trial split dff all tuned -0.453 -0.402
(a') Bayes, trial split dff heading -0.271 -0.215
(a') Bayes, trial split dff landmark -0.143 -0.069
(a') Bayes, trial split dff alignment -0.461 -0.420
(a') Bayes, trial split sp all tuned -0.473 -0.432
(a') Bayes, trial split sp heading -0.254 -0.208
(a') Bayes, trial split sp landmark -0.153 -0.083
(a') Bayes, trial split sp alignment -0.472 -0.445
(a'') Bayes, nearest frames dff all tuned -0.176 -0.163
(a'') Bayes, nearest frames dff heading -0.101 -0.097
(a'') Bayes, nearest frames dff landmark -0.044 -0.042
(a'') Bayes, nearest frames dff alignment -0.172 -0.166
(a'') Bayes, nearest frames sp all tuned -0.416 -0.366
(a'') Bayes, nearest frames sp heading -0.213 -0.168
(a'') Bayes, nearest frames sp landmark -0.133 -0.055
(a'') Bayes, nearest frames sp alignment -0.416 -0.384
What the checks show. The reimplemented HeadingDecoder decodes synthetic cells well (0.908 to 0.999 within 18 degrees, better than the naive Bayes decoder's 0.663 to 0.982), so its low accuracy on the real cells (0.202 to 0.215 in light-on; 0.197 without its temporal filter) is not a coding error in the algorithm as written here. It cannot rule out a detail in which the lab's code differs. Its accuracy is close to that of the naive Bayes decoder trained and tested on different trials (0.239 and 0.255, all tuned cells, light-on); both test on trials not used for training. Removing the interpolation from the paper's split of time points (a''; 0.2% of test bins share a frame with a training bin) lowers light-on accuracy of all tuned cells from 0.708 to 0.532 (dF/F) and from 0.712 to 0.296 (spikes). So the accuracy of reading (a) depends heavily on test bins that share frames, or neighbouring frames, with training bins. Reading (a) is still the one the Methods describe (resampling to 200 bins, then even and odd time points), and it is the reading on which the Bayes outcomes below rest. Each claim below also reports the two variants.
def row(inp, dec_name, gname):
return res[(res.input == inp) & (res.decoder == dec_name) & (res.group == gname)].iloc[0]
READ = [(inp, d) for inp in ["dff", "sp"] for d in ["(a) Bayes", "(a') Bayes, trial split", "(a'') Bayes, nearest frames", "(b) HeadingDecoder"]]
def excludes(lo, hi, c):
return not (lo <= c <= hi)
fig6a-population-light-on¶
Decoding heading from all tuned RSC cells within single 100 ms frames is accurate in the light-on condition despite the symmetric cues (mean decoder accuracy 0.82, 95% CI 0.78 to 0.86).
Reproduced when a reading of the decoder gives a mean within 0.82 +- 0.04.
for inp, d in READ:
r = row(inp, d, "all tuned")
print(f"input {inp}, decoder {d}: light-on mean {r.on_mean:.3f} [{r.on_lo:.3f}, {r.on_hi:.3f}]; within 0.82 +- 0.04: {abs(r.on_mean - 0.82) <= 0.04}")
input dff, decoder (a) Bayes: light-on mean 0.708 [0.665, 0.747]; within 0.82 +- 0.04: False input dff, decoder (a') Bayes, trial split: light-on mean 0.255 [0.210, 0.301]; within 0.82 +- 0.04: False input dff, decoder (a'') Bayes, nearest frames: light-on mean 0.532 [0.487, 0.579]; within 0.82 +- 0.04: False input dff, decoder (b) HeadingDecoder: light-on mean 0.202 [0.158, 0.251]; within 0.82 +- 0.04: False input sp, decoder (a) Bayes: light-on mean 0.712 [0.678, 0.746]; within 0.82 +- 0.04: False input sp, decoder (a') Bayes, trial split: light-on mean 0.239 [0.207, 0.274]; within 0.82 +- 0.04: False input sp, decoder (a'') Bayes, nearest frames: light-on mean 0.296 [0.263, 0.332]; within 0.82 +- 0.04: False input sp, decoder (b) HeadingDecoder: light-on mean 0.215 [0.178, 0.257]; within 0.82 +- 0.04: False
Result: not reproduced. No decoder reading gives a light-on mean within 0.82 +- 0.04: the naive Bayes decoder gives 0.708 [0.665, 0.747] (dF/F) and 0.712 [0.678, 0.746] (spikes), the variants without shared frames 0.239 to 0.532, and the HeadingDecoder 0.202 and 0.215. Chance with each decoder is 0.10 (means 0.101 to 0.102), the 0.10 reading of the claims file, so every reading is above chance but less accurate than the paper.
fig6a-population-light-off¶
In the light-off condition the decoder's accuracy from all tuned cells drops (mean 0.44, 95% CI 0.34 to 0.54; p = 7.6e-3 against light-on) but stays above chance.
Reproduced when, under a decoder reading, the mean is within 0.44 +- 0.10, the interval excludes chance (0.18 or 0.10) and the bootstrapped t-test gives p < 0.05 with light-off lower.
for inp, d in READ:
r = row(inp, d, "all tuned")
ok = abs(r.off_mean - 0.44) <= 0.10 and excludes(r.off_lo, r.off_hi, 0.18) and excludes(r.off_lo, r.off_hi, 0.10) and r.off_lower
print(f"input {inp}, decoder {d}: light-off mean {r.off_mean:.3f} [{r.off_lo:.3f}, {r.off_hi:.3f}]; p paired {r.p_paired:.3f}, "
f"unpaired {r.p_unpaired:.3f}; within 0.44 +- 0.10: {abs(r.off_mean - 0.44) <= 0.10}; above both chance levels: "
f"{r.off_lo > 0.18}; criterion met (paired / unpaired): {ok and r.p_paired < 0.05} / {ok and r.p_unpaired < 0.05}")
input dff, decoder (a) Bayes: light-off mean 0.569 [0.525, 0.612]; p paired 0.016, unpaired 0.061; within 0.44 +- 0.10: False; above both chance levels: True; criterion met (paired / unpaired): False / False input dff, decoder (a') Bayes, trial split: light-off mean 0.167 [0.137, 0.202]; p paired 0.348, unpaired 0.284; within 0.44 +- 0.10: False; above both chance levels: False; criterion met (paired / unpaired): False / False input dff, decoder (a'') Bayes, nearest frames: light-off mean 0.406 [0.361, 0.452]; p paired 0.105, unpaired 0.175; within 0.44 +- 0.10: True; above both chance levels: True; criterion met (paired / unpaired): False / False input dff, decoder (b) HeadingDecoder: light-off mean 0.137 [0.101, 0.178]; p paired 0.346, unpaired 0.361; within 0.44 +- 0.10: False; above both chance levels: False; criterion met (paired / unpaired): False / False input sp, decoder (a) Bayes: light-off mean 0.572 [0.535, 0.608]; p paired 0.000, unpaired 0.002; within 0.44 +- 0.10: False; above both chance levels: True; criterion met (paired / unpaired): False / False input sp, decoder (a') Bayes, trial split: light-off mean 0.140 [0.116, 0.165]; p paired 0.047, unpaired 0.055; within 0.44 +- 0.10: False; above both chance levels: False; criterion met (paired / unpaired): False / False input sp, decoder (a'') Bayes, nearest frames: light-off mean 0.206 [0.177, 0.233]; p paired 0.062, unpaired 0.049; within 0.44 +- 0.10: False; above both chance levels: False; criterion met (paired / unpaired): False / False input sp, decoder (b) HeadingDecoder: light-off mean 0.144 [0.112, 0.176]; p paired 0.231, unpaired 0.261; within 0.44 +- 0.10: False; above both chance levels: False; criterion met (paired / unpaired): False / False
Result: partially reproduced, under the naive Bayes decoder (an assumption). Light-off accuracy drops and stays above chance: 0.569 [0.525, 0.612] (dF/F) and 0.572 [0.535, 0.608] (spikes), both intervals above 0.18 and 0.10, light-off lower with p = 0.016 and p < 0.001 (paired; unpaired 0.061 and 0.002). The mean is outside 0.44 +- 0.10 under reading (a), so the criterion is not met in full. Without shared frames, (a'') gives 0.406 [0.361, 0.452] for dF/F, within the tolerance and above both chance levels, but the drop is not significant (p = 0.105 paired); for spikes 0.206 [0.177, 0.233]. The trial split (a') gives 0.140 and 0.167, with intervals that include 0.18 (dF/F) or lie between 0.10 and 0.18 (spikes). The HeadingDecoder gives 0.137 [0.101, 0.178] and 0.144 [0.112, 0.176], above 0.10 and below 0.18, with no significant drop (p 0.231 to 0.361).
fig6b-heading-cells¶
Heading cells alone decode heading accurately in both light-on (0.61, 95% CI 0.56 to 0.66) and light-off (0.44, 95% CI 0.36 to 0.52) conditions, with no significant difference between them (p = 0.11).
Reproduced when a decoder reading puts light-on within 0.61 +- 0.05 and light-off within 0.44 +- 0.08, both intervals exclude chance, and the t-test gives p >= 0.05.
for inp, d in READ:
r = row(inp, d, "heading")
vals = abs(r.on_mean - 0.61) <= 0.05 and abs(r.off_mean - 0.44) <= 0.08
print(f"input {inp}, decoder {d}: on {r.on_mean:.3f} [{r.on_lo:.3f}, {r.on_hi:.3f}], off {r.off_mean:.3f} [{r.off_lo:.3f}, {r.off_hi:.3f}]; "
f"p paired {r.p_paired:.3f}, unpaired {r.p_unpaired:.3f}; both within tolerance: {vals}; "
f"both intervals above 0.18: {r.on_lo > 0.18 and r.off_lo > 0.18}, above 0.10: {r.on_lo > 0.10 and r.off_lo > 0.10}")
input dff, decoder (a) Bayes: on 0.429 [0.399, 0.457], off 0.377 [0.350, 0.403]; p paired 0.697, unpaired 0.964; both within tolerance: False; both intervals above 0.18: True, above 0.10: True input dff, decoder (a') Bayes, trial split: on 0.158 [0.131, 0.187], off 0.162 [0.139, 0.184]; p paired 0.991, unpaired 0.988; both within tolerance: False; both intervals above 0.18: False, above 0.10: True input dff, decoder (a'') Bayes, nearest frames: on 0.328 [0.297, 0.357], off 0.280 [0.253, 0.305]; p paired 0.806, unpaired 0.917; both within tolerance: False; both intervals above 0.18: True, above 0.10: True input dff, decoder (b) HeadingDecoder: on 0.130 [0.113, 0.153], off 0.109 [0.091, 0.131]; p paired 0.939, unpaired 0.915; both within tolerance: False; both intervals above 0.18: False, above 0.10: False input sp, decoder (a) Bayes: on 0.415 [0.390, 0.441], off 0.341 [0.317, 0.364]; p paired 0.170, unpaired 0.534; both within tolerance: False; both intervals above 0.18: True, above 0.10: True input sp, decoder (a') Bayes, trial split: on 0.161 [0.139, 0.182], off 0.133 [0.116, 0.151]; p paired 0.857, unpaired 0.955; both within tolerance: False; both intervals above 0.18: False, above 0.10: True input sp, decoder (a'') Bayes, nearest frames: on 0.202 [0.177, 0.225], off 0.173 [0.151, 0.195]; p paired 0.837, unpaired 0.817; both within tolerance: False; both intervals above 0.18: False, above 0.10: True input sp, decoder (b) HeadingDecoder: on 0.146 [0.130, 0.164], off 0.124 [0.103, 0.146]; p paired 0.907, unpaired 0.934; both within tolerance: False; both intervals above 0.18: False, above 0.10: True
Result: partially reproduced, under the naive Bayes decoder (an assumption), on a heading-like class that is not the paper's. Heading-like cells (179 used) decode above chance in both conditions (dF/F 0.429 and 0.377; spikes 0.415 and 0.341; every interval above 0.18) with no significant difference (p 0.170 to 0.964), as the claim says. The accuracies are below the paper's under every reading (outside 0.61 +- 0.05 and 0.44 +- 0.08). The same pattern, lower, holds without shared frames: (a'') dF/F 0.328 and 0.280, both intervals above 0.18, p = 0.806; the trial split and the HeadingDecoder give 0.130 to 0.161 in light-on, above 0.10, with p 0.857 to 0.991.
fig6b-landmark-cells¶
Landmark cells decode heading very well in light-on (0.80, 95% CI 0.75 to 0.84) but fall to chance in light-off (0.20, 95% CI 0.14 to 0.27; p = 8.0e-4).
Reproduced when, under a decoder reading, light-on is within 0.80 +- 0.05, light-off within 0.20 +- 0.07, p < 0.05, and the light-off interval includes chance; the outcome under each chance level is reported.
for inp, d in READ:
r = row(inp, d, "landmark")
vals = abs(r.on_mean - 0.80) <= 0.05 and abs(r.off_mean - 0.20) <= 0.07
print(f"input {inp}, decoder {d}: on {r.on_mean:.3f} [{r.on_lo:.3f}, {r.on_hi:.3f}], off {r.off_mean:.3f} [{r.off_lo:.3f}, {r.off_hi:.3f}]; "
f"p paired {r.p_paired:.3f}, unpaired {r.p_unpaired:.3f}; both within tolerance: {vals}; "
f"light-off interval includes 0.18: {r.off_lo <= 0.18 <= r.off_hi}, includes 0.10: {r.off_lo <= 0.10 <= r.off_hi}")
input dff, decoder (a) Bayes: on 0.396 [0.377, 0.414], off 0.207 [0.192, 0.220]; p paired 0.000, unpaired 0.000; both within tolerance: False; light-off interval includes 0.18: False, includes 0.10: False input dff, decoder (a') Bayes, trial split: on 0.253 [0.228, 0.278], off 0.138 [0.126, 0.150]; p paired 0.014, unpaired 0.000; both within tolerance: False; light-off interval includes 0.18: False, includes 0.10: False input dff, decoder (a'') Bayes, nearest frames: on 0.352 [0.333, 0.371], off 0.165 [0.148, 0.184]; p paired 0.000, unpaired 0.000; both within tolerance: False; light-off interval includes 0.18: True, includes 0.10: False input dff, decoder (b) HeadingDecoder: on 0.201 [0.171, 0.234], off 0.103 [0.086, 0.120]; p paired 0.003, unpaired 0.001; both within tolerance: False; light-off interval includes 0.18: False, includes 0.10: True input sp, decoder (a) Bayes: on 0.383 [0.361, 0.401], off 0.187 [0.175, 0.199]; p paired 0.000, unpaired 0.000; both within tolerance: False; light-off interval includes 0.18: True, includes 0.10: False input sp, decoder (a') Bayes, trial split: on 0.230 [0.209, 0.251], off 0.104 [0.093, 0.115]; p paired 0.000, unpaired 0.000; both within tolerance: False; light-off interval includes 0.18: False, includes 0.10: True input sp, decoder (a'') Bayes, nearest frames: on 0.250 [0.224, 0.275], off 0.132 [0.117, 0.147]; p paired 0.000, unpaired 0.000; both within tolerance: False; light-off interval includes 0.18: False, includes 0.10: False input sp, decoder (b) HeadingDecoder: on 0.247 [0.225, 0.267], off 0.123 [0.115, 0.131]; p paired 0.000, unpaired 0.000; both within tolerance: False; light-off interval includes 0.18: False, includes 0.10: False
Result: partially reproduced, under the naive Bayes decoder (an assumption), on a landmark-like class that is not the paper's. Landmark-like cells (50 used, resampled to 332) fall towards chance in light-off, 0.207 [0.192, 0.220] (dF/F) and 0.187 [0.175, 0.199] (spikes), within 0.20 +- 0.07, with p < 0.001 for the drop. Their light-on accuracy, 0.396 and 0.383, is half the paper's 0.80, so the criterion is not met. The spikes light-off interval includes 0.18; the dF/F interval excludes it, and neither includes 0.10. Every other reading also gives a significant drop (p 0.014 or less) and light-on accuracy of 0.201 to 0.352; the light-off interval includes 0.18 under (a'') dF/F and 0.10 under the trial split for spikes and the HeadingDecoder for dF/F.
fig6b-alignment-cells¶
Alignment cells decode between heading and landmark cells, better in light-on (0.71, 95% CI 0.66 to 0.76) and fairly accurately in light-off (0.51, 95% CI 0.43 to 0.59), with no significant difference (p = 0.54).
Reproduced when, under a decoder reading, light-on is within 0.71 +- 0.05 and light-off within 0.51 +- 0.08, the light-on mean lies between those of heading and landmark cells, the light-off interval excludes chance, and p >= 0.05.
for inp, d in READ:
r, rh, rl = row(inp, d, "alignment"), row(inp, d, "heading"), row(inp, d, "landmark")
vals = abs(r.on_mean - 0.71) <= 0.05 and abs(r.off_mean - 0.51) <= 0.08
between = min(rh.on_mean, rl.on_mean) < r.on_mean < max(rh.on_mean, rl.on_mean)
print(f"input {inp}, decoder {d}: on {r.on_mean:.3f} [{r.on_lo:.3f}, {r.on_hi:.3f}], off {r.off_mean:.3f} [{r.off_lo:.3f}, {r.off_hi:.3f}]; "
f"p paired {r.p_paired:.3f}, unpaired {r.p_unpaired:.3f}; both within tolerance: {vals}; light-on between heading "
f"({rh.on_mean:.3f}) and landmark ({rl.on_mean:.3f}): {between}; light-off above 0.18: {r.off_lo > 0.18}, above 0.10: {r.off_lo > 0.10}")
input dff, decoder (a) Bayes: on 0.720 [0.680, 0.761], off 0.580 [0.538, 0.619]; p paired 0.011, unpaired 0.053; both within tolerance: True; light-on between heading (0.429) and landmark (0.396): False; light-off above 0.18: True, above 0.10: True input dff, decoder (a') Bayes, trial split: on 0.259 [0.215, 0.308], off 0.160 [0.129, 0.192]; p paired 0.250, unpaired 0.158; both within tolerance: False; light-on between heading (0.158) and landmark (0.253): False; light-off above 0.18: False, above 0.10: True input dff, decoder (a'') Bayes, nearest frames: on 0.548 [0.503, 0.596], off 0.414 [0.367, 0.464]; p paired 0.065, unpaired 0.140; both within tolerance: False; light-on between heading (0.328) and landmark (0.352): False; light-off above 0.18: True, above 0.10: True input dff, decoder (b) HeadingDecoder: on 0.203 [0.160, 0.251], off 0.151 [0.110, 0.201]; p paired 0.530, unpaired 0.527; both within tolerance: False; light-on between heading (0.130) and landmark (0.201): False; light-off above 0.18: False, above 0.10: True input sp, decoder (a) Bayes: on 0.711 [0.675, 0.745], off 0.581 [0.545, 0.613]; p paired 0.001, unpaired 0.001; both within tolerance: True; light-on between heading (0.415) and landmark (0.383): False; light-off above 0.18: True, above 0.10: True input sp, decoder (a') Bayes, trial split: on 0.239 [0.206, 0.272], off 0.136 [0.114, 0.159]; p paired 0.023, unpaired 0.023; both within tolerance: False; light-on between heading (0.161) and landmark (0.230): False; light-off above 0.18: False, above 0.10: True input sp, decoder (a'') Bayes, nearest frames: on 0.295 [0.261, 0.331], off 0.197 [0.171, 0.229]; p paired 0.022, unpaired 0.021; both within tolerance: False; light-on between heading (0.202) and landmark (0.250): False; light-off above 0.18: False, above 0.10: True input sp, decoder (b) HeadingDecoder: on 0.219 [0.181, 0.259], off 0.149 [0.116, 0.183]; p paired 0.237, unpaired 0.233; both within tolerance: False; light-on between heading (0.146) and landmark (0.247): True; light-off above 0.18: False, above 0.10: True
Result: not reproduced. Under the naive Bayes decoder both means of alignment-like cells are within tolerance (dF/F 0.720 and 0.580; spikes 0.711 and 0.581) and the light-off intervals exclude chance, but their light-on accuracy is the highest of the three classes, not between heading-like (0.429, 0.415) and landmark-like (0.396, 0.383) cells, and light-off is significantly lower by the paired test (p = 0.011 and 0.001; unpaired 0.053 and 0.001). No other reading puts both means within tolerance (light-on 0.203 to 0.548). Under the HeadingDecoder with spikes the light-on mean lies between the other classes (0.219, between 0.146 and 0.247), but far below 0.71.
fig6b-180-errors¶
Decoding from landmark cells in light-on gives an overrepresentation of 180 degree errors from the symmetric cues, while heading cells show no preponderance of 180 degree errors.
The share of light-on bins with an error within 18 degrees of 180, against 0.10 under uniform errors. Reproduced when it exceeds 0.10 for landmark cells and is no more than 0.10 for heading cells.
for inp, d in READ:
rl, rh, ra, rt = row(inp, d, "landmark"), row(inp, d, "heading"), row(inp, d, "alignment"), row(inp, d, "all tuned")
print(f"input {inp}, decoder {d}: share of errors >= 162 deg: landmark {rl.err180_on:.3f}, heading {rh.err180_on:.3f}, "
f"alignment {ra.err180_on:.3f}, all tuned {rt.err180_on:.3f}; criterion met: {rl.err180_on > 0.10 and rh.err180_on <= 0.10}")
input dff, decoder (a) Bayes: share of errors >= 162 deg: landmark 0.120, heading 0.070, alignment 0.033, all tuned 0.034; criterion met: True input dff, decoder (a') Bayes, trial split: share of errors >= 162 deg: landmark 0.179, heading 0.157, alignment 0.168, all tuned 0.170; criterion met: False input dff, decoder (a'') Bayes, nearest frames: share of errors >= 162 deg: landmark 0.135, heading 0.078, alignment 0.048, all tuned 0.050; criterion met: True input dff, decoder (b) HeadingDecoder: share of errors >= 162 deg: landmark 0.115, heading 0.104, alignment 0.148, all tuned 0.138; criterion met: False input sp, decoder (a) Bayes: share of errors >= 162 deg: landmark 0.136, heading 0.091, alignment 0.057, all tuned 0.055; criterion met: True input sp, decoder (a') Bayes, trial split: share of errors >= 162 deg: landmark 0.194, heading 0.124, alignment 0.202, all tuned 0.201; criterion met: False input sp, decoder (a'') Bayes, nearest frames: share of errors >= 162 deg: landmark 0.171, heading 0.137, alignment 0.133, all tuned 0.136; criterion met: False input sp, decoder (b) HeadingDecoder: share of errors >= 162 deg: landmark 0.151, heading 0.106, alignment 0.159, all tuned 0.149; criterion met: False
Result: reproduced, under the naive Bayes decoder (an assumption), on classes that are not the paper's. In light-on, 12.0% (dF/F) and 13.6% (spikes) of the landmark-like decoder's errors lie within 18 degrees of 180, above 0.10, and 7.0% and 9.1% of the heading-like decoder's, no more than 0.10. The result is fragile. Without shared frames it holds for dF/F (a'': landmark 0.135, heading 0.078) but not for spikes (heading 0.137). Under the trial split and the HeadingDecoder it fails (heading-like 0.104 to 0.157); those readings decode at 0.13 to 0.26, so their errors are spread widely: all their groups have 0.104 to 0.202.
fig6-heading-below-population¶
In light-on, heading cells alone decode worse than the whole tuned population, so alone they may not give an optimal representation of instantaneous heading.
Reading (a): the 95% bootstrap intervals of light-on accuracy for heading cells and for all tuned cells do not overlap. Reading (b), a mixed model of per-trial decoder error with recording or mouse as random intercept, cannot be formed here: each pseudopopulation trial pools cells of many recordings and mice, so no trial belongs to one recording or mouse. As a variant, the share of bootstrap iterations in which heading cells decode worse than all tuned cells is given (iterations paired by index).
for inp, d in READ:
rh, rt = row(inp, d, "heading"), row(inp, d, "all tuned")
worse = np.mean(per_iter[(inp, d, "heading")]["on"] < per_iter[(inp, d, "all tuned")]["on"])
print(f"input {inp}, decoder {d}: heading {rh.on_mean:.3f} [{rh.on_lo:.3f}, {rh.on_hi:.3f}], all tuned {rt.on_mean:.3f} "
f"[{rt.on_lo:.3f}, {rt.on_hi:.3f}]; intervals separate with heading lower: {rh.on_hi < rt.on_lo}; "
f"iterations with heading lower: {worse:.3f}")
input dff, decoder (a) Bayes: heading 0.429 [0.399, 0.457], all tuned 0.708 [0.665, 0.747]; intervals separate with heading lower: True; iterations with heading lower: 1.000 input dff, decoder (a') Bayes, trial split: heading 0.158 [0.131, 0.187], all tuned 0.255 [0.210, 0.301]; intervals separate with heading lower: True; iterations with heading lower: 1.000 input dff, decoder (a'') Bayes, nearest frames: heading 0.328 [0.297, 0.357], all tuned 0.532 [0.487, 0.579]; intervals separate with heading lower: True; iterations with heading lower: 1.000 input dff, decoder (b) HeadingDecoder: heading 0.130 [0.113, 0.153], all tuned 0.202 [0.158, 0.251]; intervals separate with heading lower: True; iterations with heading lower: 0.998 input sp, decoder (a) Bayes: heading 0.415 [0.390, 0.441], all tuned 0.712 [0.678, 0.746]; intervals separate with heading lower: True; iterations with heading lower: 1.000 input sp, decoder (a') Bayes, trial split: heading 0.161 [0.139, 0.182], all tuned 0.239 [0.207, 0.274]; intervals separate with heading lower: True; iterations with heading lower: 1.000 input sp, decoder (a'') Bayes, nearest frames: heading 0.202 [0.177, 0.225], all tuned 0.296 [0.263, 0.332]; intervals separate with heading lower: True; iterations with heading lower: 1.000 input sp, decoder (b) HeadingDecoder: heading 0.146 [0.130, 0.164], all tuned 0.215 [0.178, 0.257]; intervals separate with heading lower: True; iterations with heading lower: 1.000
Result: reproduced. Under every decoder and input reading the 95% interval of light-on accuracy for heading-like cells lies below that of all tuned cells (Bayes: 0.429 [0.399, 0.457] against 0.708 [0.665, 0.747] for dF/F, 0.415 [0.390, 0.441] against 0.712 [0.678, 0.746] for spikes; the variants without shared frames 0.158 to 0.328 against 0.239 to 0.532; HeadingDecoder 0.130 and 0.146 against 0.202 and 0.215), and heading-like cells decode worse in 99.8 to 100% of iterations. Reading (b), the mixed model, cannot be formed for pseudopopulation trials. The heading-like class is not the paper's.
Summary¶
| Claim | Outcome |
|---|---|
| fig6a-population-light-on | not reproduced: 0.71 (Bayes), 0.24 to 0.53 (Bayes without shared frames) or 0.20 to 0.22 (HeadingDecoder) against 0.82 |
| fig6a-population-light-off | partially reproduced, assumed (Bayes): drop significant and above chance; 0.57 against 0.44 |
| fig6b-heading-cells | partially reproduced, assumed (Bayes): above chance in both, no significant difference; accuracies below the paper's |
| fig6b-landmark-cells | partially reproduced, assumed (Bayes): falls to 0.19 to 0.21 in light-off; light-on 0.38 to 0.40 against 0.80 |
| fig6b-alignment-cells | not reproduced: values within tolerance, but not between the other classes and a significant drop |
| fig6b-180-errors | reproduced, assumed (Bayes); fails for spikes without shared frames and under the other readings |
| fig6-heading-below-population | reproduced under every decoder reading |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's: the classes are those of somata.ipynb, not the
paper's; spikes and dF/F are computed from raw fluorescence without neuropil correction; only the
7550-frame recordings with 5 trials on each side of every transition enter the
pseudopopulations; trials are resampled on heading, not time; the naive Bayes decoder is written
here after Etter et al. (2020), with a binarization threshold that is a standard choice; the
HeadingDecoder is reimplemented in numpy. The accuracies of the naive Bayes decoder depend heavily on
test bins sharing frames with training bins, as the paper's own resampling and split of time
points do; two variants without shared frames are shown.
What was found about the dataset: the files do not mark light-on and light-off; the segments are the code's hard-coded ones.