Reproducing the human iEEG factor and cluster analyses of Kauvar, Richman, Liu et al. (2025) from dandiset 001347¶

  • Paper: Kauvar, Richman, Liu et al., "Conserved brain-wide emergence of emotional response from sensory experience in humans and mice", Science 388, eadt3971 (2025), doi:10.1126/science.adt3971. Text read: the version of record with its supplement, as the claims file says.
  • Dataset: DANDI 001347, version 0.250528.0702 (the only published version, the paper's ref. 102), streamed. The claims file places the human iEEG data there (data: elsewhere).
  • Added in the first revision, after the verifier found that the derived human files (sub-H) hold what three of the blocked claims need.

The sub-H files are the paper's own intermediate products: per-channel spectrograms already z-scored and averaged over the train and test halves of each session (sgramPerChannel), per-trial spectrograms averaged within each Yeo7 network (sgramPerYeoPerTrial), and the paper's NMF factors and loadings. The spectrogram step itself (multitaper spectrograms of the raw iEEG in sub-SD*_ecephys, z-scoring, channel selection) is not repeated here; the tests start from the deposited spectrograms.

Claim Tested here
fig2g-six-factors yes, NMF of the deposited train-half spectrograms
fig2i-factors-1-2-ketamine yes, loadings from that NMF and from the deposit
fig3b-cluster-test-networks yes, the Yeo7 networks, from the deposited per-trial spectrograms
fig2d, figs14g, fig1f, figs12, fig6c, fig6e, fig6g, fig6i no (see the attempt record)
In [1]:
import warnings
import numpy as np
import pandas as pd
import scipy.stats
from sklearn.metrics import r2_score

import behavior as B

warnings.filterwarnings("ignore")
urls = B.asset_urls()
S, ch, ff, tt = B.load_sgram_per_channel(urls)
L, ch_dep, F_dep = B.load_nmf_deposit(urls)
yeo, yff, ytt = B.load_yeo_trials(urls)
print(f"dandiset {B.DANDISET}, version {B.VERSION}")

dandiset 001347, version 0.250528.0702

Data check¶

What the derived files hold: shapes, channels and patients, units and ranges, and whether the deposited loadings name the same channels.

In [2]:
pts = np.array([c.split("_")[1] for c in ch])
print("sgramPerChannel:", {k: v.shape for k, v in S.items()}, "| channels per patient:", pd.Series(pts).value_counts().sort_index().to_dict())
print(f"frequencies {ff[0]:.2f} to {ff[-1]:.2f} Hz ({len(ff)}), times {tt[0]:.3f} to {tt[-1]:.3f} s ({len(tt)}; the Methods' 0 to 1.4 s)")
print("value range (z):", {k: (round(float(v.min()), 2), round(float(v.max()), 2)) for k, v in S.items()}, "| NaN:", sum(int(np.isnan(v).sum()) for v in S.values()))
print("the loadings file lists the same 458 channels in the same order:", bool((ch == ch_dep).all()))
print("deposited loadings, channels with a value per factor:", (~np.isnan(L["pre"])).sum(0).tolist())
print("\nsgramPerYeoPerTrial:")
for r, d in yeo.items():
    common = set(d["preinfusion"]) & set(d["infusion"])
    print(f"  {r}: preinfusion {len(d['preinfusion'])} trials, infusion {len(d['infusion'])}, paired {len(common)}, patients {sorted({c[0] for c in common})}")
print(f"  frequencies {yff[0]:.2f} to {yff[-1]:.2f} Hz ({len(yff)}), times {ytt[0]:.3f} to {ytt[-1]:.3f} s ({len(ytt)})")

sgramPerChannel: {'inf_test': (458, 65, 28), 'inf_train': (458, 65, 28), 'post_test': (458, 65, 28), 'post_train': (458, 65, 28), 'pre_test': (458, 65, 28), 'pre_train': (458, 65, 28)} | channels per patient: {'SD168': 86, 'SD172': 97, 'SD181': 58, 'SD185': 43, 'SD186': 47, 'SD211': 44, 'SD245': 83}
frequencies 0.00 to 100.00 Hz (65), times 0.000 to 1.367 s (28; the Methods' 0 to 1.4 s)
value range (z): {'inf_test': (-1.96, 3.94), 'inf_train': (-2.11, 3.9), 'post_test': (-1.91, 4.09), 'post_train': (-2.16, 3.94), 'pre_test': (-1.96, 3.92), 'pre_train': (-2.04, 3.98)} | NaN: 0
the loadings file lists the same 458 channels in the same order: True
deposited loadings, channels with a value per factor: [312, 260, 185, 158, 370, 378]

sgramPerYeoPerTrial:
  Frontoparietal: preinfusion 178 trials, infusion 178, paired 178, patients [np.str_('SD172'), np.str_('SD181'), np.str_('SD185'), np.str_('SD211'), np.str_('SD245')]
  Default: preinfusion 214 trials, infusion 214, paired 214, patients [np.str_('SD168'), np.str_('SD172'), np.str_('SD181'), np.str_('SD185'), np.str_('SD186'), np.str_('SD211'), np.str_('SD245')]
  Dorsal Attention: preinfusion 153 trials, infusion 153, paired 153, patients [np.str_('SD168'), np.str_('SD172'), np.str_('SD186'), np.str_('SD211'), np.str_('SD245')]
  Limbic: preinfusion 214 trials, infusion 214, paired 214, patients [np.str_('SD168'), np.str_('SD172'), np.str_('SD181'), np.str_('SD185'), np.str_('SD186'), np.str_('SD211'), np.str_('SD245')]
  Ventral Attention: preinfusion 214 trials, infusion 214, paired 214, patients [np.str_('SD168'), np.str_('SD172'), np.str_('SD181'), np.str_('SD185'), np.str_('SD186'), np.str_('SD211'), np.str_('SD245')]
  Somatomotor: preinfusion 214 trials, infusion 214, paired 214, patients [np.str_('SD168'), np.str_('SD172'), np.str_('SD181'), np.str_('SD185'), np.str_('SD186'), np.str_('SD211'), np.str_('SD245')]
  frequencies 0.00 to 100.00 Hz (65), times -0.531 to 1.444 s (40)

fig2g-six-factors¶

Nonnegative matrix factorization of the preinfusion puff-triggered spectrograms of 458 channels from seven participants is best described by six factors, which explain 0.47 of the variance on held-out trials.

Criterion: as the Methods (Factorization) and table S3. Train-half spectrograms (every other pre- and postinfusion trial), positive and negative parts each factored with k; the other half is the test set. Reproduced when the held-out explained variance at six factors is 0.47 +- 0.047 and the test-set curve shows its elbow at six.

Method details from the paper's code (Zenodo 10.5281/zenodo.14920570, human_ieeg/scripts/analyze_airpuff_ieeg_sgram_load_allpt.py and ketman/analyze/decode_utils.py): the training array is the mean of the pre- and postinfusion train halves (X_prepost_train), flattened to channels x (frequency x time); sklearn NMF with solver 'mu' and random_state 1 (get_nmf); k from 1 to 19 per sign (nmf_k_determination, 2 to 38 factors); the held-out explained variance is r2_score of the test half against the train reconstruction (do_refit_W_for_test = False). The paper picks the elbow by eye; here the elbow is the factor count with the largest drop in marginal gain (largest second difference), with the point farthest from the chord of the curve as a second reading.

In [3]:
shape = (len(ff), len(tt))
train = ((S["pre_train"] + S["post_train"]) / 2).reshape(len(ch), -1)
test = ((S["pre_test"] + S["post_test"]) / 2).reshape(len(ch), -1)
rows = []
models = {}
for k in range(1, 20):
    m = B.fit_nmf(train, k)
    models[k] = m
    rec = m[0].transform(B.split_pos_neg(train)[0]) @ m[0].components_ - m[1].transform(B.split_pos_neg(train)[1]) @ m[1].components_
    rows.append({"factors": 2 * k, "R2 train": r2_score(train.ravel(), rec.ravel()), "R2 test (held out)": r2_score(test.ravel(), rec.ravel())})
curve = pd.DataFrame(rows)
print(curve.round(4).to_string(index=False))
y = curve["R2 test (held out)"].to_numpy(); x = curve["factors"].to_numpy()
gain = np.diff(y)
elbow_d2 = x[1:-1][np.argmax(gain[:-1] - gain[1:])]
d = np.abs((x[-1] - x[0]) * (y[0] - y) - (x[0] - x) * (y[-1] - y[0])) / np.hypot(x[-1] - x[0], y[-1] - y[0])
elbow_chord = x[np.argmax(d)]
ve6 = float(curve.loc[curve.factors == 6, "R2 test (held out)"].iloc[0])
print(f"\nheld-out explained variance at six factors: {ve6:.3f} (expected 0.47 +- 0.047; within: {abs(ve6 - 0.47) <= 0.047})")
print(f"marginal gains in held-out R2 (2 to 12 factors): {np.round(gain[:5], 4).tolist()}")
print(f"elbow by the largest drop in marginal gain: {elbow_d2} factors; by the farthest point from the chord (2 to 38): {elbow_chord} factors")
for seed in (0, 2, 3):
    m = B.fit_nmf(train, 3, random_state=seed)
    rec = m[0].transform(B.split_pos_neg(train)[0]) @ m[0].components_ - m[1].transform(B.split_pos_neg(train)[1]) @ m[1].components_
    print(f"  six factors with NMF random_state {seed}: held-out R2 {r2_score(test.ravel(), rec.ravel()):.4f}")

 factors  R2 train  R2 test (held out)
       2    0.3535              0.3291
       4    0.4698              0.4097
       6    0.5416              0.4655
       8    0.5669              0.4800
      10    0.5913              0.4848
      12    0.6101              0.4903
      14    0.6284              0.4958
      16    0.6406              0.4992
      18    0.6523              0.5013
      20    0.6639              0.5057
      22    0.6754              0.5069
      24    0.6836              0.5054
      26    0.6905              0.5064
      28    0.6960              0.5058
      30    0.7027              0.5038
      32    0.7100              0.5019
      34    0.7175              0.5012
      36    0.7233              0.4994
      38    0.7281              0.4991

held-out explained variance at six factors: 0.466 (expected 0.47 +- 0.047; within: True)
marginal gains in held-out R2 (2 to 12 factors): [0.0805, 0.0559, 0.0145, 0.0048, 0.0056]
elbow by the largest drop in marginal gain: 6 factors; by the farthest point from the chord (2 to 38): 8 factors

  six factors with NMF random_state 0: held-out R2 0.4655

  six factors with NMF random_state 2: held-out R2 0.4655

  six factors with NMF random_state 3: held-out R2 0.4655

Result: reproduced, resting on the assumption that the elbow is where the marginal gain drops most. The six-factor NMF of the train half explains 0.466 of the variance of the held-out half (0.4655 under three other NMF seeds), within 0.47 +- 0.047. The largest drop in marginal gain comes after six factors (gains 0.081 from 2 to 4 factors, 0.056 from 4 to 6, 0.015 from 6 to 8), so that rule puts the elbow at six; the farthest point from the chord of the 2-to-38-factor curve puts it at eight. The paper picks the elbow by eye.

fig2i-factors-1-2-ketamine¶

Factors 1 and 2 (a fast high-frequency rise and a delayed low-frequency fall) were less represented during ketamine than the average of pre- and postinfusion.

Criterion: channels with a mean loading of at least 0.1 on the factor over the pre- and postinfusion test trials; infusion loading against the mean of pre- and postinfusion loadings by a mixed-effects linear model grouped by participant, FDR across the six factors (table S3). Reproduced when factors 1 and 2 both have a negative coefficient with corrected p < 0.05. Expected loadings: factor 1 preinfusion 0.51 +- 0.04, infusion 0.32 +- 0.04; factor 2 preinfusion 0.35 +- 0.04, infusion 0.21 +- 0.02.

Loadings are the NMF transform of each session's test half with the six-factor fit above (code: nmf_pos.transform, nmf_neg.transform); the 0.1 threshold applies to the loading of the mean of the pre- and postinfusion test halves (W_prepost, compute_factors_using = 'prepost'). The code orders factors by the time of their largest absolute value; the figure numbers follow the deposited factor files, so each recomputed factor is matched to the deposited factor it correlates with. The deposited loadings are tested too.

In [4]:
m6 = models[3]
fac = B.nmf_factors(m6, shape)
corr = np.array([[np.corrcoef(F_dep[i].ravel(), fac[j].ravel())[0, 1] for j in range(6)] for i in range(6)])
match = corr.argmax(1)
print("correlation of each deposited factor (rows, figure numbers 0 to 5) with the recomputed factors (columns):")
print(pd.DataFrame(corr, index=[f"deposit {i}" for i in range(6)]).round(3).to_string())
print("recomputed factor matched to each figure number:", match.tolist(), "| one-to-one:", len(set(match)) == 6)
peak_t = tt[np.argmax(np.abs(fac).max(1), axis=1)]
print("time of the largest |value| of each factor in figure order (s):", np.round(peak_t[match], 3).tolist(),
      "| sign:", ["rise" if j < 3 else "fall" for j in match])
W = {s: B.nmf_loadings(m6, S[f"{s}_test"].reshape(len(ch), -1))[:, match] for s in ("pre", "inf", "post")}
Wpp = B.nmf_loadings(m6, ((S["pre_test"] + S["post_test"]) / 2).reshape(len(ch), -1))[:, match]
for s in W:
    W[s][Wpp < 0.1] = np.nan
print("recomputed loadings equal the deposited ones (NaN pattern and values):",
      all(np.allclose(W[s], L[s], equal_nan=True, atol=1e-6) for s in W))
exp = {1: (0.51, 0.04, 0.32, 0.04), 2: (0.35, 0.04, 0.21, 0.02)}
for lab, (a, b, c) in (("recomputed", (W["pre"], W["inf"], W["post"])), ("deposited", (L["pre"], L["inf"], L["post"]))):
    t2i = B.mixed_loading_test(a, b, c, ch)
    print(f"\n{lab} loadings, mixed-effects model (infusion against the mean of pre and post), FDR over six factors:")
    print(t2i.round(4).to_string(index=False))
    for f_, (ep, tp, ei, ti) in exp.items():
        r = t2i.loc[f_]
        print(f"  factor {f_}: pre {r.pre:.3f} (within {abs(r.pre - ep) <= tp}), infusion {r.infusion:.3f} (within {abs(r.infusion - ei) <= ti}), "
              f"coefficient {r.coefficient:.4f}, corrected p {r['corrected p']:.2g}")

correlation of each deposited factor (rows, figure numbers 0 to 5) with the recomputed factors (columns):
               0      1      2      3      4      5
deposit 0  1.000 -0.124  0.092  0.393  0.171  0.015
deposit 1  0.092 -0.635  1.000 -0.147  0.521 -0.120
deposit 2  0.171 -0.836  0.521 -0.492  1.000 -0.110
deposit 3  0.015  0.031 -0.120 -0.386 -0.110  1.000
deposit 4 -0.124  1.000 -0.635  0.521 -0.836  0.031
deposit 5  0.393  0.521 -0.147  1.000 -0.492 -0.386
recomputed factor matched to each figure number: [0, 2, 4, 5, 1, 3] | one-to-one: True
time of the largest |value| of each factor in figure order (s): [0.152, 0.203, 0.405, 0.304, 0.456, 0.506] | sign: ['rise', 'rise', 'fall', 'fall', 'rise', 'fall']
recomputed loadings equal the deposited ones (NaN pattern and values): True

recomputed loadings, mixed-effects model (infusion against the mean of pre and post), FDR over six factors:
 factor  channels  patients    pre  infusion   post  coefficient      p  converged  corrected p
      0       312         7 0.4825    0.4173 0.4500      -0.0489 0.0808       True       0.1213
      1       260         7 0.5076    0.3180 0.3785      -0.1250 0.0000       True       0.0000
      2       185         7 0.3542    0.2137 0.3025      -0.1147 0.0000       True       0.0000
      3       158         7 0.2817    0.2141 0.2280      -0.0407 0.0079       True       0.0159
      4       370         7 0.3237    0.3866 0.4232       0.0132 0.4060       True       0.4060
      5       378         7 0.3920    0.4299 0.4413       0.0132 0.2845       True       0.3414
  factor 1: pre 0.508 (within True), infusion 0.318 (within True), coefficient -0.1250, corrected p 2.6e-05
  factor 2: pre 0.354 (within True), infusion 0.214 (within True), coefficient -0.1147, corrected p 5.5e-08

deposited loadings, mixed-effects model (infusion against the mean of pre and post), FDR over six factors:
 factor  channels  patients    pre  infusion   post  coefficient      p  converged  corrected p
      0       312         7 0.4825    0.4173 0.4500      -0.0489 0.0808       True       0.1213
      1       260         7 0.5076    0.3180 0.3785      -0.1250 0.0000       True       0.0000
      2       185         7 0.3542    0.2137 0.3025      -0.1147 0.0000       True       0.0000
      3       158         7 0.2817    0.2141 0.2280      -0.0407 0.0079      False       0.0159
      4       370         7 0.3237    0.3866 0.4232       0.0132 0.4060       True       0.4060
      5       378         7 0.3920    0.4299 0.4413       0.0132 0.2845       True       0.3414
  factor 1: pre 0.508 (within True), infusion 0.318 (within True), coefficient -0.1250, corrected p 2.6e-05
  factor 2: pre 0.354 (within True), infusion 0.214 (within True), coefficient -0.1147, corrected p 5.5e-08

Result: reproduced. The NMF fitted here reproduces the deposited loadings exactly (same values and the same channels below the 0.1 threshold). Factor 1 (a rise peaking 0.20 s after the puff) has loadings 0.508 before and 0.318 during the infusion, and factor 2 (a fall peaking at 0.41 s) 0.354 and 0.214, all within tolerance. The mixed-effects model gives negative coefficients for both (-0.125 and -0.115, corrected p = 2.6e-05 and 5.5e-08, 260 and 185 channels from 7 patients). Factor 3 also falls (corrected p = 0.016), and factors 0, 4 and 5 do not change.

fig3b-cluster-test-networks¶

Permutation cluster tests found ketamine-altered spectrotemporal clusters in the limbic, ventral attention and somatomotor networks, and in regions including insula, cingulate, orbitofrontal and posteromedial cortex, supramarginal gyrus, motor cortex, hippocampus, basal ganglia and thalamus, while many other regions were not changed.

Criterion: as the Methods (Permutation cluster test) and table S3. Regional (or network) spectrograms per trial, trials paired by position in the puff sequence, a one-sample t-test of infusion minus preinfusion with mne permutation_cluster_1samp_test, cluster-forming threshold 0.01, the lowest cluster p of each network FDR-corrected across networks, significance at corrected p < 0.01; networks sampled in at least two participants. Reproduced when limbic, ventral attention and somatomotor networks are significant, and partially reproduced when one or two of them are.

Trials are paired by patient and trial number. The number of permutations is the code's default (512, permutation_cluster_tests); with it the smallest attainable p is 1/512. The permutation seed is a standard choice; three seeds are shown. The regional test (35 regions, 624 MB) is not run; the criterion is decided by the networks.

In [5]:
rows = []
for seed in (0, 1, 2):
    pmin = {}
    for r, d in yeo.items():
        n, npt, pv = B.paired_cluster_test(d, seed=seed)
        pmin[r] = pv.min() if len(pv) else 1.0
        if seed == 0:
            rows.append({"network": r, "pairs": n, "patients": npt, "clusters": len(pv), "clusters p < 0.01": int((pv < 0.01).sum())})
    q = scipy.stats.false_discovery_control(list(pmin.values()))
    for i, r in enumerate(pmin):
        if seed == 0:
            rows[i].update({"lowest p": pmin[r], "corrected p": q[i], "significant": q[i] < 0.01})
        else:
            rows[i][f"corrected p, seed {seed}"] = q[i]
t3b = pd.DataFrame(rows)
print(t3b.round(4).to_string(index=False))
need = ["Limbic", "Ventral Attention", "Somatomotor"]
print("\nlimbic, ventral attention and somatomotor significant (seed 0):", t3b.set_index("network").loc[need, "significant"].to_dict())
with B._open(urls["sub-H/sub-H_ses-ketChangePermutationClusterTest.nwb"]) as f:
    dep = {r: f[f"scratch/{r}_contour"][:] for r in yeo if f"{r}_contour" in f["scratch"]}
print("deposited cluster masks (ketChangePermutationClusterTest), cells inside a displayed cluster per network:",
      {r: int(np.nansum(np.asarray(v) != 0)) for r, v in dep.items()})

          network  pairs  patients  clusters  clusters p < 0.01  lowest p  corrected p  significant  corrected p, seed 1  corrected p, seed 2
   Frontoparietal    178         5        27                  0    0.0254       0.0381        False               0.0410               0.0381
          Default    214         7        17                  0    0.0566       0.0680        False               0.0539               0.0750
 Dorsal Attention    153         5        20                  0    0.0703       0.0703        False               0.1016               0.1172
           Limbic    214         7        19                  2    0.0020       0.0039         True               0.0039               0.0039
Ventral Attention    214         7        20                  3    0.0020       0.0039         True               0.0039               0.0039
      Somatomotor    214         7        23                  1    0.0020       0.0039         True               0.0039               0.0039

limbic, ventral attention and somatomotor significant (seed 0): {np.str_('Limbic'): True, np.str_('Ventral Attention'): True, np.str_('Somatomotor'): True}

deposited cluster masks (ketChangePermutationClusterTest), cells inside a displayed cluster per network: {np.str_('Frontoparietal'): 0, np.str_('Default'): 0, np.str_('Dorsal Attention'): 0, np.str_('Limbic'): 192, np.str_('Ventral Attention'): 219, np.str_('Somatomotor'): 100}

Result: reproduced. Of the six Yeo7 networks in the deposit, the limbic, ventral attention and somatomotor networks each have a cluster at the smallest attainable p (1/512), and their corrected p is 0.0039, below 0.01, under all three permutation seeds. The frontoparietal, default and dorsal attention networks are not significant (corrected p 0.038 to 0.12). The deposited cluster masks mark the same three networks and no other. The regional part of the statement (insula, cingulate and other regions) is not tested here.

Summary¶

Claim Outcome
fig2g-six-factors reproduced (assumes the elbow is the largest drop in marginal gain)
fig2i-factors-1-2-ketamine reproduced
fig3b-cluster-test-networks reproduced

How the analysis differs from the paper's: it starts from the spectrograms the paper deposited (sub-H), already z-scored and averaged, and does not recompute them from the raw iEEG. The elbow of fig. 2G is chosen by a rule where the paper chose by eye. The regional cluster tests of Fig. 3C are not run.

What was found about the dataset: the derived files carry what these three analyses need, and the NMF recomputed from them gives back the deposited loadings exactly. The factor order of the deposit is not the code's order by peak time (deposited factors 2 and 3 are the code's 3 and 2), so the figure numbers come only from the deposited factor files.