More population decoders for Tseng et al. (2022), dandiset 000579¶

Paper. Tseng, Chettih, Arlt, Barroso-Luque & Harvey (2022), Neuron 110, 2484-2502, doi:10.1016/j.neuron.2022.05.012.

Dataset. DANDI:000579, version 0.230728.1727, streamed (see reproduce_tseng_2022.ipynb for why this version).

Claims. claims/000579/10.1016-j.neuron.2022.05.012.yaml.

Why this notebook. The main notebook fits the population decoders of the Methods (Trial phase-specific logistic decoders, helper decode_session) on a random subsample of 18 sessions, as fitting them all does not fit its 30-minute limit. The second verification round found two places where the subsample leaves a claim's outcome open:

  • fig3q-cue-decoding: on the subsample PM is above V1 by 0.001 to 0.027 in 10 of 12 readings, on 12 to 14 PM decoders, so V1 against PM is not decided.
  • fig8cd-position-specific-decoders: the Figure 8C legend uses only sessions with the cue offset at 0.76 of the maze, and the 6 such sessions in the pool were not in the subsample.

This notebook fits the decoders of every session in reading b's pool (multi-plane layer 2/3 sessions above 0.65 correct over all trials, the main notebook's adopted selection) whose imaged field is V1 or PM, and of the 6 sessions with the cue offset at 0.76. The decoders, readings and statistics are exactly those of the main notebook (same helper, same seeds inside decode_session). Its results bear only on these two claims. The fig3q result was proposed as a change of outcome and the third verification round agreed with it; the fig8cd result is a check of an assumption.

Claim Tested here
fig3q-cue-decoding V1 against PM, all decoders of the V1 and PM fields
fig8cd-position-specific-decoders cue and dynamic choice, the 6 sessions with the cue offset at 0.76
every other claim in reproduce_tseng_2022.ipynb
In [1]:
import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
import time
import warnings
warnings.filterwarnings("ignore")   # empty-slice warnings; dropped fits are counted in the cells
import numpy as np
import pandas as pd
import tseng
from tseng import AREAS, hboot, p_two, p_unpaired

T0 = time.time()
print("dandiset", tseng.DANDISET, "version", tseng.VERSION)

def pb(boot):
    p = p_two(boot)
    return f"< {1 / len(boot):g}" if p == 0 else f"= {p:.3g}"

def pu(b1, b2):
    p = p_unpaired(b1, b2)
    return f"< {1 / len(b1):g}" if p == 0 else f"= {p:.3g}"

dandiset 000579 version 0.230728.1727

Data check¶

In [2]:
assets = tseng.list_assets()
cand = assets[(assets.kind == "session") & (assets.layer == "L23") & (assets.imaging == "multi")]
S = tseng.summaries(list(cand.asset_id), workers=12)
S = S.merge(cand[["asset_id", "subject", "fov_area", "date"]], on="asset_id").rename(columns={"subject": "mouse", "fov_area": "fov"})
print(f"multi-plane layer 2/3 session files: {len(S)}; mice: {sorted(S.mouse.unique(), key=int)}")
print("fraction correct over all trials:", S.fc_all.min().round(3), "to", S.fc_all.max().round(3),
      "| without visually guided trials:", S.fc.min().round(3), "to", S.fc.max().round(3))
print("cue offset (fraction of the maze):", S.cue_off.round(2).value_counts().sort_index().to_dict())
print("stem as a fraction of the maze:", S.stem_frac.round(3).value_counts().to_dict())
POOL = S[S.fc_all > 0.65]
print(f"pool (reading b, above 0.65 over all trials): {len(POOL)} sessions, {POOL.n_neurons.sum():,} neurons")
print(pd.crosstab(POOL.mouse, POOL.fov, margins=True).to_string())
off76 = POOL[POOL.cue_off.round(2) == 0.76]
print(f"pool sessions with the cue offset at 0.76: {len(off76)}, by field {off76.fov.value_counts().to_dict()}, mice {sorted(off76.mouse.unique(), key=int)}")

