Sit & Goard (2023) on dandiset 000483: RSC somata with dual cues, light and dark (Figure 5)¶
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 |
|---|---|
| fig5-three-classes | tested |
| fig5-class-proportions | tested |
| fig5e-alignment-peaks-coregistered | tested |
| fig5f-heading-flip-unchanged | tested |
| fig5gh-landmark-alignment-flip-drop | tested |
| fig5i-axons-in-clusters | tested |
The axon claims of Figure 4 are in axons.ipynb; Figures S5 and S6 are in directions.ipynb;
Figure 5K is in directions.ipynb (blocked) and the decoding of Figure 6 in decoding.ipynb. The axons are loaded here too, for fig5i-axons-in-clusters,
and tested with the same seed as in axons.ipynb. Activity is inferred spikes computed from the raw fluorescence (no
neuropil correction; see rsc.py). Light and dark segments are the code's hard-coded ones,
as the files do not mark them.
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
import statsmodels.formula.api as smf
from scipy import stats
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.metrics import silhouette_score
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.isin(["soma_dual", "adn_axons", "vis_axons"])].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: 59 files loaded and deconvolved in 236 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"])
step = rsc.wrap180(np.diff(rsc.clean_heading(r["heading"])))
# activity per segment, to see whether the hard-coded light and dark segments differ
act = r["sp"].sum(1)
rows.append(dict(experiment=m.experiment, subject=m.subject, identifier=m.identifier[:6] + " " + m.identifier.split("_")[-1],
frames=m.n_frames, rois=m.n_rois, median_step=round(float(np.median(step)), 2),
trials_all=len(r["seg_all"]), trials_light=len(r["seg_l"]), trials_dark=len(r["seg_d"]),
spikes_light=round(float(np.mean([act[s].mean() for s in L])), 2),
spikes_dark=round(float(np.mean([act[s].mean() for s in D])), 2)))
info = pd.DataFrame(rows)
print(info.to_string(index=False))
print(info.groupby("experiment").agg(files=("rois", "size"), mice=("subject", "nunique"), rois=("rois", "sum")).to_string())
experiment subject identifier frames rois median_step trials_all trials_light trials_dark spikes_light spikes_dark
adn_axons SKKS092 210217 SKKS092-axons-001 7550 27 1.93 38 17 19 0.15 0.09
adn_axons SKKS092 210217 SKKS092-axons-002 7550 60 1.93 38 19 18 0.45 0.24
adn_axons SKKS092 210304 SKKS092-ADNAxons-3rpm-001 7550 25 1.90 38 18 15 0.61 0.74
adn_axons SKKS092 210304 SKKS092-ADNAxons-3rpm-002 7550 26 1.89 39 18 19 0.47 0.65
adn_axons SKKS092 210304 SKKS092-ADNAxons-3rpm-003 7550 20 1.87 37 17 18 0.39 0.56
adn_axons SKKS093 210217 SKKS093-axons-001 7550 76 1.95 38 19 17 0.64 0.48
adn_axons SKKS093 210217 SKKS093-axons-002 7550 60 1.94 39 18 19 0.50 0.20
adn_axons SKKS093 210304 SKKS093-ADNAxons-3rpm-001 7550 67 1.89 39 18 20 0.90 0.64
adn_axons SKKS093 210304 SKKS093-ADNAxons-3rpm-002 7550 58 1.90 37 18 19 0.78 0.41
adn_axons SKKS093 210304 SKKS093-ADNAxons-3rpm-003 7550 79 1.89 39 18 18 0.82 0.72
adn_axons SKKS094 210217 SKKS094-axons-001 7550 66 1.88 37 17 20 0.81 0.56
adn_axons SKKS094 210217 SKKS094-axons-002 7550 75 1.88 38 17 18 0.82 0.47
adn_axons SKKS094 210304 SKKS094-ADNAxons-3rpm-001 7550 112 1.88 37 18 17 1.53 1.70
adn_axons SKKS094 210304 SKKS094-ADNAxons-3rpm-002 7550 110 1.84 35 18 17 1.55 1.53
adn_axons SKKS094 210304 SKKS094-ADNAxons-3rpm-003 7550 120 1.62 23 13 9 1.64 1.13
adn_axons SKKS095 210217 SKKS095-axons-001 7550 44 1.90 39 20 18 0.47 0.28
adn_axons SKKS095 210217 SKKS095-axons-002 7550 59 1.89 37 15 18 0.41 0.35
soma_dual SKKS080 201123 SKKS080-DualCue-003 1950 112 2.35 12 5 6 0.88 0.56
soma_dual SKKS080 201123 SKKS080-DualCue-004 1950 110 2.55 12 6 6 0.76 0.40
soma_dual SKKS088 210118 SKKS088-DualCue-001 1950 329 2.34 9 6 2 4.92 3.11
soma_dual SKKS088 210118 SKKS088-DualCue-002 1950 240 0.04 4 0 3 4.43 1.93
soma_dual SKKS088 210118 SKKS088-DualCue-003 1950 286 0.03 4 1 2 5.56 3.22
soma_dual SKKS088 210216 SKKS088-DualCue-001 7550 411 1.92 40 18 19 2.92 1.75
soma_dual SKKS088 210216 SKKS088-DualCue-002 7550 289 1.94 39 19 18 2.31 0.86
soma_dual SKKS089 201205 SKKS089-DualCue-003 1950 115 4.09 20 10 10 1.68 1.76
soma_dual SKKS089 201205 SKKS089-DualCue-004 1950 82 4.37 19 9 10 1.04 0.69
soma_dual SKKS089 210118 SKKS089-DualCue-002 1950 307 0.21 6 3 2 5.13 3.17
soma_dual SKKS089 210216 SKKS089-DualCue-001 7550 435 1.94 38 17 19 5.61 3.77
soma_dual SKKS089 210216 SKKS089-DualCue-002 7550 459 1.94 39 18 19 5.06 2.91
soma_dual SKKS089 210312 SKKS089-RSCSoma-3rpm-001 7550 498 1.88 38 17 18 6.81 5.90
soma_dual SKKS089 210312 SKKS089-RSCSoma-3rpm-002 7550 430 1.90 39 19 16 5.65 4.61
soma_dual SKKS089 210312 SKKS089-RSCSoma-3rpm-003 7550 518 1.88 36 18 17 6.56 4.55
soma_dual SKKS090 201205 SKKS090-DualCue-003 1950 172 4.53 22 11 11 3.35 1.75
soma_dual SKKS090 201205 SKKS090-DualCue-004 1950 145 4.66 21 10 11 2.80 1.57
soma_dual SKKS090 210203 SKKS090-DualCue-001 7550 164 1.95 38 19 18 1.92 1.66
soma_dual SKKS090 210216 SKKS090-DualCue-001 7550 368 1.82 37 17 17 2.81 3.59
soma_dual SKKS090 210216 SKKS090-DualCue-002 7550 549 1.84 36 16 18 4.87 4.26
soma_dual SKKS090 210312 SKKS090-RSCSoma-3rpm-001 7550 479 1.90 39 18 18 4.65 5.61
soma_dual SKKS090 210312 SKKS090-RSCSoma-3rpm-002 7550 549 1.92 41 19 18 4.06 5.46
soma_dual SKKS090 210312 SKKS090-RSCSoma-3rpm-003 7550 586 1.92 38 16 18 5.07 4.98
soma_dual SKKS091 210203 SKKS091-DualCue-001 7550 162 1.95 39 19 19 1.80 0.68
soma_dual SKKS091 210203 SKKS091-DualCue-002 7550 158 1.96 38 19 18 1.62 0.79
soma_dual SKKS091 210216 SKKS091-DualCue-001 7550 326 1.99 40 19 19 2.32 1.68
soma_dual SKKS091 210216 SKKS091-DualCue-002 7550 317 1.98 38 19 18 2.50 2.09
soma_dual SKKS091 210311 SKKS091-RSCSoma-3rpm-001 7550 346 1.90 38 18 15 3.83 2.68
soma_dual SKKS091 210311 SKKS091-RSCSoma-3rpm-002 7550 290 1.90 39 18 16 2.86 1.84
soma_dual SKKS091 210311 SKKS091-RSCSoma-3rpm-003 7550 307 1.91 39 18 18 2.87 2.24
soma_dual SKKS100 210216 SKKS100-DualCue-001 7550 496 1.98 40 18 19 5.58 3.93
soma_dual SKKS100 210216 SKKS100-DualCue-002 7550 440 1.97 39 19 18 4.34 2.89
soma_dual SKKS100 210311 SKKS100-RSCSoma-3rpm-001 7550 395 1.91 38 18 16 3.83 4.19
soma_dual SKKS100 210311 SKKS100-RSCSoma-3rpm-002 7550 373 1.89 37 17 19 4.63 3.42
soma_dual SKKS100 210311 SKKS100-RSCSoma-3rpm-003 7550 326 1.89 38 16 18 3.98 2.47
vis_axons SKKS097 210217 SKKS097-axons-001 7550 68 1.73 34 15 16 1.66 0.25
vis_axons SKKS097 210217 SKKS097-axons-002 7550 47 1.87 37 17 17 0.96 0.12
vis_axons SKKS098 210217 SKKS098-axons-001 7550 36 1.92 39 21 18 0.62 0.25
vis_axons SKKS098 210217 SKKS098-axons-002 7550 26 1.93 39 17 18 0.43 0.13
vis_axons SKKS103 210304 SKKS103-VISAxons-3rpm-001 7550 37 0.11 11 8 3 0.76 0.45
vis_axons SKKS103 210304 SKKS103-VISAxons-3rpm-002 7550 28 1.87 37 17 15 0.57 0.49
vis_axons SKKS103 210304 SKKS103-VISAxons-3rpm-003 7550 5 1.89 38 18 17 0.03 0.01
files mice rois
experiment
adn_axons 17 4 1084
soma_dual 35 6 11569
vis_axons 7 3 247
The files do not say when the lights were on. The code hard-codes the segments (Separator.m):
950-frame light and dark halves for 1950-frame recordings, three light/dark repeats of 1250
frames for 7550-frame ones. The summed activity per segment is shown above as a check that is
consistent with, but does not prove, the segmentation. Several 1950-frame recordings of
2021-01-18 stop rotating early and have few trials.
Clock check. As in rotation_remapping.ipynb, files that rotate one way at constant speed
cannot show a time offset between heading and activity; the clock is checked only on the files
that rotate both ways (single_cue.ipynb), and is not established for these files.
Heading-selection test (two-term Gaussian)¶
On all trials, light and dark together (the code's default, claims file), with the centred fit
adopted in single_cue.ipynb. Axons are tested as in axons.ipynb (same seed, so the same
axons pass).
t0 = time.time()
trials_all = [rsc.trial_tuning(r["sp"], r["heading"], r["seg_all"]) for r in recs]
ax_idx = [i for i in range(len(recs)) if "axons" in sub.experiment[i]]
with Pool(4) as pool:
soma_files = [i for i in range(len(recs)) if i not in ax_idx]
sel_s = rsc.selection_many([trials_all[i] for i in soma_files], "double", 306, pool, centre=True)
sel_a = rsc.selection_many([trials_all[i] for i in ax_idx], "double", 301, pool, centre=True)
sel = [None] * len(recs)
for i, s in zip(soma_files, sel_s):
sel[i] = s
for i, s in zip(ax_idx, sel_a):
sel[i] = s
info["selective"] = [s.is_hd.sum() for s in sel]
print(f"selection tests in {time.time() - t0:.0f} s")
print(info.groupby("experiment")[["rois", "selective"]].sum().to_string())
selection tests in 506 s
rois selective
experiment
adn_axons 1084 135
soma_dual 11569 1433
vis_axons 247 106
# Light and dark tuning, aligned with the light-on cross-validated shifts (Methods).
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)
# light-off trials take the shift of the mean light-on curve (all light trials)
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))]
print("Recordings without light-on trials (left out of the light/dark analyses):",
[info.identifier[i] for i in range(len(recs)) if aligned[i][0] is None])
print("Recordings without light-off trials:", [info.identifier[i] for i in range(len(recs)) if aligned[i][0] is not None and np.isnan(aligned[i][1]).all()])
Recordings without light-on trials (left out of the light/dark analyses): ['210118 SKKS088-DualCue-002'] Recordings without light-off trials: []
RSC somata: Gaussian fits and clustering¶
For each soma, light-on and light-off curves aligned with the light-on shifts are rescaled
together to 0 to 1 and each fitted with the two-term Gaussian (code
LightDarkExperiment.calculateTuningFits2); features are a1, a2, c1, c2 of both fits. Cells are
(a) the heading-selective somata, (b) all somata.
soma_idx = [i for i in range(len(recs)) if sub.experiment[i] == "soma_dual" and aligned[i][0] is not None]
def soma_features(mask_fn):
L, D, mice, recid = [], [], [], []
for i in soma_idx:
m = mask_fn(i) & ~np.isnan(aligned[i][1]).any(1)
L.append(aligned[i][0][m]); D.append(aligned[i][1][m])
mice += [sub.subject[i]] * int(m.sum()); recid += [i] * int(m.sum())
L, D = np.vstack(L), np.vstack(D)
both = rsc.rescale(np.hstack([L, D]))
Lr, Dr = both[:, :60], both[:, 60:]
pl, r2l = rsc.fit_gauss(Lr, "double", do_rescale=False)
pd_, r2d = rsc.fit_gauss(Dr, "double", do_rescale=False)
feats = np.column_stack([pl[:, [1, 2, 5, 6]], pd_[:, [1, 2, 5, 6]]])
return dict(L=Lr, D=Dr, Lraw=L, Draw=D, feats=feats, pl=pl, pd=pd_, mice=np.array(mice), rec=np.array(recid), r2l=r2l, r2d=r2d)
t0 = time.time()
soma = {"(a) selective": soma_features(lambda i: sel[i].is_hd.values),
"(b) all": soma_features(lambda i: np.ones(len(sel[i]), bool))}
for k, v in soma.items():
print(f"{k}: {len(v['feats'])} somata from {len(np.unique(v['rec']))} recordings, {len(np.unique(v['mice']))} mice; "
f"median r2 light {np.nanmedian(v['r2l']):.2f}, dark {np.nanmedian(v['r2d']):.2f}")
print(f"fits in {time.time() - t0:.0f} s")
(a) selective: 1419 somata from 34 recordings, 6 mice; median r2 light 0.54, dark 0.26 (b) all: 11329 somata from 34 recordings, 6 mice; median r2 light 0.26, dark 0.14 fits in 6 s
fig5-three-classes¶
Clustering the Gaussian-fit coefficients of RSC somatic tuning curves in light-on and light-off conditions gives an optimum of three functional classes (heading, landmark and alignment cells), though the responses span a continuum.
Gaussian mixture models (full covariance, 20 replicates, up to 1000 iterations, as the reviewer code) for k = 2 to 6, silhouette value of the hard assignments. Criterion: silhouette at k = 3 above k = 2 and k = 4.
gmms = {}
for k_name, v in soma.items():
X = v["feats"]
sil = {}
for k in range(2, 7):
g = GaussianMixture(k, covariance_type="full", n_init=20, max_iter=1000, random_state=304).fit(X)
lab = g.predict(X)
sil[k] = silhouette_score(X, lab) if len(np.unique(lab)) > 1 else np.nan
gmms[(k_name, k)] = g
print(f"{k_name}: silhouette by k:", {k: round(float(s), 4) for k, s in sil.items()},
"| local maximum at 3:", bool(sil[3] > sil[2] and sil[3] > sil[4]))
# variant: the reviewer code's fit on unaligned curves centred on their peak
for k_name, v in soma.items():
both = rsc.rescale(np.hstack([rsc.center(v["Lraw"]), rsc.center(v["Draw"])]))
pl, _ = rsc.fit_gauss(both[:, :60], "double", do_rescale=False)
pd2, _ = rsc.fit_gauss(both[:, 60:], "double", do_rescale=False)
X2 = np.column_stack([pl[:, [1, 2, 5, 6]], pd2[:, [1, 2, 5, 6]]])
sil2 = {}
for k in range(2, 5):
lab = GaussianMixture(k, covariance_type="full", n_init=20, max_iter=1000, random_state=305).fit(X2).predict(X2)
sil2[k] = silhouette_score(X2, lab)
print(f"variant, each curve centred on its own peak, {k_name}:", {k: round(float(s), 4) for k, s in sil2.items()})
(a) selective: silhouette by k: {2: 0.2723, 3: 0.0191, 4: 0.0321, 5: 0.0317, 6: 0.0384} | local maximum at 3: False
(b) all: silhouette by k: {2: 0.1897, 3: 0.1198, 4: -0.0518, 5: 0.0117, 6: -0.0034} | local maximum at 3: False
variant, each curve centred on its own peak, (a) selective: {2: 0.1136, 3: 0.3926, 4: 0.0043}
variant, each curve centred on its own peak, (b) all: {2: 0.2598, 3: -0.0148, 4: 0.0162}
Result: not reproduced. The silhouette value is highest at k = 2 under both readings of the criterion: selective cells give 0.272, 0.019 and 0.032 for k = 2, 3 and 4, and all cells 0.190, 0.120 and -0.052. Neither has a local maximum at 3. The reviewer-code variant, with each curve centred on its own peak, does have one for selective cells (0.114, 0.393 and 0.004), but not for all cells (0.260, -0.015 and 0.016). The criterion decides on readings (a) and (b) and reports the variant only, so the outcome stays not reproduced; the variant's peak at 3 is noted as a reading under which the claim would hold.
fig5-class-proportions¶
12.1% of RSC neurons are tuned in the dual-cue experiment; of these, landmark cells are the plurality (33%), followed by alignment cells (28%) and heading cells (22%).
Three-component GMM on the features of reading (a); outlier rule (a) squared Mahalanobis distance to the nearest component beyond the chi-square 95th percentile (8 df) goes to a null cluster, or (b) none; archetypal light and dark curves from each component's mean coefficients (b1 = pi/2, b2 = 3 pi/2, a0 = 0); each cell to the archetype pair with the smallest summed cosine distance. Classes are named from the archetypes: heading is the one whose light-on archetype is the least bimodal (a2/a1), and of the other two, landmark is the one with the weaker light-off archetype. Shares are given (i) as the mean over recordings, (ii) the mean over mice, (iii) pooled.
def archetype(coef):
a1, a2, c1, c2 = coef
x = rsc.X
return a1 * np.exp(-((x - np.pi / 2) / c1) ** 2) + a2 * np.exp(-((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))
v = soma["(a) selective"]
g = gmms[("(a) selective", 3)]
X = v["feats"]
arch_L = [archetype(mu[:4]) for mu in g.means_]
arch_D = [archetype(mu[4:]) for mu in g.means_]
dist = np.column_stack([cosdist(v["L"], aL) + cosdist(v["D"], aD) for aL, aD in zip(arch_L, arch_D)])
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]))]
align_c = [c for c in rest if c != landmark_c][0]
names = {heading_c: "heading", landmark_c: "landmark", align_c: "alignment"}
print("component means (a1, a2, c1, c2 light; a1, a2, c1, c2 dark):")
for c in range(3):
print(f" {names[c]:9s}", np.round(g.means_[c], 2).tolist(), f"light a2/a1 ratio {bimod[c]:.2f}, dark peak {dark_amp[c]:.2f}")
# outlier rule (a)
md2 = np.min(np.column_stack([np.einsum("ni,ij,nj->n", X - g.means_[c], np.linalg.inv(g.covariances_[c]), X - g.means_[c]) for c in range(3)]), 1)
outlier = md2 > stats.chi2.ppf(0.95, X.shape[1])
labels = {"(a) chi-square outliers to null": np.where(outlier, -1, cls), "(b) no outlier step": cls}
total = {"rec": [], "mice": {}}
n_rois_rec = {i: len(sel[i]) for i in soma_idx}
tuned_rec = np.array([sel[i].is_hd.sum() / len(sel[i]) for i in soma_idx])
mice_rec = np.array([sub.subject[i] for i in soma_idx])
print(f"\ntuned somata: pooled {100 * sum(sel[i].is_hd.sum() for i in soma_idx) / sum(n_rois_rec.values()):.2f}%, "
f"mean over recordings {100 * tuned_rec.mean():.2f}% (s.e.m. {100 * tuned_rec.std(ddof=1) / np.sqrt(len(tuned_rec)):.2f}), "
f"mean over mice {100 * pd.Series(tuned_rec).groupby(mice_rec).mean().mean():.2f}% (paper 12.1 +- 1.1)")
all_soma = [i for i in range(len(recs)) if sub.experiment[i] == "soma_dual"]
print(f"over all {len(all_soma)} soma recordings, including those without light-on trials: "
f"{100 * sum(sel[i].is_hd.sum() for i in all_soma) / sum(len(sel[i]) for i in all_soma):.2f}%")
share_rows = []
for lname, lab in labels.items():
df = pd.DataFrame({"cls": lab, "rec": v["rec"], "mouse": v["mice"]})
kept = df[df.cls >= 0]
for avg in ["(i) recordings", "(ii) mice", "(iii) pooled"]:
if avg == "(iii) pooled":
sh = kept.cls.value_counts(normalize=True)
sem = pd.Series(np.nan, index=sh.index)
else:
key = "rec" if avg.startswith("(i)") else "mouse"
tab = pd.crosstab(df[key], df.cls, normalize="index") # share of all tuned cells of the unit, null included
tab = tab.drop(columns=[-1], errors="ignore")
sh, sem = tab.mean(), tab.std(ddof=1) / np.sqrt(len(tab))
share_rows.append({"outliers": lname, "average": avg, "null": int((lab < 0).sum()),
**{f"{names[c]} %": round(100 * float(sh.get(c, 0)), 1) for c in range(3)},
**{f"{names[c]} sem": round(100 * float(sem.get(c, np.nan)), 1) for c in range(3)}})
shares = pd.DataFrame(share_rows)
print(shares.to_string(index=False))
print("\ncounts (no outlier step):", {names[c]: int((cls == c).sum()) for c in range(3)},
"| paper (from F-test degrees of freedom): heading 332, landmark 522, alignment 409")
print("counts by mouse (no outlier step):")
print(pd.crosstab(v["mice"], pd.Series(cls).map(names).values).to_string())
component means (a1, a2, c1, c2 light; a1, a2, c1, c2 dark):
heading [0.44, 0.3, 0.71, 1.88, 0.37, 0.21, 0.64, 0.88] light a2/a1 ratio 0.70, dark peak 0.37
landmark [0.58, 0.46, 0.79, 0.7, 0.05, 0.02, 0.65, 0.66] light a2/a1 ratio 0.79, dark peak 0.05
alignment [0.36, 0.28, 2.03, 0.91, 0.22, 0.26, 1.84, 1.87] light a2/a1 ratio 0.80, dark peak 0.26
tuned somata: pooled 12.53%, mean over recordings 12.10% (s.e.m. 1.44), mean over mice 11.51% (paper 12.1 +- 1.1)
over all 35 soma recordings, including those without light-on trials: 12.39%
outliers average null heading % landmark % alignment % heading sem landmark sem alignment sem
(a) chi-square outliers to null (i) recordings 58 13.3 7.4 75.3 2.3 1.7 3.0
(a) chi-square outliers to null (ii) mice 58 15.3 9.5 72.2 4.2 4.5 7.5
(a) chi-square outliers to null (iii) pooled 58 14.5 4.8 80.6 NaN NaN NaN
(b) no outlier step (i) recordings 0 13.5 7.6 78.9 2.3 1.7 3.3
(b) no outlier step (ii) mice 0 15.5 9.6 74.9 4.2 4.5 8.1
(b) no outlier step (iii) pooled 0 14.2 4.7 81.0 NaN NaN NaN
counts (no outlier step): {'heading': 202, 'landmark': 67, 'alignment': 1150} | paper (from F-test degrees of freedom): heading 332, landmark 522, alignment 409
counts by mouse (no outlier step):
col_0 alignment heading landmark
row_0
SKKS080 7 6 6
SKKS088 98 22 13
SKKS089 416 121 15
SKKS090 181 17 9
SKKS091 203 27 12
SKKS100 245 9 12
Result: not reproduced. 12.10% of somata are tuned (mean over recordings, s.e.m. 1.44; pooled 12.53%), matching 12.1 +- 1.1%. The three classes do not match. Over recordings, alignment-like cells make up 78.9%, heading 13.5% and landmark 7.6% (no outlier step); with the outlier step the shares are 75.3, 13.3 and 7.4%. Every reading gives alignment > heading > landmark, not landmark > alignment > heading. The counts are 1150, 202 and 67, against 409, 332 and 522 implied by the paper's degrees of freedom. All three light-on archetypes are bimodal (a2/a1 0.70 to 0.80), so the class called heading here is only the least bimodal of three bimodal classes.
fig5e-alignment-peaks-coregistered¶
In alignment cells the single light-off peak is always at the location of one of the two light-on peaks, so landmark and heading responses are coregistered.
Per alignment cell (no outlier step), offset of the light-off peak bin from the nearer of the two light-on peaks (the peak bin, and the peak bin of the other half circle), V-test against 0.
Taking the nearer of two peaks 180 degrees apart folds the offsets into +-90 degrees, so the V-test against 0 is positive even for light-off peaks unrelated to the light-on curve. The V of each class is therefore compared with two nulls: light-off peaks paired at random with other cells' light-on curves within the class, and light-off peaks drawn uniformly.
def peak_offsets(L, D):
out = []
for l, d in zip(L, D):
p1 = int(np.argmax(l))
other = [(p1 + 30 + k) % 60 for k in range(-15, 16)]
p2 = other[int(np.argmax(l[other]))]
pd_ = int(np.argmax(d))
o = [rsc.wrap180((pd_ - p) * 6) for p in (p1, p2)]
out.append(o[int(np.argmin(np.abs(o)))])
return np.array(out)
for cname in ["alignment", "heading", "landmark"]:
c = [k for k, n in names.items() if n == cname][0]
m = cls == c
off = peak_offsets(v["L"][m], v["D"][m])
V, p, n = rsc.vtest(off)
print(f"{cname}: n={n}, V={V:.1f}, p={p:.2g}, share within 18 deg {np.mean(np.abs(off) <= 18):.3f}, "
f"median |offset| {np.median(np.abs(off)):.0f} deg; under uniform peaks the share within 18 deg of the nearer of two peaks is about 0.21")
# Null for V: taking the nearer of two peaks 180 degrees apart folds every offset into +-90, which
# by itself gives V > 0. (1) light-off peaks of the class paired with the light-on curves of other
# cells of the class at random; (2) light-off peaks drawn uniformly over the 60 bins.
def nearer(p1, p2, pdk):
o1, o2 = rsc.wrap180((pdk - p1) * 6), rsc.wrap180((pdk - p2) * 6)
return np.where(np.abs(o1) <= np.abs(o2), o1, o2)
rng_e = np.random.default_rng(307)
for cname in ["alignment", "heading", "landmark"]:
c = [k for k, n in names.items() if n == cname][0]
m = cls == c
Lm, Dm = v["L"][m], v["D"][m]
p1 = Lm.argmax(1)
other = (p1[:, None] + 30 + np.arange(-15, 16)[None]) % 60
p2 = np.take_along_axis(other, np.take_along_axis(Lm, other, 1).argmax(1)[:, None], 1)[:, 0]
pdk = Dm.argmax(1)
off = nearer(p1, p2, pdk)
V, p, n = rsc.vtest(off)
null_perm = np.array([rsc.vtest(nearer(p1, p2, pdk[rng_e.permutation(n)]))[0] for _ in range(1000)])
null_unif = np.array([rsc.vtest(nearer(p1, p2, rng_e.integers(0, 60, n)))[0] for _ in range(1000)])
w18_perm = np.mean([np.mean(np.abs(nearer(p1, p2, pdk[rng_e.permutation(n)])) <= 18) for _ in range(200)])
print(f"{cname}: V = {V:.1f} (V/n {V / n:.3f}); random pairing V mean {null_perm.mean():.1f}, 95th percentile "
f"{np.percentile(null_perm, 95):.1f}, share at or above the real V {np.mean(null_perm >= V):.3f}; uniform peaks V mean "
f"{null_unif.mean():.1f} (V/n {null_unif.mean() / n:.3f}), 95th percentile {np.percentile(null_unif, 95):.1f}; "
f"share within 18 deg real {np.mean(np.abs(off) <= 18):.3f}, random pairing {w18_perm:.3f}")
alignment: n=1150, V=743.5, p=0, share within 18 deg 0.292, median |offset| 36 deg; under uniform peaks the share within 18 deg of the nearer of two peaks is about 0.21 heading: n=202, V=138.7, p=0, share within 18 deg 0.475, median |offset| 24 deg; under uniform peaks the share within 18 deg of the nearer of two peaks is about 0.21 landmark: n=67, V=58.2, p=0, share within 18 deg 0.537, median |offset| 18 deg; under uniform peaks the share within 18 deg of the nearer of two peaks is about 0.21
alignment: V = 743.5 (V/n 0.646); random pairing V mean 707.2, 95th percentile 728.4, share at or above the real V 0.002; uniform peaks V mean 669.5 (V/n 0.582), 95th percentile 692.3; share within 18 deg real 0.292, random pairing 0.262
heading: V = 138.7 (V/n 0.686); random pairing V mean 140.7, 95th percentile 147.8, share at or above the real V 0.677; uniform peaks V mean 116.6 (V/n 0.577), 95th percentile 125.8; share within 18 deg real 0.475, random pairing 0.433
landmark: V = 58.2 (V/n 0.869); random pairing V mean 57.1, 95th percentile 59.3, share at or above the real V 0.230; uniform peaks V mean 40.2 (V/n 0.600), 95th percentile 44.9; share within 18 deg real 0.537, random pairing 0.532
Result: partially reproduced, and the classes are not the paper's (see fig5-class-proportions). For the 1150 cells of the alignment-like class, the V-test of the criterion passes (V = 743.5, p < 1e-300), but most of that V comes from taking the nearer of two peaks: light-off peaks drawn uniformly give V = 669.5 on average (V/n 0.582 against 0.646), and light-off peaks paired at random with other alignment-like cells give 707.2 (95th percentile 728.4). The real V exceeds that null (share of random pairings at or above it 0.002), so there is a small coregistration. Only 29.2% of the cells have the light-off peak within 18 degrees of a light-on peak, against 26.2% under random pairing, with a median offset of 36 degrees, so "always" is not supported. In the heading-like and landmark-like classes the real V does not exceed random pairing (shares at or above it 0.677 and 0.230).
fig5f-heading-flip-unchanged and fig5gh-landmark-alignment-flip-drop¶
Heading cells show no significant difference in flip score between light-on and light-off conditions. Landmark and alignment cells have significantly lower flip scores in light-off than in light-on.
Flip score of each cell's aligned light-on and light-off curves; linear mixed model flip ~ light
- (1 | mouse) (statsmodels MixedLM, REML), F = z squared of the light effect with the paper's residual df 2n - 2; a model with a random intercept per cell as well is a variant.
flip_rows = []
for lname, lab in labels.items():
for cname in ["heading", "landmark", "alignment"]:
c = [k for k, n in names.items() if n == cname][0]
m = lab == c
fl, fd = rsc.flip_score(v["L"][m]), rsc.flip_score(v["D"][m])
df = pd.DataFrame({"flip": np.r_[fl, fd], "light": np.r_[np.ones(m.sum()), np.zeros(m.sum())],
"mouse": np.r_[v["mice"][m], v["mice"][m]], "cell": np.r_[np.arange(m.sum()), np.arange(m.sum())]}).dropna()
fit = smf.mixedlm("flip ~ light", df, groups=df["mouse"]).fit(reml=True)
z, p = fit.tvalues["light"], fit.pvalues["light"]
fit2 = smf.mixedlm("flip ~ light", df, groups=df["mouse"], vc_formula={"cell": "0 + C(cell)"}).fit(reml=True)
flip_rows.append(dict(outliers=lname, cls=cname, n=int(m.sum()), mean_on=round(float(np.nanmean(fl)), 3),
mean_off=round(float(np.nanmean(fd)), 3), diff_off_minus_on=round(float(fit.params["light"] * -1), 3),
F=round(float(z ** 2), 1), df=f"1,{2 * int(m.sum()) - 2}", p=float(f"{p:.3g}"),
p_cell_variant=float(f"{fit2.pvalues['light']:.3g}"),
paired_t_p=float(f"{stats.ttest_rel(fl, fd, nan_policy='omit').pvalue:.3g}")))
flips = pd.DataFrame(flip_rows)
print(flips.to_string(index=False))
print("paper: heading F(1,662) = 0.18, p = 0.67; landmark F(1,1042) = 1062.5, p = 3.1e-161; alignment F(1,816) = 209.5, p = 2.0e-42")
outliers cls n mean_on mean_off diff_off_minus_on F df p p_cell_variant paired_t_p
(a) chi-square outliers to null heading 198 0.381 0.281 -0.102 2.2 1,394 1.380000e-01 1.140000e-01 1.480000e-01
(a) chi-square outliers to null landmark 66 0.490 -0.238 -0.728 45.0 1,130 2.010000e-11 1.500000e-12 1.280000e-09
(a) chi-square outliers to null alignment 1097 0.381 0.014 -0.366 140.8 1,2192 1.750000e-32 1.550000e-32 1.080000e-30
(b) no outlier step heading 202 0.378 0.279 -0.100 2.2 1,402 1.410000e-01 1.160000e-01 1.500000e-01
(b) no outlier step landmark 67 0.478 -0.239 -0.717 44.3 1,132 2.800000e-11 2.180000e-12 1.480000e-09
(b) no outlier step alignment 1150 0.365 0.014 -0.351 135.4 1,2298 2.680000e-31 2.780000e-31 1.190000e-29
paper: heading F(1,662) = 0.18, p = 0.67; landmark F(1,1042) = 1062.5, p = 3.1e-161; alignment F(1,816) = 209.5, p = 2.0e-42
Result: fig5f-heading-flip-unchanged and fig5gh-landmark-alignment-flip-drop partially reproduced: each criterion is met, on classes that are not the paper's. The flip score of the 202 heading-like cells does not differ significantly between light-on and light-off: 0.378 against 0.279, F(1,402) = 2.2, p = 0.141; with the outlier step p = 0.138; the variant with a random intercept per cell gives p = 0.116. The class is not the paper's heading class: it has 202 cells where the paper's degrees of freedom imply 332, and its light-on archetype is bimodal (a2/a1 0.70). Flip scores fall in light-off for landmark-like cells (0.478 to -0.239, F(1,132) = 44.3, p = 2.8e-11, n = 67) and alignment-like cells (0.365 to 0.014, F(1,2298) = 135.4, p = 2.7e-31, n = 1150).
fig5i-axons-in-clusters¶
Projected into the soma clusters' discriminant space, ADN axons fall mostly within the heading cluster and visual axons mostly within the landmark cluster, with some overlap with the alignment cluster.
LDA trained on the soma features and classes (no outlier step); axon features from the same fits (light-on and light-off curves aligned with the light-on shifts, rescaled together), axons as in fig4 readings (a) and (b); 95% ellipse of each soma class: squared Mahalanobis distance in the first two discriminants below the chi-square 95th percentile with 2 df.
lda = LinearDiscriminantAnalysis(n_components=2).fit(X, cls)
Z = lda.transform(X)
ell = {}
for c in range(3):
zc = Z[cls == c]
ell[c] = (zc.mean(0), np.linalg.inv(np.cov(zc.T)))
thr = stats.chi2.ppf(0.95, 2)
def inside(z, c):
mu, icov = ell[c]
return np.einsum("ni,ij,nj->n", z - mu, icov, z - mu) <= thr
def axon_sets(selection_by_rec):
out = {}
for kind in ["adn_axons", "vis_axons"]:
la, da = [], []
for i in ax_idx:
if sub.experiment[i] != kind or aligned[i][0] is None:
continue
m = selection_by_rec(i)
la.append(aligned[i][0][m]); da.append(aligned[i][1][m])
out[kind] = (np.vstack(la), np.vstack(da))
return out
readings = {"(a) selective, centred": axon_sets(lambda i: sel[i].is_hd.values),
"(b) all axons": axon_sets(lambda i: np.ones(len(sel[i]), bool))}
for name in ["(a) selective, centred", "(b) all axons"]:
s = readings[name]
for kind in ["adn_axons", "vis_axons"]:
la, da = s[kind][0], s[kind][1]
ok = ~np.isnan(da).any(1)
both = rsc.rescale(np.hstack([la[ok], da[ok]]))
pl, _ = rsc.fit_gauss(both[:, :60], "double", do_rescale=False)
pdk, _ = rsc.fit_gauss(both[:, 60:], "double", do_rescale=False)
za = lda.transform(np.column_stack([pl[:, [1, 2, 5, 6]], pdk[:, [1, 2, 5, 6]]]))
shares_in = {names[c]: round(float(inside(za, c).mean()), 3) for c in range(3)}
nearest = pd.Series(lda.predict(np.column_stack([pl[:, [1, 2, 5, 6]], pdk[:, [1, 2, 5, 6]]]))).map(names).value_counts(normalize=True).round(3).to_dict()
print(f"{name}, {kind}: n={len(za)}; share inside each soma ellipse {shares_in}; LDA class assigned {nearest}")
fig, ax = plt.subplots(figsize=(4, 4))
for c in range(3):
ax.scatter(*Z[cls == c].T, s=3, alpha=0.3, label=names[c])
ax.set_xlabel("LD1"); ax.set_ylabel("LD2"); ax.legend(); plt.tight_layout(); plt.show()
print("LD1 range", np.round([Z[:, 0].min(), Z[:, 0].max()], 2), "LD2 range", np.round([Z[:, 1].min(), Z[:, 1].max()], 2))
(a) selective, centred, adn_axons: n=135; share inside each soma ellipse {'heading': 0.881, 'landmark': 0.444, 'alignment': 0.859}; LDA class assigned {'alignment': 0.852, 'heading': 0.148}
(a) selective, centred, vis_axons: n=106; share inside each soma ellipse {'heading': 0.925, 'landmark': 0.83, 'alignment': 0.962}; LDA class assigned {'alignment': 0.972, 'heading': 0.028}
(b) all axons, adn_axons: n=1084; share inside each soma ellipse {'heading': 0.724, 'landmark': 0.28, 'alignment': 0.843}; LDA class assigned {'alignment': 0.919, 'heading': 0.081}
(b) all axons, vis_axons: n=247; share inside each soma ellipse {'heading': 0.81, 'landmark': 0.555, 'alignment': 0.927}; LDA class assigned {'alignment': 0.98, 'heading': 0.02}
LD1 range [-4.31 4.2 ] LD2 range [-2.06 3.99]
Result: partially reproduced: the criterion is met, but the ellipses do not separate the axon groups. Under reading (a), 88.1% of ADN axons lie inside the heading-class ellipse and 83.0% of visual axons inside the landmark-class ellipse. Under (b), 72.4% and 55.5%. The soma ellipses overlap heavily: the alignment ellipse holds 84.3 to 96.2% of both axon groups, the heading ellipse holds 81.0 to 92.5% of visual axons (more than the landmark ellipse does), and the LDA assigns 85.2 to 98.0% of the axons to the alignment class. The classes are those of fig5-class-proportions, not the paper's.
Summary¶
| Claim | Outcome |
|---|---|
| fig5-three-classes | not reproduced: silhouette highest at k = 2 under both readings (the centred variant peaks at 3 for selective cells) |
| fig5-class-proportions | not reproduced: 12.1% tuned matches; class shares do not |
| fig5e-alignment-peaks-coregistered | partially reproduced: V above a random-pairing null by a small margin; 29.2% within 18 degrees (26.2% at random); classes not the paper's |
| fig5f-heading-flip-unchanged | partially reproduced: p = 0.141, on a heading-like class that is not the paper's |
| fig5gh-landmark-alignment-flip-drop | partially reproduced: both drops significant, on classes that are not the paper's |
| fig5i-axons-in-clusters | partially reproduced: criterion met; the ellipses overlap and do not separate the axon groups |
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: inferred spikes from raw fluorescence without neuropil correction; selection with the centred fit; Gaussian fits reimplemented in numpy; light-off curves take the shift of the mean light-on curve (leaveOneOutAligner is not public); the GMM, LDA and mixed models are scikit-learn and statsmodels in place of MATLAB; the class names come from a rule on the archetypes stated above. The classes differ from the paper's, so the claims that use them test other cell groups.
What was found about the dataset: 35 dual-cue soma recordings where the paper reports 36; no light-on and light-off labels; no matching of ROIs between the standard and the visual-cue-only recordings of Figure 5K.