Daume et al. (2024) on dandiset 000673: hippocampal PAC controls

Paper. Daume J, Kaminski J, Schjetnan AGP, Salimpour Y, Khan U, Kyzar M, Reed CM, Anderson WS, Valiante TA, Mamelak AN, Rutishauser U (2024). Control of working memory by phase-amplitude coupling of human hippocampal neurons. Nature 629, 393-401. doi:10.1038/s41586-024-07309-z.

Dandiset. DANDI:000673, version 0.240118.2135, the version the paper cites (reference 59); see the header of reproduce_daume_2024.ipynb for the later versions.

Text the claims were written from. The version of record (PMC11078732) with its Methods, legends and Supplementary Information, and the paper's code (SBCAT-release-NWB at a25e8cc3).

This notebook tests the controls of Extended Data Fig. 3 on the hippocampal PAC channels. It repeats, for the 586 hippocampal channels only, the channel processing of reproduce_daume_2024_lfp.ipynb (each channel's random draws are seeded by session, channel and artefact reading, so the PAC channels are the same 144 and 125 as there), and adds the theta cycle features of bycycle, the phase-phase coupling of each phase and amplitude centre pair, and a comodulogram with 0.5 Hz phase steps. It is kept apart so that neither notebook exceeds the time limit.

Claim	Tested here
edf3-hippocampus-pac-controls	yes

import os
for _k in ["OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"]:
    os.environ[_k] = "1"   # many small computations run in parallel processes
import time
import warnings

import numpy as np
import pandas as pd
from scipy import stats

import sbcat

pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
warnings.filterwarnings("ignore")

started = time.time()
S = sbcat.load_all(sbcat.list_assets())
print(f"loaded trials and units of {len(S)} files in {time.time() - started:.0f} s")
t0 = time.time()
R = sbcat.run_all_lfp(S, seed=201, workers=6, controls=True, areas=["hippocampus"])
print(f"LFP of {len(R)} sessions, {sum(len(r['channels']) for r in R)} hippocampal channels, processed in {(time.time() - t0) / 60:.1f} min")

loaded trials and units of 44 files in 133 s

LFP of 38 sessions, 586 hippocampal channels, processed in 20.7 min

HG = np.where((sbcat.AMP_CENTRES >= 70) & (sbcat.AMP_CENTRES <= 140))[0]
READINGS = ["a", "b"]
rows = []
for r in R:
    for c in r["channels"]:
        for rd in READINGS:
            o = c["readings"][rd]
            base = {"session": c["session"], "subject": c["subject"], "channel": c["channel"], "reading": rd,
                    "kept_trials": int(c["masks"][rd].sum()), "n_trials": len(c["masks"][rd])}
            if o is None:
                rows.append({**base, "n_per_load": 0})
                continue
            rows.append({**base, "n_per_load": o["n_per_load"],
                         "z_thetaHG": float(np.mean([o["z_all"][p][HG].mean() for p in sbcat.THETA])),
                         "z1_thetaHG": float(np.mean([o["z_l1"][p][HG].mean() for p in sbcat.THETA])),
                         "z3_thetaHG": float(np.mean([o["z_l3"][p][HG].mean() for p in sbcat.THETA]))})
ch = pd.DataFrame(rows)
chan_index = {(c["session"], c["channel"]): c for r in R for c in r["channels"]}


def hpc_sig(rd):
    return ch[(ch.reading == rd) & (ch.z_thetaHG > 1.64)]

Data check

print("hippocampal LFP channels:", int((ch.reading == "a").sum()), "(paper 586) in", ch.session.nunique(), "sessions of", ch.subject.nunique(), "patients")
for rd in READINGS:
    d = hpc_sig(rd)
    print(f"reading {rd}: {len(d)} channels with significant theta-high-gamma PAC (as in reproduce_daume_2024_lfp.ipynb), "
          f"load 1 mean z {d.z1_thetaHG.mean():.2f}, load 3 {d.z3_thetaHG.mean():.2f}")
cyc = np.concatenate([c["cycles"] for c in chan_index.values()])
print("theta cycles (bycycle, 3-7 Hz) with their peak in 0 to 2.5 s, per correct trial: all cycles median %.0f (range %d to %d); "
      "burst cycles median %.0f (range %d to %d); trials with no burst cycle %.1f%%" % (
          np.median(cyc[:, 5]), cyc[:, 5].min(), cyc[:, 5].max(), np.median(cyc[:, 2]), cyc[:, 2].min(), cyc[:, 2].max(),
          100 * (cyc[:, 2] == 0).mean()))
allc = cyc[:, 5] > 0
print("range of the mean asymmetry per trial (all cycles): rise-decay %.3f to %.3f, peak-trough %.3f to %.3f" % (
    (cyc[allc, 3] / cyc[allc, 5]).min(), (cyc[allc, 3] / cyc[allc, 5]).max(), (cyc[allc, 4] / cyc[allc, 5]).min(), (cyc[allc, 4] / cyc[allc, 5]).max()))
print("phase-phase pairs (phase centre, amplitude centre, n):", [(p, int(a), int(round(a / p))) for p, a in sbcat.PPC_PAIRS])

hippocampal LFP channels: 586 (paper 586) in 38 sessions of 30 patients
reading a: 144 channels with significant theta-high-gamma PAC (as in reproduce_daume_2024_lfp.ipynb), load 1 mean z 4.62, load 3 3.80
reading b: 125 channels with significant theta-high-gamma PAC (as in reproduce_daume_2024_lfp.ipynb), load 1 mean z 4.14, load 3 3.46
theta cycles (bycycle, 3-7 Hz) with their peak in 0 to 2.5 s, per correct trial: all cycles median 12 (range 6 to 16); burst cycles median 0 (range 0 to 6); trials with no burst cycle 99.9%
range of the mean asymmetry per trial (all cycles): rise-decay 0.107 to 0.840, peak-trough 0.283 to 0.746
phase-phase pairs (phase centre, amplitude centre, n): [(4, 70, 18), (4, 75, 19), (4, 80, 20), (4, 85, 21), (4, 90, 22), (4, 95, 24), (4, 100, 25), (4, 105, 26), (4, 110, 28), (4, 115, 29), (4, 120, 30), (4, 125, 31), (4, 130, 32), (4, 135, 34), (4, 140, 35), (6, 70, 12), (6, 75, 12), (6, 80, 13), (6, 85, 14), (6, 90, 15), (6, 95, 16), (6, 100, 17), (6, 105, 18), (6, 110, 18), (6, 115, 19), (6, 120, 20), (6, 125, 21), (6, 130, 22), (6, 135, 22), (6, 140, 23)]

edf3-hippocampus-pac-controls

The hippocampal load effect on PAC was not explained by theta or gamma power, theta waveform asymmetry, theta-gamma phase-phase coupling, preferred phase or PAC peak frequency.

Every test is a two-sided paired permutation t-test between loads over the PAC channels, on the trials of the PAC analysis (the equal-count subsamples of each load), as the criterion and the Extended Data Fig. 3 legend say.

• Power, theta (3-7 Hz) and high gamma (70-140 Hz): readings (a) mean squared Hilbert amplitude of the PAC filters; (b) the Welch spectrum of 0 to 2.5 s averaged over the band (log10).
• Waveform asymmetry: bycycle (version in the lock file) with f_range 3-7 Hz and its default burst thresholds, on each trial's -0.5 to 3 s epoch, cycles with their peak in 0 to 2.5 s; rise-decay (time_rdsym) and peak-trough (time_ptsym) symmetry averaged over cycles per load ("averaged estimates ... across cycles during the maintenance period from the same trials used for our PAC analysis within each load", Methods, PAC). The Methods do not say whether only burst cycles count; burst cycles (bycycle's default use) and all cycles are both shown. Each is also tested against 0.5 (symmetric), as the criterion asks.
• Phase-phase coupling, n:m phase locking value: readings (a) the ratio of the band centres, 1:21 (5 and 105 Hz); (b) each pair of a theta phase centre (4, 6 Hz) and a high-gamma amplitude centre (70 to 140 Hz) of O3, n the ratio rounded to an integer, averaged over pairs.
• Preferred phase: the load difference of each channel's preferred phase, readings (a) V test against 0; (b) permutation t-test of the wrapped differences against 0.
• PAC peak frequency (Extended Data Fig. 3h): raw MI with phase centres 3 to 7 Hz in 0.5 Hz steps (2 Hz bandwidth) and amplitude centres 70 to 140 Hz (bandwidth twice the phase centre), per load; "the frequency bin for which PAC was maximal for the theta and the gamma band separately" is read as the maximum of the mean over the other axis (main), and as the coordinates of the single largest bin (variant).