multi-plane layer 2/3 session files: 172; mice: ['3', '4', '5', '6', '7', '8', '9', '10']
fraction correct over all trials: 0.5 to 0.894 | without visually guided trials: 0.471 to 0.882
cue offset (fraction of the maze): {0.76: 8, 0.94: 1, 0.95: 4, 0.96: 103, 1.0: 56}
stem as a fraction of the maze: {0.8: 168, 0.75: 4}
pool (reading b, above 0.65 over all trials): 142 sessions, 95,988 neurons
fov    AM  MM  PM  RSC  V1  visA  All
mouse                                
10      2   3   5    2   4     3   19
3       1   0   2    4   1     1    9
4       2   4   0    2   3     2   13
5       2   3   3    2   3     2   15
6       3   1   2    9   6     3   24
7       2   3   0    3   2     5   15
8       6   2   2    3   3    11   27
9       3   2   3    3   2     7   20
All    21  18  17   28  24    34  142
pool sessions with the cue offset at 0.76: 6, by field {'RSC': 2, 'visA': 2, 'PM': 1, 'MM': 1}, mice ['3', '5', '8']

fig3q-cue-decoding: V1 against PM on every decoder of the V1 and PM fields¶

Cue decoded from populations of about 100 nearby neurons follows the same gradient as single-neuron cue encoding.

The criterion asks that under every reading V1 has the highest area mean and A the lowest, with V1 above A at p < 0.05. On the main notebook's subsample, A is the lowest area and V1 is above A in all 12 readings, and every area other than PM is below V1; only PM is above V1, in 10 readings. Here every pool session whose field is V1 or PM is fitted under the same three grouping and fold readings, and the V1 and PM decoders (area of the neuron nearest the centroid, as in the main notebook) are compared under the 12 readings (2 groupings with random-frame folds and k-means with trial folds; 2 penalties; stem bins or all bins). Hierarchical bootstrap over mice then decoders, 1,000 samples, V1 against PM as two unpaired samples, two-tailed. Decoders from sessions in other fields are not fitted here, so the V1 and PM means are those of the V1 and PM fields.

In [3]:
dsess = POOL[POOL.fov.isin(["V1", "PM"])]
PAIRS = [("kmeans", "frames"), ("tiles", "frames"), ("kmeans", "trials")]
t = time.time()
dec_res = tseng.decode_many(list(dsess.asset_id), PAIRS, workers=os.cpu_count())
print(f"{len(dsess)} sessions ({dsess.n_neurons.sum():,} neurons) fit in {time.time() - t:.0f} s on {os.cpu_count()} processes")
dfail = {a: r for a, r in dec_res.items() if isinstance(r, Exception)}
print("sessions that failed:", len(dfail), [repr(e)[:200] for e in dfail.values()])
mouse_of = S.set_index("asset_id").mouse
DEC = pd.DataFrame([r for v in dec_res.values() if not isinstance(v, Exception) for r in v if r["lam"] != "cross"])
DEC["mouse"] = DEC.asset_id.map(mouse_of)
NB = tseng.DEC_NBINS
STEM_BINS = [b for b in range(NB) if (b + 0.5) / NB < S.stem_frac.iloc[0]]
DEC["cue_stem"] = DEC[[f"cue_ll_{b}" for b in STEM_BINS]].mean(axis=1)
DEC["cue_all"] = DEC[[f"cue_ll_{b}" for b in range(NB)]].mean(axis=1)
print("stem bins:", STEM_BINS, "| cue bin fits with too few frames or one class:", int(DEC[[f'cue_ll_{b}' for b in range(NB)]].isna().sum().sum()))
for (g_, f_), x in DEC[DEC.lam == "min"].groupby(["grouping", "folds"]):
    print(f"{g_} groups, {f_} folds: {len(x)} decoders, neurons per decoder {x.n_neurons.min()} to {x.n_neurons.max()} "
          f"(median {x.n_neurons.median():.0f}); by area {x.area.value_counts().reindex(AREAS).fillna(0).astype(int).to_dict()}")
