Claar, Rembado et al. (2023) on dandiset 000458: EEG claims¶
Paper. Claar LD, Rembado I, Kuyat JR, Russo S, Marks LC, Olsen SR, Koch C (2023). Cortico-thalamo-cortical interactions modulate electrically evoked EEG responses in mice. eLife 12:RP84630. doi:10.7554/eLife.84630. The claims file was written from the version of record (full text from Europe PMC), with its figure supplement legends.
Dandiset. DANDI:000458, version 0.230317.0039. The paper's Data availability cites this version, and it is the only published version, so it is pinned. The draft lists the same 24 assets with the same total size.
What the dandiset holds. 24 sessions of 23 mice: skull EEG (30 electrodes), Neuropixels LFP and sorted units in
18 sessions, wheel speed, isoflurane epochs, and every stimulus pulse with its current, target and state. This
notebook tests the EEG claims (Figures 1, 7C, 7D and 8); claar_spikes.ipynb tests the spiking claims.
| Claim | Here |
|---|---|
| stim-depths | blocked: the stimulation depth in mm is not in the dandiset |
| eeg-good-electrodes | tested |
| fig1d-two-erp-peaks | tested |
| fig1-supp1-example | tested |
| fig1e-erp-duration | tested |
| fig1f-erp-magnitude | tested |
| fig1-supp2-first-peak-flips | blocked: the stimulation site coordinates are not in the dandiset |
| state-trial-counts | tested |
| fig7c-erp-duration-state | tested |
| fig7d-erp-magnitude-state | tested |
| fig8b-pcist-awake-anesthesia | tested |
| fig8c-pcist-depth | tested |
| fig8d-pcist-three-states | tested |
| fig8e-pcist-recovery | tested |
Readings run for every claim (claims file, header): the stimulus intensity ("high", the largest of the session's
currents, and "medium", the middle one; sessions with one current give the same trials under both), and, where the
paper does not say, all awake trials or quiet trials only. Only trials the file marks is_valid are used.
import time
import warnings
from concurrent.futures import ThreadPoolExecutor
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import claar
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
urls = claar.asset_urls()
metas = [claar.read_meta(p, u) for p, u in urls]
print(f"{len(metas)} sessions from dandiset {claar.DANDISET_ID} version {claar.DANDISET_VERSION}, metadata read in {time.time() - started:.0f} s")
24 sessions from dandiset 000458 version 0.230317.0039, metadata read in 38 s
Data check¶
rows = []
for m in metas:
el = m["estim"]
eeg_el = m["electrodes"].query("group_name == 'EEG array'")
rows.append({"subject": m["subject"], "session": m["session"], "stimulus note": m["stimulus"].replace("single pulse electrical stimuli targeted to ", ""),
"sites": ", ".join(sorted(el.site.unique())), "EEG good": int(eeg_el.is_data_valid.sum()),
"EEG hours": round(m["eeg_shape"][0] / claar.FS_EEG / 3600, 2), "probes": len(m["probes"]),
"units": 0 if m["units"] is None else len(m["units"]),
"pulses": len(el), "valid": int(el.is_valid.sum()),
"visual trials": int((m["trials"].stimulus_type == "visual").sum()),
"epochs": "; ".join(f"{r.tag} {r.start:.0f}-{r.stop:.0f} s" for r in m["epochs"].itertuples())})
info = pd.DataFrame(rows)
print(info.to_string(index=False))
print("\nsubjects:", info.subject.nunique(), " sessions:", len(info))
sites = [(i, s) for i, m in enumerate(metas) for s in sorted(m["estim"].site.unique())]
print("stimulation sites (session x target):", len(sites))
no_valid = [(metas[i]["subject"], s) for i, s in sites if not metas[i]["estim"].query("site == @s").is_valid.any()]
print("sites whose pulses are all marked invalid (left out of the analyses):", no_valid)
print(pd.Series([s for _, s in sites]).value_counts().to_string())
print("\nTrial-table columns:", list(metas[0]["trials"].columns))
print("Values of estim_target_depth:", sorted({d for m in metas for d in m["estim"].estim_target_depth.unique()}))
print("\nElectrode-table columns:", metas[0]["electrode_columns"])
eeg_xyz = pd.concat([m["electrodes"].query("group_name == 'EEG array'")[["x", "y", "z"]] for m in metas])
print("EEG electrode x, y, z over sessions: ", {c: (round(float(eeg_xyz[c].min()), 2), round(float(eeg_xyz[c].max()), 2)) for c in "xyz"},
"; distinct (x, y) positions:", len(eeg_xyz[["x", "y"]].drop_duplicates()))
probe_loc = pd.concat([m["electrodes"].query("group_name != 'EEG array'") for m in metas])
print("probe channels:", len(probe_loc), "; their location values (CCF acronyms) include:", sorted(probe_loc.location.unique())[:40])
subject session stimulus note sites EEG good EEG hours probes units pulses valid visual trials epochs
521885 20200709 MOs (superficial layers) MOs-superficial 17 0.69 0 0 360 354 0 isoflurane_induction 1309-1706 s; isoflurane_anesthesia 1706-2425 s
521886 20200716 MOs (superficial layers) MOs-superficial 20 1.17 0 0 900 893 0 isoflurane_induction 1453-1669 s; isoflurane_anesthesia 1669-2979 s
521887 20200730 MOs (superficial layers) MOs-superficial 21 1.03 0 0 800 799 0 isoflurane_induction 962-1071 s; isoflurane_anesthesia 1071-2462 s
543393 20200820 MOs (superficial layers) MOs-superficial 24 1.19 0 0 900 900 0 isoflurane_induction 1415-1684 s; isoflurane_anesthesia 1684-3035 s
543394 20200827 MOs (superficial layers) MOs-superficial 21 1.30 0 0 900 898 0 isoflurane_induction 1417-1832 s; isoflurane_anesthesia 1832-3328 s
546655 20201023 MOs (deep layers) MOs-deep 22 1.94 2 265 1159 943 300 isoflurane_induction 2006-2362 s; isoflurane_anesthesia 2362-4524 s
551397 20210211 MOs (deep layers) MOs-deep 27 2.38 2 607 1440 1440 360 isoflurane_induction 2415-2655 s; isoflurane_anesthesia 2655-4951 s
551399 20210128 MOs (deep layers) MOs-deep 22 1.91 2 593 1440 1436 0 isoflurane_induction 1787-2060 s; isoflurane_anesthesia 2060-3720 s
569062 20210218 MOs (deep layers) MOs-deep 28 2.15 3 902 1080 1079 360 isoflurane_induction 2453-2711 s; isoflurane_anesthesia 2711-5633 s
569064 20210408 MOs (deep layers) MOs-deep 30 2.82 3 481 1440 1437 480 isoflurane_induction 2477-2676 s; isoflurane_anesthesia 2676-5989 s
569068 20210304 MOs (deep layers) MOs-deep 30 2.51 3 338 1440 1435 360 isoflurane_induction 2370-2908 s; isoflurane_anesthesia 2908-5486 s
569069 20210312 MOs (deep layers) MOs-deep 22 2.42 3 520 1440 1415 360 isoflurane_induction 2341-2620 s; isoflurane_anesthesia 2620-5174 s
569070 20210401 MOs (deep layers) MOs-deep 26 2.45 0 0 1440 693 0 isoflurane_induction 1904-2194 s; isoflurane_anesthesia 2194-5610 s
569072 20210422 SSp (deep layers) SSp-deep 21 1.41 3 995 720 719 240 isoflurane_induction 2432-2655 s; isoflurane_anesthesia 2655-5051 s
569073 20210415 MOs (deep layers) MOs-deep 27 2.98 3 631 1440 1415 600 isoflurane_induction 2630-3140 s; isoflurane_anesthesia 3140-6719 s
569073 20210416 SSp (deep layers) SSp-deep 18 1.14 3 705 720 715 0 isoflurane_induction 1824-1906 s; isoflurane_anesthesia 1906-4109 s
571619 20210319 MOs (superficial and deep layers) MOs-deep, MOs-superficial 24 2.21 3 859 1440 1413 0 isoflurane_induction 3974-4107 s; isoflurane_anesthesia 4107-7953 s
571620 20210513 SSp (superficial and deep layers) SSp-deep, SSp-superficial 17 2.62 3 771 1440 1440 240 isoflurane_induction 4291-4517 s; isoflurane_anesthesia 4517-9435 s
575102 20210603 MOs (superficial layers) MOs-superficial 26 2.59 3 528 1440 1439 442 isoflurane_induction 2521-2736 s; isoflurane_anesthesia 2736-5373 s
586466 20210729 SSp (superficial and deep layers) SSp-deep, SSp-superficial 17 2.38 2 470 1440 1429 240 isoflurane_induction 4269-4480 s; isoflurane_anesthesia 4480-8548 s
586468 20210819 SSp (superficial and deep layers) SSp-deep, SSp-superficial 28 2.38 3 977 1440 720 240 isoflurane_induction 4219-4345 s; isoflurane_anesthesia 4345-8598 s
590479 20210921 SSp (superficial and deep layers) SSp-deep, SSp-superficial 18 2.06 2 616 1440 1438 0 isoflurane_induction 3489-3700 s; isoflurane_anesthesia 3700-7452 s
590480 20211111 SSp (superficial and deep layers) SSp-deep, SSp-superficial 26 2.29 3 948 1440 1437 240 isoflurane_induction 3931-4066 s; isoflurane_anesthesia 4066-8240 s
599017 20220210 SSp (superficial and deep layers) SSp-deep, SSp-superficial 21 2.40 3 943 1440 1438 240 isoflurane_induction 4106-4343 s; isoflurane_anesthesia 4343-8662 s
subjects: 23 sessions: 24
stimulation sites (session x target): 31
sites whose pulses are all marked invalid (left out of the analyses): [('586468', 'SSp-superficial')]
MOs-deep 10
SSp-deep 8
MOs-superficial 7
SSp-superficial 6
Trial-table columns: ['behavioral_epoch', 'estim_current', 'estim_target_depth', 'estim_target_region', 'is_running', 'is_valid', 'start_time', 'stimulus_description', 'stimulus_type', 'stop_time', 'trial']
Values of estim_target_depth: ['deep', 'superficial']
Electrode-table columns: ['group', 'group_name', 'id', 'is_data_valid', 'location', 'probe_horizontal_position', 'probe_vertical_position', 'reference', 'x', 'y', 'z']
EEG electrode x, y, z over sessions: {'x': (120.0, 377.0), 'y': (7.0, 64.0), 'z': (79.0, 376.0)} ; distinct (x, y) positions: 17
probe channels: 18816 ; their location values (CCF acronyms) include: ['ACAd', 'ACAv', 'AD', 'AV', 'CA1', 'CA2', 'CA3', 'CL', 'CP', 'DG-mo', 'DG-po', 'DG-sg', 'DP', 'Eth', 'HPF', 'HY', 'IAD', 'IG', 'IGL', 'ILA', 'LD', 'LGd-co', 'LGd-ip', 'LGd-sh', 'LGv', 'LP', 'LSr', 'MB', 'MD', 'MGd', 'MGm', 'MGv', 'MOp', 'MOs', 'PIL', 'PL', 'PO', 'PP', 'PoT', 'ProS']
The dandiset holds 23 subjects and 24 sessions (subject 569073 has two, MOs and SSp), matching the 23 included mice of the Methods. There are 31 stimulation sites, 13 superficial and 18 deep, the N of Figure 8B (31 sessions) and 8C (13 and 18). The trial table gives the target depth only as "superficial" or "deep"; the depth of the stimulating electrode in mm, and its coordinates, are not in the files.
# currents, stage sizes and states per site
rows = []
for i, s in sites:
el = metas[i]["estim"].query("site == @s")
lv = claar.intensity_levels(el.current)
rows.append({"subject": metas[i]["subject"], "site": s, "currents uA": "/".join(str(int(c)) for c in sorted(el.current.unique())),
"high": lv["high"], "medium": lv["medium"],
**{st: int((el.state == st).sum()) for st in ["quiet", "running", "isoflurane", "recovery1", "recovery2", "recovery3"]},
"invalid": int((~el.is_valid).sum())})
site_table = pd.DataFrame(rows)
print(site_table.to_string(index=False))
# running speed and the file's is_running flag
aw = pd.concat([m["estim"].query("behavioral_epoch == 'awake'") for m in metas])
print("\nwheel speed unit:", metas[0]["speed_unit"], "; range over sessions:",
round(min(m["speed"].min() for m in metas), 2), "to", round(max(m["speed"].max() for m in metas), 1))
print("awake pulses:", len(aw), "; is_running agrees with speed > 0 somewhere in -0.5 to +0.5 s:",
round((aw.is_running == (aw.speed_max > 0)).mean(), 4), "; disagreements:", int((aw.is_running != (aw.speed_max > 0)).sum()))
print("units' waveform_duration range (file says s):",
round(min(m["units"].waveform_duration.min() for m in metas if m["units"] is not None), 3), "to",
round(max(m["units"].waveform_duration.max() for m in metas if m["units"] is not None), 3))
subject site currents uA high medium quiet running isoflurane recovery1 recovery2 recovery3 invalid 521885 MOs-superficial 20/50/100 100.0 50.0 154 26 180 0 0 0 6 521886 MOs-superficial 20/50/100 100.0 50.0 262 38 300 300 0 0 7 521887 MOs-superficial 20 20.0 20.0 147 53 300 300 0 0 1 543393 MOs-superficial 70 70.0 70.0 246 54 300 300 0 0 0 543394 MOs-superficial 50 50.0 50.0 192 108 300 300 0 0 2 546655 MOs-deep 50/70/90 90.0 70.0 118 242 360 360 79 0 216 551397 MOs-deep 10/30/50 50.0 30.0 277 83 360 360 356 4 0 551399 MOs-deep 40/60/80 80.0 60.0 325 35 360 403 317 0 4 569062 MOs-deep 15/35/60 60.0 35.0 286 74 360 360 0 0 1 569064 MOs-deep 15/40/60 60.0 40.0 350 10 360 360 360 0 3 569068 MOs-deep 10/20/40 40.0 20.0 296 64 360 360 360 0 5 569069 MOs-deep 20/40/70 70.0 40.0 241 119 360 360 360 0 25 569070 MOs-deep 20/40/60 60.0 40.0 275 85 720 74 286 0 747 569072 SSp-deep 30/50/70 70.0 50.0 321 39 360 0 0 0 1 569073 MOs-deep 20/40/70 70.0 40.0 283 77 360 360 360 0 25 569073 SSp-deep 60/80/100 100.0 80.0 270 90 360 0 0 0 5 571619 MOs-deep 20/40/60 60.0 40.0 307 53 360 0 0 0 3 571619 MOs-superficial 20/40/60 60.0 40.0 290 70 360 0 0 0 24 571620 SSp-deep 60/80/100 100.0 80.0 312 48 360 0 0 0 0 571620 SSp-superficial 60/80/100 100.0 80.0 274 86 360 0 0 0 0 575102 MOs-superficial 50/70/90 90.0 70.0 257 103 360 360 360 0 1 586466 SSp-deep 20/40/60 60.0 40.0 320 40 360 0 0 0 4 586466 SSp-superficial 20/40/60 60.0 40.0 347 13 360 0 0 0 7 586468 SSp-deep 25/35/45 45.0 35.0 290 70 360 0 0 0 0 586468 SSp-superficial 25/35/45 45.0 35.0 236 124 360 0 0 0 720 590479 SSp-deep 40/60/80 80.0 60.0 302 58 360 0 0 0 0 590479 SSp-superficial 40/60/80 80.0 60.0 291 69 360 0 0 0 2 590480 SSp-deep 30/50/70 70.0 50.0 211 149 360 0 0 0 3 590480 SSp-superficial 30/50/70 70.0 50.0 97 263 360 0 0 0 0 599017 SSp-deep 30/50/70 70.0 50.0 59 301 360 0 0 0 1 599017 SSp-superficial 30/50/70 70.0 50.0 48 312 360 0 0 0 1 wheel speed unit: cm/s ; range over sessions: 0.0 to 88.1 awake pulses: 10640 ; is_running agrees with speed > 0 somewhere in -0.5 to +0.5 s: 0.9994 ; disagreements: 6 units' waveform_duration range (file says s): 0.096 to 2.074
Each stage has 120 pulses per current in most sessions (three currents, as in the Methods); five early superficial
sessions differ (one or three currents, 60 to 300 pulses). The recovery stages split at 30 min after the end of
the isoflurane epoch give a first and, in 11 sites, a second stage. The file's is_running agrees with the wheel
speed by the paper's rule. The wheel speed is in cm/s and never negative, as described.
# Clock check: the stimulus artifact in the raw EEG must sit at the trial onsets of the trial table.
# Contrast: raw amplitude at true onsets against onsets shifted by 1 s.
m = metas[[k for k, mm in enumerate(metas) if mm["subject"] == "571619"][0]]
f = claar.open_nwb(m["url"])
ts = f["acquisition/ElectricalSeriesEEG/timestamps"][:]
d = f["acquisition/ElectricalSeriesEEG/data"]
print("EEG timestamps: first", round(ts[0], 3), "s, median interval", round(np.median(np.diff(ts)) * 1e3, 4), "ms ->",
round(1 / np.median(np.diff(ts))), "Hz; conversion", m["eeg_conversion"], "V per unit")
on = m["estim"].start_time.values[::20]
amp = []
for shift in (0.0, 1.0):
a = []
for t in on:
k = np.searchsorted(ts, t + shift)
x = d[k - 25:k + 25, :].astype(float) * m["eeg_conversion"] * 1e6
x -= x[:20].mean(0)
a.append(np.abs(x).max(1))
amp.append(np.median(a, axis=0))
print("median over pulses and electrodes-max of |raw EEG| (uV) from -10 ms to +10 ms, true onsets:")
print(np.round(amp[0]).astype(int))
print("onsets shifted by +1 s:")
print(np.round(amp[1]).astype(int))
f.close()
EEG timestamps: first 12.534 s, median interval 0.4 ms -> 2500 Hz; conversion 1.9499999284744263e-07 V per unit
median over pulses and electrodes-max of |raw EEG| (uV) from -10 ms to +10 ms, true onsets: [ 67 58 54 50 50 46 44 41 39 36 42 35 40 42 49 49 50 61 53 59 64 68 70 92 656 696 198 127 106 107 105 104 129 112 113 118 110 113 127 115 126 127 126 129 125 135 129 121 137 135] onsets shifted by +1 s: [ 60 58 54 53 47 48 42 39 36 37 34 42 39 37 41 46 51 60 50 59 74 73 70 73 78 82 84 85 86 96 93 97 103 96 98 86 96 97 85 89 93 79 79 85 77 71 71 70 66 63]
In this session the artifact spans indices 24 and 25 of the window (656 and 696 uV), that is k - 1 and k, where k is the first EEG sample at or after the trial onset (index 25), and it is absent one second later, so the EEG and the trial table share a clock. The mask therefore starts at k - 1 (claar.MASK_START): samples k - 1 to k + 4 are replaced by k - 7 to k - 2, reversed. The same profile is printed below for every site, from the pulses of each site as it is read. The wheel speed's clock (starting at 0 s, 100 Hz) is checked by the state contrast of Figure 7, where running trials give a different ERP.
EEG of every site¶
sites = [(i, s) for i, s in sites if (metas[i]["subject"], s) not in no_valid]
def summary(args):
i, s = args
return claar.eeg_site_summary(metas[i], s)
t0 = time.time()
with ThreadPoolExecutor(max_workers=2) as ex:
res = list(ex.map(summary, sites))
print(f"{len(res)} sites summarised in {time.time() - t0:.0f} s")
for r in res:
r["depth"] = r["site"].split("-")[1]
r["region"] = r["site"].split("-")[0]
ex_site = [k for k, r in enumerate(res) if r["subject"] == "571619" and r["site"] == "MOs-deep"][0]
chk = [claar.pcist_check(r["erp"][("high", "awake")], r["times"]) for r in res]
diff = np.array([abs(a - b) for a, b in chk])
print(f"PCIst of the fast count against the original code on decimated ERPs (step 10), all {len(chk)} sites: "
f"largest absolute difference {diff.max():.2e}; values {min(min(x) for x in chk):.2f} to {max(max(x) for x in chk):.2f}")
30 sites summarised in 868 s
PCIst of the fast count against the original code on decimated ERPs (step 10), all 30 sites: largest absolute difference 0.00e+00; values 6.29 to 87.06
Artifact profile of every site (median over its pulses of the largest raw |EEG| over good electrodes, uV, relative to the mean of k - 25 to k - 11), and what remains of it in the trial-averaged ERP: the largest |z| of the z-scored GFP from -20 to -2 ms and from +2 to +20 ms (all awake trials, high current). A residue of the artifact would show as a peak just before the onset.
rows = []
for r in res:
a = r["artifact"]
z = claar.gfp_z(r["erp"][("high", "awake")], r["times"])
t = r["times"]
rows.append({"subject": r["subject"], "site": r["site"], "peak at k +": int(np.argmax(a)) - 10,
**{f"k{d:+d}": int(round(a[10 + d])) for d in (-3, -2, -1, 0, 1, 2, 5, 10)},
"max |z| -20..-2 ms": round(float(np.abs(z[(t >= -0.020) & (t <= -0.002)]).max()), 2),
"max z +2..+20 ms": round(float(z[(t >= 0.002) & (t <= 0.020)].max()), 2)})
art = pd.DataFrame(rows)
print(art.to_string(index=False))
print("\nsites whose artifact peaks at k - 1:", int((art["peak at k +"] == -1).sum()), "of", len(art),
"; at k:", int((art["peak at k +"] == 0).sum()))
fig, ax = plt.subplots(1, 2, figsize=(10, 3))
t = res[0]["times"]
w = (t >= -0.020) & (t <= 0.020)
zs = np.array([claar.gfp_z(r["erp"][("high", "awake")], t)[w] for r in res])
for k, dep in enumerate(("superficial", "deep")):
sel = [j for j, r in enumerate(res) if r["depth"] == dep]
ax[k].plot(1e3 * t[w], zs[sel].T, lw=0.6, color="0.5")
ax[k].axvspan(-0.4, 1.6, color="r", alpha=0.15)
ax[k].set(title=f"{dep}: z-scored GFP, high current, all awake", xlabel="ms", ylabel="z")
plt.tight_layout()
print("z range plotted:", round(float(zs.min()), 1), "to", round(float(zs.max()), 1), "(no points outside the axes)")
subject site peak at k + k-3 k-2 k-1 k+0 k+1 k+2 k+5 k+10 max |z| -20..-2 ms max z +2..+20 ms 521885 MOs-superficial -1 38 44 431 323 58 57 64 63 0.97 100.91 521886 MOs-superficial -1 27 34 1357 474 76 63 55 59 18.57 129.07 521887 MOs-superficial -1 34 42 2780 1179 104 88 99 98 10.86 85.77 543393 MOs-superficial -1 46 64 772 404 57 56 60 63 4.04 106.76 543394 MOs-superficial -1 50 73 1423 637 65 62 63 69 13.63 98.36 546655 MOs-deep -1 44 61 2900 2034 236 144 129 140 9.35 212.76 551397 MOs-deep 0 50 58 362 503 148 93 81 187 10.42 232.20 551399 MOs-deep -1 37 52 2043 1007 154 107 77 77 20.17 69.97 569062 MOs-deep 0 71 85 583 721 225 135 117 122 4.35 206.48 569064 MOs-deep 0 60 78 802 954 112 82 80 151 4.67 116.70 569068 MOs-deep -1 44 57 370 357 88 60 69 101 4.64 207.98 569069 MOs-deep -1 58 72 1213 694 176 124 95 103 6.14 94.75 569070 MOs-deep -1 44 63 727 630 72 65 56 72 2.53 103.96 569072 SSp-deep -1 31 40 1665 1380 87 67 61 105 4.73 57.34 569073 MOs-deep 0 38 47 365 417 79 56 57 113 1.93 80.75 569073 SSp-deep -1 35 48 2760 2577 92 70 70 76 10.26 44.15 571619 MOs-deep -1 37 49 346 318 65 51 49 56 1.55 171.34 571619 MOs-superficial -1 38 54 934 615 99 62 48 58 10.97 63.30 571620 SSp-deep -1 50 65 1226 1028 219 104 101 123 12.77 42.37 571620 SSp-superficial -1 47 69 2275 1203 111 85 84 92 16.78 28.84 575102 MOs-superficial 0 144 230 531 630 586 256 154 117 2.37 183.05 586466 SSp-deep -1 62 81 2235 909 326 175 129 165 5.74 49.00 586466 SSp-superficial -1 56 69 2160 784 311 145 83 89 4.98 9.70 586468 SSp-deep -1 41 58 1188 666 75 64 75 73 12.58 81.97 590479 SSp-deep -1 41 54 812 164 53 52 54 53 6.71 28.30 590479 SSp-superficial -1 43 56 1583 787 139 86 69 68 11.79 13.19 590480 SSp-deep -1 54 70 1427 1277 178 133 164 181 8.10 106.03 590480 SSp-superficial -1 52 82 1730 1090 161 86 70 70 21.78 24.62 599017 SSp-deep -1 61 72 721 717 129 81 96 137 3.17 67.88 599017 SSp-superficial 0 58 69 751 1228 151 83 74 81 2.81 20.89 sites whose artifact peaks at k - 1: 24 of 30 ; at k: 6 z range plotted: -1.9 to 232.2 (no points outside the axes)
The artifact is largest at k - 1 in 24 of 30 sites and at k in 6, and in every site it is far above the pre-pulse level at both samples, so the EEG of every session shares the trial table's clock and the mask covers the artifact peak. In several sites the EEG stays above its pre-pulse level after the masked samples (for example 586466 SSp-deep, 129 uV at k + 5 against 62 at k - 3; 551397 MOs-deep, 187 uV at k + 10 against 50), and the zero-phase filter spreads this step to both sides of the onset: the z-scored GFP reaches |z| = 21.78 between -20 and -2 ms in the worst site (590480 SSp-superficial) and exceeds 10 in 12 of 30 sites. After the onset this residue is not separated from the response. The -20 to 0 ms window lies inside the baselines of the GFP z-score (-2 to 0 s) and of PCIst (-0.8 to -0.002 s), so the residue raises their SD. Its effect is measured by repeating the GFP metrics, and PCIst at the medium current, with baselines that end 20 ms before the pulse (rows "-2 to -0.02 s" and matchings ending in "-excl20" below). This variant is outside the claims file's criteria and decides no outcome.
def metric_table(level, state, matching="all"):
rows = []
for r in res:
key = (level, state, matching)
if key not in r["metrics"]:
continue
rows.append({"subject": r["subject"], "site": r["site"], "depth": r["depth"], "region": r["region"],
"n": r["n"][(level, state)], **r["metrics"][key]})
return pd.DataFrame(rows)
def mean_sem(x):
x = np.asarray(x, float)
return f"{x.mean():.2f} +/- {stats.sem(x):.2f} (n={len(x)})"
eeg-good-electrodes¶
After artifact masking and visual inspection, an average of 23 +/- 4 of the 30 EEG electrodes per subject were kept as artifact-free.
good = info.groupby("subject")["EEG good"].mean()
print("good (is_data_valid) EEG electrodes per session:", info["EEG good"].tolist())
print(f"per subject (23 subjects, the two sessions of 569073 averaged): mean {good.mean():.2f}, SD {good.std():.2f}, range {good.min()}-{good.max()}")
print(f"per session (24): mean {info['EEG good'].mean():.2f}, SD {info['EEG good'].std():.2f}")
print("expected 23 +/- 4; within tolerance:", abs(good.mean() - 23) <= 4)
good (is_data_valid) EEG electrodes per session: [17, 20, 21, 24, 21, 22, 27, 22, 28, 30, 30, 22, 26, 21, 27, 18, 24, 17, 26, 17, 28, 18, 26, 21] per subject (23 subjects, the two sessions of 569073 averaged): mean 23.07, SD 4.01, range 17.0-30.0 per session (24): mean 23.04, SD 4.14 expected 23 +/- 4; within tolerance: True
Result: reproduced. The files mark good EEG electrodes (is_data_valid); their mean per subject is 23.07 (SD 4.01, range 17 to 30) against 23 +/- 4. That these marks are the paper's visual selection is consistent with the data, not stated in the files.
fig1-supp1-example¶
In the example deep MOs session, the ERP lasts 1.01 s and has a magnitude of 11.13 SD.s.
The example is "the deep MOs subject of Figure 1C and 1D (bottom left)", whose superficial MOs stimulation is shown top left "in the same subject" (Figure 1D legend). Subject 571619 is the only one with both superficial and deep MOs stimulation, so it is the example.
r = res[ex_site]
print("example:", r["subject"], r["site"])
rows = []
for lvl in claar.LEVELS:
for st in ("awake", "quiet"):
mt = r["metrics"][(lvl, st, "all")]
rows.append({"intensity": lvl, "trials": st, "n": r["n"][(lvl, st)], **{k: round(float(v), 3) for k, v in mt.items()}})
ex_tab = pd.DataFrame(rows)
print(ex_tab.to_string(index=False))
print("currents (uA):", site_table.query("subject == '571619' and site == 'MOs-deep'")[["high", "medium"]].to_dict("records"))
ok = ex_tab[(abs(ex_tab.duration_a - 1.01) <= 0.101) | (abs(ex_tab.duration_b - 1.01) <= 0.101)]
ok = ok[(abs(ok.magnitude_a - 11.13) <= 1.113) | (abs(ok.magnitude_b - 11.13) <= 1.113)]
print("\nreadings with both values within tolerance (duration 1.01 +/- 0.101 s, magnitude 11.13 +/- 1.113 SD.s):")
print(ok[["intensity", "trials", "duration_a", "duration_b", "magnitude_a", "magnitude_b"]].to_string(index=False))
fig, ax = plt.subplots(1, 2, figsize=(11, 3.2))
t = r["times"]
w = (t >= -0.3) & (t <= 1.1)
erp = r["erp"][("high", "awake")]
ax[0].plot(t[w], erp[:, w].T, color="k", lw=0.5)
ax[0].set(xlabel="time from pulse (s)", ylabel="uV", title="571619 deep MOs, awake, high current")
z = claar.gfp_z(erp, t)
ax[1].plot(t[w], z[w], color="m")
ax[1].axhline(3, ls="--", color="k")
ax[1].set(xlabel="time from pulse (s)", ylabel="GFP z")
plt.tight_layout()
print("plotted ranges: ERP", round(float(erp[:, w].min()), 1), "to", round(float(erp[:, w].max()), 1), "uV; GFP z", round(float(z[w].min()), 1), "to", round(float(z[w].max()), 1))
example: 571619 MOs-deep
intensity trials n duration_a duration_b magnitude_a magnitude_b baseline_signed_area
high awake 120 1.169 1.621 12.556 7.861 -0.0
high quiet 103 1.298 1.876 12.691 7.922 0.0
medium awake 118 1.011 1.159 11.133 7.275 0.0
medium quiet 97 1.040 1.161 9.056 5.938 0.0
currents (uA): [{'high': 60.0, 'medium': 40.0}]
readings with both values within tolerance (duration 1.01 +/- 0.101 s, magnitude 11.13 +/- 1.113 SD.s):
intensity trials duration_a duration_b magnitude_a magnitude_b
medium awake 1.011 1.159 11.133 7.275
plotted ranges: ERP -315.5 to 211.6 uV; GFP z -1.8 to 171.3
Result: reproduced. In subject 571619, deep MOs, all awake trials at the medium current (40 uA) give a duration of 1.011 s and a magnitude of 11.133 SD.s, against 1.01 s and 11.13 SD.s, with duration and magnitude read as the total time and the area above z = 3. This is the only reading with both values within tolerance; the high current gives 1.169 s and 12.556 SD.s, both above.
fig1-supp2-first-peak-flips (blocked)¶
The amplitude of the first ERP peak decays with distance from the stimulation site and eventually flips sign.
The criterion needs the distance of each electrode from the stimulation site. The EEG electrodes have coordinates (data check), the stimulating electrode has none in the files, so the claim is blocked. What can be shown without that distance is whether both signs occur at the time of the early GFP peak; this does not test the claim.
rows = []
for r in res:
erp, t = r["erp"][("high", "awake")], r["times"]
z = claar.gfp_z(erp, t)
w = np.where((t >= 0.002) & (t <= 0.050))[0]
j = w[np.argmax(z[w])]
v = erp[:, j]
rows.append({"subject": r["subject"], "site": r["site"], "peak ms": round(1e3 * t[j], 1), "electrodes": len(v),
"positive": int((v > 0).sum()), "negative": int((v < 0).sum())})
sg = pd.DataFrame(rows)
print(sg.to_string(index=False))
print("\nsites with both signs at the early GFP peak:", int(((sg.positive > 0) & (sg.negative > 0)).sum()), "of", len(sg),
"(the ERP is average-referenced, so both signs are expected at any time point)")
subject site peak ms electrodes positive negative 521885 MOs-superficial 10.8 17 10 7 521886 MOs-superficial 11.2 20 14 6 521887 MOs-superficial 8.4 21 12 9 543393 MOs-superficial 10.4 24 14 10 543394 MOs-superficial 10.0 21 13 8 546655 MOs-deep 8.8 22 10 12 551397 MOs-deep 8.0 27 12 15 551399 MOs-deep 11.6 22 13 9 569062 MOs-deep 9.2 28 16 12 569064 MOs-deep 9.6 30 16 14 569068 MOs-deep 9.2 30 16 14 569069 MOs-deep 11.2 22 14 8 569070 MOs-deep 13.2 26 14 12 569072 SSp-deep 7.2 21 11 10 569073 MOs-deep 8.8 27 13 14 569073 SSp-deep 8.4 18 7 11 571619 MOs-deep 13.6 24 12 12 571619 MOs-superficial 11.2 24 11 13 571620 SSp-deep 8.4 17 11 6 571620 SSp-superficial 6.4 17 7 10 575102 MOs-superficial 9.2 26 13 13 586466 SSp-deep 8.4 17 8 9 586466 SSp-superficial 12.4 17 7 10 586468 SSp-deep 7.6 28 15 13 590479 SSp-deep 14.0 18 9 9 590479 SSp-superficial 16.0 18 8 10 590480 SSp-deep 6.0 26 17 9 590480 SSp-superficial 8.8 26 14 12 599017 SSp-deep 6.8 21 10 11 599017 SSp-superficial 30.8 21 8 13 sites with both signs at the early GFP peak: 30 of 30 (the ERP is average-referenced, so both signs are expected at any time point)
fig1d-two-erp-peaks¶
Deep (layer 5/6) stimulation in the awake state usually evokes an ERP with two prominent peaks, near 25 ms and near 180 ms; the late component is not evident with superficial stimulation, in MOs or SSp.
Readings (claims file): (a) a local maximum of the z-scored GFP above 3 in 2 to 50 ms and another in 150 to 300 ms; (b) the same with 2 to 25 ms and 100 to 300 ms.
rows = []
for lvl in claar.LEVELS:
for st in ("awake", "quiet"):
for rd in ("a", "b"):
d = pd.DataFrame([{"depth": r["depth"], "region": r["region"], **r["peaks"][(lvl, st)][rd]} for r in res if (lvl, st) in r["peaks"]])
deep, sup = d[d.depth == "deep"], d[d.depth == "superficial"]
rows.append({"intensity": lvl, "trials": st, "reading": rd,
"deep both peaks": f"{int((deep.early & deep.late).sum())}/{len(deep)}",
"superficial no late": f"{int((~sup.late).sum())}/{len(sup)}",
"deep MOs both": f"{int((deep[deep.region=='MOs'].early & deep[deep.region=='MOs'].late).sum())}/{(deep.region=='MOs').sum()}",
"deep SSp both": f"{int((deep[deep.region=='SSp'].early & deep[deep.region=='SSp'].late).sum())}/{(deep.region=='SSp').sum()}",
"sup MOs no late": f"{int((~sup[sup.region=='MOs'].late).sum())}/{(sup.region=='MOs').sum()}",
"sup SSp no late": f"{int((~sup[sup.region=='SSp'].late).sum())}/{(sup.region=='SSp').sum()}",
"holds": ((deep.early & deep.late).mean() > 0.5) and ((~sup.late).mean() > 0.5)})
peaks_tab = pd.DataFrame(rows)
print(peaks_tab.to_string(index=False))
d = pd.DataFrame([{"subject": r["subject"], "site": r["site"], **{k: round(v, 3) if isinstance(v, float) else v for k, v in r["peaks"][("high", "awake")]["a"].items()}} for r in res])
print("\nper site, high current, all awake trials, reading (a); late_max is the largest GFP z in 150-300 ms:")
print(d.to_string(index=False))
intensity trials reading deep both peaks superficial no late deep MOs both deep SSp both sup MOs no late sup SSp no late holds
high awake a 18/18 0/12 10/10 8/8 0/7 0/5 False
high awake b 18/18 0/12 10/10 8/8 0/7 0/5 False
high quiet a 18/18 1/12 10/10 8/8 0/7 1/5 False
high quiet b 18/18 0/12 10/10 8/8 0/7 0/5 False
medium awake a 18/18 0/12 10/10 8/8 0/7 0/5 False
medium awake b 18/18 0/12 10/10 8/8 0/7 0/5 False
medium quiet a 18/18 2/12 10/10 8/8 0/7 2/5 False
medium quiet b 18/18 1/12 10/10 8/8 0/7 1/5 False
per site, high current, all awake trials, reading (a); late_max is the largest GFP z in 150-300 ms:
subject site early early_max early_t late late_max late_t
521885 MOs-superficial True 100.911 0.011 True 88.382 0.158
521886 MOs-superficial True 129.069 0.011 True 96.642 0.167
521887 MOs-superficial True 85.766 0.008 True 61.583 0.156
543393 MOs-superficial True 106.763 0.010 True 17.240 0.237
543394 MOs-superficial True 98.361 0.010 True 4.372 0.194
546655 MOs-deep True 212.762 0.009 True 85.878 0.168
551397 MOs-deep True 232.198 0.008 True 147.178 0.164
551399 MOs-deep True 69.970 0.012 True 62.833 0.183
569062 MOs-deep True 206.485 0.009 True 99.803 0.176
569064 MOs-deep True 116.696 0.010 True 89.582 0.175
569068 MOs-deep True 207.985 0.009 True 127.111 0.175
569069 MOs-deep True 94.753 0.011 True 107.174 0.189
569070 MOs-deep True 103.957 0.013 True 103.515 0.185
569072 SSp-deep True 57.341 0.007 True 42.355 0.176
569073 MOs-deep True 80.753 0.009 True 71.087 0.164
569073 SSp-deep True 44.154 0.008 True 20.616 0.170
571619 MOs-deep True 171.344 0.014 True 103.096 0.173
571619 MOs-superficial True 63.303 0.011 True 19.035 0.182
571620 SSp-deep True 42.365 0.008 True 34.234 0.178
571620 SSp-superficial True 28.845 0.006 True 6.868 0.204
575102 MOs-superficial True 183.051 0.009 True 94.100 0.180
586466 SSp-deep True 48.997 0.008 True 77.029 0.186
586466 SSp-superficial True 9.700 0.012 True 8.155 0.196
586468 SSp-deep True 81.967 0.008 True 14.045 0.160
590479 SSp-deep True 28.299 0.014 True 19.893 0.189
590479 SSp-superficial True 13.193 0.016 True 5.828 0.204
590480 SSp-deep True 106.030 0.006 True 64.165 0.202
590480 SSp-superficial True 24.623 0.009 True 4.796 0.150
599017 SSp-deep True 67.881 0.007 True 11.323 0.244
599017 SSp-superficial True 21.292 0.031 True 10.144 0.150
Result: not reproduced. Every deep site (18 of 18) has GFP peaks above z = 3 in both the early and the late window under every reading, but the late peak is also present at most superficial sites: only 0 to 2 of the 12 superficial sites lack it, under every reading. The superficial late peaks are smaller (largest late z in the per-site table) but above the threshold of the criterion.
fig1e-erp-duration¶
ERPs evoked by superficial stimulation were shorter than those evoked by deep stimulation (0.4 +/- 0.0 against 0.5 +/- 0.1 s), but the difference was not significant.
Per site, awake ERP; duration readings (a) total time above z = 3, (b) onset to the last sample above z = 3.
Claim of absence, judged on the paper's test (scipy ttest_ind, two-sided); the paper's N is 12 superficial and 18
deep, the dandiset has 13 and 18 sites.
rows = []
for lvl in claar.LEVELS:
for st, base in (("awake", "all"), ("quiet", "all"), ("awake", "excl20"), ("quiet", "excl20")):
d = metric_table(lvl, st, base)
for rd in ("duration_a", "duration_b"):
s, dp = d[d.depth == "superficial"][rd], d[d.depth == "deep"][rd]
tt = stats.ttest_ind(s, dp)
rows.append({"intensity": lvl, "trials": st, "z baseline": "-2 to 0 s" if base == "all" else "-2 to -0.02 s", "reading": rd, "superficial": mean_sem(s), "deep": mean_sem(dp),
"t": round(tt.statistic, 2), "p": round(tt.pvalue, 4),
"within tol": bool(abs(s.mean() - 0.4) <= 0.05 and abs(dp.mean() - 0.5) <= 0.1),
"p>=0.05 and sup<deep": bool(tt.pvalue >= 0.05 and s.mean() < dp.mean())})
dur_tab = pd.DataFrame(rows)
print(dur_tab.to_string(index=False))
d = metric_table("high", "awake")
print("\nper site, high current, all awake trials:")
print(d[["subject", "site", "n", "duration_a", "duration_b", "magnitude_a", "magnitude_b"]].round(3).to_string(index=False))
intensity trials z baseline reading superficial deep t p within tol p>=0.05 and sup<deep
high awake -2 to 0 s duration_a 0.56 +/- 0.14 (n=12) 0.65 +/- 0.07 (n=18) -0.63 0.5323 False True
high awake -2 to 0 s duration_b 1.44 +/- 0.17 (n=12) 1.57 +/- 0.11 (n=18) -0.68 0.4999 False True
high quiet -2 to 0 s duration_a 0.53 +/- 0.15 (n=12) 0.64 +/- 0.06 (n=18) -0.75 0.4599 False True
high quiet -2 to 0 s duration_b 1.28 +/- 0.21 (n=12) 1.69 +/- 0.09 (n=18) -2.06 0.0488 False False
high awake -2 to -0.02 s duration_a 0.61 +/- 0.14 (n=12) 0.69 +/- 0.07 (n=18) -0.58 0.5678 False True
high awake -2 to -0.02 s duration_b 1.74 +/- 0.08 (n=12) 1.70 +/- 0.08 (n=18) 0.32 0.7521 False False
high quiet -2 to -0.02 s duration_a 0.56 +/- 0.15 (n=12) 0.68 +/- 0.06 (n=18) -0.79 0.4374 False True
high quiet -2 to -0.02 s duration_b 1.37 +/- 0.20 (n=12) 1.77 +/- 0.06 (n=18) -2.27 0.0314 False False
medium awake -2 to 0 s duration_a 0.33 +/- 0.04 (n=12) 0.49 +/- 0.06 (n=18) -1.81 0.0813 False True
medium awake -2 to 0 s duration_b 1.23 +/- 0.13 (n=12) 1.36 +/- 0.15 (n=18) -0.67 0.5086 False True
medium quiet -2 to 0 s duration_a 0.33 +/- 0.06 (n=12) 0.52 +/- 0.06 (n=18) -2.03 0.0522 False True
medium quiet -2 to 0 s duration_b 1.45 +/- 0.13 (n=12) 1.35 +/- 0.14 (n=18) 0.54 0.5931 False False
medium awake -2 to -0.02 s duration_a 0.37 +/- 0.04 (n=12) 0.50 +/- 0.06 (n=18) -1.48 0.1510 True True
medium awake -2 to -0.02 s duration_b 1.48 +/- 0.11 (n=12) 1.37 +/- 0.15 (n=18) 0.59 0.5625 False False
medium quiet -2 to -0.02 s duration_a 0.37 +/- 0.06 (n=12) 0.53 +/- 0.06 (n=18) -1.74 0.0924 True True
medium quiet -2 to -0.02 s duration_b 1.57 +/- 0.12 (n=12) 1.43 +/- 0.13 (n=18) 0.78 0.4442 False False
per site, high current, all awake trials:
subject site n duration_a duration_b magnitude_a magnitude_b
521885 MOs-superficial 60 1.924 1.975 17.364 10.938
521886 MOs-superficial 100 0.402 0.629 9.707 8.722
521887 MOs-superficial 199 0.397 0.814 8.366 6.260
543393 MOs-superficial 300 0.550 1.998 4.839 2.983
543394 MOs-superficial 300 0.188 1.551 2.457 1.818
546655 MOs-deep 119 0.470 0.470 12.301 8.050
551397 MOs-deep 120 0.546 1.606 15.851 10.770
551399 MOs-deep 120 0.513 0.578 6.788 7.897
569062 MOs-deep 120 0.508 1.782 14.939 9.603
569064 MOs-deep 120 0.902 1.150 9.457 6.420
569068 MOs-deep 120 1.032 1.524 12.923 7.981
569069 MOs-deep 118 0.898 1.765 14.048 9.173
569070 MOs-deep 109 0.468 1.956 12.105 7.802
569072 SSp-deep 120 0.586 1.417 6.755 4.402
569073 MOs-deep 113 0.244 1.975 4.000 2.496
569073 SSp-deep 118 0.634 1.965 4.525 3.009
571619 MOs-deep 120 1.169 1.621 12.556 7.861
571619 MOs-superficial 115 0.475 1.546 4.015 2.861
571620 SSp-deep 120 0.535 1.924 4.902 3.468
571620 SSp-superficial 120 0.396 1.667 1.561 1.305
575102 MOs-superficial 120 1.157 2.000 14.403 9.163
586466 SSp-deep 120 0.958 1.860 14.272 9.170
586466 SSp-superficial 120 0.322 1.715 1.033 0.694
586468 SSp-deep 120 0.328 1.784 2.764 1.943
590479 SSp-deep 120 0.438 1.781 2.603 1.707
590479 SSp-superficial 120 0.314 1.750 1.081 0.735
590480 SSp-deep 119 1.048 1.993 13.852 9.632
590480 SSp-superficial 120 0.202 0.202 1.145 1.519
599017 SSp-deep 119 0.363 1.125 3.012 1.990
599017 SSp-superficial 119 0.359 1.437 2.452 1.534
Result: partially reproduced. Superficial ERPs are shorter than deep ERPs and the difference is not significant under every reading of duration (a). At the medium current with all awake trials, the reading that matches the example of fig1-supp1, superficial ERPs last 0.33 +/- 0.04 s and deep ERPs 0.49 +/- 0.06 s (n = 12 and 18), p = 0.081 (paper p = 0.0843); the deep mean is within tolerance, the superficial mean is 0.02 s below 0.4 +/- 0.05. At the high current both means are above tolerance (0.56 and 0.65 s, p = 0.53). No reading of the criterion puts both means within tolerance. With the baseline ending 20 ms before the pulse (the artifact check above), the medium current gives 0.37 +/- 0.04 and 0.50 +/- 0.06 s, p = 0.151 (quiet trials 0.37 and 0.53 s, p = 0.092), both within tolerance; the outcome on the criterion's -2 to 0 s baseline therefore turns on the artifact residue in the last 20 ms of the baseline.
fig1f-erp-magnitude¶
ERPs evoked by superficial stimulation had a significantly smaller normalized magnitude than those evoked by deep stimulation (3.3 +/- 0.8 against 6.2 +/- 0.9 SD.s).
Readings (a) the area above z = 3 (SD.s); (b) that area divided by the area under the baseline. The area under the signed baseline z-score is zero by construction (printed as baseline_signed_area), so (b) cannot be computed as written. The table also shows (b) with the area under |z| in the baseline as divisor, for comparison only: it is close to a rescaling of (a), and the paper prints the magnitude in SD.s, the unit of (a). The outcome rests on (a).
rows = []
for lvl in claar.LEVELS:
for st, base in (("awake", "all"), ("quiet", "all"), ("awake", "excl20"), ("quiet", "excl20")):
d = metric_table(lvl, st, base)
for rd in ("magnitude_a", "magnitude_b"):
s, dp = d[d.depth == "superficial"][rd], d[d.depth == "deep"][rd]
tt = stats.ttest_ind(s, dp)
rows.append({"intensity": lvl, "trials": st, "z baseline": "-2 to 0 s" if base == "all" else "-2 to -0.02 s", "reading": rd, "superficial": mean_sem(s), "deep": mean_sem(dp),
"t": round(tt.statistic, 2), "p": round(tt.pvalue, 4),
"within tol": bool(abs(s.mean() - 3.3) <= 0.8 and abs(dp.mean() - 6.2) <= 0.9),
"p<0.05 and sup<deep": bool(tt.pvalue < 0.05 and s.mean() < dp.mean())})
mag_tab = pd.DataFrame(rows)
print(mag_tab.to_string(index=False))
print("max |baseline_signed_area| over sites:", float(np.max(np.abs(metric_table('high', 'awake').baseline_signed_area))))
intensity trials z baseline reading superficial deep t p within tol p<0.05 and sup<deep
high awake -2 to 0 s magnitude_a 5.70 +/- 1.61 (n=12) 9.31 +/- 1.13 (n=18) -1.89 0.0687 False False
high awake -2 to 0 s magnitude_b 4.04 +/- 1.07 (n=12) 6.30 +/- 0.74 (n=18) -1.80 0.0833 True False
high quiet -2 to 0 s magnitude_a 4.89 +/- 1.52 (n=12) 8.81 +/- 1.13 (n=18) -2.12 0.0432 False True
high quiet -2 to 0 s magnitude_b 3.41 +/- 1.03 (n=12) 5.92 +/- 0.73 (n=18) -2.05 0.0496 True True
high awake -2 to -0.02 s magnitude_a 6.75 +/- 1.69 (n=12) 10.43 +/- 1.24 (n=18) -1.80 0.0832 False False
high awake -2 to -0.02 s magnitude_b 4.36 +/- 1.09 (n=12) 6.66 +/- 0.79 (n=18) -1.76 0.0896 False False
high quiet -2 to -0.02 s magnitude_a 5.47 +/- 1.62 (n=12) 9.69 +/- 1.21 (n=18) -2.13 0.0420 False True
high quiet -2 to -0.02 s magnitude_b 3.58 +/- 1.06 (n=12) 6.22 +/- 0.77 (n=18) -2.06 0.0486 True True
medium awake -2 to 0 s magnitude_a 3.03 +/- 0.78 (n=12) 5.86 +/- 0.90 (n=18) -2.21 0.0352 True True
medium awake -2 to 0 s magnitude_b 2.04 +/- 0.54 (n=12) 3.81 +/- 0.59 (n=18) -2.10 0.0453 False True
medium quiet -2 to 0 s magnitude_a 2.85 +/- 0.80 (n=12) 5.76 +/- 0.87 (n=18) -2.32 0.0278 True True
medium quiet -2 to 0 s magnitude_b 1.90 +/- 0.54 (n=12) 3.71 +/- 0.56 (n=18) -2.22 0.0346 False True
medium awake -2 to -0.02 s magnitude_a 3.38 +/- 0.85 (n=12) 6.22 +/- 0.95 (n=18) -2.08 0.0466 True True
medium awake -2 to -0.02 s magnitude_b 2.17 +/- 0.56 (n=12) 3.95 +/- 0.61 (n=18) -2.02 0.0528 False False
medium quiet -2 to -0.02 s magnitude_a 3.09 +/- 0.86 (n=12) 6.02 +/- 0.91 (n=18) -2.22 0.0347 True True
medium quiet -2 to -0.02 s magnitude_b 2.00 +/- 0.56 (n=12) 3.82 +/- 0.58 (n=18) -2.15 0.0400 False True
max |baseline_signed_area| over sites: 4.150390395807335e-07
Result: reproduced. With reading (a), the area above z = 3 in SD.s, the medium current gives a superficial magnitude of 3.03 +/- 0.78 and a deep magnitude of 5.86 +/- 0.90 with all awake trials, both within tolerance, t-test p = 0.035 (paper 3.3 +/- 0.8, 6.2 +/- 0.9, p = 0.0354); quiet trials give 2.85 and 5.76, p = 0.028. At the high current the means are above tolerance (5.70 and 9.31, p = 0.069 with all awake trials). With the baseline ending 20 ms before the pulse (variant, decides nothing) the medium current gives 3.38 and 6.22, p = 0.047, also within tolerance.
state-trial-counts¶
Per animal, the state comparisons averaged a median of 97 quiet-wakefulness, 21 running and 120 anesthetized trials, across 17 animals.
Session sets (claims file): (a) the deep sites, one value per animal (the two deep sites of 569073 averaged); (b) every site, one value per animal (sites of one animal averaged).
def per_animal(df, cols):
return df.groupby("subject")[cols].mean()
rows = []
for lvl in claar.LEVELS:
d = pd.DataFrame([{"subject": r["subject"], "depth": r["depth"], "quiet": r["n"][(lvl, "quiet")],
"running": r["n"][(lvl, "running")], "isoflurane": r["n"][(lvl, "isoflurane")]} for r in res])
for name, sub in (("a: deep sites", d[d.depth == "deep"]), ("b: all sites", d)):
pa = per_animal(sub, ["quiet", "running", "isoflurane"])
q = pa.quantile([0.25, 0.5, 0.75])
rows.append({"intensity": lvl, "set": name, "animals": len(pa),
**{c: f"{q.loc[0.5, c]:.1f} [{q.loc[0.25, c]:.1f}-{q.loc[0.75, c]:.1f}]" for c in pa.columns},
"within tol": bool(abs(q.loc[0.5, "quiet"] - 97) <= 9.7 and abs(q.loc[0.5, "running"] - 21) <= 2.1 and abs(q.loc[0.5, "isoflurane"] - 120) <= 12)})
counts_tab = pd.DataFrame(rows)
print(counts_tab.to_string(index=False))
print("\nvalid trials only; the paper's medians are 97 [86-104], 21 [16-27], 120 [120-120] over 17 animals")
intensity set animals quiet running isoflurane within tol
high a: deep sites 17 99.0 [90.0-104.0] 21.0 [16.0-30.0] 120.0 [119.0-120.0] True
high b: all sites 23 96.0 [88.5-105.0] 24.0 [17.5-45.5] 120.0 [119.0-120.0] False
medium a: deep sites 17 97.0 [86.0-104.0] 21.0 [16.0-27.0] 120.0 [120.0-120.0] True
medium b: all sites 23 96.0 [84.5-108.5] 24.0 [15.5-45.0] 120.0 [119.2-120.0] False
valid trials only; the paper's medians are 97 [86-104], 21 [16-27], 120 [120-120] over 17 animals
Result: reproduced. Over the 17 animals with deep stimulation, the median counts per animal are 99, 21 and 120 at the high current and 97, 21 and 120 at the medium current, against 97, 21 and 120. Over all 23 animals the running median is 24, out of tolerance.
fig7c-erp-duration-state and fig7d-erp-magnitude-state¶
ERP duration decreased from quiet wakefulness (0.5 +/- 0.1 s) to running (0.3 +/- 0.0 s) and anesthesia (0.3 +/- 0.0 s); the late component is diminished when running and absent under anesthesia.
ERP magnitude decreased from quiet wakefulness (6.1 +/- 0.9) to running (1.6 +/- 0.3) and anesthesia (1.6 +/- 0.3).
Per animal (session sets as above), states unmatched or matched by random subsampling to the smallest state;
Friedman test (scipy friedmanchisquare).
exp = {"duration": {"quiet": (0.5, 0.1), "running": (0.3, 0.05), "isoflurane": (0.3, 0.05)},
"magnitude": {"quiet": (6.1, 0.9), "running": (1.6, 0.3), "isoflurane": (1.6, 0.3)}}
state_rows = {"duration": [], "magnitude": []}
for lvl in claar.LEVELS:
for matching in ("all", "matched", "excl20"):
tabs = []
for st in ("quiet", "running", "isoflurane"):
d = metric_table(lvl, st, matching)
d["state"] = st
tabs.append(d)
d = pd.concat(tabs)
for setname, sub in (("a: deep", d[d.depth == "deep"]), ("b: all", d)):
for measure in ("duration", "magnitude"):
for rd in ("a", "b"):
col = f"{measure}_{rd}"
w = sub.pivot_table(index="subject", columns="state", values=col, aggfunc="mean").dropna()
fr = stats.friedmanchisquare(w.quiet, w.running, w.isoflurane)
within = all(abs(w[s].mean() - exp[measure][s][0]) <= exp[measure][s][1] for s in exp[measure])
top = (w.quiet.mean() > w.running.mean()) and (w.quiet.mean() > w.isoflurane.mean())
state_rows[measure].append({"intensity": lvl, "matching": matching, "set": setname, "reading": rd, "animals": len(w),
**{s: mean_sem(w[s]) for s in ("quiet", "running", "isoflurane")},
"Q": round(fr.statistic, 1), "p": f"{fr.pvalue:.2g}", "within tol": within,
"quiet largest": top, "all hold": bool(within and fr.pvalue < 0.05 and top)})
dstate = pd.DataFrame(state_rows["duration"])
print(dstate.to_string(index=False))
intensity matching set reading animals quiet running isoflurane Q p within tol quiet largest all hold
high all a: deep a 17 0.65 +/- 0.07 (n=17) 0.35 +/- 0.04 (n=17) 0.33 +/- 0.05 (n=17) 18.4 0.0001 False True False
high all a: deep b 17 1.68 +/- 0.09 (n=17) 1.27 +/- 0.12 (n=17) 1.19 +/- 0.14 (n=17) 7.9 0.019 False True False
high all b: all a 23 0.63 +/- 0.08 (n=23) 0.32 +/- 0.03 (n=23) 0.38 +/- 0.08 (n=23) 23.6 7.6e-06 False True False
high all b: all b 23 1.53 +/- 0.10 (n=23) 1.20 +/- 0.09 (n=23) 1.13 +/- 0.12 (n=23) 10.8 0.0046 False True False
high matched a: deep a 17 0.39 +/- 0.04 (n=17) 0.34 +/- 0.04 (n=17) 0.23 +/- 0.05 (n=17) 12.8 0.0016 False True False
high matched a: deep b 17 1.43 +/- 0.11 (n=17) 1.22 +/- 0.11 (n=17) 1.04 +/- 0.17 (n=17) 4.2 0.12 False True False
high matched b: all a 23 0.34 +/- 0.03 (n=23) 0.31 +/- 0.03 (n=23) 0.23 +/- 0.03 (n=23) 15.2 0.0005 False True False
high matched b: all b 23 1.25 +/- 0.11 (n=23) 1.21 +/- 0.09 (n=23) 1.19 +/- 0.12 (n=23) 0.6 0.74 False True False
high excl20 a: deep a 17 0.68 +/- 0.07 (n=17) 0.36 +/- 0.04 (n=17) 0.38 +/- 0.05 (n=17) 16.4 0.00028 False True False
high excl20 a: deep b 17 1.76 +/- 0.06 (n=17) 1.42 +/- 0.11 (n=17) 1.25 +/- 0.13 (n=17) 6.7 0.035 False True False
high excl20 b: all a 23 0.66 +/- 0.08 (n=23) 0.33 +/- 0.03 (n=23) 0.43 +/- 0.08 (n=23) 23.6 7.6e-06 False True False
high excl20 b: all b 23 1.63 +/- 0.09 (n=23) 1.36 +/- 0.09 (n=23) 1.25 +/- 0.11 (n=23) 5.6 0.062 False True False
medium all a: deep a 17 0.52 +/- 0.07 (n=17) 0.25 +/- 0.04 (n=17) 0.22 +/- 0.04 (n=17) 20.1 4.3e-05 False True False
medium all a: deep b 17 1.31 +/- 0.14 (n=17) 1.47 +/- 0.17 (n=17) 1.11 +/- 0.17 (n=17) 7.4 0.025 False False False
medium all b: all a 23 0.46 +/- 0.05 (n=23) 0.24 +/- 0.03 (n=23) 0.25 +/- 0.05 (n=23) 23.6 7.6e-06 False True False
medium all b: all b 23 1.34 +/- 0.11 (n=23) 1.37 +/- 0.14 (n=23) 1.00 +/- 0.14 (n=23) 11.0 0.004 False False False
medium matched a: deep a 17 0.33 +/- 0.04 (n=17) 0.24 +/- 0.04 (n=17) 0.11 +/- 0.03 (n=17) 13.1 0.0015 False True False
medium matched a: deep b 17 1.32 +/- 0.13 (n=17) 1.58 +/- 0.15 (n=17) 1.21 +/- 0.16 (n=17) 4.6 0.1 False False False
medium matched b: all a 23 0.31 +/- 0.04 (n=23) 0.24 +/- 0.03 (n=23) 0.12 +/- 0.03 (n=23) 16.4 0.00027 False True False
medium matched b: all b 23 1.35 +/- 0.11 (n=23) 1.40 +/- 0.13 (n=23) 1.08 +/- 0.15 (n=23) 3.1 0.21 False False False
medium excl20 a: deep a 17 0.53 +/- 0.07 (n=17) 0.26 +/- 0.04 (n=17) 0.24 +/- 0.04 (n=17) 19.9 4.8e-05 False True False
medium excl20 a: deep b 17 1.40 +/- 0.13 (n=17) 1.50 +/- 0.17 (n=17) 1.25 +/- 0.15 (n=17) 5.1 0.08 False False False
medium excl20 b: all a 23 0.48 +/- 0.05 (n=23) 0.24 +/- 0.03 (n=23) 0.26 +/- 0.05 (n=23) 23.7 7e-06 False True False
medium excl20 b: all b 23 1.45 +/- 0.11 (n=23) 1.44 +/- 0.13 (n=23) 1.22 +/- 0.12 (n=23) 8.1 0.018 False True False
Result: partially reproduced. Duration falls from quiet wakefulness to running and isoflurane in every reading except reading (b) at the medium current, and the Friedman test is significant in most (17 deep animals, medium current, unmatched, reading a: 0.52, 0.25 and 0.22 s, Q = 20.1, p = 4.3e-05; paper Q = 15.6). No reading puts the three means within tolerance: at the medium current the isoflurane mean (0.22 s) is below 0.3 +/- 0.05, at the high current the quiet mean (0.65 s) is above 0.5 +/- 0.1, and matched trials shorten the quiet mean to 0.31 to 0.39 s. With the baseline ending 20 ms before the pulse (variant) the medium current gives 0.53, 0.26 and 0.24 s, the isoflurane mean still below tolerance.
mstate = pd.DataFrame(state_rows["magnitude"])
print(mstate.to_string(index=False))
intensity matching set reading animals quiet running isoflurane Q p within tol quiet largest all hold
high all a: deep a 17 8.96 +/- 1.17 (n=17) 3.40 +/- 0.72 (n=17) 2.44 +/- 0.43 (n=17) 21.299999 2.4e-05 False True False
high all a: deep b 17 6.03 +/- 0.75 (n=17) 2.17 +/- 0.47 (n=17) 1.69 +/- 0.30 (n=17) 17.799999 0.00014 False True False
high all b: all a 23 7.96 +/- 0.97 (n=23) 2.85 +/- 0.49 (n=23) 2.92 +/- 0.60 (n=23) 26.400000 1.8e-06 False True False
high all b: all b 23 5.44 +/- 0.64 (n=23) 1.83 +/- 0.32 (n=23) 2.03 +/- 0.38 (n=23) 23.700001 7e-06 False True False
high matched a: deep a 17 3.94 +/- 0.87 (n=17) 3.08 +/- 0.65 (n=17) 0.89 +/- 0.19 (n=17) 18.400000 0.0001 False True False
high matched a: deep b 17 2.57 +/- 0.57 (n=17) 1.95 +/- 0.42 (n=17) 0.59 +/- 0.13 (n=17) 15.600000 0.0004 False True False
high matched b: all a 23 3.22 +/- 0.55 (n=23) 2.61 +/- 0.42 (n=23) 0.94 +/- 0.15 (n=23) 22.500000 1.3e-05 False True False
high matched b: all b 23 2.13 +/- 0.36 (n=23) 1.66 +/- 0.27 (n=23) 0.62 +/- 0.10 (n=23) 19.900000 4.7e-05 False True False
high excl20 a: deep a 17 9.88 +/- 1.25 (n=17) 3.52 +/- 0.76 (n=17) 2.83 +/- 0.49 (n=17) 20.600000 3.4e-05 False True False
high excl20 a: deep b 17 6.34 +/- 0.80 (n=17) 2.24 +/- 0.49 (n=17) 1.85 +/- 0.34 (n=17) 17.400000 0.00017 False True False
high excl20 b: all a 23 8.86 +/- 1.05 (n=23) 2.99 +/- 0.52 (n=23) 3.42 +/- 0.61 (n=23) 26.400000 1.8e-06 False True False
high excl20 b: all b 23 5.73 +/- 0.68 (n=23) 1.90 +/- 0.34 (n=23) 2.21 +/- 0.40 (n=23) 24.600000 4.5e-06 False True False
medium all a: deep a 17 5.89 +/- 0.91 (n=17) 1.51 +/- 0.32 (n=17) 1.23 +/- 0.29 (n=17) 21.299999 2.4e-05 False True False
medium all a: deep b 17 3.80 +/- 0.59 (n=17) 0.98 +/- 0.21 (n=17) 0.81 +/- 0.18 (n=17) 20.600000 3.4e-05 False True False
medium all b: all a 23 4.92 +/- 0.55 (n=23) 1.45 +/- 0.26 (n=23) 1.45 +/- 0.26 (n=23) 23.600000 7.6e-06 False True False
medium all b: all b 23 3.21 +/- 0.36 (n=23) 0.93 +/- 0.17 (n=23) 0.98 +/- 0.16 (n=23) 23.700001 7e-06 False True False
medium matched a: deep a 17 2.31 +/- 0.49 (n=17) 1.35 +/- 0.23 (n=17) 0.33 +/- 0.11 (n=17) 23.600000 7.3e-06 False True False
medium matched a: deep b 17 1.48 +/- 0.31 (n=17) 0.87 +/- 0.15 (n=17) 0.21 +/- 0.07 (n=17) 24.500000 4.9e-06 False True False
medium matched b: all a 23 1.98 +/- 0.32 (n=23) 1.32 +/- 0.21 (n=23) 0.44 +/- 0.11 (n=23) 34.700001 2.9e-08 False True False
medium matched b: all b 23 1.28 +/- 0.21 (n=23) 0.85 +/- 0.14 (n=23) 0.28 +/- 0.07 (n=23) 34.700001 2.9e-08 False True False
medium excl20 a: deep a 17 6.17 +/- 0.96 (n=17) 1.56 +/- 0.32 (n=17) 1.37 +/- 0.31 (n=17) 20.600000 3.4e-05 True True True
medium excl20 a: deep b 17 3.91 +/- 0.60 (n=17) 1.01 +/- 0.21 (n=17) 0.87 +/- 0.19 (n=17) 20.600000 3.4e-05 False True False
medium excl20 b: all a 23 5.21 +/- 0.58 (n=23) 1.49 +/- 0.27 (n=23) 1.67 +/- 0.28 (n=23) 23.700001 7e-06 True True True
medium excl20 b: all b 23 3.32 +/- 0.37 (n=23) 0.96 +/- 0.18 (n=23) 1.05 +/- 0.17 (n=23) 23.700001 7e-06 False True False
Result: partially reproduced. Magnitude falls from quiet wakefulness to running and isoflurane in every reading, Friedman p <= 0.0004 and quiet the largest. With reading (a) at the medium current, unmatched, 17 deep animals, the means are 5.89, 1.51 and 1.23 (Q = 21.3, p = 2.4e-05): quiet and running are within tolerance, isoflurane is 0.07 below 1.6 +/- 0.3. At the high current quiet (8.96) and running (3.40) are above tolerance. Reading (b), with the area under |z| as divisor, is shown but does not decide the outcome (fig1f). With the baseline ending 20 ms before the pulse (variant, outside the criterion) the medium current gives 6.17, 1.56 and 1.37 (Q = 20.6, p = 3.4e-05), all within tolerance; the outcome on the criterion's baseline therefore turns on the artifact residue in the last 20 ms of the baseline.
fig8b-pcist-awake-anesthesia¶
PCIst was higher in the awake than in the isoflurane-anesthetized state in all subjects (45.5 +/- 3.6 against 16.5 +/- 2.4).
PCIst of the code the paper cites, on the ERP of the good electrodes at 2.5 kHz, trials matched within site between awake (quiet and running) and anesthetized, by random subsampling or the first trials.
def pc_frame(level, fig, how):
rows = []
for r in res:
d = {s: v for (lv, fg, s, hw), v in r["pcist"].items() if lv == level and fg == fig and hw == how}
if d:
rows.append({"subject": r["subject"], "site": r["site"], "depth": r["depth"], **d})
return pd.DataFrame(rows)
rows = []
for lvl in claar.LEVELS:
for how in ("random", "first") + (("random-excl20", "first-excl20") if lvl == "medium" else ()):
d = pc_frame(lvl, "fig8b", how)
tt = stats.ttest_rel(d.awake, d.isoflurane)
rows.append({"intensity": lvl, "matching": how, "sites": len(d), "awake": mean_sem(d.awake), "isoflurane": mean_sem(d.isoflurane),
"awake > iso": f"{int((d.awake > d.isoflurane).sum())}/{len(d)}", "t": round(tt.statistic, 2), "p": f"{tt.pvalue:.3g}",
"within tol": bool(abs(d.awake.mean() - 45.5) <= 3.6 and abs(d.isoflurane.mean() - 16.5) <= 2.4)})
pc8b = pd.DataFrame(rows)
print(pc8b.to_string(index=False))
d = pc_frame("high", "fig8b", "random").merge(pc_frame("medium", "fig8b", "random"), on=["subject", "site", "depth"], suffixes=("_high", "_medium"))
print("\nper site, random matching:")
print(d.round(1).to_string(index=False))
intensity matching sites awake isoflurane awake > iso t p within tol
high random 30 44.33 +/- 3.39 (n=30) 23.55 +/- 2.59 (n=30) 29/30 6.69 2.48e-07 False
high first 30 44.24 +/- 3.41 (n=30) 23.42 +/- 2.55 (n=30) 29/30 6.61 3.03e-07 False
medium random 30 43.26 +/- 3.72 (n=30) 14.93 +/- 2.41 (n=30) 30/30 9.83 9.77e-11 True
medium first 30 43.47 +/- 3.76 (n=30) 14.84 +/- 2.43 (n=30) 30/30 9.90 8.22e-11 True
medium random-excl20 30 44.59 +/- 3.73 (n=30) 15.87 +/- 2.30 (n=30) 30/30 9.76 1.14e-10 True
medium first-excl20 30 44.85 +/- 3.78 (n=30) 15.75 +/- 2.32 (n=30) 30/30 9.83 9.67e-11 True
per site, random matching:
subject site depth awake_high isoflurane_high awake_medium isoflurane_medium
521885 MOs-superficial superficial 55.1 26.2 65.8 18.8
521886 MOs-superficial superficial 43.1 18.1 37.2 7.8
521887 MOs-superficial superficial 51.5 8.7 51.5 8.7
543393 MOs-superficial superficial 53.0 7.6 53.0 7.6
543394 MOs-superficial superficial 25.7 13.1 25.7 13.1
546655 MOs-deep deep 58.6 34.2 53.6 14.7
551397 MOs-deep deep 57.7 34.7 76.9 26.2
551399 MOs-deep deep 45.1 25.2 38.8 15.3
569062 MOs-deep deep 65.6 20.1 57.1 11.7
569064 MOs-deep deep 55.9 46.6 58.7 12.6
569068 MOs-deep deep 35.1 31.1 22.2 3.9
569069 MOs-deep deep 52.6 40.4 56.1 9.0
569070 MOs-deep deep 34.9 22.2 42.9 28.5
569072 SSp-deep deep 44.0 27.2 41.1 8.8
569073 MOs-deep deep 31.5 56.1 33.7 15.8
569073 SSp-deep deep 49.8 10.7 35.0 5.7
571619 MOs-deep deep 53.0 26.0 50.3 18.6
571619 MOs-superficial superficial 31.9 25.3 24.9 14.4
571620 SSp-deep deep 51.0 25.5 77.2 18.0
571620 SSp-superficial superficial 18.0 9.6 21.6 6.3
575102 MOs-superficial superficial 48.8 33.2 44.1 30.0
586466 SSp-deep deep 68.8 12.9 70.8 21.2
586466 SSp-superficial superficial 20.7 7.2 20.9 17.7
586468 SSp-deep deep 34.6 28.3 24.8 7.0
590479 SSp-deep deep 49.9 7.3 37.3 7.2
590479 SSp-superficial superficial 8.3 4.6 8.6 7.5
590480 SSp-deep deep 97.2 58.1 92.5 73.0
590480 SSp-superficial superficial 25.2 12.0 32.2 15.1
599017 SSp-deep deep 53.1 25.0 32.5 3.4
599017 SSp-superficial superficial 10.2 9.7 10.5 0.0
Result: reproduced. At the medium current with random matching, PCIst is 43.26 +/- 3.72 awake and 14.93 +/- 2.41 under isoflurane, both within tolerance of 45.5 +/- 3.6 and 16.5 +/- 2.4, higher awake in 30 of 30 sites, paired t-test p = 9.8e-11; first-trial matching gives 43.47 and 14.84. At the high current the anesthetized mean (23.4 to 23.6) is above tolerance and one site of 30 is higher under isoflurane. The dandiset has 30 sites with valid trials, the paper 31. With the PCIst baseline ending 20 ms before the pulse (variant) the medium current gives 44.59 and 15.87 (random) and 44.85 and 15.75 (first), within tolerance, 30 of 30 sites higher awake.
fig8c-pcist-depth¶
In the awake state PCIst was significantly lower for superficial than for deep stimulation (35.6 +/- 4.8 against 52.7 +/- 4.4).
rows = []
for lvl in claar.LEVELS:
for how, fig in (("all", "fig8c"), ("random", "fig8b"), ("first", "fig8b")):
d = pc_frame(lvl, fig, how)
s, dp = d[d.depth == "superficial"].awake, d[d.depth == "deep"].awake
tt = stats.ttest_ind(s, dp)
rows.append({"intensity": lvl, "trials": "all awake" if how == "all" else f"matched ({how})", "superficial": mean_sem(s), "deep": mean_sem(dp),
"t": round(tt.statistic, 2), "p": f"{tt.pvalue:.3g}",
"within tol": bool(abs(s.mean() - 35.6) <= 4.8 and abs(dp.mean() - 52.7) <= 4.4),
"p<0.05, sup lower": bool(tt.pvalue < 0.05 and s.mean() < dp.mean())})
pc8c = pd.DataFrame(rows)
print(pc8c.to_string(index=False))
intensity trials superficial deep t p within tol p<0.05, sup lower
high all awake 32.34 +/- 4.84 (n=12) 52.51 +/- 3.69 (n=18) -3.36 0.00225 True True
high matched (random) 32.61 +/- 4.92 (n=12) 52.15 +/- 3.63 (n=18) -3.26 0.00291 True True
high matched (first) 32.10 +/- 4.76 (n=12) 52.33 +/- 3.68 (n=18) -3.40 0.00205 True True
medium all awake 32.90 +/- 5.05 (n=12) 51.24 +/- 4.79 (n=18) -2.55 0.0164 True True
medium matched (random) 33.01 +/- 5.11 (n=12) 50.09 +/- 4.61 (n=18) -2.43 0.0217 True True
medium matched (first) 33.10 +/- 5.17 (n=12) 50.38 +/- 4.65 (n=18) -2.44 0.0215 True True
Result: reproduced. Under every reading superficial PCIst is lower than deep, within tolerance of 35.6 +/- 4.8 and 52.7 +/- 4.4: 32.34 +/- 4.84 against 52.51 +/- 3.69 at the high current with all awake trials (p = 0.0023), and 32.90 against 51.24 at the medium current (p = 0.016); n = 12 and 18.
fig8d-pcist-three-states¶
PCIst differed between all states; highest in quiet wakefulness (55.2 +/- 10.3), lower when running (39.9 +/- 7.7), and lowest under anesthesia (14.0 +/- 5.8), with every subject higher when running than under anesthesia.
Sites with at least 30 running trials at the intensity used; three states matched by random subsampling (the
criterion's method; the first trials are shown as a variant). One-way RM ANOVA (statsmodels AnovaRM), post hoc
paired t-tests with Benjamini-Hochberg.
from statsmodels.stats.anova import AnovaRM
rows = []
for lvl in claar.LEVELS:
for how in ("random", "first") + (("random-excl20", "first-excl20") if lvl == "medium" else ()):
for per in ("site", "animal"):
d = pc_frame(lvl, "fig8d", how)
if per == "animal":
d = d.groupby("subject")[["quiet", "running", "isoflurane"]].mean().reset_index().assign(site="mean")
if len(d) < 3:
rows.append({"intensity": lvl, "matching": how, "unit": per, "n": len(d)})
continue
d["unit"] = d.subject + d.site
long = d.melt(id_vars=["unit"], value_vars=["quiet", "running", "isoflurane"], var_name="state", value_name="pcist")
an = AnovaRM(long, "pcist", "unit", within=["state"]).fit().anova_table
pairs = [("quiet", "running"), ("quiet", "isoflurane"), ("running", "isoflurane")]
pp = [stats.ttest_rel(d[a], d[b]).pvalue for a, b in pairs]
pc_ = claar.bh(pp)
rows.append({"intensity": lvl, "matching": how, "unit": per, "n": len(d), "animals": d.subject.nunique(),
**{s: mean_sem(d[s]) for s in ("quiet", "running", "isoflurane")},
"F": f"F[{int(an['Num DF'].iloc[0])},{int(an['Den DF'].iloc[0])}]={an['F Value'].iloc[0]:.1f}", "p": f"{an['Pr > F'].iloc[0]:.3g}",
"post hoc q (q-r, q-i, r-i)": ", ".join(f"{x:.3g}" for x in pc_),
"running > iso": f"{int((d.running > d.isoflurane).sum())}/{len(d)}",
"within tol": bool(abs(d.quiet.mean() - 55.2) <= 10.3 and abs(d.running.mean() - 39.9) <= 7.7 and abs(d.isoflurane.mean() - 14.0) <= 5.8)})
pc8d = pd.DataFrame(rows)
print(pc8d.to_string(index=False))
print("\nsites with >= 30 running trials, high current:")
print(pc_frame("high", "fig8d", "random").round(1).to_string(index=False))
print("\nmatched trials per state and site (the smallest of quiet, running, isoflurane; the paper's average is 51):")
for lvl in claar.LEVELS:
nm = [r["n_matched"][(lvl, "fig8d")] for r in res if (lvl, "fig8d") in r["n_matched"]]
print(f" {lvl}: {nm}, mean {np.mean(nm):.1f}")
intensity matching unit n animals quiet running isoflurane F p post hoc q (q-r, q-i, r-i) running > iso within tol
high random site 13 11 45.65 +/- 8.73 (n=13) 38.69 +/- 7.33 (n=13) 15.01 +/- 4.21 (n=13) F[2,24]=14.9 6.12e-05 0.239, 0.00105, 0.00105 13/13 True
high random animal 11 11 45.30 +/- 6.87 (n=11) 39.19 +/- 6.94 (n=11) 15.61 +/- 4.03 (n=11) F[2,20]=13.9 0.000164 0.336, 0.0022, 0.0022 11/11 True
high first site 13 11 45.43 +/- 8.26 (n=13) 39.31 +/- 7.22 (n=13) 19.85 +/- 6.58 (n=13) F[2,24]=10.9 0.000421 0.225, 0.0054, 0.0054 12/13 False
high first animal 11 11 45.05 +/- 6.00 (n=11) 39.21 +/- 6.58 (n=11) 19.17 +/- 6.70 (n=11) F[2,20]=10.6 0.000719 0.227, 0.00843, 0.00843 10/11 True
medium random site 10 8 44.65 +/- 10.14 (n=10) 32.11 +/- 8.35 (n=10) 9.69 +/- 3.36 (n=10) F[2,18]=12.9 0.000331 0.0572, 0.00744, 0.0105 9/10 False
medium random animal 8 8 45.54 +/- 8.41 (n=8) 33.06 +/- 7.53 (n=8) 10.03 +/- 3.31 (n=8) F[2,14]=13.0 0.000642 0.0949, 0.00899, 0.0156 8/8 True
medium first site 10 8 44.93 +/- 10.03 (n=10) 32.77 +/- 7.39 (n=10) 10.78 +/- 5.15 (n=10) F[2,18]=23.2 1.05e-05 0.0209, 0.000934, 0.00177 10/10 True
medium first animal 8 8 45.47 +/- 7.30 (n=8) 32.94 +/- 6.20 (n=8) 9.86 +/- 3.69 (n=8) F[2,14]=26.0 1.94e-05 0.0235, 0.000777, 0.00538 8/8 True
medium random-excl20 site 10 8 45.19 +/- 10.26 (n=10) 32.38 +/- 8.36 (n=10) 10.20 +/- 3.37 (n=10) F[2,18]=12.4 0.000404 0.0585, 0.00836, 0.0116 9/10 True
medium random-excl20 animal 8 8 46.10 +/- 8.54 (n=8) 33.35 +/- 7.51 (n=8) 10.60 +/- 3.36 (n=8) F[2,14]=12.5 0.000762 0.0982, 0.00984, 0.0167 8/8 True
medium first-excl20 site 10 8 45.48 +/- 10.20 (n=10) 33.23 +/- 7.46 (n=10) 10.05 +/- 4.50 (n=10) F[2,18]=21.3 1.81e-05 0.0265, 0.00143, 0.00199 10/10 True
medium first-excl20 animal 8 8 45.98 +/- 7.42 (n=8) 33.42 +/- 6.35 (n=8) 9.44 +/- 3.43 (n=8) F[2,14]=25.2 2.31e-05 0.0301, 0.00082, 0.00536 8/8 True
sites with >= 30 running trials, high current:
subject site depth quiet running isoflurane
521887 MOs-superficial superficial 75.3 32.8 10.0
543393 MOs-superficial superficial 40.8 15.1 0.0
543394 MOs-superficial superficial 9.2 25.3 12.5
546655 MOs-deep deep 52.9 72.8 30.0
569062 MOs-deep deep 57.0 72.4 8.6
569069 MOs-deep deep 58.4 49.4 17.3
569073 SSp-deep deep 29.6 10.9 10.3
571620 SSp-superficial superficial 16.9 16.8 14.1
575102 MOs-superficial superficial 63.2 63.6 45.7
590480 SSp-deep deep 123.2 82.7 41.9
590480 SSp-superficial superficial 21.0 10.9 4.8
599017 SSp-deep deep 34.3 36.3 0.0
599017 SSp-superficial superficial 11.7 14.0 0.0
matched trials per state and site (the smallest of quiet, running, isoflurane; the paper's average is 51):
high: [53, 54, 108, 38, 30, 38, 37, 32, 30, 52, 30, 19, 16], mean 41.3
medium: [53, 54, 108, 40, 37, 37, 49, 35, 22, 18], mean 45.3
Result: partially reproduced. With random matching (the criterion's method) the state effect is significant and the order is quiet > running > isoflurane in every row. At the high current, per site, the three means are within tolerance (45.65, 38.69, 15.01; F[2,24] = 14.9, p = 6e-05; running above isoflurane in 13 of 13), but quiet against running is not significant after correction (q = 0.24); per animal (11) the same holds (q = 0.34). At the medium current the quiet against running pair has q = 0.057 per site and 0.095 per animal. Only first-trial matching at the medium current makes every pair significant with the means within tolerance (q = 0.021, 0.0009, 0.0018 per site; 8 animals, the paper's 8). Matched trials average 41 per state at the high current and 45 at the medium current, against the paper's 51. With the PCIst baseline ending 20 ms before the pulse (variant) at the medium current, random matching puts the three means within tolerance per site (45.19, 32.38, 10.20) with quiet against running q = 0.059 (0.098 per animal), and first-trial matching keeps every pair significant (q = 0.027, 0.0014, 0.0020 per site).
fig8e-pcist-recovery¶
After isoflurane was stopped, PCIst rose in all but one subject within the first 30 min (32.7 +/- 4.9 against 14.3 +/- 2.3 anesthetized and 48.0 +/- 4.0 awake) and did not increase further in the next 30 min (36.4 +/- 5.0).
Readings for the sample: (a) every site with a first recovery stage (means over them; the RM ANOVA on those with all four stages); (b) only sites with all four stages. Trials matched across the stages of a site.
rows = []
for lvl in claar.LEVELS:
for how in ("random", "first") + (("random-excl20", "first-excl20") if lvl == "medium" else ()):
d = pc_frame(lvl, "fig8e", how)
full = d.dropna(subset=["recovery2"]) if "recovery2" in d else d.iloc[:0]
full = full.assign(unit=full.subject + full.site)
long = full.melt(id_vars=["unit"], value_vars=["awake", "isoflurane", "recovery1", "recovery2"], var_name="state", value_name="pcist")
an = AnovaRM(long, "pcist", "unit", within=["state"]).fit().anova_table if len(full) >= 3 else None
pairs = [("awake", "isoflurane"), ("awake", "recovery1"), ("awake", "recovery2"), ("isoflurane", "recovery1"), ("isoflurane", "recovery2"), ("recovery1", "recovery2")]
q = claar.bh([stats.ttest_rel(full[a], full[b]).pvalue for a, b in pairs]) if len(full) >= 3 else [np.nan] * 6
for sample, dd in (("a: all with recovery1", d), ("b: all four stages", full)):
rows.append({"intensity": lvl, "matching": how, "sample": sample, "sites": len(dd),
**{s: mean_sem(dd[s].dropna()) for s in ("awake", "isoflurane", "recovery1", "recovery2")},
"rec1 > iso": f"{int((dd.recovery1 > dd.isoflurane).sum())}/{len(dd)}",
"RM ANOVA (four stages)": "" if an is None else f"F[{int(an['Num DF'].iloc[0])},{int(an['Den DF'].iloc[0])}]={an['F Value'].iloc[0]:.1f}, p={an['Pr > F'].iloc[0]:.3g}",
"q rec1-rec2": round(q[5], 3), "q awake-rec1": round(q[1], 3), "q iso-rec1": round(q[3], 3),
"within tol": bool(abs(dd.awake.mean() - 48.0) <= 4.0 and abs(dd.isoflurane.mean() - 14.3) <= 2.3 and abs(dd.recovery1.mean() - 32.7) <= 4.9 and abs(dd.recovery2.dropna().mean() - 36.4) <= 5.0)})
pc8e = pd.DataFrame(rows)
print(pc8e.to_string(index=False))
print("\nper site, high current, random matching:")
print(pc_frame("high", "fig8e", "random").round(1).to_string(index=False))
print("\nper site, high current, first trials matched:")
print(pc_frame("high", "fig8e", "first").round(1).to_string(index=False))
print("\nper site, medium current, random matching:")
print(pc_frame("medium", "fig8e", "random").round(1).to_string(index=False))
print("\nper site, medium current, first trials matched:")
print(pc_frame("medium", "fig8e", "first").round(1).to_string(index=False))
for lvl in claar.LEVELS:
print(f"\nmatched trials per stage and site (fig8e), {lvl} current:",
{f"{r['subject']} {r['site']}": r["n_matched"][(lvl, "fig8e")] for r in res if (lvl, "fig8e") in r["n_matched"]})
# a minimum of matched trials per stage (a check; the paper states none)
rows = []
for lvl in claar.LEVELS:
ok = {(r["subject"], r["site"]) for r in res if r["n_matched"].get((lvl, "fig8e"), 0) >= 20}
for how in ("random", "first"):
d = pc_frame(lvl, "fig8e", how)
d = d[[(a, b) in ok for a, b in zip(d.subject, d.site)]]
full = d.dropna(subset=["recovery2"]) if "recovery2" in d else d.iloc[:0]
for sample, dd in (("a: all with recovery1", d), ("b: all four stages", full)):
rows.append({"intensity": lvl, "matching": how, "sample": sample, "sites": len(dd),
**{s: mean_sem(dd[s].dropna()) for s in ("awake", "isoflurane", "recovery1", "recovery2")},
"rec1 > iso": f"{int((dd.recovery1 > dd.isoflurane).sum())}/{len(dd)}",
"within tol": bool(len(dd) and abs(dd.awake.mean() - 48.0) <= 4.0 and abs(dd.isoflurane.mean() - 14.3) <= 2.3 and abs(dd.recovery1.mean() - 32.7) <= 4.9 and abs(dd.recovery2.dropna().mean() - 36.4) <= 5.0)})
print("\nsites with at least 20 matched trials per stage only:")
print(pd.DataFrame(rows).to_string(index=False))
print("\nERPs under isoflurane, high current, every site with a first recovery stage: largest |z| of the GFP from -20 to -2 ms, peak z 2 to 50 ms and 100 to 800 ms, against awake:")
rows = []
for r in res:
if ("high", "fig8e") not in r["n_matched"]:
continue
t = r["times"]
row = {"subject": r["subject"], "site": r["site"]}
for st in ("awake", "isoflurane"):
z = claar.gfp_z(r["erp"][("high", st)], t)
row[f"{st} n"] = r["n"][("high", st)]
row[f"{st} pre"] = round(float(np.abs(z[(t >= -0.020) & (t <= -0.002)]).max()), 1)
row[f"{st} early"] = round(float(z[(t >= 0.002) & (t <= 0.050)].max()), 1)
row[f"{st} late"] = round(float(z[(t >= 0.100) & (t <= 0.800)].max()), 1)
rows.append(row)
print(pd.DataFrame(rows).to_string(index=False))
intensity matching sample sites awake isoflurane recovery1 recovery2 rec1 > iso RM ANOVA (four stages) q rec1-rec2 q awake-rec1 q iso-rec1 within tol
high random a: all with recovery1 13 43.58 +/- 4.77 (n=13) 25.04 +/- 4.19 (n=13) 34.69 +/- 5.86 (n=13) 44.74 +/- 5.91 (n=8) 9/13 F[3,21]=1.8, p=0.185 0.279 0.446 0.446 False
high random b: all four stages 8 40.98 +/- 6.73 (n=8) 31.72 +/- 5.54 (n=8) 35.93 +/- 6.43 (n=8) 44.74 +/- 5.91 (n=8) 6/8 F[3,21]=1.8, p=0.185 0.279 0.446 0.446 False
high first a: all with recovery1 13 44.63 +/- 3.72 (n=13) 49.67 +/- 21.78 (n=13) 33.26 +/- 5.66 (n=13) 44.76 +/- 5.93 (n=8) 7/13 F[3,21]=0.8, p=0.501 0.188 0.188 0.543 False
high first b: all four stages 8 42.81 +/- 4.63 (n=8) 72.26 +/- 33.62 (n=8) 33.65 +/- 5.99 (n=8) 44.76 +/- 5.93 (n=8) 4/8 F[3,21]=0.8, p=0.501 0.188 0.188 0.543 False
medium random a: all with recovery1 13 45.73 +/- 4.24 (n=13) 16.10 +/- 2.52 (n=13) 29.14 +/- 4.46 (n=13) 33.82 +/- 4.46 (n=8) 8/13 F[3,21]=9.7, p=0.000322 0.604 0.025 0.093 True
medium random b: all four stages 8 46.24 +/- 6.14 (n=8) 19.48 +/- 3.59 (n=8) 31.72 +/- 4.29 (n=8) 33.82 +/- 4.46 (n=8) 5/8 F[3,21]=9.7, p=0.000322 0.604 0.025 0.093 False
medium first a: all with recovery1 13 45.80 +/- 4.23 (n=13) 12.88 +/- 2.28 (n=13) 28.51 +/- 4.59 (n=13) 34.06 +/- 4.49 (n=8) 11/13 F[3,21]=16.9, p=8.32e-06 0.475 0.014 0.013 True
medium first b: all four stages 8 46.47 +/- 6.08 (n=8) 14.89 +/- 3.53 (n=8) 30.74 +/- 4.76 (n=8) 34.06 +/- 4.49 (n=8) 8/8 F[3,21]=16.9, p=8.32e-06 0.475 0.014 0.013 True
medium random-excl20 a: all with recovery1 13 47.01 +/- 4.20 (n=13) 17.47 +/- 3.01 (n=13) 30.61 +/- 4.49 (n=13) 36.31 +/- 5.10 (n=8) 8/13 F[3,21]=7.4, p=0.0015 0.593 0.027 0.121 False
medium random-excl20 b: all four stages 8 47.32 +/- 6.17 (n=8) 21.08 +/- 4.41 (n=8) 33.22 +/- 3.94 (n=8) 36.31 +/- 5.10 (n=8) 5/8 F[3,21]=7.4, p=0.0015 0.593 0.027 0.121 False
medium first-excl20 a: all with recovery1 13 47.02 +/- 4.20 (n=13) 15.16 +/- 3.03 (n=13) 29.79 +/- 4.71 (n=13) 36.52 +/- 5.12 (n=8) 10/13 F[3,21]=9.3, p=0.000423 0.490 0.030 0.042 True
medium first-excl20 b: all four stages 8 47.44 +/- 6.11 (n=8) 18.11 +/- 4.64 (n=8) 31.87 +/- 4.68 (n=8) 36.52 +/- 5.12 (n=8) 7/8 F[3,21]=9.3, p=0.000423 0.490 0.030 0.042 False
per site, high current, random matching:
subject site depth awake isoflurane recovery1 recovery2
521886 MOs-superficial superficial 43.0 17.5 30.6 NaN
521887 MOs-superficial superficial 51.5 13.1 53.4 NaN
543393 MOs-superficial superficial 53.0 7.6 4.3 NaN
543394 MOs-superficial superficial 25.6 13.0 8.3 NaN
546655 MOs-deep deep 0.0 1.8 4.2 41.7
551397 MOs-deep deep 57.3 35.4 54.9 70.4
551399 MOs-deep deep 43.7 22.5 25.9 43.6
569062 MOs-deep deep 65.7 20.5 66.9 NaN
569064 MOs-deep deep 56.7 46.6 57.8 59.0
569068 MOs-deep deep 36.1 31.1 48.9 45.1
569069 MOs-deep deep 53.8 29.8 22.1 13.5
569073 MOs-deep deep 31.5 54.2 39.2 49.4
575102 MOs-superficial superficial 48.7 32.3 34.5 35.1
per site, high current, first trials matched:
subject site depth awake isoflurane recovery1 recovery2
521886 MOs-superficial superficial 42.9 18.0 30.6 NaN
521887 MOs-superficial superficial 51.5 8.4 53.0 NaN
543393 MOs-superficial superficial 53.0 7.6 4.3 NaN
543394 MOs-superficial superficial 25.0 13.1 8.3 NaN
546655 MOs-deep deep 20.3 306.6 4.2 41.7
551397 MOs-deep deep 58.8 33.9 54.3 70.4
551399 MOs-deep deep 44.4 33.2 26.2 43.6
569062 MOs-deep deep 65.3 20.5 66.9 NaN
569064 MOs-deep deep 55.5 46.6 56.9 59.7
569068 MOs-deep deep 34.1 31.1 31.7 45.4
569069 MOs-deep deep 49.8 38.0 23.7 13.5
569073 MOs-deep deep 31.5 56.2 37.3 48.5
575102 MOs-superficial superficial 48.1 32.5 34.9 35.1
per site, medium current, random matching:
subject site depth awake isoflurane recovery1 recovery2
521886 MOs-superficial superficial 37.3 8.0 15.7 NaN
521887 MOs-superficial superficial 51.5 13.1 53.4 NaN
543393 MOs-superficial superficial 53.0 7.6 4.3 NaN
543394 MOs-superficial superficial 25.6 13.0 8.3 NaN
546655 MOs-deep deep 37.6 26.9 21.9 35.3
551397 MOs-deep deep 76.9 26.2 44.2 40.5
551399 MOs-deep deep 38.3 18.2 16.3 32.6
569062 MOs-deep deep 57.1 11.7 43.4 NaN
569064 MOs-deep deep 58.9 12.6 51.5 50.5
569068 MOs-deep deep 22.3 3.9 26.2 7.7
569069 MOs-deep deep 58.1 8.5 35.3 34.0
569073 MOs-deep deep 33.7 29.5 22.5 27.9
575102 MOs-superficial superficial 44.1 30.0 35.9 42.1
per site, medium current, first trials matched:
subject site depth awake isoflurane recovery1 recovery2
521886 MOs-superficial superficial 37.0 7.5 15.7 NaN
521887 MOs-superficial superficial 51.5 8.4 53.0 NaN
543393 MOs-superficial superficial 53.0 7.6 4.3 NaN
543394 MOs-superficial superficial 25.0 13.1 8.3 NaN
546655 MOs-deep deep 43.1 2.1 10.1 35.3
551397 MOs-deep deep 76.9 26.2 44.2 40.5
551399 MOs-deep deep 36.9 10.6 20.3 32.6
569062 MOs-deep deep 57.1 11.7 43.4 NaN
569064 MOs-deep deep 60.1 12.6 51.3 51.5
569068 MOs-deep deep 21.8 3.9 25.4 7.9
569069 MOs-deep deep 55.2 13.3 35.3 34.0
569073 MOs-deep deep 33.7 20.4 23.4 28.6
575102 MOs-superficial superficial 44.1 30.0 35.9 42.1
matched trials per stage and site (fig8e), high current: {'521886 MOs-superficial': 99, '521887 MOs-superficial': 199, '543393 MOs-superficial': 300, '543394 MOs-superficial': 298, '546655 MOs-deep': 1, '551397 MOs-deep': 117, '551399 MOs-deep': 104, '569062 MOs-deep': 119, '569064 MOs-deep': 119, '569068 MOs-deep': 118, '569069 MOs-deep': 109, '569073 MOs-deep': 113, '575102 MOs-superficial': 119}
matched trials per stage and site (fig8e), medium current: {'521886 MOs-superficial': 98, '521887 MOs-superficial': 199, '543393 MOs-superficial': 300, '543394 MOs-superficial': 298, '546655 MOs-deep': 23, '551397 MOs-deep': 120, '551399 MOs-deep': 106, '569062 MOs-deep': 120, '569064 MOs-deep': 118, '569068 MOs-deep': 119, '569069 MOs-deep': 118, '569073 MOs-deep': 111, '575102 MOs-superficial': 120}
sites with at least 20 matched trials per stage only:
intensity matching sample sites awake isoflurane recovery1 recovery2 rec1 > iso within tol
high random a: all with recovery1 12 47.22 +/- 3.37 (n=12) 26.98 +/- 4.05 (n=12) 37.23 +/- 5.74 (n=12) 45.17 +/- 6.81 (n=7) 8/12 False
high random b: all four stages 7 46.84 +/- 3.83 (n=7) 35.99 +/- 4.08 (n=7) 40.46 +/- 5.27 (n=7) 45.17 +/- 6.81 (n=7) 5/7 False
high first a: all with recovery1 12 46.65 +/- 3.39 (n=12) 28.26 +/- 4.36 (n=12) 35.68 +/- 5.56 (n=12) 45.19 +/- 6.83 (n=7) 7/12 False
high first b: all four stages 7 46.02 +/- 3.86 (n=7) 38.78 +/- 3.52 (n=7) 37.86 +/- 4.92 (n=7) 45.19 +/- 6.83 (n=7) 4/7 False
medium random a: all with recovery1 13 45.73 +/- 4.24 (n=13) 16.10 +/- 2.52 (n=13) 29.14 +/- 4.46 (n=13) 33.82 +/- 4.46 (n=8) 8/13 True
medium random b: all four stages 8 46.24 +/- 6.14 (n=8) 19.48 +/- 3.59 (n=8) 31.72 +/- 4.29 (n=8) 33.82 +/- 4.46 (n=8) 5/8 False
medium first a: all with recovery1 13 45.80 +/- 4.23 (n=13) 12.88 +/- 2.28 (n=13) 28.51 +/- 4.59 (n=13) 34.06 +/- 4.49 (n=8) 11/13 True
medium first b: all four stages 8 46.47 +/- 6.08 (n=8) 14.89 +/- 3.53 (n=8) 30.74 +/- 4.76 (n=8) 34.06 +/- 4.49 (n=8) 8/8 True
ERPs under isoflurane, high current, every site with a first recovery stage: largest |z| of the GFP from -20 to -2 ms, peak z 2 to 50 ms and 100 to 800 ms, against awake:
subject site awake n awake pre awake early awake late isoflurane n isoflurane pre isoflurane early isoflurane late
521886 MOs-superficial 100 18.6 129.1 96.6 100 19.9 83.3 7.0
521887 MOs-superficial 199 10.9 85.8 61.6 300 15.4 51.9 6.2
543393 MOs-superficial 300 4.0 106.8 23.2 300 1.7 43.1 3.3
543394 MOs-superficial 300 13.6 98.4 8.9 300 15.1 118.1 5.9
546655 MOs-deep 119 9.3 212.8 85.9 75 2.9 60.3 26.8
551397 MOs-deep 120 10.4 232.2 147.2 120 10.9 173.4 19.3
551399 MOs-deep 120 20.2 70.0 62.8 120 22.3 24.3 7.4
569062 MOs-deep 120 4.3 206.5 99.8 120 6.0 123.5 6.0
569064 MOs-deep 120 4.7 116.7 89.6 119 2.8 95.1 8.8
569068 MOs-deep 120 4.6 208.0 127.1 118 1.1 94.3 6.0
569069 MOs-deep 118 6.1 94.8 107.2 119 3.3 29.7 13.1
569073 MOs-deep 113 1.9 80.8 71.1 120 2.1 53.9 3.8
575102 MOs-superficial 120 2.4 183.1 94.1 120 1.8 76.1 17.7
Result: reproduced. At the medium current with the first trials matched, in the 8 sites with all four stages: awake 46.47, isoflurane 14.89, first recovery 30.74 and second recovery 34.06, all within tolerance; RM ANOVA F[3,21] = 16.9, p = 8.3e-06; first recovery above isoflurane in 8 of 8; recovery 2 against recovery 1 q = 0.475. In the 13 sites with a first recovery stage the means are also within tolerance but 2 of 13 sites do not rise. Random matching at the medium current puts the means of the 13 sites within tolerance, but the first recovery exceeds isoflurane in only 8 of 13. At the medium current every site has at least 23 matched trials per stage (546655 MOs-deep the fewest), so requiring 20 changes nothing; at the high current 546655 MOs-deep is matched on a single trial (isoflurane PCIst 306.6 with first trials), which drives the isoflurane mean of 49.67 +/- 21.78, and without it the high current still puts isoflurane above tolerance (26.98 to 38.78). The outcome rests on one combination of readings. With the PCIst baseline ending 20 ms before the pulse (variant, outside the criterion) that combination gives an isoflurane mean of 18.11, above 14.3 +/- 2.3, while the 13 sites with a first recovery stage (first trials) stay within tolerance (47.02, 15.16, 29.79, 36.52) with 10 of 13 rising.
Summary¶
| Claim | Outcome |
|---|---|
| stim-depths | blocked: data missing (no stimulation depth in mm) |
| eeg-good-electrodes | reproduced |
| fig1d-two-erp-peaks | not reproduced |
| fig1-supp1-example | reproduced (medium current) |
| fig1e-erp-duration | partially reproduced |
| fig1f-erp-magnitude | reproduced (medium current) |
| fig1-supp2-first-peak-flips | blocked: data missing (no stimulation coordinates) |
| state-trial-counts | reproduced |
| fig7c-erp-duration-state | partially reproduced |
| fig7d-erp-magnitude-state | partially reproduced |
| fig8b-pcist-awake-anesthesia | reproduced (medium current) |
| fig8c-pcist-depth | reproduced |
| fig8d-pcist-three-states | partially reproduced |
| fig8e-pcist-recovery | reproduced (one reading) |
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.
How the analysis differs from the paper's: the paper does not say which of the three currents it analysed, so the
highest and the middle current are both run; the middle current with all awake trials gives the paper's example
(1.011 s and 11.133 SD.s against 1.01 and 11.13) and the paper's p-values for Figures 1E and 1F (0.081 and 0.035
against 0.0843 and 0.0354), which supports it as the paper's choice. Good EEG electrodes are the file's
is_data_valid marks; trials the file marks invalid are left out, which removes one superficial site (586468, SSp)
entirely; the stimulus artifact peaks one sample before the first sample at or after the onset, so the 2 ms mask
starts one sample earlier than the onset sample; the 0.1 to 100 Hz filter is zero phase; the magnitude's division by
"the area under the baseline" cannot be computed (the signed area of a z-scored baseline is zero), so the magnitude
is the area above z = 3; PCIst uses the cited code with an exact, memory-saving count in its threshold search
(checked against the original above); recovery stages are split at 30 and 60 min after the end of the isoflurane
epoch.
The zero-phase filter spreads a residue of the stimulus artifact into the last 20 ms before the pulse (|z| of the GFP above 10 in 12 of 30 sites), inside the z-score and PCIst baselines. With baselines that end 20 ms before the pulse (a variant outside the criteria), fig1e-erp-duration and fig7d-erp-magnitude-state would have all means within tolerance at the medium current, and the one passing reading of fig8e-pcist-recovery would not; fig1f-erp-magnitude, fig8b-pcist-awake-anesthesia and fig8d-pcist-three-states would keep their results. The outcomes stay on the criteria's baselines.
What was found about the dataset: 23 subjects and 31 stimulation sites, as in the paper; the trial table gives the current, target, depth label, state and validity of every pulse, and the EEG electrode table marks good electrodes (mean 23 per subject, as in the paper); the stimulation depth in mm and the coordinates of the stimulating electrode are not in the files; all 720 pulses of one site are marked invalid without saying why.
print(f"{(time.time() - started) / 60:.1f} min")
15.7 min