Reproducing Tingley et al. (2021) from dandiset 000233¶
Paper. Tingley D, McClain K, Kaya E, Carpenter J, Buzsáki G (2021). A metabolic function of the hippocampal sharp wave-ripple. Nature 597, 82-86. doi:10.1038/s41586-021-03811-w. The claims file was written from the accepted author manuscript in PubMed Central (PMC9214835) with its Methods, Extended Data legends and Supplementary Information.
Dandiset. DANDI:000233, version 0.250624.0414. The paper cites no DANDI version (its data statement points to the lab's server). This is the most recent published version; the earlier one, 0.230223.0815, lists the same 345 files with the same sizes, so the choice does not change the data.
What the dandiset holds. 345 NWB files from 25 rats: CA1 LFP at 1,250 Hz, tables of detected SPW-Rs, the iPro2 glucose sensor current, a head accelerometer, and sleep states for six rats. It holds no spike units, no optogenetic stimulation times, no CNO or vehicle labels and injection times, no meal times, and no cohort labels, so the optogenetic, DREADD, fasting and unit claims cannot be tested.
Random numbers are drawn only for circular shuffles, with a fixed seed in each claim's section.
| Claim | Status |
|---|---|
| fig2-ripple-rate-predicts-glucose | tested |
| fig2ef-correlation-beyond-shuffle | tested |
| fig2g-clustered-ripples | tested |
| edfig4-ripple-features | tested |
| fig2h-phase-reset | tested |
| edfig5-cross-frequency-coupling | tested |
| edfig1-glucose-rhythm | tested |
| edfig1-glucose-circadian-meals | blocked: meal times not in the dandiset (circadian part run) |
| edfig1-sensor-decay | tested |
| supp-ripple-rate | tested |
| supp-ripples-stable-across-days | tested |
| edfig3-ripple-rate-rhythms | tested |
| edfig2-refractory-period | tested |
| edfig4-nrem-wake | tested |
| edfig4-no-hour-differs | tested |
| edfig9-dorsal-ventral | blocked: glucose with the electrophysiology in one of three dorsal and ventral rats (run) |
| edfig3e-ca1-firing-decrease, edfig8-lateral-septum-rectified | blocked: no spike units |
| edfig4-fasting-fed | blocked: no meal times |
| edfig6-ripple-rate-best-predictor, edfig6d-pairwise-predictors, edfig6ef-latent-variables | blocked: too large to compute the predictors (timed below) |
| fig3d, fig3-induced-effect-size, fig3ef, fig3f, edfig7ab, edfig7jm, edfig7-ppc | blocked: no stimulation times |
| fig4-cno-controls-no-change, fig4-ls-dreadd-decouples, fig4-ls-dreadd-chance-level, fig4e-dreadd-time-course, fig4c-medial-septum-no-effect, edfig10-expression-correlates-effect | blocked: no CNO or vehicle labels and no cohort labels |
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import signal, stats
import sprglu as S
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()
assets = S.list_assets()
meta = S.version_metadata()
files = S.read_all(assets)
# the glucose series is stored in every file; all are read to check that they agree
glu_all = S.read_all_glucose(assets)
print(f"dandiset {S.DANDISET} version {S.VERSION}: {len(assets)} assets, "
f"{assets['size'].sum() / 1e12:.2f} TB; read in {time.time() - started:.0f} s")
dandiset 000233 version 0.250624.0414: 345 assets, 12.32 TB; read in 57 s
Data check¶
The archive's metadata and the files are checked against each other and against the paper's description. Each check is stated under its output.
print("Archive summary:", meta["assetsSummary"]["numberOfSubjects"], "subjects,",
meta["assetsSummary"]["numberOfFiles"], "files;", meta["assetsSummary"]["variableMeasured"])
fi = pd.DataFrame([{
"subject": f["subject"], "path": f["path"].split("/")[1], "experiment": f["experiment"][:40],
"lfp_t0": f["lfp_t0"], "hours": f["lfp_n"] / f["lfp_rate"] / 3600, "lfp_rate": f["lfp_rate"],
"channels": f["lfp_channels"], "acquisition": "+".join(f["acquisition"]),
"tables": "+".join(sorted(f["tables"])), "units": f["units"],
"sleep_states": f["sleep"] is not None and len(f["sleep"]) > 0} for f in files])
print("LFP rates:", sorted(fi.lfp_rate.unique()), " units tables in any file:", fi.units.any())
print("experiment_description values:", fi.experiment.value_counts().to_dict())
print("acquisition series:", fi.acquisition.value_counts().to_dict())
print("SPW-R tables per file:", fi.tables.value_counts().to_dict())
print("Full description strings (distinct values, number of files):")
for key in ("experiment", "session_description", "subject_description"):
for v, n in pd.Series([f[key] for f in files]).value_counts().items():
print(f" {key} ({n} files): {v!r}")
gp = pd.DataFrame([(f["subject"], g) for f in files for g in f["groups"]], columns=["subject", "group"])
print("Groups and members under /intervals, /stimulus, /analysis, /processing, /general (two levels), number of subjects and files:")
print(gp.groupby("group").agg(subjects=("subject", "nunique"), files=("subject", "size")).to_string())
cond = gp[gp.group.str.match(r"/(intervals|stimulus/.+/|analysis/)")]
print("members inside /intervals, /stimulus/presentation, /stimulus/templates or /analysis:", len(cond))
print("electrode table columns:", sorted({c for f in files for c in f["electrode_columns"]}),
"; locations:", sorted({c for f in files for c in f["electrode_locations"]}),
"; groups:", sorted({c for f in files for c in f["electrode_groups"]}))
cn = sorted({c for f in files for c in f["channel_names"]})
print(f"{len(cn)} distinct channel names, e.g. {cn[:5]} ... {cn[-5:]}; names that are not a number:",
[c for c in cn if not c.isdigit()])
print("attributes of the SPW-R tables:", sorted({a for f in files for t in f["table_attrs"].values() for a in t}))
Archive summary: 25 subjects, 345 files; ['ElectricalSeries', 'ProcessingModule', 'LFP', 'ElectrodeGroup']
LFP rates: [np.float64(1250.0)] units tables in any file: False
experiment_description values: {'Consult Supplementary Table 1 from the p': 324, "Experiment condition 'ripple/glucose rec": 21}
acquisition series: {'Accelerometer+ElectricalSeriesRaw+GlucoseLevel': 332, 'ElectricalSeriesRaw+GlucoseLevel': 10, 'Accelerometer+GlucoseLevel': 3}
SPW-R tables per file: {'Ca1Ripples': 221, '': 63, 'Ripples': 22, 'Dorsalripples+Globalripples+Ventralripples': 12, 'Ca1Ripples+Nonoise_Ca1Ripples': 12, 'Globalripples+Ventralripples': 4, 'Ventralripples': 3, 'Dorsalripples+Globalripples': 3, 'Nonoise_Ca1Ripples': 3, '_Remove_Ca1Ripples': 1, 'Globalripples': 1}
Full description strings (distinct values, number of files):
experiment (324 files): 'Consult Supplementary Table 1 from the publication for more information about this session.'
experiment (21 files): "Experiment condition 'ripple/glucose recording' with surgery condition 'rCA1'."
session_description (296 files): 'Consult Supplementary Table 1 from the publication for more information about this session.'
session_description (49 files): 'no description'
subject_description (345 files): 'A total of 45 adult Long Evans rats were used in this study. Five were used for the hippocampal and lateral septal in vivo recordings26. Ten rats were used for chronic electrophysiological recordings paired with glucose monitoring. Eight rats were used for optogenetic induction of SPW-Rs and glucose monitoring. Six rats were used for the DREADD experiment, and seven rats were used as controls for this experiment. Three additional rats were used for simultaneous dorsal and ventral hippocampus recordings. Six additional rats were used for the posterior parietal cortex (PPC) optogenetic stimulation control experiments. Sample sizes were selected to match cohort sizes where applicable.\n'
Groups and members under /intervals, /stimulus, /analysis, /processing, /general (two levels), number of subjects and files:
subjects files
group
/general/devices 25 345
/general/devices/IntanDevice 25 345
/general/devices/Medtronic iPro2 CGM 25 345
/general/experiment_description 25 345
/general/experimenter 25 345
/general/extracellular_ephys 25 345
/general/institution 25 345
/general/keywords 25 345
/general/lab 25 345
/general/related_publications 25 345
/general/session_id 25 345
/general/subject 25 345
/general/subject/description 25 345
/general/subject/sex 25 345
/general/subject/species 25 345
/general/subject/strain 25 345
/general/subject/subject_id 25 345
/processing/behavior 6 86
/processing/behavior/sleep_states 6 86
/processing/ecephys 25 345
/processing/ecephys/Ca1Ripples 18 233
/processing/ecephys/Dorsalripples 3 15
/processing/ecephys/Globalripples 3 20
/processing/ecephys/LFP 25 345
/processing/ecephys/Nonoise_Ca1Ripples 1 15
/processing/ecephys/Ripples 1 22
/processing/ecephys/Ventralripples 3 19
/processing/ecephys/_Remove_Ca1Ripples 1 1
/stimulus/presentation 25 345
/stimulus/templates 25 345
members inside /intervals, /stimulus/presentation, /stimulus/templates or /analysis: 0
electrode table columns: ['channel_name', 'filtering', 'gain_to_uV', 'group', 'group_name', 'id', 'imp', 'location', 'rel_x', 'rel_y', 'shank_electrode_number', 'x', 'y', 'z'] ; locations: ['unknown'] ; groups: ['Group1', 'Group2', 'Group3', 'Group4', 'Group5', 'Group6', 'Group7', 'Group8']
256 distinct channel names, e.g. ['0', '1', '10', '100', '101'] ... ['95', '96', '97', '98', '99']; names that are not a number: []
attributes of the SPW-R tables: ['colnames', 'description', 'neurodata_type']
No file holds spike units, optogenetic stimulation times, injection times or labels (CNO or
vehicle), meal times, or a cohort label: only the files of CGM1 name their experiment
("ripple/glucose recording", surgery rCA1), the rest point to Supplementary Table 1, whose rats
are numbered 1 to 45 while the dandiset's subjects are named CGM1 to CGM60, bruce, dt15, flex1,
ros and Vanessa, with no mapping between the two. SPW-R tables come under seven names; the
dorsal and ventral tables identify the dorsal and ventral rats. The full description strings say
no more: session_description is "no description" or points to Supplementary Table 1, and
subject/description is the paper's count of 45 rats over all cohorts, the same in every file.
/intervals is absent, /stimulus/presentation, /stimulus/templates and /analysis are empty
in every file, /processing holds only the LFP, the SPW-R tables and (six subjects) sleep states,
and the electrode tables name channels by number with location "unknown", so no file holds a
stimulation, TTL or condition channel. The SPW-R tables carry no attributes beyond their column
names and description: the detection parameters the Methods say are stored with each session's
ripple events (in the lab's MATLAB files) are not in the NWB files.
dups = S.find_duplicates(files)
drop = {i for g in dups for i in g[1:]}
print(f"{len(dups)} recordings are deposited twice (same subject, same LFP start and sample count):")
dd = pd.DataFrame([{"subject": files[g[0]]["subject"], "file_1": files[g[0]]["path"].split("/")[1],
"file_2": files[g[1]]["path"].split("/")[1],
"hours": files[g[0]]["lfp_n"] / 1250 / 3600,
"same_events": all(np.array_equal(files[g[0]]["tables"][k].peak.values,
files[g[1]]["tables"][k].peak.values)
for k in files[g[0]]["tables"] if k in files[g[1]]["tables"])}
for g in dups])
print(dd.round(2).to_string())
kept = [f for i, f in enumerate(files) if i not in drop]
per_subject = fi.drop(index=list(drop)).groupby("subject").agg(
files=("path", "size"), hours=("hours", "sum"), with_sleep_states=("sleep_states", "sum"))
print(per_subject.round(1).to_string())
print("subjects:", len(per_subject), " files counted once:", len(kept), " hours of LFP:", round(per_subject.hours.sum(), 1))
27 recordings are deposited twice (same subject, same LFP start and sample count):
subject file_1 file_2 hours same_events
0 CGM1 sub-CGM1_ses-CGM1-0um-181130-112307_ecephys.nwb sub-CGM1_ses-CGM1_0um_181130_112307_ecephys.nwb 0.89 True
1 CGM1 sub-CGM1_ses-CGM1-144um-181130-121628_ecephys.nwb sub-CGM1_ses-CGM1_144um_181130_121628_ecephys.nwb 0.37 True
2 CGM1 sub-CGM1_ses-CGM1-288um-181130-123859_ecephys.nwb sub-CGM1_ses-CGM1_288um_181130_123859_ecephys.nwb 1.99 True
3 CGM1 sub-CGM1_ses-CGM1-432um-181130-143939_ecephys.nwb sub-CGM1_ses-CGM1_432um_181130_143939_ecephys.nwb 2.59 True
4 CGM1 sub-CGM1_ses-CGM1-504um-181130-180733_ecephys.nwb sub-CGM1_ses-CGM1_504um_181130_180733_ecephys.nwb 0.95 True
5 CGM1 sub-CGM1_ses-CGM1-720um-181130-194600_ecephys.nwb sub-CGM1_ses-CGM1_720um_181130_194600_ecephys.nwb 15.25 True
6 CGM1 sub-CGM1_ses-CGM1-792um-181201-110104_ecephys.nwb sub-CGM1_ses-CGM1_792um_181201_110104_ecephys.nwb 23.77 True
7 CGM3 sub-CGM3_ses-CGM3-0um-190301-122633_ecephys.nwb sub-CGM3_ses-CGM3_0um_190301_122633_ecephys.nwb 1.39 True
8 CGM3 sub-CGM3_ses-CGM3-288um-190301-142051_ecephys.nwb sub-CGM3_ses-CGM3_288um_190301_142051_ecephys.nwb 0.80 True
9 CGM3 sub-CGM3_ses-CGM3-537um-190301-150904_ecephys.nwb sub-CGM3_ses-CGM3_537um_190301_150904_ecephys.nwb 1.80 True
10 CGM3 sub-CGM3_ses-CGM3-609um-190301-183547_ecephys.nwb sub-CGM3_ses-CGM3_609um_190301_183547_ecephys.nwb 0.01 True
11 CGM3 sub-CGM3_ses-CGM3-681um-190302-212204_ecephys.nwb sub-CGM3_ses-CGM3_681um_190302_212204_ecephys.nwb 13.22 True
12 CGM3 sub-CGM3_ses-CGM3-72um-190301-134959_ecephys.nwb sub-CGM3_ses-CGM3_72um_190301_134959_ecephys.nwb 0.51 True
13 CGM30 sub-CGM30_ses-CGM30-0um-210123-103256_ecephys.nwb sub-CGM30_ses-CGM30_0um_210123_103256_ecephys.nwb 0.52 True
14 CGM30 sub-CGM30_ses-CGM30-0um-210131-101140_ecephys.nwb sub-CGM30_ses-CGM30_0um_210131_101140_ecephys.nwb 22.88 True
15 CGM30 sub-CGM30_ses-CGM30-0um-210202-081800_ecephys.nwb sub-CGM30_ses-CGM30_0um_210202_081800_ecephys.nwb 0.50 True
16 CGM31 sub-CGM31_ses-CGM31-0um-210121-143156_ecephys.nwb sub-CGM31_ses-CGM31_0um_210121_143156_ecephys.nwb 0.75 True
17 CGM31 sub-CGM31_ses-CGM31-288um-210121-151659_ecephys.nwb sub-CGM31_ses-CGM31_288um_210121_151659_ecephys.nwb 0.44 True
18 CGM31 sub-CGM31_ses-CGM31-720um-210121-154347_ecephys.nwb sub-CGM31_ses-CGM31_720um_210121_154347_ecephys.nwb 1.17 True
19 CGM31 sub-CGM31_ses-CGM31-720um-210122-074732_ecephys.nwb sub-CGM31_ses-CGM31_720um_210122_074732_ecephys.nwb 1.21 True
20 CGM31 sub-CGM31_ses-CGM31-720um-210122-090353_ecephys.nwb sub-CGM31_ses-CGM31_720um_210122_090353_ecephys.nwb 3.64 True
21 CGM32 sub-CGM32_ses-CGM32-0um-210121-144055_ecephys.nwb sub-CGM32_ses-CGM32_0um_210121_144055_ecephys.nwb 0.61 True
22 CGM32 sub-CGM32_ses-CGM32-0um-210121-151751_ecephys.nwb sub-CGM32_ses-CGM32_0um_210121_151751_ecephys.nwb 0.09 True
23 CGM32 sub-CGM32_ses-CGM32-288um-210121-152327_ecephys.nwb sub-CGM32_ses-CGM32_288um_210121_152327_ecephys.nwb 0.44 True
24 CGM32 sub-CGM32_ses-CGM32-576um-210121-155015_ecephys.nwb sub-CGM32_ses-CGM32_576um_210121_155015_ecephys.nwb 1.80 True
25 CGM32 sub-CGM32_ses-CGM32-576um-210121-181445_ecephys.nwb sub-CGM32_ses-CGM32_576um_210121_181445_ecephys.nwb 14.84 True
26 CGM32 sub-CGM32_ses-CGM32-576um-210123-091140_ecephys.nwb sub-CGM32_ses-CGM32_576um_210123_091140_ecephys.nwb 0.74 True
files hours with_sleep_states
subject
CGM1 21 259.2 21
CGM2 12 59.0 12
CGM3 29 351.1 28
CGM30 19 174.4 0
CGM31 28 330.5 0
CGM32 36 307.4 0
CGM36 17 156.8 0
CGM37 15 36.9 0
CGM39 19 132.0 0
CGM40 7 75.1 0
CGM41 28 297.3 0
CGM46 4 49.9 0
CGM48 6 3.9 0
CGM49 14 104.7 0
CGM51 7 25.0 0
CGM52 6 43.1 0
CGM55 6 34.2 0
CGM57 4 18.9 0
CGM58 12 113.5 0
CGM60 3 11.3 0
Vanessa 4 42.9 0
bruce 9 48.6 4
dt15 8 28.9 8
flex1 1 16.7 1
ros 3 22.9 0
subjects: 25 files counted once: 318 hours of LFP: 2744.1
27 recordings are deposited twice, under a file name with dashes and one with underscores (subjects CGM1, CGM3, CGM30, CGM31, CGM32); both copies have the same SPW-R events. Each is counted once from here on.
rows = []
for f in kept:
named = S.file_start_from_name(f["path"])
on_clock = S.reference_time(f["reference"]) + pd.Timedelta(seconds=f["lfp_t0"])
row = {"subject": f["subject"], "name_minus_clock_s": (on_clock - named).total_seconds() if named else np.nan}
for k, t in f["tables"].items():
if len(t):
row["start_minus_t0_plus_peak_s"] = float(np.median(t.start_time - (f["lfp_t0"] + t.peak)))
row["peak_max_minus_duration_s"] = float(t.peak.max() - f["lfp_n"] / f["lfp_rate"])
rows.append(row)
ck = pd.DataFrame(rows)
print("files whose name has no date:", ck.name_minus_clock_s.isna().sum(),
"; |file-name time - session clock| max (s):", ck.name_minus_clock_s.abs().max())
print("SPW-R start_time - (LFP start + peak), s: ", ck.start_minus_t0_plus_peak_s.describe()[["min", "50%", "max"]].round(4).to_dict())
print("last SPW-R peak - file duration, s (<= 0 means peaks are within the file):",
round(ck.peak_max_minus_duration_s.max(), 3))
ev_all = pd.concat([t.assign(subject=f["subject"], table=k) for f in kept for k, t in f["tables"].items()])
print("share of stored CA1 events below 5 s.d. per subject:",
{s_: round(float((g.peak_normed_power < 5).mean()), 3) for s_, g in ev_all[ev_all.table.isin(S.CA1_TABLES)].groupby("subject")})
print(ev_all.groupby("table").agg(events=("peak", "size"), min_normed_power=("peak_normed_power", "min"),
min_duration_ms=("duration", lambda x: 1000 * x.min()), max_duration_ms=("duration", lambda x: 1000 * x.max()),
median_frequency=("peak_frequency", "median"), min_frequency=("peak_frequency", "min"),
max_frequency=("peak_frequency", "max")).round(2).to_string())
files whose name has no date: 1 ; |file-name time - session clock| max (s): 0.0
SPW-R start_time - (LFP start + peak), s: {'min': -0.0272, '50%': -0.0152, 'max': -0.0096}
last SPW-R peak - file duration, s (<= 0 means peaks are within the file): -0.028
share of stored CA1 events below 5 s.d. per subject: {'CGM1': 0.0, 'CGM2': 0.0, 'CGM3': 0.068, 'CGM30': 0.289, 'CGM31': 0.187, 'CGM32': 0.002, 'CGM37': 0.209, 'CGM39': 0.16, 'CGM40': 0.0, 'CGM41': 0.001, 'CGM46': 0.0, 'CGM48': 0.0, 'CGM49': 0.0, 'CGM51': 0.0, 'CGM52': 0.0, 'CGM55': 0.091, 'CGM57': 0.0, 'CGM58': 0.162, 'CGM60': 0.111, 'Vanessa': 0.119, 'dt15': 0.0, 'flex1': 0.0}
events min_normed_power min_duration_ms max_duration_ms median_frequency min_frequency max_frequency table Ca1Ripples 1303128 0.23 20.0 200.0 153.75 0.00 238.69 Dorsalripples 31409 4.00 20.0 193.6 155.04 132.65 193.35 Globalripples 1648 4.00 20.0 157.6 155.19 136.30 191.60 Nonoise_Ca1Ripples 55692 4.00 20.0 199.2 158.14 127.60 206.42 Ripples 165926 5.00 20.0 149.6 158.54 130.52 204.14 Ventralripples 53145 4.00 20.0 197.6 154.29 130.80 202.77 _Remove_Ca1Ripples 243 5.00 20.0 117.6 150.42 134.21 179.86
The date and time in each file name equals the session reference time plus the LFP starting time,
so the session clock is the lab's wall clock (the reference is labelled UTC but holds local time).
In the SPW-R tables, start_time is on the session clock and peak is relative to the start of the
file (start_time is LFP start + peak minus 10 to 45 ms, the time from start to peak); the column
descriptions do not say this. SPW-R times used here are LFP start + peak. Durations lie in 20 to 200
ms, inside the Methods' 15 to 250 ms. The Methods set the peak threshold at 5 s.d.; the stored
peak normed power is at least 5 in the tables of 12 subjects but goes down to 4 in nine others
(7 to 29% of their events) and below 1 in CGM30, so those tables were made with a lower threshold.
The analysis keeps events of at least 5 s.d. (Methods) and shows the main results with all stored
events as a variant.
ss = pd.concat([f["sleep"].assign(subject=f["subject"]) for f in kept if f["sleep"] is not None])
print("sleep-state labels:", ss.label.value_counts().to_dict(), " subjects:", sorted(ss.subject.unique()))
# clock check between SPW-Rs and sleep states: SPW-R rate per state, true and with states shifted by 30 min
rows = []
for f in kept:
if f["sleep"] is None or not len(f["sleep"]) or S.ca1_table(f) is None:
continue
e = S.ripple_events(f).t.values
for shift in (0, 1800):
for lab, g in f["sleep"].groupby("label"):
st, sp = g.start.values + shift, g.stop.values + shift
n = sum(((e >= a) & (e < b)).sum() for a, b in zip(st, sp))
rows.append({"subject": f["subject"], "shift_s": shift, "state": lab, "events": n, "s": (sp - st).sum()})
r = pd.DataFrame(rows).groupby(["subject", "shift_s", "state"]).sum()
r["per_min"] = 60 * r.events / r.s
print(r.per_min.unstack("state").round(2).to_string())
sleep-state labels: {'Awake': 11349, 'Non-REM': 10437, 'REM': 3770} subjects: ['CGM1', 'CGM2', 'CGM3', 'bruce', 'dt15', 'flex1']
state Awake Non-REM REM
subject shift_s
CGM1 0 3.10 18.40 0.38
1800 10.33 10.46 10.16
CGM2 0 3.41 19.79 0.54
1800 8.55 9.19 9.50
CGM3 0 11.97 20.58 1.52
1800 13.89 14.71 14.70
dt15 0 5.51 5.93 0.65
1800 8.54 3.96 7.14
flex1 0 6.54 20.96 0.26
1800 9.94 12.55 9.99
Sleep states are present for six subjects. In CGM1, CGM2, CGM3 and flex1, SPW-Rs occur at 18 to 21 per minute in NREM sleep against 3 to 12 when awake with the states as stored, and the contrast shrinks when the states are shifted by 30 min, so the SPW-R tables and the sleep states share a clock there. In dt15 the contrast is small (5.9 against 5.5 per minute) and the check does not establish it; bruce has no SPW-R table.
rows, GLU_FIRST = [], {}
for s, idx in assets.groupby("subject").groups.items():
stored = {(glu_all[i]["t"].tobytes(), glu_all[i]["nA"].tobytes()) for i in idx}
wall = {tuple(sorted(zip(S.clock(files[i]["reference"], glu_all[i]["t"]), glu_all[i]["nA"]))) for i in idx}
refs = {files[i]["reference"] for i in idx}
rows.append({"subject": s, "files": len(idx), "distinct_series_as_stored": len(stored),
"distinct_series_on_wall_clock": len(wall), "session_references": len(refs)})
print(pd.DataFrame(rows).to_string())
for f, g in zip(files, glu_all):
GLU_FIRST.setdefault(f["subject"], g)
rows = []
for s in sorted(GLU_FIRST):
g = GLU_FIRST[s]
ref = next(f["reference"] for f in files if f["subject"] == s)
t = np.sort(g["t"])
runs, swapped, b_stored, b_swapped = S.repaired_glucose(ref, g)
rows.append({"subject": s, "readings": len(g["t"]),
"unit": g["unit"], "nA_min": g["nA"].min(), "nA_median": np.median(g["nA"]), "nA_max": g["nA"].max(),
"share_steps_300s": np.mean(np.abs(np.diff(t) - 300) < 10), "order_breaks_as_stored": int((np.abs(np.diff(g["t"]) - 300) > 10).sum()),
"breaks_sorted_as_stored": b_stored, "breaks_day_month_swapped": b_swapped, "repaired": swapped,
"first_reading": str(S.clock(ref, [runs.t.min()])[0])[:16], "last_reading": str(S.clock(ref, [runs.t.max()])[0])[:16],
"readings_over_30_days_from_median": int((np.abs(runs.t - np.median(runs.t)) > 30 * 86400).sum())})
gq = pd.DataFrame(rows)
print(gq.round(3).to_string())
print("glucose description:", GLU_FIRST["CGM1"]["description"])
subject files distinct_series_as_stored distinct_series_on_wall_clock session_references 0 CGM1 28 2 1 1 1 CGM2 12 2 1 2 2 CGM3 35 2 1 2 3 CGM30 22 1 1 1 4 CGM31 33 1 1 1 5 CGM32 42 1 1 1 6 CGM36 17 1 1 1 7 CGM37 15 1 1 1 8 CGM39 19 1 1 1 9 CGM40 7 1 1 1 10 CGM41 28 1 1 1 11 CGM46 4 1 1 1 12 CGM48 6 1 1 1 13 CGM49 14 1 1 1 14 CGM51 7 1 1 1 15 CGM52 6 1 1 1 16 CGM55 6 1 1 1 17 CGM57 4 1 1 1 18 CGM58 12 1 1 1 19 CGM60 3 1 1 1 20 Vanessa 4 1 1 1 21 bruce 9 1 1 1 22 dt15 8 1 1 1 23 flex1 1 1 1 1 24 ros 3 1 1 1
subject readings unit nA_min nA_median nA_max share_steps_300s order_breaks_as_stored breaks_sorted_as_stored breaks_day_month_swapped repaired first_reading last_reading readings_over_30_days_from_median 0 CGM1 2962 nA 0.00 33.990 106.18 0.998 7 7 17 False 2018-11-29 18:15 2018-12-12 10:56 0 1 CGM2 3865 nA 0.00 29.440 106.55 0.998 7 7 9 False 2018-12-14 15:52 2019-01-06 11:04 0 2 CGM3 3637 nA 0.13 33.330 119.94 0.998 7 7 18 False 2019-03-01 18:38 2019-03-17 09:40 0 3 CGM30 1935 nA -0.29 34.260 57.79 0.995 10 10 7 True 2021-01-24 11:02 2021-02-05 08:53 0 4 CGM31 4020 nA -2.42 33.265 62.37 0.997 13 13 6 True 2021-01-21 17:38 2021-02-12 08:34 0 5 CGM32 2860 nA -1.32 32.810 50.22 0.997 10 10 10 False 2021-01-21 18:17 2021-12-02 23:58 86 6 CGM36 1484 nA -0.72 32.535 49.38 0.997 4 4 4 False 2021-02-22 14:49 2021-03-03 13:58 0 7 CGM37 720 nA -0.63 27.015 53.68 0.999 1 1 1 False 2021-01-02 09:50 2021-01-30 23:04 0 8 CGM39 1203 nA -0.63 32.160 56.33 0.996 5 5 1 True 2021-02-01 15:57 2021-02-08 09:08 0 9 CGM40 1578 nA 24.61 34.560 58.55 0.996 7 7 1 True 2021-02-01 16:33 2021-02-08 09:13 0 10 CGM41 3316 nA -2.35 31.990 64.38 0.995 15 15 3 True 2021-02-04 17:53 2021-03-10 17:07 0 11 CGM46 2020 nA 22.82 34.765 66.34 0.996 8 8 1 True 2021-03-04 15:14 2021-03-13 08:52 0 12 CGM48 3173 nA -1.33 29.410 54.55 0.997 11 11 4 True 2021-03-01 12:01 2021-03-21 13:07 0 13 CGM49 3166 nA -1.23 32.705 61.80 0.999 3 3 3 False 2021-02-13 12:31 2021-02-28 09:16 0 14 CGM51 1679 nA -0.91 29.870 54.38 0.997 5 5 2 True 2021-02-19 17:37 2021-03-10 17:18 0 15 CGM52 1648 nA -0.69 29.545 54.49 0.999 1 1 1 False 2021-02-16 14:43 2021-02-28 09:20 0 16 CGM55 1241 nA -0.99 38.510 48.75 0.998 2 2 2 False 2021-03-13 15:56 2021-03-19 15:34 0 17 CGM57 1795 nA -1.89 33.100 53.37 0.996 7 7 1 True 2021-03-01 14:06 2021-03-10 09:11 0 18 CGM58 1958 nA -0.98 33.185 51.52 0.995 10 10 6 True 2021-03-08 15:45 2021-03-17 09:19 0 19 CGM60 815 nA 1.12 40.970 53.08 0.999 1 1 1 False 2021-03-20 14:37 2021-03-24 15:17 0 20 Vanessa 2020 nA 0.23 0.250 50.56 1.000 1 1 1 False 2021-02-03 11:32 2021-02-26 16:11 0 21 bruce 866 nA 0.29 28.125 61.01 0.995 4 4 2 True 2021-01-10 16:59 2021-01-14 15:47 0 22 dt15 1377 nA 26.65 38.620 52.97 0.999 1 1 1 False 2020-02-14 16:27 2020-02-19 11:24 0 23 flex1 2979 nA 23.29 29.870 66.22 0.997 9 9 4 True 2019-11-25 18:16 2019-12-08 10:52 0 24 ros 1747 nA 23.10 36.880 71.30 0.999 1 1 1 False 2021-03-19 16:32 2021-03-25 18:13 0 glucose description: Raw current from Medtronic iPro2 ISIG tracking.
Each subject's glucose series is stored whole in every one of its files. In CGM1, CGM2 and CGM3
one file stores it in a different order or relative to a different session reference time (CGM2
and CGM3 have two reference times), but on the wall clock every file holds the same readings. Times
are therefore put on one clock per subject, that of its first file (helper clock_offset); the
series of the first file is used. Readings are 5 min apart (300 s; a few at 301 and 302 s). The values are
called nA, and their range (about 0 to 120) is consistent with the sensor current. In twelve
subjects the timestamps are out of order and the series jumps by about a month at each midnight:
dates whose day is 12 or less were stored with day and month swapped (a reading of 4 March 2021 is
stored as 3 April). Swapping them back joins the pieces into runs of consecutive 5-min readings
(fewer breaks, table above) that overlap the electrophysiology. The repair is applied to a subject
only when it reduces the breaks. In CGM32 it does not (10 breaks either way): 86 readings are dated
more than 30 days from the subject's median date (up to 2021-12-02), which looks like 1 to 12
February with day and month swapped, and they are left as stored. No CGM32 recording falls on
those days, and the glucose-only analyses below use runs of at least 24 h, so they do not use them; the effect of the repair on the SPW-R and glucose correlation is
shown below.
SUBJ = {}
rows = []
for s in sorted(GLU_FIRST):
fs = [f for f in kept if f["subject"] == s]
ref = next(f["reference"] for f in files if f["subject"] == s)
ev, iv = S.subject_events(fs, reference=ref)
ev_all, _ = S.subject_events(fs, min_peak_sd=0, reference=ref)
glu, swapped, _, _ = S.repaired_glucose(ref, GLU_FIRST[s])
b = S.binned(glu, ev, iv) if ev is not None else None
b_all = S.binned(glu, ev_all, iv) if ev is not None else None
stored = S.binned(S.glucose_runs(GLU_FIRST[s]["t"], GLU_FIRST[s]["nA"]), ev, iv) if ev is not None else None
if b is not None:
b["clock"] = S.clock(ref, b.t)
SUBJ[s] = {"files": fs, "ref": ref, "ev": ev, "ev_all": ev_all, "iv": iv, "glu": glu, "b": b, "b_all": b_all, "stored": stored, "swapped": swapped}
ok = 0 if b is None else int((np.isfinite(b.rate) & np.isfinite(b.delta)).sum())
rows.append({"subject": s, "ca1_files": len(iv), "ca1_hours": sum(y - x for x, y in iv) / 3600,
"spw_r_stored": 0 if ev is None else len(ev_all), "spw_r_5sd": 0 if ev is None else len(ev), "glucose_hours": len(glu) * 5 / 60,
"hours_with_both": ok * 5 / 60,
"hours_with_both_as_stored": 0 if stored is None else int((np.isfinite(stored.rate) & np.isfinite(stored.delta)).sum()) * 5 / 60})
cov = pd.DataFrame(rows)
print(cov.round(1).to_string())
RATS = [s for s in cov[cov.hours_with_both > 0].subject]
print(len(RATS), "rats with CA1 SPW-Rs and glucose recorded together:", RATS,
f"; {cov.hours_with_both.sum():.1f} h with both ({cov.hours_with_both_as_stored.sum():.1f} h with the timestamps as stored)")
subject ca1_files ca1_hours spw_r_stored spw_r_5sd glucose_hours hours_with_both hours_with_both_as_stored 0 CGM1 18 257.0 165926 165926 246.8 216.8 216.8 1 CGM2 10 58.4 32809 32809 322.1 56.3 56.3 2 CGM3 27 350.5 307187 286372 303.1 284.2 284.2 3 CGM30 18 174.2 120120 85377 161.2 134.6 97.3 4 CGM31 28 330.5 243127 197710 335.0 265.8 132.9 5 CGM32 33 297.3 154965 154717 238.3 198.2 198.2 6 CGM36 0 0.0 0 0 123.7 0.0 0.0 7 CGM37 8 15.0 9562 7560 60.0 0.0 0.0 8 CGM39 14 106.7 35851 30127 100.2 51.8 17.6 9 CGM40 7 75.1 52517 52517 131.5 72.8 23.8 10 CGM41 27 297.2 162347 162142 276.3 210.6 12.7 11 CGM46 4 49.9 38632 38632 168.3 49.4 0.0 12 CGM48 3 3.5 1504 1504 264.4 2.2 0.0 13 CGM49 7 28.5 20855 20855 263.8 6.0 6.0 14 CGM51 7 25.0 10422 10422 139.9 17.3 17.3 15 CGM52 6 43.1 26556 26556 137.3 17.9 17.9 16 CGM55 4 25.5 19815 18010 103.4 16.0 16.0 17 CGM57 4 18.9 14234 14234 149.6 18.0 0.0 18 CGM58 9 83.8 35091 29418 163.2 73.8 17.5 19 CGM60 2 11.2 7875 6998 67.9 0.0 0.0 20 Vanessa 4 42.9 20343 17929 168.3 24.8 24.8 21 bruce 0 0.0 0 0 72.2 0.0 0.0 22 dt15 7 28.9 9740 9740 114.8 0.0 0.0 23 flex1 1 16.7 10985 10985 248.2 0.0 0.0 24 ros 0 0.0 0 0 145.6 0.0 0.0 18 rats with CA1 SPW-Rs and glucose recorded together: ['CGM1', 'CGM2', 'CGM3', 'CGM30', 'CGM31', 'CGM32', 'CGM39', 'CGM40', 'CGM41', 'CGM46', 'CGM48', 'CGM49', 'CGM51', 'CGM52', 'CGM55', 'CGM57', 'CGM58', 'Vanessa'] ; 1716.6 h with both (1139.2 h with the timestamps as stored)
Eighteen subjects have CA1 SPW-Rs and glucose recorded at the same time. In CGM36, bruce and ros no file has a SPW-R table; in CGM37, CGM60, dt15 and flex1 the glucose readings and the electrophysiology cover different days. Glucose and electrophysiology share no event that could show their clocks agree independently of the claims; the repair of the dates, which brings the two into overlap, is consistent with the file names but the alignment within a day is taken as stored.
fig2-ripple-rate-predicts-glucose¶
SPW-R rate is negatively correlated with future changes in interstitial glucose, with a mean peak correlation of R = -0.21 (+- 0.08) at a lag of 10.7 (+- 4.5) min, over 30 rats.
Readings, as the claims file lists them: delta current from the raw readings (a) or after a high-pass filter (b, a centred 24-h running mean subtracted within each run, a standard choice for the cutoff the paper does not give); Pearson or Spearman correlation; one correlogram per rat (a) or one per session averaged per rat (b; a session is a file, used when it holds at least 6 h of bins, a standard choice); the lag searched over 0 to +100 min (a) or -100 to +100 min (b); all rats with SPW-Rs and glucose together (a) or the ripple and glucose cohort (b; only CGM1 is labelled as such in the files). Reading (a) leaves out CNO days and the stimulation periods of the optogenetic rats; the files label neither, so every day of every rat is used. The paper does not say how the SPW-R count of a 5-min bin is aligned with the glucose change; the main analysis pairs the count and the change over the same interval (an assumption), and the two alignments one reading earlier or later are shown as variants.
def highpass(b, hours=24):
out = b.copy()
w = int(hours * 12)
out["nA"] = b.nA - b.groupby("run").nA.transform(lambda x: x.rolling(w, center=True, min_periods=1).mean())
same = out.run.eq(out.run.shift(1))
out["delta"] = np.where(same, out.nA - out.nA.shift(1), np.nan)
return out
def per_session(b, iv, min_bins=72):
"""Rows of b inside each file, one table per file with at least min_bins bins of both."""
out = []
for x, y in iv:
sel = b[(b.t - S.STEP >= x) & (b.t <= y)]
if (np.isfinite(sel.rate) & np.isfinite(sel.delta)).sum() >= min_bins:
out.append(sel)
return out
LAG_MIN = S.LAGS * 5
# the other alignments of the SPW-R interval and the glucose change, one reading apart, are the same
# pairs at a lag one reading shorter (SPW-Rs over the 5 min after each reading) or longer (delta to the next reading)
ALIGN = {"same interval (main)": 0, "SPW-Rs after the reading": -1, "delta to the next reading": 1}
CC = {}
for s in RATS:
d = SUBJ[s]
for dname, bb in (("raw", d["b"]), ("highpass", highpass(d["b"]))):
for m in ("pearson", "spearman"):
CC[(s, dname, m, "rat")] = S.xcorr(bb, method=m)
sess = per_session(bb, d["iv"])
CC[(s, dname, m, "session")] = (np.nanmean([S.xcorr(x, method=m) for x in sess], axis=0)
if sess else np.full(len(S.LAGS), np.nan))
def peak(cc, lag_range):
sel = LAG_MIN >= 0 if lag_range == "0..+100" else np.ones(len(LAG_MIN), bool)
if np.all(np.isnan(cc[sel])):
return np.nan, np.nan
i = np.nanargmin(np.where(sel, cc, np.nan))
return cc[i], LAG_MIN[i]
rows = []
for (s, dname, m, unit), cc in CC.items():
for lr in ("0..+100", "-100..+100"):
R, lag = peak(cc, lr)
rows.append({"rat": s, "delta": dname, "corr": m, "unit": unit, "lags": lr, "R": R, "lag_min": lag,
"R_at_10": cc[LAG_MIN == 10][0]})
PEAKS = pd.DataFrame(rows)
def summarise(df):
return pd.Series({"rats": df.R.notna().sum(), "mean_R": df.R.mean(), "sd_R": df.R.std(),
"mean_lag": df.lag_min.mean(), "sd_lag": df.lag_min.std(),
"rats_R_negative": (df.R < 0).sum(), "R_ok": abs(df.R.mean() + 0.21) <= 0.08,
"lag_ok": abs(df.lag_min.mean() - 10.7) <= 4.5})
tab_a = PEAKS.groupby(["delta", "corr", "unit", "lags"]).apply(summarise)
tab_a.insert(0, "rats_set", "a: all")
tab_b = PEAKS[PEAKS.rat == "CGM1"].groupby(["delta", "corr", "unit", "lags"]).apply(summarise)
tab_b.insert(0, "rats_set", "b: CGM1")
T1 = pd.concat([tab_a, tab_b]).reset_index()
T1["both_hold_and_most_negative"] = T1.R_ok & T1.lag_ok & (T1.rats_R_negative > T1.rats / 2)
print(T1.round(3).to_string())
delta corr unit lags rats_set rats mean_R sd_R mean_lag sd_lag rats_R_negative R_ok lag_ok both_hold_and_most_negative 0 highpass pearson rat -100..+100 a: all 18 -0.180 0.143 14.444 41.969 18 True True True 1 highpass pearson rat 0..+100 a: all 18 -0.162 0.117 33.056 30.541 18 True False False 2 highpass pearson session -100..+100 a: all 16 -0.177 0.109 26.875 34.199 16 True False False 3 highpass pearson session 0..+100 a: all 16 -0.171 0.113 29.062 31.582 16 True False False 4 highpass spearman rat -100..+100 a: all 18 -0.216 0.145 15.833 29.816 18 True False False 5 highpass spearman rat 0..+100 a: all 18 -0.212 0.145 25.556 29.200 18 True False False 6 highpass spearman session -100..+100 a: all 16 -0.184 0.122 8.750 37.572 16 True True True 7 highpass spearman session 0..+100 a: all 16 -0.177 0.131 20.000 25.690 15 True False False 8 raw pearson rat -100..+100 a: all 18 -0.176 0.145 14.167 41.666 18 True True True 9 raw pearson rat 0..+100 a: all 18 -0.159 0.121 32.778 30.303 18 True False False 10 raw pearson session -100..+100 a: all 16 -0.178 0.110 27.500 31.091 16 True False False 11 raw pearson session 0..+100 a: all 16 -0.176 0.112 28.438 29.929 16 True False False 12 raw spearman rat -100..+100 a: all 18 -0.212 0.145 15.833 29.816 18 True False False 13 raw spearman rat 0..+100 a: all 18 -0.209 0.144 25.556 29.200 18 True False False 14 raw spearman session -100..+100 a: all 16 -0.199 0.118 24.062 30.617 16 True False False 15 raw spearman session 0..+100 a: all 16 -0.197 0.120 24.375 30.325 16 True False False 16 highpass pearson rat -100..+100 b: CGM1 1 -0.103 NaN 5.000 NaN 1 False False False 17 highpass pearson rat 0..+100 b: CGM1 1 -0.103 NaN 5.000 NaN 1 False False False 18 highpass pearson session -100..+100 b: CGM1 1 -0.197 NaN 5.000 NaN 1 True False False 19 highpass pearson session 0..+100 b: CGM1 1 -0.197 NaN 5.000 NaN 1 True False False 20 highpass spearman rat -100..+100 b: CGM1 1 -0.196 NaN 5.000 NaN 1 True False False 21 highpass spearman rat 0..+100 b: CGM1 1 -0.196 NaN 5.000 NaN 1 True False False 22 highpass spearman session -100..+100 b: CGM1 1 -0.223 NaN 5.000 NaN 1 True False False 23 highpass spearman session 0..+100 b: CGM1 1 -0.223 NaN 5.000 NaN 1 True False False 24 raw pearson rat -100..+100 b: CGM1 1 -0.101 NaN 5.000 NaN 1 False False False 25 raw pearson rat 0..+100 b: CGM1 1 -0.101 NaN 5.000 NaN 1 False False False 26 raw pearson session -100..+100 b: CGM1 1 -0.197 NaN 5.000 NaN 1 True False False 27 raw pearson session 0..+100 b: CGM1 1 -0.197 NaN 5.000 NaN 1 True False False 28 raw spearman rat -100..+100 b: CGM1 1 -0.194 NaN 5.000 NaN 1 True False False 29 raw spearman rat 0..+100 b: CGM1 1 -0.194 NaN 5.000 NaN 1 True False False 30 raw spearman session -100..+100 b: CGM1 1 -0.224 NaN 5.000 NaN 1 True False False 31 raw spearman session 0..+100 b: CGM1 1 -0.224 NaN 5.000 NaN 1 True False False
main = PEAKS[(PEAKS["delta"] == "raw") & (PEAKS["corr"] == "spearman") & (PEAKS.unit == "rat") & (PEAKS.lags == "-100..+100")]
print("Per rat, raw delta, Spearman, one correlogram per rat, lags -100..+100 min:")
print(main.merge(cov[["subject", "hours_with_both"]], left_on="rat", right_on="subject")[["rat", "hours_with_both", "R", "lag_min", "R_at_10"]].round(3).to_string())
for minh in (6, 12, 24):
keep = cov[cov.hours_with_both >= minh].subject
x = PEAKS[PEAKS.rat.isin(keep)].groupby(["delta", "corr", "unit", "lags"]).apply(summarise)
print(f"rats with at least {minh} h of both: combinations where both values hold: {int((x.R_ok & x.lag_ok).sum())} of {len(x)};"
f" raw/spearman/rat/-100..+100: R {x.loc[('raw', 'spearman', 'rat', '-100..+100'), 'mean_R']:.3f},"
f" lag {x.loc[('raw', 'spearman', 'rat', '-100..+100'), 'mean_lag']:.1f} min, rats {int(x.loc[('raw', 'spearman', 'rat', '-100..+100'), 'rats'])}")
keep6 = cov[cov.hours_with_both >= 6].subject
x6 = PEAKS[PEAKS.rat.isin(keep6)].groupby(["delta", "corr", "unit", "lags"]).apply(summarise)
print("with at least 6 h of both, the combinations where both values hold:", list(x6[x6.R_ok & x6.lag_ok].index))
print("median of the per-rat lags (min) and rats with the peak at +5 to +15 min, all rats:")
print(PEAKS.groupby(["delta", "corr", "unit", "lags"]).lag_min.agg(
median_lag="median", rats="count", rats_5_to_15=lambda x: int(x.between(5, 15).sum())).to_string())
passing = T1[(T1.rats_set == "a: all") & T1.both_hold_and_most_negative]
for _, r_ in passing.iterrows():
sel = PEAKS[(PEAKS["delta"] == r_["delta"]) & (PEAKS["corr"] == r_["corr"]) & (PEAKS.unit == r_.unit) & (PEAKS.lags == r_.lags)]
print(f"per-rat lags (min), {r_['delta']}/{r_['corr']}/{r_.unit}/{r_.lags}:", dict(zip(sel.rat, sel.lag_min)))
# variant: other alignments of the SPW-R interval and the glucose change, one reading apart
for name, how in (("SPW-R rate over the 5 min after each reading", "rate"), ("delta to the next reading", "delta")):
al = []
for s in RATS:
for dname in ("raw", "highpass"):
bb = SUBJ[s]["b"] if dname == "raw" else highpass(SUBJ[s]["b"])
bb = bb.copy()
bb[how] = bb.groupby("run")[how].shift(-1)
for m in ("pearson", "spearman"):
cc = S.xcorr(bb, method=m)
for lr in ("0..+100", "-100..+100"):
R, lag = peak(cc, lr)
al.append({"rat": s, "delta": dname, "corr": m, "lags": lr, "R": R, "lag_min": lag})
t_al = pd.DataFrame(al).groupby(["delta", "corr", "lags"]).apply(summarise)
print(f"Alignment variant, {name} (one correlogram per rat, all rats):")
print(t_al.round(3).to_string())
# variant: all stored events, without the Methods' 5 s.d. threshold
va = []
for s in RATS:
for m in ("pearson", "spearman"):
cc = S.xcorr(SUBJ[s]["b_all"], method=m)
for lr in ("0..+100", "-100..+100"):
R, lag = peak(cc, lr)
va.append({"rat": s, "corr": m, "lags": lr, "R": R, "lag_min": lag})
print("Variant with all stored events (raw delta, one correlogram per rat, all rats):")
print(pd.DataFrame(va).groupby(["corr", "lags"]).apply(summarise).round(3).to_string())
# the effect of the date repair, in the rats it was applied to
rep = [s for s in RATS if SUBJ[s]["swapped"]]
cmp_rows = []
for s in rep:
st = SUBJ[s]["stored"]
n_st = int((np.isfinite(st.rate) & np.isfinite(st.delta)).sum())
cmp_rows.append({"rat": s, "bins_repaired": int((np.isfinite(SUBJ[s]["b"].rate) & np.isfinite(SUBJ[s]["b"].delta)).sum()),
"R10_repaired": CC[(s, "raw", "spearman", "rat")][LAG_MIN == 10][0],
"bins_as_stored": n_st, "R10_as_stored": S.xcorr(st, method="spearman")[LAG_MIN == 10][0] if n_st else np.nan})
print("Spearman R at +10 min, with the repaired and the stored glucose timestamps:")
print(pd.DataFrame(cmp_rows).round(3).to_string())
Per rat, raw delta, Spearman, one correlogram per rat, lags -100..+100 min:
rat hours_with_both R lag_min R_at_10
0 CGM1 216.750 -0.194 5.0 -0.137
1 CGM2 56.333 -0.065 5.0 -0.014
2 CGM3 284.167 -0.044 65.0 -0.033
3 CGM30 134.583 -0.116 5.0 -0.097
4 CGM31 265.750 -0.162 5.0 -0.120
5 CGM32 198.167 -0.115 10.0 -0.115
6 CGM39 51.833 -0.189 5.0 -0.123
7 CGM40 72.833 -0.182 10.0 -0.182
8 CGM41 210.583 -0.124 60.0 -0.086
9 CGM46 49.417 -0.479 0.0 -0.327
10 CGM48 2.250 -0.296 -40.0 -0.094
11 CGM49 6.000 -0.593 15.0 -0.535
12 CGM51 17.333 -0.128 -30.0 -0.117
13 CGM52 17.917 -0.407 10.0 -0.407
14 CGM55 16.000 -0.206 60.0 0.041
15 CGM57 18.000 -0.204 0.0 -0.062
16 CGM58 73.833 -0.150 55.0 -0.021
17 Vanessa 24.833 -0.163 45.0 0.118
rats with at least 6 h of both: combinations where both values hold: 1 of 16; raw/spearman/rat/-100..+100: R -0.207, lag 19.1 min, rats 17
rats with at least 12 h of both: combinations where both values hold: 1 of 16; raw/spearman/rat/-100..+100: R -0.183, lag 19.4 min, rats 16
rats with at least 24 h of both: combinations where both values hold: 1 of 16; raw/spearman/rat/-100..+100: R -0.165, lag 22.5 min, rats 12
with at least 6 h of both, the combinations where both values hold: [('highpass', 'spearman', 'session', '-100..+100')]
median of the per-rat lags (min) and rats with the peak at +5 to +15 min, all rats:
median_lag rats rats_5_to_15
delta corr unit lags
highpass pearson rat -100..+100 15.0 18 8
0..+100 15.0 18 9
session -100..+100 10.0 16 7
0..+100 10.0 16 7
spearman rat -100..+100 7.5 18 9
0..+100 10.0 18 10
session -100..+100 5.0 16 7
0..+100 10.0 16 8
raw pearson rat -100..+100 15.0 18 8
0..+100 15.0 18 9
session -100..+100 10.0 16 8
0..+100 10.0 16 8
spearman rat -100..+100 7.5 18 9
0..+100 10.0 18 10
session -100..+100 10.0 16 7
0..+100 10.0 16 7
per-rat lags (min), highpass/pearson/rat/-100..+100: {'CGM1': 5.0, 'CGM2': -20.0, 'CGM3': 20.0, 'CGM30': 80.0, 'CGM31': 5.0, 'CGM32': 15.0, 'CGM39': 5.0, 'CGM40': 15.0, 'CGM41': 60.0, 'CGM46': 0.0, 'CGM48': -40.0, 'CGM49': 15.0, 'CGM51': 10.0, 'CGM52': 15.0, 'CGM55': 60.0, 'CGM57': 60.0, 'CGM58': 55.0, 'Vanessa': -100.0}
per-rat lags (min), highpass/spearman/session/-100..+100: {'CGM1': 5.0, 'CGM2': 60.0, 'CGM3': 5.0, 'CGM30': 75.0, 'CGM31': 5.0, 'CGM32': 15.0, 'CGM39': 0.0, 'CGM40': 10.0, 'CGM41': 10.0, 'CGM46': 0.0, 'CGM48': nan, 'CGM49': nan, 'CGM51': -30.0, 'CGM52': 15.0, 'CGM55': 65.0, 'CGM57': 0.0, 'CGM58': -10.0, 'Vanessa': -85.0}
per-rat lags (min), raw/pearson/rat/-100..+100: {'CGM1': 5.0, 'CGM2': -20.0, 'CGM3': 20.0, 'CGM30': 80.0, 'CGM31': 5.0, 'CGM32': 15.0, 'CGM39': 5.0, 'CGM40': 15.0, 'CGM41': 60.0, 'CGM46': 0.0, 'CGM48': -40.0, 'CGM49': 15.0, 'CGM51': 10.0, 'CGM52': 15.0, 'CGM55': 60.0, 'CGM57': 55.0, 'CGM58': 55.0, 'Vanessa': -100.0}
Alignment variant, SPW-R rate over the 5 min after each reading (one correlogram per rat, all rats):
rats mean_R sd_R mean_lag sd_lag rats_R_negative R_ok lag_ok
delta corr lags
highpass pearson -100..+100 18 -0.179 0.143 19.444 41.969 18 True False
0..+100 18 -0.160 0.117 37.778 29.913 18 True False
spearman -100..+100 18 -0.215 0.146 20.833 29.816 18 True False
0..+100 18 -0.211 0.145 30.556 29.200 18 True False
raw pearson -100..+100 18 -0.175 0.145 19.167 41.666 18 True False
0..+100 18 -0.157 0.120 37.500 29.667 18 True False
spearman -100..+100 18 -0.211 0.145 20.833 29.816 18 True False
0..+100 18 -0.208 0.144 30.556 29.200 18 True False
Alignment variant, delta to the next reading (one correlogram per rat, all rats):
rats mean_R sd_R mean_lag sd_lag rats_R_negative R_ok lag_ok
delta corr lags
highpass pearson -100..+100 18 -0.177 0.144 9.722 41.177 18 True True
0..+100 18 -0.160 0.115 28.333 30.244 18 True False
spearman -100..+100 18 -0.216 0.145 10.833 29.816 18 True True
0..+100 18 -0.208 0.142 21.111 28.726 18 True False
raw pearson -100..+100 18 -0.173 0.146 9.444 40.869 18 True True
0..+100 18 -0.157 0.118 28.056 30.007 18 True False
spearman -100..+100 18 -0.212 0.145 10.833 29.816 18 True True
0..+100 18 -0.204 0.141 21.111 28.726 18 True False
Variant with all stored events (raw delta, one correlogram per rat, all rats):
rats mean_R sd_R mean_lag sd_lag rats_R_negative R_ok lag_ok
corr lags
pearson -100..+100 18 -0.175 0.146 14.167 41.666 18 True True
0..+100 18 -0.158 0.121 32.778 30.303 18 True False
spearman -100..+100 18 -0.212 0.145 16.111 30.222 18 True False
0..+100 18 -0.209 0.144 25.833 29.518 18 True False
Spearman R at +10 min, with the repaired and the stored glucose timestamps:
rat bins_repaired R10_repaired bins_as_stored R10_as_stored
0 CGM30 1615 -0.097 1168 -0.081
1 CGM31 3189 -0.120 1595 -0.107
2 CGM39 622 -0.123 211 0.012
3 CGM40 874 -0.182 285 -0.187
4 CGM41 2527 -0.086 152 -0.130
5 CGM46 593 -0.327 0 NaN
6 CGM48 27 -0.094 0 NaN
7 CGM51 208 -0.117 208 -0.117
8 CGM57 216 -0.062 0 NaN
9 CGM58 886 -0.021 210 -0.117
fig, ax = plt.subplots(1, 2, figsize=(11, 3.5))
for s in RATS:
ax[0].plot(LAG_MIN, CC[(s, "raw", "spearman", "rat")], color="0.7", lw=0.8)
med = np.nanmedian([CC[(s, "raw", "spearman", "rat")] for s in RATS], axis=0)
ax[0].plot(LAG_MIN, med, color="C3", lw=2)
ax[0].axvline(0, color="k", lw=0.5); ax[0].set_xlabel("offset (min; + = glucose after SPW-Rs)"); ax[0].set_ylabel("Spearman R")
ax[0].set_title("SPW-R rate vs delta current, per rat (grey) and median (red)")
allcc = np.array([CC[(s, "raw", "spearman", "rat")] for s in RATS])
ax[1].imshow(allcc, aspect="auto", cmap="RdBu_r", vmin=-0.4, vmax=0.4, extent=[-102.5, 102.5, len(RATS) - 0.5, -0.5])
ax[1].set_yticks(range(len(RATS))); ax[1].set_yticklabels(RATS, fontsize=7); ax[1].set_xlabel("offset (min)")
plt.tight_layout(); plt.show()
print(f"plotted R range {np.nanmin(allcc):.3f} to {np.nanmax(allcc):.3f}; colour scale -0.4 to 0.4, "
f"{int((np.abs(allcc) > 0.4).sum())} values outside it; median correlogram minimum {np.nanmin(med):.3f} at {LAG_MIN[np.nanargmin(med)]} min")
plotted R range -0.593 to 0.509; colour scale -0.4 to 0.4, 15 values outside it; median correlogram minimum -0.108 at 5 min
Result: reproduced, with 18 rats where the paper has 30, under 3 of the 16 combinations of readings, resting on an assumption about the alignment. With Pearson correlation, one correlogram per rat and lags searched over -100 to +100 min, the mean peak correlation is -0.176 (raw delta current) or -0.180 (high-pass) at a mean lag of 14.2 or 14.4 min, inside both tolerances; with Spearman correlation, one correlogram per session and the high-pass delta, -0.184 at 8.8 min. The peak is negative in every rat. Mean R lies within the tolerance under all 16 combinations (-0.159 to -0.216); the mean lag does not under the other 13 (15.8 to 33.1 min), and it is unstable: the lags of single rats spread with an s.d. of 26 to 42 min, and leaving out the rats with less than 6 h of data leaves 1 combination of 16 that matches. With at least 6 h of both per rat, only the high-pass, Spearman, per-session, -100 to +100 min combination matches. The median of the per-rat lags is 5 to 15 min in every combination, but only 7 to 10 rats of 16 or 18 peak between +5 and +15 min; in the raw, Pearson, per-rat combination single-rat lags run from -100 to +80 min, and the mean lies in tolerance because peaks at -100 to -20 min offset peaks at +55 to +80 min. "R is negative in most rats" holds by construction, as R is the minimum over 41 offsets. The median correlogram has its minimum, -0.108, at +5 min. CGM1, the one rat labelled as the ripple and glucose cohort, peaks at +5 min with R -0.101 (Pearson) or -0.194 (Spearman). The outcome rests on pairing the SPW-R count and the glucose change of the same 5-min interval, an assumption: counting SPW-Rs over the 5 min after each reading moves every lag by +5 min, and then none of the 8 per-rat combinations matches (mean lags 19.2 to 37.8 min), while taking the change to the next reading moves them by -5 min and 4 of 8 match (mean lags 9.4 to 10.8 min).
fig2ef-correlation-beyond-shuffle¶
The correlation between SPW-R rate and glucose change at a 10-min lag is more negative than the correlations obtained when SPW-R rate is circularly permuted.
Per rat, the correlation at +10 min against 100 circular shuffles of SPW-R rate within each run (the claims file's shuffle), with the two-sample Kolmogorov-Smirnov test the claims file names. Both correlation types are shown.
The main analysis pairs the SPW-R count and the glucose change of the same 5-min interval, as in fig2-ripple-rate-predicts-glucose (an assumption). The two alignments one reading apart are shown as variants: counting SPW-Rs over the 5 min after each reading, or taking the change to the next reading, which is the same as moving the offset by one reading.
rng = np.random.default_rng(20210811)
lag2 = int(np.where(S.LAGS == 2)[0][0])
obs_al, shuf_al = {}, {}
obs, shuf = {"pearson": [], "spearman": []}, {"pearson": [], "spearman": []}
for s in RATS:
b = SUBJ[s]["b"]
for m in obs:
obs[m].append(CC[(s, "raw", m, "rat")][lag2])
for name, k in ALIGN.items():
obs_al.setdefault((name, m), []).append(S._corr(*S.lagged_pairs(b, 2 + k), m))
for i in range(100):
sb = S.circular_shuffle(b, rng)
for name, k in ALIGN.items():
pairs = S.lagged_pairs(sb, 2 + k)
for m in shuf:
v = S._corr(*pairs, m)
shuf_al.setdefault((name, m), []).append(v)
if k == 0:
shuf[m].append(v)
for m in obs:
o, sh = np.array(obs[m]), np.array(shuf[m])
ks = stats.ks_2samp(o, sh[np.isfinite(sh)])
print(f"{m}: observed median {np.nanmedian(o):.3f} ({len(o)} rats, {np.sum(o < 0)} negative), shuffled median {np.nanmedian(sh):.4f} "
f"(n = {np.isfinite(sh).sum()}, 2.5-97.5%: {np.nanpercentile(sh, 2.5):.3f} to {np.nanpercentile(sh, 97.5):.3f}); KS D = {ks.statistic:.3f}, p = {ks.pvalue:.2e}")
per = np.array(sh).reshape(len(RATS), 100)
print(" rats whose observed value is below all of their own 100 shuffles:", int(np.sum(o < np.nanmin(per, axis=1))))
print("Alignment variants (each with its own 100 shuffles per rat), +10 min:")
for (name, m), o in obs_al.items():
o, sh = np.array(o), np.array(shuf_al[(name, m)])
ks = stats.ks_2samp(o[np.isfinite(o)], sh[np.isfinite(sh)])
per = sh.reshape(len(RATS), 100)
print(f" {name}, {m}: observed median {np.nanmedian(o):.3f} ({np.sum(o < 0)} of {np.isfinite(o).sum()} negative), shuffled median {np.nanmedian(sh):.4f}; "
f"KS p = {ks.pvalue:.2e}; rats below all their shuffles {int(np.sum(o < np.nanmin(per, axis=1)))}")
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist(np.array(shuf["spearman"]), bins=40, density=True, color="k", alpha=0.5, label="shuffled")
ax.hist(np.array(obs["spearman"]), bins=10, density=True, color="C3", alpha=0.7, label="observed")
ax.set_xlabel("Spearman R at +10 min"); ax.legend(); plt.show()
pearson: observed median -0.076 (18 rats, 15 negative), shuffled median 0.0019 (n = 1799, 2.5-97.5%: -0.228 to 0.203); KS D = 0.494, p = 1.51e-04 rats whose observed value is below all of their own 100 shuffles: 5 spearman: observed median -0.106 (18 rats, 16 negative), shuffled median 0.0046 (n = 1799, 2.5-97.5%: -0.221 to 0.202); KS D = 0.575, p = 4.03e-06 rats whose observed value is below all of their own 100 shuffles: 7 Alignment variants (each with its own 100 shuffles per rat), +10 min: same interval (main), pearson: observed median -0.076 (15 of 18 negative), shuffled median 0.0019; KS p = 1.51e-04; rats below all their shuffles 5 SPW-Rs after the reading, pearson: observed median -0.071 (15 of 18 negative), shuffled median 0.0011; KS p = 7.58e-04; rats below all their shuffles 3 delta to the next reading, pearson: observed median -0.045 (13 of 18 negative), shuffled median 0.0039; KS p = 6.02e-03; rats below all their shuffles 4 same interval (main), spearman: observed median -0.106 (16 of 18 negative), shuffled median 0.0046; KS p = 4.03e-06; rats below all their shuffles 7 SPW-Rs after the reading, spearman: observed median -0.108 (16 of 18 negative), shuffled median 0.0049; KS p = 2.49e-07; rats below all their shuffles 7 delta to the next reading, spearman: observed median -0.047 (13 of 18 negative), shuffled median 0.0038; KS p = 2.52e-03; rats below all their shuffles 3
Result: reproduced. At +10 min the observed per-rat correlations have a median of -0.106 (Spearman; 16 of 18 rats negative) against 0.005 for 1,799 circular shuffles, KS p = 4.0e-06; with Pearson, -0.076 against 0.002, KS p = 1.5e-04. Seven rats (Spearman) lie below all 100 of their own shuffles.
The result does not depend on the alignment. Counting SPW-Rs over the 5 min after each reading, the Spearman median is -0.108 (KS p = 2.5e-07; Pearson p = 7.6e-04); taking the change to the next reading, -0.047 (KS p = 2.5e-03; Pearson p = 6.0e-03), each against its own shuffles.
fig2g-clustered-ripples¶
Isolated SPW-Rs have little or no effect on glucose, whereas clustered SPW-Rs (more than 30 per minute) are followed by a several-fold larger decrease.
As the Fig. 2g legend: the mean delta current 10 min after each 5-min epoch, grouped by the number of SPW-Rs in the epoch, pooled over rats. Clustered is more than 150 SPW-Rs per epoch (a) or more than 30 (b); isolated is 1 to 5 per epoch (claims file).
The main analysis pairs the SPW-R count and the glucose change of the same 5-min interval, as in fig2-ripple-rate-predicts-glucose (an assumption). The two alignments one reading apart are shown as variants: counting SPW-Rs over the 5 min after each reading, or taking the change to the next reading, which is the same as moving the offset by one reading.
ep = []
for s in RATS:
b = SUBJ[s]["b"].copy()
b["delta_10"] = b.groupby("run").delta.shift(-2)
b["rat"] = s
ep.append(b[np.isfinite(b.rate) & np.isfinite(b.delta_10)])
ep = pd.concat(ep)
print(f"{len(ep)} epochs from {ep.rat.nunique()} rats; SPW-Rs per epoch: median {ep['count'].median():.0f}, max {ep['count'].max()}")
iso = ep[(ep["count"] >= 1) & (ep["count"] <= 5)]
none = ep[ep["count"] == 0]
for name, thr in (("a: > 150 per epoch", 150), ("b: > 30 per epoch", 30)):
cl = ep[ep["count"] > thr]
t = stats.ttest_ind(cl.delta_10, iso.delta_10, equal_var=False)
print(f"{name}: clustered n = {len(cl)} ({cl.rat.nunique()} rats), mean delta at +10 min {cl.delta_10.mean():.4f} nA; "
f"isolated n = {len(iso)}, mean {iso.delta_10.mean():.4f} nA; ratio {cl.delta_10.mean() / iso.delta_10.mean():.2f}; "
f"Welch t = {t.statistic:.2f}, p = {t.pvalue:.2g}")
print(f"epochs with no SPW-R: n = {len(none)}, mean {none.delta_10.mean():.4f} nA; all epochs mean {ep.delta_10.mean():.4f} nA")
bins = [-1, 0, 5, 15, 30, 60, 90, 120, 150, 200, 400]
ep["group"] = pd.cut(ep["count"], bins)
curve = ep.groupby("group", observed=True).delta_10.agg(["size", "mean", "sem"])
print(curve.round(4).to_string())
per_rat = ep.assign(cls=np.select([ep["count"].between(1, 5), ep["count"] > 150, ep["count"] > 30], ["isolated", "clustered_a", "clustered_b_only"], "other"))
print(per_rat.groupby(["rat", "cls"]).delta_10.mean().unstack().round(3).to_string())
print("Alignment variants (delta current 10 min after the epoch under each alignment):")
for name, k in ALIGN.items():
ev_ = []
for s in RATS:
b = SUBJ[s]["b"].copy()
b["d10"] = b.groupby("run").delta.shift(-(2 + k))
ev_.append(b[np.isfinite(b.rate) & np.isfinite(b.d10)])
ev_ = pd.concat(ev_)
iso_ = ev_[ev_["count"].between(1, 5)].d10
for thr in (150, 30):
cl_ = ev_[ev_["count"] > thr].d10
print(f" {name}, > {thr} per epoch: clustered mean {cl_.mean():.4f} nA (n = {len(cl_)}), isolated {iso_.mean():.4f} nA (n = {len(iso_)}); "
f"ratio {cl_.mean() / iso_.mean():.2f}; Welch p = {stats.ttest_ind(cl_, iso_, equal_var=False).pvalue:.2g}")
cur = ev_.groupby(pd.cut(ev_["count"], bins), observed=True).d10.mean()
print(f" curve by count: {np.round(cur.values, 3)}")
fig, ax = plt.subplots(figsize=(5, 3))
x = np.arange(len(curve)); ax.errorbar(x, curve["mean"], curve["sem"], marker="o")
ax.set_xticks(x); ax.set_xticklabels([str(i) for i in curve.index], rotation=45, fontsize=7)
ax.axhline(0, color="k", lw=0.5); ax.set_xlabel("SPW-Rs per 5-min epoch"); ax.set_ylabel("delta current at +10 min (nA)"); plt.show()
20558 epochs from 18 rats; SPW-Rs per epoch: median 45, max 611
a: > 150 per epoch: clustered n = 585 (13 rats), mean delta at +10 min -0.1910 nA; isolated n = 1718, mean 0.0496 nA; ratio -3.85; Welch t = -1.89, p = 0.059
b: > 30 per epoch: clustered n = 12595 (18 rats), mean delta at +10 min -0.0316 nA; isolated n = 1718, mean 0.0496 nA; ratio -0.64; Welch t = -1.43, p = 0.15
epochs with no SPW-R: n = 1128, mean 0.0595 nA; all epochs mean -0.0164 nA
size mean sem
group
(-1, 0] 1128 0.0595 0.0416
(0, 5] 1718 0.0496 0.0514
(5, 15] 2294 0.0210 0.0415
(15, 30] 2823 -0.0492 0.0580
(30, 60] 4814 -0.0013 0.0309
(60, 90] 3843 -0.0446 0.0371
(90, 120] 2350 -0.0994 0.0626
(120, 150] 1003 0.1248 0.1590
(150, 200] 456 -0.2183 0.1483
(200, 400] 127 -0.1019 0.0649
cls clustered_a clustered_b_only isolated other
rat
CGM1 -0.235 -0.060 0.213 0.080
CGM2 0.133 0.016 0.129 -0.083
CGM3 -0.212 0.092 -0.278 -0.098
CGM30 NaN -0.054 0.067 0.016
CGM31 -0.068 -0.079 0.045 0.024
CGM32 -0.058 -0.044 0.059 0.044
CGM39 NaN -0.065 0.184 -0.045
CGM40 0.198 -0.035 0.043 0.061
CGM41 -0.375 -0.017 0.089 -0.083
CGM46 -0.466 -0.094 0.072 0.286
CGM48 NaN -0.002 0.010 0.128
CGM49 NaN -0.091 0.795 0.593
CGM51 NaN -0.086 0.754 0.049
CGM52 -0.647 -0.077 0.146 0.219
CGM55 -0.223 0.014 -0.121 0.207
CGM57 -0.750 0.040 0.087 -0.004
CGM58 -0.547 0.004 -0.013 -0.069
Vanessa 0.000 -0.015 0.007 -0.128
Alignment variants (delta current 10 min after the epoch under each alignment):
same interval (main), > 150 per epoch: clustered mean -0.1910 nA (n = 585), isolated 0.0496 nA (n = 1718); ratio -3.85; Welch p = 0.059
same interval (main), > 30 per epoch: clustered mean -0.0316 nA (n = 12595), isolated 0.0496 nA (n = 1718); ratio -0.64; Welch p = 0.15
curve by count: [ 0.059 0.05 0.021 -0.049 -0.001 -0.045 -0.099 0.125 -0.218 -0.102]
SPW-Rs after the reading, > 150 per epoch: clustered mean -0.1018 nA (n = 585), isolated -0.0114 nA (n = 1720); ratio 8.90; Welch p = 0.66
SPW-Rs after the reading, > 30 per epoch: clustered mean -0.0396 nA (n = 12609), isolated -0.0114 nA (n = 1720); ratio 3.46; Welch p = 0.75
curve by count: [ 0.043 -0.011 0.016 -0.038 -0.04 -0.018 -0.093 0.042 -0.088 -0.155]
delta to the next reading, > 150 per epoch: clustered mean 0.0176 nA (n = 585), isolated -0.0273 nA (n = 1716); ratio -0.65; Welch p = 0.78
delta to the next reading, > 30 per epoch: clustered mean -0.0224 nA (n = 12583), isolated -0.0273 nA (n = 1716); ratio 0.82; Welch p = 0.94
curve by count: [ 0.004 -0.027 0.021 0.003 -0.016 0.01 -0.163 0.127 0.021 -0.002]
Result: reproduced under reading (a), resting on the alignment assumption. Epochs with more than 150 SPW-Rs are followed 10 min later by a mean delta current of -0.191 nA (585 epochs, 13 rats), against +0.050 nA after epochs with 1 to 5 SPW-Rs (1,718 epochs): negative and 3.9 times the magnitude, though the difference is not significant by a Welch test (p = 0.059). Under reading (b), more than 30 per epoch, the mean is -0.032 nA, smaller in magnitude than the isolated value, so (b) does not match. The curve over epoch counts is not monotonic (the 120 to 150 group is +0.125 nA), and in 3 of the 13 rats with clustered epochs their mean is not negative. Under the alignments one reading apart, counting SPW-Rs over the 5 min after each reading gives -0.102 nA after clustered epochs against -0.011 nA after isolated ones (reading (a), Welch p = 0.66), which still matches, while taking the change to the next reading gives +0.018 against -0.027 nA under (a) and -0.022 against -0.027 nA under (b), which does not.
edfig4-ripple-features¶
Of the SPW-R features, the inter-SPW-R interval (clustering) most strongly modulates the glucose effect, then amplitude, then duration, and ripple frequency least.
For each SPW-R, the delta current of the reading 10 min after the 5-min interval holding it, as in
the Extended Data Fig. 4b legend; SPW-Rs pooled over rats and binned by the percentile of each
feature (deciles, a; 5-percentile bins, b). Features from the Methods: amplitude is the maximum
normalized squared signal (peak_normed_power), duration is stop minus start, frequency is
peak_frequency, and the inter-SPW-R interval is the time from the previous SPW-R in the same
file (the claims file does not say previous or next; previous is a standard choice). The
influence of a feature is the range of the bin means (claims file).
The main analysis pairs the SPW-R count and the glucose change of the same 5-min interval, as in fig2-ripple-rate-predicts-glucose (an assumption). The two alignments one reading apart are shown as variants: counting SPW-Rs over the 5 min after each reading, or taking the change to the next reading, which is the same as moving the offset by one reading.
evs = []
for s in RATS:
d = SUBJ[s]
e = d["ev"].copy()
e["interval"] = e.groupby("file").t.diff()
b = d["b"].copy()
b["delta_10"] = b.groupby("run").delta.shift(-2)
for name, k in ALIGN.items():
b[f"d10 {name}"] = b.groupby("run").delta.shift(-(2 + k))
b = b[np.isfinite(b.rate)]
# the reading whose 5-min interval holds the SPW-R: the first reading at or after it
idx = np.searchsorted(b.t.values, e.t.values, side="left")
ok = idx < len(b)
idx = np.where(ok, idx, 0)
inside = ok & (b.t.values[idx] - e.t.values < S.STEP)
e["delta_10"] = np.where(inside, b.delta_10.values[idx], np.nan)
for name in ALIGN:
e[f"d10 {name}"] = np.where(inside, b[f"d10 {name}"].values[idx], np.nan)
e["rat"] = s
evs.append(e[np.isfinite(e.delta_10)])
evs = pd.concat(evs)
evs["duration_ms"] = 1000 * evs.duration
print(f"{len(evs)} SPW-Rs from {evs.rat.nunique()} rats with a delta current 10 min later")
feats = {"interval": "interval", "amplitude": "peak_normed_power", "duration": "duration_ms", "frequency": "peak_frequency"}
res = []
for nb, label in ((10, "a: deciles"), (20, "b: 5-percentile bins")):
for f, col in feats.items():
x = evs[[col, "delta_10"]].dropna()
q = pd.qcut(x[col].rank(method="first"), nb, labels=False)
m = x.groupby(q).delta_10.mean()
res.append({"bins": label, "feature": f, "n": len(x), "range_nA": m.max() - m.min(),
"lowest_bin_mean": m.iloc[0], "highest_bin_mean": m.iloc[-1], "min_bin": int(m.idxmin()), "max_bin": int(m.idxmax())})
R4 = pd.DataFrame(res)
print(R4.round(4).to_string())
for label, g in R4.groupby("bins"):
print(label, "order by range:", " > ".join(g.sort_values("range_nA", ascending=False).feature))
# per rat, deciles
pr = []
for s, g in evs.groupby("rat"):
row = {"rat": s, "n": len(g)}
for f, col in feats.items():
x = g[[col, "delta_10"]].dropna()
m = x.groupby(pd.qcut(x[col].rank(method="first"), 10, labels=False)).delta_10.mean()
row[f] = m.max() - m.min()
pr.append(row)
pr = pd.DataFrame(pr)
pr["largest"] = pr[list(feats)].idxmax(axis=1)
print(pr.round(4).to_string()); print("largest range per rat:", pr.largest.value_counts().to_dict())
# uncertainty: resample rats with replacement, decile bins fixed at the pooled edges
rng = np.random.default_rng(4)
rat_ids = sorted(evs.rat.unique())
sums = {}
for f, col in feats.items():
x = evs[[col, "delta_10", "rat"]].dropna()
q = pd.qcut(x[col].rank(method="first"), 10, labels=False)
g = x.groupby(["rat", q]).delta_10.agg(["sum", "size"]).unstack(fill_value=0)
sums[f] = (g["sum"].reindex(rat_ids, fill_value=0).values, g["size"].reindex(rat_ids, fill_value=0).values)
boot = []
for i in range(1000):
pick = rng.integers(0, len(rat_ids), len(rat_ids))
row = {}
for f, (sm, n) in sums.items():
m = sm[pick].sum(axis=0) / np.maximum(n[pick].sum(axis=0), 1)
row[f] = m.max() - m.min()
boot.append(row)
boot = pd.DataFrame(boot)
print("bootstrap over rats (1000), range of decile means, 2.5 to 97.5%:",
{f: (round(boot[f].quantile(0.025), 3), round(boot[f].quantile(0.975), 3)) for f in feats})
print("share of bootstraps in which each feature has the largest range:", boot.idxmax(axis=1).value_counts(normalize=True).round(3).to_dict())
print("share of bootstraps with the full order interval > amplitude > duration > frequency:",
round(float(((boot.interval > boot.amplitude) & (boot.amplitude > boot.duration) & (boot.duration > boot.frequency)).mean()), 3))
al = []
for name in ALIGN:
for nb_ in (10, 20):
row = {"alignment": name, "bins": nb_}
for f, col in feats.items():
x = evs[[col, f"d10 {name}"]].dropna()
m = x.groupby(pd.qcut(x[col].rank(method="first"), nb_, labels=False))[f"d10 {name}"].mean()
row[f] = m.max() - m.min()
row["order"] = " > ".join(sorted(feats, key=lambda f: -row[f]))
al.append(row)
print("Alignment variants, range of bin means (nA):")
print(pd.DataFrame(al).round(4).to_string())
1083233 SPW-Rs from 18 rats with a delta current 10 min later
bins feature n range_nA lowest_bin_mean highest_bin_mean min_bin max_bin 0 a: deciles interval 1083173 0.0703 0.0040 -0.0274 4 0 1 a: deciles amplitude 1083233 0.0400 -0.0401 -0.0262 3 8 2 a: deciles duration 1083233 0.0687 -0.0502 0.0060 6 9 3 a: deciles frequency 1083233 0.0527 -0.0353 -0.0671 9 3 4 b: 5-percentile bins interval 1083173 0.0966 0.0153 -0.0183 9 0 5 b: 5-percentile bins amplitude 1083233 0.0609 -0.0587 -0.0247 0 17 6 b: 5-percentile bins duration 1083233 0.0967 -0.0462 0.0084 13 19 7 b: 5-percentile bins frequency 1083233 0.0601 -0.0340 -0.0680 19 4 a: deciles order by range: interval > duration > frequency > amplitude b: 5-percentile bins order by range: duration > interval > amplitude > frequency
rat n interval amplitude duration frequency largest
0 CGM1 141750 0.0326 0.1246 0.0985 0.1112 amplitude
1 CGM2 30905 0.1789 0.1056 0.2061 0.1882 duration
2 CGM3 266538 0.2146 0.0751 0.1721 0.0895 interval
3 CGM30 71007 0.0330 0.0273 0.0199 0.2198 frequency
4 CGM31 158246 0.0408 0.0709 0.0441 0.0608 amplitude
5 CGM32 100343 0.0602 0.0304 0.0248 0.0620 frequency
6 CGM39 15742 0.0506 0.0744 0.0244 0.0409 amplitude
7 CGM40 51486 0.0835 0.1550 0.0604 0.0530 amplitude
8 CGM41 119905 0.0393 0.0468 0.0295 0.0417 amplitude
9 CGM46 38473 0.1538 0.0687 0.0450 0.0842 interval
10 CGM48 948 0.0959 0.0585 0.0553 0.0509 interval
11 CGM49 3958 0.4211 0.3058 0.2546 0.1739 interval
12 CGM51 7570 0.1864 0.1120 0.1689 0.1243 interval
13 CGM52 12651 0.2699 0.1123 0.0991 0.0798 interval
14 CGM55 11829 0.1040 0.1091 0.0593 0.0713 amplitude
15 CGM57 13672 0.1076 0.0715 0.1208 0.1052 duration
16 CGM58 27487 0.0682 0.0784 0.0531 0.0789 frequency
17 Vanessa 10723 0.0967 0.1613 0.1017 0.1233 amplitude
largest range per rat: {'amplitude': 7, 'interval': 6, 'frequency': 3, 'duration': 2}
bootstrap over rats (1000), range of decile means, 2.5 to 97.5%: {'interval': (np.float64(0.028), np.float64(0.131)), 'amplitude': (np.float64(0.019), np.float64(0.074)), 'duration': (np.float64(0.014), np.float64(0.1)), 'frequency': (np.float64(0.018), np.float64(0.085))}
share of bootstraps in which each feature has the largest range: {'interval': 0.647, 'duration': 0.14, 'amplitude': 0.136, 'frequency': 0.077}
share of bootstraps with the full order interval > amplitude > duration > frequency: 0.015
Alignment variants, range of bin means (nA):
alignment bins interval amplitude duration frequency order
0 same interval (main) 10 0.0703 0.0400 0.0687 0.0527 interval > duration > frequency > amplitude
1 same interval (main) 20 0.0966 0.0609 0.0967 0.0601 duration > interval > amplitude > frequency
2 SPW-Rs after the reading 10 0.0307 0.0465 0.0470 0.0596 frequency > duration > amplitude > interval
3 SPW-Rs after the reading 20 0.0456 0.0672 0.0860 0.0851 duration > frequency > amplitude > interval
4 delta to the next reading 10 0.0563 0.0381 0.0433 0.0392 interval > duration > frequency > amplitude
5 delta to the next reading 20 0.0671 0.0534 0.0762 0.0588 duration > interval > frequency > amplitude
Result: partially reproduced, resting on the alignment assumption. With deciles (a), the inter-SPW-R interval has the largest range of the +10-min effect across bins (0.070 nA), then duration (0.069), frequency (0.053) and amplitude (0.040), so the interval leads but the rest of the order differs. With 5-percentile bins (b), duration (0.0967) and interval (0.0966) are tied ahead of amplitude (0.061) and frequency (0.060). Per rat, amplitude has the largest range in 7 of 18 rats and the interval in 6. Resampling rats (1,000 bootstraps, decile edges fixed), the interval has the largest range in 64.7% of resamples and the full order interval > amplitude > duration > frequency occurs in 1.5%. The decile with the shortest intervals has the least negative effect (+0.004 nA, against -0.027 nA for the longest), the opposite direction to fig2g-clustered-ripples. Under the alignments one reading apart, taking the change to the next reading keeps the interval first with deciles (0.056 nA, then duration 0.043), so the outcome is the same; counting SPW-Rs over the 5 min after each reading puts the interval last under both bin widths (deciles: frequency 0.060, duration 0.047, amplitude 0.047, interval 0.031), which matches neither case of the criterion.
fig2h-phase-reset¶
The phase dispersion of the glucose rhythm decreases after SPW-Rs, more so for higher SPW-R rates.
As the claims file: glucose phase is the Hilbert phase of delta current band-passed at periods of 15 to 75 min (a) or 15 to 90 min (b) (second-order Butterworth, zero phase, within each run of at least 24 h, a standard choice); peaks are local maxima of SPW-R rate; peaks are pooled over rats and grouped into deciles of peak rate; for each decile the circular s.d. of glucose phase is taken at offsets -100 to +100 min; the decrease is the baseline (offset 0, a; mean over -100 to 0 min, b) minus the minimum from 0 to +30 min. As the minimum over seven offsets favours a decrease, the baselines minus the mean over 0 to +30 min are shown as a variant.
FS = 1 / S.STEP
def band_phase(b, lo_min, hi_min, col="delta", min_len=288):
"""Hilbert phase and envelope of `col` band-passed between periods lo_min and hi_min, per run."""
ph = np.full(len(b), np.nan); env = np.full(len(b), np.nan)
sos = signal.butter(2, [1 / (hi_min * 60), 1 / (lo_min * 60)], btype="band", fs=FS, output="sos")
vals = b[col].values
for _, idx in b.groupby("run").indices.items():
x = vals[idx]
if len(idx) < min_len:
continue
x = np.where(np.isfinite(x), x, np.nanmean(x))
z = signal.hilbert(signal.sosfiltfilt(sos, x - x.mean()))
ph[idx] = np.angle(z); env[idx] = np.abs(z)
return ph, env
def circ_sd(a):
a = a[np.isfinite(a)]
if len(a) < 3:
return np.nan
return np.sqrt(-2 * np.log(np.abs(np.mean(np.exp(1j * a)))))
OFFS = np.arange(-20, 21)
for band in ((15, 75), (15, 90)):
pk = []
for s in RATS:
b = SUBJ[s]["b"]
ph, _ = band_phase(b, *band)
r = b.rate.values
for _, idx in b.groupby("run").indices.items():
rr = np.where(np.isfinite(r[idx]), r[idx], -1)
peaks, _ = signal.find_peaks(rr)
for p in peaks:
if not np.isfinite(r[idx][p]):
continue
row = np.full(len(OFFS), np.nan)
for j, o in enumerate(OFFS):
if 0 <= p + o < len(idx):
row[j] = ph[idx][p + o]
if np.isfinite(row[OFFS == 0][0]):
pk.append((s, r[idx][p], row))
rates = np.array([x[1] for x in pk]); P = np.array([x[2] for x in pk])
dec = pd.qcut(pd.Series(rates).rank(method="first"), 10, labels=False).values
tab = []
for k in range(10):
csd = np.array([circ_sd(P[dec == k, j]) for j in range(len(OFFS))])
post = csd[(OFFS >= 0) & (OFFS <= 6)].min()
tab.append({"decile": k + 1, "peaks": int((dec == k).sum()), "median_rate": np.median(rates[dec == k]),
"csd_0": csd[OFFS == 0][0], "csd_pre_mean": csd[OFFS <= 0].mean(), "min_0_30": post,
"decrease_a": csd[OFFS == 0][0] - post, "decrease_b": csd[OFFS <= 0].mean() - post,
"decrease_a_mean": csd[OFFS == 0][0] - csd[(OFFS >= 0) & (OFFS <= 6)].mean(),
"decrease_b_mean": csd[OFFS <= 0].mean() - csd[(OFFS >= 0) & (OFFS <= 6)].mean()})
tab = pd.DataFrame(tab)
print(f"band {band[0]}-{band[1]} min: {len(rates)} peaks of SPW-R rate from {len(set(x[0] for x in pk))} rats")
print(tab.round(3).to_string())
for col in ("decrease_a", "decrease_b", "decrease_a_mean", "decrease_b_mean"):
rho = stats.spearmanr(tab[col], tab.median_rate)
print(f" {col}: top decile {tab[col].iloc[-1]:.3f}; Spearman with median rate rho = {rho.statistic:.3f}, p = {rho.pvalue:.3g}")
band 15-75 min: 5839 peaks of SPW-R rate from 18 rats decile peaks median_rate csd_0 csd_pre_mean min_0_30 decrease_a decrease_b decrease_a_mean decrease_b_mean 0 1 584 1.6 2.117 2.360 1.886 0.231 0.474 0.152 0.396 1 2 584 4.8 2.353 2.343 2.239 0.115 0.104 -0.033 -0.044 2 3 584 8.2 2.351 2.405 2.044 0.307 0.362 0.040 0.094 3 4 584 11.2 2.562 2.757 2.268 0.294 0.489 0.085 0.280 4 5 584 13.8 2.743 2.710 2.575 0.167 0.135 -0.098 -0.131 5 6 583 16.2 2.598 2.621 2.308 0.289 0.313 0.183 0.206 6 7 584 18.6 2.429 2.350 2.242 0.187 0.108 0.091 0.012 7 8 584 21.2 2.486 2.610 2.104 0.381 0.506 0.213 0.337 8 9 584 25.1 1.973 2.519 1.834 0.139 0.685 -0.027 0.519 9 10 584 32.0 2.133 2.589 2.083 0.051 0.506 -0.089 0.366 decrease_a: top decile 0.051; Spearman with median rate rho = -0.236, p = 0.511 decrease_b: top decile 0.506; Spearman with median rate rho = 0.552, p = 0.0984 decrease_a_mean: top decile -0.089; Spearman with median rate rho = -0.067, p = 0.855 decrease_b_mean: top decile 0.366; Spearman with median rate rho = 0.345, p = 0.328
band 15-90 min: 5839 peaks of SPW-R rate from 18 rats decile peaks median_rate csd_0 csd_pre_mean min_0_30 decrease_a decrease_b decrease_a_mean decrease_b_mean 0 1 584 1.6 2.135 2.336 1.899 0.235 0.437 0.158 0.359 1 2 584 4.8 2.476 2.422 2.268 0.208 0.154 0.054 -0.001 2 3 584 8.2 2.325 2.421 2.096 0.230 0.326 -0.036 0.060 3 4 584 11.2 2.910 2.746 2.314 0.596 0.432 0.322 0.158 4 5 584 13.8 2.592 2.727 2.530 0.062 0.197 -0.114 0.022 5 6 583 16.2 2.733 2.601 2.267 0.466 0.334 0.309 0.178 6 7 584 18.6 2.432 2.367 2.144 0.288 0.223 0.179 0.114 7 8 584 21.2 2.451 2.602 2.096 0.354 0.506 0.199 0.350 8 9 584 25.1 1.931 2.455 1.829 0.102 0.626 -0.040 0.484 9 10 584 32.0 2.147 2.618 2.125 0.022 0.494 -0.079 0.393 decrease_a: top decile 0.022; Spearman with median rate rho = -0.236, p = 0.511 decrease_b: top decile 0.494; Spearman with median rate rho = 0.552, p = 0.0984 decrease_a_mean: top decile -0.079; Spearman with median rate rho = -0.188, p = 0.603 decrease_b_mean: top decile 0.393; Spearman with median rate rho = 0.564, p = 0.0897
Result: not reproduced. Under neither band nor baseline is the decrease of phase dispersion correlated with SPW-R rate at p < 0.05: Spearman rho -0.24 (p = 0.51) for baseline (a) and 0.55 (p = 0.098) for baseline (b), the same in both bands. The decrease in the top decile is positive under all four readings (0.022 to 0.506), but it is positive in every decile, including the lowest (0.23 to 0.47). With the mean over 0 to +30 min in place of the minimum, rho is -0.19 to 0.56 (p = 0.09 to 0.86), and the top decile's decrease is negative with baseline (a).
edfig5-cross-frequency-coupling¶
SPW-R rate and glucose show power-power and phase-power coupling at periods of 15 to 75 min, with the faster SPW-R rhythm (about 40 min) modulating the slower glucose rhythm (about 60 min).
As the claims file: band-pass at log-spaced periods (here 12 to 120 min: a band centred at 10 min would sit at the Nyquist period of 5-min readings), Hilbert envelopes; power-power is the correlation of the SPW-R-rate and glucose envelopes for each pair of periods per rat, against 100 circular shuffles of SPW-R rate per rat, with a paired t-test across rats on the mean over 15 to 60 min. Phase-power is SPW-R rate (z-scored per rat) by glucose phase (the 15 to 75 min band), with a Rayleigh test of the rate-weighted phases pooled over rats (weights: z-scored rate minus its minimum, so that they are positive; a choice). Both signals are filtered only over stretches where both are recorded. The shortest stretch that is filtered is a choice: 12 h (six cycles of the longest period, 120 min) is the main reading, with 24 h and 6 h as variants. Bins are autocorrelated, so the pooled Rayleigh p is also compared with circular shifts of SPW-R rate.
rng = np.random.default_rng(5)
PER = np.round(np.geomspace(12, 120, 12), 1)
def comod(eg, er):
out = np.full((len(PER), len(PER)), np.nan)
for i in range(len(PER)):
for j in range(len(PER)):
ok = np.isfinite(eg[i]) & np.isfinite(er[j])
if ok.sum() > 50:
out[i, j] = np.corrcoef(eg[i][ok], er[j][ok])[0, 1]
return out
def both_segments(b):
"""Rows with both SPW-R rate and delta current, renumbered into runs of consecutive readings."""
ok = np.isfinite(b.rate.values) & np.isfinite(b.delta.values)
c = b[ok].copy()
pos = np.flatnonzero(ok)
c["run"] = np.cumsum(np.r_[True, (np.diff(pos) != 1) | (c.run.values[1:] != c.run.values[:-1])])
return c
patch = (PER >= 15) & (PER <= 60)
def rayleigh_w(ph, z):
"""Resultant length of phases weighted by z minus its minimum, and the Rayleigh p for n phases."""
w = z - z.min()
Rw = np.abs(np.sum(w * np.exp(1j * ph))) / w.sum()
n = len(ph)
return Rw, np.exp(np.sqrt(1 + 4 * n + 4 * (n ** 2 - (n * Rw) ** 2)) - (1 + 2 * n))
def coupling(min_len, show_matrix=False):
"""Power-power and phase-power coupling over stretches of at least min_len bins with both signals."""
obs_m, sh_m, rats5, mats, hrs5 = [], [], [], [], {}
for s in RATS:
b = both_segments(SUBJ[s]["b"])
eg = np.array([band_phase(b, max(p / 1.25, 10.5), p * 1.25, col="delta", min_len=min_len)[1] for p in PER])
if np.all(np.isnan(eg)):
continue
er = np.array([band_phase(b, max(p / 1.25, 10.5), p * 1.25, col="rate", min_len=min_len)[1] for p in PER])
hrs5[s] = round(float(np.isfinite(eg[0]).sum()) * 5 / 60, 1)
M = comod(eg, er)
sh = []
for i in range(100):
sb = S.circular_shuffle(b, rng)
esh = np.array([band_phase(sb, max(p / 1.25, 10.5), p * 1.25, col="rate", min_len=min_len)[1] for p in PER])
sh.append(np.nanmean(comod(eg, esh)[np.ix_(patch, patch)]))
rats5.append(s); mats.append(M); obs_m.append(np.nanmean(M[np.ix_(patch, patch)])); sh_m.append(np.mean(sh))
obs_m, sh_m = np.array(obs_m), np.array(sh_m)
t = stats.ttest_rel(obs_m, sh_m)
print(f"--- stretches of at least {min_len * 5 / 60:.0f} h with both SPW-R rate and delta current: {len(rats5)} rats; hours used per rat: {hrs5}")
print(f"power-power, mean over 15-60 min: observed {obs_m.mean():.3f} (per rat {np.round(obs_m, 3)}), shuffled {sh_m.mean():.3f}; paired t = {t.statistic:.2f}, p = {t.pvalue:.2g}")
MM = np.nanmean(mats, axis=0)
i, j = np.unravel_index(np.nanargmax(MM), MM.shape)
off = np.nanmean([MM[a, c] for a in range(len(PER)) for c in range(len(PER)) if a > c])
below = np.nanmean([MM[a, c] for a in range(len(PER)) for c in range(len(PER)) if a < c])
print(f"mean comodulogram peak R = {MM[i, j]:.3f} at glucose period {PER[i]} min and SPW-R period {PER[j]} min; "
f"mean R with glucose period longer than SPW-R period {off:.3f}, shorter {below:.3f}, diagonal {np.nanmean(np.diag(MM)):.3f}")
if show_matrix:
print(pd.DataFrame(MM, index=[f"glu {p}" for p in PER], columns=[f"spw {p}" for p in PER]).round(2).to_string())
# phase-power: glucose phase (15-75 min band) and z-scored SPW-R rate, pooled over rats
ph_all, z_all, segs5, per_rat_p = [], [], [], {}
for s in rats5:
b = both_segments(SUBJ[s]["b"])
ph, _ = band_phase(b, 15, 75, min_len=min_len)
ok = np.isfinite(ph)
b = b[ok].assign(ph=ph[ok])
b["z"] = (b.rate - b.rate.mean()) / b.rate.std()
ph_all.append(b.ph.values); z_all.append(b.z.values); segs5.append(b)
per_rat_p[s] = rayleigh_w(b.ph.values, b.z.values)[1]
ph_all = np.concatenate(ph_all); z_all = np.concatenate(z_all); w = z_all - z_all.min()
Rw, p_ray = rayleigh_w(ph_all, z_all)
mean_phase = np.degrees(np.angle(np.sum(w * np.exp(1j * ph_all))))
print(f"phase-power: {len(ph_all)} bins; weighted resultant length {Rw:.4f}, Rayleigh p = {p_ray:.3g}, preferred glucose phase {mean_phase:.0f} deg")
# null for the pooled resultant length: SPW-R rate circularly shifted within each stretch (bins are autocorrelated)
null = np.array([rayleigh_w(ph_all, np.concatenate([S.circular_shuffle(b, rng, x="z").z.values for b in segs5]))[0] for i in range(200)])
print(f" circular-shift null of the resultant length (200): 95th percentile {np.percentile(null, 95):.4f}; share >= observed {np.mean(null >= Rw):.3f}")
print(f" per-rat weighted Rayleigh p < 0.05 in {sum(p < 0.05 for p in per_rat_p.values())} of {len(per_rat_p)} rats")
pb = pd.cut(ph_all, np.linspace(-np.pi, np.pi, 13))
print(" mean z-scored SPW-R rate by glucose phase:", pd.Series(z_all).groupby(pb, observed=True).mean().round(3).values)
for min_len, show in ((288, False), (144, True), (72, False)):
coupling(min_len, show)
--- stretches of at least 24 h with both SPW-R rate and delta current: 3 rats; hours used per rat: {'CGM1': 49.1, 'CGM3': 49.9, 'CGM32': 24.0}
power-power, mean over 15-60 min: observed 0.173 (per rat [ 0.218 -0.039 0.341]), shuffled -0.014; paired t = 1.64, p = 0.24
mean comodulogram peak R = 0.487 at glucose period 79.0 min and SPW-R period 64.0 min; mean R with glucose period longer than SPW-R period 0.156, shorter 0.130, diagonal 0.203
phase-power: 1476 bins; weighted resultant length 0.0118, Rayleigh p = 0.815, preferred glucose phase -154 deg
circular-shift null of the resultant length (200): 95th percentile 0.0620; share >= observed 0.955 per-rat weighted Rayleigh p < 0.05 in 2 of 3 rats mean z-scored SPW-R rate by glucose phase: [-0.062 0.014 0. -0.172 -0.166 0.037 -0.085 0.028 0.057 0.211 0.083 0.09 ]
--- stretches of at least 12 h with both SPW-R rate and delta current: 16 rats; hours used per rat: {'CGM1': 196.2, 'CGM2': 42.9, 'CGM3': 252.8, 'CGM30': 108.7, 'CGM31': 233.8, 'CGM32': 175.2, 'CGM39': 31.5, 'CGM40': 63.1, 'CGM41': 172.9, 'CGM46': 41.1, 'CGM51': 15.4, 'CGM52': 17.9, 'CGM55': 16.0, 'CGM57': 17.9, 'CGM58': 55.6, 'Vanessa': 15.5}
power-power, mean over 15-60 min: observed 0.064 (per rat [ 0.11 0.242 -0.062 -0.127 -0.022 -0.047 0.181 -0.028 0.016 0.08
0.113 0.073 -0.105 0.012 0.252 0.328]), shuffled -0.007; paired t = 2.57, p = 0.021
mean comodulogram peak R = 0.216 at glucose period 79.0 min and SPW-R period 79.0 min; mean R with glucose period longer than SPW-R period 0.072, shorter 0.047, diagonal 0.096
spw 12.0 spw 14.8 spw 18.2 spw 22.5 spw 27.7 spw 34.2 spw 42.1 spw 51.9 spw 64.0 spw 79.0 spw 97.3 spw 120.0
glu 12.0 -0.03 0.01 0.05 0.05 -0.01 -0.01 0.03 -0.02 -0.05 0.01 0.03 -0.00
glu 14.8 -0.03 0.00 0.07 0.06 0.00 -0.00 0.03 -0.00 -0.03 0.03 0.03 -0.00
glu 18.2 -0.03 0.01 0.09 0.07 -0.00 -0.01 0.03 0.01 -0.01 0.05 0.04 0.00
glu 22.5 -0.01 0.03 0.11 0.09 0.02 0.01 0.04 0.02 0.01 0.08 0.07 0.03
glu 27.7 -0.01 0.04 0.11 0.09 0.05 0.04 0.04 0.01 0.02 0.09 0.07 0.04
glu 34.2 0.01 0.04 0.11 0.10 0.06 0.09 0.09 0.04 0.03 0.10 0.08 0.05
glu 42.1 0.03 0.04 0.08 0.08 0.05 0.05 0.13 0.10 0.09 0.15 0.11 0.06
glu 51.9 0.03 0.03 0.06 0.11 0.08 0.04 0.10 0.13 0.13 0.16 0.09 0.03
glu 64.0 0.04 0.05 0.10 0.12 0.09 0.04 0.05 0.06 0.14 0.20 0.14 0.07
glu 79.0 0.06 0.06 0.14 0.16 0.10 0.06 0.05 0.03 0.09 0.22 0.15 0.08
glu 97.3 0.06 0.06 0.15 0.17 0.11 0.06 0.04 0.01 0.06 0.18 0.15 0.10
glu 120.0 0.07 0.06 0.14 0.19 0.13 0.07 0.04 0.02 0.07 0.18 0.14 0.10
phase-power: 17477 bins; weighted resultant length 0.0334, Rayleigh p = 3.22e-09, preferred glucose phase 130 deg
circular-shift null of the resultant length (200): 95th percentile 0.0129; share >= observed 0.000 per-rat weighted Rayleigh p < 0.05 in 10 of 16 rats mean z-scored SPW-R rate by glucose phase: [ 0.041 -0.014 -0.095 -0.12 -0.104 -0.116 -0.078 0.009 0.067 0.175 0.159 0.069]
--- stretches of at least 6 h with both SPW-R rate and delta current: 16 rats; hours used per rat: {'CGM1': 207.0, 'CGM2': 49.8, 'CGM3': 269.5, 'CGM30': 131.8, 'CGM31': 259.8, 'CGM32': 175.2, 'CGM39': 50.7, 'CGM40': 71.8, 'CGM41': 200.0, 'CGM46': 49.4, 'CGM51': 15.4, 'CGM52': 17.9, 'CGM55': 16.0, 'CGM57': 17.9, 'CGM58': 71.1, 'Vanessa': 24.1}
power-power, mean over 15-60 min: observed 0.048 (per rat [ 0.111 0.244 -0.057 -0.115 -0.019 -0.047 0.131 -0.006 0.014 0.078
0.113 0.073 -0.105 0.012 0.257 0.09 ]), shuffled -0.005; paired t = 2.48, p = 0.025
mean comodulogram peak R = 0.220 at glucose period 79.0 min and SPW-R period 79.0 min; mean R with glucose period longer than SPW-R period 0.058, shorter 0.043, diagonal 0.083
phase-power: 19529 bins; weighted resultant length 0.0306, Rayleigh p = 1.1e-08, preferred glucose phase 135 deg
circular-shift null of the resultant length (200): 95th percentile 0.0144; share >= observed 0.000 per-rat weighted Rayleigh p < 0.05 in 10 of 16 rats mean z-scored SPW-R rate by glucose phase: [ 0.037 -0.007 -0.089 -0.109 -0.092 -0.118 -0.086 0.012 0.065 0.152 0.157 0.074]
Result: partially reproduced. Over stretches of at least 12 h with both signals (16 rats), power-power coupling over 15 to 60 min exceeds the shuffles (mean R 0.064 against -0.007, paired t = 2.57, p = 0.021), and the rate-weighted glucose phases are non-uniform (Rayleigh p = 3.2e-09, resultant length 0.033, above all 200 circular shifts, whose 95th percentile is 0.013; p < 0.05 in 10 of 16 rats singly). The comodulogram peak (R 0.216) lies on the diagonal, at 79.0 min for both signals, so the third condition, a peak at a longer glucose period, fails; the mean correlation where the glucose period is longer than the SPW-R period (below the diagonal of the printed matrix) is 0.072 against 0.047 where it is shorter. With 6-h stretches the result is the same (p = 0.025 and 1.1e-08, peak at 79.0 min for both). With 24-h stretches only 3 rats remain and neither test reaches p < 0.05 (p = 0.24 and 0.82). The paper's p bounds (1e-24, 1e-18) are not reached.
edfig1-glucose-rhythm¶
Interstitial glucose fluctuations show an approximately 60-min oscillation, with power elevated above a linear best fit to the spectrum at periods of 15 to 90 min.
As Extended Data Fig. 1g-i: power spectral density of delta current per rat (Welch, 24-h segments within runs, a standard choice), a line fitted on log-log axes (the claims file's standard choice) and subtracted, averaged across rats; and the mean autocorrelogram of delta current. All subjects with at least one 24-h run of readings are used (the glucose series do not need the electrophysiology). The legend says "linear best fit" without the axes, so the fit is also run on semi-log and linear axes (power normalized per subject), over the whole spectrum or over 10 to 300 min, and with 12-h segments.
psd_rows, ac_rows, used = [], [], []
for s in sorted(SUBJ):
glu = SUBJ[s]["glu"].copy()
glu["delta"] = np.where(glu.run.eq(glu.run.shift(1)), glu.nA.diff(), np.nan)
segs = [g.delta.values[1:] for _, g in glu.groupby("run") if len(g) >= 289]
if not segs:
continue
ps = []
for x in segs:
f, p = signal.welch(np.nan_to_num(x), fs=FS, nperseg=288)
ps.append(p * len(x))
psd_rows.append(np.sum(ps, axis=0) / sum(len(x) for x in segs))
ac = []
for x in segs:
x = np.nan_to_num(x - np.nanmean(x))
c = np.correlate(x, x, "full")[len(x) - 1:len(x) - 1 + 37] / np.dot(x, x)
ac.append(c)
ac_rows.append(np.mean(ac, axis=0)); used.append(s)
F1 = f[1:]; period = 1 / F1 / 60
resid = []
for p in psd_rows:
lp = np.log10(p[1:]); coef = np.polyfit(np.log10(F1), lp, 1)
resid.append(lp - np.polyval(coef, np.log10(F1)))
resid = np.array(resid); mres = resid.mean(axis=0)
print(f"{len(used)} subjects: {used}")
tbl = pd.DataFrame({"period_min": period, "mean_residual_log10": mres, "rats_positive": (resid > 0).sum(axis=0)})
print(tbl[(tbl.period_min >= 10) & (tbl.period_min <= 300)].round(3).to_string())
band = (period >= 15) & (period <= 90)
print(f"residual over 15-90 min: mean {mres[band].mean():.3f}, positive at {np.sum(mres[band] > 0)} of {band.sum()} frequencies")
ipk = np.argmax(np.where((period >= 10) & (period <= 300), mres, -np.inf))
print(f"peak of the mean residual at {period[ipk]:.1f} min")
mac = np.mean(ac_rows, axis=0); lags_min = np.arange(37) * 5
pk, _ = signal.find_peaks(mac)
print("mean autocorrelogram:", dict(zip(lags_min[:25], np.round(mac[:25], 3))))
print("first side peak of the mean autocorrelogram at", lags_min[pk[0]] if len(pk) else None, "min",
f"(value {mac[pk[0]]:.3f})" if len(pk) else "")
per_rat_pk = [period[np.argmax(np.where((period >= 10) & (period <= 300), r, -np.inf))] for r in resid]
print("per subject peak of residual (min):", dict(zip(used, np.round(per_rat_pk, 1))))
fig, ax = plt.subplots(1, 2, figsize=(10, 3))
ax[0].semilogx(period, mres); ax[0].axhline(0, color="k", lw=0.5); ax[0].set_xlabel("period (min)"); ax[0].set_ylabel("log10 power - fit")
ax[1].plot(lags_min, mac); ax[1].set_xlabel("lag (min)"); plt.show()
# variants: axes of the line, the range it is fitted over, and the segment length
def spectra(nper):
out = []
for s in sorted(SUBJ):
glu = SUBJ[s]["glu"].copy()
glu["delta"] = np.where(glu.run.eq(glu.run.shift(1)), glu.nA.diff(), np.nan)
segs = [g.delta.values[1:] for _, g in glu.groupby("run") if len(g) >= nper + 1]
if segs:
ps = [signal.welch(np.nan_to_num(x), fs=FS, nperseg=nper)[1] * len(x) for x in segs]
out.append(np.sum(ps, axis=0) / sum(len(x) for x in segs))
return signal.welch(np.zeros(nper), fs=FS, nperseg=nper)[0], out
rows = []
for nper, seg in ((288, "24 h"), (144, "12 h")):
fq, P = spectra(nper)
F, per = fq[1:], 1 / fq[1:] / 60
for axes in ("log-log", "semi-log (log power)", "linear"):
for frange in ("all", "10-300 min"):
sel = np.ones(len(F), bool) if frange == "all" else (per >= 10) & (per <= 300)
res = []
for p in P:
y = p[1:] / p[1:].sum()
X = np.log10(F) if axes == "log-log" else F
Y = y if axes == "linear" else np.log10(y)
res.append(Y - np.polyval(np.polyfit(X[sel], Y[sel], 1), X))
m = np.mean(res, axis=0)
win = (per >= 10) & (per <= 300)
b15 = (per >= 15) & (per <= 90)
rows.append({"segment": seg, "axes": axes, "fit_over": frange, "subjects": len(P),
"peak_period_min": per[np.argmax(np.where(win, m, -np.inf))],
"positive_15_90": f"{int((m[b15] > 0).sum())} of {int(b15.sum())}",
"shortest_positive_period": per[win & (m > 0)].min() if np.any(win & (m > 0)) else np.nan,
"peak_in_54_66": 54 <= per[np.argmax(np.where(win, m, -np.inf))] <= 66})
print(pd.DataFrame(rows).round(1).to_string())
25 subjects: ['CGM1', 'CGM2', 'CGM3', 'CGM30', 'CGM31', 'CGM32', 'CGM36', 'CGM37', 'CGM39', 'CGM40', 'CGM41', 'CGM46', 'CGM48', 'CGM49', 'CGM51', 'CGM52', 'CGM55', 'CGM57', 'CGM58', 'CGM60', 'Vanessa', 'bruce', 'dt15', 'flex1', 'ros']
period_min mean_residual_log10 rats_positive
4 288.000 -1.376 0
5 240.000 -1.035 0
6 205.714 -0.669 0
7 180.000 -0.546 0
8 160.000 -0.341 1
9 144.000 -0.114 7
10 130.909 0.084 19
11 120.000 0.232 20
12 110.769 0.361 23
13 102.857 0.475 24
14 96.000 0.623 23
15 90.000 0.773 23
16 84.706 0.879 24
17 80.000 0.948 24
18 75.789 1.074 25
19 72.000 1.158 25
20 68.571 1.284 25
21 65.455 1.282 25
22 62.609 1.254 25
23 60.000 1.363 25
24 57.600 1.414 25
25 55.385 1.316 25
26 53.333 1.355 25
27 51.429 1.459 24
28 49.655 1.447 25
29 48.000 1.406 25
30 46.452 1.386 25
31 45.000 1.399 24
32 43.636 1.397 25
33 42.353 1.383 25
34 41.143 1.359 24
35 40.000 1.329 25
36 38.919 1.289 25
37 37.895 1.269 24
38 36.923 1.234 25
39 36.000 1.273 25
40 35.122 1.227 25
41 34.286 1.170 24
42 33.488 1.109 24
43 32.727 1.056 25
44 32.000 0.972 25
45 31.304 0.929 24
46 30.638 0.873 25
47 30.000 0.853 25
48 29.388 0.793 24
49 28.800 0.820 25
50 28.235 0.754 25
51 27.692 0.725 24
52 27.170 0.688 25
53 26.667 0.593 25
54 26.182 0.512 24
55 25.714 0.443 23
56 25.263 0.437 25
57 24.828 0.404 23
58 24.407 0.389 25
59 24.000 0.365 23
60 23.607 0.319 23
61 23.226 0.244 22
62 22.857 0.152 17
63 22.500 0.106 16
64 22.154 0.072 13
65 21.818 0.042 11
66 21.493 -0.045 8
67 21.176 -0.070 9
68 20.870 -0.164 7
69 20.571 -0.254 5
70 20.282 -0.253 4
71 20.000 -0.281 4
72 19.726 -0.337 5
73 19.459 -0.357 6
74 19.200 -0.377 5
75 18.947 -0.471 4
76 18.701 -0.509 3
77 18.462 -0.519 4
78 18.228 -0.567 5
79 18.000 -0.620 4
80 17.778 -0.585 3
81 17.561 -0.592 3
82 17.349 -0.677 2
83 17.143 -0.679 3
84 16.941 -0.700 2
85 16.744 -0.709 3
86 16.552 -0.721 2
87 16.364 -0.726 3
88 16.180 -0.692 3
89 16.000 -0.791 2
90 15.824 -0.790 2
91 15.652 -0.744 2
92 15.484 -0.728 1
93 15.319 -0.709 2
94 15.158 -0.668 2
95 15.000 -0.690 2
96 14.845 -0.748 0
97 14.694 -0.739 1
98 14.545 -0.707 2
99 14.400 -0.679 0
100 14.257 -0.636 1
101 14.118 -0.619 2
102 13.981 -0.567 1
103 13.846 -0.548 1
104 13.714 -0.550 2
105 13.585 -0.564 1
106 13.458 -0.533 1
107 13.333 -0.489 1
108 13.211 -0.477 0
109 13.091 -0.526 0
110 12.973 -0.478 0
111 12.857 -0.413 2
112 12.743 -0.431 0
113 12.632 -0.428 0
114 12.522 -0.441 1
115 12.414 -0.429 1
116 12.308 -0.442 0
117 12.203 -0.394 0
118 12.101 -0.422 0
119 12.000 -0.381 0
120 11.901 -0.385 0
121 11.803 -0.336 2
122 11.707 -0.330 1
123 11.613 -0.349 0
124 11.520 -0.292 1
125 11.429 -0.281 3
126 11.339 -0.257 3
127 11.250 -0.273 1
128 11.163 -0.225 3
129 11.077 -0.239 0
130 10.992 -0.263 1
131 10.909 -0.251 2
132 10.827 -0.210 3
133 10.746 -0.188 4
134 10.667 -0.158 7
135 10.588 -0.164 4
136 10.511 -0.179 5
137 10.435 -0.187 3
138 10.360 -0.189 5
139 10.286 -0.150 9
140 10.213 -0.125 8
141 10.141 -0.129 8
142 10.070 -0.106 10
143 10.000 -0.427 0
residual over 15-90 min: mean 0.392, positive at 51 of 81 frequencies
peak of the mean residual at 51.4 min
mean autocorrelogram: {np.int64(0): np.float64(1.0), np.int64(5): np.float64(0.682), np.int64(10): np.float64(0.447), np.int64(15): np.float64(0.172), np.int64(20): np.float64(-0.077), np.int64(25): np.float64(-0.222), np.int64(30): np.float64(-0.27), np.int64(35): np.float64(-0.247), np.int64(40): np.float64(-0.186), np.int64(45): np.float64(-0.119), np.int64(50): np.float64(-0.06), np.int64(55): np.float64(-0.02), np.int64(60): np.float64(0.002), np.int64(65): np.float64(0.012), np.int64(70): np.float64(0.013), np.int64(75): np.float64(0.01), np.int64(80): np.float64(0.005), np.int64(85): np.float64(0.0), np.int64(90): np.float64(-0.003), np.int64(95): np.float64(-0.005), np.int64(100): np.float64(-0.004), np.int64(105): np.float64(-0.002), np.int64(110): np.float64(0.002), np.int64(115): np.float64(0.003), np.int64(120): np.float64(0.001)}
first side peak of the mean autocorrelogram at 70 min (value 0.013)
per subject peak of residual (min): {'CGM1': np.float64(33.5), 'CGM2': np.float64(32.0), 'CGM3': np.float64(57.6), 'CGM30': np.float64(46.5), 'CGM31': np.float64(36.0), 'CGM32': np.float64(51.4), 'CGM36': np.float64(41.1), 'CGM37': np.float64(57.6), 'CGM39': np.float64(49.7), 'CGM40': np.float64(51.4), 'CGM41': np.float64(51.4), 'CGM46': np.float64(55.4), 'CGM48': np.float64(37.9), 'CGM49': np.float64(55.4), 'CGM51': np.float64(60.0), 'CGM52': np.float64(57.6), 'CGM55': np.float64(41.1), 'CGM57': np.float64(46.5), 'CGM58': np.float64(43.6), 'CGM60': np.float64(49.7), 'Vanessa': np.float64(51.4), 'bruce': np.float64(45.0), 'dt15': np.float64(45.0), 'flex1': np.float64(55.4), 'ros': np.float64(42.4)}
segment axes fit_over subjects peak_period_min positive_15_90 shortest_positive_period peak_in_54_66 0 24 h log-log all 25 51.4 51 of 81 21.8 False 1 24 h log-log 10-300 min 25 43.6 51 of 81 10.1 False 2 24 h semi-log (log power) all 25 10.1 38 of 81 10.0 False 3 24 h semi-log (log power) 10-300 min 25 10.1 37 of 81 10.0 False 4 24 h linear all 25 68.6 16 of 81 10.0 False 5 24 h linear 10-300 min 25 68.6 15 of 81 10.0 False 6 12 h log-log all 25 48.0 27 of 41 21.2 False 7 12 h log-log 10-300 min 25 45.0 25 of 41 10.1 False 8 12 h semi-log (log power) all 25 10.1 19 of 41 10.0 False 9 12 h semi-log (log power) 10-300 min 25 10.1 19 of 41 10.0 False 10 12 h linear all 25 72.0 8 of 41 10.0 False 11 12 h linear 10-300 min 25 72.0 8 of 41 10.0 False
Result: not reproduced. In 25 subjects the power of delta current above a log-log linear fit is elevated from about 22 to 131 min, with its peak at 51.4 min (per subject 32 to 60 min), outside 60 +- 6 min; below 22 min the residual is negative, so it is positive at only 51 of the 81 frequencies between 15 and 90 min. The mean autocorrelogram has its first side peak at 70 min with a value of 0.013, also outside the tolerance. Neither condition holds under any of the 12 variants of axes, fit range and segment length: the residual peak lies at 43.6 to 51.4 min (log-log), 10.1 min (semi-log) or 68.6 to 72.0 min (linear), and the residual is negative somewhere between 15 and 90 min in every variant. The spectrum has a broad hump of excess power from about 22 to 131 min, but the autocorrelogram shows no rhythm at its peak: it is -0.06 at a 50-min lag, and its first side peak, at 70 min, is 0.013. Neither the expected period nor the 15 to 90 min range is matched.
edfig1-glucose-circadian-meals¶
Glucose concentrations are modulated over the circadian cycle and rise after meals.
The meal times of the fasting experiment are not in the dandiset, so the meal part and its expected values (36 meals, 4 rats) cannot be tested. The circadian part of the criterion is run: current by hour of day averaged per subject, one-way ANOVA of hour on the hourly means pooled over subjects.
hourly = []
for s in sorted(SUBJ):
glu = SUBJ[s]["glu"]
h = S.clock(SUBJ[s]["ref"], glu.t).hour
m = glu.nA.groupby(h).mean()
if len(m) == 24:
hourly.append(pd.DataFrame({"subject": s, "hour": m.index, "nA": m.values}))
hourly = pd.concat(hourly)
an = stats.f_oneway(*[g.nA.values for _, g in hourly.groupby("hour")])
print(f"{hourly.subject.nunique()} subjects; one-way ANOVA of hour of day: F = {an.statistic:.2f}, p = {an.pvalue:.3g}")
print("mean current by hour (nA):", hourly.groupby("hour").nA.mean().round(1).to_dict())
# variant within rats: the linear decay of each run of at least 24 h removed, then each subject's mean
hw = []
for s in sorted(SUBJ):
glu = SUBJ[s]["glu"].copy()
glu = glu[glu.groupby("run").t.transform("size") >= 288].copy()
if not len(glu):
continue
glu["resid"] = glu.groupby("run").nA.transform(lambda x: x - np.polyval(np.polyfit(np.arange(len(x)), x.values, 1), np.arange(len(x))))
m = glu.resid.groupby(S.clock(SUBJ[s]["ref"], glu.t).hour).mean()
if len(m) == 24:
hw.append(pd.DataFrame({"subject": s, "hour": m.index, "nA": m.values - m.mean()}))
hw = pd.concat(hw)
an2 = stats.f_oneway(*[g.nA.values for _, g in hw.groupby("hour")])
fr = stats.friedmanchisquare(*[g.sort_values("subject").nA.values for _, g in hw.groupby("hour")])
print(f"within subjects, decay removed: {hw.subject.nunique()} subjects; one-way ANOVA F = {an2.statistic:.2f}, p = {an2.pvalue:.3g}; "
f"Friedman chi2 = {fr.statistic:.1f}, p = {fr.pvalue:.3g}")
print("mean deviation by hour (nA):", hw.groupby("hour").nA.mean().round(2).to_dict())
25 subjects; one-way ANOVA of hour of day: F = 0.46, p = 0.986
mean current by hour (nA): {0: 30.7, 1: 30.6, 2: 30.1, 3: 29.3, 4: 28.9, 5: 28.7, 6: 28.6, 7: 28.6, 8: 28.9, 9: 29.2, 10: 28.7, 11: 29.6, 12: 30.1, 13: 30.1, 14: 29.7, 15: 29.9, 16: 30.0, 17: 31.0, 18: 31.0, 19: 31.1, 20: 31.3, 21: 30.8, 22: 30.8, 23: 30.9}
within subjects, decay removed: 25 subjects; one-way ANOVA F = 1.69, p = 0.024; Friedman chi2 = 59.1, p = 5.17e-05
mean deviation by hour (nA): {0: 0.22, 1: 0.35, 2: 0.04, 3: -0.5, 4: -0.64, 5: -0.83, 6: -0.84, 7: -0.65, 8: -0.33, 9: 0.09, 10: -0.48, 11: 0.26, 12: 0.6, 13: 0.65, 14: 0.04, 15: 0.16, 16: -0.17, 17: 0.33, 18: 0.02, 19: 0.42, 20: 0.47, 21: 0.17, 22: 0.3, 23: 0.34}
Result: blocked (meal times are not in the dandiset). The circadian part of the criterion does not find a modulation: one-way ANOVA of hour of day on the hourly means of 25 subjects, F = 0.46, p = 0.99; the mean current ranges from 28.6 nA (06:00 to 07:00) to 31.3 nA (20:00). The criterion's test pools hourly means across subjects whose offsets differ by tens of nA. Within subjects, with each run's linear decay and each subject's mean removed, hour of day is significant (one-way ANOVA p = 0.024, Friedman p = 5.2e-05; lowest at 06:00 to 07:00, -0.84 nA, and highest at 13:00 to 14:00, +0.65 nA), so a circadian modulation is present in that variant, which the criterion does not use.
edfig1-sensor-decay¶
The raw sensor current decays slowly over days after implantation, while the magnitude of its fluctuations (delta current) stays stable.
Per placement, daily mean current and daily mean absolute delta current against day since the first reading (claims file). A placement is a run of consecutive readings lasting at least 3 days (a standard choice); slopes are compared as fractions of the first-day value.
rows = []
for s in sorted(SUBJ):
glu = SUBJ[s]["glu"]
for r, g in glu.groupby("run"):
if len(g) < 3 * 288:
continue
day = ((g.t - g.t.iloc[0]) // 86400).astype(int)
ad = g.nA.diff().abs()
dm = g.nA.groupby(day).mean(); da = ad.groupby(day).mean()
full = g.groupby(day).size() >= 144
dm, da = dm[full], da[full]
sl_c = np.polyfit(dm.index, dm.values, 1)[0]; sl_d = np.polyfit(da.index, da.values, 1)[0]
rows.append({"subject": s, "run": r, "days": len(dm), "first_day_nA": dm.iloc[0], "slope_nA_per_day": sl_c,
"rel_slope_current": sl_c / dm.iloc[0], "first_day_abs_delta": da.iloc[0], "rel_slope_abs_delta": sl_d / da.iloc[0]})
dec = pd.DataFrame(rows)
print(dec.round(3).to_string())
print(f"placements: {len(dec)} from {dec.subject.nunique()} subjects; current slope negative in {int((dec.slope_nA_per_day < 0).sum())}; "
f"|relative slope of |delta|| smaller than that of the current in {int((dec.rel_slope_abs_delta.abs() < dec.rel_slope_current.abs()).sum())}")
print(f"median relative slope per day: current {dec.rel_slope_current.median():.3f}, |delta| {dec.rel_slope_abs_delta.median():.3f}")
subject run days first_day_nA slope_nA_per_day rel_slope_current first_day_abs_delta rel_slope_abs_delta 0 CGM1 0 4 46.672 -4.338 -0.093 1.579 -0.188 1 CGM1 2 3 36.310 -1.386 -0.038 1.411 -0.304 2 CGM2 0 3 39.819 -2.489 -0.063 1.556 -0.262 3 CGM2 4 4 30.573 -5.467 -0.179 1.291 -0.081 4 CGM2 6 5 33.531 -7.725 -0.230 1.614 -0.106 5 CGM3 2 6 36.579 -7.978 -0.218 1.295 -0.037 6 CGM30 2 4 37.023 -6.999 -0.189 0.323 0.218 7 CGM31 0 7 42.638 -4.573 -0.107 0.469 0.020 8 CGM31 2 4 43.976 -13.529 -0.308 0.268 0.674 9 CGM31 4 3 42.523 -7.023 -0.165 0.383 0.003 10 CGM32 8 7 38.716 -0.800 -0.021 0.344 0.079 11 CGM36 0 5 29.727 -3.454 -0.116 0.350 0.016 12 CGM39 0 4 42.881 -4.318 -0.101 0.345 -0.084 13 CGM40 0 5 43.573 -2.096 -0.048 0.420 -0.047 14 CGM41 0 4 44.116 -6.767 -0.153 0.448 -0.123 15 CGM41 2 7 45.282 -7.353 -0.162 0.413 0.052 16 CGM46 0 7 43.828 -2.351 -0.054 0.540 -0.058 17 CGM48 0 5 37.689 -5.295 -0.140 0.382 -0.088 18 CGM48 3 4 39.322 -11.183 -0.284 0.277 0.742 19 CGM49 0 7 40.327 -2.552 -0.063 0.425 -0.007 20 CGM49 2 4 32.282 -1.613 -0.050 0.421 1.101 21 CGM52 0 6 39.869 -7.488 -0.188 0.413 -0.100 22 CGM55 0 4 42.672 -3.383 -0.079 0.307 0.455 23 CGM57 0 6 40.342 -2.156 -0.053 0.329 0.130 24 CGM58 0 7 40.339 -7.977 -0.198 0.339 -0.152 25 Vanessa 1 7 35.659 -6.207 -0.174 0.597 -0.158 26 bruce 0 3 40.724 -8.970 -0.220 0.355 0.014 27 dt15 0 5 44.095 -2.468 -0.056 0.425 -0.055 28 flex1 0 6 36.028 -1.797 -0.050 0.409 -0.036 29 flex1 2 4 38.341 -2.938 -0.077 0.283 -0.087 30 ros 0 6 46.933 -3.419 -0.073 0.626 -0.080 placements: 31 from 22 subjects; current slope negative in 31; |relative slope of |delta|| smaller than that of the current in 18 median relative slope per day: current -0.107, |delta| -0.047
Result: reproduced. The daily mean current falls in all 31 placements (median -10.7% of the first-day value per day), and the daily mean absolute delta current changes less in relative terms (median -4.7% per day; smaller in magnitude than the current's in 18 of 31 placements).
supp-ripple-rate¶
CA1 SPW-Rs occur at 9.5 (+- 8.2) per minute on average, about 14,000 per day, ranging from none for several minutes to more than 40 in a minute.
SPW-R count in 1-min bins over all CA1 recording (minutes fully covered by a file with a CA1 SPW-R table). Readings for the rats: (a) every rat with CA1 SPW-Rs and glucose together, (b) the ripple and glucose cohort, of which only CGM1 is labelled in the files. All subjects with a CA1 table are also shown.
def minute_counts(d):
out = []
for x, y in d["iv"]:
edges = np.arange(np.ceil(x / 60) * 60, y + 1e-9, 60)
if len(edges) < 2:
continue
e = d["ev"].t.values
out.append(np.diff(np.searchsorted(e, edges)))
return np.concatenate(out) if out else np.array([], int)
MC = {s: minute_counts(SUBJ[s]) for s in SUBJ if SUBJ[s]["ev"] is not None}
rows = [{"subject": s, "minutes": len(c), "mean_per_min": c.mean(), "sd_per_min": c.std(), "max_per_min": c.max(),
"per_day": c.mean() * 1440, "zero_runs_ge_5min": int(np.sum(np.diff(np.flatnonzero(np.diff(np.r_[1, (c == 0).astype(int), 1]))) [::2] >= 5))}
for s, c in MC.items()]
mr = pd.DataFrame(rows)
print(mr.round(2).to_string())
for name, keep in (("a: rats with SPW-Rs and glucose", RATS), ("b: CGM1", ["CGM1"]), ("all subjects with a CA1 table", list(MC))):
c = np.concatenate([MC[s] for s in keep])
print(f"{name}: {len(keep)} rats, {len(c) / 60:.0f} h; mean {c.mean():.2f} (sd {c.std():.2f}) per minute, {c.mean() * 1440:.0f} per day, "
f"max {c.max()} per minute; mean of rat means {mr[mr.subject.isin(keep)].mean_per_min.mean():.2f}")
c5 = np.concatenate([minute_counts({"iv": SUBJ[s]["iv"], "ev": SUBJ[s]["ev_all"]}) for s in keep])
print(f" variant, all stored events (no 5 s.d. threshold): mean {c5.mean():.2f} per minute, {c5.mean() * 1440:.0f} per day, max {c5.max()}")
print("minutes with more than 40 SPW-Rs:", {s: int((c > 40).sum()) for s, c in MC.items()})
subject minutes mean_per_min sd_per_min max_per_min per_day zero_runs_ge_5min
0 CGM1 15404 10.76 10.77 57 15494.40 678
1 CGM2 3493 9.38 12.79 70 13509.53 148
2 CGM3 21002 15.55 16.07 110 22386.03 839
3 CGM30 10435 8.18 7.59 50 11775.98 441
4 CGM31 19803 9.98 11.25 232 14367.57 739
5 CGM32 17808 8.67 9.51 258 12483.80 794
6 CGM37 890 8.43 8.36 66 12144.54 31
7 CGM39 6387 4.71 5.14 42 6776.36 296
8 CGM40 4497 11.67 9.63 81 16810.25 145
9 CGM41 17808 9.09 8.12 104 13093.26 698
10 CGM46 2989 12.92 10.18 48 18598.11 133
11 CGM48 206 7.25 6.90 38 10443.50 5
12 CGM49 1701 12.23 9.89 45 17615.24 62
13 CGM51 1496 6.94 8.17 51 10000.11 46
14 CGM52 2579 10.28 11.86 86 14799.22 76
15 CGM55 1525 11.80 9.57 43 16998.61 65
16 CGM57 1129 12.59 9.29 88 18132.01 33
17 CGM58 5021 5.85 6.79 45 8430.35 212
18 CGM60 669 10.45 7.82 40 15050.04 22
19 Vanessa 2569 6.97 6.56 45 10033.48 126
20 dt15 1726 5.62 8.28 38 8090.20 25
21 flex1 1003 10.95 8.57 35 15763.91 43
a: rats with SPW-Rs and glucose: 18 rats, 2264 h; mean 10.19 (sd 11.15) per minute, 14671 per day, max 258 per minute; mean of rat means 9.71
variant, all stored events (no 5 s.d. threshold): mean 11.07 per minute, 15941 per day, max 352
b: CGM1: 1 rats, 257 h; mean 10.76 (sd 10.77) per minute, 15494 per day, max 57 per minute; mean of rat means 10.76
variant, all stored events (no 5 s.d. threshold): mean 10.76 per minute, 15494 per day, max 57
all subjects with a CA1 table: 22 rats, 2336 h; mean 10.13 (sd 11.08) per minute, 14584 per day, max 258 per minute; mean of rat means 9.56
variant, all stored events (no 5 s.d. threshold): mean 11.00 per minute, 15844 per day, max 352
minutes with more than 40 SPW-Rs: {'CGM1': 101, 'CGM2': 115, 'CGM3': 1453, 'CGM30': 3, 'CGM31': 139, 'CGM32': 47, 'CGM37': 6, 'CGM39': 1, 'CGM40': 29, 'CGM41': 26, 'CGM46': 17, 'CGM48': 0, 'CGM49': 4, 'CGM51': 8, 'CGM52': 64, 'CGM55': 5, 'CGM57': 9, 'CGM58': 4, 'CGM60': 0, 'Vanessa': 1, 'dt15': 0, 'flex1': 0}
Result: reproduced. Over the 18 rats with SPW-Rs and glucose (2,264 h), SPW-Rs with a peak of at least 5 s.d. occur at 10.19 (s.d. 11.15) per minute, 14,671 per day, inside both tolerances, and up to 258 in a minute; minutes with more than 40 occur in 18 of 22 subjects. With all stored events the mean is 11.07 per minute and 15,941 per day, outside the tolerances; CGM1 alone gives 10.76 per minute and 15,494 per day.
supp-ripples-stable-across-days¶
In each rat the shape, amplitude and rate of SPW-Rs are consistent across days.
Per rat and calendar day (days with at least 6 h of CA1 recording, a standard choice), the median
amplitude (peak_normed_power, the Methods' maximum normalized squared signal; peak_amplitude
shown as a variant), median duration and mean rate; Spearman correlation with day, and whether
every daily value stays within 20% of the rat's overall median (claims file).
from itertools import permutations
def spearman_exact_p(y):
"""Two-sided exact p of Spearman's rho with day order, over all orderings (up to 8 days)."""
n = len(y)
if n > 8:
return np.nan
r = stats.rankdata(y); r = r - r.mean()
d = np.arange(n) - (n - 1) / 2
null = np.abs(r[np.array(list(permutations(range(n))))] @ d)
return float(np.mean(null >= abs(r @ d) - 1e-9))
rows = []
for s, d in SUBJ.items():
if d["ev"] is None:
continue
e = d["ev"].copy(); e["day"] = S.clock(d["ref"], e.t).normalize()
secs = pd.Series(np.concatenate([np.arange(x, y, 60) for x, y in d["iv"]]))
mins_per_day = S.clock(d["ref"], secs).normalize().value_counts()
good = mins_per_day[mins_per_day >= 360].index
g = e[e.day.isin(good)].groupby("day")
daily = pd.DataFrame({"amp": g.peak_normed_power.median(), "amp_uV": g.peak_amplitude.median(),
"dur": g.duration.median(), "rate": g.size() / mins_per_day[good].reindex(g.size().index)})
if len(daily) < 3:
continue
row = {"subject": s, "days": len(daily)}
for k in daily:
rho = stats.spearmanr(np.arange(len(daily)), daily[k])
med = e[k if k not in ("amp", "amp_uV", "dur") else {"amp": "peak_normed_power", "amp_uV": "peak_amplitude", "dur": "duration"}[k]].median() if k != "rate" else daily[k].median()
row[f"{k}_rho"] = rho.statistic; row[f"{k}_p"] = rho.pvalue; row[f"{k}_p_exact"] = spearman_exact_p(daily[k].values)
row[f"{k}_within20"] = bool(np.all(np.abs(daily[k] / med - 1) <= 0.2))
row[f"{k}_maxdev"] = float(np.max(np.abs(daily[k] / med - 1)))
rows.append(row)
st = pd.DataFrame(rows)
print(st.round(3).to_string())
for k in ("amp", "amp_uV", "dur", "rate"):
ok = (st[f"{k}_p"] >= 0.05) | st[f"{k}_within20"]
print(f"{k}: consistent (Spearman p >= 0.05 or within 20%) in {int(ok.sum())} of {len(st)} rats; "
f"within 20% in {int(st[f'{k}_within20'].sum())}; Spearman p < 0.05 in {int((st[f'{k}_p'] < 0.05).sum())} "
f"(exact p < 0.05 in {int((st[f'{k}_p_exact'] < 0.05).sum())} of {int(st[f'{k}_p_exact'].notna().sum())} rats with up to 8 days)")
print("exact two-sided p of a perfectly monotone series, by number of days:", {n: float(f"{spearman_exact_p(np.arange(n)):.2g}") for n in range(3, 9)})
subject days amp_rho amp_p amp_p_exact amp_within20 amp_maxdev amp_uV_rho amp_uV_p amp_uV_p_exact amp_uV_within20 amp_uV_maxdev dur_rho dur_p dur_p_exact dur_within20 dur_maxdev rate_rho rate_p rate_p_exact rate_within20 rate_maxdev
0 CGM1 12 0.420 0.175 NaN False 0.326 0.755 0.005 NaN False 0.848 0.436 0.157 NaN True 0.089 0.000 1.000 NaN False 0.933
1 CGM2 5 0.700 0.188 0.233 False 0.558 0.600 0.285 0.350 False 0.401 0.667 0.219 0.267 False 0.341 0.900 0.037 0.083 False 0.974
2 CGM3 17 0.086 0.743 NaN True 0.137 0.316 0.216 NaN True 0.083 0.180 0.490 NaN True 0.078 -0.081 0.758 NaN False 0.560
3 CGM30 11 -0.164 0.631 NaN False 0.492 -0.500 0.117 NaN False inf 0.436 0.180 NaN True 0.178 -0.318 0.340 NaN False 0.992
4 CGM31 18 -0.123 0.627 NaN False 0.465 -0.748 0.000 NaN True 0.195 -0.422 0.081 NaN True 0.136 -0.292 0.240 NaN False 0.340
5 CGM32 15 0.239 0.390 NaN False 0.318 0.342 0.212 NaN False 1.000 -0.121 0.669 NaN True 0.065 -0.468 0.079 NaN False 0.268
6 CGM39 7 -0.321 0.482 0.498 False 0.967 -0.571 0.180 0.200 False 1.528 -0.714 0.071 0.088 False 0.225 -0.571 0.180 0.200 False 0.964
7 CGM40 5 0.700 0.188 0.233 False 0.372 0.300 0.624 0.683 False 0.587 -0.300 0.624 0.683 True 0.133 0.600 0.285 0.350 True 0.174
8 CGM41 17 -0.113 0.667 NaN False 0.420 -0.345 0.176 NaN False 1.000 0.079 0.764 NaN True 0.093 -0.458 0.064 NaN False 0.538
9 CGM46 3 0.500 0.667 1.000 True 0.162 1.000 0.000 0.333 False 0.205 0.500 0.667 1.000 True 0.067 0.500 0.667 1.000 True 0.019
10 CGM49 3 0.500 0.667 1.000 False 0.348 0.500 0.667 1.000 True 0.050 0.500 0.667 1.000 True 0.067 1.000 0.000 0.333 True 0.094
11 CGM51 3 -0.500 0.667 1.000 True 0.170 -1.000 0.000 0.333 False 0.351 0.500 0.667 1.000 True 0.070 -1.000 0.000 0.333 True 0.118
12 CGM52 4 0.600 0.400 0.417 True 0.193 -0.600 0.400 0.417 False 0.203 0.632 0.368 0.500 True 0.044 -1.000 0.000 0.083 False 0.308
13 CGM58 7 0.214 0.645 0.662 False 0.352 0.286 0.535 0.556 False 1.185 0.072 0.878 0.889 True 0.119 -0.036 0.939 0.963 False 1.048
14 Vanessa 4 -0.600 0.400 0.417 True 0.095 -1.000 0.000 0.083 False 0.202 0.800 0.200 0.333 True 0.077 -0.800 0.200 0.333 False 0.291
amp: consistent (Spearman p >= 0.05 or within 20%) in 15 of 15 rats; within 20% in 5; Spearman p < 0.05 in 0 (exact p < 0.05 in 0 of 9 rats with up to 8 days)
amp_uV: consistent (Spearman p >= 0.05 or within 20%) in 11 of 15 rats; within 20% in 3; Spearman p < 0.05 in 5 (exact p < 0.05 in 0 of 9 rats with up to 8 days)
dur: consistent (Spearman p >= 0.05 or within 20%) in 15 of 15 rats; within 20% in 13; Spearman p < 0.05 in 0 (exact p < 0.05 in 0 of 9 rats with up to 8 days)
rate: consistent (Spearman p >= 0.05 or within 20%) in 13 of 15 rats; within 20% in 4; Spearman p < 0.05 in 4 (exact p < 0.05 in 0 of 9 rats with up to 8 days)
exact two-sided p of a perfectly monotone series, by number of days: {3: 0.33, 4: 0.083, 5: 0.017, 6: 0.0028, 7: 0.0004, 8: 5e-05}
Result: reproduced. Daily median amplitude (peak normed power) and duration show no trend
with day (no Spearman p < 0.05) in all 15 rats with at least 3 days, and the daily rate in 11 of 15
(13 counting rats whose rate stays within 20%). Amplitude stays within 20% of the median in only 5
rats, and the rate in 4, so the consistency rests mostly on the absence of trends. Amplitude in
microvolts trends with day in 5 rats. These p values are scipy's asymptotic ones, which are too small with
few days: every perfectly monotone series of 3 or 4 days gets p = 0.000, where the exact two-sided
p is 0.33 for 3 days and 0.083 for 4, so with 3 or 4 days no series can reach p < 0.05. With exact
p no rat with up to 8 days has a trend in any measure: the two rate trends that fail the criterion
(CGM2, 5 days, rho 0.9, exact p 0.083; CGM52, 4 days, exact p 0.083) and the microvolt trends of
CGM46, CGM51 and Vanessa are artefacts of the asymptotic p. With exact p the rate is consistent in
all 15 rats, and amplitude in microvolts trends only in CGM1 and CGM31 (12 and 18 days). In CGM30 the median peak_amplitude over all
events is 0, which makes its deviation infinite; this variant does not enter the criterion.
edfig3-ripple-rate-rhythms¶
SPW-R rate varies with circadian and ultradian rhythms, with spectral peaks at about 24 h and about 4 h, in 1,568.5 recording hours from 8 rats.
Power spectral density of SPW-R rate in 5-min bins per rat, averaged over rats (claims file). The recordings have gaps, so the spectrum is a Lomb-Scargle periodogram of the covered bins (a standard choice for unevenly sampled data). The 8 rats with the most CA1 recording are used, and all rats with at least 48 h as a variant.
per_h = np.linspace(1, 48, 400)
freqs = 2 * np.pi / (per_h * 3600)
LS, hrs = {}, {}
for s, d in SUBJ.items():
if d["ev"] is None:
continue
t_all, r_all = [], []
for x, y in d["iv"]:
edges = np.arange(x, y - 299, 300)
if len(edges) < 2:
continue
c = np.diff(np.searchsorted(d["ev"].t.values, np.r_[edges, edges[-1] + 300]))
t_all.append(edges + 150); r_all.append(c / 5)
t_all = np.concatenate(t_all); r_all = np.concatenate(r_all)
hrs[s] = len(t_all) * 5 / 60
if hrs[s] >= 48:
p = signal.lombscargle(t_all, r_all - r_all.mean(), freqs, normalize=True)
LS[s] = p
top8 = sorted(hrs, key=hrs.get, reverse=True)[:8]
print("CA1 recording hours per subject:", {k: round(v, 1) for k, v in sorted(hrs.items(), key=lambda x: -x[1])})
for name, keep in (("8 rats with the most recording", top8), ("all rats with >= 48 h", list(LS))):
m = np.mean([LS[s] for s in keep], axis=0)
pk, _ = signal.find_peaks(m)
order = pk[np.argsort(m[pk])[::-1]]
print(f"{name}: {len(keep)} rats, {sum(hrs[s] for s in keep):.1f} h; local peaks (period h, power): "
f"{[(round(per_h[i], 1), round(m[i], 4)) for i in order[:8]]}")
prom = signal.peak_prominences(m, pk)[0]
rank = pd.Series(prom).rank(ascending=False).values
for target in (24, 4):
k = np.argmin(np.abs(per_h[pk] - target))
print(f" peak nearest {target} h at {per_h[pk][k]:.1f} h: power {m[pk][k]:.4f}, prominence {prom[k]:.4f}, "
f"rank {int(rank[k])} of {len(pk)} peaks by prominence (largest prominence {prom.max():.4f} at {per_h[pk][np.argmax(prom)]:.1f} h)")
print(f" peak nearest 24 h: {per_h[pk][np.argmin(np.abs(per_h[pk] - 24))]:.1f} h; nearest 4 h: {per_h[pk][np.argmin(np.abs(per_h[pk] - 4))]:.1f} h")
fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(per_h, np.mean([LS[s] for s in top8], axis=0)); ax.set_xlabel("period (h)"); ax.set_ylabel("normalized power"); plt.show()
CA1 recording hours per subject: {'CGM3': 349.2, 'CGM31': 329.2, 'CGM32': 295.9, 'CGM41': 295.8, 'CGM1': 256.2, 'CGM30': 173.6, 'CGM39': 106.2, 'CGM58': 83.5, 'CGM40': 74.7, 'CGM2': 58.0, 'CGM46': 49.8, 'Vanessa': 42.8, 'CGM52': 42.6, 'dt15': 28.6, 'CGM49': 28.2, 'CGM55': 25.2, 'CGM51': 24.7, 'CGM57': 18.7, 'flex1': 16.7, 'CGM37': 14.6, 'CGM60': 11.1, 'CGM48': 3.4}
8 rats with the most recording: 8 rats, 1889.8 h; local peaks (period h, power): [(np.float64(42.3), np.float64(0.0225)), (np.float64(19.5), np.float64(0.0185)), (np.float64(24.0), np.float64(0.016)), (np.float64(20.0), np.float64(0.0152)), (np.float64(16.3), np.float64(0.015)), (np.float64(34.7), np.float64(0.0142)), (np.float64(28.8), np.float64(0.0117)), (np.float64(29.9), np.float64(0.0113))]
peak nearest 24 h at 24.0 h: power 0.0160, prominence 0.0107, rank 4 of 45 peaks by prominence (largest prominence 0.0175 at 42.3 h)
peak nearest 4 h at 3.8 h: power 0.0076, prominence 0.0034, rank 12 of 45 peaks by prominence (largest prominence 0.0175 at 42.3 h)
peak nearest 24 h: 24.0 h; nearest 4 h: 3.8 h
all rats with >= 48 h: 11 rats, 2072.2 h; local peaks (period h, power): [(np.float64(19.5), np.float64(0.0218)), (np.float64(42.0), np.float64(0.0215)), (np.float64(40.9), np.float64(0.0214)), (np.float64(20.0), np.float64(0.0193)), (np.float64(35.9), np.float64(0.0188)), (np.float64(38.0), np.float64(0.0176)), (np.float64(24.0), np.float64(0.0166)), (np.float64(16.3), np.float64(0.0141))]
peak nearest 24 h at 24.0 h: power 0.0166, prominence 0.0084, rank 7 of 45 peaks by prominence (largest prominence 0.0143 at 19.5 h)
peak nearest 4 h at 4.1 h: power 0.0096, prominence 0.0014, rank 21 of 45 peaks by prominence (largest prominence 0.0143 at 19.5 h)
peak nearest 24 h: 24.0 h; nearest 4 h: 4.1 h
Result: partially reproduced. The mean Lomb-Scargle spectrum of the 8 rats with the most recording has local peaks at 24.0 h and 3.8 h, both inside the tolerances, but they are not its largest peaks (42.3 and 19.5 h are larger), and the 8 rats hold 1,889.8 h of recording, outside 1,568.5 +- 156.9 h. By prominence the 24.0-h peak ranks 4th of 45 local peaks and the 3.8-h peak 12th (prominence 0.0034 against 0.0175 for the 42.3-h peak), so on this 400-point grid a nearby local peak is weak evidence of a 4-h rhythm.
edfig2-refractory-period¶
Inter-ripple intervals drop to zero below about 20 ms, a refractory period between ripples, with peaks reflecting bursts of ripples.
Histogram of intervals between consecutive SPW-R peaks within a file, per rat, 1-ms bins (claims file); the shortest interval with counts. The legend does not name the time point of each event, so start-to-start and stop-to-stop intervals and the gap from one event's stop to the next one's start are shown as variants. The Methods' detector merges neighbouring events ("identify peaks above a threshold merging neighbouring events", SPW-R detection), which sets a floor on the gap.
rows = []
for s, d in SUBJ.items():
if d["ev"] is None:
continue
iri = d["ev"].groupby("file").t.diff().dropna().values * 1000
h, _ = np.histogram(iri, bins=np.arange(0, 501))
nz = np.flatnonzero(h)
gap = d["ev"].t.values[1:] - d["ev"].t.values[:-1]
rows.append({"subject": s, "intervals": len(iri), "shortest_ms": iri.min(), "first_bin_with_counts_ms": nz[0] if len(nz) else np.nan,
"count_below_20ms": int((iri < 20).sum()), "count_20_30ms": int(((iri >= 20) & (iri < 30)).sum()),
"histogram_peak_ms": int(np.argmax(h)), "within_18_22": 18 <= iri.min() <= 22})
rf = pd.DataFrame(rows)
print(rf.round(2).to_string())
print(f"shortest interval within 18-22 ms in {int(rf.within_18_22.sum())} of {len(rf)} rats; median shortest {rf.shortest_ms.median():.1f} ms")
# other time points of each event, and the gap from one event's stop to the next one's start
rows = []
for s, d in SUBJ.items():
if d["ev"] is None:
continue
e = d["ev"].sort_values(["file", "start_time"])
st = e.groupby("file").start_time.diff().dropna().values * 1000
gap = (e.start_time - e.groupby("file").stop_time.shift(1)).dropna().values * 1000
sp = e.groupby("file").stop_time.diff().dropna().values * 1000
rows.append({"subject": s, "shortest_start_to_start_ms": st.min(), "shortest_stop_to_stop_ms": sp.min(),
"shortest_stop_to_next_start_ms": gap.min(), "gap_1st_percentile_ms": np.percentile(gap, 1),
"gaps_below_10ms": int((gap < 10).sum()), "gaps_below_20ms": int((gap < 20).sum()),
"start_to_start_18_22": 18 <= st.min() <= 22, "gap_18_22": 18 <= gap.min() <= 22})
rg = pd.DataFrame(rows)
print(rg.round(1).to_string())
print(f"shortest start-to-start interval within 18-22 ms in {int(rg.start_to_start_18_22.sum())} of {len(rg)} rats "
f"(median {rg.shortest_start_to_start_ms.median():.1f} ms); shortest stop-to-next-start gap within 18-22 ms in {int(rg.gap_18_22.sum())} "
f"(median {rg.shortest_stop_to_next_start_ms.median():.1f} ms, range {rg.shortest_stop_to_next_start_ms.min():.1f} to {rg.shortest_stop_to_next_start_ms.max():.1f})")
subject intervals shortest_ms first_bin_with_counts_ms count_below_20ms count_20_30ms histogram_peak_ms within_18_22 0 CGM1 165908 30.4 30 0 0 76 False 1 CGM2 32799 30.4 30 0 0 87 False 2 CGM3 286345 28.8 28 0 5 76 False 3 CGM30 85361 35.2 35 0 0 87 False 4 CGM31 197682 31.2 31 0 0 87 False 5 CGM32 154684 34.4 34 0 0 87 False 6 CGM37 7552 28.8 28 0 1 103 False 7 CGM39 30113 36.8 36 0 0 87 False 8 CGM40 52510 36.8 36 0 0 80 False 9 CGM41 162115 34.4 34 0 0 80 False 10 CGM46 38628 32.8 32 0 0 87 False 11 CGM48 1501 32.0 32 0 0 71 False 12 CGM49 20848 38.4 38 0 0 87 False 13 CGM51 10415 29.6 29 0 1 75 False 14 CGM52 26550 30.4 30 0 0 72 False 15 CGM55 18006 39.2 39 0 0 76 False 16 CGM57 14230 31.2 31 0 0 80 False 17 CGM58 29409 36.0 36 0 0 92 False 18 CGM60 6996 36.8 36 0 0 108 False 19 Vanessa 17925 34.4 34 0 0 67 False 20 dt15 9733 36.8 36 0 0 91 False 21 flex1 10984 35.2 35 0 0 87 False shortest interval within 18-22 ms in 0 of 22 rats; median shortest 34.4 ms
subject shortest_start_to_start_ms shortest_stop_to_stop_ms shortest_stop_to_next_start_ms gap_1st_percentile_ms gaps_below_10ms gaps_below_20ms start_to_start_18_22 gap_18_22 0 CGM1 36.0 35.2 15.2 22.4 0 1080 False False 1 CGM2 37.6 35.2 15.2 20.8 0 279 False False 2 CGM3 35.2 35.2 15.2 16.8 0 7476 False False 3 CGM30 40.8 40.0 20.0 28.0 0 38 False True 4 CGM31 40.0 40.0 20.0 25.6 0 62 False True 5 CGM32 40.0 40.0 20.0 23.2 0 53 False True 6 CGM37 40.8 42.4 20.0 27.2 0 1 False True 7 CGM39 42.4 40.8 20.0 31.2 0 1 False True 8 CGM40 40.0 40.8 20.0 23.2 0 3 False True 9 CGM41 40.0 40.0 20.0 27.2 0 33 False True 10 CGM46 40.8 40.0 20.0 25.6 0 1 False True 11 CGM48 40.8 44.0 20.0 24.0 0 0 False True 12 CGM49 40.8 42.4 20.0 27.2 0 7 False True 13 CGM51 42.4 40.8 20.0 26.4 0 1 False True 14 CGM52 40.8 40.8 20.0 24.0 0 10 False True 15 CGM55 42.4 41.6 20.0 25.6 0 6 False True 16 CGM57 40.8 40.0 20.0 24.0 0 25 False True 17 CGM58 40.0 40.8 20.0 30.4 0 15 False True 18 CGM60 43.2 40.8 20.0 27.2 0 3 False True 19 Vanessa 40.8 40.8 20.0 25.6 0 7 False True 20 dt15 40.8 40.8 20.0 26.4 0 3 False True 21 flex1 40.8 40.8 20.0 24.8 0 2 False True shortest start-to-start interval within 18-22 ms in 0 of 22 rats (median 40.8 ms); shortest stop-to-next-start gap within 18-22 ms in 19 (median 20.0 ms, range 15.2 to 20.0)
Result: not reproduced. Between SPW-R peaks, the reading the criterion uses, the shortest interval is 28.8 to 39.2 ms (median 34.4 ms) and no interval below 20 ms occurs in any of 22 rats; no rat has its shortest interval in 18 to 22 ms. Start-to-start (median shortest 40.8 ms) does not match either. The gap from one event's stop to the next one's start has its minimum at exactly 20.0 ms in 19 rats and at 15.2 ms in CGM1, CGM2 and CGM3, so it lies in 18 to 22 ms in 19 of 22 rats. A floor at the same value in nearly every rat matches the merging of neighbouring events in the detection, so the 20-ms drop the legend describes may be this rule, measured as stop to start.
edfig4-nrem-wake¶
The SPW-R and glucose cross-correlogram does not differ between NREM sleep and waking, and its negative peak does not change across state transitions.
Sleep states are stored for six subjects; three of them (CGM1, CGM2, CGM3) have SPW-Rs and glucose together. As the claims file: blocks of 4 h; within each, the correlation at +10 min (Spearman, the correlation type under which the main claim matched) between the SPW-R rate of bins mostly in NREM (or mostly awake) and the delta current 10 min later; one value per block and state, unpaired two-sided t-test (per-rat means, paired, as a variant). Transitions: for each change from an awake to a NREM episode (and the reverse) where both episodes last at least 10 min (a choice), the most negative correlation from +5 to +20 min in the 2 h before and the 2 h after, paired t-test.
The main analysis pairs the SPW-R count and the glucose change of the same 5-min interval, as in fig2-ripple-rate-predicts-glucose (an assumption). The two alignments one reading apart are shown as variants: counting SPW-Rs over the 5 min after each reading, or taking the change to the next reading, which is the same as moving the offset by one reading.
def state_of_bins(b, ss):
lab = np.full(len(b), "", dtype=object)
for state in ("Non-REM", "Awake"):
g = ss[ss.label == state]
cov_ = S.coverage(list(zip(g.start, g.stop)), b.t.values - S.STEP, b.t.values)
lab[cov_ > 0.5] = state
return lab
def nrem_wake(k):
blocks, trans = [], []
for s in RATS:
fs = SUBJ[s]["files"]
if not any(f["sleep"] is not None and len(f["sleep"]) for f in fs):
continue
ss = pd.concat([f["sleep"] + [S.clock_offset(f, SUBJ[s]["ref"]), S.clock_offset(f, SUBJ[s]["ref"]), ""]
for f in fs if f["sleep"] is not None]).sort_values("start").reset_index(drop=True)
b = SUBJ[s]["b"].copy()
b["state"] = state_of_bins(b, ss)
b["block"] = ((b.t - b.t.min()) // (4 * 3600)).astype(int)
for blk, g in b.groupby("block"):
for state in ("Non-REM", "Awake"):
gg = g.copy(); gg.loc[gg.state != state, "rate"] = np.nan
x, y = S.lagged_pairs(gg, 2 + k)
if len(x) >= 10:
blocks.append({"rat": s, "block": blk, "state": state, "pairs": len(x), "R10": S._corr(x, y, "spearman")})
# transitions between episodes of at least 10 min
ss2 = ss[ss.label.isin(["Non-REM", "Awake"])].reset_index(drop=True)
for i in range(1, len(ss2)):
a, c = ss2.iloc[i - 1], ss2.iloc[i]
if a.label == c.label or a.stop - a.start < 600 or c.stop - c.start < 600 or abs(c.start - a.stop) > 60:
continue
t0 = c.start
vals = {}
for side, (lo, hi) in (("before", (t0 - 7200, t0)), ("after", (t0, t0 + 7200))):
g = b[(b.t > lo) & (b.t <= hi)]
cc = [S._corr(*S.lagged_pairs(g, j + k), "spearman") for j in (1, 2, 3, 4)]
vals[side] = np.nanmin(cc) if np.any(np.isfinite(cc)) else np.nan
trans.append({"rat": s, "type": f"{a.label}->{c.label}", **vals})
return pd.DataFrame(blocks), pd.DataFrame(trans)
blocks, trans = nrem_wake(0)
nrem = blocks[blocks.state == "Non-REM"].R10; wake = blocks[blocks.state == "Awake"].R10
nrem, wake = nrem.dropna(), wake.dropna()
t = stats.ttest_ind(nrem, wake)
print(f"blocks: NREM n = {len(nrem)}, mean R10 {nrem.mean():.3f}; awake n = {len(wake)}, mean {wake.mean():.3f}; unpaired t = {t.statistic:.2f}, p = {t.pvalue:.3g}")
pm = blocks.groupby(["rat", "state"]).R10.mean().unstack()
print(pm.round(3).to_string())
if len(pm) > 1:
tp = stats.ttest_rel(pm["Non-REM"], pm["Awake"], nan_policy="omit"); print(f"per-rat means, paired: t = {tp.statistic:.2f}, p = {tp.pvalue:.3g}")
for typ, g in trans.dropna().groupby("type"):
tt = stats.ttest_rel(g.before, g.after)
print(f"{typ}: {len(g)} transitions from {g.rat.nunique()} rats; most negative R (+5 to +20 min) before {g.before.mean():.3f}, after {g.after.mean():.3f}; paired t = {tt.statistic:.2f}, p = {tt.pvalue:.3g}")
print("Alignment variants:")
for name, k in ALIGN.items():
if k == 0:
continue
bl, tr = nrem_wake(k)
n_, w_ = bl[bl.state == "Non-REM"].R10.dropna(), bl[bl.state == "Awake"].R10.dropna()
txt = "; ".join(f"{typ}: {len(g)}, p = {stats.ttest_rel(g.before, g.after).pvalue:.3g}" for typ, g in tr.dropna().groupby("type"))
print(f" {name}: blocks NREM mean {n_.mean():.3f} (n = {len(n_)}), awake {w_.mean():.3f} (n = {len(w_)}), unpaired p = {stats.ttest_ind(n_, w_).pvalue:.3g}; transitions {txt}")
blocks: NREM n = 120, mean R10 -0.017; awake n = 112, mean -0.083; unpaired t = 1.92, p = 0.0567 state Awake Non-REM rat CGM1 -0.160 -0.050 CGM2 -0.136 0.159 CGM3 -0.022 -0.028 per-rat means, paired: t = 1.52, p = 0.268 Awake->Non-REM: 13 transitions from 2 rats; most negative R (+5 to +20 min) before -0.261, after -0.151; paired t = -1.55, p = 0.148 Non-REM->Awake: 23 transitions from 1 rats; most negative R (+5 to +20 min) before -0.234, after -0.165; paired t = -1.23, p = 0.23 Alignment variants:
SPW-Rs after the reading: blocks NREM mean -0.043 (n = 120), awake -0.109 (n = 114), unpaired p = 0.0378; transitions Awake->Non-REM: 13, p = 0.0843; Non-REM->Awake: 23, p = 0.28
delta to the next reading: blocks NREM mean -0.001 (n = 120), awake -0.040 (n = 112), unpaired p = 0.212; transitions Awake->Non-REM: 13, p = 0.271; Non-REM->Awake: 23, p = 0.658
Result: reproduced, on 3 rats where the paper has 8, resting on the alignment assumption. NREM and awake blocks do not differ at p < 0.05 (mean R at +10 min -0.017 against -0.083, unpaired t p = 0.057; per-rat means paired p = 0.27), and the most negative correlation does not change across awake-to-NREM (p = 0.15, 13 transitions) or NREM-to-awake transitions (p = 0.23, 23 transitions from one rat). The awake blocks are more negative, near the threshold, and the test has little power with 3 rats. Under the alignments one reading apart, counting SPW-Rs over the 5 min after each reading makes the awake blocks differ from the NREM blocks (mean -0.109 against -0.043, unpaired p = 0.038), so the claim fails there; taking the change to the next reading, no test reaches p < 0.05 (blocks p = 0.21, transitions p = 0.27 and 0.66).
edfig4-no-hour-differs¶
No hour of the day has a significantly different SPW-R and glucose correlation than any other hour (smallest P = 0.53).
As the claims file: per rat, cross-correlograms in 4-h windows of the day centred on each hour (all days pooled), integrated over offsets +5 to +20 min (the sum of R at +5, +10, +15, +20 min; Spearman, as for the main claim); Student's t-test across rats for every pair of hours. Rats with at least 24 h of both signals (a standard choice, so that each hour has data). A rat-hour in which any of the four correlations cannot be computed is left out.
The main analysis pairs the SPW-R count and the glucose change of the same 5-min interval, as in fig2-ripple-rate-predicts-glucose (an assumption). The two alignments one reading apart are shown as variants: counting SPW-Rs over the 5 min after each reading, or taking the change to the next reading, which is the same as moving the offset by one reading.
rats24 = list(cov[cov.hours_with_both >= 24].subject)
def hour_matrix(k):
"""Integrated R (+5 to +20 min) per rat and hour of day, with the lags moved by k readings (the alignment)."""
H = np.full((len(rats24), 24), np.nan)
for i, s in enumerate(rats24):
b = SUBJ[s]["b"]
hod = (b.clock.dt.hour + b.clock.dt.minute / 60).values
for hcent in range(24):
dist = np.abs(((hod - hcent) + 12) % 24 - 12)
g = b.copy(); g.loc[dist > 2, "rate"] = np.nan
v = [S._corr(*S.lagged_pairs(g, j + k), "spearman") for j in (1, 2, 3, 4)]
H[i, hcent] = np.sum(v) if np.all(np.isfinite(v)) else np.nan
return H
def hour_pairs(H):
ps = []
for a in range(24):
for c in range(a + 1, 24):
ps.append((a, c, stats.ttest_ind(H[:, a], H[:, c], nan_policy="omit").pvalue,
stats.ttest_rel(H[:, a], H[:, c], nan_policy="omit").pvalue))
return pd.DataFrame(ps, columns=["hour_a", "hour_b", "p_student", "p_paired"])
H = hour_matrix(0)
ps = hour_pairs(H)
print(f"{len(rats24)} rats: {rats24}; rat-hours with no valid integrated correlation (left out): {int(np.isnan(H).sum())} of {H.size}")
print("integrated R (+5 to +20 min) by hour, mean over rats:", np.round(np.nanmean(H, axis=0), 3))
print(f"{len(ps)} pairs of hours; Student's t-test: smallest p = {ps.p_student.min():.3f}, pairs with p < 0.05: {int((ps.p_student < 0.05).sum())}; "
f"paired variant: smallest p = {ps.p_paired.min():.3g}, pairs with p < 0.05: {int((ps.p_paired < 0.05).sum())}; "
f"expected by chance at 0.05: {0.05 * len(ps):.1f}")
for name, k in ALIGN.items():
if k:
p_ = hour_pairs(hour_matrix(k))
print(f"alignment variant, {name}: Student's t smallest p = {p_.p_student.min():.3f}, pairs with p < 0.05: {int((p_.p_student < 0.05).sum())}")
12 rats: ['CGM1', 'CGM2', 'CGM3', 'CGM30', 'CGM31', 'CGM32', 'CGM39', 'CGM40', 'CGM41', 'CGM46', 'CGM58', 'Vanessa']; rat-hours with no valid integrated correlation (left out): 3 of 288 integrated R (+5 to +20 min) by hour, mean over rats: [-0.145 -0.199 -0.031 -0.036 0.017 0.048 -0.132 -0.198 -0.13 -0.223 -0.109 -0.089 -0.267 -0.241 -0.235 -0.298 -0.283 -0.374 -0.453 -0.417 -0.508 -0.357 -0.329 -0.245] 276 pairs of hours; Student's t-test: smallest p = 0.010, pairs with p < 0.05: 10; paired variant: smallest p = 0.0098, pairs with p < 0.05: 25; expected by chance at 0.05: 13.8
alignment variant, SPW-Rs after the reading: Student's t smallest p = 0.032, pairs with p < 0.05: 2
alignment variant, delta to the next reading: Student's t smallest p = 0.015, pairs with p < 0.05: 13
Result: not reproduced. Student's t-test across 12 rats finds 10 of 276 pairs of hours at p < 0.05 (smallest p = 0.010), where the paper reports a smallest P of 0.53; the paired variant finds 25. 3 of 288 rat-hours have no valid correlation and are left out. The integrated correlation is most negative from 17:00 to 21:00 (-0.37 to -0.51) and near zero from 02:00 to 05:00 (-0.04 to 0.05). Ten pairs is fewer than the 13.8 expected by chance among 276 uncorrected tests, so the difference from the paper lies in the smallest p, not in a clear effect of hour. The two other alignments give the same outcome: smallest p 0.032 (2 pairs at p < 0.05) counting SPW-Rs after each reading, and 0.015 (13 pairs) with the change to the next reading.
edfig9-dorsal-ventral¶
Dorsal and ventral CA1 SPW-Rs are only weakly correlated, and ventral SPW-Rs are significantly less correlated with glucose fluctuations than dorsal SPW-Rs, which match the earlier dorsal cohorts.
The dorsal and ventral rats are the three subjects with Dorsalripples and Ventralripples
tables (CGM37, CGM48, Vanessa). The correlation of dorsal and ventral 5-min SPW-R rates is taken
over the files holding both tables. The glucose comparison needs glucose recorded with the
electrophysiology: CGM37 has none, CGM48 has 2.2 h (too little for a day of at least 4 h, printed
below), and Vanessa has three such days, two of which give a correlation; per-day correlations at +10 min (Spearman) for dorsal
and ventral SPW-Rs, paired (a) and unpaired (b) t-tests, and the dorsal values against the other
rats' per-rat values (unpaired t-test).
rows = []
for s in ("CGM37", "CGM48", "Vanessa"):
for f in SUBJ[s]["files"]:
if "Dorsalripples" in f["tables"] and "Ventralripples" in f["tables"]:
dur = f["lfp_n"] / f["lfp_rate"]
edges = np.arange(0, dur + 1e-9, 300)
if len(edges) < 3:
continue
dc = np.histogram(f["tables"]["Dorsalripples"].peak, edges)[0]; vc = np.histogram(f["tables"]["Ventralripples"].peak, edges)[0]
rows.append({"subject": s, "file": f["path"].split("/")[1][-30:], "bins": len(dc), "dorsal": int(dc.sum()), "ventral": int(vc.sum()),
"r_5min": np.corrcoef(dc, vc)[0, 1] if dc.std() and vc.std() else np.nan, "_d": dc, "_v": vc})
dv = pd.DataFrame(rows)
print(dv.drop(columns=["_d", "_v"]).round(3).to_string())
for s, g in dv.groupby("subject"):
d_, v_ = np.concatenate(list(g._d)), np.concatenate(list(g._v))
print(f"{s}: pooled 5-min rate correlation dorsal vs ventral r = {np.corrcoef(d_, v_)[0, 1]:.3f} over {len(d_)} bins")
for s in ("CGM37", "CGM48", "Vanessa"):
ev_, iv_ = S.subject_events(SUBJ[s]["files"], tables=("Dorsalripples",), reference=SUBJ[s]["ref"])
bb_ = S.binned(SUBJ[s]["glu"], ev_, iv_)
ok_ = np.isfinite(bb_.rate) & np.isfinite(bb_.delta)
per_day = S.clock(SUBJ[s]["ref"], bb_.t[ok_]).normalize().value_counts()
print(f"{s}: {ok_.sum() * 5 / 60:.1f} h of glucose with the dorsal and ventral tables; days with at least 4 h: {int((per_day >= 48).sum())}")
d = SUBJ["Vanessa"]
res = {}
for name in ("Dorsalripples", "Ventralripples"):
ev, iv = S.subject_events(d["files"], tables=(name,), reference=d["ref"])
bb = S.binned(d["glu"], ev, iv); bb["day"] = S.clock(d["ref"], bb.t).normalize()
res[name] = {day: S._corr(*S.lagged_pairs(g, 2), "spearman") for day, g in bb.groupby("day") if (np.isfinite(g.rate) & np.isfinite(g.delta)).sum() >= 48}
pdays = pd.DataFrame(res)
print("Vanessa, per-day Spearman R at +10 min (days with at least 4 h of both):"); print(pdays.round(3).to_string())
pdays = pdays.dropna()
if len(pdays) >= 2:
print("paired t:", stats.ttest_rel(pdays.Dorsalripples, pdays.Ventralripples), " unpaired t:", stats.ttest_ind(pdays.Dorsalripples, pdays.Ventralripples))
others = [CC[(s, "raw", "spearman", "rat")][LAG_MIN == 10][0] for s in RATS if s != "Vanessa"]
print(f"other rats' R at +10 min: mean {np.mean(others):.3f} (n = {len(others)}); Vanessa dorsal per-day values {np.round(pdays.Dorsalripples.values, 3)}; "
f"unpaired t = {stats.ttest_ind(pdays.Dorsalripples, others).statistic:.2f}, p = {stats.ttest_ind(pdays.Dorsalripples, others).pvalue:.3g}")
subject file bins dorsal ventral r_5min
0 CGM37 34um_210121_113448_ecephys.nwb 28 829 949 0.237
1 CGM37 34um_210121_135742_ecephys.nwb 52 3369 2920 0.606
2 CGM37 34um_210124_090816_ecephys.nwb 23 967 1356 0.548
3 CGM37 34um_210120_163540_ecephys.nwb 26 1634 1580 0.602
4 CGM37 34um_210121_101926_ecephys.nwb 3 112 218 0.754
5 CGM37 34um_210121_103917_ecephys.nwb 10 296 683 0.804
6 CGM37 34um_210126_144503_ecephys.nwb 33 2073 2011 0.026
7 CGM48 _0um_210301_100159_ecephys.nwb 7 274 473 0.549
8 CGM48 _0um_210302_090041_ecephys.nwb 27 945 1644 0.864
9 Vanessa essa_210219_164638_ecephys.nwb 10 604 375 0.443
10 Vanessa essa_210219_174006_ecephys.nwb 186 7438 7011 0.348
11 Vanessa essa_210226_073342_ecephys.nwb 312 11713 19797 0.134
CGM37: pooled 5-min rate correlation dorsal vs ventral r = 0.439 over 175 bins
CGM48: pooled 5-min rate correlation dorsal vs ventral r = 0.834 over 34 bins
Vanessa: pooled 5-min rate correlation dorsal vs ventral r = 0.133 over 508 bins
CGM37: 0.0 h of glucose with the dorsal and ventral tables; days with at least 4 h: 0
CGM48: 2.2 h of glucose with the dorsal and ventral tables; days with at least 4 h: 0
Vanessa: 24.8 h of glucose with the dorsal and ventral tables; days with at least 4 h: 3
Vanessa, per-day Spearman R at +10 min (days with at least 4 h of both):
Dorsalripples Ventralripples
2021-02-19 0.115 -0.088
2021-02-20 0.239 0.096
2021-02-26 NaN NaN
paired t: TtestResult(statistic=np.float64(5.838727256078682), pvalue=np.float64(0.10798626111660188), df=np.int64(1)) unpaired t: TtestResult(statistic=np.float64(1.554798920766594), pvalue=np.float64(0.26023993528172173), df=np.float64(2.0))
other rats' R at +10 min: mean -0.143 (n = 17); Vanessa dorsal per-day values [0.115 0.239]; unpaired t = 2.94, p = 0.00918
Result: blocked (data missing: enough glucose with the electrophysiology in one of three rats, on two days). Dorsal and ventral 5-min SPW-R rates are correlated with r = 0.13 (Vanessa), 0.44 (CGM37) and 0.83 (CGM48), below 0.5 in two of three rats. Only Vanessa has days with at least 4 h of glucose and both tables (three days, two of which give a correlation; CGM48 has 2.2 h in all, CGM37 none): dorsal 0.115 and 0.239, ventral -0.088 and 0.096 at +10 min, the ventral values more negative (paired p = 0.11), and Vanessa's dorsal values differ from the other rats' (p = 0.009). Vanessa's sensor reads a median of 0.25 nA, so its glucose series is doubtful.
Size of the predictor models (edfig6-ripple-rate-best-predictor, edfig6d-pairwise-predictors, edfig6ef-latent-variables)¶
Six of the seven predictors of the Methods (power spectrum slope, theta residual, theta/delta ratio, high-pass coherence between shanks, movement, and NREM or wake where no sleep states are stored) have to be computed from the LFP and the accelerometer over every hour with glucose. The LFP is stored in chunks of 707 samples by 32 channels (gzip), so even one channel means reading all 32. The cell below times the reading of 30 min of LFP and of accelerometer from each of 16 files, in 4 processes (the cores of the machine that runs the notebook) and in 16 (the reading waits on the network), and scales the wall time to all rats, to the 8 rats with the most hours (the legend's n is 8) and to the 3 rats whose brain state is stored (in the others it would first have to be scored).
a = assets[assets.path == "sub-CGM31/sub-CGM31_ses-CGM31-720um-210209-085741_ecephys.nwb"].iloc[0]
with S.open_h5(a.asset_id) as h:
lfp = h["processing/ecephys/LFP/ElectricalSeriesLFP/data"]; acc = h["acquisition/Accelerometer/data"]
print("LFP chunks", lfp.chunks, lfp.compression, "; accelerometer chunks", acc.chunks, acc.compression,
f"at {h['acquisition/Accelerometer/starting_time'].attrs['rate']} Hz")
# 30 min of LFP and accelerometer from each of 16 files, read in 4 and then in 16 processes
from concurrent.futures import ProcessPoolExecutor
pool = [f for f in kept if f["lfp_n"] / f["lfp_rate"] > 1800 and "Accelerometer" in f["acquisition"]]
pick = pool[:: max(1, len(pool) // 16)][:16]
ids = [assets.set_index("path").asset_id[f["path"]] for f in pick]
top8_both = cov.sort_values("hours_with_both", ascending=False).subject[:8].tolist()
sets = {"all 18 rats": RATS, "the 8 rats with the most hours of both": top8_both,
"the 3 rats with stored sleep states and glucose (CGM1, CGM2, CGM3)": ["CGM1", "CGM2", "CGM3"]}
for procs in (4, 16):
t0 = time.time()
with ProcessPoolExecutor(procs) as ex:
res = np.array(list(ex.map(S.time_read, ids)))
wall = time.time() - t0
sec_per_h = wall / (len(ids) * 0.5)
print(f"16 x 30 min in {procs} processes: {wall:.1f} s of wall time, {sec_per_h:.1f} s per hour of recording; "
f"per piece LFP median {np.median(res[:, 0]):.1f} s (range {res[:, 0].min():.1f} to {res[:, 0].max():.1f}), "
f"accelerometer median {np.median(res[:, 1]):.1f} s")
for name, rs in sets.items():
h_ = cov[cov.subject.isin(rs)].hours_with_both.sum()
print(f" {name}: {h_:.0f} h, reading alone about {h_ * sec_per_h / 60:.0f} min, against a notebook limit of 30 min")
LFP chunks (707, 32) gzip ; accelerometer chunks (104209, 1) gzip at 5000.0 Hz
16 x 30 min in 4 processes: 38.0 s of wall time, 4.8 s per hour of recording; per piece LFP median 4.5 s (range 1.2 to 16.8), accelerometer median 1.6 s all 18 rats: 1717 h, reading alone about 136 min, against a notebook limit of 30 min the 8 rats with the most hours of both: 1457 h, reading alone about 115 min, against a notebook limit of 30 min the 3 rats with stored sleep states and glucose (CGM1, CGM2, CGM3): 557 h, reading alone about 44 min, against a notebook limit of 30 min
16 x 30 min in 16 processes: 23.2 s of wall time, 2.9 s per hour of recording; per piece LFP median 8.7 s (range 1.7 to 21.5), accelerometer median 2.4 s all 18 rats: 1717 h, reading alone about 83 min, against a notebook limit of 30 min the 8 rats with the most hours of both: 1457 h, reading alone about 71 min, against a notebook limit of 30 min the 3 rats with stored sleep states and glucose (CGM1, CGM2, CGM3): 557 h, reading alone about 27 min, against a notebook limit of 30 min
Blocked claims¶
- edfig3e-ca1-firing-decrease, edfig8-lateral-septum-rectified. No file holds spike units (data check). The lateral septum recordings of ref. 26 are on the lab's server, which is not a pinned public source.
- edfig4-fasting-fed, edfig1-glucose-circadian-meals (meal part). Meal times are not in the dandiset.
- edfig6-ripple-rate-best-predictor, edfig6d-pairwise-predictors, edfig6ef-latent-variables. The predictors have to be computed from the LFP and the accelerometer over every hour with glucose. The reading speed varies with the network: in the runs made for this attempt a single 30-min piece of LFP took from about 1 to 22 s, and once 106 s, so the printed values change from run to run. In 4 processes the reading alone came to 4.4 to 4.6 s per hour of recording in three runs: 126 to 130 min for all 18 rats, 107 to 110 min for the 8 rats with the most hours and 41 to 42 min for the 3 rats with stored brain states. In 16 processes it came to 2.8 to 3.0 s per hour: 81 to 85 min for all rats, 69 to 72 min for the 8 rats and 26 to 28 min for the 3 rats, against the 30-min limit of a notebook. For the 3 rats the computation of the seven predictors and the 100 splits of each model at 41 offsets would come on top of the reading, and the result would cover 3 of the paper's 8 rats; a test on them is an analysis not run in this attempt. The brain state of the other 15 rats would first have to be scored.
- fig3d-induced-ripples-lower-glucose, fig3-induced-effect-size, fig3ef-stimulation-dose-response, fig3f-stimulation-vs-shuffle, edfig7ab-stimulation-protocol, edfig7jm-induced-ripples-decoupled, edfig7-ppc-stimulation-no-effect. No file holds stimulation times or a stimulation channel, and no subject is labelled as optogenetic or PPC.
- fig4-cno-controls-no-change, fig4-ls-dreadd-decouples, fig4-ls-dreadd-chance-level, fig4e-dreadd-time-course, fig4c-medial-septum-no-effect, edfig10-expression-correlates-effect. The recordings of the DREADD rats may be among the files, but nothing says which subjects had hM4Di in the lateral or medial septum or no virus, which days were CNO or vehicle days, or when the injections were given. The histology scores of Extended Data Fig. 10e are printed only in the figure.
print(f"notebook ran in {(time.time() - started) / 60:.1f} min")
notebook ran in 6.6 min
Summary¶
| Claim | Outcome |
|---|---|
| fig2-ripple-rate-predicts-glucose | reproduced under 3 of 16 combinations of readings (18 rats), resting on an assumption |
| fig2ef-correlation-beyond-shuffle | reproduced |
| fig2g-clustered-ripples | reproduced under reading (a), resting on an assumption |
| edfig4-ripple-features | partially reproduced, resting on an assumption |
| fig2h-phase-reset | not reproduced |
| edfig5-cross-frequency-coupling | partially reproduced |
| edfig1-glucose-rhythm | not reproduced |
| edfig1-glucose-circadian-meals | blocked: data missing |
| edfig1-sensor-decay | reproduced |
| supp-ripple-rate | reproduced |
| supp-ripples-stable-across-days | reproduced |
| edfig3-ripple-rate-rhythms | partially reproduced |
| edfig2-refractory-period | not reproduced |
| edfig4-nrem-wake | reproduced (3 rats), resting on an assumption |
| edfig4-no-hour-differs | not reproduced |
| edfig9-dorsal-ventral | blocked: data missing (one rat with enough glucose) |
| 19 other claims | blocked (see above) |
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: 18 rats in place of 30; the glucose dates of twelve subjects were repaired; SPW-Rs are the stored events with a peak of at least 5 s.d.; the correlation type, the delta current filter, the unit of analysis and the lag range are open in the paper and were run under every reading the claims file lists; reading (a) of the rats cannot leave out CNO days or stimulation periods, which are not labelled; the coupling analysis filters stretches of at least 12 h with both signals. The SPW-R count of a 5-min interval is paired with the glucose change over the same interval, an assumption: under the alignments one reading apart the outcomes of fig2-ripple-rate-predicts-glucose, fig2g-clustered-ripples, edfig4-ripple-features and edfig4-nrem-wake change under one of the two, while those of fig2ef-correlation-beyond-shuffle and edfig4-no-hour-differs do not.
What was found about the dataset: the glucose timestamps of twelve subjects carry dates with day
and month swapped; 27 recordings are deposited twice; the SPW-R peak column is relative to the
file while start_time is on the session clock; several SPW-R tables were made with a 4 s.d.
threshold where the Methods give 5; the session reference times are local wall-clock times labelled
UTC; there is no mapping from subjects to the paper's rats or cohorts, and no stimulation, injection,
meal or unit data; in four subjects the glucose and the electrophysiology cover different days;
the detection parameters are not stored with the SPW-R tables, and the gap from one stored event's
end to the next one's start is at least 20.0 ms in 19 rats and 15.2 ms in three.