print("paper (Figure 8D): V1 151, PM 106 decoders")

41 sessions (23,375 neurons) fit in 736 s on 4 processes
sessions that failed: 0 []
stem bins: [0, 1, 2, 3, 4, 5, 6, 7, 8] | cue bin fits with too few frames or one class: 0
kmeans groups, frames folds: 234 decoders, neurons per decoder 41 to 219 (median 99); by area {'V1': 149, 'PM': 71, 'AM': 2, 'MM': 12, 'RSC': 0, 'A': 0}
kmeans groups, trials folds: 234 decoders, neurons per decoder 41 to 219 (median 99); by area {'V1': 149, 'PM': 71, 'AM': 2, 'MM': 12, 'RSC': 0, 'A': 0}
tiles groups, frames folds: 234 decoders, neurons per decoder 84 to 114 (median 99); by area {'V1': 142, 'PM': 80, 'AM': 3, 'MM': 9, 'RSC': 0, 'A': 0}
paper (Figure 8D): V1 151, PM 106 decoders
In [4]:
rng = np.random.default_rng(211)
rows = []
for (g_, f_, lam_), sub in DEC.groupby(["grouping", "folds", "lam"]):
    for bname, col in (("stem bins", "cue_stem"), ("all bins", "cue_all")):
        B = {a: hboot([x[col].values for _, x in sub[sub.area == a].groupby("mouse")], 1000, rng) for a in ("V1", "PM")}
        rows.append(dict(groups=g_, folds=f_, penalty=lam_, bins=bname,
                         n_V1=int((sub.area == "V1").sum()), n_PM=int((sub.area == "PM").sum()),
                         V1=round(B["V1"].mean(), 3), PM=round(B["PM"].mean(), 3),
                         V1_minus_PM=round(B["V1"].mean() - B["PM"].mean(), 3), sd_diff=round(np.sqrt(B["V1"].var() + B["PM"].var()), 3),
                         p_V1_PM=pu(B["V1"], B["PM"])))
F3 = pd.DataFrame(rows)
print(F3.to_string(index=False))
print("readings with V1 above PM:", int((F3.V1_minus_PM > 0).sum()), "of", len(F3),
      "| V1 - PM from", F3.V1_minus_PM.min(), "to", F3.V1_minus_PM.max())
print("decoder values (stem bins) range:", round(DEC.cue_stem.min(), 3), "to", round(DEC.cue_stem.max(), 3), "; chance -1")
k0 = DEC[(DEC.grouping == "kmeans") & (DEC.folds == "frames") & (DEC.lam == "min") & DEC.area.isin(["V1", "PM"])]
print("per mouse, k-means / frames / minimum penalty, stem bins (mean, decoders):")
print(k0.groupby(["mouse", "area"]).cue_stem.agg(["mean", "size"]).round(3).unstack("area").to_string())
print("by field, the same reading (mean of V1 and PM decoders):")
k0f = k0.assign(fov=k0.asset_id.map(S.set_index("asset_id").fov))
print(k0f.groupby(["fov", "area"]).cue_stem.agg(["mean", "size"]).round(3).to_string())

