Daume et al. (2024) on dandiset 000673: PAC neurons¶
Paper. Daume J et al. (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.
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 selects PAC neurons (Methods, Selection of PAC neurons; sbcat.py says where each
parameter comes from) and tests the claims on them that need only spikes and the LFP of their
bundle. Artefact rejection uses reading (a) of O1 (on the epochs), the reading adopted in
reproduce_daume_2024_lfp.ipynb.
| Claim | Tested here |
|---|---|
| pac-neuron-counts | yes |
| pac-category-overlap | yes |
| fig4c-pac-neuron-firing-single-trial-pac | yes |
| fig4e-hippocampus-pac-neurons-accuracy | yes |
| fig4g-hippocampus-pac-neurons-load | yes |
| fig4-amygdala-pac-neurons | yes |
| SFC, noise correlation and decoding claims on PAC neurons | blocked (attempt record) |
In [1]:
import os
for _k in ["OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"]:
os.environ[_k] = "1"
import time
import warnings
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()
cat_table, _ = sbcat.category_table(S, seed=103) # the same selection and seed as reproduce_daume_2024.ipynb
print(f"category selection in {time.time() - t0:.0f} s")
t0 = time.time()
P = sbcat.run_all_pac_neurons(S, seed=401, workers=6)
pairs = pd.DataFrame([r for p in P for r in p["rows"]])
counts = np.concatenate([p["counts"] for p in P if len(p["counts"])])
print(f"spike counts for {len(pairs)} neuron-channel-load combinations x {counts.shape[1]} pairings in {(time.time() - t0) / 60:.1f} min")
loaded trials and units of 44 files in 131 s
category selection in 115 s
spike counts for 9640 neuron-channel-load combinations x 201 pairings in 3.4 min
Data check¶
In [2]:
mtl_units = pd.DataFrame([{"session": p["session"], "unit": u} for p in P for u in p["mtl_units"]])
mtl_units["area"] = [S[r.session]["units"].area[r.unit] for r in mtl_units.itertuples()]
with_lfp = pairs.groupby(["session", "unit"]).size().index
print("MTL neurons in sessions with LFP:", mtl_units.area.value_counts().to_dict())
avail = pairs.drop_duplicates(["session", "unit"]).area.value_counts()
print("MTL neurons with at least one LFP channel in their bundle:", avail.to_dict(), "(paper: 212 hippocampal, 358 amygdala)")
print("channels per neuron (bundle):", pairs[pairs.load == 1].groupby(["session", "unit"]).size().value_counts().sort_index().to_dict())
print("trials per combination: median %d, range %d to %d; spikes per combination: median %d, combinations with no spike: %d" % (
pairs.n_trials.median(), pairs.n_trials.min(), pairs.n_trials.max(), pairs.n_spikes.median(), int((pairs.n_spikes == 0).sum())))
print("total MTL neurons (all sessions):", int(cat_table.area.isin(sbcat.MTL).sum()))
MTL neurons in sessions with LFP: {'amygdala': 490, 'hippocampus': 338}
MTL neurons with at least one LFP channel in their bundle: {'amygdala': 422, 'hippocampus': 269} (paper: 212 hippocampal, 358 amygdala)
channels per neuron (bundle): {3: 14, 6: 63, 7: 512, 8: 102}
trials per combination: median 58, range 15 to 70; spikes per combination: median 183, combinations with no spike: 57
total MTL neurons (all sessions): 856
pac-neuron-counts¶
79 of 212 hippocampal neurons (37.29%) and 163 of 358 amygdala neurons (45.53%) qualified as PAC neurons, more than expected by chance.
In [3]:
t0 = time.time()
K, R1, _ = counts.shape
p12 = np.empty((K, R1)); p13 = np.empty((K, R1)); r2 = np.empty((K, R1))
flat = counts.reshape(-1, 20)
chunk = 100000
p12f, p13f, r2f = [], [], []
for i in range(0, len(flat), chunk):
o = sbcat.pac_neuron_tests(flat[i:i + chunk])
p12f.append(o["p12"]); p13f.append(o["p13"]); r2f.append(o["r2"])
p12 = np.concatenate(p12f).reshape(K, R1); p13 = np.concatenate(p13f).reshape(K, R1); r2 = np.concatenate(r2f).reshape(K, R1)
print(f"{3 * flat.shape[0]} Poisson GLMs in {time.time() - t0:.0f} s")
empty = counts.sum(2) == 0
p12[empty] = 1; p13[empty] = 1
def fdr_matrix(p, groups):
"""Benjamini-Hochberg within each group of rows, separately for each column (pairing)."""
q = np.empty_like(p)
for g in np.unique(groups):
m = groups == g
sub = p[m]
n = sub.shape[0]
order = np.argsort(sub, axis=0)
ranked = np.take_along_axis(sub, order, 0) * n / np.arange(1, n + 1)[:, None]
ranked = np.minimum.accumulate(ranked[::-1], axis=0)[::-1]
out = np.empty_like(sub)
np.put_along_axis(out, order, np.minimum(ranked, 1), 0)
q[m] = out
return q
readings = {
"a": (pairs.area + "|" + pairs.load.astype(str)).to_numpy(), # all combinations of an area and load
"b": (pairs.session.astype(str) + "|" + pairs.unit.astype(str) + "|" + pairs.load.astype(str)).to_numpy(), # a neuron's channels
}
neuron_key = pairs.session.astype(str) + "|" + pairs.unit.astype(str)
sel = {}
rows = []
for rd, groups in readings.items():
q12, q13 = fdr_matrix(p12, groups), fdr_matrix(p13, groups)
ok = (q12 < sbcat.PAC_NEURON_ALPHA) & (q13 < sbcat.PAC_NEURON_ALPHA) # K x pairings
per_neuron = pd.DataFrame(ok).groupby(neuron_key.to_numpy()).any() # neurons x pairings
area_of = pairs.drop_duplicates(["session", "unit"]).set_index(neuron_key.drop_duplicates())["area"]
cnt = per_neuron.groupby(area_of.reindex(per_neuron.index).to_numpy()).sum() # areas x pairings
sel[rd] = (ok[:, 0], per_neuron[0])
for a in sbcat.MTL:
obs = int(cnt.loc[a, 0]); null = cnt.loc[a].to_numpy()[1:]
n_av = int(avail[a])
rows.append({"FDR": rd, "area": a, "available": n_av, "pac_neurons": obs, "percent": 100 * obs / n_av,
"paper": {"hippocampus": 79, "amygdala": 163}[a], "null_mean": null.mean(), "null_max": null.max(),
"p_shuffle": (np.sum(null >= obs) + 1) / (len(null) + 1)})
pac_counts = pd.DataFrame(rows)
print(pac_counts.round(4).to_string(index=False))
print("paper: hippocampus 79 of 212 (37.29%), amygdala 163 of 358 (45.53%), P < 0.005 each")
# how the count depends on the threshold (reading a)
for thr in [0.001, 0.01, 0.05]:
q12, q13 = fdr_matrix(p12[:, :1], readings["a"]), fdr_matrix(p13[:, :1], readings["a"])
okt = pd.Series(((q12 < thr) & (q13 < thr))[:, 0]).groupby(neuron_key.to_numpy()).any()
print(f"FDR threshold {thr}:", okt.groupby(area_of.reindex(okt.index).to_numpy()).sum().to_dict())
print("uncorrected share of combinations with p12 < 0.01 and p13 < 0.01: observed %.3f, shuffled %.3f" % (
((p12[:, 0] < 0.01) & (p13[:, 0] < 0.01)).mean(), ((p12[:, 1:] < 0.01) & (p13[:, 1:] < 0.01)).mean()))
# Uncorrected counts per area: neurons with at least one combination (channel, load) passing both tests at p < 0.01
okr = pd.DataFrame((p12 < 0.01) & (p13 < 0.01)).groupby(neuron_key.to_numpy()).any()
cu = okr.groupby(area_of.reindex(okr.index).to_numpy()).sum()
for a in sbcat.MTL:
null = cu.loc[a].to_numpy()[1:]
print(f"{a}, uncorrected: {int(cu.loc[a, 0])} of {int(avail[a])} neurons; shuffles mean {null.mean():.1f}, max {null.max()}")
print("combinations per neuron (channels x loads): median %d, range %d to %d" % (
pairs.groupby(["session", "unit"]).size().median(), pairs.groupby(["session", "unit"]).size().min(), pairs.groupby(["session", "unit"]).size().max()))
5812920 Poisson GLMs in 244 s
FDR area available pac_neurons percent paper null_mean null_max p_shuffle
a hippocampus 269 8 2.9740 79 0.035 1 0.005
a amygdala 422 5 1.1848 163 0.005 1 0.005
b hippocampus 269 12 4.4610 79 1.295 5 0.005
b amygdala 422 24 5.6872 163 1.100 4 0.005
paper: hippocampus 79 of 212 (37.29%), amygdala 163 of 358 (45.53%), P < 0.005 each
FDR threshold 0.001: {'amygdala': 4, 'hippocampus': 4}
FDR threshold 0.01: {'amygdala': 5, 'hippocampus': 8}
FDR threshold 0.05: {'amygdala': 6, 'hippocampus': 12}
uncorrected share of combinations with p12 < 0.01 and p13 < 0.01: observed 0.018, shuffled 0.001
hippocampus, uncorrected: 41 of 269 neurons; shuffles mean 6.2, max 13
amygdala, uncorrected: 66 of 422 neurons; shuffles mean 6.9, max 14
combinations per neuron (channels x loads): median 14, range 6 to 16
Result: partially reproduced. With FDR over the combinations of an area and load (O9 a), 8 of
269 hippocampal (2.97%) and 5 of 422 amygdala neurons (1.18%) qualify, against 79 of 212 and 163
of 358; with FDR over each neuron's channels (b), 12 and 24. Both counts exceed every shuffled
count (p = 0.005 in both areas under both readings), so the selection finds more than
chance, but the counts are far out of tolerance. Even without any correction, 41 hippocampal and
66 amygdala neurons have a combination passing both tests at p < 0.01 (shuffles at most 13 and
14), fewer than the paper's counts after FDR. The paper's plot_PAC_neuron.m fits the three GLMs to
the 20 aggregated bin counts of one neuron-channel pair, with no random effect (fitglme with
fixed effects only), as sbcat.py does, so neither a per-trial layout nor a random effect is the
difference.
In [4]:
# PAC neurons under FDR reading (a), with the channel of highest model 1 R2 among their qualifying combinations
ok_a = sel["a"][0]
q = pairs.assign(ok=ok_a, r2=r2[:, 0])
pac = (q[q.ok].sort_values("r2", ascending=False).drop_duplicates(["session", "unit"])
[["session", "subject", "unit", "area", "channel", "lfp_electrode", "load", "r2"]].reset_index(drop=True))
print("PAC neurons by area:", pac.area.value_counts().to_dict())
cat_key = set(zip(cat_table.session[cat_table.is_category], cat_table.unit[cat_table.is_category]))
pac["is_category"] = [(r.session, r.unit) in cat_key for r in pac.itertuples()]
PAC neurons by area: {'hippocampus': 8, 'amygdala': 5}
pac-category-overlap¶
The overlap of PAC and category neurons (28 of 79 in the hippocampus, 68 of 163 in the amygdala) did not exceed that expected of independent populations.
In [5]:
analysed = pairs.drop_duplicates(["session", "unit"])[["session", "unit", "area"]]
analysed["is_category"] = [(r.session, r.unit) in cat_key for r in analysed.itertuples()]
rows = []
for a in sbcat.MTL:
an = analysed[analysed.area == a]
N, Kc = len(an), int(an.is_category.sum())
pa = pac[pac.area == a]
n, k = len(pa), int(pa.is_category.sum())
p_bin = stats.binom.sf(k - 1, N, (n / N) * (Kc / N))
p_hyp = stats.hypergeom.sf(k - 1, N, Kc, n)
rows.append({"area": a, "analysed": N, "category": Kc, "pac": n, "overlap": k, "percent_of_pac": 100 * k / n if n else np.nan,
"expected": n * Kc / N, "p_binomial": p_bin, "p_hypergeometric": p_hyp})
print(pd.DataFrame(rows).round(4).to_string(index=False))
print("paper: 28 of 79 (35.44%) and 68 of 163 (41.72%), both P > 0.05")
area analysed category pac overlap percent_of_pac expected p_binomial p_hypergeometric hippocampus 269 70 8 4 50.0 2.0818 0.1573 0.1251 amygdala 422 162 5 2 40.0 1.9194 0.5723 0.6358 paper: 28 of 79 (35.44%) and 68 of 163 (41.72%), both P > 0.05
Result: blocked. With 8 hippocampal and 5 amygdala PAC neurons (paper 79 and 163), the overlap (4 of 8 and 2 of 5) does not test the paper's populations; the tests are printed above.
Firing of PAC neurons¶
Rates as in reproduce_daume_2024.ipynb: maintenance 0 to 2.5 s, baseline -0.9 to -0.3 s before
picture 1 (mean over correct trials); readings of O6 (a) Hz, (b) % of baseline; neurons with
fewer incorrect trials than 1, 2 or 3 are left out of the correct against incorrect tests (O7).
In [6]:
pac_rates = sbcat.trial_rates(S, pac.assign(preferred=0))
def nmean(r, mask, reading):
base = r["base"][r["correct"]].mean()
v = r["maint"][mask].mean() if mask.sum() else np.nan
if reading == "b":
return 100 * (v - base) / base if base > 0 else np.nan
return v
def bmean(r, reading):
base = r["base"][r["correct"]].mean()
return (0.0 if base > 0 else np.nan) if reading == "b" else base
def rate_tests(idx, rng):
rows = []
for reading in ["a", "b"]:
rs = [pac_rates[i] for i in idx]
m_all = np.array([nmean(r, r["correct"], reading) for r in rs])
b = np.array([bmean(r, reading) for r in rs])
t = sbcat.perm_ttest_paired(m_all, b, rng)
rows.append({"reading": reading, "test": "maintenance vs baseline (correct)", "min_inc": "", "n": t["n"], "x": t["mean_x"], "y": t["mean_y"], "t": t["t"], "p": t["p"]})
for L in (1, 3):
mL = np.array([nmean(r, r["correct"] & (r["load"] == L), reading) for r in rs])
t = sbcat.perm_ttest_paired(mL, b, rng)
rows.append({"reading": reading, "test": f"load {L} vs baseline", "min_inc": "", "n": t["n"], "x": t["mean_x"], "y": t["mean_y"], "t": t["t"], "p": t["p"]})
m1 = np.array([nmean(r, r["correct"] & (r["load"] == 1), reading) for r in rs])
m3 = np.array([nmean(r, r["correct"] & (r["load"] == 3), reading) for r in rs])
t = sbcat.perm_ttest_paired(m3, m1, rng)
rows.append({"reading": reading, "test": "load 3 vs load 1", "min_inc": "", "n": t["n"], "x": t["mean_x"], "y": t["mean_y"], "t": t["t"], "p": t["p"]})
fv, sv = [], []
for r in rs:
fast, slow = sbcat.rt_split(r["load"], r["rt"], r["correct"])
fv.append(np.nanmean([nmean(r, fast & (r["load"] == L), reading) for L in (1, 3)]))
sv.append(np.nanmean([nmean(r, slow & (r["load"] == L), reading) for L in (1, 3)]))
t = sbcat.perm_ttest_paired(fv, sv, rng)
rows.append({"reading": reading, "test": "fast vs slow RT", "min_inc": "", "n": t["n"], "x": t["mean_x"], "y": t["mean_y"], "t": t["t"], "p": t["p"]})
for mi in [1, 2, 3]:
rk = [r for r in rs if (~r["correct"]).sum() >= mi]
c = np.array([nmean(r, r["correct"], reading) for r in rk])
inc = np.array([nmean(r, ~r["correct"], reading) for r in rk])
bb = np.array([bmean(r, reading) for r in rk])
for name, x, y in [("correct vs incorrect", c, inc), ("correct vs baseline", c, bb), ("incorrect vs baseline", inc, bb)]:
t = sbcat.perm_ttest_paired(x, y, rng)
rows.append({"reading": reading, "test": name, "min_inc": mi, "n": t["n"], "x": t["mean_x"], "y": t["mean_y"], "t": t["t"], "p": t["p"]})
return pd.DataFrame(rows)
fig4e-hippocampus-pac-neurons-accuracy and fig4g-hippocampus-pac-neurons-load¶
Hippocampal PAC neurons fired above baseline during maintenance, more in correct than in incorrect trials, and above baseline in correct trials only.
Hippocampal PAC neuron rates did not differ between loads but were above baseline in each load; they also did not differ between fast and slow RT trials.
In [7]:
rng = np.random.default_rng(402)
hpc_tests = rate_tests(np.where(pac.area == "hippocampus")[0], rng)
print(f"hippocampal PAC neurons: {int((pac.area == 'hippocampus').sum())} (paper 79)")
print(hpc_tests.round(4).to_string(index=False))
print("paper: maintenance vs baseline t78 = 2.43, p = 0.01; correct vs incorrect t62 = 3.82, p < 0.001 (n = 63); "
"correct vs baseline t62 = 2.67, p = 0.01; incorrect vs baseline t62 = -0.98, p = 0.33")
hippocampal PAC neurons: 8 (paper 79)
reading test min_inc n x y t p
a maintenance vs baseline (correct) 8 6.9347 6.6807 0.5708 0.5084
a load 1 vs baseline 8 7.6635 6.6807 1.5957 0.1435
a load 3 vs baseline 8 6.1466 6.6807 -1.6213 0.1627
a load 3 vs load 1 8 6.1466 7.6635 -3.3453 0.0086
a fast vs slow RT 8 6.5838 7.2353 -2.3625 0.0579
a correct vs incorrect 1 8 6.9347 5.6923 2.9819 0.0236
a correct vs baseline 1 8 6.9347 6.6807 0.5708 0.5081
a incorrect vs baseline 1 8 5.6923 6.6807 -1.8352 0.1050
a correct vs incorrect 2 7 7.9118 6.5055 3.1795 0.0337
a correct vs baseline 2 7 7.9118 7.6316 0.5461 0.5239
a incorrect vs baseline 2 7 6.5055 7.6316 -1.8732 0.1084
a correct vs incorrect 3 7 7.9118 6.5055 3.1795 0.0317
a correct vs baseline 3 7 7.9118 7.6316 0.5461 0.5152
a incorrect vs baseline 3 7 6.5055 7.6316 -1.8732 0.1079
b maintenance vs baseline (correct) 8 45.7638 0.0000 1.2310 0.2602
b load 1 vs baseline 8 85.0530 0.0000 1.3331 0.0823
b load 3 vs baseline 8 5.1402 0.0000 0.4830 0.7326
b load 3 vs load 1 8 5.1402 85.0530 -1.4855 0.0089
b fast vs slow RT 8 49.4449 41.1279 0.3649 0.9856
b correct vs incorrect 1 8 45.7638 -19.5272 1.3772 0.0170
b correct vs baseline 1 8 45.7638 0.0000 1.2310 0.2627
b incorrect vs baseline 1 8 -19.5272 0.0000 -1.2792 0.2646
b correct vs incorrect 2 7 10.0157 -8.0311 3.9421 0.0332
b correct vs baseline 2 7 10.0157 0.0000 0.8503 0.5176
b incorrect vs baseline 2 7 -8.0311 0.0000 -0.6926 0.5434
b correct vs incorrect 3 7 10.0157 -8.0311 3.9421 0.0310
b correct vs baseline 3 7 10.0157 0.0000 0.8503 0.5143
b incorrect vs baseline 3 7 -8.0311 0.0000 -0.6926 0.5447
paper: maintenance vs baseline t78 = 2.43, p = 0.01; correct vs incorrect t62 = 3.82, p < 0.001 (n = 63); correct vs baseline t62 = 2.67, p = 0.01; incorrect vs baseline t62 = -0.98, p = 0.33
In [8]:
print(hpc_tests[hpc_tests.test.isin(["load 1 vs baseline", "load 3 vs baseline", "load 3 vs load 1", "fast vs slow RT"])].round(4).to_string(index=False))
print("paper: load 3 - load 1 t78 = 1.38, p = 0.20; load 1 vs baseline p = 0.03; load 3 vs baseline p = 0.01; fast vs slow t78 = -0.07, p = 0.94")
reading test min_inc n x y t p
a load 1 vs baseline 8 7.6635 6.6807 1.5957 0.1435
a load 3 vs baseline 8 6.1466 6.6807 -1.6213 0.1627
a load 3 vs load 1 8 6.1466 7.6635 -3.3453 0.0086
a fast vs slow RT 8 6.5838 7.2353 -2.3625 0.0579
b load 1 vs baseline 8 85.0530 0.0000 1.3331 0.0823
b load 3 vs baseline 8 5.1402 0.0000 0.4830 0.7326
b load 3 vs load 1 8 5.1402 85.0530 -1.4855 0.0089
b fast vs slow RT 8 49.4449 41.1279 0.3649 0.9856
paper: load 3 - load 1 t78 = 1.38, p = 0.20; load 1 vs baseline p = 0.03; load 3 vs baseline p = 0.01; fast vs slow t78 = -0.07, p = 0.94
Result: blocked. The tests above are on 8 hippocampal PAC neurons, against the paper's 79, and are not a test of the paper's population.
fig4-amygdala-pac-neurons¶
Amygdala PAC neurons fired above baseline during maintenance but did not differ between correct and incorrect trials, loads, or fast and slow RT trials.
In [9]:
rng = np.random.default_rng(403)
amy_tests = rate_tests(np.where(pac.area == "amygdala")[0], rng)
print(f"amygdala PAC neurons: {int((pac.area == 'amygdala').sum())} (paper 163)")
print(amy_tests.round(4).to_string(index=False))
print("paper: maintenance vs baseline t162 = 6.40, p < 0.001; correct vs incorrect t156 = -0.77, p = 0.45 (n = 156); "
"load 3 - load 1 t162 = 0.26, p = 0.80; fast vs slow t162 = 1.87, p = 0.07")
amygdala PAC neurons: 5 (paper 163)
reading test min_inc n x y t p
a maintenance vs baseline (correct) 5 9.6658 8.3455 2.7921 0.1232
a load 1 vs baseline 5 9.5508 8.3455 3.7229 0.0619
a load 3 vs baseline 5 9.7875 8.3455 2.1874 0.1822
a load 3 vs load 1 5 9.7875 9.5508 0.5761 0.5628
a fast vs slow RT 5 9.8373 9.4976 1.5550 0.3759
a correct vs incorrect 1 5 9.6658 10.0848 -0.4855 0.7501
a correct vs baseline 1 5 9.6658 8.3455 2.7921 0.1337
a incorrect vs baseline 1 5 10.0848 8.3455 1.4582 0.2431
a correct vs incorrect 2 4 11.1053 11.5061 -0.3598 0.8731
a correct vs baseline 2 4 11.1053 9.6376 2.5305 0.2547
a incorrect vs baseline 2 4 11.5061 9.6376 1.2206 0.3824
a correct vs incorrect 3 4 11.1053 11.5061 -0.3598 0.8796
a correct vs baseline 3 4 11.1053 9.6376 2.5305 0.2575
a incorrect vs baseline 3 4 11.5061 9.6376 1.2206 0.3703
b maintenance vs baseline (correct) 5 20.4708 0.0000 2.6660 0.1283
b load 1 vs baseline 5 18.7601 0.0000 3.1624 0.0639
b load 3 vs baseline 5 22.3148 0.0000 2.2517 0.1278
b load 3 vs load 1 5 22.3148 18.7601 0.6891 0.5619
b fast vs slow RT 5 22.8754 18.1517 1.5601 0.2502
b correct vs incorrect 1 5 20.4708 25.2042 -0.5568 0.5513
b correct vs baseline 1 5 20.4708 0.0000 2.6660 0.1305
b incorrect vs baseline 1 5 25.2042 0.0000 2.1888 0.1270
b correct vs incorrect 2 4 19.8413 21.8864 -0.1964 0.8703
b correct vs baseline 2 4 19.8413 0.0000 2.0083 0.2461
b incorrect vs baseline 2 4 21.8864 0.0000 1.5374 0.2503
b correct vs incorrect 3 4 19.8413 21.8864 -0.1964 0.8775
b correct vs baseline 3 4 19.8413 0.0000 2.0083 0.2567
b incorrect vs baseline 3 4 21.8864 0.0000 1.5374 0.2577
paper: maintenance vs baseline t162 = 6.40, p < 0.001; correct vs incorrect t156 = -0.77, p = 0.45 (n = 156); load 3 - load 1 t162 = 0.26, p = 0.80; fast vs slow t162 = 1.87, p = 0.07
Result: blocked. The tests above are on 5 amygdala PAC neurons, against the paper's 163.
fig4c-pac-neuron-firing-single-trial-pac¶
Spike counts of hippocampal PAC neurons were positively related to single-trial theta-high-gamma PAC on their selected channel; amygdala PAC neurons showed no relation.
Correct trials, clean for the selected channel; PAC ~ 1 + SC * load + (1 | patient) + (1 | neuron-channel pair nested in patient) (Supplementary Table S4 and Methods), SC the spike count in 0 to 2.5 s; single-trial PAC (O5) (a) the trial's MI, raw, (b) z-scored across trials within the pair.
In [10]:
st_index = {(p["session"], ci): (p["st_mi"][ci], p["clean"][ci]) for p in P for ci in p["st_mi"]}
glm_rows = []
for a in sbcat.MTL:
parts = []
for r in pac[pac.area == a].itertuples():
s = S[r.session]
tr = s["trials"]
st, clean = st_index[(r.session, r.channel)]
keep = tr.correct.to_numpy() & clean
t0_ = tr.timestamps_Maintenance.to_numpy()[keep]
sc = sbcat.count_in(s["spikes"][r.unit], t0_, t0_ + 2.5, closed=True)
v = st[keep]
parts.append(pd.DataFrame({"pac_a": v, "pac_b": (v - np.nanmean(v)) / np.nanstd(v), "sc": sc,
"load": tr.load.to_numpy()[keep].astype(str), "patient": r.subject, "pair": f"{r.session}_{r.unit}_{r.channel}"}))
if not parts:
glm_rows.append({"area": a, "pairs": 0})
continue
df = pd.concat(parts).dropna().reset_index(drop=True)
for o5 in ["a", "b"]:
m = smf.mixedlm(f"pac_{o5} ~ sc * load", df, groups="patient", vc_formula={"pair": "0 + C(pair)"}).fit(reml=True)
glm_rows.append({"area": a, "O5": o5, "pairs": df.pair.nunique(), "trials": len(df), "estimate": m.params["sc"],
"t": m.tvalues["sc"], "p": m.pvalues["sc"], "converged": m.converged})
print(pd.DataFrame(glm_rows).to_string(index=False))
print("paper (Table S4): hippocampus estimate 0.0032, t = 2.19, p = 0.028 (n = 79); amygdala estimate -3e-5, p = 0.98 (n = 163)")
area O5 pairs trials estimate t p converged hippocampus a 8 983 8.755584e-07 0.183206 0.854636 True hippocampus b 8 983 -4.554168e-03 -1.352057 0.176357 True amygdala a 5 420 -3.455215e-05 -1.962637 0.049688 True amygdala b 5 420 -1.177520e-02 -2.038363 0.041514 True paper (Table S4): hippocampus estimate 0.0032, t = 2.19, p = 0.028 (n = 79); amygdala estimate -3e-5, p = 0.98 (n = 163)
Result: blocked. The fits above use 8 hippocampal and 5 amygdala neuron-channel pairs, against the paper's 79 and 163.
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 |
|---|---|
| pac-neuron-counts | partially reproduced (8 and 5 PAC neurons against 79 and 163) |
| pac-category-overlap | blocked (too few PAC neurons) |
| fig4c-pac-neuron-firing-single-trial-pac | blocked (too few PAC neurons) |
| fig4e-hippocampus-pac-neurons-accuracy | blocked (too few PAC neurons) |
| fig4g-hippocampus-pac-neurons-load | blocked (too few PAC neurons) |
| fig4-amygdala-pac-neurons | blocked (too few PAC neurons) |
How the analysis differs from the paper's: the selection follows the paper's plot_PAC_neuron.m (20 bins, three Poisson GLMs, likelihood-ratio tests); the artefact rejection is the automated part only (reading a), and the FDR family is not stated (two readings shown). What was found about the dataset: 269 hippocampal and 422 amygdala neurons have LFP in their bundle, more than the paper's 212 and 358, so the files do not say which LFPs the paper removed.
In [11]:
print(f"notebook ran in {(time.time() - started) / 60:.1f} min")
notebook ran in 11.6 min