Siegle, Jia et al. (2021): the cross-correlogram claims¶
Paper. Siegle JH, Jia X, et al. (2021). Survey of spiking in the mouse visual system reveals functional hierarchy. Nature 592, 86-92. doi:10.1038/s41586-020-03171-x. Claims file written from the PMC author manuscript (PMC10399640).
Data. The target dandiset, DANDI:000253
(version 0.240923.1441), holds a different experiment (see reproduce_siegle_2021.ipynb). The CCG
analyses use the Functional Connectivity sessions, which the paper deposited as
DANDI:000022, version 0.251116.2247, its
only published version. The paper's per-unit layer table (data/layer_info.csv of
AllenInstitute/neuropixels_platform_paper at commit 4f8bcc26) labels layers for ed10ce; its link to
the dandisets' unit ids is tested in the other notebook (link-unit-tables).
| Claim | Here |
|---|---|
| fig2-sharp-peak-pairs | tested |
| fig2d-v1-leads-lm-example | tested |
| fig2-neighbour-offsets-across-mice | tested |
| ed7c-v1-leads-am-follows | tested |
| fig2f-directionality-vs-anatomy | tested |
| ed7e-distance-does-not-explain | blocked: the left-right CCF coordinate is not in the session or probe files |
| ed7fg-absent-in-spontaneous | tested |
| ed8-total-hierarchy-score | tested |
| ed10ce-superficial-lead-deep | tested |
Method (siegle.py, "cross-correlograms"): 1 ms bins; per orientation (0, 45, 90, 135 degrees, contrast 0.8, 75 trials) the CCG of every pair of units of V1, LM, RL, AL, PM and AM with a mean rate above 2 Hz from 50 to 500 ms after onset, normalized by the triangular overlap and the geometric mean rate, minus its expected value under jitter of both trains within 25 ms windows, averaged over the four orientations; a sharp peak when the maximum within +/-10 ms exceeds 7 times the SD of the flanks at 50 to 100 ms. Readings named by the claims file: window 2 s (the presentation) or 3 s (presentation and grey); units (a) all non-noise, (b) with the Extended Data Fig. 1j filters (ISI violations < 0.5, amplitude cutoff < 0.1, presence ratio > 0.95), which only removes pairs. The claims that use the sharp peaks are judged under the reading of fig2-sharp-peak-pairs that matches the paper, or, when none does, the one nearest it (printed there); the other readings are shown as variants. The spontaneous period (30 min grey screen) is cut into consecutive 2 s or 3 s segments, with the same units.
Revision after the first verification. The expected CCG under jitter was wrong at negative lags in the first version of the helper (it used the lag +k window product where the lag -k one belongs); it is fixed, and a check against a direct enumeration of spike pairs is added. A within-area reference distribution now counts each pair once.
import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import io
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
from scipy import stats
import siegle as S
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
assets = S.list_assets("000022", S.SOURCES["000022"])
fc = assets[~assets.path.str.contains("probe")].reset_index(drop=True)
res = S.load_ccg_all(fc)
layer_csv = pd.read_csv(io.StringIO(requests.get(f"{S.CODE_RAW}/data/layer_info.csv", timeout=60).text))
print(f"{len(res)} sessions in {(time.time() - started) / 60:.1f} min")
26 sessions in 15.5 min
Data check¶
info = pd.DataFrame([{"subject": r["subject"], "session": r["session"], "units_above_2Hz": len(r["units"]),
"dg_trials_contrast_0.8": r["n_dg_trials"], "orientations": len(r["orientations"]),
"dg_onset_interval_s": round(r["dg_onset_interval"], 3),
"spont_s": round(r.get("spont_duration", np.nan), 1),
"segments_2s": r.get("n_spont_segments_2000"), "seconds": round(r["seconds"])} for r in res])
print(info.to_string(index=False))
print("\nmice:", info.subject.nunique(), " units:", info.units_above_2Hz.sum())
allu = pd.concat([r["units"].assign(subject=r["subject"]) for r in res])
print("units above 2 Hz by area:", allu.location.value_counts().reindex(S.CORTEX).to_dict())
print("rate from 50 to 500 ms (Hz), range:", round(allu.rate_dg.min(), 2), "to", round(allu.rate_dg.max(), 2))
subject session units_above_2Hz dg_trials_contrast_0.8 orientations dg_onset_interval_s spont_s segments_2s seconds
744912845 766640955 357 300 4 6.005 1802.5 901 104
753795604 767871931 400 300 4 6.005 1802.5 901 154
754477355 768515987 416 300 4 6.005 1802.5 901 152
754488976 771160300 377 300 4 6.005 1802.5 901 133
756578432 771990200 325 300 4 3.003 1802.5 901 116
759674767 778998620 484 300 4 6.005 1802.5 901 174
759711149 774875821 425 300 4 6.005 1802.5 901 132
760938794 778240327 513 300 4 6.005 1802.5 901 197
760946810 781842082 403 300 4 3.003 1802.5 901 119
760960650 779839471 307 300 4 3.003 1802.5 901 102
763236006 787025148 280 300 4 3.003 1802.5 901 78
763808588 789848216 192 300 4 3.003 1802.5 901 52
763884096 786091066 523 300 4 3.003 1802.5 901 249
769319621 793224716 309 300 4 3.003 1802.5 901 87
774672354 794812542 518 300 4 6.005 1802.5 901 203
791857600 816200189 320 300 4 6.005 1802.5 901 110
795770029 819701982 349 300 4 3.003 1802.5 901 131
800249582 819186360 390 300 4 3.003 1802.6 901 128
800250054 821695405 345 300 4 6.005 1802.5 901 116
803390283 831882777 513 300 4 3.003 1802.5 901 224
811322604 829720705 408 300 4 6.005 1802.5 901 158
813701555 835479236 359 300 4 6.005 1802.5 901 138
817060743 839068429 436 300 4 3.019 1802.5 901 148
820866111 840012044 440 300 4 3.003 1802.5 901 148
821469659 839557629 376 300 4 6.005 1802.5 901 107
827809879 847657808 518 300 4 6.005 1802.5 901 109
mice: 26 units: 10283
units above 2 Hz by area: {'VISp': 2196, 'VISl': 1448, 'VISrl': 1645, 'VISal': 1850, 'VISpm': 1248, 'VISam': 1896}
rate from 50 to 500 ms (Hz), range: 2.01 to 151.55
Every session has 300 trials at contrast 0.8 over the four orientations (75 each); contrast 0.8 and 0.1 trials alternate, so successive 0.8 onsets are 6 s apart and each trial has 2 s of grating and at least 1 s of grey after it. The 30 min spontaneous block is present in each session. The paper reports 25 mice for the CCG analyses; the dandiset holds 26 Functional Connectivity sessions and nothing identifies the one left out, so all 26 are used.
Check of the jitter correction¶
S.ccg_group against a direct enumeration of spike pairs on synthetic trains: the raw CCG and its
expected value under independent uniform jitter of both trains within 25 ms windows, normalized the
same way. Unit 1's rate rises later in the trial than unit 0's, so an error in the expected value at
negative lags would show. The round-1 version of the helper failed this check at negative lags.
def brute_force_ccg(spikes, onsets, T):
"""Jitter-corrected, normalized CCG of unit 1 relative to unit 0, by direct enumeration of spike pairs."""
lags = np.arange(-S.MAXLAG, S.MAXLAG + 1)
raw, expv = np.zeros(len(lags)), np.zeros(len(lags))
w = np.arange(S.JW)
diff = (w[None, :] - w[:, None]).ravel() # q - p for every pair of positions in two windows
n_sp = [0, 0]
for on in onsets:
b = [np.floor((st[(st >= on) & (st < on + T / 1000)] - on) * 1000).astype(int) for st in spikes[:2]]
n_sp[0] += len(b[0])
n_sp[1] += len(b[1])
for s in b[0]:
for s2 in b[1]:
if abs(s2 - s) <= S.MAXLAG:
raw[s2 - s + S.MAXLAG] += 1
dd = (s2 // S.JW - s // S.JW) * S.JW + diff # jittered differences, each with probability 1/625
dd = dd[np.abs(dd) <= S.MAXLAG]
np.add.at(expv, dd + S.MAXLAG, 1 / S.JW ** 2)
lam = np.array(n_sp) / (len(onsets) * T)
return (raw - expv) / (len(onsets) * (T - np.abs(lags))) / np.sqrt(lam[0] * lam[1])
rng_check = np.random.default_rng(7)
on_check = np.arange(20) * 1.0
sp_check = []
for rise in (0.02, 0.08):
st = [on + rng_check.uniform(0, 0.2, rng_check.poisson(8)) for on in on_check]
st += [on + rise + rng_check.uniform(0, 0.05, rng_check.poisson(6)) for on in on_check]
sp_check.append(np.sort(np.concatenate(st)))
c_check = S.ccg_group(sp_check, on_check, 200)[0, 1]
bf_check = brute_force_ccg(sp_check, on_check, 200)
lags_check = np.arange(-S.MAXLAG, S.MAXLAG + 1)
used = np.isfinite(c_check)
err_neg = np.max(np.abs(c_check - bf_check)[used & (lags_check < 0)])
err_pos = np.max(np.abs(c_check - bf_check)[used & (lags_check >= 0)])
print(f"max |helper - brute force|: negative lags {err_neg:.2e}, zero and positive lags {err_pos:.2e}")
print("expected CCG under jitter at -10 and +10 ms, brute force:", np.round(bf_check[[S.MAXLAG - 10, S.MAXLAG + 10]], 4),
" helper:", np.round(c_check[[S.MAXLAG - 10, S.MAXLAG + 10]], 4))
assert err_neg < 1e-6 and err_pos < 1e-6
max |helper - brute force|: negative lags 5.48e-09, zero and positive lags 4.89e-09 expected CCG under jitter at -10 and +10 ms, brute force: [ 0.0273 -0.0035] helper: [ 0.0273 -0.0035]
AREA_ORDER = ["VISp", "VISl", "VISrl", "VISal", "VISpm", "VISam"] # anatomical order of the cortex
NAMES = {a: S.NAMES[a] for a in AREA_ORDER}
HSC = np.array([S.HIERARCHY[a] for a in AREA_ORDER])
def ed1j(u):
return ((u.isi_violations < 0.5) & (u.amplitude_cutoff < 0.1) & (u.presence_ratio > 0.95)).values
def peak_table(reading, quality, seed=253):
# One row per significant pair and direction: source unit, target unit, offset (source leads when > 0).
# A within-area pair appears twice, with its offset and the negative. `ref` marks one of the two
# copies, chosen at random (a standard choice: the sign of a within-area offset has no meaning), so a
# within-area reference distribution counts each pair once; `first` marks the copy in file order, a
# variant. The two copies of a between-area pair are different area pairs, and both are marked.
rng = np.random.default_rng(seed)
rows, npairs = [], []
for r in res:
u = r["units"]
ok = np.ones(len(u), bool) if quality == "all" else ed1j(u)
p = r["peaks"][reading]
p = p[ok[p.i.values] & ok[p.j.values]]
n_ok = ok.sum()
npairs.append({"subject": r["subject"], "pairs": n_ok * (n_ok - 1) // 2, "peaks": len(p)})
loc = u.location.values
same = loc[p.i.values] == loc[p.j.values]
flip = rng.random(len(p)) < 0.5
for sgn, a, b in ((1, "i", "j"), (-1, "j", "i")):
rows.append(pd.DataFrame({"subject": r["subject"], "src": loc[p[a].values], "tgt": loc[p[b].values],
"src_id": u.unit_id.values[p[a].values], "tgt_id": u.unit_id.values[p[b].values],
"offset": sgn * p.offset.values,
"ref": ~same | (flip if sgn < 0 else ~flip), "first": ~same | (sgn > 0)}))
return pd.concat(rows, ignore_index=True), pd.DataFrame(npairs)
def ds(offsets, zero="exclude"):
pos = np.sum(offsets > 0) + (np.sum(offsets == 0) if zero == "positive" else 0)
neg = np.sum(offsets < 0)
return (pos - neg) / (pos + neg) if pos + neg else np.nan
def ds_matrix(t, mode="pooled", zero="exclude"):
# rows: target, columns: source (as Figure2/ccg_distribution_model_final.ipynb)
M = np.full((6, 6), np.nan)
for i, tg in enumerate(AREA_ORDER):
for j, sr in enumerate(AREA_ORDER):
d = t[(t.src == sr) & (t.tgt == tg)]
if mode == "pooled":
M[i, j] = ds(d.offset.values, zero)
else:
M[i, j] = np.nanmean([ds(g.offset.values, zero) for _, g in d.groupby("subject")]) if len(d) else np.nan
return M
def anatomy_corr(M):
A = HSC[:, None] - HSC[None, :] # target score minus source score
il = np.tril_indices(6)
x, y = A[il], M[il]
ok = np.isfinite(y)
r, p = stats.pearsonr(x[ok], y[ok])
return r, p, int(ok.sum())
READINGS = [(w, q) for w in ("grating_2000", "grating_3000") for q in ("all", "ed1j")]
tables = {rq: peak_table(*rq) for rq in READINGS}
fig2-sharp-peak-pairs¶
Jitter-corrected CCGs during drifting gratings showed a significant sharp peak for 16,119 of 2,089,890 possible cortical unit pairs, 0.96% +/- 0.13% per mouse in 25 mice.
rows = []
for (w, q), (t, n) in tables.items():
frac = 100 * n.peaks / n.pairs
rows.append({"window": w, "units": q, "peaks": int(n.peaks.sum()), "pairs": int(n.pairs.sum()),
"per_mouse_mean_%": frac.mean(), "per_mouse_sd_%": frac.std(ddof=0),
"per_mouse_sem_%": frac.std() / np.sqrt(len(frac)), "mice": len(n),
"pooled_%": 100 * n.peaks.sum() / n.pairs.sum()})
summ = pd.DataFrame(rows)
summ["all_within_tolerance"] = (((summ.peaks - 16119).abs() <= 1612) & ((summ.pairs - 2089890).abs() <= 208989)
& ((summ["per_mouse_mean_%"] - 0.96).abs() <= 0.13))
print(summ.round(4).to_string(index=False))
print("expected: 16119 +/- 1612 peaks, 2089890 +/- 208989 pairs, 0.96 +/- 0.13 % per mouse")
MATCH = [(w, q) for w, q, ok in zip(summ.window, summ.units, summ.all_within_tolerance) if ok]
if MATCH:
JUDGED = MATCH
else:
# no reading reproduces the parent claim: the dependent claims are judged under the reading whose pair
# count is within tolerance and whose peak count is nearest the paper's (an assumption)
near = summ[(summ.pairs - 2089890).abs() <= 208989]
near = near.loc[[(near.peaks - 16119).abs().idxmin()]]
JUDGED = list(zip(near.window, near.units))
J = JUDGED[0]
print("readings with all three values within tolerance:", MATCH)
print("readings the dependent claims are judged under:", JUDGED, "; the others are variants")
n0 = tables[J][1]
print(f"\nper mouse, {J}:")
print(n0.assign(percent=(100 * n0.peaks / n0.pairs).round(3)).to_string(index=False))
t0 = tables[J][0]
print("\noffset distribution of all sharp peaks (ms, one direction per pair):",
t0[t0.offset >= 0].offset.describe().round(2).to_dict())
loo = []
for s in n0.subject:
m = n0[n0.subject != s]
loo.append({"left_out": s, "peaks": int(m.peaks.sum()), "pairs": int(m.pairs.sum()), "per_mouse_%": (100 * m.peaks / m.pairs).mean()})
loo = pd.DataFrame(loo)
print("\n25 of the 26 mice, each one left out in turn: peaks", loo.peaks.min(), "to", loo.peaks.max(),
"; pairs", loo.pairs.min(), "to", loo.pairs.max(), "; per mouse %.3f to %.3f %%" % (loo["per_mouse_%"].min(), loo["per_mouse_%"].max()))
window units peaks pairs per_mouse_mean_% per_mouse_sd_% per_mouse_sem_% mice pooled_% all_within_tolerance
grating_2000 all 21962 2116259 1.1751 0.6907 0.1381 26 1.0378 False
grating_2000 ed1j 6460 544044 1.3264 0.7177 0.1435 26 1.1874 False
grating_3000 all 23474 2116259 1.2366 0.7088 0.1418 26 1.1092 False
grating_3000 ed1j 6885 544044 1.3997 0.7791 0.1558 26 1.2655 False
expected: 16119 +/- 1612 peaks, 2089890 +/- 208989 pairs, 0.96 +/- 0.13 % per mouse
readings with all three values within tolerance: []
readings the dependent claims are judged under: [('grating_2000', 'all')] ; the others are variants
per mouse, ('grating_2000', 'all'):
subject pairs peaks percent
744912845 63546 721 1.135
753795604 79800 1488 1.865
754477355 86320 1795 2.079
754488976 70876 664 0.937
756578432 52650 488 0.927
759674767 116886 368 0.315
759711149 90100 1917 2.128
760938794 131328 927 0.706
760946810 81003 844 1.042
760960650 46971 537 1.143
763236006 39060 528 1.352
763808588 18336 562 3.065
763884096 136503 853 0.625
769319621 47586 1184 2.488
774672354 133903 232 0.173
791857600 51040 698 1.368
795770029 60726 389 0.641
800249582 75855 407 0.537
800250054 59340 673 1.134
803390283 131328 955 0.727
811322604 83028 642 0.773
813701555 64261 554 0.862
817060743 94830 718 0.757
820866111 96580 2025 2.097
821469659 70500 506 0.718
827809879 133903 1287 0.961
offset distribution of all sharp peaks (ms, one direction per pair): {'count': 25630.0, 'mean': 2.64, 'std': 2.82, 'min': 0.0, '25%': 0.0, '50%': 2.0, '75%': 4.0, 'max': 10.0}
25 of the 26 mice, each one left out in turn: peaks 19937 to 21730 ; pairs 1979756 to 2097923 ; per mouse 1.100 to 1.215 %
Result: not reproduced. No reading gives all three values within tolerance. With 2 s windows
and all non-noise units the possible pairs match (2,116,259; paper 2,089,890 +/- 208,989), but
there are 21,962 sharp peaks (16,119 +/- 1,612) and the per mouse mean is 1.18% (0.96% +/- 0.13%).
The 3 s window gives more peaks (23,474), and the Extended Data Fig. 1j filters give 544,044
possible pairs, a quarter of the paper's count, so the paper's CCG units were not filtered that
way. With any one of the 26 mice left out, 19,937 to 21,730 peaks and 1.10% to 1.22% per mouse,
still outside. The s.e.m. across mice (0.14%) is near the paper's 0.13%, the SD 0.69%. The claims
that use the sharp peaks are judged under 2 s with all units, the reading whose pair count matches
and whose peak count is nearest the paper's. The paper's code subtracts one random jitter
resample, where this analysis subtracts the expected value as the Methods say; that variant was
not run, because the jitter function its code imports is not in the repositories the paper
links.
fig2d-v1-leads-lm-example¶
In an example mouse, the CCG peak offsets between V1 and LM were shifted positive against the V1-V1 distribution (median 3.9 ms), so V1 units on average lead LM units.
def per_mouse_test(t, src, tgt, ref_src, ref_tgt, ref="ref"):
# ref: "ref" each within-area pair once in a random orientation; "first" once in file order;
# "both" each within-area pair twice (offset and its negative), the round-1 version
out = []
for subj, g in t.groupby("subject"):
a = g[(g.src == src) & (g.tgt == tgt)].offset.values
bb = g[(g.src == ref_src) & (g.tgt == ref_tgt)]
b = (bb if ref == "both" else bb[bb[ref]]).offset.values
if len(a) and len(b):
out.append({"subject": subj, "n": len(a), "n_ref": len(b), "median": np.median(a),
"p": stats.ranksums(a, b)[1], "higher": np.median(a) > np.median(b)})
return pd.DataFrame(out)
def label(rq):
return "judged" if rq in JUDGED else "variant"
for rq in READINGS:
for ref in (("ref", "both") if rq == J else ("ref",)):
d = per_mouse_test(tables[rq][0], "VISp", "VISl", "VISp", "VISp", ref)
hit = d[(d["median"].sub(3.9).abs() <= 0.39) & (d.p < 0.05) & d.higher]
print(f"{rq} [{label(rq)}], V1-V1 reference {ref}: mice with both distributions {len(d)}; "
f"matching mice (median 3.51 to 4.29 ms, p < 0.05, higher): {len(hit)}; positive medians {int((d['median'] > 0).sum())}; "
f"p < 0.05 and higher {int(((d.p < 0.05) & d.higher).sum())}")
if rq == J and ref == "ref":
print(d.round(5).to_string(index=False))
if len(hit):
print("matching mice:")
print(hit.round(6).to_string(index=False))
('grating_2000', 'all') [judged], V1-V1 reference ref: mice with both distributions 15; matching mice (median 3.51 to 4.29 ms, p < 0.05, higher): 4; positive medians 14; p < 0.05 and higher 13
subject n n_ref median p higher
744912845 88 357 3.0 0.00000 True
754477355 273 635 4.0 0.00000 True
759711149 301 683 5.0 0.00000 True
760938794 75 306 3.0 0.00000 True
763808588 124 107 6.0 0.00000 True
774672354 1 106 7.0 0.09227 True
791857600 77 247 4.0 0.00000 True
800250054 47 70 5.0 0.00039 True
803390283 62 170 2.5 0.00037 True
811322604 16 53 4.5 0.00063 True
813701555 39 87 3.0 0.00370 True
817060743 41 442 4.0 0.00000 True
820866111 210 333 4.0 0.00000 True
821469659 1 23 -6.0 0.16994 False
827809879 141 386 2.0 0.00000 True
matching mice:
subject n n_ref median p higher
754477355 273 635 4.0 0.0 True
791857600 77 247 4.0 0.0 True
817060743 41 442 4.0 0.0 True
820866111 210 333 4.0 0.0 True
('grating_2000', 'all') [judged], V1-V1 reference both: mice with both distributions 15; matching mice (median 3.51 to 4.29 ms, p < 0.05, higher): 4; positive medians 14; p < 0.05 and higher 13
('grating_2000', 'ed1j') [variant], V1-V1 reference ref: mice with both distributions 13; matching mice (median 3.51 to 4.29 ms, p < 0.05, higher): 4; positive medians 13; p < 0.05 and higher 9
('grating_3000', 'all') [variant], V1-V1 reference ref: mice with both distributions 14; matching mice (median 3.51 to 4.29 ms, p < 0.05, higher): 7; positive medians 14; p < 0.05 and higher 13
('grating_3000', 'ed1j') [variant], V1-V1 reference ref: mice with both distributions 13; matching mice (median 3.51 to 4.29 ms, p < 0.05, higher): 4; positive medians 12; p < 0.05 and higher 10
Result: reproduced. Under the judged reading (2 s, all units), in 4 of the 15 mice with both distributions the median V1 to LM offset is 4.0 ms and the V1-LM offsets lie above the V1-V1 offsets (rank-sum p < 1e-6); 14 of the 15 medians are positive and 13 of 15 tests significant with V1-LM higher. The V1-V1 pairs counted twice (offset and its negative) also give 4 mice, and every other reading has at least 4.
fig2-neighbour-offsets-across-mice¶
Across 25 mice, the CCG peak offset distributions between neighbouring cortical areas differ from the within-area distributions for V1-LM, LM-RL, RL-AL and PM-AM, but not for AL-PM.
Per mouse two-sided rank-sum of the source-to-target offsets (lower area as source) against the lower area's within-area offsets, combined over mice with Fisher's method.
NEIGH = [("VISp", "VISl"), ("VISl", "VISrl"), ("VISrl", "VISal"), ("VISal", "VISpm"), ("VISpm", "VISam")]
expect = {("VISal", "VISpm"): False}
for rq in READINGS:
for ref in (("ref", "first", "both") if rq == J else ("ref",)):
rows = []
for a, b in NEIGH:
d = per_mouse_test(tables[rq][0], a, b, a, a, ref)
stat, p = stats.combine_pvalues(d.p.values, method="fisher") if len(d) else (np.nan, np.nan)
rows.append({"pair": f"{NAMES[a]}-{NAMES[b]}", "mice": len(d), "fisher_p": p, "significant": p < 0.05,
"expected": expect.get((a, b), True), "median_of_mouse_medians": d["median"].median() if len(d) else np.nan,
"units_pairs": int(d.n.sum()) if len(d) else 0, "ref_pairs": int(d.n_ref.sum()) if len(d) else 0})
t = pd.DataFrame(rows)
print(f"{rq} [{label(rq)}], within-area reference {ref}: agreeing pairs", int((t.significant == t.expected).sum()), "of 5")
print(t.to_string(index=False))
('grating_2000', 'all') [judged], within-area reference ref: agreeing pairs 4 of 5
pair mice fisher_p significant expected median_of_mouse_medians units_pairs ref_pairs
V1-LM 15 2.958540e-148 True True 4.00 1496 4005
LM-RL 13 3.524057e-11 True True -2.00 623 1443
RL-AL 12 5.818144e-15 True True 5.00 211 805
AL-PM 8 5.016023e-04 True False 2.25 48 768
PM-AM 14 2.689800e-06 True True 0.00 151 699
('grating_2000', 'all') [judged], within-area reference first: agreeing pairs 4 of 5
pair mice fisher_p significant expected median_of_mouse_medians units_pairs ref_pairs
V1-LM 15 9.764731e-166 True True 4.00 1496 4005
LM-RL 13 1.032661e-09 True True -2.00 623 1443
RL-AL 12 1.262686e-16 True True 5.00 211 805
AL-PM 8 1.278431e-04 True False 2.25 48 768
PM-AM 14 3.670675e-07 True True 0.00 151 699
('grating_2000', 'all') [judged], within-area reference both: agreeing pairs 4 of 5
pair mice fisher_p significant expected median_of_mouse_medians units_pairs ref_pairs
V1-LM 15 5.224761e-176 True True 4.00 1496 8010
LM-RL 13 1.585297e-13 True True -2.00 623 2886
RL-AL 12 2.528854e-17 True True 5.00 211 1610
AL-PM 8 1.144893e-04 True False 2.25 48 1536
PM-AM 14 3.796843e-07 True True 0.00 151 1398
('grating_2000', 'ed1j') [variant], within-area reference ref: agreeing pairs 4 of 5
pair mice fisher_p significant expected median_of_mouse_medians units_pairs ref_pairs
V1-LM 13 3.820296e-44 True True 4.0 397 1124
LM-RL 9 2.291793e-03 True True 0.0 175 252
RL-AL 7 6.977282e-09 True True 6.0 54 186
AL-PM 5 4.793574e-02 True False 1.0 22 152
PM-AM 11 1.223894e-03 True True 0.0 52 158
('grating_3000', 'all') [variant], within-area reference ref: agreeing pairs 4 of 5
pair mice fisher_p significant expected median_of_mouse_medians units_pairs ref_pairs
V1-LM 14 6.395974e-140 True True 4.0 1555 4181
LM-RL 12 1.406009e-14 True True -3.0 580 1481
RL-AL 12 2.446938e-17 True True 5.0 197 919
AL-PM 7 1.976632e-07 True False 5.0 49 726
PM-AM 14 6.786644e-08 True True 0.0 165 844
('grating_3000', 'ed1j') [variant], within-area reference ref: agreeing pairs 4 of 5
pair mice fisher_p significant expected median_of_mouse_medians units_pairs ref_pairs
V1-LM 13 2.636749e-37 True True 4.0 407 1160
LM-RL 8 2.716789e-03 True True -3.5 160 260
RL-AL 9 5.245877e-06 True True 4.0 41 242
AL-PM 4 2.902251e-04 True False 1.0 19 151
PM-AM 11 1.599251e-05 True True 0.0 56 180
Result: partially reproduced. Under the judged reading, V1-LM (Fisher p = 3e-148), LM-RL (3.5e-11), RL-AL (5.8e-15) and PM-AM (2.7e-06) differ from the within-area distributions, as stated, but AL-PM also differs (p = 5.0e-04, 8 mice): 4 of 5 pairs agree. The same 4 of 5 holds with the within-area pairs in file order (AL-PM p = 1.3e-04) or counted twice (1.1e-04), and under the three other readings (AL-PM p from 2.0e-07 to 0.048).
ed7c-v1-leads-am-follows¶
On average V1 units lead the units of the other cortical areas, while AM follows the other areas.
Per mouse and ordered pair, the median offset; averaged across mice.
def mean_median_matrix(t):
M = pd.DataFrame(np.nan, index=[NAMES[a] for a in AREA_ORDER], columns=[NAMES[a] for a in AREA_ORDER])
for sr in AREA_ORDER:
for tg in AREA_ORDER:
d = t[(t.src == sr) & (t.tgt == tg)]
if len(d):
M.loc[NAMES[sr], NAMES[tg]] = d.groupby("subject").offset.median().mean()
return M # rows: source, columns: target
for rq in READINGS:
M = mean_median_matrix(tables[rq][0])
v1 = M.loc["V1", ["LM", "RL", "AL", "PM", "AM"]]
am = M.loc[["V1", "LM", "RL", "AL", "PM"], "AM"]
print(f"{rq} [{label(rq)}]: V1 as source positive for {int((v1 > 0).sum())} of 5; AM as target positive for {int((am > 0).sum())} of 5; "
f"V1 row {v1.round(2).to_dict()}; AM column {am.round(2).to_dict()}")
if rq == J:
print("mean over mice of the median offset (ms), rows source, columns target:")
print(M.round(3).to_string())
('grating_2000', 'all') [judged]: V1 as source positive for 5 of 5; AM as target positive for 4 of 5; V1 row {'LM': 3.4, 'RL': 0.68, 'AL': 3.31, 'PM': 3.08, 'AM': 5.22}; AM column {'V1': 5.22, 'LM': 3.73, 'RL': 1.92, 'AL': 4.7, 'PM': -1.11}
mean over mice of the median offset (ms), rows source, columns target:
V1 LM RL AL PM AM
V1 0.000 3.400 0.676 3.306 3.077 5.225
LM -3.400 0.000 -1.962 2.409 1.889 3.733
RL -0.676 1.962 0.000 3.792 2.917 1.917
AL -3.306 -2.409 -3.792 0.000 1.375 4.700
PM -3.077 -1.889 -2.917 -1.375 0.000 -1.107
AM -5.225 -3.733 -1.917 -4.700 1.107 0.000
('grating_2000', 'ed1j') [variant]: V1 as source positive for 5 of 5; AM as target positive for 4 of 5; V1 row {'LM': 3.77, 'RL': 1.12, 'AL': 3.38, 'PM': 3.36, 'AM': 4.03}; AM column {'V1': 4.03, 'LM': 3.31, 'RL': 5.07, 'AL': 3.0, 'PM': -0.27}
('grating_3000', 'all') [variant]: V1 as source positive for 4 of 5; AM as target positive for 4 of 5; V1 row {'LM': 4.18, 'RL': -0.2, 'AL': 4.44, 'PM': 5.68, 'AM': 5.55}; AM column {'V1': 5.55, 'LM': 4.18, 'RL': 2.91, 'AL': 4.53, 'PM': -1.11}
('grating_3000', 'ed1j') [variant]: V1 as source positive for 5 of 5; AM as target positive for 4 of 5; V1 row {'LM': 3.42, 'RL': 0.95, 'AL': 3.61, 'PM': 4.11, 'AM': 4.79}; AM column {'V1': 4.79, 'LM': 4.0, 'RL': 5.69, 'AL': 4.23, 'PM': -0.64}
Result: partially reproduced. Under the judged reading the mean median offset is positive from V1 to each of the five other areas (0.68 to 5.22 ms) and into AM from V1, LM, RL and AL (1.92 to 5.22 ms); from PM to AM it is -1.11 ms. 9 of 10 pairs hold. The Extended Data Fig. 1j readings also give 9 of 10 (PM to AM -0.27 and -0.64 ms), and 3 s with all units 8 of 10 (V1 to RL -0.20 ms, PM to AM -1.11 ms), all within the criterion's partial range.
fig2f-directionality-vs-anatomy¶
The matrix of pairwise directionality scores between cortical areas is correlated with the matrix of anatomical hierarchy score differences (Pearson r = 0.74 over 21 pairs).
DS from offsets pooled over mice (zero lags in neither count), rows target and columns source;
anatomy target minus source; the 15 lower-triangle and 6 diagonal entries. Variants: DS averaged
across mice; zero lags counted positive; the paper's precomputed matrix FFscore_grating_10.npy.
fig2f = {}
for rq in READINGS:
t = tables[rq][0]
for mode in ("pooled", "mice"):
for zero in ("exclude", "positive"):
M = ds_matrix(t, mode, zero)
r, p, n = anatomy_corr(M)
fig2f[(rq, mode, zero)] = (r, p, n, M)
print(pd.DataFrame([{"window": k[0][0], "units": k[0][1], "DS": k[1], "zero lag": k[2], "r": v[0], "p": v[1], "n": v[2],
"": "judged" if (k[0] in JUDGED and k[1] == "pooled" and k[2] == "exclude") else "variant"}
for k, v in fig2f.items()]).round(5).to_string(index=False))
print("expected r = 0.74 +/- 0.074, p < 0.05")
M0 = fig2f[(J, "pooled", "exclude")][3]
print(f"\nDS matrix, {J}, pooled (rows target, columns source):")
print(pd.DataFrame(M0, index=[NAMES[a] for a in AREA_ORDER], columns=[NAMES[a] for a in AREA_ORDER]).round(3).to_string())
tJ = tables[J][0]
loo_f = pd.DataFrame([{"left_out": s, "r": anatomy_corr(ds_matrix(tJ[tJ.subject != s]))[0]} for s in sorted(tJ.subject.unique())])
loo_f["within_tolerance"] = (loo_f.r - 0.74).abs() <= 0.074
print(f"\n25 of the 26 mice, each one left out in turn: r {loo_f.r.min():.4f} to {loo_f.r.max():.4f}; "
f"within the tolerance (0.666 to 0.814) for {int(loo_f.within_tolerance.sum())} of {len(loo_f)}")
print(loo_f.sort_values("r").round(4).to_string(index=False))
ff = np.load(io.BytesIO(requests.get(f"{S.CODE_RAW}/data/processed_data/FFscore_grating_10.npy", timeout=60).content))[0]
code_order = ["VISp", "VISrl", "VISl", "VISal", "VISpm", "VISam"]
hs_code = np.array([S.HIERARCHY[a] for a in code_order])
Acode = hs_code[:, None] - hs_code[None, :]
il = np.tril_indices(6)
print("\nvariant, the paper's precomputed DS matrix (its own area order, V1 RL LM AL PM AM): r = %.4f, p = %.3g"
% stats.pearsonr(Acode[il], ff[il]))
window units DS zero lag r p n
grating_2000 all pooled exclude 0.82693 0.00000 21 judged
grating_2000 all pooled positive 0.75134 0.00009 21 variant
grating_2000 all mice exclude 0.77658 0.00003 21 variant
grating_2000 all mice positive 0.62136 0.00264 21 variant
grating_2000 ed1j pooled exclude 0.80997 0.00001 21 variant
grating_2000 ed1j pooled positive 0.76431 0.00005 21 variant
grating_2000 ed1j mice exclude 0.80579 0.00001 21 variant
grating_2000 ed1j mice positive 0.67416 0.00080 21 variant
grating_3000 all pooled exclude 0.84182 0.00000 21 variant
grating_3000 all pooled positive 0.78620 0.00002 21 variant
grating_3000 all mice exclude 0.79659 0.00002 21 variant
grating_3000 all mice positive 0.70292 0.00038 21 variant
grating_3000 ed1j pooled exclude 0.80037 0.00001 21 variant
grating_3000 ed1j pooled positive 0.76236 0.00006 21 variant
grating_3000 ed1j mice exclude 0.81732 0.00001 21 variant
grating_3000 ed1j mice positive 0.74309 0.00011 21 variant
expected r = 0.74 +/- 0.074, p < 0.05
DS matrix, ('grating_2000', 'all'), pooled (rows target, columns source):
V1 LM RL AL PM AM
V1 0.000 -0.650 -0.193 -0.570 -0.576 -0.693
LM 0.650 0.000 0.254 -0.064 -0.543 -0.611
RL 0.193 -0.254 0.000 -0.277 -0.867 -0.548
AL 0.570 0.064 0.277 0.000 -0.174 -0.602
PM 0.576 0.543 0.867 0.174 0.000 0.105
AM 0.693 0.611 0.548 0.602 -0.105 0.000
25 of the 26 mice, each one left out in turn: r 0.8024 to 0.8394; within the tolerance (0.666 to 0.814) for 3 of 26 left_out r within_tolerance 820866111 0.8024 True 800250054 0.8057 True 763884096 0.8112 True 827809879 0.8154 False 754477355 0.8223 False 800249582 0.8244 False 744912845 0.8247 False 817060743 0.8251 False 821469659 0.8253 False 763236006 0.8258 False 803390283 0.8260 False 813701555 0.8263 False 756578432 0.8264 False 760960650 0.8274 False 774672354 0.8277 False 759674767 0.8282 False 769319621 0.8283 False 795770029 0.8288 False 763808588 0.8294 False 759711149 0.8303 False 811322604 0.8306 False 760946810 0.8330 False 760938794 0.8348 False 754488976 0.8352 False 753795604 0.8359 False 791857600 0.8394 False variant, the paper's precomputed DS matrix (its own area order, V1 RL LM AL PM AM): r = 0.7351, p = 0.000147
Result: partially reproduced, resting on an assumption. Under the judged reading, with offsets pooled over mice, r = 0.827 (p < 1e-5) over the 21 entries: the correlation is there and significant, but r is above the tolerance (0.666 to 0.814). The criterion names no partial outcome; the correlation holds and only its size is outside the tolerance, so the claim holds in part. Leaving out any one mouse gives r = 0.802 to 0.839, and 3 of the 26 exclusions (mice 820866111, 800250054 and 763884096) put r within the tolerance (0.802 to 0.811). The paper used 25 of the 26 mice and the one it left out is not identified, so the outcome also rests on which 25 mice are used. Under the other readings r = 0.800 to 0.842 with pooled DS: the Extended Data Fig. 1j readings (0.810 with 2 s, 0.800 with 3 s) fall within the tolerance and 3 s with all units (0.842) does not. DS averaged across mice gives 0.777 to 0.817, and zero lags counted positive 0.621 to 0.786. The paper's precomputed DS matrix gives r = 0.735 with the same anatomy differences, which confirms their orientation.
ed7fg-absent-in-spontaneous¶
The functional hierarchy is absent during spontaneous activity (30 min grey screen).
Sharp peaks over 2 s or 3 s segments of the spontaneous block, the same units; DS averaged across mice; Pearson against the hierarchy differences over the 21 entries, at p < 0.05.
for w in ("spont_2000", "spont_3000"):
for q in ("all", "ed1j"):
t, n = peak_table(w, q)
M = ds_matrix(t, "mice")
r, p, k = anatomy_corr(M)
Mp = ds_matrix(t, "pooled")
rp, pp, kp = anatomy_corr(Mp)
print(f"{w}, {q} [{'judged' if q == J[1] else 'variant'}]: sharp peaks {int(n.peaks.sum())} of {int(n.pairs.sum())} pairs; "
f"per mouse {n.peaks.min()} to {n.peaks.max()} (median {n.peaks.median():.0f}); "
f"DS averaged across mice r = {r:.4f}, p = {p:.3g} (n = {k}); pooled r = {rp:.4f}, p = {pp:.3g}")
if q == J[1]:
lo = pd.DataFrame([dict(zip(("r", "p"), anatomy_corr(ds_matrix(t[t.subject != s], "mice"))[:2]))
for s in sorted(t.subject.unique())])
print(f" 25 of the 26 mice, each one left out in turn, DS averaged across mice: r {lo.r.min():.4f} to {lo.r.max():.4f}, "
f"p {lo.p.min():.3g} to {lo.p.max():.3g}; p >= 0.05 for {int((lo.p >= 0.05).sum())} of {len(lo)}")
spont_2000, all [judged]: sharp peaks 7326 of 2116259 pairs; per mouse 94 to 547 (median 294); DS averaged across mice r = 0.4965, p = 0.0221 (n = 21); pooled r = 0.4351, p = 0.0487
25 of the 26 mice, each one left out in turn, DS averaged across mice: r 0.4463 to 0.5549, p 0.00903 to 0.0426; p >= 0.05 for 0 of 26 spont_2000, ed1j [variant]: sharp peaks 2355 of 544044 pairs; per mouse 34 to 159 (median 88); DS averaged across mice r = 0.1424, p = 0.538 (n = 21); pooled r = 0.1914, p = 0.406
spont_3000, all [judged]: sharp peaks 7455 of 2116259 pairs; per mouse 96 to 571 (median 298); DS averaged across mice r = 0.4836, p = 0.0264 (n = 21); pooled r = 0.4448, p = 0.0433
25 of the 26 mice, each one left out in turn, DS averaged across mice: r 0.4185 to 0.5515, p 0.00956 to 0.059; p >= 0.05 for 2 of 26 spont_3000, ed1j [variant]: sharp peaks 2404 of 544044 pairs; per mouse 34 to 161 (median 92); DS averaged across mice r = 0.1696, p = 0.462 (n = 21); pooled r = 0.2330, p = 0.309
Result: not reproduced, resting on an assumption. With the judged unit set (all non-noise units) the DS matrix of the spontaneous block, averaged across mice, is correlated with the hierarchy differences: r = 0.497, p = 0.022 with 2 s segments and r = 0.484, p = 0.026 with 3 s segments (pooled DS r = 0.435 and 0.445, p = 0.049 and 0.043). The correlation is weaker than during gratings (0.827) and there are fewer sharp peaks (7,326 against 21,962 of the same 2,116,259 pairs). With the Extended Data Fig. 1j filters it is not significant (r = 0.142, p = 0.54 with 2 s; r = 0.170, p = 0.46 with 3 s), so the outcome depends on the unit set, which the paper does not state; the unit set whose pair count matches the paper's is adopted. With any one mouse left out, the 2 s correlation stays significant (p 0.009 to 0.043, 26 of 26) and the 3 s one does in 24 of 26 (p up to 0.059). The result is borderline: it is measured under a reading with 36% more sharp peaks during gratings than the paper reports (21,962 against 16,119).
ed8-total-hierarchy-score¶
The measured functional network has a total hierarchy score of 0.89.
Mean DS of each source over all targets including its own diagonal entry (np.mean(FF_score, axis=0)); THS the largest minus the smallest.
for rq in READINGS:
for mode in ("pooled", "mice"):
M = fig2f[(rq, mode, "exclude")][3]
col = np.nanmean(M, axis=0)
ths = np.nanmax(col) - np.nanmin(col)
Mn = M.copy()
np.fill_diagonal(Mn, np.nan)
coln = np.nanmean(Mn, axis=0)
print(f"{rq}, DS {mode} [{'judged' if (rq in JUDGED and mode == 'pooled') else 'variant'}]: THS = {ths:.4f} "
f"(without diagonal {np.nanmax(coln) - np.nanmin(coln):.4f}); "
f"source means {dict(zip([NAMES[a] for a in AREA_ORDER], np.round(col, 3)))}")
print("expected 0.89 +/- 0.089, between 0.15 and 1.54")
('grating_2000', 'all'), DS pooled [judged]: THS = 0.8385 (without diagonal 1.0062); source means {'V1': np.float64(0.447), 'LM': np.float64(0.053), 'RL': np.float64(0.292), 'AL': np.float64(-0.023), 'PM': np.float64(-0.377), 'AM': np.float64(-0.392)}
('grating_2000', 'all'), DS mice [variant]: THS = 0.7270 (without diagonal 0.8724); source means {'V1': np.float64(0.373), 'LM': np.float64(0.059), 'RL': np.float64(0.245), 'AL': np.float64(-0.087), 'PM': np.float64(-0.236), 'AM': np.float64(-0.354)}
('grating_2000', 'ed1j'), DS pooled [variant]: THS = 0.8433 (without diagonal 1.0120); source means {'V1': np.float64(0.467), 'LM': np.float64(0.051), 'RL': np.float64(0.295), 'AL': np.float64(-0.079), 'PM': np.float64(-0.376), 'AM': np.float64(-0.357)}
('grating_2000', 'ed1j'), DS mice [variant]: THS = 0.7979 (without diagonal 0.9575); source means {'V1': np.float64(0.432), 'LM': np.float64(0.031), 'RL': np.float64(0.312), 'AL': np.float64(-0.102), 'PM': np.float64(-0.306), 'AM': np.float64(-0.366)}
('grating_3000', 'all'), DS pooled [variant]: THS = 0.9218 (without diagonal 1.1062); source means {'V1': np.float64(0.47), 'LM': np.float64(0.081), 'RL': np.float64(0.349), 'AL': np.float64(-0.007), 'PM': np.float64(-0.441), 'AM': np.float64(-0.452)}
('grating_3000', 'all'), DS mice [variant]: THS = 0.8662 (without diagonal 1.0395); source means {'V1': np.float64(0.447), 'LM': np.float64(0.058), 'RL': np.float64(0.31), 'AL': np.float64(-0.048), 'PM': np.float64(-0.349), 'AM': np.float64(-0.419)}
('grating_3000', 'ed1j'), DS pooled [variant]: THS = 0.9320 (without diagonal 1.1184); source means {'V1': np.float64(0.48), 'LM': np.float64(0.058), 'RL': np.float64(0.369), 'AL': np.float64(-0.012), 'PM': np.float64(-0.444), 'AM': np.float64(-0.452)}
('grating_3000', 'ed1j'), DS mice [variant]: THS = 0.9794 (without diagonal 1.1753); source means {'V1': np.float64(0.476), 'LM': np.float64(0.052), 'RL': np.float64(0.296), 'AL': np.float64(-0.057), 'PM': np.float64(-0.263), 'AM': np.float64(-0.504)}
expected 0.89 +/- 0.089, between 0.15 and 1.54
Result: reproduced. Under the judged reading, with pooled offsets, THS = 0.839 (paper 0.89 +/- 0.089), between the one-to-all 0.15 and the feedforward 1.54. V1 has the highest mean DS as a source (0.447) and AM the lowest (-0.392). DS averaged across mice gives 0.727, outside the tolerance; the other pooled readings give 0.843 to 0.932, all within the tolerance, so the outcome does not rest on the adopted reading.
ed10ce-superficial-lead-deep¶
Within an area, superficial layers (2 to 4) are hierarchically lower than deep layers (5 and 6), with superficial layers tending to drive deep layers in the CCG analysis.
Layers from the paper's table (codes 2 and 4 superficial, 5 and 6 deep); pairs within one area, superficial unit as source. (a) mean over the six areas of each area's DS averaged across mice; (b) DS of the offsets pooled over the six areas, averaged across mice.
lay = layer_csv.set_index("ecephys_unit_id").cortical_layer
for rq in READINGS:
t = tables[rq][0].copy()
t["layer_src"] = t.src_id.map(lay)
t["layer_tgt"] = t.tgt_id.map(lay)
sd = t[(t.src == t.tgt) & t.layer_src.isin([2, 4]) & t.layer_tgt.isin([5, 6])]
per_area = {}
for a in AREA_ORDER:
g = sd[sd.src == a]
per_area[NAMES[a]] = np.nanmean([ds(x.offset.values) for _, x in g.groupby("subject")]) if len(g) else np.nan
pooled = np.nanmean([ds(x.offset.values) for _, x in sd.groupby("subject")])
print(f"{rq} [{label(rq)}]: pairs superficial-deep within an area {len(sd)}; (a) mean over areas {np.nanmean(list(per_area.values())):.4f}; "
f"(b) pooled over areas, mean over mice {pooled:.4f}; per area {dict((k, round(v, 3)) for k, v in per_area.items())}")
print("unit ids without a layer in the table:", int(pd.concat([r["units"] for r in res]).unit_id.map(lay).isna().sum()))
('grating_2000', 'all') [judged]: pairs superficial-deep within an area 4974; (a) mean over areas 0.3210; (b) pooled over areas, mean over mice 0.3177; per area {'V1': np.float64(0.286), 'LM': np.float64(0.316), 'RL': np.float64(0.107), 'AL': np.float64(0.199), 'PM': np.float64(0.461), 'AM': np.float64(0.557)}
('grating_2000', 'ed1j') [variant]: pairs superficial-deep within an area 1513; (a) mean over areas 0.2196; (b) pooled over areas, mean over mice 0.2636; per area {'V1': np.float64(0.283), 'LM': np.float64(0.326), 'RL': np.float64(0.165), 'AL': np.float64(0.236), 'PM': np.float64(0.238), 'AM': np.float64(0.069)}
('grating_3000', 'all') [variant]: pairs superficial-deep within an area 5327; (a) mean over areas 0.3196; (b) pooled over areas, mean over mice 0.2996; per area {'V1': np.float64(0.326), 'LM': np.float64(0.239), 'RL': np.float64(0.125), 'AL': np.float64(0.253), 'PM': np.float64(0.403), 'AM': np.float64(0.572)}
('grating_3000', 'ed1j') [variant]: pairs superficial-deep within an area 1615; (a) mean over areas 0.3199; (b) pooled over areas, mean over mice 0.2726; per area {'V1': np.float64(0.235), 'LM': np.float64(0.347), 'RL': np.float64(0.448), 'AL': np.float64(0.374), 'PM': np.float64(0.274), 'AM': np.float64(0.241)}
unit ids without a layer in the table: 49
Result: reproduced. Under the judged reading the superficial-to-deep DS within an area is (a) 0.321 and (b) 0.318, from 4,974 pairs, and positive in each of the six areas (0.11 to 0.56). Under every reading both are positive ((a) 0.220 to 0.321, (b) 0.264 to 0.318).
ed7e-distance-does-not-explain¶
The physical layout of the areas does not account for the correlation between directionality and hierarchy in terms of average distance.
Blocked. The claim needs the 3D CCF distance between units. In every session file of 000021 and
000022 the electrodes' y and z columns are identical (reproduce_siegle_2021.ipynb, data
check). The cell below reads the electrodes table of every probe (LFP) file of the two dandisets,
which carry their own, to see whether the left-right coordinate is there. A distance from the two
coordinates that are present (anterior-posterior x and dorsal-ventral y) is shown as a variant, not
judged: the six areas lie side by side across the cortex, so it leaves out the axis along which
most of them are separated.
import concurrent.futures as cf
def electrodes_yz(asset_id):
f = S.open_asset(asset_id)
e = f["general/extracellular_ephys/electrodes"]
y, z = e["y"][:], e["z"][:]
f.close()
same = (y == z) | (np.isnan(y) & np.isnan(z))
return len(y), int(same.sum())
probe_files = pd.concat([S.list_assets(d, v) for d, v in S.SOURCES.items()])
probe_files = probe_files[probe_files.path.str.contains("probe")]
with cf.ThreadPoolExecutor(16) as ex:
yz = pd.DataFrame(list(ex.map(electrodes_yz, probe_files.asset_id)), columns=["channels", "y_equals_z"])
print("probe files:", len(yz), dict(probe_files.dandiset.value_counts()), "; channels:", int(yz.channels.sum()),
"; channels with y equal to z (or both missing):", int(yz.y_equals_z.sum()),
"; files where y equals z on every channel:", int((yz.channels == yz.y_equals_z).sum()))
# variant, not judged: mean 2D (x, y) distance between simultaneously recorded units of each pair of areas
D = np.full((6, 6, len(res)), np.nan)
for s, r in enumerate(res):
u = r["units"]
for i, a in enumerate(AREA_ORDER):
for j, b in enumerate(AREA_ORDER[:i + 1]):
pa, pb = u[u.location == a][["x", "y"]].values, u[u.location == b][["x", "y"]].values
if len(pa) and len(pb):
D[i, j, s] = np.nanmean(np.hypot(pa[:, None, 0] - pb[None, :, 0], pa[:, None, 1] - pb[None, :, 1]))
Dm = np.nanmean(D, axis=2)
A = HSC[:, None] - HSC[None, :]
il = np.tril_indices(6)
ok = np.isfinite(Dm[il])
r, p = stats.pearsonr(A[il][ok], Dm[il][ok])
print(f"variant, mean 2D distance against the hierarchy difference: r = {r:.3f}, p = {p:.3g} ({int(ok.sum())} entries; "
"every lower-triangle difference is >= 0 in this area order, so signed and absolute are the same)")
print("mean 2D distance (um), rows and columns V1 LM RL AL PM AM:")
print(pd.DataFrame(Dm, index=[NAMES[a] for a in AREA_ORDER], columns=[NAMES[a] for a in AREA_ORDER]).round(0).to_string())
probe files: 325 {'000021': np.int64(182), '000022': np.int64(143)} ; channels: 27329 ; channels with y equal to z (or both missing): 27329 ; files where y equals z on every channel: 325
variant, mean 2D distance against the hierarchy difference: r = 0.581, p = 0.00574 (21 entries; every lower-triangle difference is >= 0 in this area order, so signed and absolute are the same)
mean 2D distance (um), rows and columns V1 LM RL AL PM AM:
V1 LM RL AL PM AM
V1 218.0 NaN NaN NaN NaN NaN
LM 1066.0 176.0 NaN NaN NaN NaN
RL 1420.0 1504.0 183.0 NaN NaN NaN
AL 1134.0 645.0 1168.0 173.0 NaN NaN
PM 984.0 1668.0 675.0 1346.0 194.0 NaN
AM 1301.0 1798.0 554.0 1462.0 373.0 191.0
All 325 probe files (182 of 000021, 143 of 000022, 27,329 channels) have y equal to z on
every channel, as the session files do, so the claim stays blocked. The 2D variant, from the two
coordinates present, gives r = 0.581 (p = 0.006) against the hierarchy difference; it is shown, not
judged, because it leaves out the left-right axis along which most of the six areas are
separated.
print(f"total run time {(time.time() - started) / 60:.1f} min")
total run time 18.0 min
Summary¶
| Claim | Outcome |
|---|---|
| fig2-sharp-peak-pairs | not reproduced: the pair count matches, the sharp peaks (21,962) and the per mouse fraction (1.18%) are above tolerance |
| fig2d-v1-leads-lm-example | reproduced |
| fig2-neighbour-offsets-across-mice | partially reproduced (AL-PM also differs), under every reading |
| ed7c-v1-leads-am-follows | partially reproduced (9 of 10), under every reading |
| fig2f-directionality-vs-anatomy | partially reproduced, r = 0.827 just above the tolerance; rests on the adopted reading and on which 25 mice are used (3 of 26 single-mouse exclusions give r within tolerance) |
| ed7e-distance-does-not-explain | blocked: the left-right CCF coordinate is missing from all session and probe files |
| ed7fg-absent-in-spontaneous | not reproduced with all units (p = 0.022 and 0.026), a borderline result; with the Extended Data Fig. 1j filters the correlation is absent; rests on the adopted unit set |
| ed8-total-hierarchy-score | reproduced, THS 0.839 |
| ed10ce-superficial-lead-deep | reproduced |
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 quality filter of the paper's CCG units and the CCG window are not stated. All non-noise units are adopted because they give the paper's number of possible pairs, and the 2 s window because it gives fewer excess peaks than 3 s; every reading is shown. The jitter correction uses the exact expected value under jitter of both trains within 25 ms windows, as the Methods describe it; the paper's code subtracts one random jitter resample, which adds noise to the flanks and so could give fewer sharp peaks, but its jitter function is not public and that variant was not run. The dandiset has 26 Functional Connectivity mice, against the paper's 25. No single mouse left out changes the outcome of fig2-sharp-peak-pairs, but leaving out any of 3 of the 26 mice puts the fig2f-directionality-vs-anatomy correlation within tolerance. The spontaneous block is cut into consecutive 2 s or 3 s segments and uses the units selected during gratings, choices the paper does not describe.
What was found about the data. The electrodes tables of all 58 session files and all 325 probe
files of 000021 and 000022 have identical y and z columns, so the left-right CCF coordinate is
missing.