groups  folds penalty      bins  n_V1  n_PM     V1     PM  V1_minus_PM  sd_diff  p_V1_PM
kmeans frames     1se stem bins   149    71 -0.156 -0.194        0.038    0.037  = 0.305
kmeans frames     1se  all bins   149    71 -0.219 -0.249        0.030    0.039  = 0.445
kmeans frames     min stem bins   149    71 -0.152 -0.187        0.036    0.034  = 0.287
kmeans frames     min  all bins   149    71 -0.212 -0.239        0.027    0.037  = 0.474
kmeans trials     1se stem bins   149    71 -0.177 -0.219        0.042    0.036  = 0.233
kmeans trials     1se  all bins   149    71 -0.246 -0.282        0.036    0.038  = 0.347
kmeans trials     min stem bins   149    71 -0.168 -0.207        0.039    0.037  = 0.277
kmeans trials     min  all bins   149    71 -0.236 -0.268        0.032    0.037    = 0.4
 tiles frames     1se stem bins   142    80 -0.145 -0.207        0.062    0.032 = 0.0419
 tiles frames     1se  all bins   142    80 -0.207 -0.260        0.053    0.034  = 0.105
 tiles frames     min stem bins   142    80 -0.139 -0.198        0.059    0.030 = 0.0387
 tiles frames     min  all bins   142    80 -0.200 -0.251        0.052    0.032 = 0.0977
readings with V1 above PM: 12 of 12 | V1 - PM from 0.027 to 0.062
decoder values (stem bins) range: -0.689 to -0.033 ; chance -1
per mouse, k-means / frames / minimum penalty, stem bins (mean, decoders):
        mean         size      
area      PM     V1    PM    V1
mouse                          
10    -0.176 -0.141  13.0  16.0
3     -0.178 -0.166  13.0   8.0
4        NaN -0.115   NaN  25.0
5     -0.126 -0.109  17.0  23.0
6     -0.397 -0.238   3.0  28.0
7        NaN -0.155   NaN  12.0
8     -0.280 -0.158  13.0  20.0
9     -0.129 -0.111  12.0  17.0
by field, the same reading (mean of V1 and PM decoders):
           mean  size
fov area             
PM  PM   -0.184    59
    V1   -0.145    20
V1  PM   -0.189    12
    V1   -0.152   129

Result: V1 is above PM under all 12 readings, so fig3q is reproduced (agreed in the third verification round). On every decoder of the 41 pool sessions with a V1 or PM field (23,375 neurons; 149 V1 and 71 PM decoders with k-means groups, 142 and 80 with tiles; the paper has 151 and 106), V1 minus PM is 0.027 to 0.062 in log2 likelihood. The difference is not significant in 10 readings (p = 0.098 to 0.474); it is p = 0.039 and 0.042 for tiles over the stem bins. For k-means groups, frame folds, the minimum penalty and stem bins (the only reading printed per mouse and per field), V1 is above PM in each of the 6 mice with both areas (mice 4 and 7 have no PM decoders; mouse 3 is close) and within each field (V1 field: V1 -0.152, PM -0.189; PM field: V1 -0.145, PM -0.184). The subsample's PM lead came from few decoders: its PM means (-0.153 to -0.263) are above the PM means here (-0.187 to -0.282). In the main notebook AM, MM, RSC and A are at or below -0.276 under every reading, and V1 here is -0.139 to -0.246, so with V1 and PM taken from this notebook V1 is the highest area under every reading; A is the lowest and below V1 at p < 0.001 under every reading (main notebook). The criterion ranks the area means and tests only V1 against A, so this is what it calls reproduced. The comparison takes the other four areas from the main notebook's subsample, where the highest of them (AM) is below this notebook's V1 by 0.118 to 0.157 under the 12 readings; only the V1 against PM ranking, the one the subsample left open, rests on every V1 and PM field decoder. Decoders labelled V1 or PM in sessions of other fields are not fitted here, so PM has 71 to 80 decoders against the paper's 106; in the subsample such PM decoders do not average higher (main notebook).

fig8cd-position-specific-decoders: the sessions with the cue offset at 0.76¶

Cue and choice decoders trained at one maze position lose accuracy when tested at other positions.

The Figure 8C legend: "We obtained 698 decoders for Figure 8C ... using only sessions with cue offset at 0.76 of maze length". The pool has 6 such sessions (printed in the data check). Their decoders (k-means groups, trial folds, minimum-CV penalty), each trained in one position bin and applied to every bin of the held-out trials, as in the main notebook: cue over the bins with the cue present (centre before 0.76 of the maze), log2 likelihood; dynamic choice over all bins, Spearman correlation. Per decoder, the mean over training bins at each train-test distance; the test is the mean at the largest distance minus that at zero distance, per area with decoders, hierarchical bootstrap over mice then decoders (1,000 samples), two-tailed. No V1 or AM field is among the 6 sessions, so areas with no decoders are said to have none.