rng = np.random.default_rng(309)
freqs = R[0]["freqs"]
centres = np.radians(np.linspace(-180, 180, sbcat.N_BINS + 1)[:-1] + 180 / sbcat.N_BINS)
n_samp = sbcat.EPOCH_N - 2 * sbcat.CUT
edf3_tables = {}
for rd in READINGS:
    d = hpc_sig(rd)
    V = {k: {1: [], 3: []} for k in ["theta power, Hilbert", "high-gamma power, Hilbert", "theta power, Welch",
                                     "high-gamma power, Welch", "rise-decay symmetry, burst cycles",
                                     "peak-trough symmetry, burst cycles", "rise-decay symmetry, all cycles",
                                     "peak-trough symmetry, all cycles", "phase locking 1:21 (a)",
                                     "phase locking, each centre pair (b)", "theta peak frequency (Hz)",
                                     "gamma peak frequency (Hz)", "theta peak frequency, largest bin (Hz)",
                                     "gamma peak frequency, largest bin (Hz)"]}
    pref = {1: [], 3: []}
    for r in d.itertuples():
        c = chan_index[(r.session, r.channel)]
        o = c["readings"][rd]
        for L in (1, 3):
            pos = np.searchsorted(c["trial_index"], o[f"trials_l{L}"])
            V["theta power, Hilbert"][L].append(np.log10(c["theta_power"][pos].mean()))
            V["high-gamma power, Hilbert"][L].append(np.log10(c["hg_power"][pos].mean()))
            V["theta power, Welch"][L].append(np.log10(o[f"psd_l{L}"][(freqs >= 3) & (freqs <= 7)].mean()))
            V["high-gamma power, Welch"][L].append(np.log10(o[f"psd_l{L}"][(freqs >= 70) & (freqs <= 140)].mean()))
            cy = c["cycles"][pos].sum(0)
            V["rise-decay symmetry, burst cycles"][L].append(cy[0] / cy[2] if cy[2] else np.nan)
            V["peak-trough symmetry, burst cycles"][L].append(cy[1] / cy[2] if cy[2] else np.nan)
            V["rise-decay symmetry, all cycles"][L].append(cy[3] / cy[5] if cy[5] else np.nan)
            V["peak-trough symmetry, all cycles"][L].append(cy[4] / cy[5] if cy[5] else np.nan)
            V["phase locking 1:21 (a)"][L].append(o[f"ppc_l{L}"])
            V["phase locking, each centre pair (b)"][L].append(np.mean(np.abs(c["ppc_pairs"][pos].sum(0)) / (len(pos) * n_samp)))
            fm = o[f"fine_mi_l{L}"]                           # phases x amplitudes
            V["theta peak frequency (Hz)"][L].append(sbcat.FINE_PHASE[np.argmax(fm.mean(1))])
            V["gamma peak frequency (Hz)"][L].append(sbcat.FINE_AMPS[np.argmax(fm.mean(0))])
            ip, ia = np.unravel_index(np.argmax(fm), fm.shape)
            V["theta peak frequency, largest bin (Hz)"][L].append(sbcat.FINE_PHASE[ip])
            V["gamma peak frequency, largest bin (Hz)"][L].append(sbcat.FINE_AMPS[ia])
            P = o[f"P_thetaHG_l{L}"]
            pref[L].append(np.angle(np.sum(P * np.exp(1j * centres))))
    res = []
    for k, v in V.items():
        t = sbcat.perm_ttest_paired(v[3], v[1], rng)
        res.append({"test": k, "n": t["n"], "load1": t["mean_y"], "load3": t["mean_x"], "t(3-1)": t["t"], "p": t["p"]})
    dphi = np.angle(np.exp(1j * (np.array(pref[3]) - np.array(pref[1]))))
    n = len(dphi)
    m = np.mean(np.exp(1j * dphi))
    u = n * np.abs(m) * np.cos(np.angle(m)) * np.sqrt(2 / n)
    res.append({"test": "preferred phase 3-1, V test against 0 (a)", "n": n, "t(3-1)": u, "p": 1 - stats.norm.cdf(u)})
    t = sbcat.perm_ttest_1samp(np.degrees(dphi), rng)
    res.append({"test": "preferred phase 3-1, wrapped, t-test against 0 (b) (deg)", "n": t["n"], "load3": t["mean"], "t(3-1)": t["t"], "p": t["p"]})
    tab = pd.DataFrame(res)
    edf3_tables[rd] = tab
    print(f"artefact reading {rd}: {len(d)} PAC channels")
    print(tab.round(4).to_string(index=False))
    print("mean resultant length of the preferred phase differences: %.3f; circular mean %.1f deg; |difference| < 20 deg in %d of %d channels" % (
        np.abs(m), np.degrees(np.angle(m)), int((np.abs(np.degrees(dphi)) < 20).sum()), n))
    sym = []
    for k in ["rise-decay symmetry, burst cycles", "peak-trough symmetry, burst cycles", "rise-decay symmetry, all cycles", "peak-trough symmetry, all cycles"]:
        both = np.nanmean([V[k][1], V[k][3]], 0)
        t = sbcat.perm_ttest_1samp(both - 0.5, rng)
        sym.append({"asymmetry against 0.5 (mean of loads)": k, "n": t["n"], "mean": t["mean"] + 0.5, "t": t["t"], "p": t["p"]})
    print(pd.DataFrame(sym).round(4).to_string(index=False))
    print("channels without a burst cycle in a load:", int(np.sum(np.isnan(V["rise-decay symmetry, burst cycles"][1]) | np.isnan(V["rise-decay symmetry, burst cycles"][3]))))
    for k in ["theta peak frequency (Hz)", "gamma peak frequency (Hz)"]:
        print(f"{k}: range load 1 {min(V[k][1])} to {max(V[k][1])}, load 3 {min(V[k][3])} to {max(V[k][3])}")
    print()
