Reproducing Tang, Shin & Jadhav (2023) from dandiset 000447¶
Paper. Tang W, Shin JD, Jadhav SP (2023). Geometric transformation of cognitive maps for
generalization across hippocampal-prefrontal circuits. Cell Reports 42, 112246.
doi:10.1016/j.celrep.2023.112246. The claims file
(claims/000447/10.1016-j.celrep.2023.112246.yaml) was written from the NIH author manuscript in
PubMed Central (PMC10124109) with its STAR Methods and supplement; the version of record was not
read. Method details the claims file lacks were taken from the same text and from the paper's code,
github.com/JadhavLab/GeometricTransformation at commit 5e55e69 (cited where used, in the cells and
in nfn.py).
Dataset. DANDI:000447, version 0.250724.2224. The paper cites no version (its data statement says the data are on Brandeis servers, on request), so the most recent published version is pinned. The other published version, 0.230316.2133, holds the same five files; their spike times, unit electrodes, trial tables, epochs and positions are identical byte for byte (checked while preparing, not re-run here), and the file sizes differ by about 169 kB of metadata.
What the dandiset holds. One file per rat (AM2, JS17, JS21, ZT2, JS34) for the N-F-N' day: three W-track epochs, a trials table (start and end well, inbound or outbound, correct), 2D position with speed at 30 frames per second, sorted units (spike times and a tetrode number) and the 1.5 kHz LFP. It does not hold the de novo learning sessions, the three F-N-F'-N' rats, linearized positions, ripple or theta times, waveforms, or cell labels.
Outside inputs. None. The paper's code is read for method details only.
| Claim | Tested or blocked |
|---|---|
| fig1b-faster-novel-learning | blocked, data missing (de novo sessions) |
| fig1-denovo-final-performance | blocked, data missing (de novo sessions) |
| fig1c-session-performance | blocked, data missing (de novo sessions) |
| fig1c-familiar-above-novel | tested on N, F, N' (the F-N-F'-N' rats are missing) |
| figs1b-first-session-inbound-outbound | blocked, data missing (de novo sessions) |
| cells-per-session | tested |
| figs3a-ca1-global-remapping | tested |
| figs3a-pfc-similar-ca1-lower | tested |
| fig2d-pv-similarity | tested |
| fig2e-ca1-error-familiar | tested |
| fig2e-pfc-error-stable | tested |
| fig2f-switch-dynamics | tested |
| figs3c-pfc-assimilation | tested |
| figs3d-pfc-similarity-over-trials | tested |
| figs4-ramping-fraction | tested |
| figs4-ramping-generalization | run; blocked for too few ramping cells per population |
| figs4d-reward-location-similarity | run; blocked, methods underspecified (the outcome changes with the smoothing width) |
| figs5d-path-equivalence | tested |
| figs5e-path-equivalence-stable | run; blocked, methods underspecified (smoothing width) |
| figs5f-path-equivalence-vs-selectivity | run; blocked, methods underspecified (signed or absolute SI) |
| fig3d-selectivity-index | run; blocked, methods underspecified (smoothing width) |
| fig3e-selectivity-stable | run; blocked, methods underspecified (smoothing width) |
| manifold-cell-counts | tested |
| fig3fg-manifold-shapes | tested (through 3H and 3I) |
| fig3h-inseq-outseq-distance | tested |
| fig3i-phase-location-distance | tested |
| fig3j-familiar-novel-manifold | tested |
| figs6-high-dimensional-distances | tested |
| fig4d-pfc-ccgp | run; blocked, methods underspecified (smoothing width of the single-trial maps) |
| fig4d-ca1-ccgp-chance | run; blocked, as it depends on fig4d-pfc-ccgp |
| figs7b-geometric-null | run; blocked, as it depends on fig4d-pfc-ccgp |
| figs7d-dichotomies-decodable | run (50 shuffles in place of 1,000); blocked, as it depends on fig4d-pfc-ccgp |
The blocked behaviour claims need the de novo learning sessions, which the dandiset does not hold and which the paper says are on Brandeis servers, available on request.
The smoothing width. The width of the 1D Gaussian that smooths the linearized rate maps is garbled in the paper's text (Methods, Spatial fields and linearization, "s = σ cm"), and the Methods give no smoothing for the single-trial maps of the manifolds. Every analysis on these maps is run at 8 cm (the 2D width given in the sentence before, the claims file's reading a) and at 4 and 2 cm. The claims file's rule for fig2d-pv-similarity (the outcome is that of 8 cm when the widths agree, and methods underspecified otherwise) is applied to every claim that uses the maps: where the widths give different outcomes, the claim is blocked as methods underspecified, with the result at every width. A claim that depends on a claim blocked this way is blocked for the same reason, with its results shown.
import os
import time
import warnings
# UMAP's numba code is compiled for the host CPU, and its optimization amplifies the rounding
# differences, so the embeddings differed between machines. Compile for a generic x86-64 target
# and fix the BLAS kernel and numpy's dispatch, before numpy and numba are imported.
os.environ["NUMBA_CPU_NAME"] = "generic"
os.environ["OPENBLAS_CORETYPE"] = "Haswell"
os.environ["NPY_DISABLE_CPU_FEATURES"] = "AVX512F AVX512CD AVX512_SKX AVX512_CLX AVX512_CNL AVX512_ICL AVX512_SPR"
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import nfn
warnings.filterwarnings("ignore")
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
started = time.time()
A = nfn.load_all()
for d in A:
nfn.linearize(d)
BY = {d["animal"]: d for d in A}
SMOOTHS = (8.0, 4.0, 2.0)
MAPS = {(d["animal"], s, ep): nfn.rate_maps(d, ep, np.arange(d["n_units"]), smooth_cm=s)
for d in A for s in SMOOTHS for ep in range(3)}
print(f"loaded {len(A)} files from dandiset {nfn.DANDISET_ID} version {nfn.DANDISET_VERSION}, "
f"linearized and mapped in {time.time() - started:.0f} s")
loaded 5 files from dandiset 000447 version 0.250724.2224, linearized and mapped in 16 s
Data check¶
Each check of the metadata against the values is stated under its output.
urls = nfn.asset_urls()
rows = []
for d in A:
n, rate = nfn.epoch_counts(d)
tr = d["trials"]
rows.append({"animal": d["animal"], "file": urls[d["animal"]][0].split("/")[-1], "units": d["n_units"],
"electrode rows": len(d["electrode_labels"]),
"tetrode labels": len(set(d["electrode_labels"])),
"max units/electrodes value": int(d["unit_electrode_value"].max()),
"CA1 units (tetrode reading)": int((d["unit_region"] == "CA1").sum()),
"PFC units (tetrode reading)": int((d["unit_region"] == "PFC").sum()),
"CA1 units (row reading)": int((d["unit_region_rowreading"] == "CA1").sum()),
"PFC units (row reading)": int((d["unit_region_rowreading"] == "PFC").sum()),
"epoch durations (s)": np.round(d["epochs"][:, 1] - d["epochs"][:, 0]).astype(int).tolist(),
"trials per epoch": tr.groupby("epoch").size().tolist(),
"units with no spike in N": int((n[:, 0] == 0).sum())})
print(pd.DataFrame(rows).set_index("animal").T.to_string())
spk = np.concatenate([np.concatenate(d["spikes"]) for d in A])
inside = np.concatenate([np.concatenate([((s[:, None] >= d["epochs"][:, 0]) & (s[:, None] <= d["epochs"][:, 1])).any(1)
for s in d["spikes"]]) for d in A])
print(f"\nspikes inside the three epochs: {inside.mean():.4f} of {len(spk)}")
print("position unit:", A[0]["pos_unit"], "| comments:", A[0]["pos_comments"])
for d in A:
sp = d["speed"]
dt = np.diff(d["t"])
v_from_xy = np.hypot(np.diff(d["x"]), np.diff(d["y"])) / dt
ok = np.isfinite(v_from_xy) & (dt < 0.05)
print(f"{d['animal']}: x {np.nanmin(d['x']):.0f} to {np.nanmax(d['x']):.0f}, y {np.nanmin(d['y']):.0f} to {np.nanmax(d['y']):.0f} cm; "
f"frame interval median {np.median(dt) * 1000:.1f} ms; speed column median {np.nanmedian(sp):.1f}, 99th pct {np.nanpercentile(sp, 99):.1f}; "
f"median speed from x, y {np.nanmedian(v_from_xy[ok]):.1f} cm/s, correlation with the speed column {np.corrcoef(sp[1:][ok], v_from_xy[ok])[0, 1]:.2f}")
for d in A:
nfn.linearize(d) if "wells" not in d else None
print(f"{d['animal']}: session_start_time {d['session_start']}; position reference frame {d['reference_frame']!r}; center-well position (cm) per epoch "
+ ", ".join(f"{nfn.EPOCHS[ep]} ({w[1][0]:.0f}, {w[1][1]:.0f})" for ep, w in d["wells"].items()))
print("\ntrajectory types per epoch (track 0 = start and end wells not one of the four types):")
for d in A:
print(d["animal"], d["trials"].groupby(["epoch", "track"]).size().unstack(fill_value=0).rename(columns={0: "other", **nfn.TRACKS}).to_dict("index"))
animal AM2 JS17 JS21 ZT2 JS34
file sub-JDS-NFN-AM2_behavior+ecephys.nwb sub-JDS-NFN-JS17_behavior+ecephys.nwb sub-JDS-NFN-JS21_behavior+ecephys.nwb sub-JDS-NFN-ZT2_behavior+ecephys.nwb sub-JDS-NFN-JS34_behavior+ecephys.nwb
units 57 97 61 106 67
electrode rows 62 64 64 128 64
tetrode labels 31 32 32 64 32
max units/electrodes value 31 32 30 64 32
CA1 units (tetrode reading) 32 31 28 48 54
PFC units (tetrode reading) 25 66 33 58 13
CA1 units (row reading) 32 50 21 80 32
PFC units (row reading) 25 47 40 26 35
epoch durations (s) [878, 899, 946] [944, 1022, 1012] [949, 983, 950] [1324, 1554, 1458] [1105, 1021, 1138]
trials per epoch [25, 49, 30] [33, 65, 57] [30, 48, 38] [42, 24, 60] [69, 82, 88]
units with no spike in N 1 7 7 19 17
spikes inside the three epochs: 1.0000 of 3705112
position unit: centimeters; centimeters/second | comments: Video frames taken at 30fps - row1:xpos row2:ypos row3:velocity
AM2: x 15 to 131, y 0 to 106 cm; frame interval median 30.4 ms; speed column median 2.8, 99th pct 42.6; median speed from x, y 0.0 cm/s, correlation with the speed column 0.71
JS17: x 8 to 139, y 0 to 110 cm; frame interval median 30.5 ms; speed column median 2.7, 99th pct 49.5; median speed from x, y 7.7 cm/s, correlation with the speed column 0.75
JS21: x 17 to 133, y 0 to 115 cm; frame interval median 30.4 ms; speed column median 2.9, 99th pct 46.3; median speed from x, y 7.7 cm/s, correlation with the speed column 0.71
ZT2: x 10 to 140, y 1 to 116 cm; frame interval median 30.5 ms; speed column median 3.3, 99th pct 41.4; median speed from x, y 7.3 cm/s, correlation with the speed column 0.74
JS34: x 17 to 138, y 1 to 113 cm; frame interval median 30.4 ms; speed column median 8.3, 99th pct 48.3; median speed from x, y 12.8 cm/s, correlation with the speed column 0.76
AM2: session_start_time 2023-03-14T17:31:07.596000-04:00; position reference frame 'center well'; center-well position (cm) per epoch N (32, 52), F (73, 86), N' (38, 50)
JS17: session_start_time 2023-03-14T17:23:13.190000-04:00; position reference frame 'center well'; center-well position (cm) per epoch N (75, 96), F (80, 96), N' (66, 87)
JS21: session_start_time 2023-03-14T17:27:16.913000-04:00; position reference frame 'center well'; center-well position (cm) per epoch N (38, 55), F (77, 93), N' (41, 56)
ZT2: session_start_time 2023-03-14T17:34:53.914000-04:00; position reference frame 'center well'; center-well position (cm) per epoch N (33, 64), F (88, 93), N' (28, 60)
JS34: session_start_time 2023-03-14T17:19:13.648000-04:00; position reference frame 'center well'; center-well position (cm) per epoch N (37, 57), F (82, 98), N' (32, 56)
trajectory types per epoch (track 0 = start and end wells not one of the four types):
AM2 {0: {'other': 10, 'OUT_L': 1, 'IN_L': 4, 'OUT_R': 7, 'IN_R': 3}, 1: {'other': 1, 'OUT_L': 13, 'IN_L': 14, 'OUT_R': 11, 'IN_R': 10}, 2: {'other': 5, 'OUT_L': 5, 'IN_L': 6, 'OUT_R': 8, 'IN_R': 6}}
JS17 {0: {'other': 8, 'OUT_L': 4, 'IN_L': 8, 'OUT_R': 9, 'IN_R': 4}, 1: {'other': 2, 'OUT_L': 14, 'IN_L': 13, 'OUT_R': 18, 'IN_R': 18}, 2: {'other': 1, 'OUT_L': 14, 'IN_L': 13, 'OUT_R': 14, 'IN_R': 15}}
JS21 {0: {'other': 8, 'OUT_L': 7, 'IN_L': 7, 'OUT_R': 4, 'IN_R': 4}, 1: {'other': 3, 'OUT_L': 11, 'IN_L': 11, 'OUT_R': 12, 'IN_R': 11}, 2: {'other': 3, 'OUT_L': 7, 'IN_L': 8, 'OUT_R': 11, 'IN_R': 9}}
ZT2 {0: {'other': 13, 'OUT_L': 10, 'IN_L': 2, 'OUT_R': 5, 'IN_R': 12}, 1: {'other': 7, 'OUT_L': 3, 'IN_L': 4, 'OUT_R': 6, 'IN_R': 4}, 2: {'other': 3, 'OUT_L': 16, 'IN_L': 14, 'OUT_R': 13, 'IN_R': 14}}
JS34 {0: {'other': 14, 'OUT_L': 13, 'IN_L': 14, 'OUT_R': 15, 'IN_R': 13}, 1: {'other': 7, 'OUT_L': 15, 'IN_L': 18, 'OUT_R': 23, 'IN_R': 19}, 2: {'other': 9, 'OUT_L': 22, 'IN_L': 20, 'OUT_R': 18, 'IN_R': 19}}
The units table has only spike times and an electrodes column. By the NWB format that column
holds row indices of the electrodes table, but its values run up to 64 in ZT2, whose table has 128
rows (two per tetrode), and up to 32 in the others (62 to 64 rows): the values are tetrode numbers.
Read as tetrode numbers through the label column, the regions give CA1 and PFC counts close to the
paper's (next cells, and claims cells-per-session and manifold-cell-counts); read as row indices, the
counts differ by up to 32 units per animal (ZT2). The tetrode reading is used throughout; the
format error is recorded as a finding. There are no waveforms, cell numbers within a tetrode, or
putative interneuron labels.
Each file has three epochs, read as N, F and N' in that order, as the dandiset's description says. All spikes lie inside the three epochs. Positions are in cm (the W-maze is about 80 by 80 cm in the Methods); the third position column is the speed in cm/s, as the comment says. Some units are listed after the tetrode-ordered block and have no spike in N; they were clustered only in later epochs.
# The paper's code excludes putative interneurons by hard-coded (tetrode, cell) lists. The files
# carry no cell numbers; numbering the units of a tetrode in their order in the file shows what the
# lists would pick.
EXCLUDE = {"AM2": [(6, 9)], "JS21": [(6, 3), (6, 6), (6, 7), (6, 8), (7, 1), (21, 1), (22, 1), (25, 1)],
"JS17": [(6, 1), (7, 1)], "ZT2": [(31, 1), (31, 5), (32, 2), (27, 1), (14, 10)], "JS34": []}
for d in A:
n, rate = nfn.epoch_counts(d)
ue = d["unit_electrode_value"]
order = np.array([(ue[:i] == v).sum() + 1 for i, v in enumerate(ue)])
for tet, cell in EXCLUDE[d["animal"]]:
hit = np.where((ue == tet) & (order == cell))[0]
desc = ", ".join(f"{d['unit_region'][i]} unit, rate N/F/N' {np.round(rate[i], 1).tolist()} Hz" for i in hit) or "no such unit"
print(f"{d['animal']} listed ({tet},{cell}): {desc}")
hi = np.where((rate >= nfn.CA1_MAX_RATE).any(1) & (d["unit_region"] == "CA1"))[0]
print(f" {d['animal']} CA1 units at 7 Hz or more in some epoch: "
+ "; ".join(f"tetrode {ue[i]} #{order[i]} {np.round(rate[i], 1).tolist()}" for i in hi))
AM2 listed (6,9): CA1 unit, rate N/F/N' [0.3, 0.4, 0.2] Hz AM2 CA1 units at 7 Hz or more in some epoch: JS17 listed (6,1): CA1 unit, rate N/F/N' [33.0, 40.2, 45.2] Hz JS17 listed (7,1): no such unit JS17 CA1 units at 7 Hz or more in some epoch: tetrode 6 #1 [33.0, 40.2, 45.2] JS21 listed (6,3): CA1 unit, rate N/F/N' [14.5, 7.7, 13.3] Hz JS21 listed (6,6): CA1 unit, rate N/F/N' [39.7, 28.3, 26.6] Hz JS21 listed (6,7): CA1 unit, rate N/F/N' [0.0, 0.0, 0.2] Hz JS21 listed (6,8): CA1 unit, rate N/F/N' [0.0, 0.0, 1.7] Hz JS21 listed (7,1): CA1 unit, rate N/F/N' [24.4, 24.5, 24.9] Hz JS21 listed (21,1): no such unit JS21 listed (22,1): no such unit JS21 listed (25,1): CA1 unit, rate N/F/N' [30.1, 31.7, 28.9] Hz JS21 CA1 units at 7 Hz or more in some epoch: tetrode 6 #3 [14.5, 7.7, 13.3]; tetrode 6 #4 [26.6, 25.4, 26.5]; tetrode 6 #5 [15.7, 20.3, 19.7]; tetrode 6 #6 [39.7, 28.3, 26.6]; tetrode 7 #1 [24.4, 24.5, 24.9]; tetrode 25 #1 [30.1, 31.7, 28.9] ZT2 listed (31,1): no such unit ZT2 listed (31,5): no such unit ZT2 listed (32,2): CA1 unit, rate N/F/N' [26.7, 33.4, 24.0] Hz ZT2 listed (27,1): CA1 unit, rate N/F/N' [0.0, 14.7, 1.4] Hz ZT2 listed (14,10): CA1 unit, rate N/F/N' [0.0, 0.0, 0.1] Hz ZT2 CA1 units at 7 Hz or more in some epoch: tetrode 14 #8 [34.0, 40.1, 36.1]; tetrode 32 #2 [26.7, 33.4, 24.0]; tetrode 27 #1 [0.0, 14.7, 1.4] JS34 CA1 units at 7 Hz or more in some epoch:
Where the code's lists name a unit that exists, it is mostly a CA1 unit firing at 14 to 45 Hz (JS17 (6,1); JS21 (6,3), (6,6), (7,1), (25,1); ZT2 (32,2) and (27,1), the latter 14.7 Hz in F only), which supports the tetrode reading. Other entries name tetrodes without units (JS21 21 and 22, ZT2 31) or low-rate units, so the numbering by order is not exact. The code's own rule, CA1 mean rate below 7 Hz in the epoch, removes every listed CA1 unit in the epochs where it fires at 7 Hz or more; ZT2 (27,1) stays in N and N', where it fires at 0 and 1.4 Hz. The lists are not applied. PFC units are kept whatever their rate, as the code keeps every unit tagged PFC.
# The acquisition ElectricalSeries and the processed LFP: are they the same recording?
for d in A:
import h5py, remfile
f = h5py.File(remfile.File(urls[d["animal"]][1]), "r")
a = f["acquisition/ElectricalSeries/data"]
b = f["processing/ecephys/LFP/ElectricalSeries/data"]
ta = f["acquisition/ElectricalSeries/timestamps"]
same = all(np.array_equal(a[i:i + 200, :], b[i:i + 200, :]) for i in (0, a.shape[0] // 2, a.shape[0] - 200))
rate = 1 / np.median(np.diff(ta[:20000]))
print(f"{d['animal']}: shapes {a.shape} and {b.shape}, sampling {rate:.0f} Hz, identical in 3 slices: {same}, "
f"stored at different offsets: {a.id.get_offset() != b.id.get_offset()}")
AM2: shapes (5005900, 62) and (5005900, 62), sampling 1500 Hz, identical in 3 slices: True, stored at different offsets: True
JS17: shapes (5432378, 64) and (5432378, 64), sampling 1500 Hz, identical in 3 slices: True, stored at different offsets: True
JS21: shapes (4881671, 64) and (4881671, 64), sampling 1500 Hz, identical in 3 slices: True, stored at different offsets: True
ZT2: shapes (6834880, 128) and (6834880, 128), sampling 1500 Hz, identical in 3 slices: True, stored at different offsets: True
JS34: shapes (5339810, 64) and (5339810, 64), sampling 1500 Hz, identical in 3 slices: True, stored at different offsets: True
The acquisition/ElectricalSeries is a second, byte-identical copy of the 1.5 kHz LFP (in every
file, at all three slices tested), not the 30 kHz spike-band recording. The LFP is not used here:
the analysis needs spikes and position only, apart from the exclusion of sharp-wave ripples, which is
not done (see the summary).
# Linearization: idealized W paths from the well positions (nfn.build_skeleton, an assumption,
# as the paths are not in the dandiset). Paper: "2-m-long linearized trajectories" (Figure 2B).
rows = []
for d in A:
for ep in range(3):
m = (d["epoch_of_sample"] == ep) & (d["speed"] > nfn.MIN_SPEED) & np.isfinite(d["lin_offset"])
rows.append({"animal": d["animal"], "epoch": nfn.EPOCHS[ep],
"OUT_L/IN_L length (cm)": round(d["track_length"][(ep, 1)]),
"OUT_R/IN_R length (cm)": round(d["track_length"][(ep, 3)]),
"run samples on a trajectory": int(m.sum()),
"distance to path, median (cm)": round(np.median(d["lin_offset"][m]), 1),
"90th pct": round(np.percentile(d["lin_offset"][m], 90), 1),
"99th pct": round(np.percentile(d["lin_offset"][m], 99), 1)})
print(pd.DataFrame(rows).to_string(index=False))
fig, axes = plt.subplots(3, 5, figsize=(16, 9))
for i, d in enumerate(A):
for ep in range(3):
ax = axes[ep, i]
m = (d["epoch_of_sample"] == ep) & (d["speed"] > nfn.MIN_SPEED)
ax.plot(d["x"][m], d["y"][m], ",", alpha=0.3)
for t in (1, 3):
p = d["templates"][(ep, t)]
ax.plot(p[:, 0], p[:, 1], "r-")
ax.set_aspect("equal"); ax.set_title(f"{d['animal']} {nfn.EPOCHS[ep]}", fontsize=9)
plt.tight_layout(); plt.show()
print("All run positions are plotted (none outside the axes).")
animal epoch OUT_L/IN_L length (cm) OUT_R/IN_R length (cm) run samples on a trajectory distance to path, median (cm) 90th pct 99th pct AM2 N 195 197 5878 5.6 12.8 30.4 AM2 F 183 184 9803 4.4 12.4 29.3 AM2 N' 184 188 8823 5.4 13.5 22.4 JS17 N 212 220 8271 5.5 13.9 20.3 JS17 F 199 195 11291 4.5 10.6 17.1 JS17 N' 187 193 12110 6.6 12.4 19.4 JS21 N 193 199 8716 6.4 14.7 37.5 JS21 F 195 192 8939 4.9 11.2 17.0 JS21 N' 199 202 11097 5.8 12.8 19.9 ZT2 N 195 195 12835 5.7 13.5 18.7 ZT2 F 207 201 10078 6.0 15.0 45.2 ZT2 N' 199 199 16925 5.6 12.6 18.0 JS34 N 202 204 14591 6.5 13.1 29.7 JS34 F 207 205 14904 4.1 11.6 16.9 JS34 N' 191 193 16056 6.1 11.3 16.2
All run positions are plotted (none outside the axes).
The idealized paths follow the arms the rats ran (figure). Every trajectory is 180 to 220 cm long, consistent with the paper's 2-m trajectories, and the median distance from a run position to its path is about 5 cm (the tracks are about 7 cm wide in the Methods).
# Clock check: decode N' positions from N' spikes with N' templates (in-sample) after shifting all
# spike times. If spikes and position share a clock the error is smallest at no shift.
rows = []
for d in A:
for reg in ("CA1", "PFC"):
cells = nfn.included(d, reg, 2)
m = MAPS[(d["animal"], 8.0, 2)]
cells = cells[nfn.peak_rate({t: (m[t][0], m[t][1][cells]) for t in m}) >= nfn.PEAK_HZ]
tm = {t: (m[t][0], m[t][1][cells]) for t in m}
row = {"animal": d["animal"], "region": reg, "cells": len(cells)}
for sh in (-2.0, -0.5, 0.0, 0.5, 2.0):
d2 = dict(d); d2["spikes"] = [s + sh for s in d["spikes"]]
df, _ = nfn.decode_session(d2, 2, cells, tm, n_shuffles=1, restrict_track=True)
row[f"shift {sh:+.1f} s"] = round(df.error.median(), 1)
rows.append(row)
clock = pd.DataFrame(rows)
print("Median in-sample decoding error (cm, decode restricted to the current trajectory):")
print(clock.to_string(index=False))
print("Rows whose smallest error is at zero shift:",
int((clock.filter(like="shift").idxmin(axis=1) == "shift +0.0 s").sum()), "of", len(clock))
Median in-sample decoding error (cm, decode restricted to the current trajectory): animal region cells shift -2.0 s shift -0.5 s shift +0.0 s shift +0.5 s shift +2.0 s AM2 CA1 23 44.2 14.5 6.3 14.0 45.2 AM2 PFC 23 42.3 14.4 6.1 13.3 43.8 JS17 CA1 15 63.7 23.0 6.4 23.2 58.7 JS17 PFC 48 63.1 19.9 5.1 21.0 74.0 JS21 CA1 16 44.7 18.4 6.6 16.0 46.3 JS21 PFC 20 43.1 20.8 10.3 20.1 49.2 ZT2 CA1 32 56.2 18.8 4.5 16.8 50.6 ZT2 PFC 39 53.3 18.0 5.1 17.2 63.2 JS34 CA1 38 76.3 24.1 3.8 22.8 69.1 JS34 PFC 10 62.9 30.0 20.1 30.8 81.8 Rows whose smallest error is at zero shift: 10 of 10
In every animal and region the in-sample decoding error is smallest with no shift and grows when the spikes are shifted by 0.5 or 2 s, so spikes and position share the clock in N'. Spikes and the trial table are checked the same way by the trial-by-trial decoding of Figure 2F below, which uses the trial boundaries.
Behaviour¶
fig1c-familiar-above-novel¶
In both groups of rats, performance in the familiar environment is higher than in the novel one whatever the order of exposure (N-F-N' group in Figure 1C, F-N-F'-N' group in Figure S1C).
The de novo sessions and the F-N-F'-N' rats are not in the dandiset. The test is run on N, F and N' alone: a repeated-measures ANOVA over three sessions in place of the five of Figure 1C, with Tukey's test on its error term. The novel curve is estimated on N and N' run together and the familiar curve on F alone (cal_behavperform.m runs the novel sessions together; its familiar run would include the de novo days, which are missing). The state-space model's background probability is 1/3 (cal_behavperform.m).
perf, perf_sep = [], []
for d in A:
tr = d["trials"]
nov = tr[tr.epoch.isin([0, 2])]
p_nov = nfn.state_space_correct(nov.correct.values, background=1 / 3)
p_f = nfn.state_space_correct(tr[tr.epoch == 1].correct.values, background=1 / 3)
perf.append([p_nov[nov.epoch.values == 0].mean(), p_f.mean(), p_nov[nov.epoch.values == 2].mean()])
perf_sep.append([nfn.state_space_correct(tr[tr.epoch == e].correct.values, 1 / 3).mean() for e in range(3)])
raw = pd.DataFrame([[d["trials"][d["trials"].epoch == e].correct.mean() for e in range(3)] for d in A],
index=nfn.ANIMALS, columns=["N", "F", "N'"])
P = pd.DataFrame(perf, index=nfn.ANIMALS, columns=["N", "F", "N'"])
print("Fraction of correct trials:"); print(raw.round(3))
print("\nMean state-space estimate per session (novel sessions estimated together):"); print(P.round(3))
a, t = nfn.rm_anova_tukey(P.values, list(P.columns))
print(f"\nrepeated-measures ANOVA F({a['df'][0]},{a['df'][1]}) = {a['F']:.2f}, p = {a['p']:.3g}"); print(t.round(4).to_string(index=False))
P2 = pd.DataFrame(perf_sep, index=nfn.ANIMALS, columns=["N", "F", "N'"])
a2, t2 = nfn.rm_anova_tukey(P2.values, list(P2.columns))
print(f"\nVariant, each session estimated alone: F = {a2['F']:.2f}, p = {a2['p']:.3g}"); print(t2.round(4).to_string(index=False))
print("\nPer rat, F minus N:", (P["F"] - P["N"]).round(3).to_dict())
Fraction of correct trials:
N F N'
AM2 0.560 0.816 0.767
JS17 0.697 0.831 0.860
JS21 0.600 0.854 0.816
ZT2 0.643 0.500 0.817
JS34 0.710 0.780 0.830
Mean state-space estimate per session (novel sessions estimated together):
N F N'
AM2 0.647 0.816 0.688
JS17 0.718 0.832 0.835
JS21 0.631 0.854 0.784
ZT2 0.683 0.401 0.784
JS34 0.700 0.778 0.828
repeated-measures ANOVA F(2,8) = 1.08, p = 0.383
A B mean_B_minus_A q p_tukey
N F 0.0604 1.1614 0.7012
N N' 0.1080 2.0785 0.3540
F N' 0.0477 0.9170 0.7984
Variant, each session estimated alone: F = 1.66, p = 0.25
A B mean_B_minus_A q p_tukey
N F 0.0965 1.7527 0.4649
N N' 0.1381 2.5088 0.2382
F N' 0.0416 0.7561 0.8569
Per rat, F minus N: {'AM2': 0.169, 'JS17': 0.114, 'JS21': 0.222, 'ZT2': -0.282, 'JS34': 0.078}
# Behaviour check: well-to-well transits counted from position, next to the trials table
rows = []
for d in A:
tr = d["trials"]
for ep in range(3):
a, b = d["epochs"][ep]
m = (d["t"] >= a) & (d["t"] <= b)
xy = np.c_[d["x"][m], d["y"][m]]
W = nfn.well_positions(d, ep)
dist = np.stack([np.hypot(*(xy - W[w]).T) for w in (1, 2, 3)], 1)
at = np.where(dist.min(1) < 15, dist.argmin(1) + 1, 0)
seq = at[at > 0]
visits = seq[np.r_[True, seq[1:] != seq[:-1]]] if len(seq) else seq
t_ep = tr[tr.epoch == ep]
gaps = (t_ep.start_time.values[1:] - t_ep.stop_time.values[:-1])
rows.append({"animal": d["animal"], "session": nfn.EPOCHS[ep], "duration s": round(b - a),
"trials in table": len(t_ep), "transits from position": max(len(visits) - 1, 0),
"time in trials s": round((t_ep.stop_time - t_ep.start_time).sum()),
"first trial starts s": round(t_ep.start_time.min() - a), "last trial ends s": round(b - t_ep.stop_time.max()),
"longest gap s": round(gaps.max()) if len(gaps) else None,
"time above 5 cm/s, fraction": round(float(np.mean(d["speed"][m] > 5)), 2)})
TRANS = pd.DataFrame(rows)
print("Well-to-well transits from position (a well visit = within 15 cm of the well's median end position):")
print(TRANS.to_string(index=False))
Well-to-well transits from position (a well visit = within 15 cm of the well's median end position): animal session duration s trials in table transits from position time in trials s first trial starts s last trial ends s longest gap s time above 5 cm/s, fraction AM2 N 878 25 22 878 0 0 0 0.44 AM2 F 899 49 50 899 0 0 0 0.39 AM2 N' 946 30 28 946 0 0 0 0.40 JS17 N 944 33 31 944 0 0 0 0.42 JS17 F 1022 65 64 1022 0 0 0 0.41 JS17 N' 1012 57 62 1012 0 0 0 0.41 JS21 N 949 30 28 949 0 0 0 0.50 JS21 F 983 48 48 983 0 0 0 0.34 JS21 N' 950 38 37 950 0 0 0 0.43 ZT2 N 1324 42 40 1324 0 0 0 0.47 ZT2 F 1554 24 22 1554 0 0 0 0.33 ZT2 N' 1458 60 60 1458 0 0 0 0.42 JS34 N 1105 69 66 1105 0 0 0 0.61 JS34 F 1021 82 81 1021 0 0 0 0.59 JS34 N' 1138 88 89 1138 0 0 0 0.55
Result: not reproduced (N-F-N' rats only). On N, F and N' the repeated-measures ANOVA gives F(2,8) = 1.08, p = 0.383, and Tukey's test does not put F above N (difference 0.060, p = 0.70). F is above N in four rats (0.078 to 0.222) and below it in ZT2 (-0.282), whose F session has 24 trials at 0.50 correct. Estimating each session alone gives the same (p = 0.46). ZT2's trials table is complete: counted from position, ZT2 made 22 well-to-well transits in F against 24 rows, the trials tile the whole 1554-s epoch, and ZT2 spent the smallest fraction of any session above 5 cm/s there (0.33), so the low count is the rat's behaviour, not a gap in the table. The F-N-F'-N' rats of Figure S1C are not in the dandiset, and the de novo sessions of the paper's five-group ANOVA are missing.
Cells¶
cells-per-session¶
Cells were tracked across environments, with 30.5 (SD 9.8) CA1 and 37.4 (SD 19.3) PFC cells per session per animal.
Cells are units with more than 100 spikes in the session, CA1 units with a mean rate below 7 Hz in it (claims file; the Methods say "at least 100", which is shown as a variant).
rows = []
for d in A:
for ep in range(3):
for reg in ("CA1", "PFC"):
rows.append({"animal": d["animal"], "session": nfn.EPOCHS[ep], "region": reg,
"n": len(nfn.included(d, reg, ep)), "n_ge100": len(nfn.included(d, reg, ep, strict=False))})
C = pd.DataFrame(rows)
tab = C.pivot_table(index=["animal", "session"], columns="region", values="n")
print(tab)
print("\nmean and SD over the 15 animal-sessions:")
print(pd.DataFrame({"mean": tab.mean(), "SD": tab.std(ddof=1)}).round(2))
tab2 = C.pivot_table(index=["animal", "session"], columns="region", values="n_ge100")
print("\nVariant, at least 100 spikes (Methods):", tab2.mean().round(2).to_dict(), "SD", tab2.std(ddof=1).round(2).to_dict())
once = {reg: [len(np.unique(np.concatenate([nfn.included(d, reg, ep) for ep in range(3)]))) for d in A] for reg in ("CA1", "PFC")}
print("Variant, each cell counted once per animal:", {r: (round(np.mean(v), 2), round(np.std(v, ddof=1), 2), v) for r, v in once.items()})
region CA1 PFC
animal session
AM2 F 29.0 25.0
N 29.0 25.0
N' 29.0 25.0
JS17 F 25.0 60.0
N 26.0 59.0
N' 23.0 63.0
JS21 F 13.0 32.0
N 12.0 32.0
N' 18.0 30.0
JS34 F 36.0 12.0
N 29.0 13.0
N' 43.0 12.0
ZT2 F 37.0 54.0
N 28.0 54.0
N' 36.0 52.0
mean and SD over the 15 animal-sessions:
mean SD
region
CA1 27.53 8.68
PFC 36.53 18.60
Variant, at least 100 spikes (Methods): {'CA1': 27.67, 'PFC': 36.53} SD {'CA1': 8.57, 'PFC': 18.6}
Variant, each cell counted once per animal: {'CA1': (np.float64(34.8), np.float64(11.58), [31, 29, 21, 43, 50]), 'PFC': (np.float64(38.6), np.float64(22.12), [25, 65, 32, 58, 13])}
Result: reproduced. Over the 15 animal-sessions there are 27.53 (SD 8.68) CA1 and 36.53 (SD 18.60) PFC cells per session per animal, within the stated 30.5 (SD 9.8) and 37.4 (SD 19.3). The Methods' "at least 100 spikes" gives 27.67 and 36.53; counting each cell once per animal gives 34.8 and 38.6.
Remapping (Figure 2, S3, S4)¶
fig2d-pv-similarity¶
PFC population-vector similarity across environments (F-N') is as high as across repeated exposures to the novel environment (N-N'), while CA1 PV similarity is lower across environments than within the novel environment, where the CA1 map is reinstated.
As cal_PVsimilarity.m: cells passing the inclusion rules in either session of the pair and spatially tuned (peak of the linearized map at least 3 Hz in either session); each track's map on 101 bins of normalized length, interior NaN runs filled, the four tracks concatenated, NaN set to 0; Pearson r per bin over cells; per animal the mean over bins. Smoothing 8 cm (reading a), 4 and 2 cm (variants).
def pv_cells(d, reg, e1, e2, s=8.0, rule=">="):
cells = np.union1d(nfn.included(d, reg, e1), nfn.included(d, reg, e2))
m1, m2 = MAPS[(d["animal"], s, e1)], MAPS[(d["animal"], s, e2)]
R1 = nfn.normalized({t: (m1[t][0], m1[t][1][cells]) for t in m1})
R2 = nfn.normalized({t: (m2[t][0], m2[t][1][cells]) for t in m2})
if rule == ">= native": # variant: peak on the native 2-cm maps, not on the 101 resampled bins
p1 = nfn.peak_rate({t: (m1[t][0], m1[t][1][cells]) for t in m1}); p2 = nfn.peak_rate({t: (m2[t][0], m2[t][1][cells]) for t in m2})
else:
p1, p2 = np.nanmax(R1, 1), np.nanmax(R2, 1)
tuned = ((p1 >= 3) | (p2 >= 3)) if rule.startswith(">=") else ((p1 > 3) | (p2 > 3))
return cells[tuned], R1[tuned], R2[tuned]
def pv_bins(d, reg, e1, e2, s=8.0, rule=">=", cells_mask=None):
cells, R1, R2 = pv_cells(d, reg, e1, e2, s, rule)
R1, R2 = np.nan_to_num(R1), np.nan_to_num(R2)
return np.array([np.corrcoef(R1[:, i], R2[:, i])[0, 1] if R1[:, i].std() > 0 and R2[:, i].std() > 0 else np.nan
for i in range(R1.shape[1])])
GROUPS = ["CA1 N-N'", "CA1 F-N'", "PFC N-N'", "PFC F-N'"]
PAIRS = {"N-N'": (0, 2), "F-N'": (1, 2), "F-N": (1, 0)}
for s, rule in [(8.0, ">="), (4.0, ">="), (2.0, ">="), (8.0, ">"), (8.0, ">= native")]:
M = pd.DataFrame([[np.nanmean(pv_bins(d, g.split()[0], *PAIRS[g.split()[1]], s=s, rule=rule)) for g in GROUPS] for d in A],
index=nfn.ANIMALS, columns=GROUPS)
ncell = pd.DataFrame([[len(pv_cells(d, g.split()[0], *PAIRS[g.split()[1]], s, rule)[0]) for g in GROUPS] for d in A],
index=nfn.ANIMALS, columns=GROUPS)
a, t = nfn.rm_anova_tukey(M.values, GROUPS)
print(f"--- smoothing {s:g} cm, peak rule {rule} 3 Hz ---")
print("mean PV similarity per animal:"); print(M.round(3))
print("tuned cells per animal and pair:", ncell.to_dict("list"))
print(f"repeated-measures ANOVA F({a['df'][0]},{a['df'][1]}) = {a['F']:.2f}, p = {a['p']:.2g}")
print(t[((t.A == "CA1 N-N'") & (t.B == "CA1 F-N'")) | ((t.A == "PFC N-N'") & (t.B == "PFC F-N'"))].round(4).to_string(index=False))
--- smoothing 8 cm, peak rule >= 3 Hz ---
mean PV similarity per animal:
CA1 N-N' CA1 F-N' PFC N-N' PFC F-N'
AM2 0.544 0.049 0.790 0.813
JS17 0.532 0.239 0.926 0.929
JS21 0.652 -0.138 0.727 0.647
ZT2 0.419 0.062 0.880 0.848
JS34 0.441 0.012 0.952 0.943
tuned cells per animal and pair: {"CA1 N-N'": [26, 22, 17, 34, 42], "CA1 F-N'": [27, 22, 18, 38, 44], "PFC N-N'": [25, 56, 21, 47, 11], "PFC F-N'": [23, 53, 24, 45, 10]}
repeated-measures ANOVA F(3,12) = 67.53, p = 8.8e-08
A B mean_B_minus_A q p_tukey
CA1 N-N' CA1 F-N' -0.4726 10.2519 0.0001
PFC N-N' PFC F-N' -0.0190 0.4118 0.9910
--- smoothing 4 cm, peak rule >= 3 Hz ---
mean PV similarity per animal:
CA1 N-N' CA1 F-N' PFC N-N' PFC F-N'
AM2 0.508 0.046 0.757 0.791
JS17 0.527 0.221 0.914 0.919
JS21 0.593 -0.112 0.727 0.652
ZT2 0.411 0.065 0.860 0.834
JS34 0.439 0.007 0.939 0.930
tuned cells per animal and pair: {"CA1 N-N'": [27, 25, 18, 35, 44], "CA1 F-N'": [27, 23, 19, 40, 44], "PFC N-N'": [25, 57, 27, 50, 11], "PFC F-N'": [23, 56, 29, 51, 10]}
repeated-measures ANOVA F(3,12) = 88.81, p = 1.8e-08
A B mean_B_minus_A q p_tukey
CA1 N-N' CA1 F-N' -0.4506 11.3833 0.0000
PFC N-N' PFC F-N' -0.0140 0.3543 0.9942
--- smoothing 2 cm, peak rule >= 3 Hz ---
mean PV similarity per animal:
CA1 N-N' CA1 F-N' PFC N-N' PFC F-N'
AM2 0.463 0.050 0.713 0.774
JS17 0.490 0.214 0.897 0.907
JS21 0.554 -0.091 0.699 0.637
ZT2 0.391 0.062 0.836 0.809
JS34 0.431 0.006 0.921 0.913
tuned cells per animal and pair: {"CA1 N-N'": [27, 26, 19, 35, 46], "CA1 F-N'": [28, 26, 20, 40, 46], "PFC N-N'": [25, 61, 31, 53, 11], "PFC F-N'": [25, 58, 32, 52, 10]}
repeated-measures ANOVA F(3,12) = 92.89, p = 1.4e-08
A B mean_B_minus_A q p_tukey
CA1 N-N' CA1 F-N' -0.4176 11.1113 0.0000
PFC N-N' PFC F-N' -0.0054 0.1438 0.9996
--- smoothing 8 cm, peak rule > 3 Hz ---
mean PV similarity per animal:
CA1 N-N' CA1 F-N' PFC N-N' PFC F-N'
AM2 0.544 0.049 0.790 0.813
JS17 0.532 0.239 0.926 0.929
JS21 0.652 -0.138 0.727 0.647
ZT2 0.419 0.062 0.880 0.848
JS34 0.441 0.012 0.952 0.943
tuned cells per animal and pair: {"CA1 N-N'": [26, 22, 17, 34, 42], "CA1 F-N'": [27, 22, 18, 38, 44], "PFC N-N'": [25, 56, 21, 47, 11], "PFC F-N'": [23, 53, 24, 45, 10]}
repeated-measures ANOVA F(3,12) = 67.53, p = 8.8e-08
A B mean_B_minus_A q p_tukey
CA1 N-N' CA1 F-N' -0.4726 10.2519 0.0001
PFC N-N' PFC F-N' -0.0190 0.4118 0.9910
--- smoothing 8 cm, peak rule >= native 3 Hz ---
mean PV similarity per animal:
CA1 N-N' CA1 F-N' PFC N-N' PFC F-N'
AM2 0.544 0.049 0.790 0.813
JS17 0.532 0.239 0.926 0.929
JS21 0.652 -0.138 0.727 0.647
ZT2 0.419 0.062 0.880 0.848
JS34 0.441 0.012 0.952 0.943
tuned cells per animal and pair: {"CA1 N-N'": [26, 22, 17, 34, 42], "CA1 F-N'": [27, 22, 18, 38, 44], "PFC N-N'": [25, 56, 21, 47, 11], "PFC F-N'": [23, 53, 24, 45, 10]}
repeated-measures ANOVA F(3,12) = 67.53, p = 8.8e-08
A B mean_B_minus_A q p_tukey
CA1 N-N' CA1 F-N' -0.4726 10.2519 0.0001
PFC N-N' PFC F-N' -0.0190 0.4118 0.9910
Result: reproduced. With 8-cm smoothing, CA1 PV similarity is 0.419 to 0.652 within the novel room (N-N') and -0.138 to 0.239 across rooms (F-N') in the five rats; PFC is 0.727 to 0.952 (N-N') and 0.647 to 0.943 (F-N'). The repeated-measures ANOVA gives F(3,12) = 67.53, p = 8.8e-08; Tukey's test puts CA1 F-N' below N-N' (difference 0.473, p = 0.0001) and finds no PFC difference (p = 0.99). The 4-cm and 2-cm widths, the code's "> 3 Hz" rule and the peak taken on the native 2-cm maps give the same outcome (CA1 p <= 0.0001, PFC p >= 0.99).
figs3a-ca1-global-remapping and figs3a-pfc-similar-ca1-lower¶
CA1 spatial similarity between the two environments is no different from that of cell-ID shuffles, that is, CA1 remaps globally.
Single-cell rate overlap and spatial correlation of PFC cells across environments (F-N') are as large as across repeated exposures to the novel environment (N-N'), whereas in CA1 they are lower across environments.
Cells of the PV analysis (fig2d). Spatial similarity is the Pearson r of a cell's two concatenated normalized maps over bins valid in both. Rate overlap is the lower over the higher mean rate, with the mean rate read as (a) the mean of the linearized map or (b) spikes over run time (speed above 5 cm/s on a trajectory). The shuffle pairs each cell's F map with another cell's N' map of the same animal (a random derangement within the animal; the Methods say "a different cell"). One-way ANOVA with Tukey's test over N-N', F-N' and shuffle, per region and measure. The shuffle is one draw; the result over 100 draws is shown too.
def remap_measures(d, reg, e1, e2, rng=None, s=8.0):
cells, R1, R2 = pv_cells(d, reg, e1, e2, s)
dt = nfn._sample_dt(d)
run = []
for ep in (e1, e2):
m = nfn._run_mask(d, ep); T = dt[m].sum(); tt = d["t"]
run.append(np.array([m[np.clip(np.searchsorted(tt, d["spikes"][u]), 0, len(tt) - 1)].sum() / T for u in cells]))
if rng is not None:
k = len(cells); perm = np.arange(k)
while k > 1 and np.any(perm == np.arange(k)):
perm = rng.permutation(k)
R2, run[1] = R2[perm], run[1][perm]
ma, mb = np.nanmean(R1, 1), np.nanmean(R2, 1)
sc = []
for i in range(len(cells)):
ok = np.isfinite(R1[i]) & np.isfinite(R2[i])
sc.append(stats.pearsonr(R1[i][ok], R2[i][ok])[0] if ok.sum() > 2 and R1[i][ok].std() > 0 and R2[i][ok].std() > 0 else np.nan)
return pd.DataFrame({"overlap_map": np.minimum(ma, mb) / np.maximum(ma, mb),
"overlap_run": np.minimum(*run) / np.maximum(*run), "spatial_r": sc})
rng_s3a = np.random.default_rng(301)
REMAP = {}
for reg in ("CA1", "PFC"):
REMAP[reg] = {"N-N'": pd.concat([remap_measures(d, reg, 0, 2) for d in A]),
"F-N'": pd.concat([remap_measures(d, reg, 1, 2) for d in A]),
"shuffle": pd.concat([remap_measures(d, reg, 1, 2, rng_s3a) for d in A])}
S3A = {}
for reg in ("CA1", "PFC"):
for meas in ("spatial_r", "overlap_map", "overlap_run"):
a, t = nfn.anova_tukey([REMAP[reg][g][meas] for g in REMAP[reg]], list(REMAP[reg]))
S3A[(reg, meas)] = t
med = {g: round(float(REMAP[reg][g][meas].median()), 3) for g in REMAP[reg]}
n = {g: int(REMAP[reg][g][meas].notna().sum()) for g in REMAP[reg]}
print(f"{reg} {meas}: medians {med}, n {n}, ANOVA F = {a['F']:.2f}, p = {a['p']:.2g}")
print(t[["group1", "group2", "meandiff", "p-adj", "lower", "upper"]].to_string(index=False))
CA1 spatial_r: medians {"N-N'": 0.709, "F-N'": -0.059, 'shuffle': -0.064}, n {"N-N'": 109, "F-N'": 125, 'shuffle': 125}, ANOVA F = 170.52, p = 1.1e-52
group1 group2 meandiff p-adj lower upper
F-N' N-N' 0.5240 0.0000 0.4418 0.6062
F-N' shuffle -0.0736 0.0752 -0.1530 0.0057
N-N' shuffle -0.5976 0.0000 -0.6798 -0.5155
CA1 overlap_map: medians {"N-N'": 0.437, "F-N'": 0.298, 'shuffle': 0.298}, n {"N-N'": 141, "F-N'": 149, 'shuffle': 149}, ANOVA F = 3.09, p = 0.046
group1 group2 meandiff p-adj lower upper
F-N' N-N' 0.0530 0.3171 -0.0330 0.1389
F-N' shuffle -0.0377 0.5487 -0.1225 0.0471
N-N' shuffle -0.0906 0.0360 -0.1766 -0.0047
CA1 overlap_run: medians {"N-N'": 0.486, "F-N'": 0.322, 'shuffle': 0.286}, n {"N-N'": 141, "F-N'": 149, 'shuffle': 149}, ANOVA F = 3.59, p = 0.028
group1 group2 meandiff p-adj lower upper
F-N' N-N' 0.0453 0.4470 -0.0426 0.1332
F-N' shuffle -0.0546 0.3006 -0.1413 0.0321
N-N' shuffle -0.0999 0.0213 -0.1878 -0.0120
PFC spatial_r: medians {"N-N'": 0.546, "F-N'": 0.502, 'shuffle': 0.016}, n {"N-N'": 150, "F-N'": 150, 'shuffle': 150}, ANOVA F = 125.18, p = 6.8e-44
group1 group2 meandiff p-adj lower upper
F-N' N-N' 0.0361 0.522 -0.0420 0.1142
F-N' shuffle -0.4359 0.000 -0.5140 -0.3578
N-N' shuffle -0.4720 0.000 -0.5501 -0.3939
PFC overlap_map: medians {"N-N'": 0.697, "F-N'": 0.707, 'shuffle': 0.217}, n {"N-N'": 160, "F-N'": 155, 'shuffle': 155}, ANOVA F = 88.92, p = 1.9e-33
group1 group2 meandiff p-adj lower upper
F-N' N-N' -0.0099 0.9414 -0.0804 0.0606
F-N' shuffle -0.3531 0.0000 -0.4242 -0.2821
N-N' shuffle -0.3432 0.0000 -0.4137 -0.2727
PFC overlap_run: medians {"N-N'": 0.706, "F-N'": 0.726, 'shuffle': 0.252}, n {"N-N'": 160, "F-N'": 155, 'shuffle': 155}, ANOVA F = 83.44, p = 1e-31
group1 group2 meandiff p-adj lower upper
F-N' N-N' -0.0177 0.8217 -0.0876 0.0521
F-N' shuffle -0.3428 0.0000 -0.4132 -0.2724
N-N' shuffle -0.3251 0.0000 -0.3949 -0.2552
# figs3a-ca1-global-remapping: CA1 F-N' against shuffle, spatial similarity (the criterion's measure)
t = S3A[("CA1", "spatial_r")]
row = t[(t.group1 == "F-N'") & (t.group2 == "shuffle")].iloc[0]
print(f"CA1 spatial similarity, F-N' against shuffle: mean difference {-row['meandiff']:.3f} (F-N' minus shuffle), Tukey p = {row['p-adj']:.4f}")
for meas in ("overlap_map", "overlap_run"):
r2 = S3A[("CA1", meas)]; r2 = r2[(r2.group1 == "F-N'") & (r2.group2 == "shuffle")].iloc[0]
print(f" {meas}: difference {-r2['meandiff']:.3f}, p = {r2['p-adj']:.3f}")
ps, diffs = [], []
rng_draws = np.random.default_rng(302)
for k in range(100):
sh = pd.concat([remap_measures(d, "CA1", 1, 2, rng_draws) for d in A])
a, t = nfn.anova_tukey([REMAP["CA1"]["N-N'"].spatial_r, REMAP["CA1"]["F-N'"].spatial_r, sh.spatial_r], ["N-N'", "F-N'", "shuffle"])
r = t[(t.group1 == "F-N'") & (t.group2 == "shuffle")].iloc[0]
ps.append(r["p-adj"]); diffs.append(-r["meandiff"])
ps, diffs = np.array(ps), np.array(diffs)
print(f"over 100 shuffle draws: Tukey p < 0.05 in {np.mean(ps < 0.05):.2f} of draws; median p {np.median(ps):.3f}; "
f"F-N' minus shuffle mean difference median {np.median(diffs):.3f} (range {diffs.min():.3f} to {diffs.max():.3f})")
print("per animal median spatial r, CA1 F-N':",
{d["animal"]: round(float(remap_measures(d, "CA1", 1, 2).spatial_r.median()), 3) for d in A})
CA1 spatial similarity, F-N' against shuffle: mean difference 0.074 (F-N' minus shuffle), Tukey p = 0.0752 overlap_map: difference 0.038, p = 0.549 overlap_run: difference 0.055, p = 0.301
over 100 shuffle draws: Tukey p < 0.05 in 0.13 of draws; median p 0.291; F-N' minus shuffle mean difference median 0.051 (range 0.014 to 0.105)
per animal median spatial r, CA1 F-N': {'AM2': -0.054, 'JS17': 0.312, 'JS21': -0.11, 'ZT2': -0.088, 'JS34': -0.079}
Result: reproduced. For the CA1 spatial correlation, Tukey's test finds no difference between F-N' and the cell-ID shuffle (F-N' higher by 0.074, p = 0.075). The margin is thin: over 100 shuffle draws p falls below 0.05 in 13 of them (median p 0.29), and F-N' is always higher than the shuffle (by 0.014 to 0.105). The 4-cm and 2-cm widths give the same outcome (p = 0.21 and 0.40). One of five rats (JS17, median r 0.31) keeps some CA1 spatial correlation across rooms.
for reg in ("PFC", "CA1"):
for meas in ("overlap_map", "overlap_run", "spatial_r"):
t = S3A[(reg, meas)]
r = t[(t.group1 == "F-N'") & (t.group2 == "N-N'")].iloc[0]
print(f"{reg} {meas}: N-N' minus F-N' = {r['meandiff']:.3f}, Tukey p = {r['p-adj']:.4f}")
for s in (4.0, 2.0):
out = []
for reg in ("PFC", "CA1"):
g = {"N-N'": pd.concat([remap_measures(d, reg, 0, 2, s=s) for d in A]), "F-N'": pd.concat([remap_measures(d, reg, 1, 2, s=s) for d in A]),
"shuffle": pd.concat([remap_measures(d, reg, 1, 2, np.random.default_rng(303), s=s) for d in A])}
for meas in ("overlap_map", "overlap_run", "spatial_r"):
a, t = nfn.anova_tukey([g[k][meas] for k in g], list(g))
r = t[(t.group1 == "F-N'") & (t.group2 == "N-N'")].iloc[0]
out.append(f"{reg} {meas} diff {r['meandiff']:.3f} p {r['p-adj']:.3g}")
if reg == "CA1" and meas == "spatial_r":
r = t[(t.group1 == "F-N'") & (t.group2 == "shuffle")].iloc[0]
out.append(f"CA1 spatial_r F-N' minus shuffle {-r['meandiff']:.3f} p {r['p-adj']:.3g}")
print(f"smoothing {s:g} cm:", "; ".join(out))
PFC overlap_map: N-N' minus F-N' = -0.010, Tukey p = 0.9414 PFC overlap_run: N-N' minus F-N' = -0.018, Tukey p = 0.8217 PFC spatial_r: N-N' minus F-N' = 0.036, Tukey p = 0.5220 CA1 overlap_map: N-N' minus F-N' = 0.053, Tukey p = 0.3171 CA1 overlap_run: N-N' minus F-N' = 0.045, Tukey p = 0.4470 CA1 spatial_r: N-N' minus F-N' = 0.524, Tukey p = 0.0000
smoothing 4 cm: PFC overlap_map diff -0.010 p 0.941; PFC overlap_run diff -0.014 p 0.88; PFC spatial_r diff 0.022 p 0.735; CA1 overlap_map diff 0.069 p 0.12; CA1 overlap_run diff 0.056 p 0.272; CA1 spatial_r diff 0.468 p 0; CA1 spatial_r F-N' minus shuffle 0.053 p 0.208
smoothing 2 cm: PFC overlap_map diff -0.010 p 0.94; PFC overlap_run diff -0.013 p 0.895; PFC spatial_r diff -0.001 p 1; CA1 overlap_map diff 0.058 p 0.237; CA1 overlap_run diff 0.048 p 0.388; CA1 spatial_r diff 0.418 p 0; CA1 spatial_r F-N' minus shuffle 0.036 p 0.404
Result: partially reproduced. PFC F-N' does not differ from N-N' in rate overlap (map mean, p = 0.94; run rate, p = 0.82) or spatial correlation (p = 0.52). CA1 F-N' is below N-N' in spatial correlation (difference 0.524, p < 0.0001) but not in rate overlap (p = 0.32 and 0.45), so the claim holds for one measure. The 4-cm and 2-cm widths give the same pattern (CA1 rate overlap p 0.12 to 0.39).
figs3c-pfc-assimilation¶
PFC PV similarity across environments grows with experience (F-N' higher than F-N), while CA1 does not change (F-N against F-N' not different) and stays below N-N'.
Per-bin PV similarity (fig2d), bins pooled over the five rats, one-way ANOVA per region with Tukey's test.
for s in SMOOTHS:
for reg in ("CA1", "PFC"):
G = {nm: np.concatenate([pv_bins(d, reg, *PAIRS[nm], s=s) for d in A]) for nm in ("F-N", "F-N'", "N-N'")}
a, t = nfn.anova_tukey(list(G.values()), list(G))
print(f"smoothing {s:g} cm, {reg}: means {({k: round(float(np.nanmean(v)), 3) for k, v in G.items()})}, "
f"bins {a['n']}, ANOVA F = {a['F']:.1f}, p = {a['p']:.2g}")
print(" " + t[["group1", "group2", "meandiff", "p-adj"]].round(4).to_string(index=False).replace("\n", "\n "))
smoothing 8 cm, CA1: means {'F-N': 0.062, "F-N'": 0.045, "N-N'": 0.517}, bins [2020, 2020, 2020], ANOVA F = 2372.2, p = 0
group1 group2 meandiff p-adj
F-N F-N' -0.0174 0.0649
F-N N-N' 0.4551 0.0000
F-N' N-N' 0.4726 0.0000
smoothing 8 cm, PFC: means {'F-N': 0.774, "F-N'": 0.836, "N-N'": 0.855}, bins [2020, 2020, 2020], ANOVA F = 173.3, p = 6.4e-74
group1 group2 meandiff p-adj
F-N F-N' 0.0623 0.0000
F-N N-N' 0.0813 0.0000
F-N' N-N' 0.0190 0.0001
smoothing 4 cm, CA1: means {'F-N': 0.061, "F-N'": 0.045, "N-N'": 0.496}, bins [2020, 2020, 2020], ANOVA F = 2110.3, p = 0
group1 group2 meandiff p-adj
F-N F-N' -0.0162 0.098
F-N N-N' 0.4343 0.000
F-N' N-N' 0.4506 0.000
smoothing 4 cm, PFC: means {'F-N': 0.76, "F-N'": 0.825, "N-N'": 0.839}, bins [2020, 2020, 2020], ANOVA F = 165.5, p = 1.1e-70
group1 group2 meandiff p-adj
F-N F-N' 0.0654 0.0000
F-N N-N' 0.0795 0.0000
F-N' N-N' 0.0140 0.0075
smoothing 2 cm, CA1: means {'F-N': 0.061, "F-N'": 0.048, "N-N'": 0.466}, bins [2019, 2019, 2020], ANOVA F = 1789.8, p = 0
group1 group2 meandiff p-adj
F-N F-N' -0.0124 0.2646
F-N N-N' 0.4052 0.0000
F-N' N-N' 0.4176 0.0000
smoothing 2 cm, PFC: means {'F-N': 0.738, "F-N'": 0.808, "N-N'": 0.813}, bins [2020, 2020, 2020], ANOVA F = 141.3, p = 1.1e-60
group1 group2 meandiff p-adj
F-N F-N' 0.0701 0.0000
F-N N-N' 0.0755 0.0000
F-N' N-N' 0.0054 0.5276
Result: reproduced. With bins pooled over rats (8 cm), PFC F-N' exceeds F-N (0.836 against 0.774, Tukey p < 0.0001), CA1 F-N and F-N' do not differ (0.062 and 0.045, p = 0.065), and CA1 F-N' is below N-N' (0.517, p < 0.0001). PFC F-N' is below N-N' (difference 0.019, p = 0.0001; the paper reports p = 0.0013). The 4-cm and 2-cm widths give the same outcome (CA1 F-N against F-N' p = 0.098 and 0.26).
figs3d-pfc-similarity-over-trials¶
Over trials in the novel environment, PFC similarity to the familiar map rises from the first to the last trial, while CA1 similarity does not change.
Single-trial maps (2-cm bins, smoothed at 8, 4 and 2 cm as the session maps) of each trajectory pass in the novel environment against the session-averaged F maps of the same trajectory, over the tuned cells (peak at least 3 Hz in F or in the novel sessions); per-bin PV similarity on 101 normalized bins. A trial is kept when at least 80% of these cells fire during its run samples (legend). The paper does not say how trials are numbered or which is "last". Here the kept trials of each animal are ranked in order; trial ranks 1 to K enter the test, K being the smallest number of kept trials over the animals that have any, and "last" is (i) rank K or (ii) each animal's own final kept trial (ranks 1 to K - 1 and the final trials as the groups). Readings for the novel trials: (a) N then N', (b) N, (c) N'. One-way ANOVA over ranks on per-bin values, Tukey's test of first against last. The session each kept trial comes from is shown for reading (a).
SORTED = {d["animal"]: [np.sort(x) for x in d["spikes"]] for d in A}
rows = []
for s in SMOOTHS:
for d in A:
tr = d["trials"]
mF = MAPS[(d["animal"], s, 1)]
for reg in ("CA1", "PFC"):
for name, eps in (("N then N'", [0, 2]), ("N", [0]), ("N'", [2])):
cells = np.union1d(nfn.included(d, reg, 1), np.concatenate([nfn.included(d, reg, e) for e in eps]))
pkF = nfn.peak_rate({t: (mF[t][0], mF[t][1][cells]) for t in mF})
pkN = np.max([nfn.peak_rate({t: (MAPS[(d["animal"], s, e)][t][0], MAPS[(d["animal"], s, e)][t][1][cells]) for t in mF}) for e in eps], 0)
cells = cells[(pkF >= 3) | (pkN >= 3)]
order = 0
for e in eps:
for i in tr[(tr.epoch == e) & (tr.track > 0)].index:
order += 1
track = int(tr.loc[i, "track"])
run = nfn._run_mask(d, e, track, i)
tt = d["t"][run]
if len(tt) == 0:
continue
ss_ = SORTED[d["animal"]]
active = np.mean([np.searchsorted(ss_[u], tt[-1], "right") > np.searchsorted(ss_[u], tt[0], "left") for u in cells])
c, r = nfn.rate_maps(d, e, cells, smooth_cm=s, trial=i, track_only=track)[track]
pos = np.linspace(0, 1, 101)
X = np.nan_to_num(r[:, np.abs((c / c.max())[None, :] - pos[:, None]).argmin(1)])
cF, rF = mF[track]
Y = np.nan_to_num(rF[cells][:, np.abs((cF / cF.max())[None, :] - pos[:, None]).argmin(1)])
pv = np.array([np.corrcoef(X[:, b], Y[:, b])[0, 1] if X[:, b].std() > 0 and Y[:, b].std() > 0 else np.nan for b in range(101)])
rows.append({"smooth": s, "animal": d["animal"], "region": reg, "reading": name, "order": order, "epoch": nfn.EPOCHS[e],
"active": active, "pv": pv})
TRI = pd.DataFrame(rows)
T8 = TRI[TRI.smooth == 8.0]
print("8 cm: passes and passes kept (at least 80% of tuned cells active):")
print(T8.groupby(["region", "reading", "animal"]).agg(passes=("order", "size"), kept=("active", lambda x: int((x >= 0.8).sum()))).unstack("animal"))
print("\n8 cm, N then N': kept passes per animal by session:")
print(T8[(T8.reading == "N then N'") & (T8.active >= 0.8)].groupby(["region", "animal", "epoch"]).size().unstack("epoch", fill_value=0))
S3D = []
for (s, reg, name), x in TRI[TRI.active >= 0.8].groupby(["smooth", "region", "reading"]):
x = x.copy(); x["rank"] = x.groupby("animal").order.rank().astype(int)
x["n_kept"] = x.groupby("animal")["rank"].transform("max")
K = int(x.n_kept.min())
# (i) ranks 1 to K, the last trial rank K (K the smallest number of kept passes over the animals)
xi = x[x["rank"] <= K]
groups = [np.concatenate(xi[xi["rank"] == k].pv.values) for k in range(1, K + 1)]
a, t = nfn.anova_tukey(groups, [f"{k:02d}" for k in range(1, K + 1)])
r = t[(t.group1 == "01") & (t.group2 == f"{K:02d}")]
# (ii) ranks 1 to K - 1, and each animal's own final kept pass as the last trial
groups2 = [np.concatenate(x[x["rank"] == k].pv.values) for k in range(1, K)] + [np.concatenate(x[x["rank"] == x.n_kept].pv.values)]
names2 = [f"{k:02d}" for k in range(1, K)] + ["last"]
if len(groups2) >= 2:
a2, t2 = nfn.anova_tukey(groups2, names2)
r2 = t2[((t2.group1 == "01") & (t2.group2 == "last")) | ((t2.group1 == "last") & (t2.group2 == "01"))]
sign2 = 1 if r2.group1.iloc[0] == "01" else -1
diff2, p2 = sign2 * r2.meandiff.iloc[0], r2["p-adj"].iloc[0]
else:
diff2, p2 = np.nan, np.nan
S3D.append({"smooth": s, "region": reg, "reading": name, "animals": x.animal.nunique(), "K": K,
"kept per animal": x.groupby("animal")["rank"].max().to_dict(),
"trial 1 mean": np.nanmean(groups[0]), "trial K mean": np.nanmean(groups[-1]), "ANOVA p": a["p"],
"(i) last - first": r.meandiff.iloc[0], "(i) Tukey p": r["p-adj"].iloc[0],
"own last mean": np.nanmean(groups2[-1]), "(ii) last - first": diff2, "(ii) Tukey p": p2})
S3D = pd.DataFrame(S3D)
with pd.option_context("display.float_format", lambda v: f"{v:.3g}", "display.max_colwidth", 60):
print(S3D.to_string(index=False))
8 cm: passes and passes kept (at least 80% of tuned cells active):
passes kept
animal AM2 JS17 JS21 JS34 ZT2 AM2 JS17 JS21 JS34 ZT2
region reading
CA1 N 15 25 22 55 29 15 24 20 0 0
N then N' 40 81 57 134 86 38 45 26 9 36
N' 25 56 35 79 57 25 39 26 19 38
PFC N 15 25 22 55 29 15 25 22 55 29
N then N' 40 81 57 134 86 40 81 57 132 86
N' 25 56 35 79 57 25 56 35 79 57
8 cm, N then N': kept passes per animal by session:
epoch N N'
region animal
CA1 AM2 14 24
JS17 24 21
JS21 0 26
JS34 0 9
ZT2 0 36
PFC AM2 15 25
JS17 25 56
JS21 22 35
JS34 55 77
ZT2 29 57
smooth region reading animals K kept per animal trial 1 mean trial K mean ANOVA p (i) last - first (i) Tukey p own last mean (ii) last - first (ii) Tukey p
2 CA1 N 3 15 {'AM2': 15, 'JS17': 23, 'JS21': 16} 0.115 0.041 1.02e-11 -0.0744 0.0633 0.0301 -0.0853 0.0098
2 CA1 N then N' 5 10 {'AM2': 39, 'JS17': 20, 'JS21': 23, 'JS34': 10, 'ZT2': 36} 0.0495 0.0922 8.25e-06 0.0427 0.0539 0.0395 -0.01 0.999
2 CA1 N' 5 23 {'AM2': 24, 'JS17': 36, 'JS21': 23, 'JS34': 26, 'ZT2': 36} 0.0801 0.0158 6.93e-19 -0.0644 0.0003 0.0221 -0.058 0.0021
2 PFC N 5 15 {'AM2': 15, 'JS17': 25, 'JS21': 22, 'JS34': 55, 'ZT2': 29} 0.587 0.666 1.45e-24 0.0787 0 0.663 0.0759 0
2 PFC N then N' 5 40 {'AM2': 40, 'JS17': 80, 'JS21': 57, 'JS34': 132, 'ZT2': 86} 0.588 0.661 1.49e-115 0.0726 0 0.667 0.0788 0
2 PFC N' 5 25 {'AM2': 25, 'JS17': 56, 'JS21': 35, 'JS34': 79, 'ZT2': 57} 0.657 0.693 1.56e-29 0.0363 0.437 0.661 0.0037 1
4 CA1 N 3 15 {'AM2': 15, 'JS17': 24, 'JS21': 16} 0.124 0.0328 3.32e-14 -0.0913 0.0055 0.0267 -0.0974 0.0014
4 CA1 N then N' 5 9 {'AM2': 38, 'JS17': 28, 'JS21': 20, 'JS34': 9, 'ZT2': 36} 0.0503 0.0456 0.0144 -0.0048 1 0.0496 -0.0008 1
4 CA1 N' 5 19 {'AM2': 25, 'JS17': 24, 'JS21': 20, 'JS34': 19, 'ZT2': 36} 0.0718 0.0338 1e-19 -0.038 0.279 0.0195 -0.0524 0.0082
4 PFC N 5 15 {'AM2': 15, 'JS17': 25, 'JS21': 22, 'JS34': 55, 'ZT2': 29} 0.639 0.717 1.56e-31 0.0777 0 0.709 0.0696 0
4 PFC N then N' 5 40 {'AM2': 40, 'JS17': 81, 'JS21': 57, 'JS34': 132, 'ZT2': 86} 0.642 0.711 1.82e-147 0.0692 0 0.713 0.0712 0
4 PFC N' 5 25 {'AM2': 25, 'JS17': 56, 'JS21': 35, 'JS34': 79, 'ZT2': 57} 0.704 0.748 1.19e-47 0.0445 0.0231 0.705 0.0017 1
8 CA1 N 3 15 {'AM2': 15, 'JS17': 24, 'JS21': 20} 0.136 0.0331 6.72e-13 -0.103 0.0014 0.0206 -0.115 0.0001
8 CA1 N then N' 5 9 {'AM2': 38, 'JS17': 45, 'JS21': 26, 'JS34': 9, 'ZT2': 36} 0.0527 0.0482 0.016 -0.0046 1 0.0248 -0.0279 0.528
8 CA1 N' 5 19 {'AM2': 25, 'JS17': 39, 'JS21': 26, 'JS34': 19, 'ZT2': 38} 0.114 0.00197 6.13e-29 -0.112 0 -0.000128 -0.114 0
8 PFC N 5 15 {'AM2': 15, 'JS17': 25, 'JS21': 22, 'JS34': 55, 'ZT2': 29} 0.684 0.757 4.28e-31 0.0725 0 0.736 0.0519 0.0015
8 PFC N then N' 5 40 {'AM2': 40, 'JS17': 81, 'JS21': 57, 'JS34': 132, 'ZT2': 86} 0.688 0.741 4.95e-170 0.0524 0.0018 0.738 0.0493 0.0046
8 PFC N' 5 25 {'AM2': 25, 'JS17': 56, 'JS21': 35, 'JS34': 79, 'ZT2': 57} 0.738 0.782 1.28e-55 0.0443 0.0074 0.73 -0.0082 1
Result: partially reproduced. The criterion asks for the result under each reading of which novel trials. PFC similarity rises from the first to the last trial under all three readings at 8 cm (N then N' 0.688 to 0.741, Tukey p = 0.0018; N 0.684 to 0.757, p < 0.0001; N' 0.738 to 0.782, p = 0.0074, with the last trial at rank K). CA1 does not change under N then N' (0.053 and 0.048, p = 1.0), but falls under N (0.136 to 0.033, p = 0.0014; three rats, as JS34 and ZT2 have no N pass with 80% of CA1 cells active) and under N' (0.114 to 0.002, p < 0.0001), so the CA1 part holds under one reading of three. Taking each animal's own final kept pass as the last trial gives the same for CA1 (N then N' p = 0.53, N p = 0.0001, N' p < 0.0001); for PFC it holds under N then N' and N (p = 0.0046 and 0.0015) but not N' (p = 1). Under N then N', CA1's kept passes of JS21, JS34 and ZT2 all come from N', so the first trials are N' passes for three of five rats. The 4-cm and 2-cm widths give the same outcome: CA1 falls under N and under N' with at least one of the two definitions of the last trial (p <= 0.0098), and PFC rises under N then N' and N (p < 0.0001), but not under N' at 2 cm (p = 0.44 and 1).
figs4d-reward-location-similarity¶
PV similarity between the environments at the reward wells is the same as at other locations, in both regions.
Per-bin PV similarity (fig2d, 8 cm) between F and the novel session, mean over the five rats, per region; Grubbs' test (two-sided, alpha 0.05) over (a) all 404 bins of a region and (b) the 101 bins of each trajectory. Novel session N' (the panels' pair) or N.
for nm in ("F-N'", "F-N"):
for reg in ("CA1", "PFC"):
M = np.array([pv_bins(d, reg, *PAIRS[nm]) for d in A]); m = np.nanmean(M, 0)
for s_ in (4.0, 2.0):
m_ = np.nanmean([pv_bins(d, reg, *PAIRS[nm], s=s_) for d in A], 0)
print(f" {nm} {reg}, {s_:g} cm: all bins p = {nfn.grubbs(m_)[1]:.3g}; smallest per-trajectory p = {min(nfn.grubbs(m_[k * 101:(k + 1) * 101])[1] for k in range(4)):.3g}")
G, p, i = nfn.grubbs(m)
per = [nfn.grubbs(m[k * 101:(k + 1) * 101]) for k in range(4)]
ends = {nfn.TRACKS[k + 1]: (round(m[k * 101], 2), round(m[k * 101 + 100], 2)) for k in range(4)}
print(f"{nm} {reg}: all bins G = {G:.2f}, p = {p:.3g} (most extreme bin {i}); per trajectory p = {[round(x[1], 3) for x in per]}; "
f"mean {np.nanmean(m):.3f}, range {np.nanmin(m):.3f} to {np.nanmax(m):.3f}; values at the center-well and side-well ends {ends}")
F-N' CA1, 4 cm: all bins p = 1; smallest per-trajectory p = 0.741
F-N' CA1, 2 cm: all bins p = 1; smallest per-trajectory p = 0.292
F-N' CA1: all bins G = 2.65, p = 1 (most extreme bin 365); per trajectory p = [1.0, 1.0, 1.0, 1.0]; mean 0.045, range -0.100 to 0.224; values at the center-well and side-well ends {'OUT_L': (np.float64(0.16), np.float64(-0.03)), 'IN_L': (np.float64(0.04), np.float64(0.19)), 'OUT_R': (np.float64(0.01), np.float64(-0.06)), 'IN_R': (np.float64(0.07), np.float64(0.09))}
F-N' PFC, 4 cm: all bins p = 1; smallest per-trajectory p = 0.186
F-N' PFC, 2 cm: all bins p = 0.00704; smallest per-trajectory p = 0.00117
F-N' PFC: all bins G = 2.86, p = 1 (most extreme bin 366); per trajectory p = [1.0, 1.0, np.float64(0.941), 1.0]; mean 0.836, range 0.723 to 0.901; values at the center-well and side-well ends {'OUT_L': (np.float64(0.82), np.float64(0.8)), 'IN_L': (np.float64(0.8), np.float64(0.89)), 'OUT_R': (np.float64(0.84), np.float64(0.88)), 'IN_R': (np.float64(0.87), np.float64(0.8))}
F-N CA1, 4 cm: all bins p = 1; smallest per-trajectory p = 0.16
F-N CA1, 2 cm: all bins p = 0.753; smallest per-trajectory p = 0.0797
F-N CA1: all bins G = 2.80, p = 1 (most extreme bin 1); per trajectory p = [np.float64(0.859), 1.0, np.float64(0.399), 1.0]; mean 0.062, range -0.132 to 0.351; values at the center-well and side-well ends {'OUT_L': (np.float64(0.35), np.float64(0.08)), 'IN_L': (np.float64(0.05), np.float64(0.17)), 'OUT_R': (np.float64(0.3), np.float64(0.04)), 'IN_R': (np.float64(0.27), np.float64(0.26))}
F-N PFC, 4 cm: all bins p = 1; smallest per-trajectory p = 0.116
F-N PFC, 2 cm: all bins p = 1; smallest per-trajectory p = 0.0661
F-N PFC: all bins G = 2.96, p = 1 (most extreme bin 365); per trajectory p = [1.0, 1.0, 1.0, 1.0]; mean 0.774, range 0.617 to 0.889; values at the center-well and side-well ends {'OUT_L': (np.float64(0.79), np.float64(0.76)), 'IN_L': (np.float64(0.78), np.float64(0.81)), 'OUT_R': (np.float64(0.85), np.float64(0.83)), 'IN_R': (np.float64(0.81), np.float64(0.8))}
Result: blocked, methods underspecified (the smoothing width). At 8 cm Grubbs' test finds no outlier among the bins of either region, for F-N' and F-N, over all bins (p = 1 in every case) or per trajectory (smallest p 0.40), and the values at the well ends lie inside the range of the other bins (CA1 F-N' -0.100 to 0.224, PFC 0.723 to 0.901). The 4-cm width gives the same (smallest p 0.116). With 2 cm, Grubbs' test finds an outlier among the PFC F-N' bins (all bins p = 0.0070, per trajectory p = 0.0012), so the claim would not reproduce. The outcome changes with the width, which the paper's text does not fix, and the claims file's rule for that case (fig2d-pv-similarity) makes the claim methods underspecified.
Decoding (Figure 2E, 2F)¶
fig2e-ca1-error-familiar and fig2e-pfc-error-stable¶
Decoding with N' rate maps, CA1 position error is higher in F than in N, and in F no better than with cell-ID shuffled templates.
Decoding with N' rate maps, PFC position error is similar in N and F, and in both well above chance (below the cell-ID shuffles).
Decoder of the Methods (nfn.decode_session): templates are the N' maps (8 cm) of cells included in N' with a peak of at least 3 Hz; 200-ms windows every 150 ms (the Methods' "50 ms overlap"); bins on one trajectory, mean speed above 5 cm/s, more than 15 cm from a well, at least one cell firing; error along the four trajectories laid end to end. Shuffles permute the cell identities of the templates: (a) one shuffle, (b) the median over 100. Wilcoxon paired test of each bin against its shuffle; Kruskal-Wallis with Dunn's test over (a) the four groups, (b) N and F per region. Variants: a 50-ms step (the code), 4 and 2 cm smoothing, and the decode restricted to the current trajectory. The theta restriction of the code cannot be run (no theta times in the dandiset).
def decode_all(s=8.0, seed=401, **kw):
rng = np.random.default_rng(seed); out = []
for d in A:
for reg in ("CA1", "PFC"):
cells = nfn.included(d, reg, 2); m = MAPS[(d["animal"], s, 2)]
cells = cells[nfn.peak_rate({t: (m[t][0], m[t][1][cells]) for t in m}) >= nfn.PEAK_HZ]
tm = {t: (m[t][0], m[t][1][cells]) for t in m}
for ep in (0, 1):
df, _ = nfn.decode_session(d, ep, cells, tm, rng=rng, **kw)
df["animal"], df["region"], df["env"], df["cells"] = d["animal"], reg, nfn.EPOCHS[ep], len(cells)
out.append(df)
return pd.concat(out, ignore_index=True)
def decode_stats(D):
g = {(r, e): D[(D.region == r) & (D.env == e)] for r in ("CA1", "PFC") for e in ("N", "F")}
a4, dn4 = nfn.kruskal_dunn([g[k].error for k in g], [f"{r} {e}" for r, e in g])
out = {"KW4 p": a4["p"]}
for r in ("CA1", "PFC"):
out[f"{r} N median"] = g[(r, "N")].error.median(); out[f"{r} F median"] = g[(r, "F")].error.median()
out[f"{r} Dunn4 N-F p"] = dn4.loc[f"{r} N", f"{r} F"]
out[f"{r} KW2 N-F p"] = stats.kruskal(g[(r, "N")].error, g[(r, "F")].error).pvalue
for e in ("N", "F"):
x = g[(r, e)]
out[f"{r} {e} shuffle median"] = x.error_shuffle_median.median()
out[f"{r} {e} W(a) p"] = stats.wilcoxon(x.error, x.error_shuffle1).pvalue
out[f"{r} {e} W(b) p"] = stats.wilcoxon(x.error, x.error_shuffle_median).pvalue
out[f"{r} {e} mean error minus shuffle(b)"] = (x.error - x.error_shuffle_median).mean()
return out
DEC = decode_all()
print(DEC.groupby(["region", "env", "animal"]).agg(bins=("error", "size"), cells=("cells", "first"), median_error=("error", "median"),
shuffle_median=("error_shuffle_median", "median")).round(1))
VAR = {"main (Methods, 8 cm)": decode_stats(DEC),
"step 50 ms (code)": decode_stats(decode_all(step=0.05)),
"smoothing 4 cm": decode_stats(decode_all(s=4.0)),
"smoothing 2 cm": decode_stats(decode_all(s=2.0)),
"decode restricted to the current trajectory": decode_stats(decode_all(restrict_track=True))}
VARDF = pd.DataFrame(VAR)
print("\nError in cm; p-values of the tests:")
with pd.option_context("display.float_format", lambda v: f"{v:.3g}"):
print(VARDF)
print("\nError range (cm), main:", DEC.error.min(), "to", DEC.error.max())
bins cells median_error shuffle_median
region env animal
CA1 F AM2 1138 23 274.5 228.9
JS17 1227 15 179.8 212.0
JS21 621 16 278.6 254.0
JS34 1882 38 243.1 235.1
ZT2 1423 32 222.3 199.0
N AM2 804 23 171.6 205.0
JS17 1016 15 164.2 207.8
JS21 1029 16 210.5 242.9
JS34 1866 38 130.2 224.0
ZT2 1490 32 177.7 218.0
PFC F AM2 1233 23 171.3 226.0
JS17 1443 48 95.3 227.0
JS21 1147 20 202.6 249.5
JS34 2006 10 225.2 198.0
ZT2 1564 39 169.5 203.0
N AM2 856 23 87.3 233.6
JS17 1112 48 137.5 218.5
JS21 1218 20 185.9 239.8
JS34 1965 10 171.8 187.7
ZT2 1586 39 158.5 209.0
Error in cm; p-values of the tests:
main (Methods, 8 cm) step 50 ms (code) smoothing 4 cm smoothing 2 cm decode restricted to the current trajectory
KW4 p 1.3e-148 0 1.03e-115 3.54e-94 0
CA1 N median 172 172 175 182 15
CA1 F median 230 230 216 212 56.8
CA1 Dunn4 N-F p 1.83e-140 0 2.97e-109 1.93e-88 0
CA1 KW2 N-F p 1.04e-139 0 1.19e-107 9.33e-87 0
CA1 N shuffle median 221 221 220 224 48.8
CA1 N W(a) p 6.26e-94 1.3e-282 1.94e-80 4.96e-72 0
CA1 N W(b) p 1.12e-129 0 8.51e-121 3.39e-114 0
CA1 N mean error minus shuffle(b) -56.3 -56.9 -53 -52.2 -21.8
CA1 F shuffle median 226 226 223 222 54.1
CA1 F W(a) p 2.36e-08 1.84e-23 3.44e-10 1.78e-08 0.0841
CA1 F W(b) p 0.0421 0.000562 0.793 0.00673 2.48e-09
CA1 F mean error minus shuffle(b) 8.75 8.48 4.18 -3.49 4.86
PFC N median 149 148 154 165 28.5
PFC F median 186 187 188 193 23.5
PFC Dunn4 N-F p 7.04e-06 1.4e-17 9.64e-06 1.56e-06 6.19e-05
PFC KW2 N-F p 3.85e-06 1.08e-16 4.7e-06 5.73e-07 1.36e-05
PFC N shuffle median 218 218 230 231 62
PFC N W(a) p 1.27e-91 1.05e-242 3.52e-136 3.29e-46 8.91e-294
PFC N W(b) p 6.43e-67 1.13e-192 4.09e-91 3.76e-83 2.61e-299
PFC N mean error minus shuffle(b) -41.6 -41.2 -48.6 -46.4 -24.8
PFC F shuffle median 223 223 225 212 60.7
PFC F W(a) p 7.76e-32 8.38e-93 1.51e-49 3.74e-21 6.42e-276
PFC F W(b) p 2.39e-35 6.55e-104 4.45e-50 2.41e-17 0
PFC F mean error minus shuffle(b) -27.9 -27.7 -31.6 -16 -27.2
Error range (cm), main: 0.0011821387079180568 to 784.208562520374
v = VARDF["main (Methods, 8 cm)"]
print(f"CA1: F median {v['CA1 F median']:.1f} cm against N {v['CA1 N median']:.1f}; Dunn (4 groups) p = {v['CA1 Dunn4 N-F p']:.2g}, "
f"(N and F only) p = {v['CA1 KW2 N-F p']:.2g}")
print(f"CA1 F against its shuffle: median {v['CA1 F median']:.1f} against {v['CA1 F shuffle median']:.1f} cm; "
f"mean paired difference {v['CA1 F mean error minus shuffle(b)']:.1f} cm; Wilcoxon (a) p = {v['CA1 F W(a) p']:.3g}, (b) p = {v['CA1 F W(b) p']:.3g}")
for name in VARDF.columns[1:]:
w = VARDF[name]
print(f" {name}: CA1 F > N {w['CA1 F median'] > w['CA1 N median']} (p {w['CA1 Dunn4 N-F p']:.2g}); "
f"CA1 F against shuffle p (a) {w['CA1 F W(a) p']:.3g}, (b) {w['CA1 F W(b) p']:.3g}, mean difference {w['CA1 F mean error minus shuffle(b)']:.1f} cm")
CA1: F median 230.1 cm against N 172.4; Dunn (4 groups) p = 1.8e-140, (N and F only) p = 1e-139 CA1 F against its shuffle: median 230.1 against 226.2 cm; mean paired difference 8.8 cm; Wilcoxon (a) p = 2.36e-08, (b) p = 0.0421 step 50 ms (code): CA1 F > N True (p 0); CA1 F against shuffle p (a) 1.84e-23, (b) 0.000562, mean difference 8.5 cm smoothing 4 cm: CA1 F > N True (p 3e-109); CA1 F against shuffle p (a) 3.44e-10, (b) 0.793, mean difference 4.2 cm smoothing 2 cm: CA1 F > N True (p 1.9e-88); CA1 F against shuffle p (a) 1.78e-08, (b) 0.00673, mean difference -3.5 cm decode restricted to the current trajectory: CA1 F > N True (p 0); CA1 F against shuffle p (a) 0.0841, (b) 2.48e-09, mean difference 4.9 cm
Decoder check¶
A check that the decoder works: N' is decoded with templates from the other half of the N' passes (alternate passes of each trajectory type), and for every decoded bin the decoded trajectory is compared with the current one, and with the current arm of the maze (OUT_L and IN_L run over the left arm, OUT_R and IN_R over the right arm).
# Decoder check: cross-validated decoding of N' (templates from half of the N' trials, decoding
# the other half), and the fraction of bins decoded onto the current trajectory and onto the
# current arm (OUT_L and IN_L share the left arm, OUT_R and IN_R the right arm).
ARM = {1: "L", 2: "L", 3: "R", 4: "R"}
def cv_decode(s=8.0, seed=404, **kw):
rng = np.random.default_rng(seed); out = []
for d in A:
tbt = nfn.trials_of_track(d, 2) # alternate passes of each trajectory type go to the two halves
folds = tuple(np.concatenate([np.array(tbt[t][k::2], int) for t in nfn.TRACKS]) for k in (0, 1))
for reg in ("CA1", "PFC"):
cells = nfn.included(d, reg, 2); m = MAPS[(d["animal"], s, 2)]
cells = cells[nfn.peak_rate({t: (m[t][0], m[t][1][cells]) for t in m}) >= nfn.PEAK_HZ]
for k in (0, 1):
mt = nfn.rate_maps(d, 2, cells, smooth_cm=s, trial=folds[1 - k])
df, _ = nfn.decode_session(d, 2, cells, mt, rng=rng, n_shuffles=1, trials=folds[k], **kw)
df["animal"], df["region"], df["env"] = d["animal"], reg, "N' (cross-validated)"
out.append(df)
return pd.concat(out, ignore_index=True)
CV = cv_decode()
CVR = cv_decode(restrict_track=True)
allD = pd.concat([DEC, CV], ignore_index=True)
allD["on_track"] = allD.decoded_track == allD.track
allD["on_arm"] = allD.decoded_track.map(ARM) == allD.track.map(ARM)
summ = allD.groupby(["region", "env"]).agg(bins=("error", "size"), median_error=("error", "median"),
on_trajectory=("on_track", "mean"), on_arm=("on_arm", "mean"))
onr = allD[allD.on_track].groupby(["region", "env"]).error.median().rename("median error when on the trajectory")
summ = summ.join(onr)
summ["N' cross-validated, restricted decoder median error"] = [CVR[CVR.region == r].error.median() if e.startswith("N'") else np.nan for r, e in summ.index]
print("Main decoder (Methods, 8 cm); chance on trajectory 0.25, on arm 0.5:")
print(summ.round(3).to_string())
print("\nPer animal, fraction on the current trajectory:")
print(allD.groupby(["region", "env", "animal"]).on_track.mean().unstack("animal").round(2).to_string())
Main decoder (Methods, 8 cm); chance on trajectory 0.25, on arm 0.5:
bins median_error on_trajectory on_arm median error when on the trajectory N' cross-validated, restricted decoder median error
region env
CA1 F 6291 230.072 0.249 0.481 44.984 NaN
N 6205 172.444 0.427 0.652 10.905 NaN
N' (cross-validated) 8109 12.193 0.640 0.775 5.367 6.758
PFC F 7393 186.075 0.375 0.541 16.351 NaN
N 6737 149.303 0.423 0.594 21.081 NaN
N' (cross-validated) 8431 73.839 0.517 0.647 8.822 12.574
Per animal, fraction on the current trajectory:
animal AM2 JS17 JS21 JS34 ZT2
region env
CA1 F 0.23 0.37 0.20 0.18 0.27
N 0.43 0.44 0.38 0.46 0.41
N' (cross-validated) 0.54 0.55 0.50 0.83 0.67
PFC F 0.38 0.47 0.34 0.31 0.39
N 0.50 0.45 0.41 0.35 0.46
N' (cross-validated) 0.51 0.55 0.44 0.36 0.69
Result: partially reproduced. CA1 error is higher in F than in N (median 230 against 172 cm, Dunn p = 1.8e-140), in every variant. In F the CA1 error is not below the shuffle; it is above it (mean paired difference 8.8 cm), and the Wilcoxon test finds the difference under both shuffle readings (p = 2.4e-08 and 0.042), so the claim of no difference fails by the paper's test. The other widths give the same outcome: 4 cm differs under reading (a) (p = 3.4e-10; reading (b) p = 0.79), and 2 cm puts F slightly below the shuffle (-3.5 cm, p = 1.8e-08 and 0.0067). The decoder works: decoding N' with templates from the other half of its passes gives a median error of 12 cm (64% of bins on the current trajectory, against 0.25 by chance). In N, 43% of CA1 bins are decoded onto the current trajectory and 65% onto the current arm; the large median errors come from decodes onto the other trajectory of the same arm, which the criterion counts as far.
v = VARDF["main (Methods, 8 cm)"]
print(f"PFC: N median {v['PFC N median']:.1f}, F median {v['PFC F median']:.1f} cm; Dunn (4 groups) p = {v['PFC Dunn4 N-F p']:.2g}, "
f"(N and F only) p = {v['PFC KW2 N-F p']:.2g}")
for e in ("N", "F"):
print(f"PFC {e} against shuffle: median {v[f'PFC {e} median']:.1f} against {v[f'PFC {e} shuffle median']:.1f}; mean difference "
f"{v[f'PFC {e} mean error minus shuffle(b)']:.1f} cm; Wilcoxon (a) p = {v[f'PFC {e} W(a) p']:.2g}, (b) p = {v[f'PFC {e} W(b) p']:.2g}")
for name in VARDF.columns[1:]:
w = VARDF[name]
print(f" {name}: PFC N {w['PFC N median']:.1f}, F {w['PFC F median']:.1f} cm, Dunn p {w['PFC Dunn4 N-F p']:.2g}; "
f"below shuffle in N and F: {w['PFC N mean error minus shuffle(b)'] < 0 and w['PFC F mean error minus shuffle(b)'] < 0} "
f"(max p {max(w['PFC N W(a) p'], w['PFC N W(b) p'], w['PFC F W(a) p'], w['PFC F W(b) p']):.2g})")
print("\nPer animal, PFC median error N and F (cm):")
print(DEC[DEC.region == "PFC"].groupby(["animal", "env"]).error.median().unstack().round(1))
PFC: N median 149.3, F median 186.1 cm; Dunn (4 groups) p = 7e-06, (N and F only) p = 3.9e-06 PFC N against shuffle: median 149.3 against 218.5; mean difference -41.6 cm; Wilcoxon (a) p = 1.3e-91, (b) p = 6.4e-67 PFC F against shuffle: median 186.1 against 222.9; mean difference -27.9 cm; Wilcoxon (a) p = 7.8e-32, (b) p = 2.4e-35 step 50 ms (code): PFC N 148.3, F 186.7 cm, Dunn p 1.4e-17; below shuffle in N and F: True (max p 8.4e-93) smoothing 4 cm: PFC N 154.0, F 188.1 cm, Dunn p 9.6e-06; below shuffle in N and F: True (max p 1.5e-49) smoothing 2 cm: PFC N 165.1, F 192.7 cm, Dunn p 1.6e-06; below shuffle in N and F: True (max p 2.4e-17) decode restricted to the current trajectory: PFC N 28.5, F 23.5 cm, Dunn p 6.2e-05; below shuffle in N and F: True (max p 6.4e-276) Per animal, PFC median error N and F (cm): env F N animal AM2 171.3 87.3 JS17 95.3 137.5 JS21 202.6 185.9 JS34 225.2 171.8 ZT2 169.5 158.5
Result: partially reproduced. PFC error is below the shuffle in N (149 against 218 cm, p = 1.3e-91 and 6.4e-67) and in F (186 against 223 cm, p = 7.8e-32 and 2.4e-35). It is not similar in N and F: F is higher (Dunn p = 7.0e-06; N and F alone, p = 3.9e-06), in every variant (p 6.2e-05 to 1.4e-17). With the decode restricted to the current trajectory, F is lower than N (23.5 against 28.5 cm, p = 6.2e-05), so the sign of the difference depends on the decoder, but every variant finds one. Four of five rats have higher PFC median error in F (JS17 is lower). PFC bins are decoded onto the current trajectory in 42% of N bins and 38% of F bins (cross-validated N', 52%).
fig2f-switch-dynamics¶
After the switch from N to F, CA1 decoding error rises within the first trial, while PFC decoding error stays stable.
Bins of the decoder above, grouped by trial (a row of the trials table; bins pooled over rats). Readings of trial -1: (a) the last N trial, so that N trials -10 to -1 are the last ten; (b) the code's batch script as it runs, which takes the first ten N trials. Trials 1 to 10 are the first ten F trials. Kruskal-Wallis over the 20 trials with Dunn's test (Bonferroni over all pairs). Variant decoder of the code: 200-ms bins without overlap, no speed filter, at least 2 active cells, decode restricted to the current trajectory.
def by_trial(D):
D = D.copy()
D["k"] = -1
for d in A:
tr = d["trials"]
for ep in (0, 1):
ids = list(tr[tr.epoch == ep].index); pos = {t: i for i, t in enumerate(ids)}
m = (D.animal == d["animal"]) & (D.env == nfn.EPOCHS[ep])
D.loc[m, "k"] = [pos.get(t, -999) for t in D.loc[m, "trial"]]
D.loc[m, "nk"] = len(ids)
return D
def trial_labels(D, reading):
"""Trial label: F trials 1, 2, ... from the switch; N trials (a) -1 for the last N trial and
back from it, (b) the code's first ten N trials labelled -10 to -1 (later N trials dropped)."""
if reading == "a":
n_lab = D.k - D.nk
else:
n_lab = np.where(D.k < 10, D.k - 10, -999)
return np.where(D.env == "N", n_lab, D.k + 1)
def trial_tests(D, reading):
D = D.copy()
D["lab"] = trial_labels(D, reading)
D = D[(D.lab >= -10) & (D.lab <= 10) & (D.lab != 0)]
out = {}
for r in ("CA1", "PFC"):
x = D[D.region == r]; labs = sorted(x.lab.unique())
a, dn = nfn.kruskal_dunn([x[x.lab == l].error for l in labs], labs)
med = x.groupby("lab").error.median()
out[r] = {"KW p": a["p"], "trial -1 median": med[-1], "trial 1 median": med[1], "n -1": int((x.lab == -1).sum()),
"n 1": int((x.lab == 1).sum()), "n trials -10 to -1": int((x.lab < 0).sum()), "Dunn -1 vs 1 p": dn.loc[-1, 1]}
return pd.DataFrame(out)
TR_M = by_trial(decode_all(seed=402, n_shuffles=1))
TR_C = by_trial(decode_all(seed=403, n_shuffles=1, step=0.2, min_speed=-1, min_active=2, restrict_track=True))
TR_4 = by_trial(decode_all(s=4.0, seed=404, n_shuffles=1))
TR_2 = by_trial(decode_all(s=2.0, seed=405, n_shuffles=1))
for nm, D in (("Methods decoder, 8 cm", TR_M), ("Methods decoder, 4 cm", TR_4), ("Methods decoder, 2 cm", TR_2), ("code decoder (variant), 8 cm", TR_C)):
for rd in ("a", "b"):
print(f"--- {nm}, reading ({rd}) ---")
with pd.option_context("display.float_format", lambda v: f"{v:.3g}"):
print(trial_tests(D, rd))
for rd in ("a", "b"):
x = TR_M.copy(); x["lab"] = trial_labels(x, rd)
print(f"\nMethods decoder, 8 cm, reading ({rd}): median error by trial (cm)")
print(x[(x.lab >= -10) & (x.lab <= 10) & (x.lab != 0)].groupby(["region", "lab"]).error.median().unstack().round(0).to_string())
print("\nN trials per animal (reading b uses the first ten):", {d["animal"]: int((d["trials"].epoch == 0).sum()) for d in A})
--- Methods decoder, 8 cm, reading (a) ---
CA1 PFC
KW p 2.59e-67 8.16e-43
trial -1 median 182 149
trial 1 median 223 200
n -1 258 281
n 1 191 269
n trials -10 to -1 1.64e+03 1.73e+03
Dunn -1 vs 1 p 0.00206 0.124
--- Methods decoder, 8 cm, reading (b) ---
CA1 PFC
KW p 4.97e-45 1.4e-28
trial -1 median 259 141
trial 1 median 223 200
n -1 114 124
n 1 191 269
n trials -10 to -1 1.54e+03 1.76e+03
Dunn -1 vs 1 p 1 1
--- Methods decoder, 4 cm, reading (a) ---
CA1 PFC
KW p 7.5e-62 1.2e-52
trial -1 median 184 150
trial 1 median 192 229
n -1 258 281
n 1 198 269
n trials -10 to -1 1.64e+03 1.73e+03
Dunn -1 vs 1 p 0.00793 0.000397
--- Methods decoder, 4 cm, reading (b) ---
CA1 PFC
KW p 1.34e-44 3.43e-26
trial -1 median 243 218
trial 1 median 192 229
n -1 114 124
n 1 198 269
n trials -10 to -1 1.54e+03 1.76e+03
Dunn -1 vs 1 p 1 1
--- Methods decoder, 2 cm, reading (a) ---
CA1 PFC
KW p 7.06e-45 6.13e-43
trial -1 median 188 162
trial 1 median 168 214
n -1 258 281
n 1 199 269
n trials -10 to -1 1.65e+03 1.73e+03
Dunn -1 vs 1 p 0.302 0.0763
--- Methods decoder, 2 cm, reading (b) ---
CA1 PFC
KW p 4.26e-35 3.42e-19
trial -1 median 240 154
trial 1 median 168 214
n -1 114 124
n 1 199 269
n trials -10 to -1 1.56e+03 1.76e+03
Dunn -1 vs 1 p 1 1
--- code decoder (variant), 8 cm, reading (a) ---
CA1 PFC
KW p 3.9e-109 1.61e-81
trial -1 median 40.4 70.5
trial 1 median 46.1 48.2
n -1 378 627
n 1 105 230
n trials -10 to -1 1.39e+03 2.02e+03
Dunn -1 vs 1 p 1 1
--- code decoder (variant), 8 cm, reading (b) ---
CA1 PFC
KW p 8.51e-50 8.59e-33
trial -1 median 20.3 18.4
trial 1 median 46.1 48.2
n -1 83 110
n 1 105 230
n trials -10 to -1 1e+03 1.59e+03
Dunn -1 vs 1 p 0.000266 0.00126
Methods decoder, 8 cm, reading (a): median error by trial (cm)
lab -10.0 -9.0 -8.0 -7.0 -6.0 -5.0 -4.0 -3.0 -2.0 -1.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
region
CA1 133.0 20.0 111.0 72.0 159.0 187.0 107.0 186.0 55.0 182.0 223.0 192.0 225.0 304.0 195.0 286.0 173.0 177.0 298.0 341.0
PFC 270.0 70.0 99.0 74.0 136.0 62.0 49.0 69.0 64.0 149.0 200.0 197.0 211.0 385.0 182.0 194.0 92.0 230.0 321.0 265.0
Methods decoder, 8 cm, reading (b): median error by trial (cm)
lab -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 1 2 3 4 5 6 7 8 9 10
region
CA1 173.0 210.0 228.0 192.0 170.0 203.0 201.0 178.0 40.0 259.0 223.0 192.0 225.0 304.0 195.0 286.0 173.0 177.0 298.0 341.0
PFC 110.0 133.0 234.0 194.0 288.0 140.0 153.0 154.0 166.0 141.0 200.0 197.0 211.0 385.0 182.0 194.0 92.0 230.0 321.0 265.0
N trials per animal (reading b uses the first ten): {'AM2': 25, 'JS17': 33, 'JS21': 30, 'ZT2': 42, 'JS34': 69}
Result: partially reproduced. The criterion asks for the result under each reading of trial -1. Under reading (a), trial -1 the last N trial, CA1 error rises from 182 to 223 cm between trial -1 and trial 1 (Dunn p = 0.0021), and PFC does not differ (149 and 200 cm, p = 0.12). Under reading (b), the code's first ten N trials labelled -10 to -1 (later N trials dropped), CA1 trial 1 is not above trial -1 (223 against 259 cm, p = 1), and PFC does not differ (p = 1). So the PFC part holds under both readings and the CA1 part under one. The other widths do not give more: at 4 cm CA1 rises under (a) (p = 0.0079) but PFC rises too (p = 0.0004), and at 2 cm CA1 does not rise under either reading (p = 0.30 and 1) while PFC is stable (p = 0.076 and 1). The code's decoder (variant) finds no CA1 rise under (a) (p = 1); under (b) both regions rise (p = 0.00027 and 0.0013). Calling the outcome partial at every width is a judgement: at 8 and 4 cm one part holds under each reading or one reading holds in full, but at 2 cm only the PFC claim of absence holds and CA1 rises under neither reading, which alone would be closer to not reproduced. Since the outcome is partial or worse at every width, and reproduced at none, the width rule does not block it.
Reward ramping (Figure S4)¶
figs4-ramping-fraction¶
Ramping cells make up 16.4% of CA1 cell samples (81 of 495) and 18.1% of PFC cell samples (115 of 636), a cell sample being one cell on one trajectory type.
A session map (8 cm) of a cell on a trajectory ramps when the regression of rate on the position traversed toward the reward, over the whole trajectory, has Pearson r > 0.5 and p < 0.05, its peak is at least 3 Hz and lies within 20 cm of the trajectory's end (Methods). Denominators: (a) cell samples of the three sessions pooled, (b) N' only, each over included cells with a valid map (at least 3 visited bins). Variants: 4 and 2 cm smoothing, and denominators restricted to tuned samples (peak at least 3 Hz).
def ramp_table(s):
rows = []
for d in A:
for reg in ("CA1", "PFC"):
for ep in range(3):
m = MAPS[(d["animal"], s, ep)]
for c in nfn.included(d, reg, ep):
for tr in nfn.TRACKS:
ok, rr, p, valid = nfn.ramping(m, c, tr)
pk = np.nanmax(m[tr][1][c]) if np.isfinite(m[tr][1][c]).any() else np.nan
rows.append({"animal": d["animal"], "region": reg, "ep": ep, "unit": c, "track": tr, "ramp": ok,
"valid": valid, "r": rr, "tuned": pk >= 3})
return pd.DataFrame(rows)
RAMP = {s: ramp_table(s) for s in SMOOTHS}
for s in SMOOTHS:
R = RAMP[s]
for reg in ("CA1", "PFC"):
for nm, eps in (("(a) all sessions", [0, 1, 2]), ("(b) N'", [2])):
x = R[(R.region == reg) & R.ep.isin(eps) & R.valid]
xt = x[x.tuned]
print(f"smoothing {s:g} cm, {reg} {nm}: {x.ramp.sum()} of {len(x)} = {100 * x.ramp.mean():.1f}%; "
f"of tuned samples {xt.ramp.sum()} of {len(xt)} = {100 * xt.ramp.mean():.1f}%")
R = RAMP[8.0]
x = R[(R.ep == 2) & R.valid]
print("\nN', 8 cm, ramping samples per animal and region:")
print(x.groupby(["region", "animal"]).agg(samples=("ramp", "size"), ramping=("ramp", "sum")).unstack("region"))
print("distribution of r among tuned samples (N', 8 cm):", np.nanpercentile(x[x.tuned].r, [10, 50, 90]).round(2))
smoothing 8 cm, CA1 (a) all sessions: 49 of 1652 = 3.0%; of tuned samples 49 of 813 = 6.0%
smoothing 8 cm, CA1 (b) N': 22 of 596 = 3.7%; of tuned samples 22 of 291 = 7.6%
smoothing 8 cm, PFC (a) all sessions: 67 of 2192 = 3.1%; of tuned samples 67 of 1430 = 4.7%
smoothing 8 cm, PFC (b) N': 22 of 728 = 3.0%; of tuned samples 22 of 484 = 4.5%
smoothing 4 cm, CA1 (a) all sessions: 55 of 1652 = 3.3%; of tuned samples 55 of 944 = 5.8%
smoothing 4 cm, CA1 (b) N': 21 of 596 = 3.5%; of tuned samples 21 of 335 = 6.3%
smoothing 4 cm, PFC (a) all sessions: 109 of 2192 = 5.0%; of tuned samples 109 of 1537 = 7.1%
smoothing 4 cm, PFC (b) N': 35 of 728 = 4.8%; of tuned samples 35 of 519 = 6.7%
smoothing 2 cm, CA1 (a) all sessions: 53 of 1652 = 3.2%; of tuned samples 53 of 1080 = 4.9%
smoothing 2 cm, CA1 (b) N': 22 of 596 = 3.7%; of tuned samples 22 of 376 = 5.9%
smoothing 2 cm, PFC (a) all sessions: 92 of 2192 = 4.2%; of tuned samples 92 of 1710 = 5.4%
smoothing 2 cm, PFC (b) N': 37 of 728 = 5.1%; of tuned samples 37 of 572 = 6.5%
N', 8 cm, ramping samples per animal and region:
samples ramping
region CA1 PFC CA1 PFC
animal
AM2 116 100 9 5
JS17 92 252 5 4
JS21 72 120 0 6
JS34 172 48 3 0
ZT2 144 208 5 7
distribution of r among tuned samples (N', 8 cm): [-0.7 -0.03 0.65]
Result: not reproduced. Ramping samples are 3.0% of CA1 (49 of 1652) and 3.1% of PFC (67 of 2192) cell samples over all sessions, and 3.7% (22 of 596) and 3.0% (22 of 728) in N', against 16.4% and 18.1%. No smoothing width or denominator gives more than 7.6%. The N' denominators (596 and 728) are above the stated 495 and 636 beyond tolerance.
figs4-ramping-generalization¶
PFC ramping cells generalize their reward representation across environments, whereas CA1 ramping cells remap, and the generalization of ramping cells is similar to that of other cells in each region.
Ramping cell samples identified in N' (8 cm). Per animal and trajectory, the PV over its ramping samples and the PV over the other tuned cells of the PV analysis on that trajectory, per bin (101 normalized bins), for N-N' and F-N'; a PV needs at least 3 cells (an assumption: a correlation over 2 cells is always plus or minus 1). Bins pooled over rats and trajectories; Kruskal-Wallis over the 8 groups with Dunn's test. Variant: the ramping samples pooled over animals into one population per trajectory.
RN = RAMP[8.0]; RN = RN[(RN.ep == 2) & RN.ramp]
pos101 = np.linspace(0, 1, 101)
def track_norm(m, track, cells):
c, r = m[track]; r = r[cells]
out = nfn.normalized({t: (m[t][0], m[t][1][cells]) for t in m})
k = list(nfn.TRACKS).index(track)
return out[:, k * 101:(k + 1) * 101]
def bin_corr(X, Y):
X, Y = np.nan_to_num(X), np.nan_to_num(Y)
return np.array([np.corrcoef(X[:, b], Y[:, b])[0, 1] if X[:, b].std() > 0 and Y[:, b].std() > 0 else np.nan for b in range(X.shape[1])])
G8 = {}; counts = []
pooled = {}
for reg in ("CA1", "PFC"):
for pair in ("N-N'", "F-N'"):
e1, e2 = PAIRS[pair]
for kind in ("ramping", "other"):
G8[(reg, kind, pair)] = []
for d in A:
cells_pv, _, _ = pv_cells(d, reg, e1, e2)
m1, m2 = MAPS[(d["animal"], 8.0, e1)], MAPS[(d["animal"], 8.0, e2)]
for track in nfn.TRACKS:
ramp = RN[(RN.animal == d["animal"]) & (RN.region == reg) & (RN.track == track)].unit.values
other = np.setdiff1d(cells_pv, ramp)
counts.append({"region": reg, "pair": pair, "animal": d["animal"], "track": track, "ramping": len(ramp), "other": len(other)})
for kind, cs in (("ramping", ramp), ("other", other)):
X, Y = track_norm(m1, track, cs), track_norm(m2, track, cs)
pooled.setdefault((reg, kind, pair, track), []).append((X, Y))
if len(cs) >= 3:
G8[(reg, kind, pair)].append(bin_corr(X, Y))
CNT = pd.DataFrame(counts)
print("ramping samples (N') per region and trajectory, summed over animals:",
CNT[CNT.pair == "N-N'"].groupby(["region", "track"]).ramping.sum().to_dict())
print("animal-trajectory populations with at least 3 ramping cells:",
CNT[(CNT.pair == "N-N'") & (CNT.ramping >= 3)].groupby("region").size().to_dict())
names = [f"{r} {k} {p}" for r in ("CA1", "PFC") for k in ("ramping", "other") for p in ("N-N'", "F-N'")]
vals = [np.concatenate(G8[(r, k, p)]) if G8[(r, k, p)] else np.array([]) for r in ("CA1", "PFC") for k in ("ramping", "other") for p in ("N-N'", "F-N'")]
print("bins per group:", dict(zip(names, [int(np.isfinite(v).sum()) for v in vals])))
print("medians:", dict(zip(names, [round(float(np.nanmedian(v)), 3) if len(v) else None for v in vals])))
def report(vals, names):
keep = [i for i, v in enumerate(vals) if np.isfinite(v).sum() > 0]
a, dn = nfn.kruskal_dunn([vals[i] for i in keep], [names[i] for i in keep])
print(f"Kruskal-Wallis H = {a['H']:.1f}, p = {a['p']:.2g}")
for r in ("CA1", "PFC"):
pairs = [(f"{r} ramping N-N'", f"{r} ramping F-N'"), (f"{r} ramping N-N'", f"{r} other N-N'"), (f"{r} ramping F-N'", f"{r} other F-N'")]
for a_, b_ in pairs:
if a_ in dn.index and b_ in dn.index:
print(f" Dunn {a_} vs {b_}: p = {dn.loc[a_, b_]:.3g}")
else:
print(f" Dunn {a_} vs {b_}: not testable (no population of 3 ramping cells)")
report(vals, names)
print("\nVariant, ramping and other samples pooled over animals per trajectory:")
vals2 = []
for r in ("CA1", "PFC"):
for k in ("ramping", "other"):
for p in ("N-N'", "F-N'"):
v = []
for track in nfn.TRACKS:
X = np.vstack([x for x, y in pooled[(r, k, p, track)]]); Y = np.vstack([y for x, y in pooled[(r, k, p, track)]])
if len(X) >= 3:
v.append(bin_corr(X, Y))
vals2.append(np.concatenate(v) if v else np.array([]))
print("medians:", dict(zip(names, [round(float(np.nanmedian(v)), 3) if len(v) else None for v in vals2])))
report(vals2, names)
ramping samples (N') per region and trajectory, summed over animals: {('CA1', 1): 6, ('CA1', 2): 4, ('CA1', 3): 5, ('CA1', 4): 7, ('PFC', 1): 5, ('PFC', 2): 5, ('PFC', 3): 7, ('PFC', 4): 5}
animal-trajectory populations with at least 3 ramping cells: {'CA1': 1, 'PFC': 1}
bins per group: {"CA1 ramping N-N'": 101, "CA1 ramping F-N'": 101, "CA1 other N-N'": 2020, "CA1 other F-N'": 2020, "PFC ramping N-N'": 101, "PFC ramping F-N'": 101, "PFC other N-N'": 2020, "PFC other F-N'": 2020}
medians: {"CA1 ramping N-N'": 0.498, "CA1 ramping F-N'": -0.303, "CA1 other N-N'": 0.562, "CA1 other F-N'": -0.01, "PFC ramping N-N'": 0.92, "PFC ramping F-N'": 0.846, "PFC other N-N'": 0.898, "PFC other F-N'": 0.875}
Kruskal-Wallis H = 5387.8, p = 0
Dunn CA1 ramping N-N' vs CA1 ramping F-N': p = 1.65e-22
Dunn CA1 ramping N-N' vs CA1 other N-N': p = 1
Dunn CA1 ramping F-N' vs CA1 other F-N': p = 0.014
Dunn PFC ramping N-N' vs PFC ramping F-N': p = 0.338
Dunn PFC ramping N-N' vs PFC other N-N': p = 1
Dunn PFC ramping F-N' vs PFC other F-N': p = 0.114
Variant, ramping and other samples pooled over animals per trajectory:
medians: {"CA1 ramping N-N'": 0.014, "CA1 ramping F-N'": -0.038, "CA1 other N-N'": 0.508, "CA1 other F-N'": 0.012, "PFC ramping N-N'": 0.905, "PFC ramping F-N'": 0.847, "PFC other N-N'": 0.895, "PFC other F-N'": 0.888}
Kruskal-Wallis H = 1766.5, p = 0
Dunn CA1 ramping N-N' vs CA1 ramping F-N': p = 1
Dunn CA1 ramping N-N' vs CA1 other N-N': p = 1.15e-07
Dunn CA1 ramping F-N' vs CA1 other F-N': p = 0.0954
Dunn PFC ramping N-N' vs PFC ramping F-N': p = 0.375
Dunn PFC ramping N-N' vs PFC other N-N': p = 0.304
Dunn PFC ramping F-N' vs PFC other F-N': p = 0.0164
Result: blocked (too few ramping cells). Ramping samples identified in N' are 4 to 7 per region and trajectory summed over the five rats, and only one animal-trajectory population per region has 3 or more, so the ramping groups are one population of 101 bins each. On it, CA1 ramping F-N' is below N-N' (p = 1.7e-22), PFC ramping does not change (p = 0.34), and ramping and other cells differ only for CA1 F-N' (p = 0.014). With the ramping cells pooled over animals, the CA1 ramping drop is not found (p = 1) and PFC ramping and other cells differ for F-N' (p = 0.016). The two ways of forming the populations disagree, and the ramping selection itself did not reproduce (figs4-ramping-fraction).
Single-cell geometry (Figures 3D, 3E, S5)¶
Cells found in all three sessions (included in N, F and N'; matchidx_acrossep.m). Selectivity index and path-equivalent coefficient per cell and session (nfn.selectivity, nfn.path_equivalence).
def sc_table(s):
rows = []
for d in A:
for reg in ("CA1", "PFC"):
cells = np.intersect1d(np.intersect1d(nfn.included(d, reg, 0), nfn.included(d, reg, 1)), nfn.included(d, reg, 2))
for ep in range(3):
m = MAPS[(d["animal"], s, ep)]
for c in cells:
si_in, si_out = nfn.selectivity(m, c)
si_in_g, si_out_g = nfn.selectivity(m, c, ">")
rows.append({"animal": d["animal"], "region": reg, "unit": c, "ep": ep, "si_in": si_in, "si_out": si_out,
"si_in_gt": si_in_g, "si_out_gt": si_out_g, "R": nfn.path_equivalence(m, c)})
return pd.DataFrame(rows)
SC = {s: sc_table(s) for s in SMOOTHS}
print("cells in all three sessions:", SC[8.0][SC[8.0].ep == 0].groupby("region").size().to_dict())
print("per animal:", SC[8.0][SC[8.0].ep == 0].groupby(["region", "animal"]).size().unstack(0).to_dict())
def pear(x, y):
ok = np.isfinite(x) & np.isfinite(y)
r, p = stats.pearsonr(x[ok], y[ok]); return r, p, int(ok.sum())
cells in all three sessions: {'CA1': 99, 'PFC': 171}
per animal: {'CA1': {'AM2': 27, 'JS17': 20, 'JS21': 8, 'JS34': 21, 'ZT2': 23}, 'PFC': {'AM2': 25, 'JS17': 56, 'JS21': 30, 'JS34': 12, 'ZT2': 48}}
fig3d-selectivity-index¶
Inbound and outbound trajectory selectivity indices are negatively correlated in PFC (r = -0.5), reflecting task-sequence selectivity, and positively correlated in CA1 (r = 0.56), reflecting spatial coverage.
Readings for the sessions: (a) N, F and N' pooled, (b) N', (c) F. Variants: the code's peak rule (> 3 Hz), 4 and 2 cm smoothing.
rows = []
for s in SMOOTHS:
S = SC[s]
for reading, eps in (("(a) pooled", [0, 1, 2]), ("(b) N'", [2]), ("(c) F", [1]), ("N (not a reading)", [0]), ("N and N' (not a reading)", [0, 2])):
for reg in ("PFC", "CA1"):
x = S[(S.region == reg) & S.ep.isin(eps)]
r, p, n = pear(x.si_in.values, x.si_out.values)
rg, pg, ng = pear(x.si_in_gt.values, x.si_out_gt.values)
rows.append({"smoothing": s, "reading": reading, "region": reg, "r": r, "p": p, "n": n, "r (> 3 Hz)": rg, "n (> 3 Hz)": ng})
T3D = pd.DataFrame(rows)
with pd.option_context("display.float_format", lambda v: f"{v:.3g}"):
print(T3D.to_string(index=False))
print("\nTolerance: PFC -0.55 to -0.45, CA1 0.504 to 0.616.")
for reg in ("PFC", "CA1"):
x = SC[8.0][(SC[8.0].region == reg) & (SC[8.0].ep == 2)]
print(f"{reg} N', per animal r:", {a: (round(pear(g.si_in.values, g.si_out.values)[0], 2), pear(g.si_in.values, g.si_out.values)[2])
for a, g in x.groupby("animal") if np.isfinite(g.si_in).sum() > 3})
smoothing reading region r p n r (> 3 Hz) n (> 3 Hz)
8 (a) pooled PFC -0.335 1.15e-10 351 -0.335 351
8 (a) pooled CA1 0.284 0.000215 165 0.284 165
8 (b) N' PFC -0.466 6.96e-08 121 -0.466 121
8 (b) N' CA1 0.28 0.0404 54 0.28 54
8 (c) F PFC -0.417 5.78e-06 110 -0.417 110
8 (c) F CA1 0.127 0.412 44 0.127 44
8 N (not a reading) PFC -0.143 0.12 120 -0.143 120
8 N (not a reading) CA1 0.361 0.00271 67 0.361 67
8 N and N' (not a reading) PFC -0.312 7.49e-07 241 -0.312 241
8 N and N' (not a reading) CA1 0.335 0.000172 121 0.335 121
4 (a) pooled PFC -0.319 2.57e-10 376 -0.319 376
4 (a) pooled CA1 0.308 1.71e-05 188 0.308 188
4 (b) N' PFC -0.414 1.16e-06 128 -0.414 128
4 (b) N' CA1 0.272 0.0326 62 0.272 62
4 (c) F PFC -0.454 3.56e-07 115 -0.454 115
4 (c) F CA1 0.237 0.0849 54 0.237 54
4 N (not a reading) PFC -0.165 0.0579 133 -0.165 133
4 N (not a reading) CA1 0.37 0.00139 72 0.37 72
4 N and N' (not a reading) PFC -0.283 3.49e-06 261 -0.283 261
4 N and N' (not a reading) CA1 0.334 7.81e-05 134 0.334 134
2 (a) pooled PFC -0.309 1.14e-10 418 -0.309 418
2 (a) pooled CA1 0.335 4.96e-07 215 0.335 215
2 (b) N' PFC -0.436 7.19e-08 140 -0.436 140
2 (b) N' CA1 0.249 0.0388 69 0.249 69
2 (c) F PFC -0.425 3.87e-07 132 -0.425 132
2 (c) F CA1 0.336 0.00617 65 0.336 65
2 N (not a reading) PFC -0.15 0.0701 146 -0.15 146
2 N (not a reading) CA1 0.4 0.000218 81 0.4 81
2 N and N' (not a reading) PFC -0.283 1.1e-06 286 -0.283 286
2 N and N' (not a reading) CA1 0.335 2.87e-05 150 0.335 150
Tolerance: PFC -0.55 to -0.45, CA1 0.504 to 0.616.
PFC N', per animal r: {'AM2': (np.float64(-0.52), 22), 'JS17': (np.float64(-0.41), 38), 'JS21': (np.float64(-0.42), 17), 'JS34': (np.float64(-0.87), 9), 'ZT2': (np.float64(-0.5), 35)}
CA1 N', per animal r: {'AM2': (np.float64(0.64), 15), 'JS17': (np.float64(-0.4), 7), 'JS21': (np.float64(0.89), 8), 'JS34': (np.float64(-0.04), 9), 'ZT2': (np.float64(0.06), 15)}
Result: blocked, methods underspecified (the smoothing width). PFC SI(IN) and SI(OUT) are negatively correlated under every reading and width, and in N' within tolerance at 8 cm (r = -0.466, p = 7.0e-08). CA1 is positively correlated but weaker than stated: at 8 cm r = 0.284 pooled, 0.280 in N' (p = 0.040) and 0.127 in F (p = 0.41), against 0.56. At 8 cm and 4 cm reading (c) is not significant for CA1 (p = 0.41 and 0.085), so the partial clause fails and the claim would not reproduce; at 2 cm every reading has CA1 r > 0 at p < 0.05 (0.249 to 0.336) and PFC r < 0, which would make it partially reproduced. The outcome changes with the width, and the claims file's rule for that case (fig2d-pv-similarity) makes the claim methods underspecified.
fig3e-selectivity-stable¶
Trajectory selectivity is stable across environments for PFC cells (r = 0.49 between familiar and novel) but not for CA1 cells (r = 0.12, not significant).
SI(IN) and SI(OUT) pooled in one correlation per region, F against the novel session: (a) N', (b) N.
rows = []
for s in SMOOTHS:
S = SC[s]
for reg in ("PFC", "CA1"):
x = S[S.region == reg].pivot_table(index=["animal", "unit"], columns="ep", values=["si_in", "si_out"])
for nov in (2, 0):
a = np.r_[x[("si_in", 1)], x[("si_out", 1)]]; b = np.r_[x[("si_in", nov)], x[("si_out", nov)]]
r, p, n = pear(a, b)
rows.append({"smoothing": s, "reading": f"F against {nfn.EPOCHS[nov]}", "region": reg, "r": r, "p": p, "n": n})
with pd.option_context("display.float_format", lambda v: f"{v:.3g}"):
print(pd.DataFrame(rows).to_string(index=False))
print("\nTolerance: PFC 0.441 to 0.539, CA1 0.108 to 0.132 with p > 0.05.")
smoothing reading region r p n
8 F against N' PFC 0.437 5.19e-11 206
8 F against N PFC 0.363 2.01e-07 194
8 F against N' CA1 0.318 0.0117 62
8 F against N CA1 0.384 0.000856 72
4 F against N' PFC 0.417 1.43e-10 218
4 F against N PFC 0.363 3.5e-08 218
4 F against N' CA1 0.231 0.0348 84
4 F against N CA1 0.301 0.00392 90
2 F against N' PFC 0.406 2.47e-11 250
2 F against N PFC 0.26 3.75e-05 246
2 F against N' CA1 0.128 0.196 104
2 F against N CA1 0.199 0.0324 116
Tolerance: PFC 0.441 to 0.539, CA1 0.108 to 0.132 with p > 0.05.
Result: blocked, methods underspecified (the smoothing width). PFC selectivity is correlated between F and the novel room at every width (F against N' r = 0.406 to 0.437, p <= 1.4e-10). At 8 cm and 4 cm CA1 is correlated too (F against N' r = 0.318 and 0.231, p = 0.012 and 0.035; F against N p <= 0.0039), so the claim would not reproduce. At 2 cm CA1 F against N' is not significant (r = 0.128, p = 0.20) and PFC is 0.406, below the tolerance, which would make it partially reproduced. The outcome changes with the width, and the claims file's rule for that case (fig2d-pv-similarity) makes the claim methods underspecified.
figs5d-path-equivalence¶
PFC cells show stronger path equivalence (similar firing at the same trajectory phase across trajectory types) than CA1 cells.
One-way ANOVA with Tukey's test over the six region by session groups.
for s in SMOOTHS:
S = SC[s]
names = [f"{r} {nfn.EPOCHS[e]}" for r in ("CA1", "PFC") for e in range(3)]
G = [S[(S.region == r) & (S.ep == e)].R.dropna().values for r in ("CA1", "PFC") for e in range(3)]
a, t = nfn.anova_tukey(G, names)
print(f"smoothing {s:g} cm: medians {dict(zip(names, [round(float(np.median(g)), 3) for g in G]))}; n {a['n']}; F = {a['F']:.2f}, p = {a['p']:.2g}")
for e in range(3):
r = t[(t.group1 == f"CA1 {nfn.EPOCHS[e]}") & (t.group2 == f"PFC {nfn.EPOCHS[e]}")].iloc[0]
print(f" {nfn.EPOCHS[e]}: PFC minus CA1 = {r['meandiff']:.3f}, Tukey p = {r['p-adj']:.4g}")
print("R range:", round(float(SC[8.0].R.min()), 3), "to", round(float(SC[8.0].R.max()), 3), "; missing:", int(SC[8.0].R.isna().sum()))
smoothing 8 cm: medians {'CA1 N': 0.056, 'CA1 F': 0.169, "CA1 N'": -0.009, 'PFC N': 0.304, 'PFC F': 0.405, "PFC N'": 0.381}; n [99, 99, 99, 171, 171, 171]; F = 22.97, p = 1.3e-21
N: PFC minus CA1 = 0.173, Tukey p = 0.0002
F: PFC minus CA1 = 0.225, Tukey p = 0
N': PFC minus CA1 = 0.303, Tukey p = 0
smoothing 4 cm: medians {'CA1 N': 0.02, 'CA1 F': 0.096, "CA1 N'": -0.003, 'PFC N': 0.231, 'PFC F': 0.324, "PFC N'": 0.338}; n [99, 99, 99, 171, 171, 171]; F = 25.90, p = 2.7e-24
N: PFC minus CA1 = 0.147, Tukey p = 0.0004
F: PFC minus CA1 = 0.217, Tukey p = 0
N': PFC minus CA1 = 0.277, Tukey p = 0
smoothing 2 cm: medians {'CA1 N': -0.003, 'CA1 F': 0.037, "CA1 N'": -0.005, 'PFC N': 0.159, 'PFC F': 0.242, "PFC N'": 0.273}; n [99, 99, 99, 171, 171, 171]; F = 26.55, p = 7e-25
N: PFC minus CA1 = 0.117, Tukey p = 0.002
F: PFC minus CA1 = 0.197, Tukey p = 0
N': PFC minus CA1 = 0.246, Tukey p = 0
R range: -0.496 to 0.996 ; missing: 0
Result: reproduced. PFC R exceeds CA1 R in every session (N: difference 0.173, Tukey p = 0.0002; F and N': p < 0.0001), with ANOVA F = 22.97, p = 1.3e-21, over 99 CA1 and 171 PFC cells. The 4-cm and 2-cm widths give the same (largest p 0.002).
figs5e-path-equivalence-stable¶
The path-equivalent coefficient of PFC cells is strongly correlated between environments (r = 0.64), more so than that of CA1 cells (r = 0.21).
Pearson r of R in F against R in the novel session ((a) N', (b) N), per region; Fisher's z for the difference.
for s in SMOOTHS:
S = SC[s]
for nov in (2, 0):
res = {}
for reg in ("PFC", "CA1"):
x = S[S.region == reg].pivot_table(index=["animal", "unit"], columns="ep", values="R")
res[reg] = pear(x[1].values, x[nov].values)
z, pz = nfn.fisher_z_compare(res["PFC"][0], res["PFC"][2], res["CA1"][0], res["CA1"][2])
print(f"smoothing {s:g} cm, F against {nfn.EPOCHS[nov]}: PFC r = {res['PFC'][0]:.3f} (p {res['PFC'][1]:.2g}, n {res['PFC'][2]}), "
f"CA1 r = {res['CA1'][0]:.3f} (p {res['CA1'][1]:.2g}, n {res['CA1'][2]}); Fisher z = {z:.2f}, p = {pz:.2g}")
print("Tolerance: PFC 0.576 to 0.704, CA1 0.189 to 0.231.")
smoothing 8 cm, F against N': PFC r = 0.572 (p 3e-16, n 171), CA1 r = 0.231 (p 0.021, n 99); Fisher z = 3.24, p = 0.0012 smoothing 8 cm, F against N: PFC r = 0.552 (p 4.9e-15, n 171), CA1 r = 0.095 (p 0.35, n 99); Fisher z = 4.11, p = 3.9e-05 smoothing 4 cm, F against N': PFC r = 0.613 (p 5e-19, n 171), CA1 r = 0.201 (p 0.046, n 99); Fisher z = 3.98, p = 6.8e-05 smoothing 4 cm, F against N: PFC r = 0.521 (p 2.7e-13, n 171), CA1 r = 0.127 (p 0.21, n 99); Fisher z = 3.52, p = 0.00044 smoothing 2 cm, F against N': PFC r = 0.647 (p 1.1e-21, n 171), CA1 r = 0.139 (p 0.17, n 99); Fisher z = 4.93, p = 8.2e-07 smoothing 2 cm, F against N: PFC r = 0.516 (p 5.1e-13, n 171), CA1 r = 0.153 (p 0.13, n 99); Fisher z = 3.25, p = 0.0011 Tolerance: PFC 0.576 to 0.704, CA1 0.189 to 0.231.
Result: blocked, methods underspecified (the smoothing width). F against N' at 8 cm, PFC r = 0.572 and CA1 r = 0.231, and Fisher's z = 3.24, p = 0.0012, puts PFC above CA1; CA1 is within tolerance and PFC 0.004 below it, which would make it partially reproduced. With 4-cm smoothing both values are within tolerance (0.613 and 0.201, z p = 6.8e-05), which would make it reproduced; 2 cm gives 0.647 and 0.139 (partially). F against N gives PFC 0.516 to 0.552 and CA1 0.095 to 0.153. The outcome changes with the width, and the claims file's rule for that case (fig2d-pv-similarity) makes the claim methods underspecified.
figs5f-path-equivalence-vs-selectivity¶
Path equivalence and task-sequence selectivity are uncorrelated across cells.
Pearson r per region of R against SI(IN) and against SI(OUT), cell samples of the three sessions pooled; signed and absolute SI.
for s in SMOOTHS:
S = SC[s]
for reg in ("CA1", "PFC"):
for col in ("si_in", "si_out"):
x = S[S.region == reg]
r, p, n = pear(x.R.values, x[col].values); ra, pa, na = pear(x.R.values, np.abs(x[col].values))
print(f"smoothing {s:g} cm, {reg} R against {col}: signed r = {r:.3f}, p = {p:.3g}; absolute r = {ra:.3f}, p = {pa:.3g} (n {n})")
smoothing 8 cm, CA1 R against si_in: signed r = -0.009, p = 0.904; absolute r = -0.174, p = 0.0254 (n 165) smoothing 8 cm, CA1 R against si_out: signed r = 0.004, p = 0.964; absolute r = -0.134, p = 0.0859 (n 165) smoothing 8 cm, PFC R against si_in: signed r = 0.011, p = 0.844; absolute r = -0.241, p = 5.18e-06 (n 351) smoothing 8 cm, PFC R against si_out: signed r = 0.060, p = 0.262; absolute r = -0.208, p = 8.71e-05 (n 351) smoothing 4 cm, CA1 R against si_in: signed r = -0.026, p = 0.726; absolute r = -0.153, p = 0.0365 (n 188) smoothing 4 cm, CA1 R against si_out: signed r = 0.013, p = 0.858; absolute r = -0.161, p = 0.0272 (n 188) smoothing 4 cm, PFC R against si_in: signed r = 0.020, p = 0.695; absolute r = -0.261, p = 3e-07 (n 376) smoothing 4 cm, PFC R against si_out: signed r = 0.065, p = 0.211; absolute r = -0.228, p = 7.64e-06 (n 376) smoothing 2 cm, CA1 R against si_in: signed r = -0.025, p = 0.71; absolute r = -0.118, p = 0.0831 (n 215) smoothing 2 cm, CA1 R against si_out: signed r = 0.043, p = 0.533; absolute r = -0.163, p = 0.017 (n 215) smoothing 2 cm, PFC R against si_in: signed r = 0.012, p = 0.806; absolute r = -0.252, p = 1.81e-07 (n 418) smoothing 2 cm, PFC R against si_out: signed r = 0.097, p = 0.0473; absolute r = -0.241, p = 5.92e-07 (n 418)
Result: blocked, methods underspecified. With the signed SI no correlation reaches p < 0.05 at 8 cm (smallest p 0.26). With the absolute SI three of four do (CA1 SI(IN) r = -0.174, p = 0.025; PFC SI(IN) r = -0.241 and SI(OUT) r = -0.208, p < 0.0001). The criterion makes the claim methods underspecified when only one reading finds a correlation. At 4 cm the same holds (signed: none, smallest p 0.21; absolute: four of four, p <= 0.037). At 2 cm both readings find one: the signed SI gives PFC SI(OUT) r = 0.097, p = 0.047, and the absolute SI three of four (CA1 SI(OUT) p = 0.017, PFC p <= 5.9e-07), so at 2 cm the claim would not reproduce. The outcome changes with the width, so the claim stays methods underspecified under the claims file's width rule as well.
Manifolds (Figure 3F to 3J, S6)¶
manifold-cell-counts¶
The pseudo-populations pooled across animals for the manifold analyses hold 98 CA1 and 171 PFC cells.
Readings: (a) the Methods', units included in N, F and N' with a linearized peak of at least 3 Hz in at least one session; (b) the code's, every unit included in all three sessions.
POP = {}
for reg in ("CA1", "PFC"):
b = {d["animal"]: np.intersect1d(np.intersect1d(nfn.included(d, reg, 0), nfn.included(d, reg, 1)), nfn.included(d, reg, 2)) for d in A}
a = {}
for d in A:
cs = b[d["animal"]]
pk = np.max([nfn.peak_rate({t: (MAPS[(d["animal"], 8.0, e)][t][0], MAPS[(d["animal"], 8.0, e)][t][1][cs]) for t in nfn.TRACKS}) for e in range(3)], 0)
a[d["animal"]] = cs[pk >= 3]
POP[(reg, "a")], POP[(reg, "b")] = a, b
print(f"{reg}: reading (a) {sum(len(v) for v in a.values())} cells {({k: len(v) for k, v in a.items()})}; "
f"reading (b) {sum(len(v) for v in b.values())} cells {({k: len(v) for k, v in b.items()})}")
print("Expected: CA1 98 (88.2 to 107.8), PFC 171 (153.9 to 188.1).")
CA1: reading (a) 95 cells {'AM2': 26, 'JS17': 20, 'JS21': 8, 'ZT2': 23, 'JS34': 18}; reading (b) 99 cells {'AM2': 27, 'JS17': 20, 'JS21': 8, 'ZT2': 23, 'JS34': 21}
PFC: reading (a) 151 cells {'AM2': 25, 'JS17': 51, 'JS21': 23, 'ZT2': 42, 'JS34': 10}; reading (b) 171 cells {'AM2': 25, 'JS17': 56, 'JS21': 30, 'ZT2': 48, 'JS34': 12}
Expected: CA1 98 (88.2 to 107.8), PFC 171 (153.9 to 188.1).
Result: reproduced. The code's reading (units included in N, F and N') gives 99 CA1 and 171 PFC cells, against 98 and 171. The Methods' reading (also tuned) gives 95 and 151; PFC is below the tolerance. The code's reading is used for the manifolds.
The embeddings¶
Single-trial maps of the pseudo-population in N' (and F), as cal_UMAP.m: for each trajectory type the
smallest number of passes over the animals, from the first; maps at 1-cm bins (Methods; smoothed at 8 cm
as the session maps, as the Methods give no width for them) on 201 bins of normalized length (cal_UMAP.m; about 1 cm on a 2-m trajectory), empty bins
inside a pass filled (an assumption, see nfn.single_trial_matrix); states with fewer than 5 active
cells dropped. UMAP (umap-learn) with 3 dimensions, cosine metric, 50 neighbours, min_dist 0.6, seed
0; F is embedded with the N' fit (transform). Seeds 1 and 2 are variants.
UMAP's numba code is compiled for the host CPU, and the optimization amplifies rounding differences, so before this round the embeddings, and every number derived from them, differed between machines (the verifier's copy of round 2 differed from the author's). The header cell now compiles for a generic x86-64 target and fixes the BLAS kernel and numpy's dispatch. On one machine the embedding repeats exactly: a second fit with seed 0 below gives the same coordinates (and, a check made outside the notebook, two full executions gave identical outputs in every code cell). Seeds 1 and 2 show how far the results move with the embedding.
import umap
POPSEL = {reg: POP[(reg, "b")] for reg in ("CA1", "PFC")}
ST = {}
for reg in ("CA1", "PFC"):
for ep in (2, 1):
tbt = {d["animal"]: nfn.trials_of_track(d, ep) for d in A}
tn = {t: min(len(tbt[a][t]) for a in tbt) for t in nfn.TRACKS}
Xs = []
for d in A:
X, lab = nfn.single_trial_matrix(d, ep, POPSEL[reg][d["animal"]], {t: tbt[d["animal"]][t][:tn[t]] for t in nfn.TRACKS})
Xs.append(X)
X = np.hstack(Xs); v = (X > 0).sum(1) >= 5
ST[(reg, ep)] = (X[v], lab[v], tn, int((~v).sum()))
print(f"{reg} {nfn.EPOCHS[ep]}: passes per trajectory {({nfn.TRACKS[k]: n for k, n in tn.items()})}; states {v.sum()} kept, {(~v).sum()} dropped (fewer than 5 active cells); cells {X.shape[1]}")
EMB = {}
for reg in ("CA1", "PFC"):
for seed in (0, 1, 2):
t0 = time.time()
um = umap.UMAP(n_components=3, metric="cosine", n_neighbors=50, min_dist=0.6, random_state=seed)
E = um.fit_transform(ST[(reg, 2)][0]); EF = um.transform(ST[(reg, 1)][0])
EMB[(reg, seed)] = (E, EF, um)
print(f"{reg} seed {seed}: fit and transform in {time.time() - t0:.0f} s")
for reg in ("CA1", "PFC"): # a second fit with seed 0 in the same run, to check that the embedding repeats
um2 = umap.UMAP(n_components=3, metric="cosine", n_neighbors=50, min_dist=0.6, random_state=0)
E2 = um2.fit_transform(ST[(reg, 2)][0]); EF2 = um2.transform(ST[(reg, 1)][0])
print(f"{reg} seed 0 refit: largest difference from the first fit {np.abs(E2 - EMB[(reg, 0)][0]).max():.3g} (N'), "
f"{np.abs(EF2 - EMB[(reg, 0)][1]).max():.3g} (F)")
fig = plt.figure(figsize=(14, 7))
for i, reg in enumerate(("CA1", "PFC")):
E, EF, _ = EMB[(reg, 0)]; lab = ST[(reg, 2)][1]
for j, (a_, b_) in enumerate(((0, 1), (0, 2), (1, 2))):
ax = fig.add_subplot(2, 4, i * 4 + j + 1)
for t, c in zip(nfn.TRACKS, ("tab:red", "tab:blue", "tab:green", "tab:purple")):
m = lab[:, 0] == t; ax.scatter(E[m, a_], E[m, b_], s=1, c=c, alpha=0.3, label=nfn.TRACKS[t])
ax.set_title(f"{reg} N', dims {a_ + 1}-{b_ + 1}", fontsize=9)
ax = fig.add_subplot(2, 4, i * 4 + 4)
sc_ = ax.scatter(E[:, 0], E[:, 1], s=1, c=lab[:, 2], cmap="viridis"); ax.set_title(f"{reg} N', colour = position", fontsize=9)
fig.axes[0].legend(markerscale=8, fontsize=7); plt.tight_layout(); plt.show()
print("All states are plotted (axes span the data).")
CA1 N': passes per trajectory {'OUT_L': 5, 'IN_L': 6, 'OUT_R': 8, 'IN_R': 6}; states 5020 kept, 5 dropped (fewer than 5 active cells); cells 99
CA1 F: passes per trajectory {'OUT_L': 3, 'IN_L': 4, 'OUT_R': 6, 'IN_R': 4}; states 3417 kept, 0 dropped (fewer than 5 active cells); cells 99
PFC N': passes per trajectory {'OUT_L': 5, 'IN_L': 6, 'OUT_R': 8, 'IN_R': 6}; states 5023 kept, 2 dropped (fewer than 5 active cells); cells 171
PFC F: passes per trajectory {'OUT_L': 3, 'IN_L': 4, 'OUT_R': 6, 'IN_R': 4}; states 3417 kept, 0 dropped (fewer than 5 active cells); cells 171
CA1 seed 0: fit and transform in 91 s
CA1 seed 1: fit and transform in 9 s
CA1 seed 2: fit and transform in 9 s
PFC seed 0: fit and transform in 9 s
PFC seed 1: fit and transform in 9 s
PFC seed 2: fit and transform in 9 s
CA1 seed 0 refit: largest difference from the first fit 0 (N'), 0 (F)
PFC seed 0 refit: largest difference from the first fit 0 (N'), 0 (F)
All states are plotted (axes span the data).
# Distance measures of Figure 3H to 3J and S6, used by the claims below and by the variants.
from scipy.spatial.distance import cdist
def seq_dist(E, lab, tn, align, nb=201):
idx = {(int(t), int(k), int(b)): i for i, (t, k, b) in enumerate(lab)}
out = {}
for name, pairs in (("InSeq", [(1, 4), (3, 2)]), ("OutSeq", [(1, 2), (3, 4)])):
vals = []
for a_, b_ in pairs:
for t1 in range(tn[a_]):
for t2 in range(tn[b_]):
ds = [np.linalg.norm(E[idx[(a_, t1, i)]] - E[idx[(b_, t2, j)]])
for i in range(nb) for j in [i if align == "location" else nb - 1 - i]
if (a_, t1, i) in idx and (b_, t2, j) in idx]
if ds:
vals.append(np.mean(ds))
out[name] = np.array(vals)
return out
def h_test(E, reg, align):
X, lab, tn, _ = ST[(reg, 2)]
sd = seq_dist(E, lab, tn, align)
return np.median(sd["InSeq"]), np.median(sd["OutSeq"]), stats.ranksums(sd["InSeq"], sd["OutSeq"]).pvalue, len(sd["InSeq"]), len(sd["OutSeq"])
CODE_ALIGN = {"CA1": "location", "PFC": "phase"}
def phase_loc(E, lab, mode, nb=201):
out = {"phase": [], "location": []}
sel = {(t, b): np.where((lab[:, 0] == t) & (lab[:, 2] == b))[0] for t in nfn.TRACKS for b in range(nb)}
for i in range(nb):
for a_, b_ in ((1, 2), (3, 4)):
p, q = sel[(a_, i)], sel[(b_, i)]
if len(p) and len(q): out["location"].append(cdist(E[p], E[q]).ravel())
for a_, b_ in (((1, 4), (3, 2)) if mode == "methods" else ((1, 2), (3, 4))):
p, q = sel[(a_, i)], sel[(b_, nb - 1 - i)]
if len(p) and len(q): out["phase"].append(cdist(E[p], E[q]).ravel())
return {k: np.concatenate(v) for k, v in out.items()}
def fn_dist(E, lab, EF, labF):
kE, kF = lab[:, 0] * 1000 + lab[:, 2], labF[:, 0] * 1000 + labF[:, 2]
out = np.full(len(E), np.nan)
for k in np.unique(kE):
p, q = np.where(kE == k)[0], np.where(kF == k)[0]
if len(q): out[p] = cdist(E[p], EF[q]).mean(1)
return out
def shuffled_activity(X, lab, rng):
Xs = X.copy()
for t, k in {(int(a_), int(b_)) for a_, b_, _ in lab}:
rows = np.where((lab[:, 0] == t) & (lab[:, 1] == k))[0]
Xs[rows] = X[rows][:, rng.permutation(X.shape[1])]
return Xs
Manifold variants: gap filling and smoothing of the single-trial maps¶
The manifold results are repeated without the linear filling of empty 1-cm bins inside a pass (an
assumption, nfn.single_trial_matrix). The Methods give no smoothing for the single-trial
maps; the session maps' 1D width is garbled in the text (fig2d-pv-similarity), so 4 and 2 cm are
run too. Each variant rebuilds the states, fits UMAP (seed 0) and repeats the tests of 3H, 3I and 3J,
the same distances in the neural space (Figure S6), and the 16 CCGP decoders (100 shuffles each, for
time, so the smallest p is 0.01; p_ge counts ties, and a decoder above every shuffle has p_ge = 0).
# Variants of the single-trial maps behind the manifolds: empty 1-cm bins inside a pass left empty
# (set to 0, as cal_UMAP.m does with NaN) instead of filled, and other smoothing widths. Each variant
# rebuilds the states, fits UMAP (seed 0) and repeats the tests of 3H, 3I, 3J and the 16 CCGP
# decoders (100 shuffles each, for time).
def build_states(reg, ep, fill, s):
tbt = {d["animal"]: nfn.trials_of_track(d, ep) for d in A}
tn = {t: min(len(tbt[a][t]) for a in tbt) for t in nfn.TRACKS}
Xs = []
for d in A:
X, lab = nfn.single_trial_matrix(d, ep, POPSEL[reg][d["animal"]], {t: tbt[d["animal"]][t][:tn[t]] for t in nfn.TRACKS},
smooth_cm=s, fill_gaps=fill)
Xs.append(X)
X = np.hstack(Xs); v = (X > 0).sum(1) >= 5
return X[v], lab[v], tn, int((~v).sum())
def h_test2(E, lab, tn, align):
sd = seq_dist(E, lab, tn, align)
return np.median(sd["InSeq"]), np.median(sd["OutSeq"]), stats.ranksums(sd["InSeq"], sd["OutSeq"]).pvalue
MANVAR = []
for vname, fill, s in (("main: gaps filled, 8 cm", True, 8.0), ("gaps not filled, 8 cm", False, 8.0),
("gaps filled, 4 cm", True, 4.0), ("gaps filled, 2 cm", True, 2.0), ("gaps not filled, 2 cm", False, 2.0)):
for reg in ("CA1", "PFC"):
if fill and s == 8.0:
X, lab, tn, ndrop = ST[(reg, 2)]; XF, labF, _, _ = ST[(reg, 1)]
E, EF, um = EMB[(reg, 0)]
else:
X, lab, tn, ndrop = build_states(reg, 2, fill, s); XF, labF, _, _ = build_states(reg, 1, fill, s)
um = umap.UMAP(n_components=3, metric="cosine", n_neighbors=50, min_dist=0.6, random_state=0)
E = um.fit_transform(X); EF = um.transform(XF)
i_, o_, ph = h_test2(E, lab, tn, CODE_ALIGN[reg])
pl = phase_loc(E, lab, "methods")
Xs = shuffled_activity(X, lab, np.random.default_rng([501, 0]))
ES = um.transform(Xs)
real, shuf = fn_dist(E, lab, EF, labF), fn_dist(E, lab, ES, lab)
rr, ss = real[np.isfinite(real)], shuf[np.isfinite(shuf)]
hi_, ho_, hp_ = h_test2(X, lab, tn, CODE_ALIGN[reg]) # Figure S6: the same tests in the neural space
pl_ = phase_loc(X, lab, "methods")
nr, ns = fn_dist(X, lab, XF, labF), fn_dist(X, lab, Xs, lab); nr, ns = nr[np.isfinite(nr)], ns[np.isfinite(ns)]
cc = nfn.run_ccgp(nfn.clusters(E, lab), nfn.clusters(EF, labF), n_perm=100, seed=901)
sig = (cc.p_ge < 0.05)
MANVAR.append({"variant": vname, "region": reg, "states dropped": ndrop,
"3H InSeq": i_, "3H OutSeq": o_, "3H p": ph,
"3I phase": np.median(pl["phase"]), "3I location": np.median(pl["location"]),
"3I p": stats.ranksums(pl["phase"], pl["location"]).pvalue,
"3J F-N' vs shuffle %": 100 * (np.median(rr) / np.median(ss) - 1), "3J p": stats.ranksums(rr, ss).pvalue,
"CCGP phase mean": cc[cc.variable == "phase"].ccgp.mean(), "CCGP sequence mean": cc[cc.variable == "sequence"].ccgp.mean(),
"CCGP environment mean": cc[cc.variable == "environment"].ccgp.mean(),
"decoders p_ge < 0.05": int(sig.sum()), "decoders above every shuffle": int((cc.p_ge == 0).sum()),
"S6 3H InSeq - OutSeq": hi_ - ho_, "S6 3H p": hp_,
"S6 3I phase - location": np.median(pl_["phase"]) - np.median(pl_["location"]),
"S6 3I p": stats.ranksums(pl_["phase"], pl_["location"]).pvalue,
"S6 3J F-N' vs shuffle %": 100 * (np.median(nr) / np.median(ns) - 1), "S6 3J p": stats.ranksums(nr, ns).pvalue})
MV = pd.DataFrame(MANVAR)
with pd.option_context("display.float_format", lambda v: f"{v:.3g}"):
print(MV.iloc[:, :15].to_string(index=False))
print(MV[["variant", "region"] + list(MV.columns[15:])].to_string(index=False))
variant region states dropped 3H InSeq 3H OutSeq 3H p 3I phase 3I location 3I p 3J F-N' vs shuffle % 3J p CCGP phase mean CCGP sequence mean CCGP environment mean decoders p_ge < 0.05
main: gaps filled, 8 cm CA1 5 17.3 13.7 0.0005 17.5 14.3 0 -4.37 1.38e-30 0.5 0.371 0.568 8
main: gaps filled, 8 cm PFC 2 14 17.5 8.95e-06 14.1 17.7 0 -22.4 0 0.549 0.745 0.681 12
gaps not filled, 8 cm CA1 259 6.13 3.78 4.18e-27 6.05 3.65 0 -1.87 4.63e-19 0.669 0.161 0.497 6
gaps not filled, 8 cm PFC 218 11 11.1 0.0216 10.6 10.7 9.16e-43 -7.76 5.73e-29 0.507 0.538 0.526 4
gaps filled, 4 cm CA1 9 17.9 10.2 3.15e-14 16.6 12.3 0 -10.8 7.21e-75 0.621 0.675 0.608 11
gaps filled, 4 cm PFC 2 9.18 11.4 3.24e-11 9.49 13 0 -40.8 0 0.729 0.813 0.804 14
gaps filled, 2 cm CA1 15 6.7 4.11 4.18e-27 8.17 4.14 0 -5.27 6.16e-28 0.536 0.142 0.519 4
gaps filled, 2 cm PFC 2 2.47 3.7 8.68e-24 2.42 6.85 0 -60.2 0 0.99 0.786 0.818 15
gaps not filled, 2 cm CA1 630 4.51 3.62 4.18e-27 4.32 3.32 0 -5.33 2.98e-24 0.595 0.29 0.488 5
gaps not filled, 2 cm PFC 221 16.4 16.7 0.339 6.78 7.82 2.58e-47 0.821 0.0287 0.582 0.527 0.512 5
variant region decoders above every shuffle S6 3H InSeq - OutSeq S6 3H p S6 3I phase - location S6 3I p S6 3J F-N' vs shuffle % S6 3J p
main: gaps filled, 8 cm CA1 6 7.11 2.49e-12 3.34 9.33e-129 -9.21 7.5e-78
main: gaps filled, 8 cm PFC 10 -6.34 4.79e-12 -20.6 0 -56.6 0
gaps not filled, 8 cm CA1 5 3.73 1.16e-06 2.09 2.38e-47 -8.26 1.55e-31
gaps not filled, 8 cm PFC 4 -1.61 3.15e-06 -10.7 1.95e-154 -29.1 0
gaps filled, 4 cm CA1 9 7.52 3.95e-10 1.45 3.01e-23 -8.76 4.88e-68
gaps filled, 4 cm PFC 14 -5.64 4.53e-11 -19.6 0 -50.6 0
gaps filled, 2 cm CA1 4 6.74 5.1e-08 0.987 7.8e-06 -8.93 5.19e-63
gaps filled, 2 cm PFC 15 -4.95 4.43e-11 -18.1 0 -43.4 0
gaps not filled, 2 cm CA1 5 4.04 1.24e-05 0.856 4.02e-05 -7.3 2.74e-21
gaps not filled, 2 cm PFC 4 -1.38 5.17e-05 -9.24 8.89e-130 -23.2 0
fig3h-inseq-outseq-distance¶
In the PFC manifold, trajectories of the same task sequence (InSeq) are closer than trajectories out of sequence (OutSeq); in the CA1 manifold, OutSeq pairs with the same spatial coverage are closer than InSeq pairs.
For every pair of passes of the two types, the mean Euclidean distance over matched bins (cal_UMAP_INSeq_distance.m): CA1 matched by location, PFC by phase (the code); both regions by phase and both by location are variants. Rank-sum test per region.
for seed in (0, 1, 2):
for reg in ("CA1", "PFC"):
E = EMB[(reg, seed)][0]
for align in ("location", "phase"):
i_, o_, p, ni, no = h_test(E, reg, align)
tag = " (main)" if (align == CODE_ALIGN[reg] and seed == 0) else ""
print(f"seed {seed}, {reg}, aligned by {align}{tag}: InSeq median {i_:.3f} (n {ni}), OutSeq median {o_:.3f} (n {no}), rank-sum p = {p:.2g}")
seed 0, CA1, aligned by location (main): InSeq median 17.270 (n 78), OutSeq median 13.665 (n 78), rank-sum p = 0.0005
seed 0, CA1, aligned by phase: InSeq median 17.545 (n 78), OutSeq median 13.282 (n 78), rank-sum p = 0.00056 seed 0, PFC, aligned by location: InSeq median 15.227 (n 78), OutSeq median 17.343 (n 78), rank-sum p = 0.00084
seed 0, PFC, aligned by phase (main): InSeq median 14.044 (n 78), OutSeq median 17.510 (n 78), rank-sum p = 8.9e-06 seed 1, CA1, aligned by location: InSeq median 18.073 (n 78), OutSeq median 14.503 (n 78), rank-sum p = 1.4e-05 seed 1, CA1, aligned by phase: InSeq median 19.050 (n 78), OutSeq median 14.170 (n 78), rank-sum p = 2.1e-05
seed 1, PFC, aligned by location: InSeq median 15.164 (n 78), OutSeq median 18.005 (n 78), rank-sum p = 0.0098 seed 1, PFC, aligned by phase: InSeq median 13.931 (n 78), OutSeq median 17.283 (n 78), rank-sum p = 0.00033 seed 2, CA1, aligned by location: InSeq median 16.481 (n 78), OutSeq median 14.090 (n 78), rank-sum p = 0.0026
seed 2, CA1, aligned by phase: InSeq median 17.649 (n 78), OutSeq median 14.664 (n 78), rank-sum p = 0.00093 seed 2, PFC, aligned by location: InSeq median 14.674 (n 78), OutSeq median 17.004 (n 78), rank-sum p = 0.069 seed 2, PFC, aligned by phase: InSeq median 14.421 (n 78), OutSeq median 15.983 (n 78), rank-sum p = 0.042
Result: reproduced, resting on the filling of empty bins in the single-trial maps. In CA1 (aligned by location), OutSeq pairs are closer than InSeq pairs (13.67 against 17.27, rank-sum p = 0.0005); in PFC (aligned by phase), InSeq pairs are closer (14.04 against 17.51, p = 8.9e-06). With the main alignment both regions keep the direction under seeds 1 and 2 (CA1 p = 1.4e-05 and 0.0026, PFC p = 0.00033 and 0.042). Of the other alignments, PFC by location is weaker under seed 2 (p = 0.069); the rest give p < 0.01. The 4-cm and 2-cm widths agree (p <= 3.2e-11), and so do the maps without gap filling at 8 cm (CA1 p = 4.2e-27, PFC 11.0 against 11.1, p = 0.022); without gap filling at 2 cm PFC InSeq and OutSeq do not differ (16.4 and 16.7, p = 0.34), so the outcome rests on the filling. The result also moves with the embedding: in the verifier's round-2 copy, embedded on another machine before the CPU target was fixed, CA1 with the main alignment gave p = 0.096 under seed 1. Seed 0, the standard choice the criterion names, passed in every execution.
fig3i-phase-location-distance¶
PFC neural states at the same trajectory phase are closer than those at the same spatial location, while CA1 shows the opposite.
Main reading (Methods): same-phase distances from InSeq pairs aligned by phase, same-location distances from same-side pairs aligned by location; all states at matched bins pooled. Variant (cal_UMAP_TrajPhase_vs_Spatial_distance.m): both from the same-side pairs.
for seed in (0, 1, 2):
for reg in ("CA1", "PFC"):
lab = ST[(reg, 2)][1]
for mode in ("methods", "code"):
pl = phase_loc(EMB[(reg, seed)][0], lab, mode)
print(f"seed {seed}, {reg}, {mode}: phase median {np.median(pl['phase']):.3f}, location median {np.median(pl['location']):.3f}, "
f"rank-sum p = {stats.ranksums(pl['phase'], pl['location']).pvalue:.2g} (n {len(pl['phase'])}, {len(pl['location'])})")
seed 0, CA1, methods: phase median 17.451, location median 14.326, rank-sum p = 0 (n 15644, 15647) seed 0, CA1, code: phase median 13.994, location median 14.326, rank-sum p = 2.4e-09 (n 15647, 15647) seed 0, PFC, methods: phase median 14.077, location median 17.721, rank-sum p = 0 (n 15668, 15662) seed 0, PFC, code: phase median 17.623, location median 17.721, rank-sum p = 0.68 (n 15662, 15662) seed 1, CA1, methods: phase median 18.624, location median 15.284, rank-sum p = 0 (n 15644, 15647) seed 1, CA1, code: phase median 14.623, location median 15.284, rank-sum p = 0.03 (n 15647, 15647) seed 1, PFC, methods: phase median 13.933, location median 17.404, rank-sum p = 0 (n 15668, 15662) seed 1, PFC, code: phase median 16.792, location median 17.404, rank-sum p = 1.9e-08 (n 15662, 15662) seed 2, CA1, methods: phase median 17.753, location median 14.608, rank-sum p = 0 (n 15644, 15647) seed 2, CA1, code: phase median 15.386, location median 14.608, rank-sum p = 5.7e-40 (n 15647, 15647) seed 2, PFC, methods: phase median 13.856, location median 16.870, rank-sum p = 1.1e-308 (n 15668, 15662)
seed 2, PFC, code: phase median 16.082, location median 16.870, rank-sum p = 5.9e-17 (n 15662, 15662)
Result: reproduced. In PFC, same-phase states are closer than same-location states (14.08 against 17.72, p near 0); in CA1 the opposite holds (location 14.33 against phase 17.45, p near 0). Seeds 1 and 2, the 4-cm and 2-cm widths, and the maps without gap filling all agree (p <= 9.2e-43). The code's variant (both distances from same-side pairs) fails under seed 0 (CA1 reversed, phase 13.99 against location 14.33, p = 2.4e-09; PFC p = 0.68) and for CA1 under seed 1 (reversed, p = 0.03), and holds under seed 2.
fig3j-familiar-novel-manifold¶
The PFC manifold is preserved across environments (F states lie closer to N' states than cell-ID shuffles do), while the distortion of the CA1 manifold is similar to that of random shuffles.
For each N' state, the mean Euclidean distance to the F states of the same trajectory type and bin (cal_UMAP_FNdistance.m). Chance: the N' activity with cell identities permuted, a new permutation for each pass, embedded with the same fit (cal_UMAP_FNdistance_shuffle.m). Rank-sum test per region.
FNJ = {}
for seed in (0, 1, 2):
for reg in ("CA1", "PFC"):
X, lab, _, _ = ST[(reg, 2)]; XF, labF, _, _ = ST[(reg, 1)]
E, EF, um = EMB[(reg, seed)]
Xs = shuffled_activity(X, lab, np.random.default_rng([501, seed]))
ES = um.transform(Xs)
real, shuf = fn_dist(E, lab, EF, labF), fn_dist(E, lab, ES, lab)
FNJ[(reg, seed)] = (X, lab, XF, labF, Xs)
rr, ss = real[np.isfinite(real)], shuf[np.isfinite(shuf)]
print(f"seed {seed}, {reg}: F-N' median {np.median(rr):.3f}, shuffle median {np.median(ss):.3f} "
f"(F-N' {100 * (np.median(rr) / np.median(ss) - 1):+.1f}%), rank-sum p = {stats.ranksums(rr, ss).pvalue:.2g} (n {len(rr)}, {len(ss)})")
seed 0, CA1: F-N' median 13.437, shuffle median 14.050 (F-N' -4.4%), rank-sum p = 1.4e-30 (n 5020, 5020)
seed 0, PFC: F-N' median 11.858, shuffle median 15.281 (F-N' -22.4%), rank-sum p = 0 (n 5023, 5023)
seed 1, CA1: F-N' median 13.073, shuffle median 14.033 (F-N' -6.8%), rank-sum p = 1.4e-62 (n 5020, 5020)
seed 1, PFC: F-N' median 11.055, shuffle median 16.131 (F-N' -31.5%), rank-sum p = 0 (n 5023, 5023)
seed 2, CA1: F-N' median 13.116, shuffle median 13.672 (F-N' -4.1%), rank-sum p = 1.5e-38 (n 5020, 5020)
seed 2, PFC: F-N' median 11.883, shuffle median 13.457 (F-N' -11.7%), rank-sum p = 1.5e-114 (n 5023, 5023)
Result: partially reproduced, resting on the filling of empty bins. PFC F states lie closer to N' states than the shuffle does (11.86 against 15.28, 22.4% closer, p near 0). CA1 F states are also closer than the shuffle, by less (13.44 against 14.05, 4.4% closer, p = 1.4e-30), so the claim that CA1 does not differ from the shuffle fails by the paper's test. Seeds 1 and 2 give the same (CA1 6.8% and 4.1% closer, PFC 31.5% and 11.7%). CA1 is closer than the shuffle in every width and gap-filling variant (1.9% to 10.8%, p <= 4.6e-19), and PFC is in all but one (7.8% to 60.2% closer): without gap filling at 2 cm PFC F-N' is not closer than the shuffle (0.8% farther, p = 0.029), so there neither part holds and the outcome rests on the filling. That variant moves with the embedding: in the verifier's round-2 copy, embedded on another machine before the CPU target was fixed, it gave PFC 4.9% closer (p = 6.2e-12).
fig3fg-manifold-shapes¶
UMAP embeddings show that PFC population activity lies on a ring-shaped manifold, whereas the CA1 manifold has a double V shape resembling the W-maze.
The embeddings are shown above. The paper gives no measure of the shapes; by the criterion the outcome is that of fig3h-inseq-outseq-distance and fig3i-phase-location-distance.
for reg in ("CA1", "PFC"):
E = EMB[(reg, 0)][0]; lab = ST[(reg, 2)][1]
cent = np.array([[E[(lab[:, 0] == t) & (lab[:, 2] == b)].mean(0) for b in range(0, 201, 25)] for t in nfn.TRACKS])
print(f"{reg}: mean embedded state per trajectory at normalized positions 0, 0.125, ..., 1 (dimension 1):")
print(pd.DataFrame(cent[:, :, 0], index=[nfn.TRACKS[t] for t in nfn.TRACKS], columns=[f"{b / 200:.3f}" for b in range(0, 201, 25)]).round(2))
CA1: mean embedded state per trajectory at normalized positions 0, 0.125, ..., 1 (dimension 1):
0.000 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000
OUT_L 1.78 -0.04 0.06 0.98 2.58 4.080000 4.860000 4.82 4.64
IN_L 7.49 9.52 11.75 13.55 15.56 17.190001 19.049999 21.01 14.33
OUT_R -1.30 -2.23 -2.51 -2.27 -1.25 0.830000 3.350000 5.98 7.77
IN_R 7.17 7.20 6.42 4.42 2.29 1.890000 2.920000 4.43 5.63
PFC: mean embedded state per trajectory at normalized positions 0, 0.125, ..., 1 (dimension 1):
0.000 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000
OUT_L 7.49 9.41 11.26 13.01 14.31 15.44 16.60 17.77 18.969999
IN_L 10.41 10.49 10.20 9.38 8.54 7.67 6.78 6.06 5.580000
OUT_R 4.57 5.43 6.03 6.31 5.88 5.23 4.72 4.09 3.950000
IN_R 12.73 12.33 11.54 10.68 9.44 8.31 7.78 7.77 7.980000
Result: reproduced (through 3H and 3I). Both distance tests that the criterion takes as the measure of the shapes reproduce. The embeddings and the mean state per trajectory are shown; no measure of a ring or a double V is computed, and the plotted embeddings are made of one thread per pass, not a clean ring or double V. The outcome carries the dependence of fig3h-inseq-outseq-distance on the gap filling (and its seed sensitivity, stated there).
figs6-high-dimensional-distances¶
The distance relationships of Figure 3H to 3J hold in the original high-dimensional neural space.
The same tests with Euclidean distances between population vectors (cal_PV_*.m), with the main alignment of each test.
for reg in ("CA1", "PFC"):
X, lab, tn, _ = ST[(reg, 2)]
i_, o_, p, ni, no = h_test(X, reg, CODE_ALIGN[reg])
print(f"{reg} 3H in neural space (aligned by {CODE_ALIGN[reg]}): InSeq {i_:.2f}, OutSeq {o_:.2f}, p = {p:.2g}")
for align in ("location", "phase"):
if align != CODE_ALIGN[reg]:
i2, o2, p2, _, _ = h_test(X, reg, align); print(f" variant aligned by {align}: InSeq {i2:.2f}, OutSeq {o2:.2f}, p = {p2:.2g}")
pl = phase_loc(X, lab, "methods")
print(f"{reg} 3I in neural space: phase {np.median(pl['phase']):.2f}, location {np.median(pl['location']):.2f}, p = {stats.ranksums(pl['phase'], pl['location']).pvalue:.2g}")
X, lab, XF, labF, Xs = FNJ[(reg, 0)]
real, shuf = fn_dist(X, lab, XF, labF), fn_dist(X, lab, Xs, lab)
rr, ss = real[np.isfinite(real)], shuf[np.isfinite(shuf)]
print(f"{reg} 3J in neural space: F-N' {np.median(rr):.2f}, shuffle {np.median(ss):.2f} ({100 * (np.median(rr) / np.median(ss) - 1):+.1f}%), p = {stats.ranksums(rr, ss).pvalue:.2g}")
CA1 3H in neural space (aligned by location): InSeq 52.15, OutSeq 45.04, p = 2.5e-12
variant aligned by phase: InSeq 47.95, OutSeq 47.72, p = 0.45 CA1 3I in neural space: phase 48.37, location 45.02, p = 9.3e-129 CA1 3J in neural space: F-N' 48.47, shuffle 53.39 (-9.2%), p = 7.5e-78
PFC 3H in neural space (aligned by phase): InSeq 64.27, OutSeq 70.61, p = 4.8e-12
variant aligned by location: InSeq 81.92, OutSeq 86.19, p = 1e-12 PFC 3I in neural space: phase 65.47, location 86.09, p = 0 PFC 3J in neural space: F-N' 73.10, shuffle 168.59 (-56.6%), p = 0
Result: partially reproduced. In the neural space, 3H holds (CA1 OutSeq closer, 45.04 against 52.15, p = 2.5e-12; PFC InSeq closer, 64.27 against 70.61, p = 4.8e-12) and 3I holds (PFC phase 65.47 against location 86.09; CA1 location 45.02 against phase 48.37, p = 9.3e-129). For 3J, PFC F-N' is far below the shuffle (56.6% closer) but CA1 F-N' is below it too (9.2% closer, p = 7.5e-78), so two of the three relationships hold. Every width and gap-filling variant gives the same three directions (CA1 3J 7.3% to 9.2% closer, p <= 2.7e-21).
Abstraction (Figure 4, S7)¶
fig4d-pfc-ccgp and fig4d-ca1-ccgp-chance¶
In PFC, the cross-condition generalization performance (CCGP) for trajectory phase, task sequence and environment is high and above trial-label shuffles for every decoder.
In CA1, CCGP for trajectory phase, task sequence and environment is at or below chance by the trial-label permutation test.
Clusters are the embedded states (seed 0) at 20-30% and 70-80% of the normalized position of each
trajectory (cal_all_dichotomies.m), in N' and, for the environment, in F. Decoders and permutation
test as nfn.ccgp_sets and nfn.svm_ccgp (1000 shuffles; p is the fraction of shuffles above the real
CCGP, p_ge counts ties too). The code's task-sequence test pairs take the two clusters at different
segments (CCGP_taskSeq.m); taking them at the same segment is a variant (100 shuffles, for time). The code's environment
task-sequence decoders come in identical pairs (rows 1 and 3, 2 and 4 share their training set and
F test set), as here.
CL = {}
CC = {}
for reg in ("CA1", "PFC"):
E, EF, _ = EMB[(reg, 0)]
CN, CF = nfn.clusters(E, ST[(reg, 2)][1]), nfn.clusters(EF, ST[(reg, 1)][1])
CL[reg] = (CN, CF)
print(reg, "states per cluster, N':", {f"{nfn.TRACKS[t]} {s}": len(v) for (t, s), v in CN.items()})
print(reg, "states per cluster, F:", {f"{nfn.TRACKS[t]} {s}": len(v) for (t, s), v in CF.items()})
t0 = time.time()
CC[reg] = nfn.run_ccgp(CN, CF, n_perm=1000, seed=601 if reg == "PFC" else 602)
print(f"{reg}: 16 decoders with 1000 shuffles each in {time.time() - t0:.0f} s")
print(CC[reg].round(3).to_string(index=False))
CCV = {reg: nfn.run_ccgp(*CL[reg], n_perm=100, seed=603, seq_same_phase=True) for reg in ("CA1", "PFC")}
for reg in ("CA1", "PFC"):
x = CCV[reg][CCV[reg].variable != "phase"]
print(f"variant, task-sequence clusters at the same segment, {reg}: CCGP {x.ccgp.round(3).tolist()}, p {x.p.round(3).tolist()}")
CA1 states per cluster, N': {'OUT_L a': 100, 'OUT_L b': 100, 'IN_L a': 120, 'IN_L b': 120, 'OUT_R a': 160, 'OUT_R b': 160, 'IN_R a': 120, 'IN_R b': 120}
CA1 states per cluster, F: {'OUT_L a': 60, 'OUT_L b': 60, 'IN_L a': 80, 'IN_L b': 80, 'OUT_R a': 120, 'OUT_R b': 120, 'IN_R a': 80, 'IN_R b': 80}
CA1: 16 decoders with 1000 shuffles each in 116 s
variable decoder ccgp null_95 null_median p p_ge frac_pred_1
phase test IN_R 0.500 0.521 0.500 0.135 0.846 0.000
phase test OUT_R 0.500 0.528 0.500 0.179 0.842 1.000
phase test IN_L 0.500 0.508 0.500 0.072 0.955 0.000
phase test OUT_L 0.500 0.510 0.500 0.081 0.920 1.000
sequence test OUT_L/OUT_R a 0.547 0.538 0.500 0.022 0.022 0.942
sequence test IN_R/IN_L a 0.275 0.500 0.500 1.000 1.000 0.242
sequence test OUT_L/OUT_R b 0.578 0.536 0.500 0.002 0.002 0.904
sequence test IN_R/IN_L b 0.083 0.500 0.500 1.000 1.000 0.417
environment phase, train OUT_L+IN_L+OUT_R 0.446 0.533 0.500 0.997 0.998 0.346
environment phase, train OUT_L+IN_L+IN_R 0.780 0.539 0.502 0.000 0.000 0.466
environment phase, train OUT_L+OUT_R+IN_R 0.435 0.535 0.498 1.000 1.000 0.308
environment phase, train IN_L+OUT_R+IN_R 0.568 0.532 0.498 0.000 0.002 0.707
environment sequence, held out OUT_L/OUT_R a 0.581 0.541 0.500 0.002 0.003 0.712
environment sequence, held out IN_R/IN_L a 0.577 0.500 0.500 0.000 0.000 0.389
environment sequence, held out OUT_L/OUT_R b 0.581 0.544 0.500 0.000 0.001 0.712
environment sequence, held out IN_R/IN_L b 0.577 0.500 0.500 0.000 0.000 0.389
PFC states per cluster, N': {'OUT_L a': 100, 'OUT_L b': 100, 'IN_L a': 120, 'IN_L b': 120, 'OUT_R a': 160, 'OUT_R b': 160, 'IN_R a': 120, 'IN_R b': 120}
PFC states per cluster, F: {'OUT_L a': 60, 'OUT_L b': 60, 'IN_L a': 80, 'IN_L b': 80, 'OUT_R a': 120, 'OUT_R b': 120, 'IN_R a': 80, 'IN_R b': 80}
PFC: 16 decoders with 1000 shuffles each in 154 s
variable decoder ccgp null_95 null_median p p_ge frac_pred_1
phase test IN_R 0.500 0.550 0.500 0.482 0.565 0.667
phase test OUT_R 0.500 0.541 0.500 0.458 0.537 0.750
phase test IN_L 0.696 0.546 0.500 0.000 0.000 0.696
phase test OUT_L 0.500 0.535 0.500 0.195 0.795 0.200
sequence test OUT_L/OUT_R a 0.934 0.554 0.501 0.000 0.000 0.465
sequence test IN_R/IN_L a 0.758 0.500 0.500 0.000 0.000 0.258
sequence test OUT_L/OUT_R b 0.788 0.552 0.500 0.000 0.000 0.646
sequence test IN_R/IN_L b 0.500 0.500 0.500 0.000 1.000 0.000
environment phase, train OUT_L+IN_L+OUT_R 0.560 0.533 0.500 0.002 0.004 0.633
environment phase, train OUT_L+IN_L+IN_R 0.680 0.539 0.500 0.000 0.000 0.598
environment phase, train OUT_L+OUT_R+IN_R 0.560 0.533 0.500 0.002 0.004 0.540
environment phase, train IN_L+OUT_R+IN_R 0.541 0.534 0.498 0.019 0.020 0.638
environment sequence, held out OUT_L/OUT_R a 0.675 0.547 0.500 0.000 0.000 0.594
environment sequence, held out IN_R/IN_L a 0.877 0.500 0.500 0.000 0.000 0.381
environment sequence, held out OUT_L/OUT_R b 0.675 0.547 0.500 0.000 0.000 0.594
environment sequence, held out IN_R/IN_L b 0.877 0.500 0.500 0.000 0.000 0.381
variant, task-sequence clusters at the same segment, CA1: CCGP [0.578, 0.092, 0.547, 0.267, 0.446, 0.78, 0.435, 0.568, 0.581, 0.577, 0.581, 0.577], p [0.0, 1.0, 0.06, 1.0, 0.99, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0] variant, task-sequence clusters at the same segment, PFC: CCGP [0.788, 0.758, 0.934, 0.5, 0.56, 0.68, 0.56, 0.541, 0.675, 0.877, 0.675, 0.877], p [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.02, 0.0, 0.0, 0.0, 0.0]
x = CC["PFC"]
x = x.assign(constant=(x.frac_pred_1 == 0) | (x.frac_pred_1 == 1))
print(f"PFC decoders with p < 0.001: {(x.p < 0.001).sum()} of {len(x)}; with p < 0.05: {(x.p < 0.05).sum()}")
print(f"decoders predicting one class for every test point (CCGP 0.5, every shuffle ties, p = 0 by the strict count): {int(x.constant.sum())}")
print(f"counting ties (p_ge): p < 0.001 for {(x.p_ge < 0.001).sum()} of {len(x)}; p < 0.05 for {(x.p_ge < 0.05).sum()}")
print(x.groupby("variable").agg(mean_ccgp=("ccgp", "mean"), min_ccgp=("ccgp", "min"), n_p_ge_below_0001=("p_ge", lambda v: int((v < 0.001).sum())), n=("p", "size")).round(3))
print("decoders failing (p_ge >= 0.001):", x[x.p_ge >= 0.001][["variable", "decoder", "ccgp", "null_95", "p", "p_ge", "frac_pred_1"]].round(3).to_string(index=False))
PFC decoders with p < 0.001: 10 of 16; with p < 0.05: 13
decoders predicting one class for every test point (CCGP 0.5, every shuffle ties, p = 0 by the strict count): 1
counting ties (p_ge): p < 0.001 for 9 of 16; p < 0.05 for 12
mean_ccgp min_ccgp n_p_ge_below_0001 n
variable
environment 0.681 0.541 5 8
phase 0.549 0.500 1 4
sequence 0.745 0.500 3 4
decoders failing (p_ge >= 0.001): variable decoder ccgp null_95 p p_ge frac_pred_1
phase test IN_R 0.500 0.550 0.482 0.565 0.667
phase test OUT_R 0.500 0.541 0.458 0.537 0.750
phase test OUT_L 0.500 0.535 0.195 0.795 0.200
sequence test IN_R/IN_L b 0.500 0.500 0.000 1.000 0.000
environment phase, train OUT_L+IN_L+OUT_R 0.560 0.533 0.002 0.004 0.633
environment phase, train OUT_L+OUT_R+IN_R 0.560 0.533 0.002 0.004 0.540
environment phase, train IN_L+OUT_R+IN_R 0.541 0.534 0.019 0.020 0.638
Result: blocked, methods underspecified (the smoothing width of the single-trial maps). With 8-cm maps, 10 of the 16 PFC decoders have p < 0.001 by the code's strict count, but one of them (the task-sequence decoder tested on the inbound pair at the late segment) predicts one class for every test point, so its CCGP is 0.500, every shuffle ties and p = 1 counting ties. Counting ties, 9 of 16 have p < 0.001. The failing decoders are the phase decoders tested on IN_R, OUT_R and OUT_L (CCGP 0.500 each, p >= 0.54 counting ties), that sequence decoder, and three environment phase decoders (0.560, 0.560 and 0.541; p = 0.004, 0.004 and 0.020), so the claim does not reproduce at 8 cm. Above every one of 100 shuffles (cell 79) are 14 of 16 decoders at 4 cm and 15 of 16 at 2 cm, and 4 of 16 without gap filling at 8 cm and at 2 cm. In this run no width gives all 16. In the verifier's round-2 copy, embedded on another machine before the CPU target was fixed, all 16 were above every one of 200 shuffles at 2 cm, so at the narrow widths the result sits at the edge of the criterion and changes with the embedding. The Methods give no smoothing for the single-trial maps and the session maps' width is garbled; at 8 cm the claim fails clearly, at 2 cm it can pass, so the claim stays methods underspecified, as the claims file's rule for the width (fig2d-pv-similarity) has it.
x = CC["CA1"]
x = x.assign(constant=(x.frac_pred_1 == 0) | (x.frac_pred_1 == 1))
print(f"CA1 decoders with p < 0.05: {(x.p < 0.05).sum()} of {len(x)}; counting ties (p_ge): {(x.p_ge < 0.05).sum()}; "
f"decoders predicting one class for every test point: {int(x.constant.sum())}")
print(x.groupby("variable").agg(mean_ccgp=("ccgp", "mean"), max_ccgp=("ccgp", "max"), n_p_ge_below_005=("p_ge", lambda v: int((v < 0.05).sum())), n=("p", "size")).round(3))
print("decoders above chance (p_ge < 0.05):", x[x.p_ge < 0.05][["variable", "decoder", "ccgp", "null_95", "p", "p_ge", "frac_pred_1"]].round(3).to_string(index=False))
CA1 decoders with p < 0.05: 8 of 16; counting ties (p_ge): 8; decoders predicting one class for every test point: 4
mean_ccgp max_ccgp n_p_ge_below_005 n
variable
environment 0.568 0.780 6 8
phase 0.500 0.500 0 4
sequence 0.371 0.578 2 4
decoders above chance (p_ge < 0.05): variable decoder ccgp null_95 p p_ge frac_pred_1
sequence test OUT_L/OUT_R a 0.547 0.538 0.022 0.022 0.942
sequence test OUT_L/OUT_R b 0.578 0.536 0.002 0.002 0.904
environment phase, train OUT_L+IN_L+IN_R 0.780 0.539 0.000 0.000 0.466
environment phase, train IN_L+OUT_R+IN_R 0.568 0.532 0.000 0.002 0.707
environment sequence, held out OUT_L/OUT_R a 0.581 0.541 0.002 0.003 0.712
environment sequence, held out IN_R/IN_L a 0.577 0.500 0.000 0.000 0.389
environment sequence, held out OUT_L/OUT_R b 0.581 0.544 0.000 0.001 0.712
environment sequence, held out IN_R/IN_L b 0.577 0.500 0.000 0.000 0.389
Result: blocked, methods underspecified (depends on fig4d-pfc-ccgp). The claims file lists
fig4d-pfc-ccgp under this claim's depends_on, and the protocol (Author, step 1.5) blocks a claim that
depends on a blocked claim for the same reason. Its own result does not change with the width: at
8 cm 4 of the 16 CA1 decoders predict one class for every test point (the four phase decoders, CCGP
0.500), and 8 decoders have p < 0.05 counting ties: two sequence decoders (CCGP 0.547 and 0.578) and
six environment decoders (0.568 to 0.780), against 95th percentiles of 0.500 to 0.544. The width and
gap-filling variants give 4 to 11 (100 shuffles). Tested alone, the claim would not reproduce.
figs7b-geometric-null¶
Against a geometric random model that keeps within-condition noise but destroys the geometry, PFC CCGP is significant while CA1 CCGP is not.
The noise cloud of each cluster (its states about the centroid) has its axes permuted, separately per cluster, 100 times, and CCGP is recomputed. Readings: (a) the 3 UMAP coordinates; (b) the neuron axes of the cluster's population vectors, embedded with the N' fit. The CCGP of a variable is the mean over its decoders (an assumption; the legend plots the decoders with their median), compared with the 95th percentile of the same mean under the null.
from joblib import Parallel, delayed
def null_b_batch(um, nCN, nCF, seed, n):
"""Reading (b): permute the neuron axes of each cluster's noise cloud, embed with the N' fit."""
rng = np.random.default_rng(seed); keys = list(nCN); out = []
for _ in range(n):
pN, pF = nfn.permute_noise(nCN, rng), nfn.permute_noise(nCF, rng)
allN = um.transform(np.vstack([pN[k] for k in keys])); allF = um.transform(np.vstack([pF[k] for k in keys]))
sN = np.cumsum([0] + [len(pN[k]) for k in keys]); sF = np.cumsum([0] + [len(pF[k]) for k in keys])
eN = {k: allN[sN[i]:sN[i + 1]] for i, k in enumerate(keys)}; eF = {k: allF[sF[i]:sF[i + 1]] for i, k in enumerate(keys)}
out.append(nfn.ccgp_values(eN, eF).groupby("variable").ccgp.mean())
return out
def cluster_rows(lab):
return {(t, s): (lab[:, 0] == t) & np.isin(lab[:, 2], b) for t in nfn.TRACKS for s, b in (("a", nfn.EARLY_BINS), ("b", nfn.LATE_BINS))}
GN = {}
for reg in ("CA1", "PFC"):
CN, CF = CL[reg]
real = nfn.ccgp_values(CN, CF).groupby("variable").ccgp.mean()
rng = np.random.default_rng(701 if reg == "CA1" else 702)
draws = [(nfn.permute_noise(CN, rng), nfn.permute_noise(CF, rng)) for _ in range(100)] # drawn in order, fitted in parallel
null_a = pd.DataFrame(Parallel(n_jobs=4)(delayed(lambda pN, pF: nfn.ccgp_values(pN, pF).groupby("variable").ccgp.mean())(pN, pF) for pN, pF in draws))
X, lab, _, _ = ST[(reg, 2)]; XF, labF, _, _ = ST[(reg, 1)]
um = EMB[(reg, 0)][2]
rN, rF = cluster_rows(lab), cluster_rows(labF)
nCN = {k: X[m] for k, m in rN.items()}; nCF = {k: XF[m] for k, m in rF.items()}
keys = list(nCN)
t0 = time.time()
batches = Parallel(n_jobs=4)(delayed(null_b_batch)(um, nCN, nCF, [703, reg == "PFC", i], 25) for i in range(4))
null_b = [row for b in batches for row in b]
null_b = pd.DataFrame(null_b)
GN[reg] = (real, null_a, null_b)
print(f"{reg} (reading b in {time.time() - t0:.0f} s)")
print(pd.DataFrame({"CCGP (mean of decoders)": real, "null (a) 95th pct": null_a.quantile(0.95), "p (a)": (null_a >= real).mean(),
"null (b) 95th pct": null_b.quantile(0.95), "p (b)": (null_b >= real).mean()}).round(3))
CA1 (reading b in 78 s)
CCGP (mean of decoders) null (a) 95th pct p (a) null (b) 95th pct p (b)
variable
environment 0.568 0.586 0.21 0.597 0.26
phase 0.500 0.500 1.00 0.505 0.56
sequence 0.371 0.592 0.54 0.616 0.54
PFC (reading b in 97 s)
CCGP (mean of decoders) null (a) 95th pct p (a) null (b) 95th pct p (b)
variable
environment 0.681 0.715 0.29 0.722 0.48
phase 0.549 0.616 0.59 0.679 0.76
sequence 0.745 0.755 0.08 0.684 0.00
Result: blocked, methods underspecified (depends on fig4d-pfc-ccgp). The claim depends on fig4d-pfc-ccgp, which is blocked. At 8 cm, under reading (a) no PFC variable's CCGP exceeds the 95th percentile of its null (p 0.08 to 0.59); under reading (b) the PFC task-sequence CCGP does (0.745 against 0.684, no null draw above it), and environment and phase do not (p = 0.48 and 0.76). No CA1 variable exceeds its null (p 0.21 to 1.00 in reading a, 0.26 to 0.56 in reading b). The null was not run at the other widths. Tested alone, the claim would not reproduce.
figs7d-dichotomies-decodable¶
With traditional cross-validated decoders, most of the 35 dichotomies of the 8 clusters are decoded accurately in both CA1 and PFC, with similar accuracy in the two regions.
The 35 balanced dichotomies of the 8 N' clusters, each decoded by a linear C-SVM with 4-fold cross-validation on the 3 UMAP coordinates. The trial-label permutation test uses 50 shuffles in place of the paper's 1,000 (a deviation, for time: 100 took about 3 minutes per region here, and the notebook must run in 30 minutes); p is resolved to 0.02, enough to decide p < 0.05 (at most 2 of 50 shuffles above). Rank-sum test of PFC against CA1 accuracies.
from joblib import Parallel, delayed
DI = {}
for reg in ("CA1", "PFC"):
CN = CL[reg][0]
t0 = time.time()
res = Parallel(n_jobs=4)(delayed(nfn.decode_dichotomy)(CN, g, 50, [801, i, reg == "PFC"]) for i, g in enumerate(nfn.dichotomies()))
DI[reg] = pd.DataFrame({"accuracy": [r[0] for r in res], "p": [r[1] for r in res], "p_ge": [r[3] for r in res], "null_95": [np.percentile(r[2], 95) for r in res]})
x = DI[reg]
print(f"{reg} ({time.time() - t0:.0f} s): accuracy median {x.accuracy.median():.3f} (range {x.accuracy.min():.3f} to {x.accuracy.max():.3f}); "
f"p < 0.05 for {(x.p < 0.05).sum()} of 35 ({(x.p_ge < 0.05).sum()} counting ties); null 95th percentile median {x.null_95.median():.3f}")
print(f"rank-sum PFC against CA1 accuracies: p = {stats.ranksums(DI['PFC'].accuracy, DI['CA1'].accuracy).pvalue:.3g}; "
f"median difference {DI['PFC'].accuracy.median() - DI['CA1'].accuracy.median():+.3f}")
CA1 (73 s): accuracy median 0.830 (range 0.655 to 1.000); p < 0.05 for 35 of 35 (35 counting ties); null 95th percentile median 0.537
PFC (94 s): accuracy median 0.631 (range 0.485 to 0.868); p < 0.05 for 34 of 35 (34 counting ties); null 95th percentile median 0.534 rank-sum PFC against CA1 accuracies: p = 8.63e-09; median difference -0.199
Result: blocked, methods underspecified (depends on fig4d-pfc-ccgp). The claim depends on fig4d-pfc-ccgp, which is blocked. At 8 cm most dichotomies are decoded above the shuffles in both regions (CA1 35 of 35, PFC 34 of 35, the same counting ties), but the accuracies differ: CA1 is higher (median 0.830 against 0.631, rank-sum p = 8.6e-09), so the similarity part would fail; the other widths were not run for this test.
print(f"notebook ran in {(time.time() - started) / 60:.1f} min")
notebook ran in 24.9 min
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome |
|---|---|
| fig1b-faster-novel-learning | blocked: data missing |
| fig1-denovo-final-performance | blocked: data missing |
| fig1c-session-performance | blocked: data missing |
| fig1c-familiar-above-novel | not reproduced (N-F-N' rats only) |
| figs1b-first-session-inbound-outbound | blocked: data missing |
| cells-per-session | reproduced |
| figs3a-ca1-global-remapping | reproduced |
| figs3a-pfc-similar-ca1-lower | partially reproduced |
| fig2d-pv-similarity | reproduced |
| fig2e-ca1-error-familiar | partially reproduced |
| fig2e-pfc-error-stable | partially reproduced |
| fig2f-switch-dynamics | partially reproduced (CA1 rise under one reading of trial -1 of two) |
| figs3c-pfc-assimilation | reproduced |
| figs3d-pfc-similarity-over-trials | partially reproduced (CA1 stable under one reading of three) |
| figs4-ramping-fraction | not reproduced |
| figs4-ramping-generalization | blocked: too few ramping cells |
| figs4d-reward-location-similarity | blocked: methods underspecified (smoothing width) |
| figs5d-path-equivalence | reproduced |
| figs5e-path-equivalence-stable | blocked: methods underspecified (smoothing width) |
| figs5f-path-equivalence-vs-selectivity | blocked: methods underspecified |
| fig3d-selectivity-index | blocked: methods underspecified (smoothing width) |
| fig3e-selectivity-stable | blocked: methods underspecified (smoothing width) |
| manifold-cell-counts | reproduced |
| fig3fg-manifold-shapes | reproduced |
| fig3h-inseq-outseq-distance | reproduced (assumption: empty bins filled in the single-trial maps) |
| fig3i-phase-location-distance | reproduced |
| fig3j-familiar-novel-manifold | partially reproduced (assumption: empty bins filled) |
| figs6-high-dimensional-distances | partially reproduced |
| fig4d-pfc-ccgp | blocked: methods underspecified (smoothing width of the single-trial maps) |
| fig4d-ca1-ccgp-chance | blocked: methods underspecified (depends on fig4d-pfc-ccgp) |
| figs7b-geometric-null | blocked: methods underspecified (depends on fig4d-pfc-ccgp) |
| figs7d-dichotomies-decodable | blocked: methods underspecified (depends on fig4d-pfc-ccgp) |
How the analysis differs from the paper's¶
- Linearization. The dandiset has no linearized positions; positions are projected onto W-shaped paths built from the well positions and the run positions (an assumption; 180 to 220 cm per trajectory, median 5 cm from the path).
- The 1D smoothing width is garbled in the paper's text, and no smoothing is given for the single-trial maps; 8, 4 and 2 cm are run. Where the outcome changes with the width (figs4d, figs5e, fig3d, fig3e, fig4d-pfc-ccgp), the claim is blocked as methods underspecified, by the claims file's rule for fig2d-pv-similarity, and the claims that depend on fig4d-pfc-ccgp are blocked with it.
- Sharp-wave ripple times are not excluded from the rate maps (no ripple times in the dandiset), and the decoder's theta restriction (a variant in the code) cannot be applied (no theta times).
- The code's per-animal interneuron lists cannot be applied (no cell numbers); the code's 7 Hz rule for CA1 is.
- Single-trial maps for the manifolds have the empty 1-cm bins inside a pass filled linearly (an assumption, as cal_PVsimilarity.m fills session maps). Without the filling, 3H, 3I and 3J keep their outcomes at 8 cm, but at 2 cm the PFC parts of 3H and 3J fail (in this run), and the CCGP changes.
- A CCGP decoder that predicts one class for every test point has CCGP 0.5 and ties with every shuffle; the code's strict count of shuffles above gives it p = 0. Here such decoders are counted as not significant (ties counted).
- The UMAP implementation is umap-learn with fixed seeds, not the MATLAB toolbox with the paper's saved templates; the manifold results vary with the seed (seeds 0 to 2 shown). The embeddings also differed between machines until numba's CPU target was fixed (header cell); the verifier's round-2 copy and this run differ in the decimals and, at the edges, in single tests (stated per claim).
- Figure 1 and S1 behaviour is tested only on N, F and N' of the five rats; the de novo sessions and the F-N-F'-N' rats are not in the dandiset.
- The dichotomy decoders use 50 label shuffles in place of 1,000, and the manifold variants' CCGP 100, for the 30-minute limit.
What was found about the dataset¶
- The units table's
electrodescolumn holds tetrode numbers, not row indices of the electrodes table; read as row indices, units are put in the wrong region (up to 32 units per animal). - Units have no cell number within a tetrode, no waveform and no interneuron label, so the paper's interneuron exclusion cannot be applied; units clustered only in later epochs are appended to the table with no spike in N, without a note.
acquisition/ElectricalSeriesis a byte-identical second copy of the 1.5 kHz LFP; there is no 30 kHz recording.- The position's reference frame is given as "center well", but positions are camera coordinates in which the center well sits at different places in each room. The five files' session start times fall within 16 minutes of each other on 2023-03-14, two days before the first published version: the time of conversion, not of the recordings. Epochs carry no names; N, F, N' are known only from the description's order. Linearized positions, the idealized track paths, ripple and theta times are not included.
- The de novo learning sessions and the F-N-F'-N' rats behind Figure 1C, 1B and S1 are not in the dandiset.