Claar, Rembado et al. (2023) on dandiset 000458: LFP and CSD magnitudes (Figure 3)¶
Paper. Claar et al. (2023), eLife 12:RP84630, doi:10.7554/eLife.84630;
claims file written from the version of record. Dandiset.
DANDI:000458, version 0.230317.0039, the version the
paper cites. See claar_eeg.ipynb for the EEG claims and the data check of the trial tables, and
claar_spikes.ipynb for the spiking claims and the spiking magnitudes of every site.
This notebook was added after the second verification round, which found that the LFP and CSD parts of Figure 3 were blocked for lack of effort, not of metadata. It runs the LFP and CSD magnitudes of the example session, and times the computation that Figure 3D would need. The third verification round agreed with the outcomes this analysis supports, and the attempt record now gives them: fig3abc partially reproduced, fig3d blocked as too large.
| Claim | Here |
|---|---|
| fig3abc-decay-with-distance | tested (LFP, CSD and spiking of the example session) |
| fig3d-magnitudes | not run: the LFP of every cortical probe exceeds the notebook limit (timed below) |
Readings: the stimulus intensity ("high" and "medium" current), the high-pass cutoff (0.1 and 1 Hz), the baseline of the z-score (-2 to 0 s and -0.5 to 0 s), the onset shuffle (within -2 to +2 s of each onset, "epoch"; or uniform over the recording, run here over the span of the awake pulses, "span"), and the layer-1 exclusion (none, or the top 100 um of the area). For spiking, also the kernel (5 and 10 ms) and the integrated signal (signed z, as the Methods describe it, or |z|).
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import claar
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
urls = claar.asset_urls()
metas = [claar.read_meta(p, u) for p, u in urls]
ex = [m for m in metas if m["subject"] == "571619"][0]
print(f"{len(metas)} sessions from dandiset {claar.DANDISET_ID} version {claar.DANDISET_VERSION}, metadata read in {time.time() - started:.0f} s")
24 sessions from dandiset 000458 version 0.230317.0039, metadata read in 39 s
Data check¶
f = claar.open_nwb(ex["url"])
rows = []
for pr in [p.replace("LFP", "") for p in ex["probes"]]:
ch = claar.lfp_channels(f, pr)
es = f[f"acquisition/LFP{pr}/ElectricalSeries{pr}"]
ts = es["timestamps"][:200000]
for loc, g in ch.groupby("location", sort=False):
rows.append({"probe": pr, "location": loc, "group": g.group.iloc[0], "channels": len(g), "valid": int(g.valid.sum()),
"vertical um": f"{g.vertical.min()}-{g.vertical.max()}"})
print(pr, "LFP shape", es["data"].shape, es["data"].dtype, "; conversion", float(es["data"].attrs["conversion"]), es["data"].attrs["unit"],
"; median interval", round(float(np.median(np.diff(ts))) * 1e3, 4), "ms; chunks", es["data"].chunks)
f.close()
lfp_tab = pd.DataFrame(rows)
print(lfp_tab.to_string(index=False))
print("\nexample session", ex["subject"], ex["session"], "; sites", sorted(ex["estim"].site.unique()))
probeB LFP shape (19863586, 384) int16 ; conversion 1.9499999999999999e-07 volts ; median interval 0.4 ms; chunks (77593, 2)
probeC LFP shape (19863628, 384) int16 ; conversion 1.9499999999999999e-07 volts ; median interval 0.4 ms; chunks (77593, 2)
probeF LFP shape (19863581, 384) int16 ; conversion 1.9499999999999999e-07 volts ; median interval 0.4 ms; chunks (77593, 2) probe location group channels valid vertical um probeB VPL SM-TH 34 34 20-340 probeB LD NaN 58 54 360-920 probeB root NaN 4 4 940-960 probeB fi NaN 32 32 980-1280 probeB alv NaN 12 11 1300-1940 probeB CA3 NaN 54 54 1380-1900 probeB ccb NaN 10 10 1960-2040 probeB scwm NaN 8 8 2060-2120 probeB SSp-ll SS 84 32 2140-2960 probeB SSp-tr SS 46 37 2980-3420 probeB none NaN 42 0 3440-3840 probeC ZI NaN 6 6 20-60 probeC TH NaN 14 14 80-1140 probeC MGm NaN 18 18 160-320 probeC MGv NaN 42 42 340-740 probeC MGd NaN 22 22 760-960 probeC SGN NaN 4 4 980-1000 probeC LP NaN 8 8 1060-1120 probeC root NaN 8 8 1160-1220 probeC DG-mo NaN 38 38 1240-1900 probeC DG-sg NaN 18 18 1420-1700 probeC DG-po NaN 12 12 1520-1620 probeC CA1 NaN 38 37 1920-2280 probeC ProS NaN 18 18 2300-2460 probeC alv NaN 2 2 2480-2480 probeC fp NaN 10 10 2500-2580 probeC or NaN 6 6 2600-2640 probeC VISp VIS 88 87 2660-3520 probeC none NaN 32 0 3540-3840 probeF ccg NaN 60 60 20-600 probeF ACAv ACA 72 71 620-1320 probeF ACAd ACA 34 34 1340-1660 probeF MOs MO 120 119 1680-2860 probeF none NaN 98 0 2880-3840 example session 571619 20210319 ; sites ['MOs-deep', 'MOs-superficial']
The example session (571619, the deep MOs subject of Figure 1D) has three probes; MOs (120 channels, 119 marked valid) on probe F and VISp (88, 87 valid) on probe C, each with channels in the ACSF above the brain ("none", 98 and 32), which the file marks invalid. The LFP is int16 at 2.5 kHz (0.4 ms intervals) with a conversion to volts. Channel locations are CCF areas without layers; layer 1 is approximated below by the top 100 um of each area, the channels nearest the ACSF, as an assumption, and the result is shown with and without it. The clock of the LFP is checked below with the stimulus artifact at the trial onsets.
fig3abc-decay-with-distance¶
In the example deep MOs session, evoked LFP, CSD and spiking magnitudes are large near the stimulation site (54.0, 51.4 and 30.9 SD.s in MO) and smaller in distant VIS (18.0, 4.8 and 1.8 SD.s).
Criterion: magnitudes as in fig3d-magnitudes for the session of fig1-supp1-example; reproduced when MO exceeds VIS for all three signals and all six values are within tolerance under one combination of readings; partially reproduced when the ordering holds but a value is out of tolerance.
LFP and CSD by claar.lfp_magnitudes (all awake trials, 1000 shuffles; the steps and where each comes from are in
its docstring). The spiking magnitudes are those of claar.spike_site_summary, as in the spiking notebook.
t0 = time.time()
timing, parts = {}, []
for k, (probe, group) in enumerate((("probeF", "MO"), ("probeC", "VIS"))):
tm = {}
parts.append(claar.lfp_magnitudes(ex, "MOs-deep", probe, (group,), np.random.default_rng([304, 571619, k]), timing=tm))
timing[(probe, group)] = tm
LFP = pd.concat(parts, ignore_index=True)
print(f"LFP and CSD magnitudes in {time.time() - t0:.0f} s")
for key, tm in timing.items():
print(key, {k: (round(v, 1) if isinstance(v, float) else v) for k, v in tm.items()})
print(LFP[["group", "level", "channels", "rejected", "trials"]].drop_duplicates().to_string(index=False))
LFP and CSD magnitudes in 335 s
('probeF', 'MO') {'read s': 68.5, 'total s': 190.7, 'samples': 3611371, 'columns': 217, 'artifact uV, k-3..k+5': array([ 28.7, 30.4, 31.6, 31.7, 36.3, 255.5, 141.2, 59. , 40. ])}
('probeC', 'VIS') {'read s': 54.1, 'total s': 144.8, 'samples': 3611380, 'columns': 119, 'artifact uV, k-3..k+5': array([36.9, 38.1, 40.6, 40.7, 49. , 79.7, 51.8, 45.6, 45.2])}
group level channels rejected trials
MO high 119 0 120
MO medium 119 0 118
VIS high 87 0 120
VIS medium 87 0 118
Clock check for the LFP: the median raw artifact over the area channels (uV above the mean of k - 25 to k - 11, k the first LFP sample at or after the onset) is printed for k - 3 to k + 5. It peaks at k + 2 and k + 3, inside the mask (k - 1 to k + 4), and is far above the samples before the onset, so the LFP shares the trial table's clock; its artifact lags the EEG's by about 1 ms. The rejected count is of channels replaced by the channel-level Chebyshev rule.
t0 = time.time()
tab, sp = claar.read_units(ex)
rsum = claar.spike_site_summary(ex, "MOs-deep", tab, sp)
SPK = rsum["magnitude"].query("group in ['MO', 'VIS']")
print(f"spiking summary in {time.time() - t0:.0f} s; RS units: MO {int(SPK.query('group == \"MO\"').units.iloc[0])}, VIS {int(SPK.query('group == \"VIS\"').units.iloc[0])}")
print("\nLFP and CSD magnitudes (expected MO 54.0 and 51.4, VIS 18.0 and 4.8; shuffle 'none' is the AUC itself in SD.s):")
print(LFP.pivot_table(index=["level", "highpass Hz", "layer1 um", "baseline"], columns=["signal", "shuffle", "group"], values="magnitude").round(2).to_string())
print("\nspiking magnitudes (expected MO 30.9, VIS 1.8):")
print(SPK.pivot_table(index=["level", "kernel ms", "baseline"], columns=["signal", "shuffle", "group"], values="magnitude").round(2).to_string())
spiking summary in 42 s; RS units: MO 95, VIS 47
LFP and CSD magnitudes (expected MO 54.0 and 51.4, VIS 18.0 and 4.8; shuffle 'none' is the AUC itself in SD.s):
signal CSD LFP
shuffle epoch none span epoch none span
group MO VIS MO VIS MO VIS MO VIS MO VIS MO VIS
level highpass Hz layer1 um baseline
high 0.1 0 -0.5-0 59.19 9.94 21.88 3.43 54.10 9.46 71.56 23.77 24.19 8.50 64.84 23.25
-2-0 49.87 11.68 15.21 2.88 45.29 11.58 77.11 21.42 17.84 5.15 76.77 20.43
100 -0.5-0 58.79 9.27 21.63 3.20 53.73 8.80 71.66 23.79 24.19 8.53 64.86 23.19
-2-0 48.79 11.21 14.90 2.76 44.32 11.07 77.16 21.27 17.84 5.12 76.84 20.27
1.0 0 -0.5-0 19.76 6.81 4.61 1.79 18.26 6.58 8.49 9.76 2.42 2.78 7.59 9.71
-2-0 45.72 11.35 7.77 2.30 41.77 11.10 21.02 19.58 4.38 4.01 20.29 18.42
100 -0.5-0 20.06 6.48 4.64 1.70 18.65 6.29 8.46 9.68 2.41 2.78 7.57 9.64
-2-0 46.14 10.83 7.75 2.18 42.15 10.62 20.96 19.37 4.37 3.99 20.23 18.23
medium 0.1 0 -0.5-0 52.34 5.36 20.16 1.90 49.92 5.32 66.56 11.85 24.41 4.34 70.05 11.26
-2-0 49.93 5.95 14.95 1.46 46.64 5.82 61.05 22.06 14.43 5.32 62.23 21.94
100 -0.5-0 49.50 4.72 19.06 1.68 47.34 4.70 66.53 11.88 24.45 4.37 70.21 11.34
-2-0 48.77 5.21 14.63 1.28 45.46 5.09 60.83 21.75 14.38 5.26 61.99 21.67
1.0 0 -0.5-0 21.85 4.54 5.53 1.28 22.73 4.86 9.05 6.35 2.88 1.86 9.15 6.16
-2-0 48.85 6.53 8.91 1.30 48.13 6.33 22.59 14.41 4.94 3.04 22.26 14.00
100 -0.5-0 22.01 4.07 5.51 1.15 22.92 4.33 9.01 6.29 2.87 1.85 9.11 6.10
-2-0 49.21 5.81 8.87 1.15 48.44 5.62 22.52 14.23 4.93 3.02 22.19 13.83
spiking magnitudes (expected MO 30.9, VIS 1.8):
signal signed |z|
shuffle epoch none recording epoch none recording
group MO VIS MO VIS MO VIS MO VIS MO VIS MO VIS
level kernel ms baseline
high 5 -0.5-0 1.46 -1.29 0.34 -0.35 1.29 -1.34 11.22 1.30 1.75 0.60 11.19 1.27
-2-0 3.96 -1.21 0.70 -0.33 3.09 -1.35 28.31 2.11 2.96 0.60 29.19 2.08
10 -0.5-0 0.59 -1.18 0.14 -0.43 0.44 -1.18 3.57 0.52 1.08 0.59 3.45 0.43
-2-0 2.30 -1.20 0.42 -0.41 1.58 -1.30 13.02 1.37 2.02 0.61 13.70 1.26
medium 5 -0.5-0 1.18 -1.24 0.40 -0.36 1.45 -1.33 13.14 0.97 2.11 0.56 15.21 0.95
-2-0 2.95 -2.14 0.71 -0.52 3.13 -2.16 28.75 2.76 3.04 0.64 31.37 2.65
10 -0.5-0 0.31 -1.62 0.21 -0.61 0.58 -1.64 5.04 2.09 1.41 0.84 5.91 1.92
-2-0 1.74 -2.31 0.53 -0.70 1.95 -2.25 16.50 3.17 2.41 0.80 17.76 2.91
exp = {("LFP", "MO"): (54.0, 5.4), ("LFP", "VIS"): (18.0, 1.8), ("CSD", "MO"): (51.4, 5.14), ("CSD", "VIS"): (4.8, 0.48),
("spiking", "MO"): (30.9, 3.09), ("spiking", "VIS"): (1.8, 0.18)}
L = LFP[LFP.shuffle != "none"].pivot_table(index=["level", "highpass Hz", "layer1 um", "baseline", "shuffle"], columns=["signal", "group"], values="magnitude")
S = SPK[SPK.shuffle != "none"].replace({"shuffle": {"recording": "span"}}).pivot_table(
index=["level", "baseline", "shuffle", "kernel ms", "signal"], columns="group", values="magnitude")
rows = []
for (lvl, hp, l1, base, sh), lr in L.iterrows():
for (kern, sig), sr in S.loc[(lvl, base, sh)].iterrows():
v = {("LFP", "MO"): lr[("LFP", "MO")], ("LFP", "VIS"): lr[("LFP", "VIS")], ("CSD", "MO"): lr[("CSD", "MO")],
("CSD", "VIS"): lr[("CSD", "VIS")], ("spiking", "MO"): sr["MO"], ("spiking", "VIS"): sr["VIS"]}
inside = [k for k in exp if abs(v[k] - exp[k][0]) <= exp[k][1]]
rows.append({"level": lvl, "highpass Hz": hp, "layer1 um": l1, "baseline": base, "shuffle": sh, "kernel ms": kern,
"spiking signal": sig, **{f"{a} {b}": round(x, 2) for (a, b), x in v.items()},
"MO > VIS all three": all(v[(s, "MO")] > v[(s, "VIS")] for s in ("LFP", "CSD", "spiking")),
"within tolerance": len(inside), "inside": ", ".join(f"{a} {b}" for a, b in inside)})
crit = pd.DataFrame(rows)
print(f"{len(crit)} combinations of readings; MO > VIS for all three signals in {int(crit['MO > VIS all three'].sum())}; "
f"all six within tolerance in {int((crit['within tolerance'] == 6).sum())}; largest number within tolerance {crit['within tolerance'].max()}")
for c in ("LFP MO", "LFP VIS", "CSD MO", "CSD VIS", "spiking MO", "spiking VIS"):
print(f" {c}: {crit[c].min():.2f} to {crit[c].max():.2f}; within tolerance in {int(crit.inside.str.contains(c).sum())} of {len(crit)}")
print("\nordering by signal, over the LFP and CSD readings:", {s: f"{int((L[(s, 'MO')] > L[(s, 'VIS')]).sum())}/{len(L)}" for s in ("LFP", "CSD")},
"; spiking:", f"{int((S.MO > S.VIS).sum())}/{len(S)}")
print("\ncombinations with the most values within tolerance:")
print(crit.sort_values("within tolerance", ascending=False).head(12).to_string(index=False))
print("\nmedium current, 0.1 Hz, -2 to 0 s, epoch shuffle, 5 ms kernel:")
print(crit.query("level == 'medium' and `highpass Hz` == 0.1 and baseline == '-2-0' and shuffle == 'epoch' and `kernel ms` == 5").to_string(index=False))
128 combinations of readings; MO > VIS for all three signals in 112; all six within tolerance in 0; largest number within tolerance 3
LFP MO: 7.57 to 77.16; within tolerance in 0 of 128
LFP VIS: 6.10 to 23.79; within tolerance in 16 of 128
CSD MO: 18.26 to 59.19; within tolerance in 60 of 128
CSD VIS: 4.07 to 11.68; within tolerance in 28 of 128
spiking MO: 0.31 to 31.37; within tolerance in 16 of 128
spiking VIS: -2.31 to 3.17; within tolerance in 4 of 128
ordering by signal, over the LFP and CSD readings: {'LFP': '28/32', 'CSD': '32/32'} ; spiking: 32/32
combinations with the most values within tolerance:
level highpass Hz layer1 um baseline shuffle kernel ms spiking signal LFP MO LFP VIS CSD MO CSD VIS spiking MO spiking VIS MO > VIS all three within tolerance inside
medium 0.1 100 -0.5-0 span 10 |z| 70.21 11.34 47.34 4.70 5.91 1.92 True 3 CSD MO, CSD VIS, spiking VIS
medium 0.1 100 -2-0 epoch 5 |z| 60.83 21.75 48.77 5.21 28.75 2.76 True 3 CSD MO, CSD VIS, spiking MO
high 0.1 0 -2-0 epoch 5 |z| 77.11 21.42 49.87 11.68 28.31 2.11 True 2 CSD MO, spiking MO
high 0.1 100 -2-0 epoch 5 |z| 77.16 21.27 48.79 11.21 28.31 2.11 True 2 CSD MO, spiking MO
high 1.0 100 -2-0 span 5 |z| 20.23 18.23 42.15 10.62 29.19 2.08 True 2 LFP VIS, spiking MO
high 1.0 100 -2-0 epoch 5 |z| 20.96 19.37 46.14 10.83 28.31 2.11 True 2 LFP VIS, spiking MO
high 1.0 0 -2-0 epoch 5 |z| 21.02 19.58 45.72 11.35 28.31 2.11 True 2 LFP VIS, spiking MO
high 1.0 0 -2-0 span 5 |z| 20.29 18.42 41.77 11.10 29.19 2.08 True 2 LFP VIS, spiking MO
medium 0.1 100 -2-0 span 5 |z| 61.99 21.67 45.46 5.09 31.37 2.65 True 2 CSD VIS, spiking MO
medium 0.1 0 -2-0 span 5 |z| 62.23 21.94 46.64 5.82 31.37 2.65 True 2 CSD MO, spiking MO
medium 1.0 0 -0.5-0 span 10 |z| 9.15 6.16 22.73 4.86 5.91 1.92 True 2 CSD VIS, spiking VIS
medium 1.0 0 -2-0 epoch 5 |z| 22.59 14.41 48.85 6.53 28.75 2.76 True 2 CSD MO, spiking MO
medium current, 0.1 Hz, -2 to 0 s, epoch shuffle, 5 ms kernel:
level highpass Hz layer1 um baseline shuffle kernel ms spiking signal LFP MO LFP VIS CSD MO CSD VIS spiking MO spiking VIS MO > VIS all three within tolerance inside
medium 0.1 0 -2-0 epoch 5 signed 61.05 22.06 49.93 5.95 2.95 -2.14 True 1 CSD MO
medium 0.1 0 -2-0 epoch 5 |z| 61.05 22.06 49.93 5.95 28.75 2.76 True 2 CSD MO, spiking MO
medium 0.1 100 -2-0 epoch 5 signed 60.83 21.75 48.77 5.21 2.95 -2.14 True 2 CSD MO, CSD VIS
medium 0.1 100 -2-0 epoch 5 |z| 60.83 21.75 48.77 5.21 28.75 2.76 True 3 CSD MO, CSD VIS, spiking MO
fig, ax = plt.subplots(1, 3, figsize=(12, 3.2))
for k, s in enumerate(("LFP", "CSD")):
d = LFP[(LFP.signal == s) & (LFP.shuffle == "epoch")]
for j, g in enumerate(("MO", "VIS")):
x = d[d.group == g].magnitude
ax[k].scatter(np.full(len(x), j) + np.random.default_rng(1).uniform(-0.15, 0.15, len(x)), x, s=8, color="0.4")
ax[k].hlines(exp[(s, g)][0], j - 0.3, j + 0.3, color="r")
ax[k].set(xticks=[0, 1], xticklabels=["MO", "VIS"], title=f"{s}, epoch shuffle, all readings", ylabel="magnitude")
d = SPK[SPK.shuffle == "epoch"]
for j, g in enumerate(("MO", "VIS")):
for sig, col in (("signed", "0.6"), ("|z|", "k")):
x = d[(d.group == g) & (d.signal == sig)].magnitude
ax[2].scatter(np.full(len(x), j) + (0.12 if sig == "|z|" else -0.12), x, s=8, color=col)
ax[2].hlines(exp[("spiking", g)][0], j - 0.3, j + 0.3, color="r")
ax[2].set(xticks=[0, 1], xticklabels=["MO", "VIS"], title="spiking, epoch shuffle (grey signed, black |z|)")
plt.tight_layout()
print("plotted ranges: LFP", round(LFP[(LFP.signal == "LFP") & (LFP.shuffle == "epoch")].magnitude.min(), 1), "to",
round(LFP[(LFP.signal == "LFP") & (LFP.shuffle == "epoch")].magnitude.max(), 1), "; CSD",
round(LFP[(LFP.signal == "CSD") & (LFP.shuffle == "epoch")].magnitude.min(), 1), "to",
round(LFP[(LFP.signal == "CSD") & (LFP.shuffle == "epoch")].magnitude.max(), 1), "; spiking",
round(d.magnitude.min(), 1), "to", round(d.magnitude.max(), 1), "(no points outside the axes; red, the paper's values)")
plotted ranges: LFP 6.3 to 77.2 ; CSD 4.1 to 59.2 ; spiking -2.3 to 28.8 (no points outside the axes; red, the paper's values)
Result: partially reproduced (agreed in the third verification round). Only the ordering is reproduced. MO exceeds VIS for all three signals in 112 of 128 combinations of readings; the LFP order reverses in 4 of 32 LFP readings, all at the high current with the 1 Hz cutoff and the -0.5 to 0 s baseline. No combination puts all six values within tolerance, and at most three are (medium current, 0.1 Hz, top 100 um left out, -2 to 0 s, epoch shuffle, 5 ms kernel, |z|: CSD 48.77 and 5.21, spiking MO 28.75; LFP 60.83 and 21.75 against 54.0 and 18.0, spiking VIS 2.76 against 1.8). The LFP magnitude in MO is never within tolerance (7.57 to 77.16; 60.83 and above with the 0.1 Hz cutoff). The CSD magnitude in MO is within tolerance in 60 of 128 combinations. The spiking magnitude as the Methods describe it (signed z) is 0.31 to 3.96 in MO; integrating |z| gives 3.57 to 31.37, and 28.31 to 31.37 with the 5 ms kernel and the -2 to 0 s baseline, near the paper's 30.9, so the paper's spiking magnitude probably counts falls below baseline as well as rises ("our magnitude metric accounts for both reductions as well as increases"). Leaving out the top 100 um changes no ordering and moves the CSD values by up to 2.84 (52.34 against 49.50, medium current, 0.1 Hz, -0.5 to 0 s, epoch shuffle). The layer-1 exclusion is an assumption (the files give no layers), and the onset shuffle over the whole recording is replaced by one over the span of the awake pulses.
fig3d-magnitudes: the size of the computation¶
Superficial stimulation modulated local LFP, CSD and spiking and non-local LFP and CSD, but evoked no spiking outside the stimulated area; deep stimulation evoked LFP, CSD and spiking both locally and in non-stimulated cortical regions.
The test needs the LFP and CSD magnitudes of every cortical area on every probe of every stimulation site. The cost of one probe was measured above. It grows with the number of channels read (the area's valid channels and the ACSF channels) and the samples of the awake span. Below, every site with valid awake trials and every probe area with at least 6 valid channels are counted, with the time each would take at the measured rate. The spiking part is in the spiking notebook.
rate = np.mean([tm["total s"] / (tm["columns"] * tm["samples"]) for tm in timing.values()])
read_rate = np.mean([tm["read s"] / (tm["columns"] * tm["samples"]) for tm in timing.values()])
rows = []
for m in metas:
el = m["electrodes"]
probes = [p.replace("LFP", "") for p in m["probes"]]
for site in sorted(m["estim"].site.unique()):
on = m["estim"].query("site == @site and is_valid and state in ['quiet', 'running']").start_time.values
if len(on) == 0 or not probes:
continue
samples = (on.max() - on.min() + 7.5) * claar.FS_LFP
for pr in probes:
e = el[el.group_name == pr]
groups = e[e.is_data_valid].location.map(claar.area_group)
n_acsf = int((e.location == "none").sum())
for g, n in groups.value_counts().items():
if g in claar.CORTEX and n >= 6:
rows.append({"subject": m["subject"], "site": site, "probe": pr, "group": g, "channels": int(n), "acsf": n_acsf,
"samples": samples, "minutes": (n + n_acsf) * samples * rate / 60,
"read minutes": (n + n_acsf) * samples * read_rate / 60})
est = pd.DataFrame(rows)
once = est.groupby(["subject", "site", "probe"]).agg(ch=("channels", "sum"), acsf=("acsf", "first"), samples=("samples", "first"))
read_once = ((once.ch + once.acsf) * once.samples * read_rate / 60).sum()
print(est.groupby(["site"]).agg(areas=("group", "size"), minutes=("minutes", "sum")).round(1).to_string())
print(f"\n{est[['subject', 'site']].drop_duplicates().shape[0]} sites with probes and valid awake trials, {len(est)} probe areas; "
f"estimated {est.minutes.sum():.0f} min at the measured rate ({rate * 1e9:.2f} s per 1e9 channel-samples, all readings), "
f"of which reading the LFP {est['read minutes'].sum():.0f} min ({read_once:.0f} min if the areas of one probe share a read); "
f"the notebook limit is 30 min")
print(f"elapsed in this notebook so far: {(time.time() - started) / 60:.1f} min")
areas minutes site MOs-deep 31 101.5 MOs-superficial 8 25.7 SSp-deep 30 82.8 SSp-superficial 17 45.3 24 sites with probes and valid awake trials, 86 probe areas; estimated 255 min at the measured rate (290.10 s per 1e9 channel-samples, all readings), of which reading the LFP 94 min (78 min if the areas of one probe share a read); the notebook limit is 30 min elapsed in this notebook so far: 7.0 min
Result: not run (too large). Figure 3D needs 86 probe areas over 24 stimulation sites. At the rate measured on the example's two probes (printed above; it varies with the network from run to run), their LFP and CSD magnitudes would take several hours with all readings, and reading their LFP alone more than an hour, also when the areas of one probe share a read. Both are over the 30 min notebook limit. Layer 1 is approximated as in fig3abc and does not block the claim. Split across notebooks of under 30 min, the computation would need about eight; that is left to a later attempt. The third verification round agreed with this outcome. The estimate is rough: it scales the shuffles, which grow with trials and channels, by samples as well, and it counts the ACSF channels once per area; the read-once figure is still more than twice the limit.
Summary¶
| Claim | Outcome |
|---|---|
| fig3abc-decay-with-distance | partially reproduced (assumption: layer 1 is the top 100 um; only the ordering holds) |
| fig3d-magnitudes | blocked: too large |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's: the files give no cortical layers, so layer 1 is approximated by the top 100 um of each area (an assumption; the result is shown without the exclusion too); the Chebyshev channel rejection is applied per channel, once, where the Methods apply it per sample and epoch; the CSD is computed from the trial-averaged LFP (all its steps are linear, so this equals averaging single-trial CSDs) with rows of two channels, 13 then 2 rows of smoothing for the 26 then 4 channels of the Methods; the decimation and high-pass filters are zero phase; the shuffle "over the recording" draws onsets over the span of the awake pulses; the mask is the EEG's (k - 1 to k + 4), which covers the LFP artifact at k + 2 and k + 3. The current is not given; the high and the medium current are both run.
What was found about the dataset: the LFP of every probe is int16 at 2.5 kHz with timestamps that share the trial table's clock (the artifact peaks 2 samples after the onset sample); the channels in the ACSF above the brain are labelled "none" and marked invalid; channel locations carry no layers.
print(f"{(time.time() - started) / 60:.1f} min")
7.0 min