print("paper: no load difference in any of these (Extended Data Fig. 3b, c, d, g, h; n = 137)")

artefact reading a: 144 PAC channels
                                                    test   n   load1   load3  t(3-1)      p
                                    theta power, Hilbert 144  2.6732  2.6815  1.2017 0.2294
                               high-gamma power, Hilbert 144  1.0749  1.0816  2.0863 0.0353
                                      theta power, Welch 144  1.8345  1.8456  1.4606 0.1553
                                 high-gamma power, Welch 144 -0.6816 -0.6726  2.6910 0.0060
                       rise-decay symmetry, burst cycles   8  0.4045  0.4566  1.6443 0.1736
                      peak-trough symmetry, burst cycles   8  0.4744  0.4656 -0.3031 0.7928
                         rise-decay symmetry, all cycles 144  0.4995  0.4987 -0.9558 0.3404
                        peak-trough symmetry, all cycles 144  0.4977  0.4973 -0.6712 0.4982
                                  phase locking 1:21 (a) 144  0.0075  0.0074 -0.1268 0.9014
                     phase locking, each centre pair (b) 144  0.0096  0.0096 -0.2272 0.8204
                               theta peak frequency (Hz) 144  3.5764  3.5139 -0.7399 0.4875
                               gamma peak frequency (Hz) 144 86.8056 92.9514  3.0360 0.0030
                  theta peak frequency, largest bin (Hz) 144  3.6840  3.5139 -1.7850 0.0795
                  gamma peak frequency, largest bin (Hz) 144 88.6111 90.7639  0.9705 0.3407
               preferred phase 3-1, V test against 0 (a) 144     NaN     NaN 15.8916 0.0000
