Reproducing Hristova, Fasol et al. (2025) from dandiset 001888¶
Paper. Hristova K, Fasol MCM, McLaughlin N, Nawaz MS, Taskiran M, Buller-Peralta I, Harris AP, Sutherland A, Bassi A, Ocampo-Garces A, Escudero J, Kind PC, Gonzalez-Sulser A (2025). Absence seizures and sleep-wake abnormalities in a rat model of GRIN2B neurodevelopmental disorder. Epilepsia 66:4996-5013. doi:10.1111/epi.18606. This attempt used the open-access full text (CC BY, PMC12779324) and its Supporting Information. After the first verification round, the criterion of figs5b-transition-percent was corrected to the Supplementary Methods' rule for SWDs initiated in NREM, and variants asked for by the verifier were added.
Dandiset. DANDI:001888, version draft, as
modified on 2026-07-31 09:37:44 UTC. The paper cites no DANDI version (its data statement says raw
data are available on request), and the dandiset has no published version, so the draft is the
only version there is. A draft can change: grin2b.check_version() stops the notebook if the
draft's modification time or file count differs from the state read here.
What the dandiset holds. One NWB file per rat for 17 rats: the full 73 to 99 h wireless recording (14 skull EEG channels, 2 neck EMG channels, 250.4 Hz), the two 24-h baseline windows (BL1, BL2) in the epochs table, and per window the lab's automated outputs: sleep scores in 5-s epochs, SWD events, an SWD flag per epoch, per-hour SWD totals and per-state spectra. The western blots, development measures, visual scoring and pharmacology recordings of the paper are not in it.
Analysis code. The text of the paper and its Supporting Information was read from Europe PMC (PMC12779324). The Supplementary Methods judge the state of an SWD in wake and REM on the epoch before it, and SWD initiation in NREM and wake-NREM transitions on the six epochs before it.
Animals. The files give 14 Grin2b+/- (Het), 2 wild type (WT) and one rat without genotype or sex (GRIN2B-132). The paper analysed 15 Het and 9 WT, so every genotype comparison here rests on two WT rats. No data from outside the dandiset are used.
Code read. The lab's sleep scorer (Gonzalez-Sulser-Team/AUTOMATIC-SLEEP-SCORER, commit
95b25fb; R scorer, MATLAB channel selection and cepstral SWD detection) and SWD detector
(Gonzalez-Sulser-Team/SWD-Automatic-Identification, commit 3fd91f5; MATLAB cepstral detection).
The scorer shows how it puts SWD epochs into the sleep scores (the per-epoch SWD file value 6
becomes state 4) and that the stored spectra are its own (log10 multitaper amplitude referenced to
each frequency's 5th percentile), not the Welch spectra of the Methods. Here SWD epochs are the
epochs SWD events overlap, as the Methods have it, because the scorer's flags in the dandiset each
overlap an event, mostly as its last epoch, and cover about half of the epochs events overlap (data
check); the flags are a variant. Neither repository holds
the statistics or figure code.
| Claim | Source | Status |
|---|---|---|
| fig2b-swd-hourly: more SWDs per hour in Het, no hour effect | Fig. 2B | tested |
| fig2c-swd-number: more SWDs over 24 h | Fig. 2C | tested |
| fig2c-swd-mean-duration: longer SWDs | Fig. 2C | tested |
| fig2c-swd-total-duration: more SWD time | Fig. 2C | tested |
| disc-het-swd-minutes: Het 59.97 min, WT under 1 min | Discussion | tested |
| figs5a-swd-light-dark: no light-dark difference | Fig. S5A | tested |
| figs8-swd-sex: no sex effect on SWDs | Fig. S8B | tested |
| fig2d-swd-spectrum-peak: theta peak with harmonic | Fig. 2D | tested |
| fig2e-swd-delta-power: higher SWD delta power in Het | Fig. 2E | blocked (data missing: one WT with SWDs) |
| fig3d-rem-hourly: less REM per hour | Fig. 3D | tested |
| fig3g-rem-total, fig3g-rem-bouts, fig3g-rem-bout-duration | Fig. 3G | tested |
| fig3h-nrem-total, fig3h-nrem-bouts | Fig. 3H | tested |
| fig3i-wake-total, fig3i-wake-bouts | Fig. 3I | tested |
| fig4a-rem-spectrum, fig4b-nrem-beta, fig4c-wake-power | Fig. 4 | tested |
| fig5a-swd-state-distribution: 94.8% wake, 4.4% NREM, 0.5% REM | Fig. 5A | tested |
| fig5a-het-more-nrem-swds: larger NREM share in Het | Fig. 5A | blocked (data missing: one WT with SWDs) |
| fig5c-transition-swds: all Het have transitional SWDs | Fig. 5C | tested |
| figs5b-transition-percent: 24.2% transitional | Fig. S5B | tested |
| fig5f-swd-wake-correlation: R = -0.51 | Fig. 5F | tested |
| fig5def-no-other-correlations | Fig. 5D-F | tested |
| fig1c, fig1df: western blots | Fig. 1, S2 | blocked (data missing) |
| figs4-normal-development | Fig. S4 | blocked (data missing) |
| figs1-scorer-agreement | Table S1 | blocked (data missing) |
| fig6b, fig6c, fig6d, fig7c: pharmacology | Fig. 6, 7 | blocked (data missing) |
In [1]:
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from concurrent.futures import ThreadPoolExecutor
from scipy import signal, stats
from statsmodels.formula.api import ols
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import grin2b
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
version = grin2b.check_version()
sessions = grin2b.load_all()
by_id = {s["subject_id"]: s for s in sessions}
print(f"dandiset {grin2b.DANDISET_ID} version {grin2b.DANDISET_VERSION} (modified {version['modified']}, "
f"{version['asset_count']} assets): loaded {len(sessions)} files in {time.time() - started:.0f} s")
dandiset 001888 version draft (modified 2026-07-31T09:37:44.644146Z, 17 assets): loaded 17 files in 39 s
Data check¶
Each check of the description against the values is stated under its output.
In [2]:
rows = []
for s in sessions:
rows.append({
"subject": s["subject_id"], "sex": s["sex"], "genotype_stored": s["genotype_stored"],
"genotype": s["genotype"], "age": s["age"], "strain": s["strain"],
"session_start": s["session_start_time"][:19], "hours": round(s["eeg_shape"][0] / s["rate"] / 3600, 1),
"eeg_ch": s["eeg_shape"][1], "emg_ch": s["emg_shape"][1], "rate_hz": s["rate"],
"BL1_start_h": round(s["days"][0]["t0"] / 3600, 2), "BL2_start_h": round(s["days"][1]["t0"] / 3600, 2),
"BL_lengths_h": sorted({round((d["t1"] - d["t0"]) / 3600, 4) for d in s["days"]}),
"epoch_tags": list(s["epochs"].tags)})
subjects = pd.DataFrame(rows)
print(subjects.to_string(index=False))
print("\nSubject fields in the files:", sorted({tuple(s["subject_fields"]) for s in sessions}))
print("Files whose subject holds a genotype field:", sum("genotype" in s["subject_fields"] for s in sessions))
print("\nGenotype by sex:")
print(pd.crosstab(subjects.genotype.fillna("none"), subjects.sex).to_string())
print("Paper (Supplementary Methods): Het 8 male + 7 female = 15, WT 3 male + 6 female = 9 analysed;"
" 4 Het and 9 WT more were recorded and excluded.")
print("\nRaw data: unit", sorted({s["unit"] for s in sessions}), "| conversion (V per code)",
sorted({s["conversion"] for s in sessions}), "| offset", sorted({s["offset"] for s in sessions}),
"| starting_time", sorted({s["starting_time"] for s in sessions}))
print("13 mV / 4096 codes =", 13e-3 / 4096, "V per code")
print("EEG channels:", sessions[0]["eeg_channels"], "| EMG channels:", sessions[0]["emg_channels"])
print("Same channel order in every file:", len({tuple(s["eeg_channels"]) for s in sessions}) == 1)
print("Storage chunk of the EEG (samples x channels):", sorted({s["eeg_chunks"] for s in sessions}))
subject sex genotype_stored genotype age strain session_start hours eeg_ch emg_ch rate_hz BL1_start_h BL2_start_h BL_lengths_h epoch_tags
GRIN2B-129 M Het Het P16W/P19W Long-Evans 2021-03-26T10:55:59 77.0 14 2 250.4 20.07 44.07 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-130 M Het Het P16W/P19W Long-Evans 2021-03-26T10:55:59 77.0 14 2 250.4 20.07 44.07 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-131 M Het Het P16W/P19W Long-Evans 2021-03-26T10:55:59 77.0 14 2 250.4 20.07 44.07 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-132 U NaN NaN P16W/P19W Long-Evans 2021-03-26T10:55:59 77.1 14 2 250.4 20.07 44.07 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-137 F Het Het P16W/P19W Long-Evans 2021-03-31T10:21:59 73.5 14 2 250.4 20.63 44.63 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-138 F Het Het P16W/P19W Long-Evans 2021-03-31T10:21:59 73.5 14 2 250.4 20.63 44.63 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-139 F Het Het P16W/P19W Long-Evans 2021-03-31T10:21:59 73.5 14 2 250.4 20.63 44.63 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-227 M Het Het P16W/P19W Long-Evans 2021-07-28T10:07:59 99.5 14 2 250.4 20.87 44.87 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-228 M Het Het P16W/P19W Long-Evans 2021-07-28T10:07:59 99.5 14 2 250.4 20.87 44.87 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-229 M Het Het P16W/P19W Long-Evans 2021-07-28T10:07:59 99.5 14 2 250.4 20.87 44.87 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-236 F Het Het P16W/P19W Long-Evans 2021-08-29T08:54:59 77.1 14 2 250.4 22.08 46.08 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-237 F Het Het P16W/P19W Long-Evans 2021-08-16T10:34:59 76.0 14 2 250.4 20.42 44.42 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-238 F Wt WT P16W/P19W Long-Evans 2021-08-16T10:34:59 76.1 14 2 250.4 20.42 44.42 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-239 F Het Het P16W/P19W Long-Evans 2021-08-16T10:34:59 76.0 14 2 250.4 20.42 44.42 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-240 F Wt WT P16W/P19W Long-Evans 2021-08-16T10:34:59 76.0 14 2 250.4 20.42 44.42 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-241 F Het Het P16W/P19W Long-Evans 2021-08-20T09:52:59 75.0 14 2 250.4 21.12 45.12 [24.0] [baseline_window_1, baseline_window_2]
GRIN2B-373 M Het Het P16W/P19W Long-Evans 2022-04-28T14:08:59 92.9 14 2 250.4 16.85 40.85 [24.0] [baseline_window_1, baseline_window_2]
Subject fields in the files: [('age', 'description', 'genotype', 'sex', 'species', 'strain', 'subject_id'), ('age', 'description', 'sex', 'species', 'strain', 'subject_id')]
Files whose subject holds a genotype field: 16
Genotype by sex:
sex F M U
genotype
Het 7 7 0
WT 2 0 0
none 0 0 1
Paper (Supplementary Methods): Het 8 male + 7 female = 15, WT 3 male + 6 female = 9 analysed; 4 Het and 9 WT more were recorded and excluded.
Raw data: unit ['volts'] | conversion (V per code) [3.173828125e-06] | offset [0.0] | starting_time [0.0]
13 mV / 4096 codes = 3.173828125e-06 V per code
EEG channels: ['S1_Tr_R', 'M2_Fra_R', 'M2_anterior_R', 'M1_anterior_R', 'V2_ML_R', 'V1_M_R', 'S1Hl_S1Fl_R', 'V1_M_L', 'V2_ML_L', 'S1Hl_S1Fl_L', 'M1_anterior_L', 'M2_anterior_L', 'M2_Fra_L', 'S1_Tr_L'] | EMG channels: ['EMG_R', 'EMG_L']
Same channel order in every file: True
Storage chunk of the EEG (samples x channels): [(357142, 14)]
Sixteen files give the genotype (Het or Wt) in the subject; GRIN2B-132 has no genotype and
sex U. With 14 Het (7 male, 7 female) and 2 WT (both female), the dandiset fits within the
paper's analysed groups (Het 8 male and 7 female, WT 3 male and 6 female) but holds 2 of its 9 WT;
it does not say which rats the paper analysed. All files are Long-Evans, aged P16W/P19W (a range, the same for all), recorded at
250.4 Hz with 14 EEG and 2 EMG channels, and each baseline window lasts exactly 24 h. The stated
conversion is 13 mV over 4096 codes, as the description says, with an offset of 0 (see the raw
values below).
In [3]:
rows = []
for s in sessions:
for d in s["days"]:
st, fl = d["states"], d["flags"]
ev = grin2b.event_epochs(d)
onset_ep = (d["sz_start"] // grin2b.EPOCH_S).astype(int)
z = s["seizure_totals_by_zt"]; z = z[z.baseline_window == d["bl"]].sort_values("zeitgeber_hour")
hourly = np.bincount(np.floor(d["sz_start"] / 3600).astype(int), minlength=24)[:24]
last = np.zeros(len(st), bool)
last[np.minimum(np.ceil(d["sz_stop"] / grin2b.EPOCH_S).astype(int) - 1, len(st) - 1)] = True
rows.append({
"subject": s["subject_id"], "day": d["bl"][-1], "epochs": len(st),
"grid_ok": bool(np.allclose(np.diff(d["epoch_start"]), 5) and abs(d["epoch_start"][0]) < 0.01
and np.allclose(d["flag_start"], d["epoch_start"])),
"state_codes": "".join(map(str, sorted(set(st)))), "flag_values": "/".join(map(str, sorted(set(fl)))),
"flag6_epochs": int((fl == 6).sum()), "swd_events": len(d["sz_start"]),
"events_in_day": bool(np.all((d["sz_start"] >= 0) & (d["sz_stop"] <= 86400))),
"dur_eq_stop_minus_start": bool(np.allclose(d["sz_stop"] - d["sz_start"], d["sz_dur"])) if len(d["sz_dur"]) else None,
"dur_range_s": f"{d['sz_dur'].min():.1f}-{d['sz_dur'].max():.1f}" if len(d["sz_dur"]) else "",
"onsets_in_flag6": round(float(np.mean(fl[onset_ep] == 6)), 3) if len(onset_ep) else np.nan,
"flag6_with_event": round(float(np.mean(ev[fl == 6])), 3) if (fl == 6).any() else np.nan,
"flag6_is_event_last_epoch": round(float(np.mean(last[fl == 6])), 3) if (fl == 6).any() else np.nan,
"event_epochs_flagged": round(float(np.mean(fl[ev] == 6)), 3) if ev.any() else np.nan,
"zt_swd_total_eq_events": int(z.spike_wave_discharge.sum()) == len(d["sz_start"]),
"zt_swd_max_hour_diff": int(np.abs(z.spike_wave_discharge.to_numpy() - hourly).max()),
"zt_n_events_row0_total": bool(z.number_of_seizure_events.iloc[0] == len(d["sz_start"])
and (z.number_of_seizure_events.iloc[1:] == 0).all()),
"zt_Day_row0_total": bool(z.Day.iloc[0] == len(d["sz_start"]) and (z.Day.iloc[1:] == 0).all()),
})
tables = pd.DataFrame(rows)
print(tables.to_string(index=False))
print("\nTable descriptions:")
for k, v in sessions[0]["descriptions"].items():
print(f" {k}: {v}")
subject day epochs grid_ok state_codes flag_values flag6_epochs swd_events events_in_day dur_eq_stop_minus_start dur_range_s onsets_in_flag6 flag6_with_event flag6_is_event_last_epoch event_epochs_flagged zt_swd_total_eq_events zt_swd_max_hour_diff zt_n_events_row0_total zt_Day_row0_total GRIN2B-129 1 17280 True 012 4/6 686 753 True True 0.8-21.2 0.127 1.0 0.743 0.513 True 1 True True GRIN2B-129 2 17280 True 012 4/6 855 816 True True 0.8-40.8 0.137 1.0 0.674 0.550 True 1 True True GRIN2B-130 1 17280 True 012 4/6 459 483 True True 0.8-15.8 0.093 1.0 0.789 0.512 True 1 True True GRIN2B-130 2 17280 True 012 4/6 497 531 True True 0.8-21.4 0.087 1.0 0.777 0.507 True 0 True True GRIN2B-131 1 17280 True 012 4/6 724 769 True True 0.8-50.6 0.105 1.0 0.742 0.514 True 0 True True GRIN2B-131 2 17280 True 012 4/6 848 863 True True 0.8-28.6 0.102 1.0 0.742 0.524 True 0 True True GRIN2B-132 1 17280 True 012 4/6 35 106 True True 0.8-5.0 0.038 1.0 1.000 0.255 True 0 True True GRIN2B-132 2 17280 True 012 4/6 39 108 True True 0.8-4.6 0.046 1.0 1.000 0.275 True 0 True True GRIN2B-137 1 17280 True 012 4/6 283 476 True True 0.8-11.4 0.095 1.0 0.947 0.397 True 1 True True GRIN2B-137 2 17280 True 012 4/6 464 697 True True 0.8-15.6 0.082 1.0 0.912 0.420 True 1 True True GRIN2B-138 1 17280 True 012 4/6 893 1019 True True 0.8-25.0 0.103 1.0 0.800 0.494 True 0 True True GRIN2B-138 2 17280 True 012 4/6 683 768 True True 0.8-30.0 0.098 1.0 0.773 0.497 True 0 True True GRIN2B-139 1 17280 True 012 4/6 1124 749 True True 0.8-30.6 0.091 1.0 0.566 0.623 True 1 True True GRIN2B-139 2 17280 True 012 4/6 1684 1144 True True 0.8-32.6 0.113 1.0 0.580 0.625 True 1 True True GRIN2B-227 1 17280 True 012 4/6 965 938 True True 0.8-40.6 0.117 1.0 0.731 0.540 True 1 True True GRIN2B-227 2 17280 True 012 4/6 1035 977 True True 0.8-49.0 0.141 1.0 0.704 0.554 True 0 True True GRIN2B-228 1 17280 True 012 4 0 0 True None NaN NaN NaN NaN True 0 True True GRIN2B-228 2 17280 True 012 4 0 2 True True 0.8-0.8 0.000 NaN NaN 0.000 True 0 True True GRIN2B-229 1 17280 True 012 4/6 336 544 True True 0.8-17.8 0.094 1.0 0.893 0.409 True 1 True True GRIN2B-229 2 17280 True 012 4/6 180 343 True True 0.8-8.8 0.067 1.0 0.939 0.362 True 1 True True GRIN2B-236 1 17280 True 012 4/6 936 896 True True 0.8-31.0 0.156 1.0 0.682 0.557 True 1 True True GRIN2B-236 2 17280 True 012 4/6 826 620 True True 0.8-55.4 0.194 1.0 0.602 0.625 True 0 True True GRIN2B-237 1 17280 True 012 4/6 503 636 True True 0.8-17.8 0.129 1.0 0.821 0.480 True 0 True True GRIN2B-237 2 17280 True 012 4/6 415 516 True True 0.8-20.8 0.128 1.0 0.798 0.480 True 0 True True GRIN2B-238 1 17280 True 012 4 0 0 True None NaN NaN NaN NaN True 0 True True GRIN2B-238 2 17280 True 012 4 0 0 True None NaN NaN NaN NaN True 0 True True GRIN2B-239 1 17280 True 012 4/6 671 765 True True 0.8-33.2 0.095 1.0 0.803 0.494 True 0 True True GRIN2B-239 2 17280 True 012 4/6 960 1003 True True 0.8-22.0 0.120 1.0 0.756 0.522 True 0 True True GRIN2B-240 1 17280 True 012 4/6 42 106 True True 0.8-4.6 0.038 1.0 1.000 0.292 True 0 True True GRIN2B-240 2 17280 True 012 4/6 35 119 True True 0.8-3.0 0.034 1.0 1.000 0.233 True 0 True True GRIN2B-241 1 17280 True 012 4/6 447 630 True True 0.8-19.4 0.100 1.0 0.848 0.443 True 1 True True GRIN2B-241 2 17280 True 012 4/6 701 767 True True 0.8-39.8 0.126 1.0 0.739 0.514 True 1 True True GRIN2B-373 1 17280 True 012 4/6 1292 867 True True 0.8-64.6 0.178 1.0 0.523 0.648 True 1 True True GRIN2B-373 2 17280 True 012 4/6 1759 765 True True 0.8-60.0 0.152 1.0 0.379 0.731 True 1 True True Table descriptions: seizure_events: Seizure events detected by the Gonzalez-Sulser lab automated scoring pipeline. start_time/stop_time are in seconds relative to session_start_time. seizure_totals_by_zt: Per-Zeitgeber-hour seizure event totals. ZT (Zeitgeber Time) = hours since lights-on. TODO: confirm column definitions (N_event, mean_dur, SWDs, Day) with lab. sleep_states: Automated sleep state classifications in 5-second epochs produced by the Gonzalez-Sulser lab scoring pipeline (https://github.com/Gonzalez-Sulser-Team/AUTOMATIC-SLEEP-SCORER). State codes: 0=Wake, 1=NREM, 2=REM, 4=SWD (spike-wave discharge epochs, classed separately from the three sleep states). swd_epoch_counts: Number of spike-wave discharge (SWD) events per 5-second epoch, computed by the Gonzalez-Sulser lab automated scoring pipeline. start_time/stop_time are in seconds relative to session_start_time. sleep_state_power_spectra: Average power spectral density per brain state, computed by the Gonzalez-Sulser lab automated scoring pipeline (https://github.com/Gonzalez-Sulser-Team/AUTOMATIC-SLEEP-SCORER). Frequency bins from 0 to 125.2 Hz (Nyquist of 250.4 Hz) in 0.2 Hz steps. State columns: Wake, NREM, REM, SWD (spike-wave discharge epochs). Spectra are baseline-corrected by normalizing to the average spectral power across REM, NREM, and wake for each animal. TODO: confirm absolute units (µV^2/Hz or normalised).
The scores are on a 5-s grid of 17,280 epochs per day starting at the baseline window, and every
SWD event falls inside its day with duration equal to stop minus start. Three descriptions do
not match the values:
sleep_statesis said to use code 4 for SWD epochs; no file holds code 4. The stored scores are the scorer'sdge_okoutput, before SWD epochs are set.swd_epoch_countsis said to hold the number of SWD events per epoch; it holds only 4 and 6. In the lab's scorer these are flags (6 = SWD epoch, 4 = none). Every flagged epoch overlaps an SWD event, and most flagged epochs are the last epoch an event overlaps, while only about a tenth of SWD onsets fall in a flagged epoch (columnsflag6_with_event,flag6_is_event_last_epoch,onsets_in_flag6). The flags cover about half of the epochs SWD events overlap (event_epochs_flagged), so removing flagged epochs from the sleep states, as the scorer does, leaves the other SWD epochs in the states. The Methods exclude "SWD times"; SWD epochs are therefore the epochs an SWD event overlaps here, and the flags are a variant.- In
seizure_totals_by_ztthe per-hour SWD counts are inspike_wave_discharge(their total equals the number of events, and per hour they equal the events binned by start time to within one event at an hour boundary), whilenumber_of_seizure_eventsandDayhold the day's total in the first row and 0 elsewhere. The description itself marks these columns as unconfirmed (TODO). Per-hour counts are taken from the events here.
GRIN2B-228 (Het) has 0 and 2 SWDs, and GRIN2B-238 (WT) has none; their files hold no
seizure_events table, or an empty one.
In [4]:
rows = []
for s in sessions:
for d in s["days"]:
h = np.arange(len(d["states"])) * grin2b.EPOCH_S / 3600
sleep = np.isin(d["states"], [1, 2])
rows.append({"subject": s["subject_id"], "day": d["bl"][-1],
"sleep_frac_ZT0-12": sleep[h < 12].mean(), "sleep_frac_ZT12-24": sleep[h >= 12].mean()})
lightdark = pd.DataFrame(rows)
lightdark["light_minus_dark"] = lightdark["sleep_frac_ZT0-12"] - lightdark["sleep_frac_ZT12-24"]
print(lightdark.round(3).to_string(index=False))
print("\nDays with more sleep in ZT 0-12 than ZT 12-24:", int((lightdark.light_minus_dark > 0).sum()), "of", len(lightdark))
subject day sleep_frac_ZT0-12 sleep_frac_ZT12-24 light_minus_dark GRIN2B-129 1 0.531 0.249 0.282 GRIN2B-129 2 0.491 0.103 0.389 GRIN2B-130 1 0.503 0.200 0.303 GRIN2B-130 2 0.516 0.173 0.343 GRIN2B-131 1 0.540 0.212 0.328 GRIN2B-131 2 0.579 0.234 0.345 GRIN2B-132 1 0.353 0.187 0.166 GRIN2B-132 2 0.420 0.164 0.256 GRIN2B-137 1 0.490 0.144 0.346 GRIN2B-137 2 0.496 0.159 0.337 GRIN2B-138 1 0.410 0.266 0.144 GRIN2B-138 2 0.470 0.244 0.226 GRIN2B-139 1 0.538 0.158 0.380 GRIN2B-139 2 0.443 0.144 0.299 GRIN2B-227 1 0.520 0.208 0.312 GRIN2B-227 2 0.535 0.208 0.327 GRIN2B-228 1 0.464 0.209 0.255 GRIN2B-228 2 0.419 0.216 0.203 GRIN2B-229 1 0.488 0.182 0.306 GRIN2B-229 2 0.513 0.215 0.299 GRIN2B-236 1 0.424 0.251 0.173 GRIN2B-236 2 0.443 0.223 0.220 GRIN2B-237 1 0.421 0.303 0.118 GRIN2B-237 2 0.500 0.317 0.183 GRIN2B-238 1 0.416 0.178 0.237 GRIN2B-238 2 0.327 0.259 0.068 GRIN2B-239 1 0.473 0.194 0.279 GRIN2B-239 2 0.500 0.159 0.341 GRIN2B-240 1 0.532 0.273 0.260 GRIN2B-240 2 0.473 0.259 0.213 GRIN2B-241 1 0.553 0.182 0.372 GRIN2B-241 2 0.448 0.155 0.293 GRIN2B-373 1 0.476 0.277 0.198 GRIN2B-373 2 0.425 0.245 0.181 Days with more sleep in ZT 0-12 than ZT 12-24: 34 of 34
Rats sleep mostly in the light phase. Every day has more sleep in its first 12 h than in its
second, which is consistent with each baseline window starting at lights on (ZT 0), as the
description says. The files derive session_start_time from that assumption, so the absolute
clock time is not an independent check.
In [5]:
rows = []
for s in sessions:
sp = s["spectra"]
for bl, g in sp.groupby("baseline_window"):
f = g.frequency_hz.to_numpy()
row = {"subject": s["subject_id"], "day": bl[-1], "rows": len(g), "f_min": f.min(), "f_max": f.max(),
"f_step": np.unique(np.round(np.diff(f), 3)).tolist()}
for c in ("Wake", "NREM", "REM", "SWD"):
row[c] = f"{g[c].min():.3f}..{g[c].max():.3f}" if c in g else "absent"
rows.append(row)
print(pd.DataFrame(rows).to_string(index=False))
print("\nDescription:", sessions[0]["descriptions"]["sleep_state_power_spectra"])
subject day rows f_min f_max f_step Wake NREM REM SWD GRIN2B-129 1 627 0.0 125.2 [0.2] 0.042..0.520 0.004..0.718 0.005..0.286 0.069..0.817 GRIN2B-129 2 627 0.0 125.2 [0.2] 0.046..1.237 0.007..1.527 0.006..0.962 0.093..1.486 GRIN2B-130 1 627 0.0 125.2 [0.2] 0.032..0.299 0.003..0.773 0.009..0.318 0.061..0.800 GRIN2B-130 2 627 0.0 125.2 [0.2] 0.035..0.301 0.004..0.795 0.008..0.341 0.064..0.817 GRIN2B-131 1 627 0.0 125.2 [0.2] 0.042..0.366 0.004..0.764 0.007..0.394 0.056..0.926 GRIN2B-131 2 627 0.0 125.2 [0.2] 0.047..1.114 0.007..1.612 0.010..1.076 0.080..1.518 GRIN2B-132 1 627 0.0 125.2 [0.2] 0.034..1.124 0.006..1.612 0.009..1.205 0.055..1.450 GRIN2B-132 2 627 0.0 125.2 [0.2] 0.035..1.112 0.006..1.633 0.009..1.198 0.054..1.424 GRIN2B-137 1 627 0.0 125.2 [0.2] 0.057..0.337 0.006..0.806 0.009..0.371 0.108..0.805 GRIN2B-137 2 627 0.0 125.2 [0.2] 0.053..0.325 0.006..0.751 0.008..0.338 0.080..0.756 GRIN2B-138 1 627 0.0 125.2 [0.2] 0.048..0.355 0.006..0.816 0.009..0.266 0.083..0.764 GRIN2B-138 2 627 0.0 125.2 [0.2] 0.043..0.354 0.007..0.849 0.010..0.310 0.069..0.814 GRIN2B-139 1 627 0.0 125.2 [0.2] 0.036..0.323 0.005..0.741 0.010..0.329 0.035..0.885 GRIN2B-139 2 627 0.0 125.2 [0.2] 0.042..0.335 0.005..0.729 0.010..0.305 0.039..0.897 GRIN2B-227 1 627 0.0 125.2 [0.2] 0.036..0.332 0.004..0.726 0.005..0.240 0.036..0.916 GRIN2B-227 2 627 0.0 125.2 [0.2] 0.037..0.394 0.008..0.849 0.011..0.365 0.036..1.093 GRIN2B-228 1 627 0.0 125.2 [0.2] 0.032..0.369 0.004..0.636 0.005..0.253 absent GRIN2B-228 2 627 0.0 125.2 [0.2] 0.030..0.380 0.003..0.699 0.003..0.234 absent GRIN2B-229 1 627 0.0 125.2 [0.2] 0.046..0.405 0.004..0.798 0.005..0.461 0.051..0.941 GRIN2B-229 2 627 0.0 125.2 [0.2] 0.061..0.359 0.003..0.728 0.004..0.372 0.045..0.864 GRIN2B-236 1 627 0.0 125.2 [0.2] 0.043..0.550 0.009..0.886 0.008..0.650 0.034..1.413 GRIN2B-236 2 627 0.0 125.2 [0.2] 0.044..1.136 0.009..1.600 0.011..1.075 0.038..1.654 GRIN2B-237 1 627 0.0 125.2 [0.2] 0.034..0.339 0.006..0.680 0.009..0.366 0.027..0.952 GRIN2B-237 2 627 0.0 125.2 [0.2] 0.037..0.323 0.004..0.542 0.007..0.321 0.023..0.884 GRIN2B-238 1 627 0.0 125.2 [0.2] 0.043..0.698 0.023..0.997 0.026..0.787 absent GRIN2B-238 2 627 0.0 125.2 [0.2] 0.041..0.625 0.015..0.972 0.017..0.722 absent GRIN2B-239 1 627 0.0 125.2 [0.2] 0.058..0.397 0.006..0.742 0.007..0.325 0.037..0.921 GRIN2B-239 2 627 0.0 125.2 [0.2] 0.058..0.434 0.006..0.782 0.007..0.370 0.039..0.948 GRIN2B-240 1 627 0.0 125.2 [0.2] 0.040..0.400 0.004..0.623 0.006..0.319 0.025..0.615 GRIN2B-240 2 627 0.0 125.2 [0.2] 0.043..0.408 0.003..0.653 0.006..0.340 0.027..0.536 GRIN2B-241 1 627 0.0 125.2 [0.2] 0.040..0.367 0.008..0.782 0.009..0.398 0.051..0.940 GRIN2B-241 2 627 0.0 125.2 [0.2] 0.037..1.148 0.013..1.593 0.015..1.112 0.054..1.489 GRIN2B-373 1 627 0.0 125.2 [0.2] 0.036..0.395 0.006..0.766 0.007..0.339 0.041..1.110 GRIN2B-373 2 627 0.0 125.2 [0.2] 0.038..0.509 0.006..0.791 0.008..0.389 0.043..1.118 Description: Average power spectral density per brain state, computed by the Gonzalez-Sulser lab automated scoring pipeline (https://github.com/Gonzalez-Sulser-Team/AUTOMATIC-SLEEP-SCORER). Frequency bins from 0 to 125.2 Hz (Nyquist of 250.4 Hz) in 0.2 Hz steps. State columns: Wake, NREM, REM, SWD (spike-wave discharge epochs). Spectra are baseline-corrected by normalizing to the average spectral power across REM, NREM, and wake for each animal. TODO: confirm absolute units (µV^2/Hz or normalised).
The stored spectra run from 0 to 125.2 Hz in 0.2 Hz steps for each day and state, with no SWD column for the two rats without SWDs in a day. Their values (0 to about 1.6) are not powers: the lab's scorer writes the mean over a state's epochs of the log10 multitaper amplitude, referenced per frequency to its 5th percentile over the day. The description (baseline-corrected PSD normalized to the average over states, units to be confirmed) describes the Methods' Welch pipeline, not these values. The spectral claims below use spectra computed from the raw EEG as the Methods describe; the stored spectra are a variant.
In [6]:
t_raw = time.time()
def raw_epochs(s):
# Raw EEG (SWD channel) and EMG in one storage chunk at ZT 0, 4, 8, 12, 16, 20 of each day,
# reduced at once to per-epoch spectra and EMG power, and 0.5-s SWD-band power.
days = {d["bl"]: d for d in s["days"]}
eps, eps_alt, swd_bins, codes = [], [], [], []
for w in grin2b.raw_windows(s):
d = days[w["bl"]]
e = grin2b.window_epochs(w, s, d)
e["bl"] = w["bl"]
e["zt"] = w["zt"]
eps.append(e)
ea = grin2b.window_epochs(w, s, d, "eeg_alt")
ea["bl"] = w["bl"]
eps_alt.append(ea)
t, p, ok = grin2b.swd_band_series(w, s)
swd_bins.append({"bl": w["bl"], "t": t, "p": p, "ok": ok})
x = w["eeg"]
valid = x[(x >= 0) & (x <= grin2b.ADC_MAX)]
codes.append((valid.min() if valid.size else np.nan, np.median(valid) if valid.size else np.nan,
valid.max() if valid.size else np.nan, (x == grin2b.DROPOUT).mean(),
((x < 0) | ((x > grin2b.ADC_MAX) & (x != grin2b.DROPOUT))).mean()))
e = pd.concat(eps, ignore_index=True)
e["subject"] = s["subject_id"]
return e, swd_bins, np.array(codes, float), pd.concat(eps_alt, ignore_index=True)
with ThreadPoolExecutor(6) as ex:
raw = dict(zip(by_id, ex.map(raw_epochs, sessions)))
FREQ = grin2b.welch_freq()
rows = []
for sid, (e, _, c, _) in raw.items():
rows.append({"subject": sid, "epochs": len(e), "epochs_with_bad_samples": round(e.dropout.mean(), 3),
"artefact_over_3000uV": round((e.artefact & ~e.dropout).mean(), 3),
"usable": int((~e.dropout & ~e.artefact).sum()),
"valid_code_min": np.nanmin(c[:, 0]), "valid_code_median": np.nanmedian(c[:, 1]),
"valid_code_max": np.nanmax(c[:, 2]), "samples_32767": round(c[:, 3].mean(), 4),
"samples_other_invalid": round(c[:, 4].mean(), 6)})
raw_summary = pd.DataFrame(rows)
print(f"read 12 raw windows per animal in {time.time() - t_raw:.0f} s")
print(raw_summary.to_string(index=False))
read 12 raw windows per animal in 267 s subject epochs epochs_with_bad_samples artefact_over_3000uV usable valid_code_min valid_code_median valid_code_max samples_32767 samples_other_invalid GRIN2B-129 3411 0.322 0.000 2313 0.0 2658.5 3122.0 0.0887 0.000001 GRIN2B-130 3411 0.206 0.009 2679 0.0 2705.0 4095.0 0.0325 0.000114 GRIN2B-131 3411 0.239 0.004 2583 0.0 2628.0 4095.0 0.0492 0.000103 GRIN2B-132 3411 0.352 0.000 2212 0.0 2604.0 4000.0 0.1312 0.000000 GRIN2B-137 3412 0.364 0.000 2170 0.0 2641.5 3442.0 0.0375 0.000000 GRIN2B-138 3412 0.482 0.000 1769 0.0 2585.0 4095.0 0.0782 0.000000 GRIN2B-139 3412 0.234 0.029 2517 0.0 2700.0 4095.0 0.0391 0.000178 GRIN2B-227 3412 0.264 0.000 2511 2335.0 2669.5 3797.0 0.0365 0.000000 GRIN2B-228 3412 0.309 0.000 2359 2260.0 2617.5 2922.0 0.0665 0.000000 GRIN2B-229 3412 0.240 0.000 2592 2126.0 2604.5 3098.0 0.0562 0.000000 GRIN2B-236 3410 0.416 0.000 1990 2291.0 2594.0 3043.0 0.1058 0.000000 GRIN2B-237 3411 0.221 0.000 2658 0.0 2618.5 3313.0 0.0297 0.000001 GRIN2B-238 3411 0.725 0.000 937 0.0 2614.0 3005.0 0.2148 0.000002 GRIN2B-239 3411 0.317 0.000 2329 0.0 2669.5 3005.0 0.0764 0.000001 GRIN2B-240 3411 0.218 0.001 2663 0.0 2648.5 3001.0 0.0643 0.000000 GRIN2B-241 3410 0.389 0.000 2083 2255.0 2610.0 3505.0 0.1494 0.000000 GRIN2B-373 3412 0.441 0.000 1908 0.0 2587.5 3600.0 0.1042 0.000000
The raw EEG is stored as the amplifier's unsigned 12-bit codes: valid samples sit around a median
near 2600, not around 0, so with the stated offset of 0 every channel reads about +8 mV of DC.
Signal loss is written as 32767 on every channel at once, a value the description does not
mention, and some files also hold other values outside 0 to 4095 (negative, or between 4096 and
32766), which a 12-bit amplifier cannot produce. Epochs touching any of these are dropped here;
the fraction varies widely between rats (column epochs_with_bad_samples; 0.73 for GRIN2B-238 in
the windows read). The DC level is removed before any spectrum.
In [7]:
# Clock check: do the sleep scores and SWD events line up with the raw signals they were scored from?
emg_shifts = np.arange(-60, 61)
swd_shifts = np.arange(-10, 10.5, 0.5)
rows, emg_curves, swd_curves = [], {}, {}
for s in sessions:
sid = s["subject_id"]
e, bins = raw[sid][:2]
days = {d["bl"]: d for d in s["days"]}
ok = ~e.dropout
curve = []
for k in emg_shifts:
lab = np.concatenate([days[bl]["states"][np.clip(g.epoch.to_numpy() - k, 0, grin2b.EPOCHS_PER_DAY - 1)]
for bl, g in e[ok].groupby("bl")])
v = np.concatenate([g.emg_logpow.to_numpy() for bl, g in e[ok].groupby("bl")])
curve.append(v[lab == 0].mean() - v[lab == 1].mean())
emg_curves[sid] = np.array(curve)
row = {"subject": sid, "emg_wake_minus_nrem_at_0": round(curve[60], 3),
"emg_best_shift_epochs": int(emg_shifts[np.argmax(curve)]), "emg_at_+-60": round((curve[0] + curve[-1]) / 2, 3)}
n_in = sum(grin2b.in_events(b["t"], days[b["bl"]]["sz_start"], days[b["bl"]]["sz_stop"]).sum() for b in bins)
if n_in >= 20:
c2 = []
for sh in swd_shifts:
inside, outside = [], []
for b in bins:
m = grin2b.in_events(b["t"], days[b["bl"]]["sz_start"], days[b["bl"]]["sz_stop"], sh)
inside.append(b["p"][m & b["ok"]]); outside.append(b["p"][~m & b["ok"]])
c2.append(np.concatenate(inside).mean() - np.concatenate(outside).mean())
swd_curves[sid] = np.array(c2)
row.update({"swd_bins_inside": int(n_in), "swd_contrast_at_0": round(c2[20], 3),
"swd_best_shift_s": float(swd_shifts[np.argmax(c2)]), "swd_contrast_at_+-10s": round((c2[0] + c2[-1]) / 2, 3)})
else:
row.update({"swd_bins_inside": int(n_in)})
# Per-epoch 5-10 Hz power, against the epochs SWD events overlap and the epochs the scorer flags,
# each shifted by k epochs.
ep_p = []
for b in bins:
k = np.floor(b["t"] / grin2b.EPOCH_S).astype(int)
df_ = pd.DataFrame({"epoch": k[b["ok"]], "p": b["p"][b["ok"]]}).groupby("epoch").p.mean()
ep_p.append((b["bl"], df_))
if n_in >= 20:
for name, lab_fn in (("events", grin2b.event_epochs), ("flags", lambda d: d["flags"] == grin2b.SWD_FLAG)):
best = []
for k in range(-3, 4):
ins, outs = [], []
for bl, pe in ep_p:
lab = lab_fn(days[bl])
idx = np.clip(pe.index.to_numpy() - k, 0, grin2b.EPOCHS_PER_DAY - 1)
ins.append(pe.to_numpy()[lab[idx]]); outs.append(pe.to_numpy()[~lab[idx]])
ins = np.concatenate(ins)
best.append(ins.mean() - np.concatenate(outs).mean() if ins.size else np.nan)
row[f"{name}_epoch_best_shift"] = int(np.arange(-3, 4)[np.nanargmax(best)])
row[f"{name}_epoch_contrast_at_0"] = round(best[3], 3)
rows.append(row)
clock = pd.DataFrame(rows)
print("EMG: log10 60-90 Hz power (right minus left EMG), mean in wake minus mean in NREM epochs, with the scores shifted by k epochs.")
print("SWD: log10 5-10 Hz EEG power in 0.5-s bins, mean inside minus outside SWD events, with the events shifted by s seconds.")
print("events/flags_epoch: the same power per 5-s epoch, in epochs overlapped by an event or flagged 6 by the scorer, minus the others,"
" with the labels shifted by k epochs (best of -3..3).")
print(clock.to_string(index=False))
fig, ax = plt.subplots(1, 2, figsize=(11, 3.5))
for sid, c in emg_curves.items():
ax[0].plot(emg_shifts * 5, c, lw=0.8)
for sid, c in swd_curves.items():
ax[1].plot(swd_shifts, c, lw=0.8)
ax[0].set(xlabel="shift of the sleep scores (s)", ylabel="EMG wake - NREM (log10)")
ax[1].set(xlabel="shift of the SWD events (s)", ylabel="5-10 Hz inside - outside (log10)")
plt.tight_layout(); plt.show()
print("plotted ranges: EMG", np.round([min(c.min() for c in emg_curves.values()), max(c.max() for c in emg_curves.values())], 3),
"| SWD", np.round([min(c.min() for c in swd_curves.values()), max(c.max() for c in swd_curves.values())], 3))
EMG: log10 60-90 Hz power (right minus left EMG), mean in wake minus mean in NREM epochs, with the scores shifted by k epochs. SWD: log10 5-10 Hz EEG power in 0.5-s bins, mean inside minus outside SWD events, with the events shifted by s seconds. events/flags_epoch: the same power per 5-s epoch, in epochs overlapped by an event or flagged 6 by the scorer, minus the others, with the labels shifted by k epochs (best of -3..3). subject emg_wake_minus_nrem_at_0 emg_best_shift_epochs emg_at_+-60 swd_bins_inside swd_contrast_at_0 swd_best_shift_s swd_contrast_at_+-10s events_epoch_best_shift events_epoch_contrast_at_0 flags_epoch_best_shift flags_epoch_contrast_at_0 GRIN2B-129 0.625 1 0.490 616 1.054 0.0 0.129 0.0 0.726 0.0 0.783 GRIN2B-130 1.234 0 0.841 1072 1.258 0.0 0.217 0.0 0.809 0.0 0.802 GRIN2B-131 0.804 0 0.555 1226 0.912 0.0 0.222 0.0 0.645 0.0 0.717 GRIN2B-132 0.851 0 0.668 78 0.985 0.0 -0.023 0.0 0.480 0.0 0.376 GRIN2B-137 0.801 0 0.517 384 1.392 0.0 0.322 0.0 0.802 0.0 0.855 GRIN2B-138 0.818 0 0.483 1722 1.348 0.0 0.400 0.0 0.771 0.0 0.762 GRIN2B-139 0.644 20 0.360 2996 0.671 0.0 0.153 0.0 0.452 0.0 0.478 GRIN2B-227 0.798 0 0.496 2212 1.250 0.0 0.510 0.0 0.920 0.0 0.917 GRIN2B-228 0.770 0 0.588 0 NaN NaN NaN NaN NaN NaN NaN GRIN2B-229 0.775 0 0.509 426 1.167 0.0 0.130 0.0 0.749 -1.0 0.658 GRIN2B-236 0.611 0 0.420 1576 1.260 0.0 0.771 0.0 1.097 -1.0 1.094 GRIN2B-237 0.679 0 0.487 960 0.932 -0.5 0.437 0.0 0.729 -1.0 0.651 GRIN2B-238 0.327 0 0.148 0 NaN NaN NaN NaN NaN NaN NaN GRIN2B-239 0.856 0 0.662 1162 1.205 0.0 0.400 0.0 0.897 0.0 0.925 GRIN2B-240 0.727 0 0.544 68 1.128 0.5 -0.054 0.0 0.511 0.0 0.455 GRIN2B-241 0.746 0 0.509 474 1.238 0.0 0.490 0.0 0.961 -1.0 0.989 GRIN2B-373 0.608 0 0.372 1997 1.282 0.0 0.758 0.0 1.068 0.0 1.108
plotted ranges: EMG [0.04 1.234] | SWD [-0.071 1.392]
For 15 of the 17 rats the EMG separates wake from NREM best with the scores where they are; for
GRIN2B-129 the best shift is one epoch and for GRIN2B-139 it is 20 epochs, where the curve is flat
(column emg_best_shift_epochs; bouts are long, so the contrast changes little over a few epochs,
and it is lower at 60 epochs in every rat). For every rat with enough SWD time in
the raw windows, the 5-10 Hz power inside SWD events is highest with the events where they are
and falls off within seconds. The scores and events share the clock of the raw recording in the
windows read. Rats with too few SWD bins in the windows (swd_bins_inside below 20) are not
checked for SWD timing.
At the 5-s scale, the epochs overlapped by SWD events contrast best with the labels unshifted in
every rat checked, and the flagged epochs at 0 or -1 epochs (flags_epoch_best_shift), with a
similar contrast at 0; the flags are on the same clock, consistent with their marking the end of
an SWD.
In [8]:
ids = list(by_id)
agree = []
for i, a in enumerate(ids):
for b in ids[i + 1:]:
for k in range(2):
sa, sb = by_id[a]["days"][k]["states"], by_id[b]["days"][k]["states"]
ta, tb = by_id[a]["days"][k]["sz_start"], by_id[b]["days"][k]["sz_start"]
agree.append({"pair": f"{a} {b}", "day": k + 1, "score_agreement": (sa == sb).mean(),
"identical_swd_times": len(np.intersect1d(np.round(ta, 2), np.round(tb, 2)))})
agree = pd.DataFrame(agree)
print("Highest agreement of sleep scores between two rats on the same day index:")
print(agree.sort_values("score_agreement", ascending=False).head(5).round(3).to_string(index=False))
n_ev = {sid: [len(d["sz_start"]) for d in by_id[sid]["days"]] for sid in ids}
agree["expected_by_chance"] = [n_ev[p.split()[0]][d - 1] * n_ev[p.split()[1]][d - 1] / (86400 / 0.2)
for p, d in zip(agree.pair, agree.day)]
print("Shared SWD start times (on the 0.2-s detection grid): largest", int(agree.identical_swd_times.max()),
"| largest expected by chance", round(agree.expected_by_chance.max(), 1),
"| pairs above 3 times chance and above 5:", int(((agree.identical_swd_times > 3 * agree.expected_by_chance) & (agree.identical_swd_times > 5)).sum()))
same_date = subjects.groupby("session_start").subject.apply(list)
rows = []
for date, group in same_date.items():
for i, a in enumerate(group):
for b in group[i + 1:]:
ea, eb = raw[a][0], raw[b][0]
m = ea.merge(eb, on=["bl", "epoch"], suffixes=("_a", "_b"))
m = m[~m.dropout_a & ~m.dropout_b]
rows.append({"pair": f"{a} {b}", "epochs": len(m),
"corr_emg_logpow": np.corrcoef(m.emg_logpow_a, m.emg_logpow_b)[0, 1] if len(m) > 10 else np.nan})
print("\nRats recorded together (same session start): correlation of per-epoch EMG power")
print(pd.DataFrame(rows).round(3).to_string(index=False))
Highest agreement of sleep scores between two rats on the same day index:
pair day score_agreement identical_swd_times
GRIN2B-129 GRIN2B-229 2 0.652 0
GRIN2B-137 GRIN2B-139 2 0.650 3
GRIN2B-129 GRIN2B-130 2 0.640 2
GRIN2B-129 GRIN2B-132 2 0.640 0
GRIN2B-132 GRIN2B-241 2 0.637 0
Shared SWD start times (on the 0.2-s detection grid): largest 5 | largest expected by chance 2.7 | pairs above 3 times chance and above 5: 0
Rats recorded together (same session start): correlation of per-epoch EMG power
pair epochs corr_emg_logpow
GRIN2B-129 GRIN2B-130 1828 0.340
GRIN2B-129 GRIN2B-131 1776 0.552
GRIN2B-129 GRIN2B-132 1546 0.267
GRIN2B-130 GRIN2B-131 2053 0.415
GRIN2B-130 GRIN2B-132 1791 0.441
GRIN2B-131 GRIN2B-132 1703 0.295
GRIN2B-137 GRIN2B-138 1197 0.497
GRIN2B-137 GRIN2B-139 1664 0.189
GRIN2B-138 GRIN2B-139 1384 0.286
GRIN2B-227 GRIN2B-228 1773 0.400
GRIN2B-227 GRIN2B-229 1907 0.350
GRIN2B-228 GRIN2B-229 1789 0.321
GRIN2B-237 GRIN2B-238 598 -0.193
GRIN2B-237 GRIN2B-239 1868 0.155
GRIN2B-237 GRIN2B-240 2109 0.354
GRIN2B-238 GRIN2B-239 582 0.222
GRIN2B-238 GRIN2B-240 703 0.076
GRIN2B-239 GRIN2B-240 1801 0.326
No two files share their scores beyond what two rats in the same light cycle do by chance, SWD start times coincide between files only as often as chance predicts on the 0.2-s detection grid, and the EMG power of rats recorded together is correlated (they share a room and light cycle) but far from identical. No recording appears to be deposited twice. GRIN2B-227 and GRIN2B-229 have the same number of samples, which follows from being recorded together.
Per animal, the measures below are the means over the two baseline days. Every genotype test is also shown without GRIN2B-228 (Het, 1 SWD per day), without GRIN2B-238 (WT, no SWDs) and without both, as the dandiset does not say whether they are among the paper's analysed rats. With GRIN2B-238 left out one WT remains: a rank-sum test then cannot reach p < 0.05 and a Welch test is not defined. Two more variants concern SWD epochs: removing only the SWD time from the states (event start to stop) where the main choice removes whole epochs, and forming bouts with SWD epochs dropped, so that an SWD does not split a bout.
In [9]:
def per_animal(mode="flag", bridge=False, exact=False):
# bridge: SWD epochs do not split bouts; exact: state minutes less the SWD time itself.
rows = []
for s in sessions:
days = []
for d in s["days"]:
st = grin2b.scored_states(d, mode)
m = grin2b.state_measures(st, bridge)
if exact:
m.update(grin2b.exact_minutes(d))
m.update({k: v for k, v in grin2b.swd_measures(d).items() if k != "swd_hourly"})
m["swd_epochs"] = int((st == 4).sum())
days.append(m)
r = pd.DataFrame(days).mean().to_dict()
r["subject"] = s["subject_id"]
r["sex"] = s["sex"]
r["genotype"] = s["genotype"]
rows.append(r)
return pd.DataFrame(rows).set_index("subject")
animals = {m: per_animal(m) for m in ("events", "flag", "stored")}
animals["bridged"] = per_animal("events", bridge=True)
animals["exact"] = per_animal("events", exact=True)
A = animals["events"]
cols = ["genotype", "sex", "swd_n", "swd_mean_s", "swd_total_min", "swd_epochs", "wake_min", "nrem_min", "rem_min",
"wake_bouts", "nrem_bouts", "rem_bouts", "wake_bout_s", "nrem_bout_s", "rem_bout_s"]
print(A[cols].sort_values(["genotype", "sex"]).round(2).to_string())
print("\nCounts:", A.genotype.fillna("none").value_counts().to_dict())
het = A[A.genotype == "Het"]
print(f"Het mean SWD time {het.swd_total_min.mean():.2f} min per day; SWD epochs removed from the states "
f"{het.swd_epochs.mean() * grin2b.EPOCH_S / 60:.2f} min per day; wake minutes with whole SWD epochs removed "
f"{het.wake_min.mean():.1f}, with only the SWD time removed {animals['exact'][animals['exact'].genotype == 'Het'].wake_min.mean():.1f}")
# Sensitivity: GRIN2B-228 (Het, 1 SWD per day) and GRIN2B-238 (WT, no SWDs) may not be among the
# paper's analysed rats (the Supplementary Methods excluded 4 Het and 9 WT of one founder line).
DROPS = {"without 228": ["GRIN2B-228"], "without 238": ["GRIN2B-238"],
"without 228 and 238": ["GRIN2B-228", "GRIN2B-238"]}
def variants(measure, table=None):
# The groups under the main choice and its variants: GRIN2B-132 as WT or as Het.
t = A if table is None else table
out = {"main (132 left out)": t[t.genotype.notna()].assign(g=t.genotype)}
for as_ in ("WT", "Het"):
out[f"132 as {as_}"] = t.assign(g=t.genotype.fillna(as_))
for name, ids in DROPS.items():
out[name] = out["main (132 left out)"][~out["main (132 left out)"].index.isin(ids)]
return out
def compare(measure, test, direction, extra_tables=()):
# The paper's two-group test under the main choice and each variant; direction "higher"/"lower" for Het.
rows = []
tabs = [("", A)] + list(extra_tables)
for tname, t in tabs:
for vname, d in variants(measure, t).items():
if tname and vname != "main (132 left out)":
continue
d = d[d[measure].notna()]
wt, het = d[d.g == "WT"][measure], d[d.g == "Het"][measure]
r = grin2b.two_group(wt, het, test)
diff = het.mean() - wt.mean()
rows.append({"variant": (tname + " " + vname).strip(), "n_WT": len(wt), "n_Het": len(het),
"WT_mean": wt.mean(), "Het_mean": het.mean(), "Het_minus_WT": diff,
"test": r["test"], "stat": r["stat"], "df": r["df"], "p": r["p"],
"paper_direction": (diff > 0) if direction == "higher" else (diff < 0)})
return pd.DataFrame(rows)
genotype sex swd_n swd_mean_s swd_total_min swd_epochs wake_min nrem_min rem_min wake_bouts nrem_bouts rem_bouts wake_bout_s nrem_bout_s rem_bout_s
subject
GRIN2B-137 Het F 586.5 3.00 29.69 908.0 900.88 433.46 30.00 628.0 210.0 22.5 88.94 124.90 80.64
GRIN2B-138 Het F 893.5 4.20 62.40 1590.5 808.17 456.46 42.83 692.5 205.5 41.0 70.18 134.65 63.84
GRIN2B-139 Het F 946.5 7.21 113.50 2249.0 797.75 417.33 37.50 747.0 256.0 36.0 64.53 98.40 62.74
GRIN2B-236 Het F 758.0 5.61 69.47 1500.5 841.00 444.83 29.12 507.5 225.5 36.5 100.33 120.05 47.77
GRIN2B-237 Het F 576.0 3.81 36.53 956.5 812.54 502.46 45.29 458.5 258.0 44.0 106.38 116.75 62.12
GRIN2B-239 Het F 884.0 4.44 65.80 1598.5 839.04 412.58 55.17 604.5 222.0 44.5 83.73 111.50 74.45
GRIN2B-241 Het F 698.5 3.83 45.04 1186.5 861.92 424.33 54.88 580.0 202.0 44.5 89.13 125.57 73.96
GRIN2B-129 Het M 784.5 4.70 61.59 1445.5 828.12 476.33 15.08 577.5 143.0 14.0 86.70 202.99 64.39
GRIN2B-130 Het M 507.0 4.44 37.48 938.5 862.54 444.88 54.38 533.0 218.5 39.0 97.19 122.16 83.58
GRIN2B-131 Het M 816.0 4.60 62.65 1513.0 753.29 495.08 65.54 639.0 205.0 50.0 70.71 144.98 83.02
GRIN2B-227 Het M 957.5 4.97 79.39 1827.5 765.75 453.75 68.21 656.0 241.5 55.5 70.11 112.86 74.44
GRIN2B-228 Het M 1.0 0.80 0.01 1.0 969.12 428.25 42.54 206.0 214.5 31.0 282.53 120.40 89.37
GRIN2B-229 Het M 443.5 2.69 20.37 659.5 883.08 452.54 49.42 515.0 202.5 41.0 104.44 134.44 72.23
GRIN2B-373 Het M 816.0 9.31 124.87 2200.5 754.00 464.75 37.88 487.0 255.0 30.5 92.98 109.29 74.52
GRIN2B-238 WT F 0.0 NaN 0.00 0.0 1015.33 370.42 54.25 230.5 253.0 60.5 265.89 88.83 55.44
GRIN2B-240 WT F 112.5 1.47 2.75 147.0 874.42 468.12 85.21 358.5 293.0 87.5 146.47 95.71 58.17
GRIN2B-132 NaN U 107.0 1.70 3.04 139.5 1023.38 367.75 37.25 272.5 171.5 23.5 226.74 129.56 95.52
Counts: {'Het': 14, 'WT': 2, 'none': 1}
Het mean SWD time 57.77 min per day; SWD epochs removed from the states 110.57 min per day; wake minutes with whole SWD epochs removed 834.1, with only the SWD time removed 884.1
fig2b-swd-hourly¶
Grin2b+/- rats have more SWDs than WT across the 24-h day; hour of day has no effect and does not interact with genotype (linear mixed model, genotype F = 55.36, p < .00001; hour p = .16; interaction p = .49).
In [10]:
def hourly_table(fn):
rows = []
for s in sessions:
v = np.mean([fn(d) for d in s["days"]], axis=0)
for h in range(len(v)):
rows.append({"subject": s["subject_id"], "hour": h, "value": v[h], "genotype": s["genotype"]})
return pd.DataFrame(rows)
def hourly_tests(H, within="hour"):
out = []
for vname, g in (("main (132 left out)", H[H.genotype.notna()]),
("132 as WT", H.assign(genotype=H.genotype.fillna("WT"))),
("132 as Het", H.assign(genotype=H.genotype.fillna("Het"))),
*((name, H[H.genotype.notna() & ~H.subject.isin(ids)]) for name, ids in DROPS.items())):
a = grin2b.split_plot_anova(g, "value", "subject", "genotype", within)
a.insert(0, "variant", vname)
out.append(a)
return pd.concat(out, ignore_index=True)
def mixed_check(H, within="hour"):
# The paper's model: a linear mixed model with a random intercept per animal (REML), Wald tests
# of each term with sum contrasts, on the main animals; compared with the split-plot F tests.
g = H[H.genotype.notna()]
fit = sm.MixedLM.from_formula(f"value ~ C(genotype, Sum) * C({within}, Sum)", data=g, groups=g.subject).fit(reml=True)
w = fit.wald_test_terms(scalar=True).table
w = w[w.index != "Intercept"]
return pd.DataFrame({"term": w.index, "chi2_over_df": w.statistic.to_numpy() / w.df_constraint.to_numpy(),
"df": w.df_constraint.to_numpy(), "p": w.pvalue.to_numpy(), "converged": fit.converged})
H_swd = hourly_table(lambda d: grin2b.swd_measures(d)["swd_hourly"])
print("Mean SWDs per hour by genotype (main):")
print(H_swd[H_swd.genotype.notna()].groupby(["genotype", "hour"]).value.mean().unstack().round(1).to_string())
print(hourly_tests(H_swd).round(5).to_string(index=False))
print("Mixed model (MixedLM, random intercept), main animals:")
print(mixed_check(H_swd).round(5).to_string(index=False))
Mean SWDs per hour by genotype (main):
hour 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
genotype
Het 30.1 29.7 32.0 24.4 21.1 29.1 25.9 26.1 19.2 21.5 26.3 31.2 31.4 35.2 33.1 37.3 29.2 37.8 33.8 35.9 21.8 21.2 26.2 31.0
WT 1.2 2.0 0.5 0.5 4.0 1.2 0.0 0.2 1.0 0.8 0.0 1.0 0.5 1.0 1.8 0.0 2.2 3.2 7.8 8.2 7.2 7.8 1.0 3.0
variant effect F df df_error p
main (132 left out) genotype 11.45361 1 14 0.00445
main (132 left out) hour 1.49658 23 322 0.06880
main (132 left out) genotype x hour 0.26633 23 322 0.99979
132 as WT genotype 16.38041 1 15 0.00105
132 as WT hour 1.60237 23 345 0.04062
132 as WT genotype x hour 0.27760 23 345 0.99970
132 as Het genotype 7.96273 1 15 0.01288
132 as Het hour 1.60156 23 345 0.04078
132 as Het genotype x hour 0.26995 23 345 0.99976
without 228 genotype 30.54028 1 13 0.00010
without 228 hour 1.49553 23 299 0.06987
without 228 genotype x hour 0.27887 23 299 0.99968
without 238 genotype 4.74580 1 13 0.04837
without 238 hour 1.48863 23 299 0.07216
without 238 genotype x hour 0.21056 23 299 0.99997
without 228 and 238 genotype 12.96719 1 12 0.00364
without 228 and 238 hour 1.48544 23 276 0.07410
without 228 and 238 genotype x hour 0.21097 23 276 0.99997
Mixed model (MixedLM, random intercept), main animals:
term chi2_over_df df p converged
C(genotype, Sum) 11.45242 1 0.00071 True
C(hour, Sum) 0.26326 23 0.99984 True
C(genotype, Sum):C(hour, Sum) 0.26633 23 0.99983 True
Result: reproduced. Het have more SWDs per hour than WT (genotype F = 11.45, p = 0.0045), hour has no effect at p < 0.05 (p = 0.069), and it does not interact with genotype (p = 0.9998). Het means per hour range from 19.2 to 37.8, WT means from 0 to 8.2. With GRIN2B-132 as WT or as Het the genotype effect holds (p = 0.0011 and 0.013), and hour reaches p = 0.041 in both. Without GRIN2B-228, without GRIN2B-238 or without both, genotype p = 0.0001, 0.048 and 0.0036, hour p = 0.070 to 0.074 and interaction p >= 0.9997. A random-intercept mixed model (MixedLM) gives the same genotype and interaction statistics (genotype 11.45, Wald p = 0.0007; interaction p = 0.9998); its hour term (p = 0.9998) differs from the split-plot hour term, as it weighs the two groups equally. The comparison rests on two WT rats.
fig2c-swd-number¶
Grin2b+/- rats have more SWDs over 24 h than WT (Wilcoxon rank-sum, W = 4.5, p = .0002).
In [11]:
print(A[["genotype", "swd_n"]].sort_values(["genotype", "swd_n"]).round(3).to_string())
print("Animals without a value (no SWDs):", A[A["swd_n"].isna()].index.tolist())
print(compare("swd_n", "wilcoxon", "higher").round(5).to_string(index=False))
genotype swd_n
subject
GRIN2B-228 Het 1.0
GRIN2B-229 Het 443.5
GRIN2B-130 Het 507.0
GRIN2B-237 Het 576.0
GRIN2B-137 Het 586.5
GRIN2B-241 Het 698.5
GRIN2B-236 Het 758.0
GRIN2B-129 Het 784.5
GRIN2B-131 Het 816.0
GRIN2B-373 Het 816.0
GRIN2B-239 Het 884.0
GRIN2B-138 Het 893.5
GRIN2B-139 Het 946.5
GRIN2B-227 Het 957.5
GRIN2B-238 WT 0.0
GRIN2B-240 WT 112.5
GRIN2B-132 NaN 107.0
Animals without a value (no SWDs): []
variant n_WT n_Het WT_mean Het_mean Het_minus_WT test stat df p paper_direction
main (132 left out) 2 14 56.25000 690.60714 634.35714 Wilcoxon rank-sum (W) 1.0 NaN 0.04701 True
132 as WT 3 14 73.16667 690.60714 617.44048 Wilcoxon rank-sum (W) 2.0 NaN 0.01969 True
132 as Het 2 15 56.25000 651.70000 595.45000 Wilcoxon rank-sum (W) 2.0 NaN 0.06225 True
without 228 2 13 56.25000 743.65385 687.40385 Wilcoxon rank-sum (W) 0.0 NaN 0.03359 True
without 238 1 14 112.50000 690.60714 578.10714 Wilcoxon rank-sum (W) 1.0 NaN 0.20261 True
without 228 and 238 1 13 112.50000 743.65385 631.15385 Wilcoxon rank-sum (W) 0.0 NaN 0.13621 True
Result: reproduced. Het have 690.6 SWDs per day on average against 56.3 in WT (W = 1, p = 0.047). GRIN2B-238 has none; GRIN2B-240 (112.5 per day) ranks above one Het, GRIN2B-228 (1 per day). With GRIN2B-132 as WT p = 0.020, as Het p = 0.062. Without GRIN2B-228 p = 0.034 (W = 0); with GRIN2B-238 left out one WT remains and p = 0.20 (0.14 without both).
fig2c-swd-mean-duration¶
SWDs are longer on average in Grin2b+/- rats than in WT (two-sample t test, df = 22, T = -5.26, p = .00003).
In [12]:
print(A[["genotype", "swd_mean_s"]].sort_values(["genotype", "swd_mean_s"]).round(3).to_string())
print("Animals without a value (no SWDs):", A[A["swd_mean_s"].isna()].index.tolist())
print(compare("swd_mean_s", "t", "higher").round(5).to_string(index=False))
genotype swd_mean_s
subject
GRIN2B-228 Het 0.800
GRIN2B-229 Het 2.687
GRIN2B-137 Het 3.001
GRIN2B-237 Het 3.814
GRIN2B-241 Het 3.829
GRIN2B-138 Het 4.200
GRIN2B-130 Het 4.438
GRIN2B-239 Het 4.443
GRIN2B-131 Het 4.601
GRIN2B-129 Het 4.699
GRIN2B-227 Het 4.973
GRIN2B-236 Het 5.614
GRIN2B-139 Het 7.206
GRIN2B-373 Het 9.308
GRIN2B-240 WT 1.472
GRIN2B-238 WT NaN
GRIN2B-132 NaN 1.703
Animals without a value (no SWDs): ['GRIN2B-238']
variant n_WT n_Het WT_mean Het_mean Het_minus_WT test stat df p paper_direction
main (132 left out) 1 14 1.47217 4.54377 3.07160 t -1.48132 13 0.16235 True
132 as WT 2 14 1.58766 4.54377 2.95610 t -2.02528 14 0.06234 True
132 as Het 1 15 1.47217 4.35440 2.88223 t -1.35142 14 0.19800 True
without 228 1 13 1.47217 4.83175 3.35959 t -1.84180 12 0.09034 True
without 238 1 14 1.47217 4.54377 3.07160 t -1.48132 13 0.16235 True
without 228 and 238 1 13 1.47217 4.83175 3.35959 t -1.84180 12 0.09034 True
Result: partially reproduced. Het SWDs last 4.54 s on average against 1.47 s in the only WT rat with SWDs (GRIN2B-240), the paper's direction, but the t test with one WT gives p = 0.16. With GRIN2B-132 (1.70 s) as WT p = 0.062, and without GRIN2B-228 p = 0.090. A t test against one WT value says little either way.
fig2c-swd-total-duration¶
The total time in SWDs over 24 h is higher in Grin2b+/- rats than in WT (Wilcoxon rank-sum, W = 5.5, p = .0002).
In [13]:
print(A[["genotype", "swd_total_min"]].sort_values(["genotype", "swd_total_min"]).round(3).to_string())
print("Animals without a value (no SWDs):", A[A["swd_total_min"].isna()].index.tolist())
print(compare("swd_total_min", "wilcoxon", "higher").round(5).to_string(index=False))
genotype swd_total_min
subject
GRIN2B-228 Het 0.013
GRIN2B-229 Het 20.372
GRIN2B-137 Het 29.690
GRIN2B-237 Het 36.535
GRIN2B-130 Het 37.477
GRIN2B-241 Het 45.045
GRIN2B-129 Het 61.587
GRIN2B-138 Het 62.403
GRIN2B-131 Het 62.648
GRIN2B-239 Het 65.797
GRIN2B-236 Het 69.467
GRIN2B-227 Het 79.393
GRIN2B-139 Het 113.502
GRIN2B-373 Het 124.870
GRIN2B-238 WT 0.000
GRIN2B-240 WT 2.747
GRIN2B-132 NaN 3.037
Animals without a value (no SWDs): []
variant n_WT n_Het WT_mean Het_mean Het_minus_WT test stat df p paper_direction
main (132 left out) 2 14 1.37333 57.77131 56.39798 Wilcoxon rank-sum (W) 1.0 NaN 0.03333 True
132 as WT 3 14 1.92778 57.77131 55.84353 Wilcoxon rank-sum (W) 2.0 NaN 0.01176 True
132 as Het 2 15 1.37333 54.12233 52.74900 Wilcoxon rank-sum (W) 1.0 NaN 0.02941 True
without 228 2 13 1.37333 62.21423 60.84090 Wilcoxon rank-sum (W) 0.0 NaN 0.01905 True
without 238 1 14 2.74667 57.77131 55.02464 Wilcoxon rank-sum (W) 1.0 NaN 0.26667 True
without 228 and 238 1 13 2.74667 62.21423 59.46756 Wilcoxon rank-sum (W) 0.0 NaN 0.14286 True
Result: reproduced. Het spend 57.77 min per day in SWDs against 1.37 min in WT (W = 1, p = 0.033). With GRIN2B-132 as WT p = 0.012, as Het p = 0.029. Without GRIN2B-228 p = 0.019; with GRIN2B-238 left out one WT remains and p = 0.27 (0.14 without both).
disc-het-swd-minutes¶
Over 24 h, Grin2b+/- rats spend on average 59.97 min in SWDs, while WT littermates average less than a minute.
In [14]:
for vname, d in variants("swd_total_min").items():
g = d.groupby("g").swd_total_min
print(f"{vname}: Het mean {g.mean()['Het']:.2f} min (n {g.size()['Het']}, SEM {g.sem()['Het']:.2f}), "
f"WT mean {g.mean()['WT']:.2f} min (n {g.size()['WT']}; values {sorted(d[d.g == 'WT'].swd_total_min.round(2))})")
print("Expected Het 59.97 +- 6.0 min; WT below 1 min.")
print("Het mean without GRIN2B-228:", round(A[(A.genotype == "Het") & (A.index != "GRIN2B-228")].swd_total_min.mean(), 2))
print("SWD time from the epochs SWD events overlap (5 s each), Het mean:", round(A[A.genotype == "Het"].swd_epochs.mean() * 5 / 60, 2), "min")
F = animals["flag"]
print("SWD time from the epochs the scorer flags (5 s each), Het mean:", round(F[F.genotype == "Het"].swd_epochs.mean() * 5 / 60, 2), "min")
main (132 left out): Het mean 57.77 min (n 14, SEM 9.04), WT mean 1.37 min (n 2; values [0.0, 2.75]) 132 as WT: Het mean 57.77 min (n 14, SEM 9.04), WT mean 1.93 min (n 3; values [0.0, 2.75, 3.04]) 132 as Het: Het mean 54.12 min (n 15, SEM 9.18), WT mean 1.37 min (n 2; values [0.0, 2.75]) without 228: Het mean 62.21 min (n 13, SEM 8.51), WT mean 1.37 min (n 2; values [0.0, 2.75]) without 238: Het mean 57.77 min (n 14, SEM 9.04), WT mean 2.75 min (n 1; values [2.75]) without 228 and 238: Het mean 62.21 min (n 13, SEM 8.51), WT mean 2.75 min (n 1; values [2.75]) Expected Het 59.97 +- 6.0 min; WT below 1 min. Het mean without GRIN2B-228: 62.21 SWD time from the epochs SWD events overlap (5 s each), Het mean: 110.57 min SWD time from the epochs the scorer flags (5 s each), Het mean: 60.2 min
Result: partially reproduced. The Het mean is 57.77 min per day (SEM 9.04, 14 rats), within the tolerance of 59.97 +- 6.0 min. The WT mean is 1.37 min (0 and 2.75 min), above the paper's less than a minute. Without GRIN2B-228 the Het mean is 62.21 min, also within tolerance. Without GRIN2B-238 the WT mean is 2.75 min, so the WT part fails with either WT set. Counting each epoch that an SWD event overlaps as 5 s of SWD gives 110.57 min, and counting each epoch the scorer flags gives 60.2 min, so the paper's number fits both the event durations and the flags.
figs5a-swd-light-dark¶
The number of SWDs does not differ between the light and dark phases, and Grin2b+/- rats have more SWDs than WT in both (linear mixed model, genotype p < .00001, phase p = .28, interaction p = .53).
In [15]:
H_ld = hourly_table(lambda d: np.array([grin2b.swd_measures(d)["swd_light"], grin2b.swd_measures(d)["swd_dark"]]))
H_ld["phase"] = np.where(H_ld.hour == 0, "light", "dark")
print(H_ld[H_ld.genotype.notna()].pivot_table(index="subject", columns="phase", values="value").join(A.genotype).round(1).to_string())
print(H_ld[H_ld.genotype.notna()].groupby(["genotype", "phase"]).value.mean().round(1).to_string())
print(hourly_tests(H_ld.drop(columns="hour"), within="phase").round(5).to_string(index=False))
het = H_ld[H_ld.genotype == "Het"].pivot_table(index="subject", columns="phase", values="value")
print("Het light minus dark: median", round((het.light - het.dark).median(), 1), "| Wilcoxon signed-rank p",
round(stats.wilcoxon(het.light, het.dark).pvalue, 4))
dark light genotype
subject
GRIN2B-129 532.0 252.5 Het
GRIN2B-130 281.0 226.0 Het
GRIN2B-131 584.0 232.0 Het
GRIN2B-137 300.5 286.0 Het
GRIN2B-138 582.0 311.5 Het
GRIN2B-139 440.0 506.5 Het
GRIN2B-227 574.0 383.5 Het
GRIN2B-228 1.0 0.0 Het
GRIN2B-229 274.0 169.5 Het
GRIN2B-236 443.5 314.5 Het
GRIN2B-237 107.5 468.5 Het
GRIN2B-238 0.0 0.0 WT
GRIN2B-239 286.5 597.5 Het
GRIN2B-240 87.5 25.0 WT
GRIN2B-241 418.5 280.0 Het
GRIN2B-373 410.5 405.5 Het
genotype phase
Het dark 373.9
light 316.7
WT dark 43.8
light 12.5
variant effect F df df_error p
main (132 left out) genotype 11.45361 1 14 0.00445
main (132 left out) phase 1.18972 1 14 0.29380
main (132 left out) genotype x phase 0.03017 1 14 0.86460
132 as WT genotype 16.38041 1 15 0.00105
132 as WT phase 1.33627 1 15 0.26577
132 as WT genotype x phase 0.02813 1 15 0.86905
132 as Het genotype 7.96273 1 15 0.01288
132 as Het phase 1.33653 1 15 0.26573
132 as Het genotype x phase 0.03106 1 15 0.86246
without 228 genotype 30.54028 1 13 0.00010
without 228 phase 1.18301 1 13 0.29650
without 228 genotype x phase 0.03798 1 13 0.84849
without 238 genotype 4.74580 1 13 0.04837
without 238 phase 1.18260 1 13 0.29658
without 238 genotype x phase 0.00061 1 13 0.98065
without 228 and 238 genotype 12.96719 1 12 0.00364
without 228 and 238 phase 1.17421 1 12 0.29983
without 228 and 238 genotype x phase 0.00002 1 12 0.99673
Het light minus dark: median -79.8 | Wilcoxon signed-rank p 0.1937
Result: reproduced. Phase has no effect (F = 1.19, p = 0.29) and genotype has (F = 11.45, p = 0.0045), with Het higher in both phases (316.7 SWDs in the light and 373.9 in the dark, against 12.5 and 43.8 in WT) and no interaction (p = 0.86). Within Het, light minus dark has a median of -79.8 SWDs (signed-rank p = 0.19). The variants with GRIN2B-132 give phase p = 0.27. Without GRIN2B-228, GRIN2B-238 or both, phase p = 0.30 and genotype p = 0.0001, 0.048 and 0.0036.
figs8-swd-sex¶
Sex does not affect the number, mean duration or total duration of SWDs (two-way ANOVA sex x genotype, sex p = .74, .99 and .99).
In [16]:
d = A[A.genotype.notna() & A.sex.isin(["M", "F"])].copy()
print(pd.crosstab(d.genotype, d.sex).to_string())
rows = []
for measure in ("swd_n", "swd_mean_s", "swd_total_min"):
dd = d[d[measure].notna()]
fit = ols(f"{measure} ~ C(genotype) + C(sex)", data=dd).fit()
a = sm.stats.anova_lm(fit, typ=2)
het = dd[dd.genotype == "Het"]
t = stats.ttest_ind(het[het.sex == "M"][measure], het[het.sex == "F"][measure])
rows.append({"measure": measure, "n": len(dd), "sex_F": a.loc["C(sex)", "F"], "sex_p": a.loc["C(sex)", "PR(>F)"],
"Het_male_mean": het[het.sex == "M"][measure].mean(), "Het_female_mean": het[het.sex == "F"][measure].mean(),
"Het_only_t_p": t.pvalue})
print(pd.DataFrame(rows).round(4).to_string(index=False))
for name, ids in DROPS.items():
dd0 = d[~d.index.isin(ids)]
ps = []
for measure in ("swd_n", "swd_mean_s", "swd_total_min"):
dd = dd0[dd0[measure].notna()]
ps.append(round(sm.stats.anova_lm(ols(f"{measure} ~ C(genotype) + C(sex)", data=dd).fit(), typ=2).loc["C(sex)", "PR(>F)"], 3))
print(f"{name}: sex p for number, mean duration, total duration = {ps}")
sex F M genotype Het 7 7 WT 2 0
measure n sex_F sex_p Het_male_mean Het_female_mean Het_only_t_p
swd_n 16 1.2218 0.2890 617.9286 763.2857 0.3073
swd_mean_s 15 0.0059 0.9398 4.5008 4.5867 0.9398
swd_total_min 16 0.0817 0.7795 55.1943 60.3483 0.7883
without 228: sex p for number, mean duration, total duration = [np.float64(0.659), np.float64(0.609), np.float64(0.816)]
without 238: sex p for number, mean duration, total duration = [np.float64(0.307), np.float64(0.94), np.float64(0.788)]
without 228 and 238: sex p for number, mean duration, total duration = [np.float64(0.67), np.float64(0.609), np.float64(0.824)]
Result: reproduced. The sex term has p = 0.29, 0.94 and 0.78 for the number, mean duration and total duration of SWDs. Within Het, females have 763.3 SWDs per day against 617.9 in males (t test p = 0.31), 4.59 against 4.50 s (p = 0.94) and 60.3 against 55.2 min (p = 0.79). The ANOVA has no interaction term, as there are no WT males. Without GRIN2B-228, GRIN2B-238 or both the sex terms have p from 0.31 to 0.94.
Spectra from the raw EEG¶
The Methods' spectra, on the epochs of the raw windows: usable epochs only, each epoch's Welch spectrum divided by the animal's mean spectrum over its usable wake, NREM and REM epochs, then averaged per state (SWD epochs are those SWD events overlap). Absolute spectra are kept for the SWD peak. Variants: SWD epochs from the scorer's flags; the spectral-slope rejection of the Methods (slope of log10 power against log10 frequency over 1 to 48 Hz below -5, or power at 1 Hz above 500 mV^2/Hz, a reading of the stated units chosen here); and, per animal, whichever S1 channel (right or left) has fewer epochs above the artefact threshold, as the Methods chose the cleaner channel by eye for sleep scoring.
In [17]:
def slope_bad(e):
return ((e.slope < grin2b.SLOPE_MIN) | (10 ** e.offset > grin2b.OFFSET_MAX)).to_numpy()
def artefact_rate(e):
ok = ~e.dropout
return (e.artefact[ok]).mean()
# Per animal, the S1 channel with fewer artefact epochs (ties keep S1_Tr_R).
cleaner = {sid: ("S1_Tr_L" if artefact_rate(raw[sid][3]) < artefact_rate(raw[sid][0]) else "S1_Tr_R") for sid in raw}
def state_spectra(mode="events", slope_reject=False, channel="main"):
out = {}
for s in sessions:
sid = s["subject_id"]
e = raw[sid][3] if channel == "cleaner" and cleaner[sid] == "S1_Tr_L" else raw[sid][0]
days = {d["bl"]: d for d in s["days"]}
sc = {bl: grin2b.scored_states(d, mode) for bl, d in days.items()}
st = np.array([sc[bl][k] for bl, k in zip(e.bl, e.epoch)])
use = (~e.dropout & ~e.artefact).to_numpy()
if slope_reject:
use = use & ~slope_bad(e)
P = np.vstack(e.psd.to_numpy())
base = P[use & np.isin(st, [0, 1, 2])].mean(0)
r = {"n": {}}
for c, name in grin2b.STATES.items():
m = use & (st == c)
r["n"][name] = int(m.sum())
r[name] = (P[m] / base).mean(0) if m.any() else None
r[name + "_abs"] = P[m].mean(0) if m.any() else None
out[sid] = r
return out
spectra = state_spectra("events")
spectra_flag = state_spectra("flag")
spectra_slope = state_spectra("events", slope_reject=True)
spectra_clean = state_spectra("events", channel="cleaner")
print("Usable epochs per state in the raw windows:")
print(pd.DataFrame({sid: r["n"] for sid, r in spectra.items()}).T.join(A.genotype).to_string())
# Coverage: usable epochs of each state in the windows, as a fraction of that state's epochs over both days.
rows = []
for s in sessions:
sid = s["subject_id"]
e = raw[sid][0]
tot = {c: sum(int((grin2b.scored_states(d, "events") == c).sum()) for d in s["days"]) for c in (0, 1, 2, 4)}
use = ~e.dropout & ~e.artefact
row = {"subject": sid, "genotype": s["genotype"], "window_epochs_frac_of_2_days": round(len(e) / (2 * grin2b.EPOCHS_PER_DAY), 3)}
for c, name in ((0, "Wake"), (1, "NREM"), (2, "REM"), (4, "SWD")):
row[f"{name}_covered"] = round(spectra[sid]["n"][name] / tot[c], 3) if tot[c] else np.nan
row["slope_rule_removes"] = int((use & slope_bad(e)).sum())
row["slope_min"] = round(e.slope[use].min(), 2)
row["offset_1Hz_max_mV2_per_Hz"] = float(f"{(10 ** e.offset[use]).max():.3g}")
row["artefact_rate_R"] = round(artefact_rate(raw[sid][0]), 3)
row["artefact_rate_L"] = round(artefact_rate(raw[sid][3]), 3)
row["cleaner_channel"] = cleaner[sid]
rows.append(row)
coverage = pd.DataFrame(rows)
print("\nCoverage of each state by the usable epochs of the windows, the slope rule and the cleaner S1 channel:")
print(coverage.to_string(index=False))
print("Median coverage:", coverage[[c for c in coverage if c.endswith("_covered")]].median().round(3).to_dict(),
"| epochs the slope rule removes, all animals:", int(coverage.slope_rule_removes.sum()),
"| animals where S1_Tr_L is cleaner:", coverage[coverage.cleaner_channel == "S1_Tr_L"].subject.tolist())
def band_table(state, bands=grin2b.BANDS, spec=None):
spec = spectra if spec is None else spec
rows = []
for sid, r in spec.items():
if r[state] is None:
continue
for b, v in grin2b.band_power(FREQ, r[state], bands).items():
rows.append({"subject": sid, "band": b, "value": v, "genotype": by_id[sid]["genotype"]})
return pd.DataFrame(rows)
def stored_band_table(state, bands=grin2b.BANDS):
rows = []
for s in sessions:
sp = s["spectra"]
if state not in sp:
continue
g = sp.groupby("frequency_hz")[state].mean()
f = g.index.to_numpy()
for b, (lo, hi) in bands.items():
rows.append({"subject": s["subject_id"], "band": b, "value": g[(f >= lo) & (f < hi)].mean(),
"genotype": s["genotype"]})
return pd.DataFrame(rows)
def band_anova(T, name):
rows = []
for vname, d in (("main (132 left out)", T[T.genotype.notna()]),
("132 as WT", T.assign(genotype=T.genotype.fillna("WT"))),
("132 as Het", T.assign(genotype=T.genotype.fillna("Het"))),
*((name, T[T.genotype.notna() & ~T.subject.isin(ids)]) for name, ids in DROPS.items())):
d = d.dropna(subset=["value"])
fit = ols("value ~ C(genotype) * C(band)", data=d).fit()
a = sm.stats.anova_lm(fit, typ=2)
tk = pairwise_tukeyhsd(d.value, d.genotype + ":" + d.band)
tk = pd.DataFrame(tk.summary().data[1:], columns=tk.summary().data[0])
row = {"spectra": name, "variant": vname, "n_WT": d[d.genotype == "WT"].subject.nunique(),
"n_Het": d[d.genotype == "Het"].subject.nunique(),
"genotype_F": a.loc["C(genotype)", "F"], "genotype_p": a.loc["C(genotype)", "PR(>F)"],
"interaction_F": a.loc["C(genotype):C(band)", "F"], "interaction_p": a.loc["C(genotype):C(band)", "PR(>F)"]}
means = d.groupby(["genotype", "band"]).value.mean()
row["mean_Het-WT"] = np.mean([means.get(("Het", b), np.nan) - means.get(("WT", b), np.nan) for b in grin2b.BANDS])
for b in grin2b.BANDS:
pair = tk[((tk.group1 == f"Het:{b}") & (tk.group2 == f"WT:{b}")) | ((tk.group1 == f"WT:{b}") & (tk.group2 == f"Het:{b}"))]
row[f"{b}_Het-WT"] = means.get(("Het", b), np.nan) - means.get(("WT", b), np.nan)
row[f"{b}_tukey_p"] = float(pair["p-adj"].iloc[0]) if len(pair) else np.nan
rows.append(row)
return pd.DataFrame(rows)
def spectral_claim(state):
res = [band_anova(band_table(state), "Methods, raw"),
band_anova(band_table(state, grin2b.BANDS_ALT), "Methods, raw, contiguous bands"),
band_anova(band_table(state, spec=spectra_flag), "Methods, raw, SWD epochs from scorer flags"),
band_anova(band_table(state, spec=spectra_slope), "Methods, raw, slope rejection"),
band_anova(band_table(state, spec=spectra_clean), "Methods, raw, cleaner S1 channel"),
band_anova(stored_band_table(state), "stored (scorer)")]
return pd.concat(res, ignore_index=True)
Usable epochs per state in the raw windows:
Wake NREM REM SWD genotype
GRIN2B-129 975 1248 4 86 Het
GRIN2B-130 1644 675 107 253 Het
GRIN2B-131 842 1379 177 185 Het
GRIN2B-132 1221 906 65 20 NaN
GRIN2B-137 990 1027 74 79 Het
GRIN2B-138 755 585 137 292 Het
GRIN2B-139 914 1089 117 397 Het
GRIN2B-227 1309 659 116 427 Het
GRIN2B-228 1341 892 126 0 Het
GRIN2B-229 1153 1257 83 99 Het
GRIN2B-236 1088 594 24 284 Het
GRIN2B-237 1013 1301 161 183 Het
GRIN2B-238 370 470 97 0 WT
GRIN2B-239 954 960 213 202 Het
GRIN2B-240 1436 1069 140 18 WT
GRIN2B-241 977 843 166 97 Het
GRIN2B-373 967 597 69 275 Het
Coverage of each state by the usable epochs of the windows, the slope rule and the cleaner S1 channel:
subject genotype window_epochs_frac_of_2_days Wake_covered NREM_covered REM_covered SWD_covered slope_rule_removes slope_min offset_1Hz_max_mV2_per_Hz artefact_rate_R artefact_rate_L cleaner_channel
GRIN2B-129 Het 0.099 0.049 0.109 0.011 0.030 0 -2.72 0.0588 0.000 0.000 S1_Tr_R
GRIN2B-130 Het 0.099 0.079 0.063 0.082 0.135 0 -3.11 0.7640 0.011 0.000 S1_Tr_L
GRIN2B-131 Het 0.099 0.047 0.116 0.113 0.061 0 -3.62 4.0000 0.005 0.003 S1_Tr_L
GRIN2B-132 NaN 0.099 0.050 0.103 0.073 0.072 0 -2.64 0.0754 0.000 0.000 S1_Tr_R
GRIN2B-137 Het 0.099 0.046 0.099 0.103 0.044 0 -2.90 0.1830 0.000 0.000 S1_Tr_L
GRIN2B-138 Het 0.099 0.039 0.053 0.133 0.092 0 -2.85 0.1030 0.000 0.000 S1_Tr_R
GRIN2B-139 Het 0.099 0.048 0.109 0.130 0.088 0 -3.76 5.3400 0.037 0.009 S1_Tr_L
GRIN2B-227 Het 0.099 0.071 0.061 0.071 0.117 0 -2.72 0.0752 0.000 0.000 S1_Tr_R
GRIN2B-228 Het 0.099 0.058 0.087 0.123 0.000 0 -2.63 0.0186 0.000 0.000 S1_Tr_R
GRIN2B-229 Het 0.099 0.054 0.116 0.070 0.075 0 -2.83 0.0956 0.000 0.000 S1_Tr_R
GRIN2B-236 Het 0.099 0.054 0.056 0.034 0.095 0 -2.74 0.0922 0.000 0.000 S1_Tr_R
GRIN2B-237 Het 0.099 0.052 0.108 0.148 0.096 0 -2.73 0.0974 0.000 0.000 S1_Tr_R
GRIN2B-238 WT 0.099 0.015 0.053 0.075 NaN 0 -2.76 0.0615 0.001 0.000 S1_Tr_L
GRIN2B-239 Het 0.099 0.047 0.097 0.161 0.063 0 -2.73 0.0809 0.000 0.000 S1_Tr_L
GRIN2B-240 WT 0.099 0.068 0.095 0.068 0.061 0 -2.68 0.0658 0.002 0.000 S1_Tr_L
GRIN2B-241 Het 0.099 0.047 0.083 0.126 0.041 0 -2.72 0.0520 0.000 0.000 S1_Tr_R
GRIN2B-373 Het 0.099 0.053 0.054 0.076 0.062 0 -2.71 0.0735 0.000 0.000 S1_Tr_R
Median coverage: {'Wake_covered': 0.05, 'NREM_covered': 0.095, 'REM_covered': 0.082, 'SWD_covered': 0.068} | epochs the slope rule removes, all animals: 0 | animals where S1_Tr_L is cleaner: ['GRIN2B-130', 'GRIN2B-131', 'GRIN2B-137', 'GRIN2B-139', 'GRIN2B-238', 'GRIN2B-239', 'GRIN2B-240']
The usable epochs of the windows cover a median of 5.0% of each rat's wake epochs over the two days, 9.5% of NREM, 8.2% of REM and 6.8% of SWD epochs; the lowest are GRIN2B-129's REM (1.1%, 4 epochs) and GRIN2B-238's wake (1.5%). The slope rule, read as above, removes no epoch: the steepest slope is -3.76 and the largest power at 1 Hz is 5.34 mV^2/Hz. S1_Tr_L has fewer artefact epochs than S1_Tr_R in 7 rats, among them both WT; the artefact rate is at most 3.7% on either channel. The band ANOVAs below treat the five bands of an animal as separate observations, as the paper's two-way ANOVA does, so their p-values are smaller than a test on animals would give.
fig2d-swd-spectrum-peak¶
The EEG spectrum during SWDs has a clear peak in the theta range with a harmonic, in both genotypes.
In [18]:
rows = []
band = (FREQ >= 3) & (FREQ <= 30)
for sid, r in spectra.items():
n = r["n"]["SWD"]
row = {"subject": sid, "genotype": by_id[sid]["genotype"], "swd_epochs": n}
if n >= 5:
lp = np.log10(r["SWD_abs"][band])
f = FREQ[band]
pk, prop = signal.find_peaks(lp, prominence=0.1)
peaks = f[pk]
# any peak at 5-10 Hz with another peak within 1 Hz of twice its frequency
pairs = [(p, q) for p in peaks if 5 <= p <= 10 for q in peaks if abs(q - 2 * p) <= 1]
row["peaks_hz"] = ", ".join(f"{p:.1f}" for p in peaks)
row["f0_hz"] = pairs[0][0] if pairs else np.nan
row["harmonic_hz"] = pairs[0][1] if pairs else np.nan
row["meets"] = bool(pairs)
stored = by_id[sid]["spectra"]
if "SWD" in stored:
g = stored.groupby("frequency_hz").SWD.mean()
gb = g[(g.index >= 3) & (g.index <= 30)]
spk, _ = signal.find_peaks(gb.to_numpy(), prominence=0.02)
row["stored_peaks_hz"] = ", ".join(f"{p:.1f}" for p in gb.index[spk][:4])
rows.append(row)
peak = pd.DataFrame(rows)
print(peak.to_string(index=False))
ok = peak[peak.swd_epochs >= 5]
print("\nAnimals with >= 5 usable SWD epochs:", len(ok), "| meeting the criterion:", int(ok.meets.sum()),
"| by genotype:", ok.groupby(ok.genotype.fillna("none")).meets.agg(["sum", "size"]).to_dict("index"))
fig, ax = plt.subplots(figsize=(6, 3.5))
for sid, r in spectra.items():
if r["n"]["SWD"] >= 5:
ax.semilogy(FREQ[band], r["SWD_abs"][band], lw=0.8, color="C3" if by_id[sid]["genotype"] == "WT" else "C0")
ax.set(xlabel="Hz", ylabel="SWD power (uV^2/Hz)", title="SWD epochs; WT red, Het and unlabelled blue")
plt.tight_layout(); plt.show()
vals = np.concatenate([r["SWD_abs"][band] for r in spectra.values() if r["n"]["SWD"] >= 5])
print("plotted range (uV^2/Hz):", np.round([vals.min(), vals.max()], 2))
subject genotype swd_epochs peaks_hz f0_hz harmonic_hz meets stored_peaks_hz
GRIN2B-129 Het 86 8.0, 16.0, 23.5 8.0128 16.0256 True 5.0, 8.0, 15.8, 23.8
GRIN2B-130 Het 253 7.5, 15.5, 23.0 7.5120 15.5248 True 7.8, 15.6, 23.0
GRIN2B-131 Het 185 4.0, 8.0, 16.0 8.0128 16.0256 True 8.2, 16.2, 24.2
GRIN2B-132 NaN 20 5.0, 7.5, 15.5, 24.5 7.5120 15.5248 True 7.4, 14.8
GRIN2B-137 Het 79 7.0, 14.5, 21.0 7.0112 14.5232 True 7.6, 14.4, 21.6
GRIN2B-138 Het 292 7.5, 15.5, 22.5 7.5120 15.5248 True 7.8, 15.4, 23.0
GRIN2B-139 Het 397 7.5, 14.5, 21.5 7.5120 14.5232 True 7.6, 14.6, 22.0
GRIN2B-227 Het 427 7.5, 15.0, 22.5 7.5120 15.0240 True 7.6, 14.8, 22.2
GRIN2B-228 Het 0 NaN NaN NaN NaN NaN
GRIN2B-229 Het 99 8.0, 10.5, 16.0, 20.5 8.0128 16.0256 True 7.8, 15.8, 23.0
GRIN2B-236 Het 284 8.0, 16.0, 24.0 8.0128 16.0256 True 4.2, 8.0, 16.0, 23.8
GRIN2B-237 Het 183 5.5, 8.0, 15.5, 23.5 8.0128 15.5248 True 7.6, 15.0, 22.8
GRIN2B-238 WT 0 NaN NaN NaN NaN NaN
GRIN2B-239 Het 202 7.5, 15.0, 22.5 7.5120 15.0240 True 7.8, 15.0, 22.8
GRIN2B-240 WT 18 3.5, 7.0, 14.0, 21.0, 28.0 7.0112 14.0224 True 7.6, 14.4, 21.0
GRIN2B-241 Het 97 7.5, 14.5, 21.5 7.5120 14.5232 True 7.8, 15.0, 22.2
GRIN2B-373 Het 275 8.0, 15.5, 23.0 8.0128 15.5248 True 4.4, 8.2, 16.0, 24.0
Animals with >= 5 usable SWD epochs: 15 | meeting the criterion: 15 | by genotype: {'Het': {'sum': 13, 'size': 13}, 'WT': {'sum': True, 'size': 1}, 'none': {'sum': True, 'size': 1}}
plotted range (uV^2/Hz): [ 18.44 7664.53]
Result: reproduced. All 15 rats with at least 5 usable SWD epochs in the raw windows, 13 Het, the WT GRIN2B-240 and GRIN2B-132, have a spectral peak at 7.0 to 8.0 Hz and a second peak at 14.0 to 16.0 Hz, twice its frequency; 12 of them also have a third at 21 to 24 Hz. The stored spectra show peaks near the same frequencies (7.4 to 8.2 and 14.4 to 16.2 Hz). GRIN2B-240 has 18 usable SWD epochs, the fewest.
fig3d-rem-hourly¶
Grin2b+/- rats have less REM sleep per hour than WT across the day, and the reduction is not specific to certain hours (linear mixed model, genotype F = 9.98, p = .0045; hour x genotype p = .53).
In [19]:
H_rem = hourly_table(lambda d: grin2b.hourly_minutes(grin2b.scored_states(d, "events"), 2))
print("Mean REM minutes per hour by genotype (main):")
print(H_rem[H_rem.genotype.notna()].groupby(["genotype", "hour"]).value.mean().unstack().round(2).to_string())
print(hourly_tests(H_rem).round(5).to_string(index=False))
print("Mixed model (MixedLM, random intercept), main animals:")
print(mixed_check(H_rem).round(5).to_string(index=False))
for mode in ("flag", "stored"):
H_v = hourly_table(lambda d: grin2b.hourly_minutes(grin2b.scored_states(d, mode), 2))
print(f"Variant, SWD epochs mode {mode}:")
print(hourly_tests(H_v).round(5).head(3).to_string(index=False))
Mean REM minutes per hour by genotype (main):
hour 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
genotype
Het 1.61 2.59 2.52 3.57 3.71 2.26 3.08 3.66 3.15 3.2 3.22 2.10 1.79 1.30 0.62 0.19 0.26 0.85 0.28 0.32 0.60 1.45 1.40 1.10
WT 4.40 2.90 2.50 5.56 2.77 3.60 4.62 4.92 3.81 4.0 5.65 2.88 3.35 2.69 2.56 5.40 1.33 1.79 0.46 0.23 1.38 0.27 1.75 0.92
variant effect F df df_error p
main (132 left out) genotype 4.68577 1 14 0.04818
main (132 left out) hour 10.75616 23 322 0.00000
main (132 left out) genotype x hour 1.42338 23 322 0.09645
132 as WT genotype 1.85831 1 15 0.19293
132 as WT hour 11.06965 23 345 0.00000
132 as WT genotype x hour 1.02544 23 345 0.43189
132 as Het genotype 5.18457 1 15 0.03787
132 as Het hour 11.38941 23 345 0.00000
132 as Het genotype x hour 1.48837 23 345 0.07088
without 228 genotype 4.25599 1 13 0.05968
without 228 hour 10.74550 23 299 0.00000
without 228 genotype x hour 1.46112 23 299 0.08196
without 238 genotype 7.16628 1 13 0.01901
without 238 hour 10.84440 23 299 0.00000
without 238 genotype x hour 2.64154 23 299 0.00010
without 228 and 238 genotype 6.53716 1 12 0.02515
without 228 and 238 hour 10.97554 23 276 0.00000
without 228 and 238 genotype x hour 2.75893 23 276 0.00005
Mixed model (MixedLM, random intercept), main animals:
term chi2_over_df df p converged
C(genotype, Sum) 4.68575 1 0.03041 True
C(hour, Sum) 5.58742 23 0.00000 True
C(genotype, Sum):C(hour, Sum) 1.42338 23 0.08575 True
Variant, SWD epochs mode flag:
variant effect F df df_error p
main (132 left out) genotype 4.63449 1 14 0.04926
main (132 left out) hour 10.75152 23 322 0.00000
main (132 left out) genotype x hour 1.41779 23 322 0.09890
Variant, SWD epochs mode stored:
variant effect F df df_error p
main (132 left out) genotype 4.61310 1 14 0.04972
main (132 left out) hour 10.75393 23 322 0.00000
main (132 left out) genotype x hour 1.41656 23 322 0.09945
Result: reproduced. Het have less REM per hour than WT (genotype F = 4.69, p = 0.048), and the reduction does not depend on the hour (interaction p = 0.096). The effect is at the threshold: with GRIN2B-132 as WT p = 0.19, as Het p = 0.038 (interaction p = 0.071), and with the scorer's flags or the stored scores p = 0.049 and 0.050. Without GRIN2B-228 genotype p = 0.060; without GRIN2B-238 or both genotype p = 0.019 and 0.025, but the interaction then has p = 0.0001 and 0.00005. The outcome holds with all 16 rats and with GRIN2B-132 counted as Het, and not without GRIN2B-228 or GRIN2B-238. The mixed model gives genotype 4.69 (Wald p = 0.030) and interaction p = 0.086.
fig3g-rem-total¶
Grin2b+/- rats spend less total time in REM sleep over 24 h than WT (two-sample t test, df = 22, T = 3.16, p = .0045).
In [20]:
print(A[["genotype", "rem_min"]].sort_values(["genotype", "rem_min"]).round(2).to_string())
print(compare("rem_min", "t", "lower",
[("scorer flags", animals["flag"]), ("stored scores", animals["stored"]), ("SWD time only", animals["exact"])]).round(5).to_string(index=False))
genotype rem_min
subject
GRIN2B-129 Het 15.08
GRIN2B-236 Het 29.12
GRIN2B-137 Het 30.00
GRIN2B-139 Het 37.50
GRIN2B-373 Het 37.88
GRIN2B-228 Het 42.54
GRIN2B-138 Het 42.83
GRIN2B-237 Het 45.29
GRIN2B-229 Het 49.42
GRIN2B-130 Het 54.38
GRIN2B-241 Het 54.88
GRIN2B-239 Het 55.17
GRIN2B-131 Het 65.54
GRIN2B-227 Het 68.21
GRIN2B-238 WT 54.25
GRIN2B-240 WT 85.21
GRIN2B-132 NaN 37.25
variant n_WT n_Het WT_mean Het_mean Het_minus_WT test stat df p paper_direction
main (132 left out) 2 14 69.72917 44.84524 -24.88393 t 2.16466 14 0.04818 True
132 as WT 3 14 58.90278 44.84524 -14.05754 t 1.36320 15 0.19293 True
132 as Het 2 15 69.72917 44.33889 -25.39028 t 2.27697 15 0.03787 True
without 228 2 13 69.72917 45.02244 -24.70673 t 2.06301 13 0.05968 True
without 238 1 14 85.20833 44.84524 -40.36310 t 2.67699 13 0.01901 True
without 228 and 238 1 13 85.20833 45.02244 -40.18590 t 2.55679 12 0.02515 True
scorer flags main (132 left out) 2 14 69.72917 44.94048 -24.78869 t 2.15279 14 0.04926 True
stored scores main (132 left out) 2 14 69.72917 44.96726 -24.76190 t 2.14781 14 0.04972 True
SWD time only main (132 left out) 2 14 69.72917 44.93185 -24.79732 t 2.15325 14 0.04922 True
Result: reproduced. Het have 44.8 min of REM per day against 69.7 in WT (t = 2.16, df = 14, p = 0.048). The result is at the threshold: with GRIN2B-132 as WT p = 0.19, as Het p = 0.038, and the two SWD-epoch variants give p = 0.049 and 0.050. Both WT (54.3 and 85.2 min) are above the Het median, and 5 Het rats have more REM than GRIN2B-238. Without GRIN2B-228 p = 0.060; without GRIN2B-238 p = 0.019, without both 0.025. Removing only the SWD time gives p = 0.049.
fig3g-rem-bouts¶
Grin2b+/- rats have fewer REM bouts than WT (two-sample t test, df = 22, T = 3.94, p = .00069).
In [21]:
print(A[["genotype", "rem_bouts"]].sort_values(["genotype", "rem_bouts"]).round(2).to_string())
print(compare("rem_bouts", "t", "lower",
[("scorer flags", animals["flag"]), ("stored scores", animals["stored"]), ("bridged bouts", animals["bridged"])]).round(5).to_string(index=False))
genotype rem_bouts
subject
GRIN2B-129 Het 14.0
GRIN2B-137 Het 22.5
GRIN2B-373 Het 30.5
GRIN2B-228 Het 31.0
GRIN2B-139 Het 36.0
GRIN2B-236 Het 36.5
GRIN2B-130 Het 39.0
GRIN2B-138 Het 41.0
GRIN2B-229 Het 41.0
GRIN2B-237 Het 44.0
GRIN2B-239 Het 44.5
GRIN2B-241 Het 44.5
GRIN2B-131 Het 50.0
GRIN2B-227 Het 55.5
GRIN2B-238 WT 60.5
GRIN2B-240 WT 87.5
GRIN2B-132 NaN 23.5
variant n_WT n_Het WT_mean Het_mean Het_minus_WT test stat df p paper_direction
main (132 left out) 2 14 74.00000 37.85714 -36.14286 t 4.12410 14 0.00103 True
132 as WT 3 14 57.16667 37.85714 -19.30952 t 1.96409 15 0.06833 True
132 as Het 2 15 74.00000 36.90000 -37.10000 t 4.19122 15 0.00079 True
without 228 2 13 74.00000 38.38462 -35.61538 t 3.95091 13 0.00166 True
without 238 1 14 87.50000 37.85714 -49.64286 t 4.43939 13 0.00067 True
without 228 and 238 1 13 87.50000 38.38462 -49.11538 t 4.28118 12 0.00107 True
scorer flags main (132 left out) 2 14 74.00000 37.21429 -36.78571 t 4.26936 14 0.00078 True
stored scores main (132 left out) 2 14 74.00000 37.03571 -36.96429 t 4.31222 14 0.00072 True
bridged bouts main (132 left out) 2 14 74.00000 37.03571 -36.96429 t 4.31222 14 0.00072 True
Result: reproduced. Het have 37.9 REM bouts per day against 74.0 in WT (t = 4.12, df = 14, p = 0.0010); both WT (60.5 and 87.5) have more than every Het (at most 55.5). With GRIN2B-132 as WT p = 0.068, as Het p = 0.0008; the SWD-epoch variants give p = 0.0008 and 0.0007. Without GRIN2B-228, GRIN2B-238 or both p = 0.0017, 0.0007 and 0.0011; with bouts not split by SWD epochs p = 0.0007.
fig3g-rem-bout-duration¶
REM bouts are longer on average in Grin2b+/- rats than in WT (Wilcoxon rank-sum, W = 33, p = .043).
In [22]:
print(A[["genotype", "rem_bout_s"]].sort_values(["genotype", "rem_bout_s"]).round(2).to_string())
print(compare("rem_bout_s", "wilcoxon", "higher",
[("scorer flags", animals["flag"]), ("stored scores", animals["stored"]), ("bridged bouts", animals["bridged"])]).round(5).to_string(index=False))
genotype rem_bout_s
subject
GRIN2B-236 Het 47.77
GRIN2B-237 Het 62.12
GRIN2B-139 Het 62.74
GRIN2B-138 Het 63.84
GRIN2B-129 Het 64.39
GRIN2B-229 Het 72.23
GRIN2B-241 Het 73.96
GRIN2B-227 Het 74.44
GRIN2B-239 Het 74.45
GRIN2B-373 Het 74.52
GRIN2B-137 Het 80.64
GRIN2B-131 Het 83.02
GRIN2B-130 Het 83.58
GRIN2B-228 Het 89.37
GRIN2B-238 WT 55.44
GRIN2B-240 WT 58.17
GRIN2B-132 NaN 95.52
variant n_WT n_Het WT_mean Het_mean Het_minus_WT test stat df p paper_direction
main (132 left out) 2 14 56.80815 71.93441 15.12626 Wilcoxon rank-sum (W) 2.0 NaN 0.06667 True
132 as WT 3 14 69.71237 71.93441 2.22203 Wilcoxon rank-sum (W) 16.0 NaN 0.59118 True
132 as Het 2 15 56.80815 73.50684 16.69869 Wilcoxon rank-sum (W) 2.0 NaN 0.05882 True
without 228 2 13 56.80815 70.59303 13.78488 Wilcoxon rank-sum (W) 2.0 NaN 0.07619 True
without 238 1 14 58.17270 71.93441 13.76171 Wilcoxon rank-sum (W) 1.0 NaN 0.26667 True
without 228 and 238 1 13 58.17270 70.59303 12.42033 Wilcoxon rank-sum (W) 1.0 NaN 0.28571 True
scorer flags main (132 left out) 2 14 56.80815 73.06874 16.26059 Wilcoxon rank-sum (W) 2.0 NaN 0.06667 True
stored scores main (132 left out) 2 14 56.80815 73.37888 16.57074 Wilcoxon rank-sum (W) 2.0 NaN 0.06667 True
bridged bouts main (132 left out) 2 14 56.80815 73.20604 16.39789 Wilcoxon rank-sum (W) 2.0 NaN 0.06667 True
Result: partially reproduced. Het REM bouts last 71.9 s on average against 56.8 s in WT, the paper's direction, with W = 2 and p = 0.067. With GRIN2B-132 (95.5 s) as WT the difference is 2.2 s (p = 0.59). Without GRIN2B-228 p = 0.076; with one WT p = 0.27 and 0.29. With bouts not split by SWD epochs p = 0.067.
fig3h-nrem-total¶
Total NREM sleep over 24 h does not differ between genotypes (two-sample t test, df = 22, T = -.09, p = .93).
In [23]:
print(A[["genotype", "nrem_min"]].sort_values(["genotype", "nrem_min"]).round(2).to_string())
print(compare("nrem_min", "t", "none",
[("scorer flags", animals["flag"]), ("stored scores", animals["stored"]), ("SWD time only", animals["exact"])]).round(5).to_string(index=False))
genotype nrem_min
subject
GRIN2B-239 Het 412.58
GRIN2B-139 Het 417.33
GRIN2B-241 Het 424.33
GRIN2B-228 Het 428.25
GRIN2B-137 Het 433.46
GRIN2B-236 Het 444.83
GRIN2B-130 Het 444.88
GRIN2B-229 Het 452.54
GRIN2B-227 Het 453.75
GRIN2B-138 Het 456.46
GRIN2B-373 Het 464.75
GRIN2B-129 Het 476.33
GRIN2B-131 Het 495.08
GRIN2B-237 Het 502.46
GRIN2B-238 WT 370.42
GRIN2B-240 WT 468.12
GRIN2B-132 NaN 367.75
variant n_WT n_Het WT_mean Het_mean Het_minus_WT test stat df p paper_direction
main (132 left out) 2 14 419.27083 450.50298 31.23214 t -1.28589 14 0.21934 False
132 as WT 3 14 402.09722 450.50298 48.40575 t -2.31356 15 0.03529 False
132 as Het 2 15 419.27083 444.98611 25.71528 t -0.91638 15 0.37397 False
without 228 2 13 419.27083 452.21474 32.94391 t -1.32547 13 0.20784 False
without 238 1 14 468.12500 450.50298 -17.62202 t 0.62390 13 0.54349 True
without 228 and 238 1 13 468.12500 452.21474 -15.91026 t 0.55533 12 0.58887 True
scorer flags main (132 left out) 2 14 419.31250 453.58333 34.27083 t -1.40345 14 0.18227 False
stored scores main (132 left out) 2 14 419.33333 454.87202 35.53869 t -1.45241 14 0.16843 False
SWD time only main (132 left out) 2 14 419.32333 453.18095 33.85762 t -1.38717 14 0.18708 False
Result: reproduced. The t test finds no difference in NREM time (Het 450.5 against WT 419.3 min per day, difference 31.2 min, t = -1.29, p = 0.22). With GRIN2B-132 as WT it does (48.4 min, p = 0.035); as Het p = 0.37. Without GRIN2B-228 p = 0.21; without GRIN2B-238 the Het have 17.6 min less (p = 0.54), and 15.9 min less without both (p = 0.59). Removing only the SWD time gives p = 0.19.
fig3i-wake-total¶
Grin2b+/- rats spend less total time awake over 24 h than WT (Welch t test, df = 10.25, T = 2.26, p = .047).
In [24]:
print(A[["genotype", "wake_min"]].sort_values(["genotype", "wake_min"]).round(2).to_string())
print(compare("wake_min", "welch", "lower",
[("scorer flags", animals["flag"]), ("stored scores", animals["stored"]), ("SWD time only", animals["exact"])]).round(5).to_string(index=False))
genotype wake_min
subject
GRIN2B-131 Het 753.29
GRIN2B-373 Het 754.00
GRIN2B-227 Het 765.75
GRIN2B-139 Het 797.75
GRIN2B-138 Het 808.17
GRIN2B-237 Het 812.54
GRIN2B-129 Het 828.12
GRIN2B-239 Het 839.04
GRIN2B-236 Het 841.00
GRIN2B-241 Het 861.92
GRIN2B-130 Het 862.54
GRIN2B-229 Het 883.08
GRIN2B-137 Het 900.88
GRIN2B-228 Het 969.12
GRIN2B-240 WT 874.42
GRIN2B-238 WT 1015.33
GRIN2B-132 NaN 1023.38
variant n_WT n_Het WT_mean Het_mean Het_minus_WT test stat df p paper_direction
main (132 left out) 2 14 944.87500 834.08631 -110.78869 Welch t 1.53315 1.10618 0.35056 True
132 as WT 3 14 971.04167 834.08631 -136.95536 Welch t 2.68750 2.45983 0.09233 True
132 as Het 2 15 944.87500 846.70556 -98.16944 Welch t 1.34255 1.15949 0.38408 True
without 228 2 13 944.87500 823.69872 -121.17628 Welch t 1.69038 1.07142 0.32761 True
without 238 1 14 874.41667 834.08631 -40.33036 Welch t NaN 1.00000 NaN True
without 228 and 238 1 13 874.41667 823.69872 -50.71795 Welch t NaN 1.00000 NaN True
scorer flags main (132 left out) 2 14 949.35417 881.27976 -68.07440 Welch t 1.01301 1.07600 0.48623 True
stored scores main (132 left out) 2 14 950.93750 940.16071 -10.77679 Welch t 0.16594 1.03436 0.89462 True
SWD time only main (132 left out) 2 14 949.57417 884.11589 -65.45827 Welch t 0.97793 1.07339 0.49804 True
Result: partially reproduced. Het are awake 834.1 min per day against 944.9 in WT (difference -110.8 min), the paper's direction, but the Welch test gives p = 0.35 with two WT. The difference depends on how SWDs are taken out of wake: it is -65.5 min (p = 0.50) when only the SWD time is removed, -68.1 min (p = 0.49) with the scorer's flags and -10.8 min (p = 0.89) with the stored scores. Without GRIN2B-228 it is -121.2 min (p = 0.33); with one WT the Welch test is not defined.
fig3h-nrem-bouts¶
Grin2b+/- rats have fewer NREM bouts, which are longer, than WT (Wilcoxon rank-sum, bouts W = 29.5, p = .025; mean bout duration W = 26, p = .015).
In [25]:
print(A[["genotype", "nrem_bouts", "nrem_bout_s"]].sort_values("genotype").round(2).to_string())
print("\nnrem_bouts:")
print(compare("nrem_bouts", "wilcoxon", "lower", [("scorer flags", animals["flag"]), ("stored scores", animals["stored"]), ("bridged bouts", animals["bridged"])]).round(5).to_string(index=False))
print("\nnrem_bout_s:")
print(compare("nrem_bout_s", "wilcoxon", "higher", [("scorer flags", animals["flag"]), ("stored scores", animals["stored"]), ("bridged bouts", animals["bridged"])]).round(5).to_string(index=False))
genotype nrem_bouts nrem_bout_s
subject
GRIN2B-129 Het 143.0 202.99
GRIN2B-130 Het 218.5 122.16
GRIN2B-131 Het 205.0 144.98
GRIN2B-137 Het 210.0 124.90
GRIN2B-138 Het 205.5 134.65
GRIN2B-139 Het 256.0 98.40
GRIN2B-227 Het 241.5 112.86
GRIN2B-228 Het 214.5 120.40
GRIN2B-229 Het 202.5 134.44
GRIN2B-236 Het 225.5 120.05
GRIN2B-237 Het 258.0 116.75
GRIN2B-239 Het 222.0 111.50
GRIN2B-241 Het 202.0 125.57
GRIN2B-373 Het 255.0 109.29
GRIN2B-238 WT 253.0 88.83
GRIN2B-240 WT 293.0 95.71
GRIN2B-132 NaN 171.5 129.56
nrem_bouts:
variant n_WT n_Het WT_mean Het_mean Het_minus_WT test stat df p paper_direction
main (132 left out) 2 14 273.00000 218.50000 -54.50000 Wilcoxon rank-sum (W) 25.0 NaN 0.10000 True
132 as WT 3 14 239.16667 218.50000 -20.66667 Wilcoxon rank-sum (W) 26.0 NaN 0.59118 True
132 as Het 2 15 273.00000 215.36667 -57.63333 Wilcoxon rank-sum (W) 27.0 NaN 0.08824 True
without 228 2 13 273.00000 218.80769 -54.19231 Wilcoxon rank-sum (W) 23.0 NaN 0.11429 True
without 238 1 14 293.00000 218.50000 -74.50000 Wilcoxon rank-sum (W) 14.0 NaN 0.13333 True
without 228 and 238 1 13 293.00000 218.80769 -74.19231 Wilcoxon rank-sum (W) 13.0 NaN 0.14286 True
scorer flags main (132 left out) 2 14 273.25000 211.03571 -62.21429 Wilcoxon rank-sum (W) 28.0 NaN 0.03182 True
stored scores main (132 left out) 2 14 273.00000 198.85714 -74.14286 Wilcoxon rank-sum (W) 28.0 NaN 0.03182 True
bridged bouts main (132 left out) 2 14 273.00000 194.42857 -78.57143 Wilcoxon rank-sum (W) 28.0 NaN 0.03195 True
nrem_bout_s:
variant n_WT n_Het WT_mean Het_mean Het_minus_WT test stat df p paper_direction
main (132 left out) 2 14 92.26769 127.06881 34.80113 Wilcoxon rank-sum (W) 0.0 NaN 0.01667 True
132 as WT 3 14 104.69968 127.06881 22.36914 Wilcoxon rank-sum (W) 10.0 NaN 0.19706 True
132 as Het 2 15 92.26769 127.23514 34.96745 Wilcoxon rank-sum (W) 0.0 NaN 0.01471 True
without 228 2 13 92.26769 127.58190 35.31421 Wilcoxon rank-sum (W) 0.0 NaN 0.01905 True
without 238 1 14 95.70738 127.06881 31.36143 Wilcoxon rank-sum (W) 0.0 NaN 0.13333 True
without 228 and 238 1 13 95.70738 127.58190 31.87451 Wilcoxon rank-sum (W) 0.0 NaN 0.14286 True
scorer flags main (132 left out) 2 14 92.19622 132.11738 39.92116 Wilcoxon rank-sum (W) 0.0 NaN 0.01667 True
stored scores main (132 left out) 2 14 92.28094 140.72102 48.44008 Wilcoxon rank-sum (W) 0.0 NaN 0.01667 True
bridged bouts main (132 left out) 2 14 92.26769 142.41370 50.14601 Wilcoxon rank-sum (W) 0.0 NaN 0.01667 True
Result: partially reproduced. Het NREM bouts are longer than in WT (127.1 against 92.3 s, W = 0, p = 0.017; both WT below every Het). Het have fewer NREM bouts (218.5 against 273.0 per day), the paper's direction, but p = 0.10. With the scorer's flags, the stored scores or bouts not split by SWD epochs both parts have p < 0.05 (bouts p = 0.032, duration p = 0.017). Without GRIN2B-228 bouts p = 0.11 and duration p = 0.019; with one WT neither test can reach p < 0.05.
fig3i-wake-bouts¶
Grin2b+/- rats have more wake bouts of shorter duration than WT (Wilcoxon rank-sum, W = 18, p = .0035; Welch t test on mean bout duration, df = 8.91, T = 2.97, p = .016).
In [26]:
print(A[["genotype", "wake_bouts", "wake_bout_s"]].sort_values("genotype").round(2).to_string())
print("\nwake_bouts:")
print(compare("wake_bouts", "wilcoxon", "higher", [("scorer flags", animals["flag"]), ("stored scores", animals["stored"]), ("bridged bouts", animals["bridged"])]).round(5).to_string(index=False))
print("\nwake_bout_s:")
print(compare("wake_bout_s", "welch", "lower", [("scorer flags", animals["flag"]), ("stored scores", animals["stored"]), ("bridged bouts", animals["bridged"])]).round(5).to_string(index=False))
genotype wake_bouts wake_bout_s
subject
GRIN2B-129 Het 577.5 86.70
GRIN2B-130 Het 533.0 97.19
GRIN2B-131 Het 639.0 70.71
GRIN2B-137 Het 628.0 88.94
GRIN2B-138 Het 692.5 70.18
GRIN2B-139 Het 747.0 64.53
GRIN2B-227 Het 656.0 70.11
GRIN2B-228 Het 206.0 282.53
GRIN2B-229 Het 515.0 104.44
GRIN2B-236 Het 507.5 100.33
GRIN2B-237 Het 458.5 106.38
GRIN2B-239 Het 604.5 83.73
GRIN2B-241 Het 580.0 89.13
GRIN2B-373 Het 487.0 92.98
GRIN2B-238 WT 230.5 265.89
GRIN2B-240 WT 358.5 146.47
GRIN2B-132 NaN 272.5 226.74
wake_bouts:
variant n_WT n_Het WT_mean Het_mean Het_minus_WT test stat df p paper_direction
main (132 left out) 2 14 294.50000 559.39286 264.89286 Wilcoxon rank-sum (W) 2.0 NaN 0.06667 True
132 as WT 3 14 287.16667 559.39286 272.22619 Wilcoxon rank-sum (W) 3.0 NaN 0.02059 True
132 as Het 2 15 294.50000 540.26667 245.76667 Wilcoxon rank-sum (W) 3.0 NaN 0.08824 True
without 228 2 13 294.50000 586.57692 292.07692 Wilcoxon rank-sum (W) 0.0 NaN 0.01905 True
without 238 1 14 358.50000 559.39286 200.89286 Wilcoxon rank-sum (W) 1.0 NaN 0.26667 True
without 228 and 238 1 13 358.50000 586.57692 228.07692 Wilcoxon rank-sum (W) 0.0 NaN 0.14286 True
scorer flags main (132 left out) 2 14 259.50000 578.78571 319.28571 Wilcoxon rank-sum (W) 2.0 NaN 0.06667 True
stored scores main (132 left out) 2 14 241.00000 190.78571 -50.21429 Wilcoxon rank-sum (W) 28.0 NaN 0.01667 False
bridged bouts main (132 left out) 2 14 241.00000 186.32143 -54.67857 Wilcoxon rank-sum (W) 28.0 NaN 0.01667 False
wake_bout_s:
variant n_WT n_Het WT_mean Het_mean Het_minus_WT test stat df p paper_direction
main (132 left out) 2 14 206.18100 100.56307 -105.61794 Welch t 1.71928 1.12031 0.31490 True
132 as WT 3 14 213.03406 100.56307 -112.47100 Welch t 2.95961 2.72150 0.06722 True
132 as Het 2 15 206.18100 108.97487 -97.20613 Welch t 1.57343 1.14582 0.33688 True
without 228 2 13 206.18100 86.56544 -119.61556 Welch t 1.99913 1.00847 0.29365 True
without 238 1 14 146.47344 100.56307 -45.91037 Welch t NaN 1.00000 NaN True
without 228 and 238 1 13 146.47344 86.56544 -59.90800 Welch t NaN 1.00000 NaN True
scorer flags main (132 left out) 2 14 225.26335 104.94756 -120.31579 Welch t 2.77979 1.28663 0.17278 True
stored scores main (132 left out) 2 14 239.16919 304.22780 65.05861 Welch t -2.12591 1.70805 0.18874 False
bridged bouts main (132 left out) 2 14 237.71789 275.40886 37.69097 Welch t -1.21259 1.47684 0.38426 False
Result: partially reproduced. Het have more wake bouts (559.4 against 294.5 per day, W = 2, p = 0.067) and shorter ones (100.6 against 206.2 s, Welch p = 0.31), both in the paper's direction without p < 0.05. With the stored scores, in which SWD epochs count as wake, both reverse (Het 190.8 against 241.0 bouts, p = 0.017; 304.2 against 239.2 s, p = 0.19), and with bouts not split by SWD epochs they reverse too (Het 186.3 against 241.0 bouts, p = 0.017; 275.4 against 237.7 s, p = 0.38), so the result depends on SWD epochs ending wake bouts. Without GRIN2B-228 the bout count has p = 0.019 (W = 0) and the bout duration p = 0.29.
fig4a-rem-spectrum¶
EEG spectral power during REM sleep does not differ between genotypes (two-way ANOVA genotype x band, genotype p = .14, interaction p = .69).
In [27]:
T = band_table("REM")
print(T[T.genotype.notna()].pivot_table(index=["genotype", "subject"], columns="band", values="value")[list(grin2b.BANDS)].round(3).to_string())
print("Het minus WT, mean over bands:", round(T[T.genotype == "Het"].value.mean() - T[T.genotype == "WT"].value.mean(), 4))
print(spectral_claim("REM").round(4).to_string(index=False))
band delta theta sigma beta gamma
genotype subject
Het GRIN2B-129 0.034 0.223 0.157 0.610 0.763
GRIN2B-130 0.084 0.501 0.209 0.267 0.421
GRIN2B-131 0.244 0.350 0.403 0.575 0.799
GRIN2B-137 0.099 0.657 0.543 1.066 1.311
GRIN2B-138 0.118 0.417 0.443 0.783 1.135
GRIN2B-139 0.060 0.173 0.207 0.397 0.488
GRIN2B-227 0.195 0.605 0.530 0.916 1.211
GRIN2B-228 0.178 0.471 0.487 0.432 0.637
GRIN2B-229 0.092 0.451 0.332 0.443 0.742
GRIN2B-236 0.174 0.633 0.668 0.751 1.076
GRIN2B-237 0.164 0.454 0.367 0.664 1.068
GRIN2B-239 0.127 0.488 0.410 0.623 0.890
GRIN2B-241 0.136 0.391 0.414 0.636 0.855
GRIN2B-373 0.204 0.465 0.576 2.453 2.112
WT GRIN2B-238 0.335 0.615 0.564 0.876 1.243
GRIN2B-240 0.360 0.649 0.625 0.895 1.157
Het minus WT, mean over bands: -0.1883
spectra variant n_WT n_Het genotype_F genotype_p interaction_F interaction_p mean_Het-WT delta_Het-WT delta_tukey_p theta_Het-WT theta_tukey_p sigma_Het-WT sigma_tukey_p beta_Het-WT beta_tukey_p gamma_Het-WT gamma_tukey_p
Methods, raw main (132 left out) 2 14 3.3043 0.0734 0.0302 0.9982 -0.1883 -0.2111 0.9956 -0.1839 0.9985 -0.1838 0.9985 -0.1273 0.9999 -0.2352 0.9903
Methods, raw 132 as WT 3 14 3.0469 0.0850 0.0980 0.9828 -0.1477 -0.1549 0.9981 -0.1681 0.9964 -0.1175 0.9998 -0.0695 1.0000 -0.2286 0.9688
Methods, raw 132 as Het 2 15 3.3702 0.0703 0.0245 0.9988 -0.1838 -0.2083 0.9949 -0.1748 0.9987 -0.1848 0.9979 -0.1304 0.9999 -0.2208 0.9922
Methods, raw without 228 2 13 2.8943 0.0937 0.0367 0.9974 -0.1803 -0.2143 0.9958 -0.1856 0.9986 -0.1897 0.9984 -0.1022 1.0000 -0.2100 0.9964
Methods, raw without 238 1 14 1.7311 0.1929 0.0106 0.9998 -0.1935 -0.2233 0.9996 -0.2010 0.9998 -0.2143 0.9997 -0.1369 1.0000 -0.1921 0.9999
Methods, raw without 228 and 238 1 13 1.5178 0.2228 0.0197 0.9992 -0.1856 -0.2265 0.9996 -0.2028 0.9998 -0.2201 0.9997 -0.1118 1.0000 -0.1669 1.0000
Methods, raw, contiguous bands main (132 left out) 2 14 4.6936 0.0337 0.1764 0.9498 -0.1994 -0.2169 0.9875 -0.1839 0.9962 -0.1304 0.9997 -0.1273 0.9998 -0.3387 0.8210
Methods, raw, contiguous bands 132 as WT 3 14 4.0778 0.0470 0.2503 0.9086 -0.1523 -0.1473 0.9968 -0.1681 0.9916 -0.0914 0.9999 -0.0695 1.0000 -0.2850 0.7971
Methods, raw, contiguous bands 132 as Het 2 15 4.8310 0.0310 0.1687 0.9537 -0.1956 -0.2164 0.9845 -0.1748 0.9967 -0.1295 0.9997 -0.1304 0.9997 -0.3268 0.8227
Methods, raw, contiguous bands without 228 2 13 4.0724 0.0477 0.1734 0.9513 -0.1908 -0.2165 0.9896 -0.1856 0.9966 -0.1244 0.9999 -0.1022 1.0000 -0.3252 0.8715
Methods, raw, contiguous bands without 238 1 14 1.9949 0.1626 0.0151 0.9995 -0.1839 -0.2286 0.9986 -0.2010 0.9995 -0.1607 0.9999 -0.1369 1.0000 -0.1925 0.9996
Methods, raw, contiguous bands without 228 and 238 1 13 1.7156 0.1952 0.0224 0.9990 -0.1753 -0.2283 0.9989 -0.2028 0.9996 -0.1547 1.0000 -0.1118 1.0000 -0.1790 0.9998
Methods, raw, SWD epochs from scorer flags main (132 left out) 2 14 21.7051 0.0000 0.6124 0.6551 -0.2968 -0.2094 0.8994 -0.2291 0.8397 -0.2206 0.8674 -0.3630 0.2621 -0.4621 0.0534
Methods, raw, SWD epochs from scorer flags 132 as WT 3 14 23.3613 0.0000 1.1515 0.3391 -0.2540 -0.1530 0.9502 -0.2106 0.7377 -0.1532 0.9498 -0.3014 0.2528 -0.4518 0.0089
Methods, raw, SWD epochs from scorer flags 132 as Het 2 15 19.8510 0.0000 0.5032 0.7335 -0.2856 -0.2067 0.9097 -0.2175 0.8809 -0.2194 0.8754 -0.3512 0.3125 -0.4333 0.0924
Methods, raw, SWD epochs from scorer flags without 228 2 13 20.3895 0.0000 0.5059 0.7315 -0.2976 -0.2124 0.9092 -0.2347 0.8469 -0.2297 0.8626 -0.3566 0.3315 -0.4545 0.0817
Methods, raw, SWD epochs from scorer flags without 238 1 14 10.8786 0.0016 0.1835 0.9462 -0.2982 -0.2218 0.9834 -0.2417 0.9706 -0.2467 0.9664 -0.3653 0.7285 -0.4157 0.5641
Methods, raw, SWD epochs from scorer flags without 228 and 238 1 13 10.2043 0.0022 0.1459 0.9641 -0.2990 -0.2248 0.9854 -0.2472 0.9725 -0.2558 0.9658 -0.3589 0.7827 -0.4082 0.6356
Methods, raw, slope rejection main (132 left out) 2 14 3.3043 0.0734 0.0302 0.9982 -0.1883 -0.2111 0.9956 -0.1839 0.9985 -0.1838 0.9985 -0.1273 0.9999 -0.2352 0.9903
Methods, raw, slope rejection 132 as WT 3 14 3.0469 0.0850 0.0980 0.9828 -0.1477 -0.1549 0.9981 -0.1681 0.9964 -0.1175 0.9998 -0.0695 1.0000 -0.2286 0.9688
Methods, raw, slope rejection 132 as Het 2 15 3.3702 0.0703 0.0245 0.9988 -0.1838 -0.2083 0.9949 -0.1748 0.9987 -0.1848 0.9979 -0.1304 0.9999 -0.2208 0.9922
Methods, raw, slope rejection without 228 2 13 2.8943 0.0937 0.0367 0.9974 -0.1803 -0.2143 0.9958 -0.1856 0.9986 -0.1897 0.9984 -0.1022 1.0000 -0.2100 0.9964
Methods, raw, slope rejection without 238 1 14 1.7311 0.1929 0.0106 0.9998 -0.1935 -0.2233 0.9996 -0.2010 0.9998 -0.2143 0.9997 -0.1369 1.0000 -0.1921 0.9999
Methods, raw, slope rejection without 228 and 238 1 13 1.5178 0.2228 0.0197 0.9992 -0.1856 -0.2265 0.9996 -0.2028 0.9998 -0.2201 0.9997 -0.1118 1.0000 -0.1669 1.0000
Methods, raw, cleaner S1 channel main (132 left out) 2 14 3.4558 0.0672 0.0192 0.9993 -0.1765 -0.2109 0.9917 -0.1754 0.9979 -0.1885 0.9964 -0.1303 0.9998 -0.1774 0.9977
Methods, raw, cleaner S1 channel 132 as WT 3 14 2.7243 0.1030 0.0610 0.9930 -0.1284 -0.1515 0.9969 -0.1575 0.9959 -0.1144 0.9997 -0.0590 1.0000 -0.1599 0.9954
Methods, raw, cleaner S1 channel 132 as Het 2 15 3.6190 0.0610 0.0179 0.9994 -0.1743 -0.2088 0.9902 -0.1673 0.9981 -0.1907 0.9949 -0.1359 0.9996 -0.1690 0.9980
Methods, raw, cleaner S1 channel without 228 2 13 2.9915 0.0884 0.0411 0.9967 -0.1659 -0.2133 0.9916 -0.1760 0.9980 -0.1929 0.9960 -0.1023 1.0000 -0.1452 0.9996
Methods, raw, cleaner S1 channel without 238 1 14 1.8324 0.1805 0.0060 0.9999 -0.1825 -0.1917 0.9997 -0.2151 0.9993 -0.1858 0.9998 -0.1604 0.9999 -0.1595 0.9999
Methods, raw, cleaner S1 channel without 228 and 238 1 13 1.5906 0.2121 0.0169 0.9994 -0.1720 -0.1941 0.9997 -0.2157 0.9994 -0.1902 0.9998 -0.1324 1.0000 -0.1273 1.0000
stored (scorer) main (132 left out) 2 14 1.9050 0.1719 0.1246 0.9731 -0.0746 -0.0459 1.0000 -0.1046 0.9970 -0.1300 0.9855 -0.0670 0.9999 -0.0254 1.0000
stored (scorer) 132 as WT 3 14 15.3576 0.0002 0.5272 0.7160 -0.1979 -0.2659 0.3680 -0.2813 0.2905 -0.2135 0.6752 -0.1388 0.9653 -0.0903 0.9984
stored (scorer) 132 as Het 2 15 0.4631 0.4983 0.0823 0.9876 -0.0449 0.0012 1.0000 -0.0622 1.0000 -0.1047 0.9994 -0.0482 1.0000 -0.0107 1.0000
stored (scorer) without 228 2 13 1.3812 0.2442 0.1223 0.9740 -0.0644 -0.0347 1.0000 -0.0934 0.9989 -0.1212 0.9920 -0.0570 1.0000 -0.0159 1.0000
stored (scorer) without 238 1 14 1.0766 0.3033 0.0179 0.9994 0.0740 0.0643 1.0000 0.1076 0.9996 0.0625 1.0000 0.0817 1.0000 0.0539 1.0000
stored (scorer) without 228 and 238 1 13 1.3609 0.2480 0.0185 0.9993 0.0841 0.0755 1.0000 0.1188 0.9991 0.0713 1.0000 0.0917 0.9999 0.0633 1.0000
Result: reproduced. For normalized REM band power neither genotype (F = 3.30, p = 0.073) nor the interaction (p = 0.998) has p < 0.05; Het are lower by 0.19 on average across bands. The genotype effect reaches p = 0.034 with contiguous bands and p < 0.0001 with SWD epochs from the scorer's flags, and has p = 0.17 with the stored spectra. The slope rule removes no epoch and leaves the result unchanged; with the cleaner S1 channel genotype p = 0.067 and interaction p = 0.9993. Without GRIN2B-228 p = 0.094, without GRIN2B-238 or both p = 0.19 and 0.22.
fig4b-nrem-beta¶
During NREM sleep Grin2b+/- rats have higher beta power than WT (two-way ANOVA genotype x band, interaction p = .038; Tukey beta p = .0008, other bands p >= .82).
In [28]:
T = band_table("NREM")
print(T[T.genotype.notna()].pivot_table(index=["genotype", "subject"], columns="band", values="value")[list(grin2b.BANDS)].round(3).to_string())
print("Het minus WT, mean over bands:", round(T[T.genotype == "Het"].value.mean() - T[T.genotype == "WT"].value.mean(), 4))
print(spectral_claim("NREM").round(4).to_string(index=False))
band delta theta sigma beta gamma
genotype subject
Het GRIN2B-129 1.512 1.315 1.241 0.967 0.742
GRIN2B-130 1.998 1.619 1.489 0.781 0.423
GRIN2B-131 0.743 0.645 0.594 0.422 0.409
GRIN2B-137 1.768 1.460 1.344 1.038 0.731
GRIN2B-138 1.987 1.681 1.538 1.126 0.724
GRIN2B-139 0.915 0.708 0.615 0.588 0.564
GRIN2B-227 2.260 1.659 1.464 1.057 0.716
GRIN2B-228 1.689 1.388 1.199 0.581 0.450
GRIN2B-229 1.478 1.227 1.182 0.841 0.617
GRIN2B-236 2.039 1.467 1.682 1.209 0.837
GRIN2B-237 1.580 1.284 1.250 0.959 0.749
GRIN2B-239 1.821 1.439 1.401 1.040 0.715
GRIN2B-241 1.765 1.255 1.164 0.872 0.654
GRIN2B-373 2.043 1.540 1.609 1.235 0.858
WT GRIN2B-238 1.490 1.301 1.306 1.042 0.808
GRIN2B-240 1.577 1.217 1.165 0.930 0.762
Het minus WT, mean over bands: 0.011
spectra variant n_WT n_Het genotype_F genotype_p interaction_F interaction_p mean_Het-WT delta_Het-WT delta_tukey_p theta_Het-WT theta_tukey_p sigma_Het-WT sigma_tukey_p beta_Het-WT beta_tukey_p gamma_Het-WT gamma_tukey_p
Methods, raw main (132 left out) 2 14 0.0122 0.9123 0.2605 0.9023 0.0110 0.1518 0.9996 0.0758 1.0000 0.0340 1.0000 -0.0780 1.0000 -0.1284 0.9999
Methods, raw 132 as WT 3 14 0.0941 0.7599 0.3075 0.8721 0.0249 0.1331 0.9992 0.1094 0.9998 0.0354 1.0000 -0.0537 1.0000 -0.0995 0.9999
Methods, raw 132 as Het 2 15 0.0061 0.9380 0.2573 0.9043 0.0075 0.1454 0.9996 0.0640 1.0000 0.0315 1.0000 -0.0777 1.0000 -0.1256 0.9999
Methods, raw without 228 2 13 0.0359 0.8503 0.2052 0.9347 0.0195 0.1515 0.9996 0.0717 1.0000 0.0394 1.0000 -0.0529 1.0000 -0.1125 1.0000
Methods, raw without 238 1 14 0.0830 0.7741 0.0995 0.9822 0.0408 0.1083 1.0000 0.1179 1.0000 0.1048 1.0000 -0.0219 1.0000 -0.1052 1.0000
Methods, raw without 228 and 238 1 13 0.1142 0.7366 0.0770 0.9890 0.0492 0.1081 1.0000 0.1138 1.0000 0.1102 1.0000 0.0032 1.0000 -0.0893 1.0000
Methods, raw, contiguous bands main (132 left out) 2 14 0.0504 0.8229 0.3941 0.8122 -0.0212 0.1149 0.9999 0.0758 1.0000 -0.0005 1.0000 -0.0780 1.0000 -0.2182 0.9891
Methods, raw, contiguous bands 132 as WT 3 14 0.0000 0.9999 0.4385 0.7804 0.0000 0.1057 0.9998 0.1094 0.9997 -0.0000 1.0000 -0.0537 1.0000 -0.1613 0.9947
Methods, raw, contiguous bands 132 as Het 2 15 0.0695 0.7927 0.3927 0.8133 -0.0241 0.1090 0.9999 0.0640 1.0000 -0.0005 1.0000 -0.0777 1.0000 -0.2151 0.9875
Methods, raw, contiguous bands without 228 2 13 0.0125 0.9113 0.3485 0.8442 -0.0108 0.1187 0.9999 0.0717 1.0000 0.0170 1.0000 -0.0529 1.0000 -0.2087 0.9933
Methods, raw, contiguous bands without 238 1 14 0.0136 0.9076 0.1413 0.9662 0.0156 0.0863 1.0000 0.1179 1.0000 0.0591 1.0000 -0.0219 1.0000 -0.1632 0.9999
Methods, raw, contiguous bands without 228 and 238 1 13 0.0359 0.8503 0.1252 0.9728 0.0260 0.0902 1.0000 0.1138 1.0000 0.0766 1.0000 0.0032 1.0000 -0.1537 1.0000
Methods, raw, SWD epochs from scorer flags main (132 left out) 2 14 1.0048 0.3196 0.9405 0.4458 -0.0875 0.1692 0.9970 0.0030 1.0000 -0.0683 1.0000 -0.3000 0.8724 -0.2415 0.9635
Methods, raw, SWD epochs from scorer flags 132 as WT 3 14 0.9645 0.3292 1.1972 0.3191 -0.0699 0.1517 0.9939 0.0409 1.0000 -0.0616 1.0000 -0.2699 0.7940 -0.2105 0.9451
Methods, raw, SWD epochs from scorer flags 132 as Het 2 15 1.0099 0.3182 0.9006 0.4681 -0.0852 0.1615 0.9974 -0.0048 1.0000 -0.0651 1.0000 -0.2860 0.8845 -0.2316 0.9665
Methods, raw, SWD epochs from scorer flags without 228 2 13 0.9213 0.3407 0.8328 0.5092 -0.0870 0.1703 0.9976 -0.0071 1.0000 -0.0713 1.0000 -0.2925 0.9086 -0.2345 0.9762
Methods, raw, SWD epochs from scorer flags without 238 1 14 0.1869 0.6670 0.3479 0.8445 -0.0535 0.1251 1.0000 0.0499 1.0000 0.0099 1.0000 -0.2360 0.9973 -0.2162 0.9986
Methods, raw, SWD epochs from scorer flags without 228 and 238 1 13 0.1698 0.6817 0.3015 0.8758 -0.0529 0.1262 1.0000 0.0399 1.0000 0.0069 1.0000 -0.2285 0.9984 -0.2092 0.9992
Methods, raw, slope rejection main (132 left out) 2 14 0.0122 0.9123 0.2605 0.9023 0.0110 0.1518 0.9996 0.0758 1.0000 0.0340 1.0000 -0.0780 1.0000 -0.1284 0.9999
Methods, raw, slope rejection 132 as WT 3 14 0.0941 0.7599 0.3075 0.8721 0.0249 0.1331 0.9992 0.1094 0.9998 0.0354 1.0000 -0.0537 1.0000 -0.0995 0.9999
Methods, raw, slope rejection 132 as Het 2 15 0.0061 0.9380 0.2573 0.9043 0.0075 0.1454 0.9996 0.0640 1.0000 0.0315 1.0000 -0.0777 1.0000 -0.1256 0.9999
Methods, raw, slope rejection without 228 2 13 0.0359 0.8503 0.2052 0.9347 0.0195 0.1515 0.9996 0.0717 1.0000 0.0394 1.0000 -0.0529 1.0000 -0.1125 1.0000
Methods, raw, slope rejection without 238 1 14 0.0830 0.7741 0.0995 0.9822 0.0408 0.1083 1.0000 0.1179 1.0000 0.1048 1.0000 -0.0219 1.0000 -0.1052 1.0000
Methods, raw, slope rejection without 228 and 238 1 13 0.1142 0.7366 0.0770 0.9890 0.0492 0.1081 1.0000 0.1138 1.0000 0.1102 1.0000 0.0032 1.0000 -0.0893 1.0000
Methods, raw, cleaner S1 channel main (132 left out) 2 14 0.3617 0.5495 0.4268 0.7888 0.0490 0.2203 0.9683 0.0867 1.0000 0.0601 1.0000 -0.0348 1.0000 -0.0875 1.0000
Methods, raw, cleaner S1 channel 132 as WT 3 14 0.9937 0.3221 0.5114 0.7275 0.0661 0.2064 0.9263 0.1276 0.9972 0.0663 1.0000 -0.0102 1.0000 -0.0593 1.0000
Methods, raw, cleaner S1 channel 132 as Het 2 15 0.2856 0.5946 0.4147 0.7975 0.0423 0.2085 0.9734 0.0727 1.0000 0.0549 1.0000 -0.0374 1.0000 -0.0873 1.0000
Methods, raw, cleaner S1 channel without 228 2 13 0.5488 0.4615 0.3633 0.8338 0.0610 0.2264 0.9648 0.0851 1.0000 0.0687 1.0000 -0.0063 1.0000 -0.0687 1.0000
Methods, raw, cleaner S1 channel without 238 1 14 0.2687 0.6060 0.1096 0.9787 0.0598 0.1501 0.9999 0.1023 1.0000 0.1001 1.0000 0.0097 1.0000 -0.0631 1.0000
Methods, raw, cleaner S1 channel without 228 and 238 1 13 0.3784 0.5408 0.0875 0.9860 0.0719 0.1562 0.9998 0.1008 1.0000 0.1086 1.0000 0.0383 1.0000 -0.0443 1.0000
stored (scorer) main (132 left out) 2 14 0.0680 0.7951 0.2219 0.9253 -0.0139 0.0804 0.9996 -0.0266 1.0000 -0.0712 0.9998 -0.0203 1.0000 -0.0320 1.0000
stored (scorer) 132 as WT 3 14 8.1591 0.0055 0.2410 0.9143 -0.1476 -0.1820 0.8557 -0.2082 0.7318 -0.1726 0.8906 -0.0975 0.9976 -0.0778 0.9996
stored (scorer) 132 as Het 2 15 0.0454 0.8318 0.2175 0.9279 0.0137 0.1275 0.9965 0.0115 1.0000 -0.0462 1.0000 -0.0035 1.0000 -0.0208 1.0000
stored (scorer) without 228 2 13 0.0006 0.9811 0.2289 0.9212 -0.0013 0.0951 0.9984 -0.0122 1.0000 -0.0560 1.0000 -0.0084 1.0000 -0.0247 1.0000
stored (scorer) without 238 1 14 4.6239 0.0353 0.1976 0.9387 0.1491 0.2205 0.9157 0.2025 0.9486 0.1419 0.9954 0.1353 0.9968 0.0452 1.0000
stored (scorer) without 228 and 238 1 13 5.5913 0.0213 0.2202 0.9261 0.1618 0.2352 0.8714 0.2169 0.9168 0.1572 0.9893 0.1471 0.9933 0.0525 1.0000
Result: not reproduced. NREM beta power is lower in Het than in WT (difference -0.078 in normalized power, Tukey p = 1.0), and the genotype by band interaction has p = 0.90. In every variant with both WT, Het beta is lower, by 0.0035 to 0.30. With GRIN2B-238 left out, Het beta is higher in 7 of the 12 variants, by 0.0032 to 0.147, with Tukey p >= 0.99 and interaction p >= 0.92.
fig4c-wake-power¶
During wake, spectral power is reduced in Grin2b+/- rats overall, without a difference in single bands (two-way ANOVA genotype x band, genotype p = .01, interaction p = .86).
In [29]:
T = band_table("Wake")
print(T[T.genotype.notna()].pivot_table(index=["genotype", "subject"], columns="band", values="value")[list(grin2b.BANDS)].round(3).to_string())
print("Het minus WT, mean over bands:", round(T[T.genotype == "Het"].value.mean() - T[T.genotype == "WT"].value.mean(), 4))
print(spectral_claim("Wake").round(4).to_string(index=False))
band delta theta sigma beta gamma
genotype subject
Het GRIN2B-129 0.348 0.600 0.695 1.044 1.332
GRIN2B-130 0.650 0.778 0.851 1.138 1.274
GRIN2B-131 1.579 1.719 1.790 2.037 2.011
GRIN2B-137 0.271 0.548 0.677 0.955 1.255
GRIN2B-138 0.396 0.578 0.684 0.942 1.190
GRIN2B-139 1.221 1.454 1.560 1.568 1.585
GRIN2B-227 0.437 0.703 0.808 0.979 1.124
GRIN2B-228 0.619 0.792 0.916 1.332 1.400
GRIN2B-229 0.545 0.792 0.850 1.214 1.436
GRIN2B-236 0.451 0.753 0.635 0.891 1.087
GRIN2B-237 0.388 0.722 0.780 1.106 1.312
GRIN2B-239 0.368 0.672 0.728 1.044 1.311
GRIN2B-241 0.487 0.884 0.958 1.172 1.323
GRIN2B-373 0.413 0.705 0.654 0.751 1.008
WT GRIN2B-238 0.551 0.718 0.726 0.979 1.180
GRIN2B-240 0.633 0.873 0.914 1.062 1.162
Het minus WT, mean over bands: 0.0813
spectra variant n_WT n_Het genotype_F genotype_p interaction_F interaction_p mean_Het-WT delta_Het-WT delta_tukey_p theta_Het-WT theta_tukey_p sigma_Het-WT sigma_tukey_p beta_Het-WT beta_tukey_p gamma_Het-WT gamma_tukey_p
Methods, raw main (132 left out) 2 14 0.5835 0.4475 0.0834 0.9873 0.0813 -0.0083 1.0000 0.0402 1.0000 0.0793 1.0000 0.1348 0.9999 0.1608 0.9996
Methods, raw 132 as WT 3 14 0.5384 0.4654 0.1263 0.9725 0.0636 -0.0130 1.0000 0.0037 1.0000 0.0657 1.0000 0.1150 0.9999 0.1464 0.9990
Methods, raw 132 as Het 2 15 0.5997 0.4411 0.0788 0.9886 0.0795 -0.0068 1.0000 0.0448 1.0000 0.0767 1.0000 0.1297 0.9999 0.1530 0.9996
Methods, raw without 228 2 13 0.4896 0.4866 0.0694 0.9910 0.0775 -0.0110 1.0000 0.0435 1.0000 0.0780 1.0000 0.1212 1.0000 0.1556 0.9998
Methods, raw without 238 1 14 0.0460 0.8309 0.0796 0.9883 0.0324 -0.0490 1.0000 -0.0371 1.0000 -0.0150 1.0000 0.0930 1.0000 0.1699 1.0000
Methods, raw without 228 and 238 1 13 0.0329 0.8567 0.0678 0.9914 0.0285 -0.0516 1.0000 -0.0337 1.0000 -0.0163 1.0000 0.0794 1.0000 0.1647 1.0000
Methods, raw, contiguous bands main (132 left out) 2 14 0.9673 0.3287 0.1317 0.9703 0.1020 0.0120 1.0000 0.0402 1.0000 0.0975 1.0000 0.1348 0.9999 0.2258 0.9929
Methods, raw, contiguous bands 132 as WT 3 14 0.9225 0.3399 0.1900 0.9429 0.0811 0.0027 1.0000 0.0037 1.0000 0.0854 1.0000 0.1150 0.9998 0.1986 0.9877
Methods, raw, contiguous bands 132 as Het 2 15 0.9872 0.3236 0.1251 0.9730 0.0994 0.0130 1.0000 0.0448 1.0000 0.0934 1.0000 0.1297 0.9999 0.2162 0.9933
Methods, raw, contiguous bands without 228 2 13 0.8178 0.3692 0.1241 0.9733 0.0974 0.0071 1.0000 0.0435 1.0000 0.0870 1.0000 0.1212 1.0000 0.2280 0.9941
Methods, raw, contiguous bands without 238 1 14 0.1288 0.7208 0.1100 0.9786 0.0528 -0.0281 1.0000 -0.0371 1.0000 0.0112 1.0000 0.0930 1.0000 0.2248 0.9995
Methods, raw, contiguous bands without 228 and 238 1 13 0.0992 0.7538 0.1040 0.9807 0.0481 -0.0329 1.0000 -0.0337 1.0000 0.0007 1.0000 0.0794 1.0000 0.2270 0.9996
Methods, raw, SWD epochs from scorer flags main (132 left out) 2 14 3.3446 0.0717 0.3646 0.8330 0.1654 0.0010 1.0000 0.1035 1.0000 0.1604 0.9985 0.3129 0.8681 0.2493 0.9644
Methods, raw, SWD epochs from scorer flags 132 as WT 3 14 3.8573 0.0532 0.5341 0.7110 0.1447 -0.0057 1.0000 0.0636 1.0000 0.1425 0.9971 0.2890 0.7605 0.2339 0.9171
Methods, raw, SWD epochs from scorer flags 132 as Het 2 15 3.2518 0.0754 0.3391 0.8507 0.1585 0.0023 1.0000 0.1045 0.9999 0.1533 0.9987 0.2968 0.8842 0.2358 0.9703
Methods, raw, SWD epochs from scorer flags without 228 2 13 3.1968 0.0784 0.3369 0.8521 0.1683 -0.0009 1.0000 0.1121 0.9999 0.1659 0.9985 0.3134 0.8916 0.2510 0.9710
Methods, raw, SWD epochs from scorer flags without 238 1 14 0.7746 0.3820 0.2393 0.9151 0.1126 -0.0403 1.0000 0.0218 1.0000 0.0599 1.0000 0.2648 0.9950 0.2570 0.9960
Methods, raw, SWD epochs from scorer flags without 228 and 238 1 13 0.7511 0.3896 0.2184 0.9272 0.1155 -0.0423 1.0000 0.0304 1.0000 0.0654 1.0000 0.2653 0.9962 0.2587 0.9969
Methods, raw, slope rejection main (132 left out) 2 14 0.5835 0.4475 0.0834 0.9873 0.0813 -0.0083 1.0000 0.0402 1.0000 0.0793 1.0000 0.1348 0.9999 0.1608 0.9996
Methods, raw, slope rejection 132 as WT 3 14 0.5384 0.4654 0.1263 0.9725 0.0636 -0.0130 1.0000 0.0037 1.0000 0.0657 1.0000 0.1150 0.9999 0.1464 0.9990
Methods, raw, slope rejection 132 as Het 2 15 0.5997 0.4411 0.0788 0.9886 0.0795 -0.0068 1.0000 0.0448 1.0000 0.0767 1.0000 0.1297 0.9999 0.1530 0.9996
Methods, raw, slope rejection without 228 2 13 0.4896 0.4866 0.0694 0.9910 0.0775 -0.0110 1.0000 0.0435 1.0000 0.0780 1.0000 0.1212 1.0000 0.1556 0.9998
Methods, raw, slope rejection without 238 1 14 0.0460 0.8309 0.0796 0.9883 0.0324 -0.0490 1.0000 -0.0371 1.0000 -0.0150 1.0000 0.0930 1.0000 0.1699 1.0000
Methods, raw, slope rejection without 228 and 238 1 13 0.0329 0.8567 0.0678 0.9914 0.0285 -0.0516 1.0000 -0.0337 1.0000 -0.0163 1.0000 0.0794 1.0000 0.1647 1.0000
Methods, raw, cleaner S1 channel main (132 left out) 2 14 0.1956 0.6596 0.2754 0.8929 0.0273 -0.0767 0.9999 -0.0038 1.0000 0.0302 1.0000 0.0780 0.9999 0.1090 0.9985
Methods, raw, cleaner S1 channel 132 as WT 3 14 0.0167 0.8974 0.4198 0.7939 0.0065 -0.0853 0.9990 -0.0451 1.0000 0.0120 1.0000 0.0571 1.0000 0.0940 0.9978
Methods, raw, cleaner S1 channel 132 as Het 2 15 0.2460 0.6214 0.2570 0.9045 0.0297 -0.0698 1.0000 0.0047 1.0000 0.0318 1.0000 0.0770 0.9999 0.1047 0.9986
Methods, raw, cleaner S1 channel without 228 2 13 0.0869 0.7692 0.2475 0.9102 0.0186 -0.0856 0.9998 -0.0049 1.0000 0.0240 1.0000 0.0597 1.0000 0.0996 0.9994
Methods, raw, cleaner S1 channel without 238 1 14 0.0068 0.9346 0.2085 0.9328 -0.0072 -0.1013 1.0000 -0.0570 1.0000 -0.0451 1.0000 0.0479 1.0000 0.1195 0.9998
Methods, raw, cleaner S1 channel without 228 and 238 1 13 0.0319 0.8588 0.1875 0.9440 -0.0159 -0.1102 0.9999 -0.0581 1.0000 -0.0513 1.0000 0.0297 1.0000 0.1102 0.9999
stored (scorer) main (132 left out) 2 14 0.2877 0.5934 0.1343 0.9692 -0.0266 0.0284 1.0000 -0.0646 0.9999 -0.0638 0.9999 -0.0326 1.0000 -0.0002 1.0000
stored (scorer) 132 as WT 3 14 7.6448 0.0072 0.3771 0.8243 -0.1248 -0.1626 0.8387 -0.1971 0.6333 -0.1355 0.9400 -0.0897 0.9964 -0.0390 1.0000
stored (scorer) 132 as Het 2 15 0.0084 0.9273 0.1206 0.9747 -0.0052 0.0647 1.0000 -0.0338 1.0000 -0.0453 1.0000 -0.0191 1.0000 0.0076 1.0000
stored (scorer) without 228 2 13 0.1552 0.6949 0.1264 0.9724 -0.0201 0.0361 1.0000 -0.0573 1.0000 -0.0568 1.0000 -0.0275 1.0000 0.0051 1.0000
stored (scorer) without 238 1 14 2.2319 0.1400 0.0561 0.9940 0.0987 0.1401 0.9941 0.1209 0.9980 0.1004 0.9996 0.0828 0.9999 0.0493 1.0000
stored (scorer) without 228 and 238 1 13 2.3995 0.1266 0.0565 0.9939 0.1052 0.1477 0.9928 0.1283 0.9975 0.1074 0.9994 0.0880 0.9999 0.0545 1.0000
Result: not reproduced. Normalized wake power is higher in Het than in WT on average across
bands (+0.081; genotype F = 0.58, p = 0.45). Averaged over the bands, Het are higher in every
variant on the raw EEG with both WT (+0.0065 to +0.168, column mean_Het-WT). Het are lower on
average with the stored spectra (-0.027, genotype p = 0.59; p = 0.0072 only when GRIN2B-132 is
counted as WT), and with the cleaner S1 channel without GRIN2B-238 (-0.0072, p = 0.93; -0.016
without both, p = 0.86).
fig5a-swd-state-distribution¶
Over both genotypes, SWDs occur mostly in wake (94.8%), then NREM (4.4%), and rarely in REM sleep (0.5%).
In [30]:
def onsets(rule):
rows = []
for s in sessions:
for d in s["days"]:
o = grin2b.swd_onsets(d, rule)
o["subject"] = s["subject_id"]
rows.append(o)
return pd.concat(rows, ignore_index=True)
on = {rule: onsets(rule) for rule in ("prev", "any6")}
rows = []
for rule, o in on.items():
per = pd.crosstab(o.subject, o.state, normalize="index") * 100
per = per.reindex(columns=[0, 1, 2], fill_value=0)
pooled = o.state.value_counts(normalize=True).reindex([0, 1, 2], fill_value=0) * 100
rows.append({"reading": f"({'a' if rule == 'prev' else 'c'}) mean over animals, NREM rule {rule}", "animals": len(per),
"wake": per[0].mean(), "NREM": per[1].mean(), "REM": per[2].mean()})
rows.append({"reading": f"({'b' if rule == 'prev' else 'd'}) pooled SWDs, NREM rule {rule}", "animals": len(per),
"wake": pooled[0], "NREM": pooled[1], "REM": pooled[2]})
dist = pd.DataFrame(rows).sort_values("reading")
dist["all_within_tolerance"] = ((dist.wake - 94.8).abs() <= 9.48) & ((dist.NREM - 4.4).abs() <= 0.44) & ((dist.REM - 0.5).abs() <= 0.05)
print(dist.round(3).to_string(index=False))
print("Expected 94.8 +- 9.48, 4.4 +- 0.44, 0.5 +- 0.05. SWDs:", len(on["prev"]))
print("Readings (c) and (d) use the Supplementary Methods' rule: wake and REM from the epoch before the SWD,"
" NREM when any of the six epochs before it is NREM.")
per_c = (pd.crosstab(on["any6"].subject, on["any6"].state, normalize="index") * 100).reindex(columns=[0, 1, 2], fill_value=0)
for name, ids in DROPS.items():
for lab, pp in (("(a)", None), ("(c)", per_c)):
q = (pd.crosstab(on["prev"].subject, on["prev"].state, normalize="index") * 100).reindex(columns=[0, 1, 2], fill_value=0) if pp is None else pp
q = q[~q.index.isin(ids)]
print(f"{name}, reading {lab}: wake {q[0].mean():.2f}, NREM {q[1].mean():.2f}, REM {q[2].mean():.2f} ({len(q)} animals)")
per = (pd.crosstab(on["prev"].subject, on["prev"].state, normalize="index") * 100).reindex(columns=[0, 1, 2], fill_value=0)
print(per.join(A.genotype).join(A.swd_n).round(2).to_string())
reading animals wake NREM REM all_within_tolerance
(a) mean over animals, NREM rule prev 16 94.184 5.546 0.271 False
(b) pooled SWDs, NREM rule prev 16 92.865 6.811 0.324 False
(c) mean over animals, NREM rule any6 16 91.172 8.636 0.192 False
(d) pooled SWDs, NREM rule any6 16 89.270 10.503 0.228 False
Expected 94.8 +- 9.48, 4.4 +- 0.44, 0.5 +- 0.05. SWDs: 19776
Readings (c) and (d) use the Supplementary Methods' rule: wake and REM from the epoch before the SWD, NREM when any of the six epochs before it is NREM.
without 228, reading (a): wake 93.80, NREM 5.92, REM 0.29 (15 animals)
without 228, reading (c): wake 90.58, NREM 9.21, REM 0.21 (15 animals)
without 238, reading (a): wake 94.18, NREM 5.55, REM 0.27 (16 animals)
without 238, reading (c): wake 91.17, NREM 8.64, REM 0.19 (16 animals)
without 228 and 238, reading (a): wake 93.80, NREM 5.92, REM 0.29 (15 animals)
without 228 and 238, reading (c): wake 90.58, NREM 9.21, REM 0.21 (15 animals)
0 1 2 genotype swd_n
subject
GRIN2B-129 95.54 4.46 0.00 Het 784.5
GRIN2B-130 95.07 4.54 0.39 Het 507.0
GRIN2B-131 95.89 3.55 0.55 Het 816.0
GRIN2B-132 99.53 0.47 0.00 NaN 107.0
GRIN2B-137 98.81 1.19 0.00 Het 586.5
GRIN2B-138 98.77 1.06 0.17 Het 893.5
GRIN2B-139 92.02 7.55 0.42 Het 946.5
GRIN2B-227 91.64 8.04 0.31 Het 957.5
GRIN2B-228 100.00 0.00 0.00 Het 1.0
GRIN2B-229 95.72 4.06 0.23 Het 443.5
GRIN2B-236 89.91 9.70 0.40 Het 758.0
GRIN2B-237 85.16 13.98 0.87 Het 576.0
GRIN2B-239 87.90 11.60 0.51 Het 884.0
GRIN2B-240 100.00 0.00 0.00 WT 112.5
GRIN2B-241 95.63 4.01 0.36 Het 698.5
GRIN2B-373 85.36 14.52 0.12 Het 816.0
Result: partially reproduced. Under reading (a), 94.2% of SWDs follow wake, 5.5% NREM and 0.27% REM: wake is within tolerance and the order holds, but NREM is above 4.4 +- 0.44 and REM below 0.5 +- 0.05. Reading (c), the Supplementary Methods' rule (wake and REM from the epoch before, NREM from the six epochs before), gives 91.2%, 8.6% and 0.19%, with the same pattern. No reading puts all three within tolerance (pooled, reading (b): 92.9%, 6.8% and 0.32%). Without GRIN2B-228 the values move by at most 0.6 percentage points.
fig5c-transition-swds¶
All Grin2b+/- rats have SWDs that start at wake-NREM transitions, against 3 of 9 WT rats (Fisher's exact test, p = .00062).
In [31]:
o = on["prev"]
tr = o.groupby("subject").agg(swds=("state", "size"), nrem_swds=("state", lambda x: int((x == 1).sum())),
transitional=("transition", "sum"))
tr = A[["genotype", "swd_n"]].join(tr).fillna({"swds": 0, "nrem_swds": 0, "transitional": 0})
tr["has_transitional"] = tr.transitional > 0
print(tr.sort_values("genotype").to_string())
for vname, d in variants("swd_n", tr).items():
tab = pd.crosstab(d.g, d.has_transitional).reindex(index=["Het", "WT"], columns=[True, False], fill_value=0)
p = stats.fisher_exact(tab.to_numpy())[1]
print(f"{vname}: Het {tab.loc['Het', True]}/{tab.loc['Het'].sum()}, WT {tab.loc['WT', True]}/{tab.loc['WT'].sum()}, Fisher p = {p:.4f}")
genotype swd_n swds nrem_swds transitional has_transitional subject GRIN2B-129 Het 784.5 1569.0 70.0 24 True GRIN2B-130 Het 507.0 1014.0 46.0 15 True GRIN2B-131 Het 816.0 1632.0 58.0 17 True GRIN2B-137 Het 586.5 1173.0 14.0 7 True GRIN2B-138 Het 893.5 1787.0 19.0 7 True GRIN2B-139 Het 946.5 1893.0 143.0 59 True GRIN2B-227 Het 957.5 1915.0 154.0 65 True GRIN2B-228 Het 1.0 2.0 0.0 0 False GRIN2B-229 Het 443.5 887.0 36.0 7 True GRIN2B-236 Het 758.0 1516.0 147.0 74 True GRIN2B-237 Het 576.0 1152.0 161.0 49 True GRIN2B-239 Het 884.0 1768.0 205.0 71 True GRIN2B-241 Het 698.5 1397.0 56.0 24 True GRIN2B-373 Het 816.0 1632.0 237.0 66 True GRIN2B-238 WT 0.0 0.0 0.0 0 False GRIN2B-240 WT 112.5 225.0 0.0 0 False GRIN2B-132 NaN 107.0 214.0 1.0 0 False main (132 left out): Het 13/14, WT 0/2, Fisher p = 0.0250 132 as WT: Het 13/14, WT 0/3, Fisher p = 0.0059 132 as Het: Het 13/15, WT 0/2, Fisher p = 0.0441 without 228: Het 13/13, WT 0/2, Fisher p = 0.0095 without 238: Het 13/14, WT 0/1, Fisher p = 0.1333 without 228 and 238: Het 13/13, WT 0/1, Fisher p = 0.0714
Result: partially reproduced. 13 of 14 Het rats have at least one SWD at a wake-NREM transition, against 0 of 2 WT (Fisher p = 0.025). The exception is GRIN2B-228, which has 2 SWDs in two days, none after NREM. GRIN2B-132 has none. Without GRIN2B-228 every Het counts (13 of 13, Fisher p = 0.0095), which would meet the criterion; with GRIN2B-238 left out one WT remains and p = 0.13 (0.071 without both). The outcome is judged on all 16 rats.
figs5b-transition-percent¶
In Grin2b+/- rats 24.2% of SWDs occur during the wake-NREM transition, similar to the 22.39% of the WT rats that have such SWDs (two-sample t test, p = .74).
In [32]:
# Reading (i) follows the Supplementary Methods: an SWD is initiated in NREM when any of the six epochs
# before it is NREM (rule any6), and the share of these that are at a wake-NREM transition is taken.
# Every transitional SWD has an NREM epoch among its six epochs before, so it is in this denominator.
# The first draft of the criterion took NREM from the epoch before (rule prev); that is shown as (i-prev),
# with all transitional SWDs in the numerator, as first run, and with only those whose epoch before is NREM.
o6 = on["any6"]
tr["nrem_swds_any6"] = o6.groupby("subject").state.apply(lambda x: int((x == 1).sum())).reindex(tr.index).fillna(0)
tr["trans_prev_nrem"] = o.groupby("subject").apply(lambda g: int((g.transition & (g.state == 1)).sum())).reindex(tr.index).fillna(0)
tr["pct_i_methods"] = 100 * tr.transitional / tr.nrem_swds_any6.where(tr.nrem_swds_any6 > 0)
tr["pct_i_prev_as_run"] = 100 * tr.transitional / tr.nrem_swds.where(tr.nrem_swds > 0)
tr["pct_i_prev_restricted"] = 100 * tr.trans_prev_nrem / tr.nrem_swds.where(tr.nrem_swds > 0)
tr["pct_ii_all"] = 100 * tr.transitional / tr.swds.where(tr.swds > 0)
cols = ["pct_i_methods", "pct_i_prev_as_run", "pct_i_prev_restricted", "pct_ii_all"]
print(tr[["genotype", "swds", "nrem_swds_any6", "nrem_swds", "transitional", "trans_prev_nrem"] + cols].sort_values("genotype").round(2).to_string())
names = {"pct_i_methods": "(i) of NREM-initiated SWDs, Methods rule (six epochs)",
"pct_i_prev_as_run": "(i-prev) epoch-before NREM rule, as first run",
"pct_i_prev_restricted": "(i-prev) epoch-before NREM rule, numerator restricted to these SWDs",
"pct_ii_all": "(ii) of all SWDs"}
for col, name in names.items():
het = tr[(tr.genotype == "Het")][col].dropna()
wt = tr[(tr.genotype == "WT") & tr.has_transitional][col].dropna()
line = f"{name}: Het mean {het.mean():.2f}% (n {len(het)}, SEM {het.sem():.2f}); WT with transitional SWDs: {wt.round(2).tolist()}"
for dname, ids in DROPS.items():
line += f" | {dname}: {het[~het.index.isin(ids)].mean():.2f}"
print(line)
print("Expected Het 24.2 +- 2.42%")
genotype swds nrem_swds_any6 nrem_swds transitional trans_prev_nrem pct_i_methods pct_i_prev_as_run pct_i_prev_restricted pct_ii_all subject GRIN2B-129 Het 1569.0 90.0 70.0 24 24.0 26.666667 34.285714 34.29 1.529637 GRIN2B-130 Het 1014.0 83.0 46.0 15 15.0 18.072289 32.608696 32.61 1.47929 GRIN2B-131 Het 1632.0 101.0 58.0 17 17.0 16.831683 29.310345 29.31 1.041667 GRIN2B-137 Het 1173.0 26.0 14.0 7 7.0 26.923077 50.0 50.00 0.59676 GRIN2B-138 Het 1787.0 42.0 19.0 7 7.0 16.666667 36.842105 36.84 0.391718 GRIN2B-139 Het 1893.0 222.0 143.0 59 59.0 26.576577 41.258741 41.26 3.116746 GRIN2B-227 Het 1915.0 207.0 154.0 65 65.0 31.400966 42.207792 42.21 3.394256 GRIN2B-228 Het 2.0 0.0 0.0 0 0.0 NaN NaN NaN 0.0 GRIN2B-229 Het 887.0 65.0 36.0 7 7.0 10.769231 19.444444 19.44 0.789177 GRIN2B-236 Het 1516.0 208.0 147.0 74 74.0 35.576923 50.340136 50.34 4.881266 GRIN2B-237 Het 1152.0 251.0 161.0 49 48.0 19.521912 30.434783 29.81 4.253472 GRIN2B-239 Het 1768.0 305.0 205.0 71 71.0 23.278689 34.634146 34.63 4.015837 GRIN2B-241 Het 1397.0 83.0 56.0 24 24.0 28.915663 42.857143 42.86 1.717967 GRIN2B-373 Het 1632.0 392.0 237.0 66 66.0 16.836735 27.848101 27.85 4.044118 GRIN2B-238 WT 0.0 0.0 0.0 0 0.0 NaN NaN NaN NaN GRIN2B-240 WT 225.0 1.0 0.0 0 0.0 0.0 NaN NaN 0.0 GRIN2B-132 NaN 214.0 1.0 1.0 0 0.0 0.0 0.0 0.00 0.0 (i) of NREM-initiated SWDs, Methods rule (six epochs): Het mean 22.93% (n 13, SEM 1.97); WT with transitional SWDs: [] | without 228: 22.93 | without 238: 22.93 | without 228 and 238: 22.93 (i-prev) epoch-before NREM rule, as first run: Het mean 36.31% (n 13, SEM 2.46); WT with transitional SWDs: [] | without 228: 36.31 | without 238: 36.31 | without 228 and 238: 36.31 (i-prev) epoch-before NREM rule, numerator restricted to these SWDs: Het mean 36.27% (n 13, SEM 2.47); WT with transitional SWDs: [] | without 228: 36.27 | without 238: 36.27 | without 228 and 238: 36.27 (ii) of all SWDs: Het mean 2.23% (n 14, SEM 0.44); WT with transitional SWDs: [] | without 228: 2.40 | without 238: 2.23 | without 228 and 238: 2.40 Expected Het 24.2 +- 2.42%
Result: reproduced. The Supplementary Methods judge SWD initiation in NREM on the six epochs before the SWD. Under that rule, 22.93% (SEM 1.97, 13 Het) of the Het SWDs initiated in NREM are at a wake-NREM transition (reading (i)), within 24.2 +- 2.42%. The first draft of the criterion took NREM from the epoch before the SWD, which gives 36.3%; it was corrected after the first verification round. Reading (ii), the share of all SWDs, gives 2.2%. Leaving out GRIN2B-228, which has no SWD initiated in NREM, does not change reading (i). No WT rat has transitional SWDs, so the WT value and the t test cannot be computed.
fig5f-swd-wake-correlation¶
In Grin2b+/- rats, total SWD duration is negatively correlated with total time awake (Pearson, R = -.51, p = .048).
In [33]:
rows = []
for tname, t in (("main", A), ("scorer flags", animals["flag"]), ("stored scores", animals["stored"]),
("SWD time only", animals["exact"])):
for vname, sel in (("Het", t.genotype == "Het"), ("Het and 132", t.genotype.fillna("Het") == "Het"),
("Het without 228", (t.genotype == "Het") & (t.index != "GRIN2B-228"))):
d = t[sel]
r, p = stats.pearsonr(d.swd_total_min, d.wake_min)
rows.append({"scores": tname, "animals": vname, "n": len(d), "R": r, "p": p})
print(A[A.genotype == "Het"][["swd_total_min", "wake_min"]].round(2).sort_values("swd_total_min").to_string())
print(pd.DataFrame(rows).round(4).to_string(index=False))
swd_total_min wake_min
subject
GRIN2B-228 0.01 969.12
GRIN2B-229 20.37 883.08
GRIN2B-137 29.69 900.88
GRIN2B-237 36.53 812.54
GRIN2B-130 37.48 862.54
GRIN2B-241 45.04 861.92
GRIN2B-129 61.59 828.12
GRIN2B-138 62.40 808.17
GRIN2B-131 62.65 753.29
GRIN2B-239 65.80 839.04
GRIN2B-236 69.47 841.00
GRIN2B-227 79.39 765.75
GRIN2B-139 113.50 797.75
GRIN2B-373 124.87 754.00
scores animals n R p
main Het 14 -0.8056 0.0005
main Het and 132 15 -0.8212 0.0002
main Het without 228 13 -0.7344 0.0043
scorer flags Het 14 -0.7523 0.0019
scorer flags Het and 132 15 -0.7839 0.0005
scorer flags Het without 228 13 -0.6660 0.0129
stored scores Het 14 -0.0616 0.8344
stored scores Het and 132 15 -0.2938 0.2878
stored scores Het without 228 13 0.0822 0.7895
SWD time only Het 14 -0.7419 0.0024
SWD time only Het and 132 15 -0.7764 0.0007
SWD time only Het without 228 13 -0.6531 0.0155
Result: reproduced. Over the 14 Het, total SWD time and wake time correlate at R = -0.81 (p = 0.0005), stronger than the paper's -0.51 and inside the tolerance. With the scorer's flags R = -0.75 (p = 0.0019). With only the SWD time removed from wake R = -0.74 (p = 0.0024). With the stored scores, in which SWD epochs count as wake, R = -0.06 (p = 0.83). The correlation therefore comes from SWD time being taken out of wake: in Het, 110.6 min of SWD epochs per day against 57.8 min of SWD time. Without GRIN2B-228 R = -0.73 (p = 0.0043).
fig5def-no-other-correlations¶
In Grin2b+/- rats, total SWD duration is not correlated with total REM or NREM time, the number of REM, NREM or wake bouts, or the mean REM, NREM or wake bout duration (Pearson, all p >= .066).
In [34]:
measures = ["rem_min", "nrem_min", "rem_bouts", "nrem_bouts", "wake_bouts", "rem_bout_s", "nrem_bout_s", "wake_bout_s"]
paper_r = {"rem_min": -.015, "nrem_min": -.13, "rem_bouts": -.19, "nrem_bouts": .49, "wake_bouts": .17,
"rem_bout_s": .17, "nrem_bout_s": -.37, "wake_bout_s": -.25}
rows = []
for tname, t in (("main", A), ("scorer flags", animals["flag"]), ("stored scores", animals["stored"]),
("bridged bouts", animals["bridged"]), ("SWD time only", animals["exact"]),
("main without 228", A[A.index != "GRIN2B-228"])):
d = t[t.genotype == "Het"]
for m in measures:
r, p = stats.pearsonr(d.swd_total_min, d[m])
rows.append({"scores": tname, "measure": m, "n": len(d), "R": r, "p": p, "paper_R": paper_r[m]})
corr = pd.DataFrame(rows)
print(corr.round(4).to_string(index=False))
print("\nTests with p < 0.05:", corr.groupby("scores").p.apply(lambda p: int((p < 0.05).sum())).to_dict())
scores measure n R p paper_R
main rem_min 14 -0.0697 0.8129 -0.015
main nrem_min 14 0.0152 0.9589 -0.130
main rem_bouts 14 0.0581 0.8436 -0.190
main nrem_bouts 14 0.3943 0.1630 0.490
main wake_bouts 14 0.5409 0.0458 0.170
main rem_bout_s 14 -0.3906 0.1673 0.170
main nrem_bout_s 14 -0.2183 0.4533 -0.370
main wake_bout_s 14 -0.5840 0.0283 -0.250
scorer flags rem_min 14 -0.0667 0.8207 -0.015
scorer flags nrem_min 14 0.0829 0.7782 -0.130
scorer flags rem_bouts 14 0.0335 0.9094 -0.190
scorer flags nrem_bouts 14 0.3326 0.2453 0.490
scorer flags wake_bouts 14 0.8243 0.0003 0.170
scorer flags rem_bout_s 14 -0.3494 0.2207 0.170
scorer flags nrem_bout_s 14 -0.1540 0.5993 -0.370
scorer flags wake_bout_s 14 -0.7243 0.0034 -0.250
stored scores rem_min 14 -0.0660 0.8227 -0.015
stored scores nrem_min 14 0.1052 0.7204 -0.130
stored scores rem_bouts 14 0.0264 0.9285 -0.190
stored scores nrem_bouts 14 0.1139 0.6983 0.490
stored scores wake_bouts 14 0.1280 0.6629 0.170
stored scores rem_bout_s 14 -0.3360 0.2402 0.170
stored scores nrem_bout_s 14 -0.0189 0.9488 -0.370
stored scores wake_bout_s 14 -0.0672 0.8195 -0.250
bridged bouts rem_min 14 -0.0697 0.8129 -0.015
bridged bouts nrem_min 14 0.0152 0.9589 -0.130
bridged bouts rem_bouts 14 0.0264 0.9285 -0.190
bridged bouts nrem_bouts 14 -0.1022 0.7280 0.490
bridged bouts wake_bouts 14 -0.0925 0.7532 0.170
bridged bouts rem_bout_s 14 -0.3427 0.2303 0.170
bridged bouts nrem_bout_s 14 0.0723 0.8060 -0.370
bridged bouts wake_bout_s 14 -0.2276 0.4339 -0.250
SWD time only rem_min 14 -0.0669 0.8202 -0.015
SWD time only nrem_min 14 0.0663 0.8218 -0.130
SWD time only rem_bouts 14 0.0581 0.8436 -0.190
SWD time only nrem_bouts 14 0.3943 0.1630 0.490
SWD time only wake_bouts 14 0.5409 0.0458 0.170
SWD time only rem_bout_s 14 -0.3906 0.1673 0.170
SWD time only nrem_bout_s 14 -0.2183 0.4533 -0.370
SWD time only wake_bout_s 14 -0.5840 0.0283 -0.250
main without 228 rem_min 13 -0.1058 0.7309 -0.015
main without 228 nrem_min 13 -0.1183 0.7004 -0.130
main without 228 rem_bouts 13 -0.0370 0.9046 -0.190
main without 228 nrem_bouts 13 0.4312 0.1412 0.490
main without 228 wake_bouts 13 0.2898 0.3369 0.170
main without 228 rem_bout_s 13 -0.2130 0.4848 0.170
main without 228 nrem_bout_s 13 -0.2952 0.3275 -0.370
main without 228 wake_bout_s 13 -0.4990 0.0825 -0.250
Tests with p < 0.05: {'SWD time only': 2, 'bridged bouts': 0, 'main': 2, 'main without 228': 0, 'scorer flags': 2, 'stored scores': 0}
Result: partially reproduced. Six of the eight measures are not correlated with total SWD time (p from 0.16 to 0.96), but the number of wake bouts (R = 0.54, p = 0.046) and the mean wake bout duration (R = -0.58, p = 0.028) are: two of eight. With the scorer's flags the same two correlate (R = 0.82 and -0.72); with the stored scores, with bouts not split by SWD epochs, or without GRIN2B-228, none of the eight does (wake bouts R = 0.13, -0.09 and 0.29).
Blocked claims¶
Two genotype comparisons need at least two animals with SWDs in each genotype, and only one WT rat (GRIN2B-240) has SWDs; the paper's other seven WT are not in the dandiset. The cell below also runs both analyses with GRIN2B-132 counted as WT, as an illustration only, since its genotype is unknown:
In [35]:
with_swd = A[A.swd_n > 0]
print("Animals with SWDs by genotype:", with_swd.genotype.fillna("none").value_counts().to_dict())
print("Animals with >= 5 usable SWD epochs in the raw windows:",
{g: int(n) for g, n in pd.Series({sid: r["n"]["SWD"] for sid, r in spectra.items()}).ge(5).groupby(
pd.Series({sid: by_id[sid]["genotype"] or "none" for sid in spectra})).sum().items()})
# Not a test of either claim: GRIN2B-132 has no genotype. With 132 counted as WT, two WT have SWDs.
print("\nfig2e-swd-delta-power, normalized SWD band power, animals with >= 5 usable SWD epochs:")
n_swd = {sid: r["n"]["SWD"] for sid, r in spectra.items()}
Tswd = band_table("SWD")
Tswd = Tswd[Tswd.subject.map(n_swd) >= 5]
bw = band_anova(Tswd, "Methods, raw, SWD")
print(bw[bw.variant.isin(["main (132 left out)", "132 as WT"])].round(4).to_string(index=False))
print("\nfig5a-het-more-nrem-swds, per-animal SWD state percentages (reading (a)), 132 as WT:")
pa = (pd.crosstab(on["prev"].subject, on["prev"].state, normalize="index") * 100).reindex(columns=[0, 1, 2], fill_value=0)
pa.columns = ["Wake", "NREM", "REM"]
pa = pa.join(A.genotype.fillna("WT").rename("g"))
print(pa.groupby("g").agg(["mean", "size"]).round(2).to_string())
long = pa.reset_index().melt(id_vars=["subject", "g"], var_name="state", value_name="pct")
a = sm.stats.anova_lm(ols("pct ~ C(g) * C(state)", data=long).fit(), typ=2)
print(a.round(4).to_string())
tk = pairwise_tukeyhsd(long.pct, long.g + ":" + long.state)
tk = pd.DataFrame(tk.summary().data[1:], columns=tk.summary().data[0])
for st_ in ("Wake", "NREM", "REM"):
r = tk[tk.group1.str.endswith(":" + st_) & tk.group2.str.endswith(":" + st_)]
print(f"Tukey {st_} Het-WT: p = {float(r['p-adj'].iloc[0]):.4f}")
print(f"\ntotal run time {time.time() - started:.0f} s")
Animals with SWDs by genotype: {'Het': 14, 'none': 1, 'WT': 1}
Animals with >= 5 usable SWD epochs in the raw windows: {'Het': 13, 'WT': 1, 'none': 1}
fig2e-swd-delta-power, normalized SWD band power, animals with >= 5 usable SWD epochs:
spectra variant n_WT n_Het genotype_F genotype_p interaction_F interaction_p mean_Het-WT delta_Het-WT delta_tukey_p theta_Het-WT theta_tukey_p sigma_Het-WT sigma_tukey_p beta_Het-WT beta_tukey_p gamma_Het-WT gamma_tukey_p
Methods, raw, SWD main (132 left out) 1 13 5.1120 0.0274 1.5717 0.1935 2.7371 -0.1860 1.0 0.0410 1.0 1.0085 1.0000 6.7786 0.2874 6.0436 0.4469
Methods, raw, SWD 132 as WT 2 13 10.5508 0.0018 2.9425 0.0268 2.7673 -0.3423 1.0 0.3121 1.0 1.3405 0.9994 6.7825 0.0227 5.7435 0.0967
fig5a-het-more-nrem-swds, per-animal SWD state percentages (reading (a)), 132 as WT:
Wake NREM REM
mean size mean size mean size
g
Het 93.39 14 6.30 14 0.31 14
WT 99.77 2 0.23 2 0.00 2
sum_sq df F PR(>F)
C(g) 0.0000 1.0 0.0000 1.0000
C(state) 89089.7172 2.0 3144.9194 0.0000
C(g):C(state) 135.8917 2.0 4.7971 0.0133
Residual 594.8909 42.0 NaN NaN
Tukey Wake Het-WT: p = 0.2407 Tukey NREM Het-WT: p = 0.2903 Tukey REM Het-WT: p = 1.0000 total run time 369 s
With GRIN2B-132 as WT, SWD delta power is lower in Het than in the two WT (-0.34, Tukey p = 1.0), against the paper's higher delta in Het; the genotype by band interaction has p = 0.027, from beta (+6.8, Tukey p = 0.023). The share of SWDs after NREM is 6.30% in Het against 0.23% in the two WT (interaction p = 0.013, Tukey NREM p = 0.29). Neither is a test of the claims.
Summary¶
| Claim | Outcome |
|---|---|
| fig2b-swd-hourly | reproduced |
| fig2c-swd-number | reproduced |
| fig2c-swd-mean-duration | partially reproduced |
| fig2c-swd-total-duration | reproduced |
| disc-het-swd-minutes | partially reproduced |
| figs5a-swd-light-dark | reproduced |
| figs8-swd-sex | reproduced |
| fig2d-swd-spectrum-peak | reproduced |
| fig2e-swd-delta-power | blocked (data missing: one WT with SWDs) |
| fig3d-rem-hourly | reproduced |
| fig3g-rem-total | reproduced |
| fig3g-rem-bouts | reproduced |
| fig3g-rem-bout-duration | partially reproduced |
| fig3h-nrem-total | reproduced |
| fig3h-nrem-bouts | partially reproduced |
| fig3i-wake-total | partially reproduced |
| fig3i-wake-bouts | partially reproduced |
| fig4a-rem-spectrum | reproduced |
| fig4b-nrem-beta | not reproduced |
| fig4c-wake-power | not reproduced |
| fig5a-swd-state-distribution | partially reproduced |
| fig5a-het-more-nrem-swds | blocked (data missing: one WT with SWDs) |
| fig5c-transition-swds | partially reproduced |
| figs5b-transition-percent | reproduced |
| fig5f-swd-wake-correlation | reproduced |
| fig5def-no-other-correlations | partially reproduced |
| fig1c-glun2b-reduced, fig1df-other-proteins-unchanged | blocked (data missing) |
| figs4-normal-development, figs1-scorer-agreement | blocked (data missing) |
| fig6b, fig6c, fig6d, fig7c (pharmacology) | blocked (data missing) |
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.
- Two of the paper's nine WT rats are in the dandiset, so every genotype test has two WT, and one of them has no SWDs. Several outcomes are partial because the direction holds without p < 0.05.
- The dandiset does not say which rats the paper analysed; all 16 with a genotype are used, and GRIN2B-132 (no genotype) is added to each group as a variant. Without GRIN2B-228 (Het, 1 SWD per day) fig5c-transition-swds and fig5def-no-other-correlations would meet their criteria, and fig3d-rem-hourly and fig3g-rem-total would not (p = 0.060). Without GRIN2B-238 one WT remains, and the genotype by hour interaction of fig3d-rem-hourly has p = 0.0001.
- The criterion of figs5b-transition-percent first took the state of an SWD from the epoch before it; it was corrected after the first verification round to the Supplementary Methods, which judge initiation in NREM on the six epochs before. The first rule gives 36.3%, outside tolerance.
- SWD epochs are the epochs SWD events overlap (the Methods); the lab's scorer uses its per-epoch flags, which cover about half of these epochs. The wake bout results depend on this choice: with bouts not split by SWD epochs, or with the stored scores, Het have fewer and longer wake bouts, and the correlations of SWD time with wake bouts disappear.
- Spectra are computed from about 4.8 h of raw EEG per rat (a median of 5 to 10% of each state's epochs) on the SWD-detection channel, with band limits and a Welch segment length chosen here. The spectral-slope rejection, as read here, removes no epoch; the per-animal cleaner S1 channel is a variant.
- Mixed models are tested as split-plot ANOVAs, which give the same genotype and interaction F as a fitted random-intercept model; the sex ANOVA has no interaction term.
- Bouts, the state before an SWD and wake-NREM transitions follow definitions chosen here.
What was found about the dandiset.
- Genotype and sex are missing for GRIN2B-132. The dandiset holds 2 of the paper's 9 WT and 14 or 15 of its 15 Het, and does not say which rats the paper analysed.
- The
sleep_statesdescription promises code 4 for SWD epochs, which no file holds. swd_epoch_countsholds flags (4 and 6), not counts; the flags cover about half of the epochs SWD events overlap, mostly the last epoch of each SWD.- The
seizure_totals_by_ztcolumns are shifted: per-hour counts are inspike_wave_discharge, and the day's total is in the first row ofnumber_of_seizure_eventsandDay. - The stored spectra are the scorer's log amplitude referenced to the 5th percentile, not the baseline-corrected PSD of their description.
- The raw EEG is stored as unsigned 12-bit codes with an offset of 0 (about +8 mV of DC), and signal loss is written as 32767 on all channels, which the description does not mention; up to 73% of epochs in the windows read are affected.
- The pharmacology, western blot, development and visual-scoring data are not in the dandiset, and there is no published version.