Daume et al. (2024) on dandiset 000673: LFP phase-amplitude coupling¶
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 claims on LFP phase-amplitude coupling (PAC). Each channel's LFP is
streamed, cut into -0.5 to 3 s epochs around maintenance onset, filtered trial by trial as the
paper's code does, and Tort's modulation index (MI) is z-scored against 200 trial-shuffled
surrogates (sbcat.py says where each parameter comes from). Every result is shown under both
readings of the artefact rejection (O1): (a) on the epochs, (b) on whole trials. Where they
disagree, (a) is adopted: on the one channel whose rejected trials the paper's code ships
(artif_trials.mat, sub-26 session 1, channel 33), every trial (a) removes is among the paper's,
while (b) removes 5 trials the paper kept (data check below).
| Claim | Tested here |
|---|---|
| fig2a-theta-gamma-pac | yes |
| fig2b-pac-channel-proportions | yes |
| fig2c-hippocampus-load | yes |
| fig2c-no-load-effect-other-areas | yes |
| edf2b-hippocampus-load-patients | yes |
| edf2c-low-gamma-no-load-effect | yes |
| tables1-slow-fast-theta | yes |
| fig2f-rt-single-trial-pac | yes |
| edf3-hippocampus-pac-controls | in reproduce_daume_2024_lfp_controls.ipynb |
| edf3a-kappa-load | yes |
| edf3e-nesting | yes |
| edf3f-raw-mi-load | yes |
| methods-spectral-slope | yes |
| fig3h-pac-category-firing | yes |
import os
for _k in ["OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"]:
os.environ[_k] = "1" # many small fits run in parallel processes
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy import stats
import sbcat
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
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)
print(f"LFP of {len(R)} sessions, {sum(len(r['channels']) for r in R)} channels, processed in {(time.time() - t0) / 60:.1f} min")
loaded trials and units of 44 files in 133 s
LFP of 43 sessions, 1922 channels, processed in 39.2 min
HG = np.where((sbcat.AMP_CENTRES >= 70) & (sbcat.AMP_CENTRES <= 140))[0]
LG = np.where((sbcat.AMP_CENTRES >= 30) & (sbcat.AMP_CENTRES <= 55))[0]
READINGS = ["a", "b"]
def band_z(zdict, phases, amp_idx):
"""Mean z over the given phase centres and amplitude indices (O3)."""
return float(np.mean([zdict[p][amp_idx].mean() for p in phases]))
rows = []
for r in R:
s = S[r["session"]]
for c in r["channels"]:
for rd in READINGS:
o = c["readings"][rd]
base = {"session": c["session"], "subject": c["subject"], "path": c["path"], "channel": c["channel"],
"area": c["area"], "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": band_z(o["z_all"], sbcat.THETA, HG), "z_thetaLG": band_z(o["z_all"], sbcat.THETA, LG),
"z1_thetaHG": band_z(o["z_l1"], sbcat.THETA, HG), "z3_thetaHG": band_z(o["z_l3"], sbcat.THETA, HG),
"z1_thetaLG": band_z(o["z_l1"], sbcat.THETA, LG), "z3_thetaLG": band_z(o["z_l3"], sbcat.THETA, LG),
"mi1_thetaHG": band_z(o["mi_l1"], sbcat.THETA, HG), "mi3_thetaHG": band_z(o["mi_l3"], sbcat.THETA, HG),
**{f"z_{nm}_{g}": band_z(o["z_all"], ph, ix) for nm, ph in [("slowA", (2, 4)), ("fastA", (6, 8))] for g, ix in [("HG", HG), ("LG", LG)]},
**{f"z{L}_{nm}_{g}": band_z(o[f"z_l{L}"], ph, ix) for L in (1, 3) for nm, ph in [("slowA", (2, 4)), ("fastA", (6, 8))] for g, ix in [("HG", HG), ("LG", LG)]},
**{f"zb_{k}_{g}_{cn}": v for (k, g, cn), v in o["bands"].items()}})
ch = pd.DataFrame(rows)
chan_index = {(c["session"], c["channel"]): c for r in R for c in r["channels"]}
print(ch.groupby("reading").size().to_dict())
{'a': 1922, 'b': 1922}
Data check¶
one = ch[ch.reading == "a"]
print("LFP channels by area:", one.area.value_counts().reindex(sbcat.AREAS).to_dict(), "; total", len(one))
info = pd.DataFrame([{"path": r["path"], "lfp_samples": r["lfp_shape"][0], "channels": r["lfp_shape"][1],
"lfp_s": r["lfp_shape"][0] / 400, "task_end_s": S[r["session"]]["trials"].stop_time.max()} for r in R])
print("LFP length minus the end of the last trial (s): min %.1f, max %.1f" % ((info.lfp_s - info.task_end_s).min(), (info.lfp_s - info.task_end_s).max()))
print("LFP rate and unit as stored:", {(s["lfp_info"]["rate"], s["lfp_info"]["unit"], s["lfp_info"]["conversion"]) for s in S if s["lfp_info"]})
sd = pd.Series([c["sd_uV"] for r in R for c in r["channels"]])
p99 = pd.Series([c["p99_abs_uV"] for r in R for c in r["channels"]])
print("per-channel s.d. of the epochs (stored units): median %.1f, range %.2f to %.1f; 99th percentile of |LFP|: median %.1f, max %.1f" % (
sd.median(), sd.min(), sd.max(), p99.median(), p99.max()))
for rd in READINGS:
d = ch[ch.reading == rd]
rem = d.n_trials - d.kept_trials
print(f"reading {rd}: trials removed per channel {rem.mean():.1f} +/- {rem.std():.1f} ({100 * rem.sum() / d.n_trials.sum():.1f}%); "
f"trials per load after subsampling: median {d.n_per_load.median():.0f}, range {d.n_per_load.min()} to {d.n_per_load.max()}; "
f"channels with fewer than 3 per load (left out): {int((d.n_per_load < 3).sum())}")
print("paper: 20.4 +/- 13.9 trials (14.6%) removed per wire, with visual inspection")
fp = pd.Series({r["path"]: r["lfp_fingerprint"] for r in R})
print("files sharing the sum of their first 4,000 LFP samples (duplicates):", fp[fp.duplicated(keep=False)].to_dict() or "none")
for rd in READINGS:
d = ch[ch.reading == rd]
rem = d.n_trials - d.kept_trials
print(f"reading {rd}: trials removed per channel, percentiles 0/25/50/75/90/100: {np.percentile(rem, [0, 25, 50, 75, 90, 100]).round(1).tolist()}; "
f"channels with none removed {int((rem == 0).sum())} of {len(rem)}; histogram (0-4, 5-9, ..., 50+): "
f"{np.bincount(np.minimum(rem.to_numpy().astype(int) // 5, 10), minlength=11).tolist()}")
LFP channels by area: {'hippocampus': 586, 'amygdala': 421, 'preSMA': 283, 'dACC': 325, 'vmPFC': 307} ; total 1922
LFP length minus the end of the last trial (s): min 1.1, max 211.7
LFP rate and unit as stored: {(400.0, 'microvolts', 1.0)}
per-channel s.d. of the epochs (stored units): median 34.4, range 4.26 to 525.2; 99th percentile of |LFP|: median 104.8, max 3000.9
reading a: trials removed per channel 16.0 +/- 12.0 (11.4%); trials per load after subsampling: median 58, range 15 to 70; channels with fewer than 3 per load (left out): 0
reading b: trials removed per channel 36.2 +/- 22.8 (25.9%); trials per load after subsampling: median 48, range 0 to 68; channels with fewer than 3 per load (left out): 4
paper: 20.4 +/- 13.9 trials (14.6%) removed per wire, with visual inspection
files sharing the sum of their first 4,000 LFP samples (duplicates): none
reading a: trials removed per channel, percentiles 0/25/50/75/90/100: [0.0, 9.0, 13.0, 19.0, 29.9, 106.0]; channels with none removed 3 of 1922; histogram (0-4, 5-9, ..., 50+): [123, 436, 520, 391, 186, 73, 61, 35, 26, 18, 53]
reading b: trials removed per channel, percentiles 0/25/50/75/90/100: [0.0, 20.0, 30.0, 47.0, 69.0, 132.0]; channels with none removed 2 of 1922; histogram (0-4, 5-9, ..., 50+): [18, 57, 154, 239, 227, 237, 225, 135, 113, 81, 436]
# The paper's code (NWB_samplePAC_SU.m, SBCAT-release-NWB at a25e8cc3) ships the artefact trials of
# its example: sub-26, session 1, LFP channel 33 (1-based), in artif_trials.mat. Compare them with the
# trials each reading of O1 removes on that channel.
import io
import urllib.request
from scipy.io import loadmat
url = ("https://raw.githubusercontent.com/rutishauserlab/SBCAT-release-NWB/a25e8cc3f61630a428f190381745257d4506f148/"
"NWB_SBCAT_analysis/helpers/internal/PAC/artif_trials.mat")
with urllib.request.urlopen(url) as r_:
paper_artif = set(int(v) - 1 for v in loadmat(io.BytesIO(r_.read()))["artif_trials"].ravel())
si26 = [i for i, s in enumerate(S) if "sub-26_ses-1_" in s["path"]]
print("sub-26 session 1:", [S[i]["path"] for i in si26])
c26 = chan_index[(si26[0], 32)]
print(f"channel 33 ({c26['location']}): the paper's code removes {len(paper_artif)} trials")
for rd in READINGS:
rem = set(np.where(~c26["masks"][rd])[0])
print(f" reading {rd} removes {len(rem)}; {len(rem & paper_artif)} of them among the paper's {len(paper_artif)}, "
f"{len(rem - paper_artif)} that the paper kept")
sub-26 session 1: ['sub-26/sub-26_ses-1_ecephys+image.nwb'] channel 33 (hippocampus_right): the paper's code removes 17 trials reading a removes 4; 4 of them among the paper's 17, 0 that the paper kept reading b removes 9; 4 of them among the paper's 17, 5 that the paper kept
# Clock check by contrast: RMS of the probe-locked trial-averaged LFP 0.1 to 0.6 s after onset over
# that of -0.6 to -0.1 s, at the true probe onsets and at onsets moved 1 or 2 s back (into maintenance).
ev = pd.DataFrame([{**{f"shift_{k}": v for k, v in c["evoked"].items()}, "path": c["path"]} for r in R for c in r["channels"]])
print(ev.drop(columns="path").median().round(2).to_dict())
per_sess = ev.groupby("path").median()
print("sessions whose median ratio is larger at the true onset than at both shifts:",
int(((per_sess["shift_0.0"] > per_sess["shift_-1.0"]) & (per_sess["shift_0.0"] > per_sess["shift_-2.0"])).sum()), "of", len(per_sess))
print("sessions where it is not:", per_sess.index[~((per_sess["shift_0.0"] > per_sess["shift_-1.0"]) & (per_sess["shift_0.0"] > per_sess["shift_-2.0"]))].tolist())
{'shift_-2.0': 1.58, 'shift_-1.0': 1.4, 'shift_0.0': 2.49}
sessions whose median ratio is larger at the true onset than at both shifts: 37 of 43
sessions where it is not: ['sub-1/sub-1_ses-1_ecephys+image.nwb', 'sub-11/sub-11_ses-1_ecephys+image.nwb', 'sub-12/sub-12_ses-2_ecephys+image.nwb', 'sub-16/sub-16_ses-1_ecephys+image.nwb', 'sub-23/sub-23_ses-1_ecephys+image.nwb', 'sub-28/sub-28_ses-1_ecephys+image.nwb']
fig2a-theta-gamma-pac¶
Across all recorded channels, PAC during maintenance was strongest between theta phase (3-7 Hz) and the amplitude of a low (30-55 Hz) and a high (70-140 Hz) gamma band.
for rd in READINGS:
PH = sbcat.PHASE_CENTRES if rd == "a" else np.array([2, 4, 6, 8]) # reading (b) computes phase 2-8 Hz only
Z = np.array([[c["readings"][rd]["z_all"][int(p)] for p in PH] for r in R for c in r["channels"] if c["readings"][rd] is not None])
mean_z = Z.mean(0)
print(f"reading {rd}: {len(Z)} channels; mean z comodulogram (rows phase {PH[0]}-{PH[-1]} Hz, columns amplitude 30-150 Hz):")
print(pd.DataFrame(mean_z, index=PH, columns=sbcat.AMP_CENTRES).round(2).to_string())
ip, ia = np.unravel_index(np.argmax(mean_z), mean_z.shape)
print(f"largest mean z {mean_z[ip, ia]:.2f} at phase {PH[ip]} Hz, amplitude {sbcat.AMP_CENTRES[ia]} Hz")
for p in sbcat.THETA:
prof = mean_z[list(PH).index(p)]
peaks = [int(sbcat.AMP_CENTRES[i]) for i in range(1, len(prof) - 1) if prof[i] > prof[i - 1] and prof[i] > prof[i + 1]]
print(f"phase {p} Hz: local maxima of the amplitude profile at {peaks} Hz")
by_area = {a: Z[[ (c['area'] == a) for r in R for c in r['channels'] if c['readings'][rd] is not None]].mean(0) for a in sbcat.AREAS}
print("largest mean z per area (phase, amplitude):", {a: (int(PH[np.unravel_index(m.argmax(), m.shape)[0]]), int(sbcat.AMP_CENTRES[np.unravel_index(m.argmax(), m.shape)[1]]), round(float(m.max()), 2)) for a, m in by_area.items()})
fig, ax = plt.subplots(figsize=(6, 3))
im = ax.imshow(mean_z, origin="lower", aspect="auto", extent=[27.5, 152.5, PH[0] - 1, PH[-1] + 1], cmap="viridis")
ax.set_xlabel("amplitude frequency (Hz)"); ax.set_ylabel("phase frequency (Hz)"); plt.colorbar(im, label="mean z (MI)")
ax.set_title(f"all channels, reading {rd}; range {mean_z.min():.2f} to {mean_z.max():.2f}")
plt.show()
reading a: 1922 channels; mean z comodulogram (rows phase 2-14 Hz, columns amplitude 30-150 Hz):
30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150
2 1.69 1.59 1.61 1.75 1.83 1.79 1.80 1.79 1.87 1.93 1.95 1.90 1.84 1.77 1.71 1.64 1.59 1.57 1.56 1.56 1.54 1.52 1.49 1.46 1.44
4 1.85 1.84 1.93 2.05 2.06 1.95 1.86 1.88 2.02 2.13 2.16 2.13 2.05 1.95 1.85 1.75 1.69 1.65 1.62 1.58 1.54 1.52 1.50 1.48 1.44
6 1.45 1.44 1.46 1.45 1.42 1.42 1.39 1.42 1.48 1.52 1.53 1.52 1.50 1.46 1.40 1.33 1.26 1.20 1.16 1.11 1.08 1.05 1.02 0.99 0.97
8 1.57 1.33 1.28 1.25 1.24 1.21 1.19 1.16 1.14 1.14 1.14 1.14 1.12 1.11 1.07 1.02 0.99 0.96 0.92 0.88 0.84 0.81 0.79 0.76 0.75
10 1.56 1.28 1.16 1.09 1.07 1.07 1.03 1.03 1.07 1.07 1.04 1.03 1.03 1.02 0.99 0.96 0.92 0.89 0.87 0.83 0.79 0.75 0.72 0.69 0.69
12 1.67 1.37 1.22 1.12 1.09 1.04 1.03 1.02 1.05 1.07 1.04 1.01 1.07 1.09 1.04 1.00 0.98 0.98 0.98 0.97 0.95 0.94 0.94 0.92 0.87
14 1.44 1.52 1.30 1.19 1.23 1.18 1.15 1.17 1.22 1.24 1.21 1.20 1.22 1.23 1.20 1.19 1.20 1.19 1.16 1.13 1.12 1.11 1.10 1.07 1.01
largest mean z 2.16 at phase 4 Hz, amplitude 80 Hz
phase 4 Hz: local maxima of the amplitude profile at [50, 80] Hz
phase 6 Hz: local maxima of the amplitude profile at [40, 80] Hz
largest mean z per area (phase, amplitude): {'hippocampus': (4, 75, 3.14), 'amygdala': (4, 50, 4.0), 'preSMA': (2, 60, 2.0), 'dACC': (2, 30, 1.05), 'vmPFC': (2, 80, 2.03)}
reading b: 1918 channels; mean z comodulogram (rows phase 2-8 Hz, columns amplitude 30-150 Hz):
30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150
2 1.34 1.28 1.26 1.36 1.38 1.37 1.38 1.37 1.43 1.49 1.50 1.47 1.42 1.38 1.34 1.28 1.24 1.24 1.23 1.22 1.21 1.20 1.18 1.15 1.11
4 1.39 1.45 1.55 1.64 1.61 1.50 1.44 1.49 1.60 1.71 1.73 1.69 1.63 1.55 1.48 1.39 1.33 1.30 1.29 1.26 1.24 1.22 1.22 1.20 1.16
6 1.09 1.11 1.14 1.11 1.10 1.13 1.14 1.14 1.18 1.22 1.23 1.20 1.17 1.15 1.11 1.06 1.04 1.01 0.99 0.95 0.93 0.90 0.88 0.84 0.81
8 1.17 0.98 0.97 0.94 0.98 1.04 1.02 0.96 0.92 0.91 0.88 0.87 0.87 0.88 0.83 0.77 0.74 0.72 0.70 0.69 0.67 0.64 0.61 0.59 0.57
largest mean z 1.73 at phase 4 Hz, amplitude 80 Hz
phase 4 Hz: local maxima of the amplitude profile at [45, 80] Hz
phase 6 Hz: local maxima of the amplitude profile at [40, 60, 80] Hz
largest mean z per area (phase, amplitude): {'hippocampus': (4, 75, 2.58), 'amygdala': (4, 45, 3.24), 'preSMA': (2, 60, 1.65), 'dACC': (2, 30, 0.7), 'vmPFC': (2, 80, 1.54)}
Result: reproduced. Under reading (a), over 1,922 channels, the largest mean z (2.16) is at phase 4 Hz and amplitude 80 Hz; at 4 Hz the amplitude profile has local maxima at 50 Hz (within 30 to 55) and 80 Hz (within 70 to 140). Under reading (b), computed for phase centres 2 to 8 Hz over 1,918 channels, the largest mean z (1.73) is again at 4 Hz and 80 Hz, with local maxima at 45 and 80 Hz. Both readings meet the criterion.
fig2b-pac-channel-proportions¶
137 of 586 hippocampal channels showed significant theta-high-gamma PAC; the proportions of such channels in the hippocampus, amygdala and vmPFC exceeded the trial-shuffled null, while pre-SMA and dACC had few such channels.
paper_counts = {"hippocampus": 137, "amygdala": 130, "preSMA": 5, "dACC": 14, "vmPFC": 40}
fig2b_rows = []
for rd in READINGS:
for a in sbcat.AREAS:
d = ch[(ch.reading == rd) & (ch.area == a) & (ch.n_per_load >= 3)]
sig = d.z_thetaHG > 1.64
surr = np.array([chan_index[(r.session, r.channel)]["readings"][rd]["surr_thetaHG"] for r in d.itertuples()])
null = (surr > 1.64).mean(0) # proportion of channels per surrogate (O4)
prop = sig.mean()
fig2b_rows.append({"reading": rd, "area": a, "channels": len(d), "significant": int(sig.sum()), "paper": paper_counts[a],
"within_10pct": abs(sig.sum() - paper_counts[a]) <= 0.1 * paper_counts[a],
"proportion": prop, "null_mean": null.mean(), "null_p99": np.percentile(null, 99),
"p": (np.sum(null >= prop) + 1) / (len(null) + 1)})
fig2b = pd.DataFrame(fig2b_rows)
print(fig2b.round(4).to_string(index=False))
sig_a = ch[(ch.reading == "a") & (ch.z_thetaHG > 1.64)]
print("\nsignificant hippocampal channels per patient (reading a):", sig_a[sig_a.area == "hippocampus"].groupby("subject").size().to_dict())
for thr in [1.28, 1.64, 1.96, 2.33]:
print(f"threshold z > {thr}:", ch[(ch.reading == 'a') & (ch.z_thetaHG > thr)].area.value_counts().reindex(sbcat.AREAS).to_dict())
reading area channels significant paper within_10pct proportion null_mean null_p99 p
a hippocampus 586 144 137 True 0.2457 0.0004 0.0034 0.005
a amygdala 421 130 130 True 0.3088 0.0004 0.0024 0.005
a preSMA 283 17 5 False 0.0601 0.0008 0.0071 0.005
a dACC 325 14 14 True 0.0431 0.0003 0.0031 0.005
a vmPFC 307 41 40 True 0.1336 0.0004 0.0033 0.005
b hippocampus 586 125 137 True 0.2133 0.0004 0.0034 0.005
b amygdala 417 112 130 False 0.2686 0.0006 0.0024 0.005
b preSMA 283 12 5 False 0.0424 0.0007 0.0071 0.005
b dACC 325 8 14 False 0.0246 0.0003 0.0031 0.005
b vmPFC 307 27 40 False 0.0879 0.0003 0.0033 0.005
significant hippocampal channels per patient (reading a): {'1': 9, '11': 1, '12': 2, '17': 4, '18': 3, '19': 3, '21': 4, '22': 29, '23': 1, '26': 10, '27': 1, '28': 8, '29': 10, '3': 1, '30': 7, '31': 11, '32': 4, '33': 8, '34': 2, '36': 6, '4': 8, '5': 2, '9': 10}
threshold z > 1.28: {'hippocampus': 162, 'amygdala': 156, 'preSMA': 24, 'dACC': 23, 'vmPFC': 47}
threshold z > 1.64: {'hippocampus': 144, 'amygdala': 130, 'preSMA': 17, 'dACC': 14, 'vmPFC': 41}
threshold z > 1.96: {'hippocampus': 126, 'amygdala': 116, 'preSMA': 14, 'dACC': 7, 'vmPFC': 33}
threshold z > 2.33: {'hippocampus': 112, 'amygdala': 102, 'preSMA': 11, 'dACC': 2, 'vmPFC': 25}
Result: partially reproduced. Under reading (a), 144 hippocampal (paper 137), 130 amygdala (130), 41 vmPFC (40) and 14 dACC (14) channels are significant, all within tolerance, but the pre-SMA has 17 against 5. Every area's proportion is above the 99th percentile of its null (p = 0.005, the smallest possible with 200 surrogates), the pre-SMA and dACC too. Under reading (b), the amygdala (112), dACC (8), vmPFC (27) and pre-SMA (12) counts are out of tolerance.
fig2c-hippocampus-load and fig2c-no-load-effect-other-areas¶
In the 137 hippocampal channels with significant PAC, theta-high-gamma PAC was weaker in load 3 than in load 1 trials.
Theta-high-gamma PAC did not differ between loads in the amygdala, vmPFC, pre-SMA or dACC channels with significant PAC.
rng = np.random.default_rng(301)
fig2c_rows = []
for rd in READINGS:
for a in sbcat.AREAS:
d = ch[(ch.reading == rd) & (ch.area == a) & (ch.z_thetaHG > 1.64)]
if len(d) < 2:
fig2c_rows.append({"reading": rd, "area": a, "n": len(d)})
continue
t = sbcat.perm_ttest_paired(d.z3_thetaHG, d.z1_thetaHG, rng)
fig2c_rows.append({"reading": rd, "area": a, "n": t["n"], "load1": t["mean_y"], "load3": t["mean_x"],
"t(3-1)": t["t"], "p": t["p"], "channels_load1_higher": int((d.z1_thetaHG > d.z3_thetaHG).sum())})
fig2c = pd.DataFrame(fig2c_rows)
fig2c["p_fdr"] = np.nan
for rd in READINGS:
m = (fig2c.reading == rd) & fig2c.p.notna()
fig2c.loc[m, "p_fdr"] = sbcat.fdr_bh(fig2c.loc[m, "p"])
print(fig2c.round(4).to_string(index=False))
d = ch[(ch.reading == "a") & (ch.area == "hippocampus") & (ch.z_thetaHG > 1.64)]
pp = d.assign(diff=d.z1_thetaHG - d.z3_thetaHG).groupby("subject")["diff"].mean()
print("\nhippocampus, reading a: patients with load 1 above load 3 (mean over their channels):", int((pp > 0).sum()), "of", len(pp))
print("range of the channel z values (reading a, hippocampus): load 1 %.2f to %.2f, load 3 %.2f to %.2f" % (d.z1_thetaHG.min(), d.z1_thetaHG.max(), d.z3_thetaHG.min(), d.z3_thetaHG.max()))
reading area n load1 load3 t(3-1) p channels_load1_higher p_fdr
a hippocampus 144 4.6203 3.7952 -4.1598 0.0002 99 0.0010
a amygdala 130 3.8586 4.2040 2.2201 0.0285 54 0.0475
a preSMA 17 1.8072 6.8738 3.1053 0.0007 3 0.0017
a dACC 14 1.3177 0.9822 -1.0976 0.3003 8 0.3753
a vmPFC 41 2.0866 2.1040 0.0828 0.9346 22 0.9346
b hippocampus 125 4.1373 3.4582 -3.0478 0.0031 76 0.0140
b amygdala 112 3.4510 3.7812 1.8286 0.0690 46 0.0862
b preSMA 12 2.2441 8.7685 2.6825 0.0084 2 0.0140
b dACC 8 1.8761 0.7531 -4.5003 0.0077 8 0.0140
b vmPFC 27 2.2654 2.2427 -0.0803 0.9411 14 0.9411
hippocampus, reading a: patients with load 1 above load 3 (mean over their channels): 18 of 23
range of the channel z values (reading a, hippocampus): load 1 -0.04 to 43.04, load 3 -0.09 to 41.96
print(fig2c[fig2c.area != "hippocampus"].round(4).to_string(index=False))
print("paper: amygdala t129 = 1.43, p = 0.38; vmPFC t39 = 0.16, p = 0.87; pre-SMA t4 = -0.16, p = 0.87; dACC t13 = -0.82, p = 0.73")
for rd in READINGS:
d = ch[(ch.reading == rd) & (ch.area == "preSMA") & (ch.z_thetaHG > 1.64)]
print(f"\npre-SMA PAC channels, reading {rd} (z of theta-high-gamma MI per load):")
print(d[["subject", "session", "channel", "z1_thetaHG", "z3_thetaHG"]].assign(diff=d.z3_thetaHG - d.z1_thetaHG).round(2).to_string(index=False))
r_ = sbcat.perm_ttest_paired(d.z3_thetaHG, d.z1_thetaHG, np.random.default_rng(310))
w = stats.wilcoxon(d.z3_thetaHG, d.z1_thetaHG)
top = d.assign(diff=d.z3_thetaHG - d.z1_thetaHG).sort_values("diff").iloc[:-3]
r3 = sbcat.perm_ttest_paired(top.z3_thetaHG, top.z1_thetaHG, np.random.default_rng(311))
print(f"load 3 higher in {int((d.z3_thetaHG > d.z1_thetaHG).sum())} of {len(d)}; paired permutation t-test p = {r_['p']:.4f}; "
f"Wilcoxon signed-rank p = {w.pvalue:.4f}; without the 3 channels of largest load 3 excess, t = {r3['t']:.2f}, p = {r3['p']:.4f}")
reading area n load1 load3 t(3-1) p channels_load1_higher p_fdr
a amygdala 130 3.8586 4.2040 2.2201 0.0285 54 0.0475
a preSMA 17 1.8072 6.8738 3.1053 0.0007 3 0.0017
a dACC 14 1.3177 0.9822 -1.0976 0.3003 8 0.3753
a vmPFC 41 2.0866 2.1040 0.0828 0.9346 22 0.9346
b amygdala 112 3.4510 3.7812 1.8286 0.0690 46 0.0862
b preSMA 12 2.2441 8.7685 2.6825 0.0084 2 0.0140
b dACC 8 1.8761 0.7531 -4.5003 0.0077 8 0.0140
b vmPFC 27 2.2654 2.2427 -0.0803 0.9411 14 0.9411
paper: amygdala t129 = 1.43, p = 0.38; vmPFC t39 = 0.16, p = 0.87; pre-SMA t4 = -0.16, p = 0.87; dACC t13 = -0.82, p = 0.73
pre-SMA PAC channels, reading a (z of theta-high-gamma MI per load):
subject session channel z1_thetaHG z3_thetaHG diff
10 2 21 0.82 1.17 0.34
11 3 7 -0.26 2.78 3.04
11 3 8 -0.26 2.25 2.51
11 3 9 0.55 3.04 2.49
11 3 10 -0.32 3.97 4.29
11 3 11 0.68 2.38 1.71
11 3 12 0.25 4.97 4.72
11 3 13 1.63 5.54 3.91
15 8 9 1.26 0.73 -0.53
15 8 11 0.78 1.19 0.40
15 8 12 0.75 1.25 0.50
17 10 0 1.84 0.67 -1.17
4 37 37 6.83 19.27 12.43
4 37 38 6.05 19.10 13.05
4 37 39 5.80 22.71 16.91
4 37 40 3.52 25.12 21.61
6 40 14 0.80 0.72 -0.08
load 3 higher in 14 of 17; paired permutation t-test p = 0.0010; Wilcoxon signed-rank p = 0.0011; without the 3 channels of largest load 3 excess, t = 2.70, p = 0.0045
pre-SMA PAC channels, reading b (z of theta-high-gamma MI per load):
subject session channel z1_thetaHG z3_thetaHG diff
10 2 21 0.84 0.88 0.03
11 3 7 -0.29 3.73 4.02
11 3 8 -0.05 1.64 1.69
11 3 9 0.40 1.26 0.86
11 3 10 -0.37 3.67 4.04
11 3 13 0.98 3.73 2.75
17 10 0 2.60 0.53 -2.07
4 37 37 5.61 18.27 12.66
4 37 38 6.06 18.84 12.77
4 37 39 6.37 22.64 16.27
4 37 40 3.63 29.24 25.60
6 40 14 1.13 0.79 -0.34
load 3 higher in 10 of 12; paired permutation t-test p = 0.0069; Wilcoxon signed-rank p = 0.0093; without the 3 channels of largest load 3 excess, t = 1.84, p = 0.0593
Result: fig2c-hippocampus-load reproduced. Hippocampal PAC is lower in load 3 than in load 1 under both readings (a: 144 channels, t = -4.16, FDR p = 0.0010; b: 125 channels, t = -3.05, FDR p = 0.014).
Result: fig2c-no-load-effect-other-areas not reproduced. Under reading (a) the pre-SMA (17 channels, load 3 higher, FDR p = 0.0017) and the amygdala (130 channels, load 3 higher, FDR p = 0.0475) differ between loads; under (b) the pre-SMA (FDR p = 0.014) and the dACC (8 channels, load 1 higher, FDR p = 0.014) do. The vmPFC does not differ under either (p = 0.93, 0.94). The pre-SMA channels come from six patients; four channels of sub-4 have z values of 19 to 25 in load 3 against 3.5 to 6.8 in load 1, and seven of sub-11 are also higher in load 3. Load 3 is higher in 14 of the 17 channels (Wilcoxon p = 0.0011), and without the three channels of largest excess the test still gives p = 0.0045 under reading (a); under (b) it gives p = 0.059.
edf2b-hippocampus-load-patients¶
Averaging significant PAC channels within each patient, hippocampal PAC was again stronger in load 1 than in load 3, and no other area showed a load difference.
rng = np.random.default_rng(302)
rows = []
for rd in READINGS:
for a in sbcat.AREAS:
d = ch[(ch.reading == rd) & (ch.area == a) & (ch.z_thetaHG > 1.64)]
pm = d.groupby("subject")[["z1_thetaHG", "z3_thetaHG"]].mean()
if len(pm) < 2:
rows.append({"reading": rd, "area": a, "patients": len(pm)})
continue
t = sbcat.perm_ttest_paired(pm.z1_thetaHG, pm.z3_thetaHG, rng)
rows.append({"reading": rd, "area": a, "patients": len(pm), "load1": t["mean_x"], "load3": t["mean_y"], "t(1-3)": t["t"], "p": t["p"],
"patients_load1_higher": int((pm.z1_thetaHG > pm.z3_thetaHG).sum())})
print(pd.DataFrame(rows).round(4).to_string(index=False))
print("paper: hippocampus n = 23 patients, p = 0.0049")
reading area patients load1 load3 t(1-3) p patients_load1_higher
a hippocampus 23 3.9687 3.2301 2.2616 0.0095 18
a amygdala 27 2.8140 3.1784 -1.6857 0.1021 11
a preSMA 6 1.7115 4.7877 -1.1594 0.3120 2
a dACC 8 0.9797 0.9948 -0.0417 0.9921 3
a vmPFC 12 2.1086 2.1599 -0.1211 0.9233 7
b hippocampus 20 3.6051 3.1090 1.5014 0.1127 14
b amygdala 25 2.6486 3.0270 -1.6992 0.1086 9
b preSMA 5 2.0252 5.4497 -0.9968 0.4411 2
b dACC 4 1.7814 0.5825 3.6894 0.1240 4
b vmPFC 10 1.9590 2.2559 -0.6765 0.5230 4
paper: hippocampus n = 23 patients, p = 0.0049
Result: reproduced, resting on an assumption. Under reading (a), 23 patients have hippocampal PAC channels (paper 23), their means are higher in load 1 (3.97 against 3.23, p = 0.0095; 18 of 23 patients), and no other area differs (p >= 0.10). Under reading (b) the hippocampus gives p = 0.11 over 20 patients. Reading (a) is adopted (header and data check: the paper's own rejected trials on its example channel include every trial (a) removes).
edf2c-low-gamma-no-load-effect¶
Theta-low-gamma (30-55 Hz) PAC did not differ between loads in any region.
rng = np.random.default_rng(303)
rows = []
for rd in READINGS:
for a in sbcat.AREAS:
d = ch[(ch.reading == rd) & (ch.area == a) & (ch.z_thetaLG > 1.64)]
if len(d) < 2:
rows.append({"reading": rd, "area": a, "n": len(d)})
continue
t = sbcat.perm_ttest_paired(d.z3_thetaLG, d.z1_thetaLG, rng)
rows.append({"reading": rd, "area": a, "n": t["n"], "load1": t["mean_y"], "load3": t["mean_x"], "t(3-1)": t["t"], "p": t["p"]})
lg = pd.DataFrame(rows)
lg["p_fdr"] = np.nan
for rd in READINGS:
m = (lg.reading == rd) & lg.p.notna()
lg.loc[m, "p_fdr"] = sbcat.fdr_bh(lg.loc[m, "p"])
print(lg.round(4).to_string(index=False))
print("paper n: 148, 155, 10, 22 and 48 channels (hippocampus, amygdala, pre-SMA, dACC, vmPFC)")
reading area n load1 load3 t(3-1) p p_fdr
a hippocampus 152 4.0861 3.7691 -1.5870 0.1127 0.2817
a amygdala 159 3.8603 3.7910 -0.4083 0.6840 0.6840
a preSMA 22 1.7798 4.6915 2.6656 0.0029 0.0145
a dACC 27 1.4122 1.7603 1.3539 0.1842 0.3070
a vmPFC 48 2.0555 2.2778 0.8793 0.3937 0.4921
b hippocampus 117 3.9291 3.6972 -0.9395 0.3408 0.5679
b amygdala 134 3.4373 3.4556 0.1121 0.9104 0.9104
b preSMA 13 2.3657 7.0907 2.6703 0.0033 0.0165
b dACC 18 1.6698 1.2210 -1.1637 0.2601 0.5679
b vmPFC 37 2.0083 2.1568 0.4757 0.6480 0.8100
paper n: 148, 155, 10, 22 and 48 channels (hippocampus, amygdala, pre-SMA, dACC, vmPFC)
Result: not reproduced. The pre-SMA's low-gamma PAC differs between loads under both readings (a: 22 channels, p = 0.0029; b: 13 channels, p = 0.0033, load 3 higher); no other area does (p >= 0.11). The counts of selected channels are 152, 159, 22, 27 and 48 under (a), against the paper's 148, 155, 10, 22 and 48.
tables1-slow-fast-theta¶
Split into slow (2-5 Hz) and fast (5-9 Hz) theta, PAC differed between loads only in the hippocampus and only with high gamma, most strongly for fast theta (t116 = -3.05, p = 0.0015; slow theta, t147 = -2.31, p = 0.0163).
Readings: (a) phase centres 2 and 4 Hz (slow) and 6 and 8 Hz (fast) of the comodulogram; (b) bands filtered at 2-5 and 5-9 Hz. Each under both artefact readings.
rng = np.random.default_rng(304)
rows = []
for rd in READINGS:
for band_reading in ["A", "B"]:
for th in ["slow", "fast"]:
for g in ["HG", "LG"]:
if band_reading == "A":
sel_col, c1, c3 = f"z_{th}A_{g}", f"z1_{th}A_{g}", f"z3_{th}A_{g}"
else:
sel_col, c1, c3 = f"zb_{th}_{g}_all", f"zb_{th}_{g}_l1", f"zb_{th}_{g}_l3"
for a in sbcat.AREAS:
d = ch[(ch.reading == rd) & (ch.area == a) & (ch[sel_col] > 1.64)]
if len(d) < 2:
rows.append({"artefact": rd, "bands": band_reading, "theta": th, "gamma": g, "area": a, "n": len(d)})
continue
t = sbcat.perm_ttest_paired(d[c3], d[c1], rng)
rows.append({"artefact": rd, "bands": band_reading, "theta": th, "gamma": g, "area": a, "n": t["n"],
"load1": t["mean_y"], "load3": t["mean_x"], "t(3-1)": t["t"], "p": t["p"]})
ts1 = pd.DataFrame(rows)
print(ts1.round(4).to_string(index=False))
summ = ts1.assign(sig=ts1.p < 0.05).groupby(["artefact", "bands"]).apply(
lambda g: pd.Series({"tests": int(g.p.notna().sum()), "significant": int(g.sig.sum()),
"significant_tests": "; ".join(f"{r.area} {r.theta} {r.gamma} t={r['t(3-1)']:.2f}" for _, r in g[g.sig].iterrows())}))
print(summ.to_string())
artefact bands theta gamma area n load1 load3 t(3-1) p
a A slow HG hippocampus 165 3.7467 3.4088 -2.4262 0.0169
a A slow HG amygdala 146 3.6087 3.7797 0.9749 0.3308
a A slow HG preSMA 47 1.8879 3.8419 2.8906 0.0009
a A slow HG dACC 24 1.7431 1.2795 -1.7265 0.0946
a A slow HG vmPFC 55 3.4108 3.2350 -0.8639 0.3900
a A slow LG hippocampus 161 3.9069 3.3689 -3.2119 0.0019
a A slow LG amygdala 162 3.5600 3.5795 0.1154 0.9081
a A slow LG preSMA 38 2.1721 3.8424 2.6204 0.0020
a A slow LG dACC 41 1.5722 2.1777 2.3768 0.0138
a A slow LG vmPFC 73 2.4291 2.6320 1.0579 0.2946
a A fast HG hippocampus 120 4.2900 3.4777 -3.4422 0.0010
a A fast HG amygdala 93 3.7751 3.9453 0.7144 0.4790
a A fast HG preSMA 9 2.6354 8.6338 2.4548 0.0469
a A fast HG dACC 11 1.3229 1.8863 1.5785 0.1326
a A fast HG vmPFC 30 1.9295 1.6585 -0.9639 0.3348
a A fast LG hippocampus 131 3.9171 3.9173 0.0009 0.9996
a A fast LG amygdala 143 3.6546 3.4361 -1.3237 0.1877
a A fast LG preSMA 14 2.3000 5.6921 2.2854 0.0384
a A fast LG dACC 21 1.5283 1.3624 -0.5837 0.5660
a A fast LG vmPFC 35 1.6837 1.7191 0.1152 0.9064
a B slow HG hippocampus 190 4.3955 4.0575 -2.2055 0.0280
a B slow HG amygdala 172 4.0442 4.2214 0.9193 0.3533
a B slow HG preSMA 46 2.3221 4.2663 2.7587 0.0045
a B slow HG dACC 26 2.1693 1.5868 -1.5599 0.1339
a B slow HG vmPFC 70 3.3035 3.3453 0.1826 0.8541
a B slow LG hippocampus 182 4.6042 4.3591 -1.2778 0.2080
a B slow LG amygdala 187 3.9588 4.0406 0.4230 0.6714
a B slow LG preSMA 39 2.4833 4.3512 2.6546 0.0061
a B slow LG dACC 50 1.9379 2.4940 1.9796 0.0523
a B slow LG vmPFC 89 2.5772 2.9088 1.6230 0.1148
a B fast HG hippocampus 152 4.9972 4.1736 -3.2656 0.0009
a B fast HG amygdala 118 4.1430 4.4140 1.0076 0.3110
a B fast HG preSMA 21 2.0255 5.3526 2.1789 0.0329
a B fast HG dACC 28 1.1934 1.5554 1.5995 0.1223
a B fast HG vmPFC 38 2.3468 2.0790 -0.7867 0.4322
a B fast LG hippocampus 156 4.4528 4.3714 -0.3171 0.7562
a B fast LG amygdala 161 4.0632 4.0178 -0.2293 0.8220
a B fast LG preSMA 22 2.5131 4.7692 1.6716 0.1123
a B fast LG dACC 29 1.2879 1.3164 0.0907 0.9278
a B fast LG vmPFC 52 1.9309 2.0002 0.2630 0.7988
b A slow HG hippocampus 146 3.2966 2.9252 -2.4597 0.0163
b A slow HG amygdala 120 3.2496 3.3778 0.7066 0.4904
b A slow HG preSMA 41 1.7484 3.9207 2.5564 0.0031
b A slow HG dACC 16 2.0117 1.1565 -2.4988 0.0286
b A slow HG vmPFC 46 3.0804 2.9688 -0.4694 0.6459
b A slow LG hippocampus 133 3.3423 3.0924 -1.4604 0.1475
b A slow LG amygdala 140 3.0857 3.0983 0.0673 0.9474
b A slow LG preSMA 32 2.0566 3.9793 2.2724 0.0209
b A slow LG dACC 34 1.4697 1.7693 1.1998 0.2408
b A slow LG vmPFC 56 2.3671 2.4598 0.4600 0.6524
b A fast HG hippocampus 103 3.8405 3.0784 -3.0260 0.0019
b A fast HG amygdala 82 3.2295 3.3931 0.6748 0.5041
b A fast HG preSMA 7 3.0115 10.3351 2.2134 0.1302
b A fast HG dACC 8 1.5667 1.7839 0.5169 0.6144
b A fast HG vmPFC 19 1.8833 1.8572 -0.0653 0.9519
b A fast LG hippocampus 110 3.6401 3.4325 -0.8036 0.4342
b A fast LG amygdala 114 3.4255 3.1973 -1.1487 0.2518
b A fast LG preSMA 13 2.0915 5.8178 2.5283 0.0182
b A fast LG dACC 21 1.3162 1.1110 -0.7389 0.4726
b A fast LG vmPFC 23 1.5483 1.7845 0.6128 0.5438
b B slow HG hippocampus 160 3.9950 3.5957 -2.2178 0.0276
b B slow HG amygdala 146 3.6351 3.7277 0.4782 0.6427
b B slow HG preSMA 47 1.9127 3.8406 2.5064 0.0101
b B slow HG dACC 21 2.0595 1.5965 -1.0065 0.3307
b B slow HG vmPFC 57 3.0344 3.2149 0.7356 0.4646
b B slow LG hippocampus 160 3.8897 3.9003 0.0537 0.9589
b B slow LG amygdala 159 3.5963 3.6681 0.3321 0.7421
b B slow LG preSMA 32 2.3763 4.5965 2.4492 0.0116
b B slow LG dACC 40 1.8351 1.9926 0.5236 0.6068
b B slow LG vmPFC 68 2.5563 2.6218 0.3357 0.7359
b B fast HG hippocampus 137 4.3840 3.6901 -2.5048 0.0118
b B fast HG amygdala 103 3.6853 3.8947 0.7632 0.4551
b B fast HG preSMA 10 3.0146 10.2855 2.4272 0.0283
b B fast HG dACC 15 1.5445 1.5256 -0.0580 0.9544
b B fast HG vmPFC 30 2.1165 2.0149 -0.2640 0.8002
b B fast LG hippocampus 142 3.8603 3.6368 -0.8161 0.4289
b B fast LG amygdala 130 3.6921 3.7606 0.3096 0.7592
b B fast LG preSMA 22 2.1182 4.8561 2.1769 0.0430
b B fast LG dACC 21 1.4821 1.2115 -0.7519 0.4657
b B fast LG vmPFC 31 2.0910 2.0517 -0.0990 0.9306
tests significant significant_tests
artefact bands
a A 20 8 hippocampus slow HG t=-2.43; preSMA slow HG t=2.89; hippocampus slow LG t=-3.21; preSMA slow LG t=2.62; dACC slow LG t=2.38; hippocampus fast HG t=-3.44; preSMA fast HG t=2.45; preSMA fast LG t=2.29
B 20 5 hippocampus slow HG t=-2.21; preSMA slow HG t=2.76; preSMA slow LG t=2.65; hippocampus fast HG t=-3.27; preSMA fast HG t=2.18
b A 20 6 hippocampus slow HG t=-2.46; preSMA slow HG t=2.56; dACC slow HG t=-2.50; preSMA slow LG t=2.27; hippocampus fast HG t=-3.03; preSMA fast LG t=2.53
B 20 6 hippocampus slow HG t=-2.22; preSMA slow HG t=2.51; preSMA slow LG t=2.45; hippocampus fast HG t=-2.50; preSMA fast HG t=2.43; preSMA fast LG t=2.18
Result: not reproduced. Under every combination of readings both hippocampal high-gamma tests give p < 0.05 with load 1 higher, but other tests do too: 8, 5, 6 and 6 of 20 tests are significant, always including pre-SMA tests with load 3 higher, and under (a, A) also the hippocampal slow theta to low gamma test.
fig2f-rt-single-trial-pac¶
In the hippocampus, stronger single-trial theta-high-gamma PAC went with faster RTs, while no such relation held in the amygdala or vmPFC.
Single-trial PAC (O5): (a) the MI of each trial's 0 to 2.5 s window, averaged over the theta and high-gamma pairs, raw; (b) the same z-scored across trials within each channel. Correct, clean trials of the significant channels; RT ~ 1 + PAC * load + (1 | patient) + (1 | channel nested in patient), normal response (Supplementary Table S2), fitted by REML with statsmodels MixedLM. Load is coded with load 3 as the reference, as Table S2 names the terms "Load1" and "PAC : Load1" and says "Load 1 are tested against load 3 trials"; the PAC term is then the slope in load 3 trials. The slopes in load 1 trials and averaged over loads are also printed. Table S2 gives no unit or transform for RT, and its Load1 term (-0.48) cannot be a difference of raw RTs in s or ms, so RT is fitted in s, as log10 RT, and z-scored over each session's correct trials.
def slopes(m, term, ind="load1"):
"""The term's slope in load 3 trials, in load 1 trials and averaged over the loads, with Wald z and p.
load1 is an indicator of load 1 trials, so load 3 is the reference: the paper's coding, as Tables
S2 to S4 name the terms "Load1" and "PAC : Load1" and say "Load 1 are tested against load 3 trials".
The model's own PAC (or SC) term is then the load 3 slope."""
V = m.cov_params()
b, i = term, f"{term}:{ind}"
out = {}
for name, w in [("L3", 0.0), ("L1", 1.0), ("mean", 0.5)]:
est = m.params[b] + w * m.params[i]
z = est / np.sqrt(V.loc[b, b] + w * w * V.loc[i, i] + 2 * w * V.loc[b, i])
out.update({f"{name} est": est, f"{name} t": z, f"{name} p": 2 * stats.norm.sf(abs(z))})
out.update({"Load1 est": m.params[ind], "int est": m.params[i], "int p": m.pvalues[i], "converged": m.converged})
return out
def rt_frame(rd, a):
"""Correct, clean trials of the area's significant PAC channels under artefact reading rd. RT in s,
log10 RT, and RT z-scored over the session's correct trials (the paper gives no unit or transform)."""
d = ch[(ch.reading == rd) & (ch.area == a) & (ch.z_thetaHG > 1.64)]
parts = []
for r in d.itertuples():
c = chan_index[(r.session, r.channel)]
tr = S[r.session]["trials"]
rt_c = tr.rt[tr.correct]
idx = c["trial_index"]
keep = c["masks"][rd][idx]
pac = c["st_thetaHG"][keep]
rt = tr.rt.to_numpy()[idx][keep]
parts.append(pd.DataFrame({"rt": rt, "log_rt": np.log10(rt), "z_rt": (rt - rt_c.mean()) / rt_c.std(),
"load1": (tr.load.to_numpy()[idx][keep] == 1).astype(int),
"pac_a": pac, "pac_b": (pac - np.nanmean(pac)) / np.nanstd(pac),
"patient": r.subject, "channel": f"{r.session}_{r.channel}"}))
return d.shape[0], pd.concat(parts).dropna().reset_index(drop=True)
glm_rows, frames = [], {}
for rd in READINGS:
for a in ["hippocampus", "amygdala", "vmPFC"]:
n_ch, df = rt_frame(rd, a)
frames[(rd, a)] = df
for resp in ["rt", "log_rt", "z_rt"]:
for o5 in ["a", "b"]:
base = {"artefact": rd, "O5": o5, "response": resp, "area": a, "channels": n_ch, "trials": len(df)}
try:
m = smf.mixedlm(f"{resp} ~ pac_{o5} * load1", df, groups="patient",
vc_formula={"channel": "0 + C(channel)"}).fit(reml=True)
glm_rows.append({**base, "intercept": m.params["Intercept"], **slopes(m, f"pac_{o5}")})
except Exception as e: # report a failed fit
glm_rows.append({**base, "error": str(e)[:80]})
glm = pd.DataFrame(glm_rows)
fmt = lambda v: f"{v:.4g}"
print("PAC slope in load 3 trials (L3, the paper's PAC term), in load 1 trials (L1) and averaged (mean);")
print("Load1 est is the load 1 minus load 3 difference at PAC = 0, int the PAC x Load1 term; intercept is load 3 at PAC = 0.")
print(glm.to_string(index=False, float_format=fmt))
print("paper (Table S2): hippocampus PAC -49.07, t = -4.01, p = 6e-5, Load1 -0.48, PAC:Load1 26.30 (p = 0.11), intercept 0.38;")
print("amygdala PAC -5.15, p = 0.48, Load1 -0.48; vmPFC PAC 24.51, p = 0.24, Load1 -0.34")
dh = frames[("a", "hippocampus")]
print("RT of the hippocampal trials under artefact reading (a): %d trials, median %.2f s, max %.1f s, %d above 10 s"
% (len(dh), dh.rt.median(), dh.rt.max(), (dh.rt >= 10).sum()))
for resp in ["rt", "log_rt", "z_rt"]:
print(f" {resp}: load 1 mean {dh[dh.load1 == 1][resp].mean():.3f}, load 3 mean {dh[dh.load1 == 0][resp].mean():.3f}")
# Variant on artefact reading (a): RTs below 10 s only (the paper states no trimming)
var_rows = []
for a in ["hippocampus", "amygdala", "vmPFC"]:
dv = frames[("a", a)]
dv = dv[dv.rt < 10].reset_index(drop=True)
for o5 in ["a", "b"]:
m = smf.mixedlm(f"rt ~ pac_{o5} * load1", dv, groups="patient", vc_formula={"channel": "0 + C(channel)"}).fit(reml=True)
var_rows.append({"variant": "RT < 10 s", "O5": o5, "area": a, "trials": len(dv), **slopes(m, f"pac_{o5}")})
print("\nartefact reading (a), RT below 10 s:")
print(pd.DataFrame(var_rows).to_string(index=False, float_format=fmt))
PAC slope in load 3 trials (L3, the paper's PAC term), in load 1 trials (L1) and averaged (mean);
Load1 est is the load 1 minus load 3 difference at PAC = 0, int the PAC x Load1 term; intercept is load 3 at PAC = 0.
artefact O5 response area channels trials intercept L3 est L3 t L3 p L1 est L1 t L1 p mean est mean t mean p Load1 est int est int p converged
a a rt hippocampus 144 16103 1.697 -53.02 -2.429 0.01513 23.77 1.149 0.2507 -14.63 -0.9319 0.3514 -0.285 76.79 0.007467 True
a b rt hippocampus 144 16103 1.49 -0.03468 -1.651 0.0987 0.009429 0.4657 0.6414 -0.01262 -0.8672 0.3858 0.01588 0.04411 0.1313 True
a a log_rt hippocampus 144 16103 0.1588 -10.5 -4.988 6.092e-07 -3.025 -1.515 0.1297 -6.763 -4.458 8.281e-06 -0.0736 7.477 0.006781 True
a b log_rt hippocampus 144 16103 0.1177 -0.005735 -2.838 0.004536 -0.005292 -2.716 0.006605 -0.005513 -3.937 8.261e-05 -0.04433 0.0004436 0.8747 True
a a z_rt hippocampus 144 16103 0.2791 -49.39 -4.189 2.799e-05 -13.11 -1.172 0.241 -31.25 -3.675 0.0002374 -0.315 36.28 0.01888 True
a b z_rt hippocampus 144 16103 0.08562 -0.02462 -2.178 0.0294 -0.02552 -2.342 0.01919 -0.02507 -3.2 0.001373 -0.1729 -0.0009 0.9544 True
a a rt amygdala 130 14469 1.683 -39.96 -2.206 0.0274 12.41 0.6099 0.5419 -13.77 -0.9238 0.3556 -0.1756 52.36 0.0318 True
a b rt amygdala 130 14469 1.515 0.001402 0.05383 0.9571 -0.0523 -1.996 0.04599 -0.02545 -1.381 0.1673 0.04259 -0.0537 0.1471 True
a a log_rt amygdala 130 14469 0.1402 -3.371 -2.093 0.03633 -0.3424 -0.1892 0.85 -1.857 -1.391 0.1643 -0.07287 3.029 0.1584 True
a b log_rt amygdala 130 14469 0.1261 -0.002785 -1.215 0.2244 -0.004207 -1.824 0.06816 -0.003496 -2.156 0.03109 -0.06024 -0.001422 0.6628 True
a a z_rt amygdala 130 14469 -1.695e-11 -19.57 -2.428 0.01518 -8.408 -0.9274 0.3537 -13.99 -2.087 0.03686 -0.2339 11.17 0.2971 True
a b z_rt amygdala 130 14469 0.08823 -0.01462 -1.279 0.2009 -0.01818 -1.581 0.1139 -0.0164 -2.029 0.04251 -0.1871 -0.003562 0.8266 True
a a rt vmPFC 41 4774 1.578 7.062 0.2521 0.8009 2.986 0.09942 0.9208 5.024 0.239 0.8111 -0.08777 -4.076 0.919 True
a b rt vmPFC 41 4774 1.605 -0.01137 -0.4282 0.6685 0.02697 1.029 0.3034 0.007803 0.4188 0.6754 -0.103 0.03834 0.3047 True
a a log_rt vmPFC 41 4774 0.1466 2.1 0.5487 0.5832 -1.536 -0.3743 0.7082 0.2821 0.09818 0.9218 -0.05757 -3.637 0.5065 True
a b log_rt vmPFC 41 4774 0.1545 -0.000597 -0.1645 0.8693 0.002026 0.5654 0.5718 0.0007143 0.2805 0.7791 -0.07125 0.002623 0.6075 True
a a z_rt vmPFC 41 4774 0.1641 -2.071 -0.09396 0.9251 7.01 0.2966 0.7668 2.47 0.1492 0.8814 -0.3265 9.081 0.7732 True
a b z_rt vmPFC 41 4774 0.1563 -0.01165 -0.5578 0.577 0.01293 0.6274 0.5304 0.0006417 0.04379 0.9651 -0.2923 0.02458 0.4028 True
b a rt hippocampus 125 12181 1.644 -58.59 -3.975 7.033e-05 -18.51 -1.302 0.1928 -38.55 -3.583 0.0003397 -0.195 40.08 0.03857 True
b b rt hippocampus 125 12181 1.413 -0.04016 -2.782 0.005406 -0.01709 -1.274 0.2028 -0.02862 -2.911 0.003603 -0.03784 0.02307 0.2429 True
b a log_rt hippocampus 125 12181 0.1511 -11.74 -5.235 1.647e-07 -5.929 -2.741 0.006122 -8.837 -5.391 7.01e-08 -0.06476 5.815 0.04832 True
b b log_rt hippocampus 125 12181 0.1049 -0.006701 -3.053 0.002263 -0.00668 -3.275 0.001058 -0.006691 -4.476 7.613e-06 -0.04197 2.1e-05 0.9944 True
b a z_rt hippocampus 125 12181 0.22 -46.95 -4.024 5.713e-05 -21.43 -1.905 0.05682 -34.19 -4.228 2.362e-05 -0.2599 25.52 0.1161 True
b b z_rt hippocampus 125 12181 0.03539 -0.02609 -2.153 0.03132 -0.03053 -2.71 0.006733 -0.02831 -3.422 0.0006219 -0.1592 -0.004441 0.7884 True
b a rt amygdala 112 10750 1.486 -2.226 -0.1801 0.8571 -3.665 -0.2888 0.7727 -2.946 -0.3031 0.7618 -0.08975 -1.438 0.9275 True
b b rt amygdala 112 10750 1.476 -0.00481 -0.267 0.7895 -0.02375 -1.387 0.1654 -0.01428 -1.153 0.2491 -0.09595 -0.01894 0.4475 True
b a log_rt amygdala 112 10750 0.1164 0.4187 0.2443 0.807 -1.069 -0.6063 0.5443 -0.325 -0.2399 0.8104 -0.05214 -1.488 0.4948 True
b b log_rt amygdala 112 10750 0.1182 -0.001141 -0.4595 0.6459 -0.003685 -1.562 0.1183 -0.002413 -1.414 0.1575 -0.05846 -0.002544 0.459 True
b a z_rt amygdala 112 10750 0.02387 -3.656 -0.5297 0.5963 -12.01 -1.731 0.08341 -7.834 -1.601 0.1094 -0.1322 -8.357 0.3932 True
b b z_rt amygdala 112 10750 0.008418 -0.007932 -0.7085 0.4787 -0.01756 -1.65 0.09894 -0.01275 -1.65 0.09888 -0.1674 -0.009627 0.5331 True
b a rt vmPFC 27 2667 1.714 -11.25 -0.388 0.698 -19.09 -0.8251 0.4093 -15.17 -0.8003 0.4235 -0.1336 -7.837 0.8289 True
b b rt vmPFC 27 2667 1.672 -0.03065 -1.134 0.257 -0.006331 -0.2604 0.7946 -0.01849 -1.018 0.3087 -0.1635 0.02432 0.5041 True
b a log_rt vmPFC 27 2667 0.1572 0.9489 0.1917 0.8479 -2.643 -0.6695 0.5032 -0.8472 -0.2619 0.7934 -0.05907 -3.592 0.5616 True
b b log_rt vmPFC 27 2667 0.1609 -0.001465 -0.3175 0.7509 0.001345 0.3238 0.7461 -6.047e-05 -0.0195 0.9844 -0.07281 0.00281 0.6511 True
b a z_rt vmPFC 27 2667 0.0761 1.253 0.04982 0.9603 -12.45 -0.6205 0.5349 -5.599 -0.3405 0.7335 -0.2321 -13.7 0.663 True
b b z_rt vmPFC 27 2667 0.08106 -0.00443 -0.1889 0.8502 -0.002699 -0.1279 0.8982 -0.003564 -0.2262 0.821 -0.2845 0.001731 0.9563 True
paper (Table S2): hippocampus PAC -49.07, t = -4.01, p = 6e-5, Load1 -0.48, PAC:Load1 26.30 (p = 0.11), intercept 0.38;
amygdala PAC -5.15, p = 0.48, Load1 -0.48; vmPFC PAC 24.51, p = 0.24, Load1 -0.34
RT of the hippocampal trials under artefact reading (a): 16103 trials, median 1.15 s, max 56.3 s, 87 above 10 s
rt: load 1 mean 1.490, load 3 mean 1.492
log_rt: load 1 mean 0.070, load 3 mean 0.117
z_rt: load 1 mean -0.088, load 3 mean 0.085
artefact reading (a), RT below 10 s: variant O5 area trials L3 est L3 t L3 p L1 est L1 t L1 p mean est mean t mean p Load1 est int est int p converged RT < 10 s a hippocampus 16016 -49.23 -5.135 2.822e-07 -12.99 -1.422 0.1551 -31.11 -4.494 6.99e-06 -0.2736 36.25 0.004047 True RT < 10 s b hippocampus 16016 -0.03224 -3.505 0.0004561 -0.01494 -1.676 0.0937 -0.02359 -3.692 0.0002229 -0.1318 0.01731 0.1774 True RT < 10 s a amygdala 14375 -11.93 -1.66 0.09699 0.1017 0.01259 0.99 -5.916 -0.9948 0.3198 -0.2283 12.03 0.2107 True RT < 10 s b amygdala 14375 -0.01699 -1.659 0.09715 -0.005758 -0.5573 0.5773 -0.01137 -1.568 0.1169 -0.1782 0.01123 0.4415 True RT < 10 s a vmPFC 4733 0.1917 0.01134 0.9909 -2.731 -0.1505 0.8804 -1.27 -0.1 0.9203 -0.2075 -2.923 0.9039 True RT < 10 s b vmPFC 4733 -0.008106 -0.5052 0.6134 0.004782 0.3014 0.7631 -0.001662 -0.1475 0.8827 -0.2185 0.01289 0.5684 True
Result: blocked (methods underspecified), resting on artefact reading (a). With load 3 as the reference level, as Table S2 codes load, the PAC term is the slope in load 3 trials. With RT in s under artefact reading (a), neither single-trial reading meets the criterion: under O5 (a) the hippocampal term is negative (-53.0, t = -2.43, p = 0.015) but so is the amygdala term (p = 0.027), and under O5 (b) the hippocampal term has p = 0.099. RT in s does not fit Table S2: the intercept is 1.70 (O5 a) against the paper's 0.38. With RT z-scored over each session's correct trials the intercept is 0.28 and the Load1 term -0.32 (paper 0.38 and -0.48), and under O5 (a) the hippocampal term is -49.39, t = -4.19, p = 2.8e-5 (paper -49.07, t = -4.01, p = 6e-5), but the amygdala term is also negative (p = 0.015, paper p = 0.48); under O5 (b) the hippocampal term is negative (p = 0.029) and the amygdala (p = 0.20) and vmPFC (p = 0.58) terms are not significant, which meets the criterion. log10 RT gives the same split (O5 a, amygdala p = 0.036; O5 b, hippocampus p = 0.0045, amygdala 0.22, vmPFC 0.87). So on the scale that fits Table S2 the two readings of O5 disagree. With RTs below 10 s (in s), both readings meet the criterion (hippocampus p = 2.8e-7 and 0.00046, amygdala 0.097 and 0.097, vmPFC 0.99 and 0.61). Under artefact reading (b), every response scale and both readings meet it (hippocampus p <= 0.031, amygdala p >= 0.479, vmPFC p >= 0.26). The slope in load 1 trials, which the earlier coding tested, is not significant in the hippocampus with RT in s (p = 0.25 and 0.64). The criterion makes a disagreement between the readings of O5 methods underspecified; this outcome replaced not reproduced after the third verification round. Table S2's PAC estimate (-49.07) can only come from the raw MI, which points to O5 (a), under which the outcome would be not reproduced; this is asked of a reader.
Hippocampal PAC channels: kappa, raw MI, nesting, spectral slope and controls¶
def hpc_sig(rd):
return ch[(ch.reading == rd) & (ch.area == "hippocampus") & (ch.z_thetaHG > 1.64)]
edf3a-kappa-load¶
The distribution of gamma amplitude over theta phase was wider (lower kappa) in load 3 than in load 1 across hippocampal PAC channels.
rng = np.random.default_rng(305)
for rd in READINGS:
d = hpc_sig(rd)
k1 = np.array([sbcat.circ_kappa_weighted(chan_index[(r.session, r.channel)]["readings"][rd]["P_thetaHG_l1"]) for r in d.itertuples()])
k3 = np.array([sbcat.circ_kappa_weighted(chan_index[(r.session, r.channel)]["readings"][rd]["P_thetaHG_l3"]) for r in d.itertuples()])
t = sbcat.perm_ttest_paired(k3, k1, rng)
print(f"reading {rd}: n = {t['n']}, kappa load 1 {k1.mean():.4f}, load 3 {k3.mean():.4f}, t(3-1) = {t['t']:.2f}, p = {t['p']:.4f}; "
f"channels with lower kappa in load 3: {int((k3 < k1).sum())}; kappa range {min(k1.min(), k3.min()):.4f} to {max(k1.max(), k3.max()):.4f}")
print("paper: t(136) = -3.7453, p = 0.0001")
reading a: n = 144, kappa load 1 0.0360, load 3 0.0326, t(3-1) = -4.54, p = 0.0001; channels with lower kappa in load 3: 96; kappa range 0.0050 to 0.1503 reading b: n = 125, kappa load 1 0.0382, load 3 0.0352, t(3-1) = -3.11, p = 0.0024; channels with lower kappa in load 3: 79; kappa range 0.0030 to 0.1750 paper: t(136) = -3.7453, p = 0.0001
Result: reproduced. Kappa is lower in load 3 under both readings (a: 144 channels, t = -4.54, p = 0.0001; b: 125 channels, t = -3.11, p = 0.0024).
edf3f-raw-mi-load¶
Raw (not z-scored) theta-gamma MI was also larger in load 1 than in load 3 in the hippocampal PAC channels.
rng = np.random.default_rng(306)
for rd in READINGS:
d = hpc_sig(rd)
t = sbcat.perm_ttest_paired(d.mi3_thetaHG, d.mi1_thetaHG, rng)
print(f"reading {rd}: n = {t['n']}, raw MI load 1 {t['mean_y']:.2e}, load 3 {t['mean_x']:.2e}, t(3-1) = {t['t']:.2f}, p = {t['p']:.4f}; "
f"channels with load 1 higher: {int((d.mi1_thetaHG > d.mi3_thetaHG).sum())}")
print("paper: n = 137; t(136) = -4.0264, p = 0.0001")
reading a: n = 144, raw MI load 1 2.43e-04, load 3 2.17e-04, t(3-1) = -4.03, p = 0.0002; channels with load 1 higher: 99 reading b: n = 125, raw MI load 1 2.81e-04, load 3 2.61e-04, t(3-1) = -1.85, p = 0.0657; channels with load 1 higher: 77 paper: n = 137; t(136) = -4.0264, p = 0.0001
Result: reproduced, resting on an assumption. Raw MI is larger in load 1 under reading (a) (t = -4.03, p = 0.0002, 99 of 144 channels) but not under (b) (p = 0.066). Reading (a) is adopted.
edf3e-nesting¶
110 of the 137 hippocampal PAC channels (80.29%) showed theta-high-gamma nesting, and PAC remained lower in load 3 than load 1 in those channels.
rng = np.random.default_rng(307)
print("Nesting readings: 'amp' aligns windows on the maximum of gamma amplitude (the Methods' words),")
print("'peak' on the maximum of the gamma-filtered signal (a gamma peak).")
for rd in READINGS:
d = hpc_sig(rd)
allh = ch[(ch.reading == rd) & (ch.area == "hippocampus") & (ch.n_per_load >= 3)]
for nr in ["amp", "peak"]:
nest = np.array([chan_index[(r.session, r.channel)]["readings"][rd]["nesting"][nr]["nested"] for r in d.itertuples()])
nmax = np.array([chan_index[(r.session, r.channel)]["readings"][rd]["nesting"][nr]["n_maxima"] for r in d.itertuples()])
dn = d[nest]
t = sbcat.perm_ttest_paired(dn.z3_thetaHG, dn.z1_thetaHG, rng) if len(dn) > 1 else {"t": np.nan, "p": np.nan, "mean_x": np.nan, "mean_y": np.nan}
nest_all = np.array([chan_index[(r.session, r.channel)]["readings"][rd]["nesting"][nr]["nested"] for r in allh.itertuples()])
print(f"artefact {rd}, nesting {nr}: nested {nest.sum()} of {len(d)} PAC channels ({100 * nest.mean():.2f}%), "
f"{100 * nest_all.mean():.1f}% of all {len(allh)} hippocampal channels; channels by number of local maxima {np.bincount(nmax).tolist()}; "
f"load test in nested channels t(3-1) = {t['t']:.2f}, p = {t['p']:.4f}, load 1 {t['mean_y']:.2f}, load 3 {t['mean_x']:.2f}")
print("paper: 110 of 137 (80.29%); t(109) = -4.10, p = 0.0001")
# How the count depends on the window for local maxima (45 ms in the Methods), reading (a)
d = hpc_sig("a")
rows_ = []
for nr in ["amp", "peak"]:
avgs = [chan_index[(r.session, r.channel)]["readings"]["a"]["nesting"][nr]["avg"] for r in d.itertuples()]
for w in [5, 7, 9, 11, 13, 17]:
nmax = np.array([sbcat.local_maxima(a_, w) for a_ in avgs])
rows_.append({"alignment": nr, "window_ms": round(2 * w * 2.5, 1), "nested (>= 3 maxima)": int((nmax >= 3).sum()), "of": len(d)})
print(pd.DataFrame(rows_).pivot(index="window_ms", columns="alignment", values="nested (>= 3 maxima)").to_string())
Nesting readings: 'amp' aligns windows on the maximum of gamma amplitude (the Methods' words), 'peak' on the maximum of the gamma-filtered signal (a gamma peak). artefact a, nesting amp: nested 9 of 144 PAC channels (6.25%), 5.3% of all 586 hippocampal channels; channels by number of local maxima [53, 58, 24, 7, 2]; load test in nested channels t(3-1) = 0.47, p = 0.7659, load 1 1.98, load 3 2.32 artefact a, nesting peak: nested 105 of 144 PAC channels (72.92%), 81.6% of all 586 hippocampal channels; channels by number of local maxima [3, 4, 32, 89, 15, 0, 1]; load test in nested channels t(3-1) = -3.61, p = 0.0002, load 1 4.29, load 3 3.48 artefact b, nesting amp: nested 12 of 125 PAC channels (9.60%), 8.7% of all 586 hippocampal channels; channels by number of local maxima [44, 48, 21, 11, 1]; load test in nested channels t(3-1) = 0.37, p = 0.7022, load 1 2.64, load 3 2.90 artefact b, nesting peak: nested 94 of 125 PAC channels (75.20%), 81.4% of all 586 hippocampal channels; channels by number of local maxima [2, 6, 23, 76, 15, 3]; load test in nested channels t(3-1) = -2.79, p = 0.0051, load 1 3.96, load 3 3.20 paper: 110 of 137 (80.29%); t(109) = -4.10, p = 0.0001
alignment amp peak window_ms 25.0 4 30 35.0 8 104 45.0 9 105 55.0 10 109 65.0 12 114 85.0 15 114
Result: reproduced, resting on an assumption. Aligned on gamma peaks, 105 of 144 PAC channels (72.92%) are nested under reading (a), within tolerance of 110, and PAC is lower in load 3 in them (t = -3.61, p = 0.0002). Aligned on gamma amplitude maxima, as the Methods' words read, only 9 are nested (6.25%) and the load test fails (p = 0.77). The count under the peak alignment is 104 to 114 for windows of 35 to 85 ms (30 for 25 ms), and 4 to 15 under the amplitude alignment for every window, so the outcome rests on the alignment and not on the window. The alignment on gamma peaks is adopted (see the attempt record).
methods-spectral-slope¶
The slope of the log-log power spectrum did not differ between load 1 and load 3 trials in hippocampal channels.
Welch spectrum (1 s Hann windows, 50% overlap; a standard choice) of the 0 to 2.5 s window of the subsampled trials, averaged per load; linear fit of log10 power on log10 frequency over (a) 2 to 150 Hz, (b) 30 to 150 Hz, leaving out 58-62 and 118-122 Hz (line-noise band-stops, a standard choice).
rng = np.random.default_rng(308)
freqs = R[0]["freqs"]
for rd in READINGS:
d = ch[(ch.reading == rd) & (ch.area == "hippocampus") & (ch.n_per_load >= 3)]
for fr_name, (f0, f1) in [("2-150 Hz", (2, 150)), ("30-150 Hz", (30, 150))]:
fm = (freqs >= f0) & (freqs <= f1) & ~((freqs >= 58) & (freqs <= 62)) & ~((freqs >= 118) & (freqs <= 122))
sl = {}
for L in (1, 3):
sl[L] = np.array([np.polyfit(np.log10(freqs[fm]), np.log10(chan_index[(r.session, r.channel)]["readings"][rd][f"psd_l{L}"][fm]), 1)[0]
for r in d.itertuples()])
t = sbcat.perm_ttest_paired(sl[1], sl[3], rng)
print(f"artefact {rd}, {fr_name}: n = {t['n']}, slope load 1 {sl[1].mean():.4f} +/- {sl[1].std(ddof=1):.4f}, "
f"load 3 {sl[3].mean():.4f} +/- {sl[3].std(ddof=1):.4f}, t(1-3) = {t['t']:.2f}, p = {t['p']:.4f}")
print("paper: n = 586; -1.7526 +/- 0.3902 against -1.7517 +/- 0.3928, t585 = -0.86, P = 0.39")
artefact a, 2-150 Hz: n = 586, slope load 1 -1.7260 +/- 0.4253, load 3 -1.7300 +/- 0.4314, t(1-3) = 2.66, p = 0.0073
artefact a, 30-150 Hz: n = 586, slope load 1 -1.8384 +/- 0.6506, load 3 -1.8416 +/- 0.6540, t(1-3) = 1.56, p = 0.1175
artefact b, 2-150 Hz: n = 586, slope load 1 -1.7236 +/- 0.4224, load 3 -1.7286 +/- 0.4292, t(1-3) = 3.14, p = 0.0018
artefact b, 30-150 Hz: n = 586, slope load 1 -1.8378 +/- 0.6518, load 3 -1.8400 +/- 0.6535, t(1-3) = 1.07, p = 0.2815 paper: n = 586; -1.7526 +/- 0.3902 against -1.7517 +/- 0.3928, t585 = -0.86, P = 0.39
Result: not reproduced, resting on an assumption. Over 586 hippocampal channels, the slope from 2 to 150 Hz is -1.7260 in load 1 and -1.7300 in load 3 (within tolerance of the paper's values), and the paired test finds a difference (artefact reading a: t = 2.66, p = 0.0073; b: p = 0.0018). From 30 to 150 Hz it does not (p = 0.12 and 0.28). The 2 to 150 Hz range, whose mean slope is closer to the paper's, is adopted (see the attempt record).
fig3h-pac-category-firing¶
Trial by trial, single-trial theta-high-gamma PAC was weakly positively related to the spike count of hippocampal category neurons, but not of amygdala category neurons.
Category neurons are selected as in reproduce_daume_2024.ipynb (same function and seed, so the
same 94 hippocampal and 181 amygdala neurons). Each is paired with every significant PAC channel
of its bundle (the electrodes of a session with the same location and coordinates; the files name
no bundle). Per pair, the correct trials clean for the channel; SC is the spike count 0 to 2.5 s
after maintenance onset; single-trial PAC under O5 (a) raw and (b) z-scored within the channel.
The model of Table S3, PAC ~ 1 + SC + load + SC x load + (1 | patient) + (1 | patient:pair),
fitted with statsmodels MixedLM (REML) in place of fitglme, with load 3 as the reference level
(Table S3, "Load 1 are tested against load 3 trials"), so the SC term is the slope in load 3
trials; the load 1 and average slopes are also printed. Under the trial readings of the
criterion: (a) all correct trials, (b) correct preferred trials. Artefact reading (a) of O1 is the
one adopted; reading (b) is also shown.
t0 = time.time()
cat_table, _ = sbcat.category_table(S, seed=103) # the same selection and seed as reproduce_daume_2024.ipynb
cat = cat_table[cat_table.is_category & cat_table.area.isin(sbcat.MTL)].reset_index(drop=True)
cat_trials = sbcat.trial_rates(S, cat)
print(f"category selection in {time.time() - t0:.0f} s:", cat.area.value_counts().to_dict())
h3_rows, pair_rows = [], []
for rd in READINGS:
parts = []
for i, u in cat.iterrows():
s = S[u.session]
keys = sbcat.bundle_keys(s["electrodes"])
on = s["trials"].timestamps_Maintenance.to_numpy()
sc_all = sbcat.count_in(s["spikes"][u.unit], on, on + 2.5)
for ci, e in enumerate(s["lfp_electrodes"]):
if keys[e] != keys[u.electrode] or (u.session, ci) not in chan_index:
continue
row = ch[(ch.reading == rd) & (ch.session == u.session) & (ch.channel == ci)].iloc[0]
if not row.z_thetaHG > 1.64:
continue
c = chan_index[(u.session, ci)]
idx = c["trial_index"]
keep = c["masks"][rd][idx]
st = c["st_thetaHG"][keep]
parts.append(pd.DataFrame({"pac_a": st, "pac_b": (st - np.nanmean(st)) / np.nanstd(st), "sc": sc_all[idx][keep],
"load1": (s["trials"].load.to_numpy()[idx][keep] == 1).astype(int),
"pref": cat_trials[i]["pref"][idx][keep], "patient": u.subject,
"pair": f"{u.session}_{u.unit}_{ci}", "area": u.area}))
df3 = pd.concat(parts).dropna().reset_index(drop=True)
for a in sbcat.MTL:
da = df3[df3.area == a]
pair_rows.append({"artefact": rd, "area": a, "pairs": da.pair.nunique(), "neurons": da.pair.str.rsplit("_", n=1).str[0].nunique(),
"patients": da.patient.nunique(), "trials": len(da), "spike count median": da.sc.median(), "max": da.sc.max()})
for trials in ["all", "preferred"]:
d = da[(trials == "all") | da.pref].reset_index(drop=True)
for o5 in ["a", "b"]:
try:
m = smf.mixedlm(f"pac_{o5} ~ sc * load1", d, groups="patient", vc_formula={"pair": "0 + C(pair)"}).fit(reml=True)
h3_rows.append({"artefact": rd, "area": a, "trials": trials, "O5": o5, "pairs": d.pair.nunique(), "n": len(d),
"intercept": m.params["Intercept"], **slopes(m, "sc")})
except Exception as e_: # report a failed fit
h3_rows.append({"artefact": rd, "area": a, "trials": trials, "O5": o5, "error": str(e_)[:80]})
print(pd.DataFrame(pair_rows).to_string(index=False))
print("paper: 151 hippocampal and 423 amygdala neuron-channel combinations")
print(pd.DataFrame(h3_rows).to_string(index=False, float_format=lambda v: f"{v:.4g}"))
print("SC slope in load 3 trials (L3, the paper's SC term, load 3 the reference as in Table S3), in load 1 trials (L1) and averaged")
print("paper (Table S3): hippocampus intercept -0.044, SC 0.002, t = 2.40, p = 0.017, Load1 0.076, SC:Load1 -0.004 (p = 0.010),")
print("so its load 1 slope is 0.002 - 0.004 = -0.002; amygdala SC -0.001, t = -0.75, p = 0.45, SC:Load1 0.002 (p = 0.14)")
category selection in 108 s: {'amygdala': 181, 'hippocampus': 94}
artefact area pairs neurons patients trials spike count median max
a hippocampus 166 52 15 20341 4.0 103
a amygdala 450 126 25 49813 5.0 140
b hippocampus 145 44 12 16319 5.0 103
b amygdala 410 115 23 40097 5.0 140
paper: 151 hippocampal and 423 amygdala neuron-channel combinations
artefact area trials O5 pairs n intercept L3 est L3 t L3 p L1 est L1 t L1 p mean est mean t mean p Load1 est int est int p converged
a hippocampus all a 166 20341 0.003858 7.826e-06 6.34 2.291e-10 1.586e-06 1.435 0.1513 4.706e-06 4.703 2.563e-06 7.759e-05 -6.239e-06 3.098e-07 True
a hippocampus all b 166 20341 -0.05654 0.004472 3.089 0.002008 -0.0003778 -0.2913 0.7708 0.002047 1.693 0.09049 0.0779 -0.00485 0.0002069 True
a hippocampus preferred a 166 8129 0.00386 6.152e-06 3.094 0.001973 -3.127e-06 -1.424 0.1546 1.512e-06 0.8378 0.4022 0.0001353 -9.28e-06 1.284e-05 True
a hippocampus preferred b 166 8129 -0.0505 0.002794 2.159 0.03083 -0.004517 -2.444 0.01452 -0.0008614 -0.758 0.4484 0.1287 -0.007311 0.001095 True
a amygdala all a 450 49813 0.004118 1.487e-07 0.1221 0.9028 4.994e-07 0.4289 0.668 3.241e-07 0.3088 0.7575 -6.426e-06 3.507e-07 0.7562 True
a amygdala all b 450 49813 -0.01132 0.001109 0.9815 0.3263 0.001346 1.245 0.2131 0.001228 1.253 0.21 0.002802 0.0002366 0.8178 True
a amygdala preferred a 450 19886 0.004117 4.374e-07 0.2547 0.7989 -2.572e-06 -1.193 0.2329 -1.068e-06 -0.6446 0.5192 4.331e-05 -3.01e-06 0.1431 True
a amygdala preferred b 450 19886 0.0001949 -9.379e-05 -0.09743 0.9224 -0.001295 -0.8273 0.4081 -0.0006943 -0.7494 0.4536 0.0274 -0.001201 0.5097 True
b hippocampus all a 145 16319 0.003883 8.115e-06 5.826 5.684e-09 2.057e-06 1.698 0.08951 5.086e-06 4.588 4.466e-06 6.837e-05 -6.058e-06 1.116e-05 True
b hippocampus all b 145 16319 -0.04535 0.0039 2.467 0.01364 -0.0005327 -0.3852 0.7001 0.001684 1.295 0.1952 0.0609 -0.004433 0.002043 True
b hippocampus preferred a 145 6402 0.003862 1.018e-05 5.683 1.32e-08 9.821e-07 0.4507 0.6522 5.584e-06 3.536 0.0004055 0.0001269 -9.203e-06 0.0001606 True
b hippocampus preferred b 145 6402 -0.03857 0.00283 1.908 0.05637 -0.004728 -2.261 0.02375 -0.0009489 -0.721 0.4709 0.1195 -0.007558 0.002433 True
b amygdala all a 410 40097 0.004155 -4.807e-07 -0.3476 0.7282 -1.069e-07 -0.08342 0.9335 -2.938e-07 -0.2508 0.802 -3.483e-06 3.738e-07 0.7691 True
b amygdala all b 410 40097 -0.01108 0.0006603 0.5329 0.5941 0.0008057 0.7 0.4839 0.000733 0.6932 0.4882 0.008942 0.0001454 0.8964 True
b amygdala preferred a 410 15573 0.004163 -1.563e-06 -0.797 0.4255 -3.43e-06 -1.432 0.1521 -2.496e-06 -1.34 0.1802 2.942e-05 -1.867e-06 0.4165 True
b amygdala preferred b 410 15573 0.00233 -0.001076 -1 0.3171 -0.00133 -0.8016 0.4228 -0.001203 -1.211 0.2258 0.02013 -0.0002539 0.8973 False
SC slope in load 3 trials (L3, the paper's SC term, load 3 the reference as in Table S3), in load 1 trials (L1) and averaged
paper (Table S3): hippocampus intercept -0.044, SC 0.002, t = 2.40, p = 0.017, Load1 0.076, SC:Load1 -0.004 (p = 0.010),
so its load 1 slope is 0.002 - 0.004 = -0.002; amygdala SC -0.001, t = -0.75, p = 0.45, SC:Load1 0.002 (p = 0.14)
Result: reproduced, resting on artefact reading (a). With load 3 as the reference level, as Table S3 codes load, the SC term is the slope in load 3 trials. Over 166 hippocampal pairs (paper 151) and 450 amygdala pairs (paper 423), under artefact reading (a), the hippocampal SC term is positive with p < 0.05 under every reading: all trials, O5 (a) t = 6.34, p = 2.3e-10, and O5 (b) 0.0045, t = 3.09, p = 0.0020 (paper 0.002, t = 2.40, p = 0.017); preferred trials, O5 (a) t = 3.09, p = 0.0020, and O5 (b) 0.0028, t = 2.16, p = 0.031. The amygdala term has p >= 0.326 under every reading (paper p = 0.45). Under O5 (b) on all trials the other terms are close to Table S3: intercept -0.057 (paper -0.044), Load1 0.078 (0.076), SC x Load1 -0.0049, p = 0.0002 (-0.004, p = 0.010). The slope in load 1 trials, which the earlier coding tested, is not positive (O5 b, all trials, p = 0.77; preferred, -0.0045, p = 0.015), as the paper's own load 1 slope (0.002 - 0.004) is not. Under artefact reading (b) (145 and 410 pairs) the hippocampal term is positive with p < 0.05 under three readings (p <= 0.014) and has p = 0.056 under the fourth (preferred, O5 b); the amygdala term has p >= 0.32 (one of its fits did not converge). This outcome replaced not reproduced after the third verification round.
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 |
|---|---|
| fig2a-theta-gamma-pac | reproduced |
| fig2b-pac-channel-proportions | partially reproduced (17 pre-SMA channels against 5) |
| fig2c-hippocampus-load | reproduced |
| fig2c-no-load-effect-other-areas | not reproduced (pre-SMA, and amygdala or dACC) |
| edf2b-hippocampus-load-patients | reproduced, resting on artefact reading (a) |
| edf2c-low-gamma-no-load-effect | not reproduced (pre-SMA) |
| tables1-slow-fast-theta | not reproduced (pre-SMA tests also significant) |
| fig2f-rt-single-trial-pac | blocked (methods underspecified: the readings of single-trial PAC disagree), resting on artefact reading (a) |
| edf3a-kappa-load | reproduced |
| edf3e-nesting | reproduced, resting on the alignment on gamma peaks |
| edf3f-raw-mi-load | reproduced, resting on artefact reading (a) |
| methods-spectral-slope | not reproduced, resting on the 2 to 150 Hz range |
| fig3h-pac-category-firing | reproduced, resting on artefact reading (a) |
| edf3-hippocampus-pac-controls | in reproduce_daume_2024_lfp_controls.ipynb |
How the analysis differs from the paper's: the artefact rejection is the automated part only (the visual inspection cannot be repeated), under two readings of the span it is computed on; reading (a) removes 11.4% of trials per channel and reading (b) 25.9%, against the paper's 14.6% with visual inspection; where the readings disagree, (a) is adopted, as the paper's own rejected trials on its example channel contain every trial (a) removes and only 4 of the 9 that (b) removes. The single-trial PAC, the nesting alignment, the slope's frequency range and the power estimate are not given by the paper. Mixed models are fitted with statsmodels in place of fitglme, with load 3 as the reference level as in Tables S2 and S3; RT is fitted in s, as log10 RT and z-scored within session, as Table S2 gives no unit or transform.
What was found about the dataset: the LFP of 1,922 channels is there for 43 of the 44 sessions (sub-20 has none), at 400 Hz in microvolts, and shares the trial clock (probe-locked contrast, 37 of 43 sessions). The files do not mark the trials or channels the paper's visual inspection removed. Without them, the 17 significant pre-SMA channels, with a mean z of 6.87 in load 3 against 1.81 in load 1 under reading (a), decide several of the absence claims.
print(f"notebook ran in {(time.time() - started) / 60:.1f} min")
notebook ran in 44.2 min