preferred phase 3-1, wrapped, t-test against 0 (b) (deg) 144     NaN -0.7176 -0.3871 0.7090
mean resultant length of the preferred phase differences: 0.937; circular mean -1.6 deg; |difference| < 20 deg in 106 of 144 channels
asymmetry against 0.5 (mean of loads)   n   mean       t      p
    rise-decay symmetry, burst cycles  14 0.4498 -1.9697 0.0723
   peak-trough symmetry, burst cycles  14 0.4953 -0.2648 0.7923
      rise-decay symmetry, all cycles 144 0.4991 -0.6198 0.5406
     peak-trough symmetry, all cycles 144 0.4975 -1.3177 0.1893
channels without a burst cycle in a load: 136
theta peak frequency (Hz): range load 1 3.0 to 7.0, load 3 3.0 to 7.0
gamma peak frequency (Hz): range load 1 70 to 140, load 3 70 to 140

artefact reading b: 125 PAC channels
                                                    test   n   load1   load3  t(3-1)      p
                                    theta power, Hilbert 125  2.7303  2.7380  0.9406 0.3630
                               high-gamma power, Hilbert 125  1.1164  1.1219  1.5503 0.1269
                                      theta power, Welch 125  1.8870  1.9025  1.7454 0.0856
                                 high-gamma power, Welch 125 -0.6395 -0.6329  1.8955 0.0584
                       rise-decay symmetry, burst cycles   5  0.4229  0.4566  0.7526 0.6247
                      peak-trough symmetry, burst cycles   5  0.4700  0.4762  0.1517 0.8800
                         rise-decay symmetry, all cycles 125  0.5001  0.4995 -0.5568 0.5786
                        peak-trough symmetry, all cycles 125  0.4969  0.4973  0.6054 0.5421
                                  phase locking 1:21 (a) 125  0.0089  0.0089 -0.0568 0.9552
                     phase locking, each centre pair (b) 125  0.0109  0.0108 -0.3857 0.6967
                               theta peak frequency (Hz) 125  3.5840  3.6200  0.4015 0.7268
                               gamma peak frequency (Hz) 125 88.2400 91.4400  1.3877 0.1696
                  theta peak frequency, largest bin (Hz) 125  3.6640  3.5760 -0.8087 0.4353
                  gamma peak frequency, largest bin (Hz) 125 86.9200 87.9200  0.4937 0.6342
               preferred phase 3-1, V test against 0 (a) 125     NaN     NaN 14.9686 0.0000
