Reproducing Fedele, Tzovara et al. (2020) from dandiset 000576¶
Paper. Fedele T, Tzovara A, Steiger B, Hilfiker P, Grunwald T, Stieglitz L, Jokeit H, Sarnthein J (2020). The relation between neuronal firing, local field potentials and hemodynamic activity in the human amygdala in response to aversive dynamic visual stimuli. NeuroImage 213, 116705. doi:10.1016/j.neuroimage.2020.116705. The claims file was written from the version of record (main text, legends, Tables 1 and 2, Methods); the supplement was not supplied and carries no finding. The data descriptor of the same recordings (Fedele et al. 2021, Sci Data, doi:10.1038/s41597-020-00790-x) was read for the layout of the data.
Dandiset. DANDI:000576, version 0.231010.1811. The paper (2020) predates the dandiset (2023) and cites no version of it. This is the only published version, so it is pinned; the draft lists the same 27 assets with the same total size.
What the dandiset holds. For each of the 9 patients one session of the video task: the bipolar
iEEG of each implanted amygdala (14 channels in all) at 2 kHz as 17 epochs of 26 s, the trial table,
and for 5 patients the spike times of 42 sorted units. Each session is deposited twice, as
_ecephys.nwb and _ecephys+image.nwb (the second adds a link to the video); the data in each pair
are identical (checked below), and each is counted once. The BOLD time courses (clinical fMRI) and the
continuous microwire signals (needed for the microwire high gamma power, muHGP) are not in it.
| Claim | Tested here |
|---|---|
| fig1b-delta-decrease, fig1c-hgp-increase-first-2s, other-bands-no-effect | tested |
| table2-mi-values, fig1d-delta-hgp-mi-correlation, nsoz-delta-hgp-mi-inverse | tested |
| table2-delta-reduced-12-of-14, table2-hgp-mi-positive-nsoz, fig1e-hgp-not-delta-distinguishes-soz | tested |
| hgp-soz-sensitivity-specificity, fig1f-hgp-more-specific-than-delta | tested |
| neuron-count, fig3a-firing-increase-aversive, fig3a-firing-gradual-sustained | tested |
| fig3b-firing-hgp-no-condition-difference, fig3c-ssi-frv-correlation, fig4b-regression-firing-hgp | tested |
| fig3d-firing-correlations-modalities | tested in part (HGP only; muHGP and BOLD are not in the dandiset) |
| fig3b-firing-muhgp-higher-aversive | blocked: the continuous microwire signals are not in the dandiset |
| fig2-participant4-bold-correlations, table2-bold-correlations-group, table2-bold-correlation-values, bold-hgp-uncorrelated-low-mi, fig4c-regression-firing-bold, fig4d-hgp-bold-covariation | blocked: the BOLD time courses are not public |
Method details the claims file leaves open are taken from the paper and cited where they are used,
here and in the docstring of the helper module amygdala.py.
import time
import warnings
from concurrent.futures import ThreadPoolExecutor
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import yaml
from scipy import stats
import amygdala as A
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore")
started = time.time()
all_assets = A.assets()
ecephys = sorted([a for a in all_assets if a["path"].endswith("_ecephys.nwb")], key=lambda a: a["path"])
with ThreadPoolExecutor(9) as ex:
sessions = list(ex.map(A.load_session, ecephys))
print(f"loaded {len(sessions)} sessions from dandiset {A.DANDISET_ID} version {A.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
claims_file = Path("../../claims/000576/10.1016-j.neuroimage.2020.116705.yaml")
claims = {c["id"]: c for c in yaml.safe_load(claims_file.read_text())["claims"]}
loaded 9 sessions from dandiset 000576 version 0.231010.1811 in 23 s
Data check¶
Each check of the description against the values is stated under its output.
kinds = pd.Series([("video" if a["path"].endswith(".avi") else
"ecephys+image.nwb" if a["path"].endswith("+image.nwb") else "ecephys.nwb")
for a in all_assets]).value_counts()
print("Assets in the version:", len(all_assets))
print(kinds.to_string())
def same_data(asset):
"""Compare the iEEG, spikes and trials of an `_ecephys+image.nwb` file with its `_ecephys.nwb` twin."""
twin = next(a for a in ecephys if a["path"] == asset["path"].replace("_ecephys+image.nwb", "_ecephys.nwb"))
h, g = A.open_asset(twin["asset_id"]), A.open_asset(asset["asset_id"])
st = lambda f: f["units/spike_times"][:] if "units" in f else np.array([])
return {"file": asset["path"].split("/")[-1],
"ieeg_identical": np.array_equal(h["processing/ecephys/LFP/ecephys.ieeg/data"][:],
g["acquisition/ecephys.ieeg/data"][:]),
"spikes_identical": np.array_equal(st(h), st(g)),
"trials_identical": np.array_equal(h["intervals/trials/start_time"][:], g["intervals/trials/start_time"][:])}
with ThreadPoolExecutor(9) as ex:
twins = pd.DataFrame(ex.map(same_data, [a for a in all_assets if a["path"].endswith("+image.nwb")]))
print("\nEach `_ecephys+image.nwb` file against its `_ecephys.nwb` twin:")
print(twins.to_string(index=False))
Assets in the version: 27 ecephys.nwb 9 ecephys+image.nwb 9 video 9
Each `_ecephys+image.nwb` file against its `_ecephys.nwb` twin:
file ieeg_identical spikes_identical trials_identical
sub-01_ses-20140828T132700_ecephys+image.nwb True True True
sub-03_ses-20150304T164400_ecephys+image.nwb True True True
sub-02_ses-20160302T164800_ecephys+image.nwb True True True
sub-06_ses-20131209T113600_ecephys+image.nwb True True True
sub-07_ses-20160212T171400_ecephys+image.nwb True True True
sub-08_ses-20161217T120200_ecephys+image.nwb True True True
sub-04_ses-20120829T175200_ecephys+image.nwb True True True
sub-09_ses-20170311T170500_ecephys+image.nwb True True True
sub-05_ses-20121207T173500_ecephys+image.nwb True True True
The nine _ecephys+image.nwb files repeat the iEEG, spikes and trials of the nine _ecephys.nwb
files exactly; each recording is counted once, from the _ecephys.nwb file.
paper_table1 = {1: (31, "m"), 2: (48, "m"), 3: (19, "f"), 4: (22, "m"), 5: (34, "m"), 6: (23, "m"),
7: (39, "m"), 8: (27, "f"), 9: (22, "m")} # Table 1 of the paper: age, gender
rows = []
for s in sessions:
for c, lab in enumerate(s["labels"]):
n = A.table2_number(s["subject"], lab)
rows.append({"subject": s["subject"], "age_file": s["age"], "age_paper": paper_table1[s["subject"]][0],
"sex_file": s["sex"], "sex_paper": paper_table1[s["subject"]][1], "channel": lab,
"table2_amygdala": n, "soz_file": s["soz"][c], "soz_table2": n in A.TABLE2_SOZ,
"location": s["locations"][c], "filtering": s["filtering"][c]})
chan = pd.DataFrame(rows)
print(chan.to_string(index=False))
print("\nChannels:", len(chan), "; SOZ flags agree with Table 2:", bool((chan.soz_file == chan.soz_table2).all()),
f"({int(chan.soz_file.sum())} SOZ, {int((~chan.soz_file).sum())} nSOZ)")
print("Ages that differ from the paper's Table 1:",
[r for r in chan.drop_duplicates("subject")[["subject", "age_file", "age_paper"]].to_dict("records") if r["age_file"] != f"P{r['age_paper']}Y"])
print("Channels whose location names the hippocampus or is unknown:",
chan[~chan.location.str.contains("Amygdala")][["subject", "channel", "location"]].to_dict("records"))
subject age_file age_paper sex_file sex_paper channel table2_amygdala soz_file soz_table2 location filtering
1 P31Y 31 M m AL2-AL1 1 False False Amyg, Left Amygdala lAmyg, lateral amygdala Passband, 0.5 to 1000 Hz
1 P31Y 31 M m AR2-AR1 2 True True Hipp, Right Hippocampus rHipp, rostral hippocampus Passband, 0.5 to 1000 Hz
2 P39Y 48 M m AL2-AL1 3 False False Amyg, Left Amygdala mAmyg, medial amygdala Passband, 0.5 to 1000 Hz
2 P39Y 48 M m AR2-AR1 4 True True Amyg, Right Amygdala mAmyg, medial amygdala Passband, 0.5 to 1000 Hz
3 P19Y 19 F f AL2-AL1 5 False False Amyg, Left Amygdala mAmyg, medial amygdala Passband, 0.5 to 1000 Hz
4 P22Y 22 M m AL2-AL1 6 False False Amyg, Left Amygdala mAmyg, medial amygdala Passband, 0.5 to 1000 Hz
5 P34Y 34 M m AR2-AR1 7 True True Hipp, Right Hippocampus rHipp, rostral hippocampus Passband, 0.5 to 1000 Hz
6 P23Y 23 M m AL2-AL1 8 False False unknown Passband, 0.5 to 1000 Hz
6 P23Y 23 M m AR2-AR1 9 True True unknown Passband, 0.5 to 1000 Hz
7 P39Y 39 M m AL2-AL1 10 False False Amyg, Left Amygdala lAmyg, lateral amygdala Passband, 0.5 to 1000 Hz
7 P39Y 39 M m AR2-AR1 11 True True Amyg, Right Amygdala mAmyg, medial amygdala Passband, 0.5 to 1000 Hz
8 P24Y 27 F f AL2-AL1 12 False False Amyg, Left Amygdala mAmyg, medial amygdala Passband, 0.5 to 1000 Hz
9 P22Y 22 M m AL2-AL1 13 True True Amyg, Left Amygdala mAmyg, medial amygdala Passband, 0.5 to 1000 Hz
9 P22Y 22 M m AR2-AR1 14 True True Hipp, Right Hippocampus rHipp, rostral hippocampus Passband, 0.5 to 1000 Hz
Channels: 14 ; SOZ flags agree with Table 2: True (7 SOZ, 7 nSOZ)
Ages that differ from the paper's Table 1: [{'subject': 2, 'age_file': 'P39Y', 'age_paper': 48}, {'subject': 8, 'age_file': 'P24Y', 'age_paper': 27}]
Channels whose location names the hippocampus or is unknown: [{'subject': 1, 'channel': 'AR2-AR1', 'location': 'Hipp, Right Hippocampus rHipp, rostral hippocampus'}, {'subject': 5, 'channel': 'AR2-AR1', 'location': 'Hipp, Right Hippocampus rHipp, rostral hippocampus'}, {'subject': 6, 'channel': 'AL2-AL1', 'location': 'unknown'}, {'subject': 6, 'channel': 'AR2-AR1', 'location': 'unknown'}, {'subject': 9, 'channel': 'AR2-AR1', 'location': 'Hipp, Right Hippocampus rHipp, rostral hippocampus'}]
The 14 bipolar channels are the 14 amygdalae of Table 2 (7 nSOZ, 7 SOZ), and the files' SOZ flags
agree with the stars of Table 2. Two ages differ from the paper's Table 1 and from the data
descriptor's (subject 2: 39 in the file, 48 in the paper; subject 8: 24 against 27). The location
of five channels does not name the amygdala: the right channels of subjects 1, 5 and 9 are labelled
"Right Hippocampus, rostral hippocampus", and both of subject 6 "unknown", while the labels (AR,
AL), the paper and the data descriptor place every one in the amygdala. The analysis follows the
paper and uses all 14.
rows = []
for s in sessions:
x = s["lfp"]
d = np.abs(np.diff(x, axis=0))
bounds = np.arange(1, s["n_epochs"]) * int(A.EPOCH_S * A.FS) - 1
rows.append({"subject": s["subject"], "samples": x.shape[0], "channels": x.shape[1],
"stated_unit": s["lfp_unit"], "t_first": s["t"][0], "t_last": s["t"][-1],
"dt_unique_ms": np.unique(np.round(np.diff(s["t"]) * 1e3, 6)).tolist(),
"sd_value": round(float(x.std()), 1), "max_abs_value": round(float(np.abs(x).max()), 1),
"median_step": round(float(np.median(d)), 2), "p99_step": round(float(np.percentile(d, 99)), 2),
"median_step_at_epoch_bounds": round(float(np.median(d[bounds])), 2)})
print(pd.DataFrame(rows).to_string(index=False))
subject samples channels stated_unit t_first t_last dt_unique_ms sd_value max_abs_value median_step p99_step median_step_at_epoch_bounds
1 884000 2 volts 0.0 441.9995 [0.5, 0.500001] 17.0 164.8 1.39 5.41 13.06
2 884000 2 volts 0.0 441.9995 [0.5, 0.500001] 15.0 109.4 1.49 5.83 12.15
3 884000 1 volts 0.0 441.9995 [0.5, 0.500001] 26.5 138.7 1.73 6.63 23.83
4 884000 1 volts 0.0 441.9995 [0.5, 0.500001] 25.3 1546.0 1.45 5.58 14.77
5 884000 1 volts 0.0 441.9995 [0.5, 0.500001] 101.0 770.0 1.92 14.28 36.30
6 884000 2 volts 0.0 441.9995 [0.5, 0.500001] 44.7 567.0 1.53 6.32 29.98
7 884000 2 volts 0.0 441.9995 [0.5, 0.500001] 18.7 163.6 1.41 5.41 12.05
8 884000 1 volts 0.0 441.9995 [0.5, 0.500001] 41.3 565.4 1.77 9.40 25.59
9 884000 2 volts 0.0 441.9995 [0.5, 0.500001] 27.1 542.7 1.55 6.49 19.86
Every file holds 884000 samples per channel at 0.5 ms (2 kHz), 17 epochs of 26 s. The data are stated in volts, but the values have a standard deviation of tens and reach hundreds to 1546: they are microvolts (the data descriptor's figures are in microvolts). The unit does not enter any result (powers are z-scored, correlations are ranks). The step between consecutive samples at the 16 epoch boundaries is a median 12 to 36 against 1.4 to 1.9 within epochs, above the 99th percentile: the epochs are cut out of the recording and laid end to end, not one continuous trace, although the timestamps run without a gap. Each epoch is therefore filtered on its own (helper).
rows = []
for s in sessions:
starts, stops = s["trial_start"], s["trial_stop"]
rows.append({"subject": s["subject"], "trials": s["n_epochs"], "aversive": int(s["aversive"].sum()),
"neutral": int((~s["aversive"]).sum()), "order": "".join("A" if a else "N" for a in s["aversive"]),
"trial1": f"{starts[0]:.1f}-{stops[0]:.1f}", "trial2": f"{starts[1]:.1f}-{stops[1]:.1f}",
"start_minus_26k": sorted(set(np.round(starts - np.arange(len(starts)) * A.EPOCH_S, 3))),
"durations": sorted(set(np.round(stops - starts, 3))),
"ts_index_start_2": s["trial_ts"][1][0]})
print(pd.DataFrame(rows).to_string(index=False))
subject trials aversive neutral order trial1 trial2 start_minus_26k durations ts_index_start_2
1 17 8 9 NANANANANANANANAN 0.0-26.0 28.0-52.0 [0.0, 2.0] [24.0, 26.0] 56001
2 17 8 9 NANANANANANANANAN 0.0-26.0 28.0-52.0 [0.0, 2.0] [24.0, 26.0] 56001
3 17 8 9 NANANANANANANANAN 0.0-26.0 28.0-52.0 [0.0, 2.0] [24.0, 26.0] 56001
4 17 8 9 NANANANANANANANAN 0.0-26.0 28.0-52.0 [0.0, 2.0] [24.0, 26.0] 56001
5 17 8 9 NANANANANANANANAN 0.0-26.0 28.0-52.0 [0.0, 2.0] [24.0, 26.0] 56001
6 17 8 9 NANANANANANANANAN 0.0-26.0 28.0-52.0 [0.0, 2.0] [24.0, 26.0] 56001
7 17 8 9 NANANANANANANANAN 0.0-26.0 28.0-52.0 [0.0, 2.0] [24.0, 26.0] 56001
8 17 8 9 NANANANANANANANAN 0.0-26.0 28.0-52.0 [0.0, 2.0] [24.0, 26.0] 56001
9 17 8 9 NANANANANANANANAN 0.0-26.0 28.0-52.0 [0.0, 2.0] [24.0, 26.0] 56001
Every session has 17 trials alternating NEUTRAL and AVERSIVE, 9 and 8, as Methods 4.4 give (eight AVERSIVE and nine NEUTRAL blocks). Trials 2 to 17 start 2 s into their 26 s epoch and last 24 s: the trial table marks the block and leaves out the 2 s baseline clip before it (Methods 4.4, data descriptor Fig. 1a). Trial 1 starts at 0 and lasts 26 s, including its baseline; this is taken as an inconsistency of the table, and block onset is taken at 2 s in every epoch. That the block starts there is consistent with the high gamma and firing responses shown below, which change after 2 s and not before; it is not proven by them.
U = A.units(sessions)
urows = []
for u in U:
sp = u["spikes"]
counts = np.histogram(sp, np.arange(18) * A.EPOCH_S)[0]
urows.append({"subject": u["subject"], "unit": u["unit"], "channel": u["label"], "amygdala": u["amygdala"],
"soz": u["amygdala"] in A.TABLE2_SOZ, "spikes": len(sp),
"rate_442s": round(len(sp) / 442, 2), "rate_observed": round(len(sp) / (u["obs"].sum() * 26), 2),
"epochs_without_spikes": np.where(counts == 0)[0].tolist(),
"sorted": bool(np.all(np.diff(sp) >= 0)), "first": round(sp.min(), 2), "last": round(sp.max(), 2)})
units_df = pd.DataFrame(urows)
print(units_df.to_string(index=False))
print("\nUnits per subject:", units_df.groupby("subject").size().to_dict())
print("Units per amygdala (Table 2 number):", units_df.groupby("amygdala").size().to_dict())
print("obs_intervals of every unit:", sorted({(len(s_u["obs"]), float(s_u["obs"][0, 0]), float(s_u["obs"][-1, 1]))
for s in sessions for s_u in s["units"]}))
subject unit channel amygdala soz spikes rate_442s rate_observed epochs_without_spikes sorted first last
1 0 AL2-AL1 1 False 683 1.55 1.55 [] True 0.44 441.76
1 1 AL2-AL1 1 False 505 1.14 1.14 [] True 2.10 441.07
1 2 AL2-AL1 1 False 672 1.52 1.52 [] True 1.26 441.90
1 3 AR2-AR1 2 True 932 2.11 2.11 [] True 1.41 441.84
2 0 AL2-AL1 3 False 5023 11.36 11.36 [] True 0.11 441.63
2 1 AL2-AL1 3 False 2317 5.24 5.24 [] True 0.14 441.96
2 2 AR2-AR1 4 True 989 2.24 2.38 [0] True 33.41 440.29
2 3 AR2-AR1 4 True 899 2.03 2.16 [0] True 33.04 441.69
2 4 AR2-AR1 4 True 6445 14.58 14.58 [] True 0.05 441.98
2 5 AR2-AR1 4 True 1575 3.56 3.56 [] True 0.13 441.97
2 6 AR2-AR1 4 True 4487 10.15 10.15 [] True 0.00 441.96
2 7 AR2-AR1 4 True 1317 2.98 2.98 [] True 0.03 441.97
2 8 AR2-AR1 4 True 2364 5.35 5.68 [0] True 32.92 441.82
2 9 AR2-AR1 4 True 1191 2.69 2.86 [0] True 33.55 441.81
2 10 AR2-AR1 4 True 2260 5.11 5.11 [] True 0.35 441.96
2 11 AL2-AL1 3 False 481 1.09 1.09 [] True 0.91 441.81
3 0 AL2-AL1 5 False 1206 2.73 2.73 [] True 0.66 441.98
3 1 AL2-AL1 5 False 787 1.78 1.78 [] True 0.14 441.48
3 2 AL2-AL1 5 False 1115 2.52 2.52 [] True 0.03 441.74
4 0 AL2-AL1 6 False 7219 16.33 16.33 [] True 0.01 441.98
4 1 AL2-AL1 6 False 1903 4.31 4.31 [] True 0.60 441.85
4 2 AL2-AL1 6 False 849 1.92 1.92 [] True 0.15 441.97
4 3 AL2-AL1 6 False 3180 7.19 7.19 [] True 0.13 441.99
4 4 AL2-AL1 6 False 687 1.55 1.55 [] True 1.24 441.54
4 5 AL2-AL1 6 False 2374 5.37 5.37 [] True 0.40 441.75
4 6 AL2-AL1 6 False 3667 8.30 8.30 [] True 0.09 441.93
7 0 AL2-AL1 10 False 937 2.12 3.28 [11, 12, 13, 14, 15, 16] True 0.74 285.90
7 1 AL2-AL1 10 False 390 0.88 1.36 [11, 12, 13, 14, 15, 16] True 1.62 280.54
7 2 AL2-AL1 10 False 1185 2.68 4.14 [11, 12, 13, 14, 15, 16] True 0.44 285.85
7 3 AL2-AL1 10 False 466 1.05 1.63 [11, 12, 13, 14, 15, 16] True 1.85 285.48
7 4 AL2-AL1 10 False 487 1.10 1.70 [11, 12, 13, 14, 15, 16] True 0.49 285.13
7 5 AL2-AL1 10 False 390 0.88 1.36 [11, 12, 13, 14, 15, 16] True 0.50 285.81
7 6 AL2-AL1 10 False 1189 2.69 4.16 [11, 12, 13, 14, 15, 16] True 0.07 285.89
7 7 AL2-AL1 10 False 1412 3.19 4.94 [11, 12, 13, 14, 15, 16] True 0.00 285.61
7 8 AL2-AL1 10 False 2517 5.69 8.80 [11, 12, 13, 14, 15, 16] True 0.44 285.99
7 9 AR2-AR1 11 True 641 1.45 2.24 [11, 12, 13, 14, 15, 16] True 0.05 285.93
7 10 AR2-AR1 11 True 626 1.42 2.19 [11, 12, 13, 14, 15, 16] True 0.43 284.75
7 11 AR2-AR1 11 True 1429 3.23 5.00 [11, 12, 13, 14, 15, 16] True 0.19 285.74
7 12 AR2-AR1 11 True 827 1.87 2.89 [11, 12, 13, 14, 15, 16] True 0.32 285.09
7 13 AR2-AR1 11 True 2142 4.85 7.49 [11, 12, 13, 14, 15, 16] True 0.19 285.90
7 14 AR2-AR1 11 True 1080 2.44 3.78 [11, 12, 13, 14, 15, 16] True 0.01 285.85
7 15 AR2-AR1 11 True 4304 9.74 15.05 [11, 12, 13, 14, 15, 16] True 0.17 285.98
Units per subject: {1: 4, 2: 12, 3: 3, 4: 7, 7: 16}
Units per amygdala (Table 2 number): {1: 3, 2: 1, 3: 3, 4: 9, 5: 3, 6: 7, 10: 9, 11: 7}
obs_intervals of every unit: [(17, 0.0, 442.0)]
42 units in 8 amygdalae of 5 patients (1 L and R, 2 L and R, 3 L, 4 L, 7 L and R). Three of these
amygdalae (1 R, 2 R, 7 R, 17 units) are SOZ amygdalae in Table 2, while Discussion 3.4 says there
were no microwire recordings from SOZ amygdalae; the neuron claims use all units, as the paper's
counts (41 neurons) do. Spike times are sorted and lie within
0 to 442 s. Every unit's obs_intervals claims all 17 epochs, but the 16 units of subject 7 have no
spike after epoch 11 (the last at 286 s), and 4 units of subject 2 have none in epoch 1, at rates of
2 to 5 Hz elsewhere. These epochs are taken as not recorded for those units (helper,
observed_epochs); only the first or last epochs are affected, so each unit's observed epochs are
contiguous.
# Do the spikes and the iEEG share a clock? Spike-triggered average of the bipolar channel of the
# unit's own amygdala, against the same with spike times shifted by a random 1 to 5 s.
rng_check = np.random.default_rng(0)
lags = np.arange(-200, 201) # samples, +-100 ms
sta_rows = []
lfp_of = {(s["subject"], lab): s["lfp"][:, c] for s in sessions for c, lab in enumerate(s["labels"])}
for u in U:
x = lfp_of[(u["subject"], u["label"])]
x = (x - x.mean()) / x.std()
idx = np.round(u["spikes"] * A.FS).astype(int)
idx = idx[(idx > 400) & (idx < len(x) - 400)]
sta = x[idx[:, None] + lags].mean(0)
sh = idx + rng_check.integers(2000, 10000)
sh = sh[sh < len(x) - 400]
sta_sh = x[sh[:, None] + lags].mean(0)
sta_rows.append({"subject": u["subject"], "unit": u["unit"], "peak_abs_sta_within_5ms": np.abs(sta[190:211]).max(),
"peak_abs_sta_shifted_within_5ms": np.abs(sta_sh[190:211]).max()})
sta_df = pd.DataFrame(sta_rows)
print(sta_df.groupby("subject")[["peak_abs_sta_within_5ms", "peak_abs_sta_shifted_within_5ms"]].median().round(3).to_string())
print("Units whose STA peak within +-5 ms exceeds the shifted one:",
int((sta_df.peak_abs_sta_within_5ms > sta_df.peak_abs_sta_shifted_within_5ms).sum()), "of", len(sta_df))
print("Wilcoxon, true against shifted:", stats.wilcoxon(sta_df.peak_abs_sta_within_5ms, sta_df.peak_abs_sta_shifted_within_5ms))
peak_abs_sta_within_5ms peak_abs_sta_shifted_within_5ms subject 1 0.019 0.052 2 0.143 0.038 3 0.045 0.024 4 0.063 0.039 7 0.030 0.056 Units whose STA peak within +-5 ms exceeds the shifted one: 24 of 42 Wilcoxon, true against shifted: WilcoxonResult(statistic=np.float64(289.0), pvalue=np.float64(0.04195783568411571))
The spike-triggered average of the macro bipolar iEEG within 5 ms of the spike is larger than at shifted spike times in 24 of 42 units (Wilcoxon p = 0.042), with medians by subject from 0.02 to 0.14 SD against 0.02 to 0.06. The lock is weak; the microwires protrude a few millimetres beyond the bipolar contacts, so a weak lock is not evidence of an offset. The shared clock of spikes, iEEG and trial table is consistent with the data (the firing and high gamma responses change at the block onset of the trial table, in the claim sections below) but not shown.
iEEG power: envelopes and the group test¶
Envelopes per band (Methods 4.8) for each of the 14 amygdalae, as the helper computes them. The
three readings of "z-scored all data within each trial" that the claims file lists decide the
outcomes: trial (over the whole 26 s epoch), baseline (mean and SD of the 2 s baseline clip of
the epoch) and amygdala (over all epochs of the amygdala). Two variants are shown beside them and
do not enter any count of readings: trialsd (the baseline-corrected envelope divided by its SD
within the trial; under trial the baseline correction of Methods 4.8 has no effect, as a constant
subtracted before a z-score) and baselinepool (the epoch's baseline mean, with the SD of the
baseline clips pooled over all epochs of the amygdala). The second is shown because the SD of
baseline rests on 20 samples of a 2 s clip smoothed over 1 s, roughly two independent values, so
baseline z-values are noisy; baselinepool shows which baseline results survive a stable SD.
The group test is run over the 24 s block (0 to 24 s from onset, the time base of Fig. 1, "locked to stimulus onset"); the baseline clip is the same in both conditions and follows a block of the other condition. This window is a choice stated here.
HGP_VARIANTS = {"hgp_60_150": (60.0, 150.0), "hgp_60_250": (60.0, 250.0)} # variants of the open upper edge
t_env = A.times()
block = t_env >= 0
tb = t_env[block]
ENV = {}
for s in sessions:
for c, lab in enumerate(s["labels"]):
n = A.table2_number(s["subject"], lab)
ep = A.epochs(s["lfp"][:, c])
for band in A.BANDS:
ENV[(n, band)] = A.envelope(ep, band)
for name, edges in HGP_VARIANTS.items():
ENV[(n, name)] = A.envelope(ep, edges)
aversive_of = {s["subject"]: s["aversive"] for s in sessions}
NSOZ = [n for n in range(1, 15) if n not in A.TABLE2_SOZ]
SOZ = sorted(A.TABLE2_SOZ)
ALL = list(range(1, 15))
ZS = ["trial", "baseline", "amygdala"] # the claims file's readings: these decide outcomes
ZV = ["trialsd", "baselinepool"] # variants, shown only
ZALL = ZS + ZV
SETS = {"nSOZ": NSOZ, "all14": ALL}
def stack(band, z, amygdalae, base="1s"):
Y, X, G, L = [], [], [], []
for n in amygdalae:
subj = A.TABLE2[n - 1][0]
y = A.normalise(ENV[(n, band)], z, base)[:, block]
Y.append(y); X.append(aversive_of[subj].astype(int)); G += [subj] * len(y); L += [n] * len(y)
return np.vstack(Y), np.concatenate(X), np.array(G), np.array(L)
print("envelopes:", len(ENV), "(14 amygdalae x 6 bands and 2 HGP variants), each", ENV[(1, "delta")].shape, "epochs x samples at",
A.STEP_S, "s; block samples", block.sum())
# Check of the fast mixed-model t against statsmodels MixedLM at a few time points.
y, x, g, l = stack("delta", "trialsd", ALL)
t_fast, df_fast = A.lmm_t(y, x, g)
cmp = []
for i in [10, 60, 120, 200]:
d = pd.DataFrame({"power": y[:, i], "condition": x, "patient": g})
fit = smf.mixedlm("power ~ condition", d, groups=d["patient"]).fit(reml=True)
cmp.append({"time_s": round(tb[i], 2), "t_fast": round(t_fast[i], 3), "t_mixedlm": round(fit.tvalues["condition"], 3),
"beta_fast": round(float(np.polyfit(x, y[:, i], 1)[0]), 3), "beta_mixedlm": round(fit.params["condition"], 3)})
print(pd.DataFrame(cmp).to_string(index=False))
# How many independent values the SD of the `baseline` reading rests on: the 2 s clip holds 20 samples
# of an envelope smoothed over 1 s. Lag-1 to lag-10 autocorrelation of the delta and HGP envelopes in the clip:
for band in ["delta", "hgp"]:
e = np.vstack([ENV[(n, band)][:, ~block] for n in ALL])
e = e - e.mean(1, keepdims=True)
ac = [np.mean((e[:, :-k] * e[:, k:]).sum(1) / (e ** 2).sum(1)) for k in (1, 5, 10)]
print(f"{band}: baseline-clip autocorrelation at 0.1, 0.5, 1 s lag {np.round(ac, 2).tolist()}; "
f"ratio of the largest to the smallest baseline SD over trials, median over amygdalae "
f"{np.median([np.ptp(np.log(ENV[(n, band)][:, ~block].std(1))) for n in ALL]):.2f} (natural log)")
envelopes: 112 (14 amygdalae x 6 bands and 2 HGP variants), each (17, 260) epochs x samples at 0.1 s; block samples 240 time_s t_fast t_mixedlm beta_fast beta_mixedlm 1.05 -2.125 -2.123 -0.378 -0.378 6.05 -2.594 -2.611 -0.488 -0.488 12.05 -3.730 -3.732 -0.618 -0.618 20.05 -2.185 -2.181 -0.403 -0.403 delta: baseline-clip autocorrelation at 0.1, 0.5, 1 s lag [0.83, 0.03, -0.32]; ratio of the largest to the smallest baseline SD over trials, median over amygdalae 1.88 (natural log) hgp: baseline-clip autocorrelation at 0.1, 0.5, 1 s lag [0.72, 0.02, -0.32]; ratio of the largest to the smallest baseline SD over trials, median over amygdalae 1.50 (natural log)
The condition effect and its t from the helper's within-patient formula match statsmodels' MixedLM
with a random intercept per patient at the time points shown (the small differences in t come from
MixedLM's variance estimate when the between-patient variance is near zero). The envelope in the 2 s
baseline clip decorrelates within about 0.5 s (autocorrelation 0.03 at 0.5 s), so the SD of the
baseline reading rests on a few independent values, and it varies across the trials of one
amygdala by a factor of about 6.5 for delta and 4.5 for HGP (median over amygdalae of the ratio of
the largest to the smallest); baselinepool replaces it with the SD pooled over the amygdala's
baseline clips.
def run_band(band, rng, zs=None):
"""All readings of the group test for one band: z-scoring x amygdala set x cluster statistic, for
the mixed model and the unpaired t-test variant."""
rows = []
for z in (zs or ZALL):
for set_name, amys in SETS.items():
y, x, g, l = stack(band, z, amys)
diff = y[x == 1].mean() - y[x == 0].mean()
for model in ("lmm", "ttest"):
pg = g if model == "lmm" else np.zeros(len(g), int)
for stat in ("mass", "size"):
_, df, res = A.cluster_test(y, x, g, rng, model=model, stat=stat, perm_groups=pg,
linked=l if model == "lmm" else None)
sig = [r for r in res if r["p"] < 0.05]
rows.append({"band": band, "z": z, "reading": "claims" if z in ZS else "variant",
"set": set_name, "model": model, "stat": stat,
"block_mean_A_minus_N": round(diff, 3),
"sig_clusters": "; ".join(f"{'A>N' if r['sign'] > 0 else 'A<N'} {tb[r['start']] - 0.05:.1f} to {tb[r['stop'] - 1] + 0.05:.1f} s p={r['p']:.3f}" for r in sig) or "none",
"n_sig_neg": sum(r["sign"] < 0 for r in sig), "n_sig_pos": sum(r["sign"] > 0 for r in sig),
"s_neg": round(sum((r["stop"] - r["start"]) * A.STEP_S for r in sig if r["sign"] < 0), 1),
"s_pos": round(sum((r["stop"] - r["start"]) * A.STEP_S for r in sig if r["sign"] > 0), 1),
"first_pos_start": min([tb[r["start"]] - 0.05 for r in sig if r["sign"] > 0], default=np.nan),
"smallest_p": round(min([r["p"] for r in res], default=1.0), 3),
"clusters": res})
return pd.DataFrame(rows)
SHOW = ["band", "z", "reading", "set", "model", "stat", "block_mean_A_minus_N", "n_sig_neg", "n_sig_pos", "s_neg", "s_pos", "smallest_p", "sig_clusters"]
fig1b-delta-decrease¶
In amygdalae outside the SOZ, delta power (1-4 Hz) is lower during AVERSIVE than during NEUTRAL sequences, and the reduction lasts through the AVERSIVE sequence.
Criterion: under the mixed model, at least one cluster with AVERSIVE below NEUTRAL at permutation p < 0.05 and the AVERSIVE block mean below the NEUTRAL one, under every reading; blocked as methods_underspecified when the readings disagree; partially reproduced when the clusters cover less than half of the 24 s block.
G = {}
G["delta"] = run_band("delta", np.random.default_rng(101))
d = G["delta"]
print(d[SHOW].to_string(index=False))
main = d[d.model == "lmm"].copy()
main["holds"] = (main.n_sig_neg > 0) & (main.block_mean_A_minus_N < 0)
mc = main[main.reading == "claims"]
print("\nMixed-model readings of the claims file where the claim holds:", int(mc.holds.sum()), "of", len(mc))
print(main.groupby(["z", "set"]).holds.all().unstack().to_string())
print("\nSeconds of the block covered by significant negative clusters (mixed model, cluster mass):")
print(main[main.stat == "mass"].pivot(index="z", columns="set", values="s_neg").to_string())
print("Abstract: 'lasting up to 18 s'; half the block is 12 s.")
band z reading set model stat block_mean_A_minus_N n_sig_neg n_sig_pos s_neg s_pos smallest_p sig_clusters delta trial claims nSOZ lmm mass -0.042 0 0 0.0 0 1.000 none delta trial claims nSOZ lmm size -0.042 0 0 0.0 0 1.000 none delta trial claims nSOZ ttest mass -0.042 0 0 0.0 0 1.000 none delta trial claims nSOZ ttest size -0.042 0 0 0.0 0 1.000 none delta trial claims all14 lmm mass -0.027 0 0 0.0 0 0.878 none delta trial claims all14 lmm size -0.027 0 0 0.0 0 0.876 none delta trial claims all14 ttest mass -0.027 0 0 0.0 0 0.838 none delta trial claims all14 ttest size -0.027 0 0 0.0 0 0.872 none delta baseline claims nSOZ lmm mass -1.389 5 0 14.3 0 0.004 A<N 4.3 to 7.9 s p=0.004; A<N 8.5 to 10.7 s p=0.012; A<N 11.6 to 15.4 s p=0.006; A<N 16.3 to 19.2 s p=0.008; A<N 21.4 to 23.2 s p=0.038 delta baseline claims nSOZ lmm size -1.389 5 0 14.3 0 0.006 A<N 4.3 to 7.9 s p=0.010; A<N 8.5 to 10.7 s p=0.028; A<N 11.6 to 15.4 s p=0.006; A<N 16.3 to 19.2 s p=0.014; A<N 21.4 to 23.2 s p=0.046 delta baseline claims nSOZ ttest mass -1.389 5 0 15.5 0 0.002 A<N 3.3 to 8.0 s p=0.002; A<N 8.4 to 10.7 s p=0.014; A<N 11.6 to 15.4 s p=0.004; A<N 16.3 to 19.2 s p=0.004; A<N 21.4 to 23.2 s p=0.036 delta baseline claims nSOZ ttest size -1.389 4 0 13.7 0 0.004 A<N 3.3 to 8.0 s p=0.004; A<N 8.4 to 10.7 s p=0.018; A<N 11.6 to 15.4 s p=0.006; A<N 16.3 to 19.2 s p=0.010 delta baseline claims all14 lmm mass -1.156 2 0 15.6 0 0.002 A<N 6.7 to 20.5 s p=0.002; A<N 21.4 to 23.2 s p=0.048 delta baseline claims all14 lmm size -1.156 1 0 13.8 0 0.002 A<N 6.7 to 20.5 s p=0.002 delta baseline claims all14 ttest mass -1.156 2 0 15.6 0 0.002 A<N 6.7 to 20.5 s p=0.002; A<N 21.4 to 23.2 s p=0.028 delta baseline claims all14 ttest size -1.156 1 0 13.8 0 0.002 A<N 6.7 to 20.5 s p=0.002 delta amygdala claims nSOZ lmm mass -0.350 0 0 0.0 0 0.072 none delta amygdala claims nSOZ lmm size -0.350 0 0 0.0 0 0.056 none delta amygdala claims nSOZ ttest mass -0.350 0 0 0.0 0 0.064 none delta amygdala claims nSOZ ttest size -0.350 0 0 0.0 0 0.072 none delta amygdala claims all14 lmm mass -0.337 2 0 17.1 0 0.002 A<N 5.9 to 20.3 s p=0.002; A<N 21.3 to 24.0 s p=0.026 delta amygdala claims all14 lmm size -0.337 2 0 17.1 0 0.002 A<N 5.9 to 20.3 s p=0.002; A<N 21.3 to 24.0 s p=0.010 delta amygdala claims all14 ttest mass -0.337 2 0 17.1 0 0.002 A<N 5.9 to 20.3 s p=0.002; A<N 21.3 to 24.0 s p=0.012 delta amygdala claims all14 ttest size -0.337 2 0 17.1 0 0.002 A<N 5.9 to 20.3 s p=0.002; A<N 21.3 to 24.0 s p=0.010 delta trialsd variant nSOZ lmm mass -0.498 4 0 7.0 0 0.034 A<N 5.6 to 7.4 s p=0.034; A<N 8.8 to 10.7 s p=0.036; A<N 11.8 to 13.4 s p=0.044; A<N 21.6 to 23.3 s p=0.042 delta trialsd variant nSOZ lmm size -0.498 2 0 3.7 0 0.040 A<N 5.6 to 7.4 s p=0.048; A<N 8.8 to 10.7 s p=0.040 delta trialsd variant nSOZ ttest mass -0.498 1 0 2.6 0 0.022 A<N 21.4 to 24.0 s p=0.022 delta trialsd variant nSOZ ttest size -0.498 3 0 6.4 0 0.018 A<N 5.6 to 7.5 s p=0.044; A<N 8.8 to 10.7 s p=0.044; A<N 21.4 to 24.0 s p=0.018 delta trialsd variant all14 lmm mass -0.445 2 0 16.1 0 0.002 A<N 5.8 to 19.3 s p=0.002; A<N 21.4 to 24.0 s p=0.030 delta trialsd variant all14 lmm size -0.445 2 0 16.1 0 0.002 A<N 5.8 to 19.3 s p=0.002; A<N 21.4 to 24.0 s p=0.022 delta trialsd variant all14 ttest mass -0.445 2 0 16.1 0 0.002 A<N 5.8 to 19.3 s p=0.002; A<N 21.4 to 24.0 s p=0.014 delta trialsd variant all14 ttest size -0.445 2 0 16.1 0 0.002 A<N 5.8 to 19.3 s p=0.002; A<N 21.4 to 24.0 s p=0.030 delta baselinepool variant nSOZ lmm mass -0.768 5 0 13.3 0 0.002 A<N 3.1 to 5.3 s p=0.026; A<N 5.6 to 7.5 s p=0.022; A<N 8.7 to 10.8 s p=0.022; A<N 11.6 to 16.1 s p=0.002; A<N 21.4 to 24.0 s p=0.008 delta baselinepool variant nSOZ lmm size -0.768 5 0 13.3 0 0.002 A<N 3.1 to 5.3 s p=0.026; A<N 5.6 to 7.5 s p=0.036; A<N 8.7 to 10.8 s p=0.030; A<N 11.6 to 16.1 s p=0.002; A<N 21.4 to 24.0 s p=0.018 delta baselinepool variant nSOZ ttest mass -0.768 5 0 13.3 0 0.002 A<N 3.2 to 5.3 s p=0.036; A<N 5.6 to 7.5 s p=0.026; A<N 8.7 to 10.8 s p=0.026; A<N 11.6 to 16.1 s p=0.002; A<N 21.3 to 24.0 s p=0.016 delta baselinepool variant nSOZ ttest size -0.768 5 0 13.3 0 0.002 A<N 3.2 to 5.3 s p=0.028; A<N 5.6 to 7.5 s p=0.040; A<N 8.7 to 10.8 s p=0.028; A<N 11.6 to 16.1 s p=0.002; A<N 21.3 to 24.0 s p=0.010 delta baselinepool variant all14 lmm mass -0.589 2 0 16.8 0 0.002 A<N 6.0 to 20.3 s p=0.002; A<N 21.3 to 23.8 s p=0.018 delta baselinepool variant all14 lmm size -0.589 2 0 16.8 0 0.002 A<N 6.0 to 20.3 s p=0.002; A<N 21.3 to 23.8 s p=0.028 delta baselinepool variant all14 ttest mass -0.589 2 0 16.9 0 0.002 A<N 6.0 to 20.3 s p=0.002; A<N 21.3 to 23.9 s p=0.012 delta baselinepool variant all14 ttest size -0.589 2 0 16.9 0 0.002 A<N 6.0 to 20.3 s p=0.002; A<N 21.3 to 23.9 s p=0.012 Mixed-model readings of the claims file where the claim holds: 6 of 12 set all14 nSOZ z amygdala True False baseline True True baselinepool True True trial False False trialsd True True Seconds of the block covered by significant negative clusters (mixed model, cluster mass): set all14 nSOZ z amygdala 17.1 0.0 baseline 15.6 14.3 baselinepool 16.8 13.3 trial 0.0 0.0 trialsd 16.1 7.0 Abstract: 'lasting up to 18 s'; half the block is 12 s.
fig, axes = plt.subplots(1, 2, figsize=(11, 3.5), sharey=True)
yrng = []
for ax, (set_name, amys) in zip(axes, SETS.items()):
for z, col in zip(ZALL, ["k", "C0", "C1", "C2", "C3"]):
y, x, g, l = stack("delta", z, amys)
curve = y[x == 1].mean(0) - y[x == 0].mean(0)
ax.plot(tb, curve, color=col, label=z, ls="-" if z in ZS else ":")
yrng.append((z, set_name, round(float(curve.min()), 2), round(float(curve.max()), 2)))
ax.axhline(0, color="grey", lw=0.5)
ax.set_title(f"delta, {set_name}: AVERSIVE minus NEUTRAL")
ax.set_xlabel("time from block onset (s)")
axes[0].set_ylabel("z (each reading's units)")
axes[0].legend(fontsize=8)
plt.tight_layout()
print("Range of each plotted curve (z, set, min, max):", yrng)
Range of each plotted curve (z, set, min, max): [('trial', 'nSOZ', -0.34, 0.39), ('baseline', 'nSOZ', -2.51, -0.18), ('amygdala', 'nSOZ', -0.59, -0.03), ('trialsd', 'nSOZ', -0.8, -0.07), ('baselinepool', 'nSOZ', -1.23, -0.12), ('trial', 'all14', -0.3, 0.3), ('baseline', 'all14', -1.89, -0.06), ('amygdala', 'all14', -0.54, -0.09), ('trialsd', 'all14', -0.72, -0.11), ('baselinepool', 'all14', -0.98, -0.06)]
Result: blocked (methods_underspecified). The readings of the claims file disagree, as the
criterion foresees: the claim holds in 6 of the 12 mixed-model readings. With the z-score over the
whole trial (trial) no reading finds a delta cluster, in either amygdala set (smallest cluster p
0.88 and 1.0). With baseline z-scoring AVERSIVE is below NEUTRAL in the 7 nSOZ amygdalae
(significant clusters covering 14.3 s of the block) and in all 14 (15.6 s). With amygdala
z-scoring the 14 amygdalae give 17.1 s and the 7 nSOZ amygdalae none (smallest p 0.072). The block
mean of AVERSIVE minus NEUTRAL is negative under every reading (-0.03 to -1.39). Where clusters are
found they cover 14 to 17 s of the 24 s block, against the abstract's "up to 18 s". The variants
agree with baseline: baselinepool 13.3 s (nSOZ) and 16.8 s (all 14), trialsd 7.0 s and 16.1 s.
The unpaired t-test variant gives the same pattern.
fig1c-hgp-increase-first-2s¶
High gamma power (above 60 Hz) is higher during AVERSIVE than during NEUTRAL sequences in the first 2 s of the sequence.
Criterion: as fig1b, for HGP; reproduced when a cluster with AVERSIVE above NEUTRAL at p < 0.05 lies within, or begins within, the first 2 s after onset, under every reading.
G["hgp"] = run_band("hgp", np.random.default_rng(102))
d = G["hgp"]
print(d[SHOW + ["first_pos_start"]].to_string(index=False))
lm = d[(d.model == "lmm") & (d.reading == "claims")]
tt_ = d[(d.model == "ttest") & (d.reading == "claims")]
print("\nMixed-model readings of the claims file with a positive cluster (p < 0.05) beginning in the first 2 s:",
int(((lm.n_sig_pos > 0) & (lm.first_pos_start < 2)).sum()), "of", len(lm))
print("Unpaired t-test readings with one:", int(((tt_.n_sig_pos > 0) & (tt_.first_pos_start < 2)).sum()), "of", len(tt_))
for set_name, amys in SETS.items():
for z in ZALL:
y, x, g, l = stack("hgp", z, amys)
tt, df = A.lmm_t(y, x, g)
crit = stats.t.ppf(0.975, df)
print(f"{set_name:6s} {z:12s} largest mixed-model t in 0-2 s: {tt[tb < 2].max():.2f} at {tb[tb < 2][tt[tb < 2].argmax()]:.2f} s "
f"(threshold {crit:.2f}); points above threshold in 0-2 s: {int((tt[tb < 2] > crit).sum())}; in 2-24 s: {int((tt[tb >= 2] > crit).sum())}")
band z reading set model stat block_mean_A_minus_N n_sig_neg n_sig_pos s_neg s_pos smallest_p sig_clusters first_pos_start hgp trial claims nSOZ lmm mass 0.019 0 0 0 0 0.499 none NaN hgp trial claims nSOZ lmm size 0.019 0 0 0 0 0.519 none NaN hgp trial claims nSOZ ttest mass 0.019 0 0 0 0 0.493 none NaN hgp trial claims nSOZ ttest size 0.019 0 0 0 0 0.515 none NaN hgp trial claims all14 lmm mass -0.004 0 0 0 0 0.210 none NaN hgp trial claims all14 lmm size -0.004 0 0 0 0 0.279 none NaN hgp trial claims all14 ttest mass -0.004 0 0 0 0 0.162 none NaN hgp trial claims all14 ttest size -0.004 0 0 0 0 0.232 none NaN hgp baseline claims nSOZ lmm mass 0.532 0 0 0 0 0.062 none NaN hgp baseline claims nSOZ lmm size 0.532 0 0 0 0 0.056 none NaN hgp baseline claims nSOZ ttest mass 0.532 0 0 0 0 0.058 none NaN hgp baseline claims nSOZ ttest size 0.532 0 0 0 0 0.064 none NaN hgp baseline claims all14 lmm mass -0.309 0 0 0 0 0.729 none NaN hgp baseline claims all14 lmm size -0.309 0 0 0 0 0.842 none NaN hgp baseline claims all14 ttest mass -0.309 0 0 0 0 0.663 none NaN hgp baseline claims all14 ttest size -0.309 0 0 0 0 0.784 none NaN hgp amygdala claims nSOZ lmm mass 0.130 0 0 0 0 0.160 none NaN hgp amygdala claims nSOZ lmm size 0.130 0 0 0 0 0.160 none NaN hgp amygdala claims nSOZ ttest mass 0.130 0 0 0 0 0.158 none NaN hgp amygdala claims nSOZ ttest size 0.130 0 0 0 0 0.198 none NaN hgp amygdala claims all14 lmm mass -0.104 0 0 0 0 0.481 none NaN hgp amygdala claims all14 lmm size -0.104 0 0 0 0 0.561 none NaN hgp amygdala claims all14 ttest mass -0.104 0 0 0 0 0.383 none NaN hgp amygdala claims all14 ttest size -0.104 0 0 0 0 0.385 none NaN hgp trialsd variant nSOZ lmm mass 0.166 0 0 0 0 0.369 none NaN hgp trialsd variant nSOZ lmm size 0.166 0 0 0 0 0.379 none NaN hgp trialsd variant nSOZ ttest mass 0.166 0 0 0 0 0.373 none NaN hgp trialsd variant nSOZ ttest size 0.166 0 0 0 0 0.399 none NaN hgp trialsd variant all14 lmm mass -0.136 0 0 0 0 0.531 none NaN hgp trialsd variant all14 lmm size -0.136 0 0 0 0 0.585 none NaN hgp trialsd variant all14 ttest mass -0.136 0 0 0 0 0.559 none NaN hgp trialsd variant all14 ttest size -0.136 0 0 0 0 0.667 none NaN hgp baselinepool variant nSOZ lmm mass 0.442 0 0 0 0 0.102 none NaN hgp baselinepool variant nSOZ lmm size 0.442 0 0 0 0 0.148 none NaN hgp baselinepool variant nSOZ ttest mass 0.442 0 0 0 0 0.118 none NaN hgp baselinepool variant nSOZ ttest size 0.442 0 0 0 0 0.148 none NaN hgp baselinepool variant all14 lmm mass -0.149 0 0 0 0 0.717 none NaN hgp baselinepool variant all14 lmm size -0.149 0 0 0 0 0.780 none NaN hgp baselinepool variant all14 ttest mass -0.149 0 0 0 0 0.651 none NaN hgp baselinepool variant all14 ttest size -0.149 0 0 0 0 0.737 none NaN Mixed-model readings of the claims file with a positive cluster (p < 0.05) beginning in the first 2 s: 0 of 12 Unpaired t-test readings with one: 0 of 12 nSOZ trial largest mixed-model t in 0-2 s: 2.05 at 0.95 s (threshold 1.98); points above threshold in 0-2 s: 1; in 2-24 s: 7 nSOZ baseline largest mixed-model t in 0-2 s: 3.01 at 0.95 s (threshold 1.98); points above threshold in 0-2 s: 13; in 2-24 s: 27 nSOZ amygdala largest mixed-model t in 0-2 s: 2.10 at 0.95 s (threshold 1.98); points above threshold in 0-2 s: 2; in 2-24 s: 11 nSOZ trialsd largest mixed-model t in 0-2 s: 1.94 at 0.95 s (threshold 1.98); points above threshold in 0-2 s: 0; in 2-24 s: 8 nSOZ baselinepool largest mixed-model t in 0-2 s: 2.78 at 0.95 s (threshold 1.98); points above threshold in 0-2 s: 9; in 2-24 s: 33 all14 trial largest mixed-model t in 0-2 s: 1.60 at 1.05 s (threshold 1.97); points above threshold in 0-2 s: 0; in 2-24 s: 8 all14 baseline largest mixed-model t in 0-2 s: 1.73 at 1.05 s (threshold 1.97); points above threshold in 0-2 s: 0; in 2-24 s: 1 all14 amygdala largest mixed-model t in 0-2 s: 0.54 at 1.05 s (threshold 1.97); points above threshold in 0-2 s: 0; in 2-24 s: 1 all14 trialsd largest mixed-model t in 0-2 s: 0.34 at 1.05 s (threshold 1.97); points above threshold in 0-2 s: 0; in 2-24 s: 0 all14 baselinepool largest mixed-model t in 0-2 s: 0.86 at 1.05 s (threshold 1.97); points above threshold in 0-2 s: 0; in 2-24 s: 0
# Variants of the open upper edge of the HGP band (the claims file's reading is a high-pass at 60 Hz):
# 60-150 Hz and 60-250 Hz band-pass. Shown, not used for the outcome.
for name in HGP_VARIANTS:
G[name] = run_band(name, np.random.default_rng(110 + list(HGP_VARIANTS).index(name)))
v = G[name][G[name].model == "lmm"]
print(f"{name}: mixed-model readings with a positive cluster starting in 0-2 s: "
f"{int(((v.n_sig_pos > 0) & (v.first_pos_start < 2)).sum())} of {len(v)}; any significant cluster: "
f"{int(((v.n_sig_pos + v.n_sig_neg) > 0).sum())}; smallest cluster p {v.smallest_p.min():.3f}")
print(v[v.sig_clusters != "none"][["z", "set", "stat", "sig_clusters"]].to_string(index=False))
hgp_60_150: mixed-model readings with a positive cluster starting in 0-2 s: 0 of 20; any significant cluster: 4; smallest cluster p 0.020
z set stat sig_clusters
baseline nSOZ mass A>N 6.2 to 7.8 s p=0.020
baseline nSOZ size A>N 6.2 to 7.8 s p=0.024
baselinepool nSOZ mass A>N 3.1 to 4.7 s p=0.040; A>N 6.2 to 7.8 s p=0.032
baselinepool nSOZ size A>N 3.1 to 4.7 s p=0.036; A>N 6.2 to 7.8 s p=0.036
hgp_60_250: mixed-model readings with a positive cluster starting in 0-2 s: 0 of 20; any significant cluster: 3; smallest cluster p 0.028
z set stat sig_clusters
baseline nSOZ mass A>N 6.4 to 7.8 s p=0.042
baselinepool nSOZ mass A>N 3.1 to 4.7 s p=0.028; A>N 6.4 to 7.8 s p=0.040
baselinepool nSOZ size A>N 3.1 to 4.7 s p=0.048
fig, axes = plt.subplots(1, 2, figsize=(11, 3.5))
for ax, (set_name, amys) in zip(axes, SETS.items()):
for z, col in zip(ZALL, ["k", "C0", "C1", "C2", "C3"]):
y, x, g, l = stack("hgp", z, amys)
tt, df = A.lmm_t(y, x, g)
ax.plot(tb, tt, color=col, label=z, ls="-" if z in ZS else ":")
ax.axhline(stats.t.ppf(0.975, df), color="grey", ls="--", lw=0.8)
ax.axhline(-stats.t.ppf(0.975, df), color="grey", ls="--", lw=0.8)
ax.axvline(2, color="r", lw=0.5)
ax.set_title(f"HGP, {set_name}: mixed-model t, AVERSIVE vs NEUTRAL")
ax.set_xlabel("time from block onset (s)")
axes[0].legend(fontsize=8)
plt.tight_layout()
print("All plotted t values are inside the axes (autoscaled).")
All plotted t values are inside the axes (autoscaled).
Result: not reproduced. None of the 12 mixed-model readings and none of the 12 unpaired t-test
readings of the claims file finds a significant HGP cluster, in the first 2 s or later, and neither
do the two z-scoring variants. The largest mixed-model t in the first 2 s is 3.01, at 0.95 s after
onset (baseline, nSOZ), above the cluster-forming threshold at 13 points, and that cluster's
permutation p is 0.062 (smallest of the readings 0.056). With the band limited to 60-150 or
60-250 Hz (variants) the only significant HGP clusters are AVERSIVE above NEUTRAL from about 6 to 8 s
(baseline, nSOZ, p 0.020 and 0.042) and, with baselinepool, 3.1 to 4.7 s, not in the first 2 s.
Per-amygdala time courses and the other baselines of Methods 4.8¶
Each amygdala's AVERSIVE minus NEUTRAL difference in 2 s bins of the block, for delta and HGP under
the three z-scorings of the claims file, to compare with the group curves of Fig. 1B and 1C. Then
the two other baselines Methods 4.8 say were tested ("taking as baseline the mean of all epochs, or
the mean of the previous epoch. Results remained quantitatively similar"), which the claims file
lists as variants. Both subtract one constant per trial, which a z-score over the trial (trial) or
over the baseline clip (baseline) removes, so they can change only amygdala and the variants.
All are run, mixed model with cluster mass.
edges2 = np.arange(0, 25, 2)
for band in ["delta", "hgp"]:
for z in ZS:
rows = []
for n in ALL:
yv = A.normalise(ENV[(n, band)], z)[:, block]
av = aversive_of[A.TABLE2[n - 1][0]]
dif = yv[av].mean(0) - yv[~av].mean(0)
rows.append({"amygdala": n, "soz": n in A.TABLE2_SOZ,
**{f"{a}-{a + 2}": round(float(dif[(tb >= a) & (tb < a + 2)].mean()), 2) for a in edges2[:-1]}})
tab = pd.DataFrame(rows)
print(f"\n{band}, z-scoring {z}: AVERSIVE minus NEUTRAL per amygdala in 2 s bins (s from onset)")
print(tab.to_string(index=False))
print(" amygdalae with a negative 0-24 s mean:", int((tab.iloc[:, 2:].mean(1) < 0).sum()), "of 14; nSOZ with a positive 0-2 s bin:",
int((tab[~tab.soz]["0-2"] > 0).sum()), "of 7; SOZ:", int((tab[tab.soz]["0-2"] > 0).sum()), "of 7")
rng_b = np.random.default_rng(120)
rows = []
for band in ["delta", "hgp"]:
for base in ["1s", "allepochs", "prevepoch"]:
for z in ZALL:
for set_name, amys in SETS.items():
y, x, g, l = stack(band, z, amys, base)
_, _, res = A.cluster_test(y, x, g, rng_b, model="lmm", stat="mass", perm_groups=g, linked=l)
sig = [r for r in res if r["p"] < 0.05]
rows.append({"band": band, "baseline": base, "z": z, "set": set_name,
"block_mean_A_minus_N": round(float(y[x == 1].mean() - y[x == 0].mean()), 3),
"smallest_p": round(min([r["p"] for r in res], default=1.0), 3),
"sig_clusters": "; ".join(f"{'A>N' if r['sign'] > 0 else 'A<N'} {tb[r['start']] - 0.05:.1f} to {tb[r['stop'] - 1] + 0.05:.1f} s p={r['p']:.3f}" for r in sig) or "none"})
rb = pd.DataFrame(rows)
print("\nBaselines of Methods 4.8 (mixed model, cluster mass):")
print(rb.to_string(index=False))
delta, z-scoring trial: AVERSIVE minus NEUTRAL per amygdala in 2 s bins (s from onset)
amygdala soz 0-2 2-4 4-6 6-8 8-10 10-12 12-14 14-16 16-18 18-20 20-22 22-24
1 False -0.18 0.55 -0.47 0.07 0.05 -0.02 -0.76 -0.35 0.02 0.47 0.44 -0.28
2 True 0.42 0.46 -0.07 -0.17 -0.05 -0.21 -0.23 -0.16 -0.53 0.64 0.55 -0.13
3 False -0.03 0.20 0.15 -0.24 -0.18 -0.27 -0.38 0.61 -0.01 -0.05 0.18 -0.52
4 True -0.03 0.30 0.46 -0.29 -0.69 0.08 -0.14 0.13 0.01 -0.32 0.02 0.14
5 False 0.16 -0.28 -0.48 -0.44 -0.15 0.50 0.25 -0.44 -0.13 -0.21 0.19 0.47
6 False -0.01 -0.51 -0.22 -0.33 -0.55 -0.21 0.05 -0.40 0.17 0.10 0.27 0.30
7 True 0.22 0.41 0.73 0.24 -0.25 0.14 -0.09 -0.52 -0.51 -0.13 -0.24 -0.46
8 False 0.46 0.04 -0.03 -0.41 0.15 -0.30 0.25 -0.25 0.06 0.49 -0.65 -0.63
9 True 0.40 0.30 0.29 0.10 0.03 -0.61 0.31 0.00 0.13 -0.48 -0.30 -0.13
10 False -0.08 -0.11 0.33 0.35 0.07 -0.01 -0.36 -0.12 -0.26 0.08 0.59 0.26
11 True -0.27 -0.37 -0.19 0.09 0.19 0.08 -0.26 0.27 0.04 -0.53 0.81 0.30
12 False 0.07 0.05 0.43 0.10 0.18 0.26 0.05 0.37 -0.58 -0.05 -0.45 -0.91
13 True -0.06 -0.40 -0.02 -0.26 0.17 -0.56 0.05 0.20 0.27 0.17 -0.09 -0.09
14 True -0.10 0.38 0.10 0.08 -0.08 -0.33 -0.46 -0.11 0.10 0.03 -0.07 0.03
amygdalae with a negative 0-24 s mean: 10 of 14; nSOZ with a positive 0-2 s bin: 3 of 7; SOZ: 3 of 7
delta, z-scoring baseline: AVERSIVE minus NEUTRAL per amygdala in 2 s bins (s from onset)
amygdala soz 0-2 2-4 4-6 6-8 8-10 10-12 12-14 14-16 16-18 18-20 20-22 22-24
1 False -0.85 1.06 -2.03 -1.45 -0.74 -1.17 -3.28 -1.34 -0.90 -0.20 0.27 -1.66
2 True 1.33 1.95 1.66 1.35 1.05 0.83 0.68 0.85 0.16 2.33 2.25 0.35
3 False -1.17 -1.05 -1.49 -2.32 -2.06 -2.20 -2.53 -0.80 -1.65 -1.86 -1.72 -2.16
4 True -0.98 0.07 0.12 -1.73 -2.28 -0.66 -1.48 -0.38 -1.02 -1.13 -0.76 0.22
5 False -0.38 -1.59 -2.60 -2.47 -2.07 0.04 -1.59 -2.54 -0.98 -1.22 -1.61 -0.26
6 False -1.77 -2.52 -1.98 -2.67 -3.62 -2.84 -1.65 -2.54 -1.40 -2.23 -1.24 -1.78
7 True -2.56 -3.36 1.98 0.49 -5.67 -3.45 -6.82 -7.95 -6.97 -5.78 -4.10 -4.02
8 False -1.41 -3.12 -3.72 -3.63 -2.61 -3.24 -1.18 -2.85 -3.59 -3.03 -4.96 -4.31
9 True 0.24 0.59 0.34 0.31 -0.71 -1.94 1.06 -0.10 -0.06 -0.24 -0.36 0.03
10 False 0.91 1.39 1.95 1.64 1.03 0.97 0.55 -0.85 0.34 1.44 2.36 2.35
11 True -0.12 -0.28 0.01 0.07 0.56 0.13 -0.41 0.72 0.33 -0.35 1.48 0.57
12 False -1.25 -0.95 0.28 -1.61 -1.08 0.02 -1.22 0.43 -2.96 -1.18 -2.42 -4.35
13 True -1.50 -2.10 -1.77 -2.20 -0.88 -2.80 -1.37 -1.69 -1.62 -0.89 -1.86 -2.02
14 True -0.90 0.04 -0.38 -0.68 -0.68 -1.53 -1.58 -2.10 -2.43 -1.73 -1.95 -1.29
amygdalae with a negative 0-24 s mean: 11 of 14; nSOZ with a positive 0-2 s bin: 1 of 7; SOZ: 2 of 7
delta, z-scoring amygdala: AVERSIVE minus NEUTRAL per amygdala in 2 s bins (s from onset)
amygdala soz 0-2 2-4 4-6 6-8 8-10 10-12 12-14 14-16 16-18 18-20 20-22 22-24
1 False 0.00 0.47 -0.43 -0.03 -0.03 -0.09 -0.80 -0.29 0.03 0.29 0.39 -0.24
2 True 0.35 0.42 0.07 -0.00 0.01 -0.10 -0.09 -0.04 -0.31 0.60 0.53 -0.01
3 False -0.20 -0.07 -0.00 -0.24 -0.25 -0.31 -0.41 0.26 -0.08 -0.13 0.03 -0.49
4 True -0.13 0.06 0.18 -0.31 -0.72 -0.15 -0.34 -0.05 -0.13 -0.37 -0.15 0.04
5 False -0.44 -0.68 -0.84 -0.85 -0.65 -0.20 -0.38 -0.87 -0.65 -0.71 -0.44 -0.25
6 False -0.76 -1.01 -0.79 -1.01 -1.26 -0.93 -0.70 -0.96 -0.60 -0.79 -0.57 -0.51
7 True -0.39 -0.27 0.14 -0.23 -0.71 -0.35 -0.64 -0.97 -1.07 -0.93 -0.98 -0.97
8 False 0.02 -0.37 -0.56 -0.86 -0.26 -0.54 -0.22 -0.58 -0.29 0.01 -0.93 -0.83
9 True 0.24 0.13 0.07 -0.01 -0.17 -0.63 0.16 -0.00 -0.02 -0.55 -0.29 -0.17
10 False 0.31 0.34 0.63 0.61 0.30 0.32 0.04 0.10 0.14 0.25 0.76 0.61
11 True -0.49 -0.51 -0.47 -0.22 -0.09 -0.19 -0.48 -0.12 -0.19 -0.63 0.24 -0.00
12 False -0.72 -0.64 -0.32 -0.64 -0.53 -0.45 -0.72 -0.47 -1.29 -0.72 -1.03 -1.35
13 True -0.67 -0.82 -0.56 -0.73 -0.42 -1.15 -0.58 -0.36 -0.33 -0.61 -0.54 -0.63
14 True -0.54 -0.27 -0.44 -0.50 -0.53 -0.81 -0.81 -0.82 -0.68 -0.70 -0.84 -0.54
amygdalae with a negative 0-24 s mean: 12 of 14; nSOZ with a positive 0-2 s bin: 2 of 7; SOZ: 2 of 7
hgp, z-scoring trial: AVERSIVE minus NEUTRAL per amygdala in 2 s bins (s from onset)
amygdala soz 0-2 2-4 4-6 6-8 8-10 10-12 12-14 14-16 16-18 18-20 20-22 22-24
1 False -0.07 0.35 -0.14 0.27 -0.06 -0.21 -0.16 -0.50 -0.32 0.63 0.92 0.16
2 True 0.03 0.48 -0.15 -0.31 -0.77 -0.77 -0.04 0.01 0.46 1.02 -0.29 -0.12
3 False 0.34 -0.18 -0.14 1.11 0.72 -0.13 0.16 0.01 0.07 -0.91 -0.68 -0.32
4 True 0.13 0.35 0.01 0.26 0.79 -0.06 -0.14 0.05 0.33 -0.78 -0.55 -0.66
5 False 0.82 0.63 -0.26 -0.05 -0.09 0.03 0.10 -0.25 0.15 0.09 0.34 -0.07
6 False 0.02 -0.39 -0.13 -0.11 -0.54 0.05 0.11 0.31 -0.06 0.04 0.56 0.45
7 True 0.04 0.32 0.38 0.18 0.10 -0.09 0.25 -0.49 -0.46 0.01 0.02 -0.51
8 False 0.19 0.00 0.19 0.18 -0.16 -0.19 0.12 0.39 -0.19 -0.44 -0.48 0.09
9 True -0.28 0.15 0.30 0.15 -0.05 -0.28 -0.44 0.13 -0.16 0.05 -0.07 0.04
10 False -0.17 0.40 0.51 -0.11 0.31 -0.24 -0.62 -0.10 0.62 -0.40 0.08 -0.37
11 True -0.66 0.17 0.26 0.17 -0.87 0.53 0.29 -0.22 -0.14 0.43 0.44 -0.28
12 False -0.39 0.15 0.52 -0.31 -0.22 0.19 -0.07 0.19 -0.11 -0.00 -0.22 -0.33
13 True -0.22 -0.03 -0.04 0.12 -0.22 -0.50 0.00 0.09 0.27 -0.01 0.09 -0.10
14 True 0.56 0.04 -0.21 -0.08 0.34 0.03 -0.04 -0.40 -0.03 -0.43 0.02 -0.13
amygdalae with a negative 0-24 s mean: 9 of 14; nSOZ with a positive 0-2 s bin: 4 of 7; SOZ: 4 of 7
hgp, z-scoring baseline: AVERSIVE minus NEUTRAL per amygdala in 2 s bins (s from onset)
amygdala soz 0-2 2-4 4-6 6-8 8-10 10-12 12-14 14-16 16-18 18-20 20-22 22-24
1 False 2.48 3.09 1.76 2.15 1.75 1.45 1.45 0.99 1.16 3.57 3.95 2.45
2 True -0.44 -0.20 -1.60 -1.85 -2.45 -2.50 -0.37 -0.41 -0.31 1.22 -0.77 -1.57
3 False 0.44 -0.60 -0.71 2.24 1.11 -0.26 0.26 -0.47 -0.02 -1.72 -1.09 -0.69
4 True -0.14 0.10 -0.14 -0.05 0.99 -0.21 -0.71 -0.19 0.57 -1.84 -1.55 -1.80
5 False 4.00 2.97 1.98 2.40 2.12 2.37 2.68 2.11 2.63 2.41 3.16 2.49
6 False 1.15 0.31 0.69 0.43 -0.69 0.85 1.05 1.62 -0.40 1.38 6.05 1.58
7 True -4.80 -3.07 1.34 -0.00 -2.27 -3.05 -7.97 -5.79 -11.46 -3.79 -11.97 -13.54
8 False -0.07 -0.81 -0.09 -0.42 -0.74 -0.99 -0.80 0.31 -0.87 -1.77 -1.78 -0.27
9 True 0.23 0.21 1.15 1.09 -0.54 -0.85 -0.65 -0.05 -0.58 -0.44 -0.44 0.09
10 False -0.85 0.28 0.34 -0.77 -0.27 -1.14 -1.55 -0.47 0.14 -1.37 -0.66 -1.43
11 True -0.87 0.45 0.87 0.69 -1.01 2.11 1.39 0.08 0.14 1.35 1.69 0.07
12 False -1.08 -0.42 -0.06 -1.23 -1.08 -0.47 -0.78 -0.45 -0.99 -0.82 -1.01 -0.97
13 True -0.78 -0.72 -0.68 -0.39 -0.84 -1.31 -0.42 -0.58 -0.35 -0.66 -0.68 -0.64
14 True 0.18 -0.85 -2.21 -1.31 -0.22 -0.71 -1.68 -1.40 -0.47 -1.22 -0.68 -1.58
amygdalae with a negative 0-24 s mean: 10 of 14; nSOZ with a positive 0-2 s bin: 4 of 7; SOZ: 2 of 7
hgp, z-scoring amygdala: AVERSIVE minus NEUTRAL per amygdala in 2 s bins (s from onset)
amygdala soz 0-2 2-4 4-6 6-8 8-10 10-12 12-14 14-16 16-18 18-20 20-22 22-24
1 False 0.60 0.93 0.59 0.88 0.61 0.51 0.57 0.32 0.45 1.19 1.42 0.86
2 True -0.39 -0.06 -0.53 -0.62 -1.00 -0.98 -0.50 -0.44 -0.03 0.35 -0.63 -0.46
3 False 0.33 -0.21 -0.13 1.08 0.64 -0.12 0.13 -0.09 -0.07 -0.79 -0.57 -0.30
4 True -0.57 -0.41 -0.64 -0.49 -0.12 -0.75 -0.73 -0.65 -0.42 -1.12 -0.96 -1.01
5 False 1.66 1.43 0.83 1.00 0.96 1.02 1.10 0.83 1.08 1.08 1.32 0.97
6 False 0.08 0.00 0.04 0.04 -0.35 0.13 0.19 0.25 0.07 0.18 1.44 0.28
7 True -0.38 -0.28 -0.05 -0.33 -0.41 -0.45 -0.38 -0.64 -0.81 -0.80 -0.87 -0.83
8 False 0.01 -0.14 0.06 0.04 -0.32 -0.30 -0.09 0.13 -0.31 -0.41 -0.49 -0.05
9 True -0.55 -0.38 -0.34 -0.24 -0.56 -0.76 -0.88 -0.43 -0.56 -0.56 -0.59 -0.45
10 False -0.54 -0.06 0.05 -0.48 -0.14 -0.60 -0.91 -0.49 0.12 -0.77 -0.38 -0.70
11 True -0.34 0.36 0.39 0.34 -0.52 0.66 0.47 0.06 0.08 0.54 0.52 -0.04
12 False -0.78 -0.33 -0.09 -0.75 -0.66 -0.37 -0.66 -0.38 -0.65 -0.51 -0.70 -0.84
13 True -0.27 -0.15 -0.13 0.04 -0.25 -0.43 -0.11 -0.02 0.11 -0.13 -0.01 -0.22
14 True -0.02 -0.09 -0.30 -0.37 0.03 -0.32 -0.43 -0.70 -0.31 -0.67 -0.12 -0.49
amygdalae with a negative 0-24 s mean: 10 of 14; nSOZ with a positive 0-2 s bin: 5 of 7; SOZ: 0 of 7
Baselines of Methods 4.8 (mixed model, cluster mass): band baseline z set block_mean_A_minus_N smallest_p sig_clusters delta 1s trial nSOZ -0.042 1.000 none delta 1s trial all14 -0.027 0.864 none delta 1s baseline nSOZ -1.389 0.010 A<N 4.3 to 7.9 s p=0.010; A<N 8.5 to 10.7 s p=0.028; A<N 11.6 to 15.4 s p=0.010; A<N 16.3 to 19.2 s p=0.014; A<N 21.4 to 23.2 s p=0.048 delta 1s baseline all14 -1.156 0.002 A<N 6.7 to 20.5 s p=0.002 delta 1s amygdala nSOZ -0.350 0.070 none delta 1s amygdala all14 -0.337 0.002 A<N 5.9 to 20.3 s p=0.002; A<N 21.3 to 24.0 s p=0.024 delta 1s trialsd nSOZ -0.498 0.058 none delta 1s trialsd all14 -0.445 0.004 A<N 5.8 to 19.3 s p=0.004; A<N 21.4 to 24.0 s p=0.024 delta 1s baselinepool nSOZ -0.768 0.004 A<N 3.1 to 5.3 s p=0.034; A<N 5.6 to 7.5 s p=0.032; A<N 8.7 to 10.8 s p=0.032; A<N 11.6 to 16.1 s p=0.004; A<N 21.4 to 24.0 s p=0.024 delta 1s baselinepool all14 -0.589 0.002 A<N 6.0 to 20.3 s p=0.002; A<N 21.3 to 23.8 s p=0.024 delta allepochs trial nSOZ -0.042 1.000 none delta allepochs trial all14 -0.027 0.850 none delta allepochs baseline nSOZ -1.389 0.008 A<N 4.3 to 7.9 s p=0.008; A<N 8.5 to 10.7 s p=0.034; A<N 11.6 to 15.4 s p=0.008; A<N 16.3 to 19.2 s p=0.010 delta allepochs baseline all14 -1.156 0.002 A<N 6.7 to 20.5 s p=0.002; A<N 21.4 to 23.2 s p=0.036 delta allepochs amygdala nSOZ -0.273 0.040 A<N 5.8 to 7.4 s p=0.040 delta allepochs amygdala all14 -0.206 0.020 A<N 9.1 to 11.0 s p=0.034; A<N 11.7 to 14.0 s p=0.020 delta allepochs trialsd nSOZ -0.321 0.018 A<N 5.7 to 7.4 s p=0.018 delta allepochs trialsd all14 -0.239 0.008 A<N 9.1 to 11.0 s p=0.044; A<N 11.3 to 14.0 s p=0.008; A<N 15.5 to 17.6 s p=0.032 delta allepochs baselinepool nSOZ -0.768 0.002 A<N 3.1 to 5.3 s p=0.044; A<N 5.6 to 7.5 s p=0.032; A<N 8.7 to 10.8 s p=0.038; A<N 11.6 to 16.1 s p=0.002; A<N 21.4 to 24.0 s p=0.014 delta allepochs baselinepool all14 -0.589 0.002 A<N 6.0 to 20.3 s p=0.002; A<N 21.3 to 23.8 s p=0.028 delta prevepoch trial nSOZ -0.042 1.000 none delta prevepoch trial all14 -0.027 0.886 none delta prevepoch baseline nSOZ -1.389 0.008 A<N 4.3 to 7.9 s p=0.008; A<N 8.5 to 10.7 s p=0.024; A<N 11.6 to 15.4 s p=0.008; A<N 16.3 to 19.2 s p=0.010; A<N 21.4 to 23.2 s p=0.046 delta prevepoch baseline all14 -1.156 0.002 A<N 6.7 to 20.5 s p=0.002; A<N 21.4 to 23.2 s p=0.048 delta prevepoch amygdala nSOZ -0.489 0.002 A<N 3.0 to 7.7 s p=0.002; A<N 8.4 to 19.3 s p=0.002; A<N 21.1 to 24.0 s p=0.004 delta prevepoch amygdala all14 -0.385 0.002 A<N 0.4 to 2.8 s p=0.020; A<N 5.7 to 20.4 s p=0.002; A<N 21.3 to 24.0 s p=0.002 delta prevepoch trialsd nSOZ -0.614 0.002 A<N 2.4 to 19.3 s p=0.002; A<N 20.8 to 24.0 s p=0.002 delta prevepoch trialsd all14 -0.495 0.002 A<N 0.3 to 20.6 s p=0.002; A<N 21.2 to 24.0 s p=0.002 delta prevepoch baselinepool nSOZ -0.768 0.002 A<N 3.1 to 5.3 s p=0.030; A<N 5.6 to 7.5 s p=0.026; A<N 8.7 to 10.8 s p=0.026; A<N 11.6 to 16.1 s p=0.002; A<N 21.4 to 24.0 s p=0.014 delta prevepoch baselinepool all14 -0.589 0.002 A<N 6.0 to 20.3 s p=0.002; A<N 21.3 to 23.8 s p=0.016 hgp 1s trial nSOZ 0.019 0.483 none hgp 1s trial all14 -0.004 0.218 none hgp 1s baseline nSOZ 0.532 0.082 none hgp 1s baseline all14 -0.309 0.741 none hgp 1s amygdala nSOZ 0.130 0.142 none hgp 1s amygdala all14 -0.104 0.477 none hgp 1s trialsd nSOZ 0.166 0.355 none hgp 1s trialsd all14 -0.136 0.599 none hgp 1s baselinepool nSOZ 0.442 0.140 none hgp 1s baselinepool all14 -0.149 0.685 none hgp allepochs trial nSOZ 0.019 0.483 none hgp allepochs trial all14 -0.004 0.228 none hgp allepochs baseline nSOZ 0.532 0.082 none hgp allepochs baseline all14 -0.309 0.723 none hgp allepochs amygdala nSOZ 0.278 0.004 A>N 3.2 to 4.6 s p=0.004; A>N 6.5 to 7.9 s p=0.018 hgp allepochs amygdala all14 0.130 0.010 A>N 3.3 to 4.6 s p=0.010 hgp allepochs trialsd nSOZ 0.317 0.006 A>N 0.6 to 1.9 s p=0.028; A>N 3.1 to 4.6 s p=0.006 hgp allepochs trialsd all14 0.116 0.060 none hgp allepochs baselinepool nSOZ 0.442 0.116 none hgp allepochs baselinepool all14 -0.149 0.713 none hgp prevepoch trial nSOZ 0.019 0.497 none hgp prevepoch trial all14 -0.004 0.190 none hgp prevepoch baseline nSOZ 0.532 0.054 none hgp prevepoch baseline all14 -0.309 0.737 none hgp prevepoch amygdala nSOZ 0.426 0.002 A>N 0.6 to 5.5 s p=0.002; A>N 6.1 to 8.1 s p=0.002; A>N 14.9 to 16.1 s p=0.014; A>N 17.0 to 18.5 s p=0.010 hgp prevepoch amygdala all14 0.191 0.002 A>N 2.2 to 4.7 s p=0.002; A>N 6.1 to 7.8 s p=0.036 hgp prevepoch trialsd nSOZ 0.645 0.002 A>N 0.5 to 8.4 s p=0.002; A>N 8.5 to 13.4 s p=0.002; A>N 14.4 to 18.7 s p=0.002; A>N 19.6 to 23.4 s p=0.002 hgp prevepoch trialsd all14 0.226 0.030 A>N 3.4 to 4.6 s p=0.030 hgp prevepoch baselinepool nSOZ 0.442 0.118 none hgp prevepoch baselinepool all14 -0.149 0.713 none
The per-amygdala tables show no amygdala whose HGP rises in the first 2 s under every z-scoring:
under baseline 4 of the 7 nSOZ amygdalae have a positive 0 to 2 s bin, under amygdala 5 and
under trial 4, against 0 to 4 of the 7 SOZ amygdalae. The delta difference is negative over the
block in 10 to 12 of 14 amygdalae. The baselines of Methods 4.8 leave trial and baseline
unchanged, as they must. Under amygdala the mean of all epochs gives delta clusters of 1.6 s
(nSOZ) and 4.2 s (all 14) and HGP clusters from 3.2 s; the previous epoch gives delta clusters over
most of the block in both sets and an HGP cluster in the nSOZ amygdalae from 0.6 to 5.5 s
(p = 0.002). With the trialsd variant two more HGP clusters begin in the first 2 s, both in the
nSOZ amygdalae: 0.6 to 1.9 s with the mean of all epochs (p = 0.028) and 0.5 to 8.4 s with the
previous epoch (p = 0.002). These three are the only readings here with an HGP cluster that begins
in the first 2 s, and all use a baseline other than the 1 s clip. As the blocks
alternate, the previous epoch is always a block of the other condition, so this baseline adds the
neighbouring block's condition difference to every trial; it is a variant of the claims file and
does not change the outcomes above, but the results are not "quantitatively similar" across these
baselines here.
other-bands-no-effect¶
No AVERSIVE against NEUTRAL difference is found in the theta, alpha, beta and gamma (30-60 Hz) bands.
A claim of absence, judged on the paper's test (mixed model, cluster permutation, p < 0.05), under the three z-scorings of the claims file and both amygdala sets; the two z-scoring variants are shown.
for band in ["theta", "alpha", "beta", "gamma"]:
G[band] = run_band(band, np.random.default_rng(103 + list(A.BANDS).index(band)))
ob = pd.concat([G[b] for b in ["theta", "alpha", "beta", "gamma"]])
print(ob[ob.model == "lmm"][SHOW].to_string(index=False))
lm = ob[(ob.model == "lmm") & (ob.stat == "mass")]
print("\nMixed model, cluster mass: bands with any significant cluster, by reading")
print(lm.assign(any_sig=(lm.n_sig_neg + lm.n_sig_pos) > 0).pivot_table(index="z", columns="set", values="any_sig", aggfunc="sum").to_string())
lmall = ob[ob.model == "lmm"].assign(any_sig=lambda v: (v.n_sig_neg + v.n_sig_pos) > 0)
dl_ = G["delta"][G["delta"].model == "lmm"].set_index(["z", "set", "stat"])
by_reading = lmall.groupby(["z", "set", "stat"]).any_sig.any().to_frame("other_band_cluster")
by_reading["delta_effect"] = [(dl_.loc[k].n_sig_neg > 0) and (dl_.loc[k].block_mean_A_minus_N < 0) for k in by_reading.index]
print("\nPer mixed-model reading: does another band have a significant cluster, and does delta show the effect of Fig. 1B?")
print(by_reading.to_string())
cl = by_reading.loc[ZS]
print("Readings of the claims file with no cluster in any other band:", int((~cl.other_band_cluster).sum()), "of", len(cl),
"; of those, readings with the delta effect:", int((~cl.other_band_cluster & cl.delta_effect).sum()))
print("Readings of the claims file with the delta effect:", int(cl.delta_effect.sum()), "; of those with another band clean:",
int((cl.delta_effect & ~cl.other_band_cluster).sum()))
band z reading set model stat block_mean_A_minus_N n_sig_neg n_sig_pos s_neg s_pos smallest_p sig_clusters
theta trial claims nSOZ lmm mass -0.022 0 0 0.0 0 0.395 none
theta trial claims nSOZ lmm size -0.022 0 0 0.0 0 0.521 none
theta trial claims all14 lmm mass -0.026 0 0 0.0 0 0.130 none
theta trial claims all14 lmm size -0.026 0 0 0.0 0 0.226 none
theta baseline claims nSOZ lmm mass -0.616 0 0 0.0 0 0.072 none
theta baseline claims nSOZ lmm size -0.616 0 0 0.0 0 0.090 none
theta baseline claims all14 lmm mass -0.908 2 0 9.4 0 0.002 A<N 10.9 to 17.4 s p=0.002; A<N 20.0 to 22.9 s p=0.012
theta baseline claims all14 lmm size -0.908 2 0 9.4 0 0.004 A<N 10.9 to 17.4 s p=0.004; A<N 20.0 to 22.9 s p=0.016
theta amygdala claims nSOZ lmm mass -0.171 0 0 0.0 0 0.112 none
theta amygdala claims nSOZ lmm size -0.171 0 0 0.0 0 0.106 none
theta amygdala claims all14 lmm mass -0.286 2 0 9.3 0 0.004 A<N 7.6 to 12.3 s p=0.010; A<N 12.7 to 17.3 s p=0.004
theta amygdala claims all14 lmm size -0.286 2 0 9.3 0 0.004 A<N 7.6 to 12.3 s p=0.004; A<N 12.7 to 17.3 s p=0.004
theta trialsd variant nSOZ lmm mass -0.216 0 0 0.0 0 0.102 none
theta trialsd variant nSOZ lmm size -0.216 0 0 0.0 0 0.098 none
theta trialsd variant all14 lmm mass -0.365 3 0 8.6 0 0.002 A<N 7.6 to 9.5 s p=0.040; A<N 10.2 to 12.4 s p=0.022; A<N 12.6 to 17.1 s p=0.002
theta trialsd variant all14 lmm size -0.365 3 0 8.6 0 0.004 A<N 7.6 to 9.5 s p=0.038; A<N 10.2 to 12.4 s p=0.026; A<N 12.6 to 17.1 s p=0.004
theta baselinepool variant nSOZ lmm mass -0.386 0 0 0.0 0 0.056 none
theta baselinepool variant nSOZ lmm size -0.386 0 0 0.0 0 0.068 none
theta baselinepool variant all14 lmm mass -0.488 2 0 6.9 0 0.004 A<N 7.6 to 10.0 s p=0.036; A<N 12.8 to 17.3 s p=0.004
theta baselinepool variant all14 lmm size -0.488 2 0 6.9 0 0.006 A<N 7.6 to 10.0 s p=0.032; A<N 12.8 to 17.3 s p=0.006
alpha trial claims nSOZ lmm mass -0.019 0 0 0.0 0 0.387 none
alpha trial claims nSOZ lmm size -0.019 0 0 0.0 0 0.529 none
alpha trial claims all14 lmm mass -0.022 0 0 0.0 0 0.315 none
alpha trial claims all14 lmm size -0.022 0 0 0.0 0 0.333 none
alpha baseline claims nSOZ lmm mass -0.550 1 0 2.6 0 0.016 A<N 13.3 to 15.9 s p=0.016
alpha baseline claims nSOZ lmm size -0.550 1 0 2.6 0 0.014 A<N 13.3 to 15.9 s p=0.014
alpha baseline claims all14 lmm mass -0.635 2 0 5.7 0 0.020 A<N 7.2 to 9.9 s p=0.028; A<N 13.2 to 16.2 s p=0.020
alpha baseline claims all14 lmm size -0.635 2 0 5.7 0 0.008 A<N 7.2 to 9.9 s p=0.014; A<N 13.2 to 16.2 s p=0.008
alpha amygdala claims nSOZ lmm mass -0.174 0 0 0.0 0 0.273 none
alpha amygdala claims nSOZ lmm size -0.174 0 0 0.0 0 0.240 none
alpha amygdala claims all14 lmm mass -0.265 3 0 12.0 0 0.004 A<N 6.7 to 12.1 s p=0.004; A<N 13.0 to 16.2 s p=0.004; A<N 20.6 to 24.0 s p=0.004
alpha amygdala claims all14 lmm size -0.265 3 0 12.0 0 0.002 A<N 6.7 to 12.1 s p=0.002; A<N 13.0 to 16.2 s p=0.012; A<N 20.6 to 24.0 s p=0.010
alpha trialsd variant nSOZ lmm mass -0.309 0 0 0.0 0 0.150 none
alpha trialsd variant nSOZ lmm size -0.309 0 0 0.0 0 0.120 none
alpha trialsd variant all14 lmm mass -0.409 3 0 15.0 0 0.004 A<N 6.7 to 12.1 s p=0.004; A<N 12.8 to 18.2 s p=0.004; A<N 19.8 to 24.0 s p=0.004
alpha trialsd variant all14 lmm size -0.409 3 0 15.0 0 0.002 A<N 6.7 to 12.1 s p=0.002; A<N 12.8 to 18.2 s p=0.002; A<N 19.8 to 24.0 s p=0.002
alpha baselinepool variant nSOZ lmm mass -0.348 0 0 0.0 0 0.188 none
alpha baselinepool variant nSOZ lmm size -0.348 0 0 0.0 0 0.140 none
alpha baselinepool variant all14 lmm mass -0.409 2 0 5.9 0 0.012 A<N 7.2 to 10.1 s p=0.014; A<N 13.1 to 16.1 s p=0.012
alpha baselinepool variant all14 lmm size -0.409 2 0 5.9 0 0.010 A<N 7.2 to 10.1 s p=0.010; A<N 13.1 to 16.1 s p=0.010
beta trial claims nSOZ lmm mass 0.004 0 0 0.0 0 0.968 none
beta trial claims nSOZ lmm size 0.004 0 0 0.0 0 0.978 none
beta trial claims all14 lmm mass -0.012 0 0 0.0 0 0.160 none
beta trial claims all14 lmm size -0.012 0 0 0.0 0 0.168 none
beta baseline claims nSOZ lmm mass 0.219 0 0 0.0 0 0.691 none
beta baseline claims nSOZ lmm size 0.219 0 0 0.0 0 0.731 none
beta baseline claims all14 lmm mass -0.541 0 0 0.0 0 0.192 none
beta baseline claims all14 lmm size -0.541 0 0 0.0 0 0.160 none
beta amygdala claims nSOZ lmm mass -0.038 0 0 0.0 0 1.000 none
beta amygdala claims nSOZ lmm size -0.038 0 0 0.0 0 1.000 none
beta amygdala claims all14 lmm mass -0.184 1 0 4.7 0 0.004 A<N 11.2 to 15.9 s p=0.004
beta amygdala claims all14 lmm size -0.184 1 0 4.7 0 0.002 A<N 11.2 to 15.9 s p=0.002
beta trialsd variant nSOZ lmm mass -0.073 0 0 0.0 0 1.000 none
beta trialsd variant nSOZ lmm size -0.073 0 0 0.0 0 1.000 none
beta trialsd variant all14 lmm mass -0.312 1 0 5.8 0 0.002 A<N 11.1 to 16.9 s p=0.002
beta trialsd variant all14 lmm size -0.312 1 0 5.8 0 0.002 A<N 11.1 to 16.9 s p=0.002
beta baselinepool variant nSOZ lmm mass 0.113 0 0 0.0 0 0.872 none
beta baselinepool variant nSOZ lmm size 0.113 0 0 0.0 0 0.920 none
beta baselinepool variant all14 lmm mass -0.234 0 0 0.0 0 0.164 none
beta baselinepool variant all14 lmm size -0.234 0 0 0.0 0 0.160 none
gamma trial claims nSOZ lmm mass -0.000 0 0 0.0 0 0.445 none
gamma trial claims nSOZ lmm size -0.000 0 0 0.0 0 0.439 none
gamma trial claims all14 lmm mass -0.012 0 0 0.0 0 0.176 none
gamma trial claims all14 lmm size -0.012 0 0 0.0 0 0.269 none
gamma baseline claims nSOZ lmm mass -0.145 0 0 0.0 0 0.880 none
gamma baseline claims nSOZ lmm size -0.145 0 0 0.0 0 0.902 none
gamma baseline claims all14 lmm mass -0.308 0 0 0.0 0 0.393 none
gamma baseline claims all14 lmm size -0.308 0 0 0.0 0 0.407 none
gamma amygdala claims nSOZ lmm mass -0.115 0 0 0.0 0 0.393 none
gamma amygdala claims nSOZ lmm size -0.115 0 0 0.0 0 0.361 none
gamma amygdala claims all14 lmm mass -0.219 2 0 3.7 0 0.024 A<N 9.7 to 11.6 s p=0.024; A<N 13.2 to 15.0 s p=0.040
gamma amygdala claims all14 lmm size -0.219 2 0 3.7 0 0.042 A<N 9.7 to 11.6 s p=0.042; A<N 13.2 to 15.0 s p=0.050
gamma trialsd variant nSOZ lmm mass -0.122 0 0 0.0 0 0.427 none
gamma trialsd variant nSOZ lmm size -0.122 0 0 0.0 0 0.463 none
gamma trialsd variant all14 lmm mass -0.275 1 0 1.9 0 0.024 A<N 9.7 to 11.6 s p=0.024
gamma trialsd variant all14 lmm size -0.275 1 0 1.9 0 0.044 A<N 9.7 to 11.6 s p=0.044
gamma baselinepool variant nSOZ lmm mass -0.028 0 0 0.0 0 0.760 none
gamma baselinepool variant nSOZ lmm size -0.028 0 0 0.0 0 0.822 none
gamma baselinepool variant all14 lmm mass -0.328 0 0 0.0 0 0.110 none
gamma baselinepool variant all14 lmm size -0.328 0 0 0.0 0 0.148 none
Mixed model, cluster mass: bands with any significant cluster, by reading
set all14 nSOZ
z
amygdala 4 0
baseline 2 1
baselinepool 2 0
trial 0 0
trialsd 4 0
Per mixed-model reading: does another band have a significant cluster, and does delta show the effect of Fig. 1B?
other_band_cluster delta_effect
z set stat
amygdala all14 mass True True
size True True
nSOZ mass False False
size False False
baseline all14 mass True True
size True True
nSOZ mass True True
size True True
baselinepool all14 mass True True
size True True
nSOZ mass False True
size False True
trial all14 mass False False
size False False
nSOZ mass False False
size False False
trialsd all14 mass True True
size True True
nSOZ mass False True
size False True
Readings of the claims file with no cluster in any other band: 6 of 12 ; of those, readings with the delta effect: 0
Readings of the claims file with the delta effect: 6 ; of those with another band clean: 0
Result: not reproduced (resting on an assumption). The readings of the claims file disagree:
under trial z-scoring no band has a cluster, while under baseline the mixed model finds clusters
with AVERSIVE below NEUTRAL in alpha in the 7 nSOZ amygdalae (13.3 to 15.9 s, p 0.016) and in theta
and alpha in all 14 (p 0.002 to 0.028), and under amygdala all four bands have clusters in all 14
amygdalae (none in the 7 nSOZ). The outcome rests on one assumption: the paper's pipeline is one under
which its own delta result of Fig. 1B appears, since the paper reports both from the same analysis.
Every reading of the claims file that shows the delta effect (baseline in both sets, amygdala in
all 14; 6 of 12) also shows a cluster in another band, and the trial reading, the only one with no
other-band cluster, shows no delta effect either. Of the variants, baselinepool in the 7 nSOZ
amygdalae shows the delta effect and no other-band cluster (the alpha cluster of baseline does not
survive a pooled baseline SD); in all 14 it has theta and alpha clusters.
Modulation indices (table2-mi-values and the claims built on them)¶
MI = AVERSIVE minus NEUTRAL time-averaged power over the window where the conditions differ
(Methods 4.11, Table 2 legend), per amygdala. Window readings: (a) the group clusters of the band
(Fig. 1B and 1C) under the same z-scoring, from the mixed model with cluster mass over the 7 nSOZ
amygdalae (a_nSOZ) or all 14 (a_all14); (b) clusters found per amygdala (unpaired t-test of its
8 AVERSIVE against 9 NEUTRAL trials, cluster mass, 500 permutations); and (c), a reading the
Results' text gives, fixed windows: "the temporal window of task-dependent modulation, which
occurred at stimulus onset in HGP and was sustained throughout the trial in DELTA" (Results 2.1,
paragraph 2) and "HGP enhancement in the first 2 s or sustained delta decrease" (Methods 4.11), read
as 0 to 2 s for HGP and 0 to 24 s for delta. (a) and (b) are the claims file's readings; (c) is taken
from the paper's text as step 1.6 of the protocol allows, and is labelled so in every table.
Where a reading gives no significant window (no group cluster under (a), no cluster in that amygdala under (b)), the MI is left undefined: the paper defines no MI without a window, and Table 2 prints a non-zero MI for each of the 14 amygdalae. A claim is judged only on readings that define every MI it needs; the number of amygdalae each reading defines is printed.
rng_mi = np.random.default_rng(201)
def group_mask(band, z, set_name):
g = G[band]
r = g[(g.z == z) & (g.set == set_name) & (g.model == "lmm") & (g.stat == "mass")].iloc[0]
m = np.zeros(block.sum(), bool)
for c in r["clusters"]:
if c["p"] < 0.05:
m[c["start"]:c["stop"]] = True
return m
per_amy_cache = {}
def per_amygdala_mask(n, band, z):
key = (n, band, z)
if key not in per_amy_cache:
y = A.normalise(ENV[(n, band)], z)[:, block]
x = aversive_of[A.TABLE2[n - 1][0]].astype(int)
_, _, res = A.cluster_test(y, x, np.zeros(len(x), int), rng_mi, model="ttest", stat="mass")
m = np.zeros(block.sum(), bool)
for c in res:
if c["p"] < 0.05:
m[c["start"]:c["stop"]] = True
per_amy_cache[key] = m
return per_amy_cache[key]
def mi(n, band, z, window):
y = A.normalise(ENV[(n, band)], z)[:, block]
x = aversive_of[A.TABLE2[n - 1][0]]
if window.startswith("a_"):
m = group_mask(band, z, "nSOZ" if window == "a_nSOZ" else "all14")
elif window == "b":
m = per_amygdala_mask(n, band, z)
else:
m = tb < 2 if band == "hgp" else np.ones(block.sum(), bool)
if not m.any():
return np.nan
return y[x][:, m].mean() - y[~x][:, m].mean()
WINDOWS = ["a_nSOZ", "a_all14", "b", "c"]
MI = pd.DataFrame([{"amygdala": n, "z": z, "window": w, "delta": mi(n, "delta", z, w), "hgp": mi(n, "hgp", z, w)}
for n in ALL for z in ZALL for w in WINDOWS])
exp = claims["table2-mi-values"]["expected"]
printed = pd.DataFrame({"amygdala": ALL, "delta_paper": [e["value"] for e in exp[0::2]],
"hgp_paper": [e["value"] for e in exp[1::2]],
"delta_tol": [e["tolerance"] for e in exp[0::2]], "hgp_tol": [e["tolerance"] for e in exp[1::2]]})
MI = MI.merge(printed, on="amygdala")
MI["soz"] = MI.amygdala.isin(A.TABLE2_SOZ)
MI["delta_ok"] = (MI.delta - MI.delta_paper).abs() <= MI.delta_tol
MI["hgp_ok"] = (MI.hgp - MI.hgp_paper).abs() <= MI.hgp_tol
MI["delta_sign_ok"] = np.sign(MI.delta) == np.sign(MI.delta_paper)
MI["hgp_sign_ok"] = np.sign(MI.hgp) == np.sign(MI.hgp_paper)
summary = MI.groupby(["z", "window"]).agg(within_tol=("delta_ok", "sum"), hgp_within=("hgp_ok", "sum"),
delta_sign=("delta_sign_ok", "sum"), hgp_sign=("hgp_sign_ok", "sum"),
delta_defined=("delta", lambda v: int(v.notna().sum())),
hgp_defined=("hgp", lambda v: int(v.notna().sum())))
summary["values_within_tol_of_28"] = summary.within_tol + summary.hgp_within
summary["signs_agree_of_28"] = summary.delta_sign + summary.hgp_sign
summary["reading"] = ["claims" if z in ZS and w != "c" else "claims z, window from text" if z in ZS else "variant"
for z, w in summary.index]
print(summary[["reading", "values_within_tol_of_28", "signs_agree_of_28", "delta_defined", "hgp_defined"]].to_string())
print("\nPer amygdala under the fixed windows (c), every z-scoring, with Table 2:")
wide = MI[MI.window == "c"].pivot(index="amygdala", columns="z", values=["delta", "hgp"]).round(2)
wide["delta_paper"] = printed.set_index("amygdala").delta_paper
wide["hgp_paper"] = printed.set_index("amygdala").hgp_paper
print(wide.to_string())
print("\nPer amygdala under (a), group windows (NaN: no group cluster):")
print(MI[MI.window.str.startswith("a_")].pivot(index="amygdala", columns=["window", "z"], values="delta").round(2).to_string())
print("\nPer amygdala under (b), per-amygdala clusters (NaN: no cluster in that amygdala):")
print(MI[MI.window == "b"].pivot(index="amygdala", columns="z", values=["delta", "hgp"]).round(2).to_string())
reading values_within_tol_of_28 signs_agree_of_28 delta_defined hgp_defined
z window
amygdala a_all14 claims 1 12 14 0
a_nSOZ claims 0 0 0 0
b claims 0 4 1 3
c claims z, window from text 4 23 14 14
baseline a_all14 claims 0 11 14 0
a_nSOZ claims 1 11 14 0
b claims 0 3 2 2
c claims z, window from text 0 19 14 14
baselinepool a_all14 variant 1 11 14 0
a_nSOZ variant 0 10 14 0
b variant 0 5 2 4
c variant 4 19 14 14
trial a_all14 claims 0 0 0 0
a_nSOZ claims 0 0 0 0
b claims 0 1 2 2
c claims z, window from text 1 16 14 14
trialsd a_all14 variant 2 11 14 0
a_nSOZ variant 3 11 14 0
b variant 1 5 3 2
c variant 6 22 14 14
Per amygdala under the fixed windows (c), every z-scoring, with Table 2:
delta hgp delta_paper hgp_paper
z amygdala baseline baselinepool trial trialsd amygdala baseline baselinepool trial trialsd
amygdala
1 -0.06 -1.02 -0.78 -0.04 -0.01 0.60 2.48 1.27 -0.07 0.75 -0.66 0.79
2 0.12 1.23 0.57 0.04 0.13 -0.39 -0.44 -0.57 0.03 -0.56 0.26 -0.57
3 -0.16 -1.75 -0.71 -0.05 -0.24 0.33 0.44 0.76 0.34 0.21 -0.67 0.35
4 -0.17 -0.83 -0.54 -0.03 -0.22 -0.57 -0.14 -0.37 0.13 -0.60 0.35 -0.66
5 -0.58 -1.44 -0.53 -0.05 -0.96 1.66 4.00 3.17 0.82 2.08 -0.91 1.75
6 -0.83 -2.19 -2.11 -0.11 -1.29 0.08 1.15 0.22 0.02 0.35 -1.22 0.50
7 -0.62 -4.02 -0.88 -0.04 -0.88 -0.38 -4.80 -1.89 0.04 -0.46 -0.33 -0.59
8 -0.45 -3.14 -1.10 -0.07 -0.69 0.01 -0.07 -0.23 0.19 0.06 -1.13 1.07
9 -0.10 -0.07 0.07 0.00 -0.13 -0.55 0.23 -0.51 -0.28 -1.02 -0.09 -0.52
10 0.37 1.17 0.92 0.06 0.75 -0.54 -0.85 -0.31 -0.17 -0.72 -0.39 0.04
11 -0.26 0.22 0.14 0.01 -0.33 -0.34 -0.87 -1.02 -0.66 -0.39 -0.17 0.25
12 -0.74 -1.36 -1.07 -0.04 -1.05 -0.78 -1.08 -1.13 -0.39 -0.97 -0.59 0.05
13 -0.62 -1.72 -1.28 -0.05 -0.63 -0.27 -0.78 -0.85 -0.22 -0.32 -0.82 -0.44
14 -0.62 -1.27 -0.95 -0.04 -0.70 -0.02 0.18 0.20 0.56 0.06 -0.68 -0.30
Per amygdala under (a), group windows (NaN: no group cluster):
window a_nSOZ a_all14 a_nSOZ a_all14 a_nSOZ a_all14 a_nSOZ a_all14 a_nSOZ a_all14
z trial trial baseline baseline amygdala amygdala trialsd trialsd baselinepool baselinepool
amygdala
1 NaN NaN -1.40 -1.28 NaN -0.12 -0.02 -0.07 -1.05 -0.89
2 NaN NaN 1.03 1.01 NaN 0.03 -0.07 -0.03 0.42 0.44
3 NaN NaN -2.17 -1.97 NaN -0.20 -0.63 -0.33 -0.75 -0.77
4 NaN NaN -1.01 -1.14 NaN -0.26 -0.35 -0.32 -0.51 -0.72
5 NaN NaN -1.74 -1.38 NaN -0.58 -0.79 -0.96 -0.55 -0.52
6 NaN NaN -2.35 -2.21 NaN -0.83 -1.43 -1.32 -2.24 -2.12
7 NaN NaN -3.93 -5.65 NaN -0.74 -0.78 -1.04 -0.65 -1.19
8 NaN NaN -3.47 -3.23 NaN -0.48 -1.10 -0.76 -1.52 -1.14
9 NaN NaN -0.10 -0.29 NaN -0.18 -0.29 -0.20 0.09 -0.07
10 NaN NaN 1.12 0.91 NaN 0.32 0.62 0.70 0.84 0.81
11 NaN NaN 0.21 0.23 NaN -0.22 -0.22 -0.25 0.19 0.21
12 NaN NaN -1.30 -1.45 NaN -0.77 -1.05 -1.06 -0.91 -1.15
13 NaN NaN -1.60 -1.64 NaN -0.62 -0.64 -0.60 -1.22 -1.28
14 NaN NaN -1.24 -1.58 NaN -0.69 -0.74 -0.75 -0.91 -1.09
Per amygdala under (b), per-amygdala clusters (NaN: no cluster in that amygdala):
delta hgp
z amygdala baseline baselinepool trial trialsd amygdala baseline baselinepool trial trialsd
amygdala
1 NaN NaN NaN NaN NaN 1.44 NaN 2.95 NaN 1.73
2 NaN NaN NaN NaN NaN NaN NaN NaN 1.10 NaN
3 NaN NaN NaN 1.35 NaN 1.31 2.66 2.47 1.38 NaN
4 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
5 NaN NaN NaN NaN NaN 1.34 2.84 2.29 NaN 1.86
6 -1.03 -2.49 -2.53 NaN -1.68 NaN NaN NaN NaN NaN
7 NaN NaN NaN NaN -1.44 NaN NaN NaN NaN NaN
8 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
9 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
10 NaN 2.54 1.61 NaN NaN NaN NaN NaN NaN NaN
11 NaN NaN NaN 1.10 NaN NaN NaN NaN NaN NaN
12 NaN NaN NaN NaN -1.83 NaN NaN NaN NaN NaN
13 NaN NaN NaN NaN NaN NaN NaN -1.75 NaN NaN
14 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
Result: not reproduced. No reading gives all 28 values within tolerance: at best 4 of 28 among
the claims file's z-scorings (amygdala, window c; 6 with the trialsd variant). No reading gives
the 28 signs of Table 2: at best 23 of 28 (amygdala, window c). Window (a) defines no HGP MI, as
there is no HGP group cluster, and window (b) defines only 1 to 3 MIs per band.
Does the dandiset's patient numbering match the paper's?¶
Two hints that it may not: the units are in subjects 1, 2, 3, 4 and 7 while Table 1 lists MUA for 1, 2, 3, 4, 5 R and 8 L, and subject 2's file gives age 39, the age Table 1 gives for patient 7 (both are gliosis cases with AL, AR implanted and AR in the SOZ, so they could be exchanged without the SOZ flags showing it). The MIs are compared with Table 2 under the file numbering and under two re-mappings: subjects 2 and 7 swapped (amygdalae 3 and 4 against 10 and 11), and subject 7's amygdalae taken as patients 8 L and 5 R (10 against 12, 11 against 7).
MAPS = {"file numbering": {}, "subjects 2 and 7 swapped": {3: 10, 4: 11, 10: 3, 11: 4},
"subject 7 as 8 L and 5 R": {10: 12, 12: 10, 11: 7, 7: 11}}
rows = []
for name, sw in MAPS.items():
target = [sw.get(n, n) for n in ALL]
pdl = printed.set_index("amygdala").loc[target]
for z in ZALL:
m = MI[(MI.z == z) & (MI.window == "c")].sort_values("amygdala")
dl, hg = m.delta.values, m.hgp.values
rows.append({"mapping": name, "z": z, "window": "c",
"delta_sign_agree": int((np.sign(dl) == np.sign(pdl.delta_paper.values)).sum()),
"hgp_sign_agree": int((np.sign(hg) == np.sign(pdl.hgp_paper.values)).sum()),
"within_tol_of_28": int((np.abs(dl - pdl.delta_paper.values) <= pdl.delta_tol.values).sum()
+ (np.abs(hg - pdl.hgp_paper.values) <= pdl.hgp_tol.values).sum()),
"spearman_delta_vs_table2": round(stats.spearmanr(dl, pdl.delta_paper.values).statistic, 2),
"spearman_hgp_vs_table2": round(stats.spearmanr(hg, pdl.hgp_paper.values).statistic, 2)})
print(pd.DataFrame(rows).to_string(index=False))
mapping z window delta_sign_agree hgp_sign_agree within_tol_of_28 spearman_delta_vs_table2 spearman_hgp_vs_table2
file numbering trial c 10 6 1 0.85 0.19
file numbering baseline c 11 8 0 0.67 0.56
file numbering amygdala c 12 11 4 0.56 0.75
file numbering trialsd c 12 10 6 0.60 0.73
file numbering baselinepool c 10 9 4 0.69 0.66
subjects 2 and 7 swapped trial c 12 8 1 0.75 0.41
subjects 2 and 7 swapped baseline c 13 8 0 0.58 0.61
subjects 2 and 7 swapped amygdala c 12 11 4 0.49 0.46
subjects 2 and 7 swapped trialsd c 12 10 5 0.51 0.52
subjects 2 and 7 swapped baselinepool c 12 9 3 0.65 0.75
subject 7 as 8 L and 5 R trial c 10 8 1 0.79 0.45
subject 7 as 8 L and 5 R baseline c 11 8 0 0.60 0.51
subject 7 as 8 L and 5 R amygdala c 12 11 4 0.49 0.73
subject 7 as 8 L and 5 R trialsd c 12 10 7 0.53 0.70
subject 7 as 8 L and 5 R baselinepool c 10 9 4 0.62 0.62
Neither re-mapping brings the MIs clearly closer to Table 2. Under window c and the claims file's z-scorings the signs agree in 16, 19 and 23 of 28 under the file numbering, 20, 21 and 23 with subjects 2 and 7 swapped, and 18, 19 and 23 with subject 7 as patients 8 L and 5 R; the values within tolerance stay at 0 to 4 of 28. The analysis keeps the file numbering, which the SOZ flags, sexes, pathologies and implanted sides of the files support (Table 1).
Claims on the MI values¶
Each claim below is computed under every z-scoring and window. A row counts towards the outcome when
its z-scoring is one of the claims file's (trial, baseline, amygdala) and it defines every MI the
claim needs; rows of the variant z-scorings and rows with undefined MIs are shown and marked. Window
(a) defines no HGP MI (there is no HGP group cluster), and window (b) defines only the 1 to 3
amygdalae per band that have a cluster of their own, so most claims rest on window (c). The printed
Table 2 values are run through the same computations as a reference.
soz_mask = np.array([n in A.TABLE2_SOZ for n in ALL])
paper_delta, paper_hgp = printed.delta_paper.values, printed.hgp_paper.values
def readings(need=("delta", "hgp"), amys=None):
"""(z, window, kind, delta MI, HGP MI, defined) for every reading; `defined` when every MI in
`need` is defined for the amygdalae `amys` (boolean mask over the 14, default all)."""
sel = np.ones(14, bool) if amys is None else amys
for z in ZALL:
for w in WINDOWS:
m = MI[(MI.z == z) & (MI.window == w)].sort_values("amygdala")
dl, hg = m.delta.values, m.hgp.values
ok = all(np.isfinite({"delta": dl, "hgp": hg}[b][sel]).all() for b in need)
kind = "variant" if z in ZV else "claims"
yield z, w, kind, dl, hg, ok
def counted(df):
"""Rows that decide the outcome: claims-file z-scoring and every needed MI defined."""
return df[(df.kind == "claims") & df.defined]
fig1d-delta-hgp-mi-correlation¶
Across the 14 amygdalae, the DELTA modulation index and the HGP modulation index are negatively correlated.
Criterion: Spearman r within 0.067 of -0.67 and p < 0.05 under every reading; partially reproduced when r is negative with p < 0.05 but outside tolerance.
rows = []
for z, w, kind, dl, hg, ok in readings():
r = stats.spearmanr(dl, hg) if ok else (np.nan, np.nan)
rows.append({"z": z, "window": w, "kind": kind, "defined": ok, "r": r[0], "p": r[1]})
rows.append({"z": "Table 2 printed", "window": "", "kind": "paper", "defined": True,
**dict(zip(["r", "p"], stats.spearmanr(paper_delta, paper_hgp)))})
r1d = pd.DataFrame(rows)
r1d["within_tol"] = (r1d.r + 0.67).abs() <= 0.067
print(r1d.round(3).to_string(index=False))
c = counted(r1d)
print(f"Counted readings: {len(c)}; r from {c.r.min():.2f} to {c.r.max():.2f}; within tolerance with p < 0.05: "
f"{int((c.within_tol & (c.p < 0.05)).sum())}; negative with p < 0.05: {int(((c.r < 0) & (c.p < 0.05)).sum())}")
z window kind defined r p within_tol
trial a_nSOZ claims False NaN NaN False
trial a_all14 claims False NaN NaN False
trial b claims False NaN NaN False
trial c claims True -0.301 0.296 False
baseline a_nSOZ claims False NaN NaN False
baseline a_all14 claims False NaN NaN False
baseline b claims False NaN NaN False
baseline c claims True -0.116 0.692 False
amygdala a_nSOZ claims False NaN NaN False
amygdala a_all14 claims False NaN NaN False
amygdala b claims False NaN NaN False
amygdala c claims True -0.121 0.681 False
trialsd a_nSOZ variant False NaN NaN False
trialsd a_all14 variant False NaN NaN False
trialsd b variant False NaN NaN False
trialsd c variant True -0.332 0.246 False
baselinepool a_nSOZ variant False NaN NaN False
baselinepool a_all14 variant False NaN NaN False
baselinepool b variant False NaN NaN False
baselinepool c variant True -0.037 0.899 False
Table 2 printed paper True -0.736 0.003 True
Counted readings: 3; r from -0.30 to -0.12; within tolerance with p < 0.05: 0; negative with p < 0.05: 0
Result: not reproduced. The three readings that define both MIs (window c under trial,
baseline and amygdala) give r = -0.30, -0.12 and -0.12 (p 0.30, 0.69, 0.68), negative but far
from -0.67 and not significant; the variants give -0.33 and -0.04. The printed Table 2 values give
r = -0.74 (p = 0.003).
This outcome rests on an assumption: window (c), fixed windows taken from the Results' text, is the only window that defines the MIs it needs, while the paper defines the MI over the interval where the two curves are significantly separated (Table 2 legend). The question is for the authors (methods_unavailable).
nsoz-delta-hgp-mi-inverse¶
Among the 7 amygdalae outside the SOZ, the stronger the HGP enhancement (HGP MI), the stronger the reduction of DELTA power (more negative DELTA MI).
Criterion: Spearman r between DELTA MI and HGP MI over the 7 nSOZ amygdalae negative under every reading.
rows = []
for z, w, kind, dl, hg, ok in readings(amys=~soz_mask):
r = stats.spearmanr(dl[~soz_mask], hg[~soz_mask]) if ok else (np.nan, np.nan)
rows.append({"z": z, "window": w, "kind": kind, "defined": ok, "r": r[0], "p": r[1]})
rows.append({"z": "Table 2 printed", "window": "", "kind": "paper", "defined": True,
**dict(zip(["r", "p"], stats.spearmanr(paper_delta[~soz_mask], paper_hgp[~soz_mask])))})
rn = pd.DataFrame(rows)
print(rn.round(3).to_string(index=False))
c = counted(rn)
print("Counted readings with r < 0:", int((c.r < 0).sum()), "of", len(c))
z window kind defined r p
trial a_nSOZ claims False NaN NaN
trial a_all14 claims False NaN NaN
trial b claims False NaN NaN
trial c claims True -0.571 0.180
baseline a_nSOZ claims False NaN NaN
baseline a_all14 claims False NaN NaN
baseline b claims False NaN NaN
baseline c claims True -0.143 0.760
amygdala a_nSOZ claims False NaN NaN
amygdala a_all14 claims False NaN NaN
amygdala b claims False NaN NaN
amygdala c claims True 0.071 0.879
trialsd a_nSOZ variant False NaN NaN
trialsd a_all14 variant False NaN NaN
trialsd b variant False NaN NaN
trialsd c variant True -0.071 0.879
baselinepool a_nSOZ variant False NaN NaN
baselinepool a_all14 variant False NaN NaN
baselinepool b variant False NaN NaN
baselinepool c variant True 0.214 0.645
Table 2 printed paper True -0.679 0.094
Counted readings with r < 0: 2 of 3
Result: not reproduced. Over the 7 nSOZ amygdalae, window c gives r = -0.57 (trial), -0.14
(baseline) and +0.07 (amygdala): negative in 2 of the 3 readings that define the MIs, so not
under every reading. Windows (a) and (b) define no HGP MI for these amygdalae. The variants give
-0.07 and +0.21. The printed Table 2 values give r = -0.68 (p = 0.094).
This outcome rests on an assumption: window (c), fixed windows taken from the Results' text, is the only window that defines the MIs it needs, while the paper defines the MI over the interval where the two curves are significantly separated (Table 2 legend). The question is for the authors (methods_unavailable).
table2-delta-reduced-12-of-14¶
DELTA power is reduced (negative MI) in 12 of the 14 amygdalae, whether inside or outside the SOZ.
Criterion: 11 to 13 negative DELTA MI under every reading, with negative MI both inside and outside the SOZ.
rows = []
for z, w, kind, dl, hg, ok in readings(need=("delta",)):
d_ = dl if ok else np.full(14, np.nan)
rows.append({"z": z, "window": w, "kind": kind, "defined": ok, "delta_defined": int(np.isfinite(dl).sum()),
"negative": int((d_ < 0).sum()) if ok else np.nan,
"negative_nSOZ": int((d_[~soz_mask] < 0).sum()) if ok else np.nan,
"negative_SOZ": int((d_[soz_mask] < 0).sum()) if ok else np.nan,
"positive_amygdalae": [n for n, v in zip(ALL, dl) if v > 0]})
rows.append({"z": "Table 2 printed", "window": "", "kind": "paper", "defined": True, "delta_defined": 14,
"negative": int((paper_delta < 0).sum()),
"negative_nSOZ": int((paper_delta[~soz_mask] < 0).sum()), "negative_SOZ": int((paper_delta[soz_mask] < 0).sum()),
"positive_amygdalae": [n for n, v in zip(ALL, paper_delta) if v > 0]})
r12 = pd.DataFrame(rows)
r12["holds"] = r12.negative.between(11, 13) & (r12.negative_nSOZ > 0) & (r12.negative_SOZ > 0)
print(r12.to_string(index=False))
c = counted(r12)
print("Counted readings where it holds:", int(c.holds.sum()), "of", len(c), "; failing:",
c[~c.holds][["z", "window", "negative", "negative_nSOZ", "negative_SOZ"]].to_dict("records"))
z window kind defined delta_defined negative negative_nSOZ negative_SOZ positive_amygdalae holds
trial a_nSOZ claims False 0 NaN NaN NaN [] False
trial a_all14 claims False 0 NaN NaN NaN [] False
trial b claims False 2 NaN NaN NaN [3, 11] False
trial c claims True 14 10.0 6.0 4.0 [2, 9, 10, 11] False
baseline a_nSOZ claims True 14 11.0 6.0 5.0 [2, 10, 11] True
baseline a_all14 claims True 14 11.0 6.0 5.0 [2, 10, 11] True
baseline b claims False 2 NaN NaN NaN [10] False
baseline c claims True 14 11.0 6.0 5.0 [2, 10, 11] True
amygdala a_nSOZ claims False 0 NaN NaN NaN [] False
amygdala a_all14 claims True 14 12.0 6.0 6.0 [2, 10] True
amygdala b claims False 1 NaN NaN NaN [] False
amygdala c claims True 14 12.0 6.0 6.0 [2, 10] True
trialsd a_nSOZ variant True 14 13.0 6.0 7.0 [10] True
trialsd a_all14 variant True 14 13.0 6.0 7.0 [10] True
trialsd b variant False 3 NaN NaN NaN [] False
trialsd c variant True 14 12.0 6.0 6.0 [2, 10] True
baselinepool a_nSOZ variant True 14 10.0 6.0 4.0 [2, 9, 10, 11] False
baselinepool a_all14 variant True 14 11.0 6.0 5.0 [2, 10, 11] True
baselinepool b variant False 2 NaN NaN NaN [10] False
baselinepool c variant True 14 10.0 6.0 4.0 [2, 9, 10, 11] False
Table 2 printed paper True 14 12.0 7.0 5.0 [2, 4] True
Counted readings where it holds: 5 of 6 ; failing: [{'z': 'trial', 'window': 'c', 'negative': 10.0, 'negative_nSOZ': 6.0, 'negative_SOZ': 4.0}]
Result: partially reproduced. In 5 of the 6 readings that define all 14 delta MIs the count is
11 or 12, with negative MIs both outside (6 of 7) and inside (5 or 6 of 7) the SOZ: baseline with
windows a (both sets) and c (11; positive in amygdalae 2, 10 and 11), and amygdala with windows
a (all 14) and c (12; positive in 2 and 10). Under trial with window c the count is 10 (positive
in 2, 9, 10 and 11), one below the range, so it does not hold under every reading. Table 2 has
amygdalae 2 and 4 positive. Window (b) defines 1 or 2 delta MIs. The variants give 10 to 13.
This outcome rests on an assumption: its only failing reading uses window (c), fixed windows taken
from the Results' text and not from the claims file. The readings of the claims file that define all
14 delta MIs (window a under baseline in both sets and under amygdala in all 14) give 11, 11
and 12, all in range. The question is for the authors (methods_unavailable).
table2-hgp-mi-positive-nsoz¶
All 7 amygdalae outside the SOZ have a positive HGP MI, while the response is absent in amygdalae within the SOZ (6 of 7 negative).
Criterion: both counts match (7 of 7, 6 of 7) and the difference has p < 0.05 under one of the tests: (a) Wilcoxon rank-sum, (b) unpaired t-test, (c) Fisher's exact test of the counts of positive MI (chi-square shown beside it).
rows = []
for z, w, kind, dl, hg, ok in list(readings(need=("hgp",))) + [("Table 2 printed", "", "paper", paper_delta, paper_hgp, True)]:
if not ok:
rows.append({"z": z, "window": w, "kind": kind, "defined": ok, "hgp_defined": int(np.isfinite(hg).sum())})
continue
pos_n, neg_s = int((hg[~soz_mask] > 0).sum()), int((hg[soz_mask] < 0).sum())
table = [[pos_n, 7 - pos_n], [int((hg[soz_mask] > 0).sum()), 7 - int((hg[soz_mask] > 0).sum())]]
rows.append({"z": z, "window": w, "kind": kind, "defined": ok, "hgp_defined": 14, "nSOZ_positive": pos_n, "SOZ_negative": neg_s,
"p_ranksum": stats.ranksums(hg[~soz_mask], hg[soz_mask]).pvalue,
"p_ttest": stats.ttest_ind(hg[~soz_mask], hg[soz_mask]).pvalue,
"p_fisher": stats.fisher_exact(table).pvalue,
"p_chi2": stats.chi2_contingency(table, correction=False)[1] if min(map(sum, table)) > 0 and min(map(sum, zip(*table))) > 0 else np.nan})
rh = pd.DataFrame(rows)
rh["counts_match"] = (rh.nSOZ_positive == 7) & (rh.SOZ_negative == 6)
print(rh.round(4).to_string(index=False))
c = counted(rh)
print("Counted readings:", len(c), "; with both counts matching:", int(c.counts_match.sum()))
z window kind defined hgp_defined nSOZ_positive SOZ_negative p_ranksum p_ttest p_fisher p_chi2 counts_match
trial a_nSOZ claims False 0 NaN NaN NaN NaN NaN NaN False
trial a_all14 claims False 0 NaN NaN NaN NaN NaN NaN False
trial b claims False 2 NaN NaN NaN NaN NaN NaN False
trial c claims True 14 4.0 3.0 0.5653 0.4495 1.0000 1.0000 False
baseline a_nSOZ claims False 0 NaN NaN NaN NaN NaN NaN False
baseline a_all14 claims False 0 NaN NaN NaN NaN NaN NaN False
baseline b claims False 2 NaN NaN NaN NaN NaN NaN False
baseline c claims True 14 4.0 5.0 0.1417 0.0836 0.5921 0.2801 False
amygdala a_nSOZ claims False 0 NaN NaN NaN NaN NaN NaN False
amygdala a_all14 claims False 0 NaN NaN NaN NaN NaN NaN False
amygdala b claims False 3 NaN NaN NaN NaN NaN NaN False
amygdala c claims True 14 5.0 7.0 0.1102 0.1011 0.0210 0.0053 False
trialsd a_nSOZ variant False 0 NaN NaN NaN NaN NaN NaN False
trialsd a_all14 variant False 0 NaN NaN NaN NaN NaN NaN False
trialsd b variant False 2 NaN NaN NaN NaN NaN NaN False
trialsd c variant True 14 5.0 6.0 0.1102 0.0973 0.1026 0.0308 False
baselinepool a_nSOZ variant False 0 NaN NaN NaN NaN NaN NaN False
baselinepool a_all14 variant False 0 NaN NaN NaN NaN NaN NaN False
baselinepool b variant False 4 NaN NaN NaN NaN NaN NaN False
baselinepool c variant True 14 4.0 6.0 0.0350 0.0522 0.2657 0.0943 False
Table 2 printed paper True 14 7.0 6.0 0.0040 0.0016 0.0047 0.0012 True
Counted readings: 3 ; with both counts matching: 0
Result: not reproduced. The readings that define every HGP MI (window c) give 4, 4 and 5 of 7
nSOZ amygdalae with a positive HGP MI and 3, 5 and 7 of 7 SOZ amygdalae with a negative one
(trial, baseline, amygdala); none gives 7 and 6. Under amygdala Fisher's test of the counts
gives p = 0.021, but the nSOZ count is 5, not 7.
This outcome rests on an assumption: window (c), fixed windows taken from the Results' text, is the only window that defines the MIs it needs, while the paper defines the MI over the interval where the two curves are significantly separated (Table 2 legend). The question is for the authors (methods_unavailable).
fig1e-hgp-not-delta-distinguishes-soz¶
The HGP MI, but not the DELTA MI, distinguishes amygdalae within the SOZ from those outside it.
Criterion: under a test, (a) Wilcoxon rank-sum or (b) unpaired t-test, uncorrected, HGP MI differs between nSOZ and SOZ with p < 0.05 and DELTA MI does not; not reproduced when HGP does not differ or DELTA does under every test reading; blocked when the test readings disagree.
rows = []
for z, w, kind, dl, hg, ok in list(readings()) + [("Table 2 printed", "", "paper", paper_delta, paper_hgp, True)]:
for test, fn in [("ranksum", lambda a, b: stats.ranksums(a, b).pvalue), ("ttest", lambda a, b: stats.ttest_ind(a, b).pvalue)]:
if not ok:
rows.append({"z": z, "window": w, "kind": kind, "defined": ok, "test": test})
continue
ph, pdl = fn(hg[~soz_mask], hg[soz_mask]), fn(dl[~soz_mask], dl[soz_mask])
rows.append({"z": z, "window": w, "kind": kind, "defined": ok, "test": test, "p_hgp": ph, "p_delta": pdl,
"hgp_nSOZ_minus_SOZ": np.mean(hg[~soz_mask]) - np.mean(hg[soz_mask]),
"delta_nSOZ_minus_SOZ": np.mean(dl[~soz_mask]) - np.mean(dl[soz_mask]),
"holds": ph < 0.05 and pdl >= 0.05})
re_ = pd.DataFrame(rows)
print(re_.round(4).to_string(index=False))
c = counted(re_)
print("Counted rows (readings x test) where it holds:", int(c.holds.fillna(False).sum()), "of", len(c),
f"; smallest HGP p {c.p_hgp.min():.3f}")
z window kind defined test p_hgp p_delta hgp_nSOZ_minus_SOZ delta_nSOZ_minus_SOZ holds
trial a_nSOZ claims False ranksum NaN NaN NaN NaN NaN
trial a_nSOZ claims False ttest NaN NaN NaN NaN NaN
trial a_all14 claims False ranksum NaN NaN NaN NaN NaN
trial a_all14 claims False ttest NaN NaN NaN NaN NaN
trial b claims False ranksum NaN NaN NaN NaN NaN
trial b claims False ttest NaN NaN NaN NaN NaN
trial c claims True ranksum 0.5653 0.0845 0.1629 -0.0291 False
trial c claims True ttest 0.4495 0.2410 0.1629 -0.0291 False
baseline a_nSOZ claims False ranksum NaN NaN NaN NaN NaN
baseline a_nSOZ claims False ttest NaN NaN NaN NaN NaN
baseline a_all14 claims False ranksum NaN NaN NaN NaN NaN
baseline a_all14 claims False ttest NaN NaN NaN NaN NaN
baseline b claims False ranksum NaN NaN NaN NaN NaN
baseline b claims False ttest NaN NaN NaN NaN NaN
baseline c claims True ranksum 0.1417 0.2774 1.8109 -0.4666 False
baseline c claims True ttest 0.0836 0.5752 1.8109 -0.4666 False
amygdala a_nSOZ claims False ranksum NaN NaN NaN NaN NaN
amygdala a_nSOZ claims False ttest NaN NaN NaN NaN NaN
amygdala a_all14 claims False ranksum NaN NaN NaN NaN NaN
amygdala a_all14 claims False ttest NaN NaN NaN NaN NaN
amygdala b claims False ranksum NaN NaN NaN NaN NaN
amygdala b claims False ttest NaN NaN NaN NaN NaN
amygdala c claims True ranksum 0.1102 0.9491 0.5514 -0.0245 False
amygdala c claims True ttest 0.1011 0.9023 0.5514 -0.0245 False
trialsd a_nSOZ variant False ranksum NaN NaN NaN NaN NaN
trialsd a_nSOZ variant False ttest NaN NaN NaN NaN NaN
trialsd a_all14 variant False ranksum NaN NaN NaN NaN NaN
trialsd a_all14 variant False ttest NaN NaN NaN NaN NaN
trialsd b variant False ranksum NaN NaN NaN NaN NaN
trialsd b variant False ttest NaN NaN NaN NaN NaN
trialsd c variant True ranksum 0.1102 0.4822 0.7202 -0.1072 False
trialsd c variant True ttest 0.0973 0.7279 0.7202 -0.1072 False
baselinepool a_nSOZ variant False ranksum NaN NaN NaN NaN NaN
baselinepool a_nSOZ variant False ttest NaN NaN NaN NaN NaN
baselinepool a_all14 variant False ranksum NaN NaN NaN NaN NaN
baselinepool a_all14 variant False ttest NaN NaN NaN NaN NaN
baselinepool b variant False ranksum NaN NaN NaN NaN NaN
baselinepool b variant False ttest NaN NaN NaN NaN NaN
baselinepool c variant True ranksum 0.0350 0.4822 1.2533 -0.3592 True
baselinepool c variant True ttest 0.0522 0.4178 1.2533 -0.3592 False
Table 2 printed paper True ranksum 0.0040 0.0350 1.0543 -0.5843 False
Table 2 printed paper True ttest 0.0016 0.0133 1.0543 -0.5843 False
Counted rows (readings x test) where it holds: 0 of 6 ; smallest HGP p 0.084
Result: not reproduced. Under the readings that define both MIs (window c) and neither test does
the HGP MI differ between nSOZ and SOZ amygdalae with p < 0.05 (smallest p 0.084, baseline
t-test); the DELTA MI does not differ either (p 0.08 to 0.95). With the baselinepool variant the
rank-sum test gives HGP p = 0.035 and DELTA p = 0.48 (t-test HGP p = 0.052). The printed Table 2
values give HGP p = 0.004 and 0.002, and DELTA p = 0.035 and 0.013, so by these tests the printed
DELTA MI also differs.
This outcome rests on an assumption: window (c), fixed windows taken from the Results' text, is the only window that defines the MIs it needs, while the paper defines the MI over the interval where the two curves are significantly separated (Table 2 legend). The question is for the authors (methods_unavailable).
hgp-soz-sensitivity-specificity¶
HGP enhancement indicates amygdalar functional integrity (outside the SOZ) with sensitivity 7/7 (100%), specificity 6/7 (86%) and accuracy 13/14 (93%).
Criterion: the effect is a positive HGP MI; counts must match; Clopper-Pearson 95% intervals.
rows = []
for z, w, kind, dl, hg, ok in list(readings(need=("hgp",))) + [("Table 2 printed", "", "paper", paper_delta, paper_hgp, True)]:
if not ok:
rows.append({"z": z, "window": w, "kind": kind, "defined": ok})
continue
tp, fn_ = int((hg[~soz_mask] > 0).sum()), int((hg[~soz_mask] <= 0).sum())
fp, tn = int((hg[soz_mask] > 0).sum()), int((hg[soz_mask] <= 0).sum())
rows.append({"z": z, "window": w, "kind": kind, "defined": ok, "TP": tp, "FN": fn_, "FP": fp, "TN": tn,
"sensitivity": f"{tp}/7", "specificity": f"{tn}/7", "accuracy": f"{tp + tn}/14",
"CI_sens": tuple(round(100 * v) for v in A.clopper_pearson(tp, 7)),
"CI_spec": tuple(round(100 * v) for v in A.clopper_pearson(tn, 7)),
"CI_acc": tuple(round(100 * v) for v in A.clopper_pearson(tp + tn, 14))})
rs = pd.DataFrame(rows)
rs["counts_match"] = (rs.TP == 7) & (rs.TN == 6)
print(rs.to_string(index=False))
c = counted(rs)
print("Counted readings:", len(c), "; with the counts matching:", int(c.counts_match.sum()))
print("Paper: 7/7 CI [66-100], 6/7 CI [42-99], 13/14 CI [66-99]")
# Other binomial intervals, to see which one the printed bounds follow.
from statsmodels.stats.proportion import proportion_confint
for k, n_ in [(7, 7), (6, 7), (13, 14)]:
cp1 = (stats.beta.ppf(0.05, k, n_ - k + 1) if k > 0 else 0.0, stats.beta.ppf(0.95, k + 1, n_ - k) if k < n_ else 1.0)
print(f"{k}/{n_}: two-sided Clopper-Pearson {tuple(round(100 * v) for v in A.clopper_pearson(k, n_))}; "
f"one-sided 95% bounds {tuple(round(100 * v) for v in cp1)}; Wilson {tuple(round(100 * v) for v in proportion_confint(k, n_, method='wilson'))}; "
f"Jeffreys {tuple(round(100 * v) for v in proportion_confint(k, n_, method='jeffreys'))}")
z window kind defined TP FN FP TN sensitivity specificity accuracy CI_sens CI_spec CI_acc counts_match
trial a_nSOZ claims False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
trial a_all14 claims False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
trial b claims False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
trial c claims True 4.0 3.0 4.0 3.0 4/7 3/7 7/14 (18, 90) (10, 82) (23, 77) False
baseline a_nSOZ claims False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
baseline a_all14 claims False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
baseline b claims False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
baseline c claims True 4.0 3.0 2.0 5.0 4/7 5/7 9/14 (18, 90) (29, 96) (35, 87) False
amygdala a_nSOZ claims False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
amygdala a_all14 claims False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
amygdala b claims False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
amygdala c claims True 5.0 2.0 0.0 7.0 5/7 7/7 12/14 (29, 96) (59, 100) (57, 98) False
trialsd a_nSOZ variant False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
trialsd a_all14 variant False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
trialsd b variant False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
trialsd c variant True 5.0 2.0 1.0 6.0 5/7 6/7 11/14 (29, 96) (42, 100) (49, 95) False
baselinepool a_nSOZ variant False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
baselinepool a_all14 variant False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
baselinepool b variant False NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN False
baselinepool c variant True 4.0 3.0 1.0 6.0 4/7 6/7 10/14 (18, 90) (42, 100) (42, 92) False
Table 2 printed paper True 7.0 0.0 1.0 6.0 7/7 6/7 13/14 (59, 100) (42, 100) (66, 100) True
Counted readings: 3 ; with the counts matching: 0
Paper: 7/7 CI [66-100], 6/7 CI [42-99], 13/14 CI [66-99]
7/7: two-sided Clopper-Pearson (59, 100); one-sided 95% bounds (65, 100); Wilson (65, 100); Jeffreys (71, 100)
6/7: two-sided Clopper-Pearson (42, 100); one-sided 95% bounds (48, 99); Wilson (49, 97); Jeffreys (50, 98)
13/14: two-sided Clopper-Pearson (66, 100); one-sided 95% bounds (70, 100); Wilson (69, 99); Jeffreys (71, 99)
Result: not reproduced. Window c gives sensitivity 4/7, 4/7 and 5/7, specificity 3/7, 5/7 and
7/7, accuracy 7/14, 9/14 and 12/14 (trial, baseline, amygdala); none gives 7/7, 6/7 and 13/14.
The two-sided Clopper-Pearson intervals of the paper's counts (last row) are 59-100, 42-100 and
66-100%, against the printed 66-100, 42-99 and 66-99%. No single interval method gives all three
printed intervals: the one-sided 95% bound gives 65 for 7/7, Wilson gives 65-100, 49-97 and 69-99.
This outcome rests on an assumption: window (c), fixed windows taken from the Results' text, is the only window that defines the MIs it needs, while the paper defines the MI over the interval where the two curves are significantly separated (Table 2 legend). The question is for the authors (methods_unavailable).
fig1f-hgp-more-specific-than-delta¶
The HGP response is more specific than the DELTA response for indicating amygdalar dysfunction within the SOZ, by a chi-square test that is not significant (p about 0.2).
Criterion: reading (a), over the 7 SOZ amygdalae, HGP effect absent (HGP MI not positive) against DELTA effect absent (DELTA MI not negative), chi-square of the 2 x 2 table; reproduced when HGP specificity exceeds DELTA specificity under (a); (b), the printed counts 5/5 against 1/5, is reported; blocked if the readings disagree on which is more specific.
rows = []
for z, w, kind, dl, hg, ok in list(readings(amys=soz_mask)) + [("Table 2 printed", "", "paper", paper_delta, paper_hgp, True)]:
if not ok:
rows.append({"z": z, "window": w, "kind": kind, "defined": ok})
continue
sh, sd = int((hg[soz_mask] <= 0).sum()), int((dl[soz_mask] >= 0).sum())
table = np.array([[sh, 7 - sh], [sd, 7 - sd]])
p = stats.chi2_contingency(table, correction=False)[1] if (table.sum(0) > 0).all() else np.nan
rows.append({"z": z, "window": w, "kind": kind, "defined": ok, "HGP_specific": f"{sh}/7", "DELTA_specific": f"{sd}/7",
"which": "HGP" if sh > sd else "DELTA" if sd > sh else "tie", "p_chi2": p})
rf = pd.DataFrame(rows)
print(rf.round(3).to_string(index=False))
pb = stats.chi2_contingency(np.array([[5, 0], [1, 4]]), correction=False)[1]
pby = stats.chi2_contingency(np.array([[5, 0], [1, 4]]), correction=True)[1]
print(f"Reading (b), the printed counts 5/5 against 1/5: chi-square p = {pb:.4f} (with Yates' correction {pby:.4f}); paper p ~ 0.2")
c = counted(rf)
print("Counted readings:", len(c), "; more specific:", c.which.value_counts().to_dict())
z window kind defined HGP_specific DELTA_specific which p_chi2
trial a_nSOZ claims False NaN NaN NaN NaN
trial a_all14 claims False NaN NaN NaN NaN
trial b claims False NaN NaN NaN NaN
trial c claims True 3/7 3/7 tie 1.000
baseline a_nSOZ claims False NaN NaN NaN NaN
baseline a_all14 claims False NaN NaN NaN NaN
baseline b claims False NaN NaN NaN NaN
baseline c claims True 5/7 2/7 HGP 0.109
amygdala a_nSOZ claims False NaN NaN NaN NaN
amygdala a_all14 claims False NaN NaN NaN NaN
amygdala b claims False NaN NaN NaN NaN
amygdala c claims True 7/7 1/7 HGP 0.001
trialsd a_nSOZ variant False NaN NaN NaN NaN
trialsd a_all14 variant False NaN NaN NaN NaN
trialsd b variant False NaN NaN NaN NaN
trialsd c variant True 6/7 1/7 HGP 0.008
baselinepool a_nSOZ variant False NaN NaN NaN NaN
baselinepool a_all14 variant False NaN NaN NaN NaN
baselinepool b variant False NaN NaN NaN NaN
baselinepool c variant True 6/7 3/7 HGP 0.094
Table 2 printed paper True 6/7 2/7 HGP 0.031
Reading (b), the printed counts 5/5 against 1/5: chi-square p = 0.0098 (with Yates' correction 0.0528); paper p ~ 0.2
Counted readings: 3 ; more specific: {'HGP': 2, 'tie': 1}
Result: blocked (methods_underspecified). The readings disagree on which response is more
specific. With window c, the only window that defines the MIs of all 7 SOZ amygdalae, HGP is more
specific under baseline (5/7 against 2/7, chi-square p = 0.11) and amygdala (7/7 against 1/7,
p = 0.001), and the two are tied under trial (3/7 each). The variants also give HGP (6/7 against
1/7 and 3/7). Under reading (b), the printed counts, HGP is more specific with p = 0.0098 (0.053 with
Yates' correction), not about 0.2; no reading gives a p near 0.2.
This outcome rests on an assumption: window (c), fixed windows taken from the Results' text, is the only window that defines the MIs it needs, while the paper defines the MI over the interval where the two curves are significantly separated (Table 2 legend). The question is for the authors (methods_unavailable).
print("Units in the dandiset:", len(units_df))
print("Units with rate > 1 Hz over all 17 trials (442 s):", int((units_df.rate_442s > 1).sum()))
print("Units with rate > 1 Hz over the epochs each unit was recorded in:", int((units_df.rate_observed > 1).sum()))
print("Units at or below 1 Hz over 442 s:", units_df[units_df.rate_442s <= 1][["subject", "unit", "rate_442s", "rate_observed"]].to_dict("records"))
print("Amygdalae with units:", units_df.amygdala.nunique(), sorted(units_df.amygdala.unique()),
"; patients with units:", units_df.subject.nunique(), sorted(units_df.subject.unique()))
print("Rates over observed epochs: mean %.2f, SD %.2f, min %.2f, max %.2f Hz (data descriptor: 1.66 +- 2.65 Hz)"
% (units_df.rate_observed.mean(), units_df.rate_observed.std(), units_df.rate_observed.min(), units_df.rate_observed.max()))
print("Paper Table 1 lists MUA in 1 L/R, 2 L/R, 3 L, 4 L, 5 R, 8 L; units here are in subjects",
units_df.groupby("subject").channel.unique().apply(list).to_dict())
Units in the dandiset: 42
Units with rate > 1 Hz over all 17 trials (442 s): 40
Units with rate > 1 Hz over the epochs each unit was recorded in: 42
Units at or below 1 Hz over 442 s: [{'subject': 7, 'unit': 1, 'rate_442s': 0.88, 'rate_observed': 1.36}, {'subject': 7, 'unit': 5, 'rate_442s': 0.88, 'rate_observed': 1.36}]
Amygdalae with units: 8 [np.int64(1), np.int64(2), np.int64(3), np.int64(4), np.int64(5), np.int64(6), np.int64(10), np.int64(11)] ; patients with units: 5 [np.int64(1), np.int64(2), np.int64(3), np.int64(4), np.int64(7)]
Rates over observed epochs: mean 4.65, SD 3.91, min 1.09, max 16.33 Hz (data descriptor: 1.66 +- 2.65 Hz)
Paper Table 1 lists MUA in 1 L/R, 2 L/R, 3 L, 4 L, 5 R, 8 L; units here are in subjects {1: ['AL2-AL1', 'AR2-AR1'], 2: ['AL2-AL1', 'AR2-AR1'], 3: ['AL2-AL1'], 4: ['AL2-AL1'], 7: ['AL2-AL1', 'AR2-AR1']}
Result: partially reproduced. The dandiset holds 42 units: 42 above 1 Hz over the epochs each was recorded in and 40 over the whole 442 s, both within 4.1 of 41. The 50 putative clusters cannot be tested, as only units above about 1 Hz were deposited (two at 0.88 Hz over the 442 s). The units come from 8 amygdalae of 5 patients, not 9 amygdalae. They are not the amygdalae of Table 1: subject 7 (L and R) has units here and none in Table 1, and subjects 5 (R) and 8 (L) have them in Table 1 and none here. Their mean rate is 4.65 +- 3.91 Hz (median 3.13), against the data descriptor's 1.66 +- 2.65 Hz; the table below gives each unit's rate in the baseline clips and in each condition's blocks.
# Per-unit rates by part of the epoch, and the population PSTH, to compare with the data descriptor's
# 1.66 +- 2.65 Hz (its Fig. 2c, "average firing rate for all neurons").
t100 = np.arange(260) * 0.1 + 0.05 - A.ONSET_S
rows = []
for u in U:
r, o, av = u["rate100"], u["obs"], u["aversive"]
rows.append({"subject": u["subject"], "unit": u["unit"], "amygdala": u["amygdala"],
"baseline_clip": r[o][:, t100 < 0].mean(), "aversive_block": r[o & av][:, t100 >= 0].mean(),
"neutral_block": r[o & ~av][:, t100 >= 0].mean(), "all_observed": r[o].mean(), "over_442s": r.mean()})
rates = pd.DataFrame(rows)
print(rates.round(2).to_string(index=False))
print("\nMean +- SD over units:", {k: f"{rates[k].mean():.2f} +- {rates[k].std(ddof=1):.2f}" for k in rates.columns[3:]})
print("Median over units:", rates[rates.columns[3:]].median().round(2).to_dict())
print("Units with rate above 1, 2 and 5 Hz over the observed epochs:", [int((rates.all_observed > k).sum()) for k in (1, 2, 5)])
fig, ax = plt.subplots(figsize=(7, 3))
for cond, col in [(True, "C3"), (False, "C0")]:
psth = np.mean([A.smooth_epochs(u["rate100"], 10)[u["obs"] & (u["aversive"] == cond)].mean(0) for u in U], axis=0)
ax.plot(t100, psth, color=col, label="AVERSIVE" if cond else "NEUTRAL")
ax.axvline(0, color="k", lw=0.5); ax.set_xlabel("time from block onset (s)"); ax.set_ylabel("mean rate over units (Hz)"); ax.legend(fontsize=8)
print("Population PSTH plotted (100 ms bins, 1 s smoothing), range within the autoscaled axes.")
subject unit amygdala baseline_clip aversive_block neutral_block all_observed over_442s
1 0 1 1.68 1.21 1.82 1.55 1.55
1 1 1 1.03 0.90 1.38 1.14 1.14
1 2 1 1.32 1.65 1.44 1.52 1.52
1 3 2 1.53 2.18 2.14 2.11 2.11
2 0 3 10.68 10.80 11.98 11.36 11.36
2 1 3 6.82 5.32 4.93 5.24 5.24
2 2 4 2.75 2.77 1.92 2.38 2.24
2 3 4 2.53 2.37 1.89 2.16 2.03
2 4 4 16.79 13.45 15.24 14.58 14.58
2 5 4 4.09 3.70 3.36 3.56 3.56
2 6 4 9.59 9.96 10.41 10.15 10.15
2 7 4 4.18 2.93 2.83 2.98 2.98
2 8 4 6.44 6.07 5.17 5.68 5.35
2 9 4 2.84 2.99 2.74 2.86 2.69
2 10 4 4.65 5.33 5.00 5.11 5.11
2 11 3 1.35 1.24 0.91 1.09 1.09
3 0 5 3.38 2.77 2.59 2.73 2.73
3 1 5 2.38 2.09 1.41 1.78 1.78
3 2 5 3.26 2.98 2.00 2.52 2.52
4 0 6 16.18 16.10 16.56 16.33 16.33
4 1 6 4.53 4.28 4.30 4.31 4.31
4 2 6 1.68 2.03 1.86 1.92 1.92
4 3 6 5.97 7.92 6.75 7.19 7.19
4 4 6 1.15 1.72 1.47 1.55 1.55
4 5 6 5.12 5.71 5.11 5.37 5.37
4 6 6 8.79 8.06 8.43 8.30 8.30
7 0 10 3.59 3.33 3.18 3.28 2.12
7 1 10 1.14 1.48 1.30 1.36 0.88
7 2 10 5.00 4.21 3.96 4.14 2.68
7 3 10 0.68 1.72 1.70 1.63 1.05
7 4 10 2.23 2.00 1.38 1.70 1.10
7 5 10 1.00 2.05 0.85 1.36 0.88
7 6 10 3.14 4.01 4.44 4.16 2.69
7 7 10 4.27 4.01 5.81 4.94 3.19
7 8 10 9.36 8.94 8.60 8.80 5.69
7 9 11 2.18 2.36 2.15 2.24 1.45
7 10 11 2.64 1.85 2.40 2.19 1.42
7 11 11 5.59 4.64 5.20 5.00 3.23
7 12 11 2.73 3.32 2.56 2.89 1.87
7 13 11 6.55 7.21 7.87 7.49 4.85
7 14 11 4.50 3.53 3.87 3.78 2.44
7 15 11 15.00 14.54 15.48 15.05 9.74
Mean +- SD over units: {'baseline_clip': '4.77 +- 4.02', 'aversive_block': '4.66 +- 3.73', 'neutral_block': '4.63 +- 4.08', 'all_observed': '4.65 +- 3.91', 'over_442s': '4.05 +- 3.68'}
Median over units: {'baseline_clip': 3.49, 'aversive_block': 3.32, 'neutral_block': 3.01, 'all_observed': 3.13, 'over_442s': 2.69}
Units with rate above 1, 2 and 5 Hz over the observed epochs: [42, 31, 13]
Population PSTH plotted (100 ms bins, 1 s smoothing), range within the autoscaled axes.
fig3a-firing-increase-aversive¶
All 41 neurons fire more during AVERSIVE than during NEUTRAL sequences (FRV 0.26 plus or minus 0.15).
Criterion: FRV (AVERSIVE minus NEUTRAL mean rate) over (a) the 24 s blocks or (b) the 26 s trials; reproduced when every neuron has FRV > 0 and the mean is within 0.15 of 0.26 under one reading; partially reproduced when the mean is within tolerance and most neurons have FRV > 0.
rows = []
for window in ["block", "trial"]:
for only in [True, False]:
f = np.array([A.frv(u, window, only) for u in U])
rows.append({"window": window, "observed_epochs_only": only, "n": len(f), "mean_FRV_Hz": f.mean(), "sd": f.std(ddof=1),
"positive": int((f > 0).sum()), "wilcoxon_p": stats.wilcoxon(f).pvalue})
rfrv = pd.DataFrame(rows)
print(rfrv.round(4).to_string(index=False))
FRV = np.array([A.frv(u, "block") for u in U])
units_df["FRV"] = FRV
print("\nFRV by subject (blocks, observed epochs):")
print(units_df.groupby("subject").FRV.agg(["count", "mean", "min", "max", lambda v: int((v > 0).sum())]).round(3).to_string())
# Variant (not the claim's measure): FRV relative to the unit's mean rate, and in units of the SD of its 100 ms rate.
rel = np.array([A.frv(u, "block") / u["rate100"][u["obs"]].mean() for u in U])
zfrv = np.array([A.frv(u, "block") / u["rate100"][u["obs"]].std() for u in U])
print(f"\nVariant: FRV / mean rate {rel.mean():.3f} +- {rel.std(ddof=1):.3f} ({(rel > 0).sum()} positive); "
f"FRV / SD of 100 ms rate {zfrv.mean():.3f} +- {zfrv.std(ddof=1):.3f}")
window observed_epochs_only n mean_FRV_Hz sd positive wilcoxon_p
block True 42 0.0329 0.6903 27 0.5035
block False 42 0.0075 0.6853 23 0.7854
trial True 42 0.0319 0.6534 26 0.5114
trial False 42 0.0069 0.6569 23 0.7665
FRV by subject (blocks, observed epochs):
count mean min max <lambda_0>
subject
1 4 -0.211 -0.606 0.206 2
2 12 0.046 -1.793 0.896 9
3 3 0.616 0.183 0.989 3
4 7 0.195 -0.456 1.171 4
7 16 -0.096 -1.804 1.203 9
Variant: FRV / mean rate 0.061 +- 0.231 (27 positive); FRV / SD of 100 ms rate 0.020 +- 0.092
Result: not reproduced. Over the blocks and the observed epochs, 27 of 42 units have FRV > 0, mean 0.033 +- 0.690 Hz (Wilcoxon against 0, p = 0.50); the other readings give 23 to 26 positive and means of 0.007 to 0.032 Hz. The mean is outside 0.26 +- 0.15 and not all units fire more in AVERSIVE. Relative to each unit's mean rate the FRV is 0.061 +- 0.231 (a variant, not the claim's measure).
fig3a-firing-gradual-sustained¶
The neurons' firing rose gradually and stayed enhanced throughout the AVERSIVE sequences, not as an immediate response at onset.
Criterion: population AVERSIVE minus NEUTRAL rate (100 ms bins, 1 s smoothing), (a) the difference over 0 to 2 s smaller than over 2 to 24 s and positive in each of the twelve 2 s bins; (b) its first peak later than 2 s and positive to 24 s; reproduced when both hold, not reproduced when both fail.
tr = np.arange(260) * 0.1 + 0.05 - A.ONSET_S
diffs = []
for u in U:
r = A.smooth_epochs(u["rate100"], 10)
o, av = u["obs"], u["aversive"]
diffs.append(r[o & av].mean(0) - r[o & ~av].mean(0))
diffs = np.array(diffs)
pop = diffs.mean(0)
pb = pop[tr >= 0]
tblk = tr[tr >= 0]
bins2 = [pb[(tblk >= k) & (tblk < k + 2)].mean() for k in range(0, 24, 2)]
early, late = pb[tblk < 2].mean(), pb[tblk >= 2].mean()
print("Population AVERSIVE - NEUTRAL (Hz) in 2 s bins:", np.round(bins2, 3).tolist())
print(f"(a) mean 0-2 s {early:.3f}, mean 2-24 s {late:.3f}; early < late: {early < late}; all twelve bins positive: {all(b > 0 for b in bins2)}")
from scipy.signal import find_peaks
pk = find_peaks(pb)[0]
first_peak = tblk[pk[0]] if len(pk) else np.nan
print(f"(b) first local peak at {first_peak:.2f} s; later than 2 s: {first_peak > 2}; positive at every point 2-24 s: {bool((pb[tblk >= 2] > 0).all())}; "
f"fraction of 0-24 s positive {np.mean(pb > 0):.2f}")
print(f"Baseline clip (-2 to 0 s) mean difference {pop[tr < 0].mean():.3f} Hz")
print("Per subject, mean 2-24 s difference:", {s: round(float(diffs[[u['subject'] == s for u in U]][:, tr >= 2].mean()), 3) for s in sorted(set(u['subject'] for u in U))})
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(tr, pop, "k")
ax.fill_between(tr, pop - diffs.std(0, ddof=1) / np.sqrt(len(U)), pop + diffs.std(0, ddof=1) / np.sqrt(len(U)), color="grey", alpha=0.3)
ax.axhline(0, color="grey", lw=0.5); ax.axvline(0, color="r", lw=0.5)
ax.set_xlabel("time from block onset (s)"); ax.set_ylabel("AVERSIVE - NEUTRAL (Hz)")
print(f"Plotted range {pop.min():.3f} to {pop.max():.3f} Hz (mean +- SEM band within the autoscaled axes)")
Population AVERSIVE - NEUTRAL (Hz) in 2 s bins: [0.31, -0.142, 0.066, 0.278, 0.459, -0.14, 0.069, 0.141, 0.195, -0.217, -0.06, -0.449]
(a) mean 0-2 s 0.310, mean 2-24 s 0.018; early < late: False; all twelve bins positive: False
(b) first local peak at 0.15 s; later than 2 s: False; positive at every point 2-24 s: False; fraction of 0-24 s positive 0.57
Baseline clip (-2 to 0 s) mean difference -0.078 Hz
Per subject, mean 2-24 s difference: {1: -0.15, 2: -0.029, 3: 0.543, 4: 0.218, 7: -0.09}
Plotted range -0.635 to 0.958 Hz (mean +- SEM band within the autoscaled axes)
Result: not reproduced. The population AVERSIVE minus NEUTRAL rate is highest at onset: 0.310 Hz over 0 to 2 s against 0.018 Hz over 2 to 24 s, its first peak at 0.15 s, and 4 of the twelve 2 s bins negative (positive at 57% of the time points). Neither reading of "gradual and sustained" holds. The difference is -0.078 Hz over the baseline clip, so the rise begins at the block onset of the trial table, consistent with that onset.
fig3b-firing-hgp-no-condition-difference¶
Across neurons, the correlation of firing rate with macro-electrode HGP does not differ between the AVERSIVE and NEUTRAL conditions.
A claim of absence: Wilcoxon signed-rank test of the per-neuron Spearman r (AVERSIVE against NEUTRAL)
at p < 0.05. HGP is the Methods 4.8 envelope of the bipolar channel of the neuron's amygdala, without
per-trial normalisation (a Spearman r within one block is unchanged by any per-trial normalisation).
The claims file's readings: time bases (a) 1 s, (b) 1.5 s with 5 points, (c) 15 s; spans (i) the
condition's blocks laid end to end, (ii) per block and averaged. In (i) each block is smoothed (or
binned into 1.5 s bins) on its own before the blocks are laid end to end, so that no sample of a
condition's series carries the other condition or a baseline clip. A variant, concat_session,
smooths over the whole session first (the first revision's form), to show what that mixing does.
HGP100 = {n: ENV[(n, "hgp")] for n in ALL}
rows = []
per_neuron = {}
for base in ["1s", "1.5s", "15s"]:
for span in ["concat", "block", "concat_session"]:
rr = np.array([A.condition_corr(u["rate100"], HGP100[u["amygdala"]], u["obs"], u["aversive"], base, span) for u in U])
ok = ~np.isnan(rr).any(1)
per_neuron[(base, span)] = rr
w = stats.wilcoxon(rr[ok, 0], rr[ok, 1])
rows.append({"base": base, "span": span, "kind": "variant" if span == "concat_session" else "claims", "n": int(ok.sum()),
"r_aversive": f"{np.nanmean(rr[:, 0]):.3f} +- {np.nanstd(rr[:, 0], ddof=1):.3f}",
"r_neutral": f"{np.nanmean(rr[:, 1]):.3f} +- {np.nanstd(rr[:, 1], ddof=1):.3f}",
"median_diff": np.median(rr[ok, 0] - rr[ok, 1]), "wilcoxon_p": w.pvalue})
r3b = pd.DataFrame(rows)
print(r3b.round(4).to_string(index=False))
c = r3b[r3b.kind == "claims"]
print("Readings of the claims file with a difference at p < 0.05:", int((c.wilcoxon_p < 0.05).sum()), "of", len(c))
print("Paper: 0.17 +- 0.28 AVERSIVE, 0.18 +- 0.18 NEUTRAL, Wilcoxon p > 0.05")
rr = per_neuron[("1s", "concat")]
print("\n1 s, blocks laid end to end, by subject, mean r AVERSIVE and NEUTRAL:")
print(pd.DataFrame({"subject": units_df.subject, "r_A": rr[:, 0], "r_N": rr[:, 1]}).groupby("subject").agg(["count", "mean"]).round(3).to_string())
base span kind n r_aversive r_neutral median_diff wilcoxon_p
1s concat claims 42 0.085 +- 0.130 0.094 +- 0.124 -0.0165 0.5357
1s block claims 42 0.086 +- 0.134 0.080 +- 0.116 0.0057 0.8334
1s concat_session variant 42 0.083 +- 0.131 0.096 +- 0.123 -0.0118 0.3718
1.5s concat claims 42 0.083 +- 0.216 0.164 +- 0.216 -0.0874 0.0122
1.5s block claims 42 0.111 +- 0.198 0.117 +- 0.180 -0.0234 0.7197
1.5s concat_session variant 42 0.083 +- 0.218 0.176 +- 0.215 -0.0852 0.0020
15s concat claims 42 0.039 +- 0.234 0.118 +- 0.226 -0.1305 0.0225
15s block claims 42 0.096 +- 0.218 0.076 +- 0.189 -0.0102 0.8045
15s concat_session variant 42 0.076 +- 0.292 0.167 +- 0.253 -0.1074 0.0090
Readings of the claims file with a difference at p < 0.05: 2 of 6
Paper: 0.17 +- 0.28 AVERSIVE, 0.18 +- 0.18 NEUTRAL, Wilcoxon p > 0.05
1 s, blocks laid end to end, by subject, mean r AVERSIVE and NEUTRAL:
r_A r_N
count mean count mean
subject
1 4 0.015 4 0.052
2 12 0.193 12 0.220
3 3 0.265 3 0.202
4 7 0.005 7 0.041
7 16 0.021 16 0.013
Result: partially reproduced. The readings of the claims file disagree. Per block and averaged,
no time base finds a difference (p 0.72 to 0.83), and neither does the 1 s time base over the blocks
laid end to end (r = 0.085 +- 0.130 in AVERSIVE and 0.094 +- 0.124 in NEUTRAL, p = 0.54). Over the
blocks laid end to end on the 1.5 s and 15 s time bases, each block smoothed on its own, the test
finds NEUTRAL higher (p = 0.012 and 0.023, median differences -0.087 and -0.131). So the claimed
absence holds in 4 of the 6 readings. Smoothing over the whole session before selecting the blocks
(variant concat_session) gives p = 0.002 and 0.009 for those two, so the mixing of neighbouring
blocks strengthened the difference but does not explain it. The means are within the stated spreads
under every reading.
The outcome is partially reproduced, not blocked: the six readings are named by the claims file, so no method detail is left to choose; the test was run under each, and the claimed absence holds under some and fails under others, that is, it holds only in part. Unlike the criteria of fig1b-delta-decrease and fig1f-hgp-more-specific-than-delta, this criterion does not say to block when its readings disagree.
fig3c-ssi-frv-correlation¶
Across neurons, the spatial synchronization index (SSI) is strongly positively correlated with the firing rate variability (FRV).
SSI per neuron is the mean Spearman correlation of its firing rate with each neuron of the same amygdala on a different microwire (Methods 4.11, "neighboring neurons recorded in the same amygdala but on different microwires"). The units table is printed first, to see whether it records the microwire. Where it does not, every other unit of the amygdala is used, an assumption (the Results say "all neighboring neurons"), and a variant drops the pairs that are likely on one wire by a spike train test: units sorted from one microwire cannot both be detected within about 1 ms of each other, so a pair with no spikes within +-1 ms where the flanks (2 to 12 ms) predict at least 3 is taken to share a wire. Time bases (a) 100 ms smoothed over 1 s, (b) 15 s, over the whole session (the epochs both units were recorded in). FRV over the blocks.
h_units = A.open_asset(next(a for a in ecephys if a["path"].startswith("sub-07"))["asset_id"])
print("Units table of subject 7, columns:", list(h_units["units"].attrs["colnames"]))
print("Electrode rows the units point to:", h_units["units/electrodes"][:].tolist(),
"; electrode table labels:", [A._s(x) for x in h_units["general/extracellular_ephys/electrodes/label"][:]],
"; electrode groups in the file:", [k for k in h_units["general/extracellular_ephys"].keys() if k != "electrodes"])
pairs = []
for i, u in enumerate(U):
for j, v in enumerate(U):
if j <= i or v["subject"] != u["subject"] or v["amygdala"] != u["amygdala"]:
continue
near, expected = A.coincidences(u["spikes"], v["spikes"])
pairs.append({"i": i, "j": j, "subject": u["subject"], "amygdala": u["amygdala"], "near_1ms": near,
"expected": round(expected, 1), "likely_same_wire": near == 0 and expected >= 3,
"testable": expected >= 3})
pairs = pd.DataFrame(pairs)
print(f"\nPairs in one amygdala: {len(pairs)}; testable (>= 3 expected within +-1 ms): {int(pairs.testable.sum())}; "
f"likely on one wire: {int(pairs.likely_same_wire.sum())}")
print(pairs.groupby("amygdala")[["testable", "likely_same_wire"]].sum().astype(int).to_string())
print("Ratio near/expected over testable pairs, quartiles:",
np.round(np.percentile((pairs.near_1ms / pairs.expected)[pairs.testable], [25, 50, 75]), 2).tolist())
same_wire = {(r.i, r.j) for r in pairs.itertuples() if r.likely_same_wire}
SSI = {}
for base in ["1s", "15s"]:
for drop in (False, True):
vals = []
for i, u in enumerate(U):
rs_ = []
for j, v in enumerate(U):
if j == i or v["subject"] != u["subject"] or v["amygdala"] != u["amygdala"]:
continue
if drop and (min(i, j), max(i, j)) in same_wire:
continue
o = u["obs"] & v["obs"]
a, _, _ = A.series(u["rate100"], o, base)
b, _, _ = A.series(v["rate100"], o, base)
rs_.append(A.spearman(a, b))
vals.append(np.mean(rs_) if rs_ else np.nan)
SSI[base + ("_xwire" if drop else "")] = np.array(vals)
units_df["SSI_1s"], units_df["SSI_15s"] = SSI["1s"], SSI["15s"]
print("\nUnits with an SSI (another unit in the same amygdala):", int(np.isfinite(SSI["1s"]).sum()), "of", len(U),
"; after dropping likely same-wire pairs:", int(np.isfinite(SSI["1s_xwire"]).sum()))
for base in SSI:
ok = np.isfinite(SSI[base])
r = stats.spearmanr(SSI[base][ok], FRV[ok])
print(f"base {base}{' (variant)' if 'xwire' in base else ''}: SSI mean {np.nanmean(SSI[base]):.3f} "
f"(range {np.nanmin(SSI[base]):.3f} to {np.nanmax(SSI[base]):.3f}); Spearman SSI vs FRV r = {r.statistic:.3f}, p = {r.pvalue:.2g}, n = {ok.sum()}")
if "xwire" in base:
continue
for s in sorted(units_df.subject.unique()):
m = ok & (units_df.subject.values == s)
if m.sum() > 2:
rr = stats.spearmanr(SSI[base][m], FRV[m])
print(f" subject {s}: n = {m.sum()}, r = {rr.statistic:.3f}, p = {rr.pvalue:.3f}")
print("Paper: r = 0.78, p < 1e-7")
Units table of subject 7, columns: ['offset', 'spike_times', 'obs_intervals', 'electrodes', 'electrode_group', 'waveform_mean', 'waveform_sd'] Electrode rows the units point to: [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1] ; electrode table labels: ['AL2-AL1', 'AR2-AR1'] ; electrode groups in the file: ['ieeg']
Pairs in one amygdala: 123; testable (>= 3 expected within +-1 ms): 116; likely on one wire: 4
testable likely_same_wire
amygdala
1 1 0
3 3 0
4 36 1
5 3 0
6 21 1
10 31 2
11 21 0
Ratio near/expected over testable pairs, quartiles: [0.94, 1.19, 1.53]
Units with an SSI (another unit in the same amygdala): 41 of 42 ; after dropping likely same-wire pairs: 41 base 1s: SSI mean 0.111 (range -0.018 to 0.337); Spearman SSI vs FRV r = 0.255, p = 0.11, n = 41 subject 1: n = 3, r = 0.500, p = 0.667 subject 2: n = 12, r = 0.140, p = 0.665 subject 3: n = 3, r = 0.500, p = 0.667 subject 4: n = 7, r = 0.214, p = 0.645 subject 7: n = 16, r = 0.018, p = 0.948 base 1s_xwire (variant): SSI mean 0.109 (range -0.018 to 0.337); Spearman SSI vs FRV r = 0.243, p = 0.13, n = 41 base 15s: SSI mean 0.133 (range -0.194 to 0.557); Spearman SSI vs FRV r = 0.403, p = 0.009, n = 41 subject 1: n = 3, r = 0.500, p = 0.667 subject 2: n = 12, r = 0.413, p = 0.183 subject 3: n = 3, r = 1.000, p = 0.000 subject 4: n = 7, r = 0.214, p = 0.645 subject 7: n = 16, r = 0.229, p = 0.393 base 15s_xwire (variant): SSI mean 0.129 (range -0.194 to 0.557); Spearman SSI vs FRV r = 0.353, p = 0.023, n = 41 Paper: r = 0.78, p < 1e-7
Result: not reproduced (resting on an assumption). The units table records no microwire: each unit points only to the bipolar macro channel of its amygdala, and the file holds one electrode group. SSI therefore uses every other unit of the amygdala (the assumption). The spike train test finds 4 of 116 testable pairs likely on one wire; the other pairs fire within 1 ms of each other about as often as the flanks predict (median ratio 1.19). SSI and FRV are positively correlated, but weakly: r = 0.26 (p = 0.11) on the 1 s time base and r = 0.40 (p = 0.009) on the 15 s one, over the 41 units with another unit in their amygdala, against r = 0.78, p < 1e-7. Dropping the 4 pairs gives 0.24 (p = 0.13) and 0.35 (p = 0.023). The criterion needs p < 0.05 under every reading for a partial result.
fig3d-firing-correlations-modalities (HGP part)¶
Across neurons, the ongoing firing rate correlates positively with muHGP (0.27), HGP (0.19) and BOLD (0.23), and the correlation is lowest with HGP.
Per neuron, Spearman r of firing rate with the HGP envelope of its amygdala over the session (observed epochs, baseline clips included), (a) 15 s smoothing (Fig. 3D legend), (b) 1.5 s bins smoothed over 5 points (Methods 4.10). muHGP and BOLD are not in the dandiset, so the means for them and the two Wilcoxon tests cannot be computed.
RH = {}
for base in ["15s", "1.5s", "1s"]:
vals = []
for u in U:
a, _, _ = A.series(u["rate100"], u["obs"], base)
b, _, _ = A.series(HGP100[u["amygdala"]], u["obs"], base)
vals.append(A.spearman(a, b))
RH[base] = np.array(vals)
print(f"base {base}: firing-HGP r = {np.nanmean(RH[base]):.3f} +- {np.nanstd(RH[base], ddof=1):.3f} "
f"(n = {np.isfinite(RH[base]).sum()}, {int((RH[base] > 0).sum())} positive); Wilcoxon against 0 p = {stats.wilcoxon(RH[base][np.isfinite(RH[base])]).pvalue:.3g}")
print("Paper: HGP 0.19 +- 0.17 (tolerance 0.17)")
units_df["r_hgp_15s"] = RH["15s"]
print(units_df.groupby("subject").r_hgp_15s.agg(["count", "mean", "min", "max"]).round(3).to_string())
base 15s: firing-HGP r = 0.136 +- 0.241 (n = 42, 29 positive); Wilcoxon against 0 p = 0.000618
base 1.5s: firing-HGP r = 0.139 +- 0.195 (n = 42, 32 positive); Wilcoxon against 0 p = 1.81e-05
base 1s: firing-HGP r = 0.095 +- 0.121 (n = 42, 32 positive); Wilcoxon against 0 p = 5.99e-06
Paper: HGP 0.19 +- 0.17 (tolerance 0.17)
count mean min max
subject
1 4 0.035 -0.532 0.295
2 12 0.331 -0.006 0.603
3 3 0.417 0.268 0.507
4 7 -0.012 -0.188 0.168
7 16 0.026 -0.218 0.273
Result: partially reproduced. The firing-HGP correlation is positive on average, 0.136 +- 0.241 on the 15 s time base and 0.139 +- 0.195 on the 1.5 s one (positive in 29 and 32 of 42 units, Wilcoxon against 0 p = 6e-4 and 2e-5), within the stated 0.19 +- 0.17. The muHGP and BOLD means, and the tests that HGP is lowest, cannot be computed from the dandiset.
fig4b-regression-firing-hgp¶
The per-neuron correlation of firing rate with HGP is described by both FRV and SSI in a multiple linear regression.
Dependent variable the per-neuron firing-HGP r (fig3d readings), independent FRV and SSI (fig3c readings, same time base where both exist). r read as (a) the multiple correlation R, (b) Spearman r of fitted against observed. Units with an SSI only.
rows = []
for rbase in ["15s", "1.5s"]:
for sbase in ["1s", "15s", "1s_xwire", "15s_xwire"]:
ok = np.isfinite(SSI[sbase]) & np.isfinite(RH[rbase])
X = sm.add_constant(np.column_stack([FRV[ok], SSI[sbase][ok]]))
fit = sm.OLS(RH[rbase][ok], X).fit()
rows.append({"r_base": rbase, "ssi_base": sbase, "kind": "variant" if "xwire" in sbase else "claims", "n": int(ok.sum()),
"R": np.sqrt(fit.rsquared), "model_p": f"{fit.f_pvalue:.2e}",
"spearman_fit_obs": stats.spearmanr(fit.fittedvalues, RH[rbase][ok]).statistic,
"b_FRV": fit.params[1], "p_FRV": fit.pvalues[1], "b_SSI": fit.params[2], "p_SSI": fit.pvalues[2]})
r4b = pd.DataFrame(rows)
print(r4b.round(4).to_string(index=False))
print("FRV and SSI across these units: Spearman r =", {b: round(stats.spearmanr(FRV[np.isfinite(SSI[b])], SSI[b][np.isfinite(SSI[b])]).statistic, 3) for b in SSI})
c = r4b[r4b.kind == "claims"]
print("Readings of the claims file with both coefficients p < 0.05:", int(((c.p_FRV < 0.05) & (c.p_SSI < 0.05)).sum()), "of", len(c))
print("Paper: r = 0.59, p < 1e-4, both FRV and SSI")
r_base ssi_base kind n R model_p spearman_fit_obs b_FRV p_FRV b_SSI p_SSI
15s 1s claims 41 0.5625 7.26e-04 0.5148 -0.0313 0.5170 1.5694 0.0002
15s 15s claims 41 0.6348 5.55e-05 0.6016 -0.0891 0.0722 1.0902 0.0000
15s 1s_xwire variant 41 0.5758 4.75e-04 0.5362 -0.0282 0.5538 1.5792 0.0001
15s 15s_xwire variant 41 0.6478 3.24e-05 0.6122 -0.0822 0.0870 1.0979 0.0000
1.5s 1s claims 41 0.6194 1.02e-04 0.5448 0.0168 0.6521 1.3336 0.0000
1.5s 15s claims 41 0.6848 5.99e-06 0.6357 -0.0314 0.4012 0.9166 0.0000
1.5s 1s_xwire variant 41 0.6249 8.21e-05 0.5650 0.0200 0.5882 1.3240 0.0000
1.5s 15s_xwire variant 41 0.6848 5.99e-06 0.6249 -0.0239 0.5147 0.9045 0.0000
FRV and SSI across these units: Spearman r = {'1s': np.float64(0.255), '1s_xwire': np.float64(0.243), '15s': np.float64(0.403), '15s_xwire': np.float64(0.353)}
Readings of the claims file with both coefficients p < 0.05: 0 of 4
Paper: r = 0.59, p < 1e-4, both FRV and SSI
Result: not reproduced (resting on an assumption). The regression describes the firing-HGP correlation: R = 0.56 to 0.68 over the four time-base readings, with model p from 6e-6 to 7e-4; on the 15 s time bases R = 0.635 (within 0.059 of 0.59) with model p = 5.6e-5. But SSI carries it alone (p <= 2e-4 in every reading) and FRV does not (p 0.07 to 0.65, coefficient negative in three readings), so "described by both FRV and SSI", the claim and the criterion's condition, fails in every reading. Dropping the pairs likely on one wire from SSI changes none of this. The SSI rests on the assumption of fig3c-ssi-frv-correlation.
Claims not tested here¶
- fig3b-firing-muhgp-higher-aversive: needs the continuous microwire signals; the dandiset holds only spike times and mean waveforms of the microwire data (as the data descriptor says).
- fig2-participant4-bold-correlations, table2-bold-correlations-group, table2-bold-correlation-values, bold-hgp-uncorrelated-low-mi, fig4c-regression-firing-bold, fig4d-hgp-bold-covariation: need the BOLD time courses of the clinical fMRI, which the paper names no public deposit for; the dandiset holds no fMRI.
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome |
|---|---|
| fig1b-delta-decrease | blocked (methods_underspecified): the z-scoring readings disagree |
| fig1c-hgp-increase-first-2s | not reproduced |
| other-bands-no-effect | not reproduced (assumption) |
| table2-mi-values | not reproduced |
| fig1d-delta-hgp-mi-correlation | not reproduced (assumption) |
| nsoz-delta-hgp-mi-inverse | not reproduced (assumption) |
| table2-delta-reduced-12-of-14 | partially reproduced (assumption) |
| table2-hgp-mi-positive-nsoz | not reproduced (assumption) |
| fig1e-hgp-not-delta-distinguishes-soz | not reproduced (assumption) |
| hgp-soz-sensitivity-specificity | not reproduced (assumption) |
| fig1f-hgp-more-specific-than-delta | blocked (methods_underspecified): the readings disagree |
| neuron-count | partially reproduced |
| fig3a-firing-increase-aversive | not reproduced |
| fig3a-firing-gradual-sustained | not reproduced |
| fig3b-firing-hgp-no-condition-difference | partially reproduced |
| fig3c-ssi-frv-correlation | not reproduced (assumption) |
| fig3d-firing-correlations-modalities | partially reproduced (HGP part only) |
| fig4b-regression-firing-hgp | not reproduced (assumption) |
| fig3b-firing-muhgp-higher-aversive | blocked (data_missing) |
| fig2-participant4-bold-correlations, table2-bold-correlations-group, table2-bold-correlation-values, bold-hgp-uncorrelated-low-mi, fig4c-regression-firing-bold, fig4d-hgp-bold-covariation | blocked (data_missing) |
How the analysis differs from the paper's. The paper links no code. Where its Methods leave a step
open, every reading named in the claims file was run, and the outcomes are judged on all of them;
fixed MI windows taken from the Results' text were run as a further reading, and two z-scoring
variants (trialsd, baselinepool) and the two other baselines of Methods 4.8 are shown without
deciding any outcome. The z-scoring decides most iEEG outcomes: z-scoring over the whole trial
removes the sustained delta decrease, while z-scoring against the baseline keeps it. The paper's MI
window, the group clusters, does not exist for HGP here, so the claims on the HGP MI rest on the fixed
windows, an assumption. The group test is a within-patient form of the paper's mixed model, which gives the same condition
estimate for this balanced design (t differs slightly). SSI uses every other unit of the amygdala, as the dandiset does not say which
microwire a unit was on. Epochs in which a unit has no spike at all are taken as not recorded.
What was found about the dandiset. Each session is deposited twice with identical data. The iEEG is stated in volts but holds microvolts. The 26 s epochs are laid end to end with jumps at their boundaries, and the trial table puts the first trial's start at its baseline and the others' at block onset. The units' observation intervals cover all 17 epochs although the 16 units of subject 7 have no spikes after 286 s and four units of subject 2 none in the first epoch. Two ages differ from the paper (subject 2's file gives 39, the age of the paper's patient 7), five channel locations name the hippocampus or are unknown for channels the paper places in the amygdala, no unit is assigned to a microwire, and the units are in 8 amygdalae of 5 patients, not the 9 amygdalae of the Results nor the amygdalae of Table 1 (subject 7 has units, 5 and 8 none). Re-mapping subjects 2 and 7, or 7 to the paper's 5 and 8, does not bring the MIs closer to Table 2. The BOLD time courses and the continuous microwire signals are not in the dandiset.
print(f"Notebook run time {(time.time() - started) / 60:.1f} min")
Notebook run time 4.7 min