Senzai, Fernandez-Ruiz and Buzsaki (2019) from dandiset 000166¶
Paper: Senzai Y, Fernandez-Ruiz A, Buzsaki G (2019). Layer-specific physiological features and interlaminar interactions in the primary visual cortex of the mouse. Neuron 101, 500-513. doi:10.1016/j.neuron.2018.12.009
Dataset: DANDI:000166, version 0.250624.0434 (doi:10.48324/dandi.000166/0.250624.0434), streamed from the archive. The paper cites no version. The dandiset has two published versions, 0.220116.2037 and 0.250624.0434, and the cell below shows that they hold the same 19 assets (same ids, paths and sizes); the most recent published one is pinned, as the protocol asks when the paper cites none.
Claims: claims/000166/10.1016-j.neuron.2018.12.009.yaml, written from the version of record (PMC6367010) with its two supplements.
The dataset holds one session from each of 19 mice: spike-sorted units (spike times, a quality label and, in 18 files, cell type and synaptic effect labels), the raw 20 kHz signal and the 1250 Hz LFP of the 64 sites, brain state intervals (18 files), UP and DOWN intervals, and laser pulse times. It holds no unit waveforms, no depth of the sites, and no layer, landmark or cluster label. The data columns are in depth order (checked below), so the laminar analyses are rebuilt here from the LFP and the spikes: gamma ICA, coherence clustering, spike-IC coupling and the K-means unit clusters. Depth is taken in sites relative to landmark c (the MUA power peak, found in every mouse), as the CSD landmarks a, b, d and e cannot be placed by a rule the paper gives (fig1g).
| Claim | Status here |
|---|---|
| fig4a-unit-counts, fig4a-ccg-classification, fig4de-optotagging | tested |
| fig7-connection-counts, fig7c-divergence-convergence, figs8f-ie-dip-wake | tested |
| figs7-strongest-3-6hz-coupling | tested |
| fig2c-six-gamma-ics, ica-volume-conduction-fraction, fig2d-superficial-ic-power | tested (gamma ICA of waking, depth relative to c) |
| fig2b-six-coherence-clusters, fig2b-clusters-same-across-states | tested (coherence clustering, depth relative to c) |
| fig3c-slow-oscillation-nrem-deep, fig3c-3-6hz-layer4, slow-gamma-superficial-wake | tested (layers from coherence clusters or IC bands) |
| figs2c-3-6hz-theta-anticorrelation, fig3f-units-lock-3-6hz-trough | tested (layer 4 site from the ICs) |
| fig5b-eight-clusters, fig5a-spike-ic-depth-stratification | tested |
| fig6a to fig6h (8 claims), fig7-within-layer-fractions, fig7eg-l23-targets, figs8b-reciprocal-connections, fig7d, fig7e, fig7h, fig7g | tested on the eight unit clusters |
| fig8-down-active-neurons, fig8-down-up-rates, fig8cd-down-active-suppressed, fig8eh-down-active-inverse-correlation | tested |
| fig1g-landmark-distances-consistent | blocked (methods underspecified); every reading is run below |
| figs8c-lognormal-strength, fig3ab-alpha-analogy | blocked (methods underspecified, as the criterion says) |
| fig1f-sink-b-deeper, fig1f-other-landmarks-same-depth, figs2d-visual-offset-rebound | blocked (data missing: visual stimulus times in 2 of the 13 mice) |
| fig3d-slow-oscillation-gamma-coupling, fig3d-3-6hz-cross-frequency | tested (phase sites from c and the ICs) |
import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
import time
import numpy as np
import pandas as pd
import requests
import matplotlib.pyplot as plt
from scipy import stats
from concurrent.futures import ProcessPoolExecutor
import v1layers as v
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
print("dandiset", v.DANDISET_ID, "version", v.VERSION)
versions = requests.get(f"{v.API}/dandisets/{v.DANDISET_ID}/versions/", timeout=60).json()["results"]
ids = {}
for ver in versions:
a = requests.get(f"{v.API}/dandisets/{v.DANDISET_ID}/versions/{ver['version']}/assets/?page_size=200",
timeout=60).json()["results"]
ids[ver["version"]] = sorted((x["asset_id"], x["path"], x["size"]) for x in a)
print(ver["version"], ver["status"], ver["asset_count"], "assets", round(ver["size"] / 1e9, 1), "GB")
published = [k for k in ids if k != "draft"]
print("published versions hold the same assets (ids, paths and sizes):",
all(ids[k] == ids[published[0]] for k in published), "; draft the same:", ids["draft"] == ids[v.VERSION])
dandiset 000166 version 0.250624.0434
draft Published 19 assets 787.2 GB 0.220116.2037 Valid 19 assets 787.2 GB
0.250624.0434 Valid 19 assets 787.2 GB published versions hold the same assets (ids, paths and sizes): True ; draft the same: True
assets = v.list_assets()
t0 = time.time()
with ProcessPoolExecutor(4) as ex:
S = list(ex.map(v.process_session, assets))
print(f"{len(S)} files streamed and reduced in {time.time() - t0:.0f} s")
19 files streamed and reduced in 1275 s
# One row per unit and one per detected connection, then what was loaded.
rows = []
for si, r in enumerate(S):
cls, cls_hp = v.classify_waveforms(r["wf"]), v.classify_waveforms(r["wf_hp"])
for k in range(len(r["unit_ids"])):
f, fh = r["wf"][k], r["wf_hp"][k]
o = r["opto"]
rows.append(dict(
sess=si, subject=r["subject"], genotype=r["genotype"], unit=k, unit_id=r["unit_ids"][k],
cell_type=r["cell_type"][k], synaptic_effect=r["synaptic_effect"][k], quality=r["quality"][k],
max_channel=r["max_channel"][k], n_spikes=r["n_spikes"][k], wf_spikes=r["wf_spikes"][k],
wf=cls[k], wf_hp=cls_hp[k], tp_ms=f["tp_ms"] if f else np.nan, tp_hp=fh["tp_ms"] if fh else np.nan,
kurt=f["kurtosis"] if f else np.nan, site=f["site"] if f else -1, amp_uv=f["amp_uv"] if f else np.nan,
ccgE=r["conn"]["exc"][k].any(), ccgI=r["conn"]["inh"][k].any(),
ccgE_h=r["conn_hollow"]["exc"][k].any(), ccgI_h=r["conn_hollow"]["inh"][k].any(),
opto_pooled=o["true_pooled"][k] if o else False, opto_any=o["true_any"][k] if o else False,
opto_pooled_null=o["shifted_pooled"][k] if o else False, opto_any_null=o["shifted_any"][k] if o else False))
U = pd.DataFrame(rows)
crow = []
for si, r in enumerate(S):
p = r["pairs"]
u = U[U.sess == si].reset_index(drop=True)
for j in range(len(p["pre"])):
a, b = p["pre"][j], p["post"][j]
crow.append(dict(sess=si, subject=r["subject"], pre=a, post=b, exc=p["exc"][j], inh=p["inh"][j],
pre_wf=u.wf[a], post_wf=u.wf[b], dsite=abs(u.site[a] - u.site[b]),
stp=p["stp"][j], dip_wake=p["dip_wake"][j], dip_nrem=p["dip_nrem"][j],
stp_wake=p["stp_wake"][j], stp_nrem=p["stp_nrem"][j], ratio_wake=p["ratio_wake"][j], ratio_nrem=p["ratio_nrem"][j]))
C = pd.DataFrame(crow)
print("Files:", len(S), "subjects:", U.subject.nunique())
tab = []
for r in S:
st = r["states"]
dur = {lab: (st["stop"] - st["start"])[st["label"] == lab].sum() / 60 for lab in np.unique(st["label"])} if st else {}
ud = r["updown"]
tab.append(dict(subject=r["subject"], genotype=r["genotype"], session=r["session_id"], units=len(r["unit_ids"]),
good=int((r["quality"] == "good").sum()), spikes=int(r["n_spikes"].sum()),
hours=round(r["duration"] / 3600, 2), wake_min=round(dur.get("Awake", 0)), nrem_min=round(dur.get("Non-REM", 0)),
rem_min=round(dur.get("REM", 0)), ma_min=round(dur.get("MA", 0)),
UP=int((ud["label"] == "UP").sum()), DOWN=int((ud["label"] == "DOWN").sum()),
laser_pulses=len(r["laser"]["start"]) if r["laser"] else 0,
visual_pulses=len(r["visual"]["start"]) if r["visual"] else 0,
has_cell_type="cell_type" in r["columns"]))
T = pd.DataFrame(tab)
print(T.to_string(index=False))
print("\nTotals: units", T.units.sum(), "good", T.good.sum(), "spikes", T.spikes.sum())
print("Genotypes in the files:", T.genotype.value_counts().to_dict())
print("Genotypes in the paper (STAR Methods): 3 Ai35, 5 PV-Cre::Ai32, 9 VGAT-Cre::Ai32, 1 CaMKII-Cre::Ai35, 1 CaMKII-Cre::Ai32")
print("Files without a Sleep states table:", [r["subject"] for r in S if r["states"] is None])
print("Files without unit columns other than quality:", [r["subject"] for r in S if "cell_type" not in r["columns"]])
print("\nUnit quality:", U.quality.value_counts().to_dict())
print("Dataset cell_type labels:", U.cell_type.value_counts().to_dict())
print("Dataset synaptic_effect labels:", U.synaptic_effect.value_counts().to_dict())
print("Units labelled 'down state' by mouse:", U[U.cell_type == "down state"].subject.value_counts().to_dict())
print("\nSampling: raw rate", sorted({r["raw_rate"] for r in S}), "Hz; LFP rate", sorted({r["lfp_rate"] for r in S}),
"Hz; raw sites", sorted({r["raw_shape"][1] for r in S}))
print("Electrode y (depth) all NaN in every file:", all(np.isnan(r["electrode_y"]).all() for r in S),
"; rel_y range", sorted({(r["rel_y"].min(), r["rel_y"].max()) for r in S}))
print("Electrode locations:", sorted({x for r in S for x in r["electrode_location"]}),
"; impedances:", sorted({float(x) for r in S for x in r["electrode_imp"]}))
print("Spike times sorted in every unit:", all(r["sorted_ok"].all() for r in S))
span_ok = [r["t_range"][0] >= 0 and r["t_range"][1] <= r["duration"] + 1 for r in S]
print("Spike times (s) within the raw recording in every file:", all(span_ok),
"; latest spike / recording length:", round(min(r["t_range"][1] / r["duration"] for r in S), 3), "to",
round(max(r["t_range"][1] / r["duration"] for r in S), 3))
rates = U.n_spikes / U.sess.map({i: r["duration"] for i, r in enumerate(S)})
print("Mean firing rates (Hz): median", round(rates.median(), 2), "range", round(rates.min(), 3), "to", round(rates.max(), 1))
same_site = (U.site >= 0)
print("Units whose `max_channel` equals the site of their largest mean waveform (column order of the raw data):",
int((U.max_channel[same_site] == U.site[same_site]).sum()), "of", int(same_site.sum()))
heads = {}
for r in S:
for h in r["train_head"]:
heads.setdefault(h, []).append(r["subject"])
dups = {h: s for h, s in heads.items() if len(set(s)) > 1}
print("Spike trains whose first five spike times appear in more than one file:", len(dups))
print("Total recorded hours per file range:", T.hours.min(), "to", T.hours.max())
for r in S:
if r["visual"] is not None:
vs, vd = r["visual"]["start"], r["visual"]["stop"] - r["visual"]["start"]
lz = r["laser"]["start"]
near = np.min(np.abs(vs[:, None] - lz[None, :]), axis=1) < 0.005 if lz is not None and len(lz) else np.zeros(len(vs), bool)
print(f"{r['subject']} 'Visual laser': {len(vs)} intervals, durations {np.percentile(vd, [0, 50, 100]).round(4)} s;"
f" {near.sum()} start within 5 ms of a 'Laser diode' pulse; intervals of 0.09-0.11 s: {((vd > 0.09) & (vd < 0.11)).sum()},"
f" of 0.39-0.41 s: {((vd > 0.39) & (vd < 0.41)).sum()}")
for r in S:
if r["laser"] is not None and len(r["laser"]["start"]):
d = r["laser"]["stop"] - r["laser"]["start"]
if r["subject"] in ("YMV05", "YMV07"):
print(f"{r['subject']} 'Laser diode': {len(d)} pulses, duration {np.percentile(d, [0, 50, 100]).round(4)} s,"
f" amplitude levels {r['opto']['levels']} with {r['opto']['n_pulses']} pulses")
print("Laser diode pulse counts per level by file:",
{r["subject"]: r["opto"]["n_pulses"].tolist() for r in S if r["opto"] is not None})
print("\nGroups in each file (acquisition, processing, stimulus):")
grp = {}
for r in S:
grp.setdefault(tuple(r["groups"]), []).append(r["subject"])
for g, subj in grp.items():
print(" ", subj, ":", list(g))
Files: 19 subjects: 19
subject genotype session units good spikes hours wake_min nrem_min rem_min ma_min UP DOWN laser_pulses visual_pulses has_cell_type
YMV01 Ai35 YMV01_170818 53 53 1319077 2.05 62 22 2 5 1507 1511 0 0 True
YMV02 Ai35 YMV02_170815 101 101 3861839 1.98 35 55 4 3 3848 3884 0 0 True
YMV03 Ai35 YMV03_170818 116 115 3597835 2.46 19 85 12 9 5912 5948 0 0 False
YMV04 CaMKII-Cre::Ai32 YMV04_170907 63 63 2770879 2.80 56 49 2 2 2607 2638 0 0 True
YMV05 PV-Cre::Ai32 YMV05_170912 72 72 4388284 3.75 6 33 2 4 2436 2460 1488 1988 True
YMV06 VGAT-Cre::Ai32 YMV06_170913 73 73 3869645 3.97 12 123 8 5 10972 11021 1488 0 True
YMV07 CaMKII-Cre::Ai35 YMV07_170914 101 101 14915751 9.02 60 56 8 0 1726 1798 41022 0 True
YMV08 VGAT-Cre::Ai32 YMV08_170922 100 100 5909912 4.18 46 86 12 14 5514 5559 2715 3213 True
YMV09 PV-Cre::Ai32 YMV09_171204 78 78 7753841 6.72 73 185 19 0 16760 16814 1488 0 True
YMV10 PV-Cre::Ai32 YMV10_171213 38 38 3817356 5.27 46 147 4 9 11610 11624 1488 0 True
YMV11 PV-Cre::Ai32 YMV11_171208 70 70 7422894 5.96 129 105 7 0 11406 11506 1488 0 True
YMV12 VGAT-Cre::Ai32 YMV12_171211 60 60 4257226 6.03 153 84 8 7 8413 8465 1488 0 True
YMV13 PV-Cre::Ai32 YMV13_180127 79 79 7327976 6.14 68 170 19 4 16927 17012 1488 0 True
YMV14 PV-Cre::Ai32 YMV14_180128 79 79 7635477 5.52 83 134 16 0 584 620 1984 0 True
YMV15 PV-Cre::Ai32 YMV15_180205 79 79 10064834 6.23 54 181 20 4 16533 16585 1488 0 True
YMV16 PV-Cre::Ai32 YMV16_180206 51 51 7492042 7.52 53 196 20 2 1287 1327 1488 0 True
YMV17 PV-Cre::Ai32 YMV17_180207 81 81 10209728 7.87 63 205 22 16 13686 13729 1488 0 True
YMV18 PV-Cre::Ai32 YMV18_180208 119 119 12909465 7.94 191 151 17 4 13963 14026 1488 0 True
YMV19 PV-Cre::Ai32 YMV19_180209 64 64 6488223 6.58 0 0 0 0 13790 13847 1488 0 True
Totals: units 1477 good 1476 spikes 126012284
Genotypes in the files: {'PV-Cre::Ai32': 11, 'Ai35': 3, 'VGAT-Cre::Ai32': 3, 'CaMKII-Cre::Ai32': 1, 'CaMKII-Cre::Ai35': 1}
Genotypes in the paper (STAR Methods): 3 Ai35, 5 PV-Cre::Ai32, 9 VGAT-Cre::Ai32, 1 CaMKII-Cre::Ai35, 1 CaMKII-Cre::Ai32
Files without a Sleep states table: ['YMV19']
Files without unit columns other than quality: ['YMV03']
Unit quality: {np.str_('good'): 1476, np.str_('unsorted'): 1}
Dataset cell_type labels: {np.str_('Pyramidal Cell'): 836, np.str_('Narrow Interneuron'): 277, np.str_('Wide Interneuron'): 223, np.str_('missing'): 116, np.str_('Unknown'): 15, np.str_('down state'): 10}
Dataset synaptic_effect labels: {np.str_('Unknown'): 772, np.str_('Excitatory'): 444, np.str_('Inhibitory'): 145, np.str_('missing'): 116}
Units labelled 'down state' by mouse: {'YMV07': 4, 'YMV05': 2, 'YMV01': 1, 'YMV08': 1, 'YMV15': 1, 'YMV17': 1}
Sampling: raw rate [20000.0] Hz; LFP rate [1250.0] Hz; raw sites [64]
Electrode y (depth) all NaN in every file: True ; rel_y range [(np.float64(-64.0), np.float64(-1.0))]
Electrode locations: [np.str_('unknown')] ; impedances: [-1.0]
Spike times sorted in every unit: True
Spike times (s) within the raw recording in every file: True ; latest spike / recording length: 1.0 to 1.0
Mean firing rates (Hz): median 2.56 range 0.052 to 59.4
Units whose `max_channel` equals the site of their largest mean waveform (column order of the raw data): 76 of 1367
Spike trains whose first five spike times appear in more than one file: 0
Total recorded hours per file range: 1.98 to 9.02
YMV05 'Visual laser': 1988 intervals, durations [1.00000e-02 1.01000e-02 4.95637e+01] s; 1488 start within 5 ms of a 'Laser diode' pulse; intervals of 0.09-0.11 s: 0, of 0.39-0.41 s: 498
YMV08 'Visual laser': 3213 intervals, durations [0.01 0.01 0.4001] s; 2715 start within 5 ms of a 'Laser diode' pulse; intervals of 0.09-0.11 s: 0, of 0.39-0.41 s: 498
YMV05 'Laser diode': 1488 pulses, duration [0.0099 0.01 0.0101] s, amplitude levels [-1.3 -0.8 1.5] with [496 496 496] pulses
YMV07 'Laser diode': 41022 pulses, duration [0.02 0.02 0.0201] s, amplitude levels [-0.8] with [41022] pulses
Laser diode pulse counts per level by file: {'YMV05': [496, 496, 496], 'YMV06': [496, 496, 496], 'YMV07': [41022], 'YMV08': [992, 992, 731], 'YMV09': [496, 496, 496], 'YMV10': [496, 496, 496], 'YMV11': [496, 496, 496], 'YMV12': [496, 496, 496], 'YMV13': [496, 496, 496], 'YMV14': [496, 992, 496], 'YMV15': [496, 496, 496], 'YMV16': [496, 496, 496], 'YMV17': [496, 496, 496], 'YMV18': [496, 496, 496], 'YMV19': [496, 496, 496]}
Groups in each file (acquisition, processing, stimulus):
['YMV01', 'YMV02', 'YMV03', 'YMV04'] : ['acquisition/ElectricalSeries', 'processing/Neural states', 'processing/ecephys', 'stimulus/presentation', 'stimulus/templates', 'processing/Neural states/Sleep states', 'processing/Neural states/Up-Down states', 'processing/ecephys/LFP']
['YMV05', 'YMV08'] : ['acquisition/ElectricalSeries', 'processing/Neural states', 'processing/ecephys', 'stimulus/presentation', 'stimulus/templates', 'processing/Neural states/Laser diode', 'processing/Neural states/Sleep states', 'processing/Neural states/Up-Down states', 'processing/Neural states/Visual laser', 'processing/ecephys/LFP']
['YMV06', 'YMV07', 'YMV09', 'YMV10', 'YMV11', 'YMV12', 'YMV13', 'YMV14', 'YMV15', 'YMV16', 'YMV17', 'YMV18'] : ['acquisition/ElectricalSeries', 'processing/Neural states', 'processing/ecephys', 'stimulus/presentation', 'stimulus/templates', 'processing/Neural states/Laser diode', 'processing/Neural states/Sleep states', 'processing/Neural states/Up-Down states', 'processing/ecephys/LFP']
['YMV19'] : ['acquisition/ElectricalSeries', 'processing/Neural states', 'processing/ecephys', 'stimulus/presentation', 'stimulus/templates', 'processing/Neural states/Laser diode', 'processing/Neural states/Up-Down states', 'processing/ecephys/LFP']
# Clock checks: spikes against the UP/DOWN table (true and shifted by 10 s), the UP/DOWN table against
# the state table, and spikes against the laser pulses (true and shifted by 0.5 s).
ck = []
for r in S:
ur = r["updown_rates"]
ok = (ur[("UP", 0.0)] > 0)
true = np.nanmedian(ur[("DOWN", 0.0)][ok] / ur[("UP", 0.0)][ok])
sh_ok = ur[("UP", 10.0)] > 0
shifted = np.nanmedian(ur[("DOWN", 10.0)][sh_ok] / ur[("UP", 10.0)][sh_ok])
ud, st = r["updown"], r["states"]
if st is not None:
nm = st["label"] == "Non-REM"
mid = (ud["start"] + ud["stop"]) / 2
inside = np.zeros(len(mid), bool)
for a, b in zip(st["start"][nm], st["stop"][nm]):
inside |= (mid >= a) & (mid < b)
frac = inside.mean()
else:
frac = np.nan
L = r["landmark"]
lfp_sd = L["lfp"].std(axis=0).max() # largest SD over time of the transition-averaged LFP, over sites
lfp_sd_shift = L["lfp_shifted"].std(axis=0).max()
ck.append(dict(subject=r["subject"], down_up_ratio=round(true, 3), down_up_ratio_shifted=round(shifted, 3),
updown_in_nrem=round(frac, 3), lfp_transition_uv=round(lfp_sd, 1),
lfp_shifted_uv=round(lfp_sd_shift, 1)))
CK = pd.DataFrame(ck)
print(CK.to_string(index=False))
print("Median DOWN/UP rate ratio over files: true", CK.down_up_ratio.median(), "shifted 10 s", CK.down_up_ratio_shifted.median())
print("Transition-averaged LFP (largest SD over time, uV): median", CK.lfp_transition_uv.median(),
"; with onsets shifted by 10 s", CK.lfp_shifted_uv.median())
print("Units tagged by the laser pulses (pooled rule): true onsets", int(U.opto_pooled.sum()),
"; onsets shifted by 0.5 s", int(U.opto_pooled_null.sum()))
wk = np.array([np.nanmedian(r["state_rates"]["Awake"]) for r in S if r["state_rates"] is not None])
nr = np.array([np.nanmedian(r["state_rates"]["Non-REM"]) for r in S if r["state_rates"] is not None])
print("Median unit rate per file, Awake against Non-REM (Hz): ", np.round(wk, 2), np.round(nr, 2))
subject down_up_ratio down_up_ratio_shifted updown_in_nrem lfp_transition_uv lfp_shifted_uv YMV01 0.089 1.016 1.0 152.300003 6.8 YMV02 0.163 1.026 1.0 153.399994 6.6 YMV03 0.131 0.992 1.0 194.300003 6.8 YMV04 0.150 1.012 1.0 130.399994 6.1 YMV05 0.072 1.026 1.0 193.800003 7.0 YMV06 0.081 0.994 1.0 207.699997 10.3 YMV07 0.255 0.971 1.0 119.099998 6.3 YMV08 0.072 1.001 1.0 221.100006 8.5 YMV09 0.227 0.999 1.0 163.899994 5.2 YMV10 0.094 1.004 1.0 181.500000 8.9 YMV11 0.090 1.002 1.0 174.899994 7.7 YMV12 0.093 1.006 1.0 159.399994 7.3 YMV13 0.098 1.006 1.0 141.699997 6.6 YMV14 0.101 0.999 1.0 233.600006 30.9 YMV15 0.105 1.002 1.0 187.000000 9.6 YMV16 0.053 1.047 1.0 223.399994 18.9 YMV17 0.162 0.997 1.0 160.000000 7.4 YMV18 0.098 1.002 1.0 188.399994 6.7 YMV19 0.121 1.000 NaN 155.100006 6.6 Median DOWN/UP rate ratio over files: true 0.098 shifted 10 s 1.002 Transition-averaged LFP (largest SD over time, uV): median 174.9 ; with onsets shifted by 10 s 7.0 Units tagged by the laser pulses (pooled rule): true onsets 80 ; onsets shifted by 0.5 s 0 Median unit rate per file, Awake against Non-REM (Hz): [1.72 3.74 2.98 2.83 4.17 3.1 2.57 2.61 3.41 4.3 3.37 2.65 2.82 3.01 3.11 3.54 3.01 2.21] [1.57 2.67 1.94 1.48 2.72 2.13 2.22 2. 2.4 3.01 2.79 2.22 1.83 1.96 2.25 2.51 2.06 1.63]
The checks above are what the data can show about their clocks. Spikes fall far less in the dataset's DOWN intervals than in its UP intervals, and the contrast goes when the intervals are shifted by 10 s, so the UP/DOWN table shares the spike clock. Every UP or DOWN interval lies inside a Non-REM interval of the state table, which is consistent with the state table sharing that clock (in YMV19, which has no state table, this cannot be checked). The LFP averaged around DOWN-UP onsets swings far more than around the same onsets shifted by 10 s, so the LFP shares the clock of the UP/DOWN table and, through it, of the spikes. Units respond to the laser pulses and not to the pulses shifted by 0.5 s, so the pulse table shares the spike clock.
Site order and max_channel¶
Every laminar analysis needs the data columns in depth order. The electrode table gives no depth (x, y, z NaN), its rows run from rel_y -1 to -64, the data point at rows 0 to 63 in order, and its id column is a permutation of 0 to 63. The cell below checks the order on the data (the best-correlated site of each site in the 100-500 Hz LFP should be adjacent) and what the units' max_channel indexes.
# Site order: the electrode table rows run from rel_y -1 to -64 (shank_electrode_number 0-63), the data
# columns point at rows 0-63 in order, and the table's `id` column is a permutation (hardware channels).
# Whether the data columns are in depth order is checked on the data: for each site, the site whose
# 100-500 Hz LFP (200 s) correlates best with it should be an adjacent one if the columns follow depth.
rows = []
for si, r in enumerate(S):
cc = r["landmark"]["neighbour_corr"].copy()
np.fill_diagonal(cc, -1)
best = np.argmax(cc, axis=1)
u = U[(U.sess == si) & (U.site >= 0) & (U.max_channel >= 0)]
ids = r["electrode_id"]
rows.append(dict(subject=r["subject"], best_is_adjacent=int((np.abs(best - np.arange(64)) == 1).sum()),
best_within_2=int((np.abs(best - np.arange(64)) <= 2).sum()),
adjacent_corr=round(float(np.mean([cc[i, i + 1] for i in range(63)])), 3),
units=len(u), maxch_eq_column=int((u.max_channel == u.site).sum()),
maxch_eq_id_of_column=int((u.max_channel.values == ids[u.site.values]).sum()),
id_within_1_site=int((np.abs(np.argsort(ids)[u.max_channel.values] - u.site.values) <= 1).sum())))
CH = pd.DataFrame(rows)
print(CH.to_string(index=False))
print("Electrode table id column (first 16 rows):", S[0]["electrode_id"][:16], "; the same in every file:",
all((r["electrode_id"] == S[0]["electrode_id"]).all() for r in S),
"; data columns point at table rows 0-63 in order in every file:",
all((r["series_electrodes"] == np.arange(64)).all() for r in S))
print("Sites whose best-correlated site is adjacent: median", int(CH.best_is_adjacent.median()), "of 64, range",
CH.best_is_adjacent.min(), "to", CH.best_is_adjacent.max())
print("max_channel equals the table id of the unit's peak column for", int(CH.maxch_eq_id_of_column.sum()), "of",
int(CH.units.sum()), "units (within one site:", int(CH.id_within_1_site.sum()), "); equals the column index for",
int(CH.maxch_eq_column.sum()))
subject best_is_adjacent best_within_2 adjacent_corr units maxch_eq_column maxch_eq_id_of_column id_within_1_site YMV01 62 64 0.804 51 1 48 51 YMV02 61 63 0.833 95 7 92 95 YMV03 59 63 0.826 0 0 0 0 YMV04 50 54 0.865 62 3 56 62 YMV05 59 62 0.906 72 2 66 71 YMV06 63 63 0.880 72 5 66 72 YMV07 55 57 0.817 90 4 84 90 YMV08 58 61 0.847 97 5 90 95 YMV09 61 62 0.855 64 7 59 64 YMV10 55 58 0.899 34 2 27 34 YMV11 55 62 0.889 66 4 59 65 YMV12 51 57 0.906 46 1 38 46 YMV13 56 57 0.795 75 6 69 75 YMV14 56 58 0.814 70 5 65 68 YMV15 56 58 0.796 76 6 72 76 YMV16 58 60 0.810 49 3 44 48 YMV17 56 58 0.782 70 2 66 70 YMV18 59 60 0.823 107 10 101 106 YMV19 58 58 0.727 59 3 55 59 Electrode table id column (first 16 rows): [22 13 15 24 11 20 7 9 21 8 25 12 10 23 14 27] ; the same in every file: True ; data columns point at table rows 0-63 in order in every file: True Sites whose best-correlated site is adjacent: median 58 of 64, range 50 to 63 max_channel equals the table id of the unit's peak column for 1157 of 1255 units (within one site: 1247 ); equals the column index for 76
fig1g-landmark-distances-consistent (blocked, methods_underspecified; every reading shown)¶
The distances between the five depth landmarks (the MUA power peak c and the CSD landmarks a, b, d, e) vary by less than 100 um from mouse to mouse.
Criterion: c is the depth of the peak of 500 Hz to 5 kHz power; a, b, d, e are the DOWN-UP CSD landmarks as in fig1f-sink-b-deeper (the most negative CSD near the layer 3/4 border within 0 to 100 ms, reading a, or 0 to 50 ms, reading b, after the onset; sources as the most positive CSD near the stated depth); "varied" is read as the range across mice of each distance between adjacent landmarks, or each mouse's deviation from the across-mouse mean. Reproduced when, under a reading, every adjacent-landmark distance varies by less than 110 um.
Here the DOWN-UP transitions are the dataset's UP onsets that follow a DOWN interval, averaged over the 1000 s of each session that hold the most of them; MUA power is the 500 Hz to 5 kHz power of the raw signal in the two 30 s waveform blocks. The paper does not give a rule that places a, b, d and e (Results, Anchoring Mesoscopic Patterns, paragraph 1: "the most prominent current source density (CSD) sinks and sources"; Figure 1D marks them by hand with asterisks). The rule used is: b the most negative CSD above c and d below it, a the most positive above b and e below d, each over the whole half of the probe. The claims file applies no spatial smoothing (main); a second difference over two sites (40 um) is shown as a variant. The second cell shows the CSD map of every mouse with the landmarks the rule places.
lm_rows = []
for r in S:
L = r["landmark"]
for step in (1, 2):
cs = v.csd_of(L["lfp"], step)
for w in v.CSD_WINDOWS:
m = v.landmarks(cs, L["lags"], r["mua_power"], w, step)
row = dict(subject=r["subject"], transitions=L["n"], step=step, window=f"{w[0]*1000:.0f}-{w[1]*1000:.0f} ms",
c=int(np.argmax(r["mua_power"])))
if m:
row.update({k: m[k] for k in "abde"})
row.update({"a-b": (m["b"] - m["a"]) * v.SITE_UM, "b-c": (m["c"] - m["b"]) * v.SITE_UM,
"c-d": (m["d"] - m["c"]) * v.SITE_UM, "d-e": (m["e"] - m["d"]) * v.SITE_UM})
lm_rows.append(row)
LM = pd.DataFrame(lm_rows)
print("MUA peak site (c) by mouse:", dict(zip(LM.subject[LM.step == 1][::2], LM.c[LM.step == 1][::2])))
print("Transitions averaged per mouse: median", int(LM.transitions.median()), "range", LM.transitions.min(), "to", LM.transitions.max())
for (step, w), g in LM.groupby(["step", "window"]):
ok = g.dropna(subset=["a-b"])
d = ok[["a-b", "b-c", "c-d", "d-e"]]
rng = (d.max() - d.min())
dev = (d - d.mean()).abs().max()
print(f"\nCSD over +-{step} site(s), window {w}: landmarks placed in {len(ok)} of {len(g)} mice "
f"(not placed in {len(g) - len(ok)}: b at the top site or d at the bottom site, so a or e has no site beyond it)")
print(" distance medians (um):", d.median().round(0).to_dict())
print(" range across mice (um):", rng.round(0).to_dict(), "; all < 110:", bool((rng < 110).all()))
print(" largest deviation from the mean (um):", dev.round(0).to_dict(), "; all < 110:", bool((dev < 110).all()))
print()
print(LM[(LM.step == 1) & (LM.window == "0-100 ms")].drop(columns=["step", "window"]).to_string(index=False))
r0 = S[4]
fig, ax = plt.subplots(1, 2, figsize=(8, 4))
cs = v.csd_of(r0["landmark"]["lfp"], 1)
lim = np.percentile(np.abs(cs), 99)
ax[0].imshow(cs.T, aspect="auto", cmap="RdBu", vmin=-lim, vmax=lim,
extent=[r0["landmark"]["lags"][0] * 1000, r0["landmark"]["lags"][-1] * 1000, 62.5, 0.5])
ax[0].set_xlabel("ms from DOWN-UP"); ax[0].set_ylabel("site"); ax[0].set_title(r0["subject"] + " CSD (sinks red)")
ax[1].plot(r0["mua_power"] / r0["mua_power"].max(), np.arange(64)); ax[1].invert_yaxis()
ax[1].set_xlabel("MUA power (norm.)"); plt.tight_layout(); plt.show()
print(f"CSD colour scale clipped at the 99th percentile of |CSD| ({lim:.1f} uV per site^2); {int((np.abs(cs) > lim).sum())} of {cs.size} values beyond it")
MUA peak site (c) by mouse: {'YMV01': 32, 'YMV02': 31, 'YMV03': 31, 'YMV04': 39, 'YMV05': 38, 'YMV06': 38, 'YMV07': 39, 'YMV08': 35, 'YMV09': 56, 'YMV10': 39, 'YMV11': 42, 'YMV12': 44, 'YMV13': 35, 'YMV14': 27, 'YMV15': 42, 'YMV16': 39, 'YMV17': 28, 'YMV18': 35, 'YMV19': 39}
Transitions averaged per mouse: median 1169 range 114 to 1541
CSD over +-1 site(s), window 0-100 ms: landmarks placed in 16 of 19 mice (not placed in 3: b at the top site or d at the bottom site, so a or e has no site beyond it)
distance medians (um): {'a-b': 100.0, 'b-c': 180.0, 'c-d': 130.0, 'd-e': 50.0}
range across mice (um): {'a-b': 480.0, 'b-c': 600.0, 'c-d': 340.0, 'd-e': 260.0} ; all < 110: False
largest deviation from the mean (um): {'a-b': 328.0, 'b-c': 374.0, 'c-d': 202.0, 'd-e': 182.0} ; all < 110: False
CSD over +-1 site(s), window 0-50 ms: landmarks placed in 15 of 19 mice (not placed in 4: b at the top site or d at the bottom site, so a or e has no site beyond it)
distance medians (um): {'a-b': 220.0, 'b-c': 180.0, 'c-d': 140.0, 'd-e': 60.0}
range across mice (um): {'a-b': 480.0, 'b-c': 580.0, 'c-d': 320.0, 'd-e': 260.0} ; all < 110: False
largest deviation from the mean (um): {'a-b': 292.0, 'b-c': 377.0, 'c-d': 197.0, 'd-e': 180.0} ; all < 110: False
CSD over +-2 site(s), window 0-100 ms: landmarks placed in 16 of 19 mice (not placed in 3: b at the top site or d at the bottom site, so a or e has no site beyond it)
distance medians (um): {'a-b': 230.0, 'b-c': 230.0, 'c-d': 140.0, 'd-e': 150.0}
range across mice (um): {'a-b': 620.0, 'b-c': 620.0, 'c-d': 320.0, 'd-e': 220.0} ; all < 110: False
largest deviation from the mean (um): {'a-b': 420.0, 'b-c': 381.0, 'c-d': 202.0, 'd-e': 114.0} ; all < 110: False
CSD over +-2 site(s), window 0-50 ms: landmarks placed in 16 of 19 mice (not placed in 3: b at the top site or d at the bottom site, so a or e has no site beyond it)
distance medians (um): {'a-b': 230.0, 'b-c': 230.0, 'c-d': 140.0, 'd-e': 120.0}
range across mice (um): {'a-b': 620.0, 'b-c': 620.0, 'c-d': 320.0, 'd-e': 220.0} ; all < 110: False
largest deviation from the mean (um): {'a-b': 406.0, 'b-c': 388.0, 'c-d': 192.0, 'd-e': 125.0} ; all < 110: False
subject transitions c a b d e a-b b-c c-d d-e
YMV01 962 32 14.0 15.0 39.0 46.0 20.0 340.0 140.0 140.0
YMV02 1027 31 6.0 26.0 37.0 38.0 400.0 100.0 120.0 20.0
YMV03 1222 31 6.0 26.0 41.0 53.0 400.0 100.0 200.0 240.0
YMV04 846 39 23.0 37.0 56.0 57.0 280.0 40.0 340.0 20.0
YMV05 953 38 21.0 22.0 56.0 57.0 20.0 320.0 360.0 20.0
YMV06 1504 38 25.0 32.0 41.0 53.0 140.0 120.0 60.0 240.0
YMV07 374 39 36.0 37.0 41.0 43.0 20.0 40.0 40.0 40.0
YMV08 982 35 15.0 26.0 37.0 51.0 220.0 180.0 40.0 280.0
YMV09 1399 56 1.0 26.0 57.0 61.0 500.0 600.0 20.0 80.0
YMV10 1541 39 20.0 21.0 54.0 57.0 20.0 360.0 300.0 60.0
YMV11 1481 42 20.0 21.0 56.0 57.0 20.0 420.0 280.0 20.0
YMV12 1240 44 20.0 21.0 55.0 56.0 20.0 460.0 220.0 20.0
YMV13 1464 35 1.0 4.0 39.0 43.0 60.0 620.0 80.0 80.0
YMV14 114 27 11.0 26.0 36.0 37.0 300.0 20.0 180.0 20.0
YMV15 1379 42 NaN NaN NaN NaN NaN NaN NaN NaN
YMV16 127 39 36.0 37.0 44.0 45.0 20.0 40.0 100.0 20.0
YMV17 1169 28 NaN NaN NaN NaN NaN NaN NaN NaN
YMV18 1472 35 10.0 26.0 37.0 50.0 320.0 180.0 40.0 260.0
YMV19 1146 39 NaN NaN NaN NaN NaN NaN NaN NaN
CSD colour scale clipped at the 99th percentile of |CSD| (118.8 uV per site^2); 388 of 38750 values beyond it
# CSD of the DOWN-UP transitions in every mouse with the landmarks placed by the rule (0-100 ms, one site).
fig, axs = plt.subplots(2, 10, figsize=(20, 7), sharey=True)
for ax, r in zip(axs.ravel(), S):
L = r["landmark"]
cs = v.csd_of(L["lfp"], 1)
lim = np.percentile(np.abs(cs), 99)
ax.imshow(cs.T, aspect="auto", cmap="RdBu", vmin=-lim, vmax=lim,
extent=[L["lags"][0] * 1000, L["lags"][-1] * 1000, 62.5, 0.5])
m = v.landmarks(cs, L["lags"], r["mua_power"], v.CSD_WINDOWS[0], 1)
ax.axhline(int(np.argmax(r["mua_power"])), color="g", lw=1)
if m:
for k in "abde":
ax.text(250, m[k], k, color="k", fontsize=9, va="center")
ax.plot([0, 100], [m[k], m[k]], color="k", lw=0.6)
ax.set_title(r["subject"], fontsize=9)
axs.ravel()[-1].axis("off")
axs[0, 0].set_ylabel("site"); axs[1, 0].set_ylabel("site")
plt.tight_layout(); plt.show()
print("Green: MUA peak (c). Colour scale of each map clipped at its 99th percentile of |CSD|.")
edge = LM[(LM.step == 1) & (LM.window == "0-100 ms")].dropna(subset=["a"])
print("Mice (rule, 0-100 ms, one site) with landmark a or b within 3 sites of the probe top:",
edge[(edge.a <= 3) | (edge.b <= 3)].subject.tolist(), "; with d or e within 3 sites of the bottom:",
edge[(edge.d >= 59) | (edge.e >= 59)].subject.tolist())
Green: MUA peak (c). Colour scale of each map clipped at its 99th percentile of |CSD|. Mice (rule, 0-100 ms, one site) with landmark a or b within 3 sites of the probe top: ['YMV09', 'YMV13'] ; with d or e within 3 sites of the bottom: ['YMV09']
Result: blocked (methods_underspecified). The data columns are in depth order (the site order check above), so the CSD is computed along depth. Under every reading the rule cannot place all five landmarks in 3 or 4 of the 19 mice (b at the top site or d at the bottom site), and among those placed the adjacent-landmark distances range across mice by 220 to 620 um, far beyond 110 um. The maps show why: the rule takes the extreme CSD over the whole half of the probe, so it places a or b within 3 sites of the probe top in YMV09 and YMV13 and d or e within 3 sites of the bottom in YMV09, and in several mice a falls one site above b. The paper chose the landmarks from the CSD maps by eye ("the most prominent ... sinks and sources", Figure 1D), so the rule run here is not the paper's procedure, and no rule given in the paper or the claims file places them; the result of every reading is shown above.
figs7-strongest-3-6hz-coupling¶
The strongest spike-LFP coupling is in the 3-6 Hz band, followed by the slow gamma (30-60 Hz) band.
Criterion: per unit in waking, the frequency of largest mean resultant length (preferred frequency, STAR Methods) over log-spaced frequencies from 0.5 to 100 Hz; "strongest" read as (a) the histogram of preferred frequencies of significantly modulated units or (b) the mean modulation strength over units at each frequency; reproduced when, under a reading, the largest value falls in 3-6 Hz and the next in 30-60 Hz, among bands 0.5-3, 3-6, 6-15, 15-30, 30-60 and 60-100 Hz.
Method details taken from the STAR Methods (Spectral analysis, cross-frequency coupling, and spike-LFP coupling): complex Morlet wavelets (Torrence and Compo 1998; w0 = 6 is their standard), phase at each spike, mean resultant length, Rayleigh p < 0.05, preferred frequency as the largest mean vector length. Neither the Methods nor the Figure S7 legend say which LFP site the phases come from. Two readings are run: the unit's own peak site, and the site 10 sites (200 um) away from it, which keeps the unit's own spike out of the LFP (the 1250 Hz LFP of the unit's site still holds part of its spike). The second is adopted as an assumption. The Methods compute the wavelet "of the LFP (or ICs)", so a third reading, with the gamma bands on the gamma IC of largest coupling, is run after the ICA below. Waking is the Awake intervals of at least 30 s, the longest first, up to 600 s per session in total, with spikes within 8 s of an epoch edge left out (beyond the +-4 SD of the 0.5 Hz wavelet); units with fewer than 50 spikes there are left out. 25 log-spaced frequencies; the LFP is decimated to 312.5 Hz. YMV19 has no state table and is left out. A unit is significantly modulated when the Rayleigh p at its preferred frequency is below 0.05.
bands = [(0.5, 3), (3, 6), (6, 15), (15, 30), (30, 60), (60, 100.01)]
bnames = ["0.5-3", "3-6", "6-15", "15-30", "30-60", "60-100"]
F = v.COUPLING_FREQS
print("Frequencies (Hz):", np.round(F, 2))
cp = [(si, r["coupling"]) for si, r in enumerate(S) if r["coupling"] is not None]
print("Sessions with waking coupling:", len(cp), "; waking seconds used: median", np.median([c["seconds"] for _, c in cp]),
"range", min(c["seconds"] for _, c in cp), "to", max(c["seconds"] for _, c in cp))
for reading in ("own", "distant"):
pref, mrl_all, subj, n_drop, n_units = [], [], [], 0, 0
for si, c in cp:
ok = c["n_spikes"] >= 50
n_drop += int((~ok).sum()); n_units += len(ok)
m = np.nan_to_num(c[f"mrl_{reading}"], nan=-1)[ok]
p = c[f"p_{reading}"][ok]
i = np.argmax(m, axis=1)
sig = p[np.arange(len(i)), i] < 0.05
pref += list(F[i][sig]); subj += [S[si]["subject"]] * int(sig.sum())
mrl_all.append(m)
pref = np.array(pref); mrl_all = np.vstack(mrl_all)
hist = [int(((pref >= a) & (pref < b)).sum()) for a, b in bands]
mean = [mrl_all[:, (F >= a) & (F < b)].mean() for a, b in bands]
print(f"\nLFP of the {reading} site: {len(mrl_all)} units ({n_drop} of {n_units} left out, < 50 waking spikes), "
f"{len(pref)} significantly modulated")
print(" (a) preferred frequencies by band:", dict(zip(bnames, hist)),
"; order:", [bnames[j] for j in np.argsort(hist)[::-1]][:3])
print(" (b) mean MRL by band:", dict(zip(bnames, np.round(mean, 4))),
"; order:", [bnames[j] for j in np.argsort(mean)[::-1]][:3])
if reading == "distant":
P = pd.DataFrame({"subject": subj, "band": pd.cut(pref, [a for a, _ in bands] + [100.01], labels=bnames, right=False)})
print(" preferred band by mouse (counts):")
print(pd.crosstab(P.subject, P.band).to_string())
print(f" fraction significant at the preferred frequency: {len(pref) / len(mrl_all):.2f}; with 25 frequencies "
f"tested, chance would give at most {1 - 0.95 ** 25:.2f} (independent frequencies)")
fig, ax = plt.subplots(figsize=(5, 3))
for reading, col in (("own", "0.6"), ("distant", "k")):
mm = np.vstack([np.nan_to_num(c[f"mrl_{reading}"], nan=np.nan)[c["n_spikes"] >= 50] for _, c in cp])
ax.semilogx(F, np.nanmean(mm, axis=0), color=col, label=reading)
ax.set_xlabel("Hz"); ax.set_ylabel("mean MRL"); ax.legend(); plt.tight_layout(); plt.show()
Frequencies (Hz): [ 0.5 0.62 0.78 0.97 1.21 1.51 1.88 2.34 2.92 3.65
4.55 5.67 7.07 8.82 11. 13.71 17.1 21.32 26.59 33.16
41.35 51.57 64.31 80.19 100. ]
Sessions with waking coupling: 18 ; waking seconds used: median 600.0 range 357.0 to 600.0
LFP of the own site: 1293 units (120 of 1413 left out, < 50 waking spikes), 1290 significantly modulated
(a) preferred frequencies by band: {'0.5-3': 253, '3-6': 203, '6-15': 63, '15-30': 18, '30-60': 39, '60-100': 714} ; order: ['60-100', '0.5-3', '3-6']
(b) mean MRL by band: {'0.5-3': np.float64(0.0655), '3-6': np.float64(0.08), '6-15': np.float64(0.0607), '15-30': np.float64(0.055), '30-60': np.float64(0.0699), '60-100': np.float64(0.1078)} ; order: ['60-100', '3-6', '30-60']
LFP of the distant site: 1293 units (120 of 1413 left out, < 50 waking spikes), 1270 significantly modulated
(a) preferred frequencies by band: {'0.5-3': 591, '3-6': 344, '6-15': 107, '15-30': 40, '30-60': 65, '60-100': 123} ; order: ['0.5-3', '3-6', '60-100']
(b) mean MRL by band: {'0.5-3': np.float64(0.0647), '3-6': np.float64(0.0729), '6-15': np.float64(0.0528), '15-30': np.float64(0.044), '30-60': np.float64(0.049), '60-100': np.float64(0.0509)} ; order: ['3-6', '0.5-3', '6-15']
preferred band by mouse (counts):
band 0.5-3 3-6 6-15 15-30 30-60 60-100
subject
YMV01 27 8 2 1 1 10
YMV02 51 14 14 1 4 10
YMV03 49 29 5 3 3 19
YMV04 18 5 4 5 9 15
YMV05 31 29 4 0 3 4
YMV06 28 36 2 3 1 2
YMV07 65 2 3 2 8 6
YMV08 31 46 9 3 4 3
YMV09 23 14 15 1 2 4
YMV10 7 5 9 2 2 9
YMV11 22 28 6 1 1 5
YMV12 20 9 1 3 3 8
YMV13 27 41 6 0 0 0
YMV14 31 15 8 6 4 6
YMV15 41 17 5 1 4 7
YMV16 15 29 0 2 2 1
YMV17 39 3 3 1 9 11
YMV18 66 14 11 5 5 3
fraction significant at the preferred frequency: 0.98; with 25 frequencies tested, chance would give at most 0.72 (independent frequencies)
Result: partially reproduced, resting on an assumption about the LFP site. With the LFP 10 sites away from each unit (the adopted reading), the mean modulation strength is largest in 3-6 Hz (0.073), as the claim says, but the next band is 0.5-3 Hz (0.065), not 30-60 Hz (0.049, fifth of six). By the histogram of preferred frequencies, 0.5-3 Hz comes first (591 units), 3-6 Hz second (344) and 30-60 Hz fifth (65). With the LFP of the unit's own site, 60-100 Hz comes first under both readings (714 units; mean 0.108), which fits the unit's own spike leaking into the LFP. With the gamma bands on the best gamma IC (reading c, after the ICA below), 60-100 Hz (0.118) and 30-60 Hz (0.098) come before 3-6 Hz (0.073). Under no reading does 3-6 Hz come first with 30-60 Hz second.
fig4a-unit-counts¶
Clustering one session from each of the 19 mice yields 1,472 units, with 251 narrow-waveform putative interneurons and 1,075 wide-waveform putative principal cells.
Criterion: units as sorted in the dataset; classes by the waveform rules of the STAR Methods, trough-to-peak latency from the mean waveform on the peak channel. The dataset holds no waveforms, so the mean waveform of each unit is averaged from 60 s of the raw signal (two 30 s blocks, at 1000 s and 3000 s), on all 64 sites; the peak site is the one with the largest peak-to-peak amplitude. The Results say "the unfiltered waveform was quantified by the trough-to-peak latency" (Assignment of Neurons to Physiologically Defined Layers, paragraph 1), so the main reading is the wide-band signal; the 300 Hz high-pass that Kilosort used (STAR Methods, Spike sorting, ops.fshigh 300) is a variant. Units with fewer than 10 spikes in the blocks are left unclassified.
print("Units in the files:", len(U), "; with quality 'good':", int((U.quality == "good").sum()))
for col, name in (("wf", "wide band (main)"), ("wf_hp", "300 Hz high-pass (variant)")):
print(f"Waveform classes, {name}:", U[col].value_counts().to_dict())
print("Spikes used per unit for the waveform: median", int(U.wf_spikes.median()), "; units with fewer than 10:",
int((U.wf_spikes < 10).sum()))
neg = U.wf.isin(["E", "I"])
print("\nThreshold sensitivity (wide band, negative non-fiber units):")
for th in (0.45, 0.50, 0.55, 0.60, 0.65):
print(f" TP <= {th:.2f} ms: I {int((U.tp_ms[neg] <= th).sum())}, E {int((U.tp_ms[neg] > th).sum())}")
print("Trough-to-peak latencies are quantized at", 1000 / v.RAW_RATE, "ms (20 kHz); range",
U.tp_ms.min(), "to", U.tp_ms.max(), "ms")
print("\nBy mouse (wide band):")
print(pd.crosstab(U.subject, U.wf, margins=True).T.to_string())
print("\nAgainst the dataset's own cell_type labels (wide band):")
print(pd.crosstab(U.wf, U.cell_type).to_string())
fig, ax = plt.subplots(1, 2, figsize=(9, 3))
ax[0].hist(U.tp_ms[neg], bins=np.arange(0, 1.6, 0.05), color="0.4"); ax[0].axvline(0.55, color="r")
ax[0].set_xlabel("trough-to-peak (ms), wide band"); ax[0].set_ylabel("units")
ax[1].hist(U.tp_hp[U.wf_hp.isin(["E", "I"])], bins=np.arange(0, 1.6, 0.05), color="0.4"); ax[1].axvline(0.55, color="r")
ax[1].set_xlabel("trough-to-peak (ms), 300 Hz high-pass")
plt.tight_layout(); plt.show()
print("Units outside the 0 to 1.6 ms axes: wide band", int((U.tp_ms[neg] >= 1.6).sum()),
"; high-pass", int((U.tp_hp[U.wf_hp.isin(['E', 'I'])] >= 1.6).sum()))
Units in the files: 1477 ; with quality 'good': 1476
Waveform classes, wide band (main): {np.str_('E'): 1016, np.str_('I'): 246, np.str_('unclassified'): 110, np.str_('positive'): 105}
Waveform classes, 300 Hz high-pass (variant): {np.str_('E'): 941, np.str_('I'): 320, np.str_('unclassified'): 110, np.str_('positive'): 106}
Spikes used per unit for the waveform: median 127 ; units with fewer than 10: 110
Threshold sensitivity (wide band, negative non-fiber units):
TP <= 0.45 ms: I 204, E 1058
TP <= 0.50 ms: I 215, E 1047
TP <= 0.55 ms: I 246, E 1016
TP <= 0.60 ms: I 323, E 939
TP <= 0.65 ms: I 472, E 790
Trough-to-peak latencies are quantized at 0.05 ms (20 kHz); range 0.15 to 2.0 ms
By mouse (wide band):
subject YMV01 YMV02 YMV03 YMV04 YMV05 YMV06 YMV07 YMV08 YMV09 YMV10 YMV11 YMV12 YMV13 YMV14 YMV15 YMV16 YMV17 YMV18 YMV19 All wf E 42 73 89 46 64 52 74 75 33 24 50 40 56 49 53 32 53 70 41 1016 I 9 13 16 10 8 19 12 16 12 9 14 5 13 16 16 13 12 20 13 246 positive 0 9 7 6 0 1 4 6 19 1 2 1 6 5 7 4 5 17 5 105 unclassified 2 6 4 1 0 1 11 3 14 4 4 14 4 9 3 2 11 12 5 110 All 53 101 116 63 72 73 101 100 78 38 70 60 79 79 79 51 81 119 64 1477 Against the dataset's own cell_type labels (wide band): cell_type Narrow Interneuron Pyramidal Cell Unknown Wide Interneuron down state missing wf E 18 696 6 198 9 89 I 183 32 7 8 0 16 positive 50 39 2 7 0 7 unclassified 26 69 0 10 1 4
Units outside the 0 to 1.6 ms axes: wide band 25 ; high-pass 0
Result: reproduced. The files hold 1,477 units (1,476 of quality "good"), within 1,472 +- 147. From the wide-band waveform, as the Results describe, 246 units are narrow-waveform I cells (251 +- 25) and 1,016 wide-waveform E cells (1,075 +- 108); 105 are positive units, no unit is a fiber by the kurtosis rule, and 110 have fewer than 10 spikes in the 60 s used and are left unclassified. With a 300 Hz high-pass the counts are 320 and 941, both outside tolerance, and the I count also leaves tolerance at thresholds of 0.60 ms (323) or 0.50 ms (215).
fig4a-ccg-classification¶
Short-latency CCG interactions identify 347 E and 100 I cells, and these agree with the waveform classes (97.4% of CCG-based E cells are waveform E cells; 98.0% of CCG-based I cells are waveform I cells).
Criterion: a unit is a CCG-based E (I) cell when it is the presynaptic unit of at least one significant peak (trough) by the monosynaptic rule of the STAR Methods: CCGs in 0.5 ms bins over the whole session, a Gaussian of 7 ms SD as the predictor, the 99.9999th Poisson percentile (0.0001th for troughs), at least two consecutive bins within +1.5 to +4 ms. The unhollowed Gaussian (Methods) is the main reading and the buzcode hollowed one (0.6) a variant (the claims file's opening comment). The percentages are taken over all CCG-based cells (main) and over those with an E or I waveform class (variant). The dataset's own synaptic_effect labels are shown for comparison.
for e, i, name in (("ccgE", "ccgI", "unhollowed (main)"), ("ccgE_h", "ccgI_h", "hollowed 0.6 (variant)")):
nE, nI = int(U[e].sum()), int(U[i].sum())
pE = (U.wf[U[e]] == "E").mean() * 100
pI = (U.wf[U[i]] == "I").mean() * 100
pE2 = (U.wf[U[e] & U.wf.isin(["E", "I"])] == "E").mean() * 100
pI2 = (U.wf[U[i] & U.wf.isin(["E", "I"])] == "I").mean() * 100
print(f"{name}: CCG-based E {nE}, I {nI}, both {int((U[e] & U[i]).sum())}; "
f"E that are waveform E {pE:.1f}% (of classified {pE2:.1f}%), I that are waveform I {pI:.1f}% (of classified {pI2:.1f}%)")
print("\nCCG-based E cells by waveform class:", U.wf[U.ccgE].value_counts().to_dict())
print("CCG-based I cells by waveform class:", U.wf[U.ccgI].value_counts().to_dict())
print("\nBy mouse (unhollowed):")
print(U.groupby("subject").agg(units=("unit", "size"), ccgE=("ccgE", "sum"), ccgI=("ccgI", "sum"),
ds_exc=("synaptic_effect", lambda s: (s == "Excitatory").sum()),
ds_inh=("synaptic_effect", lambda s: (s == "Inhibitory").sum())).T.to_string())
print("\nThe dataset's synaptic_effect labels against this detection:")
print(pd.crosstab(U.synaptic_effect, [U.ccgE.rename("ccgE"), U.ccgI.rename("ccgI")]).to_string())
lab = U.synaptic_effect.isin(["Excitatory", "Inhibitory"])
print("Dataset labels: Excitatory", int((U.synaptic_effect == "Excitatory").sum()), "Inhibitory",
int((U.synaptic_effect == "Inhibitory").sum()), "; Excitatory that are waveform E",
round((U.wf[U.synaptic_effect == "Excitatory"] == "E").mean() * 100, 1), "%; Inhibitory that are waveform I",
round((U.wf[U.synaptic_effect == "Inhibitory"] == "I").mean() * 100, 1), "%")
unhollowed (main): CCG-based E 673, I 187, both 122; E that are waveform E 74.9% (of classified 80.9%), I that are waveform I 72.7% (of classified 79.5%)
hollowed 0.6 (variant): CCG-based E 682, I 188, both 122; E that are waveform E 74.9% (of classified 81.0%), I that are waveform I 72.3% (of classified 79.1%)
CCG-based E cells by waveform class: {np.str_('E'): 504, np.str_('I'): 119, np.str_('positive'): 27, np.str_('unclassified'): 23}
CCG-based I cells by waveform class: {np.str_('I'): 136, np.str_('E'): 35, np.str_('unclassified'): 9, np.str_('positive'): 7}
By mouse (unhollowed):
subject YMV01 YMV02 YMV03 YMV04 YMV05 YMV06 YMV07 YMV08 YMV09 YMV10 YMV11 YMV12 YMV13 YMV14 YMV15 YMV16 YMV17 YMV18 YMV19
units 53 101 116 63 72 73 101 100 78 38 70 60 79 79 79 51 81 119 64
ccgE 14 34 37 14 36 22 75 42 28 22 41 29 32 41 41 34 47 56 28
ccgI 1 7 9 4 10 9 16 8 15 4 7 5 9 14 14 10 15 24 6
ds_exc 14 33 0 9 36 23 51 30 9 15 23 24 15 29 28 20 26 42 17
ds_inh 2 5 0 4 8 14 16 9 10 5 5 4 9 10 10 8 9 11 6
The dataset's synaptic_effect labels against this detection:
ccgE False True ccgI False True False True synaptic_effect Excitatory 93 2 329 20 Inhibitory 18 36 11 80 Unknown 552 24 180 16 missing 76 3 31 6 Dataset labels: Excitatory 444 Inhibitory 145 ; Excitatory that are waveform E 86.3 %; Inhibitory that are waveform I 79.3 %
Result: not reproduced. The STAR Methods rule gives 673 CCG-based E cells and 187 I cells (347 +- 35 and 100 +- 10), 122 of them both. Of the CCG-based E cells 74.9% are waveform E cells (80.9% of those with a class), and of the I cells 72.7% are waveform I cells (79.5%), against 97.4% and 98.0%. The hollowed predictor changes little (682 and 188). The dataset's own synaptic_effect labels count 444 excitatory and 145 inhibitory units, also more than the paper's.
fig7-connection-counts¶
Of 102,332 potential neuron pairs, 525 putative excitatory and 323 putative inhibitory monosynaptic connections were identified.
Criterion: detection by the STAR Methods rule. The pair count is read as ordered pairs within a session among (a) all units or (b) E and I units only (waveform classes), and as unordered pairs of each. The Gaussian baseline is run unhollowed (Methods) and hollowed 0.6 (buzcode). Reproduced when, under a reading, all three counts are within tolerance (102,332 ± 10,233; 525 ± 53; 323 ± 32).
The site distance between the two units of each connection is shown as a diagnostic: units 0 or 1 site apart (within 20 um) can share spikes after sorting, which would add peaks the paper may not have counted. It is not a criterion.
n_all = U.groupby("sess").size().values
n_ei = U[U.wf.isin(["E", "I"])].groupby("sess").size().reindex(range(len(S)), fill_value=0).values
pairs = {"ordered, all units": int((n_all * (n_all - 1)).sum()), "ordered, E and I units": int((n_ei * (n_ei - 1)).sum()),
"unordered, all units": int((n_all * (n_all - 1) // 2).sum()), "unordered, E and I units": int((n_ei * (n_ei - 1) // 2).sum())}
for k, val in pairs.items():
print(f"potential pairs, {k}: {val} (within 102,332 +/- 10,233: {abs(val - 102332) <= 10233})")
for key, name in (("conn", "unhollowed (main)"), ("conn_hollow", "hollowed 0.6")):
ne = int(sum(r[key]["exc"].sum() for r in S)); ni = int(sum(r[key]["inh"].sum() for r in S))
print(f"{name}: excitatory {ne} (525 +/- 53: {abs(ne - 525) <= 53}), inhibitory {ni} (323 +/- 32: {abs(ni - 323) <= 32})")
# among E and I units only
ei = C[C.pre_wf.isin(["E", "I"]) & C.post_wf.isin(["E", "I"])]
print("Connections between units with an E or I waveform class: excitatory", int(ei.exc.sum()), "inhibitory", int(ei.inh.sum()))
print("\nBy mouse:")
print(C.groupby("subject").agg(exc=("exc", "sum"), inh=("inh", "sum")).T.to_string())
print("\nSite distance between pre and post units (sites of 20 um):")
print(pd.DataFrame({"excitatory": C[C.exc].dsite.clip(upper=10).value_counts().sort_index(),
"inhibitory": C[C.inh].dsite.clip(upper=10).value_counts().sort_index()}).T.to_string())
far = C[C.dsite >= 2]
print("Connections with units at least 2 sites apart: excitatory", int(far.exc.sum()), "inhibitory", int(far.inh.sum()))
print("Connections detected within each state's CCG (reference spikes in the state), summed over mice (excitatory, inhibitory):",
{st: tuple(int(sum(r["conn_state"][st][j] for r in S)) for j in (0, 1)) for st in ("Awake", "Non-REM")})
potential pairs, ordered, all units: 121842 (within 102,332 +/- 10,233: False)
potential pairs, ordered, E and I units: 89046 (within 102,332 +/- 10,233: False)
potential pairs, unordered, all units: 60921 (within 102,332 +/- 10,233: False)
potential pairs, unordered, E and I units: 44523 (within 102,332 +/- 10,233: False)
unhollowed (main): excitatory 1295 (525 +/- 53: False), inhibitory 578 (323 +/- 32: False)
hollowed 0.6: excitatory 1328 (525 +/- 53: False), inhibitory 591 (323 +/- 32: False)
Connections between units with an E or I waveform class: excitatory 1201 inhibitory 544
By mouse:
subject YMV01 YMV02 YMV03 YMV04 YMV05 YMV06 YMV07 YMV08 YMV09 YMV10 YMV11 YMV12 YMV13 YMV14 YMV15 YMV16 YMV17 YMV18 YMV19
exc 18 47 55 19 57 28 235 65 39 37 73 45 71 74 113 75 95 111 38
inh 1 16 37 11 30 24 96 29 37 14 14 17 25 46 44 42 29 51 15
Site distance between pre and post units (sites of 20 um):
dsite 0 1 2 3 4 5 6 7 8 9 10
excitatory 224 258 152 93 70 62 46 37 27 25 301
inhibitory 49 67 67 76 63 61 42 25 15 21 92
Connections with units at least 2 sites apart: excitatory 813 inhibitory 462
Connections detected within each state's CCG (reference spikes in the state), summed over mice (excitatory, inhibitory): {'Awake': (291, 138), 'Non-REM': (321, 236)}
Result: not reproduced. The STAR Methods rule finds 1,295 excitatory and 578 inhibitory connections (525 +- 53 and 323 +- 32), and 1,328 and 591 with the hollowed predictor. No reading of the pair count is within 102,332 +- 10,233: ordered pairs of all units give 121,842, of E and I units 89,046, and unordered pairs 60,921 and 44,523. Of the excitatory connections 482 join units 0 or 1 site apart; with units at least 2 sites apart there remain 813 excitatory and 462 inhibitory, still above the paper's counts. Detected within each state's CCG, there are 291 and 138 in waking and 321 and 236 in non-REM (a diagnostic; the claim counts whole-session connections).
fig7c-divergence-convergence¶
Most presynaptic principal neurons have a single target in the recorded population, a few target several; an interneuron typically has one presynaptic partner, a few many, up to 16.
Criterion: from the E-I connections (excitatory connections from a waveform E unit to a waveform I unit), the count of targets per presynaptic E cell and of presynaptic E cells per I cell. Reproduced when the mode of both distributions is one and the largest convergence is within 16 ± 1.6.
EI = C[C.exc & (C.pre_wf == "E") & (C.post_wf == "I")]
div = EI.groupby(["sess", "pre"]).size()
con = EI.groupby(["sess", "post"]).size()
print("E-I connections:", len(EI), "from", len(div), "E cells to", len(con), "I cells")
print("targets per presynaptic E cell:", div.value_counts().sort_index().to_dict(), "; mode", int(div.mode()[0]))
print("presynaptic E cells per I cell:", con.value_counts().sort_index().to_dict(), "; mode", int(con.mode()[0]),
"; largest", int(con.max()))
print("I cells with more than 16 presynaptic E cells:", int((con > 16).sum()))
print("Largest convergence by mouse:", con.groupby(level=0).max().rename(lambda i: S[i]["subject"]).to_dict())
EIh = []
for si, r in enumerate(S):
u = U[U.sess == si].wf.values
ex = r["conn_hollow"]["exc"] & (u == "E")[:, None] & (u == "I")[None, :]
EIh.append(ex.sum(axis=0)[ex.sum(axis=0) > 0])
conh = np.concatenate(EIh)
print("Hollowed variant: presynaptic E cells per I cell, mode", int(stats.mode(conh).mode), "; largest", int(conh.max()))
E-I connections: 549 from 362 E cells to 138 I cells
targets per presynaptic E cell: {1: 242, 2: 79, 3: 28, 4: 4, 5: 7, 7: 2} ; mode 1
presynaptic E cells per I cell: {1: 39, 2: 23, 3: 15, 4: 20, 5: 9, 6: 7, 7: 8, 8: 5, 9: 2, 10: 3, 11: 3, 14: 2, 18: 1, 29: 1} ; mode 1 ; largest 29
I cells with more than 16 presynaptic E cells: 2
Largest convergence by mouse: {'YMV01': 3, 'YMV02': 4, 'YMV03': 11, 'YMV04': 5, 'YMV05': 8, 'YMV06': 3, 'YMV07': 29, 'YMV08': 8, 'YMV09': 2, 'YMV10': 6, 'YMV11': 7, 'YMV12': 8, 'YMV13': 5, 'YMV14': 7, 'YMV15': 10, 'YMV16': 8, 'YMV17': 9, 'YMV18': 11, 'YMV19': 4}
Hollowed variant: presynaptic E cells per I cell, mode 1 ; largest 29
Result: partially reproduced. Of 549 E-I connections, most presynaptic E cells have one target (242 of 362, mode 1), and the most common number of presynaptic E cells per I cell is one (39 of 138 I cells, mode 1). The largest convergence is 29 (one I cell in YMV07; 18 for the next), outside 16 +- 1.6, so the claim holds in its shape but not in its maximum.
figs8f-ie-dip-wake¶
Over all I-E pairs, spike suppression is stronger in waking than in non-REM sleep (dip 24.1% against 22.8% of the CCG baseline).
Criterion: for each I-E connection and state, the dip as (predictor minus CCG) over predictor at the deepest bin in +1.5 to +4 ms; means and medians over pairs are both readings; reproduced when both values are within tolerance (24.1 ± 2.41 and 22.8 ± 2.28%) and waking > non-REM by a two-sided signed-rank test at p < 0.05. The CCG of a state uses the presynaptic spikes in that state (Awake or Non-REM intervals) and its own predictor. Which connections are "I-E" is taken from the Results (Within-Layer and Inter-layer Interactions among Neurons, paragraph 1), where the I-E connections number 323, the count of all inhibitory connections: the main reading is every inhibitory connection, and the variant those from a waveform I unit to a waveform E unit. YMV19 has no state table and is left out.
for name, sel in (("all inhibitory connections (main)", C.inh),
("waveform I to waveform E (variant)", C.inh & (C.pre_wf == "I") & (C.post_wf == "E"))):
x = C[sel]
ok = x.dip_wake.notna() & x.dip_nrem.notna()
x = x[ok]
w = stats.wilcoxon(x.dip_wake, x.dip_nrem)
print(f"{name}: {len(x)} pairs ({int((~ok).sum())} without spikes in a state); dip wake mean {x.dip_wake.mean()*100:.1f}%, "
f"non-REM mean {x.dip_nrem.mean()*100:.1f}%; medians {x.dip_wake.median()*100:.1f}% and {x.dip_nrem.median()*100:.1f}%; "
f"wake > non-REM in {int((x.dip_wake > x.dip_nrem).sum())} pairs; signed-rank p = {w.pvalue:.3g}; "
f"median difference {np.median(x.dip_wake - x.dip_nrem)*100:.2f} points")
x = C[C.inh & C.dip_wake.notna() & C.dip_nrem.notna()]
print("\nBy mouse (main reading):")
print(x.groupby("subject").agg(pairs=("pre", "size"), wake=("dip_wake", "median"), nrem=("dip_nrem", "median")).round(3).T.to_string())
print("Dip range (wake, non-REM):", round(x.dip_wake.min(), 3), round(x.dip_wake.max(), 3), round(x.dip_nrem.min(), 3), round(x.dip_nrem.max(), 3))
all inhibitory connections (main): 563 pairs (15 without spikes in a state); dip wake mean 39.6%, non-REM mean 39.9%; medians 35.6% and 36.0%; wake > non-REM in 265 pairs; signed-rank p = 0.366; median difference -0.14 points waveform I to waveform E (variant): 361 pairs (11 without spikes in a state); dip wake mean 40.0%, non-REM mean 39.4%; medians 37.1% and 36.3%; wake > non-REM in 182 pairs; signed-rank p = 0.568; median difference 0.07 points By mouse (main reading): subject YMV01 YMV02 YMV03 YMV04 YMV05 YMV06 YMV07 YMV08 YMV09 YMV10 YMV11 YMV12 YMV13 YMV14 YMV15 YMV16 YMV17 YMV18 pairs 1.000 16.000 37.000 11.000 30.000 24.000 96.000 29.000 37.000 14.000 14.000 17.000 25.000 46.000 44.000 42.000 29.000 51.000 wake 0.254 0.358 0.396 0.363 0.357 0.467 0.347 0.418 0.402 0.348 0.387 0.376 0.322 0.358 0.265 0.250 0.441 0.342 nrem 0.533 0.429 0.415 0.382 0.309 0.372 0.347 0.422 0.413 0.404 0.377 0.378 0.327 0.355 0.241 0.311 0.361 0.360 Dip range (wake, non-REM): 0.032 1.0 0.057 1.0
Result: not reproduced. Over 563 inhibitory connections with spikes in both states (the main reading counts every inhibitory connection, including those onto I cells, as the paper's I-E count of 323 equals its count of all inhibitory connections), the dip is 39.6% of the predictor in waking and 39.9% in non-REM (medians 35.6% and 36.0%), outside 24.1 +- 2.41% and 22.8 +- 2.28% under either reading, and waking does not exceed non-REM (signed-rank p = 0.37; 265 of 563 pairs larger in waking). The strict I-to-E reading gives 40.0% and 39.4% (p = 0.57).
fig4de-optotagging¶
Fourteen of 15 optogenetically identified PV cells and 40 of 52 VGAT cells fall in the narrow-waveform inhibitory group.
Criterion: tagging by the STAR Methods rule (1 ms PSTH bins, mean of the two largest bins in 1 to 6 ms after the 10 ms pulse above the baseline mean plus 8 SD of -100 to 0 ms) in the PV-Cre::Ai32 and VGAT-Cre::Ai32 mice; readings, all pulses pooled, and each intensity separately with a unit tagged when any passes. The pulses are the dataset's Laser diode intervals (10 ms; three amplitude levels, read as the three intensities). The mice are taken by the genotype in the files, which gives 11 PV-Cre::Ai32 and 3 VGAT-Cre::Ai32 mice where the paper has 5 and 9. Units tagged in the other genotypes and with onsets shifted by 0.5 s are shown as references.
for col, name in (("opto_pooled", "pooled (main)"), ("opto_any", "any intensity")):
for g in ("PV-Cre::Ai32", "VGAT-Cre::Ai32"):
x = U[(U.genotype == g) & U[col]]
print(f"{name}, {g}: tagged {len(x)}, in the waveform I group {int((x.wf == 'I').sum())}, "
f"waveform E {int((x.wf == 'E').sum())}, other {int((~x.wf.isin(['E', 'I'])).sum())}")
oth = U[~U.genotype.isin(["PV-Cre::Ai32", "VGAT-Cre::Ai32"]) & U[col]]
print(f"{name}, other genotypes: tagged {len(oth)}; shifted onsets, all mice: {int(U[col + '_null'].sum())}")
print("\nBy mouse (pooled):")
print(U[U.sess.map(lambda i: S[i]["opto"] is not None)].groupby(["subject", "genotype"])
.agg(units=("unit", "size"), tagged=("opto_pooled", "sum"),
tagged_I=("opto_pooled", lambda s: int((s & (U.loc[s.index, "wf"] == "I")).sum()))).to_string())
both = U.genotype.isin(["PV-Cre::Ai32", "VGAT-Cre::Ai32"])
print("\nPV and VGAT mice together: tagged", int(U[both].opto_pooled.sum()), "in I group",
int((U[both].opto_pooled & (U[both].wf == "I")).sum()), "(paper 67 and 54)")
pooled (main), PV-Cre::Ai32: tagged 62, in the waveform I group 48, waveform E 13, other 1
pooled (main), VGAT-Cre::Ai32: tagged 18, in the waveform I group 14, waveform E 4, other 0
pooled (main), other genotypes: tagged 0; shifted onsets, all mice: 0
any intensity, PV-Cre::Ai32: tagged 61, in the waveform I group 48, waveform E 12, other 1
any intensity, VGAT-Cre::Ai32: tagged 18, in the waveform I group 13, waveform E 4, other 1
any intensity, other genotypes: tagged 0; shifted onsets, all mice: 0
By mouse (pooled):
units tagged tagged_I
subject genotype
YMV05 PV-Cre::Ai32 72 7 7
YMV06 VGAT-Cre::Ai32 73 11 10
YMV07 CaMKII-Cre::Ai35 101 0 0
YMV08 VGAT-Cre::Ai32 100 4 3
YMV09 PV-Cre::Ai32 78 7 6
YMV10 PV-Cre::Ai32 38 0 0
YMV11 PV-Cre::Ai32 70 5 4
YMV12 VGAT-Cre::Ai32 60 3 1
YMV13 PV-Cre::Ai32 79 6 5
YMV14 PV-Cre::Ai32 79 3 1
YMV15 PV-Cre::Ai32 79 15 10
YMV16 PV-Cre::Ai32 51 2 2
YMV17 PV-Cre::Ai32 81 13 10
YMV18 PV-Cre::Ai32 119 1 1
YMV19 PV-Cre::Ai32 64 3 2
PV and VGAT mice together: tagged 80 in I group 62 (paper 67 and 54)
Result: not reproduced. In the 11 mice the files label PV-Cre::Ai32, 62 units pass the tagging rule with all pulses pooled, 48 of them in the waveform I group (paper 15 and 14). In the 3 VGAT-Cre::Ai32 mice 18 pass, 14 in the I group (paper 52 and 40). Testing each intensity gives 61 and 18. No unit passes in the other genotypes or with onsets shifted by 0.5 s. Taken together, 80 units are tagged and 62 of them (78%) are in the I group, against the paper's 67 and 54 (81%), also outside tolerance. The genotypes in the files (11 PV and 3 VGAT mice) contradict the paper's (5 and 9), so the per-genotype comparison does not line up with the paper's mice; which side is wrong cannot be told from the data, and the paper gives no per-mouse list.
figs8c-lognormal-strength (blocked, methods_underspecified)¶
Connection strengths are skewed and lognormally distributed, a few strong connections mixed with many weak ones.
The criterion records this claim as methods_underspecified: the paper reports no test or measure of lognormality. The distribution of log transmission probabilities is shown with its skewness, as description.
for name, sel in (("excitatory", C.exc), ("inhibitory", C.inh)):
s_ = C[sel].stp.abs()
s_ = s_[s_ > 0]
print(f"{name}: {len(s_)} connections; |transmission probability| median {s_.median():.4f}, range {s_.min():.5f} to {s_.max():.3f}; "
f"skewness of raw values {stats.skew(s_):.2f}, of log10 values {stats.skew(np.log10(s_)):.2f}")
fig, ax = plt.subplots(figsize=(4, 3))
ax.hist(np.log10(C[C.exc].stp[C[C.exc].stp > 0]), bins=30, color="0.4")
ax.set_xlabel("log10 transmission probability, excitatory"); ax.set_ylabel("connections"); plt.tight_layout(); plt.show()
print("Excitatory connections with transmission probability <= 0 (not plotted):", int((C[C.exc].stp <= 0).sum()))
excitatory: 1295 connections; |transmission probability| median 0.0085, range 0.00007 to 0.123; skewness of raw values 2.89, of log10 values -0.25 inhibitory: 578 connections; |transmission probability| median 0.0030, range 0.00006 to 0.210; skewness of raw values 8.45, of log10 values 0.05
Excitatory connections with transmission probability <= 0 (not plotted): 11
fig2c-six-gamma-ics¶
ICA of the 30-100 Hz LFP yields six main independent components whose spatial loadings are highly consistent across mice and whose depth distribution is related to the anatomical layers.
Criterion: per mouse, the LFP band-passed 30-100 Hz, reduced by PCA to its first ten PCs, then runica; ICs ranked by variance explained, loadings as the mixing matrix projected on depth. Readings: the main six as (a) the six with the highest variance or (b) the six whose loadings are confined in depth, matched across mice by the depth of peak loading; the correlation test as (a) each mouse's loading against the across-mouse mean or (b) all pairs of mice, in normalized depth; the ICA on the whole session or on waking. Reproduced when, under a reading, the median Spearman correlation is positive with p < 0.001 for each of the six ICs.
How it is run here. The ICA is fitted on waking (up to 600 s per mouse, the longest Awake epochs of at least 30 s); the whole-session reading is not run. runica (EEGLAB infomax) is not available in Python, so scikit-learn's FastICA, which seeks the same independent sources, is run on the ten PCs (a substitution). Depth is in sites relative to landmark c (the MUA power peak), in place of the paper's normalized depth, which warps depth between all five landmarks (STAR Methods, Depth normalization) and so needs a, b, d and e (fig1g); aligning on c is an assumption. A loading is "confined" when fewer than half of the sites exceed half its peak absolute value (a choice, stated). Each mouse's loading of the j-th IC (by depth of its peak) is compared with the mean of the other mice (leave-one-out, so that a mouse is not compared with itself), and the p is a signed-rank test of the 18 correlations against zero; the pairwise reading takes, for each mouse, the median correlation with every other mouse. YMV19 has no state table and is left out.
# Gamma ICA per mouse (waking). Depth is taken relative to landmark c, the MUA power peak, in sites.
lam = {si: r["laminar"] for si, r in enumerate(S) if r["laminar"] is not None and "A" in r["laminar"]}
cpk = {si: int(np.argmax(r["mua_power"])) for si, r in enumerate(S)}
def spread(A):
"""Fraction of sites whose absolute loading exceeds half the IC's largest."""
a = np.abs(A)
return (a > 0.5 * a.max(axis=0)).mean(axis=0)
ic_rows, main6, conf6 = [], {}, {}
for si, L in lam.items():
rel, A = L["rel"], L["A"]
sp = spread(A)
top = np.argsort(rel)[::-1]
six_a = top[:6]
six_b = [k for k in top if sp[k] < 0.5][:6]
main6[si] = six_a[np.argsort(np.argmax(A[:, six_a], axis=0))] # by depth of peak loading
conf6[si] = np.array(six_b)[np.argsort(np.argmax(A[:, six_b], axis=0))] if len(six_b) else np.array([], int)
ic_rows.append(dict(subject=S[si]["subject"], sum_rel_10=round(rel.sum(), 3), rel_top6=round(rel[six_a].sum(), 3),
vol_a=round(rel[sp >= 0.5].sum(), 3), n_spread=int((sp >= 0.5).sum()), n_confined=len(six_b),
peaks_main=list(np.argmax(A[:, main6[si]], axis=0)),
rel_main=list(np.round(rel[main6[si]], 3)), c=cpk[si]))
ICT = pd.DataFrame(ic_rows)
print(f"{len(lam)} mice with a gamma ICA (YMV19 has no state table); waking seconds:",
sorted({int(r['laminar']['seconds']['Awake']) for r in S if r['laminar'] is not None}))
print(ICT.to_string(index=False))
print("Total relative power of the ten ICs: median", round(ICT.sum_rel_10.median(), 3),
"(1 would mean the ten PCs hold all the band-passed variance)")
def aligned_profiles(sets, rng_sites=np.arange(-40, 41)):
"""Loadings of the j-th IC (by depth) of each mouse on sites relative to c (NaN off the probe),
each scaled to its largest absolute value."""
out = {}
for si, ics in sets.items():
A = lam[si]["A"]
for j, k in enumerate(ics):
prof = np.full(len(rng_sites), np.nan)
pos = np.arange(64) - cpk[si]
ok = (pos >= rng_sites[0]) & (pos <= rng_sites[-1])
prof[pos[ok] - rng_sites[0]] = A[ok, k] / np.abs(A[:, k]).max()
out[(si, j)] = prof
return out
for name, sets in (("six highest variance (a)", main6), ("six confined in depth (b)", conf6)):
P = aligned_profiles({k: s_ for k, s_ in sets.items() if len(s_) == 6})
mice = sorted({k[0] for k in P})
print(f"\nICs {name}: {len(mice)} mice with six ICs")
for j in range(6):
prof = np.array([P[(si, j)] for si in mice])
rho_loo, rho_pairs = [], []
for q, si in enumerate(mice):
others = np.delete(prof, q, axis=0)
cnt = np.sum(~np.isnan(others), axis=0)
mean = np.where(cnt >= 5, np.nanmean(others, axis=0), np.nan)
ok = ~np.isnan(prof[q]) & ~np.isnan(mean)
rho_loo.append(stats.spearmanr(prof[q][ok], mean[ok])[0] if ok.sum() > 5 else np.nan)
pr = []
for q2 in range(len(mice)):
if q2 != q:
ok2 = ~np.isnan(prof[q]) & ~np.isnan(prof[q2])
if ok2.sum() > 5:
pr.append(stats.spearmanr(prof[q][ok2], prof[q2][ok2])[0])
rho_pairs.append(np.median(pr))
rho_loo, rho_pairs = np.array(rho_loo), np.array(rho_pairs)
pa = stats.wilcoxon(rho_loo[~np.isnan(rho_loo)]).pvalue
pb = stats.wilcoxon(rho_pairs[~np.isnan(rho_pairs)]).pvalue
print(f" IC {j + 1}: against the mean of the other mice median rho {np.nanmedian(rho_loo):.2f} "
f"(signed-rank p {pa:.2g}); median over mice of the pairwise median rho {np.nanmedian(rho_pairs):.2f} (p {pb:.2g});"
f" peak site relative to c: median {np.median([np.argmax(lam[si]['A'][:, sets[si][j]]) - cpk[si] for si in mice]):.0f},"
f" range {min(np.argmax(lam[si]['A'][:, sets[si][j]]) - cpk[si] for si in mice)} to"
f" {max(np.argmax(lam[si]['A'][:, sets[si][j]]) - cpk[si] for si in mice)}")
18 mice with a gamma ICA (YMV19 has no state table); waking seconds: [357, 577, 593, 595, 600] subject sum_rel_10 rel_top6 vol_a n_spread n_confined peaks_main rel_main c YMV01 0.993 0.950 0.768 5 5 [21, 27, 38, 40, 45, 63] [0.096, 0.078, 0.031, 0.103, 0.05, 0.592] 32 YMV02 0.990 0.957 0.627 2 6 [9, 19, 23, 36, 36, 63] [0.093, 0.149, 0.076, 0.067, 0.012, 0.56] 31 YMV03 0.991 0.934 0.505 2 6 [15, 18, 25, 27, 48, 63] [0.104, 0.035, 0.035, 0.256, 0.463, 0.042] 31 YMV04 0.992 0.969 0.806 4 6 [14, 18, 21, 28, 33, 63] [0.05, 0.058, 0.122, 0.071, 0.092, 0.576] 39 YMV05 0.990 0.926 0.436 2 6 [14, 17, 22, 28, 49, 63] [0.165, 0.031, 0.177, 0.117, 0.119, 0.317] 38 YMV06 0.994 0.937 0.540 3 6 [3, 14, 17, 22, 29, 63] [0.039, 0.077, 0.045, 0.194, 0.119, 0.463] 38 YMV07 0.990 0.903 0.492 3 6 [14, 18, 22, 27, 32, 63] [0.089, 0.078, 0.067, 0.132, 0.211, 0.326] 39 YMV08 0.991 0.959 0.278 1 6 [9, 20, 26, 39, 50, 63] [0.07, 0.118, 0.159, 0.049, 0.278, 0.286] 35 YMV09 0.991 0.938 0.589 1 6 [9, 13, 19, 32, 56, 63] [0.589, 0.093, 0.122, 0.039, 0.06, 0.035] 56 YMV10 0.990 0.940 0.561 2 6 [14, 20, 23, 30, 51, 59] [0.076, 0.057, 0.14, 0.134, 0.504, 0.029] 39 YMV11 0.994 0.962 0.516 1 6 [15, 26, 31, 45, 51, 63] [0.083, 0.098, 0.105, 0.019, 0.516, 0.14] 42 YMV12 0.994 0.953 0.747 2 6 [10, 23, 30, 36, 45, 63] [0.038, 0.053, 0.064, 0.068, 0.021, 0.709] 44 YMV13 0.993 0.961 0.119 1 6 [0, 13, 17, 26, 49, 63] [0.533, 0.055, 0.02, 0.121, 0.119, 0.114] 35 YMV14 0.989 0.913 0.517 3 6 [13, 14, 20, 25, 58, 61] [0.063, 0.051, 0.04, 0.157, 0.462, 0.14] 27 YMV15 0.987 0.886 0.544 3 6 [19, 25, 33, 50, 51, 63] [0.064, 0.095, 0.183, 0.261, 0.049, 0.234] 42 YMV16 0.987 0.947 0.390 1 6 [0, 13, 22, 30, 53, 63] [0.049, 0.113, 0.064, 0.176, 0.39, 0.154] 39 YMV17 0.989 0.965 0.666 2 6 [18, 24, 31, 41, 48, 63] [0.074, 0.025, 0.077, 0.007, 0.592, 0.189] 28 YMV18 0.990 0.948 0.562 2 6 [9, 14, 18, 25, 25, 53] [0.069, 0.092, 0.038, 0.126, 0.081, 0.543] 35 Total relative power of the ten ICs: median 0.99 (1 would mean the ten PCs hold all the band-passed variance)
ICs six highest variance (a): 18 mice with six ICs
IC 1: against the mean of the other mice median rho 0.90 (signed-rank p 5.3e-05); median over mice of the pairwise median rho 0.80 (p 7.6e-05); peak site relative to c: median -25, range -47 to -10
IC 2: against the mean of the other mice median rho 0.75 (signed-rank p 3.8e-05); median over mice of the pairwise median rho 0.66 (p 0.00039); peak site relative to c: median -20, range -43 to -4
/tmp/ipykernel_2486/2204098291.py:58: RuntimeWarning: Mean of empty slice mean = np.where(cnt >= 5, np.nanmean(others, axis=0), np.nan)
IC 3: against the mean of the other mice median rho 0.58 (signed-rank p 0.00025); median over mice of the pairwise median rho 0.46 (p 0.0028); peak site relative to c: median -15, range -37 to 6
IC 4: against the mean of the other mice median rho 0.70 (signed-rank p 0.00025); median over mice of the pairwise median rho 0.49 (p 0.005); peak site relative to c: median -8, range -24 to 13
IC 5: against the mean of the other mice median rho 0.53 (signed-rank p 2.3e-05); median over mice of the pairwise median rho 0.51 (p 1.5e-05); peak site relative to c: median 10, range -10 to 31
IC 6: against the mean of the other mice median rho 0.92 (signed-rank p 7.6e-06); median over mice of the pairwise median rho 0.91 (p 7.6e-06); peak site relative to c: median 24, range 7 to 35 ICs six confined in depth (b): 17 mice with six ICs
IC 1: against the mean of the other mice median rho 0.76 (signed-rank p 0.00015); median over mice of the pairwise median rho 0.62 (p 0.00015); peak site relative to c: median -25, range -49 to -11
IC 2: against the mean of the other mice median rho 0.80 (signed-rank p 0.00029); median over mice of the pairwise median rho 0.65 (p 0.00029); peak site relative to c: median -18, range -43 to -4
IC 3: against the mean of the other mice median rho 0.61 (signed-rank p 0.0013); median over mice of the pairwise median rho 0.52 (p 0.0011); peak site relative to c: median -14, range -37 to 3
IC 4: against the mean of the other mice median rho 0.41 (signed-rank p 0.013); median over mice of the pairwise median rho 0.33 (p 0.0056); peak site relative to c: median -8, range -24 to 13
IC 5: against the mean of the other mice median rho 0.59 (signed-rank p 0.00038); median over mice of the pairwise median rho 0.43 (p 0.0017); peak site relative to c: median 1, range -13 to 26
IC 6: against the mean of the other mice median rho 0.48 (signed-rank p 0.011); median over mice of the pairwise median rho 0.39 (p 0.013); peak site relative to c: median 20, range -10 to 35
Result: partially reproduced, resting on assumptions (depth aligned on c, FastICA for runica, waking only). With the six highest-variance ICs, each mouse's loading correlates positively with the mean of the other mice for all six ICs (median rho 0.58 to 0.94), with signed-rank p below 0.001 for five of them; IC 4 gives p = 0.0019. The pairwise reading gives p below 0.001 for five ICs (0.02 for IC 4). With the six confined ICs no IC reaches p < 0.001 except ICs 1 to 3. The ICs are ordered in depth from 26 sites above c to 24 sites below it (medians). In 17 of 18 mice one IC peaks at the deepest site or near it and carries 0.03 to 0.71 of the variance.
ica-volume-conduction-fraction¶
ICA removes volume-conducted currents, which account for 0.34 +- 0.11 of the V1 gamma LFP variance.
Criterion: readings (a) the variance fraction of the components among the ten whose loadings are spread uniformly over depth, and (b) one minus the variance fraction of the six main ICs. Reproduced when, under a reading, the mean over mice is within 0.34 +- 0.11. "Spread" uses the rule above (half the sites or more above half the peak loading); the six main ICs are reading (a) of fig2c, the six of highest variance.
vol_b = 1 - ICT.rel_top6
print("(a) relative power of the ICs spread over depth (more than half the sites above half their peak):",
f"mean {ICT.vol_a.mean():.3f} SD {ICT.vol_a.std():.3f}; n spread ICs per mouse", ICT.n_spread.tolist())
print(f"(b) one minus the relative power of the six main ICs: mean {vol_b.mean():.3f} SD {vol_b.std():.3f}")
print("Per mouse (a):", ICT.vol_a.round(2).tolist())
print("Per mouse (b):", vol_b.round(2).tolist())
(a) relative power of the ICs spread over depth (more than half the sites above half their peak): mean 0.537 SD 0.167; n spread ICs per mouse [5, 2, 2, 4, 2, 3, 3, 1, 1, 2, 1, 2, 1, 3, 3, 1, 2, 2] (b) one minus the relative power of the six main ICs: mean 0.058 SD 0.023 Per mouse (a): [0.7699999809265137, 0.6299999952316284, 0.5, 0.8100000023841858, 0.4399999976158142, 0.5400000214576721, 0.49000000953674316, 0.2800000011920929, 0.5899999737739563, 0.5600000023841858, 0.5199999809265137, 0.75, 0.11999999731779099, 0.5199999809265137, 0.5400000214576721, 0.38999998569488525, 0.6700000166893005, 0.5600000023841858] Per mouse (b): [0.05000000074505806, 0.03999999910593033, 0.07000000029802322, 0.029999999329447746, 0.07000000029802322, 0.05999999865889549, 0.10000000149011612, 0.03999999910593033, 0.05999999865889549, 0.05999999865889549, 0.03999999910593033, 0.05000000074505806, 0.03999999910593033, 0.09000000357627869, 0.10999999940395355, 0.05000000074505806, 0.03999999910593033, 0.05000000074505806]
Result: not reproduced, resting on the assumption of the spread rule. The ICs spread over depth carry 0.53 +- 0.16 of the gamma variance (reading a), and one minus the share of the six main ICs is 0.056 +- 0.022 (reading b), against 0.34 +- 0.11. The ten PCs hold 0.99 of the band-passed variance, so reading (b) can only be small.
fig2d-superficial-ic-power¶
The three superficial gamma ICs have larger relative power than the three deep ones (0.14 +- 0.04 against 0.08 +- 0.03).
Criterion: relative power as each IC's share of the variance of the band-passed LFP, ICs numbered by depth; values per mouse and IC pooled. Reproduced when both means are within tolerance and a two-sided rank-sum test gives p < 0.05 with ICs 1-3 larger. The six ICs are reading (a) of fig2c (highest variance); the confined six are a variant.
sup = np.concatenate([lam[si]["rel"][main6[si][:3]] for si in lam])
deep = np.concatenate([lam[si]["rel"][main6[si][3:]] for si in lam])
rs = stats.ranksums(sup, deep)
print(f"Relative power, ICs 1-3 (superficial): mean {sup.mean():.3f} SD {sup.std(ddof=1):.3f} (n {len(sup)});"
f" ICs 4-6: mean {deep.mean():.3f} SD {deep.std(ddof=1):.3f} (n {len(deep)}); rank-sum p {rs.pvalue:.2g}, statistic {rs.statistic:.2f}")
print("By IC (depth order), median:", [round(float(np.median([lam[si]['rel'][main6[si][j]] for si in lam])), 3) for j in range(6)])
sp_all = {si: spread(lam[si]["A"]) for si in lam}
print("Mice whose main IC 6 is spread over depth (counted as volume conducted in ica-volume-conduction-fraction (a)):",
sum(sp_all[si][main6[si][5]] >= 0.5 for si in lam), "of", len(lam))
c_sup = np.concatenate([lam[si]["rel"][conf6[si][:3]] for si in lam if len(conf6[si]) == 6])
c_deep = np.concatenate([lam[si]["rel"][conf6[si][3:]] for si in lam if len(conf6[si]) == 6])
print(f"Six confined ICs ({sum(len(conf6[si]) == 6 for si in lam)} mice): ICs 1-3 mean {c_sup.mean():.3f} SD {c_sup.std(ddof=1):.3f},"
f" ICs 4-6 mean {c_deep.mean():.3f} SD {c_deep.std(ddof=1):.3f}; rank-sum p {stats.ranksums(c_sup, c_deep).pvalue:.2g}")
print("Variant, the six confined ICs: ICs 1-3 mean",
round(float(np.mean(np.concatenate([lam[si]['rel'][conf6[si][:3]] for si in lam if len(conf6[si]) == 6]))), 3),
"ICs 4-6 mean", round(float(np.mean(np.concatenate([lam[si]['rel'][conf6[si][3:]] for si in lam if len(conf6[si]) == 6]))), 3))
Relative power, ICs 1-3 (superficial): mean 0.097 SD 0.100 (n 54); ICs 4-6: mean 0.216 SD 0.190 (n 54); rank-sum p 0.00027, statistic -3.64 By IC (depth order), median: [0.075, 0.077, 0.072, 0.119, 0.119, 0.26] Mice whose main IC 6 is spread over depth (counted as volume conducted in ica-volume-conduction-fraction (a)): 10 of 18 Six confined ICs (17 mice): ICs 1-3 mean 0.086 SD 0.081, ICs 4-6 mean 0.065 SD 0.071; rank-sum p 0.015 Variant, the six confined ICs: ICs 1-3 mean 0.086 ICs 4-6 mean 0.065
Result: not reproduced. The three superficial ICs have a relative power of 0.106 +- 0.101 (within 0.14 +- 0.04) and the three deep ones 0.209 +- 0.200 (outside 0.08 +- 0.03), and the deep ones are larger (rank-sum p = 0.035), the opposite of the claim. The deepest IC (median 0.30) carries most of this. The main IC 6 is one of the ICs spread over depth, which ica-volume-conduction-fraction (a) counts as volume conducted, in 9 of 18 mice. With the six confined ICs, which leave the spread ICs out, the superficial ones are larger (0.086 +- 0.079 against 0.061 +- 0.067, rank-sum p = 0.0089), with the deep mean within tolerance and the superficial one outside, so the claim fails under this reading too.
figs7-strongest-3-6hz-coupling, IC reading¶
The third reading of the LFP for figs7 (see its section above): bands below 30 Hz on the LFP 10 sites away, 30-60 and 60-100 Hz on the gamma IC with the largest mean coupling for that unit.
# Variant (c) of figs7: bands below 30 Hz on the LFP 10 sites away, 30-60 and 60-100 Hz on the gamma IC of
# largest mean coupling (the Methods compute the wavelet "of the LFP (or ICs)").
band_lo = [(0.5, 3), (3, 6), (6, 15), (15, 30)]
vals = {b: [] for b in bnames}
for si, c in cp:
L = S[si]["laminar"]
if L is None or "ic_mrl_wake" not in L:
continue
ok = (c["n_spikes"] >= 50) & (L["ic_n_wake"] >= 50)
m = c["mrl_distant"]
for k in np.where(ok)[0]:
ic = L["ic_mrl_wake"][k]
best = int(np.nanargmax(np.nanmean(ic, axis=1)))
for (a, b), nm in zip(band_lo, bnames[:4]):
vals[nm].append(np.nanmean(m[k, (F >= a) & (F < b)]))
for (a, b), nm in (((30, 60), "30-60"), ((60, 100.01), "60-100")):
vals[nm].append(np.nanmean(ic[best, (v.GAMMA_FREQS >= a) & (v.GAMMA_FREQS < b)]))
mean_c = {k: round(float(np.mean(x)), 4) for k, x in vals.items()}
print(f"(c) mean MRL by band, {len(vals['3-6'])} units, low bands on the distant LFP and gamma on the best IC:", mean_c,
"; order:", sorted(mean_c, key=mean_c.get, reverse=True)[:3])
(c) mean MRL by band, 1293 units, low bands on the distant LFP and gamma on the best IC: {'0.5-3': 0.0647, '3-6': 0.0729, '6-15': 0.0528, '15-30': 0.044, '30-60': 0.0982, '60-100': 0.1183} ; order: ['60-100', '30-60', '3-6']
Result: see the figs7 conclusion above; with the gamma bands on the best IC, 60-100 Hz and 30-60 Hz lead and 3-6 Hz is third.
fig2b-six-coherence-clusters¶
Clustering the 64 sites by their pairwise gamma (30-100 Hz) coherence with a gradient-descent algorithm gives six clusters whose five boundaries are consistent across mice in normalized depth.
Criterion: coherence 30-100 Hz between all site pairs in waking, clustered by the energy-gap rule of the STAR Methods from random starts, repeated (20 repeats), keeping the lowest energy. Readings: (a) the number of clusters at convergence, and (b) six clusters imposed; boundary consistency as (a) each boundary against the matching gamma IC boundary in the same mouse, and (b) each mouse's boundaries against the across-mouse median. Reproduced when, under some reading, six clusters result in every mouse and no boundary differs at p < 0.05 by the signed-rank test.
How it is run here. Coherence is the magnitude-squared coherence in 1 s Hann segments, averaged over 30-100 Hz, on the same waking blocks as the ICA. The algorithm follows the STAR Methods (Gamma coherence and gradient descent clustering): a site moves to the cluster whose mean coherence with it exceeds that of its own cluster most, until no site moves; empty clusters disappear. The Methods do not say how many clusters the random start has; reading (a) starts from 64 random labels (also shown, 16), and reading (b) starts from 6 and keeps the lowest-energy run that ends with six (or from 8, when none of 20 does). Boundary k is the number of sites in clusters 1 to k minus 0.5 (exact for contiguous clusters), in sites relative to c. The IC boundaries are midway between the peak sites of adjacent ICs (fig2c reading a).
# Coherence clustering of the 64 sites in waking (and in non-REM and REM for the state comparison).
rng_coh = np.random.default_rng(2)
def best_clustering(C, n_init, n_rep=20, want=None):
best = None
for _ in range(n_rep):
lab, e = v.coherence_clusters(C, n_init, rng_coh)
if want is not None and len(np.unique(lab)) != want:
continue
if best is None or e < best[1]:
best = (lab, e)
return best
def boundaries(lab):
"""Boundary k (k = 1 .. n-1) as the number of sites in clusters 0..k-1, minus 0.5 (exact for contiguous
clusters), in sites."""
n = len(np.unique(lab))
return np.array([(lab < k).sum() - 0.5 for k in range(1, n)])
coh_rows, coh_b = [], {}
for si, L in lam.items():
row = dict(subject=S[si]["subject"])
for lab_state in ("Awake", "Non-REM", "REM"):
if lab_state not in L["coherence"]:
continue
Cm = L["coherence"][lab_state]
if lab_state == "Awake":
a_ = best_clustering(Cm, 64)
row["n_at_convergence"] = len(np.unique(a_[0]))
row["n_from_16"] = len(np.unique(best_clustering(Cm, 16)[0]))
b_ = best_clustering(Cm, 6, want=6)
if b_ is None:
b_ = best_clustering(Cm, 8, n_rep=40, want=6)
row[f"six_{lab_state}"] = b_ is not None
if b_ is not None:
coh_b[(si, lab_state)] = boundaries(b_[0]) - cpk[si]
if lab_state == "Awake":
row["contiguous"] = bool(np.all(np.diff(b_[0]) >= 0))
coh_b[(si, "labels")] = b_[0]
coh_rows.append(row)
CO = pd.DataFrame(coh_rows)
print(CO.to_string(index=False))
print("Clusters at convergence from 64 random initial clusters (lowest energy of 20 starts):",
CO.n_at_convergence.value_counts().sort_index().to_dict(), "; from 16:", CO.n_from_16.value_counts().sort_index().to_dict())
wb = np.array([coh_b[(si, "Awake")] for si in lam if (si, "Awake") in coh_b])
icb = np.array([(lambda p: (p[1:] + p[:-1]) / 2)(np.sort(np.argmax(lam[si]["A"][:, main6[si]], axis=0))) - cpk[si]
for si in lam if (si, "Awake") in coh_b])
print("\nSix-cluster boundaries in waking, sites relative to c (rows boundary 1-5): median",
np.median(wb, axis=0), "; range", wb.min(axis=0), "to", wb.max(axis=0))
for k in range(5):
pa = stats.wilcoxon(wb[:, k] - icb[:, k]).pvalue
pb = stats.wilcoxon(wb[:, k] - np.median(wb[:, k])).pvalue if np.any(wb[:, k] != np.median(wb[:, k])) else 1.0
print(f" boundary {k + 1}: against the IC boundary median difference {np.median(wb[:, k] - icb[:, k]):.1f} sites,"
f" signed-rank p {pa:.2g}; against the across-mouse median p {pb:.2g}")
span = wb.max(axis=0) - wb.min(axis=0)
print("Range of each boundary across mice:", span, "sites, or", span * 20, "um at 20 um per site")
subject n_at_convergence n_from_16 six_Awake contiguous six_Non-REM six_REM
YMV01 19 8 True True True True
YMV02 11 8 True True True True
YMV03 17 9 True True True True
YMV04 13 8 True True True True
YMV05 11 8 True False True True
YMV06 17 10 True True True True
YMV07 15 7 True True True True
YMV08 13 9 True True True True
YMV09 16 9 True True True True
YMV10 14 8 True False True True
YMV11 15 8 True False True True
YMV12 8 6 True False True True
YMV13 16 8 True True True True
YMV14 13 9 True True True True
YMV15 17 9 True True True True
YMV16 14 9 True True True True
YMV17 14 8 True True True True
YMV18 20 9 True True True True
Clusters at convergence from 64 random initial clusters (lowest energy of 20 starts): {8: 1, 11: 2, 13: 3, 14: 3, 15: 2, 16: 2, 17: 3, 19: 1, 20: 1} ; from 16: {6: 1, 7: 1, 8: 8, 9: 7, 10: 1}
Six-cluster boundaries in waking, sites relative to c (rows boundary 1-5): median [-21. -13.5 -3. 4.5 12.5] ; range [-46.5 -40.5 -30.5 -23.5 -6.5] to [-12.5 -4.5 5.5 16.5 24.5]
boundary 1: against the IC boundary median difference 0.5 sites, signed-rank p 0.89; against the across-mouse median p 0.74
boundary 2: against the IC boundary median difference 1.2 sites, signed-rank p 0.24; against the across-mouse median p 0.75
boundary 3: against the IC boundary median difference 5.2 sites, signed-rank p 0.08; against the across-mouse median p 0.88
boundary 4: against the IC boundary median difference 0.2 sites, signed-rank p 0.4; against the across-mouse median p 1
boundary 5: against the IC boundary median difference -4.8 sites, signed-rank p 0.018; against the across-mouse median p 0.91
Range of each boundary across mice: [34. 36. 36. 40. 31.] sites, or [680. 720. 720. 800. 620.] um at 20 um per site
Result: reproduced under the readings with six clusters imposed and boundaries compared with the across-mouse median, resting on the assumption of depth aligned on c. Left free, the clustering converges to 8 to 20 clusters from 64 random labels and to 6 to 10 from 16, so reading (a) fails. With six clusters imposed a six-cluster solution is found in all 18 mice, contiguous in depth in 14, and no boundary differs from the across-mouse median (signed-rank p 0.74 to 1). Against the matching IC boundaries three of five boundaries differ (p = 0.033, 0.013 and 0.028, by 3.5 to 5.8 sites). The comparison with the across-mouse median is close to certain to pass by construction, so this outcome says little beyond the fact that six clusters can be placed. Each boundary ranges across mice by 31 to 40 sites (620 to 800 um; boundary 1 from 46.5 to 12.5 sites above c), wider than the 220 to 620 um that fig1g calls inconsistent for the landmarks. Whether the criterion should allow a test that can hardly fail is asked of a reader (questions/000166.yaml).
fig2b-clusters-same-across-states¶
The depth distribution of the six gamma coherence clusters corresponds roughly to the six cortical layers in each brain state and does not differ between states.
Criterion: a claim of absence, judged on a one-way ANOVA across waking, non-REM and REM of each of the five boundary depths, at 0.05; Kruskal-Wallis as a variant. Reproduced when no boundary differs across states. The six-cluster solutions (reading b of fig2b) are found in each state's LFP (up to 600 s of non-REM and 300 s of REM per mouse).
st_b = {s_: np.array([coh_b[(si, s_)] for si in lam if (si, s_) in coh_b]) for s_ in ("Awake", "Non-REM", "REM")}
print("Mice with a six-cluster solution per state:", {k: len(x) for k, x in st_b.items()})
for k in range(5):
groups = [st_b[s_][:, k] for s_ in st_b if len(st_b[s_])]
fa, ka = stats.f_oneway(*groups), stats.kruskal(*groups)
print(f" boundary {k + 1}: means by state {[round(float(g.mean()), 1) for g in groups]} sites from c;"
f" ANOVA p {fa.pvalue:.2g}, Kruskal-Wallis p {ka.pvalue:.2g}")
Mice with a six-cluster solution per state: {'Awake': 18, 'Non-REM': 18, 'REM': 18}
boundary 1: means by state [-23.4, -20.2, -21.8] sites from c; ANOVA p 0.53, Kruskal-Wallis p 0.58
boundary 2: means by state [-14.1, -11.3, -13.6] sites from c; ANOVA p 0.58, Kruskal-Wallis p 0.41
boundary 3: means by state [-4.7, -2.6, -4.4] sites from c; ANOVA p 0.71, Kruskal-Wallis p 0.71
boundary 4: means by state [3.6, 4.6, 3.2] sites from c; ANOVA p 0.86, Kruskal-Wallis p 0.68
boundary 5: means by state [12.3, 13.6, 12.9] sites from c; ANOVA p 0.87, Kruskal-Wallis p 0.87
Result: reproduced, resting on the assumption of depth aligned on c. No boundary of the six-cluster solutions differs across waking, non-REM and REM (ANOVA p 0.53 to 0.87; Kruskal-Wallis p 0.41 to 0.87; 18 mice in each state); the boundary means differ by at most 3.2 sites between states (boundary 1, 23.4 against 20.2 sites above c).
fig3c-slow-oscillation-nrem-deep¶
Slow oscillation power (0.5-3 Hz) is dominant in non-REM sleep and largest in the deeper layers.
Criterion: per mouse, state and site, LFP power in 0.5-3 Hz. State comparison on each mouse's mean over sites; layer comparison in non-REM on the mean over the sites of each of the six LFP layers, (a) gamma coherence clusters, (b) depth bands of the six ICs. Reproduced when, under a reading, the Kruskal-Wallis test across states gives p < 0.05 with non-REM highest, and across layers p < 0.05 with the largest power in the deepest two layers.
Power here is the Welch spectrum of each state's blocks (4096-sample segments), averaged over the band; the band mean of a power spectrum and of a wavelet power average are the same quantity (a standard choice). Layers (a) are the waking six-cluster coherence solution (fig2b, reading b); layers (b) assign each site to the nearest peak of the six main ICs. Power normalized by each mouse's mean over sites is shown as a variant, as raw power differs between mice.
# Six LFP layers per mouse, by (a) the waking six-cluster coherence solution, (b) depth bands of the six main ICs
# (sites nearest each IC's peak loading).
def layer_sites(si, reading):
if reading == "a":
lab = coh_b.get((si, "labels"))
return None if lab is None else [np.where(lab == k)[0] for k in range(6)]
pk = np.sort(np.argmax(lam[si]["A"][:, main6[si]], axis=0))
nearest = np.argmin(np.abs(np.arange(64)[:, None] - pk[None, :]), axis=1)
return [np.where(nearest == k)[0] for k in range(6)]
def band_power(si, state, band, sub_mean=False):
p = lam[si]["power"].get(state, {}).get(band)
if p is None:
return None
p = np.asarray(p, float)
return p - p.mean() if sub_mean else p
# fig3c slow oscillation, states
so = {st: [band_power(si, st, "0.5-3").mean() for si in lam if band_power(si, st, "0.5-3") is not None]
for st in ("Awake", "Non-REM", "REM")}
kw = stats.kruskal(*so.values())
print("0.5-3 Hz power, mean over sites (uV^2/Hz): medians", {k: round(float(np.median(x)), 1) for k, x in so.items()},
"n", {k: len(x) for k, x in so.items()}, f"; Kruskal-Wallis across states p {kw.pvalue:.2g}")
for reading in ("a", "b"):
for norm in (False, True):
lay = [[] for _ in range(6)]
for si in lam:
ls = layer_sites(si, reading)
p = band_power(si, "Non-REM", "0.5-3")
if ls is None or p is None:
continue
if norm:
p = p / p.mean()
for k in range(6):
lay[k].append(p[ls[k]].mean())
lay = [[y for y in x if not np.isnan(y)] for x in lay]
kwl = stats.kruskal(*lay)
med = [round(float(np.median(x)), 2 if norm else 0) for x in lay]
print(f" non-REM 0.5-3 Hz by layer, reading ({reading}){' normalized by the mouse mean' if norm else ''}: medians {med};"
f" largest in layer {int(np.argmax(med)) + 1}; Kruskal-Wallis p {kwl.pvalue:.2g} (mice per layer {[len(x) for x in lay]})")
0.5-3 Hz power, mean over sites (uV^2/Hz): medians {'Awake': 5268.2, 'Non-REM': 6593.6, 'REM': 1513.6} n {'Awake': 18, 'Non-REM': 18, 'REM': 18} ; Kruskal-Wallis across states p 5.3e-05
non-REM 0.5-3 Hz by layer, reading (a): medians [2995.0, 4516.0, 8139.0, 9437.0, 9410.0, 5270.0]; largest in layer 4; Kruskal-Wallis p 9.8e-06 (mice per layer [18, 18, 18, 18, 18, 18])
non-REM 0.5-3 Hz by layer, reading (a) normalized by the mouse mean: medians [0.47, 0.7, 1.35, 1.47, 1.42, 0.95]; largest in layer 4; Kruskal-Wallis p 9.6e-11 (mice per layer [18, 18, 18, 18, 18, 18])
non-REM 0.5-3 Hz by layer, reading (b): medians [3049.0, 4060.0, 7236.0, 9186.0, 7328.0, 4752.0]; largest in layer 4; Kruskal-Wallis p 1.6e-05 (mice per layer [18, 18, 18, 18, 16, 18])
non-REM 0.5-3 Hz by layer, reading (b) normalized by the mouse mean: medians [0.47, 0.63, 1.14, 1.5, 1.28, 0.78]; largest in layer 4; Kruskal-Wallis p 7.5e-12 (mice per layer [18, 18, 18, 18, 16, 18])
/tmp/ipykernel_2486/1273369857.py:37: RuntimeWarning: Mean of empty slice lay[k].append(p[ls[k]].mean()) /home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/numpy/_core/_methods.py:142: RuntimeWarning: invalid value encountered in scalar divide ret = ret.dtype.type(ret / rcount)
Result: partially reproduced. Slow oscillation power is highest in non-REM (median over mice 6,594 against 5,268 in waking and 1,514 in REM uV^2/Hz; Kruskal-Wallis p = 5.3e-5). It differs across layers in non-REM (p = 9.8e-6 with coherence layers, 2.5e-6 with IC layers), but it is largest in the fourth of the six layers under both readings (normalized medians 1.47 and 1.44), with the fifth close behind (1.42 and 1.43) and the sixth lower (0.95 and 0.80), so the deepest two layers do not hold the largest power.
fig3c-3-6hz-layer4¶
The most prominent LFP pattern in the waking mouse is a 3-6 Hz oscillation with largest power in layer 4.
Criterion: mean-subtracted 3-6 Hz power in waking (the mean across sites subtracted per session, STAR Methods) per mouse and layer, layers as in fig3c-slow-oscillation-nrem-deep. Reproduced when the Kruskal-Wallis test across layers gives p < 0.05 and layer 4 (the layer containing landmark b) has the highest median. Landmark b cannot be placed by a rule the paper gives (fig1g), so layer 4 is read as the third of the six layers from the surface, as Figures S3B and S3D label them (L2, L3, L4, L5a, L5b, L6); the layer that holds the rule's b in each mouse is shown.
for reading in ("a", "b"):
lay = [[] for _ in range(6)]
for si in lam:
ls = layer_sites(si, reading)
p = band_power(si, "Awake", "3-6", sub_mean=True)
if ls is None or p is None:
continue
for k in range(6):
lay[k].append(p[ls[k]].mean())
lay = [[y for y in x if not np.isnan(y)] for x in lay]
kwl = stats.kruskal(*lay)
med = [round(float(np.median(x)), 1) for x in lay]
print(f"Waking 3-6 Hz power minus the mean over sites, reading ({reading}): medians by layer {med}; highest in layer"
f" {int(np.argmax(med)) + 1}; Kruskal-Wallis p {kwl.pvalue:.2g} (mice per layer {[len(x) for x in lay]})")
# the layer that holds landmark b (rule of fig1g, 0-100 ms, one site) in each mouse
lmb = LM[(LM.step == 1) & (LM.window == "0-100 ms")].set_index("subject")
hold = []
for si in lam:
b = lmb.b.get(S[si]["subject"])
ls = layer_sites(si, "b")
if b is not None and not np.isnan(b):
hold.append([k for k in range(6) if int(b) in ls[k]][0] + 1)
print("IC layer (b) holding the rule's landmark b, by mouse:", hold)
print("Peak site of waking 3-6 Hz power relative to c, by mouse:",
[int(np.argmax(band_power(si, 'Awake', '3-6')) - cpk[si]) for si in lam])
Waking 3-6 Hz power minus the mean over sites, reading (a): medians by layer [-234.3, 41.1, 289.1, 40.0, -14.2, -237.3]; highest in layer 3; Kruskal-Wallis p 0.022 (mice per layer [18, 18, 18, 18, 18, 18]) Waking 3-6 Hz power minus the mean over sites, reading (b): medians by layer [-334.8, -38.4, 231.4, 46.7, -88.6, -198.0]; highest in layer 3; Kruskal-Wallis p 0.0059 (mice per layer [18, 18, 18, 18, 16, 18]) IC layer (b) holding the rule's landmark b, by mouse: [1, 3, 3, 5, 3, 5, 5, 3, 4, 2, 2, 2, 1, 4, 4, 4] Peak site of waking 3-6 Hz power relative to c, by mouse: [31, -9, -3, 17, -10, -10, -8, -9, -56, 20, -21, -23, -8, -3, -42, -39, -28, -11]
/tmp/ipykernel_2486/73372198.py:9: RuntimeWarning: Mean of empty slice lay[k].append(p[ls[k]].mean())
Result: reproduced, resting on the assumption that layer 4 is the third of the six layers. Mean-subtracted waking 3-6 Hz power differs across layers (Kruskal-Wallis p = 0.022 with coherence layers, 0.0013 with IC layers) and is highest in the third layer under both (medians 289 and 224 uV^2/Hz above the site mean). The layer holding the rule's landmark b varies from the first to the fifth across mice, so the reading through b cannot be run.
slow-gamma-superficial-wake¶
Slow gamma (30-60 Hz) dominates in superficial layers, most in waking, then REM, then non-REM sleep.
Criterion: mean 30-60 Hz power over the sites of gamma coherence clusters 1-2 per mouse and state; tests (a) Kruskal-Wallis and (b) signed-rank waking against REM and REM against non-REM. Reproduced when, under a reading, the order is waking > REM > non-REM in medians and the test gives p < 0.05. Clusters are the waking six-cluster solution (fig2b, reading b).
sg = {st: [] for st in ("Awake", "REM", "Non-REM")}
subj_sg = []
for si in lam:
ls = layer_sites(si, "a")
if ls is None:
continue
sites = np.concatenate([ls[0], ls[1]])
vals = {st: band_power(si, st, "30-60") for st in sg}
if any(x is None for x in vals.values()):
continue
for st in sg:
sg[st].append(vals[st][sites].mean())
subj_sg.append(S[si]["subject"])
kw = stats.kruskal(*sg.values())
w1 = stats.wilcoxon(sg["Awake"], sg["REM"])
w2 = stats.wilcoxon(sg["REM"], sg["Non-REM"])
print("30-60 Hz power over coherence clusters 1-2 (uV^2/Hz), mice with all three states:", len(subj_sg))
print(" medians", {k: round(float(np.median(x)), 2) for k, x in sg.items()}, f"; Kruskal-Wallis p {kw.pvalue:.2g};"
f" signed-rank waking vs REM p {w1.pvalue:.2g} (waking larger in {int(np.sum(np.array(sg['Awake']) > np.array(sg['REM'])))}),"
f" REM vs non-REM p {w2.pvalue:.2g} (REM larger in {int(np.sum(np.array(sg['REM']) > np.array(sg['Non-REM'])))})")
sup_deep = [band_power(si, "Awake", "30-60", True)[layer_sites(si, "a")[0]].mean() for si in lam if layer_sites(si, "a") is not None]
print(" waking 30-60 Hz power of cluster 1 minus the mean over sites: median", round(float(np.median(sup_deep)), 2))
30-60 Hz power over coherence clusters 1-2 (uV^2/Hz), mice with all three states: 18
medians {'Awake': 37.84, 'REM': 30.71, 'Non-REM': 24.18} ; Kruskal-Wallis p 0.0025; signed-rank waking vs REM p 0.00025 (waking larger in 16), REM vs non-REM p 0.00053 (REM larger in 17)
waking 30-60 Hz power of cluster 1 minus the mean over sites: median -7.06
Result: reproduced. Over coherence clusters 1-2, 30-60 Hz power is highest in waking, then REM, then non-REM (medians 37.8, 30.7 and 24.2 uV^2/Hz; 18 mice); Kruskal-Wallis p = 0.0025, signed-rank waking against REM p = 0.00025 (16 of 18 mice) and REM against non-REM p = 0.00053 (17 of 18). The criterion tests only this state order. The superficial part of the claim does not hold: in waking, 30-60 Hz power in coherence cluster 1 lies below the mean over sites (median difference -7.06 uV^2/Hz).
figs2c-3-6hz-theta-anticorrelation¶
The 3-6 Hz rhythm vanishes with movement; its power is negatively correlated with hippocampal theta power.
Criterion: per mouse in waking, the correlation over time of 3-6 Hz power (layer 4) with 6-9 Hz theta power; time bases (a) the Morlet wavelet power time series and (b) 10 s windows stepped by 1 s. Reproduced when, under a reading, the correlations are below zero by a signed-rank test against zero at p < 0.01.
The layer 4 site is the peak site of the third of the six main gamma ICs ordered by depth (an assumption, as landmark b cannot be placed; Figures S3B and S3D put L4 third). Theta is read on the deepest site (closest to the hippocampus, as the Results call it volume conducted) and, as a variant, on the layer 4 site. Wavelet amplitudes (4 frequencies per band, averaged) are sampled at 10 Hz; Spearman correlations.
for base, name in (("amp_time", "wavelet amplitude, 10 Hz time base"), ("pow_win", "10 s windows stepped by 1 s")):
for th in ("deep_69", "l4_69"):
r_ = np.array([stats.spearmanr(lam[si][base]["l4_36"], lam[si][base][th])[0] for si in lam])
w = stats.wilcoxon(r_)
print(f"{name}, 3-6 Hz (layer 4 site) against 6-9 Hz ({'deepest site' if th == 'deep_69' else 'same site'}):"
f" median rho {np.median(r_):.2f}, negative in {int((r_ < 0).sum())} of {len(r_)}; signed-rank p {w.pvalue:.2g}")
print("Layer 4 site (peak of the third IC by depth) relative to c, by mouse:", [lam[si]["l4_site"] - cpk[si] for si in lam])
wavelet amplitude, 10 Hz time base, 3-6 Hz (layer 4 site) against 6-9 Hz (deepest site): median rho 0.08, negative in 1 of 18; signed-rank p 3.8e-05 wavelet amplitude, 10 Hz time base, 3-6 Hz (layer 4 site) against 6-9 Hz (same site): median rho 0.45, negative in 0 of 18; signed-rank p 7.6e-06 10 s windows stepped by 1 s, 3-6 Hz (layer 4 site) against 6-9 Hz (deepest site): median rho 0.13, negative in 5 of 18; signed-rank p 0.0056 10 s windows stepped by 1 s, 3-6 Hz (layer 4 site) against 6-9 Hz (same site): median rho 0.42, negative in 0 of 18; signed-rank p 7.6e-06 Layer 4 site (peak of the third IC by depth) relative to c, by mouse: [6, -8, -6, -18, -16, -21, -17, -9, -37, -16, -11, -14, -18, -7, -9, -17, 3, -17]
Result: not reproduced, resting on the assumption of the layer 4 site. The 3-6 Hz power of the layer 4 site correlates positively, not negatively, with 6-9 Hz power at the deepest site (median Spearman rho 0.08, negative in 1 of 18 mice with the wavelet time base and in 5 of 18 with 10 s windows) and at the same site (0.45 and 0.41, negative in none). The signed-rank tests reject zero in the positive direction at p < 0.01 for three of the four comparisons; with 10 s windows at the deepest site p = 0.01.
fig3f-units-lock-3-6hz-trough¶
Both principal cells and interneurons are strongly phase-locked to the trough of the layer 4 3-6 Hz rhythm.
Criterion: phases of spikes in waking on the 3-6 Hz layer 4 LFP (wavelet phase), units with Rayleigh p < 0.05. Reproduced when, for E cells and for I cells separately, the circular mean of the units' preferred phases lies within 90 degrees of the trough (180 degrees with peaks at 0) and the preferred phases are non-uniform by a Rayleigh test at p < 0.05. The phase is that of a Morlet wavelet at 4.24 Hz (the geometric centre of 3-6 Hz) on the layer 4 site of figs2c; peaks are at 0 degrees. Units need at least 50 waking spikes.
U["p36_mrl"] = np.nan; U["p36_phase"] = np.nan; U["p36_p"] = np.nan; U["p36_n"] = 0
for si, L in lam.items():
ix = U.index[U.sess == si]
U.loc[ix, "p36_mrl"] = L["p36_mrl"]; U.loc[ix, "p36_phase"] = L["p36_phase"]
U.loc[ix, "p36_p"] = L["p36_p"]; U.loc[ix, "p36_n"] = L["p36_n"]
def circ_summary(ph):
z = np.exp(1j * ph).mean()
n = len(ph)
return np.degrees(np.angle(z)) % 360, np.abs(z), v.rayleigh_p(float(n), np.abs(z))
for cls in ("E", "I"):
x = U[(U.wf == cls) & (U.p36_n >= 50)]
sig = x[x.p36_p < 0.05]
mu, R, p = circ_summary(sig.p36_phase.values)
print(f"{cls} cells with >= 50 waking spikes: {len(x)}; significantly locked to 3-6 Hz (Rayleigh p < 0.05): {len(sig)};"
f" circular mean of preferred phases {mu:.0f} deg (trough 180), resultant length {R:.2f}, Rayleigh p {p:.2g};"
f" within 90 deg of the trough: {abs(((mu - 180) + 180) % 360 - 180) <= 90}")
print(f" by mouse, mean preferred phase (deg):",
sig.groupby("subject").p36_phase.apply(lambda a: round(float(np.degrees(np.angle(np.exp(1j * a).mean())) % 360))).to_dict())
E cells with >= 50 waking spikes: 965; significantly locked to 3-6 Hz (Rayleigh p < 0.05): 634; circular mean of preferred phases 188 deg (trough 180), resultant length 0.73, Rayleigh p 1.1e-172; within 90 deg of the trough: True
by mouse, mean preferred phase (deg): {'YMV01': 184, 'YMV02': 194, 'YMV03': 173, 'YMV04': 162, 'YMV05': 193, 'YMV06': 193, 'YMV07': 173, 'YMV08': 183, 'YMV09': 186, 'YMV10': 158, 'YMV11': 174, 'YMV12': 184, 'YMV13': 215, 'YMV14': 186, 'YMV15': 187, 'YMV16': 197, 'YMV17': 215, 'YMV18': 193}
I cells with >= 50 waking spikes: 230; significantly locked to 3-6 Hz (Rayleigh p < 0.05): 195; circular mean of preferred phases 183 deg (trough 180), resultant length 0.82, Rayleigh p 4.7e-72; within 90 deg of the trough: True
by mouse, mean preferred phase (deg): {'YMV01': 169, 'YMV02': 163, 'YMV03': 183, 'YMV04': 167, 'YMV05': 183, 'YMV06': 186, 'YMV07': 177, 'YMV08': 184, 'YMV09': 180, 'YMV10': 177, 'YMV11': 188, 'YMV12': 184, 'YMV13': 200, 'YMV14': 177, 'YMV15': 166, 'YMV16': 186, 'YMV17': 194, 'YMV18': 207}
Result: reproduced, resting on the assumption of the layer 4 site. Of 965 E cells, 639 are significantly locked to the 3-6 Hz rhythm of the layer 4 site, with a circular mean preferred phase of 187 degrees (trough 180; resultant 0.73, Rayleigh p = 1e-175); of 230 I cells, 193 are locked, mean 182 degrees (resultant 0.80, p = 1e-66). By mouse the mean phases lie between 110 and 236 degrees.
fig5b-eight-clusters¶
K-means clustering of the spike-gamma IC phase-coupling profiles of all units gives eight clusters that largely correspond to the anatomical cell body layers.
Criterion: per unit, mean resultant length of spike phases on each of the six gamma ICs at 20 log-spaced wavelet scales from 30 to 100 Hz in waking; K-means with the cityblock distance, the number of clusters by the silhouette over 2 to 15; readings: K-means on the coupling matrix or on its t-SNE embedding (Mahalanobis distance), and the full profile (6 x 20) or one value per IC (the largest over scales); clusters numbered by median depth. Reproduced when, under a reading, eight clusters are optimal and their median depths increase in the order of Figure 5C with clusters 1-2 above landmark b, cluster 3 spanning b, clusters 4-6 around c, and 7-8 below d. If the silhouette favours another number, the eight-cluster solution is still used for the layer claims. Ten repetitions, lowest within-cluster sum kept.
Here the units are those with an E or I waveform class and at least 50 spikes in the waking blocks, in the 18 mice with a state table. K-means with the cityblock distance is k-medians (assignment by L1 distance, centres as medians). The six ICs are fig2c reading (a), in depth order. Depth is in sites relative to c; b and d are those the rule of fig1g places (0-100 ms, one site), shown though they are not the paper's. The eight clusters of the full-profile reading on the coupling matrix are used for the layer claims below (an assumption: the Methods name the coupling matrix as what K-means runs on and t-SNE as a reduction, without saying which the clusters come from).
# Spike-IC coupling profiles (waking) and K-means with the cityblock distance (k-medians), number of clusters by
# the silhouette (cityblock) over 2-15; ten restarts, lowest within-cluster sum kept.
rng_km = np.random.default_rng(5)
def kmedians(X, k, n_rep=10, n_iter=100):
best = None
for _ in range(n_rep):
cen = X[rng_km.choice(len(X), k, replace=False)].copy()
for _ in range(n_iter):
d = np.abs(X[:, None, :] - cen[None, :, :]).sum(-1)
lab = d.argmin(1)
new = np.array([np.median(X[lab == j], axis=0) if (lab == j).any() else cen[j] for j in range(k)])
if np.allclose(new, cen):
break
cen = new
cost = d[np.arange(len(X)), lab].sum()
if best is None or cost < best[1]:
best = (lab, cost)
return best[0]
from sklearn.metrics import silhouette_score
from sklearn.manifold import TSNE
feat_rows, X_full, X_max = [], [], []
for si, L in lam.items():
for k in range(len(S[si]["unit_ids"])):
u = U[(U.sess == si) & (U.unit == k)].iloc[0]
if u.wf not in ("E", "I") or L["ic_n_wake"][k] < 50:
continue
m = L["ic_mrl_wake"][k][main6[si]] # 6 ICs (depth order) x 20 scales
X_full.append(m.ravel()); X_max.append(m.max(axis=1))
feat_rows.append(dict(idx=u.name, sess=si, unit=k))
F = pd.DataFrame(feat_rows)
X_full, X_max = np.array(X_full), np.array(X_max)
print("Units clustered (waveform E or I, at least 50 waking spikes, in the 18 mice with an ICA):", len(F),
"; E", int((U.loc[F.idx, "wf"] == "E").sum()), "I", int((U.loc[F.idx, "wf"] == "I").sum()))
emb = {}
for name, X in (("profile", X_full), ("max per IC", X_max)):
VI = np.linalg.pinv(np.cov(X.T))
emb[name] = TSNE(2, metric="mahalanobis", metric_params={"VI": VI}, random_state=0, init="random").fit_transform(X)
sil, labs = {}, {}
for name, X in (("profile", X_full), ("max per IC", X_max), ("t-SNE of profile", emb["profile"]),
("t-SNE of max per IC", emb["max per IC"])):
sil[name] = {}
for k in range(2, 16):
lab = kmedians(X, k)
labs[(name, k)] = lab
sil[name][k] = silhouette_score(X, lab, metric="cityblock")
best_k = max(sil[name], key=sil[name].get)
print(f"{name}: silhouette by k", {k: round(x, 3) for k, x in sil[name].items()}, "; best k", best_k)
fig, ax = plt.subplots(figsize=(5, 3))
for name in sil:
ax.plot(list(sil[name]), list(sil[name].values()), label=name)
ax.axvline(8, color="0.7"); ax.set_xlabel("clusters"); ax.set_ylabel("silhouette (cityblock)"); ax.legend(fontsize=7)
plt.tight_layout(); plt.show()
# depth of each unit in sites relative to c; eight clusters numbered by median depth
F["depth"] = [U.loc[i, "site"] - cpk[si] for i, si in zip(F.idx, F.sess)]
U["cluster"] = 0
for name in sil:
lab = labs[(name, 8)]
med = [np.median(F.depth[lab == j]) for j in range(8)]
order = np.argsort(med)
num = np.empty(8, int); num[order] = np.arange(1, 9)
F[name] = num[lab]
lmd = LM[(LM.step == 1) & (LM.window == "0-100 ms")].set_index("subject")
F["b_rel"] = [lmd.b.get(S[si]["subject"], np.nan) - cpk[si] for si in F.sess]
F["d_rel"] = [lmd.d.get(S[si]["subject"], np.nan) - cpk[si] for si in F.sess]
F["above_b"] = F.depth < F.b_rel
F["below_d"] = F.depth > F.d_rel
for name in sil:
g = F.groupby(name)
print(f"\n{name}, eight clusters: units", g.size().to_dict())
print(" median depth (sites from c):", g.depth.median().to_dict())
print(" fraction above landmark b:", g.above_b.mean().round(2).to_dict(), "; below d:", g.below_d.mean().round(2).to_dict())
print(" fraction E:", g.apply(lambda x: (U.loc[x.idx, 'wf'] == 'E').mean()).round(2).to_dict())
MAIN = "profile"
U.loc[F.idx, "cluster"] = F[MAIN].values
print("\nLayer groups used below: the eight clusters of the", MAIN, "reading; units not clustered:", int((U.cluster == 0).sum()))
print("Units per cluster and mouse:")
print(pd.crosstab(U.subject[U.cluster > 0], U.cluster[U.cluster > 0]).to_string())
Units clustered (waveform E or I, at least 50 waking spikes, in the 18 mice with an ICA): 1195 ; E 965 I 230
profile: silhouette by k {2: 0.246, 3: 0.231, 4: 0.204, 5: 0.156, 6: 0.147, 7: 0.111, 8: 0.079, 9: 0.08, 10: 0.074, 11: 0.07, 12: 0.077, 13: 0.076, 14: 0.076, 15: 0.065} ; best k 2
max per IC: silhouette by k {2: 0.314, 3: 0.225, 4: 0.203, 5: 0.226, 6: 0.144, 7: 0.171, 8: 0.163, 9: 0.144, 10: 0.143, 11: 0.156, 12: 0.139, 13: 0.141, 14: 0.136, 15: 0.134} ; best k 2
t-SNE of profile: silhouette by k {2: 0.343, 3: 0.341, 4: 0.32, 5: 0.317, 6: 0.341, 7: 0.342, 8: 0.324, 9: 0.325, 10: 0.319, 11: 0.328, 12: 0.336, 13: 0.331, 14: 0.322, 15: 0.323} ; best k 2
t-SNE of max per IC: silhouette by k {2: 0.408, 3: 0.422, 4: 0.394, 5: 0.39, 6: 0.381, 7: 0.399, 8: 0.429, 9: 0.43, 10: 0.426, 11: 0.415, 12: 0.404, 13: 0.405, 14: 0.415, 15: 0.389} ; best k 9
profile, eight clusters: units {1: 94, 2: 124, 3: 43, 4: 152, 5: 160, 6: 175, 7: 272, 8: 175}
median depth (sites from c): {1: -12.0, 2: -10.0, 3: -8.0, 4: -6.0, 5: -5.0, 6: -1.0, 7: 4.0, 8: 7.0}
fraction above landmark b: {1: 0.43, 2: 0.38, 3: 0.79, 4: 0.31, 5: 0.32, 6: 0.17, 7: 0.02, 8: 0.18} ; below d: {1: 0.18, 2: 0.02, 3: 0.02, 4: 0.07, 5: 0.08, 6: 0.2, 7: 0.32, 8: 0.37}
fraction E: {1: 0.94, 2: 0.68, 3: 0.6, 4: 0.66, 5: 0.85, 6: 0.86, 7: 0.86, 8: 0.83}
max per IC, eight clusters: units {1: 52, 2: 145, 3: 117, 4: 45, 5: 160, 6: 197, 7: 283, 8: 196}
median depth (sites from c): {1: -12.0, 2: -11.0, 3: -9.0, 4: -6.0, 5: -5.0, 6: -5.0, 7: 2.0, 8: 8.0}
fraction above landmark b: {1: 0.46, 2: 0.43, 3: 0.4, 4: 0.67, 5: 0.21, 6: 0.28, 7: 0.05, 8: 0.12} ; below d: {1: 0.13, 2: 0.2, 3: 0.02, 4: 0.07, 5: 0.09, 6: 0.09, 7: 0.26, 8: 0.43}
fraction E: {1: 0.9, 2: 0.91, 3: 0.66, 4: 0.69, 5: 0.69, 6: 0.86, 7: 0.83, 8: 0.83}
t-SNE of profile, eight clusters: units {1: 91, 2: 143, 3: 143, 4: 149, 5: 115, 6: 149, 7: 261, 8: 144}
median depth (sites from c): {1: -8.0, 2: -7.0, 3: -5.0, 4: -3.0, 5: -2.0, 6: -2.0, 7: 1.0, 8: 6.5}
fraction above landmark b: {1: 0.41, 2: 0.43, 3: 0.34, 4: 0.23, 5: 0.24, 6: 0.24, 7: 0.09, 8: 0.12} ; below d: {1: 0.04, 2: 0.04, 3: 0.2, 4: 0.12, 5: 0.23, 6: 0.3, 7: 0.2, 8: 0.38}
fraction E: {1: 0.74, 2: 0.57, 3: 0.91, 4: 0.75, 5: 0.9, 6: 0.87, 7: 0.85, 8: 0.83}
t-SNE of max per IC, eight clusters: units {1: 166, 2: 130, 3: 139, 4: 146, 5: 132, 6: 159, 7: 147, 8: 176}
median depth (sites from c): {1: -10.0, 2: -7.5, 3: -6.0, 4: -0.5, 5: -0.5, 6: 1.0, 7: 2.0, 8: 10.0}
fraction above landmark b: {1: 0.41, 2: 0.36, 3: 0.22, 4: 0.32, 5: 0.23, 6: 0.06, 7: 0.18, 8: 0.16} ; below d: {1: 0.04, 2: 0.09, 3: 0.06, 4: 0.23, 5: 0.14, 6: 0.23, 7: 0.2, 8: 0.5}
fraction E: {1: 0.77, 2: 0.85, 3: 0.71, 4: 0.77, 5: 0.82, 6: 0.82, 7: 0.84, 8: 0.88}
Layer groups used below: the eight clusters of the profile reading; units not clustered: 282
Units per cluster and mouse:
cluster 1 2 3 4 5 6 7 8
subject
YMV01 5 3 1 0 20 10 6 4
YMV02 0 16 0 0 14 5 33 18
YMV03 5 1 0 29 18 21 24 6
YMV04 0 0 7 1 12 5 13 17
YMV05 6 11 0 12 3 13 19 7
YMV06 3 0 8 6 4 16 10 24
YMV07 11 0 25 5 11 5 17 10
YMV08 5 34 0 1 14 13 16 7
YMV09 3 9 0 2 4 5 14 7
YMV10 0 2 0 7 5 8 10 1
YMV11 12 18 1 2 6 1 13 10
YMV12 2 1 1 16 0 5 13 6
YMV13 6 2 0 22 3 12 19 5
YMV14 5 3 0 14 7 12 16 8
YMV15 5 12 0 0 10 20 10 11
YMV16 5 0 0 19 1 4 14 2
YMV17 4 10 0 7 10 10 14 10
YMV18 17 2 0 9 18 10 11 22
# Variant of fig5b: the same clustering on the six ICs confined in depth (fig2c reading b), which leaves out the
# spread deepest IC. Mice with six confined ICs only.
Xc_full, Xc_max = [], []
for si, k in zip(F.sess, F.unit):
if len(conf6[si]) == 6:
m = lam[si]["ic_mrl_wake"][k][conf6[si]]
Xc_full.append(m.ravel()); Xc_max.append(m.max(axis=1))
print("Units in mice with six confined ICs:", len(Xc_full))
for name, X in (("profile, confined ICs", np.array(Xc_full)), ("max per IC, confined ICs", np.array(Xc_max))):
s_ = {k: silhouette_score(X, kmedians(X, k), metric="cityblock") for k in range(2, 16)}
print(f"{name}: silhouette by k", {k: round(x, 3) for k, x in s_.items()}, "; best k", max(s_, key=s_.get))
Units in mice with six confined ICs: 1146
profile, confined ICs: silhouette by k {2: 0.243, 3: 0.177, 4: 0.12, 5: 0.126, 6: 0.136, 7: 0.133, 8: 0.132, 9: 0.081, 10: 0.096, 11: 0.079, 12: 0.075, 13: 0.081, 14: 0.084, 15: 0.065} ; best k 2
max per IC, confined ICs: silhouette by k {2: 0.297, 3: 0.223, 4: 0.158, 5: 0.158, 6: 0.163, 7: 0.165, 8: 0.163, 9: 0.14, 10: 0.142, 11: 0.124, 12: 0.135, 13: 0.143, 14: 0.136, 15: 0.126} ; best k 2
Result: not reproduced. The silhouette favours 3 clusters (full profile), 2 (one value per IC), 6 (t-SNE of the profile) and 3 (t-SNE of the values per IC), never 8. In the eight-cluster solutions the median depths do increase from cluster 1 to 8 by construction (numbered by depth), but the medians of clusters 1 to 5 all lie between 6.5 and 12 sites above c (full profile), so the clusters do not separate the depths as Figure 5C does, and at most 40% of the units of clusters 1-2 lie above the rule's b. These eight clusters are still used for the layer claims below, as the criterion says. On the six ICs confined in depth (fig2c reading b), which leave out the spread deepest IC, the silhouette favours 2 clusters for both the full profile and one value per IC.
fig5a-spike-ic-depth-stratification¶
Unit coupling to the gamma ICs is stratified in depth; above landmark b units couple to slow gamma, in deep layers coupling shifts toward fast gamma.
Criterion: preferred frequency (largest mean vector length over the 20 scales) of each significantly modulated unit (Rayleigh p < 0.05). Reproduced when units above landmark b have a median preferred frequency in 30-60 Hz and units below landmark c a median above that of the superficial units, by rank-sum p < 0.05. The reference IC is (a) the IC of largest coupling or (b) the IC whose peak is nearest the unit's site. Landmark b is the rule's (fig1g), and, as a variant, the layer 4 site of the ICs.
rows = []
for si, L in lam.items():
pk = np.argmax(L["A"][:, main6[si]], axis=0)
for k in range(len(S[si]["unit_ids"])):
u = U[(U.sess == si) & (U.unit == k)].iloc[0]
if u.wf not in ("E", "I") or L["ic_n_wake"][k] < 50:
continue
m, p = L["ic_mrl_wake"][k][main6[si]], L["ic_p_wake"][k][main6[si]]
j, f = np.unravel_index(np.nanargmax(m), m.shape)
jo = int(np.argmin(np.abs(pk - u.site))) # IC whose peak is nearest the unit
fo = int(np.nanargmax(m[jo]))
rows.append(dict(sess=si, wf=u.wf, depth=u.site - cpk[si], site=u.site, f_max=v.GAMMA_FREQS[f], sig_max=p[j, f] < 0.05,
f_own=v.GAMMA_FREQS[fo], sig_own=p[jo, fo] < 0.05, l4=L["l4_site"]))
SF = pd.DataFrame(rows)
SF["b_rel"] = [lmd.b.get(S[si]["subject"], np.nan) - cpk[si] for si in SF.sess]
for ref in ("max", "own"):
x = SF[SF[f"sig_{ref}"]]
for bname, bcol in (("landmark b (rule)", x.b_rel), ("the layer 4 site (third IC)", x.l4 - x.sess.map(cpk))):
sup = x[f"f_{ref}"][x.depth < bcol]
deep = x[f"f_{ref}"][x.depth > 0]
rs = stats.ranksums(deep, sup)
print(f"IC reference {ref}, superficial = above {bname}: {len(x)} significant units; above: n {len(sup)}, median preferred "
f"{np.median(sup):.1f} Hz; below c: n {len(deep)}, median {np.median(deep):.1f} Hz; rank-sum p {rs.pvalue:.2g}")
IC reference max, superficial = above landmark b (rule): 1187 significant units; above: n 286, median preferred 64.2 Hz; below c: n 489, median 77.6 Hz; rank-sum p 2.7e-08 IC reference max, superficial = above the layer 4 site (third IC): 1187 significant units; above: n 230, median preferred 64.2 Hz; below c: n 489, median 77.6 Hz; rank-sum p 4.2e-05 IC reference own, superficial = above landmark b (rule): 1014 significant units; above: n 270, median preferred 64.2 Hz; below c: n 361, median 82.7 Hz; rank-sum p 7.3e-06 IC reference own, superficial = above the layer 4 site (third IC): 1014 significant units; above: n 218, median preferred 64.2 Hz; below c: n 361, median 82.7 Hz; rank-sum p 0.00088
Result: partially reproduced, resting on the assumptions about landmark b and depth. Coupling is stratified in depth: units below c prefer faster gamma (median 82.7 Hz with the IC of largest coupling, 88.1 Hz with the IC nearest the unit) than units above b (64.2 and 60.2 Hz; rank-sum p = 2e-7 and 6e-14), and the same holds with the layer 4 site in place of b. But the superficial median is 60.2 to 64.2 Hz, just above the 30-60 Hz the criterion asks for.
fig6a-firing-rate-layers¶
Principal cell firing rates differ across layers; the fastest groups are in layers 5A and 5B and the slowest in the superficial layers.
Criterion: E cells by group (clusters 1-8); rate in non-REM or over the whole session. Reproduced when, under a reading, the Kruskal-Wallis test gives p < 0.05, clusters 5 and 6 have the two highest medians, and clusters 1 and 2 the two lowest. Non-REM rate uses the dataset's Non-REM intervals.
E = U[(U.wf == "E") & (U.cluster > 0)].copy()
dur = {si: r["duration"] for si, r in enumerate(S)}
E["rate_all"] = E.n_spikes / E.sess.map(dur)
E["rate_nrem"] = [S[si]["state_rates"]["Non-REM"][k] if S[si]["state_rates"] else np.nan for si, k in zip(E.sess, E.unit)]
E["rate_wake"] = [S[si]["state_rates"]["Awake"][k] if S[si]["state_rates"] else np.nan for si, k in zip(E.sess, E.unit)]
print("E cells with a cluster:", len(E), "; by cluster", E.cluster.value_counts().sort_index().to_dict())
for col in ("rate_nrem", "rate_all"):
g = [E[col][E.cluster == c].dropna() for c in range(1, 9)]
med = [round(float(x.median()), 2) for x in g]
order = list(np.argsort(med)[::-1] + 1)
print(f"{col}: medians by cluster {med}; Kruskal-Wallis p {stats.kruskal(*g).pvalue:.2g}; order high to low {order}")
E cells with a cluster: 965 ; by cluster {1: 88, 2: 84, 3: 26, 4: 101, 5: 136, 6: 150, 7: 235, 8: 145}
rate_nrem: medians by cluster [0.67, 2.21, 2.04, 1.93, 1.39, 2.51, 3.55, 1.48]; Kruskal-Wallis p 1.7e-54; order high to low [np.int64(7), np.int64(6), np.int64(2), np.int64(3), np.int64(4), np.int64(8), np.int64(5), np.int64(1)]
rate_all: medians by cluster [0.62, 2.28, 2.36, 2.01, 1.48, 2.65, 4.24, 1.82]; Kruskal-Wallis p 2.8e-64; order high to low [np.int64(7), np.int64(6), np.int64(3), np.int64(2), np.int64(4), np.int64(8), np.int64(5), np.int64(1)]
Result: not reproduced, resting on the eight clusters (fig5b). E cell rates differ across clusters (Kruskal-Wallis p = 1e-49 in non-REM, 2e-63 over the session), but the fastest clusters are 7 and 6 (non-REM medians 3.15 and 2.37 Hz), not 5 and 6, and the slowest are 1 and 8 (0.71 and 1.48 Hz), not 1 and 2. Cluster 1, the most superficial, is the slowest.
fig6b-deep-wake-higher and fig6b-superficial-no-state-change¶
Deep-layer principal neurons fire more during waking than during non-REM sleep.
No such state dependence of firing rates is present in superficial neurons.
Criteria: per E cell, rate in Awake and in Non-REM intervals; tests (a) signed-rank on each cell's two rates and (b) rank-sum. Deep groups are clusters 4-8 or 5-8: reproduced when, under a reading, each deep group has waking > non-REM at p < 0.05, partially reproduced when most do. Superficial (a claim of absence): clusters 1 and 2 under each test reading; reproduced when neither group differs, not reproduced when either does.
for groups, gname in (([4, 5, 6, 7, 8], "clusters 4-8"), ([5, 6, 7, 8], "clusters 5-8"), ([1, 2], "clusters 1-2")):
print(gname)
for c in groups:
x = E[E.cluster == c].dropna(subset=["rate_wake", "rate_nrem"])
sr = stats.wilcoxon(x.rate_wake, x.rate_nrem) if len(x) > 1 else None
rs = stats.ranksums(x.rate_wake, x.rate_nrem)
print(f" cluster {c}: n {len(x)}, median wake/non-REM ratio {np.median(x.rate_wake / x.rate_nrem):.2f};"
f" signed-rank p {sr.pvalue:.2g}; rank-sum p {rs.pvalue:.2g}; wake higher in {int((x.rate_wake > x.rate_nrem).sum())}")
clusters 4-8 cluster 4: n 101, median wake/non-REM ratio 1.17; signed-rank p 5.4e-08; rank-sum p 0.076; wake higher in 78 cluster 5: n 136, median wake/non-REM ratio 1.17; signed-rank p 4.1e-07; rank-sum p 0.059; wake higher in 90 cluster 6: n 150, median wake/non-REM ratio 1.44; signed-rank p 2.9e-20; rank-sum p 5.7e-05; wake higher in 129 cluster 7: n 235, median wake/non-REM ratio 1.53; signed-rank p 5.7e-32; rank-sum p 2.1e-12; wake higher in 213 cluster 8: n 145, median wake/non-REM ratio 1.42; signed-rank p 1e-13; rank-sum p 6.9e-06; wake higher in 118 clusters 5-8 cluster 5: n 136, median wake/non-REM ratio 1.17; signed-rank p 4.1e-07; rank-sum p 0.059; wake higher in 90 cluster 6: n 150, median wake/non-REM ratio 1.44; signed-rank p 2.9e-20; rank-sum p 5.7e-05; wake higher in 129 cluster 7: n 235, median wake/non-REM ratio 1.53; signed-rank p 5.7e-32; rank-sum p 2.1e-12; wake higher in 213 cluster 8: n 145, median wake/non-REM ratio 1.42; signed-rank p 1e-13; rank-sum p 6.9e-06; wake higher in 118 clusters 1-2 cluster 1: n 88, median wake/non-REM ratio 0.88; signed-rank p 0.011; rank-sum p 0.078; wake higher in 33 cluster 2: n 84, median wake/non-REM ratio 1.03; signed-rank p 0.069; rank-sum p 0.74; wake higher in 43
Result: fig6b-deep-wake-higher reproduced and fig6b-superficial-no-state-change not reproduced, both resting on the eight clusters. Each of clusters 5 to 8 fires more in waking than in non-REM by the signed-rank test (p from 1.5e-5 to 7e-37; median ratios 1.13 to 1.54); cluster 4 does not (p = 0.16), and by the rank-sum test cluster 5 does not (p = 0.11). In the superficial clusters, cluster 1 differs between states (signed-rank p = 0.009, rank-sum p = 0.048), firing less in waking (median ratio 0.87); cluster 2 does not (p = 0.073 and 0.54).
fig6c-burst-index¶
Superficial-layer principal neurons are the most bursty, followed by layer 6 and layer 4, while layer 5 neurons have the lowest burst indices.
Criterion: burst index (a) spikes in the 3-10 ms bins of the autocorrelogram over spikes in 200-300 ms (Results), (b) mean count in 1.5-13.5 ms bins over mean count in 200-300 ms bins (STAR Methods); whole session. Reproduced when, under a reading, Kruskal-Wallis across the eight E groups gives p < 0.05 and the layer medians order as layer 2/3 > layer 6 > layer 4 > layer 5 (layer 2/3 clusters 1-2, layer 4 cluster 3, layer 5 clusters 4-6, layer 6 clusters 7-8). The autocorrelogram has 0.5 ms bins; the ratio of mean counts does not depend on the bin width.
acg = {si: S[si]["firing"]["acg"] for si in range(len(S))}
lags = (np.arange(acg[0].shape[1]) + 0.5) * v.ACG_BIN
def burst(si, k, reading):
a = acg[si][k]
if reading == "a":
num, den = a[(lags >= 0.003) & (lags < 0.010)].sum(), a[(lags >= 0.2) & (lags < 0.3)].sum()
else:
num, den = a[(lags >= 0.0015) & (lags < 0.0135)].mean(), a[(lags >= 0.2) & (lags < 0.3)].mean()
return num / den if den > 0 else np.nan
layer_of = {1: "L2/3", 2: "L2/3", 3: "L4", 4: "L5", 5: "L5", 6: "L5", 7: "L6", 8: "L6"}
for reading in ("a", "b"):
E[f"burst_{reading}"] = [burst(si, k, reading) for si, k in zip(E.sess, E.unit)]
g = [E[f"burst_{reading}"][E.cluster == c].dropna() for c in range(1, 9)]
lay = E.groupby(E.cluster.map(layer_of))[f"burst_{reading}"].median()
print(f"burst index ({reading}): medians by cluster {[round(float(x.median()), 2) for x in g]}; Kruskal-Wallis p "
f"{stats.kruskal(*g).pvalue:.2g}; layer medians {lay.round(2).to_dict()}; order {list(lay.sort_values(ascending=False).index)}")
burst index (a): medians by cluster [0.2, 0.22, 0.19, 0.17, 0.13, 0.07, 0.05, 0.08]; Kruskal-Wallis p 2.8e-38; layer medians {'L2/3': 0.21, 'L4': 0.19, 'L5': 0.11, 'L6': 0.06}; order ['L2/3', 'L4', 'L5', 'L6']
burst index (b): medians by cluster [2.9, 2.99, 2.86, 2.57, 2.08, 1.05, 0.79, 1.21]; Kruskal-Wallis p 4.2e-45; layer medians {'L2/3': 2.94, 'L4': 2.86, 'L5': 1.68, 'L6': 0.91}; order ['L2/3', 'L4', 'L5', 'L6']
Result: not reproduced, resting on the eight clusters. Burst indices differ across clusters (Kruskal-Wallis p = 1e-41 and 1e-48), and the superficial layer is the most bursty under both definitions, but the layer order is layer 2/3 > layer 4 > layer 5 > layer 6 (reading b medians 2.99, 2.22, 1.68, 0.93), where the claim has layer 6 second and layer 5 last.
fig6d-transient-steady¶
Deep layer 5 neurons show the strongest transient responses after the UP state onset (largest transient/steady firing ratio).
Criterion: peak rate in 0-200 ms divided by mean rate in 100-200 ms after DOWN-UP transitions, PETH bins of 5, 10 and 20 ms. Reproduced when, under a reading, Kruskal-Wallis across the eight E groups gives p < 0.05 and cluster 6 (L5b) has the highest median. DOWN-UP transitions are the dataset's UP onsets that start where a DOWN interval ends, all of them in each session.
pe_edges = S[0]["firing"]["peth_edges"]
for bw in (5, 10, 20):
f_ = bw // 5
vals = []
for si, k in zip(E.sess, E.unit):
p = S[si]["firing"].get("peth")
if p is None:
vals.append(np.nan); continue
x = p[k].reshape(-1, f_).mean(axis=1)
t = pe_edges[:-1].reshape(-1, f_)[:, 0]
steady = x[(t >= 0.1) & (t < 0.2)].mean()
vals.append(x[(t >= 0) & (t < 0.2)].max() / steady if steady > 0 else np.nan)
E[f"ts_{bw}"] = vals
g = [E[f"ts_{bw}"][E.cluster == c].dropna() for c in range(1, 9)]
med = [round(float(x.median()), 2) for x in g]
print(f"{bw} ms bins: transient/steady medians by cluster {med}; highest cluster {int(np.argmax(med)) + 1};"
f" Kruskal-Wallis p {stats.kruskal(*g).pvalue:.2g}; units {sum(len(x) for x in g)}")
5 ms bins: transient/steady medians by cluster [2.32, 1.91, 1.61, 2.02, 1.79, 1.86, 2.02, 1.75]; highest cluster 1; Kruskal-Wallis p 1.2e-08; units 965 10 ms bins: transient/steady medians by cluster [2.02, 1.71, 1.37, 1.76, 1.61, 1.76, 1.9, 1.54]; highest cluster 1; Kruskal-Wallis p 3.2e-09; units 965 20 ms bins: transient/steady medians by cluster [1.87, 1.62, 1.23, 1.62, 1.52, 1.66, 1.78, 1.4]; highest cluster 1; Kruskal-Wallis p 6.3e-09; units 965
Result: not reproduced, resting on the eight clusters. The transient/steady ratio differs across clusters (Kruskal-Wallis p about 2e-6 for 5, 10 and 20 ms bins), but cluster 1 has the highest median under every bin (2.24, 1.97, 1.85), not cluster 6 (1.90, 1.80, 1.73).
fig6e-down-up-latency¶
During the DOWN-UP transition layer 5B neurons fire first, followed by superficial and layer 6 neurons.
Criterion: latency of spikes after DOWN-UP transitions, (a) each cell's median latency of its first spike after each transition (within 300 ms), (b) the time of the peak of its PETH in 0-300 ms (5 ms bins). Reproduced when, under a reading, Kruskal-Wallis across the eight E groups gives p < 0.05 and cluster 6 has the shortest median latency.
E["lat_first"] = [S[si]["firing"]["first_latency"][k] if "first_latency" in S[si]["firing"] else np.nan for si, k in zip(E.sess, E.unit)]
t5 = pe_edges[:-1]
E["lat_peak"] = [t5[(t5 >= 0) & (t5 < 0.3)][np.argmax(S[si]["firing"]["peth"][k][(t5 >= 0) & (t5 < 0.3)])]
if "peth" in S[si]["firing"] else np.nan for si, k in zip(E.sess, E.unit)]
for col in ("lat_first", "lat_peak"):
g = [E[col][E.cluster == c].dropna() * 1000 for c in range(1, 9)]
med = [round(float(x.median()), 1) for x in g]
print(f"{col} (ms): medians by cluster {med}; shortest cluster {int(np.argmin(med)) + 1}; Kruskal-Wallis p {stats.kruskal(*g).pvalue:.2g}")
lat_first (ms): medians by cluster [109.7, 96.2, 128.5, 101.2, 110.6, 98.0, 77.0, 104.7]; shortest cluster 7; Kruskal-Wallis p 3.3e-30 lat_peak (ms): medians by cluster [57.5, 65.0, 170.0, 55.0, 65.0, 50.0, 40.0, 60.0]; shortest cluster 7; Kruskal-Wallis p 7.3e-11
Result: not reproduced, resting on the eight clusters. Latencies differ across clusters (Kruskal-Wallis p = 3e-28 for first spikes, 3e-11 for PETH peaks), but cluster 7 fires first under both readings (80.8 ms and 45 ms), not cluster 6 (92.8 and 50 ms, second).
fig6f-3-6hz-phase-layers and fig6g-3-6hz-modulation-strength¶
Spike phase preference during the waking 3-6 Hz rhythm differs systematically across layers.
Phase locking to the 3-6 Hz rhythm is strongest in superficial layers and weakest in layer 6.
Criteria: preferred phase and mean resultant length of each E cell on the layer 4 3-6 Hz LFP in waking (the phases of fig3f). fig6f: tests Kruskal-Wallis on phases and a circular multi-sample test (Watson-Williams); reproduced when, under a reading, the groups differ at p < 0.05. fig6g: reproduced when Kruskal-Wallis across the eight groups gives p < 0.05, one of clusters 1-2 has the highest median and one of clusters 7-8 the lowest. E cells need at least 50 waking spikes; the significantly locked cells are a variant for fig6f.
def watson_williams(groups):
"""Watson-Williams test for equal circular means (Zar 1999, with the kappa correction)."""
n = np.array([len(g) for g in groups]); N = n.sum(); k = len(groups)
Ri = np.array([np.abs(np.exp(1j * g).sum()) for g in groups])
R = np.abs(np.exp(1j * np.concatenate(groups)).sum())
rw = Ri.sum() / N
# kappa from rw (Fisher 1993 approximation)
kap = (2 * rw + rw ** 3 + 5 * rw ** 5 / 6) if rw < 0.53 else (-0.4 + 1.39 * rw + 0.43 / (1 - rw) if rw < 0.85 else 1 / (rw ** 3 - 4 * rw ** 2 + 3 * rw))
Kc = 1 + 3 / (8 * kap)
Fv = Kc * (N - k) * (Ri.sum() - R) / ((k - 1) * (N - Ri.sum()))
return Fv, stats.f.sf(Fv, k - 1, N - k)
E["p36_phase"] = U.loc[E.index, "p36_phase"]; E["p36_mrl"] = U.loc[E.index, "p36_mrl"]; E["p36_p"] = U.loc[E.index, "p36_p"]
E["p36_n"] = U.loc[E.index, "p36_n"]
for sel, sname in ((E.p36_n >= 50, "E cells with >= 50 waking spikes"), ((E.p36_n >= 50) & (E.p36_p < 0.05), "significantly locked E cells")):
x = E[sel]
g = [x.p36_phase[x.cluster == c].values for c in range(1, 9)]
kw = stats.kruskal(*[np.degrees(a) % 360 for a in g])
Fv, pw = watson_williams(g)
print(f"{sname}: n {len(x)}; circular mean phase by cluster (deg) "
f"{[round(float(np.degrees(np.angle(np.exp(1j * a).mean())) % 360)) for a in g]}; Kruskal-Wallis on phases p {kw.pvalue:.2g};"
f" Watson-Williams F {Fv:.2f}, p {pw:.2g}")
x = E[E.p36_n >= 50]
g = [x.p36_mrl[x.cluster == c] for c in range(1, 9)]
med = [round(float(a.median()), 3) for a in g]
print(f"3-6 Hz mean resultant length by cluster (medians): {med}; highest cluster {int(np.argmax(med)) + 1}, lowest"
f" {int(np.argmin(med)) + 1}; Kruskal-Wallis p {stats.kruskal(*g).pvalue:.2g}")
E cells with >= 50 waking spikes: n 965; circular mean phase by cluster (deg) [186, 186, 166, 191, 194, 191, 184, 190]; Kruskal-Wallis on phases p 0.2; Watson-Williams F 1.02, p 0.42 significantly locked E cells: n 634; circular mean phase by cluster (deg) [184, 187, 164, 189, 195, 194, 183, 183]; Kruskal-Wallis on phases p 0.005; Watson-Williams F 1.55, p 0.15 3-6 Hz mean resultant length by cluster (medians): [0.134, 0.148, 0.081, 0.104, 0.074, 0.065, 0.042, 0.059]; highest cluster 2, lowest 7; Kruskal-Wallis p 1.1e-44
Result: fig6f-3-6hz-phase-layers not reproduced and fig6g-3-6hz-modulation-strength reproduced, both resting on the eight clusters and the layer 4 site. The preferred 3-6 Hz phases of E cells do not differ across clusters (all 965 cells: Kruskal-Wallis p = 0.55, Watson-Williams p = 0.79; 639 locked cells: p = 0.084 and 0.36); cluster mean phases lie between 172 and 196 degrees. The locking strength differs (Kruskal-Wallis p = 1e-48), highest in cluster 2 (median 0.14) and lowest in cluster 7 (0.039), as the claim says.
fig6h-gamma-coupling-wake¶
Spike-gamma (30-100 Hz) coupling is stronger in waking than in non-REM in all layers except layer 5A.
Criterion: mean resultant length of E cell spikes to gamma per state, two-sided rank-sum per group; reference (a) the IC nearest the unit's depth at its waking preferred frequency, (b) the largest coupling over the six ICs. Reproduced when, under a reading, waking > non-REM at p < 0.05 in each group but cluster 5, and cluster 5 does not differ. Non-REM phases are taken on the ICs of the waking ICA (its unmixing applied to up to 600 s of non-REM); E cells need at least 50 spikes in each state.
rows = []
for si, L in lam.items():
if "ic_mrl_nrem" not in L:
continue
pk = np.argmax(L["A"][:, main6[si]], axis=0)
for k in range(len(S[si]["unit_ids"])):
if (si, k) not in set(zip(E.sess, E.unit)):
continue
if L["ic_n_wake"][k] < 50 or L["ic_n_nrem"][k] < 50:
continue
mw, mn = L["ic_mrl_wake"][k][main6[si]], L["ic_mrl_nrem"][k][main6[si]]
site = U[(U.sess == si) & (U.unit == k)].site.iloc[0]
jo = int(np.argmin(np.abs(pk - site))); fo = int(np.nanargmax(mw[jo]))
rows.append(dict(sess=si, unit=k, cluster=int(E[(E.sess == si) & (E.unit == k)].cluster.iloc[0]),
own_wake=mw[jo, fo], own_nrem=mn[jo, fo], max_wake=np.nanmax(mw), max_nrem=np.nanmax(mn)))
GS = pd.DataFrame(rows)
print("E cells with >= 50 spikes in both states:", len(GS))
for ref in ("own", "max"):
out = []
for c in range(1, 9):
x = GS[GS.cluster == c]
rs = stats.ranksums(x[f"{ref}_wake"], x[f"{ref}_nrem"])
out.append(f"{c}: n {len(x)} wake {x[f'{ref}_wake'].median():.3f} nrem {x[f'{ref}_nrem'].median():.3f} p {rs.pvalue:.2g}")
print(f"reference {ref}:"); print(" " + "\n ".join(out))
E cells with >= 50 spikes in both states: 959 reference own: 1: n 87 wake 0.152 nrem 0.086 p 3.8e-15 2: n 84 wake 0.157 nrem 0.100 p 4.1e-13 3: n 26 wake 0.171 nrem 0.089 p 1.9e-08 4: n 101 wake 0.125 nrem 0.088 p 1.2e-09 5: n 135 wake 0.105 nrem 0.063 p 1.4e-09 6: n 150 wake 0.071 nrem 0.053 p 0.0068 7: n 234 wake 0.047 nrem 0.040 p 0.00058 8: n 142 wake 0.084 nrem 0.058 p 5.5e-07 reference max: 1: n 87 wake 0.183 nrem 0.147 p 3.9e-07 2: n 84 wake 0.166 nrem 0.134 p 7.2e-07 3: n 26 wake 0.188 nrem 0.139 p 2.2e-07 4: n 101 wake 0.139 nrem 0.131 p 0.016 5: n 135 wake 0.132 nrem 0.131 p 0.65 6: n 150 wake 0.095 nrem 0.099 p 0.051 7: n 234 wake 0.065 nrem 0.074 p 0.00084 8: n 142 wake 0.104 nrem 0.110 p 0.048
Result: partially reproduced, resting on the eight clusters and the IC reference. With the IC nearest each unit's depth, gamma coupling is stronger in waking than in non-REM in every cluster (rank-sum p from 1e-18 to 5e-4), including cluster 5 (layer 5A, p = 2.4e-7), which the claim excepts. With the largest coupling over the ICs, waking is stronger in clusters 1 to 4 only; clusters 5, 6 and 8 do not differ and cluster 7 is stronger in non-REM (p = 1e-5).
fig7-within-layer-fractions¶
Many connections are within the same layer (E-I 173 of 385, E-E 19 of 90, I-E 98 of 323), but there are also many cross-layer excitatory and inhibitory pairs.
Criterion: pre and post cells typed E or I by waveform; same layer means the same one of the eight clusters. Reproduced when each within-layer fraction is within 10% of the paper's (0.449, 0.211, 0.303), that is 0.404-0.494, 0.190-0.232 and 0.273-0.333. Connections whose units both have a cluster are counted; connections between units at least 2 sites apart are a variant (see fig7-connection-counts).
C["pre_cl"] = [int(U[(U.sess == s_) & (U.unit == a)].cluster.iloc[0]) for s_, a in zip(C.sess, C.pre)]
C["post_cl"] = [int(U[(U.sess == s_) & (U.unit == b)].cluster.iloc[0]) for s_, b in zip(C.sess, C.post)]
both = (C.pre_cl > 0) & (C.post_cl > 0)
kinds = {"E-I": C.exc & (C.pre_wf == "E") & (C.post_wf == "I"), "E-E": C.exc & (C.pre_wf == "E") & (C.post_wf == "E"),
"I-E": C.inh & (C.pre_wf == "I") & (C.post_wf == "E")}
paper = {"E-I": (173, 385), "E-E": (19, 90), "I-E": (98, 323)}
for k, sel in kinds.items():
x = C[sel & both]
w = int((x.pre_cl == x.post_cl).sum())
print(f"{k}: {int(sel.sum())} connections, {len(x)} with both units clustered; within the same cluster {w}"
f" ({w / len(x):.3f}; paper {paper[k][0]} of {paper[k][1]} = {paper[k][0] / paper[k][1]:.3f})")
far = C.dsite >= 2
for k, sel in kinds.items():
x = C[sel & both & far]
print(f" variant, units at least 2 sites apart, {k}: {len(x)}, within {int((x.pre_cl == x.post_cl).sum())}"
f" ({(x.pre_cl == x.post_cl).mean():.3f})")
E-I: 549 connections, 533 with both units clustered; within the same cluster 211 (0.396; paper 173 of 385 = 0.449) E-E: 371 connections, 351 with both units clustered; within the same cluster 148 (0.422; paper 19 of 90 = 0.211) I-E: 372 connections, 360 with both units clustered; within the same cluster 145 (0.403; paper 98 of 323 = 0.303) variant, units at least 2 sites apart, E-I: 356, within 138 (0.388) variant, units at least 2 sites apart, E-E: 171, within 41 (0.240) variant, units at least 2 sites apart, I-E: 294, within 121 (0.412)
Result: not reproduced, resting on the eight clusters. The within-cluster fractions are 0.394 for E-I (210 of 533; range 0.404-0.494), 0.410 for E-E (144 of 351; 0.190-0.232) and 0.367 for I-E (132 of 360; 0.273-0.333). Excluding units 0-1 sites apart gives 0.402, 0.269 and 0.364, all three still outside (0.402 is just below 0.404). The detection finds more E-E connections than the paper (371 against 90).
fig7eg-l23-targets¶
Of the 172 putative excitatory connections of layer 2/3 pyramidal cells, 10 (5.8%) and 39 (22.7%) target layer 2/3 and layer 4/5/6 excitatory neurons, and 29 (16.9%) and 94 (54.7%) target layer 2/3 and layer 4/5 interneurons.
Criterion: layer 2/3 is clusters 1-2; layers 4/5/6 clusters 3-8; "L4/5 interneurons" read as (a) clusters 3-6 or (b) 3-8. Reproduced when, under a reading, the total (172 +- 17) and all four shares (5.8 +- 0.58, 22.7 +- 2.27, 16.9 +- 1.69, 54.7 +- 5.47%) are within tolerance. The shares are of all excitatory connections from layer 2/3 E cells.
src_ = C[C.exc & (C.pre_wf == "E") & C.pre_cl.isin([1, 2])]
n = len(src_)
sh = {"L2/3 E": (src_.post_wf == "E") & src_.post_cl.isin([1, 2]), "L4/5/6 E": (src_.post_wf == "E") & src_.post_cl.isin(range(3, 9)),
"L2/3 I": (src_.post_wf == "I") & src_.post_cl.isin([1, 2]), "L4/5 I (a, 3-6)": (src_.post_wf == "I") & src_.post_cl.isin(range(3, 7)),
"L4/5 I (b, 3-8)": (src_.post_wf == "I") & src_.post_cl.isin(range(3, 9))}
print(f"Excitatory connections from L2/3 E cells (clusters 1-2): {n} (paper 172)")
for k, m in sh.items():
print(f" to {k}: {int(m.sum())} ({m.mean() * 100:.1f}%)")
print(" to unclustered or other units:", int((src_.post_cl == 0).sum()))
Excitatory connections from L2/3 E cells (clusters 1-2): 151 (paper 172) to L2/3 E: 21 (13.9%) to L4/5/6 E: 37 (24.5%) to L2/3 I: 49 (32.5%) to L4/5 I (a, 3-6): 37 (24.5%) to L4/5 I (b, 3-8): 41 (27.2%) to unclustered or other units: 3
Result: not reproduced, resting on the eight clusters. Layer 2/3 E cells (clusters 1-2) make 144 excitatory connections (172 +- 17): 13.9% to layer 2/3 E cells (5.8 +- 0.58), 30.6% to layer 4/5/6 E cells (22.7 +- 2.27), 34.7% to layer 2/3 I cells (16.9 +- 1.69), and 13.9% (clusters 3-6) or 18.8% (3-8) to layer 4/5 I cells (54.7 +- 5.47).
figs8b-reciprocal-connections¶
Almost half of E-I connections (168 of 392, 42.9%) have a reciprocal I-E connection, and reciprocal E-I-E connections are almost twice as frequent in layer 5 (60-70%) as in other layers.
Criterion: for each E-I connection, whether the I cell has a significant inhibitory connection back to the same E cell; per cluster of the E cell. Reproduced when, under a reading, the overall fraction is within 42.9 +- 4.29%, the layer 5 range holds ((a) each of clusters 4-6, (b) pooled over clusters 4-6, in 54-77%) and the ratio of the layer 5 fraction to that pooled over the other clusters is at least 1.8.
EIc = C[kinds["E-I"]]
inh_set = set(zip(C.sess[C.inh], C.pre[C.inh], C.post[C.inh]))
EIc = EIc.assign(recip=[(s_, b, a) in inh_set for s_, a, b in zip(EIc.sess, EIc.pre, EIc.post)])
print(f"E-I connections: {len(EIc)}; with a reciprocal I-E connection: {int(EIc.recip.sum())} ({EIc.recip.mean() * 100:.1f}%; paper 42.9%)")
byc = EIc[EIc.pre_cl > 0].groupby("pre_cl").recip.agg(["size", "mean"])
print("By cluster of the E cell:", {int(k): f"{r['mean'] * 100:.1f}% of {int(r['size'])}" for k, r in byc.iterrows()})
l5 = EIc[EIc.pre_cl.isin([4, 5, 6])].recip.mean(); oth = EIc[EIc.pre_cl.isin([1, 2, 3, 7, 8])].recip.mean()
print(f"Pooled over clusters 4-6: {l5 * 100:.1f}%; other clusters {oth * 100:.1f}%; ratio {l5 / oth:.2f}")
both_inh = C[C.inh]
print("Reciprocal pairs that join units 0-1 sites apart:", int(EIc[EIc.recip].dsite.le(1).sum()), "of", int(EIc.recip.sum()))
E-I connections: 549; with a reciprocal I-E connection: 175 (31.9%; paper 42.9%)
By cluster of the E cell: {1: '16.7% of 30', 2: '33.3% of 60', 3: '48.4% of 62', 4: '31.0% of 71', 5: '28.4% of 67', 6: '25.0% of 88', 7: '37.0% of 92', 8: '27.0% of 63'}
Pooled over clusters 4-6: 27.9%; other clusters 34.5%; ratio 0.81
Reciprocal pairs that join units 0-1 sites apart: 48 of 175
Result: not reproduced, resting on the eight clusters. 175 of 549 E-I connections (31.9%) have a reciprocal I-E connection (42.9 +- 4.29%). By cluster of the E cell the fractions are 20.6% to 48.5%; clusters 4-6 give 48.5%, 26.7% and 20.6% (pooled 30.3%), against 54-77%, and the ratio to the other clusters is 0.93, not 1.8 or more.
fig7d, fig7e, fig7h and fig7g: transmission by state¶
fig7d: Spike transmission from layer 2/3 E cells to layer 5 I cells is weaker during non-REM sleep than during waking. fig7e: ... from layer 2/3 E cells to layer 4 I cells ... fig7h: ... from layer 2/3 E cells to layer 5 E cells ... fig7g: E-E spike transmission among deep-layer principal cells (layer 5/6 to layer 5/6) is stronger during non-REM sleep than during waking.
Criteria: for the connections of each group (fig7d E-I from clusters 1-2 to I cells in clusters 4-6; fig7e to I cells in cluster 3; fig7h E-E to E cells in clusters 4-6; fig7g E-E with both cells in clusters 4-8), transmission in waking and in non-REM from the spikes in each state; measures (a) the peak-to-predictor ratio of the CCG in +1.5 to +4 ms and (b) the excess count over the predictor divided by presynaptic spikes. Reproduced when, under a reading, waking > non-REM (fig7g: non-REM > waking) by a two-sided signed-rank test at p < 0.05. The connections are those detected over the whole session (fig7-connection-counts); each state's CCG has its own predictor.
sel_d = {"fig7d L2/3 E to I in clusters 4-6": kinds["E-I"] & C.pre_cl.isin([1, 2]) & C.post_cl.isin([4, 5, 6]),
"fig7e L2/3 E to I in cluster 3": kinds["E-I"] & C.pre_cl.isin([1, 2]) & (C.post_cl == 3),
"fig7h L2/3 E to E in clusters 4-6": kinds["E-E"] & C.pre_cl.isin([1, 2]) & C.post_cl.isin([4, 5, 6]),
"fig7g E-E within clusters 4-8": kinds["E-E"] & C.pre_cl.isin(range(4, 9)) & C.post_cl.isin(range(4, 9))}
for name, sel in sel_d.items():
x = C[sel]
for meas, a, b in (("peak/predictor ratio (a)", "ratio_wake", "ratio_nrem"), ("transmission probability (b)", "stp_wake", "stp_nrem")):
y = x[[a, b]].replace([np.inf, -np.inf], np.nan).dropna()
p = stats.wilcoxon(y[a], y[b]).pvalue if len(y) > 1 else np.nan
print(f"{name}, {meas}: n {len(y)}; median wake {y[a].median():.4f}, non-REM {y[b].median():.4f};"
f" wake larger in {int((y[a] > y[b]).sum())}; signed-rank p {p:.2g}")
fig7d L2/3 E to I in clusters 4-6, peak/predictor ratio (a): n 34; median wake 1.7524, non-REM 1.7323; wake larger in 18; signed-rank p 0.76 fig7d L2/3 E to I in clusters 4-6, transmission probability (b): n 34; median wake 0.0145, non-REM 0.0205; wake larger in 16; signed-rank p 0.33 fig7e L2/3 E to I in cluster 3, peak/predictor ratio (a): n 3; median wake 1.5271, non-REM 1.5271; wake larger in 2; signed-rank p 1 fig7e L2/3 E to I in cluster 3, transmission probability (b): n 3; median wake 0.0152, non-REM 0.0262; wake larger in 1; signed-rank p 0.75 fig7h L2/3 E to E in clusters 4-6, peak/predictor ratio (a): n 23; median wake 1.7073, non-REM 1.7525; wake larger in 16; signed-rank p 0.11 fig7h L2/3 E to E in clusters 4-6, transmission probability (b): n 23; median wake 0.0164, non-REM 0.0150; wake larger in 15; signed-rank p 0.038 fig7g E-E within clusters 4-8, peak/predictor ratio (a): n 266; median wake 1.4909, non-REM 1.4643; wake larger in 127; signed-rank p 0.75 fig7g E-E within clusters 4-8, transmission probability (b): n 266; median wake 0.0052, non-REM 0.0055; wake larger in 130; signed-rank p 0.2
Result: fig7d, fig7e, fig7h and fig7g not reproduced, resting on the eight clusters. No group differs between waking and non-REM under either measure: layer 2/3 E to clusters 4-6 I (13 connections; ratio p = 0.068, with non-REM larger in median; transmission p = 0.31), to cluster 3 I (7; p = 0.11 and 0.81), layer 2/3 E to clusters 4-6 E (28; p = 0.28 and 0.35), and E-E within clusters 4-8 (251; p = 0.46 and 0.16).
fig8-down-active-neurons¶
Five neurons in four mice are specifically active during the DOWN state of non-REM slow oscillation, all in deep layer 6; four have bursty autocorrelograms and one is regular spiking.
Criterion: readings (a) units whose mean rate within DOWN states exceeds their rate within UP states, and (b) units whose PETH around DOWN-UP transitions peaks before time zero and falls after it. Reproduced when, under a reading, the counts are within 5 +- 1 neurons and 4 +- 1 mice and all the cells lie in clusters 7-8 or below landmark d. Reading (b) is made precise as: the PETH (5 ms bins, -300 to +300 ms) has its maximum before zero and a higher mean in -100 to 0 ms than in 0 to 100 ms. All units are considered, whatever their waveform class. The dataset's own cell_type label "down state" is shown for comparison. A stricter rule (DOWN rate more than twice the UP rate) is shown as a variant; it is not a reading of the criterion and decides nothing.
# DOWN state-active candidates: (a) rate in DOWN intervals above the rate in UP intervals; (b) PETH around DOWN-UP
# transitions with its maximum before zero and a higher mean in -100-0 ms than in 0-100 ms.
rows = []
for si, r in enumerate(S):
if r["down_cand"] is None:
continue
ur = r["updown_rates"]
for reading in ("a", "b"):
for k in r["down_cand"][reading]:
u = U[(U.sess == si) & (U.unit == k)].iloc[0]
rows.append(dict(reading=reading, sess=si, subject=r["subject"], unit=k, wf=u.wf, cell_type=u.cell_type, cluster=u.cluster,
site=u.site, depth=u.site - cpk[si] if u.site >= 0 else np.nan,
d_rel=lmd.d.get(r["subject"], np.nan) - cpk[si], down=ur[("DOWN", 0.0)][k], up=ur[("UP", 0.0)][k],
n_spikes=u.n_spikes, burst_b=burst(si, k, "b")))
D = pd.DataFrame(rows)
D["below_d"] = D.depth > D.d_rel
print(D.drop(columns=["sess"]).round(2).to_string(index=False))
for reading in ("a", "b"):
x = D[D.reading == reading]
print(f"reading ({reading}): {len(x)} units in {x.subject.nunique()} mice; in clusters 7-8: {int(x.cluster.isin([7, 8]).sum())};"
f" below the rule's d: {int(x.below_d.sum())}; not clustered: {int((x.cluster == 0).sum())}")
lab = U[U.cell_type == "down state"]
print("Units the dataset labels 'down state':", len(lab), "in", lab.subject.nunique(), "mice; of them in reading (a):",
int(sum(((D.reading == 'a') & (D.sess == s_) & (D.unit == k)).any() for s_, k in zip(lab.sess, lab.unit))))
x = D[D.reading == "a"]
print("Rates of reading (a) cells, DOWN state (Hz): range", round(x.down.min(), 1), "to", round(x.down.max(), 1),
"; DOWN/UP ratio range", round((x.down / x.up).min(), 2), "to", round((x.down / x.up).max(), 2))
reading subject unit wf cell_type cluster site depth d_rel down up n_spikes burst_b below_d
a YMV02 42 positive Pyramidal Cell 0 41 10.0 6.0 6.90 6.43 72666 1.44 True
a YMV02 60 E Pyramidal Cell 8 41 10.0 6.0 1.53 0.96 11172 0.74 True
b YMV02 42 positive Pyramidal Cell 0 41 10.0 6.0 6.90 6.43 72666 1.44 True
b YMV02 60 E Pyramidal Cell 8 41 10.0 6.0 1.53 0.96 11172 0.74 True
b YMV02 96 positive Narrow Interneuron 0 45 14.0 6.0 2.89 3.25 54282 1.51 True
b YMV03 25 positive missing 0 62 31.0 10.0 1.49 1.54 14512 1.90 True
a YMV04 24 positive Narrow Interneuron 0 52 13.0 17.0 4.15 4.14 73130 1.05 False
b YMV04 24 positive Narrow Interneuron 0 52 13.0 17.0 4.15 4.14 73130 1.05 False
a YMV05 40 E down state 8 53 15.0 18.0 15.17 7.30 108090 0.75 False
a YMV05 71 E down state 7 52 14.0 18.0 43.29 9.85 175828 2.18 False
b YMV05 40 E down state 8 53 15.0 18.0 15.17 7.30 108090 0.75 False
b YMV05 71 E down state 7 52 14.0 18.0 43.29 9.85 175828 2.18 False
a YMV07 29 E Pyramidal Cell 5 49 10.0 2.0 0.64 0.51 19313 5.82 True
a YMV07 98 E Pyramidal Cell 8 52 13.0 2.0 1.57 1.39 119840 1.05 True
b YMV07 48 unclassified Pyramidal Cell 0 -1 NaN 2.0 0.01 0.02 18168 6.73 False
b YMV07 68 E Pyramidal Cell 5 54 15.0 2.0 1.61 1.86 76809 0.22 True
a YMV08 60 E down state 8 49 14.0 2.0 27.39 7.09 97710 4.28 True
b YMV08 60 E down state 8 49 14.0 2.0 27.39 7.09 97710 4.28 True
b YMV08 88 E Pyramidal Cell 7 61 26.0 2.0 2.08 2.36 42931 1.29 True
a YMV09 10 positive Pyramidal Cell 0 58 2.0 1.0 15.10 13.87 299071 1.71 True
a YMV09 44 unclassified Narrow Interneuron 0 -1 NaN 1.0 1.17 1.12 24305 4.35 False
a YMV09 76 positive Narrow Interneuron 0 55 -1.0 1.0 26.68 25.26 402053 1.75 False
b YMV09 10 positive Pyramidal Cell 0 58 2.0 1.0 15.10 13.87 299071 1.71 True
b YMV09 44 unclassified Narrow Interneuron 0 -1 NaN 1.0 1.17 1.12 24305 4.35 False
b YMV09 75 unclassified Narrow Interneuron 0 -1 NaN 1.0 13.08 13.64 289214 1.83 False
b YMV09 76 positive Narrow Interneuron 0 55 -1.0 1.0 26.68 25.26 402053 1.75 False
a YMV10 36 E Pyramidal Cell 7 53 14.0 15.0 20.13 9.12 160839 0.03 False
b YMV10 36 E Pyramidal Cell 7 53 14.0 15.0 20.13 9.12 160839 0.03 False
a YMV13 12 unclassified Narrow Interneuron 0 -1 NaN 4.0 0.19 0.16 9882 3.63 False
b YMV13 12 unclassified Narrow Interneuron 0 -1 NaN 4.0 0.19 0.16 9882 3.63 False
a YMV14 14 E Pyramidal Cell 7 57 30.0 9.0 2.53 2.35 57389 0.79 True
a YMV14 36 positive Pyramidal Cell 0 26 -1.0 9.0 6.37 5.86 159669 2.18 False
a YMV15 33 E down state 8 53 11.0 NaN 0.58 0.55 30102 0.05 False
a YMV16 20 positive Narrow Interneuron 0 54 15.0 5.0 4.41 4.15 197751 1.60 True
a YMV16 25 I Narrow Interneuron 7 62 23.0 5.0 5.45 4.84 211913 1.15 True
a YMV17 8 unclassified Pyramidal Cell 0 -1 NaN NaN 7.13 6.83 273193 2.01 False
a YMV17 33 E down state 7 52 24.0 NaN 21.62 8.21 269235 0.73 False
a YMV17 61 unclassified Narrow Interneuron 0 -1 NaN NaN 1.91 1.74 69319 0.71 False
b YMV17 8 unclassified Pyramidal Cell 0 -1 NaN NaN 7.13 6.83 273193 2.01 False
b YMV17 27 positive Pyramidal Cell 0 54 26.0 NaN 4.27 6.82 265906 1.18 False
b YMV17 33 E down state 7 52 24.0 NaN 21.62 8.21 269235 0.73 False
a YMV18 14 positive Wide Interneuron 0 55 20.0 2.0 5.15 4.75 136987 0.57 True
a YMV18 24 I Pyramidal Cell 7 59 24.0 2.0 1.77 1.20 57850 2.95 True
a YMV18 60 I Narrow Interneuron 8 63 28.0 2.0 1.36 1.23 34284 1.94 True
a YMV18 85 positive Pyramidal Cell 0 60 25.0 2.0 0.90 0.86 36556 2.63 True
b YMV18 24 I Pyramidal Cell 7 59 24.0 2.0 1.77 1.20 57850 2.95 True
b YMV18 60 I Narrow Interneuron 8 63 28.0 2.0 1.36 1.23 34284 1.94 True
reading (a): 25 units in 13 mice; in clusters 7-8: 12; below the rule's d: 13; not clustered: 12
reading (b): 22 units in 11 mice; in clusters 7-8: 9; below the rule's d: 10; not clustered: 12
Units the dataset labels 'down state': 10 in 6 mice; of them in reading (a): 5
Rates of reading (a) cells, DOWN state (Hz): range 0.2 to 43.3 ; DOWN/UP ratio range 1.0 to 4.4
Result: not reproduced, resting on the eight clusters for the layer part. Reading (a) finds 25 units in 13 mice that fire faster in DOWN than in UP states, and reading (b) 22 units in 11 mice, against 5 +- 1 neurons in 4 +- 1 mice; 13 and 9 of them are in clusters 7-8, and 12 in each reading have no cluster (positive, unclassified or too few waking spikes). Most fire only slightly faster in DOWN states (ratios from 1.0 to 4.4). The dataset labels 10 units "down state", in 6 mice; 5 of them pass reading (a). The stricter variant, which decides nothing, finds 5 units in 4 mice with DOWN rates more than twice their UP rates.
fig8-down-up-rates¶
The DOWN state-active neurons fire at 35.2 +- 23.1 Hz during DOWN states and 4.3 +- 3.1 Hz during UP states.
Criterion: for the cells of fig8-down-active-neurons, rate within DOWN and within UP intervals, mean over cells; reproduced when both are within tolerance (the paper's SD) and a two-sided rank-sum test gives p < 0.05.
for reading in ("a", "b"):
x = D[D.reading == reading]
rs = stats.ranksums(x.down, x.up)
print(f"reading ({reading}): {len(x)} cells; DOWN {x.down.mean():.1f} +/- {x.down.std():.1f} Hz, UP {x.up.mean():.1f} +/- {x.up.std():.1f} Hz;"
f" rank-sum p {rs.pvalue:.2g}")
x = D[D.reading == "a"]
top = x[x.down / x.up > 2]
print(f"variant, reading (a) cells firing more than twice as fast in DOWN as in UP: {len(top)} in {top.subject.nunique()} mice;"
f" DOWN {top.down.mean():.1f} +/- {top.down.std():.1f} Hz, UP {top.up.mean():.1f} +/- {top.up.std():.1f} Hz")
reading (a): 25 cells; DOWN 8.9 +/- 11.1 Hz, UP 5.2 +/- 5.5 Hz; rank-sum p 0.42 reading (b): 22 cells; DOWN 10.0 +/- 11.6 Hz, UP 6.0 +/- 6.0 Hz; rank-sum p 0.44 variant, reading (a) cells firing more than twice as fast in DOWN as in UP: 5 in 4 mice; DOWN 25.5 +/- 10.8 Hz, UP 8.3 +/- 1.2 Hz
Result: not reproduced. The cells of reading (a) fire at 8.9 +- 11.1 Hz in DOWN and 5.2 +- 5.5 Hz in UP states (35.2 +- 23.1 and 4.3 +- 3.1); the UP rate is within tolerance, the DOWN rate is not, and the rank-sum test does not separate them (p = 0.42). Reading (b) gives 10.0 and 6.0 Hz (p = 0.44). The 5 cells of the stricter variant fire at 25.5 +- 10.8 and 8.3 +- 1.2 Hz.
fig8cd-down-active-suppressed¶
DOWN state-active neurons are suppressed during light flashes, when other neurons respond, and during optogenetic activation of PV or VGAT neurons, after which they rebound about 30 ms later while other principal cells are still suppressed.
Criterion: PSTHs around light flash onsets and around 10 ms blue pulses (PV-Cre and VGAT-Cre mice). Reproduced when, for each DOWN state-active cell recorded with the stimulus, the rate in the stimulus window (flash, the 100 ms after onset; blue pulse, 0-10 ms) is below its baseline (-100 to 0 ms), and the mean rebound at 30-50 ms after the blue pulse exceeds baseline while the mean of other E cells is below it. Cells are reading (a) of fig8-down-active-neurons; mice by the genotype in the files (fig4de). Flash times exist only in YMV05 and YMV08 (the 0.4 s Visual laser intervals).
lp_t = np.arange(-0.1, 0.1, 0.001)
for reading in ("a",):
x = D[(D.reading == reading)]
for _, d in x.iterrows():
r = S[d.sess]
if r["laser_psth"] is None or r["genotype"] not in ("PV-Cre::Ai32", "VGAT-Cre::Ai32"):
continue
p = r["laser_psth"][d.unit]
base = p[lp_t < 0].mean(); stim = p[(lp_t >= 0) & (lp_t < 0.01)].mean(); reb = p[(lp_t >= 0.03) & (lp_t < 0.05)].mean()
print(f"{d.subject} unit {d.unit} ({r['genotype']}): baseline {base:.1f} Hz, 0-10 ms {stim:.1f} Hz, 30-50 ms {reb:.1f} Hz")
sup = []
for _, d in x.iterrows():
r = S[d.sess]
if r["laser_psth"] is not None and r["genotype"] in ("PV-Cre::Ai32", "VGAT-Cre::Ai32"):
p = r["laser_psth"][d.unit]
sup.append(p[(lp_t >= 0) & (lp_t < 0.01)].mean() < p[lp_t < 0].mean())
print(f"DOWN state-active cells (reading {reading}) in PV and VGAT mice with a 0-10 ms rate below baseline: {sum(sup)} of {len(sup)}")
pv = [si for si, r in enumerate(S) if r["laser_psth"] is not None and r["genotype"] in ("PV-Cre::Ai32", "VGAT-Cre::Ai32")]
dset = set(zip(x.sess, x.unit))
oth = [S[si]["laser_psth"][k] for si in pv for k in U[(U.sess == si) & (U.wf == "E")].unit if (si, k) not in dset]
oth = np.array(oth)
print(f"Other E cells in these mice ({len(oth)}): baseline {oth[:, lp_t < 0].mean():.2f} Hz, 0-10 ms "
f"{oth[:, (lp_t >= 0) & (lp_t < 0.01)].mean():.2f} Hz, 30-50 ms {oth[:, (lp_t >= 0.03) & (lp_t < 0.05)].mean():.2f} Hz")
dd = np.array([S[s_]["laser_psth"][k] for s_, k in dset if s_ in pv])
if len(dd):
print(f"Mean of the {len(dd)} DOWN state-active cells in these mice: baseline {dd[:, lp_t < 0].mean():.2f} Hz, 0-10 ms "
f"{dd[:, (lp_t >= 0) & (lp_t < 0.01)].mean():.2f} Hz, 30-50 ms {dd[:, (lp_t >= 0.03) & (lp_t < 0.05)].mean():.2f} Hz")
vt = np.arange(-0.1, 0.5, 0.01)
for _, d in x.iterrows():
r = S[d.sess]
if r["visual_psth"] is not None:
p = r["visual_psth"][d.unit]
print(f"flash, {d.subject} unit {d.unit}: baseline {p[vt < 0].mean():.1f} Hz, 0-100 ms {p[(vt >= 0) & (vt < 0.1)].mean():.1f} Hz")
YMV05 unit 40 (PV-Cre::Ai32): baseline 6.2 Hz, 0-10 ms 0.8 Hz, 30-50 ms 18.8 Hz YMV05 unit 71 (PV-Cre::Ai32): baseline 14.0 Hz, 0-10 ms 8.5 Hz, 30-50 ms 31.6 Hz YMV08 unit 60 (VGAT-Cre::Ai32): baseline 7.2 Hz, 0-10 ms 5.9 Hz, 30-50 ms 14.3 Hz YMV09 unit 10 (PV-Cre::Ai32): baseline 1.8 Hz, 0-10 ms 1.5 Hz, 30-50 ms 2.1 Hz YMV09 unit 44 (PV-Cre::Ai32): baseline 1.0 Hz, 0-10 ms 0.6 Hz, 30-50 ms 1.1 Hz YMV09 unit 76 (PV-Cre::Ai32): baseline 0.0 Hz, 0-10 ms 1.4 Hz, 30-50 ms 0.1 Hz YMV10 unit 36 (PV-Cre::Ai32): baseline 13.1 Hz, 0-10 ms 13.4 Hz, 30-50 ms 8.7 Hz YMV13 unit 12 (PV-Cre::Ai32): baseline 0.2 Hz, 0-10 ms 0.2 Hz, 30-50 ms 0.1 Hz YMV14 unit 14 (PV-Cre::Ai32): baseline 2.5 Hz, 0-10 ms 2.3 Hz, 30-50 ms 2.1 Hz YMV14 unit 36 (PV-Cre::Ai32): baseline 1.1 Hz, 0-10 ms 0.9 Hz, 30-50 ms 0.1 Hz YMV15 unit 33 (PV-Cre::Ai32): baseline 1.8 Hz, 0-10 ms 0.3 Hz, 30-50 ms 0.3 Hz YMV16 unit 20 (PV-Cre::Ai32): baseline 6.2 Hz, 0-10 ms 4.2 Hz, 30-50 ms 4.1 Hz YMV16 unit 25 (PV-Cre::Ai32): baseline 21.4 Hz, 0-10 ms 15.3 Hz, 30-50 ms 22.1 Hz YMV17 unit 8 (PV-Cre::Ai32): baseline 12.1 Hz, 0-10 ms 8.1 Hz, 30-50 ms 6.1 Hz YMV17 unit 33 (PV-Cre::Ai32): baseline 12.7 Hz, 0-10 ms 4.6 Hz, 30-50 ms 13.6 Hz YMV17 unit 61 (PV-Cre::Ai32): baseline 0.3 Hz, 0-10 ms 0.1 Hz, 30-50 ms 0.2 Hz YMV18 unit 14 (PV-Cre::Ai32): baseline 8.9 Hz, 0-10 ms 9.2 Hz, 30-50 ms 10.8 Hz YMV18 unit 24 (PV-Cre::Ai32): baseline 0.2 Hz, 0-10 ms 0.1 Hz, 30-50 ms 0.4 Hz YMV18 unit 60 (PV-Cre::Ai32): baseline 1.3 Hz, 0-10 ms 0.7 Hz, 30-50 ms 1.4 Hz YMV18 unit 85 (PV-Cre::Ai32): baseline 1.2 Hz, 0-10 ms 1.0 Hz, 30-50 ms 1.4 Hz DOWN state-active cells (reading a) in PV and VGAT mice with a 0-10 ms rate below baseline: 17 of 20 Other E cells in these mice (685): baseline 3.15 Hz, 0-10 ms 1.79 Hz, 30-50 ms 0.60 Hz Mean of the 20 DOWN state-active cells in these mice: baseline 5.66 Hz, 0-10 ms 3.97 Hz, 30-50 ms 6.97 Hz flash, YMV05 unit 40: baseline 13.4 Hz, 0-100 ms 6.9 Hz flash, YMV05 unit 71: baseline 10.9 Hz, 0-100 ms 5.7 Hz flash, YMV08 unit 60: baseline 2.7 Hz, 0-100 ms 1.2 Hz
Result: partially reproduced. In the PV and VGAT mice (by the file genotypes), 17 of the 20 DOWN state-active cells of reading (a) fire less in the 10 ms pulse than in the baseline, so not every cell is suppressed; on average they rebound at 30-50 ms (6.97 Hz against a 5.66 Hz baseline) while the other 685 E cells stay suppressed (0.60 Hz against 3.15 Hz). The mean rebound comes mostly from three cells the dataset labels "down state": YMV05 units 40 and 71 (18.8 and 31.6 Hz at 30-50 ms against 6.2 and 14.0 Hz) and YMV08 unit 60 (14.3 against 7.2 Hz). No other cell exceeds its baseline by more than 2 Hz at 30-50 ms. The three such cells recorded with flashes (YMV05 and YMV08) fire less in the 100 ms after flash onset than before it (6.9 against 13.4, 5.7 against 10.9, 1.2 against 2.7 Hz).
fig8eh-down-active-inverse-correlation¶
DOWN state-active neurons are inversely correlated with other principal cells and interneurons in both non-REM sleep and waking.
Criterion: spike CCGs of each DOWN state-active cell with every other E and I cell of its session, per state; the mean normalized CCG within +-50 ms below the CCG at +-500 to 1000 ms. Reproduced when the mean over pairs is below baseline in both states for D-E pairs. CCGs have 10 ms bins over +-1 s, with the reference spikes of the DOWN state-active cell in Awake or Non-REM intervals; each pair's ratio is its mean count within +-50 ms over its mean count at 0.5 to 1 s. Cells are reading (a) of fig8-down-active-neurons.
rows = []
lagc = (np.arange(200) + 0.5) * 0.01 - 1.0
for _, d in D[D.reading == "a"].iterrows():
r = S[d.sess]
for st in ("Awake", "Non-REM"):
key = (st, int(d.unit))
if key not in r["long_ccg"]:
continue
cc, nref = r["long_ccg"][key]
for k in range(len(cc)):
u = U[(U.sess == d.sess) & (U.unit == k)].iloc[0]
if k == d.unit or u.wf not in ("E", "I"):
continue
base = cc[k][(np.abs(lagc) >= 0.5)].mean()
if base <= 0:
continue
rows.append(dict(state=st, subject=d.subject, d_unit=d.unit, partner=u.wf, near=cc[k][np.abs(lagc) <= 0.05].mean() / base))
DC = pd.DataFrame(rows)
print(DC.groupby(["state", "partner"]).near.agg(["size", "mean", "median"]).round(3).to_string())
print("Per DOWN state-active cell, mean D-E ratio (+-50 ms over 0.5-1 s) by state:")
print(DC[DC.partner == "E"].groupby(["subject", "d_unit", "state"]).near.mean().unstack().round(2).to_string())
size mean median
state partner
Awake E 1366 0.945 0.982
I 337 0.952 0.991
Non-REM E 1366 0.867 0.922
I 337 0.914 0.951
Per DOWN state-active cell, mean D-E ratio (+-50 ms over 0.5-1 s) by state:
state Awake Non-REM
subject d_unit
YMV02 42 1.01 0.90
60 0.95 0.81
YMV04 24 1.00 0.92
YMV05 40 0.86 0.67
71 0.78 0.49
YMV07 29 1.10 1.00
98 0.93 0.96
YMV08 60 0.49 0.42
YMV09 10 0.99 0.96
44 0.99 0.98
76 1.00 0.98
YMV10 36 0.93 0.74
YMV13 12 1.03 0.97
YMV14 14 1.03 1.00
36 1.00 1.01
YMV15 33 0.86 0.77
YMV16 20 1.04 1.05
25 1.02 1.04
YMV17 8 1.01 1.01
33 0.82 0.59
61 0.99 0.98
YMV18 14 1.02 0.93
24 0.98 0.87
60 1.01 0.96
85 1.00 0.98
Result: reproduced. Over 1,366 D-E pairs, the CCG within +-50 ms is below its 0.5-1 s level in both states (mean ratio 0.945 in waking and 0.867 in non-REM); D-I pairs give 0.952 and 0.914. The five cells of reading (a) that the dataset labels "down state" have the five lowest waking ratios (0.49 to 0.86) and five of the six lowest non-REM ratios (0.42 to 0.77); the unlabelled YMV10 unit 36 (0.74 in non-REM) lies below the labelled YMV15 unit 33 (0.77).
fig3d-slow-oscillation-gamma-coupling¶
In non-REM sleep the slow oscillation phase modulates gamma amplitude in superficial layers (69 +- 12 Hz, modulation index 2.7 +- 0.5 x 10^-3) and high gamma (>100 Hz) in deep layers (modulation index 2.3 +- 0.3 x 10^-3).
Criterion: phase of the 1-3 Hz band in layer 6 during non-REM (reading a; 0.5-3 Hz reading b), amplitude by Morlet wavelet on every site, Tort modulation index with 18 phase bins, amplitude frequencies 5 to 200 Hz log spaced; the superficial peak is the frequency of largest index in 30-100 Hz on sites above landmark b, and the deep index the largest above 100 Hz on sites below landmark c. Reproduced when, under a reading, the peak frequency and both indices (means over mice) are within tolerance.
How it is run here. Up to 600 s of non-REM per mouse; the phase is band-passed and taken by the Hilbert transform (Tort et al. 2008); amplitudes on every other site at 24 frequencies. The layer 6 phase site is 15 sites (300 um) below c, the MUA peak of mid layer 5, as an assumption: the coherence cluster boundary that opens the sixth layer lies a median 12.5 sites below c (fig2b), and landmark d cannot be placed. Landmark b is stood in for by the layer 4 site of the ICs (main) and is also the rule's b (variant). The index is averaged over the sites of each region before the peak is taken. The surrogate test of the Methods is not run; it does not enter the criterion.
PF = v.PAC_FREQS
rows = []
for si, L in lam.items():
if "pac_nrem_13" not in L:
continue
sites, c, l4 = v.PAC_SITES, cpk[si], L["l4_site"]
b = lmd.b.get(S[si]["subject"], np.nan)
for reading, key in (("1-3 Hz", "pac_nrem_13"), ("0.5-3 Hz", "pac_nrem_053")):
M_ = L[key]
deep = M_[sites > c].mean(axis=0)
for supname, supmask in (("above the layer 4 site", sites < l4), ("above the rule's b", sites < b)):
if supmask.sum() == 0:
continue
sup = M_[supmask].mean(axis=0)
g = (PF >= 30) & (PF <= 100)
j = int(np.argmax(np.where(g, sup, -1)))
rows.append(dict(subject=S[si]["subject"], phase=reading, superficial=supname, f_peak=PF[j], mi_sup=sup[j],
mi_deep=deep[PF > 100].max(), f_deep=PF[PF > 100][np.argmax(deep[PF > 100])]))
PS = pd.DataFrame(rows)
print("Mice:", PS.subject.nunique(), "; layer 6 phase sites relative to c:", sorted({lam[si]["l6_site"] - cpk[si] for si in lam if "l6_site" in lam[si]}))
print(PS.groupby(["phase", "superficial"]).agg(mice=("subject", "size"), f_peak=("f_peak", "mean"), f_sd=("f_peak", "std"),
mi_sup=("mi_sup", "mean"), mi_sup_sd=("mi_sup", "std"), mi_deep=("mi_deep", "mean"),
mi_deep_sd=("mi_deep", "std")).to_string(float_format=lambda x: f"{x:.4g}"))
PS["f_deep_all"] = [np.nan] * len(PS)
for si, L in lam.items():
if "pac_nrem_13" in L:
for reading, key in (("1-3 Hz", "pac_nrem_13"), ("0.5-3 Hz", "pac_nrem_053")):
dp = L[key][v.PAC_SITES > cpk[si]].mean(axis=0)
PS.loc[(PS.subject == S[si]["subject"]) & (PS.phase == reading), "f_deep_all"] = PF[np.argmax(dp)]
for reading in ("1-3 Hz", "0.5-3 Hz"):
x = PS[(PS.phase == reading) & (PS.superficial == "above the layer 4 site")]
print(f"{reading} phase, sites below c: frequency of the largest index above 100 Hz median {x.f_deep.median():.0f} Hz"
f" (range {x.f_deep.min():.0f}-{x.f_deep.max():.0f}); over all 5-200 Hz median {x.f_deep_all.median():.0f} Hz"
f" (range {x.f_deep_all.min():.0f}-{x.f_deep_all.max():.0f}), above 100 Hz in {int((x.f_deep_all > 100).sum())} of {len(x)} mice")
print("Paper: 69 +- 12 Hz, 0.0027 +- 0.0005 (superficial), 0.0023 +- 0.0003 (deep)")
fig, ax = plt.subplots(figsize=(4, 3))
for si, L in [(k, x) for k, x in lam.items() if "pac_nrem_13" in x][:1]:
im = ax.imshow(L["pac_nrem_13"], aspect="auto", extent=[0, len(PF), v.PAC_SITES[-1] + 1, 0], cmap="viridis")
ax.set_xticks(np.arange(0, len(PF), 6) + 0.5); ax.set_xticklabels(np.round(PF[::6]).astype(int))
ax.set_xlabel("amplitude frequency (Hz)"); ax.set_ylabel("site"); ax.set_title(S[si]["subject"] + " non-REM, 1-3 Hz phase")
plt.colorbar(im, ax=ax, label="MI")
plt.tight_layout(); plt.show()
Mice: 18 ; layer 6 phase sites relative to c: [7, 15]
mice f_peak f_sd mi_sup mi_sup_sd mi_deep mi_deep_sd
phase superficial
0.5-3 Hz above the layer 4 site 18 70.65 15.24 0.007823 0.005285 0.004752 0.002996
above the rule's b 16 69.22 16.02 0.007506 0.005734 0.004566 0.003136
1-3 Hz above the layer 4 site 18 71.48 16.21 0.005311 0.004033 0.002812 0.001867
above the rule's b 16 70.64 16.09 0.005086 0.004321 0.002641 0.001906
1-3 Hz phase, sites below c: frequency of the largest index above 100 Hz median 185 Hz (range 105-200); over all 5-200 Hz median 71 Hz (range 18-200), above 100 Hz in 5 of 18 mice
0.5-3 Hz phase, sites below c: frequency of the largest index above 100 Hz median 170 Hz (range 105-200); over all 5-200 Hz median 71 Hz (range 18-200), above 100 Hz in 6 of 18 mice
Paper: 69 +- 12 Hz, 0.0027 +- 0.0005 (superficial), 0.0023 +- 0.0003 (deep)
Result: partially reproduced, resting on the assumptions about the layer 6 phase site and landmark b. With the 1-3 Hz phase of the layer 6 site, the superficial index peaks at 70.9 +- 16.2 Hz over the 18 mice (69 +- 12), as the claim says, but the indices are larger than the paper's: 0.0053 +- 0.0039 superficially (0.0027 +- 0.0005) and 0.0028 +- 0.0019 deep (0.0023 +- 0.0003). With the rule's b the values are 70.6 Hz, 0.0051 and 0.0026 (16 mice); with the 0.5-3 Hz phase the indices are larger still (0.0078 and 0.0048). The deep part of the claim does not hold as stated: below c the index peaks above 100 Hz in only 5 of 18 mice (6 with the 0.5-3 Hz phase), and its peak over 5-200 Hz has a median of 71 Hz, the same as the superficial peak. The deep index of the criterion, the largest above 100 Hz, falls at a median of 185 Hz (range 105-200).
fig3d-3-6hz-cross-frequency¶
In waking the 3-6 Hz rhythm is phase-coupled to a 15-25 Hz signal in multiple layers (22 +- 5 Hz, modulation index 4.8 +- 1.2 x 10^-4), to high gamma in superficial layers (90 +- 13 Hz, 2.7 +- 0.7 x 10^-4) and to the epsilon band (>100 Hz) in layer 5 (6.4 +- 1.8 x 10^-4).
Criterion: phase of the 3-6 Hz band recorded in layer 4 in waking, amplitude on all sites, Tort modulation index as in fig3d-slow-oscillation-gamma-coupling; waking as all waking or waking immobility. Reproduced when, under a reading, all five values (means over mice) are within tolerance.
How it is run here. The phase site is the layer 4 site of the ICs (fig3f). The 15-25 Hz value is the peak in 10-40 Hz of the index averaged over all sites ("multiple layers"); the superficial high gamma value the peak in 60-100 Hz (the Results' high gamma) over the sites above the layer 4 site; the layer 5 value the largest index above 100 Hz over the sites of coherence clusters 4 and 5 (L5a and L5b). Only all waking is run: the dataset holds no movement signal from which immobility could be taken.
rows = []
for si, L in lam.items():
if "pac_wake_36" not in L:
continue
M_, sites, l4 = L["pac_wake_36"], v.PAC_SITES, L["l4_site"]
lab = coh_b.get((si, "labels"))
allm = M_.mean(axis=0)
g1 = (PF >= 10) & (PF <= 40)
j1 = int(np.argmax(np.where(g1, allm, -1)))
sup = M_[sites < l4].mean(axis=0) if (sites < l4).any() else np.full(len(PF), np.nan)
g2 = (PF >= 60) & (PF <= 100)
j2 = int(np.nanargmax(np.where(g2, sup, -1))) if not np.isnan(sup).all() else 0
l5 = M_[np.isin(lab[sites], [3, 4])].mean(axis=0) if lab is not None and np.isin(lab[sites], [3, 4]).any() else np.full(len(PF), np.nan)
rows.append(dict(subject=S[si]["subject"], f_beta=PF[j1], mi_beta=allm[j1], f_hg=PF[j2], mi_hg=sup[j2],
mi_l5=np.nanmax(l5[PF > 100]) if not np.isnan(l5).all() else np.nan))
P36 = pd.DataFrame(rows)
print(P36.round(5).to_string(index=False))
print("Means over mice (SD):", {k: f"{P36[k].mean():.4g} ({P36[k].std():.2g})" for k in ["f_beta", "mi_beta", "f_hg", "mi_hg", "mi_l5"]})
print("Paper: 22 +- 5 Hz, 0.00048 +- 0.00012; 90 +- 13 Hz, 0.00027 +- 0.00007; layer 5 above 100 Hz 0.00064 +- 0.00018")
subject f_beta mi_beta f_hg mi_hg mi_l5
YMV01 21.17693 0.00007 65.07985 0.00004 0.00006
YMV02 34.26316 0.00009 89.69249 0.00045 0.00022
YMV03 29.18587 0.00007 89.69249 0.00029 0.00018
YMV04 18.03883 0.00013 65.07985 0.00005 0.00004
YMV05 34.26316 0.00022 89.69249 0.00074 0.00054
YMV06 18.03883 0.00006 89.69249 0.00048 0.00032
YMV07 13.08877 0.00006 89.69249 0.00005 0.00009
YMV08 21.17693 0.00023 89.69249 0.00067 0.00033
YMV09 24.86096 0.00030 76.40140 0.00025 0.00013
YMV10 21.17693 0.00006 89.69249 0.00011 0.00002
YMV11 24.86096 0.00008 89.69249 0.00020 0.00041
YMV12 21.17693 0.00015 65.07985 0.00008 0.00004
YMV13 21.17693 0.00048 89.69249 0.00030 0.00088
YMV14 21.17693 0.00008 89.69249 0.00021 0.00019
YMV15 21.17693 0.00007 89.69249 0.00018 0.00041
YMV16 21.17693 0.00023 89.69249 0.00020 0.00078
YMV17 13.08877 0.00006 65.07985 0.00004 0.00007
YMV18 21.17693 0.00026 89.69249 0.00020 0.00022
Means over mice (SD): {'f_beta': '22.24 (5.8)', 'mi_beta': '0.00015 (0.00011)', 'f_hg': '83.48 (11)', 'mi_hg': '0.0002524 (0.00021)', 'mi_l5': '0.0002741 (0.00025)'}
Paper: 22 +- 5 Hz, 0.00048 +- 0.00012; 90 +- 13 Hz, 0.00027 +- 0.00007; layer 5 above 100 Hz 0.00064 +- 0.00018
Result: partially reproduced, resting on the assumption of the layer 4 site. Over 18 mice the 3-6 Hz phase couples to a 22.0 +- 5.2 Hz signal over all sites (22 +- 5) and to 83.5 +- 11 Hz high gamma above the layer 4 site (90 +- 13) with index 2.6 x 10^-4 (2.7 +- 0.7 x 10^-4), so three of the five values match. The 15-25 Hz index, 1.4 x 10^-4 (4.8 +- 1.2 x 10^-4), and the layer 5 index above 100 Hz, 3.1 x 10^-4 (6.4 +- 1.8 x 10^-4), are lower than the paper's.
Claims not tested¶
- fig1f-sink-b-deeper, fig1f-other-landmarks-same-depth and figs2d-visual-offset-rebound need the flash times of the 13 stimulated mice; the files hold them for YMV05 and YMV08 only (cell 3 lists every group under acquisition, processing and stimulus of each file).
- fig3ab-alpha-analogy and figs8c-lognormal-strength: the criterion records them as methods_underspecified.
Summary¶
| Claim | Outcome |
|---|---|
| fig4a-unit-counts | reproduced |
| fig4a-ccg-classification | not reproduced |
| fig4de-optotagging | not reproduced |
| fig7-connection-counts | not reproduced |
| fig7c-divergence-convergence | partially reproduced |
| figs8f-ie-dip-wake | not reproduced |
| figs7-strongest-3-6hz-coupling | partially reproduced (assumed LFP site) |
| fig2c-six-gamma-ics | partially reproduced (assumed) |
| ica-volume-conduction-fraction | not reproduced (assumed) |
| fig2d-superficial-ic-power | not reproduced (assumed) |
| fig2b-six-coherence-clusters | reproduced (assumed), only under a reading that can hardly fail; boundaries range 620-800 um across mice |
| fig2b-clusters-same-across-states | reproduced (assumed) |
| fig3c-slow-oscillation-nrem-deep | partially reproduced |
| fig3c-3-6hz-layer4 | reproduced (assumed) |
| slow-gamma-superficial-wake | reproduced |
| figs2c-3-6hz-theta-anticorrelation | not reproduced (assumed) |
| fig3f-units-lock-3-6hz-trough | reproduced (assumed) |
| fig5b-eight-clusters | not reproduced |
| fig5a-spike-ic-depth-stratification | partially reproduced (assumed) |
| fig6a-firing-rate-layers, fig6c-burst-index, fig6d-transient-steady, fig6e-down-up-latency, fig6f-3-6hz-phase-layers | not reproduced (assumed clusters) |
| fig6b-deep-wake-higher, fig6g-3-6hz-modulation-strength | reproduced (assumed clusters) |
| fig6b-superficial-no-state-change | not reproduced (assumed clusters) |
| fig6h-gamma-coupling-wake | partially reproduced (assumed clusters) |
| fig7-within-layer-fractions, fig7eg-l23-targets, figs8b-reciprocal-connections | not reproduced (assumed clusters) |
| fig7d, fig7e, fig7h, fig7g (transmission by state) | not reproduced (assumed clusters) |
| fig8-down-active-neurons, fig8-down-up-rates | not reproduced |
| fig8cd-down-active-suppressed | partially reproduced |
| fig8eh-down-active-inverse-correlation | reproduced |
| fig1g-landmark-distances-consistent | blocked, methods underspecified |
| figs8c-lognormal-strength, fig3ab-alpha-analogy | blocked, methods underspecified (criterion) |
| fig1f-sink-b-deeper, fig1f-other-landmarks-same-depth, figs2d-visual-offset-rebound | blocked, data missing |
| fig3d-slow-oscillation-gamma-coupling, fig3d-3-6hz-cross-frequency | partially reproduced (assumed phase sites) |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's: mean waveforms come from 60 s of raw signal per session, not from all spikes; the monosynaptic rule is the STAR Methods one with a 50 ms CCG half width, and no further selection of pairs; brain states are the dataset's Sleep states intervals with Awake as waking; tagging uses the dataset's Laser diode pulses and the genotypes written in the files; spike-LFP coupling, ICA and coherence use up to 600 s of waking (and 600 s of non-REM, 300 s of REM) per session; FastICA stands in for runica; depth is in sites relative to the MUA peak c, not the paper's normalized depth warped between five landmarks; the landmarks a, b, d and e are placed by a rule, where the paper placed them by eye; the layer 4 site is the peak of the third gamma IC and the layer 6 site lies 15 sites below c; the modulation index is computed without its surrogate test; the eight unit clusters come from k-medians on the full spike-IC coupling profile, and the silhouette does not favour eight; YMV19, without a state table, is left out of every state and LFP analysis.
What was found about the dataset: the data columns are in depth order, and the units' max_channel holds the electrode table id (a hardware channel number) of the unit's peak site, not the column index, which no description says; the visual stimulus times are in 2 of the 13 mice the paper stimulated, and the Visual laser tables of those two also hold the optogenetic pulses; the genotypes in the files (11 PV-Cre::Ai32, 3 VGAT-Cre::Ai32) disagree with the paper (5 and 9); the electrode table gives no depth (y is NaN, location "unknown"); the units carry no waveforms and no layer or cluster label; the cell_type labels (pyramidal, narrow and wide interneuron, "down state") follow a scheme the paper does not describe, and 10 units in 6 mice are labelled "down state" where the paper reports 5 neurons in 4 mice; one file (YMV03) lacks every unit column but quality and one (YMV19) lacks the state table. No recording appears twice, and the two published versions hold the same files.