Reproducing Grosmark & Buzsáki (2016) from dandiset 000044¶
Paper: A. D. Grosmark and G. Buzsáki, "Diversity in neural firing dynamics supports both rigid and learned hippocampal sequences", Science 351:1440 (2016), doi:10.1126/science.aad1935.
Dandiset: DANDI:000044, version
0.261005.2101. The paper cites no version (it predates the dandiset and points to CRCNS). This
is the most recent published version. It is also the one that matches the files best: in the two
earlier versions (0.210812.1516 and 0.250624.0426) the rate of every position series holds the
sampling period (0.0256 s) in place of the rate (39.06 Hz), so the positions would span 880 to
2360 hours in place of the MAZE epoch. The cell below shows this. Otherwise the versions differ only
in subject metadata.
Text used: the full main text and the legends of Figures 1 to 4, from the NIH author manuscript in PubMed Central (PMC4919122). The supplementary Materials and Methods and figs. S1 to S22 could not be retrieved (PMC and Science serve them only behind a browser check), so every definition the main text does not give is a common choice, fixed in the claims file before the analysis. Neither the paper nor the dandiset links analysis code.
Eight claims rest on the selection of strongly contributing cells (PCC). The test fixed for it
before the analysis has a false-positive rate of 56%, far from the 5% that "significant" implies,
so these claims are blocked as methods_underspecified; their analyses are still shown, with a
best-template variant of the contribution whose false-positive rate (21%) is also far above 5%.
fig2d-shuffle-high-dpcc ranks all place cells by dPCC and does not use that selection, so it is
tested.
| Claim | Status |
|---|---|
| sessions-and-tracks: 4 rats, 8 sessions, two linear and one circular track | tested |
| place-cell-count: 491 place cells | tested |
| fig1c-post-gt-pre: sequence scores higher in POST than PRE | tested |
| fig1c-maze-gt-pre: sequence scores higher in MAZE than PRE | tested |
| fig2b-strong-contributors: 216 strongly contributing neurons | blocked, methods_underspecified (analysis shown) |
| fig2-pre-post-overlap: 73% of PRE contributors also contribute in POST | blocked, methods_underspecified (analysis shown) |
| fig2d-shuffle-high-dpcc: shuffling the top 10 to 20% dPCC cells removes replay | tested |
| fig2e-remove-plastic: removing plastic cells removes replay | blocked, methods_underspecified (analysis shown) |
| figs2-track-distribution: rigid and plastic cells similarly placed on the track | blocked, methods_underspecified (analysis shown) |
| fig3a-plastic-specificity: plastic cells more specific, fewer fields | blocked, methods_underspecified (analysis shown) |
| fig3b-plastic-improve: plastic cells sharpen during the session | blocked, methods_underspecified (analysis shown) |
| fig4a-rigid-faster: rigid cells fire faster | blocked, methods_underspecified (analysis shown) |
| fig4a-plastic-ripple-burst-coactivation: plastic cells more ripple-related, bursty and coactive | blocked, methods_underspecified (analysis shown) |
| fig4c-pre-fast-coactivation: PRE coactivation dominated by fast cells | tested |
| fig4c-slow-coactivation-increase: slow cells gain coactivation PRE to POST | tested |
| fig4c-slow-replay: slow cells show most replay of coactivation | tested |
| rate-split-medians: low and high-rate medians 0.39 and 1.12 Hz | tested |
| figs16-low-rate-info-gain: low-rate cells more informative, larger ripple-gain increase | tested |
| fig4d-slow-sequences: slow but not fast cells increase sequence content | tested |
| fig4e-gain-sequences: increase limited to high ripple-gain cells | tested |
import os
for k in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
os.environ[k] = "1"
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
from concurrent.futures import ProcessPoolExecutor
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import grosmark as g
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
print("dandiset", g.DANDISET, "version", g.VERSION)
dandiset 000044 version 0.261005.2101
Versions¶
The rate of the position series in the first published version and in the pinned one, against the length of the MAZE epoch.
OLD = dict(g.asset_urls("0.210812.1516"))
URLS = g.asset_urls()
rows = []
for path, url in URLS:
with g.open_nwb(url) as f:
beh = f["processing/behavior"]
name = [k for k in beh if "Linearized" in k][0]
ser = beh[name][list(beh[name])[0]]
n = ser["data"].shape[0]
rate_new = float(ser["starting_time"].attrs["rate"])
ep = f["intervals/epochs"]
maze = float(ep["stop_time"][1] - ep["start_time"][1])
with g.open_nwb(OLD[path]) as f:
beh = f["processing/behavior"]
ser = beh[name][list(beh[name])[0]]
rate_old = float(ser["starting_time"].attrs["rate"])
rows.append(dict(file=path.split("/")[1][:30], samples=n, maze_s=round(maze, 1),
rate_0_210812=rate_old, span_old_h=round(n / rate_old / 3600, 1),
rate_pinned=round(rate_new, 3), span_pinned_s=round(n / rate_new, 1)))
print(pd.DataFrame(rows).to_string(index=False))
file samples maze_s rate_0_210812 span_old_h rate_pinned span_pinned_s sub-Achilles_ses-Achilles-1025 80762 2067.5 0.0256 876.3 39.063 2067.5 sub-Achilles_ses-Achilles-1101 104457 2674.1 0.0256 1133.4 39.062 2674.1 sub-Buddy_ses-Buddy-06272013_b 90942 2328.1 0.0256 986.8 39.063 2328.1 sub-Cicero_ses-Cicero-09012014 217235 5561.2 0.0256 2357.2 39.063 5561.2 sub-Cicero_ses-Cicero-09102014 118360 3030.0 0.0256 1284.3 39.063 3030.0 sub-Cicero_ses-Cicero-09172014 120227 3077.8 0.0256 1304.6 39.063 3077.8 sub-Gatsby_ses-Gatsby-08022013 91215 2335.1 0.0256 989.7 39.063 2335.1 sub-Gatsby_ses-Gatsby-08282013 105180 2692.6 0.0256 1141.3 39.063 2692.6
Data check¶
Each file is streamed and reduced in parallel: units, epochs, states, positions, one LFP channel (the good channel with the most 150 to 250 Hz power in a POST non-REM chunk), the ripple envelope at the candidate events, and the checks below. This takes a few minutes.
import time
T0 = time.time()
with ProcessPoolExecutor(4) as ex:
S = list(ex.map(g.prepare_session, *zip(*URLS)))
print(f"loaded in {time.time() - T0:.0f} s")
NAMES = [s["session"] for s in S]
print(len(S), "sessions loaded:", NAMES)
loaded in 163 s 8 sessions loaded: ['Achilles_10252013', 'Achilles_11012013', 'Buddy_06272013', 'Cicero_09012014', 'Cicero_09102014', 'Cicero_09172014', 'Gatsby_08022013', 'Gatsby_08282013']
rows = []
for s in S:
ep = s["epochs"]
ct = pd.Series(s["cell_type"]).value_counts()
loc = pd.Series(s["location"][s["cell_type"] == "excitatory"]).value_counts()
t, x, v = g.lin_speed(s)
tx, sp = g.xy_speed(s)
run = np.isfinite(x) & (np.abs(np.nan_to_num(v)) > g.RUN_SPEED)
last_spike = max(q.max() for q in s["spikes"] if len(q))
rows.append(dict(
subject=s["subject"], session=s["session"], track=s["track"],
PRE_h=round((ep["PREEpoch"][1] - ep["PREEpoch"][0]) / 3600, 2),
MAZE_min=round((ep["MazeEpoch"][1] - ep["MazeEpoch"][0]) / 60, 1),
POST_h=round((ep["POSTEpoch"][1] - ep["POSTEpoch"][0]) / 3600, 2),
epochs_end_s=round(ep["POSTEpoch"][1], 1), lfp_end_s=round(s["lfp_seconds"], 1),
last_spike_s=round(last_spike, 1),
pyr=int(ct.get("excitatory", 0)), int_=int(ct.get("inhibitory", 0)),
pyr_lCA1=int(loc.get("lCA1", 0)), pyr_rCA1=int(loc.get("rCA1", 0)),
pos_rate=round(s["lin"]["rate"], 2), pos_unit=s["lin"]["unit"],
lin_max_m=round(float(np.nanmax(s["lin"]["data"])), 3),
on_track_frac=round(float(np.isfinite(x).mean()), 3), run_s=round(run.sum() / s["lin"]["rate"], 0),
lin_speed_run=round(float(np.nanmedian(np.abs(v[run]))), 1),
xy_speed_run=round(float(np.nanmedian(sp[run])), 1),
run_up_frac=round(float(np.mean(v[np.isfinite(v)] > g.RUN_SPEED)), 3),
run_down_frac=round(float(np.mean(v[np.isfinite(v)] < -g.RUN_SPEED)), 3),
directions=str(g.directions(s)),
))
DATA = pd.DataFrame(rows)
print(DATA.to_string(index=False))
print()
print("subjects", DATA.subject.nunique(), "sessions", len(DATA), "pyramidal", DATA.pyr.sum(),
"interneurons", DATA.int_.sum())
print("sessions per subject:", DATA.groupby("subject").size().to_dict())
subject session track PRE_h MAZE_min POST_h epochs_end_s lfp_end_s last_spike_s pyr int_ pyr_lCA1 pyr_rCA1 pos_rate pos_unit lin_max_m on_track_frac run_s lin_speed_run xy_speed_run run_up_frac run_down_frac directions Achilles Achilles_10252013 1.6mLinearMaze 5.02 34.5 4.09 34861.1 34861.1 34861.1 120 17 57 63 39.06 meters 1.600 0.128 244.0 49.8 52.0 0.512 0.411 [1, -1] Achilles Achilles_11012013 CircularMaze 5.37 44.6 4.06 36629.2 36629.2 36629.2 92 12 38 54 39.06 meters 2.900 0.585 984.0 24.0 25.3 0.156 0.474 [-1] Buddy Buddy_06272013 1.6mLinearMaze 2.98 38.8 2.22 21054.7 21054.7 21054.7 48 20 46 2 39.06 meters 1.600 0.072 161.0 53.7 56.9 0.457 0.502 [1, -1] Cicero Cicero_09012014 1.6mLinearMaze 4.31 92.7 3.80 34756.8 34756.8 34756.8 55 18 27 28 39.06 meters 1.600 0.218 632.0 18.9 20.5 0.265 0.255 [1, -1] Cicero Cicero_09102014 CircularMaze 4.79 50.5 4.51 36498.6 36498.6 36498.6 81 24 21 60 39.06 meters 2.842 0.287 603.0 22.5 24.7 0.066 0.627 [-1] Cicero Cicero_09172014 2mLinearMaze 4.44 51.3 4.02 33560.2 33560.2 33560.2 59 13 22 37 39.06 meters 2.000 0.321 599.0 16.7 18.6 0.293 0.313 [1, -1] Gatsby Gatsby_08022013 1.6mLinearMaze 4.74 38.9 3.06 30413.6 32002.4 32002.4 66 14 18 48 39.06 meters 1.600 0.187 340.0 36.8 37.9 0.413 0.367 [1, -1] Gatsby Gatsby_08282013 CircularMaze 4.66 44.9 4.42 35347.9 35347.9 35347.9 41 10 41 0 39.06 meters 2.841 0.363 669.0 33.5 34.3 0.101 0.583 [-1]
subjects 4 sessions 8 pyramidal 562 interneurons 128
sessions per subject: {'Achilles': 2, 'Buddy': 1, 'Cicero': 3, 'Gatsby': 2}
# sleep states, LFP channel and clock checks
rows = []
for s in S:
st = s["states"]
dur = pd.Series(st["stop"] - st["start"]).groupby(st["label"]).sum()
td0, td1 = s["theta_delta"][0.0], s["theta_delta"][600.0]
ov = s["ripple_overlap"]
ev = pd.Series(s["events"]["epoch"]).value_counts()
rows.append(dict(session=s["session"], states=",".join(sorted(set(st["label"]))),
awake_h=round(dur.get("Awake", 0) / 3600, 2), nrem_h=round(dur.get("Non-REM", 0) / 3600, 2),
rem_h=round(dur.get("REM", 0) / 3600, 2),
states_end_s=round(float(st["stop"].max()), 0),
channel=s["ripple_channel"], bad_ch=int(s["bad_electrode"].sum()),
TD_REM=round(td0["REM"], 2), TD_NREM=round(td0["Non-REM"], 2),
TD_REM_shift600=round(td1["REM"], 2), TD_NREM_shift600=round(td1["Non-REM"], 2),
burst_rip=round(ov[0.0], 3), burst_rip_shift05=round(ov[0.5], 3),
burst_rip_shift2=round(ov[2.0], 3),
ev_PRE=int(ev.get("PRE", 0)), ev_MAZE=int(ev.get("MAZE", 0)), ev_POST=int(ev.get("POST", 0))))
print(pd.DataFrame(rows).to_string(index=False))
session states awake_h nrem_h rem_h states_end_s channel bad_ch TD_REM TD_NREM TD_REM_shift600 TD_NREM_shift600 burst_rip burst_rip_shift05 burst_rip_shift2 ev_PRE ev_MAZE ev_POST Achilles_10252013 Awake,Non-REM,REM 4.65 4.41 0.58 34856.0 2 0 4.06 0.30 0.33 0.38 0.733 0.141 0.103 2839 258 3382 Achilles_11012013 Awake,Non-REM,REM 5.29 4.47 0.33 36624.0 3 0 1.05 0.26 0.24 0.25 0.766 0.161 0.110 2423 195 4138 Buddy_06272013 Awake,Non-REM,REM 3.35 2.21 0.22 20995.0 14 3 3.06 0.26 0.25 0.37 0.650 0.168 0.153 1533 121 1458 Cicero_09012014 Awake,Non-REM,REM 4.81 3.94 0.84 34751.0 78 0 3.31 0.26 0.34 0.40 0.656 0.115 0.105 1709 1009 2993 Cicero_09102014 Awake,Non-REM,REM 5.75 3.65 0.67 36493.0 78 0 6.47 0.28 0.32 0.42 0.670 0.124 0.120 1291 83 3302 Cicero_09172014 Awake,Non-REM,REM 5.47 3.19 0.59 33555.0 8 0 0.61 0.15 0.16 0.18 0.679 0.120 0.104 1524 266 2677 Gatsby_08022013 Awake,Non-REM,REM 4.94 3.47 0.42 31997.0 36 0 3.31 0.29 0.34 0.33 0.631 0.103 0.119 3421 399 2723 Gatsby_08282013 Awake,Non-REM,REM 4.78 4.50 0.46 35342.0 53 3 1.62 0.31 0.34 0.31 0.712 0.154 0.140 2670 289 4379
# spikes against positions: median error (cm) of decoding position during running, rate maps from
# even minutes and 0.5 s windows from odd minutes, with the spikes shifted by each lag (s)
rows = []
for s in S:
e = g.position_error(s)
rows.append(dict(session=s["session"], **{f"lag {k:+g}": round(v, 1) for k, v in e.items()}))
POSERR = pd.DataFrame(rows)
print(POSERR.to_string(index=False))
print("lag with the smallest error per session:",
POSERR.set_index("session").astype(float).idxmin(axis=1).to_dict())
# recordings deposited twice: identical spike trains or positions across files
sig_units = {}
for s in S:
for u, q in zip(s["unit_id"], s["spikes"]):
sig_units.setdefault((len(q), round(float(q[:50].sum()), 5)), []).append((s["session"], int(u)))
dup_units = [v for v in sig_units.values() if len({x[0] for x in v}) > 1]
sig_pos = {}
for s in S:
sig_pos.setdefault(round(float(np.nansum(s["xy"]["data"])), 4), []).append(s["session"])
print("unit spike trains found in more than one file:", len(dup_units))
print("position series found in more than one file:", [v for v in sig_pos.values() if len(v) > 1])
session lag -30 lag -5 lag -1 lag +0 lag +1 lag +5 lag +30
Achilles_10252013 41.2 34.0 11.0 6.6 14.8 43.2 47.5
Achilles_11012013 83.0 59.1 18.1 5.6 24.7 65.3 78.1
Buddy_06272013 41.4 46.1 17.5 11.7 24.8 45.8 41.5
Cicero_09012014 47.6 33.7 16.3 10.2 17.1 39.7 48.3
Cicero_09102014 71.2 43.0 10.2 5.0 12.2 46.7 69.4
Cicero_09172014 53.3 32.6 13.3 8.8 20.0 31.5 49.6
Gatsby_08022013 43.5 34.9 12.0 9.2 17.2 43.2 45.5
Gatsby_08282013 79.8 65.1 16.4 4.9 21.6 52.5 79.9
lag with the smallest error per session: {'Achilles_10252013': 'lag +0', 'Achilles_11012013': 'lag +0', 'Buddy_06272013': 'lag +0', 'Cicero_09012014': 'lag +0', 'Cicero_09102014': 'lag +0', 'Cicero_09172014': 'lag +0', 'Gatsby_08022013': 'lag +0', 'Gatsby_08282013': 'lag +0'}
unit spike trains found in more than one file: 0
position series found in more than one file: []
Checks of the description against the data, from the outputs above:
- Four subjects and eight sessions, each with PRE, MAZE and POST epochs, as the metadata state. Durations are printed; the description's "~4 hour" PRE and POST and "~45 minute" MAZE are compared there.
- Units are labelled
excitatoryorinhibitory; pyramidal cells here are the excitatory units. Their hemisphere is inlocation(lCA1 or rCA1); the electrodes'locationis "unknown". - Positions are stated in metres. The linearized position spans 0 to 1.6 m on the 1.6 m tracks and 0 to 2 m on the 2 m track, consistent with the description. On the circular track (described as 1 m in diameter, a circumference of 3.14 m) it reaches 2.84 to 2.9 m, consistent with the series description, which has it run from one edge of the reward area to the other. The running speed from the 2D positions agrees with that from the linearized position, so the 2D positions are consistent with metres.
- The sleep states are consistent with the LFP: the theta/delta ratio is far higher in REM than in non-REM in every session, and the contrast is lost when the labels are shifted by 600 s.
- Spikes and LFP share a clock: the population bursts in non-REM hold a ripple far more often with the envelope as stored than shifted by 0.5 or 2 s.
- Spikes and positions share a clock: the decoding error of position is smallest at lag 0.
- No unit spike train or position series appears in more than one file.
- In Gatsby_08022013 the POST epoch ends at 30413.6 s, while the LFP, the spikes and the sleep states run to about 32000 s. Firing rates over the session use the LFP duration.
- On the circular track the linearized position mostly decreases while the animal runs
(
run_down_fracagainstrun_up_frac), so its single template is built from running toward decreasing position. The description says the position increases clockwise, not which way the animals ran. - The description mentions neck EMG and accelerometer signals; the files hold neither, and no raw
wide-band data (
acquisitionis empty). These are not needed here.
Templates, event scores and per-cell contributions¶
For each session: rate maps on running in each direction, place cells, the event scores (rZ) of every candidate event, each active cell's contribution to its events (for PCC), and the same contributions with every rate map shifted (the false-positive rate of the PCC test). Seeds are fixed per session and per step. The contributions are computed twice: as the criterion defines them, on the template the event was scored with, and as a variant on the better of the event's templates (see fig2b-strong-contributors).
T0 = time.time()
with ProcessPoolExecutor(4) as ex:
FP = list(ex.map(g.first_pass, S, [1000 + i for i in range(8)], [2000 + i for i in range(8)],
[3000 + i for i in range(8)]))
print(f"first pass in {time.time() - T0:.0f} s")
EVENTS = pd.concat([pd.DataFrame(dict(session=s["session"], subject=s["subject"], epoch=fp["decoded"]["epoch"],
rZ=fp["decoded"]["rZ"], r=fp["decoded"]["r"],
n_active=fp["decoded"]["n_active"]))
for s, fp in zip(S, FP)], ignore_index=True)
print("candidate events scored, by session and epoch:")
print(EVENTS.groupby(["session", "epoch"]).size().unstack().to_string())
cand = pd.DataFrame([pd.Series(s["events"]["epoch"]).value_counts() for s in S], index=NAMES)
print("candidate events before the active-cell condition:")
print(cand.to_string())
first pass in 260 s
candidate events scored, by session and epoch:
epoch MAZE POST PRE
session
Achilles_10252013 255 3371 2831
Achilles_11012013 193 4079 2400
Buddy_06272013 76 1136 1246
Cicero_09012014 953 2870 1646
Cicero_09102014 78 3244 1284
Cicero_09172014 264 2661 1519
Gatsby_08022013 376 2602 3294
Gatsby_08282013 265 3838 2228
candidate events before the active-cell condition:
POST PRE MAZE
Achilles_10252013 3382 2839 258
Achilles_11012013 4138 2423 195
Buddy_06272013 1458 1533 121
Cicero_09012014 2993 1709 1009
Cicero_09102014 3302 1291 83
Cicero_09172014 2677 1524 266
Gatsby_08022013 2723 3421 399
Gatsby_08282013 4379 2670 289
sessions-and-tracks¶
Recordings were made in four rats, in eight novel exploration sessions, each a PRE sleep, a MAZE epoch on one of two linear tracks or a circular track, and a POST sleep.
order_ok = [s["epochs"]["PREEpoch"][1] <= s["epochs"]["MazeEpoch"][0] + 1e-6 and
s["epochs"]["MazeEpoch"][1] <= s["epochs"]["POSTEpoch"][0] + 1e-6 for s in S]
print(DATA[["subject", "session", "track", "PRE_h", "MAZE_min", "POST_h"]].to_string(index=False))
print("subjects:", DATA.subject.nunique(), " sessions:", len(DATA), " epochs in PRE, MAZE, POST order:",
sum(order_ok), "of", len(S))
print("tracks:", DATA.track.value_counts().to_dict())
subject session track PRE_h MAZE_min POST_h
Achilles Achilles_10252013 1.6mLinearMaze 5.02 34.5 4.09
Achilles Achilles_11012013 CircularMaze 5.37 44.6 4.06
Buddy Buddy_06272013 1.6mLinearMaze 2.98 38.8 2.22
Cicero Cicero_09012014 1.6mLinearMaze 4.31 92.7 3.80
Cicero Cicero_09102014 CircularMaze 4.79 50.5 4.51
Cicero Cicero_09172014 2mLinearMaze 4.44 51.3 4.02
Gatsby Gatsby_08022013 1.6mLinearMaze 4.74 38.9 3.06
Gatsby Gatsby_08282013 CircularMaze 4.66 44.9 4.42
subjects: 4 sessions: 8 epochs in PRE, MAZE, POST order: 8 of 8
tracks: {'1.6mLinearMaze': 4, 'CircularMaze': 3, '2mLinearMaze': 1}
Result: reproduced. The dandiset holds 4 subjects and 8 sessions, each with PRE, MAZE and POST epochs in that order (8 of 8), on a 1.6 m linear track (4 sessions), a 2 m linear track (1) and a circular track (3). The rooms named in the paper are not in the files.
place-cell-count¶
The Bayesian decoder was built from 491 place cells over the eight sessions.
A place cell is a pyramidal cell whose smoothed rate map on running (above 5 cm/s) peaks at 1 Hz or more in at least one direction. 0.5 and 2 Hz are shown as variants.
rows = []
for s, fp in zip(S, FP):
pk = fp["tpl"]["peak"]
rows.append(dict(subject=s["subject"], session=s["session"], pyramidal=len(pk),
place_05Hz=int((pk >= 0.5).sum()), place_1Hz=int((pk >= 1).sum()),
place_2Hz=int((pk >= 2).sum()),
**{f"template_{d:+d}": len(c) for d, c in fp["cells"].items()}))
PC = pd.DataFrame(rows)
print(PC.to_string(index=False))
print("totals:", PC[["pyramidal", "place_05Hz", "place_1Hz", "place_2Hz"]].sum().to_dict())
print("by subject (1 Hz):", PC.groupby("subject").place_1Hz.sum().to_dict())
print("expected 491, tolerance 49; within:", abs(PC.place_1Hz.sum() - 491) <= 49)
subject session pyramidal place_05Hz place_1Hz place_2Hz template_+1 template_-1
Achilles Achilles_10252013 120 112 99 87 81.0 84
Achilles Achilles_11012013 92 92 82 75 NaN 82
Buddy Buddy_06272013 48 35 30 24 23.0 26
Cicero Cicero_09012014 55 41 31 24 26.0 27
Cicero Cicero_09102014 81 63 54 45 NaN 54
Cicero Cicero_09172014 59 50 42 34 38.0 40
Gatsby Gatsby_08022013 66 49 44 33 36.0 38
Gatsby Gatsby_08282013 41 38 35 30 NaN 35
totals: {'pyramidal': 562, 'place_05Hz': 480, 'place_1Hz': 417, 'place_2Hz': 352}
by subject (1 Hz): {'Achilles': 181, 'Buddy': 30, 'Cicero': 127, 'Gatsby': 79}
expected 491, tolerance 49; within: False
Result: not reproduced. With a 1 Hz peak, 417 of the 562 pyramidal cells are place cells, below the paper's 491 by more than the tolerance of 49. With a 0.5 Hz peak the count is 480, within the tolerance, and with 2 Hz it is 352. The criterion fixes 1 Hz, and the other thresholds are variants that do not decide the outcome; the count depends strongly on the threshold. The paper's definition is in the Methods, which were not available.
fig1c-post-gt-pre and fig1c-maze-gt-pre¶
Sequence scores (rZ) of candidate events are higher in POST sleep than in PRE sleep (replay). Sequence content is higher in MAZE immobility events than in PRE sleep events.
Kruskal-Wallis across the three epochs on all events pooled, then Tukey's HSD on the ranks.
def ranks_tukey(groups):
allv = np.concatenate(groups)
rk = stats.rankdata(allv)
out, i = [], 0
for gr in groups:
out.append(rk[i:i + len(gr)])
i += len(gr)
return stats.tukey_hsd(*out)
ep_names = ["PRE", "MAZE", "POST"]
grp = [EVENTS.loc[EVENTS.epoch == e, "rZ"].dropna().values for e in ep_names]
for e, v in zip(ep_names, grp):
print(f"{e:5s} n={len(v):5d} mean rZ={v.mean():.3f} SE={v.std(ddof=1) / np.sqrt(len(v)):.3f} "
f"median={np.median(v):.3f} range {v.min():.2f} to {v.max():.2f}; sign-rank vs 0 p={stats.wilcoxon(v).pvalue:.2g}")
kw = stats.kruskal(*grp)
print(f"Kruskal-Wallis H={kw.statistic:.1f} p={kw.pvalue:.2g}")
tk = ranks_tukey(grp)
for i, j in [(2, 0), (1, 0), (2, 1)]:
print(f"Tukey on ranks {ep_names[i]} vs {ep_names[j]}: p={tk.pvalue[i, j]:.3g}, "
f"mean difference {grp[i].mean() - grp[j].mean():+.3f}")
per = EVENTS.groupby(["subject", "session", "epoch"]).rZ.mean().unstack()[ep_names]
per["POST>PRE"] = per.POST > per.PRE
per["MAZE>PRE"] = per.MAZE > per.PRE
print(per.round(3).to_string())
print("sessions POST>PRE:", int(per["POST>PRE"].sum()), "of 8; MAZE>PRE:", int(per["MAZE>PRE"].sum()), "of 8")
print("by subject:")
print(EVENTS.groupby(["subject", "epoch"]).rZ.mean().unstack()[ep_names].round(3).to_string())
fig, ax = plt.subplots(figsize=(5, 3.5))
for e, v in zip(ep_names, grp):
ax.plot(np.sort(v), np.linspace(0, 1, len(v)), label=e)
ax.set_xlabel("rZ"); ax.set_ylabel("cumulative fraction"); ax.legend()
ax.set_xlim(-3, 8)
print("points outside the x axis:", int(sum(((v < -3) | (v > 8)).sum() for v in grp)))
PRE n=16448 mean rZ=0.770 SE=0.010 median=0.682 range -1.85 to 6.80; sign-rank vs 0 p=0
MAZE n= 2460 mean rZ=1.270 SE=0.029 median=1.213 range -1.70 to 6.10; sign-rank vs 0 p=5.5e-254
POST n=23801 mean rZ=0.763 SE=0.008 median=0.677 range -1.86 to 5.94; sign-rank vs 0 p=0
Kruskal-Wallis H=290.2 p=9.6e-64
Tukey on ranks POST vs PRE: p=0.412, mean difference -0.007
Tukey on ranks MAZE vs PRE: p=0, mean difference +0.500
Tukey on ranks POST vs MAZE: p=0, mean difference -0.507
epoch PRE MAZE POST POST>PRE MAZE>PRE
subject session
Achilles Achilles_10252013 1.136 1.677 1.334 True True
Achilles_11012013 0.170 0.526 0.326 True True
Buddy Buddy_06272013 1.024 1.580 1.209 True True
Cicero Cicero_09012014 1.103 1.423 1.251 True True
Cicero_09102014 0.266 0.491 0.338 True True
Cicero_09172014 0.894 1.337 0.912 True True
Gatsby Gatsby_08022013 1.231 1.749 1.296 True True
Gatsby_08282013 0.090 0.261 0.125 True True
sessions POST>PRE: 8 of 8; MAZE>PRE: 8 of 8
by subject:
epoch PRE MAZE POST
subject
Achilles 0.693 1.181 0.782
Buddy 1.024 1.580 1.209
Cicero 0.790 1.349 0.811
Gatsby 0.771 1.134 0.598
points outside the x axis: 0
The weighted correlation threshold is not used here; as a variant, the fraction of events with rZ above 1.96 per epoch.
print(EVENTS.assign(sig=EVENTS.rZ > 1.96).groupby("epoch").sig.mean().round(3).to_dict())
{'MAZE': 0.316, 'POST': 0.191, 'PRE': 0.189}
Variant, within sessions: the pooled test mixes sessions whose rZ differ in level (the circular track sessions, with one template, score lower in every epoch) and in the number of events. A one-sided rank-sum test of POST against PRE within each session, a sign test over the 8 sessions, and the pooled comparison after subtracting each session's PRE mean.
rows = []
for (sb, ss), df in EVENTS.groupby(["subject", "session"]):
a = df.loc[df.epoch == "POST", "rZ"]; b = df.loc[df.epoch == "PRE", "rZ"]; m = df.loc[df.epoch == "MAZE", "rZ"]
rows.append(dict(subject=sb, session=ss, POST_minus_PRE=a.mean() - b.mean(),
p_POST_gt_PRE=stats.mannwhitneyu(a, b, alternative="greater").pvalue,
MAZE_minus_PRE=m.mean() - b.mean(),
p_MAZE_gt_PRE=stats.mannwhitneyu(m, b, alternative="greater").pvalue))
W = pd.DataFrame(rows)
print(W.round(4).to_string(index=False))
n_pos = int((W.POST_minus_PRE > 0).sum())
print(f"sessions with POST > PRE: {n_pos} of 8, sign test p = {stats.binomtest(n_pos, 8, alternative='greater').pvalue:.4f}; "
f"with one-sided p < 0.05: {int((W.p_POST_gt_PRE < 0.05).sum())}")
cen = EVENTS.copy()
cen["rZc"] = cen.rZ - cen.groupby("session").rZ.transform(lambda v: v[cen.loc[v.index, "epoch"] == "PRE"].mean())
a = cen.loc[cen.epoch == "POST", "rZc"]; b = cen.loc[cen.epoch == "PRE", "rZc"]
print(f"pooled after subtracting each session's PRE mean: POST {a.mean():+.3f}, rank-sum one-sided p = "
f"{stats.mannwhitneyu(a, b, alternative='greater').pvalue:.3g}")
print("events per epoch on circular and linear tracks:")
print(EVENTS.merge(DATA[["session", "track"]]).assign(circular=lambda d: d.track == "CircularMaze")
.groupby(["circular", "epoch"]).rZ.agg(["size", "mean"]).round(3).to_string())
subject session POST_minus_PRE p_POST_gt_PRE MAZE_minus_PRE p_MAZE_gt_PRE
Achilles Achilles_10252013 0.1984 0.0000 0.5418 0.0000
Achilles Achilles_11012013 0.1555 0.0000 0.3562 0.0005
Buddy Buddy_06272013 0.1847 0.0005 0.5561 0.0001
Cicero Cicero_09012014 0.1481 0.0002 0.3203 0.0000
Cicero Cicero_09102014 0.0713 0.0454 0.2248 0.0405
Cicero Cicero_09172014 0.0184 0.3815 0.4433 0.0000
Gatsby Gatsby_08022013 0.0650 0.0166 0.5186 0.0000
Gatsby Gatsby_08282013 0.0349 0.0770 0.1713 0.0290
sessions with POST > PRE: 8 of 8, sign test p = 0.0039; with one-sided p < 0.05: 6
pooled after subtracting each session's PRE mean: POST +0.106, rank-sum one-sided p = 2.34e-14
events per epoch on circular and linear tracks:
size mean
circular epoch
False MAZE 1924 1.515
POST 12640 1.207
PRE 10536 1.112
True MAZE 536 0.390
POST 11161 0.260
PRE 5912 0.161
Variant, by track: the criterion's test on the five linear-track sessions alone and on the three circular-track sessions alone. Then, to compare the event pipeline with the paper's, the number of scored events per session and epoch and the distribution of MAZE rZ in each session.
ev_t = EVENTS.merge(DATA[["session", "track"]])
ev_t["linear"] = ~ev_t.track.str.startswith("Circular")
for lin, df in ev_t.groupby("linear"):
grp_t = [df.loc[df.epoch == e, "rZ"].dropna().values for e in ep_names]
tk_t = ranks_tukey(grp_t)
print(("linear" if lin else "circular") + f" tracks ({df.session.nunique()} sessions): Kruskal-Wallis p={stats.kruskal(*grp_t).pvalue:.2g}; "
+ "; ".join(f"{ep_names[i]} {grp_t[i].mean():.3f} vs {ep_names[j]} {grp_t[j].mean():.3f}, Tukey on ranks p={tk_t.pvalue[i, j]:.3g}"
for i, j in [(2, 0), (1, 0)]))
print()
print("scored events per session and epoch:")
print(EVENTS.groupby(["session", "epoch"]).size().unstack()[ep_names].to_string())
print("MAZE rZ per session:")
print(EVENTS[EVENTS.epoch == "MAZE"].groupby("session").rZ.describe(percentiles=[0.1, 0.5, 0.9])
.round(3).to_string())
circular tracks (3 sessions): Kruskal-Wallis p=6.5e-07; POST 0.260 vs PRE 0.161, Tukey on ranks p=1.13e-05; MAZE 0.390 vs PRE 0.161, Tukey on ranks p=0.000695
linear tracks (5 sessions): Kruskal-Wallis p=7.1e-32; POST 1.207 vs PRE 1.112, Tukey on ranks p=5.1e-07; MAZE 1.515 vs PRE 1.112, Tukey on ranks p=0
scored events per session and epoch:
epoch PRE MAZE POST
session
Achilles_10252013 2831 255 3371
Achilles_11012013 2400 193 4079
Buddy_06272013 1246 76 1136
Cicero_09012014 1646 953 2870
Cicero_09102014 1284 78 3244
Cicero_09172014 1519 264 2661
Gatsby_08022013 3294 376 2602
Gatsby_08282013 2228 265 3838
MAZE rZ per session:
count mean std min 10% 50% 90% max
session
Achilles_10252013 255.0 1.677 1.332 -1.264 -0.171 1.587 3.518 4.742
Achilles_11012013 193.0 0.526 1.342 -1.611 -1.102 0.426 2.247 4.663
Buddy_06272013 76.0 1.580 1.183 -1.034 0.088 1.558 3.056 4.499
Cicero_09012014 953.0 1.423 1.351 -1.700 -0.315 1.388 3.189 6.099
Cicero_09102014 78.0 0.491 1.173 -1.463 -0.857 0.392 2.276 3.516
Cicero_09172014 264.0 1.337 1.361 -1.242 -0.300 1.251 3.228 5.417
Gatsby_08022013 376.0 1.749 1.396 -1.477 -0.174 1.793 3.492 5.516
Gatsby_08282013 265.0 0.261 1.226 -1.481 -1.188 0.040 1.971 4.981
Result: fig1c-maze-gt-pre reproduced; fig1c-post-gt-pre not reproduced. MAZE events score higher than PRE events (mean rZ 1.270 against 0.770; Tukey on ranks p printed as 0), and MAZE is higher in all 8 sessions. Pooled over sessions, as the criterion fixes the test, POST does not exceed PRE (0.763 against 0.770, Tukey on ranks p = 0.412). Within sessions it does: POST is higher in 8 of 8 sessions (sign test p = 0.0039, one-sided rank-sum p < 0.05 in 6), and with each session's PRE mean subtracted POST is higher by 0.106 (p = 2.3e-14). The pooled test is dominated by the circular-track sessions, which score lower in every epoch (PRE 0.161, POST 0.260, against 1.112 and 1.207 on the linear tracks) and contribute 11161 of the 23801 POST events but 5912 of the 16448 PRE events. The same test run on each kind of track finds POST above PRE on the five linear-track sessions (Tukey on ranks p = 5.1e-07) and on the three circular-track sessions (p = 1.1e-05), and MAZE above PRE on both (p printed as 0, and p = 0.0007). The scored events per session (76 to 953 in MAZE, 1246 to 3294 in PRE, 1136 to 4079 in POST) and the MAZE rZ of each session are printed for comparison with the paper's figs. S2 to S4, which could not be read.
fig2b-strong-contributors¶
216 neurons showed significant PCC scores in PRE or POST and were considered to contribute strongly to sequences.
A cell's contribution to an event is |r| minus the mean |r| over 20 decodings with that cell's rate map circularly shifted, with shifts drawn afresh for every event. Strongly contributing in an epoch: one-sided Wilcoxon signed-rank test of its contributions against 0, p < 0.05, with at least 5 events. The false-positive rate comes from the same test with every rate map shifted.
On the linear tracks the criterion takes |r| on the template the event was scored with, the one of the two with the larger rZ. That |r| was selected for being the larger, so it tends to fall when any part of the template is perturbed, whatever the cell. A variant ("best template") removes this: |r| is the larger over the event's templates both with the original maps and with the cell's map shifted.
def pcc_table(rows, s, place):
df = pd.DataFrame([(e, ep, c, r0 - np.nanmean(rs)) for e, ep, c, r0, rs in rows],
columns=["event", "epoch", "cell", "contrib"])
out = []
for j in place:
rec = dict(session=s["session"], subject=s["subject"], cell=j)
for e in ("PRE", "POST"):
v = df.loc[(df.cell == j) & (df.epoch == e), "contrib"].values
rec[f"n_{e}"] = len(v)
rec[f"PCC_{e}"] = v.mean() if len(v) else np.nan
rec[f"p_{e}"] = stats.wilcoxon(v, alternative="greater").pvalue if len(v) >= 5 and np.any(v != 0) else np.nan
out.append(rec)
return pd.DataFrame(out)
def strong_cols(df):
df["strong_PRE"] = df.p_PRE < 0.05
df["strong_POST"] = df.p_POST < 0.05
df["strong"] = df.strong_PRE | df.strong_POST
df["dPCC"] = df.PCC_POST - df.PCC_PRE
return df
tables = {}
for key in ("contrib", "contrib_null", "contrib_best", "contrib_null_best"):
frames = []
for s, fp in zip(S, FP):
place = np.flatnonzero(fp["tpl"]["peak"] >= g.PLACE_PEAK)
frames.append(pcc_table(fp[key], s, place))
tables[key] = strong_cols(pd.concat(frames, ignore_index=True))
CELLS, NULL = tables["contrib"], tables["contrib_null"]
BEST, NULL_BEST = tables["contrib_best"], tables["contrib_null_best"]
# PCC from the even-numbered candidate events only, for the cross-validated controls of fig2d and fig2e
EVEN = strong_cols(pd.concat([pcc_table([r for r in fp["contrib"] if r[0] % 2 == 0], s,
np.flatnonzero(fp["tpl"]["peak"] >= g.PLACE_PEAK))
for s, fp in zip(S, FP)], ignore_index=True))
print("place cells", len(CELLS), "; with >= 5 events in PRE:", int((CELLS.n_PRE >= 5).sum()),
"in POST:", int((CELLS.n_POST >= 5).sum()))
print("strongly contributing: PRE", int(CELLS.strong_PRE.sum()), " POST", int(CELLS.strong_POST.sum()),
" either", int(CELLS.strong.sum()), f"({CELLS.strong.mean():.1%} of place cells)")
print("false-positive rate (all maps shifted): PRE", f"{NULL.strong_PRE.mean():.3f}", " POST",
f"{NULL.strong_POST.mean():.3f}", " either", f"{NULL.strong.mean():.3f}",
f"-> expected {NULL.strong.mean() * len(CELLS):.0f} cells by chance")
print("expected 216, tolerance 21.6; within:", abs(CELLS.strong.sum() - 216) <= 21.6)
for a in (0.01, 0.001):
n = int(((CELLS.p_PRE < a) | (CELLS.p_POST < a)).sum())
nn = ((NULL.p_PRE < a) | (NULL.p_POST < a)).mean()
print(f"variant alpha={a}: {n} strongly contributing, false-positive rate {nn:.3f}")
print(CELLS.groupby(["subject", "session"]).agg(place=("cell", "size"), strong=("strong", "sum"),
strong_PRE=("strong_PRE", "sum"), strong_POST=("strong_POST", "sum"))
.join(BEST.groupby(["subject", "session"]).agg(strong_best=("strong", "sum")))
.join(NULL.groupby(["subject", "session"]).agg(fpr=("strong", "mean")))
.join(NULL_BEST.groupby(["subject", "session"]).agg(fpr_best=("strong", "mean"))).round(3).to_string())
TRACK = DATA.set_index("session").track
print("false-positive rate by track (either epoch): criterion",
NULL.groupby(NULL.session.map(TRACK)).strong.mean().round(3).to_dict(),
" best template", NULL_BEST.groupby(NULL_BEST.session.map(TRACK)).strong.mean().round(3).to_dict())
print()
print("variant, best template:")
print("strongly contributing: PRE", int(BEST.strong_PRE.sum()), " POST", int(BEST.strong_POST.sum()),
" either", int(BEST.strong.sum()), f"({BEST.strong.mean():.1%} of place cells)")
print("false-positive rate (all maps shifted): PRE", f"{NULL_BEST.strong_PRE.mean():.3f}", " POST",
f"{NULL_BEST.strong_POST.mean():.3f}", " either", f"{NULL_BEST.strong.mean():.3f}",
f"-> expected {NULL_BEST.strong.mean() * len(BEST):.0f} cells by chance")
for a in (0.01, 0.001):
n = int(((BEST.p_PRE < a) | (BEST.p_POST < a)).sum())
nn = ((NULL_BEST.p_PRE < a) | (NULL_BEST.p_POST < a)).mean()
print(f"best template, alpha={a}: {n} strongly contributing, false-positive rate {nn:.3f}")
print("dPCC of the criterion against the best-template dPCC, Spearman over place cells:",
round(stats.spearmanr(CELLS.dPCC, BEST.dPCC, nan_policy="omit").statistic, 3))
/tmp/ipykernel_2198/4204492205.py:2: RuntimeWarning: Mean of empty slice df = pd.DataFrame([(e, ep, c, r0 - np.nanmean(rs)) for e, ep, c, r0, rs in rows],
/tmp/ipykernel_2198/4204492205.py:2: RuntimeWarning: Mean of empty slice df = pd.DataFrame([(e, ep, c, r0 - np.nanmean(rs)) for e, ep, c, r0, rs in rows],
/tmp/ipykernel_2198/4204492205.py:2: RuntimeWarning: Mean of empty slice df = pd.DataFrame([(e, ep, c, r0 - np.nanmean(rs)) for e, ep, c, r0, rs in rows],
place cells 417 ; with >= 5 events in PRE: 416 in POST: 417
strongly contributing: PRE 248 POST 290 either 303 (72.7% of place cells)
false-positive rate (all maps shifted): PRE 0.448 POST 0.463 either 0.561 -> expected 234 cells by chance
expected 216, tolerance 21.6; within: False
variant alpha=0.01: 277 strongly contributing, false-positive rate 0.434
variant alpha=0.001: 243 strongly contributing, false-positive rate 0.319
place strong strong_PRE strong_POST strong_best fpr fpr_best
subject session
Achilles Achilles_10252013 99 96 93 95 50 0.889 0.182
Achilles_11012013 82 43 15 39 40 0.049 0.073
Buddy Buddy_06272013 30 28 25 26 15 0.667 0.233
Cicero Cicero_09012014 31 30 27 30 18 0.871 0.323
Cicero_09102014 54 13 6 10 14 0.241 0.241
Cicero_09172014 42 41 38 39 17 0.905 0.262
Gatsby Gatsby_08022013 44 42 39 42 22 0.841 0.295
Gatsby_08282013 35 10 5 9 12 0.200 0.229
false-positive rate by track (either epoch): criterion {'1.6mLinearMaze': 0.843, '2mLinearMaze': 0.905, 'CircularMaze': 0.14} best template {'1.6mLinearMaze': 0.235, '2mLinearMaze': 0.262, 'CircularMaze': 0.158}
variant, best template:
strongly contributing: PRE 102 POST 158 either 188 (45.1% of place cells)
false-positive rate (all maps shifted): PRE 0.132 POST 0.118 either 0.206 -> expected 86 cells by chance
best template, alpha=0.01: 140 strongly contributing, false-positive rate 0.086
best template, alpha=0.001: 95 strongly contributing, false-positive rate 0.031
dPCC of the criterion against the best-template dPCC, Spearman over place cells: 0.823
Result: blocked (methods_underspecified). The test fixed in the criterion finds 303 strongly contributing cells (72.7% of 417 place cells), against the paper's 216, but with every rate map shifted it calls 56.1% of cells significant (0.843 on the 1.6 m tracks, 0.905 on the 2 m track, 0.140 on the circular track), so about 234 would be expected by chance. The excess is on the linear tracks, where the event is scored on the template with the larger rZ. The best-template variant brings the false-positive rate down to 0.206 (0.132 in PRE, 0.118 in POST; 0.158 to 0.262 by track) and finds 188 cells, outside the tolerance; at alpha 0.001 its false-positive rate is 0.031 and it finds 95 cells. Neither the criterion's test nor the variant at alpha 0.05 has the 5% false-positive rate that a significance test implies, so neither selection can be taken as the paper's. The false-positive rates are measured on every fourth PRE and POST event; a biased test becomes more significant with more events, so on all events they would likely be higher.
The best-template variant does not give identical numbers on every machine. This run prints the values above; the verifier's executed copy of the same code and seeds printed 190 cells and a false-positive rate of 0.129 in PRE, with the dependent best-template results shifting by similar amounts (cells 32 to 59). The cause was not found. No conclusion rests on these differences.
fig2-pre-post-overlap¶
Most neurons strongly contributing to PRE sequences (73%) also contributed to POST sequences.
def overlap(T, label):
pre = T[T.strong_PRE]
frac = pre.strong_POST.mean()
print(f"{label}: PRE strongly contributing {len(pre)}; also in POST {int(pre.strong_POST.sum())} ({frac:.1%}); "
f"chance (fraction of all place cells strongly contributing in POST) {T.strong_POST.mean():.1%}")
return frac
frac = overlap(CELLS, "criterion")
tab = pd.crosstab(CELLS.strong_PRE, CELLS.strong_POST)
print(tab.to_string())
print("Fisher exact test of association p =", f"{stats.fisher_exact(tab.values, alternative='greater').pvalue:.2g}")
print("by subject:", CELLS[CELLS.strong_PRE].groupby("subject").strong_POST.agg(["size", "mean"]).round(3).to_dict("index"))
print("expected 73%, tolerance 7.3; within:", abs(frac * 100 - 73) <= 7.3)
overlap(BEST, "variant, best template")
criterion: PRE strongly contributing 248; also in POST 235 (94.8%); chance (fraction of all place cells strongly contributing in POST) 69.5%
strong_POST False True
strong_PRE
False 114 55
True 13 235
Fisher exact test of association p = 2.4e-44
by subject: {'Achilles': {'size': 108, 'mean': 0.954}, 'Buddy': {'size': 25, 'mean': 0.92}, 'Cicero': {'size': 71, 'mean': 0.93}, 'Gatsby': {'size': 44, 'mean': 0.977}}
expected 73%, tolerance 7.3; within: False
variant, best template: PRE strongly contributing 102; also in POST 72 (70.6%); chance (fraction of all place cells strongly contributing in POST) 37.9%
np.float64(0.7058823529411765)
Result: blocked (methods_underspecified), as it depends on fig2b-strong-contributors. Of 248 PRE strongly contributing cells, 235 (94.8%) also contribute in POST, against a chance level of 69.5%. With the best-template variant 72 of 102 (70.6%), within 73 +/- 7.3, against a chance level of 37.9%.
Rigid and plastic cells¶
Strongly contributing cells with a PCC in both PRE and POST, split at the median dPCC of the
session's strongly contributing cells (reading (a)); reading (b), the median over all of the
session's place cells, is a variant. A second variant (cls_v) takes the strongly contributing
cells and dPCC of the best-template contributions, split as in reading (a).
CELLS["cls"] = ""
CELLS["cls_b"] = ""
CELLS["cls_v"] = ""
CELLS["dPCC_best"] = BEST.dPCC.values
CELLS["dPCC_even"] = EVEN.dPCC.values
CELLS["cls_cv"] = ""
for sess, df in CELLS.groupby("session"):
ok = df.strong & df.dPCC.notna()
med = df.loc[ok, "dPCC"].median()
CELLS.loc[ok[ok].index, "cls"] = np.where(df.loc[ok, "dPCC"] < med, "rigid", "plastic")
med_b = df.dPCC.median()
CELLS.loc[ok[ok].index, "cls_b"] = np.where(df.loc[ok, "dPCC"] < med_b, "rigid", "plastic")
medc = df.loc[ok, "dPCC_even"].median()
CELLS.loc[ok[ok].index, "cls_cv"] = np.where(df.loc[ok, "dPCC_even"] < medc, "rigid", "plastic")
okv = BEST.loc[df.index, "strong"] & df.dPCC_best.notna()
medv = df.loc[okv, "dPCC_best"].median()
CELLS.loc[okv[okv].index, "cls_v"] = np.where(df.loc[okv, "dPCC_best"] < medv, "rigid", "plastic")
print("strong cells without a PCC in one epoch (not classified):", int((CELLS.strong & CELLS.dPCC.isna()).sum()))
print(CELLS[CELLS.cls != ""].groupby(["subject", "session", "cls"]).size().unstack().to_string())
print("totals (a):", CELLS.cls.value_counts().to_dict(), " (b):", CELLS.cls_b.value_counts().to_dict(),
" best-template variant:", CELLS.cls_v.value_counts().to_dict())
print("cells classed in both (a) and the best-template variant:", int(((CELLS.cls != "") & (CELLS.cls_v != "")).sum()),
"; same class:", int(((CELLS.cls != "") & (CELLS.cls == CELLS.cls_v)).sum()))
print("dPCC from even events against dPCC from all events, Spearman:",
round(stats.spearmanr(CELLS.dPCC_even, CELLS.dPCC, nan_policy="omit").statistic, 3))
strong cells without a PCC in one epoch (not classified): 1
cls plastic rigid
subject session
Achilles Achilles_10252013 48 48
Achilles_11012013 22 21
Buddy Buddy_06272013 14 14
Cicero Cicero_09012014 15 14
Cicero_09102014 7 6
Cicero_09172014 21 20
Gatsby Gatsby_08022013 21 21
Gatsby_08282013 5 5
totals (a): {'plastic': 153, 'rigid': 149, '': 115} (b): {'plastic': 161, 'rigid': 141, '': 115} best-template variant: {'': 237, 'plastic': 92, 'rigid': 88}
cells classed in both (a) and the best-template variant: 175 ; same class: 142
dPCC from even events against dPCC from all events, Spearman: 0.695
Per-cell measures: rates, ripple gain, bursting, spatial information, fields, the in-field ratio over the session, and the 100 ms coactivation matrices.
MEAS = []
for s, fp in zip(S, FP):
place = np.flatnonzero(fp["tpl"]["peak"] >= g.PLACE_PEAK)
MEAS.append(g.cell_measures(s, fp["tpl"], place))
cols = ["rate_session", "rate_sleep", "rate_prepost", "gain", "gain_pre", "gain_post", "burst", "info",
"n_fields", "peak_pos", "best_dir", "infield_change"]
CELLS = CELLS.join(pd.concat([pd.DataFrame({c: m[c] for c in cols}) for m in MEAS], ignore_index=True))
CELLS["gain_change"] = CELLS.gain_post - CELLS.gain_pre
CELLS["coact"] = np.concatenate([np.nanmean(np.where(np.eye(len(m["C_session"]), dtype=bool), np.nan, m["C_session"]), 1)
for m in MEAS])
print(CELLS[cols + ["gain_change", "coact"]].describe().T[["count", "min", "50%", "max"]].round(3).to_string())
count min 50% max rate_session 417.0 0.101 0.653 3.683 rate_sleep 417.0 0.121 0.671 4.200 rate_prepost 417.0 0.109 0.623 4.019 gain 417.0 0.303 4.025 11.017 gain_pre 417.0 0.000 4.155 11.790 gain_post 417.0 0.422 4.183 13.968 burst 417.0 0.014 0.151 0.733 info 417.0 0.024 0.779 2.987 n_fields 417.0 0.000 2.000 6.000 peak_pos 417.0 0.003 0.445 0.997 best_dir 417.0 -1.000 -1.000 1.000 infield_change 393.0 -3.528 0.103 3.071 gain_change 417.0 -8.615 0.209 10.036 coact 417.0 -0.015 0.025 0.067
Re-decoding for the variant templates¶
Within each session: the top or bottom 10, 20, 50 and 75% of place cells by dPCC shuffled
(fig2d), plastic or rigid cells removed (fig2e), templates of the low or high-rate half (fig4d)
and of the low or high-gain half (fig4e). The rate split is made here with reading (a), a median
split within each session; the section rate-split-medians shows both readings. 50 column-cycle
shuffles per event here, in place of 200, so the base scores differ a little from those of
fig1c. The variants no_plastic_v and no_rigid_v use the best-template classes.
Controls for the selection by dPCC, which is measured on the same PRE and POST events that are
decoded again: random{X} shuffles a random X% of the place cells with a dPCC, and no_random
removes as many randomly chosen strongly contributing cells as there are plastic cells. A
cross-validated control (_cv) ranks and classes cells by the dPCC of the even-numbered
candidate events only, and its replay table uses only the odd-numbered events, which played no
part in choosing the cells. _best variants rank cells by the best-template dPCC.
Replay is judged, as the criteria fix it, on all events pooled (present); diff_within and
p_within are a variant with each session's PRE mean subtracted from its PRE and POST scores. The
same table is then given for the five linear-track sessions alone.
LEVELS = [10, 20, 50, 75]
def variants_for(i, s, fp):
df = CELLS[CELLS.session == s["session"]]
cells = fp["cells"]
rng = np.random.default_rng(9000 + i)
def restrict(keep):
keep = set(int(k) for k in keep)
return {d: np.array([c for c in cl if c in keep], int) for d, cl in cells.items()}
def drop(rem):
rem = set(int(k) for k in rem)
return {d: np.array([c for c in cl if c not in rem], int) for d, cl in cells.items()}
v = {"base": (cells, None, 5000 + i)}
ranked = df[df.dPCC.notna()].sort_values("dPCC")
n = len(ranked)
for k, X in enumerate(LEVELS):
m = int(round(n * X / 100))
v[f"top{X}"] = (cells, ranked.cell.values[n - m:], 5100 + 10 * k + i)
v[f"bottom{X}"] = (cells, ranked.cell.values[:m], 5200 + 10 * k + i)
for k, X in enumerate((10, 20, 75)):
m = int(round(n * X / 100))
v[f"random{X}"] = (cells, rng.choice(ranked.cell.values, m, replace=False), 5300 + 10 * k + i)
ranked_e = df[df.dPCC_even.notna()].sort_values("dPCC_even")
ranked_b = df[df.dPCC_best.notna()].sort_values("dPCC_best")
for k, (name, X) in enumerate([("top", 10), ("top", 20), ("bottom", 75)]):
for tag, rk, base_seed in (("cv", ranked_e, 5400), ("best", ranked_b, 5500)):
m = int(round(len(rk) * X / 100))
sel = rk.cell.values[len(rk) - m:] if name == "top" else rk.cell.values[:m]
v[f"{name}{X}_{tag}"] = (cells, sel, base_seed + 10 * k + i)
v["no_plastic_cv"] = (drop(df.cell[df.cls_cv == "plastic"]), None, 6500 + i)
v["no_rigid_cv"] = (drop(df.cell[df.cls_cv == "rigid"]), None, 6600 + i)
v["no_plastic"] = (drop(df.cell[df.cls == "plastic"]), None, 6000 + i)
v["no_rigid"] = (drop(df.cell[df.cls == "rigid"]), None, 6100 + i)
classed = df.cell[df.cls != ""].values
v["no_random"] = (drop(rng.choice(classed, int((df.cls == "plastic").sum()), replace=False)), None, 6400 + i)
v["no_plastic_v"] = (drop(df.cell[df.cls_v == "plastic"]), None, 6200 + i)
v["no_rigid_v"] = (drop(df.cell[df.cls_v == "rigid"]), None, 6300 + i)
rs = df.rate_sleep.rank(method="first").values
low = rs <= len(df) // 2
v["low_rate"] = (restrict(df.cell[low]), None, 7000 + i)
v["high_rate"] = (restrict(df.cell[~low]), None, 7100 + i)
gr = df.gain.rank(method="first").values
hi = gr > len(df) // 2
v["high_gain"] = (restrict(df.cell[hi]), None, 8000 + i)
v["low_gain"] = (restrict(df.cell[~hi]), None, 8100 + i)
return v
VARS = [variants_for(i, s, fp) for i, (s, fp) in enumerate(zip(S, FP))]
T0 = time.time()
with ProcessPoolExecutor(4) as ex:
VRES = list(ex.map(g.variant_pass, S, [fp["tpl"] for fp in FP], VARS))
print(f"variants in {time.time() - T0:.0f} s")
LINEAR = np.array([not s["track"].startswith("Circular") for s in S])
def replay(name, keep=None, odd=False):
R = [r for r, k in zip(VRES, keep if keep is not None else [True] * len(VRES)) if k]
if odd: # only the odd-numbered candidate events
R = [{name: {q: w[r[name]["index"] % 2 == 1] for q, w in r[name].items()}} for r in R]
pre = np.concatenate([r[name]["rZ"][r[name]["epoch"] == "PRE"] for r in R])
post = np.concatenate([r[name]["rZ"][r[name]["epoch"] == "POST"] for r in R])
p = stats.mannwhitneyu(post, pre, alternative="greater").pvalue
sess = sum(np.nanmean(r[name]["rZ"][r[name]["epoch"] == "POST"]) > np.nanmean(r[name]["rZ"][r[name]["epoch"] == "PRE"])
for r in R)
# variant: each session's PRE mean subtracted from its PRE and POST scores
pre_c = np.concatenate([r[name]["rZ"][r[name]["epoch"] == "PRE"] - r[name]["rZ"][r[name]["epoch"] == "PRE"].mean() for r in R])
post_c = np.concatenate([r[name]["rZ"][r[name]["epoch"] == "POST"] - r[name]["rZ"][r[name]["epoch"] == "PRE"].mean() for r in R])
pc = stats.mannwhitneyu(post_c, pre_c, alternative="greater").pvalue
return dict(variant=name, n_PRE=len(pre), n_POST=len(post), mean_PRE=pre.mean(), mean_POST=post.mean(),
diff=post.mean() - pre.mean(), p_one_sided=p, present=p < 0.05, sessions_POST_gt_PRE=sess,
diff_within=post_c.mean(), p_within=pc)
REPLAY = pd.DataFrame([replay(k) for k in VARS[0]]).set_index("variant")
REPLAY_LIN = pd.DataFrame([replay(k, LINEAR) for k in VARS[0]]).set_index("variant")
CV_NAMES = ["base", "top10_cv", "top20_cv", "bottom75_cv", "no_plastic_cv", "no_rigid_cv"]
REPLAY_CV = pd.DataFrame([replay(k, odd=True) for k in CV_NAMES]).set_index("variant")
REPLAY_CV_LIN = pd.DataFrame([replay(k, LINEAR, odd=True) for k in CV_NAMES]).set_index("variant")
print("all 8 sessions:")
print(REPLAY.round(4).to_string())
print()
print(f"linear-track sessions only ({int(LINEAR.sum())}):")
print(REPLAY_LIN.round(4).to_string())
variants in 136 s
all 8 sessions:
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 16448 23801 0.7951 0.7868 -0.0083 0.9160 False 8 0.1067 0.0000
top10 16448 23801 0.7957 0.6956 -0.1001 1.0000 False 4 0.0113 0.3377
bottom10 16448 23801 0.7221 0.7566 0.0345 0.0412 True 8 0.1436 0.0000
top20 16448 23801 0.7748 0.6477 -0.1271 1.0000 False 4 -0.0136 0.9655
bottom20 16448 23801 0.7224 0.7392 0.0167 0.2128 False 8 0.1243 0.0000
top50 16448 23801 0.7094 0.5497 -0.1597 1.0000 False 1 -0.0571 1.0000
bottom50 16448 23801 0.6313 0.6705 0.0393 0.0061 True 7 0.1276 0.0000
top75 16448 23801 0.6511 0.4427 -0.2084 1.0000 False 1 -0.1006 1.0000
bottom75 16448 23801 0.5979 0.6197 0.0218 0.2085 False 8 0.1024 0.0000
random10 16448 23801 0.7615 0.7640 0.0025 0.6325 False 8 0.1211 0.0000
random20 16448 23801 0.7511 0.7298 -0.0213 0.9933 False 8 0.0810 0.0000
random75 16448 23801 0.5956 0.5224 -0.0732 1.0000 False 3 0.0121 0.1715
top10_cv 16448 23801 0.7879 0.7185 -0.0694 1.0000 False 6 0.0357 0.0183
top10_best 16448 23801 0.7821 0.6473 -0.1348 1.0000 False 3 -0.0149 0.9824
top20_cv 16448 23801 0.7573 0.6444 -0.1129 1.0000 False 5 -0.0081 0.8898
top20_best 16448 23801 0.7811 0.6325 -0.1486 1.0000 False 1 -0.0346 0.9994
bottom75_cv 16448 23801 0.6045 0.5974 -0.0071 0.8751 False 7 0.0690 0.0000
bottom75_best 16448 23801 0.6436 0.6567 0.0131 0.3711 False 7 0.0963 0.0000
no_plastic_cv 14191 20401 0.7140 0.5741 -0.1398 1.0000 False 3 -0.0144 0.9091
no_rigid_cv 13867 21509 0.6882 0.6906 0.0024 0.6099 False 7 0.1206 0.0000
no_plastic 14432 20805 0.7594 0.5911 -0.1683 1.0000 False 3 -0.0466 0.9999
no_rigid 13449 21035 0.6468 0.6850 0.0382 0.0194 True 8 0.1526 0.0000
no_random 13702 20514 0.6419 0.6195 -0.0225 0.9889 False 7 0.0819 0.0000
no_plastic_v 15260 22235 0.6964 0.5717 -0.1247 1.0000 False 5 -0.0045 0.7121
no_rigid_v 15537 22832 0.6935 0.7135 0.0200 0.2390 False 8 0.1402 0.0000
low_rate 8212 12549 0.6998 0.6763 -0.0235 0.9641 False 6 0.1157 0.0000
high_rate 14812 21567 0.8081 0.7594 -0.0488 0.9999 False 6 0.0583 0.0001
high_gain 12342 18081 0.7231 0.7222 -0.0009 0.7097 False 5 0.1114 0.0000
low_gain 12086 17653 0.8412 0.7304 -0.1108 1.0000 False 3 0.0155 0.1416
linear-track sessions only (5):
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 10536 12640 1.1443 1.2452 0.1008 0.0000 True 5 0.1261 0.0000
top10 10536 12640 1.1300 1.1448 0.0148 0.2289 False 3 0.0383 0.0177
bottom10 10536 12640 1.0496 1.1964 0.1468 0.0000 True 5 0.1723 0.0000
top20 10536 12640 1.1051 1.0666 -0.0385 0.9951 False 2 -0.0089 0.8169
bottom20 10536 12640 1.0499 1.1506 0.1007 0.0000 True 5 0.1194 0.0000
top50 10536 12640 1.0258 0.9478 -0.0780 1.0000 False 1 -0.0716 1.0000
bottom50 10536 12640 0.9163 1.0715 0.1552 0.0000 True 5 0.1717 0.0000
top75 10536 12640 0.9570 0.8507 -0.1063 1.0000 False 1 -0.0843 1.0000
bottom75 10536 12640 0.8755 0.9732 0.0976 0.0000 True 5 0.0936 0.0000
random10 10536 12640 1.1272 1.2267 0.0994 0.0000 True 5 0.1262 0.0000
random20 10536 12640 1.0748 1.1631 0.0883 0.0000 True 5 0.1042 0.0000
random75 10536 12640 0.8732 0.8605 -0.0126 0.8006 False 1 -0.0197 0.8887
top10_cv 10536 12640 1.1088 1.1232 0.0143 0.3044 False 4 0.0384 0.0265
top10_best 10536 12640 1.1218 1.0766 -0.0452 0.9992 False 2 -0.0084 0.8362
top20_cv 10536 12640 1.0888 1.0776 -0.0113 0.8285 False 3 0.0085 0.4225
top20_best 10536 12640 1.1077 1.0408 -0.0669 1.0000 False 1 -0.0346 0.9928
bottom75_cv 10536 12640 0.8529 0.9488 0.0958 0.0000 True 5 0.1066 0.0000
bottom75_best 10536 12640 0.8750 0.9844 0.1095 0.0000 True 5 0.1328 0.0000
no_plastic_cv 8471 9510 1.1178 1.0800 -0.0378 0.9670 False 2 -0.0378 0.9679
no_rigid_cv 8180 10710 1.0947 1.1957 0.1010 0.0000 True 5 0.1684 0.0000
no_plastic 8715 9958 1.1527 1.0932 -0.0595 0.9997 False 2 -0.0609 0.9997
no_rigid 7761 10209 1.0870 1.2232 0.1362 0.0000 True 5 0.1927 0.0000
no_random 7982 9690 1.0228 1.1170 0.0942 0.0000 True 5 0.1121 0.0000
no_plastic_v 9565 11430 1.0258 0.9891 -0.0367 0.9866 False 4 0.0110 0.2750
no_rigid_v 9960 12117 1.0513 1.1892 0.1379 0.0000 True 5 0.1696 0.0000
low_rate 5308 6268 0.9704 1.0476 0.0772 0.0046 True 4 0.1224 0.0000
high_rate 9517 11652 1.1767 1.2201 0.0435 0.0175 True 3 0.0661 0.0004
high_gain 7641 9040 1.0731 1.0965 0.0233 0.2495 False 2 0.0453 0.0288
low_gain 8169 10026 1.1554 1.1655 0.0101 0.2083 False 3 0.0497 0.0028
fig2d-shuffle-high-dpcc¶
The PRE to POST increase in sequence score is eliminated when the place fields of the 10 to 20% of cells with the highest dPCC are shuffled, and remains after shuffling the bottom 75%.
COLS = ["n_PRE", "n_POST", "mean_PRE", "mean_POST", "diff", "p_one_sided", "present", "sessions_POST_gt_PRE", "diff_within", "p_within"]
rows_2d = ["base"] + [f"top{X}" for X in LEVELS] + [f"bottom{X}" for X in LEVELS] + ["random10", "random20", "random75",
"top10_best", "top20_best", "bottom75_best"]
print(REPLAY.loc[rows_2d, COLS].round(4).to_string())
elim = (not REPLAY.loc["top10", "present"]) or (not REPLAY.loc["top20", "present"])
remain = REPLAY.loc["bottom75", "present"]
print("replay absent after shuffling top 10% or 20%:", elim, "; present after shuffling bottom 75%:", remain)
print("linear-track sessions only:")
print(REPLAY_LIN.loc[rows_2d, COLS].round(4).to_string())
print("cross-validated: cells ranked on even events, replay on odd events, all sessions:")
print(REPLAY_CV.loc[["base", "top10_cv", "top20_cv", "bottom75_cv"], COLS].round(4).to_string())
print("cross-validated, linear-track sessions only:")
print(REPLAY_CV_LIN.loc[["base", "top10_cv", "top20_cv", "bottom75_cv"], COLS].round(4).to_string())
def by_subject(name):
out = {}
for s, r in zip(S, VRES):
d = out.setdefault(s["subject"], {"PRE": [], "POST": []})
for e in ("PRE", "POST"):
d[e].append(r[name]["rZ"][r[name]["epoch"] == e])
return {k: round(np.concatenate(v["POST"]).mean() - np.concatenate(v["PRE"]).mean(), 3) for k, v in out.items()}
for name in ("base", "top10", "top20", "bottom75", "random10", "random20"):
print(name, "POST minus PRE mean rZ by subject:", by_subject(name))
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 16448 23801 0.7951 0.7868 -0.0083 0.9160 False 8 0.1067 0.0000
top10 16448 23801 0.7957 0.6956 -0.1001 1.0000 False 4 0.0113 0.3377
top20 16448 23801 0.7748 0.6477 -0.1271 1.0000 False 4 -0.0136 0.9655
top50 16448 23801 0.7094 0.5497 -0.1597 1.0000 False 1 -0.0571 1.0000
top75 16448 23801 0.6511 0.4427 -0.2084 1.0000 False 1 -0.1006 1.0000
bottom10 16448 23801 0.7221 0.7566 0.0345 0.0412 True 8 0.1436 0.0000
bottom20 16448 23801 0.7224 0.7392 0.0167 0.2128 False 8 0.1243 0.0000
bottom50 16448 23801 0.6313 0.6705 0.0393 0.0061 True 7 0.1276 0.0000
bottom75 16448 23801 0.5979 0.6197 0.0218 0.2085 False 8 0.1024 0.0000
random10 16448 23801 0.7615 0.7640 0.0025 0.6325 False 8 0.1211 0.0000
random20 16448 23801 0.7511 0.7298 -0.0213 0.9933 False 8 0.0810 0.0000
random75 16448 23801 0.5956 0.5224 -0.0732 1.0000 False 3 0.0121 0.1715
top10_best 16448 23801 0.7821 0.6473 -0.1348 1.0000 False 3 -0.0149 0.9824
top20_best 16448 23801 0.7811 0.6325 -0.1486 1.0000 False 1 -0.0346 0.9994
bottom75_best 16448 23801 0.6436 0.6567 0.0131 0.3711 False 7 0.0963 0.0000
replay absent after shuffling top 10% or 20%: True ; present after shuffling bottom 75%: False
linear-track sessions only:
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 10536 12640 1.1443 1.2452 0.1008 0.0000 True 5 0.1261 0.0000
top10 10536 12640 1.1300 1.1448 0.0148 0.2289 False 3 0.0383 0.0177
top20 10536 12640 1.1051 1.0666 -0.0385 0.9951 False 2 -0.0089 0.8169
top50 10536 12640 1.0258 0.9478 -0.0780 1.0000 False 1 -0.0716 1.0000
top75 10536 12640 0.9570 0.8507 -0.1063 1.0000 False 1 -0.0843 1.0000
bottom10 10536 12640 1.0496 1.1964 0.1468 0.0000 True 5 0.1723 0.0000
bottom20 10536 12640 1.0499 1.1506 0.1007 0.0000 True 5 0.1194 0.0000
bottom50 10536 12640 0.9163 1.0715 0.1552 0.0000 True 5 0.1717 0.0000
bottom75 10536 12640 0.8755 0.9732 0.0976 0.0000 True 5 0.0936 0.0000
random10 10536 12640 1.1272 1.2267 0.0994 0.0000 True 5 0.1262 0.0000
random20 10536 12640 1.0748 1.1631 0.0883 0.0000 True 5 0.1042 0.0000
random75 10536 12640 0.8732 0.8605 -0.0126 0.8006 False 1 -0.0197 0.8887
top10_best 10536 12640 1.1218 1.0766 -0.0452 0.9992 False 2 -0.0084 0.8362
top20_best 10536 12640 1.1077 1.0408 -0.0669 1.0000 False 1 -0.0346 0.9928
bottom75_best 10536 12640 0.8750 0.9844 0.1095 0.0000 True 5 0.1328 0.0000
cross-validated: cells ranked on even events, replay on odd events, all sessions:
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 8228 11905 0.7924 0.7875 -0.0048 0.7530 False 7 0.1072 0.0000
top10_cv 8228 11905 0.7804 0.7315 -0.0490 0.9974 False 5 0.0569 0.0036
top20_cv 8228 11905 0.7420 0.6705 -0.0715 1.0000 False 6 0.0332 0.0597
bottom75_cv 8228 11905 0.6146 0.5841 -0.0305 0.9786 False 4 0.0426 0.0095
cross-validated, linear-track sessions only:
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 5262 6310 1.1366 1.2525 0.1158 0.0000 True 5 0.1404 0.0000
top10_cv 5262 6310 1.1021 1.1484 0.0463 0.0238 True 4 0.0731 0.0011
top20_cv 5262 6310 1.0772 1.1188 0.0415 0.0378 True 4 0.0644 0.0035
bottom75_cv 5262 6310 0.8660 0.9470 0.0810 0.0002 True 3 0.0881 0.0000
base POST minus PRE mean rZ by subject: {'Achilles': np.float64(0.085), 'Buddy': np.float64(0.178), 'Cicero': np.float64(0.02), 'Gatsby': np.float64(-0.17)}
top10 POST minus PRE mean rZ by subject: {'Achilles': np.float64(0.006), 'Buddy': np.float64(0.124), 'Cicero': np.float64(-0.114), 'Gatsby': np.float64(-0.226)}
top20 POST minus PRE mean rZ by subject: {'Achilles': np.float64(-0.073), 'Buddy': np.float64(0.103), 'Cicero': np.float64(-0.107), 'Gatsby': np.float64(-0.237)}
bottom75 POST minus PRE mean rZ by subject: {'Achilles': np.float64(0.016), 'Buddy': np.float64(0.057), 'Cicero': np.float64(0.102), 'Gatsby': np.float64(-0.102)}
random10 POST minus PRE mean rZ by subject: {'Achilles': np.float64(0.078), 'Buddy': np.float64(0.172), 'Cicero': np.float64(0.029), 'Gatsby': np.float64(-0.146)}
random20 POST minus PRE mean rZ by subject: {'Achilles': np.float64(0.042), 'Buddy': np.float64(0.07), 'Cicero': np.float64(0.022), 'Gatsby': np.float64(-0.166)}
Result: partially reproduced. The criterion ranks all place cells by dPCC and does not use the selection of fig2b-strong-contributors, so the claim is tested. On the pooled test of the criterion replay is absent after shuffling the top 10% (POST minus PRE -0.100, p = 1.0) and the top 20% (-0.127, p = 1.0), but it is also absent after shuffling the bottom 75% (+0.022, p = 0.21), so one of the two conditions holds. The pooled test finds no replay before any shuffling either (-0.008, p = 0.916), so on all sessions it cannot show a removal.
On the five linear-track sessions alone, replay is present before shuffling (+0.101, p < 0.0001), absent after shuffling the top 10% (+0.015, p = 0.23) or 20% (-0.039, p = 0.995), and present after shuffling the bottom 75% (+0.098, p < 0.0001), the paper's pattern. Shuffling a random 10% or 20% of the cells leaves it (+0.099 and +0.088, p < 0.0001); ranking by the best-template dPCC gives the same pattern. dPCC is measured on the same events that are decoded again. When cells are ranked on the even-numbered events and replay is measured on the odd-numbered ones, replay on the linear tracks falls from +0.116 to +0.046 (p = 0.024) after shuffling the top 10% and to +0.042 (p = 0.038) after the top 20%, and remains after shuffling the bottom 75% (+0.081, p = 0.0002). Out of sample, shuffling the high-dPCC cells reduces replay but does not remove it.
fig2e-remove-plastic¶
The PRE to POST increase in sequence content is attributable to a small number of plastic cells; removing them, unlike removing rigid cells, removes the increase.
rows_2e = ["base", "no_plastic", "no_rigid", "no_random", "no_plastic_v", "no_rigid_v"]
print(REPLAY.loc[rows_2e, COLS].round(4).to_string())
print("linear-track sessions only:")
print(REPLAY_LIN.loc[rows_2e, COLS].round(4).to_string())
print("cross-validated: classes from even events, replay on odd events, all sessions:")
print(REPLAY_CV.loc[["base", "no_plastic_cv", "no_rigid_cv"], COLS].round(4).to_string())
print("cross-validated, linear-track sessions only:")
print(REPLAY_CV_LIN.loc[["base", "no_plastic_cv", "no_rigid_cv"], COLS].round(4).to_string())
for name in ("no_plastic", "no_rigid", "no_random"):
print(name, "POST minus PRE mean rZ by subject:", by_subject(name))
print("ratio of differences (no plastic / no rigid):",
round(REPLAY.loc["no_plastic", "diff"] / REPLAY.loc["no_rigid", "diff"], 3))
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 16448 23801 0.7951 0.7868 -0.0083 0.9160 False 8 0.1067 0.0000
no_plastic 14432 20805 0.7594 0.5911 -0.1683 1.0000 False 3 -0.0466 0.9999
no_rigid 13449 21035 0.6468 0.6850 0.0382 0.0194 True 8 0.1526 0.0000
no_random 13702 20514 0.6419 0.6195 -0.0225 0.9889 False 7 0.0819 0.0000
no_plastic_v 15260 22235 0.6964 0.5717 -0.1247 1.0000 False 5 -0.0045 0.7121
no_rigid_v 15537 22832 0.6935 0.7135 0.0200 0.2390 False 8 0.1402 0.0000
linear-track sessions only:
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 10536 12640 1.1443 1.2452 0.1008 0.0000 True 5 0.1261 0.0000
no_plastic 8715 9958 1.1527 1.0932 -0.0595 0.9997 False 2 -0.0609 0.9997
no_rigid 7761 10209 1.0870 1.2232 0.1362 0.0000 True 5 0.1927 0.0000
no_random 7982 9690 1.0228 1.1170 0.0942 0.0000 True 5 0.1121 0.0000
no_plastic_v 9565 11430 1.0258 0.9891 -0.0367 0.9866 False 4 0.0110 0.2750
no_rigid_v 9960 12117 1.0513 1.1892 0.1379 0.0000 True 5 0.1696 0.0000
cross-validated: classes from even events, replay on odd events, all sessions:
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 8228 11905 0.7924 0.7875 -0.0048 0.7530 False 7 0.1072 0.0000
no_plastic_cv 7093 10256 0.7133 0.5894 -0.1239 1.0000 False 4 -0.0032 0.6855
no_rigid_cv 6924 10772 0.7024 0.6807 -0.0218 0.9333 False 6 0.0836 0.0001
cross-validated, linear-track sessions only:
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 5262 6310 1.1366 1.2525 0.1158 0.0000 True 5 0.1404 0.0000
no_plastic_cv 4224 4793 1.1138 1.1047 -0.0091 0.5426 False 3 -0.0092 0.5363
no_rigid_cv 4081 5367 1.0984 1.1880 0.0896 0.0027 True 5 0.1428 0.0000
no_plastic POST minus PRE mean rZ by subject: {'Achilles': np.float64(-0.135), 'Buddy': np.float64(0.022), 'Cicero': np.float64(-0.116), 'Gatsby': np.float64(-0.279)}
no_rigid POST minus PRE mean rZ by subject: {'Achilles': np.float64(0.065), 'Buddy': np.float64(0.299), 'Cicero': np.float64(0.149), 'Gatsby': np.float64(-0.17)}
no_random POST minus PRE mean rZ by subject: {'Achilles': np.float64(0.029), 'Buddy': np.float64(0.013), 'Cicero': np.float64(0.032), 'Gatsby': np.float64(-0.166)}
ratio of differences (no plastic / no rigid): -4.403
Result: blocked (methods_underspecified), as it depends on fig2b-strong-contributors. On the pooled test, removing plastic cells gives POST minus PRE -0.168 (p = 1.0) and removing rigid cells +0.038 (p = 0.019), the pattern the criterion asks for; removing as many randomly chosen strongly contributing cells gives -0.023 (p = 0.99). On the linear-track sessions alone the three give -0.060 (p = 0.9997), +0.136 and +0.094 (both p < 0.0001). With classes from the even-numbered events and replay on the odd-numbered ones, the linear tracks give -0.009 (p = 0.54) without plastic cells and +0.090 (p = 0.0027) without rigid cells. With the best-template classes, pooled over all sessions, -0.125 (p = 1.0) and +0.020 (p = 0.24).
figs2-track-distribution¶
Rigid and plastic neurons were similarly distributed along the track.
R = CELLS[CELLS.cls == "rigid"]; P = CELLS[CELLS.cls == "plastic"]
ks = stats.ks_2samp(R.peak_pos, P.peak_pos)
print(f"peak position (fraction of track): rigid n={len(R)} median {R.peak_pos.median():.3f}; "
f"plastic n={len(P)} median {P.peak_pos.median():.3f}; KS D={ks.statistic:.3f} p={ks.pvalue:.3g}")
Rb = CELLS[CELLS.cls_b == "rigid"]; Pb = CELLS[CELLS.cls_b == "plastic"]
print(f"reading (b): KS p={stats.ks_2samp(Rb.peak_pos, Pb.peak_pos).pvalue:.3g}")
Rv = CELLS[CELLS.cls_v == "rigid"]; Pv = CELLS[CELLS.cls_v == "plastic"]
print(f"best-template variant: rigid n={len(Rv)} plastic n={len(Pv)}, KS D={stats.ks_2samp(Rv.peak_pos, Pv.peak_pos).statistic:.3f} "
f"p={stats.ks_2samp(Rv.peak_pos, Pv.peak_pos).pvalue:.3g}")
print("by subject (median peak position):")
print(CELLS[CELLS.cls != ""].groupby(["subject", "cls"]).peak_pos.median().unstack().round(3).to_string())
peak position (fraction of track): rigid n=149 median 0.419; plastic n=153 median 0.486; KS D=0.145 p=0.0723 reading (b): KS p=0.0787 best-template variant: rigid n=88 plastic n=92, KS D=0.173 p=0.117 by subject (median peak position): cls plastic rigid subject Achilles 0.457 0.456 Buddy 0.319 0.269 Cicero 0.494 0.461 Gatsby 0.506 0.357
Result: blocked (methods_underspecified), as it depends on fig2b-strong-contributors. The peak positions of rigid and plastic cells do not differ by the KS test fixed in the criterion (medians 0.419 and 0.486 of the track, D = 0.145, p = 0.072); with reading (b) p = 0.079 and with the best-template classes p = 0.12. The first two are close to the threshold.
fig3a-plastic-specificity¶
Plastic neurons had higher place-specific indices and fewer place fields than rigid cells.
def compare(col, a="plastic", b="rigid", cls="cls", alt="greater"):
A = CELLS.loc[CELLS[cls] == a, col].dropna(); B = CELLS.loc[CELLS[cls] == b, col].dropna()
p = stats.mannwhitneyu(A, B).pvalue
print(f"{col}: {a} n={len(A)} median {A.median():.3f} mean {A.mean():.3f}; {b} n={len(B)} median "
f"{B.median():.3f} mean {B.mean():.3f}; rank-sum p={p:.3g}")
return A, B, p
_, _, p_info = compare("info")
_, _, p_nf = compare("n_fields")
print("reading (b):"); compare("info", cls="cls_b"); compare("n_fields", cls="cls_b")
print("best-template variant:"); compare("info", cls="cls_v"); compare("n_fields", cls="cls_v")
print(CELLS[CELLS.cls != ""].groupby(["subject", "cls"])[["info", "n_fields"]].mean().unstack().round(3).to_string())
info: plastic n=153 median 0.904 mean 1.026; rigid n=149 median 0.594 mean 0.713; rank-sum p=1.89e-05
n_fields: plastic n=153 median 1.000 mean 1.484; rigid n=149 median 1.000 mean 1.685; rank-sum p=0.0732
reading (b):
info: plastic n=161 median 0.900 mean 1.014; rigid n=141 median 0.593 mean 0.709; rank-sum p=2.33e-05
n_fields: plastic n=161 median 1.000 mean 1.547; rigid n=141 median 1.000 mean 1.624; rank-sum p=0.506
best-template variant:
info: plastic n=92 median 0.976 mean 1.101; rigid n=88 median 0.760 mean 0.812; rank-sum p=0.00353
n_fields: plastic n=92 median 1.000 mean 1.576; rigid n=88 median 2.000 mean 1.784; rank-sum p=0.104
info n_fields
cls plastic rigid plastic rigid
subject
Achilles 1.161 0.831 1.471 1.783
Buddy 1.260 0.829 1.357 1.429
Cicero 0.796 0.507 1.535 1.600
Gatsby 0.917 0.655 1.500 1.692
Result: blocked (methods_underspecified), as it depends on fig2b-strong-contributors. Plastic cells carry more spatial information than rigid cells (median 0.904 against 0.594 bits per spike, p = 1.9e-05; best-template classes p = 0.0035), but their number of fields does not differ (means 1.484 and 1.685, p = 0.073; best-template classes p = 0.10).
fig3b-plastic-improve¶
Plastic but not rigid neurons increased their spatial specificity (in-field versus out-of-field firing ratio) steadily during the MAZE session.
for c in ("plastic", "rigid"):
v = CELLS.loc[CELLS.cls == c, "infield_change"]
print(f"{c}: n={v.notna().sum()} (no field: {v.isna().sum()}), median change of log10 ratio "
f"{v.median():.3f}, mean {v.mean():.3f}, cells increasing {(v > 0).sum()}; one-sided sign-rank p="
f"{stats.wilcoxon(v.dropna(), alternative='greater').pvalue:.3g}")
A = CELLS.loc[CELLS.cls == "plastic", "infield_change"].dropna(); B = CELLS.loc[CELLS.cls == "rigid", "infield_change"].dropna()
print(f"plastic vs rigid rank-sum p={stats.mannwhitneyu(A, B).pvalue:.3g}")
for c in ("plastic", "rigid"):
v = CELLS.loc[CELLS.cls_b == c, "infield_change"].dropna()
print(f"reading (b) {c}: n={len(v)} median {v.median():.3f} p={stats.wilcoxon(v, alternative='greater').pvalue:.3g}")
for c in ("plastic", "rigid"):
v = CELLS.loc[CELLS.cls_v == c, "infield_change"].dropna()
print(f"best-template variant {c}: n={len(v)} median {v.median():.3f} p={stats.wilcoxon(v, alternative='greater').pvalue:.3g}")
print(CELLS[CELLS.cls != ""].groupby(["subject", "cls"]).infield_change.median().unstack().round(3).to_string())
plastic: n=145 (no field: 8), median change of log10 ratio 0.155, mean 0.155, cells increasing 86; one-sided sign-rank p=0.00791 rigid: n=138 (no field: 11), median change of log10 ratio 0.129, mean 0.105, cells increasing 80; one-sided sign-rank p=0.0301 plastic vs rigid rank-sum p=0.405 reading (b) plastic: n=153 median 0.146 p=0.0115 reading (b) rigid: n=130 median 0.136 p=0.0204 best-template variant plastic: n=86 median 0.154 p=0.0197 best-template variant rigid: n=87 median 0.047 p=0.278 cls plastic rigid subject Achilles 0.074 0.035 Buddy 0.202 0.026 Cicero 0.220 0.152 Gatsby 0.380 0.213
Result: blocked (methods_underspecified), as it depends on fig2b-strong-contributors. Both groups increase their in-field ratio over the session (change of log10 ratio, plastic median 0.155, p = 0.0079; rigid 0.129, p = 0.030; between groups p = 0.41). With the best-template classes plastic cells increase (p = 0.020) and rigid cells do not (p = 0.28).
fig4a-rigid-faster¶
Session-wide firing rates of rigid neurons were higher than those of plastic neurons.
compare("rate_session", a="rigid", b="plastic")
print("reading (b):"); compare("rate_session", a="rigid", b="plastic", cls="cls_b")
print("best-template variant:"); compare("rate_session", a="rigid", b="plastic", cls="cls_v")
print(CELLS[CELLS.cls != ""].groupby(["subject", "cls"]).rate_session.median().unstack().round(3).to_string())
rate_session: rigid n=149 median 0.700 mean 0.887; plastic n=153 median 0.591 mean 0.785; rank-sum p=0.0616 reading (b): rate_session: rigid n=141 median 0.700 mean 0.898; plastic n=161 median 0.591 mean 0.780; rank-sum p=0.0532 best-template variant: rate_session: rigid n=88 median 0.664 mean 0.849; plastic n=92 median 0.626 mean 0.833; rank-sum p=0.743 cls plastic rigid subject Achilles 0.518 0.626 Buddy 0.471 0.415 Cicero 0.828 0.887 Gatsby 0.638 0.773
Result: blocked (methods_underspecified), as it depends on fig2b-strong-contributors. Rigid cells fire somewhat faster (median 0.700 against 0.591 Hz) without significance (p = 0.062; reading (b) p = 0.053; best-template classes p = 0.74).
fig4a-plastic-ripple-burst-coactivation¶
Ripple-related spiking, bursting and pairwise coactivation were higher in plastic than in rigid neurons.
res = {}
for col in ("gain", "burst", "coact"):
A, B, p = compare(col)
res[col] = (A.median() > B.median()) and p < 0.05
print("higher in plastic with p < 0.05:", res, "->", sum(res.values()), "of 3")
print("reading (b):")
for col in ("gain", "burst", "coact"):
compare(col, cls="cls_b")
print("best-template variant:")
for col in ("gain", "burst", "coact"):
compare(col, cls="cls_v")
print(CELLS[CELLS.cls != ""].groupby(["subject", "cls"])[["gain", "burst", "coact"]].median().unstack().round(3).to_string())
gain: plastic n=153 median 4.489 mean 4.748; rigid n=149 median 3.939 mean 4.131; rank-sum p=0.0152
burst: plastic n=153 median 0.180 mean 0.194; rigid n=149 median 0.147 mean 0.171; rank-sum p=0.0282
coact: plastic n=153 median 0.028 mean 0.027; rigid n=149 median 0.027 mean 0.027; rank-sum p=0.9
higher in plastic with p < 0.05: {'gain': np.True_, 'burst': np.True_, 'coact': np.False_} -> 2 of 3
reading (b):
gain: plastic n=161 median 4.455 mean 4.724; rigid n=141 median 3.923 mean 4.123; rank-sum p=0.0185
burst: plastic n=161 median 0.180 mean 0.193; rigid n=141 median 0.147 mean 0.171; rank-sum p=0.0262
coact: plastic n=161 median 0.027 mean 0.027; rigid n=141 median 0.027 mean 0.027; rank-sum p=0.964
best-template variant:
gain: plastic n=92 median 4.512 mean 4.767; rigid n=88 median 4.406 mean 4.258; rank-sum p=0.152
burst: plastic n=92 median 0.188 mean 0.204; rigid n=88 median 0.163 mean 0.179; rank-sum p=0.112
coact: plastic n=92 median 0.030 mean 0.027; rigid n=88 median 0.027 mean 0.026; rank-sum p=0.412
gain burst coact
cls plastic rigid plastic rigid plastic rigid
subject
Achilles 5.112 4.167 0.248 0.212 0.027 0.027
Buddy 4.553 4.077 0.138 0.133 0.029 0.028
Cicero 4.154 3.794 0.102 0.076 0.032 0.030
Gatsby 3.732 3.403 0.147 0.141 0.022 0.021
Result: blocked (methods_underspecified), as it depends on fig2b-strong-contributors. Ripple gain (4.489 against 3.939, p = 0.015) and bursting (0.180 against 0.147, p = 0.028) are higher in plastic cells; coactivation is not (0.028 against 0.027, p = 0.90). With the best-template classes none reaches p < 0.05 (gain p = 0.15, bursting p = 0.11, coactivation p = 0.41).
fig4c: coactivation by firing rate¶
Place cells of each session sorted by session firing rate; windows of 20% of the cells stepped by 1%. For each window: mean pairwise correlation (100 ms bins) in PRE, the mean POST minus PRE change, and the partial correlation across its pairs of MAZE running and POST coactivation given PRE. Curves are averaged over sessions by window position.
POS = np.arange(0, 81)
curves = {"pre": [], "change": [], "replay": []}
for s, fp, m in zip(S, FP, MEAS):
df = CELLS[CELLS.session == s["session"]].reset_index(drop=True)
order = np.argsort(df.rate_session.values)
N = len(order); w = max(2, int(round(0.2 * N)))
cp, cc, cr = [], [], []
for p in POS:
a = int(round(p / 100 * N)); idx = order[a:a + w]
iu = np.triu_indices(len(idx), 1)
pre_ = m["C_pre"][np.ix_(idx, idx)][iu]; post_ = m["C_post"][np.ix_(idx, idx)][iu]; run_ = m["C_run"][np.ix_(idx, idx)][iu]
cp.append(np.nanmean(pre_)); cc.append(np.nanmean(post_ - pre_)); cr.append(g.partial_corr(run_, post_, pre_))
curves["pre"].append(cp); curves["change"].append(cc); curves["replay"].append(cr)
CURVES = {k: np.array(v) for k, v in curves.items()}
COACT = {}
for k, v in CURVES.items():
mean = np.nanmean(v, 0)
rho = stats.spearmanr(POS, mean)
per = [stats.spearmanr(POS, row, nan_policy="omit").statistic for row in v]
COACT[k] = (rho.statistic, per)
print(f"{k}: averaged curve slowest window {mean[0]:.4f}, fastest window {mean[-1]:.4f}, range {np.nanmin(mean):.4f} "
f"to {np.nanmax(mean):.4f}; Spearman rho with window position {rho.statistic:.3f} (p={rho.pvalue:.2g}); "
f"per session rho {np.round(per, 2).tolist()}; positive in {sum(r > 0 for r in per)} of 8")
fig, axs = plt.subplots(1, 3, figsize=(11, 3))
for ax, k in zip(axs, CURVES):
ax.plot(POS, CURVES[k].T, color="0.8", lw=0.8); ax.plot(POS, np.nanmean(CURVES[k], 0), color="k")
ax.set_title(k); ax.set_xlabel("window start (rate percentile)")
print("all values are within the axes (autoscaled)")
pre: averaged curve slowest window 0.0208, fastest window 0.0486, range 0.0204 to 0.0518; Spearman rho with window position 0.861 (p=6.9e-25); per session rho [0.69, -0.07, 0.72, 0.16, 0.5, -0.08, 0.21, 0.92]; positive in 6 of 8 change: averaged curve slowest window 0.0019, fastest window -0.0007, range -0.0026 to 0.0090; Spearman rho with window position -0.068 (p=0.55); per session rho [-0.65, -0.06, 0.11, -0.15, -0.14, 0.35, 0.89, -0.62]; positive in 3 of 8 replay: averaged curve slowest window 0.1858, fastest window 0.1146, range 0.0646 to 0.2959; Spearman rho with window position 0.095 (p=0.4); per session rho [-0.22, -0.67, -0.65, 0.17, -0.74, 0.68, 0.21, 0.71]; positive in 4 of 8
all values are within the axes (autoscaled)
By subject, the averaged curves' Spearman correlations.
subj = np.array([s["subject"] for s in S])
for k, v in CURVES.items():
print(k, {sb: round(stats.spearmanr(POS, np.nanmean(v[subj == sb], 0)).statistic, 3) for sb in np.unique(subj)})
fastest = CURVES["change"][:, -1]
print("POST minus PRE change in the fastest window per session:", np.round(fastest, 4).tolist(),
"sign-rank p =", round(stats.wilcoxon(fastest).pvalue, 3))
pre {np.str_('Achilles'): np.float64(0.728), np.str_('Buddy'): np.float64(0.716), np.str_('Cicero'): np.float64(0.276), np.str_('Gatsby'): np.float64(0.811)}
change {np.str_('Achilles'): np.float64(-0.45), np.str_('Buddy'): np.float64(0.111), np.str_('Cicero'): np.float64(0.021), np.str_('Gatsby'): np.float64(0.075)}
replay {np.str_('Achilles'): np.float64(-0.588), np.str_('Buddy'): np.float64(-0.646), np.str_('Cicero'): np.float64(-0.07), np.str_('Gatsby'): np.float64(0.671)}
POST minus PRE change in the fastest window per session: [-0.0013, 0.0008, 0.0284, 0.0115, 0.0089, 0.0002, 0.0135, -0.0675] sign-rank p = 0.383
Result: reproduced (fig4c-pre-fast-coactivation). PRE coactivation rises with firing rate, from 0.0208 in the slowest window to 0.0486 in the fastest; the Spearman correlation on the averaged curve is 0.861 and is positive in 6 of 8 sessions. The change in the fastest window does not differ from 0 over sessions (sign-rank p = 0.383), consistent with the statement that these correlations remain into POST.
Result: not reproduced (fig4c-slow-coactivation-increase). The POST minus PRE change does not fall with firing rate (Spearman correlation -0.068 on the averaged curve, negative in 5 of 8 sessions); the criterion needs -0.5 or below.
Result: not reproduced (fig4c-slow-replay). The partial correlation of MAZE running and POST coactivation given PRE is 0.186 in the slowest window and 0.115 in the fastest, but the Spearman correlation over windows is 0.095 and is negative in only 4 of 8 sessions; the criterion needs -0.5 or below and 5 sessions.
rate-split-medians¶
Place cells divided into equal groups by sleep firing rate have medians of 0.39 Hz (low-rate) and 1.12 Hz (high-rate).
Sleep firing rate is the rate in non-REM sleep of PRE and POST. Reading (a): median split within each session; reading (b): median split of all place cells pooled. The rate over all of PRE and POST is a variant.
def split_medians(col):
out = {}
lo_a = CELLS.groupby("session")[col].transform(lambda x: x.rank(method="first") <= len(x) // 2)
out["a"] = (CELLS.loc[lo_a, col].median(), CELLS.loc[~lo_a, col].median())
lo_b = CELLS[col].rank(method="first") <= len(CELLS) // 2
out["b"] = (CELLS.loc[lo_b, col].median(), CELLS.loc[~lo_b, col].median())
return out, lo_a, lo_b
SPLIT, LOW_A, LOW_B = split_medians("rate_sleep")
for r, (lo, hi) in SPLIT.items():
print(f"reading ({r}): low median {lo:.3f} Hz, high median {hi:.3f} Hz; within tolerance: "
f"{abs(lo - 0.39) <= 0.039 and abs(hi - 1.12) <= 0.112}")
v, _, _ = split_medians("rate_prepost")
print("variant, rate over all of PRE and POST:", {k: tuple(round(x, 3) for x in t) for k, t in v.items()})
v, _, _ = split_medians("rate_session")
print("variant, session rate:", {k: tuple(round(x, 3) for x in t) for k, t in v.items()})
print("sleep rate of all place cells: median", round(CELLS.rate_sleep.median(), 3), "range",
round(CELLS.rate_sleep.min(), 3), "to", round(CELLS.rate_sleep.max(), 3))
print(CELLS.assign(low=LOW_A).groupby(["subject", "low"]).rate_sleep.median().unstack().round(3).to_string())
reading (a): low median 0.435 Hz, high median 1.091 Hz; within tolerance: False
reading (b): low median 0.435 Hz, high median 1.092 Hz; within tolerance: False
variant, rate over all of PRE and POST: {'a': (np.float64(0.395), np.float64(1.041)), 'b': (np.float64(0.398), np.float64(1.046))}
variant, session rate: {'a': (np.float64(0.43), np.float64(1.061)), 'b': (np.float64(0.431), np.float64(1.064))}
sleep rate of all place cells: median 0.671 range 0.121 to 4.2
low False True
subject
Achilles 0.855 0.408
Buddy 0.880 0.374
Cicero 1.373 0.533
Gatsby 1.280 0.510
Result: partially reproduced. Split by non-REM sleep rate, the low-rate median is 0.435 Hz under both readings, just outside 0.39 +/- 0.039, and the high-rate median 1.091 (a) or 1.092 (b), within 1.12 +/- 0.112. With the rate over all of PRE and POST (a variant) both are within tolerance (0.395 and 1.041). The split is of 417 place cells, not the paper's 491.
figs16-low-rate-info-gain¶
Low-rate neurons carried more spatial information per spike on the maze, and increased their within-ripple firing-rate gain from PRE to POST more than high-rate cells.
def low_high(col, low):
A = CELLS.loc[low, col].dropna(); B = CELLS.loc[~low, col].dropna()
p = stats.mannwhitneyu(A, B).pvalue
print(f"{col}: low-rate n={len(A)} median {A.median():.3f}; high-rate n={len(B)} median {B.median():.3f}; rank-sum p={p:.3g}")
return A.median() > B.median() and p < 0.05
print("reading (a), split within each session:")
ok_a = [low_high("info", LOW_A), low_high("gain_change", LOW_A)]
print("reading (b), pooled split:")
ok_b = [low_high("info", LOW_B), low_high("gain_change", LOW_B)]
print("held under (a):", ok_a, " under (b):", ok_b)
print(CELLS.assign(low=LOW_A).groupby(["subject", "low"])[["info", "gain_change"]].median().unstack().round(3).to_string())
reading (a), split within each session:
info: low-rate n=207 median 0.975; high-rate n=210 median 0.622; rank-sum p=4.06e-08
gain_change: low-rate n=207 median 0.251; high-rate n=210 median 0.155; rank-sum p=0.948
reading (b), pooled split:
info: low-rate n=208 median 0.947; high-rate n=209 median 0.582; rank-sum p=4.03e-09
gain_change: low-rate n=208 median 0.268; high-rate n=209 median 0.114; rank-sum p=0.503
held under (a): [np.True_, np.False_] under (b): [np.True_, np.False_]
info gain_change
low False True False True
subject
Achilles 0.677 1.059 0.370 0.577
Buddy 0.843 0.977 0.257 1.615
Cicero 0.438 0.701 -0.474 -0.643
Gatsby 0.527 0.863 0.884 0.317
Result: partially reproduced. Under reading (a), which the criterion selects as neither reading matched, low-rate cells carry more spatial information (median 0.975 against 0.622 bits per spike, p = 4.1e-08), but their ripple-gain change from PRE to POST is not larger than that of high-rate cells (0.251 against 0.155, p = 0.948). Reading (b) gives the same.
fig4d-slow-sequences¶
Slow-firing, but not fast-firing, neurons increased their contribution to sequences from PRE to POST sleep, in a decoding restricted to each group.
print(REPLAY.loc[["base", "low_rate", "high_rate"], ["n_PRE", "n_POST", "mean_PRE", "mean_POST", "diff",
"p_one_sided", "present", "sessions_POST_gt_PRE", "diff_within", "p_within"]].round(4).to_string())
for name in ("low_rate", "high_rate"):
print(name, "POST minus PRE mean rZ by subject:", by_subject(name))
print("linear-track sessions only:")
print(REPLAY_LIN.loc[["base", "low_rate", "high_rate"], COLS].round(4).to_string())
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 16448 23801 0.7951 0.7868 -0.0083 0.9160 False 8 0.1067 0.0000
low_rate 8212 12549 0.6998 0.6763 -0.0235 0.9641 False 6 0.1157 0.0000
high_rate 14812 21567 0.8081 0.7594 -0.0488 0.9999 False 6 0.0583 0.0001
low_rate POST minus PRE mean rZ by subject: {'Achilles': np.float64(0.054), 'Buddy': np.float64(-0.129), 'Cicero': np.float64(-0.001), 'Gatsby': np.float64(-0.114)}
high_rate POST minus PRE mean rZ by subject: {'Achilles': np.float64(0.004), 'Buddy': np.float64(0.147), 'Cicero': np.float64(0.0), 'Gatsby': np.float64(-0.228)}
linear-track sessions only:
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 10536 12640 1.1443 1.2452 0.1008 0.0000 True 5 0.1261 0.0000
low_rate 5308 6268 0.9704 1.0476 0.0772 0.0046 True 4 0.1224 0.0000
high_rate 9517 11652 1.1767 1.2201 0.0435 0.0175 True 3 0.0661 0.0004
Result: not reproduced. On the pooled test there is no replay with either half (low-rate POST minus PRE -0.024, p = 0.964; high-rate -0.049, p = 0.9999), as for all cells. In the within-session variant both halves show replay, larger for low-rate cells (+0.116 against +0.058, both p <= 0.0001), which is the partial pattern of the criterion. On the linear-track sessions alone the pooled test also finds replay in both halves, larger for low-rate cells (+0.077, p = 0.0046, against +0.044, p = 0.018).
fig4e-gain-sequences¶
PRE to POST increases in sequence content were limited to cells with high ripple-specific recruitment (gain).
print(REPLAY.loc[["base", "high_gain", "low_gain"], ["n_PRE", "n_POST", "mean_PRE", "mean_POST", "diff",
"p_one_sided", "present", "sessions_POST_gt_PRE", "diff_within", "p_within"]].round(4).to_string())
for name in ("high_gain", "low_gain"):
print(name, "POST minus PRE mean rZ by subject:", by_subject(name))
print("linear-track sessions only:")
print(REPLAY_LIN.loc[["base", "high_gain", "low_gain"], COLS].round(4).to_string())
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 16448 23801 0.7951 0.7868 -0.0083 0.9160 False 8 0.1067 0.0000
high_gain 12342 18081 0.7231 0.7222 -0.0009 0.7097 False 5 0.1114 0.0000
low_gain 12086 17653 0.8412 0.7304 -0.1108 1.0000 False 3 0.0155 0.1416
high_gain POST minus PRE mean rZ by subject: {'Achilles': np.float64(0.102), 'Buddy': np.float64(-0.025), 'Cicero': np.float64(-0.043), 'Gatsby': np.float64(-0.067)}
low_gain POST minus PRE mean rZ by subject: {'Achilles': np.float64(-0.105), 'Buddy': np.float64(0.165), 'Cicero': np.float64(0.035), 'Gatsby': np.float64(-0.371)}
linear-track sessions only:
n_PRE n_POST mean_PRE mean_POST diff p_one_sided present sessions_POST_gt_PRE diff_within p_within
variant
base 10536 12640 1.1443 1.2452 0.1008 0.0000 True 5 0.1261 0.0000
high_gain 7641 9040 1.0731 1.0965 0.0233 0.2495 False 2 0.0453 0.0288
low_gain 8169 10026 1.1554 1.1655 0.0101 0.2083 False 3 0.0497 0.0028
Result: not reproduced. On the pooled test neither half shows replay (high-gain -0.001, p = 0.71; low-gain -0.111, p = 1.0). In the within-session variant high-gain cells show replay (+0.111, p < 0.0001) and low-gain cells do not (+0.016, p = 0.14), the paper's pattern. On the linear-track sessions alone the pooled test finds replay in neither half (+0.023, p = 0.25; +0.010, p = 0.21) and the within-session variant finds it in both (+0.045, p = 0.029; +0.050, p = 0.0028), so the difference between the halves in the within-session variant over all sessions comes from the circular-track sessions.
Summary¶
| Claim | Outcome |
|---|---|
| sessions-and-tracks | reproduced |
| place-cell-count (491) | not reproduced (417; 480 at 0.5 Hz) |
| fig1c-post-gt-pre | not reproduced by the pooled test; higher within all 8 sessions |
| fig1c-maze-gt-pre | reproduced |
| fig2b-strong-contributors (216) | blocked, methods_underspecified |
| fig2-pre-post-overlap (73%) | blocked, methods_underspecified |
| fig2d-shuffle-high-dpcc | partially reproduced |
| fig2e-remove-plastic | blocked, methods_underspecified |
| figs2-track-distribution | blocked, methods_underspecified |
| fig3a-plastic-specificity | blocked, methods_underspecified |
| fig3b-plastic-improve | blocked, methods_underspecified |
| fig4a-rigid-faster | blocked, methods_underspecified |
| fig4a-plastic-ripple-burst-coactivation | blocked, methods_underspecified |
| fig4c-pre-fast-coactivation | reproduced |
| fig4c-slow-coactivation-increase | not reproduced |
| fig4c-slow-replay | not reproduced |
| rate-split-medians | partially reproduced |
| figs16-low-rate-info-gain | partially reproduced |
| fig4d-slow-sequences | not reproduced |
| fig4e-gain-sequences | not reproduced |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's:
- The supplementary Materials and Methods were not available. Place cells, candidate events (population bursts holding a ripple on one LFP channel), the decoder's details, the rZ shuffle, the PCC and its test, bursting, ripple gain and the field and in-field measures are common choices fixed before the analysis.
- Immobility for MAZE candidate events is taken from the 2D head speed, running for the rate maps from the speed of the linearized position, both at 5 cm/s.
- The sequence score is decoded with each direction's template on the linear tracks, the event taking the larger rZ, and with one template on the circular track. The circular sessions score lower in every epoch and supply a larger share of POST than of PRE events, so pooling all events, as the criteria fix it, hides the POST over PRE increase that every session shows and that the same test finds on each kind of track. The within-session variant recovers the paper's pattern for fig2d, fig2e and fig4e, and the linear-track sessions alone for fig2d and fig2e; for fig4d both give the criterion's partial pattern. The outcomes are judged on the pooled test.
- The PCC test fixed in the criterion, a signed-rank test of a cell's contributions over events, has a false-positive rate of 56%, mostly on the linear tracks, where the contribution is taken on the template the event was selected for. This blocks the eight claims that rest on the selection of strongly contributing cells. A best-template variant of the contribution does better (21%) but is still far from 5%.
- dPCC and the rigid and plastic classes are measured on the same events that are decoded again for fig2d and fig2e. A cross-validated control (cells chosen on even-numbered events, replay measured on odd-numbered ones) shows that, on the linear tracks, shuffling the top 10 to 20% of cells reduces replay out of sample without removing it, while removing plastic cells still removes it on the linear tracks and within sessions.
What was found about the dandiset:
- The data hold what the paper's replay analyses need: sorted units with a cell type, positions, epochs, sleep states and the LFP, all on one clock (shown by decoding position, ripple-burst coincidence and the theta/delta ratio of the sleep states). No candidate events, ripples, place cell lists or PCC values are stored, so the paper's selections cannot be checked against it.
- In versions 0.210812.1516 and 0.250624.0426 the position series store the sampling period (0.0256 s) as their rate; version 0.261005.2101 corrects it.
- In Gatsby_08022013 the POST epoch ends at 30413.6 s while the spikes, LFP and states run to 32002 s.
- On the circular track the animals ran toward decreasing linearized position; the description says the position increases clockwise but not which way they ran.
- The description lists neck EMG and accelerometer signals; the files hold neither.
- Electrode
locationis "unknown" for every channel; the units carry the hemisphere (lCA1 or rCA1).