Reproducing Petersen & Buzsaki (2020) from dandiset 000059¶
Paper. Petersen PC, Buzsaki G (2020). Cooling of medial septum reveals theta phase lag coordination of hippocampal cell assemblies. Neuron 107, 731-744. doi:10.1016/j.neuron.2020.05.023. This attempt used the full text of the accepted author manuscript in PubMed Central (PMC7442698). The supplementary PDF could not be obtained, so supplementary figures are tested only where the main text states their result.
Dandiset. DANDI:000059, version 0.250624.0444. The paper cites no DANDI version (its data statement points to the lab's website and CRCNS). This is the most recent published version. Version 0.230907.2101 holds the same 100 files with the same blobs. The first version, 0.210812.1514, held 54 combined files, among them four sessions without behaviour, and its files have no trials table, which every behavioural claim here needs; the data check shows this.
Code. The paper's MATLAB code is public (GitHub petersenpeter/Code_Petersen_Buzsaki_Neuron_2020,
commit aec2f7e, read for this attempt, not run). It fixes the theta measures, the per-trial
averages, the onset-aligned trial series behind Figures 2G, 2H and 6, and the lists of sessions;
cooling.py says which parameter comes from where.
What the dandiset holds. 50 behavioural sessions of five rats on the figure-8 maze, each with a processed file (trials with cooling state and choice error, position, speed, MS temperature, and units in 24 sessions) and a raw file (20 kHz wide-band and 1250 Hz LFP). The files carry no cell types and no place field definitions, and no hippocampal temperature recordings.
The notebook streams every file at the pinned version. One bootstrap draws random numbers; it has its own seeded generator.
| Claim | Status |
|---|---|
| session-structure, cooling-time-course | tested |
| fig2g-theta-freq-temperature, fig2h-theta-power-temperature, fig2j-frequency-consistent-across-animals, fig2i-theta-pre-cooling-post, disc-theta-frequency-12-percent | tested |
| fig6d-errors-increase, fig6e-errors-recover, fig6c-worst-trials-4-to-8, fig6f-errors-temperature, fig6g-errors-theta-frequency, fig6h-errors-theta-power, fig6i-errors-not-speed | tested |
| fig4e-speed-decreases-within-session, fig5-speed-control-cooling | tested |
| suppl2-place-cell-rates-decrease, disc-interneuron-rates-decrease, suppl2d-interneuron-oscillation-slows, suppl1g-waveforms-unchanged | tested, with putative cell types made here |
| fig4b-fp-flfp-decrease, fig4c-fp-flfp-correlation, fig4ac-fp-flfp-difference-preserved, fig4a-frequency-temperature-slopes, fig4d-precession-slope-decrease, fig4gh-fp-speed | tested, with place fields defined here |
| fig3c-field-size-correlated, fig3d-size-fp-negative, fig3e-model-field-size, fig3f-information-unchanged, fig3hi-population-vector-unchanged | tested, with place fields defined here |
| fig5d-ccg-period, fig5cg-compression, fig5h-time-lags-increase, fig5i-phase-lags-unchanged, fig5j-fraction-active-unchanged, fig5k-spikes-per-cycle-unchanged | tested |
| suppl1-hippocampal-temperature | blocked: the hippocampal temperature recordings are not in the dandiset |
In [1]:
import os
os.environ["OMP_NUM_THREADS"] = "1"
os.environ["OPENBLAS_NUM_THREADS"] = "1"
import time
import warnings
from concurrent.futures import ProcessPoolExecutor
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import cooling as C
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
In [2]:
started = time.time()
S = C.sessions()
with ProcessPoolExecutor(6) as ex:
beh = dict(zip(S.index, ex.map(C.behaviour_session, S.index)))
unit_keys = [k for k, r in beh.items() if S.loc[k, "in_code_list"] and r["has_state"] and r["n_good_units"] > 0]
with ProcessPoolExecutor(4) as ex:
units = dict(zip(unit_keys, ex.map(C.unit_session, unit_keys)))
In [3]:
print(f"dandiset {C.DANDISET_ID} version {C.DANDISET_VERSION}: {len(S)} sessions, "
f"loaded in {time.time() - started:.0f} s ({len(units)} with unit analyses)")
dandiset 000059 version 0.250624.0444: 50 sessions, loaded in 509 s (13 with unit analyses)
Data check¶
Each check of the description against the values is stated under its output.
In [4]:
from dandi.dandiapi import DandiAPIClient
with DandiAPIClient() as client:
vers = {}
for v in ["0.210812.1514", "0.230907.2101", "0.250624.0444"]:
assets = list(client.get_dandiset(C.DANDISET_ID, v).get_assets())
vers[v] = {a.path: (a.size, getattr(a, "blob", None)) for a in assets}
old_asset = next(a for a in client.get_dandiset(C.DANDISET_ID, "0.210812.1514").get_assets()
if "MS10-170319" in a.path)
old_url = old_asset.get_content_url(follow_redirects=1, strip_query=True)
for v, a in vers.items():
print(f"{v}: {len(a)} files, {sum(s for s, _ in a.values()) / 1e12:.2f} TB")
same = vers["0.230907.2101"].keys() == vers["0.250624.0444"].keys() and all(
vers["0.230907.2101"][p] == vers["0.250624.0444"][p] for p in vers["0.250624.0444"])
print("0.230907.2101 and 0.250624.0444 hold the same paths, sizes and blobs:", same)
old_sessions = sorted({C.session_key(p) for p in vers["0.210812.1514"]})
new_sessions = sorted({C.session_key(p) for p in vers["0.250624.0444"]})
print("sessions only in 0.210812.1514:", sorted(set(old_sessions) - set(new_sessions)))
f_old = C.open_h5(old_url)
print("0.210812.1514, file", old_asset.path.split("/")[-1], "has a trials table:", "intervals/trials" in f_old,
"; has units:", "units" in f_old, "; has LFP:", "processing/ecephys/LFP" in f_old)
f_old.close()
0.210812.1514: 54 files, 3.26 TB 0.230907.2101: 100 files, 2.94 TB 0.250624.0444: 100 files, 2.94 TB 0.230907.2101 and 0.250624.0444 hold the same paths, sizes and blobs: True sessions only in 0.210812.1514: ['MS12-170712-101120', 'MS13-171204-102904', 'MS13-171205-110831', 'MS13-171206-132039']
0.210812.1514, file sub-MS10_ses-Peter-MS10-170319-180352-concat_behavior+ecephys.nwb has a trials table: False ; has units: True ; has LFP: True
The most recent version and the one before it are the same files. The first version holds four more sessions (MS12-170712, MS13-171204, -171205, -171206), none with behaviour, and its files have no trials table, so the pinned version is the only one with the cooling states and the choice errors.
In [5]:
rows = []
for k, r in beh.items():
tr = r.get("trials")
states = tr.state.value_counts().to_dict() if tr is not None and r["has_state"] else {}
rows.append({"session": k, "subject": r["subject"], "in_code_list": S.loc[k, "in_code_list"],
"fig5_list": S.loc[k, "in_fig5_list"], "code_control": S.loc[k, "code_control"],
"trials": r["n_trials"], "has_state": r["has_state"],
"pre": states.get("Pre-Cooling", 0), "cool_on": states.get("Cooling on", 0),
"cool_off": states.get("Cooling off", 0), "post": states.get("Post-Cooling", 0),
"errors": int(tr.error.sum()) if tr is not None else 0,
"good_units": r["n_good_units"], "all_units": sum(r["unit_quality_counts"].values()),
"frame_ms": round(1000 * r["frame_dt"], 2), "theta_ch": r["theta_channel"], "n_el": r["n_electrodes"],
"T_min": round(r["temp_min"], 2), "T_max": round(r["temp_max"], 2), "T_d2_sd": r["temp_d2_sd"],
"raw_GB": round(S.loc[k, "raw_size_gb"], 1)})
info = pd.DataFrame(rows).set_index("session")
print("Sessions per subject:", info.groupby("subject").size().to_dict(), "; total", len(info))
print("Paper: 53 sessions in 5 rats (Figure 2). The code's behavioural list (BatchProcessing.m) has",
sum(len(v) for v in C.CODE_SESSIONS.values()), "sessions;", int(info.in_code_list.sum()), "of them are in the dandiset.")
missing = [f"{a}-{d}" for a, ds in C.CODE_SESSIONS.items() for d in ds if f"{a}-{d}" not in info.index]
print("Code sessions not in the dandiset:", missing)
print("Dandiset sessions not in the code's behavioural list:", info.index[~info.in_code_list].tolist())
print("Code's Figure 5 sessions in the dandiset:", int(info.fig5_list.sum()), "of", len(C.FIG5_SESSIONS))
print("\nTrial-table columns that differ from (condition, cooling state, error):",
{k: r["trial_columns"] for k, r in beh.items() if set(r["trial_columns"]) != {"condition", "cooling state", "error", "id", "start_time", "stop_time"}})
print("\nUnit qualities over all files:", pd.DataFrame([r["unit_quality_counts"] for r in beh.values()]).sum().astype(int).to_dict())
print("Sessions whose units are all 'unsorted':", [k for k, r in beh.items() if r["unit_quality_counts"] and set(r["unit_quality_counts"]) == {"unsorted"}])
print("\nSubject metadata:", {r["subject"]: (r["subject_meta"].get("age"), r["subject_meta"].get("weight"), r["subject_meta"].get("sex"))
for r in beh.values()})
print()
print(info.drop(columns=["T_d2_sd"]).to_string())
Sessions per subject: {'MS10': 9, 'MS12': 8, 'MS13': 5, 'MS21': 15, 'MS22': 13} ; total 50
Paper: 53 sessions in 5 rats (Figure 2). The code's behavioural list (BatchProcessing.m) has 53 sessions; 47 of them are in the dandiset.
Code sessions not in the dandiset: ['MS12-170715-111545', 'MS12-170712-101120', 'MS13-171206-132039', 'MS13-171205-110831', 'MS13-171204-102904', 'MS21-180806-111316']
Dandiset sessions not in the code's behavioural list: ['MS10-170307-154746', 'MS21-180719-155941', 'MS22-180719-122813']
Code's Figure 5 sessions in the dandiset: 17 of 18
Trial-table columns that differ from (condition, cooling state, error): {'MS12-170719-095305': ['condition', 'cooling state', 'error', 'id', 'start_time', 'stop_time', 'temperature'], 'MS21-180627-143449': ['condition', 'error', 'id', 'start_time', 'stop_time'], 'MS21-180712-103200': []}
Unit qualities over all files: {'unsorted': 4664, 'noise': 13188, 'good': 1269, 'mua': 2995, 'good2': 203}
Sessions whose units are all 'unsorted': ['MS10-170307-154746', 'MS10-170317-153237', 'MS10-170319-180352', 'MS10-170320-212321', 'MS10-170321-164406', 'MS10-170323-173246', 'MS12-170723-141121']
Subject metadata: {'MS10': ('P3M/P6M', '395g', 'M'), 'MS12': ('P3M/P6M', '410g', 'M'), 'MS13': ('P3M/P6M', '320g', 'M'), 'MS21': ('P3M/P6M', '250g', 'M'), 'MS22': ('P3M/P6M', '260g', 'M')}
subject in_code_list fig5_list code_control trials has_state pre cool_on cool_off post errors good_units all_units frame_ms theta_ch n_el T_min T_max raw_GB
session
MS10-170307-154746 MS10 False True False 37 True 21 10 0 6 7 0 409 25.00 46 64 37.00 37.00 41.6
MS10-170311-180956 MS10 True False False 130 True 43 31 0 56 16 0 0 25.00 46 64 17.00 37.00 19.9
MS10-170313-122631 MS10 True False False 143 True 39 47 0 57 22 0 0 25.00 46 64 16.98 37.00 11.8
MS10-170314-163038 MS10 True False False 142 True 52 47 0 43 16 0 0 25.00 46 47 17.00 37.00 2.8
MS10-170317-153237 MS10 True True False 201 True 52 55 1 93 11 0 395 25.00 46 64 17.00 37.00 12.2
MS10-170319-180352 MS10 True False False 204 True 49 92 0 63 16 0 381 25.00 52 64 17.00 37.00 14.5
MS10-170320-212321 MS10 True False False 179 True 53 60 0 66 17 0 365 33.35 52 64 17.00 37.00 10.2
MS10-170321-164406 MS10 True False False 130 True 52 54 0 24 20 0 344 8.35 52 64 17.00 37.00 14.3
MS10-170323-173246 MS10 True False False 154 True 48 67 9 30 14 0 361 8.35 52 64 17.00 37.00 3.9
MS12-170713-093400 MS12 True False False 115 True 32 51 13 19 5 0 0 8.35 8 128 19.73 37.19 117.0
MS12-170714-122034 MS12 True True False 151 True 39 57 19 36 7 0 0 8.35 8 128 18.49 50.16 72.5
MS12-170716-172307 MS12 True True False 146 True 41 52 17 36 13 36 915 8.35 19 128 18.58 37.44 80.5
MS12-170717-111614 MS12 True True False 134 True 39 54 15 26 7 0 0 8.35 67 128 18.75 36.96 107.9
MS12-170718-101720 MS12 True False False 125 True 38 60 10 17 3 0 0 8.35 19 128 19.57 36.57 119.1
MS12-170719-095305 MS12 True True False 149 True 35 47 26 41 5 48 1114 8.35 19 128 20.06 36.92 115.9
MS12-170720-093528 MS12 True False False 155 True 39 32 15 69 10 0 0 8.35 19 128 20.53 36.84 100.6
MS12-170723-141121 MS12 True False False 140 True 32 42 12 54 6 0 468 8.35 19 128 23.07 36.82 108.5
MS13-171127-133330 MS13 True False False 160 True 42 35 34 49 5 0 0 8.35 20 96 17.87 35.95 49.8
MS13-171128-113924 MS13 True True False 159 True 62 30 17 50 8 73 1014 8.35 20 128 18.50 37.82 79.7
MS13-171129-105507 MS13 True True False 121 True 42 41 9 29 5 0 0 8.35 20 119 19.10 37.57 78.2
MS13-171130-121758 MS13 True True False 174 True 33 52 30 59 13 0 0 8.35 82 119 18.54 37.55 65.2
MS13-171201-130527 MS13 True False False 133 True 44 24 9 56 7 67 928 8.35 20 128 20.27 37.62 93.5
MS21-180621-141817 MS21 True False False 125 True 33 27 8 57 5 23 584 8.35 63 128 18.24 40.40 9.5
MS21-180622-174411 MS21 True False False 131 True 36 20 17 58 8 0 0 8.35 63 128 20.65 40.81 7.8
MS21-180623-170717 MS21 True False False 122 True 39 33 14 36 9 0 0 8.35 63 128 18.18 45.48 7.9
MS21-180625-153927 MS21 True True False 128 True 35 36 14 43 23 83 1024 8.35 63 128 16.73 37.35 38.9
MS21-180627-143449 MS21 True True False 179 False 0 0 0 0 33 160 1730 8.35 79 128 17.91 37.31 41.2
MS21-180628-155921 MS21 True True False 168 True 39 29 24 76 13 109 854 8.35 63 128 20.23 37.36 23.8
MS21-180629-110332 MS21 True True False 169 True 40 50 20 59 16 172 1827 8.35 63 128 18.99 37.21 50.0
MS21-180711-112927 MS21 True False False 129 True 36 38 15 40 13 0 0 8.35 63 128 18.07 37.51 70.2
MS21-180712-103200 MS21 True True False 0 False 0 0 0 0 0 59 855 8.35 63 128 25.05 36.95 94.3
MS21-180719-155941 MS21 False False True 132 True 39 36 11 46 11 139 1384 8.35 40 128 17.00 37.00 19.0
MS21-180720-110037 MS21 True False False 187 True 38 55 14 80 17 0 0 8.35 63 128 20.74 45.80 124.0
MS21-180724-144402 MS21 True False False 192 True 45 71 22 54 13 0 0 8.35 63 128 20.17 37.76 96.9
MS21-180730-120605 MS21 True False False 138 True 36 23 16 63 5 0 0 8.35 63 128 22.50 38.46 78.9
MS21-180802-120147 MS21 True False False 152 True 41 20 14 77 15 0 0 8.35 63 128 21.80 38.23 74.6
MS21-180803-142301 MS21 True False False 121 True 33 48 10 30 13 0 0 8.35 63 128 19.92 37.76 56.1
MS22-180621-131028 MS22 True False False 115 True 31 25 17 42 14 27 476 8.35 63 128 22.38 45.28 12.6
MS22-180622-115734 MS22 True False False 80 True 29 31 5 15 16 0 0 8.35 63 128 21.60 49.94 56.7
MS22-180623-170141 MS22 True False False 113 True 28 34 8 43 23 0 0 8.35 63 128 21.24 41.95 16.5
MS22-180625-162932 MS22 True False False 116 True 27 44 17 28 18 0 0 8.35 63 128 21.85 37.75 20.0
MS22-180626-110916 MS22 True False False 134 True 30 31 14 59 22 0 0 8.35 63 128 23.34 37.80 112.7
MS22-180627-114856 MS22 True False False 162 True 33 45 13 71 19 0 0 8.35 63 128 23.49 37.32 48.5
MS22-180628-120341 MS22 True True False 159 True 36 59 7 57 24 79 1191 8.35 0 128 22.68 37.43 68.9
MS22-180629-110319 MS22 True True False 143 True 40 50 15 38 11 172 1509 8.35 6 128 23.11 37.33 51.8
MS22-180711-112912 MS22 True True False 111 True 39 29 14 29 14 68 1183 8.35 6 128 23.35 37.29 63.6
MS22-180712-102504 MS22 True False False 151 True 38 52 17 44 9 0 0 8.35 6 128 22.56 37.90 97.5
MS22-180719-122813 MS22 False False True 142 True 77 26 8 31 11 77 1664 8.35 55 128 17.00 37.00 59.9
MS22-180720-110055 MS22 True False False 141 True 45 47 14 35 6 80 1344 8.35 6 128 23.64 37.48 106.7
MS22-180724-144421 MS22 True False False 138 True 37 59 16 26 12 0 0 8.35 6 128 24.03 48.72 93.3
The dandiset holds 50 sessions; the paper reports 53 in five rats, and the code's behavioural list has 53. 47 of the code's sessions are here; six are not (MS12-170712, -170715, MS13-171204, -171205, -171206, MS21-180806). Three sessions here are not in that list: MS10-170307, and the two sessions of 2018-07-19, which the code's recording list describes as controls without cooling. MS21-180627 has a trials table without the cooling state, and MS21-180712 has no trials. Units are curated (good, good2, mua, noise) in 17 sessions; in seven (the six MS10 sessions with units and MS12-170723) every unit is 'unsorted', so no unit can be selected there. No file has a cell type, a place field, or an electrode location (every location is 'unknown'). The subject age is a range (P3M/P6M).
In [6]:
print("Frame interval of position, speed and temperature (ms):", sorted(set(info.frame_ms)),
"; Methods: OptiTrack at 120 Hz (8.33 ms)")
print("Position, speed and temperature share one timestamp array in every file:", all(r["same_clock"] for r in beh.values()))
print("Speed unit:", {r["speed_unit"] for r in beh.values()}, "; position conversion to metres:", {r["pos_conversion"] for r in beh.values()})
xr = np.array([r["x_range"] for r in beh.values()]); yr = np.array([r["y_range"] for r in beh.values()])
print(f"x range over sessions {xr[:, 0].min():.1f} to {xr[:, 1].max():.1f}, y range {yr[:, 0].min():.1f} to {yr[:, 1].max():.1f} (in cm)")
print("Speed percentiles 1, 50, 99 (cm/s), median over sessions:", np.round(np.median([r["speed_pct"] for r in beh.values()], axis=0), 1).tolist())
lin = pd.DataFrame({k: u["lin"] for k, u in units.items()}).T
print("\nMaze fitted per unit session: stem length L and circle radius R (cm)")
print(lin.round(1).to_string())
print("Diameter 2R:", round(2 * lin.R.min(), 1), "to", round(2 * lin.R.max(), 1), "cm; Methods: 110 cm diameter")
Frame interval of position, speed and temperature (ms): [8.35, 25.0, 33.35] ; Methods: OptiTrack at 120 Hz (8.33 ms)
Position, speed and temperature share one timestamp array in every file: True
Speed unit: {'cm/s'} ; position conversion to metres: {0.01}
x range over sessions -74.6 to 67.7, y range -114.5 to 67.3 (in cm)
Speed percentiles 1, 50, 99 (cm/s), median over sessions: [1.7, 14.4, 118.7]
Maze fitted per unit session: stem length L and circle radius R (cm)
L R frac_linearized
MS12-170716-172307 99.0 49.5 1.0
MS12-170719-095305 99.4 49.7 1.0
MS13-171128-113924 98.9 49.5 1.0
MS13-171201-130527 97.9 48.9 1.0
MS21-180621-141817 100.7 50.4 1.0
MS21-180625-153927 99.6 49.8 1.0
MS21-180628-155921 100.8 50.4 1.0
MS21-180629-110332 99.8 49.9 1.0
MS22-180621-131028 101.3 50.6 1.0
MS22-180628-120341 98.9 49.5 1.0
MS22-180629-110319 97.5 48.7 1.0
MS22-180711-112912 97.5 48.8 1.0
MS22-180720-110055 98.4 49.2 1.0
Diameter 2R: 97.5 to 101.3 cm; Methods: 110 cm diameter
Position is stored in centimetres (conversion 0.01 to metres) and speed in cm/s. The maze circle fitted to the arm positions has a radius of about 50 cm (diameter about 100 cm against the 110 cm of the Methods, which may be the outer edge of the track). Position is sampled at 120 Hz in the 2017-2018 MS12 to MS22 sessions and at 40 Hz in the MS10 sessions.
In [7]:
tshape = pd.DataFrame({k: {"T_min": r["temp_min"], "T_max": r["temp_max"], "unique_values": r["temp_n_unique"],
"frac_at_max": r["temp_frac_at_max"], "sd_2nd_diff_1Hz": r["temp_d2_sd"]}
for k, r in beh.items()}).T
tshape["code_has_temp_channel"] = [not (k.startswith("MS10") or k in C.CODE_CONTROLS) for k in tshape.index]
print(tshape.round(4).to_string())
grp = tshape.groupby("code_has_temp_channel")
print("\nSD of the second difference of temperature at 1 Hz (C), by whether the code's recording list has a temperature channel:")
print(grp.sd_2nd_diff_1Hz.describe().round(4).to_string())
fig, ax = plt.subplots(1, 3, figsize=(15, 3))
for a, k in zip(ax, ["MS10-170319-180352", "MS21-180719-155941", "MS21-180628-155921"]):
a.plot(beh[k]["temp_1hz_t"], beh[k]["temp_1hz"], lw=0.8)
a.set_title(k); a.set_xlabel("time (s)"); a.set_ylabel("MS temperature (C)")
plt.tight_layout(); plt.show()
print("Plotted ranges (C):", {k: (round(beh[k]["temp_1hz"].min(), 2), round(beh[k]["temp_1hz"].max(), 2))
for k in ["MS10-170319-180352", "MS21-180719-155941", "MS21-180628-155921"]})
T_min T_max unique_values frac_at_max sd_2nd_diff_1Hz code_has_temp_channel
MS10-170307-154746 37.0000 37.0000 1.0 1.0000 0.0000 False
MS10-170311-180956 17.0000 37.0000 58202.0 0.2633 0.0206 False
MS10-170313-122631 16.9777 37.0000 92783.0 0.1552 0.0246 False
MS10-170314-163038 17.0000 37.0000 50064.0 0.2570 0.0222 False
MS10-170317-153237 17.0000 37.0000 60000.0 0.3478 0.0191 False
MS10-170319-180352 17.0000 37.0000 80358.0 0.1933 0.0256 False
MS10-170320-212321 17.0000 37.0000 45001.0 0.3422 0.0136 False
MS10-170321-164406 17.0000 37.0000 165107.0 0.2908 0.0207 False
MS10-170323-173246 17.0000 37.0000 149526.0 0.2983 0.0170 False
MS12-170713-093400 19.7332 37.1922 245888.0 0.0000 0.1135 True
MS12-170714-122034 18.4883 50.1643 312653.0 0.0000 0.2549 True
MS12-170716-172307 18.5784 37.4376 276429.0 0.0000 0.0518 True
MS12-170717-111614 18.7501 36.9551 267761.0 0.0000 0.0616 True
MS12-170718-101720 19.5652 36.5724 297659.0 0.0000 0.0551 True
MS12-170719-095305 20.0605 36.9175 387733.0 0.0000 0.0140 True
MS12-170720-093528 20.5322 36.8443 292416.0 0.0000 0.0081 True
MS12-170723-141121 23.0687 36.8235 288316.0 0.0000 0.0398 True
MS13-171127-133330 17.8655 35.9522 270305.0 0.0000 0.0191 True
MS13-171128-113924 18.5004 37.8232 312838.0 0.0000 0.0221 True
MS13-171129-105507 19.0983 37.5706 266332.0 0.0000 0.1074 True
MS13-171130-121758 18.5401 37.5508 317625.0 0.0000 0.1115 True
MS13-171201-130527 20.2721 37.6232 203305.0 0.0000 0.0172 True
MS21-180621-141817 18.2381 40.3965 267810.0 0.0000 0.1149 True
MS21-180622-174411 20.6450 40.8108 220990.0 0.0000 0.0088 True
MS21-180623-170717 18.1842 45.4813 216006.0 0.0000 0.0346 True
MS21-180625-153927 16.7320 37.3512 215519.0 0.0000 0.0100 True
MS21-180627-143449 17.9093 37.3109 237891.0 0.0000 0.0727 True
MS21-180628-155921 20.2315 37.3577 210408.0 0.0000 0.0107 True
MS21-180629-110332 18.9900 37.2132 208505.0 0.0000 0.0105 True
MS21-180711-112927 18.0670 37.5082 250750.0 0.0000 0.0075 True
MS21-180712-103200 25.0488 36.9479 87788.0 0.0000 0.0505 True
MS21-180719-155941 17.0000 37.0000 180010.0 0.3391 0.0137 False
MS21-180720-110037 20.7409 45.7992 287853.0 0.0000 0.0312 True
MS21-180724-144402 20.1658 37.7594 263654.0 0.0000 0.0073 True
MS21-180730-120605 22.4984 38.4624 211180.0 0.0000 0.0097 True
MS21-180802-120147 21.8019 38.2259 236113.0 0.0000 0.0050 True
MS21-180803-142301 19.9210 37.7573 187675.0 0.0000 0.0102 True
MS22-180621-131028 22.3802 45.2808 356812.0 0.0000 0.1846 True
MS22-180622-115734 21.5986 49.9391 189203.0 0.0000 0.0420 True
MS22-180623-170141 21.2363 41.9523 195302.0 0.0000 0.0224 True
MS22-180625-162932 21.8525 37.7517 204527.0 0.0000 0.0095 True
MS22-180626-110916 23.3408 37.7988 187607.0 0.0000 0.0100 True
MS22-180627-114856 23.4917 37.3170 259029.0 0.0000 0.0081 True
MS22-180628-120341 22.6824 37.4339 255110.0 0.0000 0.0762 True
MS22-180629-110319 23.1095 37.3333 198664.0 0.0000 0.0090 True
MS22-180711-112912 23.3535 37.2892 237184.0 0.0000 0.0067 True
MS22-180712-102504 22.5606 37.9018 282412.0 0.0000 0.0071 True
MS22-180719-122813 17.0000 37.0000 126881.0 0.3333 0.0169 False
MS22-180720-110055 23.6435 37.4797 226322.0 0.0000 0.0085 True
MS22-180724-144421 24.0259 48.7159 207028.0 0.0000 0.0333 True
SD of the second difference of temperature at 1 Hz (C), by whether the code's recording list has a temperature channel:
count mean std min 25% 50% 75% max
code_has_temp_channel
False 11.0 0.0176 0.0071 0.000 0.0153 0.0191 0.0215 0.0256
True 39.0 0.0430 0.0538 0.005 0.0092 0.0191 0.0535 0.2549
Plotted ranges (C): {'MS10-170319-180352': (np.float64(17.0), np.float64(37.0)), 'MS21-180719-155941': (np.float64(17.0), np.float64(37.0)), 'MS21-180628-155921': (np.float64(20.23), np.float64(37.36))}
The code's recording list gives no temperature channel for the nine MS10 sessions and for the two
2018-07-19 control sessions (ch_temp = 0). In the dandiset these eleven sessions still have a
temperature trace, and it is not a recording. MS10-170307 is a constant 37.0 C. Each of the other
ten runs from 17.0 to 37.0 C (exactly in nine; MS10-170313 reaches 16.98 C) and sits at exactly 37.0 C for 16% to 35% of its samples
(frac_at_max), while no recorded trace takes its maximum value more than once (frac_at_max 0).
These are smooth curves of a cooling time course, consistent with a model; their roughness (SD of
the second difference) does not tell them apart from recordings. The trials of the two control
sessions are also labelled with cooling states. Recorded traces contain artefacts up to 50 C in some
sessions (T_max above).
In [8]:
series = pd.Series([p for r in beh.values() for p, _ in r["series"] + r["raw_info"]["raw_series"]])
print("Every series (NWB type ending in 'Series', and LFP) in the processed and raw files, and the number of sessions holding it:")
print(series.value_counts().to_string())
print("Files listed:", len(beh), "processed and", sum(1 for r in beh.values() if r["raw_info"].get("raw_series")), "raw")
print("Series whose path contains 'temp':", sorted({p for p in series if "temp" in p.lower()}))
ri = {k: r["raw_info"] for k, r in beh.items() if "lfp_column" in r["raw_info"]}
print("\nSessions whose LFP was read:", len(ri), "; LFP electrodes region equal to rows 0 to n-1 in order in", sum(v["region_is_identity"] for v in ri.values()))
print("LFP column read is the processed file's theta_reference row in every session:", all(v["lfp_column"] == beh[k]["theta_channel"] for k, v in ri.items()))
print("Raw file flags the same single theta_reference row as the processed file:", all(v["raw_theta_rows"] == [beh[k]["theta_channel"]] for k, v in ri.items()))
Every series (NWB type ending in 'Series', and LFP) in the processed and raw files, and the number of sessions holding it: processing/behavior/Acceleration 50 processing/behavior/SubjectPosition/SpatialSeries 50 processing/behavior/SubjectSpeed 50 processing/behavior/Temperature 50 acquisition/ElectricalSeries 50 processing/ecephys/LFP 50 processing/ecephys/LFP/LFP 50 Files listed: 50 processed and 50 raw Series whose path contains 'temp': ['processing/behavior/Temperature'] Sessions whose LFP was read: 49 ; LFP electrodes region equal to rows 0 to n-1 in order in 49 LFP column read is the processed file's theta_reference row in every session: True Raw file flags the same single theta_reference row as the processed file: True
Each of the 100 files was listed. The processed files hold four series in processing/behavior
(Acceleration, SubjectPosition, SubjectSpeed and Temperature) and the raw files hold the wide-band
ElectricalSeries and the LFP; the only temperature series is the MS Temperature, so the files
hold no hippocampal temperature. In the 49 sessions whose LFP was read, the LFP's electrodes
region lists the electrode rows in order, so the LFP column read is the row flagged
theta_reference, and the raw file flags the same row as the processed file.
In [9]:
def state_theta(tr):
m = tr.groupby("state").theta_freq.mean()
return {s: round(m.get(s, np.nan), 3) for s in ["Pre-Cooling", "Cooling on", "Cooling off", "Post-Cooling"]}
print("Mean trial theta frequency (Hz) by labelled state, in sessions the code lists as controls or that are not in its list:")
for k in sorted(C.CODE_CONTROLS) + ["MS10-170307-154746"]:
tr = beh[k]["trials"]
th = state_theta(tr)
pre = tr.theta_freq[tr.state == "Pre-Cooling"]; cool = tr.theta_freq[tr.state == "Cooling on"]
print(f" {k}: {th}; cooling minus pre {cool.mean() - pre.mean():+.3f} Hz, rank-sum P = {stats.ranksums(cool, pre).pvalue:.2g}")
ms10 = [k for k in info.index if k.startswith("MS10") and info.loc[k, "in_code_list"]]
d = [beh[k]["trials"].groupby("state").theta_freq.mean() for k in ms10]
print("\nMS10 sessions of the code's list: cooling minus pre theta (Hz):", [round(x["Cooling on"] - x["Pre-Cooling"], 2) for x in d])
others = [k for k in info.index if info.loc[k, "in_code_list"] and beh[k]["has_state"] and not k.startswith("MS10")]
d2 = [beh[k]["trials"].groupby("state").theta_freq.mean() for k in others]
print("Other sessions of the code's list: cooling minus pre theta (Hz), median", round(np.median([x["Cooling on"] - x["Pre-Cooling"] for x in d2]), 2),
", range", round(min(x["Cooling on"] - x["Pre-Cooling"] for x in d2), 2), "to", round(max(x["Cooling on"] - x["Pre-Cooling"] for x in d2), 2))
Mean trial theta frequency (Hz) by labelled state, in sessions the code lists as controls or that are not in its list:
MS21-180719-155941: {'Pre-Cooling': np.float64(7.767), 'Cooling on': np.float64(8.206), 'Cooling off': np.float64(8.309), 'Post-Cooling': np.float64(8.42)}; cooling minus pre +0.439 Hz, rank-sum P = 0.00055
MS22-180719-122813: {'Pre-Cooling': np.float64(8.731), 'Cooling on': np.float64(8.528), 'Cooling off': np.float64(8.976), 'Post-Cooling': np.float64(9.094)}; cooling minus pre -0.203 Hz, rank-sum P = 0.00013
MS10-170307-154746: {'Pre-Cooling': np.float64(8.271), 'Cooling on': np.float64(7.731), 'Cooling off': nan, 'Post-Cooling': np.float64(7.871)}; cooling minus pre -0.541 Hz, rank-sum P = 0.0002
MS10 sessions of the code's list: cooling minus pre theta (Hz): [np.float64(-1.06), np.float64(-0.76), np.float64(-1.08), np.float64(-0.86), np.float64(-0.7), np.float64(-0.37), np.float64(-0.72), np.float64(-0.65)]
Other sessions of the code's list: cooling minus pre theta (Hz), median -0.82 , range -1.23 to 0.5
In [10]:
print("Clock check: cooling minus pre theta frequency (Hz) with the theta series shifted in time against the trials.")
sh = pd.DataFrame({k: r["shift_contrast"] for k, r in beh.items() if r.get("shift_contrast")}).T
sel = [k for k in sh.index if info.loc[k, "in_code_list"]]
print("Median over the code's sessions, by shift (s):")
print(sh.loc[sel].median().round(3).to_string())
best = sh.loc[sel].idxmin(axis=1)
print("\nShift giving the most negative contrast, count of sessions:", best.value_counts().sort_index().to_dict())
far = [k for k in sel if best[k] not in (-60, -30, 0, 30, 60)]
print("Sessions where it is more than 60 s from zero:", {k: int(best[k]) for k in far})
print("Controls:", sh.loc[sorted(C.CODE_CONTROLS)].round(3).to_string())
Clock check: cooling minus pre theta frequency (Hz) with the theta series shifted in time against the trials.
Median over the code's sessions, by shift (s):
-600 0.484
-300 0.132
-120 -0.439
-60 -0.601
-30 -0.537
0 -0.809
30 -0.641
60 -0.611
120 -0.520
300 -0.249
600 NaN
Shift giving the most negative contrast, count of sessions: {-120: 2, -60: 2, 0: 29, 30: 6, 60: 4, 120: 2}
Sessions where it is more than 60 s from zero: {'MS10-170323-173246': 120, 'MS12-170718-101720': -120, 'MS12-170723-141121': -120, 'MS22-180724-144421': 120}
Controls: -600 -300 -120 -60 -30 0 30 60 120 300 600
MS21-180719-155941 0.081 0.119 -0.023 0.230 0.002 0.439 -0.158 0.120 -0.214 -0.426 -0.632
MS22-180719-122813 0.143 -0.081 -0.207 -0.148 -0.256 -0.203 -0.277 -0.195 -0.172 -0.226 NaN
Over the code's sessions the theta slowing is strongest with no shift, and in each session the best shift is within 120 s of zero (29 at zero); the cooling lasts about ten minutes, so the contrast changes slowly with the shift. The LFP and the trial table share a clock as far as this check can tell. In the eight MS10 sessions theta slows by 0.37 to 1.08 Hz in the "Cooling on" trials, so they were cooled, as the code's hand-noted onsets say, although their temperature is not a recording. In the two control sessions the "Cooling on" trials differ from the pre trials by -0.20 and +0.44 Hz, the weakest contrasts but one (MS12-170723, +0.50 Hz, has a recorded temperature drop) against a median of -0.82 Hz in the other sessions. This is consistent with the code's description of them as sessions without cooling; it does not prove it.
In [11]:
uc = pd.concat([u["units"].assign(session=k) for k, u in units.items()], ignore_index=True)
print("Good and good2 units in unit sessions:", len(uc), "; by session:", uc.groupby("session").size().to_dict())
print("\nClock check of spikes against position: spatial information (bits/spike) of every unit's map, true and with spikes shifted by 60 s")
print(uc.groupby("session")[["info_all", "info_all_shift60"]].median().round(3).to_string())
print("\nClock check of spikes against LFP: theta phase locking (mean resultant length) during runs, true and with spikes shifted by 0.5 s")
print(uc.groupby("session")[["theta_mrl", "theta_mrl_shift"]].median().round(3).to_string())
Good and good2 units in unit sessions: 1037 ; by session: {'MS12-170716-172307': 36, 'MS12-170719-095305': 48, 'MS13-171128-113924': 73, 'MS13-171201-130527': 67, 'MS21-180621-141817': 23, 'MS21-180625-153927': 83, 'MS21-180628-155921': 109, 'MS21-180629-110332': 172, 'MS22-180621-131028': 27, 'MS22-180628-120341': 79, 'MS22-180629-110319': 172, 'MS22-180711-112912': 68, 'MS22-180720-110055': 80}
Clock check of spikes against position: spatial information (bits/spike) of every unit's map, true and with spikes shifted by 60 s
info_all info_all_shift60
session
MS12-170716-172307 0.895 0.010
MS12-170719-095305 0.888 0.009
MS13-171128-113924 1.503 0.015
MS13-171201-130527 1.491 0.009
MS21-180621-141817 0.232 0.282
MS21-180625-153927 0.771 0.011
MS21-180628-155921 1.035 0.016
MS21-180629-110332 1.090 0.018
MS22-180621-131028 0.117 0.112
MS22-180628-120341 0.624 0.012
MS22-180629-110319 1.109 0.013
MS22-180711-112912 0.623 0.009
MS22-180720-110055 0.590 0.004
Clock check of spikes against LFP: theta phase locking (mean resultant length) during runs, true and with spikes shifted by 0.5 s
theta_mrl theta_mrl_shift
session
MS12-170716-172307 0.184 0.055
MS12-170719-095305 0.156 0.035
MS13-171128-113924 0.221 0.085
MS13-171201-130527 0.187 0.083
MS21-180621-141817 0.368 0.109
MS21-180625-153927 0.156 0.063
MS21-180628-155921 0.172 0.081
MS21-180629-110332 0.185 0.062
MS22-180621-131028 0.344 0.106
MS22-180628-120341 0.290 0.062
MS22-180629-110319 0.237 0.073
MS22-180711-112912 0.172 0.060
MS22-180720-110055 0.229 0.062
In eleven of the thirteen unit sessions the spatial information of the maps falls 40-fold or more when the spikes are shifted by 60 s, so spikes and position share a clock as far as this check can tell. In MS21-180621 and MS22-180621 it does not fall (0.23 against 0.28, and 0.12 against 0.11 bits/spike), and these sessions have 4 and 0 place cells. Their spikes are locked to theta (mean resultant length 0.37 and 0.34, against 0.11 when shifted), so they belong to the recording; the outputs cannot say whether the spike and position clocks differ there or the units were not spatially tuned (both are the rat's first session). The clock of spikes and position could not be established in these two sessions, and they are left out of the unit claims. Theta phase locking falls with a 0.5 s shift in every session, so spikes and LFP share a clock. Spikes and the wide-band data share it too, since mean waveforms at the spike times have the shape of action potentials (below).
Session sets used below (fixed in the claims file): BEH, the 45 sessions of the code's behavioural list that are here and have cooling states; TEMP, the 37 BEH sessions with a recorded temperature (BEH without MS10); UNIT, the 11 BEH sessions with curated units whose spikes carry spatial information (the 13 with curated units less the two 180621 sessions); FIG5, the 9 UNIT sessions of the code's Figure 5 list.
In [12]:
BEH = [k for k in info.index if info.loc[k, "in_code_list"] and info.loc[k, "has_state"]]
TEMP = [k for k in BEH if not k.startswith("MS10")]
NO_CLOCK = ["MS21-180621-141817", "MS22-180621-131028"]
UNIT = [k for k in BEH if k in units and k not in NO_CLOCK]
FIG5 = [k for k in UNIT if S.loc[k, "in_fig5_list"]]
for name, ks in [("BEH", BEH), ("TEMP", TEMP), ("UNIT", UNIT), ("FIG5", FIG5)]:
print(f"{name}: {len(ks)} sessions;", pd.Series([k.split('-')[0] for k in ks]).value_counts().sort_index().to_dict())
TR = {k: beh[k]["trials"] for k in BEH}
BEH: 45 sessions; {'MS10': 8, 'MS12': 8, 'MS13': 5, 'MS21': 12, 'MS22': 12}
TEMP: 37 sessions; {'MS12': 8, 'MS13': 5, 'MS21': 12, 'MS22': 12}
UNIT: 11 sessions; {'MS12': 2, 'MS13': 2, 'MS21': 3, 'MS22': 4}
FIG5: 9 sessions; {'MS12': 2, 'MS13': 1, 'MS21': 3, 'MS22': 3}
In [13]:
print("Putative cell types (made here; the files have none):", uc.cell_type.value_counts().to_dict())
print("By session:"); print(pd.crosstab(uc.session, uc.cell_type).to_string())
fig, ax = plt.subplots(1, 2, figsize=(11, 3))
ax[0].hist(uc.ttp.dropna(), bins=np.arange(0, 1.3, 0.025)); ax[0].axvline(C.NARROW_TTP_MS, color="k", ls="--")
ax[0].set_xlabel("trough-to-peak (ms)"); ax[0].set_ylabel("units")
ax[1].scatter(uc.ttp, uc.rate_trials, s=4); ax[1].set_yscale("log"); ax[1].axhline(C.WIDE_INT_RATE, color="k", ls="--")
ax[1].set_xlabel("trough-to-peak (ms)"); ax[1].set_ylabel("rate in trials (Hz)")
plt.tight_layout(); plt.show()
print("Trough-to-peak range (ms):", uc.ttp.min(), "to", uc.ttp.max(), "; rate range (Hz):", round(uc.rate_trials.min(), 3), "to", round(uc.rate_trials.max(), 1),
"; units with rate 0 not drawn on the log axis:", int((uc.rate_trials <= 0).sum()))
print("Units without 30 spikes in either waveform window (unclassified):", int((uc.cell_type == "unclassified").sum()))
print("Peak-to-peak amplitude of mean waveforms (uV), median by type:", uc.groupby("cell_type").amp_pre.median().round(1).to_dict())
for thr in [0.35, 0.425, 0.5]:
t_ = np.where(uc.ttp.isna(), "unclassified", np.where((uc.ttp <= thr) | (uc.rate_trials > C.WIDE_INT_RATE), "interneuron", "pyramidal"))
print(f"threshold {thr} ms:", pd.Series(t_).value_counts().to_dict())
Putative cell types (made here; the files have none): {'pyramidal': 623, 'unclassified': 271, 'interneuron': 143}
By session:
cell_type interneuron pyramidal unclassified
session
MS12-170716-172307 3 29 4
MS12-170719-095305 5 28 15
MS13-171128-113924 6 33 34
MS13-171201-130527 2 35 30
MS21-180621-141817 2 9 12
MS21-180625-153927 13 57 13
MS21-180628-155921 17 86 6
MS21-180629-110332 14 113 45
MS22-180621-131028 12 9 6
MS22-180628-120341 11 38 30
MS22-180629-110319 18 92 62
MS22-180711-112912 19 39 10
MS22-180720-110055 21 55 4
Trough-to-peak range (ms): 0.2 to 1.15 ; rate range (Hz): 0.0 to 74.3 ; units with rate 0 not drawn on the log axis: 56
Units without 30 spikes in either waveform window (unclassified): 271
Peak-to-peak amplitude of mean waveforms (uV), median by type: {'interneuron': 130.7, 'pyramidal': 136.5, 'unclassified': 153.4}
threshold 0.35 ms: {'pyramidal': 643, 'unclassified': 271, 'interneuron': 123}
threshold 0.425 ms: {'pyramidal': 623, 'unclassified': 271, 'interneuron': 143}
threshold 0.5 ms: {'pyramidal': 590, 'unclassified': 271, 'interneuron': 176}
Trough-to-peak times are bimodal around the 0.425 ms threshold, with a narrow group below it. Changing the threshold to 0.35 or 0.5 ms moves few units (printed above).
In [14]:
fields = pd.concat([units[k]["fields"].assign(session=k, subject=units[k]["subject"]) for k in UNIT if len(units[k]["fields"])], ignore_index=True)
print("Place fields of putative pyramidal cells:", len(fields), "in", fields.groupby(["session", "unit"]).ngroups, "cells;",
"by subject:", fields.groupby("subject").size().to_dict())
print("Pyramidal cells per UNIT session:", uc[(uc.cell_type == "pyramidal") & uc.session.isin(UNIT)].groupby("session").size().to_dict())
print("Place cells per UNIT session:", {k: units[k]["n_place_cells"] for k in UNIT})
print("Field length (cm): median", fields.stop_cm.sub(fields.start_cm).median(), ", range", fields.stop_cm.sub(fields.start_cm).min(), "to", fields.stop_cm.sub(fields.start_cm).max())
print("Fields by segment (0 stem, 1 right arm, 2 left arm):", fields.segment.value_counts().sort_index().to_dict())
print("Paper: n = 365 place fields (Figure 4B); 2948 pairs across 18 sessions (Figure 5C).")
Place fields of putative pyramidal cells: 643 in 374 cells; by subject: {'MS12': 95, 'MS13': 51, 'MS21': 249, 'MS22': 248}
Pyramidal cells per UNIT session: {'MS12-170716-172307': 29, 'MS12-170719-095305': 28, 'MS13-171128-113924': 33, 'MS13-171201-130527': 35, 'MS21-180625-153927': 57, 'MS21-180628-155921': 86, 'MS21-180629-110332': 113, 'MS22-180628-120341': 38, 'MS22-180629-110319': 92, 'MS22-180711-112912': 39, 'MS22-180720-110055': 55}
Place cells per UNIT session: {'MS12-170716-172307': 22, 'MS12-170719-095305': 23, 'MS13-171128-113924': 19, 'MS13-171201-130527': 17, 'MS21-180625-153927': 37, 'MS21-180628-155921': 59, 'MS21-180629-110332': 60, 'MS22-180628-120341': 12, 'MS22-180629-110319': 61, 'MS22-180711-112912': 27, 'MS22-180720-110055': 37}
Field length (cm): median 55.0 , range 20.0 to 165.0
Fields by segment (0 stem, 1 right arm, 2 left arm): {0: 154, 1: 252, 2: 237}
Paper: n = 365 place fields (Figure 4B); 2948 pairs across 18 sessions (Figure 5C).
In [15]:
print("Recordings deposited twice: identical trial start times, temperature traces or trial counts across sessions")
sig_tr = {k: tuple(np.round(r["trials"].start.values[:20], 3)) for k, r in beh.items() if r.get("trials") is not None}
sig_T = {k: (round(float(np.sum(r["temp_1hz"])), 3), len(r["temp_1hz"])) for k, r in beh.items()}
dup_tr = [(a, b) for i, a in enumerate(sig_tr) for b in list(sig_tr)[i + 1:] if sig_tr[a] == sig_tr[b]]
dup_T = [(a, b) for i, a in enumerate(sig_T) for b in list(sig_T)[i + 1:] if sig_T[a] == sig_T[b]]
print("pairs with identical first 20 trial starts:", dup_tr)
print("pairs with identical temperature traces:", dup_T)
Recordings deposited twice: identical trial start times, temperature traces or trial counts across sessions pairs with identical first 20 trial starts: [] pairs with identical temperature traces: []
No recording appears twice.
Helpers for the claims¶
Wilcoxon signed-rank and rank-sum tests are two-sided (scipy, exact where scipy uses it). Each test prints its n.
In [16]:
def wsr(a, b=None):
a = np.asarray(a, float)
d = a if b is None else a - np.asarray(b, float)
d = d[np.isfinite(d)]
if len(d) < 2 or np.all(d == 0):
return np.nan, len(d)
return float(stats.wilcoxon(d).pvalue), len(d)
def rks(a, b):
a, b = np.asarray(a, float), np.asarray(b, float)
a, b = a[np.isfinite(a)], b[np.isfinite(b)]
return float(stats.ranksums(a, b).pvalue), len(a), len(b)
def within(x, value, tol):
return bool(np.isfinite(x) and abs(x - value) <= tol + 1e-12)
def subj(k):
return k.split("-")[0]
session-structure¶
A session holds about 40 control trials, 20 to 60 cooling trials and 20 to 80 recovery trials, 80 to 200 trials in all with a mean of 150.
In [17]:
counts = pd.DataFrame({k: {"trials": len(t), "pre": (t.state == "Pre-Cooling").sum(), "cooling": (t.state == "Cooling on").sum(),
"cooling_off": (t.state == "Cooling off").sum(), "post": (t.state == "Post-Cooling").sum(),
"cooling_blocks": int(((t.state == "Cooling on") & (t.state.shift() != "Cooling on")).sum())}
for k, t in TR.items()}).T
print(counts.describe().round(1).to_string())
print("\nMean trials per session:", round(counts.trials.mean(), 1), "(expected 150 +- 15); range", counts.trials.min(), "to", counts.trials.max(), "(allowed 72 to 220)")
print("Median pre", counts.pre.median(), "(36 to 44); median cooling", counts.cooling.median(), "(18 to 66); median post", counts.post.median(), "(18 to 88)")
print("Sessions with two cooling blocks:", counts.index[counts.cooling_blocks > 1].tolist())
print("\nBy subject (mean):"); print(counts.groupby(lambda k: subj(k)).mean().round(1).to_string())
ok_mean = within(counts.trials.mean(), 150, 15)
ok_parts = [counts.trials.between(72, 220).all(), 36 <= counts.pre.median() <= 44, 18 <= counts.cooling.median() <= 66, 18 <= counts.post.median() <= 88]
print("\nmean within tolerance:", ok_mean, "; range, pre, cooling, post checks:", ok_parts)
def runs(states):
out = []
for st in states:
if out and out[-1][0] == st:
out[-1][1] += 1
else:
out.append([st, 1])
return " ".join(f"{ {'Pre-Cooling': 'pre', 'Cooling on': 'on', 'Cooling off': 'off', 'Post-Cooling': 'post'}[a]}x{n}" for a, n in out)
print("\nState sequences of the sessions with two cooling blocks (runs of trials), and the trials kept by the onset-aligned series:")
for k in counts.index[counts.cooling_blocks > 1]:
print(f" {k}: {runs(TR[k].state)}; kept {C.aligned_series(TR[k], columns=('error',))['n_kept']} of {len(TR[k])}")
trials pre cooling cooling_off post cooling_blocks
count 45.0 45.0 45.0 45.0 45.0 45.0
mean 143.8 39.3 44.1 12.9 47.4 1.3
std 25.3 7.3 15.0 7.8 18.1 0.4
min 80.0 27.0 20.0 0.0 15.0 1.0
25% 128.0 35.0 31.0 9.0 35.0 1.0
50% 141.0 39.0 47.0 14.0 44.0 1.0
75% 159.0 42.0 54.0 17.0 59.0 2.0
max 204.0 62.0 92.0 34.0 93.0 2.0
Mean trials per session: 143.8 (expected 150 +- 15); range 80 to 204 (allowed 72 to 220)
Median pre 39.0 (36 to 44); median cooling 47.0 (18 to 66); median post 44.0 (18 to 88)
Sessions with two cooling blocks: ['MS10-170313-122631', 'MS10-170319-180352', 'MS12-170713-093400', 'MS12-170714-122034', 'MS12-170716-172307', 'MS12-170717-111614', 'MS12-170718-101720', 'MS12-170719-095305', 'MS13-171127-133330', 'MS13-171128-113924', 'MS13-171130-121758', 'MS22-180621-131028']
By subject (mean):
trials pre cooling cooling_off post cooling_blocks
MS10 160.4 48.5 56.6 1.2 54.0 1.2
MS12 139.4 36.9 49.4 15.9 37.2 1.8
MS13 149.4 44.6 36.4 19.8 48.6 1.6
MS21 146.8 37.6 37.5 15.7 56.1 1.0
MS22 130.2 34.4 42.2 13.1 40.6 1.1
mean within tolerance: True ; range, pre, cooling, post checks: [np.True_, np.True_, np.True_, np.True_]
State sequences of the sessions with two cooling blocks (runs of trials), and the trials kept by the onset-aligned series:
MS10-170313-122631: prex39 onx28 postx52 onx19 postx5; kept 119 of 143
MS10-170319-180352: prex49 onx49 postx63 onx43; kept 161 of 204
MS12-170713-093400: prex32 onx27 offx13 postx19 onx24; kept 91 of 115
MS12-170714-122034: prex39 onx36 offx12 postx18 onx21 offx7 postx18; kept 105 of 151
MS12-170716-172307: prex41 onx25 offx17 postx36 onx27; kept 119 of 146
MS12-170717-111614: prex39 onx31 offx15 postx26 onx23; kept 111 of 134
MS12-170718-101720: prex38 onx35 offx10 postx17 onx25; kept 100 of 125
MS12-170719-095305: prex35 onx20 offx15 postx22 onx27 offx11 postx19; kept 92 of 149
MS13-171127-133330: prex42 onx10 offx17 postx15 onx25 offx17 postx34; kept 84 of 160
MS13-171128-113924: prex62 onx25 offx17 postx50 onx5; kept 154 of 159
MS13-171130-121758: prex33 onx31 offx15 postx39 onx21 offx15 postx20; kept 118 of 174
MS22-180621-131028: prex31 onx7 offx9 postx38 onx18 offx8 postx4; kept 85 of 115
Result: reproduced. The 45 BEH sessions hold a mean of 143.8 trials (80 to 204), with medians of 39 pre-cooling, 47 cooling and 44 post-cooling trials, all within the criterion. 12 sessions have two cooling blocks, which the paper does not mention. In each of them the second block follows post-cooling trials, and the onset-aligned series used below stops before it.
cooling-time-course¶
MS cooling peaked about 60 s after its onset and lasted about 10 to 12 minutes.
In [18]:
rows = []
for k in TEMP:
r = beh[k]; t1, T1 = r["temp_1hz_t"], r["temp_1hz"]
tr = TR[k]
base = float(np.nanmedian(tr.temperature[tr.state == "Pre-Cooling"]))
t_first_cool = float(tr.start[tr.state == "Cooling on"].iloc[0])
after = t1 >= float(tr.start.iloc[0])
cand = np.flatnonzero(after & (T1 < base - 1) & (t1 > t_first_cool - 600))
if not len(cand):
continue
i0 = cand[0]
win = (t1 >= t1[i0]) & (t1 <= t1[i0] + 300)
imin = np.flatnonzero(win)[np.argmin(T1[win])]
below = np.flatnonzero((t1 > t1[i0]) & (T1 < 32))
dur32 = np.nan
if len(below):
j = below[0]
back = np.flatnonzero((t1 > t1[j]) & (T1 >= 32))
dur32 = (t1[back[0]] - t1[j]) / 60 if len(back) else np.nan
rec = np.flatnonzero((t1 > t1[imin]) & (T1 >= base - 1))
dur_base = (t1[rec[0]] - t1[i0]) / 60 if len(rec) else np.nan
rows.append({"session": k, "onset_t": t1[i0], "baseline_C": base, "onset_to_first_cooling_trial_s": t_first_cool - t1[i0],
"time_to_min_s": t1[imin] - t1[i0], "min_C": T1[imin], "below32_min": dur32, "to_baseline_min": dur_base})
ct = pd.DataFrame(rows).set_index("session")
print(ct.drop(columns="onset_t").round(1).to_string())
print("\nTEMP sessions measured:", len(ct), "of", len(TEMP))
print("Median time to minimum:", round(ct.time_to_min_s.median(), 1), "s (expected 60 +- 6); IQR", ct.time_to_min_s.quantile([.25, .75]).round(1).tolist())
print("Median duration below 32 C:", round(ct.below32_min.median(), 2), "min; until within 1 C of baseline:", round(ct.to_baseline_min.median(), 2),
"min (expected 9 to 13.2)")
print("Median minimum temperature:", round(ct.min_C.median(), 1), "C; range", round(ct.min_C.min(), 1), "to", round(ct.min_C.max(), 1))
print("By subject (median):"); print(ct.groupby(lambda k: subj(k))[["time_to_min_s", "below32_min", "to_baseline_min", "min_C"]].median().round(1).to_string())
for thr in [0.5, 1, 2]:
vals = []
for k in TEMP:
r = beh[k]; t1, T1 = r["temp_1hz_t"], r["temp_1hz"]; tr = TR[k]
base = float(np.nanmedian(tr.temperature[tr.state == "Pre-Cooling"])); tc = float(tr.start[tr.state == "Cooling on"].iloc[0])
cand = np.flatnonzero((t1 >= float(tr.start.iloc[0])) & (T1 < base - thr) & (t1 > tc - 600))
if len(cand):
w = (t1 >= t1[cand[0]]) & (t1 <= t1[cand[0]] + 300); vals.append(t1[np.flatnonzero(w)[np.argmin(T1[w])]] - t1[cand[0]])
print(f"onset threshold {thr} C below baseline: median time to minimum {np.median(vals):.1f} s (n = {len(vals)})")
# onset-aligned mean trace and the time to 90% of the drop
grid = np.arange(-60, 1201)
trace = np.vstack([np.interp(ct.onset_t[k] + grid, beh[k]["temp_1hz_t"], beh[k]["temp_1hz"], left=np.nan, right=np.nan) for k in ct.index])
mtrace = np.nanmean(trace, axis=0)
t90 = []
for k, row in zip(ct.index, trace):
after = grid >= 0
drop = ct.baseline_C[k] - ct.min_C[k]
hit = np.flatnonzero(after & (row <= ct.baseline_C[k] - 0.9 * drop))
t90.append(grid[hit[0]] if len(hit) else np.nan)
ct["t90_s"] = t90
print("\nOnset-aligned mean MS temperature (C) at", {s_: round(float(mtrace[grid == s_][0]), 2) for s_ in [-60, 0, 30, 60, 90, 120, 180, 300, 600, 900, 1200]})
print("Minimum of the mean trace at", int(grid[np.nanargmin(mtrace)]), "s after onset; it is within 0.5 C of that minimum from",
int(grid[np.flatnonzero(mtrace <= np.nanmin(mtrace) + 0.5)[0]]), "to", int(grid[np.flatnonzero(mtrace <= np.nanmin(mtrace) + 0.5)[-1]]), "s")
print("Time to 90% of each session's drop: median", np.nanmedian(t90), "s, IQR", np.nanpercentile(t90, [25, 75]).tolist())
art = ct[(ct.baseline_C > 38) | (ct.onset_to_first_cooling_trial_s > 300)]
print("Sessions with a baseline above 38 C or an onset more than 300 s before the first cooling trial:", art.index.tolist())
cl = ct.drop(art.index)
print(f"variant, without them (n = {len(cl)}): median time to minimum {cl.time_to_min_s.median():.1f} s, to 90% {cl.t90_s.median():.1f} s, "
f"until within 1 C of baseline {cl.to_baseline_min.median():.2f} min")
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(grid, mtrace); ax.axvline(0, color="k", lw=0.5); ax.set_xlabel("time from onset (s)"); ax.set_ylabel("mean MS temperature (C)")
plt.tight_layout(); plt.show()
print("Plotted range:", round(float(np.nanmin(mtrace)), 2), "to", round(float(np.nanmax(mtrace)), 2), "C")
baseline_C onset_to_first_cooling_trial_s time_to_min_s min_C below32_min to_baseline_min
session
MS12-170713-093400 36.9 18.3 121.0 19.8 9.0 12.2
MS12-170714-122034 36.8 17.1 76.0 21.0 9.1 13.2
MS12-170716-172307 37.1 13.4 149.0 20.4 9.1 14.8
MS12-170717-111614 36.7 17.4 193.0 21.5 10.6 14.5
MS12-170718-101720 36.4 25.3 122.0 20.2 10.7 16.8
MS12-170719-095305 36.6 49.8 122.0 20.1 10.6 16.9
MS12-170720-093528 36.6 19.6 74.0 20.5 9.9 15.1
MS12-170723-141121 36.6 31.1 141.0 23.1 12.6 21.2
MS13-171127-133330 35.8 41.9 24.0 24.9 2.7 4.0
MS13-171128-113924 37.4 17.8 42.0 20.4 6.4 8.2
MS13-171129-105507 36.6 16.1 138.0 19.2 10.1 12.2
MS13-171130-121758 37.3 14.4 131.0 23.7 5.8 7.8
MS13-171201-130527 37.5 28.2 131.0 20.3 5.1 7.4
MS21-180621-141817 37.0 24.4 98.0 18.3 9.6 11.6
MS21-180622-174411 40.2 72.4 124.0 20.6 6.1 NaN
MS21-180623-170717 42.3 13.9 85.0 18.2 6.5 NaN
MS21-180625-153927 37.0 26.4 215.0 16.7 8.8 11.8
MS21-180628-155921 37.0 42.1 81.0 20.2 5.2 8.0
MS21-180629-110332 36.9 21.3 121.0 19.0 8.3 11.0
MS21-180711-112927 37.0 20.3 68.0 18.1 10.3 13.8
MS21-180720-110037 41.3 413.5 46.0 39.0 10.2 1.8
MS21-180724-144402 37.5 17.2 94.0 20.2 10.8 15.7
MS21-180730-120605 37.3 25.5 89.0 22.5 4.7 8.0
MS21-180802-120147 37.8 20.4 105.0 21.8 5.8 9.8
MS21-180803-142301 37.4 17.0 236.0 19.9 9.9 12.9
MS22-180621-131028 39.1 32.0 99.0 24.3 4.1 5.8
MS22-180622-115734 38.7 27.0 132.0 21.6 8.0 NaN
MS22-180623-170141 37.6 21.5 184.0 21.2 8.8 12.1
MS22-180625-162932 37.3 39.4 129.0 21.9 8.0 11.8
MS22-180626-110916 37.4 33.5 115.0 23.3 7.9 10.9
MS22-180627-114856 37.1 28.5 155.0 23.5 8.5 12.4
MS22-180628-120341 37.1 29.0 139.0 22.7 9.0 12.4
MS22-180629-110319 37.0 19.4 130.0 23.1 9.8 13.8
MS22-180711-112912 36.5 30.4 100.0 23.4 8.1 11.3
MS22-180712-102504 37.4 21.1 113.0 22.6 8.1 12.1
MS22-180720-110055 37.3 12.0 67.0 23.6 7.4 13.3
MS22-180724-144421 37.8 29.0 140.0 24.0 10.8 17.2
TEMP sessions measured: 37 of 37
Median time to minimum: 121.0 s (expected 60 +- 6); IQR [89.0, 138.0]
Median duration below 32 C: 8.77 min; until within 1 C of baseline: 12.16 min (expected 9 to 13.2)
Median minimum temperature: 21.2 C; range 16.7 to 39.0
By subject (median):
time_to_min_s below32_min to_baseline_min min_C
session
MS12 122.0 10.2 14.9 20.5
MS13 131.0 5.8 7.8 20.4
MS21 96.0 8.5 11.3 20.0
MS22 129.5 8.1 12.1 23.2
onset threshold 0.5 C below baseline: median time to minimum 122.0 s (n = 37)
onset threshold 1 C below baseline: median time to minimum 121.0 s (n = 37)
onset threshold 2 C below baseline: median time to minimum 119.0 s (n = 37)
Onset-aligned mean MS temperature (C) at {-60: 37.55, 0: 36.13, 30: 25.66, 60: 23.52, 90: 22.68, 120: 22.6, 180: 23.83, 300: 26.9, 600: 33.47, 900: 36.49, 1200: 35.13}
Minimum of the mean trace at 121 s after onset; it is within 0.5 C of that minimum from 70 to 152 s
Time to 90% of each session's drop: median 53.0 s, IQR [44.0, 63.0]
Sessions with a baseline above 38 C or an onset more than 300 s before the first cooling trial: ['MS21-180622-174411', 'MS21-180623-170717', 'MS21-180720-110037', 'MS22-180621-131028', 'MS22-180622-115734']
variant, without them (n = 32): median time to minimum 121.5 s, to 90% 53.5 s, until within 1 C of baseline 12.22 min
Plotted range: 22.6 to 37.66 C
Result: partially reproduced. The cooling lasted a median 12.16 minutes from onset until the temperature was again within 1 C of baseline (reading b, within 9 to 13.2 minutes; 8.77 minutes below 32 C under reading a). The minimum came a median 121 s after the onset (IQR 89 to 138 s), twice the 60 s of the Methods, and the onset threshold (0.5, 1 or 2 C) does not change this (119 to 122 s). The temperature falls fast and then stays near its minimum: 90% of each session's drop is reached a median 53 s after onset (IQR 44 to 63 s), and the onset-aligned mean trace has its minimum at 121 s and is within 0.5 C of it from 70 to 152 s. Leaving out the five sessions with an artefact baseline above 38 C or an early onset gives 121.5 s to the minimum and 12.22 minutes.
fig2g-theta-freq-temperature¶
Theta frequency decreased linearly with MS temperature (R = 0.81, P = 5.4e-48).
In [19]:
AL = {k: C.aligned_series(TR[k], columns=("error", "temperature", "speed", "theta_freq", "theta_power", "theta_freq_10", "theta_power_10")) for k in BEH}
def stack(keys, col, mask_ms10_temp=False):
rows = []
for k in keys:
x = AL[k][col].copy()
if mask_ms10_temp and col == "temperature" and k.startswith("MS10"):
x[:] = np.nan
rows.append(x)
return np.vstack(rows)
def mean_series(keys, col, mask=False):
return np.nanmean(stack(keys, col, mask), axis=0)
pos = np.arange(1, C.ALIGN_LENGTH + 1)
readings = {"(a) TEMP sessions": TEMP, "(b) BEH with MS10 curves": BEH}
res_2g = {}
for name, ks in readings.items():
r, p, n = C.pearson(mean_series(ks, "temperature"), mean_series(ks, "theta_freq"))
res_2g[name] = (r, p, n)
print(f"{name}: R = {r:.3f}, P = {p:.2g}, positions = {n}, sessions = {len(ks)}")
r, p, n = C.pearson(mean_series(TEMP, "temperature"), mean_series(BEH, "theta_freq"))
print(f"variant, temperature over TEMP and theta over BEH (as the code with MS10 temperature missing): R = {r:.3f}, P = {p:.2g}, n = {n}")
r, p, n = C.pearson(mean_series(TEMP, "temperature"), mean_series(TEMP, "theta_freq_10"))
print(f"variant, trial theta over frames above 10 cm/s, TEMP: R = {r:.3f}, P = {p:.2g}")
pooled = pd.concat([TR[k] for k in TEMP])
r, p, n = C.pearson(pooled.temperature, pooled.theta_freq)
sl = np.polyfit(pooled.temperature[pooled.theta_freq.notna()], pooled.theta_freq.dropna(), 1)[0]
print(f"variant, single trials pooled over TEMP: R = {r:.3f}, P = {p:.2g}, trials = {n}, slope {sl:.3f} Hz/C")
print("Expected R 0.81 +- 0.081 with P < 1e-10; within under:", [k for k, v in res_2g.items() if within(v[0], 0.81, 0.081) and v[1] < 1e-10])
mt, mf = mean_series(TEMP, "temperature"), mean_series(TEMP, "theta_freq")
fig, ax = plt.subplots(1, 2, figsize=(11, 3.5))
ax[0].plot(pos, mt, "b"); ax0 = ax[0].twinx(); ax0.plot(pos, mf, "m"); ax[0].set_xlabel("aligned trial position (65 = first cooling trial)")
ax[0].set_ylabel("temperature (C)"); ax0.set_ylabel("theta (Hz)")
ax[1].scatter(mt, mf, s=6); ax[1].set_xlabel("mean temperature (C)"); ax[1].set_ylabel("mean theta frequency (Hz)")
plt.tight_layout(); plt.show()
print("Plotted ranges: temperature", np.nanmin(mt).round(2), "to", np.nanmax(mt).round(2), "C; theta", np.nanmin(mf).round(2), "to", np.nanmax(mf).round(2), "Hz")
(a) TEMP sessions: R = 0.789, P = 2.8e-43, positions = 198, sessions = 37 (b) BEH with MS10 curves: R = 0.758, P = 3e-38, positions = 198, sessions = 45 variant, temperature over TEMP and theta over BEH (as the code with MS10 temperature missing): R = 0.758, P = 3.5e-38, n = 198 variant, trial theta over frames above 10 cm/s, TEMP: R = 0.783, P = 2.8e-42 variant, single trials pooled over TEMP: R = 0.672, P = 0, trials = 5187, slope 0.065 Hz/C Expected R 0.81 +- 0.081 with P < 1e-10; within under: ['(a) TEMP sessions', '(b) BEH with MS10 curves']
Plotted ranges: temperature 22.59 to 37.89 C; theta 7.67 to 9.57 Hz
Result: reproduced. With the code's onset-aligned series, theta frequency correlates with MS temperature at R = 0.789 (P = 2.8e-43, 198 positions) over the 37 sessions with a recorded temperature, and at R = 0.758 with the MS10 temperature curves included, both within 0.81 +- 0.081. Over single trials pooled the correlation is R = 0.672, with a slope of 0.065 Hz per C.
fig2h-theta-power-temperature¶
Theta power also decreased with MS temperature, less tightly than frequency (R = 0.54, P = 1.2e-16).
In [20]:
res_2h = {}
for name, ks in readings.items():
r, p, n = C.pearson(mean_series(ks, "temperature"), mean_series(ks, "theta_power"))
res_2h[name] = (r, p, n)
print(f"{name}: R = {r:.3f}, P = {p:.2g}, positions = {n}; frequency R under the same reading {res_2g[name][0]:.3f}")
r, p, n = C.pearson(mean_series(TEMP, "temperature"), mean_series(TEMP, "theta_power_10"))
print(f"variant, frames above 10 cm/s, TEMP: R = {r:.3f}, P = {p:.2g}")
pooled_p = pd.concat([TR[k].assign(tp=TR[k].theta_power - TR[k].theta_power.mean()) for k in TEMP])
r, p, n = C.pearson(pooled_p.temperature, pooled_p.tp)
print(f"variant, single trials pooled (power minus session mean), TEMP: R = {r:.3f}, P = {p:.2g}, n = {n}")
print("Within tolerance with P < 1e-5 and below the frequency R under:",
[k for k, v in res_2h.items() if within(v[0], 0.54, 0.054) and v[1] < 1e-5 and v[0] < res_2g[k][0]])
(a) TEMP sessions: R = -0.056, P = 0.44, positions = 198; frequency R under the same reading 0.789 (b) BEH with MS10 curves: R = 0.123, P = 0.085, positions = 198; frequency R under the same reading 0.758 variant, frames above 10 cm/s, TEMP: R = 0.135, P = 0.058 variant, single trials pooled (power minus session mean), TEMP: R = 0.107, P = 9.7e-15, n = 5187 Within tolerance with P < 1e-5 and below the frequency R under: []
Result: not reproduced. Theta power, computed as the code computes it, does not follow MS temperature: R = -0.056 (P = 0.44) over TEMP sessions and R = 0.123 (P = 0.085) with the MS10 curves, against R = 0.54. Over single trials pooled R = 0.107. The power measure is the code's (spectrogram magnitude of the band-passed LFP), on the electrode the files flag as theta reference.
fig2j-frequency-consistent-across-animals¶
The relationship of theta frequency with MS temperature was more consistent across animals than that of theta power.
In [21]:
rows = []
for a in sorted({subj(k) for k in TEMP}):
ks = [k for k in TEMP if subj(k) == a]
rf = C.pearson(mean_series(ks, "temperature"), mean_series(ks, "theta_freq"))
rp = C.pearson(mean_series(ks, "temperature"), mean_series(ks, "theta_power"))
rows.append({"animal": a, "sessions": len(ks), "R_freq": rf[0], "P_freq": rf[1], "R_power": rp[0], "P_power": rp[1]})
pa = pd.DataFrame(rows).set_index("animal")
print(pa.to_string(float_format=lambda x: f"{x:.3g}"))
print("\nAll frequency R positive with P < 0.05:", bool(((pa.R_freq > 0) & (pa.P_freq < 0.05)).all()))
print("Range of frequency R:", round(pa.R_freq.max() - pa.R_freq.min(), 3), "; range of power R:", round(pa.R_power.max() - pa.R_power.min(), 3))
print("MS10 (no recorded temperature) is not in this table.")
sessions R_freq P_freq R_power P_power animal MS12 8 0.685 4.81e-23 -0.346 9.3e-06 MS13 5 0.846 1.95e-43 0.122 0.131 MS21 12 0.954 2.04e-95 0.376 1.83e-07 MS22 12 0.919 1.92e-71 -0.0969 0.203 All frequency R positive with P < 0.05: True Range of frequency R: 0.269 ; range of power R: 0.722 MS10 (no recorded temperature) is not in this table.
Result: reproduced. In each of the four rats with a recorded temperature, theta frequency correlates positively with temperature (R = 0.685 to 0.954, every P below 1e-22), a range of 0.269, while the power correlations range over 0.722 (from -0.346 to 0.376). The outcome rests on power varying across rats, not on a consistent power effect (see fig2h-theta-power-temperature). MS10 has no recorded temperature and is not in the test.
fig2i-theta-pre-cooling-post¶
Theta frequency and power were lower during MS cooling than before and after it.
In [22]:
sm = pd.DataFrame({k: {f"{m}_{C.SHORT[s]}": TR[k].loc[TR[k].state == s, m].mean() for m in ["theta_freq", "theta_power"] for s in C.STATES} for k in BEH}).T
print(sm.groupby(lambda k: subj(k)).mean().round(3).to_string())
print("\nMean over sessions:", sm.mean().round(3).to_dict())
out = {}
for m in ["theta_freq", "theta_power"]:
for other in ["pre", "post"]:
p, n = wsr(sm[f"{m}_cooling"], sm[f"{m}_{other}"])
lower = int((sm[f"{m}_cooling"] < sm[f"{m}_{other}"]).sum())
out[(m, other)] = (p, lower > n - lower)
print(f"{m}: cooling against {other}, cooling lower in {lower} of {n} sessions, signed-rank P = {p:.2g}, median difference {np.median(sm[f'{m}_cooling'] - sm[f'{m}_{other}']):.3g}")
pn = sm[[c for c in sm.columns if "power" in c]].div(sm[[c for c in sm.columns if "power" in c]].mean(axis=1), axis=0)
print("\nPower relative to the session mean of the three states, mean over sessions:", pn.mean().round(3).to_dict())
theta_freq_pre theta_freq_cooling theta_freq_post theta_power_pre theta_power_cooling theta_power_post
MS10 7.811 7.038 7.897 27669.561 22897.731 27833.071
MS12 8.082 7.674 8.455 45186.609 47974.823 43815.427
MS13 9.007 8.118 8.871 47137.315 45803.194 49558.150
MS21 8.489 7.574 8.366 58297.688 54009.993 59295.749
MS22 8.806 8.038 8.827 53780.929 52551.309 50155.489
Mean over sessions: {'theta_freq_pre': 8.438, 'theta_freq_cooling': 7.681, 'theta_freq_post': 8.477, 'theta_power_pre': 48077.319, 'theta_power_cooling': 46105.156, 'theta_power_post': 47430.969}
theta_freq: cooling against pre, cooling lower in 44 of 45 sessions, signed-rank P = 1.4e-12, median difference -0.809
theta_freq: cooling against post, cooling lower in 45 of 45 sessions, signed-rank P = 5.7e-14, median difference -0.803
theta_power: cooling against pre, cooling lower in 26 of 45 sessions, signed-rank P = 0.048, median difference -992
theta_power: cooling against post, cooling lower in 19 of 45 sessions, signed-rank P = 0.94, median difference 770
Power relative to the session mean of the three states, mean over sessions: {'theta_power_pre': 1.019, 'theta_power_cooling': 0.972, 'theta_power_post': 1.009}
Result: partially reproduced. Theta frequency is lower during cooling than before (44 of 45 sessions, P = 1.4e-12) and after (45 of 45, P = 5.7e-14), by medians of 0.81 and 0.80 Hz. Theta power is lower than before at P = 0.048 (26 of 45 sessions) but not lower than after (19 of 45, P = 0.94); during cooling it is 0.972 of the session mean of the three states.
disc-theta-frequency-12-percent¶
During MS cooling LFP theta frequency decreased by 12%.
In [23]:
dec_a = 100 * (1 - sm.theta_freq_cooling / sm.theta_freq_pre)
dec_b = 100 * (1 - sm.theta_freq_cooling / ((sm.theta_freq_pre + sm.theta_freq_post) / 2))
pooled_all = pd.concat([TR[k] for k in BEH])
dec_c = 100 * (1 - pooled_all.theta_freq[pooled_all.state == "Cooling on"].mean() / pooled_all.theta_freq[pooled_all.state == "Pre-Cooling"].mean())
print(f"(a) mean over sessions of the decrease against pre: {dec_a.mean():.2f}% (median {dec_a.median():.2f}, range {dec_a.min():.2f} to {dec_a.max():.2f})")
print(f"(b) against the mean of pre and post: {dec_b.mean():.2f}%")
print(f"(c) pooled trials: {dec_c:.2f}%")
print("By subject, reading (a):", dec_a.groupby(lambda k: subj(k)).mean().round(2).to_dict())
mins = pd.Series({k: 100 * (1 - TR[k].theta_freq[TR[k].state == "Cooling on"].min() / sm.loc[k, "theta_freq_pre"]) for k in BEH})
print(f"variant, decrease at the slowest cooling trial: mean {mins.mean():.2f}%")
print("Within 12 +- 1.2 under:", [n for n, v in [("a", dec_a.mean()), ("b", dec_b.mean()), ("c", dec_c)] if within(v, 12, 1.2)])
(a) mean over sessions of the decrease against pre: 8.88% (median 9.56, range -7.19 to 13.90)
(b) against the mean of pre and post: 9.18%
(c) pooled trials: 9.23%
By subject, reading (a): {'MS10': 9.85, 'MS12': 4.76, 'MS13': 9.86, 'MS21': 10.74, 'MS22': 8.72}
variant, decrease at the slowest cooling trial: mean 16.93%
Within 12 +- 1.2 under: []
Result: partially reproduced. Theta frequency decreases during cooling under every reading, by 8.88% (a), 9.18% (b) and 9.23% (c), all outside 12 +- 1.2%. The slowest cooling trial of a session is a mean 16.93% below pre. The paper's own within-field values (8.4 and 7.7 Hz) are matched here (see fig4b-fp-flfp-decrease).
fig6d-errors-increase¶
MS cooling increased choice errors about 3-fold, within each rat (P < 0.05 for all animals) and across all animals (P = 1e-6).
In [24]:
def rates_table(off_as=None, keys=BEH):
return pd.DataFrame({k: C.error_rates(TR[k], off_as) for k in keys}).T
er = rates_table()
print(er[["pre", "cooling", "post"]].groupby(lambda k: subj(k)).mean().round(2).to_string())
print("\nMean error rate (%) over", len(er), "sessions: pre", round(er.pre.mean(), 2), ", cooling", round(er.cooling.mean(), 2), ", post", round(er.post.mean(), 2))
ratio_a = er.cooling.mean() / er.pre.mean()
ratio_b = er.cooling.mean() / ((er.pre + er.post) / 2).mean()
ok = er.pre > 0
ratio_c = (er.cooling[ok] / er.pre[ok]).median()
print(f"ratio (a) {ratio_a:.2f}; (b) {ratio_b:.2f}; (c) {ratio_c:.2f} over {int(ok.sum())} sessions with pre > 0 (expected 3.0 +- 0.3)")
p_all, n_all = wsr(er.cooling, er.pre)
print(f"cooling against pre over sessions: higher in {(er.cooling > er.pre).sum()} of {n_all}, signed-rank P = {p_all:.2g}")
rat = {}
for a in sorted({subj(k) for k in BEH}):
e = er.loc[[k for k in BEH if subj(k) == a]]
p, n = wsr(e.cooling, e.pre)
rat[a] = (p, n)
print(f" {a}: {n} sessions, pre {e.pre.mean():.1f}%, cooling {e.cooling.mean():.1f}%, ratio {e.cooling.mean() / e.pre.mean():.2f}, higher in {(e.cooling > e.pre).sum()}, P = {p:.3g}")
print("Rats with at least 6 sessions and P < 0.05:", [a for a, (p, n) in rat.items() if n >= 6 and p < 0.05], "of", [a for a, (p, n) in rat.items() if n >= 6])
for off in ["cooling", "post"]:
e2 = rates_table(off)
print(f"variant, 'Cooling off' trials counted as {off}: pre {e2.pre.mean():.2f}, cooling {e2.cooling.mean():.2f}, post {e2.post.mean():.2f}, ratio (a) {e2.cooling.mean() / e2.pre.mean():.2f}, P = {wsr(e2.cooling, e2.pre)[0]:.2g}")
e3 = rates_table(keys=TEMP)
print(f"variant, TEMP sessions only: ratio (a) {e3.cooling.mean() / e3.pre.mean():.2f}, P = {wsr(e3.cooling, e3.pre)[0]:.2g}")
e4 = pd.DataFrame({k: C.error_rates(beh[k]["trials"]) for k in BEH + sorted(C.CODE_CONTROLS) + ["MS10-170307-154746"]}).T
print(f"variant, with the two controls and MS10-170307: ratio (a) {e4.cooling.mean() / e4.pre.mean():.2f}, P = {wsr(e4.cooling, e4.pre)[0]:.2g}")
print("Controls alone:", e4.loc[sorted(C.CODE_CONTROLS), ["pre", "cooling", "post"]].round(1).to_dict("index"))
pre cooling post MS10 4.13 20.05 8.44 MS12 4.26 8.09 2.39 MS13 3.40 9.61 4.50 MS21 6.84 16.47 4.86 MS22 11.18 17.86 10.36 Mean error rate (%) over 45 sessions: pre 6.67 , cooling 15.23 , post 6.49 ratio (a) 2.28; (b) 2.31; (c) 2.53 over 42 sessions with pre > 0 (expected 3.0 +- 0.3) cooling against pre over sessions: higher in 40 of 45, signed-rank P = 3e-07 MS10: 8 sessions, pre 4.1%, cooling 20.1%, ratio 4.86, higher in 8, P = 0.00781 MS12: 8 sessions, pre 4.3%, cooling 8.1%, ratio 1.90, higher in 7, P = 0.0391 MS13: 5 sessions, pre 3.4%, cooling 9.6%, ratio 2.83, higher in 5, P = 0.0625 MS21: 12 sessions, pre 6.8%, cooling 16.5%, ratio 2.41, higher in 10, P = 0.0269 MS22: 12 sessions, pre 11.2%, cooling 17.9%, ratio 1.60, higher in 10, P = 0.064 Rats with at least 6 sessions and P < 0.05: ['MS10', 'MS12', 'MS21'] of ['MS10', 'MS12', 'MS21', 'MS22'] variant, 'Cooling off' trials counted as cooling: pre 6.67, cooling 13.23, post 6.49, ratio (a) 1.98, P = 1.6e-06 variant, 'Cooling off' trials counted as post: pre 6.67, cooling 15.23, post 6.37, ratio (a) 2.28, P = 3e-07 variant, TEMP sessions only: ratio (a) 1.96, P = 3.7e-05
variant, with the two controls and MS10-170307: ratio (a) 2.03, P = 8.3e-06
Controls alone: {'MS21-180719-155941': {'pre': 17.9, 'cooling': 5.6, 'post': 4.3}, 'MS22-180719-122813': {'pre': 6.5, 'cooling': 7.7, 'post': 9.7}}
Result: partially reproduced. Errors rise from 6.67% before cooling to 15.23% during it, higher in 40 of 45 sessions (P = 3e-07), a ratio of 2.28 (a), 2.31 (b) or 2.53 (c), below 3.0 +- 0.3 under every reading. Within rats the increase is significant for MS10 (P = 0.0078), MS12 (P = 0.039) and MS21 (P = 0.027) but not for MS22 (P = 0.064, 10 of 12 sessions higher); MS13 has 5 sessions, all higher (P = 0.0625, the smallest an exact test allows).
fig6e-errors-recover¶
The error rate after cooling did not differ from the error rate before it.
In [25]:
p, n = wsr(er.post, er.pre)
print(f"post against pre over {n} sessions: post higher in {(er.post > er.pre).sum()}, lower in {(er.post < er.pre).sum()}, "
f"signed-rank P = {p:.3g}; median difference {np.median(er.post - er.pre):.2f} points; means {er.pre.mean():.2f} and {er.post.mean():.2f}%")
p2, n2 = wsr(er.cooling, er.post)
print(f"cooling against post: P = {p2:.2g}")
post against pre over 45 sessions: post higher in 20, lower in 23, signed-rank P = 0.837; median difference -0.13 points; means 6.67 and 6.49% cooling against post: P = 8e-08
Result: reproduced. After cooling the error rate is 6.49% against 6.67% before (median difference -0.13 points, P = 0.837 over 45 sessions), and it is lower than during cooling (P = 8e-08).
fig6c-worst-trials-4-to-8¶
The worst memory performance occurred between trials 4 and 8 after the onset of cooling, when MS temperature reached its minimum.
In [26]:
me = mean_series(BEH, "error_smooth")
mt_all = mean_series(TEMP, "temperature")
mf_all = mean_series(BEH, "theta_freq")
onset0 = C.ALIGN_ONSET - 1 # 0-based index of the first cooling trial, which is trial 1 after onset
k_after = np.arange(1, 41)
idx = onset0 + k_after - 1
peak_k = int(k_after[np.nanargmax(me[idx])])
print("Peak of the mean smoothed error series over trials 1 to 40 after onset: trial", peak_k, f"(error rate {np.nanmax(me[idx]):.3f})")
print("Trial of the minimum mean temperature (TEMP):", int(k_after[np.nanargmin(mt_all[idx])]), "; of the minimum mean theta frequency (BEH):", int(k_after[np.nanargmin(mf_all[idx])]))
print("Mean smoothed error rate by trial after onset (1 to 20):", np.round(me[idx][:20], 3).tolist())
print("Mean smoothed error rate over the 40 trials before onset:", round(float(np.nanmean(me[onset0 - 40:onset0])), 3))
raw = mean_series(BEH, "error")
print("Unsmoothed mean error by trial after onset (1 to 20):", np.round(raw[idx][:20], 3).tolist())
print("Mean smoothed error rate by trial after onset (21 to 40):", np.round(me[idx][20:40], 3).tolist())
n_ses = np.isfinite(stack(BEH, "error_smooth")).sum(axis=0)
print("Sessions contributing, by trial after onset:", {int(k_): int(n_ses[onset0 + k_ - 1]) for k_ in [1, 5, 7, 10, 15, 20, 25, 30, 35, 40]})
for a in sorted({subj(k) for k in BEH}):
m_a = mean_series([k for k in BEH if subj(k) == a], "error_smooth")
print(f" {a}: peak at trial {int(k_after[np.nanargmax(m_a[idx])])}")
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(np.arange(-40, 41), me[onset0 - 41:onset0 + 40], "r"); ax.axvline(0, color="k", lw=0.5)
ax.set_xlabel("trial after onset (1 = first cooling trial)"); ax.set_ylabel("smoothed error rate")
plt.tight_layout(); plt.show()
print("Plotted range:", np.nanmin(me[onset0 - 41:onset0 + 40]).round(3), "to", np.nanmax(me[onset0 - 41:onset0 + 40]).round(3))
Peak of the mean smoothed error series over trials 1 to 40 after onset: trial 7 (error rate 0.193)
Trial of the minimum mean temperature (TEMP): 7 ; of the minimum mean theta frequency (BEH): 7
Mean smoothed error rate by trial after onset (1 to 20): [0.092, 0.116, 0.14, 0.158, 0.173, 0.185, 0.193, 0.193, 0.183, 0.167, 0.155, 0.15, 0.151, 0.153, 0.157, 0.16, 0.165, 0.174, 0.183, 0.186]
Mean smoothed error rate over the 40 trials before onset: 0.066
Unsmoothed mean error by trial after onset (1 to 20): [0.044, 0.178, 0.133, 0.2, 0.156, 0.156, 0.178, 0.333, 0.178, 0.133, 0.089, 0.156, 0.178, 0.156, 0.111, 0.222, 0.156, 0.089, 0.178, 0.333]
Mean smoothed error rate by trial after onset (21 to 40): [0.182, 0.172, 0.165, 0.164, 0.163, 0.16, 0.153, 0.143, 0.132, 0.122, 0.114, 0.108, 0.103, 0.098, 0.091, 0.084, 0.077, 0.076, 0.078, 0.083]
Sessions contributing, by trial after onset: {1: 45, 5: 45, 7: 45, 10: 45, 15: 45, 20: 45, 25: 45, 30: 45, 35: 45, 40: 45}
MS10: peak at trial 21
MS12: peak at trial 26
MS13: peak at trial 14
MS21: peak at trial 20
MS22: peak at trial 19
Plotted range: 0.045 to 0.193
Result: reproduced. The mean smoothed error rate peaks at trial 7 after the onset of cooling (0.193, against 0.066 over the 40 trials before), the same trial as the minimum of the mean temperature and of the mean theta frequency. The peak is narrow: trials 7 and 8 are at 0.193, and a second rise reaches 0.186 at trial 20 before the rate falls to 0.083 at trial 40. All 45 sessions contribute at every trial from 1 to 40. In single rats the peak of the mean falls later, at trials 14 to 26.
fig6f-errors-temperature¶
The incidence of errors correlated negatively with MS temperature (R = -0.35, P < 0.01).
In [27]:
win = slice(9, 160) # positions 10 to 160, as the code's correlation matrix
def r_err(keys_err, keys_x, col, sl=slice(None)):
return C.pearson(mean_series(keys_err, "error_smooth")[sl], mean_series(keys_x, col)[sl])
res_6f = {"(a) TEMP, positions 10 to 160": r_err(TEMP, TEMP, "temperature", win),
"(b) TEMP, all positions": r_err(TEMP, TEMP, "temperature"),
"(c) BEH with MS10 curves, positions 10 to 160": r_err(BEH, BEH, "temperature", win)}
for k, v in res_6f.items():
print(f"{k}: R = {v[0]:.3f}, P = {v[1]:.2g}, positions = {v[2]}")
v = r_err(BEH, TEMP, "temperature", win)
print(f"variant, errors over BEH and temperature over TEMP, positions 10 to 160: R = {v[0]:.3f}, P = {v[1]:.2g}")
pooled_e = pd.concat([TR[k] for k in TEMP])
v = C.pearson(pooled_e.error.astype(float), pooled_e.temperature)
print(f"variant, single trials pooled over TEMP (error 0/1): R = {v[0]:.3f}, P = {v[1]:.2g}, n = {v[2]}")
print("Per animal, positions 10 to 160:", {a: tuple(np.round(r_err([k for k in TEMP if subj(k) == a], [k for k in TEMP if subj(k) == a], "temperature", win)[:2], 4))
for a in sorted({subj(k) for k in TEMP})})
print("Within -0.35 +- 0.035 with P < 0.01 under:", [k for k, v in res_6f.items() if within(v[0], -0.35, 0.035) and v[1] < 0.01])
(a) TEMP, positions 10 to 160: R = -0.838, P = 5.9e-41, positions = 151
(b) TEMP, all positions: R = -0.789, P = 2.5e-43, positions = 198
(c) BEH with MS10 curves, positions 10 to 160: R = -0.868, P = 3.1e-47, positions = 151
variant, errors over BEH and temperature over TEMP, positions 10 to 160: R = -0.881, P = 3.4e-50
variant, single trials pooled over TEMP (error 0/1): R = -0.127, P = 5.5e-20, n = 5187
Per animal, positions 10 to 160: {'MS12': (np.float64(-0.298), np.float64(0.0004)), 'MS13': (np.float64(-0.5125), np.float64(0.0)), 'MS21': (np.float64(-0.8597), np.float64(0.0)), 'MS22': (np.float64(-0.5577), np.float64(0.0))}
Within -0.35 +- 0.035 with P < 0.01 under: []
Result: partially reproduced. Errors correlate negatively with MS temperature, and more strongly than reported: R = -0.838 (P = 5.9e-41) over positions 10 to 160, -0.789 over all positions, and -0.868 with the MS10 curves, none within -0.35 +- 0.035. Over single trials pooled R = -0.127 (P = 5.5e-20). Per rat R ranges from -0.298 to -0.860.
fig6g-errors-theta-frequency¶
The incidence of errors correlated negatively with theta frequency (R = -0.35, P < 0.01).
In [28]:
res_6g = {"(a) BEH, positions 10 to 160": r_err(BEH, BEH, "theta_freq", win), "(b) BEH, all positions": r_err(BEH, BEH, "theta_freq")}
for k, v in res_6g.items():
print(f"{k}: R = {v[0]:.3f}, P = {v[1]:.2g}, positions = {v[2]}")
v = r_err(TEMP, TEMP, "theta_freq", win); print(f"variant, TEMP sessions, positions 10 to 160: R = {v[0]:.3f}, P = {v[1]:.2g}")
v = r_err(BEH, BEH, "theta_freq_10", win); print(f"variant, theta over frames above 10 cm/s: R = {v[0]:.3f}, P = {v[1]:.2g}")
print("Per animal, positions 10 to 160:", {a: tuple(np.round(r_err([k for k in BEH if subj(k) == a], [k for k in BEH if subj(k) == a], "theta_freq", win)[:2], 4))
for a in sorted({subj(k) for k in BEH})})
print("Within -0.35 +- 0.035 with P < 0.01 under:", [k for k, v in res_6g.items() if within(v[0], -0.35, 0.035) and v[1] < 0.01])
(a) BEH, positions 10 to 160: R = -0.868, P = 3.9e-47, positions = 151
(b) BEH, all positions: R = -0.733, P = 1.2e-34, positions = 198
variant, TEMP sessions, positions 10 to 160: R = -0.887, P = 8e-52
variant, theta over frames above 10 cm/s: R = -0.868, P = 4.6e-47
Per animal, positions 10 to 160: {'MS10': (np.float64(-0.7121), np.float64(0.0)), 'MS12': (np.float64(-0.3406), np.float64(0.0)), 'MS13': (np.float64(-0.5716), np.float64(0.0)), 'MS21': (np.float64(-0.8289), np.float64(0.0)), 'MS22': (np.float64(-0.5369), np.float64(0.0))}
Within -0.35 +- 0.035 with P < 0.01 under: []
Result: partially reproduced. Errors correlate negatively with theta frequency at R = -0.868 (P = 3.9e-47) over positions 10 to 160 and R = -0.733 over all positions, stronger than -0.35 +- 0.035. Per rat R ranges from -0.341 to -0.829.
fig6h-errors-theta-power¶
The incidence of errors correlated less with theta power; only 2 of 5 animals had a significant correlation (R beyond -0.3, P < 0.01).
In [29]:
rows = []
for a in sorted({subj(k) for k in BEH}):
ks = [k for k in BEH if subj(k) == a]
for col in ["theta_power", "theta_freq", "temperature", "speed"]:
if col == "temperature" and a == "MS10":
continue
r, p, n = r_err(ks, ks, col)
rows.append({"animal": a, "measure": col, "R": r, "P": p, "positions": n})
pa6 = pd.DataFrame(rows)
print(pa6.pivot(index="animal", columns="measure", values="R").round(3).to_string())
print(pa6.pivot(index="animal", columns="measure", values="P").to_string(float_format=lambda x: f"{x:.2g}"))
pw = pa6[pa6.measure == "theta_power"]
n_sig = int(((pw.R < -0.3) & (pw.P < 0.01)).sum())
print("\nAnimals with errors and theta power at R < -0.3 and P < 0.01:", n_sig, pw[(pw.R < -0.3) & (pw.P < 0.01)].animal.tolist())
for thr in [0.0, -0.2, -0.3]:
print(f" with R < {thr} and P < 0.01: {int(((pw.R < thr) & (pw.P < 0.01)).sum())}")
measure speed temperature theta_freq theta_power animal MS10 0.073 NaN -0.641 -0.425 MS12 -0.055 -0.303 -0.307 -0.038 MS13 -0.493 -0.525 -0.590 0.063 MS21 0.216 -0.833 -0.802 -0.202 MS22 0.247 -0.550 -0.524 0.325 measure speed temperature theta_freq theta_power animal MS10 0.32 NaN 2.9e-23 1.1e-09 MS12 0.49 0.00011 9.3e-05 0.64 MS13 8e-11 2.9e-12 8e-16 0.44 MS21 0.0034 7.9e-48 6.3e-42 0.0063 MS22 0.001 3.7e-15 1.2e-13 1.2e-05 Animals with errors and theta power at R < -0.3 and P < 0.01: 1 ['MS10'] with R < 0.0 and P < 0.01: 2 with R < -0.2 and P < 0.01: 2 with R < -0.3 and P < 0.01: 1
Result: partially reproduced. One rat (MS10, R = -0.425, P = 1.1e-09) has errors correlated with theta power beyond -0.3 at P < 0.01, against two in the paper. A second (MS21, R = -0.202, P = 0.0063) is significant with a weaker correlation, and MS22 is significant with the opposite sign (R = 0.325).
fig6i-errors-not-speed¶
The incidence of errors did not correlate with the running speed of the rat (R about 0.1, P > 0.05 for all animals).
In [30]:
sp6 = pa6[pa6.measure == "speed"].set_index("animal")
print(sp6[["R", "P", "positions"]].to_string(float_format=lambda x: f"{x:.3g}"))
print("Animals with P < 0.05:", sp6.index[sp6.P < 0.05].tolist())
v = r_err(BEH, BEH, "speed"); print(f"all BEH sessions together: R = {v[0]:.3f}, P = {v[1]:.2g}")
v = r_err(BEH, BEH, "speed", win); print(f"all BEH sessions, positions 10 to 160: R = {v[0]:.3f}, P = {v[1]:.2g}")
R P positions animal MS10 0.0729 0.319 189 MS12 -0.055 0.494 157 MS13 -0.493 8.04e-11 154 MS21 0.216 0.00344 181 MS22 0.247 0.00102 174 Animals with P < 0.05: ['MS13', 'MS21', 'MS22'] all BEH sessions together: R = 0.093, P = 0.19 all BEH sessions, positions 10 to 160: R = -0.055, P = 0.5
Result: not reproduced. Errors correlate with running speed in three of five rats at P < 0.05: MS13 (R = -0.493, P = 8e-11), MS21 (R = 0.216, P = 0.0034) and MS22 (R = 0.247, P = 0.001); MS10 and MS12 do not (P = 0.32 and 0.49). The signs differ between rats, and over all sessions together R = 0.093 (P = 0.19).
fig4e-speed-decreases-within-session¶
Running speed decreased relatively monotonically within a session.
In [31]:
rho = pd.Series({k: stats.spearmanr(np.arange(len(TR[k])), TR[k].speed, nan_policy="omit")[0] for k in BEH})
p, n = wsr(rho)
print(f"Spearman rho of trial speed with trial number: negative in {(rho < 0).sum()} of {len(rho)} sessions ({100 * (rho < 0).mean():.0f}%), "
f"median {rho.median():.3f}, signed-rank P = {p:.2g}")
print("By subject: negative / sessions:", {a: f"{(rho[[k for k in BEH if subj(k) == a]] < 0).sum()}/{len([k for k in BEH if subj(k) == a])}" for a in sorted({subj(k) for k in BEH})})
ms = mean_series(BEH, "speed")
print("Mean aligned speed (cm/s) in blocks of 20 positions from 41:", [round(float(np.nanmean(ms[i:i + 20])), 1) for i in range(40, 200, 20)])
rho_run = pd.Series({k: stats.spearmanr(np.arange(len(TR[k])), TR[k].speed_run, nan_policy="omit")[0] for k in BEH})
print(f"variant, speed over frames above 10 cm/s: negative in {(rho_run < 0).sum()} of {len(rho_run)}, median {rho_run.median():.3f}")
Spearman rho of trial speed with trial number: negative in 44 of 45 sessions (98%), median -0.452, signed-rank P = 4e-13
By subject: negative / sessions: {'MS10': '8/8', 'MS12': '8/8', 'MS13': '5/5', 'MS21': '12/12', 'MS22': '11/12'}
Mean aligned speed (cm/s) in blocks of 20 positions from 41: [86.0, 80.6, 79.6, 77.4, 75.6, 75.3, 75.4, 73.7]
variant, speed over frames above 10 cm/s: negative in 44 of 45, median -0.458
Result: reproduced. Trial speed falls with trial number in 44 of 45 sessions (median Spearman rho -0.452, P = 4e-13), in every rat, and the aligned mean falls from 86.0 to 73.7 cm/s over the session.
fig5-speed-control-cooling¶
Mean running speed was 87.3 +- 8.0 cm/s in control trials and 86.2 +- 9.9 cm/s during MS cooling.
In [32]:
F5B = [k for k in S.index[S.in_fig5_list] if beh[k]["has_state"]]
print("Sessions:", len(F5B), F5B)
for col, name in [("speed", "(a) all frames"), ("speed_run", "(b) frames above 10 cm/s")]:
ctl = pd.Series({k: beh[k]["trials"].loc[beh[k]["trials"].state.isin(["Pre-Cooling", "Post-Cooling"]), col].mean() for k in F5B})
cool = pd.Series({k: beh[k]["trials"].loc[beh[k]["trials"].state == "Cooling on", col].mean() for k in F5B})
print(f"{name}: control {ctl.mean():.1f} +- {ctl.std():.1f}, cooling {cool.mean():.1f} +- {cool.std():.1f} cm/s (mean +- SD over sessions); "
f"within: {within(ctl.mean(), 87.3, 8.0)}, {within(cool.mean(), 86.2, 9.9)}; cooling minus control P = {wsr(cool, ctl)[0]:.2g}")
if col == "speed":
print(" by subject, control:", ctl.groupby(lambda k: subj(k)).mean().round(1).to_dict(), "cooling:", cool.groupby(lambda k: subj(k)).mean().round(1).to_dict())
ctl_b = pd.Series({k: TR[k].loc[TR[k].state.isin(["Pre-Cooling", "Post-Cooling"]), "speed"].mean() for k in BEH})
cool_b = pd.Series({k: TR[k].loc[TR[k].state == "Cooling on", "speed"].mean() for k in BEH})
print(f"variant, all BEH sessions, all frames: control {ctl_b.mean():.1f} +- {ctl_b.std():.1f}, cooling {cool_b.mean():.1f} +- {cool_b.std():.1f}")
Sessions: 15 ['MS10-170307-154746', 'MS10-170317-153237', 'MS12-170714-122034', 'MS12-170716-172307', 'MS12-170717-111614', 'MS12-170719-095305', 'MS13-171128-113924', 'MS13-171129-105507', 'MS13-171130-121758', 'MS21-180625-153927', 'MS21-180628-155921', 'MS21-180629-110332', 'MS22-180628-120341', 'MS22-180629-110319', 'MS22-180711-112912']
(a) all frames: control 78.9 +- 10.3, cooling 74.9 +- 14.0 cm/s (mean +- SD over sessions); within: False, False; cooling minus control P = 0.041
by subject, control: {'MS10': 65.0, 'MS12': 71.7, 'MS13': 76.8, 'MS21': 92.9, 'MS22': 85.7} cooling: {'MS10': 61.1, 'MS12': 64.6, 'MS13': 68.7, 'MS21': 89.7, 'MS22': 89.3}
(b) frames above 10 cm/s: control 79.6 +- 10.3, cooling 76.0 +- 13.6 cm/s (mean +- SD over sessions); within: True, False; cooling minus control P = 0.035
variant, all BEH sessions, all frames: control 79.4 +- 7.5, cooling 78.6 +- 10.0
Result: not reproduced. Over the 15 sessions of the code's Figure 5 list, neither reading gives both values within the paper's spreads. Over all frames the mean speed is 78.9 cm/s in control trials and 74.9 cm/s during cooling, both outside. Over frames above 10 cm/s control is 79.6 cm/s, within 87.3 +- 8.0, and cooling 76.0 cm/s, just below 86.2 +- 9.9. Speed is lower during cooling than in control (P = 0.035 to 0.041), where the paper's values are close to each other.
Unit claims: selection¶
Units of the UNIT sessions typed as above; place fields of putative pyramidal cells by the automatic criteria of the Methods. The counts below are what each unit claim starts from.
In [33]:
UC = uc[uc.session.isin(UNIT)].copy()
UC["subject"] = UC.session.map(subj)
F = fields[fields.session.isin(UNIT)].copy()
print("UNIT sessions:", len(UNIT), "; units", len(UC), UC.cell_type.value_counts().to_dict())
print("Place fields:", len(F), "in", F.groupby(["session", "unit"]).ngroups, "place cells; by subject", F.groupby("subject").size().to_dict())
F3 = F[(F.n_pre >= 50) & (F.n_cooling >= 50) & (F.n_post >= 50)].copy()
F3 = F3[F3[["fp_pre", "fp_cooling", "fp_post", "flfp_pre", "flfp_cooling", "flfp_post"]].notna().all(axis=1)]
print("Fields with at least 50 in-field spikes in each of pre, cooling and post:", int(((F.n_pre >= 50) & (F.n_cooling >= 50) & (F.n_post >= 50)).sum()),
"; of them with an autocorrelogram peak (fp) in all three states:", len(F3))
print("Dropped for spikes:", len(F) - int(((F.n_pre >= 50) & (F.n_cooling >= 50) & (F.n_post >= 50)).sum()), "; for no fp peak:",
int(((F.n_pre >= 50) & (F.n_cooling >= 50) & (F.n_post >= 50)).sum()) - len(F3))
print("F3 by subject:", F3.groupby("subject").size().to_dict(), "; by session:", F3.groupby("session").size().to_dict())
UNIT sessions: 11 ; units 987 {'pyramidal': 605, 'unclassified': 253, 'interneuron': 129}
Place fields: 643 in 374 place cells; by subject {'MS12': 95, 'MS13': 51, 'MS21': 249, 'MS22': 248}
Fields with at least 50 in-field spikes in each of pre, cooling and post: 443 ; of them with an autocorrelogram peak (fp) in all three states: 361
Dropped for spikes: 200 ; for no fp peak: 82
F3 by subject: {'MS12': 59, 'MS13': 19, 'MS21': 145, 'MS22': 138} ; by session: {'MS12-170716-172307': 35, 'MS12-170719-095305': 24, 'MS13-171128-113924': 11, 'MS13-171201-130527': 8, 'MS21-180625-153927': 34, 'MS21-180628-155921': 63, 'MS21-180629-110332': 48, 'MS22-180628-120341': 16, 'MS22-180629-110319': 59, 'MS22-180711-112912': 27, 'MS22-180720-110055': 36}
suppl2-place-cell-rates-decrease¶
Within-field and peak firing rates of place cells decreased moderately during MS cooling (11% and 12%; P < 1e-9 and P < 1e-11).
In [34]:
Fr = F[(F.time_pre > 0) & (F.time_cooling > 0) & (F.infield_rate_pre > 0) & (F.peak_rate_pre > 0)].copy()
print("Fields with time in the field in both states and a pre rate above 0:", len(Fr), "of", len(F))
res_r = {}
for m in ["infield_rate", "peak_rate"]:
pc = 100 * (1 - Fr[f"{m}_cooling"] / Fr[f"{m}_pre"])
a, b, c = pc.median(), 100 * (1 - Fr[f"{m}_cooling"].mean() / Fr[f"{m}_pre"].mean()), pc.mean()
p, n = wsr(Fr[f"{m}_cooling"], Fr[f"{m}_pre"])
res_r[m] = (a, b, c, p)
print(f"{m}: pre mean {Fr[m + '_pre'].mean():.2f} Hz, cooling {Fr[m + '_cooling'].mean():.2f} Hz; decrease (a) median {a:.1f}%, (b) ratio of means {b:.1f}%, "
f"(c) mean {c:.1f}%; signed-rank P = {p:.2g} (n = {n}); lower in {(Fr[m + '_cooling'] < Fr[m + '_pre']).sum()}")
pp, nn = wsr(Fr[f"{m}_post"], Fr[f"{m}_cooling"])
print(f" post against cooling: post mean {Fr[m + '_post'].mean():.2f} Hz, P = {pp:.2g}")
print(" by subject, reading (b):", {a_: round(100 * (1 - g[f'{m}_cooling'].mean() / g[f'{m}_pre'].mean()), 1) for a_, g in Fr.groupby("subject")})
for mins in [20, 50, 100]:
g = Fr[(Fr.n_pre >= mins) & (Fr.n_cooling >= mins)]
print(f"variant, fields with at least {mins} spikes in pre and cooling (n = {len(g)}): within-field decrease (b) {100 * (1 - g.infield_rate_cooling.mean() / g.infield_rate_pre.mean()):.1f}%, "
f"peak {100 * (1 - g.peak_rate_cooling.mean() / g.peak_rate_pre.mean()):.1f}%")
pyr = UC[UC.cell_type == "pyramidal"]
print(f"variant, all putative pyramidal cells, rate over trial frames above 10 cm/s (n = {len(pyr)}): pre {pyr.rate_pre.mean():.2f}, cooling {pyr.rate_cooling.mean():.2f} Hz, "
f"decrease {100 * (1 - pyr.rate_cooling.mean() / pyr.rate_pre.mean()):.1f}%, P = {wsr(pyr.rate_cooling, pyr.rate_pre)[0]:.2g}")
print("Readings with both decreases within tolerance (11 +- 1.1, 12 +- 1.2):",
[r for i, r in enumerate("abc") if within(res_r["infield_rate"][i], 11, 1.1) and within(res_r["peak_rate"][i], 12, 1.2)])
Fields with time in the field in both states and a pre rate above 0: 636 of 643
infield_rate: pre mean 10.94 Hz, cooling 10.09 Hz; decrease (a) median 3.7%, (b) ratio of means 7.7%, (c) mean -111.2%; signed-rank P = 0.0024 (n = 636); lower in 343
post against cooling: post mean 11.55 Hz, P = 4.5e-09
by subject, reading (b): {'MS12': np.float64(-10.1), 'MS13': np.float64(27.8), 'MS21': np.float64(12.8), 'MS22': np.float64(4.2)}
peak_rate: pre mean 19.56 Hz, cooling 18.25 Hz; decrease (a) median 5.3%, (b) ratio of means 6.7%, (c) mean -70.5%; signed-rank P = 0.0013 (n = 636); lower in 352
post against cooling: post mean 20.39 Hz, P = 3.1e-08
by subject, reading (b): {'MS12': np.float64(-1.3), 'MS13': np.float64(25.7), 'MS21': np.float64(9.8), 'MS22': np.float64(2.1)}
variant, fields with at least 20 spikes in pre and cooling (n = 559): within-field decrease (b) 10.0%, peak 8.9%
variant, fields with at least 50 spikes in pre and cooling (n = 463): within-field decrease (b) 10.8%, peak 9.9%
variant, fields with at least 100 spikes in pre and cooling (n = 332): within-field decrease (b) 10.9%, peak 10.8%
variant, all putative pyramidal cells, rate over trial frames above 10 cm/s (n = 605): pre 2.59, cooling 2.40 Hz, decrease 7.5%, P = 0.0088
Readings with both decreases within tolerance (11 +- 1.1, 12 +- 1.2): []
Result: partially reproduced. Within-field and peak rates of 636 place fields are lower during cooling (P = 0.0024 and P = 0.0013), by 7.7% and 6.7% as ratios of means and 3.7% and 5.3% as medians, short of 11% and 12% and of the P < 0.001 of the criterion. Fields with at least 50 spikes in both states come closer (10.8% and 9.9%). The decrease differs between rats (MS12 rises by 10.1%).
disc-interneuron-rates-decrease¶
Firing rates of interneurons decreased during MS cooling by about 14%.
In [35]:
inn = UC[(UC.cell_type == "interneuron") & (UC.rate_pre > 0)]
pc = 100 * (1 - inn.rate_cooling / inn.rate_pre)
p, n = wsr(inn.rate_cooling, inn.rate_pre)
print(f"Putative interneurons: {len(inn)}; rate pre {inn.rate_pre.mean():.2f}, cooling {inn.rate_cooling.mean():.2f}, post {inn.rate_post.mean():.2f} Hz")
print(f"decrease (a) median {pc.median():.1f}%, (b) ratio of means {100 * (1 - inn.rate_cooling.mean() / inn.rate_pre.mean()):.1f}%; lower in {(inn.rate_cooling < inn.rate_pre).sum()}; signed-rank P = {p:.2g} (n = {n})")
print("By subject, reading (b):", {a: round(100 * (1 - g.rate_cooling.mean() / g.rate_pre.mean()), 1) for a, g in inn.groupby("subject")})
for thr in [0.35, 0.5]:
t_ = (UC.ttp <= thr) | (UC.rate_trials > C.WIDE_INT_RATE)
g = UC[t_ & UC.ttp.notna() & (UC.rate_pre > 0)]
print(f"variant, trough-to-peak threshold {thr} ms (n = {len(g)}): decrease (a) {np.median(100 * (1 - g.rate_cooling / g.rate_pre)):.1f}%, "
f"(b) {100 * (1 - g.rate_cooling.mean() / g.rate_pre.mean()):.1f}%, P = {wsr(g.rate_cooling, g.rate_pre)[0]:.2g}")
nar = UC[(UC.ttp <= C.NARROW_TTP_MS) & (UC.rate_pre > 0)]
print(f"variant, narrow waveform only (n = {len(nar)}): decrease (b) {100 * (1 - nar.rate_cooling.mean() / nar.rate_pre.mean()):.1f}%, P = {wsr(nar.rate_cooling, nar.rate_pre)[0]:.2g}")
Putative interneurons: 127; rate pre 20.97, cooling 19.06, post 22.25 Hz
decrease (a) median 6.5%, (b) ratio of means 9.1%; lower in 75; signed-rank P = 0.0017 (n = 127)
By subject, reading (b): {'MS12': np.float64(-21.7), 'MS13': np.float64(9.9), 'MS21': np.float64(18.4), 'MS22': np.float64(4.2)}
variant, trough-to-peak threshold 0.35 ms (n = 110): decrease (a) 7.4%, (b) 9.6%, P = 0.00047
variant, trough-to-peak threshold 0.5 ms (n = 155): decrease (a) 6.5%, (b) 9.4%, P = 0.0016
variant, narrow waveform only (n = 96): decrease (b) 9.9%, P = 0.0037
Result: partially reproduced. Putative interneurons fire less during cooling (127 units, lower in 75, P = 0.0017), by 6.5% (median) or 9.1% (ratio of means), below 14 +- 1.4%; other waveform thresholds give 9.4% to 9.9%.
suppl2d-interneuron-oscillation-slows¶
The theta oscillation frequency of putative fast-spiking interneurons decreased during MS cooling, commensurate with the LFP theta.
In [36]:
io = UC[(UC.cell_type == "interneuron")][["session", "subject", "osc_pre", "osc_cooling", "osc_post"]].dropna(subset=["osc_pre", "osc_cooling"])
p, n = wsr(io.osc_cooling, io.osc_pre)
print(f"Interneurons with an autocorrelogram peak in pre and cooling: {len(io)} of {int((UC.cell_type == 'interneuron').sum())}")
print(f"oscillation frequency pre {io.osc_pre.mean():.2f}, cooling {io.osc_cooling.mean():.2f}, post {io.osc_post.mean():.2f} Hz; lower in {(io.osc_cooling < io.osc_pre).sum()}; signed-rank P = {p:.2g} (n = {n})")
print(f"decrease {100 * (1 - io.osc_cooling.mean() / io.osc_pre.mean()):.1f}%; LFP theta decrease in the same sessions (reading (a)) {dec_a[UNIT].mean():.1f}%")
print("By subject:", io.groupby("subject")[["osc_pre", "osc_cooling"]].mean().round(2).to_dict("index"))
nar = UC[(UC.ttp <= C.NARROW_TTP_MS)].dropna(subset=["osc_pre", "osc_cooling"])
print(f"variant, narrow waveform only (n = {len(nar)}): pre {nar.osc_pre.mean():.2f}, cooling {nar.osc_cooling.mean():.2f} Hz, P = {wsr(nar.osc_cooling, nar.osc_pre)[0]:.2g}")
Interneurons with an autocorrelogram peak in pre and cooling: 114 of 129
oscillation frequency pre 9.55, cooling 8.66, post 9.42 Hz; lower in 105; signed-rank P = 8.7e-16 (n = 114)
decrease 9.4%; LFP theta decrease in the same sessions (reading (a)) 8.8%
By subject: {'MS12': {'osc_pre': 9.52, 'osc_cooling': 8.23}, 'MS13': {'osc_pre': 9.12, 'osc_cooling': 8.26}, 'MS21': {'osc_pre': 9.47, 'osc_cooling': 8.33}, 'MS22': {'osc_pre': 9.66, 'osc_cooling': 8.97}}
variant, narrow waveform only (n = 83): pre 9.32, cooling 8.43 Hz, P = 3.1e-11
Result: reproduced. The oscillation frequency of putative interneurons falls from 9.55 to 8.66 Hz during cooling (lower in 105 of 114, P = 8.7e-16), a decrease of 9.4%, close to the 8.8% decrease of LFP theta in the same sessions.
suppl1g-waveforms-unchanged¶
The waveforms of hippocampal action potentials were not affected by MS cooling.
In [37]:
wv = UC[(UC.nwf_pre >= 30) & (UC.nwf_cooling >= 30)]
win_t = pd.DataFrame({k: units[k]["waveform_windows"] for k in UNIT}).T
print("Waveform windows (s) and the mean temperature of the coldest cooling trial:"); print(win_t.to_string())
print(f"\nUnits with at least 30 spikes in each window: {len(wv)} of {len(UC)}")
for m in ["ttp", "amp"]:
p, n = wsr(wv[f"{m}_cooling"], wv[f"{m}_pre"])
d = wv[f"{m}_cooling"] - wv[f"{m}_pre"]
print(f"{m}: pre median {wv[m + '_pre'].median():.3f}, cooling median {wv[m + '_cooling'].median():.3f}; median change {d.median():.4f} "
f"({100 * (d / wv[m + '_pre']).median():.1f}%), mean change {d.mean():.4f}; higher in {(d > 0).sum()}, lower in {(d < 0).sum()}, equal in {(d == 0).sum()}; "
f"signed-rank P = {p:.2g} (n = {n})")
print("Spikes averaged per unit, median (IQR): pre", wv.nwf_pre.median(), wv.nwf_pre.quantile([.25, .75]).tolist(),
"; cooling", wv.nwf_cooling.median(), wv.nwf_cooling.quantile([.25, .75]).tolist(), "; fewer in cooling than pre in", int((wv.nwf_cooling < wv.nwf_pre).sum()), "units")
cap = wv[(wv.nsp_pre >= 300) & (wv.nsp_cooling >= 300)]
print(f"units with at least 300 spikes in both windows, so both means use the 300-spike cap (n = {len(cap)}): amplitude median change "
f"{100 * ((cap.amp_cooling - cap.amp_pre) / cap.amp_pre).median():.1f}%, P = {wsr(cap.amp_cooling, cap.amp_pre)[0]:.2g}; "
f"ttp higher in {(cap.ttp_cooling > cap.ttp_pre).sum()}, lower in {(cap.ttp_cooling < cap.ttp_pre).sum()}, P = {wsr(cap.ttp_cooling, cap.ttp_pre)[0]:.2g}")
print("Units whose largest channel moved between windows:", int((wv.chan_pre != wv.chan_cooling).sum()))
for t_ in ["pyramidal", "interneuron"]:
g = wv[wv.cell_type == t_]
print(f" {t_} (n = {len(g)}): ttp P = {wsr(g.ttp_cooling, g.ttp_pre)[0]:.2g}, amplitude change {100 * ((g.amp_cooling - g.amp_pre) / g.amp_pre).median():.1f}%, P = {wsr(g.amp_cooling, g.amp_pre)[0]:.2g}")
for mins in [50, 100]:
g = UC[(UC.nwf_pre >= mins) & (UC.nwf_cooling >= mins)]
print(f"variant, at least {mins} spikes (n = {len(g)}): ttp P = {wsr(g.ttp_cooling, g.ttp_pre)[0]:.2g}, amplitude P = {wsr(g.amp_cooling, g.amp_pre)[0]:.2g}, "
f"median amplitude change {100 * ((g.amp_cooling - g.amp_pre) / g.amp_pre).median():.1f}%")
print("Trough-to-peak is measured on 50 us samples, so a change below one sample (0.05 ms) cannot be seen.")
Waveform windows (s) and the mean temperature of the coldest cooling trial:
pre cooling cooling_window_temp
MS12-170716-172307 (7022.75615, 7142.75615) (8817.10535, 8937.10535) 18.867797
MS12-170719-095305 (19127.8269, 19247.8269) (21038.440225, 21158.440225) 20.356578
MS13-171128-113924 (13651.6122, 13771.6122) (14220.337775, 14340.337775) 20.633155
MS13-171201-130527 (12333.46545, 12453.46545) (12801.729, 12921.729) 20.393548
MS21-180625-153927 (58.55275, 178.55275) (514.7125249999999, 634.7125249999999) 17.052737
MS21-180628-155921 (44.8137, 164.81369999999998) (527.768325, 647.768325) 20.245797
MS21-180629-110332 (6379.86145, 6499.86145) (6806.29755, 6926.29755) 19.010863
MS22-180628-120341 (9768.0003, 9888.0003) (10302.256525, 10422.256525) 22.734461
MS22-180629-110319 (6357.4987, 6477.4987) (6819.020875, 6939.020875) 23.136223
MS22-180711-112912 (4837.7538, 4957.7538) (5335.5362000000005, 5455.5362000000005) 23.898493
MS22-180720-110055 (13596.84975, 13716.84975) (14212.535875000001, 14332.535875000001) 24.064289
Units with at least 30 spikes in each window: 646 of 987
ttp: pre median 0.750, cooling median 0.800; median change 0.0000 (0.0%), mean change 0.0213; higher in 294, lower in 143, equal in 209; signed-rank P = 4.9e-10 (n = 646)
amp: pre median 135.199, cooling median 138.450; median change 0.9783 (0.8%), mean change 1.1003; higher in 363, lower in 283, equal in 0; signed-rank P = 0.00057 (n = 646)
Spikes averaged per unit, median (IQR): pre 247.5 [118.25, 298.75] ; cooling 200.0 [98.0, 298.0] ; fewer in cooling than pre in 369 units
units with at least 300 spikes in both windows, so both means use the 300-spike cap (n = 208): amplitude median change 0.7%, P = 0.032; ttp higher in 78, lower in 33, P = 2.2e-05
Units whose largest channel moved between windows: 61
pyramidal (n = 524): ttp P = 3.7e-09, amplitude change 0.8%, P = 0.0015
interneuron (n = 122): ttp P = 0.066, amplitude change 0.7%, P = 0.2
variant, at least 50 spikes (n = 582): ttp P = 2.8e-10, amplitude P = 0.00027, median amplitude change 0.8%
variant, at least 100 spikes (n = 445): ttp P = 3.9e-08, amplitude P = 0.0026, median amplitude change 0.7%
Trough-to-peak is measured on 50 us samples, so a change below one sample (0.05 ms) cannot be seen.
Result: not reproduced. Across 646 units the waveforms differ between the pre window and the coldest cooling window: trough-to-peak by the signed-rank test (P = 4.9e-10, median change 0, since most units do not change by a whole 0.05 ms sample) and amplitude by a median 0.8% (P = 0.00057). Fewer spikes are averaged in the cooling window (median 200 against 247.5), but in the 208 units whose two means both use 300 spikes the changes remain (amplitude 0.7%, P = 0.032; trough-to-peak P = 2.2e-05), so the spike count does not explain them. Both changes are small, and the two windows start 7.1 to 31.8 minutes apart (printed above), so electrode drift cannot be told apart from cooling here; for interneurons alone neither differs (P = 0.066 and 0.2).
fig4b-fp-flfp-decrease¶
MS cooling decreased both the oscillation frequency of place cells (fp) and the LFP theta frequency in their fields (fLFP) by about the same degree, and both recovered (fLFP 8.4, 7.7, 8.5 Hz; fp 10.6, 9.5, 10.4 Hz before, during and after; n = 365).
In [38]:
exp4b = {"flfp_pre": (8.4, 0.38), "flfp_cooling": (7.7, 0.36), "flfp_post": (8.5, 0.29), "fp_pre": (10.6, 0.95), "fp_cooling": (9.5, 0.84), "fp_post": (10.4, 0.72)}
print(f"Fields: {len(F3)}")
for c, (v, tol) in exp4b.items():
print(f"{c}: mean {F3[c].mean():.2f} +- {F3[c].std():.2f} (expected {v} +- {tol}); within: {within(F3[c].mean(), v, tol)}")
for m in ["fp", "flfp"]:
p, n = wsr(F3[f"{m}_cooling"], F3[f"{m}_pre"])
p2, _ = wsr(F3[f"{m}_post"], F3[f"{m}_cooling"])
print(f"{m}: cooling lower than pre in {(F3[m + '_cooling'] < F3[m + '_pre']).sum()} of {n}, P = {p:.2g}; post higher than cooling P = {p2:.2g}; "
f"decrease {100 * (1 - F3[m + '_cooling'].mean() / F3[m + '_pre'].mean()):.1f}%")
print("\nBy subject (means):"); print(F3.groupby("subject")[list(exp4b)].mean().round(2).to_string())
print("By session (means):"); print(F3.groupby("session")[list(exp4b)].mean().round(2).to_string())
print("(oscillation_frequency needs 50 spikes, so 50 is the smallest minimum that can be tested)")
for mins in [100, 200]:
g = F[(F.n_pre >= mins) & (F.n_cooling >= mins) & (F.n_post >= mins)].dropna(subset=list(exp4b))
print(f"variant, at least {mins} spikes per state (n = {len(g)}):", g[list(exp4b)].mean().round(2).to_dict())
Fields: 361
flfp_pre: mean 8.54 +- 0.43 (expected 8.4 +- 0.38); within: True
flfp_cooling: mean 7.77 +- 0.37 (expected 7.7 +- 0.36); within: True
flfp_post: mean 8.60 +- 0.35 (expected 8.5 +- 0.29); within: True
fp_pre: mean 10.55 +- 1.19 (expected 10.6 +- 0.95); within: True
fp_cooling: mean 9.57 +- 1.08 (expected 9.5 +- 0.84); within: True
fp_post: mean 10.39 +- 0.98 (expected 10.4 +- 0.72); within: True
fp: cooling lower than pre in 311 of 361, P = 3.5e-44; post higher than cooling P = 4.5e-37; decrease 9.2%
flfp: cooling lower than pre in 352 of 361, P = 1.2e-60; post higher than cooling P = 6.8e-61; decrease 9.0%
By subject (means):
flfp_pre flfp_cooling flfp_post fp_pre fp_cooling fp_post
subject
MS12 7.90 7.56 8.44 10.44 9.06 10.20
MS13 9.19 8.14 9.02 11.21 10.17 11.17
MS21 8.51 7.57 8.38 10.66 9.54 10.42
MS22 8.75 8.01 8.82 10.39 9.75 10.31
By session (means):
flfp_pre flfp_cooling flfp_post fp_pre fp_cooling fp_post
session
MS12-170716-172307 7.98 7.59 8.63 10.56 9.03 10.26
MS12-170719-095305 7.79 7.50 8.17 10.26 9.11 10.12
MS13-171128-113924 9.18 8.07 8.90 10.84 9.99 10.85
MS13-171201-130527 9.21 8.23 9.20 11.72 10.43 11.62
MS21-180625-153927 8.59 7.45 8.53 10.62 9.30 10.24
MS21-180628-155921 8.43 7.57 8.41 10.19 9.51 10.51
MS21-180629-110332 8.56 7.67 8.24 11.30 9.75 10.44
MS22-180628-120341 8.48 7.61 8.79 9.78 8.80 10.52
MS22-180629-110319 8.63 7.93 8.90 10.10 9.75 10.45
MS22-180711-112912 8.77 7.97 8.75 10.73 9.76 10.42
MS22-180720-110055 9.06 8.36 8.77 10.86 10.15 9.92
(oscillation_frequency needs 50 spikes, so 50 is the smallest minimum that can be tested)
variant, at least 100 spikes per state (n = 267): {'flfp_pre': 8.5, 'flfp_cooling': 7.76, 'flfp_post': 8.58, 'fp_pre': 10.54, 'fp_cooling': 9.57, 'fp_post': 10.37}
variant, at least 200 spikes per state (n = 121): {'flfp_pre': 8.48, 'flfp_cooling': 7.79, 'flfp_post': 8.6, 'fp_pre': 10.58, 'fp_cooling': 9.55, 'fp_post': 10.34}
Result: reproduced. Over 361 place fields (the paper has 365), fLFP is 8.54, 7.77 and 8.60 Hz and fp 10.55, 9.57 and 10.39 Hz before, during and after cooling, all six within the paper's spreads. Both fall during cooling (fp in 311 fields, P = 3.5e-44; fLFP in 352, P = 1.2e-60) by 9.2% and 9.0%, and both recover (P = 4.5e-37 and 6.8e-61). Each rat shows the same pattern, and every session shows the decrease; fp does not recover in MS22-180720-110055 (9.92 Hz after, 10.15 Hz during cooling). Fields with at least 100 or 200 spikes per state give the same means.
fig4c-fp-flfp-correlation¶
fp and fLFP were correlated across place fields before, during and after cooling (R = 0.61, 0.59 and 0.36; P < 1e-6 each).
In [39]:
exp4c = {"pre": (0.61, 0.061), "cooling": (0.59, 0.059), "post": (0.36, 0.036)}
for st, (v, tol) in exp4c.items():
r, p, n = C.pearson(F3[f"flfp_{st}"], F3[f"fp_{st}"])
print(f"{st}: R = {r:.3f}, P = {p:.2g}, n = {n} (expected {v} +- {tol}); within: {within(r, v, tol)}")
print(" within subjects:", {a: round(C.pearson(g[f"flfp_{st}"], g[f"fp_{st}"])[0], 2) for a, g in F3.groupby("subject")})
pre: R = 0.226, P = 1.5e-05, n = 361 (expected 0.61 +- 0.061); within: False
within subjects: {'MS12': -0.05, 'MS13': 0.41, 'MS21': 0.35, 'MS22': 0.38}
cooling: R = 0.444, P = 6.7e-19, n = 361 (expected 0.59 +- 0.059); within: False
within subjects: {'MS12': 0.13, 'MS13': 0.19, 'MS21': 0.43, 'MS22': 0.5}
post: R = 0.157, P = 0.0029, n = 361 (expected 0.36 +- 0.036); within: False
within subjects: {'MS12': 0.23, 'MS13': 0.64, 'MS21': 0.1, 'MS22': 0.13}
Result: partially reproduced. fp and fLFP correlate positively across fields before (R = 0.226, P = 1.5e-05), during (R = 0.444) and after cooling (R = 0.157, P = 0.0029), weaker than 0.61, 0.59 and 0.36.
fig4ac-fp-flfp-difference-preserved¶
During MS cooling fp and fLFP decreased in parallel, so the difference fp - fLFP remained constant.
In [40]:
dpre, dcool, dpost = F3.fp_pre - F3.flfp_pre, F3.fp_cooling - F3.flfp_cooling, F3.fp_post - F3.flfp_post
p, n = wsr(dcool, dpre)
print(f"fp - fLFP: pre {dpre.mean():.3f}, cooling {dcool.mean():.3f}, post {dpost.mean():.3f} Hz (means); median change cooling minus pre {np.median(dcool - dpre):.3f} Hz; "
f"signed-rank P = {p:.3g} (n = {n})")
print(f"post against cooling: P = {wsr(dpost, dcool)[0]:.2g}")
print("By subject, mean change:", (dcool - dpre).groupby(F3.subject).mean().round(3).to_dict())
r_pre = (F3.fp_pre / F3.flfp_pre).mean(); r_cool = (F3.fp_cooling / F3.flfp_cooling).mean()
print(f"variant, ratio fp / fLFP: pre {r_pre:.3f}, cooling {r_cool:.3f}, P = {wsr(F3.fp_cooling / F3.flfp_cooling, F3.fp_pre / F3.flfp_pre)[0]:.2g}")
fp - fLFP: pre 2.009, cooling 1.806, post 1.790 Hz (means); median change cooling minus pre -0.128 Hz; signed-rank P = 0.000132 (n = 361)
post against cooling: P = 0.98
By subject, mean change: {'MS12': -1.032, 'MS13': 0.012, 'MS21': -0.186, 'MS22': 0.103}
variant, ratio fp / fLFP: pre 1.237, cooling 1.232, P = 0.84
Result: not reproduced. fp - fLFP falls from 2.009 Hz before cooling to 1.806 Hz during it (median change -0.128 Hz, P = 0.000132), so the difference is not preserved by the paper's default test, and it stays at 1.790 Hz after cooling (P = 0.98 against cooling). The ratio fp/fLFP is preserved (1.237 and 1.232, P = 0.84). The change differs between rats (MS12 -1.03 Hz, MS13 and MS22 near 0).
fig4a-frequency-temperature-slopes¶
fLFP and fp measured within place fields decreased linearly with MS temperature (fLFP R = 0.66, slope 0.071 Hz per C; fp R = 0.38, slope 0.09 Hz per C).
In [41]:
W = pd.concat([units[k]["windows"].assign(session=k) for k in UNIT if len(units[k]["windows"])], ignore_index=True)
print("20-trial windows with more than 50 in-field spikes:", len(W), "from", W.groupby(["session", "unit", "field"]).ngroups, "fields")
out4a = {}
for m in ["flfp", "fp"]:
g = W.dropna(subset=[m, "temp"])
r, p, n = C.pearson(g.temp, g[m])
slope = np.polyfit(g.temp, g[m], 1)[0]
out4a[m] = (r, slope, p)
print(f"{m}: R = {r:.3f}, slope {slope:.4f} Hz/C, P = {p:.2g}, windows = {n}")
print("Expected fLFP R 0.66 +- 0.066 and slope 0.071 +- 0.0071; fp R 0.38 +- 0.038 and slope 0.09 +- 0.009")
print("within:", within(out4a["flfp"][0], 0.66, 0.066), within(out4a["flfp"][1], 0.071, 0.0071), within(out4a["fp"][0], 0.38, 0.038), within(out4a["fp"][1], 0.09, 0.009))
print("Temperature range of windows:", round(W.temp.min(), 1), "to", round(W.temp.max(), 1), "C")
print("By subject (R, slope):", {a: {m: (round(C.pearson(g.dropna(subset=[m]).temp, g.dropna(subset=[m])[m])[0], 3), round(np.polyfit(g.dropna(subset=[m]).temp, g.dropna(subset=[m])[m], 1)[0], 4)) for m in ["flfp", "fp"]}
for a, g in W.groupby(W.session.map(subj))})
gm = W.groupby(["session", "unit", "field"])
wc = W.assign(temp_c=W.temp - gm.temp.transform("mean"), flfp_c=W.flfp - gm.flfp.transform("mean"), fp_c=W.fp - gm.fp.transform("mean"))
for m in ["flfp", "fp"]:
g = wc.dropna(subset=[m + "_c"])
print(f"variant, within-field (each field's mean removed), {m}: R = {C.pearson(g.temp_c, g[m + '_c'])[0]:.3f}, slope {np.polyfit(g.temp_c, g[m + '_c'], 1)[0]:.4f}")
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5))
for a, m in zip(ax, ["flfp", "fp"]):
g = W.dropna(subset=[m]); a.scatter(g.temp, g[m], s=1, alpha=0.2); a.set_xlabel("temperature at spikes (C)"); a.set_ylabel(m + " (Hz)")
plt.tight_layout(); plt.show()
print("Plotted ranges: fLFP", round(W.flfp.min(), 2), "to", round(W.flfp.max(), 2), "Hz; fp", round(W.fp.min(), 2), "to", round(W.fp.max(), 2), "Hz")
20-trial windows with more than 50 in-field spikes: 52474 from 563 fields flfp: R = 0.678, slope 0.0716 Hz/C, P = 0, windows = 52474
fp: R = 0.338, slope 0.0859 Hz/C, P = 0, windows = 48355 Expected fLFP R 0.66 +- 0.066 and slope 0.071 +- 0.0071; fp R 0.38 +- 0.038 and slope 0.09 +- 0.009 within: True True False True Temperature range of windows: 17.9 to 37.6 C
By subject (R, slope): {'MS12': {'flfp': (0.58, np.float64(0.0584)), 'fp': (0.38, np.float64(0.1012))}, 'MS13': {'flfp': (0.798, np.float64(0.0748)), 'fp': (0.228, np.float64(0.0766))}, 'MS21': {'flfp': (0.732, np.float64(0.0665)), 'fp': (0.326, np.float64(0.076))}, 'MS22': {'flfp': (0.782, np.float64(0.0785)), 'fp': (0.292, np.float64(0.0758))}}
variant, within-field (each field's mean removed), flfp: R = 0.847, slope 0.0697
variant, within-field (each field's mean removed), fp: R = 0.428, slope 0.0783
Plotted ranges: fLFP 6.51 to 9.67 Hz; fp 6.71 to 16.13 Hz
Result: partially reproduced. Over 52474 windows of 20 trials, fLFP rises with temperature at R = 0.678 and 0.0716 Hz per C, and fp at a slope of 0.0859 Hz per C, all three within tolerance; the fp correlation is R = 0.338, just below 0.38 +- 0.038. The windows overlap (stepped by one trial), so the P values over windows overstate the evidence.
fig3c-field-size-correlated¶
Place field size was highly correlated before and during MS cooling (R = 0.91).
Result: reproduced. Field sizes before and during cooling correlate at R = 0.8197 over 635 fields, at the edge of 0.91 +- 0.091 (0.819); with at least 20 or 50 spikes per state R is 0.8348 and 0.8308, and with at least 100 it is 0.8136, outside. Sizes do not differ between the states (62.7 and 64.9 cm, P = 0.084).
fig3d-size-fp-negative¶
Place field size and fp were negatively correlated before, during and after cooling (R = -0.49, -0.40, -0.40; P < 1e-9).
In [43]:
exp3d = {"pre": (-0.49, 0.049), "cooling": (-0.40, 0.04), "post": (-0.40, 0.04)}
for st, (v, tol) in exp3d.items():
r, p, n = C.pearson(F3[f"size_{st}"], F3[f"fp_{st}"])
print(f"{st}: R = {r:.3f}, P = {p:.2g}, n = {n} (expected {v} +- {tol}); within: {within(r, v, tol)}")
print(" by subject:", {a: round(C.pearson(g[f"size_{st}"], g[f"fp_{st}"])[0], 2) for a, g in F3.groupby("subject")})
for st in exp3d:
print(f"variant, Spearman, {st}: rho = {stats.spearmanr(F3[f'size_{st}'], F3[f'fp_{st}'])[0]:.3f}")
pre: R = -0.311, P = 1.6e-09, n = 361 (expected -0.49 +- 0.049); within: False
by subject: {'MS12': -0.27, 'MS13': -0.51, 'MS21': -0.3, 'MS22': -0.3}
cooling: R = -0.307, P = 2.4e-09, n = 361 (expected -0.4 +- 0.04); within: False
by subject: {'MS12': -0.17, 'MS13': -0.09, 'MS21': -0.34, 'MS22': -0.32}
post: R = -0.323, P = 3.3e-10, n = 361 (expected -0.4 +- 0.04); within: False
by subject: {'MS12': -0.17, 'MS13': -0.6, 'MS21': -0.32, 'MS22': -0.34}
variant, Spearman, pre: rho = -0.267
variant, Spearman, cooling: rho = -0.279 variant, Spearman, post: rho = -0.293
Result: partially reproduced. Field size and fp correlate negatively before, during and after cooling (R = -0.311, -0.307 and -0.323, every P below 3e-09), weaker than -0.49 and -0.40.
fig3e-model-field-size¶
Measured place field size correlated with the size predicted by SPF = v / (fp - fLFP) before (R = 0.55) and during (R = 0.54) cooling, and weakly after (R = 0.11).
In [44]:
exp3e = {"pre": (0.55, 0.055), "cooling": (0.54, 0.054), "post": (0.11, 0.011)}
for st, (v, tol) in exp3e.items():
dfq = F3[f"fp_{st}"] - F3[f"flfp_{st}"]
ok = dfq > 0.2
pred = F3.loc[ok, f"v_{st}"] / dfq[ok]
r, p, n = C.pearson(F3.loc[ok, f"size_{st}"], pred)
print(f"{st}: fields {n} (of {len(F3)}), R = {r:.3f}, P = {p:.2g} (expected {v} +- {tol}); within: {within(r, v, tol)}; "
f"median predicted {pred.median():.1f} cm, measured {F3.loc[ok, f'size_{st}'].median():.1f} cm")
rs = stats.spearmanr(F3.loc[ok, f"size_{st}"], pred)[0]
print(f" variant, Spearman rho = {rs:.3f}; by subject R:", {a: round(C.pearson(F3.loc[ok & (F3.subject == a), f"size_{st}"], pred[F3.loc[ok, 'subject'] == a])[0], 2) for a in sorted(F3.subject.unique())})
pre: fields 339 (of 361), R = 0.246, P = 4.7e-06 (expected 0.55 +- 0.055); within: False; median predicted 46.7 cm, measured 65.0 cm
variant, Spearman rho = 0.240; by subject R: {'MS12': 0.23, 'MS13': 0.49, 'MS21': 0.15, 'MS22': 0.35}
cooling: fields 346 (of 361), R = 0.344, P = 4.9e-11 (expected 0.54 +- 0.054); within: False; median predicted 48.4 cm, measured 65.0 cm
variant, Spearman rho = 0.325; by subject R: {'MS12': 0.46, 'MS13': -0.23, 'MS21': 0.25, 'MS22': 0.43}
post: fields 339 (of 361), R = 0.325, P = 9.1e-10 (expected 0.11 +- 0.011); within: False; median predicted 46.7 cm, measured 65.0 cm
variant, Spearman rho = 0.343; by subject R: {'MS12': 0.16, 'MS13': -0.08, 'MS21': 0.47, 'MS22': 0.23}
Result: partially reproduced. Measured sizes correlate with v / (fp - fLFP) before (R = 0.246) and during cooling (R = 0.344), both with P below 1e-5 but below 0.55 and 0.54, and the post value (R = 0.325) is far above the paper's 0.11.
fig3f-information-unchanged¶
The information content of place cell spikes was not altered by MS cooling (P > 0.26).
In [45]:
pcells = F[["session", "unit"]].drop_duplicates()
pi = UC.merge(pcells, on=["session", "unit"])
print("Place cells:", len(pi))
for a, b in [("pre", "cooling"), ("cooling", "post"), ("pre", "post")]:
p, n = wsr(pi[f"info_{b}"], pi[f"info_{a}"])
print(f"{a} against {b}: medians {pi['info_' + a].median():.3f} and {pi['info_' + b].median():.3f} bits/spike; median change {np.nanmedian(pi['info_' + b] - pi['info_' + a]):.3f}; "
f"signed-rank P = {p:.2g} (n = {n})")
pe = pi[pi.n_eq >= 20]
print(f"\nvariant, each state's spikes subsampled to the smallest count of the three (20 draws averaged), cells with at least 20 spikes in each state (n = {len(pe)}, median {pe.n_eq.median():.0f} spikes):")
for a, b in [("pre", "cooling"), ("cooling", "post"), ("pre", "post")]:
p, n = wsr(pe[f"info_eq_{b}"], pe[f"info_eq_{a}"])
p0, _ = wsr(pe[f"info_{b}"], pe[f"info_{a}"])
print(f" {a} against {b}: medians {pe['info_eq_' + a].median():.3f} and {pe['info_eq_' + b].median():.3f}; median change {np.nanmedian(pe['info_eq_' + b] - pe['info_eq_' + a]):.3f}; "
f"P = {p:.2g} (n = {n}); without subsampling on the same cells P = {p0:.2g}")
print("By subject, median pre and cooling:", pi.groupby("subject")[["info_pre", "info_cooling"]].median().round(3).to_dict("index"))
Place cells: 374
pre against cooling: medians 1.253 and 1.230 bits/spike; median change -0.029; signed-rank P = 0.014 (n = 374)
cooling against post: medians 1.230 and 1.193 bits/spike; median change -0.006; signed-rank P = 0.31 (n = 372)
pre against post: medians 1.253 and 1.193 bits/spike; median change -0.021; signed-rank P = 0.0064 (n = 372)
variant, each state's spikes subsampled to the smallest count of the three (20 draws averaged), cells with at least 20 spikes in each state (n = 360, median 251 spikes):
pre against cooling: medians 1.247 and 1.227; median change -0.032; P = 0.0073 (n = 360); without subsampling on the same cells P = 0.002
cooling against post: medians 1.227 and 1.205; median change -0.008; P = 0.29 (n = 360); without subsampling on the same cells P = 0.13
pre against post: medians 1.247 and 1.205; median change -0.018; P = 0.0092 (n = 360); without subsampling on the same cells P = 0.00079
By subject, median pre and cooling: {'MS12': {'info_pre': 1.117, 'info_cooling': 0.905}, 'MS13': {'info_pre': 1.434, 'info_cooling': 1.394}, 'MS21': {'info_pre': 1.257, 'info_cooling': 1.246}, 'MS22': {'info_pre': 1.223, 'info_cooling': 1.216}}
Result: not reproduced. The paper's default paired test finds information content lower during cooling than before (median change -0.029 against a median of 1.253 bits/spike, P = 0.014 over 374 place cells) and lower after cooling than before (P = 0.0064); cooling and post do not differ (P = 0.31). With each state's spikes subsampled to the same count (360 cells), pre against cooling still gives P = 0.0073 and pre against post P = 0.0092, so the decrease does not come from the different numbers of spikes in the states.
fig3hi-population-vector-unchanged¶
The spatial scale of the population vector correlation of pyramidal cells was not affected by MS cooling.
In [46]:
pvt = pd.DataFrame({k: {"pre": units[k]["pv"]["pre"], "cooling": units[k]["pv"]["cooling"], "post": units[k]["pv"]["post"],
"pyramidal": int(((UC.session == k) & (UC.cell_type == "pyramidal")).sum())} for k in UNIT}).T
print("Distance (cm) at which the population vector correlation drops below 0.5:"); print(pvt.round(2).to_string())
p, n = wsr(pvt.cooling.astype(float), pvt.pre.astype(float))
print(f"cooling against pre: median change {np.nanmedian(pvt.cooling.astype(float) - pvt.pre.astype(float)):.2f} cm, signed-rank P = {p:.3g} (n = {n})")
print(f"cooling against post: P = {wsr(pvt.cooling.astype(float), pvt.post.astype(float))[0]:.3g}")
Distance (cm) at which the population vector correlation drops below 0.5:
pre cooling post pyramidal
MS12-170716-172307 22.24 22.92 22.45 29.0
MS12-170719-095305 25.66 29.97 25.43 28.0
MS13-171128-113924 21.72 24.91 21.53 33.0
MS13-171201-130527 27.03 18.30 21.25 35.0
MS21-180625-153927 27.62 24.13 28.42 57.0
MS21-180628-155921 23.54 23.07 24.22 86.0
MS21-180629-110332 23.14 23.39 22.49 113.0
MS22-180628-120341 36.88 31.54 27.92 38.0
MS22-180629-110319 25.50 26.80 25.83 92.0
MS22-180711-112912 25.54 25.18 30.32 39.0
MS22-180720-110055 27.42 28.27 27.44 55.0
cooling against pre: median change 0.25 cm, signed-rank P = 0.966 (n = 11)
cooling against post: P = 0.898
Result: reproduced. The distance at which the population vector correlation falls below 0.5 does not change with cooling (median change 0.25 cm, P = 0.966 over 11 sessions; 18.3 to 36.9 cm across sessions and states).
fig4d-precession-slope-decrease¶
The phase precession slope decreased slightly but significantly during MS cooling (P < 0.007 and P < 0.0003).
In [47]:
fsl = F3.dropna(subset=["slope_pre", "slope_cooling", "slope_post"])
print(f"Fields with a slope in every state: {len(fsl)} of {len(F3)}; slopes in rad/cm")
print("Median slope pre, cooling, post:", fsl[["slope_pre", "slope_cooling", "slope_post"]].median().round(4).to_dict(),
"; fraction negative in pre:", round(float((fsl.slope_pre < 0).mean()), 3))
for name, f_ in [("(a) magnitude", np.abs), ("(b) signed", lambda x: x)]:
p1, n1 = wsr(f_(fsl.slope_cooling), f_(fsl.slope_pre))
p2, n2 = wsr(f_(fsl.slope_cooling), f_(fsl.slope_post))
d1 = np.median(f_(fsl.slope_cooling) - f_(fsl.slope_pre)); d2 = np.median(f_(fsl.slope_cooling) - f_(fsl.slope_post))
print(f"{name}: cooling minus pre median {d1:.4f}, P = {p1:.2g}; cooling minus post median {d2:.4f}, P = {p2:.2g} (n = {n1})")
neg = fsl[(fsl.slope_pre < 0) & (fsl.slope_cooling < 0) & (fsl.slope_post < 0)]
print(f"variant, fields precessing (negative slope) in all three states (n = {len(neg)}): magnitude cooling minus pre median {np.median(neg.slope_cooling.abs() - neg.slope_pre.abs()):.4f}, "
f"P = {wsr(neg.slope_cooling.abs(), neg.slope_pre.abs())[0]:.2g}; against post P = {wsr(neg.slope_cooling.abs(), neg.slope_post.abs())[0]:.2g}")
print("Mean resultant length of the fits, median:", fsl[["rbar_pre", "rbar_cooling", "rbar_post"]].median().round(3).to_dict())
Fields with a slope in every state: 361 of 361; slopes in rad/cm
Median slope pre, cooling, post: {'slope_pre': -0.0874, 'slope_cooling': -0.0874, 'slope_post': -0.0889} ; fraction negative in pre: 0.875
(a) magnitude: cooling minus pre median 0.0000, P = 0.83; cooling minus post median -0.0020, P = 0.022 (n = 361)
(b) signed: cooling minus pre median -0.0004, P = 0.32; cooling minus post median 0.0019, P = 0.15 (n = 361)
variant, fields precessing (negative slope) in all three states (n = 286): magnitude cooling minus pre median 0.0002, P = 0.82; against post P = 0.08
Mean resultant length of the fits, median: {'rbar_pre': 0.403, 'rbar_cooling': 0.419, 'rbar_post': 0.406}
Result: partially reproduced. Over 361 fields, the magnitude of the precession slope during cooling does not differ from pre (median change 0, P = 0.83) and is smaller than post (median -0.002 rad/cm, P = 0.022), so under reading (a) one comparison holds; the signed slope does not differ in either (P = 0.32 and 0.15).
fig4gh-fp-speed¶
fp depended on running speed more strongly than fLFP, and fp during cooling was lower at all velocities.
In [48]:
rng_4gh = np.random.default_rng(4)
bins = list(range(20, 120, 20))
long = []
for _, r_ in F.iterrows():
for st in ["pre", "cooling", "post"]:
for b in bins:
fp_, fl_ = r_[f"fp_{st}_v{b}"], r_[f"flfp_{st}_v{b}"]
if np.isfinite(fp_) and np.isfinite(fl_):
long.append({"field": (r_.session, r_.unit, r_.field), "state": st, "speed": b + 10, "fp": fp_, "flfp": fl_})
L = pd.DataFrame(long)
print("Field, state and speed-bin values:", len(L), "from", L.field.nunique(), "fields; per bin:", L.groupby("speed").size().to_dict())
s_fp = stats.linregress(L.speed, L.fp); s_fl = stats.linregress(L.speed, L.flfp)
print(f"fp slope {s_fp.slope:.4f} Hz per cm/s (P = {s_fp.pvalue:.2g}); fLFP slope {s_fl.slope:.4f} (P = {s_fl.pvalue:.2g})")
codes = pd.factorize(L.field)[0]
rows_by = [np.flatnonzero(codes == i) for i in range(codes.max() + 1)]
boot = []
for _ in range(1000):
pick = rng_4gh.integers(0, len(rows_by), len(rows_by))
ii = np.concatenate([rows_by[i] for i in pick])
x = L.speed.values[ii]
boot.append(np.polyfit(x, L.fp.values[ii], 1)[0] - np.polyfit(x, L.flfp.values[ii], 1)[0])
lo, hi = np.percentile(boot, [2.5, 97.5])
print(f"difference of slopes (fp minus fLFP) {s_fp.slope - s_fl.slope:.4f}, 95% bootstrap interval over fields {lo:.4f} to {hi:.4f}")
print("\nMedian fp by speed bin, fields measured in both pre and cooling:")
part2 = []
for b in bins:
pr = L[(L.state == "pre") & (L.speed == b + 10)].set_index("field").fp
co = L[(L.state == "cooling") & (L.speed == b + 10)].set_index("field").fp
both = pr.index.intersection(co.index)
if len(both) >= 10:
part2.append(co[both].median() < pr[both].median())
print(f" {b}-{b + 20} cm/s: n = {len(both)}, pre {pr[both].median():.2f}, cooling {co[both].median():.2f} Hz, P = {wsr(co[both], pr[both])[0]:.2g}")
else:
print(f" {b}-{b + 20} cm/s: n = {len(both)} (fewer than 10)")
print("fp lower during cooling in every bin with at least 10 fields:", all(part2), f"({sum(part2)} of {len(part2)})")
print("Mean fp and fLFP by speed bin (all states):", L.groupby("speed")[["fp", "flfp"]].mean().round(2).to_dict("index"))
Field, state and speed-bin values: 1420 from 474 fields; per bin: {30: 21, 50: 67, 70: 274, 90: 502, 110: 556}
fp slope 0.0136 Hz per cm/s (P = 3.2e-12); fLFP slope 0.0061 (P = 1.5e-15)
difference of slopes (fp minus fLFP) 0.0074, 95% bootstrap interval over fields 0.0040 to 0.0112
Median fp by speed bin, fields measured in both pre and cooling:
20-40 cm/s: n = 0 (fewer than 10)
40-60 cm/s: n = 3 (fewer than 10)
60-80 cm/s: n = 38, pre 9.95, cooling 8.93 Hz, P = 4.8e-06
80-100 cm/s: n = 65, pre 10.20, cooling 9.26 Hz, P = 8.6e-06
100-120 cm/s: n = 110, pre 10.31, cooling 9.52 Hz, P = 5.7e-09
fp lower during cooling in every bin with at least 10 fields: True (3 of 3)
Mean fp and fLFP by speed bin (all states): {30: {'fp': 9.44, 'flfp': 7.98}, 50: {'fp': 9.52, 'flfp': 7.99}, 70: {'fp': 9.68, 'flfp': 8.12}, 90: {'fp': 10.06, 'flfp': 8.23}, 110: {'fp': 10.29, 'flfp': 8.38}}
Result: reproduced. fp rises with running speed at 0.0136 Hz per cm/s (P = 3.2e-12), more than fLFP at 0.0061 (difference 0.0074, bootstrap interval 0.0040 to 0.0112), and in the three speed bins with at least 10 fields fp is lower during cooling (medians 8.93 to 9.52 against 9.95 to 10.31 Hz, P below 1e-5 each). Below 60 cm/s too few fields could be measured.
fig5d-ccg-period¶
The theta-repeated cross-correlograms of place cell pairs had a period of about 85 ms before cooling, which lengthened during MS cooling (P = 0.002), while the period of the phase cross-correlograms was unaffected (P = 0.32).
In [49]:
P5 = pd.concat([units[k]["pairs"].assign(session=k) for k in FIG5 if len(units[k]["pairs"])], ignore_index=True)
print("FIG5 sessions:", len(FIG5), "; pairs of overlapping fields:", len(P5), "; by session:", P5.groupby("session").size().to_dict())
print("Pairs with a main peak (time) in pre, cooling, post:", {c: int(P5[f"tau_{c}"].notna().sum()) for c in ["pre", "cooling", "post"]},
"; with at least 50 spikes of each cell:", {c: int(((P5[f"na_{c}"] >= 50) & (P5[f"nb_{c}"] >= 50)).sum()) for c in ["pre", "cooling", "post"]})
rows = []
for k, g in P5.groupby("session"):
ctl_t = pd.concat([g.period_pre, g.period_post]); ctl_p = pd.concat([g.phase_period_pre, g.phase_period_post])
rows.append({"session": k, "control_ms": ctl_t.median(), "cooling_ms": g.period_cooling.median(), "pre_ms": g.period_pre.median(), "post_ms": g.period_post.median(),
"control_rad_pi": ctl_p.median() / np.pi, "cooling_rad_pi": g.phase_period_cooling.median() / np.pi,
"n_control": int(ctl_t.notna().sum()), "n_cooling": int(g.period_cooling.notna().sum())})
per = pd.DataFrame(rows).set_index("session")
print(per.round(2).to_string())
print(f"\nMedian control period over sessions: {per.control_ms.median():.1f} ms (expected 85 +- 8.5); cooling {per.cooling_ms.median():.1f} ms")
p_t, n_t = wsr(per.cooling_ms, per.control_ms); p_p, n_p = wsr(per.cooling_rad_pi, per.control_rad_pi)
print(f"cooling longer than control in {(per.cooling_ms > per.control_ms).sum()} of {n_t} sessions, signed-rank P = {p_t:.3g}")
print(f"phase period: control {per.control_rad_pi.median():.3f} pi, cooling {per.cooling_rad_pi.median():.3f} pi rad; longer in {(per.cooling_rad_pi > per.control_rad_pi).sum()} of {n_p}, P = {p_p:.3g}")
print("LFP theta period in the same sessions (1000 / mean trial theta), control and cooling (ms):",
round(float(np.median([1000 / ((sm.loc[k, 'theta_freq_pre'] + sm.loc[k, 'theta_freq_post']) / 2) for k in FIG5])), 1),
round(float(np.median([1000 / sm.loc[k, 'theta_freq_cooling'] for k in FIG5])), 1))
print("Pooled pairs, median period: pre", round(P5.period_pre.median(), 1), ", cooling", round(P5.period_cooling.median(), 1), ", post", round(P5.period_post.median(), 1), "ms")
# per-pair periods, and the nearest neighbouring peak in place of the largest
bins_p = np.arange(60, 210, 10)
for c in ["pre", "cooling", "post"]:
v = P5[f"period_{c}"].dropna(); vn = P5[f"period_near_{c}"].dropna()
print(f"per-pair period, {c}: histogram over {bins_p[0]} to {bins_p[-1]} ms in 10 ms bins (largest neighbour) {np.histogram(v, bins_p)[0].tolist()}; "
f"(nearest neighbour) {np.histogram(vn, bins_p)[0].tolist()}")
both = P5.dropna(subset=["period_pre", "period_near_pre"])
print(f"pre pairs where the largest neighbour is more than 1.25 times as far as the nearest: {(both.period_pre > 1.25 * both.period_near_pre).mean():.3f} of {len(both)}")
print("|Ltr| of pairs (cm), quartiles:", P5.ltr.abs().quantile([.25, .5, .75]).round(1).tolist())
for lim in [10, 20]:
g = P5[P5.ltr.abs() < lim]
print(f"pairs with |Ltr| < {lim} cm (n = {len(g)}): median control period {pd.concat([g.period_pre, g.period_post]).median():.1f} ms (largest), "
f"{pd.concat([g.period_near_pre, g.period_near_post]).median():.1f} ms (nearest); cooling {g.period_cooling.median():.1f} and {g.period_near_cooling.median():.1f} ms")
pn_ = P5.groupby("session").apply(lambda g: pd.Series({"control": pd.concat([g.period_near_pre, g.period_near_post]).median(), "cooling": g.period_near_cooling.median()}))
p_n, n_n = wsr(pn_.cooling, pn_.control)
print(f"variant, nearest neighbour: median control period over sessions {pn_.control.median():.1f} ms, cooling {pn_.cooling.median():.1f} ms; "
f"longer in {(pn_.cooling > pn_.control).sum()} of {n_n}, P = {p_n:.3g}")
vn_pre = P5.period_near_pre.dropna()
print(f"variant, nearest neighbour: pre pairs with a period below 70 ms {(vn_pre < 70).sum()} of {len(vn_pre)} (the search starts at 60 ms)")
pnp_ = P5.groupby("session").apply(lambda g: pd.Series({"control": pd.concat([g.phase_period_near_pre, g.phase_period_near_post]).median() / np.pi,
"cooling": g.phase_period_near_cooling.median() / np.pi}))
p_np, n_np = wsr(pnp_.cooling, pnp_.control)
print(f"variant, nearest neighbour, phase period: control {pnp_.control.median():.3f} pi, cooling {pnp_.cooling.median():.3f} pi rad; "
f"longer in {(pnp_.cooling > pnp_.control).sum()} of {n_np}, P = {p_np:.3g}")
fig, ax = plt.subplots(figsize=(6, 3))
ax.hist([pd.concat([P5.period_pre, P5.period_post]).dropna(), pd.concat([P5.period_near_pre, P5.period_near_post]).dropna()], bins=bins_p, label=["largest", "nearest"])
ax.set_xlabel("pair period, control (ms)"); ax.set_ylabel("pairs"); ax.legend(); plt.tight_layout(); plt.show()
print("Plotted: all control pair periods; by construction they lie in 60 to 200 ms")
# variant: period of the summed correlogram of all pairs, by spectral peak and by the fit
from scipy.ndimage import gaussian_filter1d
for k in FIG5:
cs = units[k]["ccg_sum"]
ctl = cs["pre"]["time"] + cs["post"]["time"]
sp_ = []
for c in [ctl, cs["cooling"]["time"]]:
res = c - gaussian_filter1d(c, 40, mode="nearest")
sel = np.abs(C.TIME_LAGS) <= 300
z = res[sel] * np.hanning(sel.sum()); spec = np.abs(np.fft.rfft(z, 1 << 16)); fr = np.fft.rfftfreq(1 << 16, 0.001)
band = (fr >= 5) & (fr <= 14); sp_.append(1000 / fr[band][np.argmax(spec[band])])
print(f"variant, summed correlogram spectral period {k}: control {sp_[0]:.1f}, cooling {sp_[1]:.1f} ms")
FIG5 sessions: 9 ; pairs of overlapping fields: 4845 ; by session: {'MS12-170716-172307': 286, 'MS12-170719-095305': 234, 'MS13-171128-113924': 85, 'MS21-180625-153927': 379, 'MS21-180628-155921': 1212, 'MS21-180629-110332': 956, 'MS22-180628-120341': 86, 'MS22-180629-110319': 1265, 'MS22-180711-112912': 342}
Pairs with a main peak (time) in pre, cooling, post: {'pre': 3502, 'cooling': 3863, 'post': 4100} ; with at least 50 spikes of each cell: {'pre': 3510, 'cooling': 3877, 'post': 4114}
control_ms cooling_ms pre_ms post_ms control_rad_pi cooling_rad_pi n_control n_cooling
session
MS12-170716-172307 109.50 117.75 111.00 108.00 1.74 1.77 471 264
MS12-170719-095305 114.50 119.00 113.50 115.00 1.76 1.76 367 215
MS13-171128-113924 116.50 112.00 105.00 123.75 1.73 1.84 133 61
MS21-180625-153927 107.00 112.50 105.75 108.00 1.71 1.69 595 289
MS21-180628-155921 112.50 114.50 115.50 110.00 1.72 1.72 2100 887
MS21-180629-110332 109.00 114.00 109.00 109.00 1.67 1.69 1329 804
MS22-180628-120341 109.25 115.00 109.00 109.50 1.79 1.78 144 72
MS22-180629-110319 107.00 112.00 106.50 107.00 1.76 1.74 1984 1050
MS22-180711-112912 108.00 112.50 107.00 108.50 1.74 1.79 476 221
Median control period over sessions: 109.2 ms (expected 85 +- 8.5); cooling 114.0 ms
cooling longer than control in 8 of 9 sessions, signed-rank P = 0.0195
phase period: control 1.740 pi, cooling 1.760 pi rad; longer in 4 of 9, P = 0.59
LFP theta period in the same sessions (1000 / mean trial theta), control and cooling (ms): 115.7 129.5
Pooled pairs, median period: pre 110.0 , cooling 114.0 , post 109.0 ms
per-pair period, pre: histogram over 60 to 200 ms in 10 ms bins (largest neighbour) [14, 71, 397, 641, 608, 391, 291, 365, 342, 191, 73, 62, 37, 16]; (nearest neighbour) [211, 522, 947, 980, 556, 199, 51, 20, 9, 4, 0, 0, 0, 0]
per-pair period, cooling: histogram over 60 to 200 ms in 10 ms bins (largest neighbour) [14, 43, 153, 461, 868, 846, 525, 357, 265, 180, 80, 36, 31, 4]; (nearest neighbour) [265, 462, 738, 830, 823, 507, 175, 50, 9, 2, 0, 1, 1, 0]
per-pair period, post: histogram over 60 to 200 ms in 10 ms bins (largest neighbour) [12, 69, 399, 877, 728, 438, 310, 417, 431, 202, 83, 57, 53, 24]; (nearest neighbour) [227, 506, 975, 1382, 704, 224, 65, 12, 0, 3, 1, 1, 0, 0]
pre pairs where the largest neighbour is more than 1.25 times as far as the nearest: 0.477 of 3499
|Ltr| of pairs (cm), quartiles: [15.0, 30.0, 50.0]
pairs with |Ltr| < 10 cm (n = 716): median control period 98.5 ms (largest), 91.0 ms (nearest); cooling 107.5 and 96.0 ms
pairs with |Ltr| < 20 cm (n = 1622): median control period 99.0 ms (largest), 91.0 ms (nearest); cooling 108.0 and 97.5 ms
variant, nearest neighbour: median control period over sessions 91.0 ms, cooling 95.5 ms; longer in 7 of 9, P = 0.0312
variant, nearest neighbour: pre pairs with a period below 70 ms 211 of 3499 (the search starts at 60 ms)
variant, nearest neighbour, phase period: control 1.285 pi, cooling 1.285 pi rad; longer in 5 of 9, P = 0.367
Plotted: all control pair periods; by construction they lie in 60 to 200 ms variant, summed correlogram spectral period MS12-170716-172307: control 119.8, cooling 118.7 ms variant, summed correlogram spectral period MS12-170719-095305: control 111.6, cooling 122.3 ms variant, summed correlogram spectral period MS13-171128-113924: control 108.3, cooling 100.4 ms variant, summed correlogram spectral period MS21-180625-153927: control 112.8, cooling 143.1 ms variant, summed correlogram spectral period MS21-180628-155921: control 145.6, cooling 166.8 ms variant, summed correlogram spectral period MS21-180629-110332: control 114.6, cooling 128.8 ms variant, summed correlogram spectral period MS22-180628-120341: control 114.2, cooling 99.0 ms variant, summed correlogram spectral period MS22-180629-110319: control 127.3, cooling 149.6 ms variant, summed correlogram spectral period MS22-180711-112912: control 199.8, cooling 123.0 ms
Result: partially reproduced. Under the criterion's estimator the period of the pair cross-correlograms lengthens during cooling, from a median 109.2 to 114.0 ms (longer in 8 of 9 sessions, P = 0.0195), and the phase period does not change (1.740 and 1.760 pi rad, P = 0.59), but the control period is 109.2 ms, not 85 +- 8.5 ms. That estimator takes the largest neighbouring peak within 60 to 200 ms, and in 47.7% of pre pairs this peak is more than 1.25 times as far as the nearest one, so it often picks the second neighbour. Taking the nearest neighbouring peak (a variant) gives a control period of 91.0 ms, within tolerance, which lengthens to 95.5 ms during cooling (7 of 9 sessions, P = 0.031); pairs with |Ltr| below 10 cm give 91.0 ms. Under that variant the phase period does not change either (1.285 pi rad in both, P = 0.367). The nearest-peak estimator is biased the other way: 211 of 3499 pre pairs have a period of 60 to 70 ms, below any plausible theta period, at the floor of the 60 ms search window, so its 91.0 ms is not a clean test of 85 ms either. Both control periods are shorter than the LFP theta period (115.7 ms), as the paper says of its 85 ms.
fig5cg-compression¶
Distance-time compression of place cell pairs (0.54 cm/ms before cooling) decreased during MS cooling (P = 0.02), while distance-phase compression was unchanged (P = 0.92).
In [50]:
def slope_xy(x, y):
ok = np.isfinite(x) & np.isfinite(y)
if ok.sum() < 5:
return np.nan
return float(np.polyfit(x[ok], y[ok], 1)[0])
c_pre = slope_xy(P5.tau_pre.values, P5.ltr.values)
print(f"Compression over pooled pre pairs (slope of Ltr on tau): {c_pre:.3f} cm/ms (expected 0.54 +- 0.054), pairs = {int(P5.tau_pre.notna().sum())}")
print(f" variant, inverse slope of tau on Ltr: {1 / slope_xy(P5.ltr.values, P5.tau_pre.values):.3f} cm/ms")
print(f" variant, half-period window: {slope_xy(P5.tau_half_pre.values, P5.ltr.values):.3f} cm/ms")
bt = P5.dropna(subset=["tau_pre"])
print(f" pre pairs whose lag differs between the 70 ms window and the half-period window: {(bt.tau_pre != bt.tau_half_pre).mean():.3f} of {len(bt)} "
f"(no peak in the half-period window: {bt.tau_half_pre.isna().mean():.3f})")
r_, p_, n_ = C.pearson(P5.tau_pre, P5.ltr); print(f" correlation of tau and Ltr in pre: R = {r_:.3f}, P = {p_:.2g}")
rows = []
for k, g in P5.groupby("session"):
t_ctl = np.r_[g.tau_pre.values, g.tau_post.values]; l_ctl = np.r_[g.ltr.values, g.ltr.values]
ph_ctl = np.r_[g.phase_pre.values, g.phase_post.values]
rows.append({"session": k, "c_control": slope_xy(t_ctl, l_ctl), "c_cooling": slope_xy(g.tau_cooling.values, g.ltr.values),
"c_pre": slope_xy(g.tau_pre.values, g.ltr.values),
"cphase_control": slope_xy(ph_ctl, l_ctl), "cphase_cooling": slope_xy(g.phase_cooling.values, g.ltr.values),
"inv_control": 1 / slope_xy(l_ctl, t_ctl), "inv_cooling": 1 / slope_xy(g.ltr.values, g.tau_cooling.values),
"half_control": slope_xy(np.r_[g.tau_half_pre.values, g.tau_half_post.values], l_ctl), "half_cooling": slope_xy(g.tau_half_cooling.values, g.ltr.values)})
cmp_ = pd.DataFrame(rows).set_index("session")
print(cmp_.round(3).to_string())
for a, b, name in [("c_cooling", "c_control", "time compression (cm/ms)"), ("cphase_cooling", "cphase_control", "phase compression (cm/rad)"),
("inv_cooling", "inv_control", "variant, inverse slope (cm/ms)"), ("half_cooling", "half_control", "variant, half-period window (cm/ms)")]:
p, n = wsr(cmp_[a], cmp_[b])
print(f"{name}: control median {cmp_[b].median():.3f}, cooling {cmp_[a].median():.3f}; cooling lower in {(cmp_[a] < cmp_[b]).sum()} of {n}, signed-rank P = {p:.3g}")
Compression over pooled pre pairs (slope of Ltr on tau): 0.273 cm/ms (expected 0.54 +- 0.054), pairs = 3502
variant, inverse slope of tau on Ltr: 7.175 cm/ms
variant, half-period window: 0.156 cm/ms
pre pairs whose lag differs between the 70 ms window and the half-period window: 0.105 of 3502 (no peak in the half-period window: 0.003)
correlation of tau and Ltr in pre: R = 0.195, P = 2.1e-31
c_control c_cooling c_pre cphase_control cphase_cooling inv_control inv_cooling half_control half_cooling
session
MS12-170716-172307 0.115 0.066 0.022 2.588 2.214 15.942 20.565 0.025 0.050
MS12-170719-095305 0.255 0.191 0.296 5.469 0.841 9.203 12.262 0.250 0.079
MS13-171128-113924 0.352 0.078 0.315 4.588 4.524 2.627 9.801 0.053 -0.099
MS21-180625-153927 0.478 0.274 0.427 6.459 5.365 6.510 10.420 0.340 0.259
MS21-180628-155921 0.233 0.062 0.203 3.347 3.389 9.216 31.586 0.117 0.028
MS21-180629-110332 0.298 0.204 0.338 6.090 6.078 5.079 7.921 0.162 0.118
MS22-180628-120341 0.314 0.458 0.529 5.622 7.475 9.085 4.930 0.361 0.483
MS22-180629-110319 0.323 0.274 0.370 4.943 5.312 5.362 5.641 0.133 0.130
MS22-180711-112912 0.165 0.150 -0.020 2.312 1.991 13.209 15.311 0.073 -0.018
time compression (cm/ms): control median 0.298, cooling 0.191; cooling lower in 8 of 9, signed-rank P = 0.0547
phase compression (cm/rad): control median 4.943, cooling 4.524; cooling lower in 6 of 9, signed-rank P = 0.426
variant, inverse slope (cm/ms): control median 9.085, cooling 10.420; cooling lower in 1 of 9, signed-rank P = 0.0547
variant, half-period window (cm/ms): control median 0.133, cooling 0.079; cooling lower in 7 of 9, signed-rank P = 0.129
Result: not reproduced. The distance-time compression over pooled pre pairs is 0.273 cm/ms, half of 0.54, and the per-session compression is lower during cooling in 8 of 9 sessions without reaching significance (0.298 against 0.191 cm/ms, P = 0.0547). The phase compression does not change (P = 0.426). Time lag and distance correlate weakly in these pairs (R = 0.195), so the slope depends on the fit: the inverse regression over pooled pre pairs gives 7.175 cm/ms. In 10.5% of pre pairs the lag changes when the window is narrowed to half the theta period, and that window gives 0.156 cm/ms, so the compression magnitude depends on how the lag is picked and cannot be judged from this method. The direction of the change also depends on the fit: the inverse slope is higher during cooling (9.085 against 10.420 cm/ms, lower in 1 of 9 sessions).
fig5h-time-lags-increase¶
During MS cooling the theta time lags between place cell pairs increased compared with before and after cooling (P = 8e-5 and P = 0.00032), and did not differ between before and after (P = 0.79).
In [51]:
for a, b in [("cooling", "pre"), ("cooling", "post"), ("post", "pre")]:
x, y = P5[f"tau_{a}"].abs(), P5[f"tau_{b}"].abs()
ok = x.notna() & y.notna()
p, n = wsr(x[ok], y[ok])
print(f"(a) pairs, |tau| {a} against {b}: medians {x[ok].median():.1f} and {y[ok].median():.1f} ms, mean difference {(x[ok] - y[ok]).mean():.2f} ms, "
f"{a} longer in {(x[ok] > y[ok]).sum()}, shorter in {(x[ok] < y[ok]).sum()}; signed-rank P = {p:.3g} (n = {n})")
ses = P5.assign(**{f"abs_{c}": P5[f"tau_{c}"].abs() for c in ["pre", "cooling", "post"]}).groupby("session")[["abs_pre", "abs_cooling", "abs_post"]].median()
print(ses.round(1).to_string())
for a, b in [("cooling", "pre"), ("cooling", "post"), ("post", "pre")]:
p, n = wsr(ses[f"abs_{a}"], ses[f"abs_{b}"])
print(f"(b) sessions, median |tau| {a} against {b}: longer in {(ses['abs_' + a] > ses['abs_' + b]).sum()} of {n}, P = {p:.3g}")
for a, b in [("cooling", "pre"), ("cooling", "post")]:
x, y = P5[f"tau_half_{a}"].abs(), P5[f"tau_half_{b}"].abs(); ok = x.notna() & y.notna()
print(f"variant, half-period window, {a} against {b}: medians {x[ok].median():.1f} and {y[ok].median():.1f}, P = {wsr(x[ok], y[ok])[0]:.3g} (n = {int(ok.sum())})")
PU = pd.concat([units[k]["pairs"].assign(session=k) for k in UNIT if len(units[k]["pairs"])], ignore_index=True)
x, y = PU.tau_cooling.abs(), PU.tau_pre.abs(); ok = x.notna() & y.notna()
print(f"variant, all UNIT sessions ({PU.session.nunique()}), cooling against pre: medians {x[ok].median():.1f} and {y[ok].median():.1f}, P = {wsr(x[ok], y[ok])[0]:.3g} (n = {int(ok.sum())})")
fig, ax = plt.subplots(1, 2, figsize=(10, 3.2))
for c, col in [("pre", "r"), ("cooling", "b"), ("post", "g")]:
v = np.sort(P5[f"tau_{c}"].abs().dropna()); ax[0].plot(v, np.linspace(0, 1, len(v)), col, label=c)
v = np.sort(P5[f"phase_{c}"].abs().dropna()); ax[1].plot(v, np.linspace(0, 1, len(v)), col)
ax[0].set_xlabel("|time lag| (ms)"); ax[1].set_xlabel("|phase lag| (rad)"); ax[0].legend(); ax[0].set_ylabel("cumulative fraction")
plt.tight_layout(); plt.show()
print("Plotted ranges: |tau| 0 to", P5[["tau_pre", "tau_cooling", "tau_post"]].abs().max().max(), "ms; |phase| 0 to", round(P5[["phase_pre", "phase_cooling", "phase_post"]].abs().max().max(), 3), "rad")
(a) pairs, |tau| cooling against pre: medians 30.5 and 28.5 ms, mean difference 1.51 ms, cooling longer in 1682, shorter in 1366; signed-rank P = 1.14e-06 (n = 3136)
(a) pairs, |tau| cooling against post: medians 30.5 and 27.5 ms, mean difference 1.97 ms, cooling longer in 1912, shorter in 1545; signed-rank P = 1.3e-10 (n = 3571)
(a) pairs, |tau| post against pre: medians 27.5 and 28.5 ms, mean difference -0.71 ms, post longer in 1498, shorter in 1562; signed-rank P = 0.0967 (n = 3139)
abs_pre abs_cooling abs_post
session
MS12-170716-172307 26.5 36.0 27.5
MS12-170719-095305 27.5 27.5 29.5
MS13-171128-113924 30.5 34.5 31.0
MS21-180625-153927 29.5 30.5 24.5
MS21-180628-155921 30.5 32.5 29.5
MS21-180629-110332 29.5 28.5 26.5
MS22-180628-120341 20.5 30.5 21.5
MS22-180629-110319 25.5 28.5 25.5
MS22-180711-112912 28.5 30.5 28.5
(b) sessions, median |tau| cooling against pre: longer in 7 of 9, P = 0.0234
(b) sessions, median |tau| cooling against post: longer in 8 of 9, P = 0.0156
(b) sessions, median |tau| post against pre: longer in 4 of 9, P = 0.812
variant, half-period window, cooling against pre: medians 29.5 and 26.5, P = 3.84e-22 (n = 3122)
variant, half-period window, cooling against post: medians 28.5 and 24.5, P = 1.97e-29 (n = 3548)
variant, all UNIT sessions (11), cooling against pre: medians 29.5 and 28.5, P = 6.63e-05 (n = 3521)
Plotted ranges: |tau| 0 to 68.5 ms; |phase| 0 to 3.094 rad
Result: reproduced. Over 3136 pairs, time lags are longer during cooling than before (median 30.5 against 28.5 ms, P = 1.14e-06) and after (P = 1.3e-10), and do not differ between before and after (P = 0.0967); per session (reading b) the same holds (P = 0.023, 0.016 and 0.81).
fig5i-phase-lags-unchanged¶
Theta phase lags between place cell pairs were not affected by MS cooling (P > 0.24 between all conditions).
In [52]:
for a, b in [("cooling", "pre"), ("cooling", "post"), ("post", "pre")]:
x, y = P5[f"phase_{a}"].abs(), P5[f"phase_{b}"].abs()
ok = x.notna() & y.notna()
p, n = wsr(x[ok], y[ok])
print(f"|phase lag| {a} against {b}: medians {x[ok].median():.3f} and {y[ok].median():.3f} rad, mean difference {(x[ok] - y[ok]).mean():.4f} rad; signed-rank P = {p:.3g} (n = {n})")
ses_p = P5.assign(**{f"abs_{c}": P5[f"phase_{c}"].abs() for c in ["pre", "cooling", "post"]}).groupby("session")[["abs_pre", "abs_cooling", "abs_post"]].median()
print("variant, session medians:", {f"{a} vs {b}": round(wsr(ses_p[f'abs_{a}'], ses_p[f'abs_{b}'])[0], 3) for a, b in [("cooling", "pre"), ("cooling", "post"), ("post", "pre")]})# variant: phase lag picked within +-pi/2 in place of +-pi (the phase correlogram period is about 1.74 pi)
for c in ["pre", "cooling", "post"]:
ok = P5[f"phase_{c}"].notna()
same = np.isclose(P5.loc[ok, f"phase_{c}"], P5.loc[ok, f"phase_half_{c}"])
print(f"{c}: pairs whose phase lag changes with the +-pi/2 window {1 - same.mean():.3f} of {int(ok.sum())} (none in it: {int(P5.loc[ok, f'phase_half_{c}'].isna().sum())})")
for a, b in [("cooling", "pre"), ("cooling", "post"), ("post", "pre")]:
x, y = P5[f"phase_half_{a}"].abs(), P5[f"phase_half_{b}"].abs(); ok = x.notna() & y.notna()
print(f"variant, +-pi/2 window, |phase lag| {a} against {b}: medians {x[ok].median():.3f} and {y[ok].median():.3f} rad, "
f"mean difference {(x[ok] - y[ok]).mean():.4f} rad; signed-rank P = {wsr(x[ok], y[ok])[0]:.3g} (n = {int(ok.sum())})")
|phase lag| cooling against pre: medians 1.367 and 1.304 rad, mean difference 0.0509 rad; signed-rank P = 0.00151 (n = 3137)
|phase lag| cooling against post: medians 1.335 and 1.272 rad, mean difference 0.0561 rad; signed-rank P = 0.00433 (n = 3572)
|phase lag| post against pre: medians 1.272 and 1.335 rad, mean difference -0.0202 rad; signed-rank P = 0.506 (n = 3138)
variant, session medians: {'cooling vs pre': 0.062, 'cooling vs post': 0.004, 'post vs pre': 0.961}
pre: pairs whose phase lag changes with the +-pi/2 window 0.433 of 3499 (none in it: 293)
cooling: pairs whose phase lag changes with the +-pi/2 window 0.442 of 3865 (none in it: 351)
post: pairs whose phase lag changes with the +-pi/2 window 0.418 of 4102 (none in it: 348)
variant, +-pi/2 window, |phase lag| cooling against pre: medians 0.675 and 0.644 rad, mean difference 0.0251 rad; signed-rank P = 0.0196 (n = 2679)
variant, +-pi/2 window, |phase lag| cooling against post: medians 0.675 and 0.675 rad, mean difference 0.0159 rad; signed-rank P = 0.162 (n = 3025)
variant, +-pi/2 window, |phase lag| post against pre: medians 0.675 and 0.644 rad, mean difference 0.0150 rad; signed-rank P = 0.179 (n = 2668)
Result: not reproduced. Phase lags are slightly larger during cooling than before (median 1.367 against 1.304 rad, mean difference 0.051 rad, P = 0.0015) and after (P = 0.0043), and before and after do not differ (P = 0.51). The change is small, but by the paper's test the phase lags are not preserved. The +-pi window spans about 1.15 periods of the phase correlogram, so it can hold two candidate peaks. With the lag picked within +-pi/2 (a variant), 43.3% of pre pairs get another lag or none (293 have no peak in the window); cooling is still above pre (mean difference 0.025 rad, P = 0.0196) but not above post (P = 0.162), and pre and post do not differ (P = 0.178). Part of the increase may come from how the lag is picked; the increase over pre remains.
fig5j-fraction-active-unchanged¶
The fraction of place cells active in a theta cycle was not significantly affected by MS cooling (P = 0.2).
In [53]:
CY = pd.concat([units[k]["cycles"].assign(session=k) for k in FIG5], ignore_index=True)
print("Theta cycles during runs in trials:", CY.groupby("state").size().to_dict(), "; place cells per session:", {k: units[k]["n_place_cells"] for k in FIG5})
print("Mean fraction active:", CY.groupby("state").frac_active.mean().round(4).to_dict(), "; median:", CY.groupby("state").frac_active.median().round(4).to_dict())
for a, b in [("pre", "cooling"), ("cooling", "post"), ("pre", "post")]:
p, na, nb = rks(CY.frac_active[CY.state == a], CY.frac_active[CY.state == b])
print(f"rank-sum {a} against {b}: P = {p:.3g} (n = {na}, {nb})")
sesj = CY.groupby(["session", "state"]).frac_active.mean().unstack()
print(sesj.round(4).to_string())
print("variant, session means, signed-rank:", {f"{a} vs {b}": round(wsr(sesj[a], sesj[b])[0], 3) for a, b in [("pre", "cooling"), ("cooling", "post"), ("pre", "post")]},
"; cooling lower than pre in", int((sesj.cooling < sesj.pre).sum()), "of", len(sesj))
print("Mean cycle duration (ms) by state:", CY.groupby("state").duration_ms.mean().round(1).to_dict())
print("Cycles per session and state:"); print(CY.groupby(["session", "state"]).size().unstack().to_string())
Theta cycles during runs in trials: {'cooling': 10152, 'post': 11669, 'pre': 9037} ; place cells per session: {'MS12-170716-172307': 22, 'MS12-170719-095305': 23, 'MS13-171128-113924': 19, 'MS21-180625-153927': 37, 'MS21-180628-155921': 59, 'MS21-180629-110332': 60, 'MS22-180628-120341': 12, 'MS22-180629-110319': 61, 'MS22-180711-112912': 27}
Mean fraction active: {'cooling': 0.2405, 'post': 0.2444, 'pre': 0.226} ; median: {'cooling': 0.2333, 'post': 0.2432, 'pre': 0.2174}
rank-sum pre against cooling: P = 1.53e-22 (n = 9037, 10152)
rank-sum cooling against post: P = 0.00186 (n = 10152, 11669)
rank-sum pre against post: P = 8.12e-44 (n = 9037, 11669)
state cooling post pre
session
MS12-170716-172307 0.3151 0.2842 0.2802
MS12-170719-095305 0.2788 0.2716 0.2276
MS13-171128-113924 0.1589 0.1754 0.1751
MS21-180625-153927 0.2133 0.2467 0.2385
MS21-180628-155921 0.2426 0.2543 0.2456
MS21-180629-110332 0.1710 0.1785 0.1779
MS22-180628-120341 0.2640 0.2972 0.2762
MS22-180629-110319 0.2254 0.2332 0.1990
MS22-180711-112912 0.2323 0.2756 0.2276
variant, session means, signed-rank: {'pre vs cooling': 0.734, 'cooling vs post': 0.074, 'pre vs post': 0.004} ; cooling lower than pre in 5 of 9
Mean cycle duration (ms) by state: {'cooling': 131.1, 'post': 117.2, 'pre': 117.8}
Cycles per session and state:
state cooling post pre
session
MS12-170716-172307 1616 1054 1162
MS12-170719-095305 1599 1215 882
MS13-171128-113924 968 1555 1602
MS21-180625-153927 1046 1247 708
MS21-180628-155921 633 1772 970
MS21-180629-110332 1213 1367 803
MS22-180628-120341 1279 1609 1035
MS22-180629-110319 1042 973 939
MS22-180711-112912 756 877 936
Result: not reproduced. Over 9037 pre and 10152 cooling theta cycles, the fraction of place cells active per cycle is higher during cooling (mean 0.2405 against 0.226, rank-sum P = 1.5e-22). Per session the cooling mean is lower than pre in 5 of 9 sessions (P = 0.734), and cooling does not differ from post (P = 0.074), while pre and post differ (P = 0.004). The pooled test mixes sessions whose means range from 0.159 to 0.315 and that contribute different numbers of cycles in each state, so it reflects the make-up of the pool as much as an effect. Cycles are longer during cooling (131.1 ms against 117.8 and 117.2 ms).
fig5k-spikes-per-cycle-unchanged¶
The number of spikes emitted by all neurons per theta cycle was preserved during MS cooling (P > 0.058 between all conditions).
In [54]:
print("Mean spikes per cycle (all good units):", CY.groupby("state").spikes.mean().round(2).to_dict())
for a, b in [("pre", "cooling"), ("cooling", "post"), ("pre", "post")]:
p, na, nb = rks(CY.spikes[CY.state == a], CY.spikes[CY.state == b])
print(f"rank-sum {a} against {b}: P = {p:.3g} (n = {na}, {nb})")
sesk = CY.groupby(["session", "state"]).spikes.mean().unstack()
print(sesk.round(2).to_string())
print("variant, session means, signed-rank:", {f"{a} vs {b}": round(wsr(sesk[a], sesk[b])[0], 3) for a, b in [("pre", "cooling"), ("cooling", "post"), ("pre", "post")]})
norm = CY.groupby("session").spikes.transform("mean")
print("variant, spikes per cycle over each session's mean, by state:", (CY.spikes / norm).groupby(CY.state).mean().round(3).to_dict())
Mean spikes per cycle (all good units): {'cooling': 44.98, 'post': 51.68, 'pre': 45.68}
rank-sum pre against cooling: P = 0.0758 (n = 9037, 10152)
rank-sum cooling against post: P = 3.64e-53 (n = 10152, 11669)
rank-sum pre against post: P = 9.54e-53 (n = 9037, 11669)
state cooling post pre
session
MS12-170716-172307 21.85 19.07 17.66
MS12-170719-095305 23.86 23.11 19.33
MS13-171128-113924 23.08 21.27 21.66
MS21-180625-153927 43.32 48.35 48.66
MS21-180628-155921 84.55 88.25 86.49
MS21-180629-110332 67.96 74.49 70.55
MS22-180628-120341 30.99 35.91 29.78
MS22-180629-110319 90.74 91.37 84.91
MS22-180711-112912 59.97 64.51 58.78
variant, session means, signed-rank: {'pre vs cooling': 0.57, 'cooling vs post': 0.129, 'pre vs post': 0.02}
variant, spikes per cycle over each session's mean, by state: {'cooling': 1.01, 'post': 1.027, 'pre': 0.953}
Result: not reproduced. Spikes per theta cycle do not differ between pre and cooling (45.68 and 44.98, P = 0.0758) but are higher after cooling (51.68) than during it or before it (both P below 1e-50), so not all comparisons are null. Per session, cooling differs from neither pre (P = 0.57) nor post (P = 0.129), and post differs from pre (P = 0.02). Session means range from 17.66 to 91.37 spikes, so the pooled tests depend on how many cycles each session contributes.
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 |
|---|---|
| session-structure | reproduced |
| cooling-time-course | partially reproduced (minimum after 121 s, not 60 s) |
| fig2g-theta-freq-temperature | reproduced (R = 0.789) |
| fig2h-theta-power-temperature | not reproduced (R = -0.056) |
| fig2j-frequency-consistent-across-animals | reproduced |
| fig2i-theta-pre-cooling-post | partially reproduced (frequency yes, power only against pre) |
| disc-theta-frequency-12-percent | partially reproduced (8.9 to 9.2%) |
| fig6d-errors-increase | partially reproduced (2.3-fold; MS13 and MS22 not significant) |
| fig6e-errors-recover | reproduced |
| fig6c-worst-trials-4-to-8 | reproduced (trial 7) |
| fig6f-errors-temperature | partially reproduced (R = -0.84, stronger) |
| fig6g-errors-theta-frequency | partially reproduced (R = -0.87, stronger) |
| fig6h-errors-theta-power | partially reproduced (1 rat, not 2) |
| fig6i-errors-not-speed | not reproduced (3 of 5 rats correlate) |
| fig4e-speed-decreases-within-session | reproduced |
| fig5-speed-control-cooling | not reproduced (cooling 74.9 to 76.0 cm/s, below the spread) |
| suppl2-place-cell-rates-decrease | partially reproduced (7.7% and 6.7%) |
| disc-interneuron-rates-decrease | partially reproduced (6.5 to 9.1%) |
| suppl2d-interneuron-oscillation-slows | reproduced |
| suppl1g-waveforms-unchanged | not reproduced (small changes, drift not excluded) |
| fig4b-fp-flfp-decrease | reproduced (all six means) |
| fig4c-fp-flfp-correlation | partially reproduced (weaker R) |
| fig4ac-fp-flfp-difference-preserved | not reproduced (difference falls 0.13 Hz; ratio preserved) |
| fig4a-frequency-temperature-slopes | partially reproduced (3 of 4 values) |
| fig3c-field-size-correlated | reproduced (R = 0.8197, at the edge) |
| fig3d-size-fp-negative | partially reproduced (weaker R) |
| fig3e-model-field-size | partially reproduced (weaker R) |
| fig3f-information-unchanged | not reproduced (small decrease, P = 0.014) |
| fig3hi-population-vector-unchanged | reproduced |
| fig4d-precession-slope-decrease | partially reproduced (one of two comparisons) |
| fig4gh-fp-speed | reproduced |
| fig5d-ccg-period | partially reproduced (lengthens; 109 ms, not 85; 91 ms with the nearest peak) |
| fig5cg-compression | not reproduced (0.27 cm/ms; decrease P = 0.055) |
| fig5h-time-lags-increase | reproduced |
| fig5i-phase-lags-unchanged | not reproduced (small increase, P = 0.0015) |
| fig5j-fraction-active-unchanged | not reproduced (pooled cycles differ) |
| fig5k-spikes-per-cycle-unchanged | not reproduced (post differs) |
| suppl1-hippocampal-temperature | blocked: data missing |
How the analysis differs from the paper's. 45 of the paper's 53 sessions have cooling states here (six are not deposited, one lacks states, one lacks trials); temperature analyses use the 37 with a recorded temperature. Cell types are not in the files and were made here from waveform width and rate. Place fields use the automatic criteria the Methods give for Figure 5, for all figures, where the paper defined fields by hand for Figures 3 and 4. Pair lags are the peak of each pair's correlogram and compression a regression, where the paper fits a surface to the density of all pairs. The CCG period is measured per pair, with an estimator that often picks the second neighbouring peak (the nearest-peak variant is shown). Units come from 11 sessions in 4 rats (the paper's Figure 5 has 18 sessions in 5 rats); MS10 has no curated units. Two unit sessions were left out because spikes carry no spatial information there.
What was found about the dandiset. The temperature of the nine MS10 sessions and of the two control sessions of 2018-07-19 is not a recording (the paper's code lists no temperature channel for them); it runs from 17 to 37 C (exactly in nine of ten), or is constant. The two controls also carry cooling labels although the code describes them as sessions without cooling, and their theta does not slow as in cooled sessions. Six of the paper's sessions are missing, one session has no cooling states and one no trials. No file holds cell types, place fields or electrode locations, and the units of seven sessions are all 'unsorted'. The trial state 'Cooling off' is not described. In two sessions spikes carry no spatial information. The first published version has no trials table.