Top-ranked units (Gerken, Darcher et al. 2025, dandiset 001616)¶
Paper: Gerken, Darcher et al. (2025), eLife reviewed preprint 106758 v1, doi:10.7554/eLife.106758.1, full text. Dandiset: DANDI:001616, version 0.260702.0824.
Companion to reproduce_gerken_2025.ipynb, with the same loading, samples, folds and logistic
regression (summer.py). It holds two claims:
| Claim | Status |
|---|---|
| fig7b-top500-overlap | tested |
| fig7a-500-neurons | tested, with the logistic regression in place of the LSTM |
The notebook takes about 8 minutes.
In [1]:
import os
for var in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
os.environ.setdefault(var, "1")
import time
import numpy as np
import pandas as pd
import summer
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
started_all = started = time.time()
sessions = summer.load_all()
units = summer.unit_table(sessions)
spikes = summer.movie_spikes(sessions)
print(f"loaded {len(sessions)} files, {len(units)} units, from dandiset {summer.DANDISET_ID} "
f"version {summer.DANDISET_VERSION} in {time.time() - started:.0f} s")
loaded 29 files, 2286 units, from dandiset 001616 version 0.260702.0824 in 116 s
In [2]:
edges = sessions[0]["edges_s"]
S = summer.sample_windows(spikes, edges)
n = S.shape[0]
y = summer.bin_labels(sessions[0], "summer")
print(f"patients {units.patient.nunique()}, units {len(units)}; by region {units.region.value_counts().to_dict()}")
print(f"samples x units: {S.shape}; Summer positive in {np.mean(y == 1):.3f} of the samples, values {sorted(np.unique(y).tolist())}")
for f in range(5):
print(f"test fold {f + 1}: " + "; ".join(f"validation fold {v + 1} train {len(a)}, validation {len(b)}, test {len(c)}"
for v in range(5) if v != f for a, b, c in [summer.split(n, f, val_fold=v)]))
patients 29, units 2286; by region {'H': 846, 'A': 559, 'EC': 493, 'PHC': 293, 'PIC': 68, 'other': 27}
samples x units: (59080, 2286); Summer positive in 0.527 of the samples, values [0, 1]
test fold 1: validation fold 2 train 30822, validation 7194, test 7194; validation fold 3 train 26202, validation 7194, test 7194; validation fold 4 train 26202, validation 7194, test 7194; validation fold 5 train 30822, validation 7194, test 7194
test fold 2: validation fold 1 train 30822, validation 7194, test 7194; validation fold 3 train 30822, validation 7194, test 7194; validation fold 4 train 26202, validation 7194, test 7194; validation fold 5 train 26202, validation 7194, test 7194
test fold 3: validation fold 1 train 26202, validation 7194, test 7194; validation fold 2 train 30822, validation 7194, test 7194; validation fold 4 train 30822, validation 7194, test 7194; validation fold 5 train 26202, validation 7194, test 7194
test fold 4: validation fold 1 train 26202, validation 7194, test 7194; validation fold 2 train 26202, validation 7194, test 7194; validation fold 3 train 30822, validation 7194, test 7194; validation fold 5 train 30822, validation 7194, test 7194
test fold 5: validation fold 1 train 30822, validation 7194, test 7194; validation fold 2 train 26202, validation 7194, test 7194; validation fold 3 train 26202, validation 7194, test 7194; validation fold 4 train 30822, validation 7194, test 7194
fig7b-top500-overlap¶
Across the five rankings of neurons by logistic regression weights for Summer, the top 500 neurons share 78 neurons, against about 1 expected by chance.
In [3]:
rng_rank = np.random.default_rng(37)
started = time.time()
keys, jobs = [], []
for f in range(5):
for v in [j for j in range(5) if j != f]:
keys.append((f, v))
jobs.append((summer.fit_logreg, (S, y, *summer.split(n, f, val_fold=v))))
kfits = dict(zip(keys, summer.run_jobs(jobs)))
print(f"{len(jobs)} sub-split fits in {time.time() - started:.0f} s")
rankings, zero_in_top = [], []
for f in range(5):
coefs = [kfits[(f, v)]["coef"] for v in range(5) if v != f]
order = summer.ranking(coefs, rng_rank)
rankings.append(order)
nz = (np.abs(np.array(coefs)) > 0).any(0)
zero_in_top.append(int((~nz[order[:500]]).sum()))
print(f"test fold {f + 1}: C {[kfits[(f, v)]['C'] for v in range(5) if v != f]}, units with a non-zero coefficient in any of the 4 models "
f"{int(nz.sum())}, in all 4 {int((np.abs(np.array(coefs)) > 0).all(0).sum())}; top 500 tied at zero {zero_in_top[-1]}; "
f"sub-split test kappa {np.round([kfits[(f, v)]['kappa'] for v in range(5) if v != f], 3).tolist()}")
N = S.shape[1]
rows = []
for k_ in (100, 350, 500, 750, 1000):
common = set(rankings[0][:k_])
for o in rankings[1:]:
common &= set(o[:k_])
sims = [len(set.intersection(*[set(rng_rank.choice(N, k_, replace=False)) for _ in range(5)])) for _ in range(200)]
rows.append({"top": k_, "common": len(common), "chance_formula": N * (k_ / N) ** 5, "chance_simulated": np.mean(sims),
"paper_chance": {100: 0, 350: 0, 500: 1, 750: 9, 1000: 37}[k_]})
ov = pd.DataFrame(rows)
print(ov.round(2).to_string(index=False))
common500 = set(rankings[0][:500]).intersection(*[set(o[:500]) for o in rankings[1:]])
print("common 500 by region:", units.loc[sorted(common500)].region.value_counts().to_dict(),
"; patients:", units.loc[sorted(common500)].patient.nunique())
# Variant: zero-tied units ordered by unit index instead of a random draw (same order in every ranking).
order_idx = [np.lexsort((np.arange(N), -np.abs(np.array([kfits[(f, v)]["coef"] for v in range(5) if v != f])).mean(0),
-(np.abs(np.array([kfits[(f, v)]["coef"] for v in range(5) if v != f])) > 0).sum(0))) for f in range(5)]
print("variant, ties by unit index: common in top 500:", len(set(order_idx[0][:500]).intersection(*[set(o[:500]) for o in order_idx[1:]])))
20 sub-split fits in 73 s test fold 1: C [0.003, 0.01, 0.003, 0.03], units with a non-zero coefficient in any of the 4 models 1281, in all 4 72; top 500 tied at zero 0; sub-split test kappa [0.459, 0.424, 0.488, 0.434] test fold 2: C [0.001, 0.01, 0.01, 0.03], units with a non-zero coefficient in any of the 4 models 1291, in all 4 47; top 500 tied at zero 0; sub-split test kappa [0.258, 0.299, 0.439, 0.288] test fold 3: C [0.01, 0.01, 0.003, 0.01], units with a non-zero coefficient in any of the 4 models 1195, in all 4 95; top 500 tied at zero 0; sub-split test kappa [0.302, 0.296, 0.316, 0.376] test fold 4: C [0.003, 0.003, 0.03, 0.03], units with a non-zero coefficient in any of the 4 models 1437, in all 4 78; top 500 tied at zero 0; sub-split test kappa [0.26, 0.397, 0.406, 0.108] test fold 5: C [0.03, 0.001, 0.003, 0.003], units with a non-zero coefficient in any of the 4 models 1101, in all 4 28; top 500 tied at zero 0; sub-split test kappa [0.324, 0.422, 0.51, 0.335]
top common chance_formula chance_simulated paper_chance
100 20 0.00 0.00 0
350 65 0.19 0.20 0
500 95 1.14 1.12 1
750 156 8.69 8.81 9
1000 252 36.62 35.70 37
common 500 by region: {'A': 34, 'H': 21, 'EC': 20, 'PHC': 13, 'other': 4, 'PIC': 3} ; patients: 26
variant, ties by unit index: common in top 500: 95
Result: partially reproduced. The five rankings share 95 units in their top 500, against 78 in the paper (tolerance 70.2 to 85.8) and 1.14 expected by chance (1.12 simulated), 83 times chance. No unit in a top 500 was tied at zero, so the tie-break does not matter (95 with ties by unit index). The overlaps for the top 100, 350, 750 and 1000 are 20, 65, 156 and 252, against chance values of 0.00, 0.19, 8.69 and 36.6 that match the paper's (0, 0, 9, 37). The 95 units come from 26 patients and every region.
fig7a-500-neurons¶
Decoding Summer from the 500 top-ranked neurons gives the peak performance (kappa 0.38), above that of the full population.
For each test fold, the logistic regression is refitted on each of its four sub-splits with only the top k units of that fold's ranking, and scored on the shared test fold. Kappa is averaged over the four sub-splits, then over the five test folds (Methods); the SEM is over the five fold means.
In [4]:
rng_top = np.random.default_rng(39)
SIZES = (23, 114, 150, 229, 500, 750, 1000, 1500)
started = time.time()
keys = [(k_, f, v) for k_ in SIZES for f in range(5) for v in range(5) if v != f]
tfits = dict(zip(keys, summer.run_jobs([(summer.fit_logreg, (S, y, *summer.split(n, f, val_fold=v), rankings[f][:k_]))
for k_, f, v in keys])))
tfits |= {(N, f, v): r for (f, v), r in kfits.items()}
print(f"{len(keys)} fits on top-ranked units in {time.time() - started:.0f} s")
rows = []
for k_ in SIZES + (N,):
per_fold = np.array([np.mean([tfits[(k_, f, v)]["kappa"] for v in range(5) if v != f]) for f in range(5)])
pmax = max(summer.permutation_p(tfits[(k_, f, v)], rng_top) for f in range(5) for v in range(5) if v != f)
rows.append({"units": k_, "percent": round(100 * k_ / N, 1), "kappa": per_fold.mean(), "sem": per_fold.std(ddof=1) / np.sqrt(5),
"fold_means": np.round(per_fold, 3).tolist(), "max_perm_p": pmax})
top = pd.DataFrame(rows).set_index("units")
print(top.round(4).to_string())
print("paper (LSTM): peak 0.38 at 500 units; 0.52 at 150 units with a ranking that saw the test data")
best = top.kappa.idxmax()
print(f"500 units {top.kappa[500]:.4f}; all {N} units {top.kappa[N]:.4f}; difference {top.kappa[500] - top.kappa[N]:.4f}; "
f"highest at {best} units ({top.kappa[best]:.4f})")
160 fits on top-ranked units in 162 s
percent kappa sem fold_means max_perm_p units 23 1.0 0.3410 0.0288 [0.424, 0.345, 0.29, 0.266, 0.38] 0.001 114 5.0 0.3779 0.0228 [0.428, 0.321, 0.329, 0.425, 0.387] 0.001 150 6.6 0.3892 0.0164 [0.441, 0.37, 0.345, 0.406, 0.384] 0.001 229 10.0 0.3939 0.0194 [0.455, 0.354, 0.351, 0.397, 0.413] 0.001 500 21.9 0.3751 0.0262 [0.429, 0.34, 0.35, 0.312, 0.445] 0.001 750 32.8 0.3713 0.0260 [0.43, 0.336, 0.344, 0.309, 0.437] 0.001 1000 43.7 0.3629 0.0274 [0.434, 0.325, 0.335, 0.298, 0.422] 0.001 1500 65.6 0.3621 0.0285 [0.428, 0.331, 0.334, 0.287, 0.429] 0.001 2286 100.0 0.3570 0.0293 [0.451, 0.321, 0.322, 0.293, 0.398] 0.001 paper (LSTM): peak 0.38 at 500 units; 0.52 at 150 units with a ranking that saw the test data 500 units 0.3751; all 2286 units 0.3570; difference 0.0181; highest at 229 units (0.3939)
Result: partially reproduced. With the logistic regression, the top 500 units decode Summer with kappa 0.375, above all 2286 units (0.357), by 0.018, less than the SEM over folds (0.026 and 0.029). The peak is not at 500: 229 units give 0.394 and 150 units 0.389; from 229 units upwards the kappa falls steadily, and the top 23 units still give 0.341. Every fit is above chance by the test-label shuffle (p = 0.001). The paper's LSTM peaks at 0.38 with 500 units.
In [5]:
print(f"total run time {(time.time() - started_all) / 60:.1f} min")
total run time 6.2 min