Reproducing Nair et al. (2018) from CRCNS data set bf-1¶
Paper. Nair J, Klaassen A-L, Arato J, Vyssotski AL, Harvey M, Rainer G (2018). Basal forebrain contributes to default mode network regulation. PNAS 115:1352-1357. doi:10.1073/pnas.1712431115. The claims file was written from the version of record and its supplement (SI Results, SI Materials and Methods, legends of Figs. S1 to S7).
Data set. CRCNS bf-1, "Basal forebrain LFP recordings from rats in home cage and during arena
exploration", doi:10.6080/K0MK6B2Q, read from the public bucket
crcnsarchive under bf-1/live/. Version: unversioned. CRCNS data sets have no versions; the
paper cites the data set's DOI only, and the bucket holds one copy (live, no old), so that copy
is the one read. It is pinned by the MD5 sum of its single data file, from the data set's
checksums.md5: basal_forebrain_lfp.tar.gz a8acbf0400260f8395257778ed04dd8c (22,301,364
bytes); the notebook checks the sum of what it streams against it. The data descriptor is version
0.5 (Jan 23, 2018).
What the data set holds. BF LFP at 400 Hz from six rats, as 2-s segments with artifact segments removed, one movement-sensor value per segment, and two home-cage and two arena sessions per rat. It holds no ACC or VC recordings, no single units, no object sessions, no behaviour scoring or video, no day labels, no identification of bipolar electrodes and no recording across the transfer between arena and home cage.
| Claim | Tested here |
|---|---|
| fig3a-bf-gamma-reduction | tested |
| fig3a-each-hemisphere-significant | tested |
| fig3b-arena-gamma-unrelated-to-movement | tested |
| fig3b-home-cage-gamma-increases-with-movement | tested |
| fig3b-matched-movement | tested |
| fig3c-home-cage-above-both-sessions | tested in part (day factor dropped, as the criterion says; the session order is not in the data set) |
| fig3d-gamma-rises-with-time-in-arena | tested |
| figs2-bf-low-frequency-bands | tested |
| figs3-unimodal-states | tested |
| figs5-arena-excursions | tested |
| bipolar-gamma-reduction | blocked: metadata_insufficient (no rat is identified as bipolar) |
| fig3c-consistent-over-days | blocked: metadata_insufficient (no day labels) |
| fig3d-rapid-switch-to-home-cage | blocked: data_missing (no recording across the transfer) |
| behavior-arena-time-budget, behavior-home-cage-less-exploration, behavior-movement-speed, fig3f-grooming-quiet-arena, figs6-object-more-stationary | blocked: data_missing (no video or behaviour scoring) |
| fig4-object-exploration-suppressed | blocked: data_missing (no object sessions) |
| fig5a-phase-locked-neurons, fig5cd-pli-max-at-gamma, spike-width-phase-locked | blocked: data_missing (no single units) |
| fig6a-acc-gamma-home-cage, fig6a-vc-gamma-unchanged, fig6b-coherence, fig6c-granger-bf-to-acc, fig6d-trivariate-granger, figs2-acc-vc-low-frequency-bands | blocked: data_missing (no ACC or VC) |
Method details common to the tested claims are in the docstring of bflfp.py, with the passage of
the paper or the claims file each comes from.
import time
import warnings
import diptest
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.anova import AnovaRM
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import bflfp
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
sums = bflfp.checksums()
recordings, archive = bflfp.load_recordings()
print(f"data set {bflfp.DATASET_ID} ({bflfp.DATASET_VERSION}, doi:{bflfp.DATASET_DOI}): "
f"{len(recordings)} files streamed in {time.time() - started:.0f} s")
print("checksums.md5 lists:", sums)
print("MD5 of the streamed archive:", archive["md5"], f"({archive['bytes']} bytes);",
"matches checksums.md5:", archive["md5"] == sums[bflfp.ARCHIVE])
data set bf-1 (unversioned, doi:10.6080/K0MK6B2Q): 24 files streamed in 3 s
checksums.md5 lists: {'basal_forebrain_lfp.tar.gz': 'a8acbf0400260f8395257778ed04dd8c', 'docs/crcns_bf-1_data_description.docx': '3c5ed5f3c32b0104e0800708399e0da2', 'docs/crcns_bf-1_data_description.pdf': '958e4e1c07db4af4b9a8308581d5d332', 'filelist.txt': 'f0d111ec1369a817591aacf2ebcc8e46'}
MD5 of the streamed archive: a8acbf0400260f8395257778ed04dd8c (22301364 bytes); matches checksums.md5: True
Data check¶
rows = []
for r in recordings:
starts = r["t"][0]
dt = np.diff(r["t"], axis=0)
span = starts[-1] + 2.0
rows.append({"rat": r["rat"], "condition": r["condition"], "session": r["session"],
"variables": ",".join(r["variables"]), "samples_per_segment": r["lfp"].shape[0],
"segments": r["lfp"].shape[1], "minutes_kept": r["lfp"].shape[1] * 2 / 60,
"first_start_s": starts[0], "last_end_s": span,
"removed_pct": 100 * (1 - r["lfp"].shape[1] * 2 / span),
"sample_dt_ms": f"{1000 * dt.min():.4f}-{1000 * dt.max():.4f}",
"segment_steps_s": ",".join(f"{x:g}" for x in np.unique(np.round(np.diff(starts), 3))[:5]),
"within_segment_contiguous": bool(np.allclose(r["t"][-1] - r["t"][0], 799 / 400)),
"lfp_min": r["lfp"].min(), "lfp_max": r["lfp"].max(), "lfp_sd": r["lfp"].std(),
"mov_min": r["mov"].min(), "mov_median": np.median(r["mov"]), "mov_max": r["mov"].max(),
"md5": r["md5"][:8]})
info = pd.DataFrame(rows)
print("Files by rat and condition:")
print(pd.crosstab(info.rat, info.condition).to_string())
print("\nVariables in every file:", info.variables.unique().tolist())
print("Condition labels in the file names:", sorted(info.condition.unique()))
print(info.drop(columns=["variables"]).round(2).to_string(index=False))
print("\nRemoved segments, %% of each file's span: %.1f to %.1f (SI: 1-4%% of epochs rejected)"
% (info.removed_pct.min(), info.removed_pct.max()))
Files by rat and condition: condition Arena HCage rat 1 2 2 2 2 2 3 2 2 4 2 2 5 2 2 6 2 2 Variables in every file: ['LFP_data,T_data,mov_data'] Condition labels in the file names: ['Arena', 'HCage'] rat condition session samples_per_segment segments minutes_kept first_start_s last_end_s removed_pct sample_dt_ms segment_steps_s within_segment_contiguous lfp_min lfp_max lfp_sd mov_min mov_median mov_max md5 1 Arena 1 800 283 9.43 0.0 620.0 8.71 2.5000-2.5000 2,4,6,8,10 True -350.07 315.89 69.09 19.54 441.42 887.43 3884e20b 1 Arena 2 800 281 9.37 0.0 620.0 9.35 2.5000-2.5000 2,4,6 True -301.28 282.71 66.73 190.96 558.23 945.24 950fdd71 1 HCage 1 800 948 31.60 0.0 2000.0 5.20 2.5000-2.5000 2,4,6,8,10 True -566.13 992.54 90.42 16.82 220.02 848.53 7d43c330 1 HCage 2 800 354 11.80 0.0 794.0 10.83 2.5000-2.5000 2,4,6,8,10 True -397.00 347.20 81.17 18.36 327.37 835.87 5a959f9d 2 Arena 1 800 302 10.07 0.0 620.0 2.58 2.5000-2.5000 2,4,8 True -346.23 372.96 64.61 24.68 402.51 941.77 0bce913e 2 Arena 2 800 255 8.50 0.0 526.0 3.04 2.5000-2.5000 2,4,8 True -479.75 430.21 63.14 48.48 499.75 995.71 018d34c7 2 HCage 1 800 623 20.77 0.0 1400.0 11.00 2.5000-2.5000 2,4,6,8,10 True -376.99 335.59 77.42 18.26 327.58 667.31 70a3cf3c 2 HCage 2 800 317 10.57 0.0 700.0 9.43 2.5000-2.5000 2,4,6,8 True -366.81 326.60 71.43 107.69 285.76 769.97 327145d6 3 Arena 1 800 332 11.07 0.0 688.0 3.49 2.5000-2.5000 2,4,6 True -352.93 342.38 78.04 103.32 370.09 911.24 b6ca8558 3 Arena 2 800 301 10.03 4.0 698.0 13.75 2.5000-2.5000 2,4,6,8 True -311.88 303.27 77.30 25.50 416.09 796.76 c2acda89 3 HCage 1 800 449 14.97 2.0 1000.0 10.20 2.5000-2.5000 2,4,6,8,10 True -370.54 352.23 84.04 19.63 223.57 677.81 104c6fb3 3 HCage 2 800 291 9.70 6.0 650.0 10.46 2.5000-2.5000 2,4,8,12,14 True -380.24 401.03 84.38 19.63 200.24 677.81 ebab2252 4 Arena 1 800 354 11.80 0.0 758.0 6.60 2.5000-2.5000 2,4,6,10 True -239.91 242.51 62.18 15.26 91.48 318.56 afa67f08 4 Arena 2 800 335 11.17 0.0 694.0 3.46 2.5000-2.5000 2,4,8 True -245.80 291.40 60.22 19.94 111.91 307.82 5648f61a 4 HCage 1 800 359 11.97 0.0 758.0 5.28 2.5000-2.5000 2,4,6,8 True -245.80 279.63 67.68 3.94 60.07 204.94 a8f68eae 4 HCage 2 800 289 9.63 0.0 700.0 17.43 2.5000-2.5000 2,4,6,8,14 True -214.48 271.77 68.53 5.27 47.24 198.02 54a177ac 5 Arena 1 800 312 10.40 2.0 704.0 11.36 2.5000-2.5000 2,4,6,10 True -391.54 358.53 88.39 6.98 84.54 223.36 0c270ea7 5 Arena 2 800 323 10.77 0.0 694.0 6.92 2.5000-2.5000 2,4,6,8 True -450.07 352.72 86.12 18.86 103.21 382.62 39de92d5 5 HCage 1 800 290 9.67 0.0 758.0 23.48 2.5000-2.5000 2,4,6,8,10 True -487.20 493.32 107.05 3.03 61.46 210.53 633e0bab 5 HCage 2 800 286 9.53 0.0 650.0 12.00 2.5000-2.5000 2,4,6,8,12 True -514.48 522.58 115.25 5.24 74.42 230.08 41742613 6 Arena 1 800 291 9.70 0.0 660.0 11.82 2.5000-2.5000 2,4,6,8,10 True -277.33 275.38 67.15 20.51 314.32 599.66 396870a9 6 Arena 2 800 350 11.67 0.0 780.0 10.26 2.5000-2.5000 2,4,6,10,22 True -281.30 281.30 67.09 41.09 346.79 693.43 e240c424 6 HCage 1 800 370 12.33 0.0 898.0 17.59 2.5000-2.5000 2,4,6,8,12 True -310.50 298.79 70.88 17.01 187.56 803.56 4ab28b37 6 HCage 2 800 261 8.70 0.0 600.0 13.00 2.5000-2.5000 2,4,6,10,14 True -314.53 326.17 76.95 16.37 304.30 776.14 ab2e1dbd Removed segments, % of each file's span: 2.6 to 23.5 (SI: 1-4% of epochs rejected)
Every file holds LFP_data, T_data and mov_data with 800 samples per segment, and sample
intervals of exactly 2.5 ms (400 Hz), as the descriptor says. Segment starts step by multiples of
2 s, so segments are consecutive 2-s blocks with some removed, as the descriptor says. Each file's
clock starts at or near 0, so no two files share a clock: an arena file cannot be placed in time
relative to a home-cage file. The descriptor names the home-cage condition Hcage; the files use
HCage. Arena files span 8.8 to 13 min and keep 8.5 to 11.8 min of segments; home-cage files keep 8.7 to
31.6 min. The share of removed
segments (2.6% to 23.5% of a file's span) is larger than the SI's 1-4% of rejected epochs, so the
data set's cleaning differs from the paper's rule or removed more. The LFP has no stated unit; its
range (a few hundred, standard deviations of 60 to 115) is consistent with microvolts, the unit of
the paper's Figure 4 axis. The movement values have no unit (the descriptor says dimensionless); in
rats 4 and 5 they are about a quarter of those in the other rats, so they are only compared within
a rat.
shared = bflfp.shared_segments(recordings)
print("Pairs of files that share identical LFP segments:")
print(shared.to_string(index=False) if len(shared) else "none")
by_key = {bflfp.key(r): r for r in recordings}
for _, row in shared.iterrows():
a, b = by_key[row["first"]], by_key[row["second"]]
pos = {a["lfp"][:, j].tobytes(): j for j in range(a["lfp"].shape[1])}
pairs = [(pos.get(b["lfp"][:, j].tobytes()), j) for j in range(b["lfp"].shape[1])]
matched = [(i, j) for i, j in pairs if i is not None]
same_mov = np.mean([a["mov"][i] == b["mov"][j] for i, j in matched])
print(f"{row['second']}: {len(matched)} of {b['lfp'].shape[1]} segments are segments of {row['first']} "
f"({a['lfp'].shape[1]} segments), shifted by {row['time_offsets_s']} s; they cover "
f"{a['t'][0, matched[0][0]]:.0f} to {a['t'][0, matched[-1][0]] + 2:.0f} s of it; "
f"movement values equal on {100 * same_mov:.0f}% of the shared segments")
lag = [np.mean([a["mov"][i + k] == b["mov"][j] for i, j in matched if 0 <= i + k < len(a["mov"])]) for k in (-1, 1)]
print(f" movement values equal at a shift of -1 segment: {100 * lag[0]:.0f}%, +1 segment: {100 * lag[1]:.0f}%")
DUPLICATES = {row["second"] for _, row in shared.iterrows()}
print("Left out as a second deposit of the same recording:", DUPLICATES)
Pairs of files that share identical LFP segments:
first second shared_segments time_offsets_s
(3, HCage, 1) (3, HCage, 2) 287 [100.0]
(3, 'HCage', 2): 287 of 291 segments are segments of (3, 'HCage', 1) (449 segments), shifted by [np.float64(100.0)] s; they cover 106 to 750 s of it; movement values equal on 100% of the shared segments
movement values equal at a shift of -1 segment: 0%, +1 segment: 0%
Left out as a second deposit of the same recording: {(3, 'HCage', 2)}
Rat 3's home-cage session 2 is not a separate recording: all but 4 of its 291 segments are segments of rat 3's home-cage session 1, 100 s later on that file's clock, with the same LFP. It is counted once: the analysis uses session 1 and leaves session 2 out. Rat 3 so has one home-cage recording. Where a test pairs sessions, the pair of rat 3's arena session 2 is dropped (11 pairs), and the result with the files as delivered (12 pairs) is shown as a variant.
# Spectra: line noise at 50 Hz (mains in Switzerland) lies inside the 30-80 Hz gamma band.
rows = []
for r in recordings:
f, p = bflfp.spectrum(r["lfp"])
mean = p.mean(axis=1)
near = (f >= 45) & (f <= 55) & (f != 50)
band = (f >= 20) & (f <= 100)
rows.append({"rat": r["rat"], "condition": r["condition"], "session": r["session"],
"peak_20_100_Hz": f[band][np.argmax(mean[band])],
"power_50Hz_over_45_55": mean[f == 50][0] / mean[near].mean()})
lines = pd.DataFrame(rows)
print(lines.round(2).to_string(index=False))
print("Frequency resolution of a 2-s epoch: %.2f Hz; highest frequency %.0f Hz" % (f[1] - f[0], f[-1]))
rat condition session peak_20_100_Hz power_50Hz_over_45_55 1 Arena 1 57.5 0.80 1 Arena 2 20.0 1.06 1 HCage 1 50.0 1.59 1 HCage 2 51.5 1.18 2 Arena 1 59.5 1.88 2 Arena 2 50.0 3.13 2 HCage 1 50.0 26.05 2 HCage 2 57.5 0.79 3 Arena 1 20.0 1.10 3 Arena 2 20.0 1.03 3 HCage 1 48.5 1.02 3 HCage 2 51.0 0.94 4 Arena 1 50.0 1.16 4 Arena 2 50.0 2.12 4 HCage 1 46.5 1.18 4 HCage 2 48.0 1.05 5 Arena 1 57.0 0.99 5 Arena 2 61.0 1.00 5 HCage 1 51.0 1.22 5 HCage 2 53.5 1.07 6 Arena 1 29.5 1.01 6 Arena 2 25.0 1.52 6 HCage 1 50.0 1.26 6 HCage 2 51.5 1.01 Frequency resolution of a 2-s epoch: 0.50 Hz; highest frequency 200 Hz
The spectral peak between 20 and 100 Hz lies at 46.5 to 57.5 Hz in every home-cage file, the BF gamma the paper describes; in arena files it lies at 50 to 61 Hz in half and at 20 to 30 Hz in the others. One file, rat 2 home cage session 1, has a narrow line at 50 Hz, 26 times the power of its neighbours; a few others have a smaller one (up to 3 times). The paper does not mention line noise; the 30-80 Hz band includes 50 Hz, and the effect of leaving out 49-51 Hz is shown under fig3a-bf-gamma-reduction.
# Do the movement values share the clock of the LFP? Head movement adds broadband muscle and motion
# artefact to an LFP, seen best above the LFP bands (100-200 Hz here). Within each file the
# correlation of that power with the movement values should peak at zero shift of the movement
# values (shifts in segments, by position in the file); theta power (6-10 Hz) is shown alongside.
shifts = range(-5, 6)
rows = []
for r in recordings:
f, p = bflfp.spectrum(r["lfp"])
for name, (lo, hi) in (("100-200 Hz", (100, 200)), ("theta", (6, 10))):
x, mv = np.log10(p[(f >= lo) & (f <= hi)].mean(axis=0)), r["mov"]
rho = {k: (stats.spearmanr(x[k:], mv[:len(mv) - k])[0] if k >= 0 else stats.spearmanr(x[:k], mv[-k:])[0])
for k in shifts}
rows.append({"measure": name, "rat": r["rat"], "condition": r["condition"], "session": r["session"],
**{f"s{k:+d}": v for k, v in rho.items()}, "best_shift": max(rho, key=rho.get)})
sync = pd.DataFrame(rows)
print("Spearman correlation of log band power with movement; s+k pairs the LFP of segment i+k with the movement of segment i:")
print(sync.round(2).to_string(index=False))
for name, g in sync.groupby("measure"):
print(f"\n{name}: best shift per file {g.best_shift.value_counts().sort_index().to_dict()}; median correlation "
f"at 0 {g['s+0'].median():.2f}, at +1 {g['s+1'].median():.2f}, at -1 {g['s-1'].median():.2f}, "
f"at +-5 {np.median(np.r_[g['s+5'], g['s-5']]):.2f}")
Spearman correlation of log band power with movement; s+k pairs the LFP of segment i+k with the movement of segment i:
measure rat condition session s-5 s-4 s-3 s-2 s-1 s+0 s+1 s+2 s+3 s+4 s+5 best_shift
100-200 Hz 1 Arena 1 0.31 0.34 0.35 0.39 0.42 0.60 0.57 0.42 0.38 0.36 0.32 0
theta 1 Arena 1 0.36 0.30 0.22 0.16 0.31 0.34 0.27 0.19 0.21 0.26 0.28 -5
100-200 Hz 1 Arena 2 0.30 0.31 0.35 0.33 0.33 0.44 0.48 0.37 0.28 0.30 0.39 1
theta 1 Arena 2 0.11 0.06 0.18 0.18 0.24 0.33 0.21 0.15 0.15 0.19 0.25 0
100-200 Hz 1 HCage 1 0.04 0.05 0.07 0.11 0.17 0.29 0.28 0.18 0.14 0.13 0.11 0
theta 1 HCage 1 0.32 0.32 0.32 0.31 0.32 0.37 0.34 0.26 0.24 0.23 0.20 0
100-200 Hz 1 HCage 2 0.15 0.18 0.25 0.32 0.47 0.68 0.66 0.46 0.32 0.22 0.18 0
theta 1 HCage 2 -0.05 -0.03 0.01 0.02 0.01 -0.01 -0.05 -0.06 -0.07 -0.06 -0.05 -2
100-200 Hz 2 Arena 1 0.33 0.33 0.30 0.29 0.38 0.63 0.66 0.42 0.31 0.28 0.27 1
theta 2 Arena 1 -0.11 -0.15 -0.13 -0.14 -0.12 -0.01 -0.06 -0.20 -0.19 -0.10 -0.07 0
100-200 Hz 2 Arena 2 0.21 0.26 0.24 0.29 0.33 0.50 0.53 0.31 0.19 0.14 0.11 1
theta 2 Arena 2 0.21 0.16 0.15 0.24 0.26 0.26 0.19 0.11 0.12 0.20 0.17 0
100-200 Hz 2 HCage 1 0.10 0.10 0.09 0.09 0.15 0.25 0.21 0.16 0.18 0.18 0.18 0
theta 2 HCage 1 0.12 0.16 0.18 0.17 0.23 0.32 0.25 0.11 0.10 0.10 0.11 0
100-200 Hz 2 HCage 2 -0.33 -0.35 -0.33 -0.28 -0.21 -0.09 -0.03 -0.05 -0.04 -0.05 -0.12 1
theta 2 HCage 2 0.14 0.14 0.14 0.18 0.23 0.25 0.20 0.19 0.23 0.24 0.22 0
100-200 Hz 3 Arena 1 0.09 0.12 0.11 0.07 0.13 0.28 0.28 0.14 0.08 -0.00 0.00 0
theta 3 Arena 1 -0.12 -0.13 -0.13 -0.16 -0.12 -0.12 -0.19 -0.16 -0.11 -0.09 -0.12 4
100-200 Hz 3 Arena 2 0.10 0.08 0.09 0.07 0.06 0.26 0.30 0.17 0.14 0.10 0.15 1
theta 3 Arena 2 -0.00 -0.07 -0.10 -0.13 -0.05 0.04 -0.06 -0.07 -0.11 -0.05 -0.07 0
100-200 Hz 3 HCage 1 -0.06 -0.04 -0.03 -0.01 0.02 0.11 0.13 0.01 -0.05 -0.09 -0.11 1
theta 3 HCage 1 -0.09 -0.10 -0.10 -0.06 -0.04 -0.06 -0.06 -0.04 -0.03 -0.04 -0.09 3
100-200 Hz 3 HCage 2 0.02 -0.00 -0.01 -0.00 -0.01 0.08 0.08 -0.01 -0.04 -0.08 -0.09 1
theta 3 HCage 2 0.03 0.03 -0.01 0.06 0.12 0.12 0.12 0.11 0.12 0.06 0.00 0
100-200 Hz 4 Arena 1 0.21 0.20 0.22 0.22 0.24 0.43 0.44 0.30 0.18 0.15 0.20 1
theta 4 Arena 1 -0.05 -0.09 -0.08 -0.07 -0.06 0.08 0.03 -0.10 -0.02 0.04 -0.04 0
100-200 Hz 4 Arena 2 0.05 0.08 0.13 0.19 0.19 0.30 0.46 0.29 0.21 0.20 0.17 1
theta 4 Arena 2 0.07 0.05 0.01 0.01 0.04 0.10 0.02 -0.02 0.10 0.17 0.14 4
100-200 Hz 4 HCage 1 0.10 0.08 0.08 0.16 0.27 0.45 0.53 0.35 0.25 0.19 0.13 1
theta 4 HCage 1 0.03 0.03 0.03 0.00 0.01 0.08 0.03 -0.00 -0.05 -0.04 0.01 0
100-200 Hz 4 HCage 2 0.23 0.26 0.27 0.26 0.31 0.44 0.46 0.36 0.26 0.20 0.24 1
theta 4 HCage 2 0.09 0.09 0.09 0.13 0.12 0.14 0.16 0.16 0.17 0.16 0.17 3
100-200 Hz 5 Arena 1 0.25 0.30 0.30 0.28 0.31 0.49 0.45 0.34 0.30 0.26 0.20 0
theta 5 Arena 1 0.10 0.13 0.16 0.16 0.17 0.26 0.31 0.14 0.04 0.04 0.08 1
100-200 Hz 5 Arena 2 0.04 0.10 0.06 0.06 0.05 0.22 0.32 0.11 -0.05 0.04 0.10 1
theta 5 Arena 2 -0.02 -0.03 -0.07 -0.06 -0.06 -0.03 -0.04 0.02 0.02 -0.00 -0.00 2
100-200 Hz 5 HCage 1 0.17 0.20 0.24 0.28 0.29 0.47 0.51 0.34 0.17 0.04 0.01 1
theta 5 HCage 1 0.10 0.05 0.00 -0.06 -0.06 -0.06 -0.03 -0.04 -0.02 -0.03 -0.08 -5
100-200 Hz 5 HCage 2 -0.17 -0.22 -0.28 -0.26 -0.26 -0.19 -0.20 -0.29 -0.37 -0.45 -0.43 -5
theta 5 HCage 2 -0.17 -0.18 -0.18 -0.09 -0.04 -0.10 -0.16 -0.19 -0.19 -0.23 -0.20 -1
100-200 Hz 6 Arena 1 0.07 0.10 0.10 0.17 0.35 0.62 0.57 0.32 0.14 0.12 0.14 0
theta 6 Arena 1 0.08 0.14 0.16 0.16 0.14 0.22 0.24 0.08 0.10 0.12 0.14 1
100-200 Hz 6 Arena 2 0.21 0.24 0.22 0.28 0.38 0.55 0.54 0.36 0.27 0.29 0.31 0
theta 6 Arena 2 0.15 0.16 0.18 0.16 0.19 0.29 0.21 0.04 0.06 0.12 0.13 0
100-200 Hz 6 HCage 1 0.33 0.30 0.32 0.37 0.47 0.57 0.56 0.46 0.39 0.38 0.38 0
theta 6 HCage 1 0.14 0.16 0.18 0.20 0.23 0.36 0.35 0.21 0.09 0.07 0.07 0
100-200 Hz 6 HCage 2 0.13 0.20 0.24 0.29 0.39 0.60 0.56 0.40 0.27 0.17 0.08 0
theta 6 HCage 2 -0.12 -0.12 -0.13 -0.12 -0.04 -0.03 -0.10 -0.10 -0.11 -0.10 -0.11 0
100-200 Hz: best shift per file {-5: 1, 0: 10, 1: 13}; median correlation at 0 0.44, at +1 0.46, at -1 0.28, at +-5 0.13
theta: best shift per file {-5: 2, -2: 1, -1: 1, 0: 13, 1: 2, 2: 1, 3: 2, 4: 2}; median correlation at 0 0.11, at +1 0.07, at -1 0.08, at +-5 0.07
# Which way does movement go with gamma and theta power, by condition (epochs, per file)?
tab_all = bflfp.epoch_table(recordings, 2, "asis")
rows = []
for (rat, cond, sess), g in tab_all.groupby(["rat", "condition", "session"]):
rows.append({"rat": rat, "condition": cond, "session": sess,
"rho_mov_gamma": stats.spearmanr(g.mov, g.gamma)[0],
"rho_mov_theta": stats.spearmanr(g.mov, g.theta)[0]})
print(pd.DataFrame(rows).round(2).to_string(index=False))
print("Median movement value by condition, per rat:")
print(tab_all.groupby(["rat", "condition"]).mov.median().unstack().round(0).to_string())
rat condition session rho_mov_gamma rho_mov_theta 1 Arena 1 0.01 0.34 1 Arena 2 0.13 0.33 1 HCage 1 0.42 0.37 1 HCage 2 0.49 -0.01 2 Arena 1 0.43 -0.01 2 Arena 2 0.06 0.26 2 HCage 1 0.45 0.32 2 HCage 2 0.50 0.25 3 Arena 1 0.27 -0.12 3 Arena 2 0.07 0.04 3 HCage 1 0.42 -0.06 3 HCage 2 0.45 0.12 4 Arena 1 -0.06 0.08 4 Arena 2 -0.11 0.10 4 HCage 1 0.30 0.08 4 HCage 2 0.47 0.14 5 Arena 1 -0.13 0.26 5 Arena 2 0.05 -0.03 5 HCage 1 0.35 -0.06 5 HCage 2 0.27 -0.10 6 Arena 1 0.01 0.22 6 Arena 2 0.21 0.29 6 HCage 1 0.56 0.36 6 HCage 2 0.35 -0.03 Median movement value by condition, per rat: condition Arena HCage rat 1 496.0 248.0 2 440.0 316.0 3 394.0 213.0 4 102.0 55.0 5 95.0 68.0 6 329.0 238.0
In 23 of 24 files the correlation of 100-200 Hz power with the movement values is largest at a shift of 0 or +1 segment and falls off at larger shifts (printed above), so the movement values are aligned with the LFP segments to within one segment (2 s). In 13 files the peak is at +1: the movement value of a segment goes best with the LFP of the next segment, so the movement stream may lead the LFP by up to one segment; this is consistent with the data and not shown to be an offset. Theta power gives a weaker and less regular relation. Movement is higher in the arena than in the home cage in every rat, as the conditions imply. This is consistent with the labels and the movement stream, and does not by itself prove each file's label. No two files share a clock, so alignment across files cannot be checked.
Units of analysis and readings shared by the claims¶
Gamma power and the other band powers are computed per 2-s segment (bflfp.epoch_table). The
claims file's header names readings that apply to several claims, and each is run:
- epoch length: 2-s segments (main) or 10-s epochs of five contiguous segments;
- artifacts: segments as delivered (main) or the SI's peak-to-peak rule applied again;
- unit of pairing: per rat, the two sessions of a condition averaged (rat 3's home cage is session 1 alone), n = 6; or session pairs (session 1 with session 1), n = 11 without the duplicate, and 12 with the files as delivered.
def drop_duplicates(tab):
keys = list(zip(tab.rat, tab.condition, tab.session))
return tab[[k not in DUPLICATES for k in keys]].copy()
def unit_means(tab, unit, measure="gamma"):
# Per recording, the mean over its epochs; then per rat (sessions averaged) or per session pair.
m = tab.groupby(["rat", "condition", "session"])[measure].mean().reset_index()
if unit != "session_asis":
m = m[[(a, b, c) not in DUPLICATES for a, b, c in zip(m.rat, m.condition, m.session)]]
if unit == "rat":
return m.groupby(["rat", "condition"])[measure].mean().unstack()
return m.set_index(["rat", "session", "condition"])[measure].unstack().dropna()
tables = {(ep, art): bflfp.epoch_table(recordings, ep, art) for ep in (2, 10) for art in ("asis", "ptp")}
for (ep, art), t in tables.items():
kept = t.groupby(["rat", "condition", "session"]).size()
print(f"{ep}-s epochs, artifacts {art}: {len(t)} epochs, per file {kept.min()} to {kept.max()}")
n2 = sum(r["lfp"].shape[1] for r in recordings)
print(f"The peak-to-peak rule drops {n2 - len(tables[(2, 'ptp')])} of {n2} segments "
f"({100 * (1 - len(tables[(2, 'ptp')]) / n2):.1f}%); the SI reports 1-4%")
main = drop_duplicates(tables[(2, "asis")])
2-s epochs, artifacts asis: 8556 epochs, per file 255 to 948 2-s epochs, artifacts ptp: 7118 epochs, per file 221 to 808 10-s epochs, artifacts asis: 1476 epochs, per file 41 to 180 10-s epochs, artifacts ptp: 1261 epochs, per file 34 to 154 The peak-to-peak rule drops 1438 of 8556 segments (16.8%); the SI reports 1-4%
fig3a-bf-gamma-reduction¶
BF gamma power was reduced by 35 +/- 3% when rats moved from the home cage to the novel arena, a highly significant reduction. (Expected 35 +/- 3%; criterion: reduction within tolerance under reading (a) mean over recordings of 100 x (1 - arena / home cage) or (b) from the group means, and a paired t test over rats at p < 0.05 with home cage higher; partially reproduced when the test holds and no reading is within tolerance.)
rows = []
for (ep, art), t in tables.items():
for unit in ("rat", "session", "session_asis"):
w = unit_means(t, unit)
h, a = w["HCage"], w["Arena"]
test = stats.ttest_rel(h, a)
rows.append({"epoch_s": ep, "artifacts": art, "unit": unit, "n": len(w),
"reading_a_pct": (100 * (1 - a / h)).mean(), "reading_b_pct": 100 * (1 - a.mean() / h.mean()),
"sem_a": stats.sem(100 * (1 - a / h)), "t": test.statistic, "p": test.pvalue,
"home_higher_in": f"{int((h > a).sum())}/{len(w)}"})
grid = pd.DataFrame(rows)
print(grid.round(4).to_string(index=False))
print("\nReadings with the reduction within 32-38%:",
int(((grid.reading_a_pct.between(32, 38)) | (grid.reading_b_pct.between(32, 38))).sum()), "of", len(grid))
w = unit_means(main, "rat")
print("\nPer rat (main reading: 2-s, as delivered), gamma power (uV^2/Hz if the LFP is in uV):")
print(w.assign(reduction_pct=100 * (1 - w.Arena / w.HCage)).round(1).to_string())
epoch_s artifacts unit n reading_a_pct reading_b_pct sem_a t p home_higher_in
2 asis rat 6 48.0473 48.0196 3.5832 3.8443 0.0121 6/6
2 asis session 11 47.7245 47.9953 3.4179 5.0958 0.0005 11/11
2 asis session_asis 12 47.7780 48.0134 3.1206 5.2753 0.0003 12/12
2 ptp rat 6 47.2548 47.5060 3.5077 3.6860 0.0142 6/6
2 ptp session 11 46.8431 47.4976 3.5070 4.8875 0.0006 11/11
2 ptp session_asis 12 46.9023 47.5004 3.2020 5.0768 0.0004 12/12
10 asis rat 6 48.1002 48.0612 3.6526 3.8799 0.0116 6/6
10 asis session 11 47.6699 48.0389 3.6382 5.1140 0.0005 11/11
10 asis session_asis 12 47.7852 48.0891 3.3232 5.3084 0.0002 12/12
10 ptp rat 6 48.0663 48.0006 3.4964 3.9241 0.0111 6/6
10 ptp session 11 47.6082 47.9685 3.5607 5.2056 0.0004 11/11
10 ptp session_asis 12 47.7115 48.0129 3.2521 5.4115 0.0002 12/12
Readings with the reduction within 32-38%: 0 of 12
Per rat (main reading: 2-s, as delivered), gamma power (uV^2/Hz if the LFP is in uV):
condition Arena HCage reduction_pct
rat
1 20.2 56.0 63.9
2 31.3 60.2 48.1
3 18.7 36.7 49.0
4 23.3 42.2 44.9
5 79.4 145.0 45.2
6 22.6 36.0 37.3
# Variants not named by the criterion: 49-51 Hz left out of the band (line noise), and the amplitude
# spectrum |FFT| in place of power, which the data descriptor's example uses (abs(fft(...))).
rows = []
for r in recordings:
if bflfp.key(r) in DUPLICATES:
continue
f, p = bflfp.spectrum(r["lfp"])
band = (f >= 30) & (f <= 80)
notch = band & ~((f >= 49) & (f <= 51))
rows.append({"rat": r["rat"], "condition": r["condition"], "session": r["session"],
"power": p[band].mean(), "power_no_50": p[notch].mean(), "amplitude": np.sqrt(p[band]).mean()})
var = pd.DataFrame(rows)
for measure in ("power", "power_no_50", "amplitude"):
w = var.groupby(["rat", "condition"])[measure].mean().unstack()
h, a = w["HCage"], w["Arena"]
print(f"{measure:12s} reduction (a) {(100 * (1 - a / h)).mean():.1f}% (b) {100 * (1 - a.mean() / h.mean()):.1f}% "
f"paired t p = {stats.ttest_rel(h, a).pvalue:.3g}, n = {len(w)}")
power reduction (a) 48.0% (b) 48.0% paired t p = 0.0121, n = 6 power_no_50 reduction (a) 44.7% (b) 44.5% paired t p = 0.0121, n = 6 amplitude reduction (a) 24.0% (b) 23.9% paired t p = 0.000619, n = 6
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5))
for r in recordings:
if bflfp.key(r) in DUPLICATES:
continue
f, p = bflfp.spectrum(r["lfp"])
keep = (f >= 1) & (f <= 100)
ax[0].semilogy(f[keep], p[keep].mean(axis=1), color="C0" if r["condition"] == "HCage" else "C3", lw=0.7)
ax[0].set(xlabel="frequency (Hz)", ylabel="power (LFP unit^2/Hz)", title="Mean spectrum per file (blue home cage, red arena)")
w = unit_means(main, "rat")
for rat, row in w.iterrows():
ax[1].plot([0, 1], [row.HCage, row.Arena], "o-", color="0.4")
ax[1].set(xticks=[0, 1], xticklabels=["home cage", "arena"], ylabel="gamma power", title="Per rat")
plt.tight_layout()
plt.show()
print("Plotted gamma power per rat ranges from %.1f to %.1f" % (w.values.min(), w.values.max()))
Plotted gamma power per rat ranges from 18.7 to 145.0
Result: partially reproduced. Home-cage gamma exceeds arena gamma in every rat and every session pair, and the paired t test gives p = 0.012 over 6 rats (p = 0.00047 over 11 session pairs). The reduction is 48% under both percentage readings and every epoch, artifact and pairing reading (46.8% to 48.1%), outside 35 +/- 3%. Leaving out 49-51 Hz does not change this; the amplitude spectrum, which the paper does not use, gives a smaller reduction (printed above).
fig3a-each-hemisphere-significant¶
The BF gamma reduction in the arena was significant in every one of the 22 hemispheres analyzed individually. (Criterion: per recording, a rat's arena session against its home-cage session of the same number, unpaired t test of epoch gamma, every recording p < 0.01 with home cage higher.)
rows = []
for ep_art, t in (("2-s, as delivered", tables[(2, "asis")]), ("2-s, peak-to-peak rule", tables[(2, "ptp")]),
("10-s, as delivered", tables[(10, "asis")])):
for (rat, sess), g in t[t.condition == "Arena"].groupby(["rat", "session"]):
for home_session, label in ((sess, "same number"), (1, "home session 1")):
if label == "home session 1" and not (rat == 3 and sess == 2):
continue
h = t[(t.rat == rat) & (t.condition == "HCage") & (t.session == home_session)]
test = stats.ttest_ind(h.gamma, g.gamma)
rows.append({"epochs": ep_art, "rat": rat, "session": sess, "home_cage_file": f"session {home_session}",
"n_home": len(h), "n_arena": len(g), "home_mean": h.gamma.mean(), "arena_mean": g.gamma.mean(),
"t": test.statistic, "p": test.pvalue})
each = pd.DataFrame(rows)
print(each.to_string(index=False, float_format=lambda x: f"{x:.3g}"))
for ep_art, g in each.groupby("epochs", sort=False):
ok = (g.p < 0.01) & (g.home_mean > g.arena_mean)
print(f"{ep_art}: {int(ok.sum())} of {len(g)} comparisons p < 0.01 with home cage higher; largest p = {g.p.max():.2g}")
epochs rat session home_cage_file n_home n_arena home_mean arena_mean t p
2-s, as delivered 1 1 session 1 948 283 57.7 24 17.8 3.51e-63
2-s, as delivered 1 2 session 2 354 281 54.3 16.5 23.8 1.08e-89
2-s, as delivered 2 1 session 1 623 302 68.4 34.5 27.5 2.43e-122
2-s, as delivered 2 2 session 2 317 255 52 28.1 16.4 7.13e-50
2-s, as delivered 3 1 session 1 449 332 36.7 18.6 29.6 2.08e-129
2-s, as delivered 3 2 session 2 291 301 36.6 18.9 26.9 9.89e-105
2-s, as delivered 3 2 session 1 449 301 36.7 18.9 28.2 6.68e-120
2-s, as delivered 4 1 session 1 359 354 45.1 26.8 19.8 7.34e-70
2-s, as delivered 4 2 session 2 289 335 39.4 19.7 21.1 6.67e-75
2-s, as delivered 5 1 session 1 290 312 135 84.8 20.4 2.1e-70
2-s, as delivered 5 2 session 2 286 323 155 74.1 28.5 2.91e-114
2-s, as delivered 6 1 session 1 370 291 31.7 23.5 12 2.75e-30
2-s, as delivered 6 2 session 2 261 350 40.3 21.6 20.6 6.84e-72
2-s, peak-to-peak rule 1 1 session 1 808 226 52.5 22.6 15.3 4.92e-48
2-s, peak-to-peak rule 1 2 session 2 303 230 49.2 16.1 21.2 1.71e-72
2-s, peak-to-peak rule 2 1 session 1 528 252 65 33 26.6 3.44e-111
2-s, peak-to-peak rule 2 2 session 2 260 221 46.5 26.7 15.4 8.52e-44
2-s, peak-to-peak rule 3 1 session 1 365 277 35.5 18.3 26.6 8.44e-106
2-s, peak-to-peak rule 3 2 session 2 230 234 35.5 18.6 23.3 7.07e-80
2-s, peak-to-peak rule 3 2 session 1 365 234 35.5 18.6 24.6 4.72e-93
2-s, peak-to-peak rule 4 1 session 1 307 285 43.2 25.7 18.2 6.67e-59
2-s, peak-to-peak rule 4 2 session 2 240 275 38.6 18.4 21.5 9.16e-74
2-s, peak-to-peak rule 5 1 session 1 251 248 129 80.2 19.6 6.62e-64
2-s, peak-to-peak rule 5 2 session 2 243 271 147 69.4 28.9 1.63e-109
2-s, peak-to-peak rule 6 1 session 1 299 244 29.5 23 9.65 1.92e-20
2-s, peak-to-peak rule 6 2 session 2 230 291 38.5 21.1 18.5 5.25e-59
10-s, as delivered 1 1 session 1 180 51 57.5 23.3 10.1 5.45e-20
10-s, as delivered 1 2 session 2 61 47 53.1 16.4 14 6.44e-26
10-s, as delivered 2 1 session 1 101 58 67.5 34 18.5 3.06e-41
10-s, as delivered 2 2 session 2 52 48 51 27.8 11.2 3.68e-19
10-s, as delivered 3 1 session 1 79 62 36.3 18.3 18.7 7.27e-40
10-s, as delivered 3 2 session 2 53 45 36.7 18.7 16.6 6.21e-30
10-s, as delivered 3 2 session 1 79 45 36.3 18.7 16.5 6.77e-33
10-s, as delivered 4 1 session 1 65 63 45 26.7 13.7 1.28e-26
10-s, as delivered 4 2 session 2 43 63 40.1 19.4 14 1.11e-25
10-s, as delivered 5 1 session 1 46 49 132 83 11.6 7.67e-20
10-s, as delivered 5 2 session 2 47 58 153 73 16.7 5.58e-31
10-s, as delivered 6 1 session 1 54 50 30.4 23.6 5.71 1.11e-07
10-s, as delivered 6 2 session 2 41 60 40.6 21.3 12.8 1.06e-22
2-s, as delivered: 13 of 13 comparisons p < 0.01 with home cage higher; largest p = 2.7e-30
2-s, peak-to-peak rule: 13 of 13 comparisons p < 0.01 with home cage higher; largest p = 1.9e-20
10-s, as delivered: 13 of 13 comparisons p < 0.01 with home cage higher; largest p = 1.1e-07
Result: reproduced. Each of the 12 arena recordings has lower gamma than its home-cage recording at p < 0.01 by the unpaired t test over epochs, under every epoch reading, including rat 3's arena session 2 against either of its home-cage files (largest p printed above; for that recording p < 1e-100 with 2-s epochs as delivered, 7.1e-80 and 4.7e-93 under the peak-to-peak rule, 6.2e-30 and 6.8e-33 with 10-s epochs). The unit is a rat and session, not a hemisphere.
fig3b-arena-gamma-unrelated-to-movement¶
In the arena, BF gamma power did not vary with movement-sensor quartile. (Criterion, absence: 2-s epochs in quartiles of the movement value within each condition, within each recording or over a rat's pooled sessions; one-way ANOVA over quartiles on per-rat quartile means; reproduced when p > 0.1, not reproduced when p < 0.05. Variants: the ANOVA over epochs, and a repeated-measures ANOVA.)
Result: reproduced. In the arena the per-rat quartile means of gamma differ by a few percent (Q4/Q1 near 1), and the one-way ANOVA over quartiles gives p > 0.1 with quartiles formed within each recording or within each rat, under every epoch and artifact reading. The epoch-level and repeated-measures variants also give p > 0.1 (epochs 0.60 and 0.19, repeated measures 0.38 and 0.28), as does the one-way ANOVA on quartile means divided by each rat's mean (0.12 and 0.33). The criterion's test cannot see a within-rat effect against the between-rat spread (it misses the home-cage rise below), so these variants are the informative ones here.
fig3b-home-cage-gamma-increases-with-movement¶
In the home cage, BF gamma power increased with movement-sensor quartile. (Criterion: per-rat quartile means as for the arena claim; reproduced when the one-way ANOVA over quartiles gives p < 0.05 and the fourth quartile mean exceeds the first.)
print("Home cage, 2-s epochs as delivered:")
print(quartile_tests(main, "HCage").to_string(index=False, float_format=lambda x: f"{x:.3g}"))
print("\nHome cage, other epoch and artifact readings:")
for (ep, art), t in tables.items():
if (ep, art) == (2, "asis"):
continue
r = quartile_tests(drop_duplicates(t), "HCage")
print(f" {ep}-s, {art}: one-way p = {r.oneway_p.round(3).tolist()}, rm p = {r.rm_p.round(6).tolist()}, "
f"Q4/Q1 = {r.Q4_over_Q1.round(2).tolist()}")
print("\nPer-rat quartile means, home cage (quartiles within each recording):")
print(qm[qm.condition == "HCage"].pivot(index="rat", columns="q", values="gamma").round(1).to_string())
print("\nBetween-rat spread of home-cage gamma (per-rat means): %.1f to %.1f"
% (unit_means(main, "rat").HCage.min(), unit_means(main, "rat").HCage.max()))
Home cage, 2-s epochs as delivered:
quartiles_within Q1 Q2 Q3 Q4 Q4_over_Q1 rats_Q4_gt_Q1 oneway_F oneway_p eta2 oneway_normalized_p epochs_F epochs_p rm_F rm_p
recording 50.3 62.1 68.4 72 1.43 6 0.311 0.817 0.0446 1.56e-09 73 3.89e-46 19.4 2e-05
rat 49.5 60.9 70.1 72.3 1.46 6 0.371 0.774 0.0528 2.48e-10 81.7 1.71e-51 21.1 1.23e-05
Home cage, other epoch and artifact readings:
2-s, ptp: one-way p = [0.836, 0.803], rm p = [1.1e-05, 5e-06], Q4/Q1 = [1.42, 1.43]
10-s, asis: one-way p = [0.892, 0.839], rm p = [0.00089, 8e-06], Q4/Q1 = [1.33, 1.36]
10-s, ptp: one-way p = [0.885, 0.835], rm p = [5.2e-05, 2.1e-05], Q4/Q1 = [1.34, 1.35] Per-rat quartile means, home cage (quartiles within each recording): q 1 2 3 4 rat 1 40.7 50.9 61.8 73.6 2 47.0 62.5 69.7 72.4 3 30.9 36.4 38.7 41.0 4 35.3 40.3 44.3 50.5 5 122.2 146.4 157.0 154.2 6 25.9 35.9 38.7 40.6 Between-rat spread of home-cage gamma (per-rat means): 36.0 to 145.0
# Variant: the data check found the movement value of a segment may go with the LFP of the next
# segment. Pair each LFP segment with the movement value of the segment before it in its file.
shifted = main.copy()
shifted["mov"] = shifted.groupby(["rat", "condition", "session"]).mov.shift(1)
shifted = shifted.dropna(subset=["mov"])
for cond in ("Arena", "HCage"):
r = quartile_tests(shifted, cond)
print(f"{cond}, movement shifted by one segment: one-way p = {r.oneway_p.round(3).tolist()}, "
f"rm p = {r.rm_p.round(6).tolist()}, Q4/Q1 = {r.Q4_over_Q1.round(2).tolist()}")
Arena, movement shifted by one segment: one-way p = [0.999, 0.998], rm p = [0.136259, 0.227404], Q4/Q1 = [1.05, 0.97] HCage, movement shifted by one segment: one-way p = [0.79, 0.747], rm p = [2.2e-05, 1.4e-05], Q4/Q1 = [1.46, 1.49]
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5), sharey=True)
for i, cond in enumerate(("HCage", "Arena")):
d = qm[qm.condition == cond]
for rat, g in d.groupby("rat"):
ax[i].plot(g.q, g.gamma, "o-", color="0.5", lw=0.8)
ax[i].plot(range(1, 5), d.groupby("q").gamma.mean(), "k-", lw=2)
ax[i].set(title=cond, xlabel="movement quartile", xticks=range(1, 5))
ax[0].set_ylabel("gamma power")
plt.tight_layout()
plt.show()
print("Plotted values range from %.1f to %.1f" % (qm.gamma.min(), qm.gamma.max()))
Plotted values range from 17.9 to 157.0
Result: not reproduced. Home-cage gamma rises from the first to the fourth movement quartile in
every rat (Q4/Q1 1.43 and 1.46), but the criterion's test, a one-way ANOVA over quartiles on the
per-rat quartile means, gives p = 0.82 and 0.77, because per-rat home-cage gamma ranges from 36 to
145 (eta squared 0.045). The variants that take the rat into account find the rise: the
repeated-measures ANOVA gives p = 2.0e-05, the ANOVA over epochs p = 3.9e-46, and the one-way ANOVA
on quartile means divided by each rat's mean p = 1.6e-09. Pairing each segment with the movement
value of the segment before it, a variant, gives one-way p of 0.75 to 0.79. Whether the claims
file's per-rat one-way ANOVA is the paper's test is asked of a reader (questions/crcns-bf-1.yaml).
fig3b-matched-movement¶
Gamma suppression in the arena remained highly significant when epochs matched for movement were compared, for example the fourth home-cage quartile against the third arena quartile. (Criterion: reading (a) that pair; (b) the pair of quartiles, one per condition, whose mean movement values are closest; one-way ANOVA of the per-rat gamma means of the two groups, p < 0.05 with the home cage higher.)
rows = []
for scope in ("recording", "rat"):
qm_s, _ = quartile_means(main, scope)
mov = qm_s.groupby(["condition", "q"]).mov.mean().unstack()
diff = pd.DataFrame([(h, a, abs(mov.loc["HCage", h] - mov.loc["Arena", a])) for h in range(1, 5) for a in range(1, 5)],
columns=["home_q", "arena_q", "mov_diff"])
closest = diff.sort_values("mov_diff").iloc[0]
for reading, (hq, aq) in (("(a) home Q4 vs arena Q3", (4, 3)),
("(b) closest movement", (int(closest.home_q), int(closest.arena_q)))):
h = qm_s[(qm_s.condition == "HCage") & (qm_s.q == hq)].set_index("rat")
a = qm_s[(qm_s.condition == "Arena") & (qm_s.q == aq)].set_index("rat")
oneway = stats.f_oneway(h.gamma, a.gamma)
paired = stats.ttest_rel(h.gamma, a.gamma)
# variant: each rat's two values divided by their mean
scale = (h.gamma + a.gamma) / 2
oneway_norm = stats.f_oneway(h.gamma / scale, a.gamma / scale)
rows.append({"quartiles_within": scope, "reading": reading, "home_q": hq, "arena_q": aq,
"home_mov": h.mov.mean(), "arena_mov": a.mov.mean(), "home_gamma": h.gamma.mean(),
"arena_gamma": a.gamma.mean(), "rats_home_higher": int((h.gamma > a.gamma).sum()),
"oneway_p": oneway.pvalue, "paired_t_p_variant": paired.pvalue,
"oneway_normalized_p_variant": oneway_norm.pvalue})
if scope == "recording":
print("Mean movement value per quartile (mean over rats), quartiles within each recording:")
print(mov.round(0).to_string())
matched = pd.DataFrame(rows)
print(matched.to_string(index=False, float_format=lambda x: f"{x:.3g}"))
# Variant: movement matched within each rat, the home-cage and arena quartiles closest in that rat.
rows = []
qm_r, _ = quartile_means(main, "recording")
for rat, d in qm_r.groupby("rat"):
h, a = d[d.condition == "HCage"].set_index("q"), d[d.condition == "Arena"].set_index("q")
hq, aq = min(((x, y) for x in range(1, 5) for y in range(1, 5)), key=lambda p: abs(h.mov[p[0]] - a.mov[p[1]]))
rows.append({"rat": rat, "home_q": hq, "arena_q": aq, "home_mov": h.mov[hq], "arena_mov": a.mov[aq],
"home_gamma": h.gamma[hq], "arena_gamma": a.gamma[aq]})
within = pd.DataFrame(rows)
print("\nVariant, quartiles matched within each rat:")
print(within.round(1).to_string(index=False))
print("one-way p = %.3g, paired t p = %.3g" % (stats.f_oneway(within.home_gamma, within.arena_gamma).pvalue,
stats.ttest_rel(within.home_gamma, within.arena_gamma).pvalue))
Mean movement value per quartile (mean over rats), quartiles within each recording:
q 1 2 3 4
condition
Arena 202.0 287.0 333.0 421.0
HCage 77.0 161.0 228.0 356.0
quartiles_within reading home_q arena_q home_mov arena_mov home_gamma arena_gamma rats_home_higher oneway_p paired_t_p_variant oneway_normalized_p_variant
recording (a) home Q4 vs arena Q3 4 3 356 333 72 32 6 0.07 0.00714 7.45e-08
recording (b) closest movement 4 3 356 333 72 32 6 0.07 0.00714 7.45e-08
rat (a) home Q4 vs arena Q3 4 3 360 337 72.3 31.6 6 0.0703 0.00692 5.7e-08
rat (b) closest movement 4 3 360 337 72.3 31.6 6 0.0703 0.00692 5.7e-08
Variant, quartiles matched within each rat:
rat home_q arena_q home_mov arena_mov home_gamma arena_gamma
1 3 1 297.6 324.8 61.8 19.9
2 4 3 478.7 473.9 72.4 30.6
3 3 1 279.8 268.5 38.7 18.3
4 4 3 107.5 107.4 50.5 22.8
5 4 4 124.4 123.4 154.2 78.5
6 2 1 196.8 191.4 35.9 21.5
one-way p = 0.0985, paired t p = 0.00914
Result: not reproduced. Both readings pick the home-cage fourth and arena third quartiles, and the movement-matched home-cage gamma (72.0) exceeds the arena gamma (32.0) in 6 of 6 rats, but the criterion's one-way ANOVA of the per-rat means of the two groups gives p = 0.070, above 0.05. A paired t test of the same values, a variant, gives p = 0.0071, and the one-way ANOVA on values divided by each rat's mean gives p = 7.5e-08. As for the home-cage rise, the choice of test is asked of a reader.
fig3c-home-cage-above-both-sessions¶
BF gamma was larger in the home cage than in either of the two daily arena sessions, and slightly larger (5.62 +/- 1.2%) in the second arena session than in the first, with no interaction between day and session. (Criterion: levels home cage, arena 1, arena 2; Tukey's HSD, with Bonferroni-corrected paired t tests as the variant; increase as 100 x (arena 2 / arena 1 - 1) averaged over recordings, or from the group means. Without day labels the two arena sessions are compared within rat by a paired t test and the day factor is dropped.)
# The descriptor calls sessions 1 and 2 "two separate examples" of each session type; nothing in the
# file names, variables or documents gives a day or the order of trials within a day. The arena 2
# against arena 1 comparison is printed but cannot test the claim's second part.
for label, t in (("2-s, as delivered", main), ("2-s, peak-to-peak rule", drop_duplicates(tables[(2, "ptp")])),
("10-s, as delivered", drop_duplicates(tables[(10, "asis")]))):
m = t.groupby(["rat", "condition", "session"]).gamma.mean()
home = m.xs("HCage", level="condition").groupby("rat").mean()
a1 = m.xs("Arena", level="condition").xs(1, level="session")
a2 = m.xs("Arena", level="condition").xs(2, level="session")
levels = pd.DataFrame({"home": home, "arena1": a1, "arena2": a2})
print(f"--- {label}")
if label.startswith("2-s, as"):
print(levels.round(1).to_string())
tukey = pairwise_tukeyhsd(np.r_[home, a1, a2], ["home"] * 6 + ["arena1"] * 6 + ["arena2"] * 6)
print(tukey.summary())
for x, y, name in ((home, a1, "home vs arena1"), (home, a2, "home vs arena2"), (a2, a1, "arena2 vs arena1")):
p = stats.ttest_rel(x, y).pvalue
print(f"paired t {name}: mean difference {np.mean(x - y):.2f}, p = {p:.3g}, Bonferroni (x3) p = {min(1, 3 * p):.3g}")
inc = 100 * (a2 / a1 - 1)
print("increase arena2 over arena1: (a) mean over rats %.2f%% (per rat %s); (b) group means %.2f%%"
% (inc.mean(), inc.round(1).tolist(), 100 * (a2.mean() / a1.mean() - 1)))
--- 2-s, as delivered
home arena1 arena2
rat
1 56.0 24.0 16.5
2 60.2 34.5 28.1
3 36.7 18.6 18.9
4 42.2 26.8 19.7
5 145.0 84.8 74.1
6 36.0 23.5 21.6
Multiple Comparison of Means - Tukey HSD, FWER=0.05
=====================================================
group1 group2 meandiff p-adj lower upper reject
-----------------------------------------------------
arena1 arena2 -5.5423 0.9476 -51.5418 40.4572 False
arena1 home 27.3301 0.2998 -18.6694 73.3296 False
arena2 home 32.8724 0.1857 -13.1271 78.8719 False
-----------------------------------------------------
paired t home vs arena1: mean difference 27.33, p = 0.0126, Bonferroni (x3) p = 0.0379
paired t home vs arena2: mean difference 32.87, p = 0.0118, Bonferroni (x3) p = 0.0354
paired t arena2 vs arena1: mean difference -5.54, p = 0.0201, Bonferroni (x3) p = 0.0604
increase arena2 over arena1: (a) mean over rats -15.86% (per rat [-31.3, -18.5, 1.9, -26.6, -12.6, -8.1]); (b) group means -15.68%
--- 2-s, peak-to-peak rule
Multiple Comparison of Means - Tukey HSD, FWER=0.05 ===================================================== group1 group2 meandiff p-adj lower upper reject ----------------------------------------------------- arena1 arena2 -5.4244 0.9443 -49.0225 38.1737 False arena1 home 25.4338 0.312 -18.1643 69.0319 False arena2 home 30.8582 0.1911 -12.7399 74.4562 False ----------------------------------------------------- paired t home vs arena1: mean difference 25.43, p = 0.0151, Bonferroni (x3) p = 0.0452 paired t home vs arena2: mean difference 30.86, p = 0.0137, Bonferroni (x3) p = 0.0412 paired t arena2 vs arena1: mean difference -5.42, p = 0.0211, Bonferroni (x3) p = 0.0633 increase arena2 over arena1: (a) mean over rats -16.05% (per rat [-28.8, -19.2, 1.8, -28.4, -13.5, -8.2]); (b) group means -16.04%
--- 10-s, as delivered
Multiple Comparison of Means - Tukey HSD, FWER=0.05 ===================================================== group1 group2 meandiff p-adj lower upper reject ----------------------------------------------------- arena1 arena2 -5.3891 0.9484 -50.4534 39.6752 False arena1 home 27.0522 0.2929 -18.0121 72.1166 False arena2 home 32.4414 0.1816 -12.6229 77.5057 False ----------------------------------------------------- paired t home vs arena1: mean difference 27.05, p = 0.0125, Bonferroni (x3) p = 0.0375 paired t home vs arena2: mean difference 32.44, p = 0.0111, Bonferroni (x3) p = 0.0334 paired t arena2 vs arena1: mean difference -5.39, p = 0.0169, Bonferroni (x3) p = 0.0507 increase arena2 over arena1: (a) mean over rats -15.84% (per rat [-29.8, -18.1, 2.0, -27.1, -12.0, -9.8]); (b) group means -15.47%
Result: partially reproduced. The home cage exceeds both arena sessions by Bonferroni-corrected paired t tests (p = 0.038 and 0.035), but not by Tukey's HSD on the per-rat values (p = 0.30 and 0.19), which treats the 18 values as independent groups and ignores the pairing within rat. The criterion names Tukey as the standard reading and the paired tests as a variant; the outcome rests on the variant, chosen because the rats are paired across conditions. The second part of the claim, a rise from the first to the second daily arena session, cannot be tested: the descriptor calls sessions 1 and 2 "two separate examples" of each session type, and nothing in the data set gives a day or the order of trials. Arena session 2 has 15.9% lower gamma than session 1 (paired t p = 0.020), against the paper's +5.62%; this is reported, not judged, as session 2 need not be the second trial of a day.
fig3d-gamma-rises-with-time-in-arena¶
Within arena sessions BF gamma stayed suppressed but rose with time spent in the arena (r = 0.72, n = 22). (Expected r = 0.72 +/- 0.072; criterion: reading (a) Pearson r between time and gamma averaged across recordings in 1-min bins over the 10-min session, with 2-s epochs as the variant; (b) the mean over recordings of each recording's r; (c) r over all epochs pooled; reproduced when under one reading r is within tolerance and positive at p < 0.05.)
arena = main[main.condition == "Arena"].copy()
arena["minute"] = (arena.t // 60).astype(int)
rows = []
for window, a in (("first 10 min", arena[arena.t < 600]), ("whole file", arena)):
per = a.groupby(["rat", "session", "minute"]).gamma.mean().reset_index()
curve = per.groupby("minute").gamma.mean()
r_a, p_a = stats.pearsonr(curve.index, curve.values)
per_epoch = a.groupby("t").gamma.mean()
r_a2, p_a2 = stats.pearsonr(per_epoch.index, per_epoch.values)
rs = np.array([stats.pearsonr(g.t, g.gamma)[0] for _, g in a.groupby(["rat", "session"])])
t_b = stats.ttest_1samp(rs, 0)
r_c, p_c = stats.pearsonr(a.t, a.gamma)
rows += [{"window": window, "reading": "(a) 1-min bins, mean over recordings", "r": r_a, "p": p_a, "n": len(curve)},
{"window": window, "reading": "(a) variant, 2-s epochs, mean over recordings", "r": r_a2, "p": p_a2, "n": len(per_epoch)},
{"window": window, "reading": "(b) mean of per-recording r", "r": rs.mean(), "p": t_b.pvalue, "n": len(rs)},
{"window": window, "reading": "(c) all epochs pooled", "r": r_c, "p": p_c, "n": len(a)}]
if window == "first 10 min":
print("Gamma per minute, mean over the 12 arena recordings:", curve.round(1).tolist())
print("Per-recording r:", np.round(rs, 2).tolist())
print("Recordings contributing to each minute:", per.groupby("minute").size().tolist())
else:
print("Whole file, recordings contributing to each minute:", per.groupby("minute").size().tolist())
corr = pd.DataFrame(rows)
print(corr.to_string(index=False, float_format=lambda x: f"{x:.3g}"))
print("Readings within 0.648-0.792 and positive at p < 0.05:",
corr[(corr.r.between(0.648, 0.792)) & (corr.p < 0.05)][["window", "reading"]].values.tolist())
Gamma per minute, mean over the 12 arena recordings: [29.3, 29.1, 32.2, 32.2, 31.6, 33.1, 33.9, 33.1, 32.7, 34.4]
Per-recording r: [-0.02, -0.27, 0.06, 0.03, 0.28, 0.36, 0.36, 0.16, 0.07, 0.4, 0.24, -0.13]
Recordings contributing to each minute: [12, 12, 12, 12, 12, 12, 12, 12, 12, 11]
Whole file, recordings contributing to each minute: [12, 12, 12, 12, 12, 12, 12, 12, 12, 11, 11, 7, 2]
window reading r p n
first 10 min (a) 1-min bins, mean over recordings 0.859 0.00147 10
first 10 min (a) variant, 2-s epochs, mean over recordings 0.197 0.00061 300
first 10 min (b) mean of per-recording r 0.129 0.0562 12
first 10 min (c) all epochs pooled 0.0451 0.00944 3312
whole file (a) 1-min bins, mean over recordings 0.285 0.345 13
whole file (a) variant, 2-s epochs, mean over recordings 0.1 0.0481 388
whole file (b) mean of per-recording r 0.165 0.011 12
whole file (c) all epochs pooled 0.0742 5.99e-06 3719
Readings within 0.648-0.792 and positive at p < 0.05: []
# Where is arena entry? Each file's clock starts at 0 and files span 526 to 780 s, longer than the
# SI's 10-min trial. Fig. 3D's time axis runs from -10 to 10 min with the transfer to the home cage at
# 0, so the 10-min window is also anchored at the file's end. Reading (b) is the same under any shift
# of a file's clock; (a) and (c) are not. Bin width is shown, not judged (the criterion names 1 min).
arena["from_end"] = arena.t - (arena.groupby(["rat", "session"]).t.transform("max") + 2.0)
windows = {"first 10 min (t = 0 is entry)": (arena[arena.t < 600], "t"),
"last 10 min (file end is transfer)": (arena[arena.from_end >= -600], "from_end")}
rows = []
for window, (a, col) in windows.items():
for width in (20, 30, 60, 120):
b = a.assign(bin=(a[col] // width).astype(int))
curve = b.groupby(["rat", "session", "bin"]).gamma.mean().groupby("bin").mean()
r, p = stats.pearsonr(curve.index, curve.values)
rows.append({"window": window, "reading": f"(a) {width}-s bins", "r": r, "p": p, "n": len(curve)})
rs = np.array([stats.pearsonr(g[col], g.gamma)[0] for _, g in a.groupby(["rat", "session"])])
rows.append({"window": window, "reading": "(b) mean of per-recording r", "r": rs.mean(),
"p": stats.ttest_1samp(rs, 0).pvalue, "n": len(rs)})
r, p = stats.pearsonr(a[col], a.gamma)
rows.append({"window": window, "reading": "(c) all epochs pooled", "r": r, "p": p, "n": len(a)})
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda x: f"{x:.3g}"))
# Movement and gamma in the first minute, the middle and the last minute of each arena file: a
# transfer or handling at a file's edge would show as a change of the movement value there.
rows = []
for (rat, s), g in arena.groupby(["rat", "session"]):
first, last = g[g.t < 60], g[g.from_end >= -60]
mid = g[(g.t >= 60) & (g.from_end < -60)]
rows.append({"rat": rat, "session": s, "span_s": -g.from_end.min() + g.t.min(),
"mov_first_min": first.mov.mean(), "mov_middle": mid.mov.mean(), "mov_last_min": last.mov.mean(),
"gamma_first_min": first.gamma.mean(), "gamma_middle": mid.gamma.mean(), "gamma_last_min": last.gamma.mean()})
edges = pd.DataFrame(rows)
print(edges.round(1).to_string(index=False))
print("Files whose first-minute movement exceeds the middle: %d of %d; last minute below the middle: %d of %d"
% ((edges.mov_first_min > edges.mov_middle).sum(), len(edges), (edges.mov_last_min < edges.mov_middle).sum(), len(edges)))
window reading r p n
first 10 min (t = 0 is entry) (a) 20-s bins 0.645 0.000119 30
first 10 min (t = 0 is entry) (a) 30-s bins 0.744 0.000168 20
first 10 min (t = 0 is entry) (a) 60-s bins 0.859 0.00147 10
first 10 min (t = 0 is entry) (a) 120-s bins 0.86 0.0618 5
first 10 min (t = 0 is entry) (b) mean of per-recording r 0.129 0.0562 12
first 10 min (t = 0 is entry) (c) all epochs pooled 0.0451 0.00944 3312
last 10 min (file end is transfer) (a) 20-s bins 0.662 6.86e-05 30
last 10 min (file end is transfer) (a) 30-s bins 0.747 0.000152 20
last 10 min (file end is transfer) (a) 60-s bins 0.885 0.000664 10
last 10 min (file end is transfer) (a) 120-s bins 0.969 0.00646 5
last 10 min (file end is transfer) (b) mean of per-recording r 0.13 0.0214 12
last 10 min (file end is transfer) (c) all epochs pooled 0.0429 0.0137 3296
rat session span_s mov_first_min mov_middle mov_last_min gamma_first_min gamma_middle gamma_last_min
1 1 620.0 511.2 447.8 137.4 22.6 24.1 24.3
1 2 620.0 587.3 557.7 494.6 16.2 16.7 15.0
2 1 620.0 419.1 398.7 652.4 31.3 33.5 45.3
2 2 526.0 482.9 503.6 428.7 28.1 28.6 24.0
3 1 688.0 403.9 389.2 372.4 17.4 18.6 19.1
3 2 698.0 402.9 420.4 354.3 17.2 18.8 21.4
4 1 758.0 103.2 95.2 80.7 23.2 27.3 25.3
4 2 694.0 114.7 113.6 84.5 19.0 19.7 20.6
5 1 704.0 106.7 75.9 78.5 76.5 84.8 94.1
5 2 694.0 92.8 106.6 81.0 61.5 74.7 81.3
6 1 660.0 361.0 313.2 257.5 17.4 24.2 23.6
6 2 780.0 390.0 337.6 275.7 21.1 21.6 23.0
Files whose first-minute movement exceeds the middle: 9 of 12; last minute below the middle: 10 of 12
fig, ax = plt.subplots(figsize=(5, 3.5))
per = arena[arena.t < 600].groupby(["rat", "session", "minute"]).gamma.mean().reset_index()
for _, g in per.groupby(["rat", "session"]):
ax.plot(g.minute + 0.5, g.gamma / g.gamma.mean(), color="0.7", lw=0.8)
ax.plot(np.arange(10) + 0.5, per.groupby("minute").gamma.mean() / per.gamma.mean(), "k-", lw=2)
ax.set(xlabel="time in arena (min)", ylabel="gamma / recording mean", title="Arena gamma by minute")
plt.show()
print("Plotted ratios range from %.2f to %.2f" % ((per.gamma / per.groupby(["rat", "session"]).gamma.transform("mean")).min(),
(per.gamma / per.groupby(["rat", "session"]).gamma.transform("mean")).max()))
Plotted ratios range from 0.72 to 1.45
Result: partially reproduced. Gamma rises with time in the arena. Under reading (a), 1-min bins over the first 10 min, r = 0.859 (p = 0.0015), above the tolerance of 0.72 +/- 0.072; reading (b) gives a mean per-recording r of 0.129 (p = 0.056) and reading (c) r = 0.045 (p = 0.0094). No reading the criterion names puts r within tolerance, while the positive correlation holds at p < 0.05 under (a) and (c).
Whether a file's t = 0 is arena entry is not documented. With the 10-min window anchored at the file's end, which Fig. 3D's axis makes the transfer, 1-min bins give r = 0.885 (p = 0.00066) and reading (b) r = 0.130 (p = 0.021). Reading (b) does not depend on where a file's clock starts, and over whole files it gives r = 0.165 (p = 0.011). In 9 of 12 arena files the first minute has more movement than the middle, consistent with entry near t = 0 but not proof of it. Several of these differences are small. The outcome rests on the criterion's 1-min bins. Other bin widths would put r within tolerance: 30-s bins give r = 0.744 and 0.747 under the two anchors, and 20-s bins with the end anchor give 0.662 (p = 6.9e-5); 20-s bins from the start give 0.645, just below, and 2-min bins 0.860 and 0.969 (shown, not judged).
figs2-bf-low-frequency-bands¶
In the BF, delta and beta power, like gamma, were higher in the home cage than in the arena, while theta power was strongly suppressed in the home cage. (Criterion: band power per rat and condition, paired t test per band; reproduced when all four bands differ at p < 0.05 in the stated directions; partially when theta and gamma do and one of delta and beta does not.)
expected_sign = {"delta": 1, "theta": -1, "beta": 1, "gamma": 1}
rows = []
for (ep, art), t in tables.items():
for unit in ("rat", "session", "session_asis"):
for band, sign in expected_sign.items():
w = unit_means(t, unit, band)
test = stats.ttest_rel(w.HCage, w.Arena)
rows.append({"epoch_s": ep, "artifacts": art, "unit": unit, "band": band, "n": len(w),
"home": w.HCage.mean(), "arena": w.Arena.mean(), "home_over_arena": (w.HCage / w.Arena).mean(),
"t": test.statistic, "p": test.pvalue,
"holds": bool(test.pvalue < 0.05 and np.sign(test.statistic) == sign)})
bands = pd.DataFrame(rows)
print(bands[(bands.epoch_s == 2) & (bands.artifacts == "asis")].to_string(index=False, float_format=lambda x: f"{x:.3g}"))
summary = bands.groupby(["epoch_s", "artifacts", "unit"]).holds.agg(lambda x: "".join("Y" if v else "n" for v in x))
print("\nBands holding (delta, theta, beta, gamma) under every reading:")
print(summary.to_string())
epoch_s artifacts unit band n home arena home_over_arena t p holds
2 asis rat delta 6 334 223 1.51 3.01 0.0299 True
2 asis rat theta 6 161 252 0.639 -8.17 0.000448 True
2 asis rat beta 6 46.7 35 1.33 3.2 0.0239 True
2 asis rat gamma 6 62.7 32.6 1.98 3.84 0.0121 True
2 asis session delta 11 314 202 1.55 3.89 0.00301 True
2 asis session theta 11 157 247 0.641 -9.31 3.04e-06 True
2 asis session beta 11 47.1 35.2 1.34 4.04 0.00237 True
2 asis session gamma 11 65 33.8 2.01 5.1 0.000467 True
2 asis session_asis delta 12 333 223 1.52 4.17 0.00157 True
2 asis session_asis theta 12 162 252 0.645 -10.2 5.85e-07 True
2 asis session_asis beta 12 46.6 35 1.33 4.32 0.00122 True
2 asis session_asis gamma 12 62.7 32.6 2 5.28 0.000262 True
Bands holding (delta, theta, beta, gamma) under every reading:
epoch_s artifacts unit
2 asis rat YYYY
session YYYY
session_asis YYYY
ptp rat YYYY
session YYYY
session_asis YYYY
10 asis rat YYYY
session YYYY
session_asis YYYY
ptp rat YYYY
session YYYY
session_asis YYYY
Result: reproduced. Over the six rats, delta, beta and gamma power are higher in the home cage and theta lower, each at p < 0.05 by the paired t test (theta p < 0.001), and all four hold under every epoch, artifact and pairing reading (printed above).
figs3-unimodal-states¶
Within each condition, the distributions of BF theta and gamma power over 2-s epochs were unimodal, with no evidence of substates, in each of six rats. (Criterion, absence: Hartigan dip test per rat, condition and band; reproduced when every test gives p > 0.1, not reproduced when any gives p < 0.05; readings (a) a rat's two sessions of a condition pooled, (b) each session alone.)
rows = []
for reading in ("(a) sessions pooled", "(b) each session"):
groups = main.groupby(["rat", "condition"]) if reading.startswith("(a)") else main.groupby(["rat", "condition", "session"])
for k, g in groups:
for band in ("theta", "gamma"):
for scale, x in (("raw", g[band].values), ("log10", np.log10(g[band].values))):
dip, p = diptest.diptest(x)
rows.append({"reading": reading, "scale": scale, "group": k, "band": band, "n_epochs": len(g),
"minutes": len(g) * 2 / 60, "dip": dip, "p": p})
dips = pd.DataFrame(rows)
print(dips[dips.scale == "raw"].to_string(index=False, float_format=lambda x: f"{x:.3g}"))
for (reading, scale), g in dips.groupby(["reading", "scale"]):
print(f"{reading}, {scale} power: {len(g)} tests, smallest p = {g.p.min():.3f}, tests with p <= 0.1: {int((g.p <= 0.1).sum())}, "
f"with p < 0.05: {int((g.p < 0.05).sum())}")
# Variant: the duplicate file kept (rat 3 home cage session 2).
dup = tables[(2, "asis")]
for k in DUPLICATES:
g = dup[(dup.rat == k[0]) & (dup.condition == k[1]) & (dup.session == k[2])]
print("Variant, the left-out file", k, {b: round(diptest.diptest(g[b].values)[1], 3) for b in ("theta", "gamma")})
reading scale group band n_epochs minutes dip p
(a) sessions pooled raw (1, Arena) theta 564 18.8 0.011 0.959
(a) sessions pooled raw (1, Arena) gamma 564 18.8 0.0103 0.984
(a) sessions pooled raw (1, HCage) theta 1302 43.4 0.0065 0.991
(a) sessions pooled raw (1, HCage) gamma 1302 43.4 0.00714 0.971
(a) sessions pooled raw (2, Arena) theta 557 18.6 0.0152 0.554
(a) sessions pooled raw (2, Arena) gamma 557 18.6 0.00995 0.991
(a) sessions pooled raw (2, HCage) theta 940 31.3 0.00801 0.985
(a) sessions pooled raw (2, HCage) gamma 940 31.3 0.00736 0.992
(a) sessions pooled raw (3, Arena) theta 633 21.1 0.0118 0.854
(a) sessions pooled raw (3, Arena) gamma 633 21.1 0.0117 0.866
(a) sessions pooled raw (3, HCage) theta 449 15 0.0116 0.983
(a) sessions pooled raw (3, HCage) gamma 449 15 0.0117 0.98
(a) sessions pooled raw (4, Arena) theta 689 23 0.00743 0.996
(a) sessions pooled raw (4, Arena) gamma 689 23 0.00898 0.99
(a) sessions pooled raw (4, HCage) theta 648 21.6 0.00896 0.992
(a) sessions pooled raw (4, HCage) gamma 648 21.6 0.0126 0.745
(a) sessions pooled raw (5, Arena) theta 635 21.2 0.00813 0.995
(a) sessions pooled raw (5, Arena) gamma 635 21.2 0.0117 0.859
(a) sessions pooled raw (5, HCage) theta 576 19.2 0.00856 0.995
(a) sessions pooled raw (5, HCage) gamma 576 19.2 0.0143 0.627
(a) sessions pooled raw (6, Arena) theta 641 21.4 0.00926 0.991
(a) sessions pooled raw (6, Arena) gamma 641 21.4 0.0112 0.905
(a) sessions pooled raw (6, HCage) theta 631 21 0.00881 0.992
(a) sessions pooled raw (6, HCage) gamma 631 21 0.0105 0.958
(b) each session raw (1, Arena, 1) theta 283 9.43 0.0179 0.816
(b) each session raw (1, Arena, 1) gamma 283 9.43 0.0178 0.826
(b) each session raw (1, Arena, 2) theta 281 9.37 0.0166 0.911
(b) each session raw (1, Arena, 2) gamma 281 9.37 0.0154 0.958
(b) each session raw (1, HCage, 1) theta 948 31.6 0.00743 0.992
(b) each session raw (1, HCage, 1) gamma 948 31.6 0.00738 0.992
(b) each session raw (1, HCage, 2) theta 354 11.8 0.0144 0.93
(b) each session raw (1, HCage, 2) gamma 354 11.8 0.0186 0.581
(b) each session raw (2, Arena, 1) theta 302 10.1 0.0194 0.639
(b) each session raw (2, Arena, 1) gamma 302 10.1 0.0135 0.99
(b) each session raw (2, Arena, 2) theta 255 8.5 0.021 0.649
(b) each session raw (2, Arena, 2) gamma 255 8.5 0.0181 0.868
(b) each session raw (2, HCage, 1) theta 623 20.8 0.0107 0.952
(b) each session raw (2, HCage, 1) gamma 623 20.8 0.00917 0.991
(b) each session raw (2, HCage, 2) theta 317 10.6 0.0109 0.996
(b) each session raw (2, HCage, 2) gamma 317 10.6 0.0143 0.965
(b) each session raw (3, Arena, 1) theta 332 11.1 0.0144 0.952
(b) each session raw (3, Arena, 1) gamma 332 11.1 0.0165 0.822
(b) each session raw (3, Arena, 2) theta 301 10 0.0141 0.981
(b) each session raw (3, Arena, 2) gamma 301 10 0.0165 0.884
(b) each session raw (3, HCage, 1) theta 449 15 0.0116 0.983
(b) each session raw (3, HCage, 1) gamma 449 15 0.0117 0.98
(b) each session raw (4, Arena, 1) theta 354 11.8 0.0123 0.991
(b) each session raw (4, Arena, 1) gamma 354 11.8 0.012 0.991
(b) each session raw (4, Arena, 2) theta 335 11.2 0.0147 0.935
(b) each session raw (4, Arena, 2) gamma 335 11.2 0.0131 0.986
(b) each session raw (4, HCage, 1) theta 359 12 0.0174 0.692
(b) each session raw (4, HCage, 1) gamma 359 12 0.0155 0.855
(b) each session raw (4, HCage, 2) theta 289 9.63 0.0138 0.99
(b) each session raw (4, HCage, 2) gamma 289 9.63 0.0139 0.988
(b) each session raw (5, Arena, 1) theta 312 10.4 0.0136 0.986
(b) each session raw (5, Arena, 1) gamma 312 10.4 0.0172 0.81
(b) each session raw (5, Arena, 2) theta 323 10.8 0.0176 0.751
(b) each session raw (5, Arena, 2) gamma 323 10.8 0.0124 0.992
(b) each session raw (5, HCage, 1) theta 290 9.67 0.0213 0.501
(b) each session raw (5, HCage, 1) gamma 290 9.67 0.0283 0.102
(b) each session raw (5, HCage, 2) theta 286 9.53 0.0126 0.993
(b) each session raw (5, HCage, 2) gamma 286 9.53 0.0198 0.654
(b) each session raw (6, Arena, 1) theta 291 9.7 0.0154 0.95
(b) each session raw (6, Arena, 1) gamma 291 9.7 0.015 0.962
(b) each session raw (6, Arena, 2) theta 350 11.7 0.0116 0.993
(b) each session raw (6, Arena, 2) gamma 350 11.7 0.0139 0.956
(b) each session raw (6, HCage, 1) theta 370 12.3 0.0114 0.992
(b) each session raw (6, HCage, 1) gamma 370 12.3 0.0138 0.949
(b) each session raw (6, HCage, 2) theta 261 8.7 0.0163 0.948
(b) each session raw (6, HCage, 2) gamma 261 8.7 0.0159 0.959
(a) sessions pooled, log10 power: 24 tests, smallest p = 0.462, tests with p <= 0.1: 0, with p < 0.05: 0
(a) sessions pooled, raw power: 24 tests, smallest p = 0.554, tests with p <= 0.1: 0, with p < 0.05: 0
(b) each session, log10 power: 46 tests, smallest p = 0.035, tests with p <= 0.1: 1, with p < 0.05: 1
(b) each session, raw power: 46 tests, smallest p = 0.102, tests with p <= 0.1: 0, with p < 0.05: 0
Variant, the left-out file (3, 'HCage', 2) {'theta': 0.993, 'gamma': 0.956}
Result: reproduced. Every Hartigan dip test of theta and gamma power over 2-s epochs, per rat
and condition, gives p > 0.1, with sessions pooled (24 tests, smallest p = 0.554) or taken alone (46
tests, smallest p = 0.102; rat 3 has one home-cage file). The p-values come from the diptest
package's tabulated null distribution. On log10 power, a variant the criterion does not name, the
pooled reading also gives p > 0.1 (smallest 0.462), and each session alone gives one test below 0.05:
rat 5's home-cage session 1, gamma, p = 0.035, the same test that gives the smallest p on raw power.
figs5-arena-excursions¶
During about 20 min of arena exploration, BF gamma made on average eight brief excursions into the home-cage range, lasting 8 to 44 s (16 +/- 8 s), with a clear rise in gamma power (0.56 +/- 0.27) and no change in theta power (-0.01 +/- 0.82) relative to 4 to 8 s before the peak. (Criterion: boxcar of 3, threshold at the intersection of the home-cage and arena gamma histograms, readings for the bins (20 equal bins; kernel density), the power scale (raw, log10, normalized to the rat's mean) and the arena data (two sessions together; each alone, count doubled). Reproduced when, under one reading, count and mean duration are within tolerance, durations within 8 to 44 s widened by 10%, the gamma rise is significant (p < 0.05) and the theta change is not (p > 0.1); partially reproduced when the gamma and theta tests hold and the count or duration does not.)
def scaled(x, scale, ref):
return {"raw": x, "log10": np.log10(x), "normalized": x / ref}[scale]
def run_excursions(scale, method, data):
per_rat, events = [], []
for rat in range(1, 7):
t = main[main.rat == rat]
ref_g, ref_t = t.gamma.mean(), t.theta.mean()
home, ar = t[t.condition == "HCage"], t[t.condition == "Arena"]
if data == "together":
sets = [(home, ar, 1)]
else:
sets = [(home[home.session == s] if (home.session == s).any() else home, ar[ar.session == s], 2) for s in (1, 2)]
for h, a, factor in sets:
find = bflfp.histogram_intersection if method == "20 bins" else bflfp.kde_intersection
thr = find(scaled(h.gamma.values, scale, ref_g), scaled(a.gamma.values, scale, ref_g))
ev = []
for s, g in a.groupby("session"): # files do not share a clock: one series per file
ev += bflfp.excursions(scaled(g.gamma.values, scale, ref_g), scaled(g.theta.values, scale, ref_t),
g.t.values, thr)
per_rat.append({"rat": rat, "threshold": thr, "count_20min": factor * len(ev),
"arena_minutes": len(a) * 2 / 60})
events += [dict(e, rat=rat) for e in ev]
per_rat, events = pd.DataFrame(per_rat), pd.DataFrame(events)
ok = events.dropna(subset=["gamma_before"])
by_rat = ok.groupby("rat")[["gamma_peak", "gamma_before", "theta_peak", "theta_before"]].mean()
g_test = stats.ttest_rel(by_rat.gamma_peak, by_rat.gamma_before)
t_test = stats.ttest_rel(by_rat.theta_peak, by_rat.theta_before)
return per_rat, events, {
"scale": scale, "bins": method, "arena_data": data,
"threshold_range": f"{per_rat.threshold.min():.3g}-{per_rat.threshold.max():.3g}",
"count_per_20min": per_rat.groupby("rat").count_20min.mean().mean(),
"mean_duration_s": events.duration_s.mean(), "sd_duration_s": events.duration_s.std(),
"min_duration_s": events.duration_s.min(), "max_duration_s": events.duration_s.max(),
"n_events": len(events), "events_without_before": len(events) - len(ok),
"gamma_rise": (by_rat.gamma_peak - by_rat.gamma_before).mean(), "gamma_p": g_test.pvalue,
"theta_change": (by_rat.theta_peak - by_rat.theta_before).mean(), "theta_p": t_test.pvalue,
# scale-free comparison with the paper's 0.56 / 2.6: the rise over the mean threshold, shown not judged
"rise_over_threshold": (by_rat.gamma_peak - by_rat.gamma_before).mean() / per_rat.threshold.mean()}
results = []
for scale in ("raw", "log10", "normalized"):
for method in ("20 bins", "kernel density"):
for data in ("together", "each alone"):
results.append(run_excursions(scale, method, data)[2])
exc = pd.DataFrame(results)
exc["count_ok"] = exc.count_per_20min.between(7.2, 8.8)
exc["duration_ok"] = exc.mean_duration_s.between(8, 24) & (exc.min_duration_s >= 7.2) & (exc.max_duration_s <= 48.4)
exc["gamma_ok"] = exc.gamma_p < 0.05
exc["theta_ok"] = exc.theta_p > 0.1
print(exc.to_string(index=False, float_format=lambda x: f"{x:.3g}"))
print("Paper's rise over threshold: 0.56 / 2.6 = %.3f; here, raw and normalized scales: %s"
% (0.56 / 2.6, exc[exc.scale != "log10"].rise_over_threshold.round(3).tolist()))
scale bins arena_data threshold_range count_per_20min mean_duration_s sd_duration_s min_duration_s max_duration_s n_events events_without_before gamma_rise gamma_p theta_change theta_p rise_over_threshold count_ok duration_ok gamma_ok theta_ok
raw 20 bins together 25.3-106 14.5 5.45 4.89 2 28 87 2 18.9 0.0317 11.5 0.0282 0.427 False False True False
raw 20 bins each alone 24.9-109 12.3 5.59 5.06 2 28 74 3 18.9 0.04 -16 0.415 0.429 False False True True
raw kernel density together 24.8-108 13.8 5.35 5.17 2 30 83 2 18.6 0.0269 11.1 0.244 0.417 False False True True
raw kernel density each alone 24.5-108 10.8 5.42 5.39 2 28 65 2 19.2 0.0284 -8.26 0.526 0.438 False False True True
log10 20 bins together 1.39-2.03 13.7 5.12 5.08 2 28 82 2 0.144 0.000417 0.0346 0.193 0.0913 False False True True
log10 20 bins each alone 1.38-2.03 13 5.56 5.44 2 28 78 3 0.164 0.000787 -0.00473 0.894 0.105 False False True True
log10 kernel density together 1.39-2.03 13.7 5.1 4.98 2 30 82 2 0.146 0.000488 0.0315 0.346 0.092 False False True True
log10 kernel density each alone 1.38-2.03 10.7 5.38 5.83 2 28 64 3 0.17 0.00052 -0.028 0.505 0.108 False False True True
normalized 20 bins together 0.699-0.992 14.5 5.45 4.89 2 28 87 2 0.349 0.000852 0.0593 0.04 0.387 False False True False
normalized 20 bins each alone 0.551-1.11 12.3 5.59 5.06 2 28 74 3 0.347 0.00287 -0.0805 0.398 0.391 False False True True
normalized kernel density together 0.674-0.998 13.8 5.35 5.17 2 30 83 2 0.352 0.000112 0.0558 0.256 0.39 False False True True
normalized kernel density each alone 0.541-1.11 10.8 5.42 5.39 2 28 65 2 0.366 0.00116 -0.0428 0.532 0.411 False False True True
Paper's rise over threshold: 0.56 / 2.6 = 0.215; here, raw and normalized scales: [0.427, 0.429, 0.417, 0.438, 0.387, 0.391, 0.39, 0.411]
per_rat, events, _ = run_excursions("log10", "20 bins", "together")
print("log10 scale, 20 bins, sessions together: per rat")
print(per_rat.round(3).to_string(index=False))
print("Durations (s), counts:", events.duration_s.value_counts().sort_index().to_dict())
print("Excursions lasting 8 s or more: %d of %d; their mean count per rat %.1f"
% ((events.duration_s >= 8).sum(), len(events), (events.duration_s >= 8).groupby(events.rat).sum().mean()))
fig, ax = plt.subplots(figsize=(9, 3))
t = main[(main.rat == 1) & (main.condition == "Arena") & (main.session == 1)]
ax.plot(t.t / 60, bflfp.boxcar(np.log10(t.gamma.values)), lw=0.8)
ax.axhline(per_rat[per_rat.rat == 1].threshold.iloc[0], color="k", ls=":")
ax.set(xlabel="time (min)", ylabel="log10 gamma, boxcar 3", title="Rat 1, arena session 1")
plt.show()
print("Plotted range: %.2f to %.2f" % (bflfp.boxcar(np.log10(t.gamma.values)).min(), bflfp.boxcar(np.log10(t.gamma.values)).max()))
log10 scale, 20 bins, sessions together: per rat
rat threshold count_20min arena_minutes
1 1.481 11 18.800
2 1.628 11 18.567
3 1.392 8 21.100
4 1.520 16 22.967
5 2.029 16 21.167
6 1.445 20 21.367
Durations (s), counts: {2.0: 37, 4.0: 17, 6.0: 15, 8.0: 2, 10.0: 4, 12.0: 2, 14.0: 1, 18.0: 1, 22.0: 1, 26.0: 1, 28.0: 1}
Excursions lasting 8 s or more: 13 of 82; their mean count per rat 2.2
Plotted range: 1.15 to 1.60
Result: partially reproduced. Gamma rises at the excursion peak under all 12 readings (paired t test over rats, p from 0.0001 to 0.04), and the theta change is not significant (p > 0.1) under 10 of them, among them every log10 reading; under the raw and normalized scales with 20 bins and the sessions taken together the theta test gives p = 0.028 and 0.040. No reading gives the paper's count or durations: 10.7 to 14.5 excursions per 20 min in place of 8, mean durations of 5.1 to 5.6 s in place of 16, and durations from 2 s, below the paper's shortest of 8 s. No reading puts the threshold near 2.6, so the 0.56 and -0.01 values are reported and not judged. As a scale-free comparison, the gamma rise over the mean threshold is 0.39 to 0.44 on the raw and normalized scales, against 0.56 / 2.6 = 0.215 in the paper (shown, not judged).
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 |
|---|---|
| fig3a-bf-gamma-reduction | partially reproduced (48% reduction, p = 0.012 over 6 rats) |
| fig3a-each-hemisphere-significant | reproduced (12 of 12 recordings p < 0.01) |
| fig3b-arena-gamma-unrelated-to-movement | reproduced |
| fig3b-home-cage-gamma-increases-with-movement | not reproduced by the criterion's test (rise in every rat; repeated-measures p < 0.001) |
| fig3b-matched-movement | not reproduced by the criterion's test (home higher in every rat; one-way p about 0.07) |
| fig3c-home-cage-above-both-sessions | partially reproduced (home cage above both arena sessions by Bonferroni-corrected paired t tests, not by Tukey's HSD; session order not in the data set) |
| fig3d-gamma-rises-with-time-in-arena | partially reproduced (r = 0.86 in 1-min bins; within tolerance in 30-s bins), rests on the criterion's bin width |
| figs2-bf-low-frequency-bands | reproduced |
| figs3-unimodal-states | reproduced |
| figs5-arena-excursions | partially reproduced (more and shorter excursions) |
| bipolar-gamma-reduction, fig3c-consistent-over-days | blocked, metadata insufficient |
| fig3d-rapid-switch-to-home-cage and 15 claims on behaviour, objects, units, ACC and VC | blocked, data missing |
How the analysis differs from the paper's: the data set gives one BF signal per rat for six rats, so tests run over rats (n = 6) or session pairs, where the paper used 22 hemispheres; the segments come already cleaned and down-sampled; the gamma measure is a plain periodogram mean over 30-80 Hz; the data set has no day labels and no order of sessions within a day.
What was found about the data set: rat 3's home-cage session 2 is a second deposit of 287 segments
of its session 1; the files hold more removed segments (2.6% to 23.5%) than the SI's 1-4%; each file
has its own clock, so no recording spans the transfer between arena and home cage; the LFP unit is
not stated; movement values differ in scale between rats; the descriptor names Hcage where files
say HCage, and calls the archive a file to unpack where it is a gzipped tar; the data set holds
only the BF LFP of six rats, not the ACC, VC, unit, object or behaviour data of the paper.