Reproducing Valero et al. (2022) from dandiset 000568¶
Paper. Valero M, Zutshi I, Yoon E, Buzsaki G (2022). Probing subthreshold dynamics of hippocampal neurons by pulsed optogenetics. Science 375, 570-574. doi:10.1126/science.abm1891. This attempt used the full text of the NIH author manuscript in PubMed Central (PMC9632609): the main report and its figure legends. The Materials and Methods and figures S1 to S15 are in the supplement, which could not be obtained (PMC's copy is behind a captcha and the publisher's is not open access). Claims about supplementary figures are taken from the main text's statements of them, and method details come from the code named below.
Dandiset. DANDI:000568, version 0.230705.1633, streamed. The paper cites no DANDI version; its data statement points to the Buzsaki Lab Databank, whose page lists this dandiset for the paper (its other copy is on Globus, which needs an account and has no pinned identifier, so it was not used). This is the only published version. The draft holds the same 138 assets (paths, sizes and asset ids, printed below), so nothing the analysis uses differs.
What the dandiset holds. 28 sessions from 4 mice, each deposited twice (a full file and a no-raw-data file with the same units). Units carry spike times only, with no cell type. There are epochs, ripple intervals and linearized position. Laps and uLED pulse times are in 2 sessions only (fCamk1 sess9 and sess12), so every claim that needs light responses or place fields is tested on these two sessions of one mouse.
Code read. The paper names buzcode (github.com/buzsakilab/buzcode, commit e615663): the lab's general toolbox, with place-field and ripple functions but no script for this paper's figures. Its field rule (bz_findPlaceFields1D) is ported here with its defaults. The first author's later lab code, HippoCookBook (github.com/valegarman/HippoCookBook, commit 0cbded3), is not cited by the paper. It holds the light-response test (getuLEDResponse.m) used here. Whether either matches the paper's Methods could not be checked, as the Methods were not available.
Outside inputs. None.
| Claim | Figure | Tested on |
|---|---|---|
| cohort-counts | Results | all 28 sessions |
| fig1c-pulse-protocol | 1C | 2 sessions with pulse times |
| fig1e-responsive-fraction | 1E | 2 sessions |
| fig1i-swr-decrease | 1F-I | 2 sessions |
| figs10-swr-gain-vs-baseline-rate | S10 | 24 sessions with ripples |
| fig2b-theta-trough-gain | 2A-B | 2 sessions |
| figs11-place-cell-fraction | S11 | 2 sessions |
| fig2e-infield-gain | 2C-E | 2 sessions |
| figs12-gain-rate-correlations | S12 | 2 sessions |
| fig3a-place-cells-more-excitable | 3A | 2 sessions |
| fig3c-unmasked-fields | 3B-C | 2 sessions |
| fig3e-ghost-fields | 3E | 2 sessions |
| fig3f-light-stability | 3F | 2 sessions |
| fig3g-spatial-information | 3G | 2 sessions |
| fig4c-swr-cofiring-place-pairs | 4B-C | 2 sessions |
| fig4d-light-unmasked-cofiring | 4C-D | 2 sessions |
| fig1j-light-intensity | 1J, 2F, 3D | blocked: data missing (no intensity levels) |
| fig1km-intracellular | 1K-M | blocked: data missing |
| figs9-excitatory-vs-inhibitory-swr | S9 | blocked: metadata insufficient (no cell types) |
| fig3h-position-decoding | 3H | blocked: methods underspecified |
| fig4ef-assemblies | 4E-F | blocked: methods underspecified |
| fig4gh-replay | 4G-H | blocked: methods underspecified |
The run takes about seven minutes, most of it streaming the files.
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import valero as v
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
def ci95(k, n):
"""Clopper-Pearson 95% interval of a fraction k/n, in percent."""
r = stats.binomtest(int(k), int(n)).proportion_ci()
return f"{100 * r.low:.1f} to {100 * r.high:.1f}%"
started = time.time()
assets = v.assets()
small = [a for a in assets if "no-raw-data" in a["path"]]
full = {v.session_key(a["path"]): a for a in assets if "no-raw-data" not in a["path"]}
sessions = v.load_sessions([a["url"] for a in small])
for a, s in zip(small, sessions):
s["key"] = v.session_key(a["path"])
s["full_url"] = full[s["key"]]["url"] if s["key"] in full else None
from dandi.dandiapi import DandiAPIClient
with DandiAPIClient() as client:
listing = {ver: {(a.path, a.size, a.identifier) for a in client.get_dandiset(v.DANDISET_ID, ver).get_assets()}
for ver in (v.DANDISET_VERSION, "draft")}
published = [x.identifier for x in client.get_dandiset(v.DANDISET_ID).get_versions() if x.identifier != "draft"]
print("Published versions:", published, "; assets in", v.DANDISET_VERSION, len(listing[v.DANDISET_VERSION]),
"; in the draft", len(listing["draft"]), "; same paths, sizes and asset ids:", listing[v.DANDISET_VERSION] == listing["draft"])
print(f"{len(assets)} NWB files in dandiset {v.DANDISET_ID} version {v.DANDISET_VERSION}; "
f"loaded {len(sessions)} no-raw-data files in {time.time() - started:.0f} s")
Published versions: ['0.230705.1633'] ; assets in 0.230705.1633 138 ; in the draft 138 ; same paths, sizes and asset ids: True 56 NWB files in dandiset 000568 version 0.230705.1633; loaded 28 no-raw-data files in 73 s
Data check¶
Each check of the metadata against the values is stated under its output.
# Every file, paired by session; a pair holds the same units, spikes, ripples and position
# when these fingerprints agree.
from concurrent.futures import ThreadPoolExecutor
with ThreadPoolExecutor(16) as ex:
prints = pd.DataFrame(list(ex.map(v.file_fingerprint, assets)))
prints["kind"] = np.where(prints.path.str.contains("no-raw-data"), "no_raw", "full")
prints["size_GB"] = [round(a["size"] / 1e9, 3) for a in assets]
cols = ["n_units", "n_spikes", "spike_index_sum", "spike_head_sum", "spike_tail_sum", "n_ripples",
"ripple_sum", "n_position", "position_sum"]
# Files are paired by the session_id they hold (without "_no_raw_data"), not by their names.
pairs = prints.pivot(index="pair_key", columns="kind",
values=cols + ["size_GB", "has_raw", "has_lfp", "raw_s", "lfp_s", "last_epoch_end_s"])
same = pd.Series({k: all(np.isclose(pairs.loc[k, (c, "full")], pairs.loc[k, (c, "no_raw")]) for c in cols)
for k in pairs.index}, name="identical_content")
print("Files per kind:", prints.kind.value_counts().to_dict())
print("Sessions (by session_id):", prints.pair_key.nunique(), "; sessions with both files:",
int(pairs[("n_units", "full")].notna().sum()))
print("Sessions whose two files hold the same units, spikes, ripples and position:", int(same.sum()), "of", len(same))
print("Subjects:", sorted(prints.subject.unique()))
print("Units table columns (distinct over all 56 files):", prints.units_columns.unique().tolist())
print("session_start_time values:", prints.session_start_time.unique().tolist())
print("\nFiles whose name lacks the session:",
prints[~prints.path.str.contains("ses-")][["path", "session_id"]].values.tolist())
table = pd.DataFrame({"units": pairs[("n_units", "no_raw")], "spikes": pairs[("n_spikes", "no_raw")],
"ripples": pairs[("n_ripples", "no_raw")], "full_GB": pairs[("size_GB", "full")],
"session_s": pairs[("last_epoch_end_s", "full")].round(0),
"full_raw_s": pairs[("raw_s", "full")].round(0), "full_lfp_s": pairs[("lfp_s", "full")].round(0),
"identical": same})
print(table.to_string())
short = table[(table.full_raw_s.notna() & (table.full_raw_s < 0.95 * table.session_s)) |
(table.full_lfp_s < 0.95 * table.session_s)]
print("\nFull files whose raw or LFP signal covers less than 95% of the session:", short.index.tolist())
print("Full files without a raw signal:", table[table.full_raw_s.isna()].index.tolist())
Files per kind: {'no_raw': 28, 'full': 28}
Sessions (by session_id): 28 ; sessions with both files: 28
Sessions whose two files hold the same units, spikes, ripples and position: 28 of 28
Subjects: ['fCamk1', 'fCamk2', 'fCamk3', 'fCamk5']
Units table columns (distinct over all 56 files): ['spike_times,unit_name,group_id,clu_id', 'spike_times,unit_name,clu_id', 'spike_times,unit_name,clu_id,group_id']
session_start_time values: ['2016-03-15T00:00:00-04:00']
Files whose name lacks the session: [['sub-fCamk3/sub-fCamk3_behavior+ecephys+image.nwb', 'fCamk3_201111_sess20']]
units spikes ripples full_GB session_s full_raw_s full_lfp_s identical
fCamk1_200827_sess9 57 2259343 4315 27.061 17920.0 17920.0 17920.0 True
fCamk1_200901_sess12 47 1800735 1796 22.555 14810.0 14810.0 14810.0 True
fCamk1_200902_sess13 37 2230994 7196 1.285 19646.0 NaN 19646.0 True
fCamk1_200904_sess15 59 2512149 6258 1.209 18410.0 NaN 18410.0 True
fCamk1_200908_sess16 59 2647836 8096 1.382 21079.0 NaN 21079.0 True
fCamk1_200909_sess17 58 2663673 4096 1.492 22389.0 NaN 22389.0 True
fCamk1_200910_sess18 54 3831190 13117 2.221 34028.0 NaN 34028.0 True
fCamk1_200911_sess19 38 1614478 10448 1.533 23459.0 NaN 23459.0 True
fCamk2_201012_sess1 64 725574 1110 17.113 11849.0 11849.0 11849.0 True
fCamk2_201013_sess2 30 782147 0 27.904 19760.0 19710.0 19710.0 True
fCamk2_201014_sess3 35 428376 49 21.123 14794.0 14794.0 14794.0 True
fCamk2_201015_sess4 33 690199 0 24.501 12781.0 16736.0 16736.0 True
fCamk3_201028_sess10_cleanned 43 3832688 12881 48.868 31968.0 31968.0 31968.0 True
fCamk3_201029_sess11_cleanned 40 2752937 9756 51.631 33743.0 33743.0 33743.0 True
fCamk3_201030_sess12 47 3997813 11835 54.146 35029.0 35029.0 35029.0 True
fCamk3_201102_sess13 99 4405532 2123 50.443 32645.0 32645.0 32645.0 True
fCamk3_201103_sess14 44 3570580 12910 51.329 32929.0 32929.0 32929.0 True
fCamk3_201105_sess16 48 4327891 13274 53.154 33898.0 33898.0 33898.0 True
fCamk3_201106_sess17 46 2240864 9056 0.028 34621.0 378.0 6850.0 True
fCamk3_201109_sess18 66 2858559 13753 48.746 31725.0 31725.0 31725.0 True
fCamk3_201110_sess19 61 3809039 15205 50.408 32645.0 32645.0 32645.0 True
fCamk3_201111_sess20 63 2899263 5024 48.582 31836.0 31836.0 31836.0 True
fCamk3_201113_sess22 58 2704593 8558 50.284 33039.0 33039.0 33039.0 True
fCamk5_210406_sess10 70 5752737 217 53.185 34976.0 34976.0 34976.0 True
fCamk5_210407_sess11 59 2855703 6093 47.629 31899.0 31899.0 31899.0 True
fCamk5_210408_sess12 62 3699871 6842 42.758 28876.0 28876.0 28876.0 True
fCamk5_210412_sess14 50 2369131 5718 49.972 34070.0 34070.0 34070.0 True
fCamk5_210414_sess16 54 2501927 5023 48.665 33229.0 33229.0 33229.0 True
Full files whose raw or LFP signal covers less than 95% of the session: ['fCamk3_201106_sess17']
Full files without a raw signal: ['fCamk1_200902_sess13', 'fCamk1_200904_sess15', 'fCamk1_200908_sess16', 'fCamk1_200909_sess17', 'fCamk1_200910_sess18', 'fCamk1_200911_sess19']
The dandiset holds 56 NWB files and 82 videos. The NWB files cover 28 sessions: each session is
deposited twice, as a full file and a no-raw-data file. In all 28 the two files hold the same units,
spike times, ripples and position (the fingerprints above agree), so each session is counted once
and read from its no-raw-data file. Units carry spike times, a unit name, a cluster id and, in most
files, the shank (group_id); no file gives a cell type. The full file of fCamk3 sess20 lacks the session in its name;
the session_id inside it identifies it. The full file of fCamk3 sess17 is 28 MB: its raw signal
covers 378 s and its LFP 6850 s of a 34621 s session. Six full files of fCamk1 have the LFP but no
raw signal. In fCamk2 sess4 the signals run 16736 s, past the last epoch at 12781 s.
session_start_time is 2016-03-15 in every file, while the session ids date the recordings to 2020
and 2021, so the start time is wrong.
# Epochs, laps, ripples and pulse series per session; baseline rates of the units.
rows = []
for s in sessions:
e0 = s["epochs"][0]
dur = e0["stop"] - e0["start"]
s["base_rate"] = np.array([((x >= e0["start"]) & (x < e0["stop"])).sum() / dur for x in s["spikes"]])
s["pyr"] = s["base_rate"] < v.PYR_MAX_RATE_HZ
paradigms = [e.get("behavioral_paradigm", "") for e in s["epochs"]]
manip = [e.get("manipulation", "") for e in s["epochs"]]
s["stim_protocol"] = any("uled" in m.lower() for m in manip)
rows.append({"session": s["key"].split("_ses-")[1], "subject": s["subject"], "units": len(s["spikes"]),
"putative_pyr": int(s["pyr"].sum()), "first_epoch_s": round(dur),
"epochs": len(s["epochs"]), "labelled": any(paradigms), "uLED_epochs": s["stim_protocol"],
"chemogenetic": any(("CNO" in p) or ("DMSO" in p) for p in paradigms),
"laps": len(s["laps"]), "ripples": len(s["ripples"]),
"pulse_series": len(s["pulses"]), "position_samples": len(s.get("pos_t", []))})
overview = pd.DataFrame(rows)
print(overview.to_string(index=False))
print("\nSessions per subject:", overview.groupby("subject").size().to_dict())
print("Sessions with pulse times:", overview[overview.pulse_series > 0].session.tolist())
print("Sessions with laps:", overview[overview.laps > 0].session.tolist())
print("Sessions without epoch labels:", overview[~overview.labelled].session.tolist())
print("Sessions with CNO or DMSO epochs:", int(overview.chemogenetic.sum()))
print("\nEpoch names of each kind of session:")
for s in sessions:
if s["key"].endswith(("fCamk1-200827-sess9", "fCamk3-201029-sess11-cleanned", "fCamk2-201012-sess1")):
print(" ", s["key"], [(e.get("behavioral_paradigm") or e.get("epoch_name"), e.get("manipulation", "")) for e in s["epochs"]])
session subject units putative_pyr first_epoch_s epochs labelled uLED_epochs chemogenetic laps ripples pulse_series position_samples
fCamk1-200827-sess9 fCamk1 57 51 5195 5 True True False 111 4315 6 53143
fCamk1-200901-sess12 fCamk1 47 40 4203 5 True True False 132 1796 6 45583
fCamk1-200902-sess13 fCamk1 37 32 7323 5 True True False 0 7196 0 42626
fCamk1-200904-sess15 fCamk1 59 53 931 6 True True False 0 6258 0 65004
fCamk1-200908-sess16 fCamk1 59 54 4794 6 True True False 0 8096 0 62195
fCamk1-200909-sess17 fCamk1 58 54 4636 5 True True False 0 4096 0 64885
fCamk1-200910-sess18 fCamk1 54 50 8239 4 True True False 0 13117 0 80741
fCamk1-200911-sess19 fCamk1 38 36 5375 5 True True False 0 10448 0 98527
fCamk2-201012-sess1 fCamk2 64 60 3795 5 False False False 0 1110 0 0
fCamk2-201013-sess2 fCamk2 30 29 4275 5 False False False 0 0 0 0
fCamk2-201014-sess3 fCamk2 35 35 4566 5 False False False 0 49 0 0
fCamk2-201015-sess4 fCamk2 33 33 3955 4 False False False 0 0 0 0
fCamk3-201028-sess10-cleanned fCamk3 43 37 2739 10 True True True 0 12881 0 186724
fCamk3-201029-sess11-cleanned fCamk3 40 36 2822 10 True True True 0 9756 0 216059
fCamk3-201030-sess12 fCamk3 47 42 2982 10 True True True 0 11835 0 208953
fCamk3-201102-sess13 fCamk3 99 85 2794 10 False False False 0 2123 0 125797
fCamk3-201103-sess14 fCamk3 44 40 2704 11 True True True 0 12910 0 132518
fCamk3-201105-sess16 fCamk3 48 43 2942 10 True True True 0 13274 0 160990
fCamk3-201106-sess17 fCamk3 46 41 2680 10 True True True 0 9056 0 175950
fCamk3-201109-sess18 fCamk3 66 62 2808 10 True True True 0 13753 0 110046
fCamk3-201110-sess19 fCamk3 61 56 2891 10 True True True 0 15205 0 148382
fCamk3-201111-sess20 fCamk3 63 59 3098 10 True True True 0 5024 0 187405
fCamk3-201113-sess22 fCamk3 58 54 2986 10 True True True 0 8558 0 110657
fCamk5-210406-sess10 fCamk5 70 62 3559 10 False False False 0 217 0 183576
fCamk5-210407-sess11 fCamk5 59 55 3078 10 True True True 0 6093 0 84191
fCamk5-210408-sess12 fCamk5 62 55 2771 10 True True True 0 6842 0 80493
fCamk5-210412-sess14 fCamk5 50 46 2928 10 True True True 0 5718 0 73377
fCamk5-210414-sess16 fCamk5 54 51 2921 11 True True True 0 5023 0 80984
Sessions per subject: {'fCamk1': 8, 'fCamk2': 4, 'fCamk3': 11, 'fCamk5': 5}
Sessions with pulse times: ['fCamk1-200827-sess9', 'fCamk1-200901-sess12']
Sessions with laps: ['fCamk1-200827-sess9', 'fCamk1-200901-sess12']
Sessions without epoch labels: ['fCamk2-201012-sess1', 'fCamk2-201013-sess2', 'fCamk2-201014-sess3', 'fCamk2-201015-sess4', 'fCamk3-201102-sess13', 'fCamk5-210406-sess10']
Sessions with CNO or DMSO epochs: 14
Epoch names of each kind of session:
sub-fCamk1_ses-fCamk1-200827-sess9 [('BaselinePre', 'None'), ('PreStim', 'uLED random stim'), ('Maze', 'uLED random stim'), ('PostStim', 'uLED random stim'), ('BaselinePost', 'None')]
sub-fCamk2_ses-fCamk2-201012-sess1 [('fCamk2_201012_090020', ''), ('fCamk2_201012_115808', ''), ('fCamk2_201012_130347', ''), ('fCamk2_201012_134511', ''), ('fCamk2_201012_140000', '')]
sub-fCamk3_ses-fCamk3-201029-sess11-cleanned [('BaselinePre', 'None'), ('PreStim', 'uLED random stim'), ('Maze', 'uLED random stim'), ('PostStim', 'uLED random stim'), ('BaselinePost', 'None'), ('CNOPre', 'None'), ('CNOPreStim', 'uLED random stim'), ('CNOMaze', 'uLED random stim'), ('CNOPostStim', 'uLED random stim'), ('CNOPost', 'None')]
There are four subjects, matching the paper's four mice. Pulse times are deposited for 2 sessions only (fCamk1 sess9 and sess12). The other 20 sessions whose epochs name "uLED random stim" have no pulse times, so their light responses cannot be computed. Only the same two sessions have laps; position is present in 24 sessions. Six sessions (all four of fCamk2, fCamk3 sess13 and fCamk5 sess10) name their epochs only by a time stamp, without paradigm or manipulation. 14 sessions have a second block of epochs named CNO or DMSO, a chemogenetic manipulation that the paper does not describe. fCamk2 has no position and 0 to 1110 ripples per session.
# Baseline firing rates, and the putative pyramidal threshold.
rates = np.concatenate([s["base_rate"] for s in sessions])
print("Units:", len(rates), "; baseline rate percentiles 5/25/50/75/95 (Hz):",
np.percentile(rates, [5, 25, 50, 75, 95]).round(2).tolist(), "; max", rates.max().round(1))
for thr in (2.5, 5, 10, 15):
print(f" below {thr:>4} Hz: {(rates < thr).sum()} units ({(rates < thr).mean():.1%})")
isi_viol = []
for s in sessions:
for x in s["spikes"]:
if len(x) > 100:
isi_viol.append(np.mean(np.diff(x) < 0.002))
isi_viol = np.array(isi_viol)
print("Units whose spike times are not sorted:", sum(int(np.any(np.diff(x) < 0)) for s in sessions for x in s["spikes"]))
print("Fraction of ISIs below 2 ms per unit, percentiles 50/90/99:", np.percentile(isi_viol, [50, 90, 99]).round(4).tolist())
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist(np.log10(rates + 1e-3), bins=60, color="0.4")
ax.axvline(np.log10(v.PYR_MAX_RATE_HZ), color="C3", lw=1)
ax.set_xlabel("log10 baseline rate (Hz)")
ax.set_ylabel("units")
ax.set_title(f"all {len(rates)} units; red: {v.PYR_MAX_RATE_HZ} Hz")
plt.show()
Units: 1481 ; baseline rate percentiles 5/25/50/75/95 (Hz): [0.02, 0.53, 1.05, 1.99, 7.52] ; max 25.5 below 2.5 Hz: 1205 units (81.4%) below 5 Hz: 1351 units (91.2%) below 10 Hz: 1441 units (97.3%) below 15 Hz: 1468 units (99.1%)
Units whose spike times are not sorted: 0 Fraction of ISIs below 2 ms per unit, percentiles 50/90/99: [0.0015, 0.0054, 0.0258]
The files carry no cell type. Baseline-epoch rates have a median of 1.05 Hz. 1351 of 1481 units (91.2%) are below 5 Hz, the threshold used here for putative pyramidal cells; with 2.5 or 10 Hz the count is 1205 or 1441. Spike times are sorted in every unit, and ISIs below 2 ms are rare (median 0.15% of a unit's ISIs, 99th percentile 2.6%).
# The two sessions with pulse times: the pulse series, their stated units and values, and the
# position and laps.
opto = [s for s in sessions if s["pulses"]]
for s in opto:
print(s["key"])
for site, p in s["pulses"].items():
d = (p["off"] - p["on"]) * 1000
print(f" {site}: {len(p['on'])} pulses, stored value(s) {p['values'].tolist()} '{p['unit']}', "
f"duration median {np.median(d):.2f} ms (range {d.min():.2f} to {d.max():.2f}), "
f"first {p['on'][0]:.1f} s, last {p['on'][-1]:.1f} s")
print(" stimulation epochs:", [(e["behavioral_paradigm"], round(e["start"]), round(e["stop"]))
for e in s["epochs"] if "uLED" in e["manipulation"]])
x = s["pos_x"]
print(f" linearized position: {len(x)} samples, {np.nanmin(x):.1f} to {np.nanmax(x):.1f}, stated unit "
f"{s['pos_attrs']['unit']!r} with conversion {s['pos_attrs']['conversion']}")
laps = s["laps"]
mz = v.epoch_interval(s, "Maze")[0]
print(f" laps: {len(laps)}, all inside the Maze epoch: {bool(np.all((laps[:, 0] >= mz[0]) & (laps[:, 1] <= mz[1])))}, "
f"directions {np.unique(laps[:, 2], return_counts=True)[1].tolist()}, duration median {np.median(laps[:, 1] - laps[:, 0]):.1f} s")
on = np.sort(np.concatenate([p["on"] for p in s["pulses"].values()]))
per_lap = np.array([((on >= a) & (on < b)).sum() / (b - a) for a, b, _ in laps])
no_stim = np.flatnonzero(per_lap < 1) + 1
print(f" laps with fewer than 1 pulse/s (numbered from 1): {no_stim.tolist()}")
start_x = np.interp(laps[:, 0], s["pos_t"], s["pos_x"])
stop_x = np.interp(laps[:, 1], s["pos_t"], s["pos_x"])
for d in (0, 1):
k = laps[:, 2] == d
print(f" direction {d}: position at lap start median {np.median(start_x[k]):.1f}, at lap end median {np.median(stop_x[k]):.1f}")
sub-fCamk1_ses-fCamk1-200827-sess9
Site3: 10414 pulses, stored value(s) [19.0] 'watts', duration median 19.87 ms (range 15.77 to 19.87), first 5198.0 s, last 11797.5 s
Site4: 10720 pulses, stored value(s) [19.0] 'watts', duration median 19.87 ms (range 19.83 to 19.87), first 5198.0 s, last 11797.5 s
Site5: 10451 pulses, stored value(s) [19.0] 'watts', duration median 19.87 ms (range 19.83 to 19.87), first 5198.0 s, last 11797.5 s
Site6: 129 pulses, stored value(s) [19.0] 'watts', duration median 19.87 ms (range 19.83 to 19.87), first 5198.0 s, last 11797.5 s
Site7: 216 pulses, stored value(s) [19.0] 'watts', duration median 19.87 ms (range 19.83 to 19.87), first 5198.0 s, last 11797.5 s
Site8: 10158 pulses, stored value(s) [19.0] 'watts', duration median 19.87 ms (range 19.83 to 19.87), first 5198.1 s, last 11797.5 s
stimulation epochs: [('PreStim', 5195, 8285), ('Maze', 8285, 10081), ('PostStim', 10081, 11804)]
linearized position: 53143 samples, 0.0 to 110.0, stated unit 'cm' with conversion 100.0
laps: 111, all inside the Maze epoch: True, directions [55, 56], duration median 12.4 s
laps with fewer than 1 pulse/s (numbered from 1): [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111]
direction 0: position at lap start median 0.0, at lap end median 101.1
direction 1: position at lap start median 101.1, at lap end median 0.0
sub-fCamk1_ses-fCamk1-200901-sess12
Site3: 12628 pulses, stored value(s) [62711.0] 'watts', duration median 19.63 ms (range 19.60 to 19.63), first 4208.8 s, last 11406.5 s
Site4: 12456 pulses, stored value(s) [62623.0] 'watts', duration median 19.63 ms (range 19.60 to 19.63), first 4208.8 s, last 11408.9 s
Site5: 12590 pulses, stored value(s) [62687.0] 'watts', duration median 19.63 ms (range 19.60 to 19.63), first 4208.8 s, last 11408.2 s
Site6: 120 pulses, stored value(s) [62711.0] 'watts', duration median 19.63 ms (range 19.60 to 19.63), first 4208.8 s, last 11392.8 s
Site7: 221 pulses, stored value(s) [62911.0] 'watts', duration median 19.63 ms (range 19.60 to 19.63), first 4208.8 s, last 11403.2 s
Site8: 12294 pulses, stored value(s) [63599.0] 'watts', duration median 19.63 ms (range 19.60 to 19.63), first 4208.8 s, last 11407.5 s
stimulation epochs: [('PreStim', 4203, 8093), ('Maze', 8093, 9638), ('PostStim', 9638, 11409)]
linearized position: 45583 samples, 0.0 to 110.0, stated unit 'cm' with conversion 100.0
laps: 132, all inside the Maze epoch: True, directions [66, 66], duration median 9.8 s
laps with fewer than 1 pulse/s (numbered from 1): [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132]
direction 0: position at lap start median 5.8, at lap end median 100.3
direction 1: position at lap start median 100.3, at lap end median 5.8
The pulse series are stated in watts but hold one value per uLED: 19 in sess9 and 62623 to 63599 in sess12. These cannot be watts (the paper gives 0.02 to 0.1 uW), and with one value per uLED they carry no intensity levels. Each pulse is four samples (off, on, on, off) spanning 19.87 ms (sess9) or 19.63 ms (sess12), consistent with the paper's 20 ms pulses. The files hold series for six uLEDs (Site3 to Site8) of the probe's twelve; Site6 and Site7 have only 120 to 221 pulses. The next check shows that most pulses are not in the files.
The linearized position runs from 0 to 110, the track's length in cm, but carries a conversion factor of 100, which would make it 0 to 110 m; the conversion is wrong. Laps lie inside the Maze epoch and alternate direction, starting and ending near the track ends. Three blocks of laps have no pulses (sess9: laps 1-11, 52-61, 102-111; sess12: 1-10, 61-70, 121-132), the baseline runs of the paper's design. In sess9 they are not at laps 61-70 and 121-130 as in Figure 3F's legend.
# Are the pulses of every uLED in the files? The device, site and series descriptions, the
# deposited pulse rate in PreStim, and the offsets from the end of one deposited pulse to the
# start of the next (any uLED).
desc = {s["key"]: v.describe_file(a["url"]) for a, s in zip(small, sessions) if s["pulses"]}
d0 = desc[opto[0]["key"]]
print("Device:", d0["devices"])
print("Stimulus sites in the file:", list(d0["sites"]), "; description:", sorted({x["description"] for x in d0["sites"].values()}))
print("Pulse series descriptions (distinct over both files):", sorted({x["description"] for d in desc.values() for x in d["series"].values()}))
print("Units columns:", d0["units_columns"], "; electrodes columns:", d0["electrodes_columns"])
edges = np.arange(0, 210, 10)
for s in opto:
pre = v.epoch_interval(s, "PreStim")[0]
iv = np.concatenate([np.column_stack([p["on"], p["off"]]) for p in s["pulses"].values()])
iv = iv[np.argsort(iv[:, 0])]
iv = iv[(iv[:, 0] >= pre[0]) & (iv[:, 0] < pre[1])]
gap = (iv[1:, 0] - iv[:-1, 1]) * 1000
print(f"\n{s['key'].split('_ses-')[1]}: units per shank {dict(zip(*np.unique(s['shank'], return_counts=True)))}; "
f"PreStim deposited pulses {len(iv)} in {pre[1] - pre[0]:.0f} s, {len(iv) / (pre[1] - pre[0]):.2f} per s")
print(" end-to-start offsets in 10 ms bins from 0 to 200 ms (count):", np.histogram(gap, edges)[0].tolist(),
f"; below 0: {int((gap < 0).sum())}, above 200: {int((gap >= 200).sum())}")
print(" fraction in 20-40 ms {:.3f}, 40-60 ms {:.3f}, 60-100 ms {:.3f}, 100 ms or more {:.3f}".format(
*[np.mean((gap >= a) & (gap < b)) for a, b in ((20, 40), (40, 60), (60, 100), (100, np.inf))]))
print(" fraction under 20 ms {:.3f}, under 10 ms or negative {:.3f}".format(np.mean(gap < 20), np.mean(gap < 10)))
if s is opto[0]:
h = np.histogram(gap, np.arange(0, 125, 5))[0]
print(" 5 ms bins from 0 to 120 ms (count):", h.tolist())
Device: {'N1-F21-O36 | 18': '12 µLEDs, 10 x 15 µm each, 3 per shank\nEmission Peak λ = 460 nm and FWHM = 40 nm\nTypical irradiance of 33 mW/mm² (@ max operating current of 100 µA)\n32 recording channels, 8 per shank\nElectrode impedance of 1000 - 1500 kΩ at 1 kHz\n'}
Stimulus sites in the file: ['OptogeneticStimulusSite3', 'OptogeneticStimulusSite4', 'OptogeneticStimulusSite5', 'OptogeneticStimulusSite6', 'OptogeneticStimulusSite7', 'OptogeneticStimulusSite8'] ; description: ['Microled site in Neurolight probe. Microscopic LED 10 x 15 µm each, 3 per shank. Each μLED has an emission area of 150 μm2']
Pulse series descriptions (distinct over both files): ['μLEDs were controlled with current (2-4.5 μA generating 0.02-0.1μW of total light power;ref (15)) provided by a 12-channel current generator (OSC1Lite, NeuroNex Michigan Hub)driven by an Arduino, which delivered trapezoid (1ms rise time)blue light (centered emission at 460 nm, emission surface area = 150 mm2) 20 ms pulses atrandom sites with a randomly variable (40-60ms) offset']
Units columns: ['spike_times', 'unit_name', 'group_id', 'clu_id'] ; electrodes columns: ['location', 'group', 'group_name', 'channel_name', 'rel_x', 'gain_to_uV', 'rel_y']
fCamk1-200827-sess9: units per shank {np.str_('Group1'): np.int64(16), np.str_('Group2'): np.int64(19), np.str_('Group3'): np.int64(9), np.str_('Group4'): np.int64(13)}; PreStim deposited pulses 17796 in 3090 s, 5.76 per s
end-to-start offsets in 10 ms bins from 0 to 200 ms (count): [564, 1, 5584, 0, 1503, 0, 2722, 1, 1221, 0, 1456, 0, 947, 0, 860, 0, 577, 0, 477, 0] ; below 0: 115, above 200: 1767
fraction in 20-40 ms 0.314, 40-60 ms 0.084, 60-100 ms 0.222, 100 ms or more 0.342
fraction under 20 ms 0.038, under 10 ms or negative 0.038
5 ms bins from 0 to 120 ms (count): [564, 0, 0, 1, 5582, 2, 0, 0, 1503, 0, 0, 0, 2722, 0, 1, 0, 1221, 0, 0, 0, 1456, 0, 0, 0]
fCamk1-200901-sess12: units per shank {np.str_('Group1'): np.int64(18), np.str_('Group2'): np.int64(7), np.str_('Group3'): np.int64(10), np.str_('Group4'): np.int64(12)}; PreStim deposited pulses 29186 in 3891 s, 7.50 per s
end-to-start offsets in 10 ms bins from 0 to 200 ms (count): [717, 96, 4292, 4533, 756, 850, 1398, 1768, 1715, 1053, 806, 974, 912, 901, 680, 603, 610, 603, 540, 452] ; below 0: 119, above 200: 4807
fraction in 20-40 ms 0.302, 40-60 ms 0.055, 60-100 ms 0.203, 100 ms or more 0.407
fraction under 20 ms 0.032, under 10 ms or negative 0.029
The device holds 12 uLEDs, 3 per shank, and the files hold pulse series for six of them (Site3 to Site8), all with the same description. Most pulses of the protocol are not in the files. In PreStim the deposited pulses come at 5.76 and 7.50 per s, where 20 ms pulses separated by 20 to 40 ms offsets would come at 17 to 25 per s. In sess9 the offsets between consecutive deposited pulses fall only near 20, 40, 60, 80, 100 ms and so on (5 ms bins above). An offset of 60 ms is a 20 ms offset, a 20 ms pulse and another 20 ms offset, so a pulse lies in it that is not in the file. In sess12, 30.2% of offsets fall within 20 to 40 ms, 5.5% within 40 to 60 ms, and 20.3% within 60 to 100 ms, the range of one missing pulse with its two offsets. No putative pyramidal cell on shank Group3 responds to any deposited uLED (fig1e section), so that shank's uLEDs are likely among the missing ones.
Under the stated protocol the next pulse at any site starts 20 to 40 ms after the end of the last, so the 20 ms before each deposited onset and the 20 ms after each deposited offset are dark. The files do not follow this exactly: 3.8% (sess9) and 3.2% (sess12) of the end-to-start offsets between deposited pulses in PreStim are under 20 ms, nearly all of them under 10 ms or negative. Deposited pulses are cut out of these windows, but pulses missing from the files may fall in them as often. From here on, control time for rate maps and in-field rates is restricted to these windows (about 20% of lap time, printed below). All time without a deposited pulse (88% to 90% of lap time) may hold the light of missing pulses, and is shown as a variant. The 20 ms control windows before each onset, used for fig1i and fig2b, are dark under the protocol.
# Clock checks. (1) Light responses: the population PSTH around the pulse onsets of each uLED,
# against the same onsets shifted by 250 ms. (2) Ripples: the population rate in the deposited
# ripple intervals against the same intervals shifted by 1 s, in every session with ripples.
for s in opto:
pyr_spk = [x for x, p in zip(s["spikes"], s["pyr"]) if p]
pre = v.epoch_interval(s, "PreStim")[0]
for site, p in s["pulses"].items():
on = p["on"][(p["on"] >= pre[0]) & (p["on"] < pre[1])]
if len(on) < v.MIN_PULSES:
continue
true = np.mean([v.count_in_windows(x, on, 0.02).sum() for x in pyr_spk]) / len(on) / 0.02
before = np.mean([v.count_in_windows(x, on - 0.02, 0.02).sum() for x in pyr_spk]) / len(on) / 0.02
shifted = np.mean([v.count_in_windows(x, on + 0.25, 0.02).sum() for x in pyr_spk]) / len(on) / 0.02
print(f"{s['key'].split('_ses-')[1]} {site}: mean pyramidal rate 0-20 ms after onset {true:.2f} Hz, "
f"20 ms before {before:.2f} Hz, onsets shifted by 250 ms {shifted:.2f} Hz")
rip_rows = []
for s in sessions:
r = s["ripples"]
if len(r) < 50:
continue
allspk = np.sort(np.concatenate(s["spikes"]))
dur = (r[:, 1] - r[:, 0]).sum()
inside = (np.searchsorted(allspk, r[:, 1]) - np.searchsorted(allspk, r[:, 0])).sum() / dur
sh = (np.searchsorted(allspk, r[:, 1] + 1) - np.searchsorted(allspk, r[:, 0] + 1)).sum() / dur
rip_rows.append({"session": s["key"].split("_ses-")[1], "ripples": len(r),
"duration_ms_median": round(1000 * np.median(r[:, 1] - r[:, 0]), 1),
"peak_inside": bool(np.all((r[:, 2] >= r[:, 0]) & (r[:, 2] <= r[:, 1]))),
"pop_rate_in_ripples": round(inside, 1), "pop_rate_shifted_1s": round(sh, 1)})
rip_check = pd.DataFrame(rip_rows)
print(rip_check.to_string(index=False))
print("Sessions where the population rate inside ripples is at least twice the shifted rate:",
int((rip_check.pop_rate_in_ripples > 2 * rip_check.pop_rate_shifted_1s).sum()), "of", len(rip_check))
fCamk1-200827-sess9 Site3: mean pyramidal rate 0-20 ms after onset 2.83 Hz, 20 ms before 1.58 Hz, onsets shifted by 250 ms 1.83 Hz fCamk1-200827-sess9 Site4: mean pyramidal rate 0-20 ms after onset 1.84 Hz, 20 ms before 1.58 Hz, onsets shifted by 250 ms 1.83 Hz fCamk1-200827-sess9 Site5: mean pyramidal rate 0-20 ms after onset 2.52 Hz, 20 ms before 1.52 Hz, onsets shifted by 250 ms 1.83 Hz fCamk1-200827-sess9 Site8: mean pyramidal rate 0-20 ms after onset 1.81 Hz, 20 ms before 1.62 Hz, onsets shifted by 250 ms 1.80 Hz
fCamk1-200901-sess12 Site3: mean pyramidal rate 0-20 ms after onset 1.47 Hz, 20 ms before 1.25 Hz, onsets shifted by 250 ms 1.51 Hz fCamk1-200901-sess12 Site4: mean pyramidal rate 0-20 ms after onset 1.94 Hz, 20 ms before 1.24 Hz, onsets shifted by 250 ms 1.48 Hz
fCamk1-200901-sess12 Site5: mean pyramidal rate 0-20 ms after onset 1.38 Hz, 20 ms before 1.28 Hz, onsets shifted by 250 ms 1.48 Hz fCamk1-200901-sess12 Site8: mean pyramidal rate 0-20 ms after onset 2.23 Hz, 20 ms before 1.28 Hz, onsets shifted by 250 ms 1.49 Hz
session ripples duration_ms_median peak_inside pop_rate_in_ripples pop_rate_shifted_1s
fCamk1-200827-sess9 4315 36.8 True 352.6 112.3
fCamk1-200901-sess12 1796 31.2 True 293.8 112.5
fCamk1-200902-sess13 7196 36.8 True 288.9 104.3
fCamk1-200904-sess15 6258 39.2 True 332.9 129.9
fCamk1-200908-sess16 8096 36.0 True 322.1 120.1
fCamk1-200909-sess17 4096 29.6 True 320.4 118.2
fCamk1-200910-sess18 13117 37.6 True 313.9 114.4
fCamk1-200911-sess19 10448 36.0 True 227.3 70.6
fCamk2-201012-sess1 1110 48.0 True 189.7 52.5
fCamk3-201028-sess10-cleanned 12881 36.0 True 284.6 96.3
fCamk3-201029-sess11-cleanned 9756 36.8 True 253.5 94.6
fCamk3-201030-sess12 11835 36.0 True 243.7 92.8
fCamk3-201102-sess13 2123 54.4 True 202.0 155.9
fCamk3-201103-sess14 12910 38.4 True 278.0 93.6
fCamk3-201105-sess16 13274 38.4 True 286.1 109.5
fCamk3-201106-sess17 9056 36.8 True 199.5 100.7
fCamk3-201109-sess18 13753 36.8 True 236.3 110.1
fCamk3-201110-sess19 15205 37.6 True 230.6 92.7
fCamk3-201111-sess20 5024 50.4 True 198.0 86.9
fCamk3-201113-sess22 8558 37.6 True 236.7 100.0
fCamk5-210406-sess10 217 108.8 True 277.3 226.8
fCamk5-210407-sess11 6093 36.8 True 231.6 96.0
fCamk5-210408-sess12 6842 36.8 True 299.7 118.0
fCamk5-210412-sess14 5718 36.8 True 196.1 84.9
fCamk5-210414-sess16 5023 38.4 True 177.7 84.0
Sessions where the population rate inside ripples is at least twice the shifted rate: 22 of 25
Pulse onsets share the clock of the spikes: for every uLED tested, the mean rate of the putative pyramidal cells in the 20 ms after onset is above the rate in the 20 ms before (by 0.10 to 1.25 Hz). For 6 of 8 it is also above the rate in the same window 250 ms later. The other two are sess12 Site3 and Site5, and the later window often holds pulses of other uLEDs. Responses per unit are tested under fig1e. The deposited ripples share the clock of the spikes: the population rate inside them is at least twice the rate 1 s later in 22 of 25 sessions with ripples. In fCamk3 sess13 (1.3 times), fCamk3 sess17 (1.98 times) and fCamk5 sess10 (1.2 times; 217 events with a median duration of 109 ms) the ripple intervals mark population bursts weakly. The LFP clock is checked in the fig2b section: the stored LFP rate is 1250 Hz from 0 s, the LFP and the raw signal both end with the last epoch, and ripple-band power at the deposited ripple peaks is 194 and 248 times its mean in BaselinePre and 113 and 78 times inside the Maze epoch. Position shares the clock of the laps (lap ends at track ends, above).
# Responsiveness of every unit to every uLED, in the PreStim epoch (main) and over all
# stimulation epochs (reading b of fig1e). One generator per session.
for i, s in enumerate(opto):
rng = np.random.default_rng(1000 + i)
pre = v.epoch_interval(s, "PreStim")[0]
s["resp_table"] = pd.DataFrame(v.responsive_sites(s, pre, rng))
stim = [(e["start"], e["stop"]) for e in s["epochs"] if "uLED" in e["manipulation"]]
whole = (min(a for a, _ in stim), max(b for _, b in stim))
s["resp_table_all"] = pd.DataFrame(v.responsive_sites(s, whole, rng))
t = s["resp_table"]
s["resp_sites"] = {u: sorted(g[g.responsive.fillna(False).astype(bool)].site) for u, g in t.groupby("unit")}
s["responsive"] = np.array([len(s["resp_sites"].get(u, [])) > 0 for u in range(len(s["spikes"]))])
s["tested_sites"] = sorted(t[t.tested].site.unique())
print(s["key"], "tested uLEDs:", s["tested_sites"], "; responsive units:", int(s["responsive"].sum()),
"of", len(s["spikes"]), "; responsive putative pyramidal:", int((s["responsive"] & s["pyr"]).sum()),
"of", int(s["pyr"].sum()))
sub-fCamk1_ses-fCamk1-200827-sess9 tested uLEDs: ['Site3', 'Site4', 'Site5', 'Site8'] ; responsive units: 39 of 57 ; responsive putative pyramidal: 35 of 51
sub-fCamk1_ses-fCamk1-200901-sess12 tested uLEDs: ['Site3', 'Site4', 'Site5', 'Site8'] ; responsive units: 25 of 47 ; responsive putative pyramidal: 22 of 40
Claim cohort-counts¶
CA1 pyramidal neurons were recorded and probed in four CaMKIIa-Cre::Ai32 mice, 822 pyramidal neurons in all, 43.3 plus or minus 8.37 per session.
def cohort(rows):
return {"sessions": len(rows), "mice": rows.subject.nunique(), "pyramidal": int(rows.putative_pyr.sum()),
"per_session_mean": round(rows.putative_pyr.mean(), 1), "per_session_sd": round(rows.putative_pyr.std(), 2)}
readings = {"(a) all sessions": overview, "(b) sessions naming the uLED protocol": overview[overview.uLED_epochs]}
res = pd.DataFrame({k: cohort(r) for k, r in readings.items()}).T
res["mice_ok"] = res.mice == 4
res["total_ok"] = (res.pyramidal - 822).abs() <= 82
res["mean_ok"] = (res.per_session_mean - 43.3).abs() <= 8.37
print(res.to_string())
print("\nPer subject, reading (a):")
print(overview.groupby("subject").agg(sessions=("session", "size"), units=("units", "sum"),
putative_pyr=("putative_pyr", "sum")).to_string())
print("\nVariant thresholds for putative pyramidal cells, reading (a) and (b):")
for thr in (2.5, 10):
n = [int(sum((s["base_rate"] < thr).sum() for s in sessions)),
int(sum((s["base_rate"] < thr).sum() for s in sessions if s["stim_protocol"]))]
print(f" below {thr} Hz: {n[0]} and {n[1]}")
sessions mice pyramidal per_session_mean per_session_sd mice_ok total_ok mean_ok
(a) all sessions 28.0 4.0 1351.0 48.2 12.18 True False True
(b) sessions naming the uLED protocol 22.0 3.0 1047.0 47.6 8.56 False False True
Per subject, reading (a):
sessions units putative_pyr
subject
fCamk1 8 409 370
fCamk2 4 162 157
fCamk3 11 615 555
fCamk5 5 295 269
Variant thresholds for putative pyramidal cells, reading (a) and (b):
below 2.5 Hz: 1205 and 940
below 10 Hz: 1441 and 1114
Result: partially reproduced. Under reading (a), all 28 sessions, the dandiset holds 4 mice and 48.2 putative pyramidal cells per session (SD 12.18), within 43.3 plus or minus 8.37, but 1351 in all, against 822 plus or minus 82. Under reading (b), the 22 sessions naming the uLED protocol, it holds 3 mice (fCamk2 has no epoch labels), 1047 cells and 47.6 per session. At 43.3 per session, 822 cells would come from about 19 sessions; the dandiset holds 28 and does not say which the paper used. The thresholds 2.5 and 10 Hz give 1205 and 1441 cells under reading (a).
Claim fig1c-pulse-protocol¶
uLEDs were driven with 20 ms pulses at random sites with randomly variable 20 to 40 ms offsets, so that the stimulation of each site recurred at about 0.3 to 0.6 s intervals.
rows = []
for s in opto:
stim = [(e["start"], e["stop"]) for e in s["epochs"] if "uLED" in e["manipulation"]]
allon = []
for site, p in s["pulses"].items():
keep = np.zeros(len(p["on"]), bool)
for a, b in stim:
keep |= (p["on"] >= a) & (p["on"] < b)
on, off = p["on"][keep], p["off"][keep]
allon.append(np.column_stack([on, off]))
ipi = np.diff(on)
rows.append({"session": s["key"].split("_ses-")[1], "uLED": site, "pulses": len(on),
"duration_ms": round(1000 * np.median(off - on), 2),
"interval_median_s": round(np.median(ipi), 3),
"interval_5_95_s": np.percentile(ipi, [5, 95]).round(3).tolist(),
"intervals_in_0.3_0.6": round(np.mean((ipi >= 0.3) & (ipi <= 0.6)), 3)})
iv = np.concatenate(allon)
iv = iv[np.argsort(iv[:, 0])]
gap = (iv[1:, 0] - iv[:-1, 1]) * 1000
print(f"{s['key'].split('_ses-')[1]}: offset from the end of a pulse to the start of the next (any uLED), "
f"median {np.median(gap):.1f} ms, 5th to 95th percentile {np.percentile(gap, 5):.1f} to {np.percentile(gap, 95):.1f} ms, "
f"overlapping pulses {int((gap < 0).sum())} of {len(gap)}; within 20 to 40 ms {np.mean((gap >= 20) & (gap <= 40)):.1%}, "
f"within 40 to 60 ms (file description) {np.mean((gap >= 40) & (gap <= 60)):.1%}; onset to next onset "
f"within 40 to 60 ms {np.mean((np.diff(iv[:, 0]) >= 0.04) & (np.diff(iv[:, 0]) <= 0.06)):.1%}")
s["gap_ms"] = gap
pp = pd.DataFrame(rows)
print(pp.to_string(index=False))
used = pp[pp.pulses >= v.MIN_PULSES]
print("\nMedian duration over uLEDs:", used.duration_ms.median(), "ms")
print("Median per-uLED interval of uLEDs with at least 1000 pulses, range:", used.interval_median_s.min(), "to", used.interval_median_s.max(), "s")
fig, axs = plt.subplots(1, 2, figsize=(9, 3))
for s in opto:
g = s["gap_ms"]
axs[0].hist(np.clip(g, 0, 400), bins=80, histtype="step", label=s["key"].split("_ses-")[1])
axs[0].set_xlabel("offset to next pulse, any uLED (ms; clipped at 400)")
axs[0].legend(fontsize=7)
s = opto[0]
p = s["pulses"]["Site3"]
axs[1].hist(np.clip(np.diff(p["on"]), 0, 3), bins=100, color="0.4")
axs[1].set_xlabel("interval between Site3 pulses, sess9 (s; clipped at 3)")
plt.tight_layout()
plt.show()
print("Values above the axes (clipped):", [int((s["gap_ms"] > 400).sum()) for s in opto], "offsets;",
int((np.diff(p["on"]) > 3).sum()), "intervals")
fCamk1-200827-sess9: offset from the end of a pulse to the start of the next (any uLED), median 60.4 ms, 5th to 95th percentile 20.2 to 241.1 ms, overlapping pulses 257 of 42087; within 20 to 40 ms 31.8%, within 40 to 60 ms (file description) 8.2%; onset to next onset within 40 to 60 ms 31.8%
fCamk1-200901-sess12: offset from the end of a pulse to the start of the next (any uLED), median 80.0 ms, 5th to 95th percentile 21.6 to 328.1 ms, overlapping pulses 241 of 50308; within 20 to 40 ms 30.2%, within 40 to 60 ms (file description) 5.5%; onset to next onset within 40 to 60 ms 30.2%
session uLED pulses duration_ms interval_median_s interval_5_95_s intervals_in_0.3_0.6
fCamk1-200827-sess9 Site3 10414 19.87 0.301 [0.04, 1.184] 0.264
fCamk1-200827-sess9 Site4 10720 19.87 0.281 [0.04, 1.146] 0.270
fCamk1-200827-sess9 Site5 10451 19.87 0.301 [0.04, 1.184] 0.279
fCamk1-200827-sess9 Site6 129 19.87 25.775 [1.501, 109.239] 0.000
fCamk1-200827-sess9 Site7 216 19.87 15.712 [0.811, 72.564] 0.023
fCamk1-200827-sess9 Site8 10158 19.87 0.301 [0.04, 1.244] 0.275
fCamk1-200901-sess12 Site3 12628 19.63 0.368 [0.052, 1.535] 0.259
fCamk1-200901-sess12 Site4 12456 19.63 0.380 [0.053, 1.534] 0.266
fCamk1-200901-sess12 Site5 12590 19.63 0.380 [0.052, 1.528] 0.265
fCamk1-200901-sess12 Site6 120 19.63 42.775 [4.997, 187.96] 0.000
fCamk1-200901-sess12 Site7 221 19.63 21.341 [2.116, 103.517] 0.005
fCamk1-200901-sess12 Site8 12294 19.63 0.383 [0.053, 1.581] 0.261
Median duration over uLEDs: 19.75 ms
Median per-uLED interval of uLEDs with at least 1000 pulses, range: 0.281 to 0.383 s
Values above the axes (clipped): [363, 1260] offsets; 16 intervals
Result: partially reproduced. Pulses last 19.75 ms (median over uLEDs), and the median interval between pulses of one uLED is 0.281 to 0.383 s for the four uLEDs with at least 1000 pulses, both as the paper states. The median offset from the end of one pulse to the start of the next is 60.4 ms in sess9 and 80.0 ms in sess12, outside 18 to 44 ms. Of these offsets 31.8% and 30.2% fall within 20 to 40 ms, and the 95th percentiles are 241 and 328 ms. The data check shows that most pulses are not in the files (5.76 and 7.50 deposited pulses per s in PreStim, and offsets in sess9 near multiples of 20 ms only), so the longer offsets span missing pulses, and the protocol's offsets cannot be measured from the files. Measured onset to onset, the same 31.8% and 30.2% fall within the 40 to 60 ms of the file's description, so the paper's range and the file's agree if the file's is onset to onset.
Claim fig1e-responsive-fraction¶
Of 822 pyramidal neurons, 611 responded to light from the uLEDs.
rows = []
for s in opto:
for label, t in (("(a) PreStim", s["resp_table"]), ("(b) all stimulation epochs", s["resp_table_all"])):
resp = t[t.tested].groupby("unit").responsive.any().reindex(range(len(s["spikes"])), fill_value=False).values
rows.append({"session": s["key"].split("_ses-")[1], "reading": label, "pyr": int(s["pyr"].sum()),
"pyr_responsive": int((resp & s["pyr"]).sum()), "all_units": len(resp), "all_responsive": int(resp.sum())})
fr = pd.DataFrame(rows)
print(fr.to_string(index=False))
tot = fr.groupby("reading")[["pyr", "pyr_responsive", "all_units", "all_responsive"]].sum()
tot["pct_pyr_responsive"] = (100 * tot.pyr_responsive / tot.pyr).round(1)
tot["within_tolerance"] = (tot.pct_pyr_responsive - 74.3).abs() <= 7.4
tot["ci95"] = [ci95(a, b) for a, b in zip(tot.pyr_responsive, tot.pyr)]
print(tot.to_string())
# Which part of the three-way test excludes cells: PreStim, putative pyramidal cells
parts = []
for s in opto:
t = s["resp_table"][s["resp_table"].tested].copy()
t = t[t.unit.isin(np.flatnonzero(s["pyr"]))]
parts.append(t)
parts = pd.concat(parts)
print("\nPutative pyramidal unit-uLED tests in PreStim:", len(parts),
"; above the bootstrap bound:", int((parts.rate_during > parts.boot_hi).sum()),
"; z > 1.96:", int((parts.z > v.Z_THRESHOLD).sum()), "; KS p < 0.01:", int((parts.ks_p < v.KS_ALPHA).sum()),
"; all three:", int(parts.responsive.sum()))
print("Variant, z > 1.96 alone (any uLED), PreStim:",
sum(int(((s["resp_table"][s["resp_table"].tested].groupby("unit").z.max() > v.Z_THRESHOLD)
.reindex(range(len(s["spikes"])), fill_value=False).values & s["pyr"]).sum()) for s in opto),
"of", int(sum(s["pyr"].sum() for s in opto)))# Responsive putative pyramidal cells by shank (PreStim), and the uLEDs they respond to
sh = pd.concat([pd.DataFrame({"session": s["key"].split("_ses-")[1], "shank": s["shank"][s["pyr"]],
"responsive": s["responsive"][s["pyr"]]}) for s in opto])
print("\nPutative pyramidal cells by shank, responsive (PreStim):")
print(sh.groupby(["session", "shank"]).responsive.agg(["sum", "size"]).T.to_string())
for label, key in (("(a) PreStim", "resp_table"), ("(b) all stimulation epochs", "resp_table_all")):
k_ = n_ = 0
for s in opto:
t = s[key]
resp = t[t.tested].groupby("unit").responsive.any().reindex(range(len(s["spikes"])), fill_value=False).values
keep = s["pyr"] & (s["shank"] != "Group3")
k_, n_ = k_ + int((resp & keep).sum()), n_ + int(keep.sum())
print(f"Variant, without shank Group3, {label}: {k_} of {n_} ({100 * k_ / n_:.1f}%, 95% interval {ci95(k_, n_)})")
pairs_sh = pd.concat([t.assign(shank=s["shank"][t.unit.values], session=s["key"].split("_ses-")[1])
for s in opto for t in [s["resp_table"][s["resp_table"].tested & s["resp_table"].responsive.fillna(False).astype(bool)]]],
ignore_index=True)
pairs_sh = pairs_sh[[bool(s["pyr"][u]) for u, sess in zip(pairs_sh.unit, pairs_sh.session) for s in opto if s["key"].endswith(sess)]]
print("Responsive unit-uLED pairs of putative pyramidal cells, shank by uLED:")
print(pd.crosstab([pairs_sh.session, pairs_sh.shank], pairs_sh.site).to_string())
# Null for the test: the same test on sham onsets. (1) Each PreStim onset jittered by up to
# +-200 ms; (2) the PreStim onsets moved into BaselinePre, which has no light.
rng = np.random.default_rng(1010)
sham = []
for s in opto:
pre = v.epoch_interval(s, "PreStim")[0]
base = v.epoch_interval(s, "BaselinePre")[0]
jit = {k: np.sort(p["on"] + rng.uniform(-0.2, 0.2, len(p["on"]))) for k, p in s["pulses"].items()}
moved = {k: p["on"] - pre[0] + base[0] for k, p in s["pulses"].items()}
for name, ons, win in (("jittered +-200 ms, PreStim", jit, pre),
("moved into BaselinePre", moved, (base[0], min(base[1], base[0] + pre[1] - pre[0])))):
t = pd.DataFrame(v.responsive_sites(s, win, rng, ons))
t = t[t.tested]
resp = t.groupby("unit").responsive.any().reindex(range(len(s["spikes"])), fill_value=False).values
tp = t[t.unit.isin(np.flatnonzero(s["pyr"]))]
sham.append({"session": s["key"].split("_ses-")[1], "sham": name, "uLEDs_tested": t.site.nunique(),
"pyr_unit_uLED_tests": len(tp), "tests_passed": int(tp.responsive.sum()),
"pyr": int(s["pyr"].sum()), "pyr_passing_any": int((resp & s["pyr"]).sum())})
sham = pd.DataFrame(sham)
print("\nSham onsets, the same test:")
print(sham.to_string(index=False))
st_ = sham.groupby("sham")[["pyr_unit_uLED_tests", "tests_passed", "pyr", "pyr_passing_any"]].sum()
st_["pct_pyr_passing"] = (100 * st_.pyr_passing_any / st_.pyr).round(1)
st_["pct_tests_passed"] = (100 * st_.tests_passed / st_.pyr_unit_uLED_tests).round(1)
print(st_.to_string())
# The same null for reading (b): every onset of the stimulation epochs jittered by up to +-200 ms.
rng = np.random.default_rng(1011)
sham_b = []
for s in opto:
stim = [(e["start"], e["stop"]) for e in s["epochs"] if "uLED" in e["manipulation"]]
whole = (min(a for a, _ in stim), max(b for _, b in stim))
jit = {k: np.sort(p["on"] + rng.uniform(-0.2, 0.2, len(p["on"]))) for k, p in s["pulses"].items()}
t = pd.DataFrame(v.responsive_sites(s, whole, rng, jit))
t = t[t.tested]
resp = t.groupby("unit").responsive.any().reindex(range(len(s["spikes"])), fill_value=False).values
tp = t[t.unit.isin(np.flatnonzero(s["pyr"]))]
sham_b.append({"session": s["key"].split("_ses-")[1], "sham": "jittered +-200 ms, all stimulation epochs",
"pyr_unit_uLED_tests": len(tp), "tests_passed": int(tp.responsive.sum()),
"pyr": int(s["pyr"].sum()), "pyr_passing_any": int((resp & s["pyr"]).sum())})
print("\nSham onsets for reading (b):")
print(pd.DataFrame(sham_b).to_string(index=False))
session reading pyr pyr_responsive all_units all_responsive
fCamk1-200827-sess9 (a) PreStim 51 35 57 39
fCamk1-200827-sess9 (b) all stimulation epochs 51 40 57 44
fCamk1-200901-sess12 (a) PreStim 40 22 47 25
fCamk1-200901-sess12 (b) all stimulation epochs 40 26 47 29
pyr pyr_responsive all_units all_responsive pct_pyr_responsive within_tolerance ci95
reading
(a) PreStim 91 57 104 64 62.6 False 51.9 to 72.6%
(b) all stimulation epochs 91 66 104 73 72.5 True 62.2 to 81.4%
Putative pyramidal unit-uLED tests in PreStim: 364 ; above the bootstrap bound: 91 ; z > 1.96: 71 ; KS p < 0.01: 91 ; all three: 68
Variant, z > 1.96 alone (any uLED), PreStim: 59 of 91
Putative pyramidal cells by shank, responsive (PreStim):
session fCamk1-200827-sess9 fCamk1-200901-sess12
shank Group1 Group2 Group3 Group4 Group1 Group2 Group3 Group4
sum 13 14 0 8 9 4 0 9
size 14 18 8 11 15 6 9 10
Variant, without shank Group3, (a) PreStim: 57 of 74 (77.0%, 95% interval 65.8 to 86.0%)
Variant, without shank Group3, (b) all stimulation epochs: 66 of 74 (89.2%, 95% interval 79.8 to 95.2%)
Responsive unit-uLED pairs of putative pyramidal cells, shank by uLED:
site Site3 Site4 Site5 Site8
session shank
fCamk1-200827-sess9 Group1 13 6 0 0
Group2 0 0 14 0
Group4 0 0 0 8
fCamk1-200901-sess12 Group1 5 9 0 0
Group2 0 0 4 0
Group4 0 0 0 9
Sham onsets, the same test:
session sham uLEDs_tested pyr_unit_uLED_tests tests_passed pyr pyr_passing_any
fCamk1-200827-sess9 jittered +-200 ms, PreStim 4 204 0 51 0
fCamk1-200827-sess9 moved into BaselinePre 4 204 0 51 0
fCamk1-200901-sess12 jittered +-200 ms, PreStim 4 160 0 40 0
fCamk1-200901-sess12 moved into BaselinePre 4 160 0 40 0
pyr_unit_uLED_tests tests_passed pyr pyr_passing_any pct_pyr_passing pct_tests_passed
sham
jittered +-200 ms, PreStim 364 0 91 0 0.0 0.0
moved into BaselinePre 364 0 91 0 0.0 0.0
Sham onsets for reading (b):
session sham pyr_unit_uLED_tests tests_passed pyr pyr_passing_any
fCamk1-200827-sess9 jittered +-200 ms, all stimulation epochs 204 0 51 0
fCamk1-200901-sess12 jittered +-200 ms, all stimulation epochs 160 0 40 0
Result: reproduced under reading (b). Over all stimulation epochs, 66 of 91 putative pyramidal cells (72.5%, 95% interval 62.2 to 81.4%) respond to at least one uLED, within 74.3 plus or minus 7.4%. Under reading (a), the PreStim epoch alone, 57 of 91 (62.6%, 51.9 to 72.6%) do, outside tolerance. Of 364 unit and uLED tests in PreStim, 68 pass all three parts of the test. On sham onsets, jittered by up to 200 ms or moved into BaselinePre, none of the 364 PreStim tests passes. For reading (b), onsets over all stimulation epochs jittered by up to 200 ms, none of 364 tests passes either. No putative pyramidal cell on shank Group3 responds (0 of 8 and 0 of 9), while the other shanks respond each to their own uLEDs (Site3 and Site4 on Group1, Site5 on Group2, Site8 on Group4). The uLEDs of Group3 are likely among those whose pulses are not in the files, so some cells counted as non-responsive may respond to light that cannot be tested here. Without Group3, 57 of 74 (77.0%) respond under reading (a), within tolerance, and 66 of 74 (89.2%) under reading (b), above it. These are two sessions of one mouse, not the paper's 822 cells.
# Light and control windows of each session, used by the claims below.
for s in opto:
s["control"] = v.control_onsets(s)
s["light_iv"] = v.all_light_intervals(s)
s["dark_iv"] = v.dark_intervals(s)
s["cell_onsets"] = {u: np.sort(np.concatenate([s["pulses"][k]["on"] for k in sites]))
for u, sites in s["resp_sites"].items() if sites}
s["all_tested_onsets"] = np.sort(np.concatenate([s["pulses"][k]["on"] for k in s["tested_sites"]]))
print(s["key"].split("_ses-")[1], "control windows:", len(s["control"]),
"; light-on time fraction over the stimulation epochs:",
round((s["light_iv"][:, 1] - s["light_iv"][:, 0]).sum() /
sum(e["stop"] - e["start"] for e in s["epochs"] if "uLED" in e["manipulation"]), 3),
"; dark under the protocol (20 ms before each onset and after each offset):",
round((s["dark_iv"][:, 1] - s["dark_iv"][:, 0]).sum() /
sum(e["stop"] - e["start"] for e in s["epochs"] if "uLED" in e["manipulation"]), 3))
fCamk1-200827-sess9 control windows: 40525 ; light-on time fraction over the stimulation epochs: 0.126 ; dark under the protocol (20 ms before each onset and after each offset): 0.205
fCamk1-200901-sess12 control windows: 48645 ; light-on time fraction over the stimulation epochs: 0.136 ; dark under the protocol (20 ms before each onset and after each offset): 0.25
Claim fig1i-swr-decrease¶
Light-induced spike responses of individual pyramidal cells decreased during sharp-wave ripples (n = 485 neurons, P < 1e-55, Wilcoxon paired signed-rank test).
def overlaps(starts, width, iv):
# whether [start, start + width) overlaps any interval of iv (sorted by start)
j = np.searchsorted(iv[:, 0], starts + width) - 1
ok = j >= 0
out = np.zeros(len(starts), bool)
out[ok] = iv[j[ok], 1] > starts[ok]
return out
rows = []
for s in opto:
home = [(e["start"], e["stop"]) for e in s["epochs"] if e["behavioral_paradigm"] in ("PreStim", "PostStim")]
rip = s["ripples"][:, :2]
ctrl = v.restrict(s["control"], home)
c_in = overlaps(ctrl, v.CONTROL_WINDOW_S, rip)
for u in np.flatnonzero(s["responsive"] & s["pyr"]):
on = v.restrict(s["cell_onsets"][u], home)
l_in = overlaps(on, v.PULSE_WINDOW_S, rip)
x = s["spikes"][u]
lc, cc = v.count_in_windows(x, on, 0.02), v.count_in_windows(x, ctrl, 0.02)
rows.append({"session": s["key"].split("_ses-")[1], "unit": u,
"pulses_in": int(l_in.sum()), "pulses_out": int((~l_in).sum()),
"light_in": lc[l_in].mean() / 0.02, "light_out": lc[~l_in].mean() / 0.02,
"ctrl_in": cc[c_in].mean() / 0.02, "ctrl_out": cc[~c_in].mean() / 0.02,
"ctrl_windows_in": int(c_in.sum())})
swr = pd.DataFrame(rows)
swr["d_in"] = swr.light_in - swr.ctrl_in
swr["d_out"] = swr.light_out - swr.ctrl_out
swr["g_in"] = swr.light_in / swr.ctrl_in.replace(0, np.nan)
swr["g_out"] = swr.light_out / swr.ctrl_out.replace(0, np.nan)
print("Cells:", len(swr), "; pulses overlapping a ripple per cell, median", swr.pulses_in.median(),
"(range", swr.pulses_in.min(), "to", swr.pulses_in.max(), "); control windows in ripples per session:",
swr.groupby("session").ctrl_windows_in.first().to_dict())
print(swr[["light_in", "ctrl_in", "light_out", "ctrl_out", "d_in", "d_out"]].median().round(2).to_string())
for name, a, b in (("(a) difference light - control (Hz)", "d_in", "d_out"), ("(b) ratio light / control", "g_in", "g_out")):
k = swr[[a, b]].dropna()
k = k[np.isfinite(k).all(axis=1)]
diff = k[a] - k[b]
w = stats.wilcoxon(k[a], k[b])
print(f"{name}: n = {len(k)}, median inside {k[a].median():.2f}, outside {k[b].median():.2f}, "
f"median paired difference {diff.median():.2f}, lower inside in {(diff < 0).sum()} of {len(k)}, "
f"Wilcoxon p = {w.pvalue:.2g}")
for sess, g in k.assign(session=swr.loc[k.index, "session"]).groupby("session"):
dd = g[a] - g[b]
print(f" {sess}: n = {len(g)}, median difference {dd.median():.2f}, lower inside {(dd < 0).sum()}, "
f"p = {stats.wilcoxon(g[a], g[b]).pvalue:.2g}" if len(g) > 5 else f" {sess}: n = {len(g)}")
fig, ax = plt.subplots(figsize=(4, 4))
ax.scatter(swr.d_out, swr.d_in, s=10, c="C0")
lim = [min(swr.d_out.min(), swr.d_in.min()), max(swr.d_out.max(), swr.d_in.max())]
ax.plot(lim, lim, "k:", lw=1)
ax.set_xlabel("light - control outside ripples (Hz)")
ax.set_ylabel("light - control inside ripples (Hz)")
plt.show()
print("Axes span every point: x", np.round([swr.d_out.min(), swr.d_out.max()], 1), "y", np.round([swr.d_in.min(), swr.d_in.max()], 1))
Cells: 57 ; pulses overlapping a ripple per cell, median 111.0 (range 67 to 234 ); control windows in ripples per session: {'fCamk1-200827-sess9': 421, 'fCamk1-200901-sess12': 272}
light_in 2.78
ctrl_in 2.14
light_out 4.13
ctrl_out 1.34
d_in 0.30
d_out 2.68
(a) difference light - control (Hz): n = 57, median inside 0.30, outside 2.68, median paired difference -2.13, lower inside in 49 of 57, Wilcoxon p = 6.8e-09
fCamk1-200827-sess9: n = 35, median difference -1.80, lower inside 28, p = 7.8e-06
fCamk1-200901-sess12: n = 22, median difference -2.18, lower inside 21, p = 9.5e-07
(b) ratio light / control: n = 56, median inside 1.25, outside 2.89, median paired difference -1.81, lower inside in 53 of 56, Wilcoxon p = 8.6e-10
fCamk1-200827-sess9: n = 35, median difference -1.84, lower inside 33, p = 8.6e-08
fCamk1-200901-sess12: n = 21, median difference -1.70, lower inside 20, p = 1.9e-06
Axes span every point: x [ 0.6 15.6] y [-4.4 5.1]
Result: reproduced. In 57 light-responsive putative pyramidal cells, the light minus control rate is 0.30 Hz (median) in windows overlapping ripples and 2.68 Hz outside them. It is lower inside ripples in 49 of 57 cells (Wilcoxon p = 6.8e-09), and in each session (p = 7.8e-06 and 9.5e-07). As a ratio, light over control is 1.25 inside and 2.89 outside (53 of 56 lower, p = 8.6e-10). The control windows, the 20 ms before each onset, are dark under the protocol, so the missing pulses do not enter them. The paper's n of 485 and P < 1e-55 are out of reach of two sessions.
Claim figs10-swr-gain-vs-baseline-rate¶
The firing-rate change of pyramidal neurons during SPW-Rs is negatively correlated with their baseline firing rate (rho = -0.19, P < 1e-7).
rows = []
for s in sessions:
e0 = s["epochs"][0]
r = s["ripples"]
r = r[(r[:, 0] >= e0["start"]) & (r[:, 1] < e0["stop"])]
if len(r) < 50:
continue
t_in = (r[:, 1] - r[:, 0]).sum()
t_out = (e0["stop"] - e0["start"]) - t_in
for u in np.flatnonzero(s["pyr"]):
x = s["spikes"][u]
n_all = ((x >= e0["start"]) & (x < e0["stop"])).sum()
n_in = (np.searchsorted(x, r[:, 1]) - np.searchsorted(x, r[:, 0])).sum()
rows.append({"session": s["key"].split("_ses-")[1], "subject": s["subject"], "base": n_all / (e0["stop"] - e0["start"]),
"in": n_in / t_in, "out": (n_all - n_in) / t_out, "ripples": len(r)})
g = pd.DataFrame(rows)
g["log2_ratio"] = np.log2(g["in"] / g["out"])
g["diff"] = g["in"] - g["out"]
print("Sessions:", g.session.nunique(), "; putative pyramidal cells:", len(g),
"; cells with no spikes in ripples (left out of reading a):", int((g["in"] == 0).sum()))
for name, col in (("(a) log2 ratio in/out", "log2_ratio"), ("(b) difference in - out", "diff")):
k = g[np.isfinite(g[col])]
rho = stats.spearmanr(k[col], k.base)
print(f"{name}: n = {len(k)}, rho = {rho.statistic:.3f}, p = {rho.pvalue:.2g}, "
f"within -0.19 +- 0.068: {abs(rho.statistic + 0.19) <= 0.068}")
print("\nPer subject and per session, reading (a):")
k = g[np.isfinite(g.log2_ratio)]
for key, grp in list(k.groupby("subject")) + list(k.groupby("session")):
rr = stats.spearmanr(grp.log2_ratio, grp.base)
print(f" {key}: n = {len(grp)}, rho = {rr.statistic:.2f}, p = {rr.pvalue:.2g}")
bad = rip_check[rip_check.pop_rate_in_ripples < 2 * rip_check.pop_rate_shifted_1s].session.tolist()
kk = k[~k.session.isin(bad)]
rr = stats.spearmanr(kk.log2_ratio, kk.base)
print(f"Variant, without the sessions whose ripples fail the clock check {bad}: n = {len(kk)}, rho = {rr.statistic:.3f}, p = {rr.pvalue:.2g}")
labelled = [s["key"].split("_ses-")[1] for s in sessions if s["epochs"][0].get("behavioral_paradigm") == "BaselinePre"]
kl = k[k.session.isin(labelled)]
rr = stats.spearmanr(kl.log2_ratio, kl.base)
print(f"Variant, sessions whose first epoch is labelled BaselinePre ({kl.session.nunique()} sessions): n = {len(kl)}, rho = {rr.statistic:.3f}, p = {rr.pvalue:.2g}")
kk = g[g.ripples >= 1000]
kk = kk[np.isfinite(kk.log2_ratio)]
rr = stats.spearmanr(kk.log2_ratio, kk.base)
print(f"Variant, sessions with at least 1000 ripples in the first epoch: n = {len(kk)}, rho = {rr.statistic:.3f}, p = {rr.pvalue:.2g}")
print("Variant, base rate measured outside ripples (not the same spikes as the change), reading (a):",
stats.spearmanr(k.log2_ratio, k.out).statistic.round(3))
fig, ax = plt.subplots(figsize=(4, 3.5))
ax.scatter(np.log10(k.base), k.log2_ratio, s=4, alpha=0.5)
ax.set_xlabel("log10 baseline rate (Hz)")
ax.set_ylabel("log2 rate in / out of ripples")
plt.show()
print("Range plotted: x", np.round([np.log10(k.base).min(), np.log10(k.base).max()], 2), "y", np.round([k.log2_ratio.min(), k.log2_ratio.max()], 2))
Sessions: 24 ; putative pyramidal cells: 1192 ; cells with no spikes in ripples (left out of reading a): 25 (a) log2 ratio in/out: n = 1166, rho = -0.437, p = 1.8e-55, within -0.19 +- 0.068: False (b) difference in - out: n = 1192, rho = -0.002, p = 0.95, within -0.19 +- 0.068: False Per subject and per session, reading (a): fCamk1: n = 369, rho = -0.45, p = 1.8e-19 fCamk2: n = 52, rho = -0.71, p = 2.9e-09 fCamk3: n = 540, rho = -0.42, p = 2.6e-24 fCamk5: n = 205, rho = -0.37, p = 4.9e-08 fCamk1-200827-sess9: n = 51, rho = -0.52, p = 0.0001 fCamk1-200901-sess12: n = 40, rho = -0.43, p = 0.0062 fCamk1-200902-sess13: n = 32, rho = -0.28, p = 0.12 fCamk1-200904-sess15: n = 53, rho = -0.59, p = 3e-06 fCamk1-200908-sess16: n = 54, rho = -0.39, p = 0.0037 fCamk1-200909-sess17: n = 53, rho = -0.46, p = 0.00059 fCamk1-200910-sess18: n = 50, rho = -0.49, p = 0.00028 fCamk1-200911-sess19: n = 36, rho = -0.40, p = 0.017 fCamk2-201012-sess1: n = 52, rho = -0.71, p = 2.9e-09 fCamk3-201028-sess10-cleanned: n = 37, rho = -0.52, p = 0.00092 fCamk3-201029-sess11-cleanned: n = 36, rho = -0.46, p = 0.005 fCamk3-201030-sess12: n = 42, rho = -0.75, p = 9.7e-09 fCamk3-201102-sess13: n = 70, rho = -0.13, p = 0.29 fCamk3-201103-sess14: n = 40, rho = -0.39, p = 0.012 fCamk3-201105-sess16: n = 43, rho = -0.35, p = 0.02 fCamk3-201106-sess17: n = 41, rho = -0.55, p = 0.00019 fCamk3-201109-sess18: n = 62, rho = -0.62, p = 6.9e-08 fCamk3-201110-sess19: n = 56, rho = -0.41, p = 0.0017 fCamk3-201111-sess20: n = 59, rho = -0.46, p = 0.00027 fCamk3-201113-sess22: n = 54, rho = -0.32, p = 0.019 fCamk5-210407-sess11: n = 55, rho = -0.56, p = 9.2e-06 fCamk5-210408-sess12: n = 55, rho = -0.35, p = 0.0092 fCamk5-210412-sess14: n = 45, rho = -0.04, p = 0.81 fCamk5-210414-sess16: n = 50, rho = -0.34, p = 0.016 Variant, without the sessions whose ripples fail the clock check ['fCamk3-201102-sess13', 'fCamk3-201106-sess17', 'fCamk5-210406-sess10']: n = 1055, rho = -0.468, p = 1.9e-58 Variant, sessions whose first epoch is labelled BaselinePre (22 sessions): n = 1044, rho = -0.450, p = 2.9e-53 Variant, sessions with at least 1000 ripples in the first epoch: n = 512, rho = -0.453, p = 2.8e-27 Variant, base rate measured outside ripples (not the same spikes as the change), reading (a): -0.466
Range plotted: x [-2.26 0.69] y [-5.39 5.09]
Result: partially reproduced. Over 1166 putative pyramidal cells of 24 sessions, the log2 ratio of the rate inside ripples to the rate outside is negatively correlated with the baseline rate, rho = -0.437 (p = 1.8e-55). It is negative in every subject (-0.37 to -0.71) and every session, and significant (p < 0.05) in all sessions but three (fCamk1 sess13, p = 0.12; fCamk3 sess13, p = 0.29; fCamk5 sess14, p = 0.81). Over all cells it is stronger than the paper's -0.19 plus or minus 0.068. Under the other reading, the difference inside minus outside, rho = -0.002 (p = 0.95). Without the two sessions in this analysis whose ripples fail the clock check (fCamk3 sess13 and sess17; the third, fCamk5 sess10, has under 50 ripples in its first epoch), n = 1055 and rho = -0.468; in the 22 sessions whose first epoch is labelled BaselinePre, rho = -0.450; against the rate outside ripples, which shares no spikes with the rate inside, rho = -0.466. Which measure of change the paper used is not known, and the two readings disagree.
Claim fig2b-theta-trough-gain¶
During track running, the light-induced rate gain is largest at the trough of the theta cycle (P < 1e-6, Wilcoxon test).
rows, check, lfp_checks = [], [], []
for s in opto:
b0 = v.epoch_interval(s, "BaselinePre")[0]
t0 = b0[0] + 600
lfp = v.read_lfp(s["full_url"], t0, t0 + 300)
peaks = s["ripples"][:, 2]
peaks = peaks[(peaks > t0) & (peaks < t0 + 300)]
ch = v.ripple_channel(lfp, t0, peaks)
mz = v.epoch_interval(s, "Maze")[0]
x = v.read_lfp(s["full_url"], mz[0], mz[1], [ch])[:, 0]
phase = v.theta_phase(x)
# clock check of the LFP: ripple-band power at ripple peaks against shifted peaks
from scipy import signal as sg
sos = sg.butter(3, v.RIPPLE_BAND, btype="band", fs=v.LFP_RATE, output="sos")
y = sg.sosfiltfilt(sos, lfp[:, ch].astype(float)) ** 2
ii = ((peaks - t0) * v.LFP_RATE).astype(int)
ii = ii[(ii > 2000) & (ii < len(y) - 2000)]
check.append((s["key"].split("_ses-")[1], ch, len(ii), y[ii].mean() / y.mean(), y[ii + 1250].mean() / y.mean()))
# the same check inside the Maze epoch, on the channel used for theta
ym = sg.sosfiltfilt(sos, x.astype(float)) ** 2
pm = s["ripples"][:, 2]
im = ((pm[(pm > mz[0]) & (pm < mz[1])] - mz[0]) * v.LFP_RATE).astype(int)
im = im[(im > 2000) & (im < len(ym) - 2000)]
# theta-band signal, for the phase convention and for light artifacts
th_sos = sg.butter(3, v.THETA_BAND, btype="band", fs=v.LFP_RATE, output="sos")
yt = sg.sosfiltfilt(th_sos, x.astype(float))
lfp_checks.append({"session": s["key"].split("_ses-")[1], "ripples_in_maze": len(im),
"maze_power_at_peaks_x_mean": round(ym[im].mean() / ym.mean(), 1),
"maze_power_1s_later_x_mean": round(ym[im + 1250].mean() / ym.mean(), 1),
"yt": yt, "mz0": mz[0], "x": x})
def ph(t):
i = ((t - mz[0]) * v.LFP_RATE).astype(int)
return phase[np.clip(i, 0, len(phase) - 1)]
laps = s["laps"][:, :2]
ctrl = s["control"][v.in_intervals(s["control"], laps)]
c_ph = ph(ctrl + 0.01)
for u in np.flatnonzero(s["responsive"] & s["pyr"]):
on = s["cell_onsets"][u]
on = on[v.in_intervals(on, laps)]
l_ph = ph(on + 0.01)
xs = s["spikes"][u]
lc, cc = v.count_in_windows(xs, on, 0.02) / 0.02, v.count_in_windows(xs, ctrl, 0.02) / 0.02
row = {"session": s["key"].split("_ses-")[1], "unit": u}
trough_l, trough_c = np.abs(l_ph) > np.pi / 2, np.abs(c_ph) > np.pi / 2
row["gain_trough"] = lc[trough_l].mean() - cc[trough_c].mean()
row["gain_peak"] = lc[~trough_l].mean() - cc[~trough_c].mean()
edges = np.linspace(-np.pi, np.pi, 7)
for k in range(6):
a = (l_ph >= edges[k]) & (l_ph < edges[k + 1])
b = (c_ph >= edges[k]) & (c_ph < edges[k + 1])
row[f"bin{k}"] = lc[a].mean() - cc[b].mean()
rows.append(row)
for c in check:
print(f"{c[0]}: LFP channel {c[1]}; ripple-band power at {c[2]} ripple peaks {c[3]:.1f} x the mean, "
f"1 s after them {c[4]:.1f} x")
# The stored timing of the LFP, against the raw signal of the same file
for s in opto:
dd = v.describe_file(s["full_url"])
print(f"{s['key'].split('_ses-')[1]}: LFP rate {dd['lfp']['rate']} Hz, starting time {dd['lfp']['starting_time']} s, "
f"{dd['lfp']['samples']} samples = {dd['lfp']['samples'] / dd['lfp']['rate']:.2f} s; raw rate {dd['raw']['rate']} Hz, "
f"{dd['raw']['samples'] / dd['raw']['rate']:.2f} s; last epoch ends at {s['epochs'][-1]['stop']:.2f} s")
print(pd.DataFrame([{k: c[k] for k in ("session", "ripples_in_maze", "maze_power_at_peaks_x_mean", "maze_power_1s_later_x_mean")}
for c in lfp_checks]).to_string(index=False))
# Phase convention: theta phase of spikes of putative pyramidal cells in control windows during
# laps (CA1 pyramidal cells fire most near the trough of pyramidal-layer theta); and light
# artifacts: the theta-band LFP averaged around deposited onsets during laps.
for s, c in zip(opto, lfp_checks):
laps = s["laps"][:, :2]
ctrl = s["control"][v.in_intervals(s["control"], laps)]
civ = v.merge_intervals(np.column_stack([ctrl, ctrl + v.CONTROL_WINDOW_S]))
spk = np.concatenate([x_[v.in_intervals(x_, civ)] for x_, p_ in zip(s["spikes"], s["pyr"]) if p_])
ph_s = np.angle(sg.hilbert(c["yt"]))[np.clip(((spk - c["mz0"]) * v.LFP_RATE).astype(int), 0, len(c["yt"]) - 1)]
hist = np.histogram(ph_s, np.linspace(-np.pi, np.pi, 7))[0]
mean_ph = np.degrees(np.angle(np.mean(np.exp(1j * ph_s))))
on = np.sort(np.concatenate([p["on"] for p in s["pulses"].values()]))
on = on[v.in_intervals(on, laps)]
io = ((on - c["mz0"]) * v.LFP_RATE).astype(int)
io = io[(io > 250) & (io < len(c["yt"]) - 250)]
seg = np.arange(-125, 250)
pta = c["yt"][io[:, None] + seg].mean(axis=0)
raw_pta = c["x"][io[:, None] + seg].astype(float).mean(axis=0)
raw_pta -= raw_pta[:125].mean()
print(f"{c['session']}: control-window spikes of putative pyramidal cells during laps {len(spk)}; counts in 60-degree bins "
f"from -180 (trough) through 0 (peak): {hist.tolist()}; circular mean phase {mean_ph:.0f} degrees")
print(f" {len(io)} deposited onsets during laps: theta-band LFP averaged around them, largest absolute value "
f"{np.abs(pta).max():.1f} against a theta-band SD of {c['yt'].std():.1f} (raw units); broadband average from "
f"0 to 40 ms after onset, largest deviation {np.abs(raw_pta[125:175]).max():.1f}")
th = pd.DataFrame(rows)
d = th.gain_trough - th.gain_peak
w = stats.wilcoxon(th.gain_trough, th.gain_peak)
print(f"\nCells: {len(th)}; gain trough half median {th.gain_trough.median():.2f} Hz, peak half {th.gain_peak.median():.2f} Hz; "
f"trough higher in {(d > 0).sum()} of {len(d)}; median difference {d.median():.2f} Hz; Wilcoxon p = {w.pvalue:.2g}")
for sess, grp in th.groupby("session"):
dd = grp.gain_trough - grp.gain_peak
print(f" {sess}: n = {len(grp)}, median difference {dd.median():.2f}, p = {stats.wilcoxon(grp.gain_trough, grp.gain_peak).pvalue:.2g}")
centres = np.degrees((np.linspace(-np.pi, np.pi, 7)[:-1] + np.pi / 6)).round()
binmed = th[[f"bin{k}" for k in range(6)]].median()
print("Variant, median gain (Hz) in 60-degree bins centred at", centres.tolist(), "(0 = peak, +-180 = trough):",
binmed.round(2).tolist(), "; largest at", centres[int(np.argmax(binmed.values))])
fCamk1-200827-sess9: LFP channel 26; ripple-band power at 25 ripple peaks 193.7 x the mean, 1 s after them 0.3 x fCamk1-200901-sess12: LFP channel 21; ripple-band power at 28 ripple peaks 248.3 x the mean, 1 s after them 0.4 x
fCamk1-200827-sess9: LFP rate 1250.0 Hz, starting time 0.0 s, 22400160 samples = 17920.13 s; raw rate 30000.0 Hz, 17920.13 s; last epoch ends at 17920.13 s
fCamk1-200901-sess12: LFP rate 1250.0 Hz, starting time 0.0 s, 18512640 samples = 14810.11 s; raw rate 30000.0 Hz, 14810.11 s; last epoch ends at 14810.11 s
session ripples_in_maze maze_power_at_peaks_x_mean maze_power_1s_later_x_mean
fCamk1-200827-sess9 44 112.9 0.9
fCamk1-200901-sess12 27 78.0 6.1
fCamk1-200827-sess9: control-window spikes of putative pyramidal cells during laps 18565; counts in 60-degree bins from -180 (trough) through 0 (peak): [3199, 2786, 2638, 3125, 3431, 3386]; circular mean phase 128 degrees
10547 deposited onsets during laps: theta-band LFP averaged around them, largest absolute value 27.1 against a theta-band SD of 405.6 (raw units); broadband average from 0 to 40 ms after onset, largest deviation 216.3
fCamk1-200901-sess12: control-window spikes of putative pyramidal cells during laps 9193; counts in 60-degree bins from -180 (trough) through 0 (peak): [1523, 1531, 1439, 1493, 1657, 1550]; circular mean phase 130 degrees
8013 deposited onsets during laps: theta-band LFP averaged around them, largest absolute value 8.5 against a theta-band SD of 770.1 (raw units); broadband average from 0 to 40 ms after onset, largest deviation 130.3
Cells: 57; gain trough half median 2.04 Hz, peak half 2.07 Hz; trough higher in 22 of 57; median difference -0.09 Hz; Wilcoxon p = 0.14
fCamk1-200827-sess9: n = 35, median difference -0.07, p = 0.46
fCamk1-200901-sess12: n = 22, median difference -0.14, p = 0.16
Variant, median gain (Hz) in 60-degree bins centred at [-150.0, -90.0, -30.0, 30.0, 90.0, 150.0] (0 = peak, +-180 = trough): [1.88, 1.96, 1.6, 2.13, 2.02, 1.92] ; largest at 30.0
Result: not reproduced. In 57 light-responsive putative pyramidal cells during laps, the gain (light minus control) is 2.04 Hz in the trough half of theta and 2.07 Hz in the peak half (medians). The trough is higher in 22 of 57 cells (Wilcoxon p = 0.14), and in neither session (p = 0.46 and 0.16). In 60-degree bins the median gain is largest at 30 degrees past the peak (2.13 Hz), and 1.88 and 1.92 Hz in the bins next to the trough. The checks above support the phase: the stored LFP rate is 1250 Hz from 0 s, ripple-band power at 44 and 27 deposited ripple peaks inside the Maze epoch is 113 and 78 times its mean, and spikes of putative pyramidal cells in control windows have a circular mean phase of 128 and 130 degrees, in the trough half, though weakly modulated (bin counts differ by at most 30%). Averaged around deposited onsets, the theta-band LFP deviates by at most 27.1 and 8.5 against a theta-band SD of 406 and 770, so light artifacts in the theta band are small. The LFP channel, band and phase convention are this attempt's choices, since the Methods were not available.
# Rate maps and place fields over the laps of the Maze epoch, per running direction.
for s in opto:
grid = v.lap_grid(s, s["laps"])
s["grid"] = grid
on_any = v.in_intervals(grid["t"], s["light_iv"])
# control time: windows dark under the protocol; all time without a deposited pulse is a variant
s["ctrl_mask_nolight"] = ~on_any
s["ctrl_mask"] = v.in_intervals(grid["t"], s["dark_iv"]) & ~on_any
s["ctrl_maps"] = [v.rate_map(x, grid, s["ctrl_mask"]) for x in s["spikes"]]
s["ctrl_fields"] = [[v.find_fields(m[d]) for d in (0, 1)] for m in s["ctrl_maps"]]
s["place"] = np.array([any(len(f) for f in ff) for ff in s["ctrl_fields"]])
# light masks: the cell's own uLEDs, and every tested uLED
def mask_for(onsets):
iv = np.column_stack([onsets, onsets + v.PULSE_WINDOW_S])
return v.in_intervals(grid["t"], v.merge_intervals(iv))
s["light_mask"] = {u: mask_for(on) for u, on in s["cell_onsets"].items()}
s["all_light_mask"] = mask_for(s["all_tested_onsets"])
s["light_maps"] = {u: v.rate_map(s["spikes"][u], grid, m) for u, m in s["light_mask"].items()}
s["all_light_maps"] = [v.rate_map(x, grid, s["all_light_mask"]) for x in s["spikes"]]
print(f"{s['key'].split('_ses-')[1]}: grid {len(grid['t'])} steps of {v.GRID_S * 1000:.0f} ms over {len(s['laps'])} laps; "
f"control time (dark windows) {s['ctrl_mask'].mean():.1%}; time without a deposited pulse {s['ctrl_mask_nolight'].mean():.1%}; time in the light windows of all tested uLEDs {s['all_light_mask'].mean():.1%}; "
f"place cells among putative pyramidal {int((s['place'] & s['pyr']).sum())} of {int(s['pyr'].sum())}")
fCamk1-200827-sess9: grid 860909 steps of 2 ms over 111 laps; control time (dark windows) 19.7%; time without a deposited pulse 87.9%; time in the light windows of all tested uLEDs 12.1%; place cells among putative pyramidal 38 of 51
fCamk1-200901-sess12: grid 745747 steps of 2 ms over 132 laps; control time (dark windows) 19.2%; time without a deposited pulse 89.5%; time in the light windows of all tested uLEDs 10.7%; place cells among putative pyramidal 34 of 40
Claim figs11-place-cell-fraction¶
More than half of the pyramidal neurons were place cells, 73% of light-responsive and 71% of non-responsive neurons, with no difference between them (P = 0.90, chi-squared test).
rows = []
for s in opto:
for u in np.flatnonzero(s["pyr"]):
rows.append({"session": s["key"].split("_ses-")[1], "unit": u, "responsive": bool(s["responsive"][u]),
"place": bool(s["place"][u]), "n_fields": sum(len(f) for f in s["ctrl_fields"][u])})
pc = pd.DataFrame(rows)
tab = pd.crosstab(pc.responsive, pc.place)
print(tab)
frac = pc.groupby("responsive").place.mean() * 100
print("Place cells among responsive: %.1f%% (95%% interval %s), among non-responsive: %.1f%% (%s)" % (
frac[True], ci95(pc[pc.responsive].place.sum(), pc.responsive.sum()), frac[False],
ci95(pc[~pc.responsive].place.sum(), (~pc.responsive).sum())))
chi = stats.chi2_contingency(tab.values)
print(f"Chi-squared (Yates) = {chi.statistic:.3f}, p = {chi.pvalue:.3f}; without correction p = {stats.chi2_contingency(tab.values, correction=False).pvalue:.3f}")
print("Per session:")
print(pc.groupby(["session", "responsive"]).place.agg(["mean", "size"]).round(3).to_string())
print("Place cells over all putative pyramidal cells: %.1f%%; fields per place cell, median %d" %
(100 * pc.place.mean(), pc[pc.place].n_fields.median()))
# Threshold variant: minimum peak rate of a field
for peak in (1.0, 3.0):
old = v.FIELD_MIN_PEAK_HZ
v.FIELD_MIN_PEAK_HZ = peak
pl = np.concatenate([[any(len(v.find_fields(s["ctrl_maps"][u][d])) for d in (0, 1)) for u in np.flatnonzero(s["pyr"])] for s in opto])
v.FIELD_MIN_PEAK_HZ = old
resp = pc.responsive.values
print(f"Variant, field peak at least {peak} Hz: place cells {100 * pl[resp].mean():.1f}% of responsive, {100 * pl[~resp].mean():.1f}% of non-responsive")# Variants: the field rule used before verification (no sub-field or dip rule), and all time
# without a deposited pulse as control time
resp = pc.responsive.values
var = []
for mask_name in ("ctrl_mask", "ctrl_mask_nolight"):
maps = {(k, u): (s["ctrl_maps"][u] if mask_name == "ctrl_mask" else v.rate_map(s["spikes"][u], s["grid"], s[mask_name]))
for k, s in enumerate(opto) for u in np.flatnonzero(s["pyr"])}
for rule_name, rule in (("buzcode port", v.find_fields), ("earlier rule", v.find_fields_simple)):
pl = np.array([any(len(rule(m[d])) for d in (0, 1)) for m in maps.values()])
tab_v = pd.crosstab(resp, pl).reindex(index=[False, True], columns=[False, True], fill_value=0).values
var.append({"control time": mask_name, "field rule": rule_name, "place % responsive": round(100 * pl[resp].mean(), 1),
"place % non-responsive": round(100 * pl[~resp].mean(), 1), "place % all": round(100 * pl.mean(), 1),
"chi2 p": round(stats.chi2_contingency(tab_v).pvalue, 4),
"cells changed vs main": int((pl != pc.place.values).sum())})
print(pd.DataFrame(var).to_string(index=False))
place False True
responsive
False 10 24
True 9 48
Place cells among responsive: 84.2% (95% interval 72.1 to 92.5%), among non-responsive: 70.6% (52.5 to 84.9%)
Chi-squared (Yates) = 1.639, p = 0.201; without correction p = 0.122
Per session:
mean size
session responsive
fCamk1-200827-sess9 False 0.562 16
True 0.829 35
fCamk1-200901-sess12 False 0.833 18
True 0.864 22
Place cells over all putative pyramidal cells: 79.1%; fields per place cell, median 2
Variant, field peak at least 1.0 Hz: place cells 91.2% of responsive, 79.4% of non-responsive
Variant, field peak at least 3.0 Hz: place cells 73.7% of responsive, 61.8% of non-responsive
control time field rule place % responsive place % non-responsive place % all chi2 p cells changed vs main
ctrl_mask buzcode port 84.2 70.6 79.1 0.2005 0
ctrl_mask earlier rule 84.2 67.6 78.0 0.1131 1
ctrl_mask_nolight buzcode port 82.5 64.7 75.8 0.0969 5
ctrl_mask_nolight earlier rule 78.9 52.9 69.2 0.0180 11
Result: partially reproduced. 48 of 57 light-responsive putative pyramidal cells (84.2%, 95% interval 72.1 to 92.5%) and 24 of 34 non-responsive ones (70.6%, 52.5 to 84.9%) are place cells. The paper's chi-squared test finds no difference (p = 0.20), as in the paper. The non-responsive fraction is within tolerance of 71%; the responsive one is above 73 plus or minus 7.3%. In sess9, 83% of responsive and 56% of non-responsive cells are place cells; in sess12, 86% and 83%. Field peak thresholds of 1 and 3 Hz give 91.2% against 79.4%, and 73.7% against 61.8%. Over all putative pyramidal cells, 79.1% are place cells, more than half as stated. With all time without a deposited pulse as control time and the earlier field rule, the fractions are 78.9% and 52.9% (p = 0.018); the control time alone changes 5 cells, the field rule alone 1, and both together 11. Some non-responsive cells may respond to uLEDs whose pulses are not in the files (fig1e).
Claim fig2e-infield-gain¶
The light-induced rate increase is higher within than outside the place field (n = 553 place fields, P < 1e-115, Wilcoxon test).
def rate_where(x, grid, mask):
# rate (Hz) of spikes falling in grid steps where mask holds
i = np.searchsorted(grid["t"], x, side="right") - 1
ok = (i >= 0) & (i < len(grid["t"]))
i = i[ok][(x[ok] - grid["t"][i[ok]]) < v.GRID_S]
t = mask.sum() * v.GRID_S
return mask[i].sum() / t if t > 0 else np.nan, t
rows = []
for s in opto:
grid = s["grid"]
b = v.bin_index(grid["x"])
for u in np.flatnonzero(s["responsive"] & s["pyr"] & s["place"]):
for d in (0, 1):
for (a, z, peak) in s["ctrl_fields"][u][d]:
infield = (grid["dir"] == d) & (b >= a) & (b < z)
outfield = (grid["dir"] == d) & ~((b >= a) & (b < z))
li, lt = rate_where(s["spikes"][u], grid, s["light_mask"][u] & infield)
ci, _ = rate_where(s["spikes"][u], grid, s["ctrl_mask"] & infield)
lo, _ = rate_where(s["spikes"][u], grid, s["light_mask"][u] & outfield)
co, _ = rate_where(s["spikes"][u], grid, s["ctrl_mask"] & outfield)
rows.append({"session": s["key"].split("_ses-")[1], "unit": u, "dir": d, "start": a, "stop": z,
"light_in": li, "ctrl_in": ci, "light_out": lo, "ctrl_out": co, "light_s_in_field": lt,
"home_rate": s["base_rate"][u]})
inf = pd.DataFrame(rows)
inf["gain_in"] = inf.light_in - inf.ctrl_in
inf["gain_out"] = inf.light_out - inf.ctrl_out
print("Place fields of responsive putative pyramidal cells:", len(inf), "from", inf.groupby(["session", "unit"]).ngroups,
"cells; light time in field per field, median", round(inf.light_s_in_field.median(), 1), "s")
print(inf[["light_in", "ctrl_in", "light_out", "ctrl_out", "gain_in", "gain_out"]].median().round(2).to_string())
d = inf.gain_in - inf.gain_out
w = stats.wilcoxon(inf.gain_in, inf.gain_out)
print(f"In-field gain higher in {(d > 0).sum()} of {len(d)}; median difference {d.median():.2f} Hz; Wilcoxon p = {w.pvalue:.2g}")
for sess, grp in inf.groupby("session"):
dd = grp.gain_in - grp.gain_out
print(f" {sess}: n = {len(grp)}, median difference {dd.median():.2f}, p = {stats.wilcoxon(grp.gain_in, grp.gain_out).pvalue:.2g}")
ratio_in = inf.light_in / inf.ctrl_in.replace(0, np.nan)
ratio_out = inf.light_out / inf.ctrl_out.replace(0, np.nan)
k = np.isfinite(ratio_in) & np.isfinite(ratio_out)
print(f"Variant, ratio light/control: in field median {ratio_in[k].median():.2f}, out {ratio_out[k].median():.2f}, "
f"n = {k.sum()}, p = {stats.wilcoxon(ratio_in[k], ratio_out[k]).pvalue:.2g}")
Place fields of responsive putative pyramidal cells: 80 from 48 cells; light time in field per field, median 12.0 s light_in 7.17 ctrl_in 2.86 light_out 2.11 ctrl_out 0.85 gain_in 3.95 gain_out 1.37 In-field gain higher in 61 of 80; median difference 2.42 Hz; Wilcoxon p = 1e-08 fCamk1-200827-sess9: n = 51, median difference 2.78, p = 9.5e-07 fCamk1-200901-sess12: n = 29, median difference 1.70, p = 0.0018 Variant, ratio light/control: in field median 2.19, out 2.77, n = 79, p = 4.3e-05
Result: reproduced. Over 80 control place fields of 48 light-responsive putative pyramidal cells, the gain is 3.95 Hz in the field and 1.37 Hz outside it (medians). In-field gain is higher for 61 of 80 fields (Wilcoxon p = 1e-08), and in both sessions (p = 9.5e-07 and 0.0018). As a ratio of light to control rate, the in-field ratio is lower (2.19 against 2.77, p = 4.3e-05), so the increase is larger in hertz but not relative to the control rate.
Claim figs12-gain-rate-correlations¶
The in-field gain is positively correlated with the out-of-field firing rate (rho = 0.17, P < 1e-5) and with the home-cage firing rate of the neuron (rho = 0.25, P < 1e-7).
r1 = stats.spearmanr(inf.gain_in, inf.ctrl_out)
r2 = stats.spearmanr(inf.gain_in, inf.home_rate)
print(f"n = {len(inf)} fields. In-field gain against out-of-field rate: rho = {r1.statistic:.3f}, p = {r1.pvalue:.2g}, "
f"within 0.17 +- 0.084: {abs(r1.statistic - 0.17) <= 0.084}")
print(f"In-field gain against home-cage rate: rho = {r2.statistic:.3f}, p = {r2.pvalue:.2g}, "
f"within 0.25 +- 0.084: {abs(r2.statistic - 0.25) <= 0.084}")
for sess, grp in inf.groupby("session"):
print(f" {sess}: n = {len(grp)}, rho out-of-field {stats.spearmanr(grp.gain_in, grp.ctrl_out).statistic:.2f}, "
f"rho home cage {stats.spearmanr(grp.gain_in, grp.home_rate).statistic:.2f}")
one = inf.sort_values("gain_in", ascending=False).groupby(["session", "unit"]).first()
print(f"Variant, one field per cell (n = {len(one)}): rho out-of-field {stats.spearmanr(one.gain_in, one.ctrl_out).statistic:.3f}, "
f"rho home cage {stats.spearmanr(one.gain_in, one.home_rate).statistic:.3f}")
n = 80 fields. In-field gain against out-of-field rate: rho = -0.018, p = 0.87, within 0.17 +- 0.084: False In-field gain against home-cage rate: rho = -0.063, p = 0.58, within 0.25 +- 0.084: False fCamk1-200827-sess9: n = 51, rho out-of-field 0.03, rho home cage -0.12 fCamk1-200901-sess12: n = 29, rho out-of-field -0.13, rho home cage 0.06 Variant, one field per cell (n = 48): rho out-of-field 0.131, rho home cage -0.003
Result: not reproduced. Over 80 fields, the in-field gain is correlated with neither the out-of-field rate (rho = -0.018, p = 0.87) nor the home-cage rate (rho = -0.063, p = 0.58); both are outside tolerance of 0.17 and 0.25. With one field per cell (48 cells), rho is 0.131 for the out-of-field rate, within tolerance, and -0.003 for the home-cage rate.
Claim fig3a-place-cells-more-excitable¶
Light responses in the home cage before the track were stronger in place cells than in non-place cells (P < 1e-4, Kruskal-Wallis test).
rows = []
for s in opto:
t = s["resp_table"][s["resp_table"].tested].copy()
t["resp_hz"] = t.rate_during - t.rate_before
best = t.groupby("unit").agg(resp_hz=("resp_hz", "max"), z=("z", "max"))
for u in np.flatnonzero(s["pyr"]):
rows.append({"session": s["key"].split("_ses-")[1], "unit": u, "place": bool(s["place"][u]),
"resp_hz": best.resp_hz.get(u, np.nan), "z": best.z.get(u, np.nan)})
ex = pd.DataFrame(rows).dropna()
for col in ("resp_hz", "z"):
a, b = ex[ex.place][col], ex[~ex.place][col]
kw = stats.kruskal(a, b)
print(f"{col}: place cells n = {len(a)}, median {a.median():.2f}; non-place n = {len(b)}, median {b.median():.2f}; "
f"Kruskal-Wallis p = {kw.pvalue:.2g}")
for sess, grp in ex.groupby("session"):
a, b = grp[grp.place].resp_hz, grp[~grp.place].resp_hz
print(f" {sess}: place {len(a)} median {a.median():.2f}, non-place {len(b)} median {b.median():.2f}, p = {stats.kruskal(a, b).pvalue:.2g}")
# firing rate as a confound: the comparison within baseline-rate halves
ex["base"] = [s["base_rate"][u] for sess, u in zip(ex.session, ex.unit) for s in opto if s["key"].endswith(sess)]
med = ex.base.median()
for name, part in (("below", ex[ex.base < med]), ("above", ex[ex.base >= med])):
a, b = part[part.place].resp_hz, part[~part.place].resp_hz
print(f" baseline rate {name} the median: place {len(a)} median {a.median():.2f}, non-place {len(b)} median {b.median():.2f}, "
f"p = {stats.kruskal(a, b).pvalue:.2g}")
resp_hz: place cells n = 72, median 1.90; non-place n = 19, median 1.03; Kruskal-Wallis p = 0.19 z: place cells n = 72, median 3.48; non-place n = 19, median 1.97; Kruskal-Wallis p = 0.25 fCamk1-200827-sess9: place 38 median 2.29, non-place 13 median 1.03, p = 0.027 fCamk1-200901-sess12: place 34 median 0.94, non-place 6 median 1.19, p = 0.85 baseline rate below the median: place 34 median 1.55, non-place 11 median 1.61, p = 0.6 baseline rate above the median: place 38 median 2.25, non-place 8 median 0.82, p = 0.2
Result: not reproduced. In the home cage (PreStim), the light response (rate during the pulse minus before it, best uLED) has a median of 1.90 Hz in 72 place cells and 1.03 Hz in 19 non-place cells, higher in place cells as stated, but not significantly (Kruskal-Wallis p = 0.19; z-scored, p = 0.25). It is significant in sess9 (2.29 against 1.03 Hz, p = 0.027) and not in sess12 (0.94 against 1.19 Hz, p = 0.85). Within halves of the baseline rate neither difference is significant (p = 0.6 and 0.2).
Claim fig3c-unmasked-fields¶
Optogenetic depolarization revealed place fields in the majority of non-place cells, 289 of 417 (69.3%).
rows = []
for s in opto:
for u in np.flatnonzero(s["pyr"] & ~s["place"]):
r = {"session": s["key"].split("_ses-")[1], "unit": u, "responsive": bool(s["responsive"][u])}
if r["responsive"]:
r["light_fields_a"] = [v.find_fields(s["light_maps"][u][d]) for d in (0, 1)]
r["light_fields_b"] = [v.find_fields(s["all_light_maps"][u][d]) for d in (0, 1)]
rows.append(r)
nonp = pd.DataFrame(rows)
nonp["unmasked_a"] = nonp.apply(lambda r: any(len(f) for f in r.light_fields_a) if r.responsive else False, axis=1)
nonp["unmasked_b"] = nonp.light_fields_b.apply(lambda ff: any(len(f) for f in ff))
ra = nonp[nonp.responsive]
pa = 100 * ra.unmasked_a.mean()
pb = 100 * nonp.unmasked_b.mean()
print(f"(a) light-responsive non-place cells: {int(ra.unmasked_a.sum())} of {len(ra)} with a light field ({pa:.1f}%), "
f"within 69.3 +- 6.9: {abs(pa - 69.3) <= 6.9}; 95% interval {ci95(ra.unmasked_a.sum(), len(ra))}")
print(f"(b) all non-place cells, every tested uLED: {int(nonp.unmasked_b.sum())} of {len(nonp)} ({pb:.1f}%), "
f"within 69.3 +- 6.9: {abs(pb - 69.3) <= 6.9}; 95% interval {ci95(nonp.unmasked_b.sum(), len(nonp))}")
print(nonp.groupby(["session", "responsive"])[["unmasked_a", "unmasked_b"]].agg(["sum", "size"]).to_string())
# Null: the same field rule on control maps made from as little time as the light maps
rng = np.random.default_rng(31)
null = []
for s in opto:
frac = s["all_light_mask"].mean() / s["ctrl_mask"].mean()
for u in np.flatnonzero(s["pyr"] & ~s["place"]):
keep = s["ctrl_mask"] & (rng.random(len(s["ctrl_mask"])) < frac)
m = v.rate_map(s["spikes"][u], s["grid"], keep)
null.append(any(len(v.find_fields(m[d])) for d in (0, 1)))
print(f"Reference, non-place cells with a field in a control map thinned at random to the time of the light maps: "
f"{sum(null)} of {len(null)} ({100 * np.mean(null):.1f}%)")
# Null for reading (a): for each light-responsive non-place cell, control maps thinned at random to
# the time of that cell's own light map (its responsive uLEDs only), 20 draws per cell.
rng = np.random.default_rng(34)
na = []
for s in opto:
for u in np.flatnonzero(s["pyr"] & ~s["place"] & s["responsive"]):
frac = s["light_mask"][u].mean() / s["ctrl_mask"].mean()
hits = []
for _ in range(20):
keep = s["ctrl_mask"] & (rng.random(len(s["ctrl_mask"])) < frac)
m = v.rate_map(s["spikes"][u], s["grid"], keep)
hits.append(any(len(v.find_fields(m[d])) for d in (0, 1)))
na.append({"session": s["key"].split("_ses-")[1], "unit": u, "light_time_pct_of_laps": round(100 * s["light_mask"][u].mean(), 1),
"light_field": bool(ra.set_index(["session", "unit"]).unmasked_a[(s["key"].split("_ses-")[1], u)]),
"pct_thinned_draws_with_field": 100 * np.mean(hits)})
na = pd.DataFrame(na)
print("\nReference for reading (a), per light-responsive non-place cell:")
print(na.to_string(index=False))
print(f"Expected share with a field in thinned control maps: {na.pct_thinned_draws_with_field.mean():.1f}% "
f"(against {pa:.1f}% with a light field)")
(a) light-responsive non-place cells: 7 of 9 with a light field (77.8%), within 69.3 +- 6.9: False; 95% interval 40.0 to 97.2%
(b) all non-place cells, every tested uLED: 6 of 19 (31.6%), within 69.3 +- 6.9: False; 95% interval 12.6 to 56.6%
unmasked_a unmasked_b
sum size sum size
session responsive
fCamk1-200827-sess9 False 0 7 2 7
True 4 6 2 6
fCamk1-200901-sess12 False 0 3 1 3
True 3 3 1 3
Reference, non-place cells with a field in a control map thinned at random to the time of the light maps: 3 of 19 (15.8%)
Reference for reading (a), per light-responsive non-place cell:
session unit light_time_pct_of_laps light_field pct_thinned_draws_with_field
fCamk1-200827-sess9 12 6.2 True 40.0
fCamk1-200827-sess9 19 3.1 True 60.0
fCamk1-200827-sess9 23 3.1 True 0.0
fCamk1-200827-sess9 24 3.1 True 15.0
fCamk1-200827-sess9 27 3.1 False 25.0
fCamk1-200827-sess9 38 3.1 False 0.0
fCamk1-200901-sess12 13 2.6 True 55.0
fCamk1-200901-sess12 22 2.6 True 60.0
fCamk1-200901-sess12 27 2.6 True 30.0
Expected share with a field in thinned control maps: 31.7% (against 77.8% with a light field)
Result: not reproduced. Under reading (a), 7 of 9 light-responsive non-place cells (77.8%, 95% interval 40.0 to 97.2%) have a light field; under reading (b), 6 of 19 non-place cells (31.6%, 12.6 to 56.6%). Neither is within 69.3 plus or minus 6.9%. Reading (a) misses by 1.6 points, and its interval, from 9 cells, includes the paper's value. Reading (b) is above the 15.8% (3 of 19) that the same rule finds in control maps thinned at random to the time of the light maps. The light maps of reading (a) use only each cell's own uLEDs, 2.6% to 6.2% of lap time against 10.7% and 12.1% for reading (b). Control maps thinned to each cell's own light time give a field in 31.7% of draws on average over the 9 cells, against 77.8% with a light field. The two sessions hold 19 non-place cells, against the paper's 417.
# Field groups for figures 3E to 3G: control fields of light-responsive place cells, and light
# fields (reading a) of light-responsive non-place cells, each with its direction.
groups = []
for s in opto:
occ_ctrl = {d: np.bincount(v.bin_index(s["grid"]["x"][(s["grid"]["dir"] == d) & s["ctrl_mask"]]), minlength=v.N_BINS) for d in (0, 1)}
for u in np.flatnonzero(s["responsive"] & s["pyr"]):
lm, cm = s["light_maps"][u], s["ctrl_maps"][u]
occ_l = {d: np.bincount(v.bin_index(s["grid"]["x"][(s["grid"]["dir"] == d) & s["light_mask"][u]]), minlength=v.N_BINS) for d in (0, 1)}
fields = s["ctrl_fields"][u] if s["place"][u] else [v.find_fields(lm[d]) for d in (0, 1)]
for d in (0, 1):
if not fields[d]:
continue
groups.append({"session": s["key"].split("_ses-")[1], "unit": u, "dir": d,
"group": "place" if s["place"][u] else "unmasked",
"peak_ctrl": int(np.nanargmax(cm[d])) if np.any(cm[d] > 0) else np.nan,
"peak_light": int(np.nanargmax(lm[d])) if np.any(lm[d] > 0) else np.nan,
"si_ctrl": v.spatial_info(cm[d], occ_ctrl[d]), "si_light": v.spatial_info(lm[d], occ_l[d])})
fg = pd.DataFrame(groups)
print(fg.groupby("group").size().to_string())
# Variant control: per cell, a random draw of control windows (20 ms before a pulse, no uLED on)
# as many as the cell's light windows over the laps, as the main text's "random intervals (20 ms)
# between the light pulses" may mean. Used as a variant in figures 3F and 3G.
rng = np.random.default_rng(30)
for s in opto:
laps = s["laps"][:, :2]
ctrl = s["control"][v.in_intervals(s["control"], laps)]
s["ctrl_matched_mask"] = {}
for u, on in s["cell_onsets"].items():
n = int(v.in_intervals(on, laps).sum())
pick = np.sort(rng.choice(ctrl, size=min(n, len(ctrl)), replace=False))
iv = v.merge_intervals(np.column_stack([pick, pick + v.CONTROL_WINDOW_S]))
s["ctrl_matched_mask"][u] = v.in_intervals(s["grid"]["t"], iv)
print("Matched control time as a fraction of light time, median over cells:",
round(np.median([s["ctrl_matched_mask"][u].sum() / s["light_mask"][u].sum() for s in opto for u in s["ctrl_matched_mask"]]), 3))
group place 70 unmasked 14
Matched control time as a fraction of light time, median over cells: 1.0
Claim fig3e-ghost-fields¶
The spatial location of light-induced spikes matches that of spikes without stimulation, both for place fields of place cells (P < 1e-40) and for unmasked place fields of non-place cells, whose sparse control spikes form ghost fields (P < 1e-10).
rng = np.random.default_rng(32)
for grp, k in fg.dropna(subset=["peak_ctrl", "peak_light"]).groupby("group"):
r = stats.pearsonr(k.peak_ctrl, k.peak_light)
same = np.mean(np.abs(k.peak_ctrl - k.peak_light) <= 2)
sh = np.array([np.mean(np.abs(k.peak_ctrl.values - rng.permutation(k.peak_light.values)) <= 2) for _ in range(500)])
print(f"{grp}: n = {len(k)} field directions, Pearson r = {r.statistic:.3f}, p = {r.pvalue:.2g}; "
f"peaks within 2 bins (4.4 cm): {same:.1%}, shuffled pairings {sh.mean():.1%} (95th percentile {np.percentile(sh, 95):.1%}), "
f"shuffle p = {(np.sum(sh >= same) + 1) / 501:.3f}")
for sess, kk in k.groupby("session"):
if len(kk) > 3:
print(f" {sess}: n = {len(kk)}, r = {stats.pearsonr(kk.peak_ctrl, kk.peak_light).statistic:.2f}")
print("Field directions dropped for an empty map:", int(fg[["peak_ctrl", "peak_light"]].isna().any(axis=1).sum()))
fig, axs = plt.subplots(1, 2, figsize=(7, 3.3))
for ax, (grp, k) in zip(axs, fg.dropna(subset=["peak_ctrl", "peak_light"]).groupby("group")):
ax.scatter(k.peak_ctrl * v.TRACK_CM / v.N_BINS, k.peak_light * v.TRACK_CM / v.N_BINS, s=10)
ax.set_title(grp)
ax.set_xlabel("control peak (cm)")
ax.set_ylabel("light peak (cm)")
ax.set_xlim(0, 110)
ax.set_ylim(0, 110)
plt.tight_layout()
plt.show()
place: n = 70 field directions, Pearson r = 0.806, p = 3.9e-17; peaks within 2 bins (4.4 cm): 65.7%, shuffled pairings 12.2% (95th percentile 18.6%), shuffle p = 0.002
fCamk1-200827-sess9: n = 44, r = 0.87
fCamk1-200901-sess12: n = 26, r = 0.71
unmasked: n = 14 field directions, Pearson r = 0.321, p = 0.26; peaks within 2 bins (4.4 cm): 35.7%, shuffled pairings 14.6% (95th percentile 28.6%), shuffle p = 0.014
fCamk1-200827-sess9: n = 8, r = 0.43
fCamk1-200901-sess12: n = 6, r = 0.35
Field directions dropped for an empty map: 0
Result: partially reproduced. For 70 field directions of light-responsive place cells (control fields), the light map's peak follows the control map's peak (Pearson r = 0.806, p = 3.9e-17). For 14 field directions of the 7 non-place cells with a light field, r = 0.321 (p = 0.26), so the criterion's test fails for this group. In the shuffle variant, 35.7% of these peaks are within 4.4 cm of the control peak, against 14.6% for shuffled pairings (p = 0.014). Of the 14 unmasked field directions, 8 come from sess9.
Claim fig3f-light-stability¶
Place fields during light stimulation were more stable between the first and second halves of the session than during control epochs, for both true place fields and unmasked place fields (P < 1e-9 and P < 1e-19, two-way ANOVA).
rows = []
for s in opto:
grid = s["grid"]
lap_no = grid["lap"] + 1
h1 = (lap_no >= 11) & (lap_no <= 60)
h2 = (lap_no >= 71) & (lap_no <= 120)
on_all = np.sort(np.concatenate([p["on"] for p in s["pulses"].values()]))
stim_lap = np.array([((on_all >= a) & (on_all < b)).sum() / (b - a) >= 1 for a, b, _ in s["laps"]])
h1s, h2s = h1 & stim_lap[grid["lap"]], h2 & stim_lap[grid["lap"]]
sub = fg[fg.session == s["key"].split("_ses-")[1]]
for _, r in sub.iterrows():
u, d = r.unit, r.dir
lm1 = v.rate_map(s["spikes"][u], grid, s["light_mask"][u] & h1)[d]
lm2 = v.rate_map(s["spikes"][u], grid, s["light_mask"][u] & h2)[d]
cm1 = v.rate_map(s["spikes"][u], grid, s["ctrl_mask"] & h1)[d]
cm2 = v.rate_map(s["spikes"][u], grid, s["ctrl_mask"] & h2)[d]
mm1 = v.rate_map(s["spikes"][u], grid, s["ctrl_matched_mask"][u] & h1)[d]
mm2 = v.rate_map(s["spikes"][u], grid, s["ctrl_matched_mask"][u] & h2)[d]
cs1 = v.rate_map(s["spikes"][u], grid, s["ctrl_mask"] & h1s)[d]
cs2 = v.rate_map(s["spikes"][u], grid, s["ctrl_mask"] & h2s)[d]
rows.append({"session": r.session, "unit": u, "dir": d, "group": r.group,
"stab_light": v.map_corr(lm1, lm2), "stab_ctrl": v.map_corr(cm1, cm2),
"stab_ctrl_matched": v.map_corr(mm1, mm2), "stab_ctrl_stimlaps": v.map_corr(cs1, cs2)})
st = pd.DataFrame(rows)
print("Laps per session:", {s["key"].split("_ses-")[1]: len(s["laps"]) for s in opto}, "; halves are laps 11-60 and 71-120")
for grp, k in st.groupby("group"):
k2 = k.dropna()
w = stats.wilcoxon(k2.stab_light, k2.stab_ctrl)
dd = k2.stab_light - k2.stab_ctrl
print(f"{grp}: n = {len(k2)} (dropped {len(k) - len(k2)} with a flat or empty half map); stability light median "
f"{k2.stab_light.median():.3f}, control {k2.stab_ctrl.median():.3f}; light higher in {(dd > 0).sum()}; "
f"Wilcoxon p = {w.pvalue:.2g}")
for grp, k in st.groupby("group"):
k2 = k.dropna(subset=["stab_light", "stab_ctrl_matched"])
dd = k2.stab_light - k2.stab_ctrl_matched
print(f"Variant, control from time-matched control windows, {grp}: n = {len(k2)}, light median {k2.stab_light.median():.3f}, "
f"matched control {k2.stab_ctrl_matched.median():.3f}; light higher in {(dd > 0).sum()}; Wilcoxon p = "
f"{stats.wilcoxon(k2.stab_light, k2.stab_ctrl_matched).pvalue:.2g}")
for grp, k in st.groupby("group"):
k2 = k.dropna(subset=["stab_light", "stab_ctrl_stimlaps"])
dd = k2.stab_light - k2.stab_ctrl_stimlaps
print(f"Variant, control maps from stimulated laps only, {grp}: n = {len(k2)}, light median {k2.stab_light.median():.3f}, "
f"control {k2.stab_ctrl_stimlaps.median():.3f}; light higher in {(dd > 0).sum()}; Wilcoxon p = "
f"{stats.wilcoxon(k2.stab_light, k2.stab_ctrl_stimlaps).pvalue:.2g}")
import statsmodels.formula.api as smf
import statsmodels.api as sm
long = pd.concat([st.assign(cond="light", stab=st.stab_light), st.assign(cond="control", stab=st.stab_ctrl)]).dropna(subset=["stab"])
an = sm.stats.anova_lm(smf.ols("stab ~ C(cond) * C(group)", data=long).fit(), typ=2)
print("\nVariant, two-way ANOVA of stability on condition and group:")
print(an.round(4).to_string())
Laps per session: {'fCamk1-200827-sess9': 111, 'fCamk1-200901-sess12': 132} ; halves are laps 11-60 and 71-120
place: n = 70 (dropped 0 with a flat or empty half map); stability light median 0.664, control 0.714; light higher in 28; Wilcoxon p = 0.058
unmasked: n = 14 (dropped 0 with a flat or empty half map); stability light median 0.187, control 0.527; light higher in 6; Wilcoxon p = 0.39
Variant, control from time-matched control windows, place: n = 70, light median 0.664, matched control 0.602; light higher in 35; Wilcoxon p = 0.68
Variant, control from time-matched control windows, unmasked: n = 14, light median 0.187, matched control 0.293; light higher in 6; Wilcoxon p = 0.67
Variant, control maps from stimulated laps only, place: n = 70, light median 0.664, control 0.714; light higher in 28; Wilcoxon p = 0.059
Variant, control maps from stimulated laps only, unmasked: n = 14, light median 0.187, control 0.527; light higher in 6; Wilcoxon p = 0.39
Variant, two-way ANOVA of stability on condition and group:
sum_sq df F PR(>F)
C(cond) 0.3850 1.0 4.0415 0.046
C(group) 2.6881 1.0 28.2216 0.000
C(cond):C(group) 0.0377 1.0 0.3961 0.530
Residual 15.6210 164.0 NaN NaN
Result: not reproduced. Between laps 11-60 and 71-120, light maps are not more stable than control maps. The medians are 0.664 against 0.714 for place fields (light higher in 28 of 70, p = 0.058) and 0.187 against 0.527 for unmasked fields (light higher in 6 of 14, p = 0.39). Against control maps from as many 20 ms control windows as light windows, light is higher for place fields (0.664 against 0.602) but not significantly (p = 0.68), and lower for unmasked fields (0.187 against 0.293, p = 0.67). Leaving the laps without stimulation out of the control maps gives the same medians (p = 0.059 and 0.39); this variant cannot differ much, since dark windows exist only next to deposited pulses, so laps without stimulation hold almost no control time. The two-way ANOVA finds an effect of condition (p = 0.046), with light the less stable in both groups by the medians above.
Claim fig3g-spatial-information¶
Spatial information of control and light-stimulated epochs is correlated, and the light-boosted increase is stronger for unmasked place fields of non-place cells than for place cells (P < 1e-15, Kruskal-Wallis test).
k = fg.dropna(subset=["si_ctrl", "si_light"]).copy()
k["boost"] = k.si_light - k.si_ctrl
for grp, kk in k.groupby("group"):
r = stats.spearmanr(kk.si_ctrl, kk.si_light)
print(f"{grp}: n = {len(kk)}, SI control median {kk.si_ctrl.median():.2f}, light {kk.si_light.median():.2f} bits/spike; "
f"Spearman rho = {r.statistic:.3f}, p = {r.pvalue:.2g}; increase median {kk.boost.median():.2f}")
kw = stats.kruskal(k[k.group == "unmasked"].boost, k[k.group == "place"].boost)
print(f"Increase, unmasked against place: Kruskal-Wallis p = {kw.pvalue:.2g}")
print("Field directions dropped for an empty map:", len(fg) - len(k))
si_m = []
for _, r in k.iterrows():
s = [x for x in opto if x["key"].endswith(r.session)][0]
g_ = s["grid"]
mask = s["ctrl_matched_mask"][r.unit]
occ = np.bincount(v.bin_index(g_["x"][(g_["dir"] == r.dir) & mask]), minlength=v.N_BINS)
si_m.append(v.spatial_info(v.rate_map(s["spikes"][r.unit], g_, mask)[r.dir], occ))
k["si_ctrl_matched"] = si_m
k["boost_matched"] = k.si_light - k.si_ctrl_matched
for grp, kk in k.dropna(subset=["si_ctrl_matched"]).groupby("group"):
r = stats.spearmanr(kk.si_ctrl_matched, kk.si_light)
print(f"Variant, time-matched control, {grp}: n = {len(kk)}, SI matched control median {kk.si_ctrl_matched.median():.2f}, "
f"rho = {r.statistic:.3f}, p = {r.pvalue:.2g}; increase median {kk.boost_matched.median():.2f}")
km = k.dropna(subset=["si_ctrl_matched"])
print("Variant, increase over time-matched control, unmasked against place: Kruskal-Wallis p = %.2g" %
stats.kruskal(km[km.group == "unmasked"].boost_matched, km[km.group == "place"].boost_matched).pvalue)
k["ratio"] = k.si_light / k.si_ctrl
print("Variant, ratio light/control: unmasked median %.2f, place median %.2f, p = %.2g" % (
k[k.group == "unmasked"].ratio.median(), k[k.group == "place"].ratio.median(),
stats.kruskal(k[k.group == "unmasked"].ratio, k[k.group == "place"].ratio).pvalue))
place: n = 70, SI control median 0.28, light 0.26 bits/spike; Spearman rho = 0.655, p = 7.6e-10; increase median -0.03 unmasked: n = 14, SI control median 0.22, light 0.25 bits/spike; Spearman rho = 0.574, p = 0.032; increase median 0.03 Increase, unmasked against place: Kruskal-Wallis p = 0.039 Field directions dropped for an empty map: 0
Variant, time-matched control, place: n = 70, SI matched control median 0.36, rho = 0.501, p = 1e-05; increase median -0.10 Variant, time-matched control, unmasked: n = 14, SI matched control median 0.46, rho = 0.609, p = 0.021; increase median -0.17 Variant, increase over time-matched control, unmasked against place: Kruskal-Wallis p = 0.3 Variant, ratio light/control: unmasked median 1.08, place median 0.90, p = 0.14
Result: partially reproduced. Control and light spatial information are correlated for place fields (rho = 0.655, p = 7.6e-10) and for unmasked fields (rho = 0.574, p = 0.032). The increase from control to light is larger for unmasked fields (median 0.03 bits/spike) than for place fields (-0.03; Kruskal-Wallis p = 0.039), far from the paper's P < 1e-15. Against time-matched control maps the increase is negative in both groups (unmasked -0.17, place -0.10) and does not differ (p = 0.3); as a ratio it does not differ either (1.08 against 0.90, p = 0.14). Unmasked fields are chosen for having a light field and no control field, which favours a larger increase in that group by construction; the paper defines its groups the same way. The paper's control epochs are "random intervals (20 ms) between the light pulses", which the time-matched variant follows most closely, and the dark-window control maps are themselves a revision of the criterion's control (all time with no uLED on). Under the time-matched control both correlations hold (place rho = 0.501, p = 1e-05; unmasked rho = 0.609, p = 0.021) but the increase does not differ between groups, so part 1 holds and part 2 does not. The outcome was changed from reproduced after a third verification round agreed.
Claim fig4c-swr-cofiring-place-pairs¶
Cofiring of place cell pairs during SPW-Rs in the home cage is positively correlated with their spatial overlap on the track (P < 1e-27), while pairs of a place cell and a non-place cell show no such relationship (P = 0.11).
def ripple_counts(s):
home = [(e["start"], e["stop"]) for e in s["epochs"] if e["behavioral_paradigm"] in ("BaselinePre", "BaselinePost")]
r = s["ripples"]
keep = np.zeros(len(r), bool)
for a, b in home:
keep |= (r[:, 0] >= a) & (r[:, 1] < b)
r = r[keep]
return np.array([np.searchsorted(x, r[:, 1]) - np.searchsorted(x, r[:, 0]) for x in s["spikes"]]), len(r)
pairs = []
for s in opto:
counts, nrip = ripple_counts(s)
cof = np.corrcoef(counts)
s["cofiring"] = cof
pyr = np.flatnonzero(s["pyr"])
cmaps = [np.concatenate(s["ctrl_maps"][u]) for u in range(len(s["spikes"]))]
lmaps = {u: np.concatenate(m) for u, m in s["light_maps"].items()}
for i_, a in enumerate(pyr):
for b in pyr[i_ + 1:]:
kind = {2: "PC-PC", 1: "PC-nonPC", 0: "nonPC-nonPC"}[int(s["place"][a]) + int(s["place"][b])]
pairs.append({"session": s["key"].split("_ses-")[1], "a": a, "b": b, "kind": kind, "cofiring": cof[a, b],
"overlap_ctrl": v.map_corr(cmaps[a], cmaps[b]),
"overlap_light": v.map_corr(lmaps[a], lmaps[b]) if (a in lmaps and b in lmaps) else np.nan,
"both_responsive": bool(s["responsive"][a] and s["responsive"][b])})
print(s["key"].split("_ses-")[1], "home-cage ripples without stimulation:", nrip)
pr = pd.DataFrame(pairs)
for kind, k in pr.dropna(subset=["cofiring", "overlap_ctrl"]).groupby("kind"):
r = stats.spearmanr(k.cofiring, k.overlap_ctrl)
print(f"{kind}: n = {len(k)} pairs, Spearman rho = {r.statistic:.3f}, p = {r.pvalue:.2g}")
for sess, kk in k.groupby("session"):
rr = stats.spearmanr(kk.cofiring, kk.overlap_ctrl)
print(f" {sess}: n = {len(kk)}, rho = {rr.statistic:.3f}, p = {rr.pvalue:.2g}")
# Variant: cofiring as the correlation of participation (any spike) per ripple
for s in opto:
counts, _ = ripple_counts(s)
s["cofiring_bin"] = np.corrcoef((counts > 0).astype(float))
pr["cofiring_bin"] = [[x for x in opto if x["key"].endswith(r.session)][0]["cofiring_bin"][r.a, r.b] for r in pr.itertuples()]
for kind, k in pr.dropna(subset=["cofiring_bin", "overlap_ctrl"]).groupby("kind"):
r = stats.spearmanr(k.cofiring_bin, k.overlap_ctrl)
print(f"Variant, participation correlation, {kind}: n = {len(k)}, rho = {r.statistic:.3f}, p = {r.pvalue:.2g}")
print("Pairs dropped (a cell without spikes in ripples or a flat map):", int(pr[["cofiring", "overlap_ctrl"]].isna().any(axis=1).sum()))
fCamk1-200827-sess9 home-cage ripples without stimulation: 3131
fCamk1-200901-sess12 home-cage ripples without stimulation: 1120
PC-PC: n = 1264 pairs, Spearman rho = 0.060, p = 0.032
fCamk1-200827-sess9: n = 703, rho = 0.119, p = 0.0015
fCamk1-200901-sess12: n = 561, rho = 0.019, p = 0.66
PC-nonPC: n = 698 pairs, Spearman rho = 0.062, p = 0.1
fCamk1-200827-sess9: n = 494, rho = 0.034, p = 0.45
fCamk1-200901-sess12: n = 204, rho = 0.068, p = 0.34
nonPC-nonPC: n = 93 pairs, Spearman rho = 0.094, p = 0.37
fCamk1-200827-sess9: n = 78, rho = 0.113, p = 0.33
fCamk1-200901-sess12: n = 15, rho = -0.121, p = 0.67
Variant, participation correlation, PC-PC: n = 1264, rho = 0.054, p = 0.056
Variant, participation correlation, PC-nonPC: n = 698, rho = 0.072, p = 0.058
Variant, participation correlation, nonPC-nonPC: n = 93, rho = 0.199, p = 0.056
Pairs dropped (a cell without spikes in ripples or a flat map): 0
Result: reproduced. Over 1264 place cell pairs, cofiring during home-cage ripples is correlated with spatial overlap (Spearman rho = 0.060, p = 0.032), weakly and in sess9 only (rho = 0.119, p = 0.0015; sess12 rho = 0.019, p = 0.66). Over 698 pairs of a place cell and a non-place cell, the paper's test finds no correlation (rho = 0.062, p = 0.1), as in the paper, but the coefficient is as large as for place cell pairs, so the absence comes from fewer pairs, not a weaker relation. Participation per ripple gives rho = 0.054 (p = 0.056) and 0.072 (p = 0.058). The cofiring measure is this attempt's choice; the paper's is in the supplement.
Claim fig4d-light-unmasked-cofiring¶
With light-induced spikes, the relationship between SPW-R cofiring and spatial overlap appears for pairs that include non-place cells (P < 1e-5, Spearman correlation).
k = pr[(pr.kind == "PC-nonPC") & pr.both_responsive].dropna(subset=["cofiring", "overlap_light"])
r = stats.spearmanr(k.cofiring, k.overlap_light)
print(f"PC-nonPC pairs, both light-responsive, light maps: n = {len(k)}, Spearman rho = {r.statistic:.3f}, p = {r.pvalue:.2g}")
kc = pr[(pr.kind == "PC-nonPC") & pr.both_responsive].dropna(subset=["cofiring", "overlap_ctrl"])
rc = stats.spearmanr(kc.cofiring, kc.overlap_ctrl)
print(f"The same pairs with control maps: n = {len(kc)}, rho = {rc.statistic:.3f}, p = {rc.pvalue:.2g}")
for kind in ("PC-PC", "nonPC-nonPC"):
kk = pr[(pr.kind == kind) & pr.both_responsive].dropna(subset=["cofiring", "overlap_light"])
if len(kk) > 5:
rr = stats.spearmanr(kk.cofiring, kk.overlap_light)
print(f"Variant, {kind} pairs, both responsive, light maps: n = {len(kk)}, rho = {rr.statistic:.3f}, p = {rr.pvalue:.2g}")
for sess, kk in k.groupby("session"):
if len(kk) > 5:
rr = stats.spearmanr(kk.cofiring, kk.overlap_light)
print(f" {sess}: n = {len(kk)}, rho = {rr.statistic:.3f}, p = {rr.pvalue:.2g}")
print(f"\nRun time so far {(time.time() - started) / 60:.1f} min")
PC-nonPC pairs, both light-responsive, light maps: n = 231, Spearman rho = 0.075, p = 0.25 The same pairs with control maps: n = 231, rho = 0.110, p = 0.094 Variant, PC-PC pairs, both responsive, light maps: n = 577, rho = 0.097, p = 0.02 Variant, nonPC-nonPC pairs, both responsive, light maps: n = 18, rho = -0.152, p = 0.55 fCamk1-200827-sess9: n = 174, rho = 0.108, p = 0.16 fCamk1-200901-sess12: n = 57, rho = -0.062, p = 0.65 Run time so far 6.7 min
Result: not reproduced. Over 231 pairs of a light-responsive place cell and a light-responsive non-place cell, cofiring during ripples is not significantly correlated with the overlap of their light maps (rho = 0.075, p = 0.25; sess9 rho = 0.108, p = 0.16; sess12 rho = -0.062, p = 0.65). With control maps the same pairs give rho = 0.110 (p = 0.094). Pairs of two responsive place cells give rho = 0.097 (p = 0.02) with light maps.
Blocked claims: what the files hold¶
The stored value of every pulse series is printed in the data check (one value per uLED), as are the units table's columns. Below, a proxy for figs9-excitatory-vs-inhibitory-swr, shown for the verifier and not deciding its outcome.
# Proxy for figs9: units at or above 5 Hz in the first epoch as putative interneurons, against
# the putative pyramidal cells, by the log2 ratio of the rate inside deposited ripples to the
# rate outside (first epoch, sessions with at least 50 ripples there), as in figs10.
rows = []
for s in sessions:
e0 = s["epochs"][0]
r = s["ripples"]
r = r[(r[:, 0] >= e0["start"]) & (r[:, 1] < e0["stop"])]
if len(r) < 50:
continue
t_in = (r[:, 1] - r[:, 0]).sum()
t_out = (e0["stop"] - e0["start"]) - t_in
for u, x in enumerate(s["spikes"]):
n_all = ((x >= e0["start"]) & (x < e0["stop"])).sum()
n_in = (np.searchsorted(x, r[:, 1]) - np.searchsorted(x, r[:, 0])).sum()
if n_in > 0 and n_all > n_in:
rows.append({"pyr": bool(s["pyr"][u]), "log2_ratio": np.log2((n_in / t_in) / ((n_all - n_in) / t_out))})
px = pd.DataFrame(rows)
a, b = px[px.pyr].log2_ratio, px[~px.pyr].log2_ratio
print(f"Putative pyramidal (below 5 Hz): n = {len(a)}, median log2 ratio {a.median():.2f}; units at or above 5 Hz: "
f"n = {len(b)}, median {b.median():.2f}; Mann-Whitney p = {stats.mannwhitneyu(a, b).pvalue:.2g}")
Putative pyramidal (below 5 Hz): n = 1166, median log2 ratio 1.19; units at or above 5 Hz: n = 121, median 0.75; Mann-Whitney p = 0.00065
Units at or above 5 Hz raise their rate inside ripples less than the putative pyramidal cells (median log2 ratio 0.75 against 1.19, p = 0.00065), in the direction of the claim. The split is by baseline rate, and within the putative pyramidal cells alone the ripple gain falls with baseline rate (figs10, rho = -0.437), so a split by rate would show this with no interneurons among the faster units. The claim stays blocked.
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 |
|---|---|
| cohort-counts | partially reproduced |
| fig1c-pulse-protocol | partially reproduced |
| fig1e-responsive-fraction | reproduced |
| fig1i-swr-decrease | reproduced |
| figs10-swr-gain-vs-baseline-rate | partially reproduced |
| fig2b-theta-trough-gain | not reproduced |
| figs11-place-cell-fraction | partially reproduced |
| fig2e-infield-gain | reproduced |
| figs12-gain-rate-correlations | not reproduced |
| fig3a-place-cells-more-excitable | not reproduced |
| fig3c-unmasked-fields | not reproduced |
| fig3e-ghost-fields | partially reproduced |
| fig3f-light-stability | not reproduced |
| fig3g-spatial-information | partially reproduced |
| fig4c-swr-cofiring-place-pairs | reproduced |
| fig4d-light-unmasked-cofiring | not reproduced |
| fig1j-light-intensity | blocked: data missing |
| fig1km-intracellular | blocked: data missing |
| figs9-excitatory-vs-inhibitory-swr | blocked: metadata insufficient |
| fig3h-position-decoding | blocked: methods underspecified |
| fig4ef-assemblies | blocked: methods underspecified |
| fig4gh-replay | blocked: methods underspecified |
How the analysis differs from the paper's. The Materials and Methods were not available. The light-response test is the one in the first author's later code, the place fields use a port of buzcode's bz_findPlaceFields1D with its defaults, and the theta phase, cofiring and stability measures are common choices. Pyramidal cells are units with a baseline rate below 5 Hz, as the files carry no cell type. Every claim on light responses or place fields rests on 2 sessions of one mouse (91 putative pyramidal cells), the only sessions with pulse times, against the paper's 822 cells. Most pulses are not in the files, so control time is restricted to the windows that the stated protocol keeps dark (the 20 ms before each deposited onset and after each offset); all time without a deposited pulse is shown as a variant. Cells on the shank whose uLEDs are missing cannot be tested for light responses. The dark windows assume every offset between pulses is at least 20 ms; 3% to 4% of the deposited offsets are shorter.
What was found about the dandiset.
- Pulse times are in 2 of 28 sessions; 20 more sessions name uLED stimulation in their epochs but hold no pulse times, so their light responses cannot be computed. Laps are also in these 2 only.
- The pulse series are stated in watts but hold one arbitrary value per uLED (19 in one session, about 63000 in the other). They carry no light intensity, so the intensity results cannot be tested. The files hold series for 6 of the probe's 12 uLEDs, and the deposited pulses are about a third of those the protocol implies (5.76 and 7.50 per s); no cell on one shank responds to any deposited uLED.
- Every session is deposited twice, as a full file and a no-raw-data file with the same units, spikes, ripples and position. One full file lacks the session in its name. The full file of fCamk3 sess17 holds 378 s of raw signal and 6850 s of LFP of a 34621 s session, and six fCamk1 full files have no raw signal.
- Units carry no cell type and the small files no waveforms.
session_start_timeis 2016-03-15 in every file, though the sessions are from 2020 and 2021.- The linearized position has a conversion factor of 100 with values already in cm.
- Six sessions have epochs without paradigm or manipulation. 14 sessions include CNO or DMSO epochs, which the paper does not describe. The dandiset holds 28 sessions and 1351 putative pyramidal cells, where the paper reports 822 cells (about 19 sessions), and does not say which sessions the paper used.
- In three sessions the deposited ripple intervals mark population bursts weakly; fCamk2 has 0 to 1110 ripples per session and no position.
- The file description of the pulses gives 40 to 60 ms offsets, the paper 20 to 40 ms.
- The paper reports 417 non-place cells among 822 pyramidal neurons (51%), while 73% of responsive and 71% of non-responsive neurons are place cells, which leaves about 230 non-place cells.