In [5]:
t = time.time()
x_res = tseng.decode_many(list(off76.asset_id), [("kmeans", "trials")], workers=os.cpu_count())
print(f"{len(off76)} sessions fit in {time.time() - t:.0f} s")
xfail = {a: r for a, r in x_res.items() if isinstance(r, Exception)}
print("sessions that failed:", len(xfail), [repr(e)[:200] for e in xfail.values()])
XP = [dict(r, mouse=mouse_of[r["asset_id"]]) for v in x_res.values() if not isinstance(v, Exception) for r in v if r["lam"] == "cross"]
rng = np.random.default_rng(214)
for v, label in (("cue", "cue, log2 likelihood"), ("dyn", "dynamic choice, Spearman")):
    rows = []
    for r in XP:
        Mx = r[f"x_{v}"]
        dist = np.abs(np.subtract.outer(np.arange(NB), np.arange(NB)))
        prof = np.array([np.nanmean(Mx[(dist == d_) & ~np.isnan(Mx)]) if ((dist == d_) & ~np.isnan(Mx)).any() else np.nan
                         for d_ in range(NB)])
        dmax = int(np.where(~np.isnan(prof))[0].max())
        rows.append(dict(mouse=r["mouse"], area=r["area"], cue_off=r["cue_off"], d0=prof[0], dmax=prof[dmax], max_dist=dmax,
                         **{f"dist{d_}": prof[d_] for d_ in range(NB)}))
    R8 = pd.DataFrame(rows)
    print(f"{label}: {len(R8)} decoders from {R8.mouse.nunique()} mice; cue offsets {R8.cue_off.round(2).value_counts().to_dict()}; "
          f"largest distance {R8.max_dist.min()} to {R8.max_dist.max()} bins")
    print("  mean performance by train-test distance (bins):", R8[[f"dist{d_}" for d_ in range(NB)]].mean().round(3).dropna().to_dict())
    for a in AREAS:
        sub = R8[R8.area == a]
        if len(sub) == 0:
            print(f"  {a}: no decoders")
            continue
        bdiff = hboot([(x.dmax - x.d0).values for _, x in sub.groupby("mouse")], 1000, rng)
        print(f"  {a}: {len(sub)} decoders, {sub.mouse.nunique()} mice; zero distance {sub.d0.mean():.3f}, largest distance {sub.dmax.mean():.3f}; "
              f"largest - zero {bdiff.mean():.3f} ± {bdiff.std():.3f}, p {pb(bdiff)}; decoders with a drop {(sub.dmax < sub.d0).mean():.3f}")

6 sessions fit in 110 s
sessions that failed: 0 []
cue, log2 likelihood: 53 decoders from 3 mice; cue offsets {0.76: 53}; largest distance 7 to 7 bins
  mean performance by train-test distance (bins): {'dist0': -0.451, 'dist1': -0.541, 'dist2': -0.763, 'dist3': -0.942, 'dist4': -1.099, 'dist5': -1.24, 'dist6': -1.31, 'dist7': -1.251}
  V1: 2 decoders, 1 mice; zero distance -0.227, largest distance -1.082; largest - zero -0.855 ± 0.052, p < 0.001; decoders with a drop 1.000
  PM: 8 decoders, 2 mice; zero distance -0.231, largest distance -1.322; largest - zero -1.086 ± 0.082, p < 0.001; decoders with a drop 1.000
  AM: 5 decoders, 1 mice; zero distance -0.395, largest distance -1.459; largest - zero -1.061 ± 0.106, p < 0.001; decoders with a drop 1.000
  MM: 9 decoders, 3 mice; zero distance -0.537, largest distance -1.210; largest - zero -0.718 ± 0.129, p < 0.001; decoders with a drop 1.000
  RSC: 18 decoders, 3 mice; zero distance -0.405, largest distance -1.265; largest - zero -0.852 ± 0.074, p < 0.001; decoders with a drop 1.000
  A: 11 decoders, 1 mice; zero distance -0.684, largest distance -1.148; largest - zero -0.465 ± 0.054, p < 0.001; decoders with a drop 1.000
