Reproducing Peyrache et al. (2009) from CRCNS pfc-6¶
Paper. Peyrache A, Khamassi M, Benchenane K, Wiener SI, Battaglia FP (2009). Replay of rule-learning related neural patterns in the prefrontal cortex during sleep. Nature Neuroscience 12(7):919-926. doi:10.1038/nn.2337
Dataset. CRCNS pfc-6, "Activity of neurons in rat medial prefrontal cortex during learning
and sleep" (doi:10.6080/K0KH0KH5), read anonymously from the
public bucket at https://crcnsarchive.s3-us-west-2.amazonaws.com/pfc-6/live/. The data set holds
one archive, mPFC_Data.tar.gz, with seven text files per session (spikes, cell types, behaviour,
position, awake and SWS intervals), described by docs/crcns_pfc-6_data_description.pdf
(version 0.7, March 12, 2018).
Version. CRCNS data sets are not versioned (docs/archives.md), so the version is recorded as
unversioned, and the file read is pinned by its MD5 from the data set's checksums.md5:
mPFC_Data.tar.gz 2c255b44209ece01326dd09efaabd92a, checked below against the bucket's ETag.
The data set's description is checksummed too (pdf 3464d68562a83bb15d843999ec3227c8, docx
c536d32bb1a791ac00a9a7d1056f57a6); the learning sessions are read from the docx. This is the
only published version (live/; there is no old/ folder), and the paper cites none.
Claims. claims/crcns-pfc-6/10.1038-nn.2337.yaml, written from the version of record and its
supplement, and checked by a second reader. Methods taken from the paper are cited where they
are used and in pfc6.py.
| Claim | Tested here |
|---|---|
| sessions-and-cells | tested |
| rules-learned | tested |
| learning-sessions | tested |
| fig1b-signal-components | tested |
| fig2b-post-greater-than-pre | tested |
| supp5-sws-only | tested |
| supp6-pre-post-sws-alike | tested (SWS duration and firing rates; the LFP measures are not in the data set) |
| fig2b-encoding-correlation | tested |
| fig2b-pre-weaker | tested |
| supp8-other-tetrode-pairs | blocked: no tetrode per cell |
| fig2e-skewness | tested |
| fig2f-tail-contribution | tested |
| supp10-no-coactivation | tested |
| supp11-shuffle-control | tested |
| fig5a-spwr-bins, fig5b-spwr-triggered, fig5c-peak-40ms, fig5d-population-asymmetric, fig5e-autocorrelogram | blocked: no ripple times or hippocampal LFP |
| replay-in-delta-spindle-periods, fig6a-delta-waves, fig6c-spindles | blocked: no prefrontal LFP |
| fig6b-down-up | tested (reactivation part; the SPWR part needs ripple times) |
| fig6d-reactivation-spwr-correlation | blocked: no ripple times or LFP |
| fig7b-post-spwr-after-learning, fig7b-controls, fig7c-choice-point, fig7d-two-factors | blocked: no ripple times |
| supp12-behavior-unchanged | tested |
Rest epochs. The data set gives the awake interval and the SWS intervals, not the bounds of
the pre and post rest epochs. Pre is taken as the time before the awake epoch and post as the time
after it to the end of the position record, or to the end of the last post SWS interval if that is
later (an assumption; a 30 min variant is shown for supp5-sws-only,
the claim that compares non-SWS bins). Every analysis that uses SWS bins only is unaffected by it except
through the z-transform reading epochz.
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy import stats
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from statsmodels.stats.anova import AnovaRM, anova_lm
import pfc6
pd.set_option("display.width", 200)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
T0 = time.time()
Data check¶
The whole archive mPFC_Data.tar.gz is streamed once and parsed in memory. Its S3 ETag is its
MD5 for a single-part upload, so it can be compared with checksums.md5.
fs_info = pfc6.fsspec.filesystem("https").info(pfc6.URL)
print("dataset", pfc6.DATASET, "| file", pfc6.URL.rsplit("/", 1)[1], "| size", fs_info["size"],
"| ETag", fs_info.get("ETag"), "| checksums.md5", pfc6.MD5)
S = pfc6.load_sessions()
LEARN = pfc6.learning_trials_from_docs()
print(f"loaded {len(S)} sessions in {time.time() - T0:.0f} s")
info = pd.DataFrame([{
"session": k, "rat": s["rat"], "cells": pfc6.n_cells(s), "cells_with_spikes": len(np.unique(s["spikes_id"])),
"pyr": int(np.sum(s["celltype"] == 2)), "int": int(np.sum(s["celltype"] == 1)), "unk": int(np.sum(s["celltype"] == 0)),
"pre_min": (s["wake"][0] - pfc6.rest_epochs(s)[0][0]) / 60000, "wake_min": np.diff(s["wake"])[0] / 60000,
"post_min": (pfc6.rest_epochs(s)[1][1] - s["wake"][1]) / 60000,
"sws_pre_min": pfc6.sws_minutes(s["sws_pre"]), "sws_post_min": pfc6.sws_minutes(s["sws_post"]),
"trials": len(s["behavior"]), "rules": "".join(str(r) for r in np.unique(s["behavior"][:, 2]).astype(int)),
"included": pfc6.included(s)} for k, s in S.items()])
print(info.round(1).to_string(index=False))
print("\nSessions by rat:", info.groupby("rat").size().to_dict(),
"| included by rat:", info[info.included].groupby("rat").size().to_dict())
print("Cells, all sessions:", info.cells.sum(), "| with at least one spike:", info.cells_with_spikes.sum())
dataset crcns:pfc-6/unversioned | file mPFC_Data.tar.gz | size 295335687 | ETag "2c255b44209ece01326dd09efaabd92a" | checksums.md5 2c255b44209ece01326dd09efaabd92a
loaded 86 sessions in 57 s
session rat cells cells_with_spikes pyr int unk pre_min wake_min post_min sws_pre_min sws_post_min trials rules included
150628 15 31 31 28 2 1 24.9 28.3 22.4 12.1 11.0 23 1 True
150629 15 21 21 19 2 0 36.1 24.7 26.7 9.1 5.7 17 12 True
150630 15 24 24 21 3 0 29.8 37.6 29.1 11.7 13.8 18 12 True
150701 15 21 21 18 3 0 30.8 31.1 20.6 0.0 0.0 15 2 False
150704 15 25 25 21 4 0 23.5 30.4 32.3 0.0 12.9 25 2 False
150705 15 19 19 17 2 0 29.5 33.7 24.7 13.7 14.7 29 2 True
150706 15 19 19 16 2 1 27.8 26.1 28.6 12.6 15.6 20 2 True
150707 15 13 13 11 2 0 28.4 34.9 24.8 14.4 8.9 33 2 True
150708 15 8 8 6 1 1 27.0 20.2 0.1 0.0 0.0 18 2 False
150711 15 21 21 16 5 0 26.8 25.5 31.5 10.3 6.6 28 2 True
150712 15 18 18 15 3 0 27.4 33.2 32.4 10.3 5.2 41 23 True
150713 15 25 25 17 7 1 28.9 36.1 39.2 10.9 8.0 45 3 True
150714 15 19 19 13 3 3 36.2 36.6 33.4 19.0 3.8 40 3 False
150715 15 23 23 18 4 1 27.0 39.1 30.1 8.7 6.8 51 34 True
181011 18 40 40 37 2 1 27.9 26.6 25.5 12.6 12.4 9 1 True
181012 18 44 44 41 3 0 30.1 26.0 30.5 8.0 15.4 13 1 True
181014 18 37 37 34 3 0 31.5 25.1 29.8 10.7 12.2 24 2 True
181017 18 52 52 48 3 1 25.4 37.1 31.6 7.5 7.0 37 2 True
181018 18 52 52 45 5 2 39.3 44.8 30.0 21.3 14.2 42 2 True
181019 18 57 55 49 4 4 37.2 24.8 27.1 5.9 8.4 27 2 True
181020 18 55 55 50 4 1 33.9 33.9 33.2 18.6 10.7 29 2 True
181021 18 45 45 39 6 0 33.3 26.9 28.3 8.8 8.8 30 2 True
181024 18 56 56 49 5 2 0.4 31.8 33.3 0.0 0.0 35 23 False
181025 18 39 39 31 2 6 37.4 20.3 31.4 11.4 8.4 20 3 True
181026 18 46 46 42 3 1 25.2 26.6 33.7 5.6 13.0 34 34 True
181027 18 42 42 32 8 2 38.2 30.3 23.6 3.5 0.6 34 4 False
181102 18 57 57 45 11 1 29.7 30.7 36.4 9.5 8.7 38 4 True
181103 18 44 44 34 5 5 27.2 27.8 38.0 3.7 7.0 22 4 False
190213 19 21 21 17 3 1 31.2 28.4 27.2 4.3 10.1 32 1 True
190214 19 25 25 21 2 2 36.1 26.1 35.9 12.9 19.0 29 1 True
190224 19 22 21 17 3 2 29.4 22.2 26.5 6.7 12.8 15 2 True
190226 19 27 27 23 3 1 41.9 14.0 29.8 24.8 11.2 7 2 True
190227 19 24 24 19 5 0 32.0 33.5 26.3 10.7 11.6 20 2 True
190228 19 27 27 23 3 1 54.2 50.8 25.2 12.4 15.1 15 2 True
190301 19 20 20 19 1 0 38.8 52.5 34.5 11.6 15.0 12 2 True
190302 19 22 22 18 3 1 30.7 29.3 29.0 8.0 4.7 14 2 True
190303 19 19 19 15 2 2 34.6 48.5 26.4 17.4 11.0 29 2 True
190307 19 18 18 15 3 0 37.3 73.9 41.6 9.3 7.7 25 2 True
190308 19 17 16 11 4 2 35.2 56.7 27.3 15.8 5.6 22 2 True
200102 20 20 20 16 4 0 33.2 25.3 26.9 14.7 5.6 41 4 True
200103 20 22 22 19 3 0 24.0 24.2 30.6 3.0 3.7 44 4 False
200104 20 18 18 18 0 0 32.3 23.2 32.1 5.3 9.3 54 4 True
200105 20 22 22 16 6 0 25.7 23.6 22.8 3.3 2.7 42 4 False
200106 20 20 20 17 2 1 36.3 29.2 58.7 9.9 8.6 59 4 True
200107 20 23 23 0 0 23 49.2 27.1 58.6 9.4 13.0 40 4 True
200109 20 22 22 20 1 1 39.8 48.7 58.6 2.8 13.9 68 14 False
200110 20 22 22 19 3 0 58.5 33.8 43.4 12.3 6.0 52 34 True
200111 20 28 28 21 6 1 35.1 32.6 43.7 8.6 22.0 50 23 True
200112 20 23 23 20 3 0 43.7 24.5 42.9 11.0 6.2 49 2 True
200116 20 24 24 22 2 0 34.6 37.0 28.5 0.7 13.4 45 2 False
200118 20 16 16 10 4 2 24.2 25.4 25.1 5.2 6.3 34 12 True
200119 20 19 19 13 4 2 24.5 54.0 22.3 3.8 6.2 26 1 False
200124 20 19 19 16 3 0 22.7 48.6 25.8 12.8 11.5 28 14 True
200125 20 14 14 11 3 0 25.5 27.2 32.8 7.3 11.6 35 4 True
200126 20 19 19 15 3 1 34.7 23.1 37.9 3.6 8.0 46 4 False
200127 20 13 13 9 3 1 27.5 25.7 30.6 3.4 3.5 35 4 False
200130 20 11 11 7 4 0 36.8 27.6 41.3 2.4 17.0 34 4 False
200131 20 21 21 9 0 12 37.0 45.4 39.1 14.1 18.7 36 14 True
200203 20 21 21 16 3 2 40.7 44.4 28.7 4.1 5.2 38 12 True
200206 20 21 21 16 5 0 29.7 34.1 30.1 4.7 6.1 38 2 True
200208 20 11 11 9 2 0 38.1 34.3 57.5 0.8 11.8 35 2 False
200209 20 13 13 0 13 0 33.9 40.3 25.9 9.0 9.4 34 23 True
200210 20 14 14 11 3 0 40.7 36.3 32.8 8.8 10.3 46 2 True
200213 20 16 16 12 2 2 32.8 26.3 -0.9 0.0 0.0 39 23 False
200214 20 26 26 22 3 1 48.4 30.4 36.4 0.6 6.5 50 34 False
200216 20 16 16 15 0 1 32.9 21.9 28.9 8.0 7.6 18 4 True
200222 20 15 15 12 2 1 36.5 21.9 23.4 9.5 6.4 24 4 True
200223 20 19 19 18 1 0 25.4 31.7 22.4 3.7 6.3 30 4 False
200224 20 15 15 4 0 11 33.8 36.0 15.1 1.5 5.6 40 4 False
200227 20 15 15 11 2 2 33.1 43.1 27.2 1.1 8.2 47 4 False
200228 20 21 21 18 3 0 37.8 42.7 24.6 4.7 14.3 34 4 True
200301 20 50 50 13 2 35 27.5 33.8 32.1 3.5 10.2 32 14 False
200302 20 13 13 8 3 2 25.5 30.7 26.3 1.5 1.3 49 12 False
200303 20 19 19 11 5 3 41.0 25.0 34.6 1.5 2.2 37 2 False
200308 20 17 17 12 2 3 28.3 23.2 25.3 1.4 5.1 23 2 False
200309 20 13 13 9 4 0 33.8 55.1 25.5 1.8 4.5 37 2 False
201219 20 18 18 15 1 2 25.1 20.7 21.0 7.0 0.7 17 1 False
201220 20 26 26 22 2 2 25.4 36.2 29.2 11.5 12.2 37 12 True
201221 20 37 36 29 5 3 42.5 38.5 33.2 20.4 16.0 30 2 True
201222 20 27 27 22 5 0 32.6 46.2 39.8 9.4 6.4 48 2 True
201223 20 26 26 23 3 0 29.7 27.2 43.3 9.5 12.3 27 2 True
201226 20 28 28 24 3 1 31.0 28.3 24.2 4.9 10.1 39 2 True
201227 20 23 23 20 2 1 34.1 29.8 33.4 2.2 5.5 42 2 False
201228 20 30 30 28 2 0 29.1 29.5 33.7 5.1 7.2 37 23 True
201229 20 21 21 20 1 0 30.1 21.1 25.7 9.0 3.3 26 3 False
201230 20 13 12 10 2 1 35.5 24.7 42.8 3.6 6.8 45 34 False
Sessions by rat: {15: 14, 18: 14, 19: 11, 20: 47} | included by rat: {15: 10, 18: 11, 19: 11, 20: 24}
Cells, all sessions: 2155 | with at least one spike: 2149
checks = []
for k, s in S.items():
t, c = s["spikes_t"], s["spikes_id"]
w0, w1 = s["wake"]
(p0, p1), (q0, q1) = pfc6.rest_epochs(s)
b = s["behavior"]
rule, right, went_lit = b[:, 2], b[:, 4] == 0, b[:, 4] == b[:, 5]
expected = np.where(rule == 1, right, np.where(rule == 3, ~right, np.where(rule == 2, went_lit, ~went_lit)))
checks.append({
"session": k,
"ids_in_1..n": bool(c.min() >= 1 and c.max() <= pfc6.n_cells(s)),
"sorted_by_cell_then_time": bool(np.all(np.diff(c) >= 0) and all(np.all(np.diff(t[c == i]) >= 0) for i in np.unique(c))),
"one_wake_row": s["n_wake_rows"] == 1,
"sws_pre_in_pre": bool(len(s["sws_pre"]) == 0 or (s["sws_pre"].min() >= p0 and s["sws_pre"].max() <= p1)),
"sws_post_in_post": bool(len(s["sws_post"]) == 0 or (s["sws_post"].min() >= q0 and s["sws_post"].max() <= q1)),
"trials_in_wake": bool(b[:, 0].min() >= w0 and b[:, 1].max() <= w1),
"correct_matches_rule": bool(np.all(expected == (b[:, 3] == 1))),
"spikes_after_pos_end_min": (t.max() - s["pos"][-1, 0]) / 60000,
"pos_dt_ms": float(np.median(np.diff(s["pos"][:, 0]))),
})
checks = pd.DataFrame(checks)
for col in checks.columns[1:8]:
print(f"{col}: {int(checks[col].sum())} of {len(checks)} sessions", "" if checks[col].all() else f"(fails: {checks.session[~checks[col]].tolist()})")
for k in checks.session[~checks.sws_post_in_post]:
print(f" {k}: wake ends {S[k]['wake'][1] / 1000:.0f} s, post SWS intervals (s) {np.round(S[k]['sws_post'] / 1000).tolist()}, included {pfc6.included(S[k])}")
print("position sampling interval (ms), median per session: %.1f to %.1f" % (checks.pos_dt_ms.min(), checks.pos_dt_ms.max()))
late = checks[checks.spikes_after_pos_end_min > 1]
print("sessions with spikes more than 1 min after the end of the position record:\n", late[["session", "spikes_after_pos_end_min"]].round(1).to_string(index=False))
allpos = np.vstack([s["pos"] for s in S.values()])
print("position x range %.1f to %.1f, y range %.1f to %.1f (units not stated; pixels assumed)" % (
allpos[:, 1].min(), allpos[:, 1].max(), allpos[:, 2].min(), allpos[:, 2].max()))
print("spike times range %.0f to %.0f ms" % (min(s["spikes_t"].min() for s in S.values()), max(s["spikes_t"].max() for s in S.values())))
print("behavior columns: rule values", sorted(np.unique(np.vstack([s["behavior"] for s in S.values()])[:, 2]).astype(int)),
"| correct values", sorted(np.unique(np.vstack([s["behavior"] for s in S.values()])[:, 3]).astype(int)),
"| arm values", sorted(np.unique(np.vstack([s["behavior"] for s in S.values()])[:, 4]).astype(int)),
"| light values", sorted(np.unique(np.vstack([s["behavior"] for s in S.values()])[:, 5]).astype(int)))
trial_len = np.concatenate([np.diff(s["behavior"][:, :2], axis=1).ravel() for s in S.values()]) / 1000
print("trial length (s): median %.1f, range %.1f to %.1f" % (np.median(trial_len), trial_len.min(), trial_len.max()))
# recordings deposited twice: identical spike trains of any cell across sessions, or identical position tracks
sig = {}
for k, s in S.items():
for i in np.unique(s["spikes_id"]):
tt = s["spikes_t"][s["spikes_id"] == i]
sig.setdefault((len(tt), round(float(tt[:10].sum()), 1)), set()).add(k)
pos_sig = {}
for k, s in S.items():
pos_sig.setdefault((len(s["pos"]), round(float(s["pos"][:100, 1:].sum()), 1)), set()).add(k)
print("cells with an identical spike train in two sessions:", sum(len(v) > 1 for v in sig.values()),
"| sessions with an identical position track:", sum(len(v) > 1 for v in pos_sig.values()))
print("learning sessions and trials from the data set's description:", LEARN)
ids_in_1..n: 86 of 86 sessions sorted_by_cell_then_time: 86 of 86 sessions one_wake_row: 86 of 86 sessions sws_pre_in_pre: 86 of 86 sessions sws_post_in_post: 82 of 86 sessions (fails: ['150701', '150708', '181024', '200213']) trials_in_wake: 86 of 86 sessions correct_matches_rule: 86 of 86 sessions 150701: wake ends 3713 s, post SWS intervals (s) [[0.0, 0.0]], included False 150708: wake ends 2829 s, post SWS intervals (s) [[0.0, 0.0]], included False 181024: wake ends 1928 s, post SWS intervals (s) [[0.0, 0.0]], included False 200213: wake ends 3545 s, post SWS intervals (s) [[0.0, 0.0]], included False position sampling interval (ms), median per session: 33.3 to 33.4 sessions with spikes more than 1 min after the end of the position record: session spikes_after_pos_end_min 181024 40.9 181027 11.8 200106 176.1 200213 29.6 200224 27.6 position x range 28.4 to 282.8, y range 22.6 to 238.2 (units not stated; pixels assumed) spike times range 0 to 18014354 ms behavior columns: rule values [np.int64(1), np.int64(2), np.int64(3), np.int64(4)] | correct values [np.int64(0), np.int64(1)] | arm values [np.int64(0), np.int64(1)] | light values [np.int64(0), np.int64(1)] trial length (s): median 4.2, range 1.5 to 435.1
cells with an identical spike train in two sessions: 0 | sessions with an identical position track: 0
learning sessions and trials from the data set's description: {'201222': 23, '201227': 23, '201229': 14, '181012': 5, '181020': 10, '150628': 12, '150630': 5, '150707': 18, '190214': 11, '190228': 8}
Do the streams share the clock of the trial table?¶
Two contrasts. The arm the rat chose (behavior column 5) should be decodable from the position and from the spikes at the end of each trial, and should not be when the trial times are moved by 30 s. And head speed from the position record should be lower inside the SWS intervals than in the rest of the rest epochs.
def arm_decoding(s, shift_ms, use="spikes"):
b = s["behavior"]
y = b[:, 4].astype(int)
if min(np.bincount(y, minlength=2)) < 3:
return np.nan
ends = b[:, 1] + shift_ms
if use == "pos":
p = s["pos"]
X = np.column_stack([np.interp(ends - 500, p[:, 0], p[:, 1]), np.interp(ends - 500, p[:, 0], p[:, 2])])
else:
t, c = s["spikes_t"], s["spikes_id"] - 1
X = np.zeros((len(ends), pfc6.n_cells(s)))
for k, e in enumerate(ends):
m = (t >= e - 2000) & (t < e)
np.add.at(X[k], c[m], 1)
X = np.sqrt(X)
folds = min(5, min(np.bincount(y, minlength=2)))
return cross_val_score(LogisticRegression(max_iter=2000), X, y, cv=folds).mean()
clock = []
for k, s in S.items():
r = {"session": k}
for use in ("pos", "spikes"):
for shift in (0, 30000):
r[f"{use}_{shift // 1000}s"] = arm_decoding(s, shift, use)
y = s["behavior"][:, 4]
r["chance"] = max(np.mean(y), 1 - np.mean(y))
p = s["pos"]
sp = np.hypot(np.diff(p[:, 1]), np.diff(p[:, 2])) / np.diff(p[:, 0]) * 1000
tm = p[1:, 0]
rest = np.zeros(len(tm), bool)
sws = np.zeros(len(tm), bool)
for (a, e), iv in zip(pfc6.rest_epochs(s), (s["sws_pre"], s["sws_post"])):
rest |= (tm >= a) & (tm < e)
for x0, x1 in iv:
sws |= (tm >= x0) & (tm < x1)
r["speed_sws"] = np.nanmean(sp[sws] > 5) if sws.any() else np.nan
r["speed_rest_nonsws"] = np.nanmean(sp[rest & ~sws] > 5) if (rest & ~sws).any() else np.nan
clock.append(r)
clock = pd.DataFrame(clock)
print("Arm decoding accuracy (5-fold logistic regression), median over sessions:")
print(clock[["pos_0s", "pos_30s", "spikes_0s", "spikes_30s", "chance"]].median().round(3).to_string())
print("sessions where spike decoding with true times beats the 30 s shifted times:",
int(np.sum(clock.spikes_0s > clock.spikes_30s)), "of", int(clock.spikes_0s.notna().sum()))
print("Wilcoxon, spikes true against shifted: p = %.2g" % stats.wilcoxon(clock.spikes_0s, clock.spikes_30s, nan_policy="omit").pvalue)
print("sessions where spike decoding with true times is at or below chance:",
clock.loc[clock.spikes_0s <= clock.chance, "session"].tolist())
both = clock.dropna(subset=["speed_sws", "speed_rest_nonsws"])
print("fraction of position frames with head speed above 5 units/s, median over sessions: SWS %.4f, rest outside SWS %.4f; lower in SWS in %d of %d sessions"
% (both.speed_sws.median(), both.speed_rest_nonsws.median(), int(np.sum(both.speed_sws < both.speed_rest_nonsws)), len(both)))
Arm decoding accuracy (5-fold logistic regression), median over sessions: pos_0s 1.000 pos_30s 0.786 spikes_0s 0.812 spikes_30s 0.676 chance 0.730 sessions where spike decoding with true times beats the 30 s shifted times: 55 of 71 Wilcoxon, spikes true against shifted: p = 1.8e-09 sessions where spike decoding with true times is at or below chance: ['150630', '150706', '150714', '181024', '190213', '190214', '190307', '190308', '200102', '200105', '200124', '200127', '200203', '200208', '200209', '200302', '200308', '200309'] fraction of position frames with head speed above 5 units/s, median over sessions: SWS 0.0056, rest outside SWS 0.2277; lower in SWS in 82 of 82 sessions
Data check. The ETag of mPFC_Data.tar.gz equals its MD5 in checksums.md5. The archive holds 86 sessions
(rat 15, 14; rat 18, 14; rat 19, 11; rat 20, 47), 2155 cells and 2834 trials; the description lists no session count.
Cell ids run from 1 to the number of lines of CellType.dat, and 6 cells have no spike. Spikes are sorted by cell, then
time; every session has one awake interval, and its trials lie inside it. The correct column agrees with the rule, the
arm chosen and the lit arm in every trial, so the rule codes, arm codes and light codes of the description are
consistent with the data. SWS intervals lie in the rest epochs; four sessions without post SWS store a single [0, 0]
interval, and one an empty file. No cell's spike train or session's position track appears twice. Position units are not
stated (x 28 to 283, y 23 to 238, 30 Hz). Two sessions have cell types all 0 (200107) or all 1 (200209), which is
unlikely for prefrontal recordings; cell types are not used here. In five sessions spikes continue 12 to 176 min after
the position record ends. Clocks: the chosen arm is decoded from position at trial end (median accuracy 1.00, 0.79 with
times shifted by 30 s) and from spikes (0.81 against 0.68 shifted; 55 of 71 sessions better with true times, Wilcoxon
p = 1.8e-09), so position and spikes are consistent with the trial table's clock; in the 18 sessions listed where spike
decoding is at or below chance this check does not establish it. Head movement is rarer inside SWS intervals than in
the rest of the rest epochs in all 82 sessions with both, consistent with the SWS intervals sharing the position clock.
sessions-and-cells¶
The paper recorded 1,692 mPFC cells from four rats in 63 sessions (rat 15, 16 sessions; rat 18, 11; rat 19, 12; rat 20, 24), and analyzed only sessions with at least ten cells and at least 4 min of SWS in each rest epoch.
Reading (a) counts all sessions; reading (b) the sessions that pass the inclusion rule.
paper_by_rat = {15: 16, 18: 11, 19: 12, 20: 24}
tab = pd.DataFrame({
"paper_sessions": pd.Series(paper_by_rat),
"a_sessions": info.groupby("rat").size(), "a_cells": info.groupby("rat").cells.sum(),
"b_sessions": info[info.included].groupby("rat").size(), "b_cells": info[info.included].groupby("rat").cells.sum(),
"b_cells_with_spikes": info[info.included].groupby("rat").cells_with_spikes.sum()})
tab.loc["total"] = tab.sum()
print(tab.to_string())
for name, n_s, n_c in (("(a) all sessions", tab.loc["total", "a_sessions"], tab.loc["total", "a_cells"]),
("(b) inclusion rule", tab.loc["total", "b_sessions"], tab.loc["total", "b_cells"])):
print(f"{name}: {n_s} sessions (expected 63, 57 to 69: {abs(n_s - 63) <= 6}), "
f"{n_c} cells (expected 1692, 1523 to 1861: {abs(n_c - 1692) <= 169})")
fails = info[~info.included].assign(
reason=lambda d: np.where(d.cells < 10, "fewer than 10 cells", "") + np.where(
(d.sws_pre_min < 4) | (d.sws_post_min < 4), " SWS under 4 min", ""))
print("\nExcluded sessions by rat and reason:")
print(fails.groupby(["rat", "reason"]).size().to_string())
for thr in (3.0, 3.5, 4.0, 5.0):
m = (info.cells >= 10) & (info.sws_pre_min >= thr) & (info.sws_post_min >= thr)
print(f"variant: SWS threshold {thr} min -> {int(m.sum())} sessions, {int(info.cells[m].sum())} cells, by rat {info[m].groupby('rat').size().to_dict()}")
INC = info.loc[info.included, "session"].tolist()
paper_sessions a_sessions a_cells b_sessions b_cells b_cells_with_spikes
15 16 14 287 10 214 214
18 11 14 666 11 524 522
19 12 11 242 11 242 240
20 24 47 960 24 519 518
total 63 86 2155 56 1499 1494
(a) all sessions: 86 sessions (expected 63, 57 to 69: False), 2155 cells (expected 1692, 1523 to 1861: False)
(b) inclusion rule: 56 sessions (expected 63, 57 to 69: False), 1499 cells (expected 1692, 1523 to 1861: False)
Excluded sessions by rat and reason:
rat reason
15 SWS under 4 min 3
fewer than 10 cells SWS under 4 min 1
18 SWS under 4 min 3
20 SWS under 4 min 23
variant: SWS threshold 3.0 min -> 66 sessions, 1738 cells, by rat {15: 11, 18: 12, 19: 11, 20: 32}
variant: SWS threshold 3.5 min -> 63 sessions, 1682 cells, by rat {15: 11, 18: 12, 19: 11, 20: 29}
variant: SWS threshold 4.0 min -> 56 sessions, 1499 cells, by rat {15: 10, 18: 11, 19: 11, 20: 24}
variant: SWS threshold 5.0 min -> 50 sessions, 1365 cells, by rat {15: 10, 18: 11, 19: 9, 20: 20}
Result: not reproduced. Neither reading matches. All sessions (a) give 86 sessions and 2155 cells; the inclusion rule (b) gives 56 sessions and 1499 cells, against 63 (57 to 69) and 1,692 (1523 to 1861). Under (b) rats 18 and 20 match the paper's per-rat counts exactly (11 and 24), rat 19 has 11 against 12, and rat 15 has 10 against 16; the data set holds only 14 sessions of rat 15 in all. A threshold of 3.5 min gives 63 sessions and 1682 cells in total, but per-rat counts (11, 12, 11, 29) that do not match. Two of the description's ten learning sessions (the paper's Figure 7 counts ten) fail the 4 min rule with the data set's SWS intervals (201227 with 2.2 min of pre SWS, 201229 with 3.3 min of post SWS), so the deposited SWS intervals, or the rule as applied, differ from the paper's.
rules-learned¶
All four rats learned the right-arm and lit-arm rules; two learned the left-arm rule and one the dark-arm rule.
Rule codes are from the data description (1 right, 2 lit arm, 3 left, 4 dark arm). Sessions are ordered by date within each rat; the session id holds month and day only, and rat 20's December sessions are taken to precede its January to March ones.
RULE = {1: "right", 2: "lit", 3: "left", 4: "dark"}
def session_order(rat):
ks = sorted(k for k in S if S[k]["rat"] == rat)
months = {int(k[2:4]) for k in ks}
key = lambda k: (int(k[2:4]) + (12 if (int(k[2:4]) < 7 and max(months) >= 10) else 0), int(k[4:6]))
return sorted(ks, key=key)
def rule_runs(trials):
"""Split (rule, correct) rows into runs of consecutive trials under one rule."""
out, cur = [], []
for r, c in trials:
if cur and r != cur[-1][0]:
out.append(cur); cur = []
cur.append((r, c))
if cur:
out.append(cur)
return [(run[0][0], [c for _, c in run]) for run in out]
learned = {"a_within_session": {}, "b_across_sessions": {}}
encountered = {}
for rat in sorted({s["rat"] for s in S.values()}):
ks = session_order(rat)
allrows = []
for k in ks:
rows = list(zip(S[k]["behavior"][:, 2].astype(int), S[k]["behavior"][:, 3].astype(int)))
allrows += rows
for rule, corr in rule_runs(rows):
if pfc6.reached_criterion(corr):
learned["a_within_session"].setdefault(rule, set()).add(rat)
for rule, corr in rule_runs(allrows):
encountered.setdefault(rule, set()).add(rat)
if pfc6.reached_criterion(corr):
learned["b_across_sessions"].setdefault(rule, set()).add(rat)
print(f"rat {rat}: rule sequence over sessions", [RULE[r] for r, _ in rule_runs(allrows)])
expected = {1: 4, 2: 4, 3: 2, 4: 1}
for reading, d in learned.items():
counts = {RULE[r]: len(d.get(r, set())) for r in (1, 2, 3, 4)}
print(reading, "rats per rule:", counts, "rats:", {RULE[r]: sorted(d.get(r, set())) for r in (1, 2, 3, 4)},
"| all four match:", all(len(d.get(r, set())) == v for r, v in expected.items()))
changed = {}
for rat in sorted({s["rat"] for s in S.values()}):
allrows = []
for k in session_order(rat):
allrows += list(zip(S[k]["behavior"][:, 2].astype(int), S[k]["behavior"][:, 3].astype(int)))
for rule, _ in rule_runs(allrows)[:-1]:
changed.setdefault(rule, set()).add(rat)
print("variant (not a reading of the criterion): rats whose rule was changed after it, as the task changes the rule at criterion:",
{RULE[r]: len(changed.get(r, set())) for r in (1, 2, 3, 4)})
best = []
for rat in sorted({s["rat"] for s in S.values()}):
allrows = []
for k in session_order(rat):
allrows += list(zip(S[k]["behavior"][:, 2].astype(int), S[k]["behavior"][:, 3].astype(int)))
for rule, corr in rule_runs(allrows):
c = np.asarray(corr, int)
longest = max(len(x) for x in "".join(map(str, c)).split("0"))
win12 = max((c[i:i + 12].sum() for i in range(len(c) - 11)), default=None)
best.append({"rat": rat, "rule": RULE[rule], "trials": len(c), "longest_correct_run": longest, "best_12_window": win12})
best = pd.DataFrame(best).groupby(["rat", "rule"], sort=False).agg(
blocks=("trials", "size"), trials=("trials", "sum"), longest_correct_run=("longest_correct_run", "max"), best_12_window=("best_12_window", "max"))
print("\nPer rat and rule, over blocks of consecutive trials under the rule across sessions (criterion: a run of 10, or 11 correct in 12):")
print(best.to_string())
print("rats that encountered each rule:", {RULE[r]: sorted(v) for r, v in sorted(encountered.items())})
rat 15: rule sequence over sessions ['right', 'lit', 'right', 'lit', 'left', 'dark']
rat 18: rule sequence over sessions ['right', 'lit', 'left', 'dark']
rat 19: rule sequence over sessions ['right', 'lit']
rat 20: rule sequence over sessions ['right', 'lit', 'left', 'dark', 'right', 'dark', 'left', 'lit', 'right', 'dark', 'right', 'lit', 'left', 'lit', 'left', 'dark', 'right', 'lit']
a_within_session rats per rule: {'right': 3, 'lit': 4, 'left': 3, 'dark': 1} rats: {'right': [15, 19, 20], 'lit': [15, 18, 19, 20], 'left': [15, 18, 20], 'dark': [20]} | all four match: False
b_across_sessions rats per rule: {'right': 3, 'lit': 4, 'left': 3, 'dark': 1} rats: {'right': [15, 19, 20], 'lit': [15, 18, 19, 20], 'left': [15, 18, 20], 'dark': [20]} | all four match: False
variant (not a reading of the criterion): rats whose rule was changed after it, as the task changes the rule at criterion: {'right': 4, 'lit': 3, 'left': 3, 'dark': 1}
Per rat and rule, over blocks of consecutive trials under the rule across sessions (criterion: a run of 10, or 11 correct in 12):
blocks trials longest_correct_run best_12_window
rat rule
15 right 2 38 6 11.0
lit 2 215 10 11.0
left 1 117 14 12.0
dark 1 33 2 6.0
18 right 1 22 9 10.0
lit 1 199 14 12.0
left 1 55 13 12.0
dark 1 118 4 8.0
19 right 1 61 27 12.0
lit 1 159 12 12.0
20 right 5 185 11 11.0
lit 5 668 13 12.0
left 4 150 10 11.0
dark 4 814 13 12.0
rats that encountered each rule: {'right': [15, 18, 19, 20], 'lit': [15, 18, 19, 20], 'left': [15, 18, 20], 'dark': [15, 18, 20]}
Result: not reproduced. By the criterion (ten correct in a row, or at most one error in 12), counted within sessions (a) or across consecutive sessions (b), the counts are the same: right arm 3 rats (rat 18 never reaches criterion on it in the recorded trials), lit arm 4, left arm 3 (rats 15, 18 and 20), dark arm 1. The paper's 4, 4, 2, 1 match for two of the four rules. The rule was changed after the right-arm rule for all four rats, so the recorded trials may not hold every trial on which the criterion was reached. The table of longest runs shows how near each rat came: on the right-arm rule rat 18 has 22 trials, a longest run of 9 and at most 10 correct in 12, one trial short of either form of the criterion. On the left-arm rule rats 15 and 18 both pass clearly (runs of 14 and 13), so the extra rat there is not a matter of the criterion's edge.
learning-sessions¶
Ten sessions, from four rats, contained rule learning, with a learning point defined as the first trial of a block of at least three consecutive rewarded trials after which performance to the end of the session stayed above 80%.
Readings of the criterion: performance counted (a) from the first trial of the block, (b) from the trial after it; sessions (i) with a learning point after the first trial, (ii) the same, keeping only sessions where the rule in force at the point had not reached the rule-change criterion in an earlier session. The data set's description lists ten learning sessions with their learning trials; they are compared below.
# a block that runs to the end of the session counts as passing under reading (b)
assert pfc6.learning_point([0, 0, 1, 1, 1, 1, 1], include_block=True) == 3
assert pfc6.learning_point([0, 0, 1, 1, 1, 1, 1], include_block=False) == 3
assert pfc6.learning_point([1, 1, 1, 0, 0], include_block=False) is None
reached_before = {}
lp_rows = []
for rat in sorted({s["rat"] for s in S.values()}):
done = set()
for k in session_order(rat):
b = S[k]["behavior"]
for perf in ("a", "b"):
lp = pfc6.learning_point(b[:, 3], include_block=(perf == "a"))
rule_at = int(b[lp - 1, 2]) if lp else None
lp_rows.append({"session": k, "rat": rat, "perf": perf, "lp": lp, "rule": rule_at,
"i": lp is not None and lp > 1,
"ii": lp is not None and lp > 1 and rule_at not in done,
"included": k in INC, "doc_lp": LEARN.get(k)})
for rule, corr in rule_runs(list(zip(b[:, 2].astype(int), b[:, 3].astype(int)))):
if pfc6.reached_criterion(corr):
done.add(rule)
lp_df = pd.DataFrame(lp_rows)
for perf in ("a", "b"):
for sel in ("i", "ii"):
d = lp_df[(lp_df.perf == perf) & lp_df[sel]]
print(f"reading ({perf})({sel}): {len(d)} sessions, rats {sorted(d.rat.unique())}, "
f"9 to 11 over four rats: {9 <= len(d) <= 11 and d.rat.nunique() == 4}; "
f"of the description's ten: {int(d.session.isin(list(LEARN)).sum())}")
print("\nThe description's ten sessions, with the learning point computed here:")
doc = lp_df[lp_df.session.isin(list(LEARN))].pivot(index="session", columns="perf", values="lp")
doc["description"] = [LEARN[k] for k in doc.index]
doc["rat"] = [S[k]["rat"] for k in doc.index]
doc["included"] = [k in INC for k in doc.index]
print(doc.to_string())
print("\nSessions under reading (a)(ii) that are not in the description's list:")
d = lp_df[(lp_df.perf == "a") & lp_df["ii"] & ~lp_df.session.isin(list(LEARN))]
print(d[["session", "rat", "lp", "rule", "included"]].to_string(index=False))
reading (a)(i): 34 sessions, rats [np.int64(15), np.int64(18), np.int64(19), np.int64(20)], 9 to 11 over four rats: False; of the description's ten: 10 reading (a)(ii): 16 sessions, rats [np.int64(15), np.int64(18), np.int64(19), np.int64(20)], 9 to 11 over four rats: False; of the description's ten: 8 reading (b)(i): 22 sessions, rats [np.int64(15), np.int64(18), np.int64(19), np.int64(20)], 9 to 11 over four rats: False; of the description's ten: 7 reading (b)(ii): 14 sessions, rats [np.int64(15), np.int64(18), np.int64(19), np.int64(20)], 9 to 11 over four rats: False; of the description's ten: 7 The description's ten sessions, with the learning point computed here: perf a b description rat included session 150628 12.0 12.0 12 15 True 150630 5.0 15.0 5 15 True 150707 10.0 18.0 18 15 True 181012 5.0 5.0 5 18 True 181020 4.0 4.0 10 18 True 190214 11.0 NaN 11 19 True 190228 8.0 8.0 8 19 True 201222 37.0 NaN 23 20 True 201227 31.0 40.0 23 20 False 201229 14.0 NaN 14 20 False Sessions under reading (a)(ii) that are not in the description's list: session rat lp rule included 150713 15 31.0 3.0 True 181011 18 6.0 1.0 True 181024 18 23.0 3.0 False 190213 19 7.0 1.0 True 190227 19 16.0 2.0 True 190301 19 2.0 2.0 True 201219 20 8.0 1.0 False 200109 20 62.0 4.0 False
Result: not reproduced. No reading gives 9 to 11 learning sessions over four rats: (a)(i) 34
sessions, (a)(ii) 16, (b)(i) 22 and (b)(ii) 14, each over four rats. Under reading (b) a block that runs to the end of
the session leaves no trials after it and counts as passing (checked above on a synthetic input); an earlier version of
learning_point dropped such sessions. The data set's description lists ten learning sessions from four rats, the
paper's count. Under reading (a) the definition gives the description's learning trial in 6 of those 10 sessions
(150628, 150630, 181012, 190214, 190228, 201229); under reading (b) in 4 (150628, 150707, 181012, 190228). Neither
reading recovers 181020 (trial 4 against 10), 201222 (37 or none against 23) or 201227 (31 or 40 against 23). Two of the ten
(201227, 201229) do not pass the 4 min SWS inclusion rule with the data set's SWS intervals.
Reactivation: awake templates and rest-epoch reactivation strengths¶
For every included session: awake-epoch template (100 ms bins over the awake interval,
z-transformed per cell, C = QQ'/B, eigenvectors), and R_l(t) in the pre and post rest epochs
under the two z-transform readings of the claims file: epochz, z over the whole rest
epoch, and swsz, z over the SWS bins alone. Cells silent in an epoch get Q = 0 there.
RX = {k: pfc6.session_reactivation(S[k]) for k in INC}
comp = []
for k, r in RX.items():
sig = r["lam"] > r["lmax"]
non = (r["lam"] > r["lmin"]) & ~sig
for l in range(r["n"]):
row = {"session": k, "rat": S[k]["rat"], "l": l, "lam": r["lam"][l], "phi": r["phi"][l],
"cls": "signal" if sig[l] else ("nonsignal" if non[l] else "below_lmin")}
for ep in ("pre", "post"):
e = r[ep]
row[f"{ep}_epochz"] = e["R_sws_epochz"][l].mean()
row[f"{ep}_swsz"] = e["R_sws_swsz"][l].mean()
row[f"{ep}_nonsws"] = e["R_nonsws_epochz"][l].mean() if (~e["sws"]).any() else np.nan
row[f"{ep}_skew_epochz"] = stats.skew(e["R_sws_epochz"][l])
row[f"{ep}_skew_swsz"] = stats.skew(e["R_sws_swsz"][l])
comp.append(row)
comp = pd.DataFrame(comp)
comp["group"] = pd.cut(comp.phi, [-np.inf, 1, 1.4, 1.8, np.inf], labels=["<1", "1-1.4", "1.4-1.8", ">1.8"])
SIG = comp[comp.cls == "signal"]
GROUPS = [">1.8", "1.4-1.8", "1-1.4"]
print("components by class:", comp.cls.value_counts().to_dict())
print("silent cells in the awake epoch, all sessions:", sum(r["silent_awake"] for r in RX.values()))
print("SWS bins per session, pre: %d to %d, post: %d to %d" % (
min(r["pre"]["sws"].sum() for r in RX.values()), max(r["pre"]["sws"].sum() for r in RX.values()),
min(r["post"]["sws"].sum() for r in RX.values()), max(r["post"]["sws"].sum() for r in RX.values())))
print(f"elapsed {time.time() - T0:.0f} s")
components by class: {'nonsignal': 730, 'below_lmin': 484, 'signal': 285}
silent cells in the awake epoch, all sessions: 12
SWS bins per session, pre: 2448 to 14895, post: 2799 to 13194
elapsed 95 s
fig1b-signal-components¶
Comparing the awake-epoch correlation matrix eigenvalues with the Marcenko-Pastur bound typically gives 1 to 6 signal components per session; 323 components were signal and 811 nonsignal in total.
per = comp.groupby("session").cls.value_counts().unstack(fill_value=0)
per["rat"] = [S[k]["rat"] for k in per.index]
n_sig, n_non = int(per.signal.sum()), int(per.nonsignal.sum())
print(f"included sessions ({len(per)}): signal {n_sig} (expected 323, 291 to 355: {abs(n_sig - 323) <= 32}), "
f"nonsignal {n_non} (expected 811, 730 to 892: {abs(n_non - 811) <= 81}), below lambda_min {int(per.below_lmin.sum())}")
print("signal components per session:", per.signal.value_counts().sort_index().to_dict())
print("sessions with 1 to 6 signal components: %d of %d" % (((per.signal >= 1) & (per.signal <= 6)).sum(), len(per)))
print("by rat:\n", per.groupby("rat")[["signal", "nonsignal", "below_lmin"]].sum().to_string())
# variant: every session with an awake epoch and at least two cells (reading (a) of sessions-and-cells)
va = {"signal": 0, "nonsignal": 0}
for k, s in S.items():
cw, _ = pfc6.bin_counts(s, *s["wake"])
lam, _, lmin, lmax = pfc6.template(pfc6.zscore(cw))
va["signal"] += int(np.sum(lam > lmax)); va["nonsignal"] += int(np.sum((lam > lmin) & (lam <= lmax)))
print("variant, all 86 sessions:", va)
def awake_classes(s, drop_silent):
cw, _ = pfc6.bin_counts(s, *s["wake"])
if drop_silent:
cw = cw[cw.sum(1) > 0]
lam, _, lmin, lmax = pfc6.template(pfc6.zscore(cw))
phi = lam[lam > lmax] / lmax
return int(np.sum(lam > lmax)), int(np.sum((lam > lmin) & (lam <= lmax))), [int(np.sum(phi > 1.8)), int(np.sum((phi > 1.4) & (phi <= 1.8))), int(np.sum(phi <= 1.4))]
for thr in (4.0, 3.5):
ks = info.session[(info.cells >= 10) & (info.sws_pre_min >= thr) & (info.sws_post_min >= thr)]
for drop in (False, True):
res = [awake_classes(S[k], drop) for k in ks]
nsig = [r[0] for r in res]
print(f"variant, SWS threshold {thr} min ({len(ks)} sessions), cells silent in the awake epoch {'left out' if drop else 'kept'}: "
f"signal {sum(nsig)} (291 to 355: {291 <= sum(nsig) <= 355}), nonsignal {sum(r[1] for r in res)} (730 to 892: {730 <= sum(r[1] for r in res) <= 892}), "
f"sessions with 1 to 6 signal {sum(1 <= x <= 6 for x in nsig)}, group sizes >1.8, 1.4-1.8, 1-1.4: {np.sum([r[2] for r in res], axis=0).tolist()} (expected 10, 40, 273)")
print("per session: signal %.2f, nonsignal %.2f here; paper 323/63 = %.2f, 811/63 = %.2f" % (n_sig / len(per), n_non / len(per), 323 / 63, 811 / 63))
included sessions (56): signal 285 (expected 323, 291 to 355: False), nonsignal 730 (expected 811, 730 to 892: True), below lambda_min 484
signal components per session: {2: 1, 3: 12, 4: 14, 5: 12, 6: 6, 7: 4, 9: 4, 10: 2, 11: 1}
sessions with 1 to 6 signal components: 45 of 56
by rat:
cls signal nonsignal below_lmin
rat
15 40 112 62
18 95 257 172
19 44 126 72
20 106 235 178
variant, all 86 sessions: {'signal': 411, 'nonsignal': 1050}
variant, SWS threshold 4.0 min (56 sessions), cells silent in the awake epoch kept: signal 285 (291 to 355: False), nonsignal 730 (730 to 892: True), sessions with 1 to 6 signal 45, group sizes >1.8, 1.4-1.8, 1-1.4: [9, 35, 241] (expected 10, 40, 273)
variant, SWS threshold 4.0 min (56 sessions), cells silent in the awake epoch left out: signal 287 (291 to 355: False), nonsignal 727 (730 to 892: False), sessions with 1 to 6 signal 45, group sizes >1.8, 1.4-1.8, 1-1.4: [9, 35, 243] (expected 10, 40, 273)
variant, SWS threshold 3.5 min (63 sessions), cells silent in the awake epoch kept: signal 315 (291 to 355: True), nonsignal 820 (730 to 892: True), sessions with 1 to 6 signal 51, group sizes >1.8, 1.4-1.8, 1-1.4: [11, 39, 265] (expected 10, 40, 273)
variant, SWS threshold 3.5 min (63 sessions), cells silent in the awake epoch left out: signal 317 (291 to 355: True), nonsignal 817 (730 to 892: True), sessions with 1 to 6 signal 51, group sizes >1.8, 1.4-1.8, 1-1.4: [11, 39, 267] (expected 10, 40, 273) per session: signal 5.09, nonsignal 13.04 here; paper 323/63 = 5.13, 811/63 = 12.87
Result: partially reproduced. Over the 56 included sessions there are 285 signal components (expected 291 to 355) and 730 nonsignal (expected 730 to 892, at the lower edge). Per session these are 5.09 and 13.04, against 5.13 and 12.87 for the paper's totals over 63 sessions, so the shortfall follows the missing sessions. 45 of 56 sessions have 1 to 6 signal components, and 11 have more (up to 11). Leaving out the 12 cells silent in the awake epoch, which cannot be z-transformed, changes the totals little (287 and 727). With a 3.5 min SWS threshold, a variant outside the criterion that gives the paper's total of sessions but not its counts per rat, the totals are 315 and 820, both within tolerance, with 51 of 63 sessions at 1 to 6 signal components and group sizes 11, 39 and 265 (expected 10, 40 and 273).
fig2b-post-greater-than-pre¶
For signal components, the average reactivation strength over SWS is greater in the post epoch than in the pre epoch, in each of the three groups of encoding strength.
def group_tests(d, a, b, label):
rows = []
for g in GROUPS + ["<1"]:
x = d[d.group == g]
if g == "<1":
x = x[x.cls == "nonsignal"]
if len(x) < 2:
rows.append({"group": g, "n": len(x)}); continue
rows.append({"group": g, "n": len(x), f"mean_{a}": x[a].mean(), f"mean_{b}": x[b].mean(),
"mean_diff": (x[b] - x[a]).mean(), "sd_diff": (x[b] - x[a]).std(),
"t": stats.ttest_rel(x[b], x[a]).statistic, "p_t": stats.ttest_rel(x[b], x[a]).pvalue,
"p_wilcoxon": stats.wilcoxon(x[b], x[a]).pvalue})
out = pd.DataFrame(rows)
print(label); print(out.to_string(index=False, float_format=lambda v: f"{v:.3g}")); print()
return out
print("group sizes (expected 10, 40, 273 for >1.8, 1.4-1.8, 1-1.4):", SIG.group.value_counts()[GROUPS].to_dict(),
"| within tolerance:", [abs(SIG.group.value_counts()[g] - e) <= t for g, e, t in zip(GROUPS, (10, 40, 273), (1, 4, 27))])
print("opposite order variant (273, 40, 10 for >1.8, 1.4-1.8, 1-1.4): within tolerance:",
[abs(SIG.group.value_counts()[g] - e) <= t for g, e, t in zip(GROUPS, (273, 40, 10), (27, 4, 1))])
print()
FIG2B = {z: group_tests(comp, f"pre_{z}", f"post_{z}", f"z over {'the whole rest epoch' if z == 'epochz' else 'SWS bins only'}")
for z in ("epochz", "swsz")}
print("by rat (epochz, all signal components): post minus pre, mean and paired t-test")
for rat, x in SIG.groupby("rat"):
print(f" rat {rat}: n={len(x)}, diff {np.mean(x.post_epochz - x.pre_epochz):.3f}, p = {stats.ttest_rel(x.post_epochz, x.pre_epochz).pvalue:.2g}")
fig, ax = plt.subplots(1, 2, figsize=(9, 3))
for a, z in zip(ax, ("epochz", "swsz")):
t = FIG2B[z]
a.bar(np.arange(4) - 0.2, t.iloc[:, 2], 0.4, label="pre")
a.bar(np.arange(4) + 0.2, t.iloc[:, 3], 0.4, label="post")
a.set_xticks(range(4), t.group); a.set_title(z); a.set_ylabel("mean R over SWS"); a.legend()
plt.tight_layout(); plt.show()
print("plotted range: %.3f to %.3f" % (min(FIG2B[z].iloc[:, 2:4].min().min() for z in FIG2B), max(FIG2B[z].iloc[:, 2:4].max().max() for z in FIG2B)))
group sizes (expected 10, 40, 273 for >1.8, 1.4-1.8, 1-1.4): {'>1.8': 9, '1.4-1.8': 35, '1-1.4': 241} | within tolerance: [np.True_, np.False_, np.False_]
opposite order variant (273, 40, 10 for >1.8, 1.4-1.8, 1-1.4): within tolerance: [np.False_, np.False_, np.False_]
z over the whole rest epoch
group n mean_pre_epochz mean_post_epochz mean_diff sd_diff t p_t p_wilcoxon
>1.8 9 0.2 0.471 0.271 0.172 4.74 0.00146 0.00391
1.4-1.8 35 0.146 0.29 0.144 0.198 4.3 0.000135 9.18e-05
1-1.4 241 0.0901 0.114 0.0238 0.0989 3.74 0.000233 8.65e-06
<1 730 -0.00401 -0.00668 -0.00267 0.0569 -1.27 0.206 0.287
z over SWS bins only
group n mean_pre_swsz mean_post_swsz mean_diff sd_diff t p_t p_wilcoxon
>1.8 9 0.165 0.423 0.258 0.145 5.32 0.00071 0.00391
1.4-1.8 35 0.123 0.261 0.138 0.144 5.69 2.18e-06 3.25e-06
1-1.4 241 0.0763 0.1 0.0237 0.0714 5.16 5.23e-07 1.73e-07
<1 730 -0.00284 -0.00537 -0.00253 0.0418 -1.63 0.103 0.313
by rat (epochz, all signal components): post minus pre, mean and paired t-test
rat 15: n=40, diff 0.057, p = 0.00045
rat 18: n=95, diff 0.076, p = 1.2e-06
rat 19: n=44, diff -0.003, p = 0.87
rat 20: n=106, diff 0.037, p = 0.003
plotted range: -0.007 to 0.471
Result: reproduced. Mean reactivation over SWS is higher in post than in pre in all three signal groups, under both z-transform readings and both tests: with z over the whole epoch, paired t p = 0.0015, 0.00014 and 0.00023 for Phi above 1.8, 1.4 to 1.8 and 1 to 1.4 (Wilcoxon 0.0039, 9.2e-05, 8.6e-06); with z over SWS bins, t p = 0.00071, 2.2e-06 and 5.2e-07 (Wilcoxon 0.0039, 3.3e-06, 1.7e-07). Nonsignal components do not differ (p 0.10 to 0.31). The group sizes are 9, 35 and 241 against 10, 40 and 273: the first is within tolerance, the other two are not (about 88% of the paper's, as for the sessions). By rat the effect holds in rats 15, 18 and 20 and not in rat 19 (difference -0.003, p = 0.87). Under the 3.5 min SWS variant of fig1b the group sizes are 11, 39 and 265, all within tolerance.
supp5-sws-only¶
Replay is concentrated in SWS; reactivation strength of signal components shows effects of SWS against non-SWS and of post against pre, and no post-pre effect is seen in rest periods not classified as SWS.
Non-SWS bins are the bins of the rest epochs (see the data check) outside SWS intervals, with Q z-transformed over the whole rest epoch.
x = SIG.dropna(subset=["pre_nonsws", "post_nonsws"])
long = pd.concat([pd.DataFrame({"R": x[f"{ep}_epochz" if st == "SWS" else f"{ep}_nonsws"], "state": st, "epoch": ep})
for st in ("SWS", "nonSWS") for ep in ("pre", "post")])
fit = smf.ols("R ~ C(state) * C(epoch)", data=long).fit()
print(f"signal components with non-SWS bins in both epochs: {len(x)} of {len(SIG)}")
print(anova_lm(fit, typ=2).to_string(float_format=lambda v: f"{v:.3g}"))
long["cid"] = np.tile(x.session.astype(str).values + "_" + x.l.astype(str).values, 4)
print("repeated-measures ANOVA (component as subject, state and epoch within):")
print(AnovaRM(long, "R", "cid", within=["state", "epoch"]).fit().anova_table.to_string(float_format=lambda v: f"{v:.3g}"))
print(long.groupby(["state", "epoch"]).R.agg(["mean", "sem", "count"]).to_string(float_format=lambda v: f"{v:.4f}"))
print("\nnon-SWS, post against pre over signal components: mean diff %.4f, paired t p = %.3g, Wilcoxon p = %.3g" % (
(x.post_nonsws - x.pre_nonsws).mean(), stats.ttest_rel(x.post_nonsws, x.pre_nonsws).pvalue,
stats.wilcoxon(x.post_nonsws, x.pre_nonsws).pvalue))
print("SWS, post against pre over the same components: mean diff %.4f, paired t p = %.3g" % (
(x.post_epochz - x.pre_epochz).mean(), stats.ttest_rel(x.post_epochz, x.pre_epochz).pvalue))
for g in GROUPS:
y = x[x.group == g]
print(f" group {g}: n={len(y)}, non-SWS post-pre {np.mean(y.post_nonsws - y.pre_nonsws):.4f} (t p = {stats.ttest_rel(y.post_nonsws, y.pre_nonsws).pvalue:.2g}), "
f"SWS post-pre {np.mean(y.post_epochz - y.pre_epochz):.4f}")
# variant: rest epochs limited to the 30 min next to the awake epoch (the Methods give 20 to 30 min)
v30 = []
for k, l in zip(SIG.session, SIG.l):
r = RX[k]; row = {}
for ep in ("pre", "post"):
e = r[ep]; mid = (e["edges"][:-1] + e["edges"][1:]) / 2
near = (mid >= S[k]["wake"][0] - 1.8e6) if ep == "pre" else (mid <= S[k]["wake"][1] + 1.8e6)
R = e["R_all_epochz"][l]
row[f"{ep}_sws"] = R[near & e["sws"]].mean() if (near & e["sws"]).any() else np.nan
row[f"{ep}_non"] = R[near & ~e["sws"]].mean() if (near & ~e["sws"]).any() else np.nan
v30.append(row)
v30 = pd.DataFrame(v30).dropna()
long30 = pd.concat([pd.DataFrame({"R": v30[f"{ep}_{st}"], "state": st, "epoch": ep}) for st in ("sws", "non") for ep in ("pre", "post")])
a30 = anova_lm(smf.ols("R ~ C(state) * C(epoch)", data=long30).fit(), typ=2)
print(f"\nvariant, 30 min rest epochs ({len(v30)} components): state p = {a30.loc['C(state)', 'PR(>F)']:.3g}, epoch p = {a30.loc['C(epoch)', 'PR(>F)']:.3g}, "
f"interaction p = {a30.loc['C(state):C(epoch)', 'PR(>F)']:.3g}; means SWS pre {v30.pre_sws.mean():.4f} post {v30.post_sws.mean():.4f}, non-SWS pre {v30.pre_non.mean():.4f} post {v30.post_non.mean():.4f}; "
f"non-SWS post-pre paired t p = {stats.ttest_rel(v30.post_non, v30.pre_non).pvalue:.3g}, Wilcoxon p = {stats.wilcoxon(v30.post_non, v30.pre_non).pvalue:.3g}")
nons = pd.DataFrame([{"session": k, "pre_nonsws_min": (~r["pre"]["sws"]).sum() / 600, "post_nonsws_min": (~r["post"]["sws"]).sum() / 600} for k, r in RX.items()])
print("non-SWS minutes per session: pre %.1f to %.1f, post %.1f to %.1f" % (nons.pre_nonsws_min.min(), nons.pre_nonsws_min.max(), nons.post_nonsws_min.min(), nons.post_nonsws_min.max()))
signal components with non-SWS bins in both epochs: 285 of 285
sum_sq df F PR(>F)
C(state) 0.114 1 6.08 0.0138
C(epoch) 0.229 1 12.2 0.0005
C(state):C(epoch) 0.0927 1 4.94 0.0265
Residual 21.3 1.14e+03 NaN NaN
repeated-measures ANOVA (component as subject, state and epoch within):
F Value Num DF Den DF Pr > F
state 9.83 1 284 0.0019
epoch 31.4 1 284 4.94e-08
state:epoch 19.2 1 284 1.64e-05
mean sem count
state epoch
SWS post 0.1469 0.0090 285
pre 0.1005 0.0076 285
nonSWS post 0.1088 0.0079 285
pre 0.0985 0.0078 285
non-SWS, post against pre over signal components: mean diff 0.0103, paired t p = 0.0416, Wilcoxon p = 0.0179
SWS, post against pre over the same components: mean diff 0.0464, paired t p = 5.83e-09
group >1.8: n=9, non-SWS post-pre 0.0183 (t p = 0.49), SWS post-pre 0.2714
group 1.4-1.8: n=35, non-SWS post-pre 0.0234 (t p = 0.21), SWS post-pre 0.1439
group 1-1.4: n=241, non-SWS post-pre 0.0081 (t p = 0.13), SWS post-pre 0.0238
variant, 30 min rest epochs (285 components): state p = 0.0178, epoch p = 0.000472, interaction p = 0.0261; means SWS pre 0.1010 post 0.1487, non-SWS pre 0.0998 post 0.1104; non-SWS post-pre paired t p = 0.0438, Wilcoxon p = 0.0207
non-SWS minutes per session: pre 10.0 to 46.3, post 10.0 to 50.1
Result: partially reproduced. Over the 285 signal components, the ANOVA finds the post against pre and the SWS against non-SWS effects, with SWS above non-SWS and post above pre. The four means come from the same component, so a repeated-measures ANOVA (component as subject) is the test used: state p = 0.0019, epoch p = 4.9e-08, interaction p = 1.6e-05, so both main effects pass p < 0.01. An ordinary two-way ANOVA, which treats the four means as independent, gives state p = 0.014, above 0.01, and epoch p = 0.0005. In non-SWS bins post exceeds pre slightly but significantly (difference 0.0103 against 0.0464 in SWS; paired t p = 0.042, Wilcoxon p = 0.018), so the claim's third part, no post-pre effect outside SWS, fails under both readings; within each signal group the non-SWS difference is not significant (p 0.13 to 0.49). With rest epochs limited to 30 min next to the awake epoch the non-SWS result is the same (non-SWS t p = 0.044; ordinary ANOVA state p = 0.018, epoch p = 0.00047).
supp6-pre-post-sws-alike¶
Pre and post epoch SWS do not differ in average SWS duration, average firing rates, rates of SPWRs and delta waves, or LFP power in the delta and spindle bands.
The data set holds no LFP and no SPWR or delta-wave times, so only SWS duration and firing rates can be tested. Paired t-tests over the included sessions (and over cells for one reading).
rows, cell_rows = [], []
for k in INC:
s = S[k]
r = {"session": k, "rat": s["rat"]}
for ep, iv in (("pre", s["sws_pre"]), ("post", s["sws_post"])):
dur = (iv[:, 1] - iv[:, 0]) / 1000
r[f"{ep}_total_s"] = dur.sum(); r[f"{ep}_episode_s"] = dur.mean(); r[f"{ep}_episodes"] = len(dur)
cnt = np.zeros(pfc6.n_cells(s))
for a, b in iv:
m = (s["spikes_t"] >= a) & (s["spikes_t"] < b)
np.add.at(cnt, s["spikes_id"][m] - 1, 1)
rate = cnt / dur.sum()
r[f"{ep}_rate"] = rate.mean()
for i, v in enumerate(rate):
cell_rows.append({"session": k, "cell": i, "epoch": ep, "rate": v})
rows.append(r)
sw = pd.DataFrame(rows)
cells = pd.DataFrame(cell_rows).pivot_table(index=["session", "cell"], columns="epoch", values="rate")
def paired(a, b, label):
t = stats.ttest_rel(b, a)
print(f"{label}: pre {np.mean(a):.3f}, post {np.mean(b):.3f}, mean diff {np.mean(b - a):.3f} "
f"({100 * np.mean(b - a) / np.mean(a):+.1f}%), paired t = {t.statistic:.2f}, p = {t.pvalue:.3g}, n = {len(a)} -> difference found: {t.pvalue < 0.05}")
paired(sw.pre_total_s, sw.post_total_s, "SWS duration, total per epoch (s)")
paired(sw.pre_episode_s, sw.post_episode_s, "SWS duration, mean episode (s)")
paired(sw.pre_rate, sw.post_rate, "SWS firing rate, session mean (Hz)")
paired(cells["pre"], cells["post"], "SWS firing rate, cells (Hz)")
wil = [stats.wilcoxon(sw[f"post_{m}"], sw[f"pre_{m}"]).pvalue for m in ("total_s", "episode_s", "rate")]
wil.append(stats.wilcoxon(cells["post"], cells["pre"]).pvalue)
print("Wilcoxon (variant): total %.3g, episode %.3g, session rate %.3g, cell rate %.3g" % tuple(wil))
print("by rat, session-mean rate post minus pre (Hz):", sw.assign(d=sw.post_rate - sw.pre_rate).groupby("rat").d.mean().round(3).to_dict())
print("by rat, total SWS post minus pre (s):", sw.assign(d=sw.post_total_s - sw.pre_total_s).groupby("rat").d.mean().round(0).to_dict())
SWS duration, total per epoch (s): pre 629.671, post 623.049, mean diff -6.622 (-1.1%), paired t = -0.16, p = 0.87, n = 56 -> difference found: False
SWS duration, mean episode (s): pre 202.792, post 222.227, mean diff 19.435 (+9.6%), paired t = 0.83, p = 0.412, n = 56 -> difference found: False
SWS firing rate, session mean (Hz): pre 3.583, post 3.576, mean diff -0.007 (-0.2%), paired t = -0.09, p = 0.927, n = 56 -> difference found: False
SWS firing rate, cells (Hz): pre 3.501, post 3.528, mean diff 0.027 (+0.8%), paired t = 0.46, p = 0.647, n = 1499 -> difference found: False
Wilcoxon (variant): total 0.935, episode 0.478, session rate 0.903, cell rate 0.00133
by rat, session-mean rate post minus pre (Hz): {15: 0.27, 18: 0.162, 19: -0.234, 20: -0.096}
by rat, total SWS post minus pre (s): {15: -105.0, 18: -4.0, 19: -55.0, 20: 55.0}
Result: reproduced, for the measures the data set holds. Pre and post SWS do not differ by paired t-test over the 56 sessions in total SWS time (630 and 623 s, p = 0.87), mean episode duration (203 and 222 s, p = 0.41), or mean firing rate per session (3.58 and 3.58 Hz, p = 0.93); over 1499 cells the rates do not differ either (p = 0.65), though a Wilcoxon test over cells, a variant, gives p = 0.0013 for a 0.8% difference. SPWR rate, delta-wave counts and LFP power cannot be tested: the data set holds no LFP or event times.
fig2b-encoding-correlation¶
Across signal components, post epoch reactivation strength correlates positively with encoding strength (r squared 0.61).
Reading (a), all post SWS bins. Reading (b), bins outside SPWRs, needs SPWR times, which the data set does not hold.
for z in ("epochz", "swsz"):
r, p = stats.pearsonr(SIG.phi, SIG[f"post_{z}"])
print(f"{z}: n = {len(SIG)}, r = {r:.3f}, r^2 = {r * r:.3f} (expected 0.61, 0.549 to 0.671: {abs(r * r - 0.61) <= 0.061}), p = {p:.3g} (below 1e-30: {p < 1e-30})")
sp = stats.spearmanr(SIG.phi, SIG[f"post_{z}"])
print(f" Spearman rho = {sp.statistic:.3f}; without the {int(np.sum(SIG.phi > 3))} components with Phi > 3: r^2 = {stats.pearsonr(SIG.phi[SIG.phi <= 3], SIG[f'post_{z}'][SIG.phi <= 3])[0] ** 2:.3f}")
print("by rat (epochz):", {rat: round(stats.pearsonr(x.phi, x.post_epochz)[0] ** 2, 3) for rat, x in SIG.groupby("rat")})
fig, ax = plt.subplots(figsize=(4, 3))
ax.scatter(SIG.phi, SIG.post_epochz, s=6, label="post"); ax.scatter(SIG.phi, SIG.pre_epochz, s=6, label="pre")
ax.set_xlabel("Phi"); ax.set_ylabel("mean R over SWS (epochz)"); ax.legend(); plt.tight_layout(); plt.show()
print("plotted ranges: Phi %.2f to %.2f, R %.3f to %.3f" % (SIG.phi.min(), SIG.phi.max(), min(SIG.pre_epochz.min(), SIG.post_epochz.min()), max(SIG.pre_epochz.max(), SIG.post_epochz.max())))
epochz: n = 285, r = 0.620, r^2 = 0.384 (expected 0.61, 0.549 to 0.671: False), p = 1.23e-31 (below 1e-30: True)
Spearman rho = 0.534; without the 0 components with Phi > 3: r^2 = 0.384
swsz: n = 285, r = 0.646, r^2 = 0.417 (expected 0.61, 0.549 to 0.671: False), p = 4.53e-35 (below 1e-30: True)
Spearman rho = 0.561; without the 0 components with Phi > 3: r^2 = 0.417
by rat (epochz): {15: np.float64(0.196), 18: np.float64(0.427), 19: np.float64(0.031), 20: np.float64(0.5)}
plotted ranges: Phi 1.00 to 2.14, R -0.127 to 0.889
Result: not reproduced. Mean post SWS reactivation correlates positively with Phi over the 285 signal components (r = 0.620, p = 1.2e-31 with z over the epoch; r = 0.646, p = 4.5e-35 with z over SWS), but r squared is 0.384 and 0.417, below the expected 0.549 to 0.671. The r itself (0.62 to 0.65) is close to the 0.61 the paper prints as r squared. Reading (b), bins outside SPWRs, cannot be run.
fig2b-pre-weaker¶
In pre epoch SWS the relation between reactivation strength and encoding strength is significantly weaker than in post, in slope and in correlation.
def steiger(r12, r13, r23, n):
"""Steiger (1980) Z for two correlations sharing variable 1."""
z12, z13 = np.arctanh(r12), np.arctanh(r13)
rm2 = (r12 ** 2 + r13 ** 2) / 2
f = min((1 - r23) / (2 * (1 - rm2)), 1.0)
h = (1 - f * rm2) / (1 - rm2)
Z = (z12 - z13) * np.sqrt(n - 3) / np.sqrt(2 * (1 - r23) * h)
return Z, 2 * stats.norm.sf(abs(Z))
rng = np.random.default_rng(20090701)
for z in ("epochz", "swsz"):
n = len(SIG)
long = pd.concat([pd.DataFrame({"R": SIG[f"{ep}_{z}"].values, "phi": SIG.phi.values, "post": int(ep == "post")}) for ep in ("pre", "post")])
fit = smf.ols("R ~ phi * post", data=long).fit()
b_pre = np.polyfit(SIG.phi, SIG[f"pre_{z}"], 1)[0]; b_post = np.polyfit(SIG.phi, SIG[f"post_{z}"], 1)[0]
boots = []
idx = np.arange(n)
for _ in range(2000):
i = rng.choice(idx, n)
boots.append(np.polyfit(SIG.phi.values[i], SIG[f"post_{z}"].values[i], 1)[0] - np.polyfit(SIG.phi.values[i], SIG[f"pre_{z}"].values[i], 1)[0])
boots = np.array(boots)
zb = (b_post - b_pre) / boots.std()
r_post = stats.pearsonr(SIG.phi, SIG[f"post_{z}"])[0]; r_pre = stats.pearsonr(SIG.phi, SIG[f"pre_{z}"])[0]
r_pp = stats.pearsonr(SIG[f"pre_{z}"], SIG[f"post_{z}"])[0]
zf = (np.arctanh(r_post) - np.arctanh(r_pre)) / np.sqrt(2 / (n - 3))
Zs, ps = steiger(r_post, r_pre, r_pp, n)
print(f"{z}: slope pre {b_pre:.4f}, post {b_post:.4f}")
print(f" (a) interaction in OLS: coef {fit.params['phi:post']:.4f}, p = {fit.pvalues['phi:post']:.3g}")
print(f" (b) bootstrap of the slope difference (2000 resamples of components): {np.mean(boots <= 0):.4f} of resamples <= 0; "
f"normal approximation z = {zb:.2f}, p = {2 * stats.norm.sf(abs(zb)):.3g}")
print(f" correlation pre r = {r_pre:.3f}, post r = {r_post:.3f}, r(pre, post) = {r_pp:.3f}")
print(f" (a) Fisher z, independent: z = {zf:.2f}, p = {2 * stats.norm.sf(abs(zf)):.3g}")
print(f" (b) Steiger, dependent: Z = {Zs:.2f}, p = {ps:.3g}")
epochz: slope pre 0.1676, post 0.4336 (a) interaction in OLS: coef 0.2660, p = 2.4e-08 (b) bootstrap of the slope difference (2000 resamples of components): 0.0000 of resamples <= 0; normal approximation z = 6.46, p = 1.04e-10 correlation pre r = 0.283, post r = 0.620, r(pre, post) = 0.579 (a) Fisher z, independent: z = 5.16, p = 2.52e-07 (b) Steiger, dependent: Z = 7.20, p = 6.18e-13
swsz: slope pre 0.1498, post 0.3975 (a) interaction in OLS: coef 0.2477, p = 2.58e-10 (b) bootstrap of the slope difference (2000 resamples of components): 0.0000 of resamples <= 0; normal approximation z = 7.28, p = 3.33e-13 correlation pre r = 0.319, post r = 0.646, r(pre, post) = 0.659 (a) Fisher z, independent: z = 5.21, p = 1.93e-07 (b) Steiger, dependent: Z = 7.91, p = 2.53e-15
Result: partially reproduced. In pre SWS the slope on Phi is lower (0.168 against 0.434 with z over the epoch) and so is the correlation (r = 0.283 against 0.620). The correlation difference passes p < 1e-5 under every reading (Fisher z p = 2.5e-07, Steiger p = 6.2e-13; with z over SWS 1.9e-07 and 2.5e-15). The slope difference is significant but does not reach p < 1e-20 under any reading (interaction p = 2.4e-08 and 2.6e-10; bootstrap normal approximation p = 1.0e-10 and 3.3e-13).
supp8-other-tetrode-pairs¶
Blocked: the data set does not give the tetrode of each cell, so same-tetrode pairs cannot be found (see the summary).
fig2e-skewness¶
Post epoch SWS reactivation strength distributions are heavy tailed, more skewed than pre epoch ones, most markedly for components of higher encoding strength.
for z in ("epochz", "swsz"):
print(f"{z}: skewness of R over SWS bins, post minus pre, one-sample t-test per group")
means = []
for g in GROUPS:
x = SIG[SIG.group == g]
dlt = x[f"post_skew_{z}"] - x[f"pre_skew_{z}"]
t = stats.ttest_1samp(dlt, 0)
means.append(dlt.mean())
print(f" {g}: n = {len(x)}, mean diff {dlt.mean():.3f} (s.e.m. {dlt.sem():.3f}), t = {t.statistic:.2f}, p = {t.pvalue:.3g}; "
f"skew pre {x[f'pre_skew_{z}'].mean():.2f}, post {x[f'post_skew_{z}'].mean():.2f}")
print(f" mean difference rises with Phi (>1.8 > 1.4-1.8 > 1-1.4): {means[0] > means[1] > means[2]}; "
f"(>1.8 > 1.4-1.8): {means[0] > means[1]}")
x = comp[comp.cls == "nonsignal"]
print(f" nonsignal: mean diff {(x[f'post_skew_{z}'] - x[f'pre_skew_{z}']).mean():.3f}, p = {stats.ttest_1samp(x[f'post_skew_{z}'] - x[f'pre_skew_{z}'], 0).pvalue:.3g}")
epochz: skewness of R over SWS bins, post minus pre, one-sample t-test per group >1.8: n = 9, mean diff -1.192 (s.e.m. 2.321), t = -0.51, p = 0.622; skew pre 11.45, post 10.26 1.4-1.8: n = 35, mean diff 3.443 (s.e.m. 1.246), t = 2.76, p = 0.00916; skew pre 5.23, post 8.68 1-1.4: n = 241, mean diff -0.151 (s.e.m. 0.552), t = -0.27, p = 0.785; skew pre 5.51, post 5.36 mean difference rises with Phi (>1.8 > 1.4-1.8 > 1-1.4): False; (>1.8 > 1.4-1.8): False nonsignal: mean diff 0.132, p = 0.728 swsz: skewness of R over SWS bins, post minus pre, one-sample t-test per group >1.8: n = 9, mean diff -1.516 (s.e.m. 2.339), t = -0.65, p = 0.535; skew pre 10.61, post 9.09 1.4-1.8: n = 35, mean diff 3.248 (s.e.m. 0.987), t = 3.29, p = 0.00234; skew pre 5.14, post 8.39 1-1.4: n = 241, mean diff -0.440 (s.e.m. 0.521), t = -0.84, p = 0.399; skew pre 5.42, post 4.98 mean difference rises with Phi (>1.8 > 1.4-1.8 > 1-1.4): False; (>1.8 > 1.4-1.8): False nonsignal: mean diff 0.142, p = 0.698
Result: not reproduced. The post minus pre skewness of R over SWS bins is positive and significant only for Phi 1.4 to 1.8 (3.44, p = 0.0092 with z over the epoch; 3.25, p = 0.0023 with z over SWS). For Phi above 1.8 it is negative (-1.19, p = 0.62; -1.52, p = 0.54), and for 1 to 1.4 near zero (p = 0.79 and 0.40), so neither reading of the groups has every tested group positive with p < 0.05, and the difference does not rise with Phi.
fig2f-tail-contribution¶
About 40 to 50% of the net reactivation (post minus pre) comes from events with reactivation strength beyond the 99th percentile.
Cumulative contribution C(r) = (1/N) sum of R values up to r (the integral of u P(u) up to r), post minus pre; the fraction beyond the post 99th percentile r99 is 1 - [C_post(r99) - C_pre(r99)] / [mean_post - mean_pre]. Reading (a) pools the bins of a group's components; reading (b) computes the fraction per component and averages.
def tail_fraction(post, pre, q=99):
r99 = np.percentile(post, q)
net = post.mean() - pre.mean()
below = post[post <= r99].sum() / len(post) - pre[pre <= r99].sum() / len(pre)
return 1 - below / net, net
def bins_of(z, ep, rows):
key = "R_sws_epochz" if z == "epochz" else "R_sws_swsz"
return [RX[k][ep][key][l] for k, l in zip(rows.session, rows.l)]
FIG2F = {}
for z in ("epochz", "swsz"):
for g in GROUPS + ["nonsignal"]:
rows = SIG[SIG.group == g] if g != "nonsignal" else comp[comp.cls == "nonsignal"]
post, pre = bins_of(z, "post", rows), bins_of(z, "pre", rows)
fa, net = tail_fraction(np.concatenate(post), np.concatenate(pre))
per_c = np.array([tail_fraction(a, b)[0] for a, b in zip(post, pre)])
FIG2F[(z, g)] = (fa, np.mean(per_c))
print(f"{z} {g:>9}: n = {len(rows)}, (a) pooled {100 * fa:.1f}% (net {net:.4f}), (b) per-component mean {100 * np.mean(per_c):.1f}% "
f"(median {100 * np.median(per_c):.1f}%, net negative in {int(np.sum([a.mean() < b.mean() for a, b in zip(post, pre)]))} components, largest |ratio| {100 * np.abs(per_c).max():.0f}%)")
for z in ("epochz", "swsz"):
for i, rd in enumerate(("(a)", "(b)")):
vals = [FIG2F[(z, g)][i] for g in GROUPS]
print(f"{z} {rd}: signal groups in 36 to 55%: {sum(0.36 <= v <= 0.55 for v in vals)} of 3")
epochz >1.8: n = 9, (a) pooled 52.2% (net 0.2881), (b) per-component mean -80.3% (median 51.7%, net negative in 0 components, largest |ratio| 1145%) epochz 1.4-1.8: n = 35, (a) pooled 56.1% (net 0.1412), (b) per-component mean -31.0% (median 38.8%, net negative in 6 components, largest |ratio| 2715%) epochz 1-1.4: n = 241, (a) pooled 16.1% (net 0.0261), (b) per-component mean 96.9% (median 41.7%, net negative in 99 components, largest |ratio| 41025%)
epochz nonsignal: n = 730, (a) pooled 201.7% (net -0.0035), (b) per-component mean 57.5% (median 39.3%, net negative in 377 components, largest |ratio| 14924%) swsz >1.8: n = 9, (a) pooled 51.0% (net 0.2738), (b) per-component mean 46.3% (median 47.4%, net negative in 0 components, largest |ratio| 68%) swsz 1.4-1.8: n = 35, (a) pooled 60.6% (net 0.1438), (b) per-component mean 45.8% (median 47.1%, net negative in 6 components, largest |ratio| 149%) swsz 1-1.4: n = 241, (a) pooled 36.6% (net 0.0242), (b) per-component mean 5459.1% (median 42.2%, net negative in 93 components, largest |ratio| 1319000%)
swsz nonsignal: n = 730, (a) pooled 127.7% (net -0.0032), (b) per-component mean -164.0% (median 33.9%, net negative in 370 components, largest |ratio| 131670%) epochz (a): signal groups in 36 to 55%: 1 of 3 epochz (b): signal groups in 36 to 55%: 0 of 3 swsz (a): signal groups in 36 to 55%: 2 of 3 swsz (b): signal groups in 36 to 55%: 2 of 3
Result: partially reproduced. Pooling bins by group (reading a), the share of the net post minus pre reactivation from bins above the post 99th percentile is 52.2%, 56.1% and 16.1% (Phi above 1.8, 1.4 to 1.8, 1 to 1.4) with z over the epoch, and 51.0%, 60.6% and 36.6% with z over SWS: one and two of three groups in 36 to 55%. Averaging per component (reading b) gives two of three in range with z over SWS (46.3%, 45.8%) and none with z over the epoch, as components whose net is near zero give extreme ratios (99 of 241 components in the lowest group have a negative net). Those means (for example 5459%, with one component's ratio above 10^6 %) are not well defined; the count of groups in range rests on readings (a) and on the two highest groups of (b) with z over SWS, whose per-component ratios stay within 149%.
supp10-no-coactivation¶
During sleep, reactivation strengths of different simultaneously recorded signal components tend not to peak together, with a trough at zero lag in their cross-correlograms, while in the awake epoch they coactivate more than by chance.
Cross-correlograms over lags of -1 to 1 s (100 ms bins), computed within each SWS interval so that no pair of bins spans a gap, averaged over all pairs of signal components of a session and then over pairs pooled. Readings: (a) Pearson cross-correlation (chance 0), (b) mean product divided by the product of the means (chance 1).
MAXLAG = 10
cc = {(rd, ep, z): [] for rd in "ab" for ep in ("post", "awake") for z in ("epochz", "swsz")}
for k, r in RX.items():
sig_l = np.flatnonzero(r["lam"] > r["lmax"])
if len(sig_l) < 2:
continue
mask = r["post"]["sws"]
st, sp_ = pfc6.runs(mask)
lens = sp_ - st
cst = np.r_[0, np.cumsum(lens)[:-1]]
for i in range(len(sig_l)):
for j in range(i + 1, len(sig_l)):
a, b = sig_l[i], sig_l[j]
for z, key in (("epochz", "R_sws_epochz"), ("swsz", "R_sws_swsz")):
x, y = r["post"][key][a], r["post"][key][b]
cc[("a", "post", z)].append(pfc6.xcorr_runs(x, y, cst, cst + lens, MAXLAG))
cc[("b", "post", z)].append(pfc6.xcorr_ratio_runs(x, y, cst, cst + lens, MAXLAG))
x, y = r["R_awake"][a], r["R_awake"][b]
for z in ("epochz", "swsz"):
if z == "epochz":
cc[("a", "awake", z)].append(pfc6.xcorr_runs(x, y, [0], [len(x)], MAXLAG))
cc[("b", "awake", z)].append(pfc6.xcorr_ratio_runs(x, y, [0], [len(x)], MAXLAG))
lags = np.arange(-MAXLAG, MAXLAG + 1) * 100
far = (np.abs(lags) >= 200) & (np.abs(lags) <= 500)
print("pairs of signal components:", len(cc[("a", "post", "epochz")]), "from", sum((r["lam"] > r["lmax"]).sum() >= 2 for r in RX.values()), "sessions")
for rd, chance in (("a", 0.0), ("b", 1.0)):
aw = np.nanmean(cc[(rd, "awake", "epochz")], axis=0)
for z in ("epochz", "swsz"):
m = np.nanmean(cc[(rd, "post", z)], axis=0)
sem = np.nanstd(cc[(rd, "post", z)], axis=0) / np.sqrt(len(cc[(rd, "post", z)]))
t0 = stats.ttest_1samp(np.array(cc[(rd, "post", z)])[:, MAXLAG], chance)
print(f"({rd}) post SWS {z}: zero lag {m[MAXLAG]:.4f} (s.e.m. {sem[MAXLAG]:.4f}; vs chance {chance}: t = {t0.statistic:.2f}, p = {t0.pvalue:.2g}), "
f"mean at 200-500 ms {m[far].mean():.4f}, min over lags at {lags[np.argmin(m)]} ms; "
+ (f"below chance: {m[MAXLAG] < chance}, below 200-500 ms: {m[MAXLAG] < m[far].mean()}" if rd == "a" else
"no comparison made: the ratio divides by mean R near zero for many components"))
t1 = stats.ttest_1samp(np.array(cc[(rd, "awake", "epochz")])[:, MAXLAG], chance)
print(f"({rd}) awake: zero lag {aw[MAXLAG]:.4f} (vs chance: t = {t1.statistic:.2f}, p = {t1.pvalue:.2g}), at 200-500 ms {aw[far].mean():.4f}; above chance: {aw[MAXLAG] > chance}")
fig, ax = plt.subplots(1, 2, figsize=(9, 3))
for a, rd in zip(ax, "ab"):
a.plot(lags, np.nanmean(cc[(rd, "post", "epochz")], axis=0), label="post SWS")
a.plot(lags, np.nanmean(cc[(rd, "awake", "epochz")], axis=0), label="awake")
a.set_xlabel("lag (ms)"); a.set_title(f"reading ({rd})"); a.legend()
plt.tight_layout(); plt.show()
print("plotted ranges: (a) %.3f to %.3f, (b) %.3f to %.3f" % (
min(np.nanmean(cc[("a", e, "epochz")], axis=0).min() for e in ("post", "awake")), max(np.nanmean(cc[("a", e, "epochz")], axis=0).max() for e in ("post", "awake")),
min(np.nanmean(cc[("b", e, "epochz")], axis=0).min() for e in ("post", "awake")), max(np.nanmean(cc[("b", e, "epochz")], axis=0).max() for e in ("post", "awake"))))
pairs of signal components: 704 from 56 sessions (a) post SWS epochz: zero lag -0.0035 (s.e.m. 0.0033; vs chance 0.0: t = -1.06, p = 0.29), mean at 200-500 ms 0.0017, min over lags at 0 ms; below chance: True, below 200-500 ms: True (a) post SWS swsz: zero lag -0.0099 (s.e.m. 0.0028; vs chance 0.0: t = -3.53, p = 0.00045), mean at 200-500 ms 0.0016, min over lags at 0 ms; below chance: True, below 200-500 ms: True (a) awake: zero lag 0.0329 (vs chance: t = 8.24, p = 8.5e-16), at 200-500 ms 0.0072; above chance: True (b) post SWS epochz: zero lag 0.8576 (s.e.m. 32.8811; vs chance 1.0: t = -0.00, p = 1), mean at 200-500 ms 3.3880, min over lags at -800 ms; no comparison made: the ratio divides by mean R near zero for many components (b) post SWS swsz: zero lag -337.1756 (s.e.m. 430.3435; vs chance 1.0: t = -0.79, p = 0.43), mean at 200-500 ms -11.4161, min over lags at 0 ms; no comparison made: the ratio divides by mean R near zero for many components (b) awake: zero lag 9.1428 (vs chance: t = 3.88, p = 0.00012), at 200-500 ms 2.0042; above chance: True
plotted ranges: (a) -0.004 to 0.033, (b) -70.323 to 30.992
Result: reproduced. Under reading (a), Pearson cross-correlation averaged over 704 pairs of signal components, post SWS has its minimum at zero lag, below chance (-0.0035 with z over the epoch, p = 0.29 against 0; -0.0099 with z over SWS, p = 0.00045) and below the mean at 200 to 500 ms (0.0017 and 0.0016), and in the awake epoch the zero-lag value is above chance (0.0329, p = 8.5e-16). Reading (b), divided by the product of the means, is unstable because mean R over SWS is near zero for many components (zero-lag s.e.m. 33 and 430); no comparison with chance is made under it in sleep, and the outcome rests on reading (a).
supp11-shuffle-control¶
Reactivation strength computed after permuting cell identities in each time bin, which keeps global firing rate fluctuations, does not differ between pre and post and is more than an order of magnitude smaller than for the actual patterns.
The mean of R_l(t) over all permutations of the cells within a bin has a closed form
(pfc6.shuffle_expected_R), which is the limit of the mean over shuffles. It is used for every
component; on one session it is compared with the mean over 1000 actual shuffles.
rng = np.random.default_rng(11)
sh = []
for k, l, g in zip(SIG.session, SIG.l, SIG.group):
r = RX[k]
row = {"session": k, "group": g}
for ep in ("pre", "post"):
for z, qk, rk in (("epochz", "q_sws_epochz", "R_sws_epochz"), ("swsz", "q_sws_swsz", "R_sws_swsz")):
row[f"{ep}_{z}_shuf"] = pfc6.shuffle_expected_R(r[ep][qk], r["vec"][:, [l]])[0].mean()
row[f"{ep}_{z}_actual"] = r[ep][rk][l].mean()
sh.append(row)
sh = pd.DataFrame(sh)
for z in ("epochz", "swsz"):
print(f"{z}:")
for g in GROUPS:
x = sh[sh.group == g]
a, b = x[f"pre_{z}_shuf"], x[f"post_{z}_shuf"]
act = (x[f"pre_{z}_actual"].mean() + x[f"post_{z}_actual"].mean()) / 2
shm = (a.mean() + b.mean()) / 2
print(f" {g}: n = {len(x)}, shuffled pre {a.mean():.5f}, post {b.mean():.5f}; paired t p = {stats.ttest_rel(b, a).pvalue:.3g}, "
f"Wilcoxon p = {stats.wilcoxon(b, a).pvalue:.3g}; actual mean {act:.4f}; |shuffled| / actual = {abs(shm) / act:.4f}")
# Monte Carlo check on the first included session
k0 = INC[0]; r0 = RX[k0]; l0 = int(np.flatnonzero(r0["lam"] > r0["lmax"])[0])
q0 = r0["post"]["q_sws_epochz"]
mc = pfc6.shuffle_R(q0, r0["vec"][:, [l0]], 1000, rng)[0]
ex = pfc6.shuffle_expected_R(q0, r0["vec"][:, [l0]])[0]
print(f"session {k0}, component {l0}, post SWS ({q0.shape[1]} bins): mean over 1000 shuffles {mc.mean():.5f}, closed form {ex.mean():.5f}; "
f"per-bin correlation {np.corrcoef(mc, ex)[0, 1]:.4f}")
epochz: >1.8: n = 9, shuffled pre 0.01621, post 0.02037; paired t p = 0.0307, Wilcoxon p = 0.0117; actual mean 0.3356; |shuffled| / actual = 0.0545 1.4-1.8: n = 35, shuffled pre 0.01105, post 0.01217; paired t p = 0.13, Wilcoxon p = 0.023; actual mean 0.2182; |shuffled| / actual = 0.0532 1-1.4: n = 241, shuffled pre 0.00610, post 0.00598; paired t p = 0.847, Wilcoxon p = 0.215; actual mean 0.1020; |shuffled| / actual = 0.0592 swsz: >1.8: n = 9, shuffled pre 0.01495, post 0.01807; paired t p = 0.0339, Wilcoxon p = 0.0195; actual mean 0.2939; |shuffled| / actual = 0.0562 1.4-1.8: n = 35, shuffled pre 0.00891, post 0.01083; paired t p = 0.067, Wilcoxon p = 0.0174; actual mean 0.1918; |shuffled| / actual = 0.0514 1-1.4: n = 241, shuffled pre 0.00509, post 0.00496; paired t p = 0.759, Wilcoxon p = 0.408; actual mean 0.0882; |shuffled| / actual = 0.0570
session 150628, component 0, post SWS (6606 bins): mean over 1000 shuffles -0.00860, closed form -0.00815; per-bin correlation 0.7213
Result: partially reproduced. Shuffled reactivation is 5 to 6% of the actual in every signal group (|shuffled| / actual 0.053 to 0.059 with z over the epoch, 0.051 to 0.057 with z over SWS bins), more than an order of magnitude smaller. It does not differ between pre and post for Phi 1 to 1.4 (t p = 0.85 and 0.76) or, by t-test, for 1.4 to 1.8 (p = 0.13 and 0.067), but for Phi above 1.8 post is higher (0.0204 against 0.0162, t p = 0.031; with z over SWS p = 0.034), and the Wilcoxon reading finds a difference in the two highest groups (p 0.012 to 0.023). The closed-form shuffle mean agrees with 1000 shuffles on one session (-0.0082 and -0.0086).
fig6b-down-up¶
Reactivation and SPWR occurrence relative to putative DOWN to UP transitions, with DOWN states defined as decreases of neural activity lasting at least 80 ms, are comparable with the delta wave results.
The SPWR part cannot be tested (no ripple times). DOWN states, from the spikes of all cells of the session in 10 ms bins inside SWS: (a) at least 80 ms with no spike, (b) at least 80 ms with the count below a fifth of the SWS mean count per bin. Transitions are the ends of DOWN states. Event-triggered averages of R_l over -4 to 4 s (100 ms bins of the rest epoch, Q z-transformed over the whole epoch, since bins outside SWS fall in the window), per component; each bin against the baseline (-4 to -2 s and 2 to 4 s) by a t-test over components at p < 0.001 (Figure 6a legend, as the criterion says).
def down_up_times(s, iv, reading):
out = []
t = np.sort(s["spikes_t"])
for a, b in iv:
nb = int((b - a) // 10)
if nb < 8:
continue
cnt = np.bincount(((t[(t >= a) & (t < a + nb * 10)] - a) // 10).astype(int), minlength=nb)
low = cnt == 0 if reading == "a" else cnt < cnt.mean() / 5
st, sp_ = pfc6.runs(low)
keep = (sp_ - st) >= 8
out.append(a + sp_[keep] * 10.0)
return np.concatenate(out) if out else np.zeros(0)
mean10 = [sum(np.sum((S[k]["spikes_t"] >= a) & (S[k]["spikes_t"] < b)) for a, b in iv) / (pfc6.sws_minutes(iv) * 6000)
for k in RX for iv in (S[k]["sws_pre"], S[k]["sws_post"])]
print(f"mean population spike count per 10 ms SWS bin, per session and epoch: {min(mean10):.2f} to {max(mean10):.2f}; "
f"a fifth of it is below one spike in {sum(m / 5 < 1 for m in mean10)} of {len(mean10)}, so reading (b) marks the same bins as (a)")
W = 40
lag_s = np.arange(-W, W + 1) * 0.1
base = (np.abs(lag_s) >= 2 - 1e-9)
FIG6B = {}
for reading in ("a", "b"):
eta = {"pre": [], "post": []}; eta_non = {"pre": [], "post": []}; n_ev = {"pre": 0, "post": 0}
for k, r in RX.items():
s = S[k]
for ep, iv in (("pre", s["sws_pre"]), ("post", s["sws_post"])):
ev = down_up_times(s, iv, reading)
edges = r[ep]["edges"]
kb = ((ev - edges[0]) // pfc6.BIN_MS).astype(int)
kb = kb[(kb - W >= 0) & (kb + W < len(edges) - 1)]
n_ev[ep] += len(kb)
if len(kb) == 0:
continue
idx = kb[:, None] + np.arange(-W, W + 1)[None, :]
R = r[ep]["R_all_epochz"]
for l in range(r["n"]):
if r["lam"][l] > r["lmax"]:
eta[ep].append(R[l][idx].mean(0))
elif r["lam"][l] > r["lmin"]:
eta_non[ep].append(R[l][idx].mean(0))
print(f"reading ({reading}): transitions used, pre {n_ev['pre']}, post {n_ev['post']}")
for ep in ("post", "pre"):
for lab, E in (("signal", np.array(eta[ep])), ("nonsignal", np.array(eta_non[ep]))):
d = E - E[:, base].mean(1, keepdims=True)
t, p = stats.ttest_1samp(d, 0)
sigb = (p < 0.001) & ~base
m = d.mean(0)
pk = lag_s[np.argmax(np.where(base, -np.inf, m))]
pre_sig = lag_s[sigb & (lag_s < 0) & (m > 0)]
FIG6B[(reading, ep, lab)] = (m, sigb)
print(f" {ep} {lab}: n = {len(E)}, peak above baseline {m.max():.4f} at {pk:+.1f} s (baseline {E[:, base].mean():.4f}); "
f"significant bins {int(sigb.sum())}, of which positive before the transition: {np.round(pre_sig, 1).tolist()}")
mpost, mpre = FIG6B[(reading, "post", "signal")][0], FIG6B[(reading, "pre", "signal")][0]
win = (lag_s >= -0.55) & (lag_s <= -0.25)
print(f" largest rise above baseline: post {mpost.max():.4f}, pre {mpre.max():.4f} (pre lower: {mpre.max() < mpost.max()}); "
f"mean rise over -0.5 to -0.3 s: post {mpost[win].mean():.4f}, pre {mpre[win].mean():.4f} (pre lower: {mpre[win].mean() < mpost[win].mean()})")
print(f" event-triggered average without baseline subtraction, -0.5 to -0.3 s: post {np.mean(eta['post'], axis=0)[win].mean():.4f}, pre {np.mean(eta['pre'], axis=0)[win].mean():.4f}")
Ep, Eq = np.array(eta["post"]), np.array(eta["pre"])
dp = (Ep[:, win].mean(1) - Ep[:, base].mean(1)) - (Eq[:, win].mean(1) - Eq[:, base].mean(1))
print(f" per component, rise over -0.5 to -0.3 s post minus pre: mean {dp.mean():.4f}, paired t p = {stats.ttest_1samp(dp, 0).pvalue:.3g}")
for ep in ("post", "pre"):
m, sigb = FIG6B[("a", ep, "signal")]
print(f"{ep} signal, significant bins (s) and their mean above baseline:", [(round(float(x), 1), round(float(y), 4)) for x, y in zip(lag_s[sigb], m[sigb])])
print("pre and post signal averages minus baseline, -1 to 1 s (reading a):")
print(pd.DataFrame({"lag_s": np.round(lag_s, 1), "pre": FIG6B[("a", "pre", "signal")][0], "post": FIG6B[("a", "post", "signal")][0]})
.query("-1 <= lag_s <= 1").to_string(index=False, float_format=lambda v: f"{v:.4f}"))
fig, ax = plt.subplots(figsize=(5, 3))
for ep in ("post", "pre"):
ax.plot(lag_s, FIG6B[("a", ep, "signal")][0], label=f"{ep} signal")
ax.plot(lag_s, FIG6B[("a", "post", "nonsignal")][0], label="post nonsignal")
ax.axvline(0, color="k", lw=0.5); ax.set_xlabel("time from DOWN-UP transition (s)"); ax.set_ylabel("R minus baseline"); ax.legend()
plt.tight_layout(); plt.show()
print("plotted range: %.4f to %.4f" % (min(FIG6B[("a", e, c)][0].min() for e in ("pre", "post") for c in ("signal", "nonsignal") if (e, c) != ("pre", "nonsignal")),
max(FIG6B[("a", e, c)][0].max() for e in ("pre", "post") for c in ("signal", "nonsignal") if (e, c) != ("pre", "nonsignal"))))
mean population spike count per 10 ms SWS bin, per session and epoch: 0.21 to 2.22; a fifth of it is below one spike in 112 of 112, so reading (b) marks the same bins as (a)
reading (a): transitions used, pre 14596, post 15415 post signal: n = 285, peak above baseline 0.1317 at -0.3 s (baseline 0.1650); significant bins 4, of which positive before the transition: [-0.5, -0.4, -0.3] post nonsignal: n = 730, peak above baseline 0.0126 at +1.6 s (baseline -0.0086); significant bins 0, of which positive before the transition: [] pre signal: n = 285, peak above baseline 0.1447 at +0.2 s (baseline 0.1193); significant bins 1, of which positive before the transition: [] pre nonsignal: n = 730, peak above baseline 0.0233 at -1.4 s (baseline -0.0073); significant bins 0, of which positive before the transition: [] largest rise above baseline: post 0.1317, pre 0.1447 (pre lower: False); mean rise over -0.5 to -0.3 s: post 0.1065, pre 0.0634 (pre lower: True) event-triggered average without baseline subtraction, -0.5 to -0.3 s: post 0.2715, pre 0.1826 per component, rise over -0.5 to -0.3 s post minus pre: mean 0.0432, paired t p = 0.152
reading (b): transitions used, pre 14596, post 15415 post signal: n = 285, peak above baseline 0.1317 at -0.3 s (baseline 0.1650); significant bins 4, of which positive before the transition: [-0.5, -0.4, -0.3] post nonsignal: n = 730, peak above baseline 0.0126 at +1.6 s (baseline -0.0086); significant bins 0, of which positive before the transition: [] pre signal: n = 285, peak above baseline 0.1447 at +0.2 s (baseline 0.1193); significant bins 1, of which positive before the transition: [] pre nonsignal: n = 730, peak above baseline 0.0233 at -1.4 s (baseline -0.0073); significant bins 0, of which positive before the transition: [] largest rise above baseline: post 0.1317, pre 0.1447 (pre lower: False); mean rise over -0.5 to -0.3 s: post 0.1065, pre 0.0634 (pre lower: True) event-triggered average without baseline subtraction, -0.5 to -0.3 s: post 0.2715, pre 0.1826 per component, rise over -0.5 to -0.3 s post minus pre: mean 0.0432, paired t p = 0.152 post signal, significant bins (s) and their mean above baseline: [(-0.5, 0.0952), (-0.4, 0.0928), (-0.3, 0.1317), (0.0, -0.0592)] pre signal, significant bins (s) and their mean above baseline: [(0.2, 0.1447)] pre and post signal averages minus baseline, -1 to 1 s (reading a): lag_s pre post -1.0000 -0.0102 0.0585 -0.9000 0.1281 0.0261 -0.8000 0.0082 0.0343 -0.7000 0.0525 0.0472 -0.6000 0.0597 0.0625 -0.5000 0.0672 0.0952 -0.4000 0.0534 0.0928 -0.3000 0.0695 0.1317 -0.2000 0.0252 0.0805 -0.1000 -0.0214 -0.0221 0.0000 -0.0278 -0.0592 0.1000 0.0519 0.1067 0.2000 0.1447 0.0216 0.3000 0.0488 0.0346 0.4000 0.0072 0.0353 0.5000 -0.0521 0.0794 0.6000 0.0205 0.0234 0.7000 -0.0330 0.0214 0.8000 -0.0216 0.0003 0.9000 0.0211 -0.0059 1.0000 -0.0047 0.0427
plotted range: -0.0592 to 0.1447
Result: partially reproduced, for the reactivation part, under an assumption. In post SWS the signal components' reactivation rises significantly above baseline 0.5 to 0.3 s before DOWN to UP transitions (three bins at p < 0.001, 15415 transitions); the fourth significant bin is a dip below baseline at the transition (0 s, -0.059). Nonsignal components have no significant bin. In the same window the pre rise is lower (0.063 against 0.107), but per component the post minus pre difference is not significant (0.043, p = 0.15). Pre has one significant bin, at +0.2 s after the transition, where its average (0.145) is above the post one at any lag (largest 0.132 at -0.3 s). So the post rise before the transition holds, and the pre average is lower only in the window where post is significant, without a significant difference. The assumption is that the pre average is lower only if both comparisons agree or the per-component difference is significant; neither holds, so this part is not established. Comparing by the window alone would make the claim reproduced. Both definitions of DOWN states give the same transitions: the mean population count per 10 ms SWS bin is 0.21 to 2.22 spikes, so a fifth of it is below one spike in every session and reading (b) marks exactly the bins without spikes, as (a) does. The SPWR part is not tested.
supp12-behavior-unchanged¶
Rule acquisition does not change trial duration or the running speed along the trajectory.
The ten learning sessions and their learning trials are those of the data set's description (see learning-sessions). Trial length is the trial's end minus its start. Trials before the learning trial are "before", the learning trial and those after it are "after". Speed along the trajectory is shown in ten bins of the fraction of the trial's path length travelled (the data set holds x, y positions without the maze's geometry, so the paper's linearization is not repeated).
tl = []
speed = []
for k, lp in LEARN.items():
s = S[k]
b = s["behavior"]
p = s["pos"]
for i, (a, e) in enumerate(b[:, :2]):
tl.append({"session": k, "rat": s["rat"], "after": i + 1 >= lp, "len_s": (e - a) / 1000})
m = (p[:, 0] >= a) & (p[:, 0] <= e)
if m.sum() < 5:
continue
xy = p[m, 1:]
step = np.hypot(*np.diff(xy, axis=0).T)
dt = np.diff(p[m, 0]) / 1000
frac = np.cumsum(step) / max(step.sum(), 1e-9)
bins = np.minimum((frac * 10).astype(int), 9)
v = step / dt
for bi in range(10):
if np.any(bins == bi):
speed.append({"session": k, "after": i + 1 >= lp, "bin": bi, "speed": v[bins == bi].mean()})
tl = pd.DataFrame(tl)
sm = tl.groupby(["session", "after"]).len_s.mean().unstack()
t = stats.ttest_rel(sm[True], sm[False])
print(f"(b) session means, paired: before {sm[False].mean():.2f} s, after {sm[True].mean():.2f} s, diff {np.mean(sm[True] - sm[False]):+.2f} s, "
f"t = {t.statistic:.2f}, p = {t.pvalue:.3g}, n = {len(sm)} -> difference found: {t.pvalue < 0.05}")
t2 = stats.ttest_ind(tl.len_s[tl.after], tl.len_s[~tl.after])
print(f"(a) trials pooled, unpaired: before {tl.len_s[~tl.after].mean():.2f} s (n = {int((~tl.after).sum())}), after {tl.len_s[tl.after].mean():.2f} s "
f"(n = {int(tl.after.sum())}), t = {t2.statistic:.2f}, p = {t2.pvalue:.3g} -> difference found: {t2.pvalue < 0.05}")
print("Mann-Whitney on pooled trials (variant): p = %.3g; medians before %.2f, after %.2f s" % (
stats.mannwhitneyu(tl.len_s[tl.after], tl.len_s[~tl.after]).pvalue, tl.len_s[~tl.after].median(), tl.len_s[tl.after].median()))
print("per session (mean s before, after):\n", sm.round(2).to_string())
lpa = lp_df[lp_df.perf == "a"].set_index("session").lp
v = tl.assign(after=[i >= lpa[k] for k, i in zip(tl.session, tl.groupby("session").cumcount() + 1)])
sv = v.groupby(["session", "after"]).len_s.mean().unstack()
tv = stats.ttest_rel(sv[True], sv[False])
print(f"variant, learning points of learning-sessions reading (a) in the same ten sessions: before {sv[False].mean():.2f} s, after {sv[True].mean():.2f} s, "
f"paired t p = {tv.pvalue:.3g}, n = {len(sv)}")
sp = pd.DataFrame(speed).groupby(["session", "after", "bin"]).speed.mean().unstack("after")
print("\nspeed by path-fraction bin (position units/s), session means, paired t-test per bin:")
for bi in range(10):
x = sp.xs(bi, level="bin").dropna()
p = stats.ttest_rel(x[True], x[False]).pvalue
print(f" bin {bi}: before {x[False].mean():.1f}, after {x[True].mean():.1f}, p = {p:.3g}, Bonferroni over 10 bins {min(1, 10 * p):.2g} (n = {len(x)})")
(b) session means, paired: before 5.48 s, after 8.23 s, diff +2.74 s, t = 1.54, p = 0.158, n = 10 -> difference found: False (a) trials pooled, unpaired: before 5.35 s (n = 119), after 6.81 s (n = 157), t = 1.09, p = 0.278 -> difference found: False Mann-Whitney on pooled trials (variant): p = 0.62; medians before 4.30, after 4.30 s per session (mean s before, after): after False True session 150628 6.20 4.78 150630 4.18 8.65 150707 5.31 3.92 181012 6.03 22.72 181020 3.97 4.51 190214 5.07 5.53 190228 9.60 16.92 201222 4.45 4.59 201227 6.16 5.06 201229 3.87 5.59 variant, learning points of learning-sessions reading (a) in the same ten sessions: before 5.65 s, after 8.23 s, paired t p = 0.19, n = 10 speed by path-fraction bin (position units/s), session means, paired t-test per bin: bin 0: before 23.3, after 24.9, p = 0.0342, Bonferroni over 10 bins 0.34 (n = 10) bin 1: before 33.7, after 37.1, p = 0.373, Bonferroni over 10 bins 1 (n = 10) bin 2: before 40.7, after 40.7, p = 0.996, Bonferroni over 10 bins 1 (n = 10) bin 3: before 43.2, after 41.5, p = 0.652, Bonferroni over 10 bins 1 (n = 10) bin 4: before 47.7, after 45.7, p = 0.501, Bonferroni over 10 bins 1 (n = 10) bin 5: before 55.8, after 52.8, p = 0.399, Bonferroni over 10 bins 1 (n = 10) bin 6: before 62.7, after 59.1, p = 0.498, Bonferroni over 10 bins 1 (n = 10) bin 7: before 63.4, after 63.9, p = 0.831, Bonferroni over 10 bins 1 (n = 10) bin 8: before 59.2, after 58.9, p = 0.898, Bonferroni over 10 bins 1 (n = 10) bin 9: before 46.0, after 44.8, p = 0.699, Bonferroni over 10 bins 1 (n = 10)
Result: reproduced. Trial length does not differ before and after the learning trial in the ten learning sessions: session means 5.48 against 8.23 s (paired t p = 0.16, n = 10), and over trials 5.35 against 6.81 s (unpaired t p = 0.28; medians both 4.30 s). The test has little power (a 2.7 s difference with n = 10 is not significant). The sessions and learning trials are those of the data set's description, since learning-sessions did not recover them; with the learning points of its reading (a) in the same ten sessions the result is the same (paired t p = 0.19). Speed along the path does not differ in 9 of 10 bins; the first bin is faster after learning (24.9 against 23.3 units/s, p = 0.034 uncorrected, 0.34 after Bonferroni over the ten bins). The bins are fractions of each trial's path length, not the paper's linearized maze.
Claims blocked for missing LFP and event times¶
The data set holds spikes, cell types, behaviour, position, and the awake and SWS intervals. It holds no LFP and no SPWR, delta-wave or spindle times, and no tetrode of each cell. The claims fig5a-spwr-bins, fig5b-spwr-triggered, fig5c-peak-40ms, fig5d-population-asymmetric, fig5e-autocorrelogram, replay-in-delta-spindle-periods, fig6a-delta-waves, fig6c-spindles, fig6d-reactivation-spwr-correlation, fig7b-post-spwr-after-learning, fig7b-controls, fig7c-choice-point and fig7d-two-factors need ripple times or the prefrontal LFP, and supp8-other-tetrode-pairs needs the tetrode of each cell.
print(f"notebook run time {(time.time() - T0) / 60:.1f} min")
notebook run time 2.0 min
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome |
|---|---|
| sessions-and-cells | not reproduced (56 sessions and 1499 cells under the inclusion rule) |
| rules-learned | not reproduced (right 3, lit 4, left 3, dark 1 rats) |
| learning-sessions | not reproduced by the definition (14 to 34 sessions); the description lists ten sessions from four rats |
| fig1b-signal-components | partially reproduced (285 signal, 730 nonsignal; same per session) |
| fig2b-post-greater-than-pre | reproduced |
| supp5-sws-only | partially reproduced (post exceeds pre in non-SWS bins too) |
| supp6-pre-post-sws-alike | reproduced for SWS duration and firing rates; LFP measures not testable |
| fig2b-encoding-correlation | not reproduced (r squared 0.38 to 0.42) |
| fig2b-pre-weaker | partially reproduced (slope p not below 1e-20) |
| supp8-other-tetrode-pairs | blocked: data missing |
| fig2e-skewness | not reproduced |
| fig2f-tail-contribution | partially reproduced |
| supp10-no-coactivation | reproduced |
| supp11-shuffle-control | partially reproduced |
| fig5a, fig5b, fig5c, fig5d, fig5e | blocked: data missing |
| replay-in-delta-spindle-periods, fig6a, fig6c, fig6d | blocked: data missing |
| fig6b-down-up | partially reproduced for reactivation (pre lower in the window before the transition, not at its peak, and not significantly) |
| fig7b, fig7b-controls, fig7c, fig7d | blocked: data missing |
| supp12-behavior-unchanged | reproduced |
How the analysis differs from the paper's. The rest epochs' bounds are not in the data set and are taken from the awake interval and the position record. The data set holds 56 sessions that pass the inclusion rule, where the paper analysed 63. Each claim with an open z-transform was run with z over the whole rest epoch and over its SWS bins; no outcome changes between them. The shuffle control uses the exact mean over all permutations instead of 1000 shuffles per bin, checked against 1000 shuffles on one session. The supp5 ANOVA is a repeated-measures ANOVA over components. Speed along the trajectory uses the fraction of path length, not the paper's linearized maze.
What was found about the data set. It holds spikes, cell types, behaviour, position and the awake and SWS intervals of 86 sessions, but no LFP, no SPWR, delta or spindle times and no tetrode per cell, so 14 of the 29 claims cannot be tested. It holds 14 sessions of rat 15 against the paper's 16, and the inclusion rule keeps 10 of them. Two of the ten learning sessions it lists do not pass the 4 min SWS rule with its own SWS intervals. Sessions without SWS store a [0, 0] interval or an empty file. Position units and the rest epochs' bounds are not given, two sessions have implausible cell types (all unknown, all interneurons), and in five sessions spikes run on after the position record.