Kauvar, Richman, Liu et al. (2025): human iEEG spectral response (Figs. 2G to 2J, 3B, 3C; figs. S12, S13, S15, S27)¶
Paper: Kauvar IV, Richman EB, Liu TX et al. (2025), Conserved brain-wide emergence of emotional response from sensory experience in humans and mice. Science 388, eadt3971, doi:10.1126/science.adt3971. The claims file (claims/001347/10.1126-science.adt3971.yaml) was written from the version of record with its supplementary materials.
Dataset: DANDI:001347, version 0.250528.0702, streamed. The paper cites the dandiset without a version; this is the only published version (28 May 2025), so it is pinned.
What this notebook uses: the aggregated human products in sub-H/ (per-channel trial means split into
train and test halves, per-channel, per-network and per-region per-trial spectrograms), and, to compare,
the paper's own deposited outputs in the same folder (NMF factors and loadings, cluster-test maps). The
Yeo7 network of each channel is not in the dataset; it is read from the paper's code record (Zenodo
14920570, anat/pts_ch_labels.pkl) and checked against the dataset in the link claim
link-ieeg-yeo-labels. Nothing random is drawn: NMF uses random_state 1 and the cluster tests seed 1,
as in the paper's code; the label permutation of figs27 is enumerated exactly. figs15 and figs27 use the
per-participant 20-min LFP files and, for figs27, the CADSS scores of the code record (data/cadss_across_pts.csv,
link claim cadss-scores-link in notebook 01).
| Claim | Data | Status |
|---|---|---|
| fig2g-nmf-six-factors | dataset | tested |
| link-ieeg-yeo-labels | with labels | tested (link claim) |
| fig2h-factor-network-loading | with labels | tested |
| fig2i-ketamine-reduces-factors-1-2 | dataset | tested |
| fig2j-factor-change-by-network | with labels | tested |
| figs13-auditory-control-no-response | dataset | tested |
| fig3b-cluster-test-yeo-networks | dataset | tested |
| fig3c-cluster-test-regions | dataset | tested |
| figs15-ketamine-spectral-changes-frontal | dataset | tested (approximate windows) |
| figs12-erp-not-explained-by-eyelid-closure | dataset | blocked (data_missing) |
| figs27-slow-oscillations-with-strong-dissociation | dataset, CADSS from the code record | tested (approximate windows) |
import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import time
import warnings
from collections import Counter
import numpy as np
import pandas as pd
import yaml
import matplotlib.pyplot as plt
import dandi_io
import ieeg_spectral as ie
warnings.filterwarnings("ignore") # sklearn ConvergenceWarning and statsmodels fit notes; iterations are reported
pd.set_option("display.width", 200)
np.set_printoptions(precision=3, suppress=True, linewidth=160)
T0 = time.time()
claims = {c["id"]: c for c in yaml.safe_load(open("../../claims/001347/10.1126-science.adt3971.yaml"))["claims"]}
print("dandiset", dandi_io.DANDISET, "version", dandi_io.VERSION)
dandiset 001347 version 0.250528.0702
Data check¶
Checks made below, each against what the claims file and the paper say:
- The per-channel file holds 458 channels from 7 participants with the counts of Table S3 (SD168 86, SD172 97, SD181 58, SD185 43, SD186 47, SD211 44, SD245 83), all in the 21 regions of Fig. 3C.
- Its train and test arrays are exactly the means of the odd-numbered and even-numbered trials (code split 1 and split 0) of the per-trial file, over the bins 0 < t < 1.4 s, for one channel per participant and each session. This ties the two files together and fixes what "train" and "test" are.
- Trials per participant and session (the per-trial file pads with NaN to 39), which participants have an auditory-control session, and that the per-network and per-region files pair infusion and preinfusion trials by participant and trial index (same labels row by row).
- Frequency axis 0 to 100 Hz in 65 bins; time axes; values are z-scores (no unit is stored).
- A clock check: the per-network preinfusion spectrograms show the puff response after t = 0 and not before.
- What the per-participant LFP and behaviour files hold, for figs15, figs27 and figs12.
C = ie.load_channel_means()
chs = C["channels"]; n_ch = len(chs)
pts = np.array([c.split("_")[1] for c in chs]); regs = np.array([c.split("_")[0] for c in chs])
print("per-channel file:", {k: C[k].shape for k in ie.SPLITS})
print("value range (z):", {k: (round(float(C[k].min()), 2), round(float(C[k].max()), 2)) for k in ie.SPLITS})
print("NaN values:", sum(int(np.isnan(C[k]).sum()) for k in ie.SPLITS))
print("frequencies (Hz):", C["ff"][0], "to", C["ff"][-1], len(C["ff"]), "bins, step", np.diff(C["ff"])[0])
print("time (s, from the first bin after the puff):", C["tt"][0], "to", round(C["tt"][-1], 3), len(C["tt"]), "bins")
print("channels:", n_ch, "participants:", dict(sorted(Counter(pts).items())))
print("regions:", len(set(regs)), "all in the 21 regions of Fig. 3C:", set(regs) <= set(ie.PLOT_REGIONS))
print("channels per region:", dict(Counter(regs).most_common()))
s_tr, names_tr, tt_tr = ie.open_trial_file()
pos_tr = {c: i for i, c in enumerate(names_tr)}
print("\nper-trial file: channels", len(names_tr), "by participant", dict(sorted(Counter(c.split("_")[1] for c in names_tr).items())))
print(" shape per session", s_tr["puff_sgram_preinfusion"].shape, "time", tt_tr[0], "to", tt_tr[-1], "s,", len(tt_tr), "bins")
print(" the 458 channels are all in it:", set(chs) <= set(names_tr),
"; labels outside the 21 regions (white matter, mixed labels) among the other", len(names_tr) - n_ch)
TIND = np.where((tt_tr > 0) & (tt_tr < 1.4))[0]
print(" NMF bins 0 < t < 1.4 s:", TIND[0], "to", TIND[-1], len(TIND), "bins; first bin", round(tt_tr[TIND[0]], 4),
"s, so the per-channel file's time axis is t - ", round(tt_tr[TIND[0]], 4))
rows = []
for pt in sorted(set(pts)):
i = np.where(pts == pt)[0][0]
row = dict(participant=pt, channel=chs[i])
for sess, key in [("preinfusion", "pre"), ("infusion", "inf"), ("postinfusion", "post")]:
x = ie.read_channel_trials(s_tr, pos_tr[chs[i]], sess)
row[f"n_{key}"] = x.shape[0]
d_tr = np.abs(x[1::2][:, :, TIND].mean(0) - C[f"{key}_train"][i]).max()
d_te = np.abs(x[0::2][:, :, TIND].mean(0) - C[f"{key}_test"][i]).max()
row[f"maxdiff_{key}"] = max(d_tr, d_te)
row["n_control"] = ie.read_channel_trials(s_tr, pos_tr[chs[i]], "control").shape[0]
rows.append(row)
trials = pd.DataFrame(rows).set_index("participant")
print(trials.to_string())
print("total trials per session:", trials[["n_pre", "n_inf", "n_post", "n_control"]].sum().to_dict())
CTRL_PTS = trials.index[trials.n_control > 0].tolist()
print("participants with an auditory control:", CTRL_PTS)
per-channel file: {'pre_train': (458, 65, 28), 'pre_test': (458, 65, 28), 'post_train': (458, 65, 28), 'post_test': (458, 65, 28), 'inf_train': (458, 65, 28), 'inf_test': (458, 65, 28)}
value range (z): {'pre_train': (-2.04, 3.98), 'pre_test': (-1.96, 3.92), 'post_train': (-2.16, 3.94), 'post_test': (-1.91, 4.09), 'inf_train': (-2.11, 3.9), 'inf_test': (-1.96, 3.94)}
NaN values: 0
frequencies (Hz): 0.0 to 100.0 65 bins, step 1.5625
time (s, from the first bin after the puff): 0.0 to 1.367 28 bins
channels: 458 participants: {np.str_('SD168'): 86, np.str_('SD172'): 97, np.str_('SD181'): 58, np.str_('SD185'): 43, np.str_('SD186'): 47, np.str_('SD211'): 44, np.str_('SD245'): 83}
regions: 21 all in the 21 regions of Fig. 3C: True
channels per region: {np.str_('PMC'): 53, np.str_('INS'): 52, np.str_('IFG'): 45, np.str_('MTG'): 40, np.str_('HIPP'): 33, np.str_('ORB'): 30, np.str_('SMG'): 27, np.str_('ITG'): 23, np.str_('STS'): 20, np.str_('STG'): 20, np.str_('THAL POS'): 17, np.str_('TP'): 16, np.str_('MCC'): 15, np.str_('ACC'): 12, np.str_('AMY'): 11, np.str_('THAL ANT'): 9, np.str_('SS'): 9, np.str_('BG'): 7, np.str_('MOT'): 7, np.str_('MFG'): 7, np.str_('PHG'): 5}
per-trial file: channels 795 by participant {'SD168': 152, 'SD172': 139, 'SD181': 104, 'SD185': 85, 'SD186': 71, 'SD211': 84, 'SD245': 160}
shape per session (795, 39, 65, 80) time -2.0 to 2.0 s, 80 bins
the 458 channels are all in it: True ; labels outside the 21 regions (white matter, mixed labels) among the other 337
NMF bins 0 < t < 1.4 s: 40 to 67 28 bins; first bin 0.0253 s, so the per-channel file's time axis is t - 0.0253
channel n_pre maxdiff_pre n_inf maxdiff_inf n_post maxdiff_post n_control
participant
SD168 PMC_SD168_LPC1-LPC2 24 0.0 24 0.0 24 0.0 0
SD172 THAL ANT_SD172_LaTH1-LaTH2 39 0.0 39 0.0 39 0.0 0
SD181 ACC_SD181_RaCN1-RaCN2 30 0.0 30 0.0 30 0.0 30
SD185 BG_SD185_RThM4-RThM5 31 0.0 31 0.0 31 0.0 31
SD186 BG_SD186_LaTH5-LaTH6 12 0.0 12 0.0 12 0.0 12
SD211 THAL ANT_SD211_LANT1-LANT2 39 0.0 39 0.0 39 0.0 39
SD245 MCC_SD245_LPCn1-LPCn2 39 0.0 39 0.0 39 0.0 39
total trials per session: {'n_pre': 214, 'n_inf': 214, 'n_post': 214, 'n_control': 151}
participants with an auditory control: [np.str_('SD181'), np.str_('SD185'), np.str_('SD186'), np.str_('SD211'), np.str_('SD245')]
YEO = ie.load_group_file(ie.F_YEO)
REG = ie.load_group_file(ie.F_REG)
for lab, G in [("network", YEO), ("region", REG)]:
print(f"per-{lab} file: {len(G['units'])} units; time {G['tt'][0]:.3f} to {G['tt'][-1]:.3f} s ({len(G['tt'])} bins); "
f"freq {G['ff'][0]} to {G['ff'][-1]} Hz ({len(G['ff'])} bins)")
paired = all((G["sess_pts"][(u, "preinfusion")] == G["sess_pts"][(u, "infusion")]).all() and
(G["sess_tri"][(u, "preinfusion")] == G["sess_tri"][(u, "infusion")]).all() for u in G["units"])
print(" infusion and preinfusion rows carry the same participant and trial index in every unit:", paired)
print("network file time bins against per-trial bins 29 to 68, max difference (s):", round(float(np.abs(YEO["tt"] - tt_tr[29:69]).max()), 5))
ntr = pd.DataFrame({u: {"trials": len(YEO["sess_pts"][(u, "preinfusion")]),
"participants": len(set(YEO["sess_pts"][(u, "preinfusion")]))} for u in YEO["units"]}).T
print("trials per network (preinfusion):\n", ntr.to_string())
ntr_r = pd.DataFrame({u: {"trials": len(REG["sess_pts"][(u, "preinfusion")]),
"participants": len(set(REG["sess_pts"][(u, "preinfusion")]))} for u in REG["units"]}).T
print("trials per region (preinfusion), the 21 of Fig. 3C:\n", ntr_r.loc[ie.PLOT_REGIONS].T.to_string())
print("other regions in the file:", [u for u in REG["units"] if u not in ie.PLOT_REGIONS])
# clock check: mean |z| over trials and frequencies before and after the puff, preinfusion networks
pre_y = YEO["h"]["scratch/puff_region_sgram_preinfusion"][:]
m = np.nanmean(pre_y, axis=1) # networks x freq x time
before = np.abs(m[:, :, YEO["tt"] < 0]).mean((1, 2)); after = np.abs(m[:, :, (YEO["tt"] > 0) & (YEO["tt"] < 1.4)]).mean((1, 2))
print("clock check, mean |z| of the trial mean, before / after the puff:")
for u, b, a in zip(YEO["units"], before, after):
print(f" {u:18s} {b:.3f} / {a:.3f}")
per-network file: 6 units; time -0.531 to 1.444 s (40 bins); freq 0.0 to 100.0 Hz (65 bins)
infusion and preinfusion rows carry the same participant and trial index in every unit: True
per-region file: 35 units; time -0.531 to 1.444 s (40 bins); freq 0.0 to 100.0 Hz (65 bins)
infusion and preinfusion rows carry the same participant and trial index in every unit: True
network file time bins against per-trial bins 29 to 68, max difference (s): 0.00063
trials per network (preinfusion):
trials participants
Frontoparietal 178 5
Default 214 7
Dorsal Attention 153 5
Limbic 214 7
Ventral Attention 214 7
Somatomotor 214 7
trials per region (preinfusion), the 21 of Fig. 3C:
ORB IFG MFG INS ACC MCC PMC MOT SS SMG PHG STG STS MTG ITG TP HIPP BG AMY THAL ANT THAL POS
trials 151 175 61 183 61 70 175 61 61 153 102 114 102 214 120 108 214 151 105 171 184
participants 5 6 2 6 2 2 6 2 2 5 3 4 3 7 4 3 7 5 4 5 6
other regions in the file: ['AG', 'CS', 'FG', 'FO', 'IPS', 'MFC', 'MOG', 'MSFG', 'OP', 'PRECG', 'SFG', 'SPL', 'THAL MID', 'VMPFC']
clock check, mean |z| of the trial mean, before / after the puff: Frontoparietal 0.057 / 0.103 Default 0.042 / 0.062 Dorsal Attention 0.072 / 0.104 Limbic 0.042 / 0.094 Ventral Attention 0.085 / 0.173 Somatomotor 0.093 / 0.186
# Per-participant LFP and behaviour files: what they hold (for figs15, figs27, figs12)
assets = dandi_io.list_assets()
ec = sorted(p for p in assets if p.startswith("sub-SD") and p.endswith("_ecephys.nwb"))
rows = []
for p in ec:
h = dandi_io.open_h5(p)
v = h["processing/ecephys/LFP/v_trace"]
rate = v["starting_time"].attrs["rate"]
rows.append(dict(file=p.split("/")[1], minutes=v["data"].shape[0] / rate / 60, rate=rate,
channels=v["data"].shape[1], description=h["session_description"][()].decode()[:60],
intervals="intervals" in h, epochs="intervals/epochs" in h,
top_groups=",".join(k for k in h.keys() if k not in ("specifications", "file_create_date", "general",
"identifier", "session_description",
"session_start_time", "timestamps_reference_time"))))
EC = pd.DataFrame(rows)
print(EC.to_string())
print("durations (min):", sorted(set(EC.minutes.round(2))))
b = dandi_io.open_h5("sub-SD245/sub-SD245_ses-puff_behavior.nwb")
print("\nbehaviour file sub-SD245_ses-puff_behavior.nwb, acquisition:", list(b["acquisition"].keys()))
print("top-level groups:", [k for k in b.keys() if k != "specifications"], "processing:", list(b["processing"].keys()))
print("any name with 'blink' or 'cadss' in any behaviour file:")
hits = []
for p in sorted(x for x in assets if x.startswith("sub-SD") and x.endswith("_behavior.nwb")):
h = dandi_io.open_h5(p)
h.visit(lambda name: hits.append((p, name)) if ("blink" in name.lower() or "cadss" in name.lower()) else None)
print(hits if hits else "none")
file minutes rate channels description intervals epochs top_groups 0 sub-SD168_ses-puffBaseline_ecephys.nwb 20.0 200.0 168 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 1 sub-SD168_ses-puffInfusion_ecephys.nwb 20.0 200.0 168 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 2 sub-SD168_ses-puffPostinfusion_ecephys.nwb 20.0 200.0 168 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 3 sub-SD172_ses-puffBaseline_ecephys.nwb 20.0 200.0 174 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 4 sub-SD172_ses-puffInfusion_ecephys.nwb 20.0 200.0 174 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 5 sub-SD172_ses-puffPostinfusion_ecephys.nwb 20.0 200.0 174 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 6 sub-SD181_ses-puffControl_ecephys.nwb 20.0 200.0 122 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 7 sub-SD181_ses-puffInfusion_ecephys.nwb 20.0 200.0 122 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 8 sub-SD181_ses-puffPostinfusion_ecephys.nwb 20.0 200.0 122 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 9 sub-SD181_ses-puffPreinfusion_ecephys.nwb 20.0 200.0 122 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 10 sub-SD185_ses-puffBaseline_ecephys.nwb 20.0 200.0 116 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 11 sub-SD185_ses-puffControl_ecephys.nwb 20.0 200.0 116 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 12 sub-SD185_ses-puffInfusion_ecephys.nwb 20.0 200.0 116 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 13 sub-SD185_ses-puffPostinfusion_ecephys.nwb 20.0 200.0 116 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 14 sub-SD185_ses-puffPreinfusion_ecephys.nwb 20.0 200.0 116 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 15 sub-SD186_ses-puffControl_ecephys.nwb 20.0 200.0 88 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 16 sub-SD186_ses-puffInfusion_ecephys.nwb 20.0 200.0 88 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 17 sub-SD186_ses-puffPostinfusion_ecephys.nwb 20.0 200.0 88 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 18 sub-SD186_ses-puffPreinfusion_ecephys.nwb 20.0 200.0 88 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 19 sub-SD211_ses-puffBaseline_ecephys.nwb 20.0 200.0 112 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 20 sub-SD211_ses-puffControl_ecephys.nwb 20.0 200.0 112 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 21 sub-SD211_ses-puffInfusion_ecephys.nwb 20.0 200.0 112 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 22 sub-SD211_ses-puffPostinfusion_ecephys.nwb 20.0 200.0 112 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 23 sub-SD211_ses-puffPostinfusionfast_ecephys.nwb 20.0 200.0 112 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 24 sub-SD211_ses-puffPreinfusion_ecephys.nwb 20.0 200.0 112 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 25 sub-SD245_ses-puffBaseline_ecephys.nwb 20.0 200.0 206 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 26 sub-SD245_ses-puffControl_ecephys.nwb 20.0 200.0 206 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 27 sub-SD245_ses-puffInfusion_ecephys.nwb 20.0 200.0 206 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 28 sub-SD245_ses-puffPostinfusion_ecephys.nwb 20.0 200.0 206 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus 29 sub-SD245_ses-puffPreinfusion_ecephys.nwb 20.0 200.0 206 Voltage traces of multichannel SEEG around eyepuffs, 10min p False False acquisition,analysis,processing,stimulus durations (min): [20.0]
behaviour file sub-SD245_ses-puff_behavior.nwb, acquisition: ['aff_w', 'affective-baseline_puffeye', 'affective-infusion_puffeye', 'affective-postinfusion_puffeye', 'affective_normed-baseline_puffeye', 'affective_normed-infusion_puffeye', 'affective_normed-postinfusion_puffeye', 'by_trial_closure-baseline_puffeye', 'by_trial_closure-infusion_puffeye', 'by_trial_closure-postinfusion_puffeye', 'events-baseline_puffeye-off', 'events-baseline_puffeye-on', 'events-infusion_puffeye-off', 'events-infusion_puffeye-on', 'events-postinfusion_puffeye-off', 'events-postinfusion_puffeye-on', 'eyeclosure-baseline_puffeye', 'eyeclosure-infusion_puffeye', 'eyeclosure-postinfusion_puffeye', 'fps', 'ref_w', 'reflexive-baseline_puffeye', 'reflexive-infusion_puffeye', 'reflexive-postinfusion_puffeye', 'which_trials'] top-level groups: ['acquisition', 'analysis', 'file_create_date', 'general', 'identifier', 'processing', 'session_description', 'session_start_time', 'stimulus', 'timestamps_reference_time'] processing: [] any name with 'blink' or 'cadss' in any behaviour file:
none
The per-participant LFP files are 20-min recordings (10 min before and after the first puff of each session) at 200 Hz; they hold no infusion start time, no epochs or annotations of sleep or seizure activity, and the behaviour files hold puff-aligned eye-closure traces, puff events and trial lists but no spontaneous blink onsets and no CADSS scores (printed above). The infusion file therefore covers about minutes 20 to 40 of the infusion (sections figs15 and figs27); the missing blink onsets block figs12 (end).
fig2g-nmf-six-factors¶
Non-negative matrix factorisation of trial-averaged, puff-triggered channel spectrograms yields six shared spectrotemporal factors, with explained variance 0.47 on held-out trials, from 458 channels of 7 participants.
The training set is the mean of the preinfusion and postinfusion training halves (odd trials), the test
set the same for the even trials, as in the code (compute_factors_using = 'prepost'). NMF is fitted with
k = 1 to 14 components per part (2 to 28 in all), as in fig. S17A. Held-out explained variance is read two
ways: (a) the code's, the training reconstruction compared with the test data; (b) loadings refitted
on the test data. The refit with k = 3 is compared with the paper's deposited factors and loadings.
X_tr = ((C["pre_train"] + C["post_train"]) / 2).reshape(n_ch, -1)
X_te = ((C["pre_test"] + C["post_test"]) / 2).reshape(n_ch, -1)
SHAPE = C["pre_train"].shape[1:]
t = time.time()
ev = ie.explained_variance_curve(X_tr, X_te, range(1, 15))
print(f"fits took {time.time() - t:.1f} s")
print(ev.round(3).to_string(index=False))
print("fits stopping at the 200-iteration default (not converged), positive / negative part:",
int((ev.iter_pos >= 200).sum()), "/", int((ev.iter_neg >= 200).sum()), "of", len(ev))
fits took 6.9 s k n_components r2_train r2_test_code r2_test_refit iter_pos iter_neg 1 2 0.353 0.329 0.361 10 10 2 4 0.470 0.410 0.464 40 50 3 6 0.542 0.466 0.535 110 50 4 8 0.567 0.480 0.560 120 70 5 10 0.591 0.485 0.582 140 100 6 12 0.610 0.490 0.600 170 90 7 14 0.628 0.496 0.615 200 120 8 16 0.641 0.499 0.626 200 120 9 18 0.652 0.501 0.638 200 200 10 20 0.664 0.506 0.649 200 200 11 22 0.675 0.507 0.659 200 200 12 24 0.684 0.505 0.665 200 200 13 26 0.690 0.506 0.670 200 180 14 28 0.696 0.506 0.676 200 160 fits stopping at the 200-iteration default (not converged), positive / negative part: 8 / 4 of 14
exp = claims["fig2g-nmf-six-factors"]["expected"][0]
six = ev.set_index("n_components").loc[6]
print(f"test explained variance, six components: code reading {six.r2_test_code:.3f}, refit reading {six.r2_test_refit:.3f}; "
f"expected {exp['value']} +/- {exp['tolerance']}")
gains = ev.set_index("n_components")[["r2_test_code", "r2_test_refit", "r2_train"]].diff()
print("gain from each step of two components:\n", gains.round(3).T.to_string())
for col in ["r2_test_code", "r2_test_refit", "r2_train"]:
y = ev[col].values; x = ev.n_components.values
g = np.diff(y) # gain into x[1:]
drop = g[:-1] - g[1:] # fall in gain after x[1:-1]
print(f"{col}: elbow by largest fall in gain = {x[1:-1][np.argmax(drop)]}, "
f"by largest distance from the chord = {ie.elbow(x, y):.0f}")
# the same rules over a shorter range of k, to show their dependence on it
for kmax in [8, 10, 12]:
sub = ev[ev.k <= kmax]
print(f" k up to {kmax}: chord elbow of r2_test_code = {ie.elbow(sub.n_components, sub.r2_test_code):.0f}")
mp, mn = ie.fit_nmf_pair(X_tr, 3)
cc = ie.components(mp, mn, SHAPE)
FAC, LOAD = ie.load_deposited_nmf()
print("\nrefit (k = 3 per part) against the deposited factors, in the code's order", ie.CODE_ORDER)
for i, fi in enumerate(ie.CODE_ORDER):
print(f" factor {i} = fitted {fi}: max |difference| {np.abs(cc[fi] - FAC[i]).max():.1e}, r = {np.corrcoef(cc[fi].ravel(), FAC[i].ravel())[0, 1]:.6f}")
print("deposited channel list equals the per-channel file's:", (LOAD["channels"] == chs).all())
# per-participant breakdown of the code-reading test explained variance at six components
H = np.vstack([mp.components_, -mn.components_])
rec_tr = ie.transform(mp, mn, X_tr) @ H
per_pt = {pt: round(ie.r2(X_te[pts == pt], rec_tr[pts == pt]), 3) for pt in sorted(set(pts))}
print("test explained variance (code reading) by participant:", per_pt)
fig, axs = plt.subplots(1, 7, figsize=(16, 2.4))
axs[0].plot(ev.n_components, ev.r2_test_code, ".-", label="test (code)")
axs[0].plot(ev.n_components, ev.r2_test_refit, ".-", label="test (refit)")
axs[0].plot(ev.n_components, ev.r2_train, ".-", label="train"); axs[0].legend(fontsize=6); axs[0].set_xlabel("components")
for i, fi in enumerate(ie.CODE_ORDER):
axs[i + 1].imshow(cc[fi], origin="lower", aspect="auto", cmap="PRGn", vmin=-1, vmax=1,
extent=[C["tt"][0], C["tt"][-1], C["ff"][0], C["ff"][-1]])
axs[i + 1].set_title(f"factor {i}", fontsize=8)
plt.tight_layout(); plt.show()
print("factor value range:", round(float(cc.min()), 2), round(float(cc.max()), 2), "(colour scale -1 to 1)")
test explained variance, six components: code reading 0.466, refit reading 0.535; expected 0.47 +/- 0.047 gain from each step of two components: n_components 2 4 6 8 10 12 14 16 18 20 22 24 26 28 r2_test_code NaN 0.081 0.056 0.014 0.005 0.006 0.005 0.003 0.002 0.004 0.001 -0.002 0.001 -0.001 r2_test_refit NaN 0.103 0.071 0.025 0.022 0.018 0.016 0.010 0.012 0.012 0.010 0.005 0.006 0.006 r2_train NaN 0.116 0.072 0.025 0.024 0.019 0.018 0.012 0.012 0.012 0.011 0.008 0.007 0.005 r2_test_code: elbow by largest fall in gain = 6, by largest distance from the chord = 8 r2_test_refit: elbow by largest fall in gain = 6, by largest distance from the chord = 8 r2_train: elbow by largest fall in gain = 6, by largest distance from the chord = 6 k up to 8: chord elbow of r2_test_code = 6 k up to 10: chord elbow of r2_test_code = 6 k up to 12: chord elbow of r2_test_code = 6
refit (k = 3 per part) against the deposited factors, in the code's order [0, 2, 4, 5, 1, 3]
factor 0 = fitted 0: max |difference| 6.4e-15, r = 1.000000
factor 1 = fitted 2: max |difference| 6.7e-15, r = 1.000000
factor 2 = fitted 4: max |difference| 4.0e-14, r = 1.000000
factor 3 = fitted 5: max |difference| 2.2e-13, r = 1.000000
factor 4 = fitted 1: max |difference| 6.3e-15, r = 1.000000
factor 5 = fitted 3: max |difference| 6.6e-14, r = 1.000000
deposited channel list equals the per-channel file's: True
test explained variance (code reading) by participant: {np.str_('SD168'): np.float64(0.413), np.str_('SD172'): np.float64(0.565), np.str_('SD181'): np.float64(0.34), np.str_('SD185'): np.float64(0.454), np.str_('SD186'): np.float64(0.292), np.str_('SD211'): np.float64(0.654), np.str_('SD245'): np.float64(0.525)}
factor value range: -1.04 2.13 (colour scale -1 to 1)
Open G, reading (b): factor 1 is the positive factor with the most energy above 50 Hz, factor 2 the negative factor with the most energy below 20 Hz, the others numbered within their part by peak time.
ORDER_B, desc = ie.description_order(cc, C["ff"], C["tt"])
print(pd.DataFrame(desc, index=[f"fitted {i}" for i in range(6)]).round(3).to_string())
print("reading (b) order:", ORDER_B, " code order:", ie.CODE_ORDER, " same:", ORDER_B == ie.CODE_ORDER)
share_above_50Hz share_below_20Hz peak_time_s fitted 0 0.051 0.915 0.152 fitted 1 0.186 0.571 0.456 fitted 2 0.860 0.016 0.203 fitted 3 0.646 0.087 0.506 fitted 4 0.058 0.852 0.405 fitted 5 0.099 0.065 0.304 reading (b) order: [0, 2, 4, 5, 1, 3] code order: [0, 2, 4, 5, 1, 3] same: True
Result: reproduced, resting on the choice of elbow rule. Refitting NMF (three positive and three negative components) on the pre/postinfusion training means of the 458 channels from 7 participants gives the deposited factors to within 3e-13 and the deposited loadings exactly (next sections). The held-out explained variance with six components is 0.466 under the code's reading (training reconstruction against the test data) and 0.535 with loadings refitted on the test data; the first is within 0.47 +/- 0.047. The test curve gains 0.081 and 0.056 in the steps to 4 and 6 components and 0.014 or less in every later step. The largest fall in gain is at 6 components under both readings; the point farthest from the chord is at 8 over the full range of 2 to 28 components and at 6 when the range ends at 16, 20 or 24 components. The paper names no rule; the largest fall in gain, the usual meaning of an elbow, is assumed. By participant the test explained variance runs from 0.292 (SD186) to 0.654 (SD211). Both Open G readings number the factors the same way. Of the 28 fits, 12 stop at the default 200 iterations, as the paper's code would.
link-ieeg-yeo-labels (link claim)¶
The Yeo7 network of each bipolar channel is not in the dataset. The code reads it from
anat/pts_ch_labels.pkl (key pt_elec_yeos_bipo), matching channels by participant and bipolar electrode
name. The link holds when every one of the 458 channels is found, its atlas region in the pickle equals
the region in the dataset's channel name, the dataset's own per-network electrode counts
(anatomy_data of the per-network file, by participant, network and region) agree with the pickle,
and the labels turn the deposited channel loadings into the deposited network loadings of Fig. 2H.
labels, md5, member = ie.fetch_code_labels()
print("fetched", member, "md5", md5)
yeo, anat = ie.yeo_labels(chs, labels)
print("channels found:", int((yeo != "Missing").sum()), "of", n_ch)
print("atlas region equals the channel-name region:", int(sum(a == r for a, r in zip(anat, regs))), "of", n_ch)
in_net = np.isin(yeo, ie.YEO_REGIONS)
print("channels in one of the six networks:", int(in_net.sum()), dict(Counter(yeo[in_net]).most_common()))
print("other labels (excluded from network means, as in the code):", int((~in_net).sum()), dict(Counter(yeo[~in_net]).most_common(8)), "...")
# counts by participant, network and region: dataset per-network file against the pickle
rows = []
for pt in sorted(set(pts)):
y = np.array(labels["pt_elec_yeos_bipo"][pt]); a = np.array(labels["pt_elec_anats_bipo"][pt])
for net in ie.YEO_REGIONS:
for k, v in Counter(a[y == net]).items():
rows.append(dict(patient=pt, region=net, anat=k, n_pkl=v))
cmp = YEO["anat"].merge(pd.DataFrame(rows), on=["patient", "region", "anat"], how="outer")
same = (cmp.nelec == cmp.n_pkl)
print(f"\nper-network file rows: {len(YEO['anat'])}; equal counts in the pickle: {int(same.sum())}; "
f"rows where the pickle has more electrodes: {int((cmp.n_pkl > cmp.nelec).sum() + cmp.nelec.isna().sum())}; "
f"rows where the dataset has more: {int((cmp.nelec > cmp.n_pkl).sum() + cmp.n_pkl.isna().sum())}")
print(cmp[~same].to_string(index=False))
# deposited network loadings recomputed from the deposited channel loadings with these labels
dep_pre = ie.group_means(LOAD["pre"], yeo, ie.YEO_REGIONS, normalize_row=True)
dep_d = ie.group_means(LOAD["inf"] - LOAD["pre"], yeo, ie.YEO_REGIONS)
print("deposited network order:", LOAD["yeo"])
print("max |difference|, Fig. 2H matrix:", float(np.nanmax(np.abs(dep_pre.values - LOAD["yeo_pre"]))),
"; Fig. 2J matrix:", float(np.nanmax(np.abs(dep_d.values - LOAD["yeo_inf_minus_pre"]))))
fetched kauvar-richman-liu-2025_zenodo/anat/pts_ch_labels.pkl md5 bd1fb567121aaf0291337783df57bcfc
channels found: 458 of 458
atlas region equals the channel-name region: 458 of 458
channels in one of the six networks: 221 {np.str_('Default'): 72, np.str_('Ventral Attention'): 57, np.str_('Somatomotor'): 37, np.str_('Limbic'): 28, np.str_('Frontoparietal'): 19, np.str_('Dorsal Attention'): 8}
other labels (excluded from network means, as in the code): 237 {np.str_('Depth'): 82, np.str_('Depth/Default'): 19, np.str_('Default/Depth'): 19, np.str_('Depth/Ventral Attention'): 12, np.str_('Ventral Attention/Depth'): 10, np.str_('Depth/Somatomotor'): 8, np.str_('Depth/Frontoparietal'): 7, np.str_('Frontoparietal/Depth'): 7} ...
per-network file rows: 126; equal counts in the pickle: 120; rows where the pickle has more electrodes: 14; rows where the dataset has more: 0
region patient anat nelec n_pkl
Default SD172 PMC 6.0 7
Frontoparietal SD172 IFG PMC NaN 1
Default SD181 TP 4.0 6
Frontoparietal SD181 IFG EXCLUDE NaN 1
Ventral Attention SD185 EXCLUDE NaN 1
Ventral Attention SD185 INS NaN 1
Ventral Attention SD185 INS MOT NaN 1
Ventral Attention SD185 MOT 2.0 3
Ventral Attention SD186 SMG 2.0 3
Default SD211 MTG 2.0 5
Default SD245 EXCLUDE AREA NaN 2
Default SD245 STG NaN 1
Dorsal Attention SD245 IPS 1.0 3
Dorsal Attention SD245 PMC NaN 1
deposited network order: ['Frontoparietal', 'Default', 'Dorsal Attention', 'Limbic', 'Ventral Attention', 'Somatomotor']
max |difference|, Fig. 2H matrix: 0.0 ; Fig. 2J matrix: 0.0
Result: the link holds. All 458 channels are found in the code record's label file, and their atlas region equals the region in the dataset's channel name for 458 of 458. Of the 126 participant, network and region counts in the dataset's per-network file, 120 equal the label file's; in the other 14 rows the label file has more electrodes (including labels marked EXCLUDE), never fewer, as expected if channels removed before the dataset's aggregates were built are still listed in the label file. With these labels the deposited channel loadings reproduce the deposited Fig. 2H and Fig. 2J matrices exactly (maximum difference 0.0). Yeo7 labels are missing from the dataset for every channel; only 221 of the 458 channels carry a single network label and enter the network means.
fig2h-factor-network-loading¶
Factors 0 to 2 load mostly on attentional and somatomotor networks, and factors 3 to 5 mostly on limbic, default-mode and frontoparietal networks.
Preinfusion test loadings, thresholded as in the code, averaged over each network's channels, each row divided by its sum. Both Open G readings give the same numbering (above), so one table serves both.
W_raw = {"pre": ie.transform(mp, mn, C["pre_test"].reshape(n_ch, -1)),
"inf": ie.transform(mp, mn, C["inf_test"].reshape(n_ch, -1)),
"post": ie.transform(mp, mn, C["post_test"].reshape(n_ch, -1))}
W_pp = ie.transform(mp, mn, ((C["pre_test"] + C["post_test"]) / 2).reshape(n_ch, -1))
Wt, EXCL = ie.threshold_loadings(W_raw, W_pp)
W = {k: v[:, ie.CODE_ORDER] for k, v in Wt.items()} # columns are factors 0..5
EXCL = EXCL[:, ie.CODE_ORDER]
print("refit loadings equal the deposited ones (NaN in the same places):",
all(np.allclose(W[k], LOAD[k], equal_nan=True, atol=1e-10) for k in ["pre", "inf", "post"]))
print("channels kept per factor (loading of the pre/post test mean >= 0.1):", (~EXCL).sum(0))
M = ie.group_means(W["pre"], yeo, ie.YEO_REGIONS, normalize_row=True)
M.columns = [f"f{i}" for i in range(6)]
M["share f0-f2"] = M[["f0", "f1", "f2"]].sum(1); M["share f3-f5"] = M[["f3", "f4", "f5"]].sum(1)
M["channels"] = [int((yeo == r).sum()) for r in ie.YEO_REGIONS]
print(M.round(3).to_string())
want = {"Dorsal Attention": "share f0-f2", "Ventral Attention": "share f0-f2", "Somatomotor": "share f0-f2",
"Limbic": "share f3-f5", "Default": "share f3-f5", "Frontoparietal": "share f3-f5"}
ok = {r: bool(M.loc[r, c] > 0.5) for r, c in want.items()}
print("rows meeting the criterion:", ok, "->", sum(ok.values()), "of 6")
# participant breakdown: share of factors 0-2 per network, per participant
rows = {}
for pt in sorted(set(pts)):
sel = pts == pt
Mp = ie.group_means(W["pre"][sel], yeo[sel], ie.YEO_REGIONS, normalize_row=True)
rows[pt] = Mp.iloc[:, :3].sum(1).where(Mp.notna().any(axis=1))
print("share of factors 0-2 by participant (NaN: no channel):\n", pd.DataFrame(rows).round(2).to_string())
refit loadings equal the deposited ones (NaN in the same places): True
channels kept per factor (loading of the pre/post test mean >= 0.1): [312 260 185 158 370 378]
f0 f1 f2 f3 f4 f5 share f0-f2 share f3-f5 channels
g
Frontoparietal 0.303 0.144 0.090 0.123 0.158 0.182 0.537 0.463 19
Default 0.149 0.205 0.111 0.131 0.177 0.228 0.464 0.536 72
Dorsal Attention 0.264 0.195 0.212 0.069 0.109 0.150 0.671 0.329 8
Limbic 0.202 0.208 0.095 0.150 0.155 0.190 0.505 0.495 28
Ventral Attention 0.256 0.243 0.152 0.098 0.128 0.123 0.651 0.349 57
Somatomotor 0.271 0.251 0.140 0.116 0.116 0.106 0.662 0.338 37
rows meeting the criterion: {'Dorsal Attention': True, 'Ventral Attention': True, 'Somatomotor': True, 'Limbic': False, 'Default': True, 'Frontoparietal': False} -> 4 of 6
share of factors 0-2 by participant (NaN: no channel):
SD168 SD172 SD181 SD185 SD186 SD211 SD245
g
Frontoparietal NaN 0.57 0.15 0.71 NaN 0.37 0.00
Default 0.44 0.48 0.15 0.26 0.32 0.14 0.53
Dorsal Attention 0.47 0.82 NaN NaN NaN 1.00 NaN
Limbic 0.73 0.53 0.33 0.47 0.51 0.27 0.41
Ventral Attention 0.69 0.77 0.56 0.84 0.66 0.82 0.66
Somatomotor 0.62 0.77 0.41 0.59 0.78 0.87 0.65
Result: partially reproduced. Factors 0 to 2 carry more than half of the dorsal attention (0.671), ventral attention (0.651) and somatomotor (0.662) rows, and factors 3 to 5 more than half of the default row (0.536). They do not in the limbic row (factors 3 to 5 carry 0.495) or the frontoparietal row (0.463). Four of the six rows meet the criterion, the partial range. The loadings are identical to the deposited ones and the matrix to the deposited Fig. 2H matrix, so the two failing rows are in the paper's own numbers. Dorsal attention rests on 8 channels; by participant the share of factors 0 to 2 in the limbic row runs from 0.27 to 0.73.
fig2i-ketamine-reduces-factors-1-2¶
Factors 1 and 2 are reduced during ketamine and recover after infusion; in the comparison of infusion with the mean of preinfusion and postinfusion, factors 0, 4 and 5 show no significant change.
Mixed model of loading on condition with channels grouped by participant (statsmodels mixedlm), infusion
against the mean of pre and post, Benjamini-Hochberg across the six factors; post-hoc preinfusion against
infusion and postinfusion against infusion, corrected together (12 tests, as the code's
comparison_list = ['pre', 'post'] does). Coefficients are infusion minus the comparison.
from statsmodels.stats.multitest import multipletests
main, post_hoc = [], []
for f in range(6):
r = ie.mixedlm_compare((W["pre"][:, f] + W["post"][:, f]) / 2, W["inf"][:, f], pts); r["factor"] = f; main.append(r)
for lab, ref in [("pre", W["pre"][:, f]), ("post", W["post"][:, f])]:
r = ie.mixedlm_compare(ref, W["inf"][:, f], pts); r.update(factor=f, against=lab); post_hoc.append(r)
main = pd.DataFrame(main).set_index("factor"); main["p_bh"] = multipletests(main.p, method="fdr_bh")[1]
post_hoc = pd.DataFrame(post_hoc); post_hoc["p_bh"] = multipletests(post_hoc.p, method="fdr_bh")[1]
means = pd.DataFrame({k: np.nanmean(W[k], 0) for k in ["pre", "inf", "post"]})
main = main.join(means)
fmt = lambda x: f"{x:.3g}"
print("infusion against mean of pre and post:\n", main.to_string(float_format=fmt))
print("\npost hoc (coefficient = infusion - comparison):\n", post_hoc.set_index(["factor", "against"]).to_string(float_format=fmt))
print("\nkept channels by participant and factor:\n", pd.DataFrame({f"f{f}": pd.Series(~np.isnan(W["pre"][:, f])).groupby(pts).sum() for f in range(6)}).to_string())
exp = claims["fig2i-ketamine-reduces-factors-1-2"]["expected"]
meas = [main.loc[1, "pre"], main.loc[1, "inf"], main.loc[1, "post"], main.loc[2, "pre"], main.loc[2, "inf"], main.loc[2, "post"]]
for e, m_ in zip(exp, meas):
print(f" {e['quantity']}: {m_:.3f} (expected {e['value']} +/- {e['tolerance']}) {'within' if abs(m_ - e['value']) <= e['tolerance'] + 1e-9 else 'outside'}")
print("\nmean loading by participant and session, factors 1 and 2:")
for f in [1, 2]:
tab = pd.DataFrame({k: pd.Series(W[k][:, f]).groupby(pts).mean() for k in ["pre", "inf", "post"]})
tab["n_ch"] = pd.Series(~np.isnan(W["pre"][:, f])).groupby(pts).sum()
print(f"factor {f}:\n", tab.round(3).to_string())
print("\nparticipants with kept channels per factor, here against the claims file's statistics (Table S3):")
st = {s["about"].split(",")[0]: s["n"] for s in claims["fig2i-ketamine-reduces-factors-1-2"]["statistics"] if "n" in s}
for f in range(6):
kept = pd.Series(~np.isnan(W["pre"][:, f])).groupby(pts).sum()
print(f" factor {f}: {int((kept > 0).sum())} participants, {int(kept.sum())} channels, fewest per participant "
f"{int(kept.min())} ({kept.idxmin()}); Table S3: {st[f'factor {f}']}")
infusion against mean of pre and post:
coef p n_ch n_pt converged p_bh pre inf post
factor
0 -0.0489 0.0808 312 7 True 0.121 0.482 0.417 0.45
1 -0.125 8.55e-06 260 7 True 2.56e-05 0.508 0.318 0.379
2 -0.115 9.18e-09 185 7 True 5.51e-08 0.354 0.214 0.302
3 -0.0407 0.00793 158 7 True 0.0159 0.282 0.214 0.228
4 0.0132 0.406 370 7 True 0.406 0.324 0.387 0.423
5 0.0132 0.285 378 7 True 0.341 0.392 0.43 0.441
post hoc (coefficient = infusion - comparison):
coef p n_ch n_pt converged p_bh
factor against
0 pre -0.0651 0.021 312 7 True 0.036
post -0.0326 0.253 312 7 True 0.304
1 pre -0.19 1.5e-10 260 7 True 1.8e-09
post -0.0605 0.0314 260 7 True 0.0464
2 pre -0.14 2.08e-09 185 7 True 1.25e-08
post -0.0888 1.29e-05 185 7 True 5.16e-05
3 pre -0.0675 0.000215 158 7 True 0.000516
post -0.0139 0.388 158 7 False 0.394
4 pre 0.0629 9.96e-05 370 7 True 0.000299
post -0.0366 0.0348 370 7 True 0.0464
5 pre 0.0378 0.00492 378 7 True 0.00985
post -0.0114 0.394 378 7 True 0.394
kept channels by participant and factor:
f0 f1 f2 f3 f4 f5
SD168 45 70 31 19 73 74
SD172 76 50 24 28 81 79
SD181 35 24 19 22 55 49
SD185 36 21 22 23 28 31
SD186 34 35 36 13 43 44
SD211 22 23 22 19 26 29
SD245 64 37 31 34 64 72
factor 1 loading preinfusion: 0.508 (expected 0.51 +/- 0.02) within
factor 1 loading infusion: 0.318 (expected 0.32 +/- 0.02) within
factor 1 loading postinfusion: 0.379 (expected 0.38 +/- 0.02) within
factor 2 loading preinfusion: 0.354 (expected 0.35 +/- 0.02) within
factor 2 loading infusion: 0.214 (expected 0.21 +/- 0.01) within
factor 2 loading postinfusion: 0.302 (expected 0.3 +/- 0.02) within
mean loading by participant and session, factors 1 and 2:
factor 1:
pre inf post n_ch
SD168 0.585 0.381 0.417 70
SD172 0.454 0.233 0.267 50
SD181 0.334 0.230 0.283 24
SD185 0.518 0.441 0.407 21
SD186 0.494 0.332 0.440 35
SD211 0.669 0.526 0.645 23
SD245 0.451 0.158 0.279 37
factor 2:
pre inf post n_ch
SD168 0.307 0.225 0.291 31
SD172 0.482 0.204 0.246 24
SD181 0.199 0.294 0.214 19
SD185 0.503 0.296 0.409 22
SD186 0.209 0.154 0.366 36
SD211 0.502 0.274 0.454 22
SD245 0.355 0.129 0.155 31
participants with kept channels per factor, here against the claims file's statistics (Table S3):
factor 0: 7 participants, 312 channels, fewest per participant 22 (SD211); Table S3: 312 channels, 7 participants
factor 1: 7 participants, 260 channels, fewest per participant 21 (SD185); Table S3: 260 channels, 7 participants
factor 2: 7 participants, 185 channels, fewest per participant 19 (SD181); Table S3: 185 channels, 7 participants
factor 3: 7 participants, 158 channels, fewest per participant 13 (SD186); Table S3: 158 channels, 5 participants
factor 4: 7 participants, 370 channels, fewest per participant 26 (SD211); Table S3: 370 channels, 6 participants
factor 5: 7 participants, 378 channels, fewest per participant 29 (SD211); Table S3: 378 channels, 7 participants
Result: reproduced. Against the mean of preinfusion and postinfusion, infusion lowers factor 1 (coefficient -0.125, corrected p 2.56e-05, 260 channels) and factor 2 (-0.115, 5.51e-08, 185 channels). Factors 0 (-0.049, 0.121), 4 (0.013, 0.406) and 5 (0.013, 0.341) are not significant; factor 3 is (-0.041, 0.0159). In this run the factor 3 fit reports not converged and every post-hoc fit converges; in the run of the second verification the factor 3 fit converged and the post-hoc factor 3 fit against postinfusion did not, with the same coefficients and p values, so the flag varies between runs. Post hoc, infusion is below postinfusion for factor 1 (-0.061, corrected p 0.0464) and factor 2 (-0.089, 5.16e-05), the recovery the claim states. The six means are within tolerance (factor 1: 0.508, 0.318, 0.379; factor 2: 0.354, 0.214, 0.302). Every coefficient, corrected p and channel count equals Table S3 to its precision; the signs of the post-hoc coefficients are reversed there (Table S3 reports the comparison minus infusion). Table S3 gives 5 participants for factor 3 and 6 for factor 4 with the same channel counts (158 and 370); here all 7 participants keep channels for every factor, the fewest being 13 (SD186) for factor 3 and 26 (SD211) for factor 4, so the difference is not a participant with few channels. Factor 1 falls from preinfusion to infusion in all 7 participants and factor 2 in 6 of 7 (SD181 rises from 0.199 to 0.294). As Table S3 also shows, the post-hoc tests are significant for factor 0 (pre, 0.036), factor 4 (pre, 0.000299; post, 0.0464) and factor 5 (pre, 0.00985).
fig2j-factor-change-by-network¶
During ketamine, factor 1 is suppressed in limbic, dorsal attention and ventral attention networks, and factor 2 in somatomotor, dorsal attention and ventral attention networks.
Change in loading (infusion minus preinfusion) averaged over each network's channels, channels kept as in fig2i. Readings of "suppressed": (a) negative in each named network, (b) the named networks are the three most negative.
D = ie.group_means(W["inf"] - W["pre"], yeo, ie.YEO_REGIONS)
D.columns = [f"f{i}" for i in range(6)]
nch = pd.DataFrame({f"n_f{f}": pd.Series(~np.isnan(W["pre"][:, f])).groupby(yeo).sum().reindex(ie.YEO_REGIONS) for f in [1, 2]})
print(D.round(3).join(nch).to_string())
named = {1: ["Limbic", "Dorsal Attention", "Ventral Attention"], 2: ["Somatomotor", "Dorsal Attention", "Ventral Attention"]}
for f, nets in named.items():
col = D[f"f{f}"]
a = all(col[n] < 0 for n in nets)
b = set(col.nsmallest(3).index) == set(nets)
print(f"factor {f}: (a) negative in all named: {a}; (b) named are the three most negative: {b} (three most negative: {list(col.nsmallest(3).index)})")
print("\nchange by participant, factor 1 and 2, named networks (NaN: no kept channel):")
for f, nets in named.items():
rows = {}
for pt in sorted(set(pts)):
sel = pts == pt
rows[pt] = ie.group_means((W["inf"] - W["pre"])[sel], yeo[sel], nets).iloc[:, f]
print(f"factor {f}:\n", pd.DataFrame(rows).round(3).to_string())
f0 f1 f2 f3 f4 f5 n_f1 n_f2
g
Frontoparietal -0.165 -0.166 0.010 -0.110 -0.059 0.115 7 2
Default -0.103 -0.269 -0.072 -0.071 0.098 -0.015 27 10
Dorsal Attention -0.113 -0.225 -0.190 0.136 0.126 -0.050 7 5
Limbic -0.147 -0.184 -0.023 -0.049 0.013 0.086 16 7
Ventral Attention -0.050 -0.266 -0.157 -0.076 -0.017 0.108 45 36
Somatomotor -0.041 -0.083 -0.108 0.054 -0.001 -0.037 33 31
factor 1: (a) negative in all named: True; (b) named are the three most negative: False (three most negative: ['Default', 'Ventral Attention', 'Dorsal Attention'])
factor 2: (a) negative in all named: True; (b) named are the three most negative: True (three most negative: ['Dorsal Attention', 'Ventral Attention', 'Somatomotor'])
change by participant, factor 1 and 2, named networks (NaN: no kept channel):
factor 1:
SD168 SD172 SD181 SD185 SD186 SD211 SD245
g
Limbic -0.586 -0.336 -0.078 -0.352 -0.146 -0.089 0.034
Dorsal Attention -0.206 -0.285 NaN NaN NaN -0.234 NaN
Ventral Attention -0.189 -0.273 -0.199 -0.295 -0.165 -0.115 -0.377
factor 2:
SD168 SD172 SD181 SD185 SD186 SD211 SD245
g
Somatomotor -0.005 -0.169 -0.052 -0.297 0.128 -0.141 -0.130
Dorsal Attention -0.031 NaN NaN NaN NaN -0.430 NaN
Ventral Attention -0.108 -0.263 0.038 -0.209 -0.096 -0.485 -0.206
Result: reproduced under reading (a). Infusion minus preinfusion, factor 1 is negative in the limbic (-0.184), dorsal attention (-0.225) and ventral attention (-0.266) networks, and factor 2 in the somatomotor (-0.108), dorsal attention (-0.190) and ventral attention (-0.157) networks. Under reading (b) factor 2 holds and factor 1 does not: its three most negative networks are default (-0.269), ventral attention and dorsal attention, with limbic fourth; factor 1 falls in all six networks (-0.083 to -0.269). The dorsal attention values rest on 7 channels from 3 participants (factor 1) and 5 channels from 2 participants (factor 2). The matrix equals the deposited Fig. 2J matrix (link section).
figs13-auditory-control-no-response¶
The brain-wide eyepuff-evoked response is absent in the auditory control, where the device sounds but the air is directed away from the participant.
The control session's trial mean (all trials, as the code's X_ctrl) of each channel of the participants
with a control session is read from the per-trial file and projected on the fig2g factors; channels
excluded for a factor in fig2i are excluded here too (code W_ctrl). The criterion compares the mean
loading over channels in the control with the preinfusion mean; both the preinfusion mean over the same
channels and over all channels are shown. The preinfusion loadings are those of the test half of the trials
(fig2i); as a variant, preinfusion loadings from all preinfusion trials, the same trial use as the control.
t = time.time()
ctrl_idx = np.where(np.isin(pts, CTRL_PTS))[0]
Xc = np.stack([ie.read_channel_trials(s_tr, pos_tr[chs[i]], "control")[:, :, TIND].mean(0) for i in ctrl_idx])
print(f"read control trial means of {len(ctrl_idx)} channels in {time.time() - t:.0f} s; NaN: {int(np.isnan(Xc).sum())}")
Wc = ie.transform(mp, mn, Xc.reshape(len(ctrl_idx), -1))[:, ie.CODE_ORDER]
Wc[EXCL[ctrl_idx]] = np.nan
res = pd.DataFrame({"control": np.nanmean(Wc, 0), "pre (same channels)": np.nanmean(W["pre"][ctrl_idx], 0),
"pre (all channels)": np.nanmean(W["pre"], 0), "n_ch": (~np.isnan(Wc)).sum(0)})
res["control lower than pre (same channels)"] = res["control"] < res["pre (same channels)"]
res["ratio control / pre"] = res["control"] / res["pre (same channels)"]
print(res.round(3).to_string())
print("factors with lower control loading:", int(res["control lower than pre (same channels)"].sum()), "of 6")
print("\nmean loading by participant, control / preinfusion (same channels):")
cp = pts[ctrl_idx]
tab = pd.DataFrame({f"f{f}": [f"{np.nanmean(Wc[cp == p, f]):.2f}/{np.nanmean(W['pre'][ctrl_idx][cp == p, f]):.2f}" for p in CTRL_PTS] for f in range(6)}, index=CTRL_PTS)
print(tab.to_string())
lower_pt = np.array([[np.nanmean(Wc[cp == p, f]) < np.nanmean(W["pre"][ctrl_idx][cp == p, f]) for f in range(6)] for p in CTRL_PTS])
print("participant-factor pairs with lower control loading:", int(lower_pt.sum()), "of", lower_pt.size)
# variant: preinfusion loadings from all preinfusion trials, as the control uses all of its trials
t = time.time()
Xp_all = np.stack([ie.read_channel_trials(s_tr, pos_tr[chs[i]], "preinfusion")[:, :, TIND].mean(0) for i in ctrl_idx])
Wp_all = ie.transform(mp, mn, Xp_all.reshape(len(ctrl_idx), -1))[:, ie.CODE_ORDER]
Wp_all[EXCL[ctrl_idx]] = np.nan
print(f"\nvariant, preinfusion from all preinfusion trials (read in {time.time() - t:.0f} s):")
res2 = pd.DataFrame({"control": np.nanmean(Wc, 0), "pre (all trials)": np.nanmean(Wp_all, 0),
"pre (test half)": np.nanmean(W["pre"][ctrl_idx], 0), "n_ch": (~np.isnan(Wc)).sum(0)})
res2["control lower than pre (all trials)"] = res2["control"] < res2["pre (all trials)"]
res2["ratio control / pre (all trials)"] = res2["control"] / res2["pre (all trials)"]
print(res2.round(3).to_string())
print("factors with lower control loading:", int(res2["control lower than pre (all trials)"].sum()), "of 6")
lower_pt2 = np.array([[np.nanmean(Wc[cp == p, f]) < np.nanmean(Wp_all[cp == p, f]) for f in range(6)] for p in CTRL_PTS])
print("participant-factor pairs with lower control loading:", int(lower_pt2.sum()), "of", lower_pt2.size)
print("pairs where it is not (participant, factor, control, pre all trials):",
[(str(p), f, round(float(np.nanmean(Wc[cp == p, f])), 2), round(float(np.nanmean(Wp_all[cp == p, f])), 2))
for j, p in enumerate(CTRL_PTS) for f in range(6) if not lower_pt2[j, f]])
read control trial means of 275 channels in 17 s; NaN: 0
control pre (same channels) pre (all channels) n_ch control lower than pre (same channels) ratio control / pre
0 0.157 0.460 0.482 191 True 0.341
1 0.229 0.488 0.508 140 True 0.469
2 0.174 0.342 0.354 130 True 0.509
3 0.104 0.315 0.282 111 True 0.329
4 0.211 0.325 0.324 216 True 0.648
5 0.299 0.429 0.392 225 True 0.698
factors with lower control loading: 6 of 6
mean loading by participant, control / preinfusion (same channels):
f0 f1 f2 f3 f4 f5
SD181 0.14/0.41 0.13/0.33 0.20/0.20 0.11/0.34 0.19/0.36 0.32/0.39
SD185 0.11/0.36 0.19/0.52 0.20/0.50 0.08/0.39 0.21/0.34 0.29/0.38
SD186 0.24/0.60 0.25/0.49 0.18/0.21 0.13/0.35 0.30/0.41 0.38/0.63
SD211 0.34/0.54 0.41/0.67 0.16/0.50 0.15/0.26 0.17/0.27 0.26/0.39
SD245 0.08/0.44 0.19/0.45 0.15/0.36 0.08/0.27 0.19/0.25 0.25/0.37
participant-factor pairs with lower control loading: 30 of 30
variant, preinfusion from all preinfusion trials (read in 13 s):
control pre (all trials) pre (test half) n_ch control lower than pre (all trials) ratio control / pre (all trials)
0 0.157 0.458 0.460 191 True 0.342
1 0.229 0.467 0.488 140 True 0.490
2 0.174 0.325 0.342 130 True 0.536
3 0.104 0.232 0.315 111 True 0.446
4 0.211 0.275 0.325 216 True 0.766
5 0.299 0.306 0.429 225 True 0.977
factors with lower control loading: 6 of 6
participant-factor pairs with lower control loading: 24 of 30
pairs where it is not (participant, factor, control, pre all trials): [('SD181', 2, 0.2, 0.16), ('SD181', 5, 0.32, 0.27), ('SD185', 5, 0.29, 0.27), ('SD186', 2, 0.18, 0.16), ('SD211', 5, 0.26, 0.24), ('SD245', 5, 0.25, 0.23)]
Result: reproduced. Over the 275 channels of the 5 participants with an auditory-control session (SD181, SD185, SD186, SD211, SD245), the mean control loading is below the preinfusion mean of the same channels for all six factors (ratios 0.33 to 0.70), and for 30 of 30 participant and factor pairs. With preinfusion loadings from all preinfusion trials (the control's trial use), the control mean is still lower for all six factors, by less (ratios 0.34 to 0.98; factor 5 0.299 against 0.306), and for 24 of 30 participant and factor pairs (the six exceptions are factors 2 and 5).
fig3b-cluster-test-yeo-networks¶
Permutation cluster tests find ketamine-altered spectrotemporal clusters in the limbic, ventral attention and somatomotor Yeo7 networks.
Per network, the infusion minus preinfusion difference of trials paired by participant and trial index
(the per-network file's rows), mne permutation_cluster_1samp_test, threshold the two-tailed t at
p = 0.01, 512 permutations, seed 1, -0.5 to 1.5 s. The lowest cluster p of each network (1 when no cluster
is below 0.01, as the code does) is Benjamini-Hochberg corrected across the six networks. The deposited
change maps and contours are compared with these. Cluster-forming thresholds of p = 0.05 and 0.001 (the
code's PARAMS 1 and 2) are shown as variants.
DEP = ie.load_deposited_clusters()
def run_units(G, units, thresh_p=0.01):
out = {}
for u in units:
j = G["units"].index(u)
pre = G["h"]["scratch/puff_region_sgram_preinfusion"][j]; inf = G["h"]["scratch/puff_region_sgram_infusion"][j]
v = ~np.isnan(pre[:, 0, 0])
r = ie.cluster_test(pre[v], inf[v], thresh_p=thresh_p)
r["pts"] = G["sess_pts"][(u, "preinfusion")]; r["diff"] = inf[v] - pre[v]
out[u] = r
return out
t = time.time()
RY = run_units(YEO, ie.YEO_REGIONS)
qy = ie.fdr_units(RY, ie.YEO_REGIONS)
tab = pd.DataFrame({u: dict(trials=RY[u]["n_trials"], clusters=RY[u]["n_clusters"], signif_clusters=RY[u]["n_signif"],
min_p=RY[u]["min_p"], p_bh=qy[u], significant=qy[u] < 0.01,
dep_contour_px=int(DEP[u][1].sum()), contour_px=int(RY[u]["contour"].sum()),
T_vs_deposited_maxdiff=float(np.abs(DEP[u][0] - RY[u]["T"]).max()))
for u in ie.YEO_REGIONS}).T
print(f"({time.time() - t:.0f} s)\n", tab.to_string())
print("min p is bounded below by 1/512 =", round(1 / 512, 5))
for thr in [0.05, 0.001]:
R2 = run_units(YEO, ie.YEO_REGIONS, thresh_p=thr); q2 = ie.fdr_units(R2, ie.YEO_REGIONS)
print(f"variant threshold p = {thr}: significant networks", [u for u in ie.YEO_REGIONS if q2[u] < 0.01])
print("\nmean infusion - preinfusion change inside the significant clusters, by participant (z):")
rows = {}
for u in ie.YEO_REGIONS:
if RY[u]["contour"].sum() == 0: continue
m = RY[u]["diff"][:, RY[u]["contour"] > 0].mean(1)
rows[u] = pd.Series(m).groupby(RY[u]["pts"]).mean()
print(pd.DataFrame(rows).round(3).to_string())
(2 s)
trials clusters signif_clusters min_p p_bh significant dep_contour_px contour_px T_vs_deposited_maxdiff
Frontoparietal 178 27 0 0.027344 1.0 False 0 0 0.0
Default 214 17 0 0.044922 1.0 False 0 0 0.0
Dorsal Attention 153 20 0 0.101562 1.0 False 0 0 0.0
Limbic 214 19 2 0.001953 0.003906 True 192 192 0.0
Ventral Attention 214 20 2 0.001953 0.003906 True 219 219 0.0
Somatomotor 214 23 1 0.001953 0.003906 True 100 100 0.0
min p is bounded below by 1/512 = 0.00195
variant threshold p = 0.05: significant networks ['Dorsal Attention', 'Limbic', 'Ventral Attention', 'Somatomotor']
variant threshold p = 0.001: significant networks ['Limbic', 'Ventral Attention']
mean infusion - preinfusion change inside the significant clusters, by participant (z):
Limbic Ventral Attention Somatomotor
SD168 -0.114 0.005 0.194
SD172 -0.343 0.024 0.136
SD181 -0.211 -0.101 0.204
SD185 -0.522 -0.049 0.257
SD186 -0.096 -0.030 0.214
SD211 -0.197 0.072 0.550
SD245 -0.234 -0.032 0.133
Result: reproduced. Limbic, ventral attention and somatomotor each have significant clusters (lowest cluster p 0.00195, the floor for 512 permutations; corrected 0.0039), and frontoparietal (lowest 0.027), default (0.045) and dorsal attention (0.10) have none below 0.01. The t maps and contours equal the deposited ones in all six networks (maximum difference 0). With a cluster-forming threshold of p = 0.05, dorsal attention is significant as well; with 0.001, somatomotor is not. Inside the significant clusters the limbic change is negative and the somatomotor change positive in each of the 7 participants; the ventral attention mean over its two clusters is near zero in every participant (-0.101 to 0.072).
fig3c-cluster-test-regions¶
Ketamine alters the eyepuff spectrotemporal response in insula, cingulate, orbitofrontal, posteromedial, supramarginal and motor cortex and in hippocampus, basal ganglia and thalamus, while many other regions are unchanged.
The fig3b test on the 21 regions of Fig. 3C, corrected across these 21 (the criterion). As variants, the correction across all 35 regions of the file and the deposited contours (33 regions) are shown. Reading (a) of the criterion is run; reading (b) (adding the single-channel test) is needed only if (a) fails.
t = time.time()
RR = run_units(REG, REG["units"])
q21 = ie.fdr_units(RR, ie.PLOT_REGIONS); q35 = ie.fdr_units(RR, REG["units"])
LISTED = ["INS", "ACC", "MCC", "ORB", "PMC", "SMG", "MOT", "HIPP", "BG", "THAL"]
tab = pd.DataFrame({u: dict(trials=RR[u]["n_trials"], participants=len(set(RR[u]["pts"])), min_p=RR[u]["min_p"],
p_bh_21=q21.get(u, np.nan), p_bh_35=q35[u],
dep_contour_px=int(DEP[u][1].sum()) if u in DEP else np.nan,
contour_px_21=int(RR[u]["contour"].sum()) if q21.get(u, 1) < 0.01 else 0,
T_vs_deposited_maxdiff=float(np.abs(DEP[u][0] - RR[u]["T"]).max()) if u in DEP else np.nan)
for u in REG["units"]}).T
print(f"({time.time() - t:.0f} s)")
print(tab.loc[ie.PLOT_REGIONS].to_string())
print("\nother regions (not in the 21):\n", tab.drop(ie.PLOT_REGIONS).to_string())
sig21 = {u for u in ie.PLOT_REGIONS if q21[u] < 0.01}
listed_ok = {l: (l in sig21) if l != "THAL" else ("THAL ANT" in sig21 or "THAL POS" in sig21) for l in LISTED}
others = [u for u in ie.PLOT_REGIONS if u not in ["INS", "ACC", "MCC", "ORB", "PMC", "SMG", "MOT", "HIPP", "BG", "THAL ANT", "THAL POS"]]
print("\nreading (a), listed regions significant:", listed_ok, "->", sum(listed_ok.values()), "of 10")
print("other regions of the 21:", others, "; not significant:", sum(u not in sig21 for u in others), "of", len(others))
print(" counting one thalamic region (ANT) among the listed and THAL POS among the others: not significant",
sum(u not in sig21 for u in others + ["THAL POS"]), "of", len(others) + 1)
print("significant with correction across 35:", sorted(u for u in REG["units"] if q35[u] < 0.01))
print("regions with a deposited contour:", sorted(u for u in DEP if DEP[u][1].sum() > 0 and u not in ie.YEO_REGIONS))
for thr in [0.05, 0.001]:
R2 = run_units(REG, ie.PLOT_REGIONS, thresh_p=thr); q2 = ie.fdr_units(R2, ie.PLOT_REGIONS)
print(f"variant threshold p = {thr}: significant of the 21:", sorted(u for u in ie.PLOT_REGIONS if q2[u] < 0.01))
print("\nmean change inside the significant clusters by participant (z), listed regions:")
rows = {}
for u in sorted(sig21):
m = RR[u]["diff"][:, RR[u]["contour"] > 0].mean(1)
rows[u] = pd.Series(m).groupby(RR[u]["pts"]).mean()
print(pd.DataFrame(rows).round(3).to_string())
(11 s)
trials participants min_p p_bh_21 p_bh_35 dep_contour_px contour_px_21 T_vs_deposited_maxdiff
ORB 151.0 5.0 0.001953 0.003729 0.005697 251.0 251.0 0.0
IFG 175.0 6.0 0.017578 1.000000 1.000000 0.0 0.0 0.0
MFG 61.0 2.0 0.046875 1.000000 1.000000 0.0 0.0 0.0
INS 183.0 6.0 0.001953 0.003729 0.005697 346.0 346.0 0.0
ACC 61.0 2.0 0.001953 0.003729 0.005697 35.0 35.0 0.0
MCC 70.0 2.0 0.001953 0.003729 0.005697 90.0 90.0 0.0
PMC 175.0 6.0 0.001953 0.003729 0.005697 205.0 205.0 0.0
MOT 61.0 2.0 0.001953 0.003729 0.005697 88.0 88.0 0.0
SS 61.0 2.0 0.007812 0.012620 0.017090 0.0 0.0 0.0
SMG 153.0 5.0 0.001953 0.003729 0.005697 303.0 303.0 0.0
PHG 102.0 3.0 0.031250 1.000000 1.000000 0.0 0.0 0.0
STG 114.0 4.0 0.005859 0.010254 0.015775 0.0 0.0 0.0
STS 102.0 3.0 0.083984 1.000000 1.000000 0.0 0.0 0.0
MTG 214.0 7.0 0.009766 0.014648 0.017989 0.0 0.0 0.0
ITG 120.0 4.0 0.447266 1.000000 1.000000 0.0 0.0 0.0
TP 108.0 3.0 0.269531 1.000000 1.000000 0.0 0.0 0.0
HIPP 214.0 7.0 0.001953 0.003729 0.005697 74.0 74.0 0.0
BG 151.0 5.0 0.001953 0.003729 0.005697 82.0 82.0 0.0
AMY 105.0 4.0 0.019531 1.000000 1.000000 0.0 0.0 0.0
THAL ANT 171.0 5.0 0.001953 0.003729 0.005697 71.0 71.0 0.0
THAL POS 184.0 6.0 0.001953 0.003729 0.005697 70.0 70.0 0.0
other regions (not in the 21):
trials participants min_p p_bh_21 p_bh_35 dep_contour_px contour_px_21 T_vs_deposited_maxdiff
AG 70.0 2.0 0.523438 NaN 1.000000 0.0 0.0 0.0
CS 39.0 1.0 0.359375 NaN 1.000000 0.0 0.0 0.0
FG 69.0 2.0 0.007812 NaN 0.017090 0.0 0.0 0.0
FO 30.0 1.0 0.011719 NaN 1.000000 0.0 0.0 0.0
IPS 39.0 1.0 0.152344 NaN 1.000000 NaN 0.0 NaN
MFC 31.0 1.0 0.152344 NaN 1.000000 0.0 0.0 0.0
MOG 39.0 1.0 0.294922 NaN 1.000000 0.0 0.0 0.0
MSFG 39.0 1.0 0.917969 NaN 1.000000 0.0 0.0 0.0
OP 39.0 1.0 0.009766 NaN 0.017989 0.0 0.0 0.0
PRECG 39.0 1.0 0.001953 NaN 0.005697 97.0 0.0 0.0
SFG 39.0 1.0 0.009766 NaN 0.017989 NaN 0.0 NaN
SPL 12.0 1.0 0.271484 NaN 1.000000 0.0 0.0 0.0
THAL MID 31.0 1.0 0.007812 NaN 0.017090 0.0 0.0 0.0
VMPFC 31.0 1.0 0.884766 NaN 1.000000 0.0 0.0 0.0
reading (a), listed regions significant: {'INS': True, 'ACC': True, 'MCC': True, 'ORB': True, 'PMC': True, 'SMG': True, 'MOT': True, 'HIPP': True, 'BG': True, 'THAL': True} -> 10 of 10
other regions of the 21: ['IFG', 'MFG', 'SS', 'PHG', 'STG', 'STS', 'MTG', 'ITG', 'TP', 'AMY'] ; not significant: 10 of 10
counting one thalamic region (ANT) among the listed and THAL POS among the others: not significant 10 of 11
significant with correction across 35: ['ACC', 'BG', 'HIPP', 'INS', 'MCC', 'MOT', 'ORB', 'PMC', 'PRECG', 'SMG', 'THAL ANT', 'THAL POS']
regions with a deposited contour: ['ACC', 'BG', 'HIPP', 'INS', 'MCC', 'MOT', 'ORB', 'PMC', 'PRECG', 'SMG', 'THAL ANT', 'THAL POS', 'Visual']
variant threshold p = 0.05: significant of the 21: ['ACC', 'BG', 'HIPP', 'INS', 'MCC', 'MOT', 'MTG', 'ORB', 'PMC', 'SMG', 'STG', 'THAL ANT', 'THAL POS']
variant threshold p = 0.001: significant of the 21: ['ACC', 'BG', 'HIPP', 'INS', 'MCC', 'MOT', 'ORB', 'PMC', 'SMG', 'THAL ANT']
mean change inside the significant clusters by participant (z), listed regions:
ACC BG HIPP INS MCC MOT ORB PMC SMG THAL ANT THAL POS
SD168 NaN NaN 0.447 0.082 NaN NaN NaN -0.093 -0.063 -0.225 -0.280
SD172 NaN 0.457 0.278 -0.105 NaN NaN -0.261 -0.163 0.169 -0.398 -0.193
SD181 -0.306 0.512 0.071 -0.116 NaN -0.380 -0.198 -0.781 NaN -0.343 NaN
SD185 -0.316 -0.046 -0.023 NaN -0.801 -0.042 -0.312 -0.143 NaN NaN -0.351
SD186 NaN 0.412 0.263 -0.020 NaN NaN -0.157 -0.134 -0.118 NaN 0.025
SD211 NaN 0.635 0.107 -0.069 NaN NaN NaN NaN -0.070 -0.295 -0.242
SD245 NaN NaN 0.500 -0.124 -0.159 NaN -0.176 -0.036 -0.043 -0.353 -0.319
Result: reproduced. With correction across the 21 regions, all ten listed regions are significant (INS, ACC, MCC, ORB, PMC, SMG, MOT, HIPP, BG, and both THAL ANT and THAL POS; lowest cluster p 0.00195, corrected 0.0037), and none of the other ten (IFG, MFG, SS, PHG, STG, STS, MTG, ITG, TP, AMY) is; counting THAL POS among the others, 10 of 11 are not significant, above the five required. The t maps equal the deposited ones in all 33 deposited regions, and the significant set equals the regions with a deposited contour, which add PRECG (one participant, outside the 21). ACC, MCC and MOT rest on 2 participants each (61 to 70 trials). At a cluster-forming threshold of 0.001 all listed regions but THAL POS stay significant; at 0.05, MTG and STG join. Reading (b) was not needed and not run.
figs15-ketamine-spectral-changes-frontal¶
Ketamine changes spectral power across the brain, with the largest changes in frontal cortex, led by orbitofrontal cortex.
An approximate test on the per-participant 20-min LFP files. Timing: the Methods (iEEG ketamine infusion protocol) say "The eyepuff assay was administered preceding infusion, at approximately 30 minutes after the infusion start, and 110 minutes after infusion start". Each file runs from 10 min before to 10 min after the first puff, so the infusion file covers about minutes 20 to 40 of the infusion, and its second half (first puff to the end) is the 30 to 40 min peak window of fig. S15C. Because the start is approximate and the infusion start time is not in the dataset, the first half (about 20 to 30 min), the whole file, and the second half without the eyepuff assay (from 2 s after the last puff) are run as well.
Baseline: the paper's 25-min pre-ketamine window does not fit in a 20-min file. The preinfusion eyepuff file
is used (the code's preinfusion session: Baseline for SD168, SD172 and SD245, Preinfusion for the others;
ieeg_cd.PRE_SESSION_CODE). Main: its first 10 min, before the first puff, since the Methods exclude the
eyepuff assay from preinfusion windows ("excluding periods of time in which the patient was undergoing the
eyepuff assay", Methods, phase locking); variants: the whole file, and SD245's Preinfusion file in place of
its Baseline file.
Spectrogram: multitaper with the parameters the Methods give for fig. S27 (time-half-bandwidth 30, 59 tapers,
30-s windows, 15-s stride); the paper gives none for fig. S15, so this is an assumption. Power is taken as
log10 (an assumption; linear power is a variant), z-scored per channel and frequency to the baseline
windows and averaged over the ketamine windows. Magnitude (criterion): per participant and region the mean
over channels of the mean |z| over 1 to 100 Hz (58 to 62 Hz, the notch, left out), then the median over
participants. Channels: within-shank bipolar pairs whose atlas label in the code record
(pt_elec_anats_bipo) is a single one of the 21 regions of Fig. 3C, without out-of-brain or excluded
contacts and without the code's bad channels (utils.bad_ch). The files carry no sleep or seizure
annotations, so no period is excluded for them. Cohorts: (a) the 7 included participants; (b) every
participant with a 40-min infusion in the dataset. Only the 7 included participants have LFP files, and the
dataset does not say which of them had a continuous 40-min infusion (the Methods mention target-controlled
and bolus infusions for some), so (b) is the same 7 or an unidentifiable subset of them; only (a) can be run.
import ieeg_cd as cd
import human_eye as he
t = time.time()
BAD = cd.code_bad_channels(cd.fetch_code_zip())
inf_pts = sorted(p.split("/")[0][4:] for p in assets if p.endswith("_ses-puffInfusion_ecephys.nwb"))
print("participants with an infusion LFP file:", inf_pts, "; the 7 included:", inf_pts == sorted(cd.PARTICIPANTS))
REST = ie.rest_all(cd.PARTICIPANTS, labels, BAD)
print(f"streamed {sum(len(d) for d in REST.values())} files and computed their spectrograms in {time.time() - t:.0f} s")
rows = []
for pt, d in REST.items():
for k, R in d.items():
rows.append(dict(participant=pt, role=k, session=R["session"], minutes=R["duration"] / 60, first_puff_s=R["first_puff"],
last_puff_s=round(R["last_puff"], 1), puffs=R["n_puffs"], bipolar_pairs=R["n_pairs"], kept=len(R["ch"]),
windows=len(R["t"]), log10_power_range=f"{np.log10(R['S'].min()):.1f} to {np.log10(R['S'].max()):.1f}"))
print(pd.DataFrame(rows).to_string(index=False))
R0 = REST[cd.PARTICIPANTS[0]]["inf"]
print("frequency bins:", R0["f"][0], "to", R0["f"][-1], "Hz, step", ie.FBIN, "; window starts", R0["t"][0], "to", R0["t"][-1], "s")
for w in ["before", "after", "after_no_puffs", "whole"]:
lo, hi = ie.window_bounds(R0, w)
print(f" {w:15s} {ie.REST_WINDOWS[w]}: {lo:.0f} to {hi:.0f} s in {R0['session']} of {cd.PARTICIPANTS[0]}, "
f"{int(ie.window_mask(R0['t'], lo, hi).sum())} windows")
REG_TAB = pd.crosstab(pd.Series(np.concatenate([REST[p]["inf"]["region"] for p in cd.PARTICIPANTS]), name="region"),
pd.Series(np.concatenate([[p] * len(REST[p]["inf"]["ch"]) for p in cd.PARTICIPANTS]), name="participant"))
print("\nkept channels by region and participant:\n", REG_TAB.to_string())
participants with an infusion LFP file: ['SD168', 'SD172', 'SD181', 'SD185', 'SD186', 'SD211', 'SD245'] ; the 7 included: True
streamed 22 files and computed their spectrograms in 94 s
participant role session minutes first_puff_s last_puff_s puffs bipolar_pairs kept windows log10_power_range
SD168 pre Baseline 20.0 600.005 939.2 40 154 86 79 -2.0 to 7.9
SD168 inf Infusion 20.0 600.005 939.2 40 154 86 79 -1.9 to 6.7
SD168 post Postinfusion 20.0 600.005 939.1 40 154 86 79 -2.0 to 6.8
SD172 pre Baseline 20.0 600.005 939.1 40 162 97 79 -2.6 to 7.9
SD172 inf Infusion 20.0 600.005 939.1 40 162 97 79 -2.6 to 7.3
SD172 post Postinfusion 20.0 600.005 939.1 40 162 97 79 -2.7 to 7.9
SD181 pre Preinfusion 20.0 600.005 939.3 40 113 58 79 -3.1 to 5.1
SD181 inf Infusion 20.0 600.005 939.1 40 113 58 79 -2.9 to 4.6
SD181 post Postinfusion 20.0 600.010 939.1 40 113 58 79 -3.1 to 4.7
SD185 pre Preinfusion 20.0 600.010 939.2 40 106 43 79 -3.0 to 8.0
SD185 inf Infusion 20.0 600.005 861.0 32 106 43 79 -2.9 to 7.2
SD185 post Postinfusion 20.0 600.005 939.1 40 106 43 79 -3.1 to 7.5
SD186 pre Preinfusion 20.0 600.015 939.2 40 81 47 79 -2.6 to 5.1
SD186 inf Infusion 20.0 600.055 861.1 32 81 47 79 -2.5 to 4.9
SD186 post Postinfusion 20.0 600.060 861.0 32 81 47 79 -2.6 to 5.2
SD211 pre Preinfusion 20.0 600.010 939.2 40 105 44 79 -3.5 to 5.0
SD211 inf Infusion 20.0 600.005 939.2 40 105 44 79 -3.3 to 4.8
SD211 post Postinfusion 20.0 600.010 939.2 40 105 44 79 -3.3 to 4.8
SD245 pre Baseline 20.0 600.005 939.2 40 190 83 79 -2.9 to 5.8
SD245 inf Infusion 20.0 600.010 939.2 40 190 83 79 -2.9 to 5.6
SD245 post Postinfusion 20.0 600.010 939.2 40 190 83 79 -2.9 to 6.1
SD245 pre_alt Preinfusion 20.0 600.010 939.2 40 190 83 79 -3.0 to 5.4
frequency bins: 0.5 to 99.5 Hz, step 0.5 ; window starts 0.0 to 1170.0 s
before first half of the file (0 s to the first puff): 0 to 600 s in Infusion of SD168, 39 windows
after second half (first puff to the end), about 30 to 40 min of the infusion: 600 to 1200 s in Infusion of SD168, 39 windows
after_no_puffs second half from 2 s after the last puff (fig. S15 only): 941 to 1200 s in Infusion of SD168, 16 windows
whole the whole 20-min file: 0 to 1200 s in Infusion of SD168, 79 windows
kept channels by region and participant:
participant SD168 SD172 SD181 SD185 SD186 SD211 SD245
region
ACC 0 0 8 4 0 0 0
AMY 5 1 3 0 2 0 0
BG 0 2 2 1 1 1 0
HIPP 13 3 4 3 4 4 2
IFG 3 13 11 1 7 0 10
INS 3 20 5 0 2 3 19
ITG 0 5 5 0 2 11 0
MCC 0 0 0 1 0 0 14
MFG 0 0 3 4 0 0 0
MOT 0 0 2 5 0 0 0
MTG 9 13 2 5 4 6 1
ORB 0 7 5 7 7 0 4
PHG 3 1 0 0 0 1 0
PMC 18 12 1 3 5 0 14
SMG 2 4 0 0 5 8 8
SS 0 0 1 8 0 0 0
STG 10 1 0 0 5 0 4
STS 14 5 0 0 0 0 1
THAL ANT 4 2 1 0 0 1 1
THAL POS 2 4 0 1 3 2 5
TP 0 4 5 0 0 7 0
rows, F15 = [], {}
for base_key in ["pre", "pre_alt"]:
for bw in ["before", "whole"]:
for kw in ["after", "after_no_puffs", "before", "whole"]:
for log in [True, False]:
if base_key == "pre_alt" and not log:
continue
mag = ie.s15_magnitudes(REST, base_key, bw, kw, log=log)
v = ie.s15_verdict(mag)
F15[(base_key, bw, kw, log)] = (mag, v)
med = v["median"]
rows.append(dict(baseline=("code session" if base_key == "pre" else "SD245 Preinfusion") + ", " + bw,
ketamine=kw, power="log10" if log else "linear", largest=med.index[0],
largest_value=med.iloc[0], ORB=med["ORB"], ORB_rank=v["rank"]["ORB"],
IFG_rank=v["rank"]["IFG"], MFG_rank=v["rank"]["MFG"], median_region=med.median(),
n_regions=v["n_regions"], ORB_first=v["orb_first"], frontal_above_median=v["frontal_above"],
met=v["orb_first"] and v["frontal_above"]))
V15 = pd.DataFrame(rows)
print(V15.to_string(index=False, float_format=lambda x: f"{x:.3f}"))
print("variants meeting the criterion:", int(V15.met.sum()), "of", len(V15), "; ORB first:", int(V15.ORB_first.sum()),
"; ORB, IFG and MFG above the median region:", int(V15.frontal_above_median.sum()))
print("largest region across the variants:", dict(Counter(V15.largest).most_common()))
mag, v = F15[("pre", "before", "after", True)]
print("\nmain reading (baseline: code preinfusion session before the first puff; ketamine: first puff to end; log10):")
print(pd.DataFrame({"median": v["median"], "participants": v["n_pts"].reindex(v["median"].index),
"rank": pd.Series(v["rank"])}).T.round(3).to_string())
print("\nmagnitude by participant (rows) and region:\n", mag[v["median"].index].round(2).to_string())
mean15 = mag.mean().sort_values(ascending=False)
print("\nmean +/- s.e.m. over participants (the fig. S15H legend's statistic; the criterion uses the median):\n",
pd.DataFrame({"mean": mean15, "sem": mag.sem()[mean15.index], "participants": mag.notna().sum()[mean15.index]}).T.round(3).to_string())
print("under the mean, ORB ranks", int(mean15.index.get_loc("ORB")) + 1, "of", len(mean15), "; largest", mean15.index[0])
prank = mag.rank(axis=1, ascending=False)
print("\nrank of ORB, IFG, MFG within each participant's regions (NaN: no channel), and the participant's number of regions:")
print(pd.DataFrame({r: prank[r] for r in ["ORB", "IFG", "MFG"]}).assign(n_regions=mag.notna().sum(1)).to_string())
fig, ax = plt.subplots(figsize=(8, 2.4))
ax.bar(range(len(v["median"])), v["median"].values, color=["C1" if r in ("ORB", "IFG", "MFG") else "C0" for r in v["median"].index])
ax.set_xticks(range(len(v["median"]))); ax.set_xticklabels(v["median"].index, rotation=90, fontsize=7); ax.set_ylabel("median |z|")
plt.tight_layout(); plt.show()
print("plotted range:", round(float(v["median"].min()), 3), "to", round(float(v["median"].max()), 3))
baseline ketamine power largest largest_value ORB ORB_rank IFG_rank MFG_rank median_region n_regions ORB_first frontal_above_median met
code session, before after log10 MCC 2.644 2.141 3 2 4 1.386 21 False True False
code session, before after linear MCC 3.398 2.500 3 2 4 1.541 21 False True False
code session, before after_no_puffs log10 MCC 2.847 2.237 3 2 4 1.593 21 False True False
code session, before after_no_puffs linear MCC 3.745 2.609 3 2 4 1.787 21 False True False
code session, before before log10 MCC 3.014 2.280 2 3 4 1.408 21 False True False
code session, before before linear MCC 4.318 2.692 2 3 4 1.580 21 False True False
code session, before whole log10 MCC 2.845 2.332 2 3 4 1.375 21 False True False
code session, before whole linear MCC 3.898 2.773 2 3 4 1.587 21 False True False
code session, whole after log10 MCC 1.854 1.687 3 4 2 1.297 21 False True False
code session, whole after linear MCC 2.077 1.854 3 5 2 1.343 21 False True False
code session, whole after_no_puffs log10 MCC 1.996 1.992 2 4 3 1.487 21 False True False
code session, whole after_no_puffs linear ORB 2.404 2.404 1 4 3 1.563 21 True True True
code session, whole before log10 MCC 1.984 1.669 4 3 2 1.276 21 False True False
code session, whole before linear MCC 2.252 1.810 4 3 2 1.354 21 False True False
code session, whole whole log10 MCC 1.927 1.671 3 4 2 1.285 21 False True False
code session, whole whole linear MCC 2.178 1.832 3 4 2 1.425 21 False True False
SD245 Preinfusion, before after log10 MCC 2.528 2.141 3 2 4 1.237 21 False True False
SD245 Preinfusion, before after_no_puffs log10 MCC 2.749 2.237 3 2 4 1.567 21 False True False
SD245 Preinfusion, before before log10 MCC 2.910 2.280 2 3 4 1.282 21 False True False
SD245 Preinfusion, before whole log10 MCC 2.734 2.332 2 3 4 1.278 21 False True False
SD245 Preinfusion, whole after log10 MFG 1.810 1.687 3 4 1 1.222 21 False True False
SD245 Preinfusion, whole after_no_puffs log10 ORB 1.992 1.992 1 4 2 1.442 21 True True True
SD245 Preinfusion, whole before log10 MCC 1.842 1.669 4 3 2 1.201 21 False True False
SD245 Preinfusion, whole whole log10 MFG 1.777 1.671 3 4 1 1.179 21 False True False
variants meeting the criterion: 2 of 24 ; ORB first: 2 ; ORB, IFG and MFG above the median region: 24
largest region across the variants: {'MCC': 20, 'ORB': 2, 'MFG': 2}
main reading (baseline: code preinfusion session before the first puff; ketamine: first puff to end; log10):
MCC IFG ORB MFG PMC ACC SMG MOT STS STG ITG INS PHG SS THAL ANT MTG AMY HIPP BG TP THAL POS
median 2.644 2.267 2.141 2.048 1.721 1.634 1.583 1.536 1.509 1.391 1.386 1.348 1.216 1.197 1.193 1.17 1.165 1.141 0.972 0.952 0.831
participants 2.000 6.000 5.000 2.000 6.000 2.000 5.000 2.000 3.000 4.000 4.000 6.000 3.000 2.000 5.000 7.00 4.000 7.000 5.000 3.000 6.000
rank 1.000 2.000 3.000 4.000 5.000 6.000 7.000 8.000 9.000 10.000 11.000 12.000 13.000 14.000 15.000 16.00 17.000 18.000 19.000 20.000 21.000
magnitude by participant (rows) and region:
MCC IFG ORB MFG PMC ACC SMG MOT STS STG ITG INS PHG SS THAL ANT MTG AMY HIPP BG TP THAL POS
SD168 NaN 2.15 NaN NaN 0.94 NaN 1.06 NaN 1.51 1.02 NaN 1.13 1.22 NaN 1.61 1.32 1.14 1.14 NaN NaN 0.59
SD172 NaN 1.79 2.08 NaN 1.53 NaN 1.58 NaN 0.64 1.76 1.28 1.30 0.85 NaN 0.91 1.10 1.19 1.19 0.98 0.73 0.71
SD181 NaN 2.39 2.39 2.25 1.91 1.46 NaN 1.80 NaN NaN 1.55 2.02 NaN 1.42 1.11 1.63 1.72 2.38 0.97 1.19 NaN
SD185 3.02 3.06 2.14 1.85 1.93 1.81 NaN 1.27 NaN NaN NaN NaN NaN 0.97 NaN 1.24 NaN 1.27 1.13 NaN 1.20
SD186 NaN 1.12 1.35 NaN 1.14 NaN 0.88 NaN NaN 0.78 1.27 0.88 NaN NaN NaN 0.64 0.42 0.77 0.48 NaN 0.56
SD211 NaN NaN NaN NaN NaN NaN 1.86 NaN NaN NaN 1.49 1.40 1.39 NaN 1.19 1.17 NaN 0.97 0.87 0.95 0.95
SD245 2.27 3.69 3.34 NaN 2.38 NaN 2.02 NaN 2.36 2.15 NaN 1.93 NaN NaN 4.18 1.10 NaN 0.71 NaN NaN 1.38
mean +/- s.e.m. over participants (the fig. S15H legend's statistic; the criterion uses the median):
MCC IFG ORB MFG THAL ANT PMC ACC MOT STS SMG INS STG ITG HIPP SS MTG PHG AMY TP THAL POS BG
mean 2.644 2.367 2.262 2.048 1.801 1.638 1.634 1.536 1.503 1.480 1.443 1.427 1.397 1.205 1.197 1.173 1.152 1.118 0.957 0.898 0.886
sem 0.376 0.372 0.321 0.202 0.607 0.221 0.173 0.262 0.496 0.222 0.184 0.320 0.073 0.212 0.227 0.112 0.159 0.267 0.135 0.137 0.108
participants 2.000 6.000 5.000 2.000 5.000 6.000 2.000 2.000 3.000 5.000 6.000 4.000 4.000 7.000 2.000 7.000 3.000 4.000 3.000 6.000 5.000
under the mean, ORB ranks 3 of 21 ; largest MCC
rank of ORB, IFG, MFG within each participant's regions (NaN: no channel), and the participant's number of regions:
ORB IFG MFG n_regions
SD168 NaN 1.0 NaN 12
SD172 1.0 2.0 NaN 16
SD181 1.0 2.0 4.0 15
SD185 3.0 1.0 5.0 12
SD186 1.0 4.0 NaN 12
SD211 NaN NaN NaN 10
SD245 3.0 2.0 NaN 12
plotted range: 0.831 to 2.644
Result: not reproduced, on approximate windows. In the main reading (baseline the preinfusion file before the first puff, ketamine the second half of the infusion file, about 30 to 40 min, log power), the median magnitude is largest in MCC (2.644, 2 participants), then IFG (2.267, 6 participants), ORB (2.141, 5) and MFG (2.048, 2), against a median region of 1.386 over 21 regions. The three frontal regions rank 2 to 4, above the median region, but ORB is third, so the criterion fails on its first part. Over the 24 window, baseline and power-scale variants, ORB, IFG and MFG are above the median region in all 24, and ORB is first in 2 (both with the whole preinfusion file as baseline and the infusion file after the last puff, 16 windows); MCC is first in 20 and MFG in 2. Within participants, ORB is the largest region in SD172, SD181 and SD186 and third in SD185 and SD245, the other two participants with ORB channels. MCC rests on SD185 (1 channel) and SD245 (14 channels). The criterion takes the median over participants; the fig. S15H legend gives the magnitude as mean +/- s.e.m. Under the mean MCC is still first (2.644), then IFG (2.367) and ORB (2.262), so ORB is again third. Cohort (b) cannot be separated from (a): the 7 participants with LFP files are the included ones, and the dataset does not record the infusion type. The test is approximate: the paper's 25-min baseline does not fit in the 20-min files, and the 30 to 40 min window is placed from the Methods' "approximately 30 minutes" and the 10 min of recording before the first puff.
figs27-slow-oscillations-with-strong-dissociation¶
All participants with strong dissociation (CADSS above 10) show ketamine-induced oscillations below 10 Hz, and high-CADSS participants tend to have a larger fraction of channels with such oscillations.
The detection of the Methods (Band-limited oscillation detection; criterion), on the channels and
spectrograms of figs15: segment spectra are the mean of the 30-s multitaper windows inside each 3-min
segment; FOOOF removes the aperiodic part; scipy.signal.find_peaks finds the peaks; peaks are clustered
by 80% overlap of peak width; clusters in fewer than 40% of a window's segments are dropped; a ketamine
cluster below 10 Hz counts if no baseline cluster lies within 2 Hz (new frequency) or if a one-sided
rank-sum test of its peak heights against those of a baseline cluster within 2 Hz gives p < 0.001.
Details the paper does not give, each an assumption: FOOOF 'fixed' aperiodic mode with default settings,
fitted on 1 to 40 Hz; peaks above the aperiodic fit with a prominence of at least 0.1 log10 units (0.05
and 0.2 as variants); a peak's width at half its prominence; clusters by single linkage, two peaks
joined when their overlap is at least 80% of the narrower one's width ('min'; 80% of the wider one's,
'max', as a variant); one height per segment (the highest); a cluster's frequency the median of its peak
centres, which also decides "within 2 Hz" and "below 10 Hz".
Windows: baseline the preinfusion and postinfusion files combined (Methods: "pre- and post-infusion windows were combined"; the whole 20-min files, 12 segments), with the preinfusion file alone (6 segments) and its first half (3) as variants; ketamine window the second half of the infusion file (30 to 40 min, the criterion's fallback, 3 segments), with the first half (3) and the whole file (6) as variants. No sleep or ictal annotations exist, so none are excluded. As a reference without ketamine, the same detection tests the second half of the postinfusion file against the preinfusion file.
CADSS: the infusion score per participant (maximum over the infusion time points, as the code's
utils.load_full_cadss and notebook 01) from the code record's data/cadss_across_pts.csv. High CADSS
under both readings of Open H, > 10 (Results) and >= 10 (Methods, fig. S27F). Cohorts as in figs15: only
(a) can be run. The label permutation of the Methods (not required) is enumerated over every assignment of
the high label.
CAD = he.load_cadss().loc[cd.PARTICIPANTS, ["pre", "infusion", "post"]]
print("CADSS (code record), the 7 participants:\n", CAD.to_string())
HIGH = {">10": [p for p in cd.PARTICIPANTS if CAD.loc[p, "infusion"] > 10],
">=10": [p for p in cd.PARTICIPANTS if CAD.loc[p, "infusion"] >= 10]}
print("high CADSS:", HIGH)
t = time.time()
FL = ie.s27_flats(REST, {"pre": ["whole"], "post": ["whole", "after"], "inf": ["after", "before", "whole"]})
print(f"FOOOF fits: {sum(v[1].shape[0] for v in FL.values())} channel segments in {len(FL)} segment spectra, {time.time() - t:.0f} s; "
f"non-finite values in the flattened spectra: {sum(int((~np.isfinite(v[1])).sum()) for v in FL.values())}")
import scipy.stats as sst
print("lowest one-sided rank-sum p reachable (no ties), ketamine against baseline segments:")
for nk, nb_ in [(3, 12), (3, 6), (3, 3), (6, 12), (6, 6), (6, 3)]:
print(f" {nk} against {nb_}: {sst.mannwhitneyu(np.arange(nb_, nb_ + nk), np.arange(nb_), alternative='greater').pvalue:.4g}")
CADSS (code record), the 7 participants:
pre infusion post
pt
SD168 0 9.0 3.0
SD172 0 30.0 0.0
SD181 0 10.0 0.0
SD185 0 0.0 0.0
SD186 1 33.0 2.0
SD211 3 26.0 0.0
SD245 2 39.0 1.0
high CADSS: {'>10': ['SD172', 'SD186', 'SD211', 'SD245'], '>=10': ['SD172', 'SD181', 'SD186', 'SD211', 'SD245']}
FOOOF fits: 10992 channel segments in 168 segment spectra, 35 s; non-finite values in the flattened spectra: 0 lowest one-sided rank-sum p reachable (no ties), ketamine against baseline segments: 3 against 12: 0.002198 3 against 6: 0.0119 3 against 3: 0.05 6 against 12: 5.387e-05 6 against 6: 0.001082 6 against 3: 0.0119
BASES = {"pre + post": [("pre", "whole"), ("post", "whole")], "pre": [("pre", "whole")], "pre first half": [("pre", "before")]}
def summarise(r, hi):
n_osc = {p: r[p]["n_osc"] for p in r}; n_ch = {p: r[p]["n_ch"] for p in r}
d, pp, nn = ie.label_permutation(n_osc, n_ch, hi)
return sum(n_osc[p] > 0 for p in hi), d, pp, nn
t = time.time()
rows, S27 = [], {}
for bname, base in BASES.items():
for kw in ["after", "before", "whole"]:
for prom in [0.05, 0.1, 0.2]:
for ov in ["min", "max"]:
r = ie.s27_channels(REST, FL, base, kw, prom, ov)
S27[(bname, kw, prom, ov)] = r
p0 = cd.PARTICIPANTS[0]
row = dict(baseline=bname, ketamine=kw, prominence=prom, overlap=ov,
segments=f"{r[p0]['n_base_seg']}/{r[p0]['n_ket_seg']}")
for thr, hi in HIGH.items():
k, d, pp, nn = summarise(r, hi)
row[f"high{thr} with osc"] = f"{k}/{len(hi)}"
row[f"diff{thr}"] = d; row[f"perm_p{thr}"] = pp
row.update({p: r[p]["n_osc"] for p in cd.PARTICIPANTS})
rows.append(row)
V27 = pd.DataFrame(rows)
print(f"({time.time() - t:.0f} s) channels with a ketamine-unique oscillation below 10 Hz, by participant, for each variant;"
" diff = fraction of channels in high minus low CADSS participants, perm_p its exact permutation p:")
print(V27.to_string(index=False, float_format=lambda x: f"{x:.3f}"))
for thr, hi in HIGH.items():
ok = V27[f"high{thr} with osc"] == f"{len(hi)}/{len(hi)}"
print(f"CADSS {thr}: variants in which every high-CADSS participant has at least one channel: {int(ok.sum())} of {len(V27)}")
print(f"(the label permutation has {summarise(S27[('pre + post', 'after', 0.1, 'min')], HIGH['>10'])[3]} assignments for >10 "
f"and {summarise(S27[('pre + post', 'after', 0.1, 'min')], HIGH['>=10'])[3]} for >=10)")
(22 s) channels with a ketamine-unique oscillation below 10 Hz, by participant, for each variant; diff = fraction of channels in high minus low CADSS participants, perm_p its exact permutation p:
baseline ketamine prominence overlap segments high>10 with osc diff>10 perm_p>10 high>=10 with osc diff>=10 perm_p>=10 SD168 SD172 SD181 SD185 SD186 SD211 SD245
pre + post after 0.050 min 12/3 4/4 -0.019 0.571 5/5 0.200 0.048 8 31 34 3 5 7 17
pre + post after 0.050 max 12/3 4/4 -0.042 0.657 5/5 0.128 0.143 4 22 25 3 6 1 6
pre + post after 0.100 min 12/3 4/4 -0.044 0.629 5/5 0.186 0.048 5 28 33 2 4 6 8
pre + post after 0.100 max 12/3 4/4 -0.055 0.714 5/5 0.125 0.286 5 24 26 2 6 1 2
pre + post after 0.200 min 12/3 4/4 -0.045 0.629 5/5 0.188 0.143 3 26 32 1 7 5 2
pre + post after 0.200 max 12/3 3/4 -0.075 0.714 4/5 0.133 0.238 3 20 28 1 4 2 0
pre + post before 0.050 min 12/3 4/4 -0.052 0.714 5/5 0.139 0.286 11 32 32 6 11 4 10
pre + post before 0.050 max 12/3 4/4 -0.054 0.657 5/5 0.106 0.286 7 31 26 6 8 2 1
pre + post before 0.100 min 12/3 4/4 -0.082 0.743 5/5 0.110 0.286 9 29 31 7 10 3 4
pre + post before 0.100 max 12/3 3/4 -0.072 0.714 4/5 0.112 0.286 5 26 28 6 10 1 0
pre + post before 0.200 min 12/3 4/4 -0.019 0.629 5/5 0.125 0.238 4 20 22 3 13 3 1
pre + post before 0.200 max 12/3 3/4 -0.029 0.657 4/5 0.101 0.238 3 16 19 2 9 2 0
pre + post whole 0.050 min 12/6 4/4 -0.058 0.657 5/5 0.200 0.143 12 38 40 3 10 6 10
pre + post whole 0.050 max 12/6 4/4 -0.059 0.686 5/5 0.152 0.238 5 27 31 4 10 2 3
pre + post whole 0.100 min 12/6 4/4 -0.103 0.714 5/5 0.148 0.238 10 28 38 5 10 4 7
pre + post whole 0.100 max 12/6 4/4 -0.079 0.714 5/5 0.139 0.238 5 22 31 3 9 2 2
pre + post whole 0.200 min 12/6 4/4 -0.041 0.629 5/5 0.163 0.190 6 26 30 2 13 4 1
pre + post whole 0.200 max 12/6 3/4 -0.055 0.686 4/5 0.147 0.190 3 20 28 2 10 3 0
pre after 0.050 min 6/3 4/4 -0.039 0.657 5/5 0.142 0.238 14 34 31 3 3 6 16
pre after 0.050 max 6/3 4/4 -0.044 0.686 5/5 0.126 0.238 12 28 28 2 8 5 8
pre after 0.100 min 6/3 4/4 -0.076 0.714 5/5 0.111 0.381 14 30 31 3 5 5 9
pre after 0.100 max 6/3 4/4 -0.051 0.686 5/5 0.119 0.286 11 30 28 3 9 4 4
pre after 0.200 min 6/3 4/4 -0.091 0.714 5/5 0.116 0.333 13 27 32 1 7 6 2
pre after 0.200 max 6/3 4/4 -0.078 0.714 5/5 0.133 0.238 7 22 30 1 7 4 1
pre before 0.050 min 6/3 4/4 -0.084 0.800 5/5 0.060 0.429 19 34 29 7 10 3 10
pre before 0.050 max 6/3 4/4 -0.067 0.743 5/5 0.101 0.333 14 32 29 4 11 3 4
pre before 0.100 min 6/3 4/4 -0.128 0.857 5/5 0.023 0.476 19 30 29 7 9 2 4
pre before 0.100 max 6/3 4/4 -0.104 0.743 5/5 0.092 0.333 10 26 31 6 11 1 2
pre before 0.200 min 6/3 4/4 -0.086 0.800 5/5 0.029 0.381 16 19 22 3 12 3 2
pre before 0.200 max 6/3 4/4 -0.062 0.771 5/5 0.065 0.333 10 16 21 2 10 3 2
pre whole 0.050 min 6/6 4/4 -0.069 0.800 5/5 0.100 0.381 19 37 31 3 5 5 11
pre whole 0.050 max 6/6 4/4 -0.056 0.714 5/5 0.120 0.286 14 31 30 3 10 5 7
pre whole 0.100 min 6/6 4/4 -0.124 0.829 5/5 0.048 0.429 19 29 31 5 7 3 7
pre whole 0.100 max 6/6 4/4 -0.082 0.714 5/5 0.116 0.286 12 26 31 2 9 4 4
pre whole 0.200 min 6/6 4/4 -0.084 0.771 5/5 0.069 0.333 17 27 27 2 10 5 2
pre whole 0.200 max 6/6 4/4 -0.092 0.743 5/5 0.097 0.286 11 20 29 2 9 5 2
pre first half after 0.050 min 3/3 4/4 -0.035 0.629 5/5 0.134 0.190 14 32 30 4 3 8 17
pre first half after 0.050 max 3/3 4/4 -0.010 0.571 5/5 0.136 0.095 15 30 27 2 11 10 10
pre first half after 0.100 min 3/3 4/4 -0.061 0.714 5/5 0.128 0.238 15 31 32 3 5 10 10
pre first half after 0.100 max 3/3 4/4 -0.029 0.600 5/5 0.142 0.190 12 31 29 3 11 9 5
pre first half after 0.200 min 3/3 4/4 -0.096 0.714 5/5 0.124 0.286 14 28 34 1 6 9 2
pre first half after 0.200 max 3/3 4/4 -0.095 0.714 5/5 0.115 0.238 10 20 31 1 6 7 2
pre first half before 0.050 min 3/3 4/4 -0.070 0.800 5/5 0.074 0.381 16 31 28 7 9 4 11
pre first half before 0.050 max 3/3 4/4 -0.063 0.714 5/5 0.104 0.238 13 29 29 5 9 8 5
pre first half before 0.100 min 3/3 4/4 -0.124 0.886 5/5 0.048 0.524 17 29 31 7 7 6 4
pre first half before 0.100 max 3/3 4/4 -0.085 0.714 5/5 0.107 0.238 10 26 31 6 9 7 3
pre first half before 0.200 min 3/3 4/4 -0.063 0.743 5/5 0.074 0.333 13 21 24 3 13 5 2
pre first half before 0.200 max 3/3 4/4 -0.039 0.657 5/5 0.094 0.143 7 15 21 2 9 6 3
pre first half whole 0.050 min 3/6 4/4 -0.064 0.800 5/5 0.103 0.381 16 33 30 4 4 6 12
pre first half whole 0.050 max 3/6 4/4 -0.066 0.714 5/5 0.088 0.333 19 28 29 3 10 9 9
pre first half whole 0.100 min 3/6 4/4 -0.111 0.771 5/5 0.066 0.381 19 30 32 5 7 7 7
pre first half whole 0.100 max 3/6 4/4 -0.085 0.714 5/5 0.111 0.286 15 25 32 2 8 10 5
pre first half whole 0.200 min 3/6 4/4 -0.074 0.714 5/5 0.084 0.333 17 28 28 2 10 8 2
pre first half whole 0.200 max 3/6 4/4 -0.084 0.743 5/5 0.103 0.238 11 19 29 2 7 9 3
CADSS >10: variants in which every high-CADSS participant has at least one channel: 50 of 54
CADSS >=10: variants in which every high-CADSS participant has at least one channel: 50 of 54
(the label permutation has 35 assignments for >10 and 21 for >=10)
MAIN27 = S27[("pre + post", "after", 0.1, "min")]
rows = []
for p in cd.PARTICIPANTS:
r = MAIN27[p]
fr = [u[0] for h in r["hits"] for u in h[2]]
kinds = Counter(u[1] for h in r["hits"] for u in h[2])
rows.append(dict(participant=p, CADSS=CAD.loc[p, "infusion"], high_gt10=p in HIGH[">10"], high_ge10=p in HIGH[">=10"],
channels=r["n_ch"], with_osc=r["n_osc"], fraction=r["n_osc"] / r["n_ch"],
new=kinds.get("new", 0), amplified=kinds.get("amplified", 0),
freq_range=f"{min(fr):.1f} to {max(fr):.1f}" if fr else "",
regions=dict(Counter(h[1] for h in r["hits"]).most_common(4))))
T27 = pd.DataFrame(rows).set_index("participant")
print("main reading (baseline pre + post, ketamine 30 to 40 min, prominence 0.1, overlap of the narrower peak):")
print(T27.to_string(float_format=lambda x: f"{x:.3f}"))
for thr, hi in HIGH.items():
k, d, pp, nn = summarise(MAIN27, hi)
lo = [p for p in cd.PARTICIPANTS if p not in hi]
fh = sum(MAIN27[p]["n_osc"] for p in hi) / sum(MAIN27[p]["n_ch"] for p in hi)
fl_ = sum(MAIN27[p]["n_osc"] for p in lo) / sum(MAIN27[p]["n_ch"] for p in lo)
print(f"CADSS {thr}: {k} of {len(hi)} high participants with >= 1 channel; channel fraction high {fh:.3f}, low {fl_:.3f}, "
f"difference {d:.3f}, permutation p {pp:.3f} ({nn} assignments)")
print("\nall channels with an oscillation, by region (main reading):", dict(Counter(h[1] for p in cd.PARTICIPANTS for h in MAIN27[p]["hits"]).most_common()))
# reference without ketamine: postinfusion second half against the preinfusion file
REF = ie.s27_channels(REST, FL, [("pre", "whole")], "after", 0.1, "min", ket_key="post")
print("\nreference, postinfusion (second half) against preinfusion, channels with a 'unique' oscillation below 10 Hz:",
{p: f"{REF[p]['n_osc']}/{REF[p]['n_ch']}" for p in cd.PARTICIPANTS})
same = S27[("pre", "after", 0.1, "min")]
print("the same baseline (preinfusion file) with the infusion second half:", {p: f"{same[p]['n_osc']}/{same[p]['n_ch']}" for p in cd.PARTICIPANTS})
above = {p: same[p]["n_osc"] > REF[p]["n_osc"] for p in cd.PARTICIPANTS}
for thr, hi in HIGH.items():
print(f"CADSS {thr}: high participants with more channels on infusion than in the reference: {sum(above[p] for p in hi)} of {len(hi)}")
print("all participants:", sum(above.values()), "of", len(above), "; total channels, infusion", sum(same[p]["n_osc"] for p in cd.PARTICIPANTS),
"against reference", sum(REF[p]["n_osc"] for p in cd.PARTICIPANTS))
main reading (baseline pre + post, ketamine 30 to 40 min, prominence 0.1, overlap of the narrower peak):
CADSS high_gt10 high_ge10 channels with_osc fraction new amplified freq_range regions
participant
SD168 9.000 False False 86 5 0.058 5 0 3.5 to 8.8 {'STS': 3, 'PMC': 1, 'STG': 1}
SD172 30.000 True True 97 28 0.289 28 0 4.2 to 9.8 {'INS': 8, 'PMC': 7, 'IFG': 6, 'ORB': 2}
SD181 10.000 False True 58 33 0.569 33 0 5.2 to 9.0 {'IFG': 10, 'ACC': 5, 'TP': 5, 'INS': 4}
SD185 0.000 False False 43 2 0.047 2 0 7.5 to 7.5 {'ORB': 1, 'HIPP': 1}
SD186 33.000 True True 47 4 0.085 4 0 5.0 to 8.0 {'MTG': 2, 'ORB': 1, 'PMC': 1}
SD211 26.000 True True 44 6 0.136 6 0 2.8 to 6.5 {'HIPP': 2, 'TP': 2, 'MTG': 1, 'ITG': 1}
SD245 39.000 True True 83 8 0.096 8 0 4.0 to 8.2 {'PMC': 2, 'MCC': 2, 'THAL POS': 1, 'ORB': 1}
CADSS >10: 4 of 4 high participants with >= 1 channel; channel fraction high 0.170, low 0.214, difference -0.044, permutation p 0.629 (35 assignments)
CADSS >=10: 5 of 5 high participants with >= 1 channel; channel fraction high 0.240, low 0.054, difference 0.186, permutation p 0.048 (21 assignments)
all channels with an oscillation, by region (main reading): {np.str_('IFG'): 16, np.str_('PMC'): 12, np.str_('INS'): 12, np.str_('ORB'): 8, np.str_('TP'): 7, np.str_('HIPP'): 5, np.str_('MTG'): 5, np.str_('ACC'): 5, np.str_('STS'): 3, np.str_('MFG'): 3, np.str_('ITG'): 2, np.str_('THAL POS'): 2, np.str_('MCC'): 2, np.str_('STG'): 1, np.str_('BG'): 1, np.str_('MOT'): 1, np.str_('SS'): 1}
reference, postinfusion (second half) against preinfusion, channels with a 'unique' oscillation below 10 Hz: {'SD168': '16/86', 'SD172': '24/97', 'SD181': '1/58', 'SD185': '1/43', 'SD186': '3/47', 'SD211': '6/44', 'SD245': '9/83'}
the same baseline (preinfusion file) with the infusion second half: {'SD168': '14/86', 'SD172': '30/97', 'SD181': '31/58', 'SD185': '3/43', 'SD186': '5/47', 'SD211': '5/44', 'SD245': '9/83'}
CADSS >10: high participants with more channels on infusion than in the reference: 2 of 4
CADSS >=10: high participants with more channels on infusion than in the reference: 3 of 5
all participants: 4 of 7 ; total channels, infusion 97 against reference 60
Result: reproduced as the criterion is written, on approximate windows; the detection is not specific to ketamine. In the main reading (baseline the preinfusion and postinfusion files, ketamine the second half of the infusion file, prominence 0.1), every high-CADSS participant has at least one channel with a ketamine-unique oscillation below 10 Hz: SD172 28 of 97 channels, SD186 4 of 47, SD211 6 of 44 and SD245 8 of 83 (CADSS > 10, 4 of 4), and SD181 33 of 58 (CADSS >= 10, 5 of 5). This holds in 50 of the 54 variants under each threshold; the four exceptions all use the combined baseline with the wider-peak overlap rule, and in each SD245 has none. All detected oscillations are at new frequencies: with 3 ketamine segments the one-sided rank-sum test cannot reach p < 0.001 (lowest reachable p 0.0022 against 12 segments). The low-CADSS participants also have such channels (SD168 5 of 86, SD185 2 of 43), and the same detection run on the second half of the postinfusion file against the preinfusion file finds them in all 7 participants (60 channels, against 97 for the infusion file with the same baseline); only 2 of the 4 participants with CADSS
10 have more channels on infusion than in that reference. For the second part (not required), the channel fraction is higher in high-CADSS participants under >= 10 (0.240 against 0.054, exact permutation p 0.048 over 21 assignments) and not under > 10 (0.170 against 0.214, p 0.629 over 35). The windows and the peak threshold are assumptions, and sleep and ictal periods could not be excluded.
Blocked claim¶
figs12-erp-not-explained-by-eyelid-closure needs the onsets of spontaneous blinks. The behaviour files hold
continuous eye-closure traces and puff events but no spontaneous blink onsets (printed in the data check);
the paper's code lists them by hand (for_reference/analyze_airpuff_ieeg_natural_blinks.py, lines 42 to
79, 52 onsets in the baseline sessions). The claims file says to block in that case; blocked, data_missing.
Summary¶
| Claim | Outcome | Key numbers |
|---|---|---|
| fig2g-nmf-six-factors | reproduced (assumed elbow rule) | 458 channels; held-out EV 0.466 at six components; elbow at 6 by largest fall in gain, 8 by chord distance |
| link-ieeg-yeo-labels | reproduced | 458 of 458 channels found; 120 of 126 dataset network counts equal; deposited Fig. 2H and 2J matrices recomputed exactly |
| fig2h-factor-network-loading | partially reproduced | 4 of 6 rows; limbic (0.495) and frontoparietal (0.463) fall short |
| fig2i-ketamine-reduces-factors-1-2 | reproduced | factor 1 -0.125 (p 2.56e-05), factor 2 -0.115 (p 5.51e-08); factors 0, 4, 5 not significant; means within tolerance |
| fig2j-factor-change-by-network | reproduced (reading a) | all six named changes negative |
| figs13-auditory-control-no-response | reproduced | control below preinfusion for 6 of 6 factors (test-half preinfusion, ratios 0.33 to 0.70; all preinfusion trials, 0.34 to 0.98) |
| fig3b-cluster-test-yeo-networks | reproduced | limbic, ventral attention, somatomotor significant; maps equal the deposited ones |
| fig3c-cluster-test-regions | reproduced | 10 of 10 listed significant, 10 of 10 others not |
| figs15-ketamine-spectral-changes-frontal | not reproduced (approximate windows) | MCC largest (2.644), IFG 2.267, ORB third 2.141, MFG 2.048, median region 1.386; ORB first in 2 of 24 variants |
| figs12-erp-not-explained-by-eyelid-closure | blocked, data_missing | spontaneous blink onsets only in the code |
| figs27-slow-oscillations-with-strong-dissociation | reproduced as written (approximate windows) | 4 of 4 (CADSS > 10) and 5 of 5 (>= 10) high-CADSS participants with a channel; 50 of 54 variants; the postinfusion reference finds such channels in all 7 participants |
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.
Differences from the paper's method: every test starts from the dataset's precomputed spectrograms, not from the LFP, so the multitaper, z-scoring, bipolar referencing and channel selection are taken as deposited (the train and test means were shown to equal the means of the per-trial file). The Yeo7 label of each channel comes from the code record. The elbow rule for fig2g is an assumption. For fig3c only reading (a) was run. figs15 and figs27 are approximate tests on the 20-min LFP files: the 25-min pre-ketamine baseline is replaced by the preinfusion eyepuff file, the peak ketamine window by the second half of the infusion file (placed at about 30 to 40 min from the Methods' eyepuff timing), no sleep or ictal periods are excluded, the fig. S15 spectrogram uses the fig. S27 multitaper parameters, and the FOOOF range, peak threshold, peak width and clustering rule of figs27 are assumptions; CADSS scores come from the code record.
What was found about the dataset: the aggregated products are complete and consistent with one another and with the paper's deposited outputs, which this notebook regenerates exactly (NMF factors and loadings, network and region t maps and contours). The per-channel Yeo7 labels are not in the dataset. The per-participant LFP files are 20-min windows around each puff session (10 min before and after the first puff), with no infusion start time and no sleep or seizure annotations; the behaviour files carry no spontaneous blink onsets or CADSS scores. The network file's time bins sit 0.6 ms off the per-trial file's.
print(f"notebook ran in {(time.time() - T0) / 60:.1f} min")
notebook ran in 5.0 min