dynamic choice, Spearman: 53 decoders from 3 mice; cue offsets {0.76: 53}; largest distance 10 to 10 bins
  mean performance by train-test distance (bins): {'dist0': 0.602, 'dist1': 0.56, 'dist2': 0.497, 'dist3': 0.439, 'dist4': 0.386, 'dist5': 0.345, 'dist6': 0.3, 'dist7': 0.259, 'dist8': 0.199, 'dist9': 0.121, 'dist10': 0.122}
  V1: 2 decoders, 1 mice; zero distance 0.595, largest distance 0.143; largest - zero -0.451 ± 0.044, p < 0.001; decoders with a drop 1.000
  PM: 8 decoders, 2 mice; zero distance 0.629, largest distance 0.158; largest - zero -0.491 ± 0.046, p < 0.001; decoders with a drop 1.000
  AM: 5 decoders, 1 mice; zero distance 0.561, largest distance 0.036; largest - zero -0.526 ± 0.037, p < 0.001; decoders with a drop 1.000

  MM: 9 decoders, 3 mice; zero distance 0.546, largest distance 0.078; largest - zero -0.460 ± 0.075, p < 0.001; decoders with a drop 1.000
  RSC: 18 decoders, 3 mice; zero distance 0.641, largest distance 0.142; largest - zero -0.500 ± 0.101, p < 0.001; decoders with a drop 1.000
  A: 11 decoders, 1 mice; zero distance 0.585, largest distance 0.134; largest - zero -0.450 ± 0.047, p < 0.001; decoders with a drop 1.000

Result: the claims file's own reading also gives the drop. On the 53 decoders of the 6 pool sessions with the cue offset at 0.76 (3 mice; the paper has 698 decoders from such sessions), cue decoding over the cue-present bins falls from -0.451 at the training position to -1.251 seven bins away, and dynamic choice from 0.602 to 0.122 ten bins away. In every area the largest minus zero distance is negative (cue -0.465 to -1.086, dynamic choice -0.450 to -0.526, p < 0.001 in each) and every decoder drops. V1 has 2 decoders and AM 5, each from one mouse, and A 11 from one mouse, so the per-area tests rest on few mice. The outcome of the main notebook (reproduced, on the assumption that the cue test may use sessions with other cue offsets) therefore does not depend on that assumption.

Summary¶

Claim What this notebook shows Recorded outcome
fig3q-cue-decoding on every V1 and PM field decoder, V1 is above PM under all 12 readings (by 0.027 to 0.062, mostly not significant), so V1 is the highest area reproduced (proposed after round 2, agreed in the third round)
fig8cd-position-specific-decoders on the 6 sessions with the 0.76 cue offset, cue and dynamic-choice decoding fall with train-test distance in every area reproduced (verified); the assumption is not needed

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. As in the main notebook (same decoders, readings and statistics), but on a selection of sessions: all 41 pool sessions with a V1 or PM field for fig3q, and the 6 pool sessions with the 0.76 cue offset for fig8cd, against the paper's 974 and 698 decoders. For fig3q the other four areas are taken from the main notebook's subsample of 18 sessions.

What was found about the dataset. The 6 sessions the files mark with the cue offset at 0.76 hold 53 decoders, against the 698 of the Figure 8C legend; 8 of the 172 multi-plane layer 2/3 session files have that offset.

In [6]:
print(f"notebook ran in {(time.time() - T0) / 60:.1f} min")

notebook ran in 14.6 min