preferred phase 3-1, wrapped, t-test against 0 (b) (deg) 125     NaN -0.7463 -0.4377 0.6655
mean resultant length of the preferred phase differences: 0.947; circular mean -0.7 deg; |difference| < 20 deg in 93 of 125 channels
asymmetry against 0.5 (mean of loads)   n   mean       t      p
    rise-decay symmetry, burst cycles  11 0.4521 -1.4101 0.1925
   peak-trough symmetry, burst cycles  11 0.4967 -0.1299 0.8983
      rise-decay symmetry, all cycles 125 0.4998 -0.1561 0.8814
     peak-trough symmetry, all cycles 125 0.4971 -1.4028 0.1717
channels without a burst cycle in a load: 120
theta peak frequency (Hz): range load 1 3.0 to 7.0, load 3 3.0 to 7.0
gamma peak frequency (Hz): range load 1 70 to 140, load 3 70 to 140

paper: no load difference in any of these (Extended Data Fig. 3b, c, d, g, h; n = 137)

Result: not reproduced, resting on an assumption. Under artefact reading (a), the one adopted, over 144 PAC channels, two load tests give p < 0.05: high-gamma power is higher in load 3 (Hilbert t = 2.09, p = 0.035; Welch t = 2.69, p = 0.006), and the gamma peak frequency is higher in load 3 (93.0 against 86.8 Hz in load 1, t = 3.04, p = 0.003; as the coordinate of the single largest bin, p = 0.34). Theta power (p = 0.23, 0.16), rise-decay and peak-trough symmetry (all cycles p = 0.34 and 0.50; burst cycles, found in only 8 channels in both loads, p = 0.17 and 0.79), phase-phase locking (1:21 p = 0.90; each centre pair p = 0.82), the theta peak frequency (p = 0.49; largest bin p = 0.080) and the preferred phase t-test (p = 0.71) do not differ. The V test of the phase differences against 0 gives p < 0.0001 (mean resultant length 0.94, circular mean -1.6 degrees), which shows differences concentrated at 0; whether the criterion's rule should count it is asked of a reader. No asymmetry differs from 0.5 (p >= 0.072). Under artefact reading (b) (125 channels) no load test gives p < 0.05 (high-gamma power p = 0.13 and 0.058; gamma peak frequency p = 0.17), so the outcome rests on reading (a).

Summary

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.

Claim	Outcome
edf3-hippocampus-pac-controls	not reproduced, resting on artefact reading (a) (high-gamma power and gamma peak frequency differ between loads)

How the analysis differs from the paper's: the artefact rejection is the automated part only, reading (a) adopted as in reproduce_daume_2024_lfp.ipynb; bycycle 1.2 is run with a replacement of its burst detection that applies the same default thresholds to a writable copy (it fails under pandas 3 otherwise); the power estimate, the phase-phase ratio, whether only burst cycles count, and the definition of the peak frequency are not given and are each shown under two readings.

What was found about the dataset: nothing beyond reproduce_daume_2024_lfp.ipynb. With bycycle's default burst thresholds, 99.9% of the 0 to 2.5 s windows of correct trials have no burst cycle.

print(f"notebook ran in {(time.time() - started) / 60:.1f} min")

notebook ran in 22.9 min