Reproducing Xue, Buccino et al. (2022) from dandiset 000223¶
Paper. Xue X, Buccino AP, Saseendran Kumar S, Hierlemann A, Bartram J (2022). Inferring monosynaptic connections from paired dendritic spine Ca2+ imaging and large-scale recording of extracellular spiking. Journal of Neural Engineering 19(4), 046044. doi:10.1088/1741-2552/ac8765. The claims file was written from the accepted author manuscript in PubMed Central (PMC7613561) and its supplement; this attempt read the same text. The paper's code is starquakes/mea_spineca_mapping, read at commit 99c2a478393341a24c58d6ef6ab5db5b0b528d03.
Dandiset. DANDI:000223, version 0.260528.0906. The paper cites the dandiset without a version. It has two published versions, 0.220823.0826 (August 2022, the month the paper appeared) and 0.260528.0906 (May 2026). The data check below shows that both list the same 20 assets with the same asset and blob identifiers, so their files are identical, and the most recent published version is pinned, as the protocol asks when the paper cites none.
Method. The pipeline of the claims file: the code's dF/F, demixing, convolution, network burst
detection and surrogates, with every parameter's source given in the helper spineca.py. Details
the claims file leaves open and the paper settles are cited where they are used. The claims file's
rule for the open burst parameters is followed: for each recording, the pair (p_thresh, hyperactive
strategy) whose best-R is nearest to Table 1 is carried forward to every claim.
| claim | tested here |
|---|---|
| data-scale | tested |
| table1-best-r | tested |
| fig1e-best-r-gap | tested |
| table1-optimal-thresholds | tested |
| mean-max-fscore | tested |
| acceptance-9-of-20 | tested |
| best-r-accepted-vs-rejected | tested |
| best-r-not-predictive | tested |
| fig5a-rejection-scenarios | tested |
| fig5a-threshold-order | tested |
| table1-network-measures | tested |
| fig5c-spike-contrast-vs-fscore | tested |
| fig5d-out-in-vs-fscore | tested |
In [1]:
import io
import itertools
import time
import warnings
from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor
import elephant
import matplotlib.pyplot as plt
import neo
import numpy as np
import pandas as pd
import quantities as pq
import requests
from elephant.spike_train_synchrony import spike_contrast
from scipy import stats
from scipy.optimize import linear_sum_assignment
import spineca as sc
warnings.filterwarnings("ignore")
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
started = time.time()
assets = sc.list_assets()
with ThreadPoolExecutor(10) as ex:
sessions = list(ex.map(lambda a: sc.load_session(a["asset_id"]), assets))
for a, s in zip(assets, sessions):
s["path"], s["blob"], s["size"] = a["path"], a["blob"], a["size"]
table = pd.DataFrame(sc.TABLE1, columns=sc.TABLE1_COLUMNS)
table.index = np.arange(1, len(table) + 1)
print(f"{len(sessions)} files of dandiset {sc.DANDISET_ID} version {sc.DANDISET_VERSION} loaded "
f"in {time.time() - started:.0f} s")
print("elephant", elephant.__version__, "(the paper used 0.10.0)")
20 files of dandiset 000223 version 0.260528.0906 loaded in 31 s elephant 1.2.1 (the paper used 0.10.0)
Data check¶
In [2]:
earlier = sc.list_assets("0.220823.0826")
same = {(a["asset_id"], a["blob"], a["path"]) for a in earlier} == {(a["asset_id"], a["blob"], a["path"]) for a in assets}
print(f"Version 0.220823.0826 has {len(earlier)} assets, {sc.DANDISET_VERSION} has {len(assets)}; "
f"same asset ids, blobs and paths: {same}")
info = sc.version_info()
print("Archive metadata: species", [x.get("name") for x in info["metadata"]["assetsSummary"].get("species", [])],
"| subjects", info["metadata"]["assetsSummary"].get("numberOfSubjects"),
"| files", info["metadata"]["assetsSummary"].get("numberOfFiles"))
files = pd.DataFrame([{
"file": s["path"].split("/")[1].replace("_ecephys+ophys.nwb", ""), "subject": s["subject_id"],
"start": s["session_start"][:19], "GB": round(s["size"] / 1e9, 2), "units": len(s["trains"]),
"spikes": sum(len(t) for t in s["trains"]), "channels": s["n_channels"],
"MEA_s": round(s["mea_duration"], 2), "frames": s["n_frames"],
"img_Hz_meta": round(s["imaging_rate"], 2), "img_Hz_ts": round(1 / np.median(np.diff(s["frame_times"])), 2),
"first_frame_s": round(s["frame_times"][0], 3), "last_frame_s": round(s["frame_times"][-1], 2),
"image": "x".join(map(str, s["image_shape"][1:])), "rois": s["roi_traces"].shape[1],
"roi_pixels": [int(m.sum()) for m in s["masks"]], "blob": s["blob"][:8]} for s in sessions])
print(files.sort_values(["subject", "start"]).to_string(index=False))
dups = files[files.duplicated("blob", keep=False)]
print("\nFiles that share a blob (the same bytes):")
print(dups[["file", "subject", "start", "units", "frames", "blob"]].to_string(index=False))
# keep one file per recording: the copy whose name carries the session
unique = [i for i, s in enumerate(sessions)
if not (files.blob.duplicated(keep=False)[i] and "_ses-" not in s["path"])]
print(f"\nRecordings counted once: {len(unique)}; subjects: {files.loc[unique, 'subject'].nunique()}")
# other pairs with the same unit, spike and frame counts
sig = files.loc[unique].groupby(["units", "spikes", "frames"]).size()
print("Other recordings sharing unit, spike and frame counts:", int((sig > 1).sum()))
Version 0.220823.0826 has 20 assets, 0.260528.0906 has 20; same asset ids, blobs and paths: True
Archive metadata: species ['Rattus norvegicus - Norway rat'] | subjects 10 | files 20
file subject start GB units spikes channels MEA_s frames img_Hz_meta img_Hz_ts first_frame_s last_frame_s image rois roi_pixels blob
sub-2282_ses-20190914T145458 2282 2019-09-14T14:54:58 4.97 159 139365 964 262.74 1800 6.85 6.85 0.038 262.70 512x512 2 [2, 9] aff5f64d
sub-2282_ses-20190914T151231 2282 2019-09-14T15:12:31 4.22 185 152473 964 262.72 1800 6.85 6.85 0.038 262.68 512x512 2 [4, 11] 3cc7fcd7
sub-2282_ses-20190914T162358 2282 2019-09-14T16:23:58 5.77 78 68372 996 246.43 1800 7.30 7.30 0.023 246.41 1024x1024 2 [10, 32] 57b07b7f
sub-2282_ses-20190914T164312 2282 2019-09-14T16:43:12 4.28 132 179532 955 187.07 1364 7.29 7.29 0.023 187.05 1024x1024 2 [9, 22] 5ede395d
sub-2282_ses-20190914T170145 2282 2019-09-14T17:01:45 3.79 66 42088 963 163.92 1200 7.32 7.31 0.023 163.89 1024x1024 2 [5, 28] 2796a7cd
sub-2610_ses-20201005T153308 2610 2020-10-05T15:33:08 4.41 87 41612 931 184.10 1800 9.78 9.78 0.050 184.05 1024x1024 2 [16, 25] 64257919
sub-2615_ses-20201007T162625 2615 2020-10-07T16:26:25 3.69 50 17129 939 124.08 1791 14.43 14.43 0.034 124.04 1024x1024 2 [5, 44] 6301e682
sub-2748_ses-20200928T173432 2748 2020-09-28T17:34:32 2.77 59 15244 965 105.55 1032 9.78 9.78 0.050 105.50 1024x1024 2 [12, 19] 9bd574f6
sub-2763_ses-20200828T154656 2763 2020-08-28T15:46:56 3.21 55 20899 976 122.73 1200 9.78 9.78 0.050 122.68 1024x1024 2 [4, 30] d1fb3728
sub-2763_ses-20200828T161323 2763 2020-08-28T16:13:23 2.99 50 13979 983 122.73 1200 9.78 9.78 0.050 122.68 1024x1024 2 [8, 42] 05846ee6
sub-2763_ses-20200828T163435 2763 2020-08-28T16:34:35 3.13 66 15880 978 122.73 1200 9.78 9.78 0.050 122.68 1024x1024 2 [13, 60] 8d7aaaed
sub-2765_ses-20200814T153629 2765 2020-08-14T15:36:29 3.27 59 14447 955 122.73 1200 9.78 9.78 0.050 122.68 1024x1024 2 [7, 21] 796c27a6
sub-3112_ses-20210227T185041 3112 2021-02-27T18:50:41 4.15 133 73691 964 184.10 1800 9.78 9.78 0.101 184.00 1024x1024 2 [6, 17] 6c0250a1
sub-3112 3112 2021-02-27T19:02:33 4.18 133 65536 964 184.10 1800 9.78 9.78 0.050 184.05 1024x1024 2 [6, 6] c5f793ff
sub-3112_ses-20210227T190233 3112 2021-02-27T19:02:33 4.18 133 65536 964 184.10 1800 9.78 9.78 0.050 184.05 1024x1024 2 [6, 6] c5f793ff
sub-3278_ses-20210614T151114 3278 2021-06-14T15:11:14 4.94 101 40777 993 184.10 1800 9.78 9.78 0.050 184.05 1024x1024 2 [10, 24] 654828e2
sub-3443_ses-20210814T190234 3443 2021-08-14T19:02:34 4.92 128 27542 1019 184.10 1800 9.78 9.78 0.050 184.05 1024x1024 2 [13, 44] de95efaf
sub-3443_ses-20210814T193211 3443 2021-08-14T19:32:11 5.52 114 16149 1019 200.98 1965 9.78 9.78 0.050 200.93 1024x1024 2 [7, 21] fcd37156
sub-3445_ses-20210814T223214 3445 2021-08-14T22:32:14 4.95 99 8492 1019 184.10 1800 9.78 9.78 0.050 184.05 1024x1024 2 [6, 8] 92bc19e2
sub-3445_ses-20210814T223641 3445 2021-08-14T22:36:41 4.94 179 24888 1019 184.10 1800 9.78 9.78 0.050 184.05 1024x1024 4 [9, 12, 8, 12] 9ca64bdd
Files that share a blob (the same bytes):
file subject start units frames blob
sub-3112 3112 2021-02-27T19:02:33 133 1800 c5f793ff
sub-3112_ses-20210227T190233 3112 2021-02-27T19:02:33 133 1800 c5f793ff
Recordings counted once: 19; subjects: 10
Other recordings sharing unit, spike and frame counts: 0
In [3]:
rows = []
for i in unique:
s = sessions[i]
rates = sc.unit_rates(s)
sorted_ok = all(np.all(np.diff(t) >= 0) for t in s["trains"])
allsp = np.concatenate(s["trains"])
rows.append({"file": files.file[i], "units": len(s["trains"]), "units_0_spikes": int(sum(len(t) == 0 for t in s["trains"])),
"rate_min_Hz": rates.min().round(3), "rate_median_Hz": np.median(rates).round(2),
"rate_max_Hz": rates.max().round(1), "spike_min_s": allsp.min().round(3),
"spike_max_s": allsp.max().round(2), "MEA_s": round(s["mea_duration"], 2),
"sorted": sorted_ok, "MEA_Hz": s["mea_rate"], "MEA_start": s["mea_start"],
"ROI_Hz": round(s["roi_rate"], 3), "indicator": s["indicator"]})
units_check = pd.DataFrame(rows)
print(units_check.to_string(index=False))
print("\nUnit table columns:", sorted({c for i in unique for c in sessions[i]["unit_cols"]}))
print("Imaging rate range over recordings (timestamps):", files.loc[unique, "img_Hz_ts"].min(), "to",
files.loc[unique, "img_Hz_ts"].max(), "Hz; Methods: 7 to 15 Hz")
print("Recording length:", files.loc[unique, "MEA_s"].min(), "to", files.loc[unique, "MEA_s"].max(),
"s; Methods: typically 2 to 3 min")
print("Last frame time minus MEA duration (s):",
np.round([sessions[i]["frame_times"][-1] - sessions[i]["mea_duration"] for i in unique], 3))
print("Spine ROI not larger than its shaft ROI in every pair:",
all(int(sessions[i]["masks"][a].sum()) <= int(sessions[i]["masks"][b].sum())
for i in unique for _, a, b in sc.spines(sessions[i])))
dd = []
for i in unique:
s = sessions[i]
for c in range(s["roi_traces"].shape[1]):
d = sc.dff(s["roi_traces"][:, c])
dd.append(np.nanmax(np.abs(d - s["dff_stored"][:, c])))
print(f"Stored Dff against dF/F recomputed from the ROI traces (25th percentile, 100 frames): "
f"largest absolute difference per ROI {np.min(dd):.4f} to {np.max(dd):.4f}")
file units units_0_spikes rate_min_Hz rate_median_Hz rate_max_Hz spike_min_s spike_max_s MEA_s sorted MEA_Hz MEA_start ROI_Hz indicator sub-3112_ses-20210227T190233 133 0 0.212 2.00 14.4 0.006 184.09 184.10 True 20000.0 0.0 9.777 unknown sub-2748_ses-20200928T173432 59 1 0.000 1.87 9.3 0.002 105.53 105.55 True 20000.0 0.0 9.777 unknown sub-2763_ses-20200828T161323 50 0 0.473 1.57 10.0 0.014 122.52 122.73 True 20000.0 0.0 9.777 unknown sub-2763_ses-20200828T154656 55 0 0.456 1.71 14.4 0.012 122.73 122.73 True 20000.0 0.0 9.777 unknown sub-2763_ses-20200828T163435 66 0 0.391 1.32 12.2 0.146 122.59 122.73 True 20000.0 0.0 9.777 unknown sub-2765_ses-20200814T153629 59 0 0.073 1.27 10.2 0.007 122.69 122.73 True 20000.0 0.0 9.777 unknown sub-2615_ses-20201007T162625 50 0 0.185 1.10 25.2 0.008 124.07 124.08 True 20000.0 0.0 14.435 unknown sub-3112_ses-20210227T185041 133 0 0.473 2.26 15.2 0.016 184.07 184.10 True 20000.0 0.0 9.777 unknown sub-2610_ses-20201005T153308 87 0 0.223 1.67 16.1 0.011 184.09 184.10 True 20000.0 0.0 9.777 unknown sub-3443_ses-20210814T190234 128 0 0.185 0.98 4.1 0.040 184.10 184.10 True 20000.0 0.0 9.777 unknown sub-3445_ses-20210814T223641 179 0 0.027 0.50 4.3 0.000 184.09 184.10 True 20000.0 0.0 9.777 unknown sub-3278_ses-20210614T151114 101 1 0.000 1.74 8.2 0.020 184.10 184.10 True 20000.0 0.0 9.777 unknown sub-3445_ses-20210814T223214 99 0 0.005 0.34 1.5 0.050 184.04 184.10 True 20000.0 0.0 9.777 unknown sub-3443_ses-20210814T193211 114 0 0.075 0.56 3.4 0.368 200.95 200.98 True 20000.0 0.0 9.777 unknown sub-2282_ses-20190914T170145 66 0 0.634 3.41 16.3 0.002 163.91 163.92 True 20000.0 0.0 7.321 unknown sub-2282_ses-20190914T151231 185 0 0.346 2.77 11.2 0.010 262.72 262.72 True 20000.0 0.0 6.851 unknown sub-2282_ses-20190914T164312 132 0 0.144 5.72 62.5 0.000 187.07 187.07 True 20000.0 0.0 7.291 unknown sub-2282_ses-20190914T145458 159 0 0.765 3.09 8.9 0.004 262.73 262.74 True 20000.0 0.0 6.851 unknown sub-2282_ses-20190914T162358 78 0 0.101 3.17 17.4 0.037 246.42 246.43 True 20000.0 0.0 7.304 unknown Unit table columns: ['spike_times'] Imaging rate range over recordings (timestamps): 6.85 to 14.43 Hz; Methods: 7 to 15 Hz Recording length: 105.55 to 262.74 s; Methods: typically 2 to 3 min Last frame time minus MEA duration (s): [-0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.034 -0.101 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.023 -0.038 -0.023 -0.038 -0.023] Spine ROI not larger than its shaft ROI in every pair: True Stored Dff against dF/F recomputed from the ROI traces (25th percentile, 100 frames): largest absolute difference per ROI 0.0043 to 0.2348
The two copies of the sub-3112 recording of 2021-02-27 19:02:33 are the same blob, so the
dandiset holds 19 distinct recordings from 10 subjects, as the paper reports 19 recordings from
10 cultures; each subject id is a culture. The units table holds only spike times: no
electrode, waveform or quality columns. ROIs come in pairs, the smaller first, which fits the
code's order (spine, then its shaft; notebook 1). The stored Dff series is not exactly the
dF/F of the ROI traces with the paper's settings, so it is not used; dF/F is recomputed.
In [4]:
# Clock check: the shaft trace carries the back-propagating action potentials of the imaged cell,
# which fire with the network. Its correlation with the population rate (sum of all units'
# convolved rates at the frame times) should peak near zero lag if the two streams share a clock.
lags = np.arange(-30, 31)
clock = []
for i in unique:
s = sessions[i]
s["rates"], s["U"] = sc.convolved_traces(s)
pop = s["rates"].sum(axis=0)
for k, a, b in sc.spines(s):
sh = sc.dff(s["roi_traces"][:, b])
r = [np.corrcoef(sh[max(L, 0):len(sh) + min(L, 0)], pop[max(-L, 0):len(pop) - max(L, 0)])[0, 1] for L in lags]
far = [np.corrcoef(np.roll(sh, sft), pop)[0, 1] for sft in (300, 450, 600)]
clock.append({"file": files.file[i], "spine": k, "peak_lag_frames": int(lags[np.argmax(r)]),
"peak_lag_s": round(lags[np.argmax(r)] / s["imaging_rate"], 2), "R_peak": round(max(r), 3),
"R_lag0": round(r[30], 3), "R_shift_30_to_60s_max": round(max(far), 3)})
clock = pd.DataFrame(clock)
print(clock.to_string(index=False))
print("Peak lag within 2 frames of zero:", int((clock.peak_lag_frames.abs() <= 2).sum()), "of", len(clock))
file spine peak_lag_frames peak_lag_s R_peak R_lag0 R_shift_30_to_60s_max sub-3112_ses-20210227T190233 0 0 0.0 0.904 0.904 0.023 sub-2748_ses-20200928T173432 0 0 0.0 0.697 0.697 -0.005 sub-2763_ses-20200828T161323 0 0 0.0 0.560 0.560 0.023 sub-2763_ses-20200828T154656 0 0 0.0 0.772 0.772 0.048 sub-2763_ses-20200828T163435 0 0 0.0 0.655 0.655 0.063 sub-2765_ses-20200814T153629 0 0 0.0 0.880 0.880 0.054 sub-2615_ses-20201007T162625 0 0 0.0 0.916 0.916 0.132 sub-3112_ses-20210227T185041 0 0 0.0 0.882 0.882 0.045 sub-2610_ses-20201005T153308 0 0 0.0 0.436 0.436 0.007 sub-3443_ses-20210814T190234 0 0 0.0 0.606 0.606 0.060 sub-3445_ses-20210814T223641 0 0 0.0 0.698 0.698 0.066 sub-3445_ses-20210814T223641 1 0 0.0 0.679 0.679 0.033 sub-3278_ses-20210614T151114 0 0 0.0 0.387 0.387 0.032 sub-3445_ses-20210814T223214 0 1 0.1 0.558 0.541 0.037 sub-3443_ses-20210814T193211 0 0 0.0 0.892 0.892 0.058 sub-2282_ses-20190914T170145 0 0 0.0 0.286 0.286 0.039 sub-2282_ses-20190914T151231 0 0 0.0 0.130 0.130 0.063 sub-2282_ses-20190914T164312 0 0 0.0 0.262 0.262 0.023 sub-2282_ses-20190914T145458 0 0 0.0 0.553 0.553 0.044 sub-2282_ses-20190914T162358 0 0 0.0 0.829 0.829 0.070 Peak lag within 2 frames of zero: 20 of 20
In [5]:
# The helper computes the convolution at the frame times directly. Compare with elephant's
# instantaneous_rate at 20 kHz, sampled with np.searchsorted, as the code does, for three units.
from elephant import kernels
from elephant.statistics import instantaneous_rate
s = sessions[unique[0]]
mea_times = np.arange(s["n_mea_samples"]) / sc.MEA_FS
ds_idx = np.searchsorted(mea_times, s["frame_times"])
print("Sample index of each frame equal to the helper's:",
np.array_equal(ds_idx, np.round(sc.frame_sample_times(s) * sc.MEA_FS).astype(int)))
kern = kernels.ExponentialKernel(sigma=0.5 * pq.s)
for u in [0, len(s["trains"]) // 2, len(s["trains"]) - 1]:
st = neo.SpikeTrain(s["trains"][u] * pq.s, t_stop=s["mea_duration"] * pq.s, sampling_rate=sc.MEA_FS * pq.Hz)
ifr = instantaneous_rate(st, sampling_period=st.sampling_period, center_kernel=False, kernel=kern)
e = ifr.magnitude.squeeze()[ds_idx]
print(f"{files.file[unique[0]]} unit {u}: {len(st)} spikes, max |elephant - helper| = "
f"{np.max(np.abs(e - s['rates'][u])):.3f} Hz (peak {e.max():.2f} Hz), R = {np.corrcoef(e, s['rates'][u])[0, 1]:.6f}")
del mea_times
Sample index of each frame equal to the helper's: True
sub-3112_ses-20210227T190233 unit 0: 1076 spikes, max |elephant - helper| = 0.335 Hz (peak 45.43 Hz), R = 0.999978
sub-3112_ses-20210227T190233 unit 66: 437 spikes, max |elephant - helper| = 0.242 Hz (peak 36.75 Hz), R = 0.999977
sub-3112_ses-20210227T190233 unit 132: 100 spikes, max |elephant - helper| = 0.134 Hz (peak 19.42 Hz), R = 0.999978
Files and Table 1 rows¶
The paper gives no map from Table 1 rows to files; the claims file names culture and unit count as the keys. Subjects are matched to cultures by their sorted unit counts, and files to rows by unit count within a culture. Where two assignments fit equally, the one whose spike contrast (all units, whole recording, elephant defaults) is nearer to Table 1 is taken, as an assumption, and the other is shown. In the file with two spines, the first ROI pair is taken as Table 1's first spine (row 9), also an assumption.
In [6]:
# Table 1 recordings: one per row, except rows 9 and 10, two spines of one recording
rec_key = pd.Series(np.where(table.index == 10, 9, table.index), index=table.index)
rec_rows = pd.DataFrame([{"culture": table.culture[k], "units": table.units[k], "rows": list(g.index)}
for k, g in table.groupby(rec_key)])
subj = files.loc[unique].groupby("subject")
cult = rec_rows.groupby("culture")
S_ids, C_ids = list(subj.groups), list(cult.groups)
def group_cost(su, cu):
a = sorted(files.loc[subj.groups[su], "units"])
b = sorted(rec_rows.loc[cult.groups[cu], "units"])
return sum(abs(x - y) for x, y in zip(a, b)) if len(a) == len(b) else 10 ** 6
cost = np.array([[group_cost(su, cu) for cu in C_ids] for su in S_ids])
ri, ci = linear_sum_assignment(cost)
subject_culture = {S_ids[r]: C_ids[c] for r, c in zip(ri, ci)}
ties = [(S_ids[a], S_ids[b]) for a, b in itertools.combinations(ri, 2)
if cost[a, ci[a]] + cost[b, ci[b]] == cost[a, ci[b]] + cost[b, ci[a]]]
def spike_contrast_of(s, t0=None, t1=None, units=None):
t0 = 0.0 if t0 is None else t0
t1 = s["mea_duration"] if t1 is None else t1
idx = range(len(s["trains"])) if units is None else units
sts = [neo.SpikeTrain(s["trains"][u][(s["trains"][u] >= t0) & (s["trains"][u] <= t1)] * pq.s,
t_start=t0 * pq.s, t_stop=t1 * pq.s) for u in idx]
return spike_contrast(sts)
sc_all = {i: spike_contrast_of(sessions[i]) for i in unique}
def file_rows(subject_culture):
"""Files to Table 1 recordings within each culture, by unit count; returns {file index: rows}
and the list of tied pairs."""
out, tied = {}, []
for su, cu in subject_culture.items():
fi = list(subj.groups[su])
rr = list(cult.groups[cu])
c = np.array([[abs(files.units[f] - rec_rows.units[r]) for r in rr] for f in fi])
a, b = linear_sum_assignment(c)
for x, y in zip(a, b):
out[fi[x]] = rec_rows.rows[rr[y]]
for x, y in itertools.combinations(range(len(a)), 2):
if c[a[x], b[x]] + c[a[y], b[y]] == c[a[x], b[y]] + c[a[y], b[x]]:
tied.append((fi[a[x]], fi[a[y]]))
return out, tied
mapping, tied_files = file_rows(subject_culture)
def sc_distance(m):
return sum(abs(sc_all[f] - table.spike_contrast[r[0]]) for f, r in m.items())
print("Subjects tied between two cultures:", ties)
print("Files tied between two rows:", [(files.file[a], files.file[b]) for a, b in tied_files])
alternatives = []
for su1, su2 in ties: # subjects tied for two cultures: keep the order nearer in spike contrast
alt_sc = dict(subject_culture)
alt_sc[su1], alt_sc[su2] = alt_sc[su2], alt_sc[su1]
alt = file_rows(alt_sc)[0]
print(f"Swap subjects {su1}/{su2}: summed |spike contrast - Table 1| {sc_distance(mapping):.3f} as assigned, "
f"{sc_distance(alt):.3f} swapped")
alternatives.append(("swap subjects " + su1 + "/" + su2, (su1, su2)))
if sc_distance(alt) < sc_distance(mapping):
subject_culture, mapping = alt_sc, alt
for f1, f2 in tied_files: # files tied for two rows within a culture: the same rule
alt = dict(mapping)
alt[f1], alt[f2] = alt[f2], alt[f1]
print(f"Swap files {files.file[f1]}/{files.file[f2]}: {sc_distance(mapping):.3f} as assigned, "
f"{sc_distance(alt):.3f} swapped")
alternatives.append(("swap files " + files.file[f1] + "/" + files.file[f2], (f1, f2)))
if sc_distance(alt) < sc_distance(mapping):
mapping = alt
link = pd.DataFrame([{"row": r, "culture": table.culture[r], "file": files.file[f], "subject": files.subject[f],
"units_file": files.units[f], "units_table": table.units[r],
"spike_contrast_file": round(sc_all[f], 3), "spike_contrast_table": table.spike_contrast[r]}
for f, rs in mapping.items() for r in rs]).sort_values("row")
print(link.to_string(index=False))
# one entry per tested spine (Table 1 row)
spine_rows = {}
for f, rs in mapping.items():
for (k, a, b), r in zip(sc.spines(sessions[f]), rs):
spine_rows[r] = {"file": f, "k": k, "spine_roi": a, "shaft_roi": b}
print("Spines:", len(spine_rows), "| ROI pairs in the recordings:", sum(len(sc.spines(sessions[i])) for i in unique))
Subjects tied between two cultures: [('2748', '2765')]
Files tied between two rows: [('sub-3112_ses-20210227T190233', 'sub-3112_ses-20210227T185041')]
Swap subjects 2748/2765: summed |spike contrast - Table 1| 0.389 as assigned, 0.189 swapped
Swap files sub-3112_ses-20210227T190233/sub-3112_ses-20210227T185041: 0.189 as assigned, 0.089 swapped
row culture file subject units_file units_table spike_contrast_file spike_contrast_table
1 1 sub-2282_ses-20190914T162358 2282 78 78 0.378 0.38
2 1 sub-2282_ses-20190914T164312 2282 132 132 0.172 0.17
3 1 sub-2282_ses-20190914T170145 2282 66 66 0.170 0.17
4 1 sub-2282_ses-20190914T151231 2282 185 185 0.229 0.23
5 1 sub-2282_ses-20190914T145458 2282 159 159 0.261 0.26
6 2 sub-2763_ses-20200828T161323 2763 50 50 0.403 0.41
7 2 sub-2763_ses-20200828T154656 2763 55 55 0.483 0.48
8 2 sub-2763_ses-20200828T163435 2763 66 66 0.616 0.62
9 3 sub-3445_ses-20210814T223641 3445 179 179 0.274 0.28
10 3 sub-3445_ses-20210814T223641 3445 179 179 0.274 0.28
11 3 sub-3445_ses-20210814T223214 3445 99 99 0.205 0.20
12 4 sub-3443_ses-20210814T193211 3443 114 114 0.290 0.29
13 4 sub-3443_ses-20210814T190234 3443 128 128 0.362 0.36
14 5 sub-3112_ses-20210227T185041 3112 133 133 0.623 0.64
15 5 sub-3112_ses-20210227T190233 3112 133 133 0.692 0.69
16 6 sub-2748_ses-20200928T173432 2748 59 58 0.423 0.43
17 7 sub-2615_ses-20201007T162625 2615 50 50 0.304 0.31
18 8 sub-3278_ses-20210614T151114 3278 101 100 0.241 0.25
19 9 sub-2610_ses-20201005T153308 2610 87 87 0.504 0.50
20 10 sub-2765_ses-20200814T153629 2765 59 59 0.539 0.53
Spines: 20 | ROI pairs in the recordings: 20
Smoothed counts at the burst baseline¶
A burst ends where the lowess-smoothed population count falls to the baseline (the 85th percentile of the counts, an integer). Many smoothed values equal that integer up to rounding error, so whether they pass depended on the last bits of the arithmetic, which differ between machines. The helper rounds the smoothed counts to 9 decimals before the comparison. The cell counts, for the 5-point smoothing, the bins within 1e-9 of the baseline that the unrounded values put above it.
In [7]:
import statsmodels.api as sm
near = []
for f, s in enumerate(sessions):
for strategy in sc.STRATEGIES:
_, _, n, sm_n, _ = sc._burst_counts(s, strategy, "span5")
raw = sm.nonparametric.lowess(n, np.arange(len(n)), frac=5 / len(n))[:, 1]
base = np.percentile(n, 85)
tie = np.abs(raw - base) < 1e-9
near.append({"file": files.file[f], "strategy": strategy, "bins": len(n), "baseline": base,
"within_1e-9": int(tie.sum()), "exactly_equal": int((raw == base).sum()),
"above_unrounded": int((tie & (raw > base)).sum()),
"max_abs_rounding_change": float(np.abs(raw - sm_n).max())})
near = pd.DataFrame(near)
print(near.to_string(index=False))
print("Files with a bin within 1e-9 of the baseline and above it unrounded:",
near[near.above_unrounded > 0].file.nunique(), "of", near.file.nunique())
file strategy bins baseline within_1e-9 exactly_equal above_unrounded max_abs_rounding_change
sub-3112 pc 18400 2.0 315 51 138 4.999846e-10
sub-3112 hz 18400 2.0 378 58 166 4.999982e-10
sub-3112_ses-20210227T190233 pc 18400 2.0 315 51 138 4.999846e-10
sub-3112_ses-20210227T190233 hz 18400 2.0 378 58 166 4.999982e-10
sub-2748_ses-20200928T173432 pc 10600 1.0 379 64 136 4.999881e-10
sub-2748_ses-20200928T173432 hz 10600 1.0 488 80 201 4.991128e-10
sub-2763_ses-20200828T161323 pc 12300 1.0 357 223 51 4.753322e-10
sub-2763_ses-20200828T161323 hz 12300 2.0 405 352 19 4.716941e-10
sub-2763_ses-20200828T154656 pc 12300 2.0 191 43 69 4.999927e-10
sub-2763_ses-20200828T154656 hz 12300 2.0 206 51 83 4.999687e-10
sub-2763_ses-20200828T163435 pc 12300 0.0 10592 10592 0 4.707846e-10
sub-2763_ses-20200828T163435 hz 12300 1.0 1682 1679 2 4.123537e-10
sub-2765_ses-20200814T153629 pc 12300 1.0 1747 1728 5 4.707846e-10
sub-2765_ses-20200814T153629 hz 12300 1.0 352 219 45 4.735130e-10
sub-2615_ses-20201007T162625 pc 12400 1.0 780 117 352 4.998751e-10
sub-2615_ses-20201007T162625 hz 12400 1.0 780 117 352 4.998751e-10
sub-3112_ses-20210227T185041 pc 18400 3.0 189 39 69 4.999046e-10
sub-3112_ses-20210227T185041 hz 18400 3.0 189 39 69 4.999046e-10
sub-2610_ses-20201005T153308 pc 18400 2.0 268 93 82 4.997313e-10
sub-2610_ses-20201005T153308 hz 18400 2.0 268 93 82 4.997313e-10
sub-3443_ses-20210814T190234 pc 18400 2.0 624 104 256 4.999456e-10
sub-3443_ses-20210814T190234 hz 18400 2.0 669 119 281 4.999915e-10
sub-3445_ses-20210814T223641 pc 18400 2.0 416 65 185 4.999812e-10
sub-3445_ses-20210814T223641 hz 18400 2.0 466 72 214 4.999768e-10
sub-3278_ses-20210614T151114 pc 18400 4.0 216 29 97 4.997669e-10
sub-3278_ses-20210614T151114 hz 18400 4.0 207 30 90 4.999743e-10
sub-3445_ses-20210814T223214 pc 18400 1.0 3126 3082 21 4.707846e-10
sub-3445_ses-20210814T223214 hz 18400 1.0 3113 3063 25 4.707846e-10
sub-3443_ses-20210814T193211 pc 20100 1.0 2608 2494 54 4.735130e-10
sub-3443_ses-20210814T193211 hz 20100 1.0 1279 224 536 4.999832e-10
sub-2282_ses-20190914T170145 pc 16400 4.0 301 27 142 4.999501e-10
sub-2282_ses-20190914T170145 hz 16400 4.0 301 27 142 4.999501e-10
sub-2282_ses-20190914T151231 pc 26300 11.0 286 99 102 4.999695e-10
sub-2282_ses-20190914T151231 hz 26300 11.0 290 95 84 4.999821e-10
sub-2282_ses-20190914T164312 pc 18700 14.0 45 2 21 4.999823e-10
sub-2282_ses-20190914T164312 hz 18700 12.0 60 12 23 4.999912e-10
sub-2282_ses-20190914T145458 pc 26300 11.0 322 142 96 4.999716e-10
sub-2282_ses-20190914T145458 hz 26300 11.0 310 117 96 4.999476e-10
sub-2282_ses-20190914T162358 pc 24600 4.0 382 124 132 4.999867e-10
sub-2282_ses-20190914T162358 hz 24600 4.0 382 124 132 4.999867e-10
Files with a bin within 1e-9 of the baseline and above it unrounded: 20 of 20
Claim data-scale¶
The paper made 19 paired spine Ca2+ imaging and HD-MEA recordings in 10 cultures, tested 20 spines (one per recording, two in one recording), and the recordings have a mean of 102 ± 43 sorted units each.
In [8]:
n_units = files.loc[unique, "units"].to_numpy()
print("Recordings:", len(unique), "| subjects (cultures):", files.loc[unique, "subject"].nunique(),
"| spines (ROI pairs):", sum(len(sc.spines(sessions[i])) for i in unique))
print(f"Units per recording, 19 recordings: mean {n_units.mean():.1f}, SD (ddof 0) {n_units.std():.1f}, "
f"SD (ddof 1) {n_units.std(ddof=1):.1f}")
per_spine = np.array([files.units[spine_rows[r]['file']] for r in sorted(spine_rows)])
print(f"Units over the 20 spine rows: mean {per_spine.mean():.1f}, SD (ddof 0) {per_spine.std():.1f}")
link["units_within_10pct"] = (link.units_file - link.units_table).abs() <= 0.1 * link.units_table
nonempty = {r: sum(len(t) > 0 for t in sessions[f]["trains"]) for f, rs in mapping.items() for r in rs}
print("Rows whose unit count differs from Table 1, with the units that have at least one spike:",
[(r, int(files.units[f]), nonempty[r], int(table.units[r])) for f, rs in mapping.items() for r in rs
if files.units[f] != table.units[r]])
print("Per-row unit counts within 10% of Table 1:", int(link.units_within_10pct.sum()), "of", len(link),
"| differing rows:", link.loc[link.units_file != link.units_table, ["row", "units_file", "units_table"]].to_dict("records"))
Recordings: 19 | subjects (cultures): 10 | spines (ROI pairs): 20
Units per recording, 19 recordings: mean 101.7, SD (ddof 0) 42.6, SD (ddof 1) 43.7
Units over the 20 spine rows: mean 105.6, SD (ddof 0) 44.8
Rows whose unit count differs from Table 1, with the units that have at least one spike: [(16, 59, 58, 58), (18, 101, 100, 100)]
Per-row unit counts within 10% of Table 1: 20 of 20 | differing rows: [{'row': 16, 'units_file': 59, 'units_table': 58}, {'row': 18, 'units_file': 101, 'units_table': 100}]
Result: reproduced. The dandiset holds 19 distinct recordings (20 files, two of them the same blob) from 10 subjects, with 20 spine and shaft ROI pairs, two in one recording. Units per recording are 101.7 ± 42.6 (ddof 0) over the 19 recordings, against 102 ± 43. Every per-row unit count is within 10% of Table 1; rows 16 and 18 have one unit more than Table 1, and in both files that extra unit has no spikes, so the counts of units with spikes equal Table 1's.
Claim table1-best-r¶
For each of the 20 tested spines, the best-R (the Pearson R of the demixed spine trace with the convolved trace of its best-matching unit, outside network bursts) has the value given in Table 1, from 0.148 to 0.765.
Every reading of the open burst details is run: p_thresh 96 to 99, hyperactive strategy (i)
99.75th percentile ("pc") or (ii) 15 Hz ("hz"), with the Methods' 5-point smoothing; the code's
lowess (frac 0.0002) is a variant. Units whose convolved trace is constant over the kept frames
have no R; the code's np.argsort would rank them first, so they are left out of the ranking
(their number is printed).
In [9]:
readings = {}
for r, d in spine_rows.items():
s = sessions[d["file"]]
readings[r] = sc.best_r_readings(s, d["spine_roi"], d["shaft_roi"])
allr = pd.DataFrame([{"row": r, **{k: v for k, v in x.items() if not isinstance(v, np.ndarray)}}
for r, xs in readings.items() for x in xs])
allr["table"] = allr.row.map(table.best_R)
allr["match"] = (allr.best_R - allr.table).abs() <= 0.1 * allr.table
main = allr[allr.smoothing == "span5"]
print("Best-R under each reading (5-point smoothing); columns strategy and p_thresh:")
print(main.pivot_table(index="row", columns=["strategy", "p_thresh"], values="best_R").round(3)
.assign(table=table.best_R).to_string())
print("\nUnits without an R (constant over the kept frames), largest per spine:",
main.groupby("row").n_undefined.max().to_dict())
# carried forward: within a recording, the (p_thresh, strategy) pair nearest to Table 1
chosen = {}
for f in set(d["file"] for d in spine_rows.values()):
rs = [r for r, d in spine_rows.items() if d["file"] == f]
best_pair = min(itertools.product(sc.P_THRESH, sc.STRATEGIES), key=lambda ps: sum(
abs(main[(main.row == r) & (main.p_thresh == ps[0]) & (main.strategy == ps[1])].best_R.iloc[0] - table.best_R[r]) for r in rs))
for r in rs:
chosen[r] = next(x for x in readings[r] if x["smoothing"] == "span5" and x["p_thresh"] == best_pair[0]
and x["strategy"] == best_pair[1])
summary = pd.DataFrame([{
"row": r, "file": files.file[spine_rows[r]["file"]], "p_thresh": c["p_thresh"], "strategy": c["strategy"],
"bursts": c["n_bursts"], "kept_frames": c["n_keep"], "of_frames": sessions[spine_rows[r]["file"]]["n_frames"],
"best_unit": c["best"], "best_R": round(c["best_R"], 3), "table": table.best_R[r],
"match_carried": abs(c["best_R"] - table.best_R[r]) <= 0.1 * table.best_R[r],
"match_any": bool(main[main.row == r].match.any()),
"R_97_pc": round(main[(main.row == r) & (main.p_thresh == 97) & (main.strategy == "pc")].best_R.iloc[0], 3),
"R_lowess_frac_same_pair": round(allr[(allr.row == r) & (allr.smoothing == "frac") & (allr.p_thresh == c["p_thresh"])
& (allr.strategy == c["strategy"])].best_R.iloc[0], 3)}
for r, c in sorted(chosen.items())])
print("\nCarried-forward reading per spine:")
print(summary.to_string(index=False))
print("\nSpines matching Table 1 within 10% under at least one reading:", int(summary.match_any.sum()), "of 20")
print("Under the code's default pair (97, pc):",
int(((summary.R_97_pc - summary.table).abs() <= 0.1 * summary.table).sum()), "of 20")
print("With the lowess variant (frac 0.0002), any reading:",
int(allr[allr.smoothing == "frac"].groupby("row").match.any().sum()), "of 20")
print("Best-R range over the 20 spines (carried forward):", summary.best_R.min(), "to", summary.best_R.max())
def n_match(swaps=()):
"""Spines matching under some reading, with the Table 1 values of each pair of rows swapped."""
tab = table.best_R.copy()
for a, b in swaps:
tab[a], tab[b] = table.best_R[b], table.best_R[a]
t = main.row.map(tab)
return int(((main.best_R - t).abs() <= 0.1 * t).groupby(main.row).any().sum())
print("Spines matching under some reading, tied assignments as taken:", n_match(),
"| rows 16 and 20 (subjects 2748, 2765) swapped:", n_match([(16, 20)]),
"| rows 14 and 15 (two sub-3112 files) swapped:", n_match([(14, 15)]),
"| both swapped:", n_match([(16, 20), (14, 15)]))
Best-R under each reading (5-point smoothing); columns strategy and p_thresh:
strategy hz pc table
p_thresh 96 97 98 99 96 97 98 99
row
1 0.489 0.522 0.538 0.574 0.489 0.522 0.538 0.574 0.523
2 0.514 0.494 0.503 0.509 0.506 0.515 0.513 0.510 0.516
3 0.159 0.213 0.221 0.224 0.159 0.213 0.221 0.224 0.217
4 0.655 0.643 0.635 0.615 0.657 0.643 0.635 0.615 0.644
5 0.468 0.430 0.492 0.545 0.472 0.431 0.493 0.545 0.429
6 0.599 0.669 0.639 0.653 0.631 0.594 0.698 0.659 0.647
7 0.290 0.364 0.344 0.602 0.269 0.360 0.344 0.605 0.346
8 0.462 0.441 0.442 0.363 0.477 0.457 0.454 0.370 0.340
9 0.673 0.680 0.679 0.639 0.673 0.680 0.679 0.639 0.671
10 0.762 0.766 0.768 0.729 0.762 0.766 0.768 0.729 0.765
11 0.189 0.189 0.178 0.115 0.136 0.136 0.163 0.137 0.148
12 0.585 0.584 0.583 0.432 0.580 0.584 0.584 0.436 0.597
13 0.133 0.131 0.118 0.139 0.141 0.132 0.117 0.143 0.160
14 0.354 0.466 0.525 0.549 0.354 0.466 0.525 0.549 0.394
15 0.322 0.395 0.627 0.701 0.320 0.393 0.595 0.700 0.400
16 0.366 0.355 0.454 0.494 0.366 0.363 0.480 0.537 0.676
17 0.642 0.642 0.632 0.627 0.642 0.642 0.632 0.627 0.602
18 0.318 0.326 0.319 0.344 0.311 0.317 0.317 0.341 0.311
19 0.149 0.222 0.387 0.415 0.149 0.222 0.387 0.415 0.219
20 0.283 0.390 0.379 0.494 0.293 0.370 0.383 0.481 0.294
Units without an R (constant over the kept frames), largest per spine: {1: 0, 2: 0, 3: 0, 4: 0, 5: 0, 6: 0, 7: 0, 8: 0, 9: 0, 10: 0, 11: 0, 12: 0, 13: 0, 14: 0, 15: 0, 16: 1, 17: 0, 18: 1, 19: 0, 20: 0}
Carried-forward reading per spine:
row file p_thresh strategy bursts kept_frames of_frames best_unit best_R table match_carried match_any R_97_pc R_lowess_frac_same_pair
1 sub-2282_ses-20190914T162358 97 pc 106 1210 1800 76 0.522 0.523 True True 0.522 0.522
2 sub-2282_ses-20190914T164312 97 pc 46 1080 1364 116 0.515 0.516 True True 0.515 0.515
3 sub-2282_ses-20190914T170145 97 pc 76 829 1200 37 0.213 0.217 True True 0.213 0.218
4 sub-2282_ses-20190914T151231 97 pc 201 1099 1800 183 0.643 0.644 True True 0.643 0.643
5 sub-2282_ses-20190914T145458 97 hz 237 1032 1800 60 0.430 0.429 True True 0.431 0.430
6 sub-2763_ses-20200828T161323 99 hz 43 927 1200 1 0.653 0.647 True True 0.594 0.653
7 sub-2763_ses-20200828T154656 98 pc 42 918 1200 27 0.344 0.346 True True 0.360 0.338
8 sub-2763_ses-20200828T163435 99 hz 31 973 1200 32 0.363 0.340 True True 0.457 0.363
9 sub-3445_ses-20210814T223641 96 pc 102 1124 1800 115 0.673 0.671 True True 0.680 0.671
10 sub-3445_ses-20210814T223641 96 pc 102 1124 1800 115 0.762 0.765 True True 0.766 0.762
11 sub-3445_ses-20210814T223214 99 pc 62 1406 1800 59 0.137 0.148 True True 0.136 0.137
12 sub-3443_ses-20210814T193211 96 hz 99 1288 1965 30 0.585 0.597 True True 0.584 0.585
13 sub-3443_ses-20210814T190234 99 pc 76 1261 1800 111 0.143 0.160 False False 0.132 0.146
14 sub-3112_ses-20210227T185041 96 pc 42 1435 1800 13 0.354 0.394 False False 0.466 0.359
15 sub-3112_ses-20210227T190233 97 hz 39 1460 1800 31 0.395 0.400 True True 0.393 0.393
16 sub-2748_ses-20200928T173432 99 pc 32 813 1032 3 0.537 0.676 False False 0.363 0.525
17 sub-2615_ses-20201007T162625 99 pc 25 1514 1791 27 0.627 0.602 True True 0.642 0.631
18 sub-3278_ses-20210614T151114 96 pc 124 1086 1800 53 0.311 0.311 True True 0.317 0.301
19 sub-2610_ses-20201005T153308 97 pc 81 1217 1800 69 0.222 0.219 True True 0.222 0.213
20 sub-2765_ses-20200814T153629 96 pc 58 856 1200 51 0.293 0.294 True True 0.370 0.293
Spines matching Table 1 within 10% under at least one reading: 17 of 20
Under the code's default pair (97, pc): 15 of 20
With the lowess variant (frac 0.0002), any reading: 19 of 20
Best-R range over the 20 spines (carried forward): 0.137 to 0.762
Spines matching under some reading, tied assignments as taken: 17 | rows 16 and 20 (subjects 2748, 2765) swapped: 16 | rows 14 and 15 (two sub-3112 files) swapped: 17 | both swapped: 16
In [10]:
fig, ax = plt.subplots(figsize=(6, 5))
for r in sorted(chosen):
v = main[main.row == r].best_R
ax.plot([table.best_R[r]] * len(v), v, ".", color="0.6")
ax.plot(summary.table, summary.best_R, "o", color="C0", label="carried forward")
ax.plot([0, 0.85], [0, 0.85], "k-", lw=0.8)
ax.plot([0, 0.85], [0, 0.85 * 1.1], "k:", lw=0.8)
ax.plot([0, 0.85], [0, 0.85 * 0.9], "k:", lw=0.8)
ax.set_xlabel("Table 1 best-R")
ax.set_ylabel("best-R here (grey: every reading)")
ax.set_xlim(0, 0.85)
ax.set_ylim(0, 0.85)
ax.legend()
print("Plotted best-R range:", round(main.best_R.min(), 3), "to", round(main.best_R.max(), 3), "(all inside the axes)")
Plotted best-R range: 0.115 to 0.768 (all inside the axes)
Result: partially reproduced. 17 of the 20 spines have a best-R within 10% of Table 1 under at least one reading of p_thresh and hyperactive strategy; with the code's lowess smoothing (frac 0.0002), 19 of 20. The three that miss are row 13 (nearest 0.143 against 0.160), row 14 (nearest 0.354, just below the bound of 0.355 for 0.394) and row 16 (0.537 at most against 0.676). The code's default pair (97, pc) matches 14 of 20. The row-to-file assignment of rows 14, 15, 16 and 20 rests on the spike-contrast tie-break, but the outcome does not: with rows 16 and 20 swapped 16 spines match, with rows 14 and 15 swapped 17, and with both swapped 16, all in the band of 15 to 19.
Within a spine, best-R moves widely across the eight readings (row 15, 0.320 to 0.701; row 7, 0.269 to 0.605), so a match under some reading is a weak test, and every later claim (thresholds, F grids, OUT/IN) inherits the reading chosen for being nearest to Table 1. Its outcomes are conditional on that choice. With every chain of close bursts merged into one, as Methods 2.6 says, in place of the code's single pass (claim table1-network-measures), 16 of 20 carried-forward best-R values match, against 17; row 20 moves from 0.293 to 0.340.
Claim fig1e-best-r-gap¶
In the example recording of Figure 1(e), the best-R shows the largest increase over the next ordered R value among all units.
In [11]:
gaps = []
for r, c in sorted(chosen.items()):
v = np.sort(c["r"][np.isfinite(c["r"])])[::-1]
g = -np.diff(v)
gaps.append({"row": r, "best_R": round(v[0], 3), "second_R": round(v[1], 3), "first_gap": round(g[0], 3),
"largest_other_gap": round(g[1:].max(), 3), "first_gap_largest": bool(g[0] > g[1:].max()),
"rank_of_first_gap": int(1 + np.sum(g[1:] > g[0])), "units_with_R": len(v)})
gaps = pd.DataFrame(gaps)
print(gaps.to_string(index=False))
print("Spines whose first gap is the largest:", int(gaps.first_gap_largest.sum()), "of", len(gaps))
row best_R second_R first_gap largest_other_gap first_gap_largest rank_of_first_gap units_with_R 1 0.522 0.311 0.212 0.060 True 1 78 2 0.515 0.385 0.130 0.058 True 1 132 3 0.213 0.165 0.048 0.048 True 1 66 4 0.643 0.460 0.183 0.128 True 1 185 5 0.430 0.426 0.004 0.039 False 45 159 6 0.653 0.426 0.227 0.067 True 1 50 7 0.344 0.328 0.016 0.081 False 6 55 8 0.363 0.313 0.050 0.066 False 2 66 9 0.673 0.634 0.039 0.086 False 7 179 10 0.762 0.714 0.048 0.168 False 3 179 11 0.137 0.134 0.003 0.030 False 17 99 12 0.585 0.564 0.021 0.042 False 12 114 13 0.143 0.118 0.026 0.015 True 1 128 14 0.354 0.214 0.140 0.025 True 1 133 15 0.395 0.359 0.036 0.079 False 2 133 16 0.537 0.507 0.029 0.055 False 6 58 17 0.627 0.610 0.016 0.089 False 13 50 18 0.311 0.307 0.004 0.066 False 26 100 19 0.222 0.214 0.008 0.058 False 21 87 20 0.293 0.285 0.008 0.027 False 21 59 Spines whose first gap is the largest: 7 of 20
Result: reproduced. With the carried-forward reading, the gap between best-R and the second R is the largest of all consecutive gaps for 7 of the 20 spines (rows 1, 2, 3, 4, 6, 13 and 14), so the example of Figure 1(e) can be one of them. For the other 13 it is not (rank 2 to 45). The claim is about one unnamed example, so this is a weak test.
Noise and signal estimates¶
The surrogates need the SD of the spine's pixels and of background pixels "close to the spine"
(Methods 2.7). The dandiset holds the spine and shaft ROIs but no background ROI, so, as the
claims file prescribes, three background squares near each spine are placed by a rule fixed
before any R was computed (spineca.background_squares), and each is run. The image series is
read in a 96 by 96 pixel window around each spine, over all frames.
The figure after the table shows each window's mean image with the spine ROI, the other ROIs and the three squares. In most windows the squares lie in dark gaps with no visible process. In row 17, where processes fill the window, they lie in the dimmest gaps between them. In row 10 the rule fails its aim for regions 0 and 2, which sit on bright processes: their noise SD is 126 and 96, against 39 for region 1. Row 10's decision is the same in every region (claim acceptance-9-of-20), and the 19-recording readings of the F-score claims use row 9's estimate for that recording. The rule was fixed before any R was computed and is kept.
In [12]:
t0 = time.time()
bounds = {r: sc.patch_bounds(sessions[d["file"]], d["spine_roi"]) for r, d in spine_rows.items()}
patches = sc.read_patches([(r, sessions[d["file"]]["asset_id"], sessions[d["file"]]["n_frames"], bounds[r])
for r, d in spine_rows.items()], workers=24, frames_per_block=75)
print(f"Read {len(patches)} windows in {time.time() - t0:.0f} s")
regions, estimates, rows_n, mean_imgs = {}, {}, [], {}
for r, d in sorted(spine_rows.items()):
s, P = sessions[d["file"]], patches[r]
yy, xx = np.nonzero(s["masks"][d["spine_roi"]])
mean_trace = P[:, yy - bounds[r][0], xx - bounds[r][2]].mean(axis=1)
mean_img = P.mean(axis=0)
mean_imgs[r] = mean_img
squares, cut, n_cand = sc.background_squares(s, d["spine_roi"], bounds[r], mean_img)
regions[r] = squares
estimates[r] = []
for g, q in enumerate(squares):
e = sc.noise_signal(P, s, d["spine_roi"], bounds[r], q, np.random.default_rng([223, 1, r, g]))
estimates[r].append(e)
rows_n.append({"row": r, "region": g, "square_rc": (int(q[0] + bounds[r][0]), int(q[1] + bounds[r][2])),
"dist_px": round(q[4], 1), "square_mean": round(q[3], 1), "window_p25": round(cut, 1),
"spine_mean": round(mean_img[yy - bounds[r][0], xx - bounds[r][2]].mean(), 1),
"spine_px": len(yy), "signal_amp": round(e["amp"], 2), "noise_sd": round(e["sigma"], 2),
"baseline": round(e["baseline"], 1),
"roi_trace_equal": bool(np.allclose(mean_trace, s["roi_traces"][:, d["spine_roi"]]))})
noise_table = pd.DataFrame(rows_n)
print(noise_table.to_string(index=False))
print("Spine ROI trace in the file equal to the mean of its mask pixels in every window:", noise_table.roi_trace_equal.all())
print("Image values range over the windows:", min(int(p.min()) for p in patches.values()), "to",
max(int(p.max()) for p in patches.values()), "(uint16, no unit given)")
Read 20 windows in 290 s
row region square_rc dist_px square_mean window_p25 spine_mean spine_px signal_amp noise_sd baseline roi_trace_equal 1 0 (815, 469) 12.7 520.4 524.3 754.3 10 117.02 32.57 715.5 True 1 1 (821, 493) 15.0 524.2 524.3 754.3 10 114.52 34.08 715.5 True 1 2 (827, 493) 15.2 523.7 524.3 754.3 10 115.14 34.44 715.5 True 2 0 (271, 561) 22.3 524.9 525.3 731.4 9 95.57 36.40 714.0 True 2 1 (277, 603) 26.0 521.4 525.3 731.4 9 97.01 33.79 714.0 True 2 2 (259, 609) 27.6 525.0 525.3 731.4 9 94.95 35.92 714.0 True 3 0 (24, 673) 9.5 519.2 520.7 880.6 5 139.86 31.35 861.0 True 3 1 (30, 691) 9.8 520.5 520.7 880.6 5 140.99 31.45 861.0 True 3 2 (30, 673) 10.6 518.1 520.7 880.6 5 142.71 31.47 861.0 True 4 0 (388, 430) 9.8 621.3 634.0 2595.2 4 721.15 69.21 2420.8 True 4 1 (388, 436) 12.7 629.2 634.0 2595.2 4 716.28 73.56 2420.8 True 4 2 (382, 412) 15.9 629.5 634.0 2595.2 4 716.26 73.25 2420.8 True 5 0 (110, 176) 9.7 629.3 653.2 1562.5 2 425.22 69.62 1434.8 True 5 1 (98, 200) 17.2 649.4 653.2 1562.5 2 420.72 74.12 1434.8 True 5 2 (128, 176) 23.3 646.7 653.2 1562.5 2 424.15 70.68 1434.8 True 6 0 (108, 510) 27.0 503.2 509.5 1322.5 8 452.91 22.11 1216.5 True 6 1 (108, 504) 27.1 507.5 509.5 1322.5 8 452.15 22.65 1216.5 True 6 2 (108, 516) 28.3 504.4 509.5 1322.5 8 453.11 21.92 1216.5 True 7 0 (104, 500) 27.6 508.6 511.7 1012.3 4 239.58 24.26 946.5 True 7 1 (104, 506) 27.7 503.2 511.7 1012.3 4 241.76 21.86 946.5 True 7 2 (104, 512) 29.1 504.2 511.7 1012.3 4 241.06 22.96 946.5 True 8 0 (75, 467) 21.2 510.9 515.2 1374.0 13 496.96 24.15 1275.5 True 8 1 (117, 449) 27.2 511.2 515.2 1374.0 13 495.85 24.81 1275.5 True 8 2 (117, 455) 27.2 508.6 515.2 1374.0 13 496.18 23.97 1275.5 True 9 0 (278, 629) 9.1 545.1 556.6 2090.5 9 1034.39 33.26 1688.5 True 9 1 (278, 623) 14.9 548.8 556.6 2090.5 9 1031.97 35.59 1688.5 True 9 2 (266, 629) 17.2 556.1 556.6 2090.5 9 1028.08 36.40 1688.5 True 10 0 (42, 610) 27.0 737.7 797.2 2224.7 8 955.93 126.15 1761.5 True 10 1 (36, 610) 27.6 583.6 797.2 2224.7 8 1029.61 38.82 1761.5 True 10 2 (30, 556) 29.3 736.8 797.2 2224.7 8 985.64 95.92 1761.5 True 11 0 (86, 760) 10.0 569.3 571.3 871.3 6 169.51 42.39 808.5 True 11 1 (80, 760) 13.1 569.3 571.3 871.3 6 172.83 41.55 808.5 True 11 2 (86, 754) 15.8 551.4 571.3 871.3 6 178.35 34.13 808.5 True 12 0 (177, 238) 12.9 529.1 573.6 1517.8 7 392.78 28.23 1356.0 True 12 1 (183, 244) 15.4 562.9 573.6 1517.8 7 388.77 33.91 1356.0 True 12 2 (183, 238) 17.6 535.9 573.6 1517.8 7 392.40 31.14 1356.0 True 13 0 (219, 589) 12.2 546.7 567.9 1032.5 13 235.88 34.42 983.0 True 13 1 (219, 595) 17.0 552.4 567.9 1032.5 13 233.72 35.19 983.0 True 13 2 (219, 565) 17.6 554.7 567.9 1032.5 13 229.30 36.22 983.0 True 14 0 (367, 389) 17.4 501.1 501.2 690.8 6 138.33 22.59 649.0 True 14 1 (361, 383) 17.6 499.6 501.2 690.8 6 139.30 21.82 649.0 True 14 2 (373, 395) 20.9 499.3 501.2 690.8 6 139.89 21.64 649.0 True 15 0 (724, 554) 12.4 525.6 526.1 799.0 6 200.18 33.80 723.5 True 15 1 (718, 548) 14.8 525.2 526.1 799.0 6 197.14 32.15 723.5 True 15 2 (730, 560) 15.2 523.6 526.1 799.0 6 198.83 32.46 723.5 True 16 0 (748, 590) 9.4 514.9 519.1 894.0 12 271.62 26.93 800.8 True 16 1 (766, 584) 15.2 513.5 519.1 894.0 12 273.88 25.72 800.8 True 16 2 (748, 596) 15.2 512.5 519.1 894.0 12 273.60 25.40 800.8 True 17 0 (574, 378) 9.0 569.9 636.5 1030.6 5 285.31 44.52 903.0 True 17 1 (574, 384) 9.2 573.1 636.5 1030.6 5 291.43 45.37 903.0 True 17 2 (574, 372) 12.3 586.8 636.5 1030.6 5 273.87 48.36 903.0 True 18 0 (724, 205) 12.9 513.0 517.5 856.4 10 118.94 25.12 843.5 True 18 1 (736, 199) 15.3 510.9 517.5 856.4 10 120.54 25.64 843.5 True 18 2 (730, 199) 15.3 507.8 517.5 856.4 10 122.21 24.68 843.5 True 19 0 (405, 437) 9.1 534.3 541.0 907.6 16 246.60 32.90 823.5 True 19 1 (405, 443) 9.3 521.3 541.0 907.6 16 254.45 28.71 823.5 True 19 2 (411, 449) 9.6 534.5 541.0 907.6 16 253.33 31.61 823.5 True 20 0 (720, 95) 15.7 517.9 525.4 801.3 7 196.32 27.23 753.0 True 20 1 (720, 89) 15.8 523.2 525.4 801.3 7 196.45 28.95 753.0 True 20 2 (720, 101) 17.7 515.3 525.4 801.3 7 198.16 26.97 753.0 True Spine ROI trace in the file equal to the mean of its mask pixels in every window: True Image values range over the windows: 402 to 23416 (uint16, no unit given)
In [13]:
# Mean image of each spine's window: spine ROI (red), other ROIs (white), the three background squares
fig, axs = plt.subplots(4, 5, figsize=(15, 12.5))
for ax, (r, d) in zip(axs.ravel(), sorted(spine_rows.items())):
s, (r0, r1, c0, c1) = sessions[d["file"]], bounds[r]
ax.imshow(mean_imgs[r], cmap="gray", vmax=np.percentile(mean_imgs[r], 99.5))
others = s["masks"].any(axis=0)[r0:r1, c0:c1]
ax.contour(others, levels=[0.5], colors="w", linewidths=0.5)
ax.contour(s["masks"][d["spine_roi"]][r0:r1, c0:c1], levels=[0.5], colors="r", linewidths=1)
for g, (i, j, side, _, _) in enumerate(regions[r]):
ax.add_patch(plt.Rectangle((j - 0.5, i - 0.5), side, side, fill=False, color="cyan", lw=1))
ax.text(j, i - 1, str(g), color="cyan", fontsize=8)
ax.set_title(f"row {r}", fontsize=9)
ax.axis("off")
plt.tight_layout()
plt.show()
Claim acceptance-9-of-20¶
The best-matching unit passed the acceptance test, and a monosynaptic connection was accepted, for 9 of the 20 spines (Table 1 rows 1, 2, 3, 4, 6, 9, 10, 12 and 16).
1'000 GT surrogates of the best unit and 1'000 test surrogates of the second-best unit, R with the best unit's convolved trace over the kept frames, for each of the three noise regions. Reading (a) takes Table 1's quantiles; reading (b), the optimal quantiles found in claim table1-optimal-thresholds, is run in that section, below.
In [14]:
pr_main = sc.release_probabilities(np.random.default_rng([223, 0]))
print(f"Release probabilities: {len(pr_main)}, median {np.median(pr_main):.3f}, range {pr_main.min():.4f} to {pr_main.max():.3f}")
dist = {}
for r, d in sorted(spine_rows.items()):
s, c = sessions[d["file"]], chosen[r]
for g, e in enumerate(estimates[r]):
dist[(r, g)] = sc.acceptance_distributions(s, s["U"], c["keep"], c["best"], c["second"], pr_main, e,
np.random.default_rng([223, 2, r, g]))
def decide(gt, te, q_gt, q_test, best_r):
g, t = np.nanquantile(gt, q_gt), np.nanquantile(te, q_test)
return g, t, bool(best_r > max(g, t))
acc_a = []
for (r, g), (gt, te) in dist.items():
gthr, tthr, ok = decide(gt, te, table.GT_Q[r], table.test_Q[r], chosen[r]["best_R"])
acc_a.append({"row": r, "region": g, "best_R": round(chosen[r]["best_R"], 3), "GT_thr": round(gthr, 3),
"GT_thr_table": table.GT_thr[r], "test_thr": round(tthr, 3), "test_thr_table": table.test_thr[r],
"accept": ok, "accept_table": table.accept[r]})
acc_a = pd.DataFrame(acc_a)
acc_a["agree"] = acc_a.accept == acc_a.accept_table
acc_a["GT_thr_10pct"] = (acc_a.GT_thr - acc_a.GT_thr_table).abs() <= 0.1 * acc_a.GT_thr_table
acc_a["test_thr_10pct"] = (acc_a.test_thr - acc_a.test_thr_table).abs() <= 0.1 * acc_a.test_thr_table
print(acc_a[acc_a.region == 0].drop(columns="region").to_string(index=False))
print("\nPer noise region, reading (a): accepted, decisions agreeing with Table 1, thresholds within 10%:")
print(acc_a.groupby("region").agg(accepted=("accept", "sum"), agree=("agree", "sum"),
GT_thr_10pct=("GT_thr_10pct", "sum"), test_thr_10pct=("test_thr_10pct", "sum")).to_string())
print("\nDecision per spine and region (1 = accept):")
print(acc_a.pivot_table(index="row", columns="region", values="accept").astype(int).assign(table=table.accept.astype(int)).T.to_string())
Release probabilities: 1000, median 0.123, range 0.0161 to 0.741
row best_R GT_thr GT_thr_table test_thr test_thr_table accept accept_table agree GT_thr_10pct test_thr_10pct
1 0.522 0.195 0.214 0.301 0.296 True True True True True
2 0.515 0.341 0.291 0.338 0.271 True True True False False
3 0.213 0.112 0.112 0.264 0.190 False True False True False
4 0.643 0.288 0.311 0.481 0.489 True True True True True
5 0.430 0.348 0.358 0.546 0.520 False False True True True
6 0.653 0.420 0.448 0.607 0.619 True True True True True
7 0.344 0.398 0.503 0.660 0.585 False False True False False
8 0.363 0.226 0.227 0.433 0.475 False False True True True
9 0.673 0.497 0.507 0.649 0.639 True True True True True
10 0.762 0.483 0.507 0.632 0.639 True True True True True
11 0.137 0.264 0.354 0.258 0.351 False False True False False
12 0.585 0.331 0.353 0.587 0.544 False True False True True
13 0.143 0.250 0.267 0.103 0.068 False False True True False
14 0.354 0.412 0.456 0.347 0.435 False False True True False
15 0.395 0.503 0.436 0.578 0.604 False False True False True
16 0.537 0.684 0.383 0.437 0.666 False True False False False
17 0.627 0.361 0.500 0.716 0.707 False False True False True
18 0.311 0.196 0.214 0.330 0.324 False False True True True
19 0.222 0.449 0.233 0.190 0.298 False False True False False
20 0.293 0.417 0.414 0.538 0.364 False False True True False
Per noise region, reading (a): accepted, decisions agreeing with Table 1, thresholds within 10%:
accepted agree GT_thr_10pct test_thr_10pct
region
0 6 17 13 11
1 7 18 13 11
2 7 18 12 10
Decision per spine and region (1 = accept):
row 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
region
0 1 1 0 1 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0
1 1 1 0 1 0 1 0 0 1 1 0 1 0 0 0 0 0 0 0 0
2 1 1 0 1 0 1 0 0 1 1 0 1 0 0 0 0 0 0 0 0
table 1 1 1 1 0 1 0 0 1 1 0 1 0 0 0 1 0 0 0 0
In [15]:
# Rows near a threshold in Table 1: decisions over five seeds, every region, reading (a)
seeds = []
for r in [x for x in [3, 6, 16, 18] if x in spine_rows]:
d, c = spine_rows[r], chosen[r]
s = sessions[d["file"]]
for g, e in enumerate(estimates[r]):
for seed in range(5):
pr = sc.release_probabilities(np.random.default_rng([223, 5, seed]))
gt, te = sc.acceptance_distributions(s, s["U"], c["keep"], c["best"], c["second"], pr, e,
np.random.default_rng([223, 6, r, g, seed]))
seeds.append({"row": r, "region": g, "seed": seed, "accept": decide(gt, te, table.GT_Q[r], table.test_Q[r], c["best_R"])[2]})
seeds = pd.DataFrame(seeds)
print(seeds.groupby(["row", "region"]).accept.agg(["sum", "count"]).rename(columns={"sum": "accepted"})
.assign(table=lambda x: [table.accept[i] for i, _ in x.index]).to_string())
# Variants on region 0: the code's stored release probabilities, and a Bernoulli draw per spike
url = ("https://raw.githubusercontent.com/starquakes/mea_spineca_mapping/"
"99c2a478393341a24c58d6ef6ab5db5b0b528d03/notebooks/log_release_Pr.npy")
pr_stored = np.load(io.BytesIO(requests.get(url, timeout=60).content))
print(f"\nStored release probabilities: {len(pr_stored)}, median {np.median(pr_stored):.3f}")
var = []
for r, d in sorted(spine_rows.items()):
s, c = sessions[d["file"]], chosen[r]
for name, pr, bern in [("stored Pr", pr_stored, False), ("Bernoulli", pr_main, True)]:
gt, te = sc.acceptance_distributions(s, s["U"], c["keep"], c["best"], c["second"], pr, estimates[r][0],
np.random.default_rng([223, 7, r]), bernoulli=bern)
var.append({"row": r, "variant": name, "accept": decide(gt, te, table.GT_Q[r], table.test_Q[r], c["best_R"])[2]})
var = pd.DataFrame(var)
var["agree"] = var.accept == var.row.map(table.accept)
print(var.groupby("variant").agg(accepted=("accept", "sum"), agree=("agree", "sum")).to_string())
accepted count table
row region
3 0 0 5 True
1 0 5 True
2 0 5 True
6 0 5 5 True
1 5 5 True
2 5 5 True
16 0 0 5 True
1 0 5 True
2 0 5 True
18 0 0 5 False
1 0 5 False
2 0 5 False
Stored release probabilities: 1000, median 0.133
accepted agree variant Bernoulli 6 17 stored Pr 6 17
In [16]:
# Best and second-best unit of each spine (carried-forward reading); the test surrogates are of the second
sb = pd.DataFrame([{"row": r, "best_unit": c["best"], "best_R": round(c["best_R"], 3), "second_unit": c["second"],
"second_R": round(c["second_R"], 3), "accept_table": table.accept[r]} for r, c in sorted(chosen.items())])
print(sb.to_string(index=False))
# GT-R and test-R distributions of rows 3 and 16 (accepted in Table 1, rejected here), every region
qs = [0.05, 0.25, 0.5, 0.75, 0.95]
dq = []
for r in [3, 16]:
for g in range(3):
gt, te = dist[(r, g)]
for name, x, tq, tt in [("GT-R", gt, table.GT_Q[r], table.GT_thr[r]), ("test-R", te, table.test_Q[r], table.test_thr[r])]:
dq.append({"row": r, "region": g, "dist": name, **{f"q{q:g}": round(np.nanquantile(x, q), 3) for q in qs},
"Q_table": tq, "at_Q_table": round(np.nanquantile(x, tq), 3), "thr_table": tt,
"best_R": round(chosen[r]["best_R"], 3), "best_R_table": table.best_R[r]})
print()
print(pd.DataFrame(dq).to_string(index=False))
# Rows 9 and 10, two spines of one recording: thresholds of row 10 with row 9's noise estimate and seed
if chosen[9]["best"] == chosen[10]["best"] and chosen[9]["second"] == chosen[10]["second"]:
for g in range(3):
d, c = spine_rows[10], chosen[10]
gt, te = sc.acceptance_distributions(sessions[d["file"]], sessions[d["file"]]["U"], c["keep"], c["best"], c["second"],
pr_main, estimates[9][g], np.random.default_rng([223, 2, 9, g]))
g9 = decide(*dist[(9, g)], table.GT_Q[9], table.test_Q[9], chosen[9]["best_R"])
g10 = decide(gt, te, table.GT_Q[10], table.test_Q[10], c["best_R"])
print(f"Region {g}: row 9 GT/test thr {g9[0]:.3f}/{g9[1]:.3f}; row 10 with row 9's estimate {g10[0]:.3f}/{g10[1]:.3f} "
f"(Table 1, {table.GT_thr[10]}/{table.test_thr[10]} for both)")
row best_unit best_R second_unit second_R accept_table 1 76 0.522 30 0.311 True 2 116 0.515 109 0.385 True 3 37 0.213 23 0.165 True 4 183 0.643 179 0.460 True 5 60 0.430 122 0.426 False 6 1 0.653 43 0.426 True 7 27 0.344 23 0.328 False 8 32 0.363 57 0.313 False 9 115 0.673 44 0.634 True 10 115 0.762 44 0.714 True 11 59 0.137 21 0.134 False 12 30 0.585 57 0.564 True 13 111 0.143 125 0.118 False 14 13 0.354 14 0.214 False 15 31 0.395 14 0.359 False 16 3 0.537 30 0.507 True 17 27 0.627 45 0.610 False 18 53 0.311 100 0.307 False 19 69 0.222 8 0.214 False 20 51 0.293 46 0.285 False row region dist q0.05 q0.25 q0.5 q0.75 q0.95 Q_table at_Q_table thr_table best_R best_R_table 3 0 GT-R 0.058 0.135 0.205 0.281 0.411 0.172 0.112 0.112 0.213 0.217 3 0 test-R -0.002 0.077 0.130 0.184 0.257 0.960 0.264 0.190 0.213 0.217 3 1 GT-R 0.062 0.142 0.206 0.286 0.418 0.172 0.115 0.112 0.213 0.217 3 1 test-R -0.006 0.072 0.127 0.180 0.254 0.960 0.266 0.190 0.213 0.217 3 2 GT-R 0.058 0.146 0.199 0.280 0.413 0.172 0.121 0.112 0.213 0.217 3 2 test-R -0.004 0.067 0.123 0.177 0.250 0.960 0.259 0.190 0.213 0.217 16 0 GT-R 0.475 0.616 0.711 0.784 0.862 0.424 0.684 0.383 0.537 0.676 16 0 test-R 0.189 0.344 0.425 0.484 0.562 0.545 0.437 0.666 0.537 0.676 16 1 GT-R 0.467 0.615 0.711 0.792 0.869 0.424 0.683 0.383 0.537 0.676 16 1 test-R 0.181 0.346 0.419 0.483 0.561 0.545 0.430 0.666 0.537 0.676 16 2 GT-R 0.455 0.613 0.711 0.790 0.865 0.424 0.684 0.383 0.537 0.676 16 2 test-R 0.188 0.345 0.420 0.483 0.574 0.545 0.434 0.666 0.537 0.676
Region 0: row 9 GT/test thr 0.497/0.649; row 10 with row 9's estimate 0.497/0.649 (Table 1, 0.507/0.639 for both)
Region 1: row 9 GT/test thr 0.479/0.653; row 10 with row 9's estimate 0.479/0.653 (Table 1, 0.507/0.639 for both)
Region 2: row 9 GT/test thr 0.493/0.662; row 10 with row 9's estimate 0.493/0.662 (Table 1, 0.507/0.639 for both)
Claim table1-optimal-thresholds¶
For each recording, simulated ground-truth connections of every unit give an F-score surface over the quantile pair (Q_test, Q_GT), whose maximum and optimal quantiles are those of Table 1 (max_F 0.485 to 0.870).
1'000 surrogates of every unit, for each spine's noise estimate and each of its three noise regions. Main reading: the Methods' classification of a surrogate whose best unit is not the ground truth (FP if it passes, else TN); the code's (FN) is a variant.
In [17]:
t0 = time.time()
grid = []
with ProcessPoolExecutor(4) as pool:
for r, d in sorted(spine_rows.items()):
s, c = sessions[d["file"]], chosen[r]
for g, e in enumerate(estimates[r]):
R = sc.r_tensor_parallel(s, s["U"], c["keep"], pr_main, e, seed=223000 + 10 * r + g, pool=pool)
fg = sc.fscore_grid(R)
(Fm, cm), (Fc, _) = fg["methods"], fg["code"]
mm, om = sc.grid_optimum(Fm)
mc, oc = sc.grid_optimum(Fc)
i_t, i_g = np.argwhere(Fm == mm)[0]
grid.append({"row": r, "region": g, "units": len(s["trains"]), "maxF": round(mm, 3), "table_maxF": table.max_F[r],
"Q_test": om[0][0], "Q_GT": om[0][1], "n_max": len(om), "table_Q_test": table.test_Q[r],
"table_Q_GT": table.GT_Q[r], "maxF_code": round(mc, 3), "Q_test_code": oc[0][0],
"Q_GT_code": oc[0][1], "all_max": om,
"TP": int(cm["TP"][i_t, i_g]), "FN": int(cm["FN"][i_t, i_g]), "FP": int(cm["FP"][i_t, i_g]),
"TN": int(cm["TN"][i_t, i_g]), **{k: round(v, 3) for k, v in fg["stats"].items()}})
del R
grid = pd.DataFrame(grid)
print(f"F-score grids computed in {time.time() - t0:.0f} s")
grid["maxF_10pct"] = (grid.maxF - grid.table_maxF).abs() <= 0.1 * grid.table_maxF
grid["maxF_code_10pct"] = (grid.maxF_code - grid.table_maxF).abs() <= 0.1 * grid.table_maxF
grid["Q_10pct"] = [any(abs(qt - t) <= 0.1 * t and abs(qg - gq) <= 0.1 * gq for qt, qg in m)
for m, t, gq in zip(grid.all_max, grid.table_Q_test, grid.table_Q_GT)]
print(grid.drop(columns="all_max").to_string(index=False))
recs = grid[grid.row != 10] # recording 3 counted once (row 9's estimate); row 10 shown above
print("\nRecordings (19) with max F within 10% of Table 1, per noise region; Methods reading, code reading, "
"and with Q_test and Q_GT also within 10%:")
print(recs.groupby("region").agg(maxF=("maxF_10pct", "sum"), maxF_code=("maxF_code_10pct", "sum"),
quantiles=("Q_10pct", "sum")).to_string())
F-score grids computed in 700 s
row region units maxF table_maxF Q_test Q_GT n_max table_Q_test table_Q_GT maxF_code Q_test_code Q_GT_code TP FN FP TN hit_fraction gt_r_median test_r_median maxF_10pct maxF_code_10pct Q_10pct
1 0 78 0.783 0.787 0.970 0.101 1 0.960 0.091 0.764 0.960 0.111 34040 36268 2746 82946 0.901 0.380 0.221 True True False
1 1 78 0.782 0.787 0.960 0.141 1 0.960 0.091 0.764 0.960 0.121 36987 33324 4569 81120 0.901 0.375 0.218 True True False
1 2 78 0.781 0.787 0.960 0.141 1 0.960 0.091 0.763 0.960 0.131 37135 33211 4746 80908 0.902 0.378 0.219 True True False
2 0 132 0.789 0.854 0.970 0.242 1 0.960 0.162 0.780 0.970 0.192 68391 53079 9541 132989 0.920 0.358 0.164 True True False
2 1 132 0.790 0.854 0.970 0.212 1 0.960 0.162 0.780 0.970 0.192 69129 52362 9919 132590 0.920 0.363 0.169 True True False
2 2 132 0.790 0.854 0.970 0.212 1 0.960 0.162 0.781 0.970 0.212 69378 52128 9979 132515 0.920 0.356 0.166 True True False
3 0 66 0.844 0.835 0.960 0.242 1 0.960 0.172 0.825 0.960 0.212 37189 19877 3595 71339 0.865 0.298 0.043 True True False
3 1 66 0.843 0.835 0.960 0.253 1 0.960 0.172 0.824 0.960 0.232 37120 19875 3661 71344 0.864 0.296 0.044 True True False
3 2 66 0.842 0.835 0.960 0.253 1 0.960 0.172 0.824 0.960 0.232 37432 19629 3896 71043 0.865 0.297 0.042 True True False
4 0 185 0.837 0.850 0.960 0.222 1 0.949 0.162 0.827 0.960 0.202 106208 66651 9127 188014 0.934 0.351 0.110 True True False
4 1 185 0.837 0.850 0.960 0.212 1 0.949 0.162 0.827 0.960 0.202 107169 65821 9588 187422 0.935 0.352 0.112 True True False
4 2 185 0.836 0.850 0.970 0.182 1 0.949 0.162 0.826 0.960 0.192 103182 69983 7795 189040 0.936 0.352 0.113 True True False
5 0 159 0.777 0.797 0.949 0.253 1 0.929 0.162 0.763 0.939 0.242 75471 70278 9459 162792 0.917 0.362 0.184 True True False
5 1 159 0.778 0.797 0.960 0.212 1 0.929 0.162 0.765 0.939 0.222 73472 72339 8151 164038 0.917 0.358 0.183 True True False
5 2 159 0.778 0.797 0.960 0.202 1 0.929 0.162 0.764 0.939 0.232 73449 72525 8041 163985 0.918 0.360 0.183 True True False
6 0 50 0.669 0.599 0.949 0.323 1 0.929 0.192 0.633 0.939 0.273 15926 19573 4948 59553 0.710 0.586 0.355 False True False
6 1 50 0.674 0.599 0.949 0.283 1 0.929 0.192 0.633 0.929 0.273 15904 19503 4736 59857 0.708 0.589 0.357 False True False
6 2 50 0.670 0.599 0.949 0.303 1 0.929 0.192 0.632 0.929 0.303 16045 19543 5008 59404 0.712 0.585 0.354 False True False
7 0 55 0.642 0.604 0.970 0.061 1 0.929 0.152 0.585 0.929 0.242 12950 24181 2975 69894 0.675 0.613 0.487 True True False
7 1 55 0.647 0.604 0.960 0.303 1 0.929 0.152 0.586 0.929 0.303 14442 22686 4197 68675 0.675 0.615 0.493 True True False
7 2 55 0.644 0.604 0.970 0.000 5 0.929 0.152 0.583 0.929 0.263 12767 24125 2805 70303 0.671 0.612 0.490 True True False
8 0 66 0.555 0.544 0.970 0.081 1 0.929 0.030 0.503 0.929 0.051 12828 25317 6551 87304 0.578 0.676 0.521 True True False
8 1 66 0.554 0.544 0.970 0.071 1 0.929 0.030 0.507 0.929 0.071 12764 25380 6500 87356 0.578 0.675 0.524 True True False
8 2 66 0.553 0.544 0.970 0.091 1 0.929 0.030 0.507 0.929 0.061 12662 25429 6443 87466 0.577 0.674 0.527 True True False
9 0 179 0.781 0.795 0.960 0.343 1 0.949 0.222 0.762 0.949 0.374 83574 57139 15093 202194 0.786 0.392 0.075 True True False
9 1 179 0.777 0.795 0.960 0.404 1 0.949 0.222 0.759 0.949 0.364 80547 60194 13805 203454 0.786 0.392 0.075 True True False
9 2 179 0.781 0.795 0.960 0.364 1 0.949 0.222 0.761 0.949 0.374 82062 58706 14144 203088 0.786 0.394 0.074 True True False
10 0 179 0.785 0.795 0.960 0.364 1 0.949 0.222 0.765 0.949 0.364 82923 58364 13717 202996 0.789 0.387 0.071 True True False
10 1 179 0.780 0.795 0.960 0.343 1 0.949 0.222 0.761 0.949 0.374 82432 58220 14517 202831 0.786 0.392 0.073 True True False
10 2 179 0.782 0.795 0.960 0.313 1 0.949 0.222 0.762 0.949 0.343 83069 57620 14508 202803 0.786 0.386 0.073 True True False
11 0 99 0.801 0.775 0.960 0.323 1 0.939 0.232 0.794 0.949 0.323 49900 31721 7590 108789 0.824 0.378 0.063 True True False
11 1 99 0.797 0.775 0.960 0.333 1 0.939 0.232 0.792 0.949 0.333 48798 32714 7323 109165 0.823 0.378 0.063 True True False
11 2 99 0.793 0.775 0.960 0.354 1 0.939 0.232 0.789 0.949 0.333 48314 33154 7453 109079 0.823 0.388 0.065 True True False
12 0 114 0.851 0.870 0.960 0.263 1 0.960 0.152 0.834 0.960 0.232 64749 33724 5781 123746 0.864 0.380 0.055 True True False
12 1 114 0.849 0.870 0.960 0.253 1 0.960 0.152 0.833 0.960 0.222 65196 33397 6176 123231 0.865 0.380 0.054 True True False
12 2 114 0.848 0.870 0.960 0.263 1 0.960 0.152 0.831 0.960 0.232 64595 33874 6049 123482 0.864 0.378 0.055 True True False
13 0 128 0.808 0.860 0.970 0.222 1 0.960 0.172 0.794 0.960 0.263 59773 47419 5843 142965 0.837 0.396 0.089 True True False
13 1 128 0.809 0.860 0.960 0.303 1 0.960 0.172 0.794 0.960 0.253 67208 40047 9860 138885 0.838 0.393 0.090 True True False
13 2 128 0.806 0.860 0.960 0.313 1 0.960 0.172 0.791 0.960 0.263 66822 40335 10047 138796 0.837 0.392 0.090 True True False
14 0 133 0.586 0.485 0.939 0.131 1 0.869 0.111 0.519 0.859 0.162 29544 52352 13008 171096 0.616 0.632 0.580 False True False
14 1 133 0.586 0.485 0.949 0.040 1 0.869 0.111 0.517 0.838 0.192 27654 54146 10905 173295 0.615 0.636 0.588 False True False
14 2 133 0.586 0.485 0.939 0.101 1 0.869 0.111 0.520 0.869 0.101 29419 52586 12852 171143 0.617 0.635 0.590 False True True
15 0 133 0.578 0.520 0.949 0.212 1 0.879 0.091 0.513 0.869 0.212 27794 52791 12119 173296 0.606 0.613 0.552 False True False
15 1 133 0.578 0.520 0.949 0.172 1 0.879 0.091 0.514 0.869 0.172 27797 52798 12218 173187 0.606 0.618 0.556 False True False
15 2 133 0.575 0.520 0.949 0.172 1 0.879 0.091 0.510 0.869 0.172 27783 52550 12532 173135 0.604 0.618 0.558 False True False
16 0 59 0.530 0.525 0.960 0.000 1 0.545 0.424 0.515 0.929 0.111 13656 20346 10024 73974 0.576 0.600 0.428 True True False
16 1 59 0.533 0.525 0.960 0.000 1 0.545 0.424 0.514 0.929 0.081 13784 20162 10047 74007 0.575 0.603 0.424 True True False
16 2 59 0.532 0.525 0.960 0.000 1 0.545 0.424 0.513 0.929 0.091 13628 20339 9906 74127 0.576 0.600 0.433 True True False
17 0 50 0.481 0.537 0.960 0.212 1 0.889 0.232 0.468 0.919 0.172 10288 17987 9397 62328 0.566 0.485 0.336 False False True
17 1 50 0.478 0.537 0.960 0.121 1 0.889 0.232 0.464 0.899 0.131 10218 17874 9454 62454 0.562 0.481 0.337 False False False
17 2 50 0.477 0.537 0.960 0.091 1 0.889 0.232 0.466 0.909 0.232 10416 17626 9864 62094 0.561 0.473 0.330 False False False
18 0 101 0.731 0.723 0.960 0.313 1 0.949 0.182 0.693 0.960 0.242 39422 34174 9582 118822 0.729 0.336 0.104 True True False
18 1 101 0.729 0.723 0.960 0.263 1 0.949 0.182 0.693 0.960 0.242 40163 33362 10310 118165 0.728 0.340 0.104 True True False
18 2 101 0.730 0.723 0.960 0.293 1 0.949 0.182 0.694 0.949 0.293 39527 34148 9738 118587 0.729 0.340 0.104 True True False
19 0 87 0.637 0.610 0.939 0.222 1 0.909 0.182 0.599 0.899 0.232 26491 35775 9959 101775 0.716 0.499 0.370 True True False
19 1 87 0.634 0.610 0.939 0.202 1 0.909 0.182 0.595 0.899 0.263 26203 35813 9952 102032 0.713 0.507 0.373 True True False
19 2 87 0.635 0.610 0.939 0.242 1 0.909 0.182 0.597 0.909 0.253 25855 36256 9532 102357 0.714 0.506 0.372 True True False
20 0 59 0.534 0.577 0.960 0.000 1 0.879 0.121 0.478 0.899 0.000 11027 22307 6434 78232 0.565 0.597 0.474 True False False
20 1 59 0.523 0.577 0.960 0.000 1 0.879 0.121 0.474 0.879 0.000 10763 22556 6654 78027 0.565 0.600 0.474 True False False
20 2 59 0.534 0.577 0.960 0.000 1 0.879 0.121 0.477 0.899 0.000 10940 22577 6268 78215 0.568 0.603 0.470 True False False
Recordings (19) with max F within 10% of Table 1, per noise region; Methods reading, code reading, and with Q_test and Q_GT also within 10%:
maxF maxF_code quantiles
region
0 15 17 1
1 15 17 0
2 15 17 1
In [18]:
fig, ax = plt.subplots(figsize=(5, 5))
for g, m in zip(range(3), "o^s"):
x = recs[recs.region == g]
ax.plot(x.table_maxF, x.maxF, m, label=f"region {g}, Methods", alpha=0.7)
ax.plot(x.table_maxF, x.maxF_code, m, mfc="none", color="0.5", alpha=0.7)
ax.plot([0.3, 1], [0.3, 1], "k-", lw=0.8)
ax.set_xlabel("Table 1 max F")
ax.set_ylabel("max F here (open: code reading)")
ax.set_xlim(0.3, 1)
ax.set_ylim(0.3, 1)
ax.legend(fontsize=8)
print("Plotted max F range:", grid[["maxF", "maxF_code"]].min().min(), "to", grid[["maxF", "maxF_code"]].max().max())
Plotted max F range: 0.464 to 0.851
In [19]:
# Share of surrogates whose best unit is their own unit (Methods 2.8: "in rare cases, the best-R of
# the correlation test with ground-truth unit was not from the ground-truth unit"), with the median
# over units of the median GT-R and median test-R
hit = grid.pivot_table(index="row", columns="region", values=["hit_fraction", "gt_r_median", "test_r_median"])
hit[("units", "")] = grid[grid.region == 0].set_index("row").units
hit[("maxF_methods_minus_code", 0)] = (grid[grid.region == 0].set_index("row").maxF - grid[grid.region == 0].set_index("row").maxF_code).round(3)
print(hit.round(3).to_string())
print("Misidentified share (1 - hit fraction) over rows and regions:", round(1 - grid.hit_fraction.max(), 3), "to",
round(1 - grid.hit_fraction.min(), 3))
x = grid[grid.region == 0]
print("Spearman, hit fraction against units:", round(stats.spearmanr(x.hit_fraction, x.units).statistic, 3),
"| against max F (Methods):", round(stats.spearmanr(x.hit_fraction, x.maxF).statistic, 3),
"| against max F Methods minus code:", round(stats.spearmanr(x.hit_fraction, x.maxF - x.maxF_code).statistic, 3))
# Rows 16 and 18 hold a unit with no spikes; the grid again without it, region 0
empty_var = []
with ProcessPoolExecutor(4) as pool:
for r, d in sorted(spine_rows.items()):
s, c = sessions[d["file"]], chosen[r]
live = [u for u in range(len(s["trains"])) if len(s["trains"][u]) > 0]
if len(live) == len(s["trains"]):
continue
s2 = dict(s, trains=[s["trains"][u] for u in live])
R = sc.r_tensor_parallel(s2, s["U"][live], c["keep"], pr_main, estimates[r][0], seed=223000 + 10 * r, pool=pool)
fg = sc.fscore_grid(R)
g0 = grid[(grid.row == r) & (grid.region == 0)].iloc[0]
empty_var.append({"row": r, "units": len(s["trains"]), "units_with_spikes": len(live),
"maxF": g0.maxF, "maxF_without_empty": round(sc.grid_optimum(fg["methods"][0])[0], 3),
"maxF_code": g0.maxF_code, "maxF_code_without_empty": round(sc.grid_optimum(fg["code"][0])[0], 3),
"table_maxF": table.max_F[r]})
del R
print()
print(pd.DataFrame(empty_var).to_string(index=False))
gt_r_median hit_fraction test_r_median units maxF_methods_minus_code region 0 1 2 0 1 2 0 1 2 0 row 1 0.380 0.375 0.378 0.901 0.901 0.902 0.221 0.218 0.219 78 0.019 2 0.358 0.363 0.356 0.920 0.920 0.920 0.164 0.169 0.166 132 0.009 3 0.298 0.296 0.297 0.865 0.864 0.865 0.043 0.044 0.042 66 0.019 4 0.351 0.352 0.352 0.934 0.935 0.936 0.110 0.112 0.113 185 0.010 5 0.362 0.358 0.360 0.917 0.917 0.918 0.184 0.183 0.183 159 0.014 6 0.586 0.589 0.585 0.710 0.708 0.712 0.355 0.357 0.354 50 0.036 7 0.613 0.615 0.612 0.675 0.675 0.671 0.487 0.493 0.490 55 0.057 8 0.676 0.675 0.674 0.578 0.578 0.577 0.521 0.524 0.527 66 0.052 9 0.392 0.392 0.394 0.786 0.786 0.786 0.075 0.075 0.074 179 0.019 10 0.387 0.392 0.386 0.789 0.786 0.786 0.071 0.073 0.073 179 0.020 11 0.378 0.378 0.388 0.824 0.823 0.823 0.063 0.063 0.065 99 0.007 12 0.380 0.380 0.378 0.864 0.865 0.864 0.055 0.054 0.055 114 0.017 13 0.396 0.393 0.392 0.837 0.838 0.837 0.089 0.090 0.090 128 0.014 14 0.632 0.636 0.635 0.616 0.615 0.617 0.580 0.588 0.590 133 0.067 15 0.613 0.618 0.618 0.606 0.606 0.604 0.552 0.556 0.558 133 0.065 16 0.600 0.603 0.600 0.576 0.575 0.576 0.428 0.424 0.433 59 0.015 17 0.485 0.481 0.473 0.566 0.562 0.561 0.336 0.337 0.330 50 0.013 18 0.336 0.340 0.340 0.729 0.728 0.729 0.104 0.104 0.104 101 0.038 19 0.499 0.507 0.506 0.716 0.713 0.714 0.370 0.373 0.372 87 0.038 20 0.597 0.600 0.603 0.565 0.565 0.568 0.474 0.474 0.470 59 0.056 Misidentified share (1 - hit fraction) over rows and regions: 0.064 to 0.439 Spearman, hit fraction against units: 0.551 | against max F (Methods): 0.884 | against max F Methods minus code: -0.578
row units units_with_spikes maxF maxF_without_empty maxF_code maxF_code_without_empty table_maxF 16 59 58 0.530 0.528 0.515 0.518 0.525 18 101 100 0.731 0.732 0.693 0.695 0.723
Result: partially reproduced. Under the Methods' classification (the claims file's main reading), the maximum F is within 10% of Table 1 for 15 of the 19 recordings with every one of the three noise regions; it misses for rows 6, 14 and 15 (higher here: 0.669 to 0.674, 0.586 and 0.575 to 0.578 against 0.599, 0.485 and 0.520) and row 17 (lower: 0.477 to 0.481 against 0.537). With the code's classification (a variant), 17 of 19 match. The optimal quantiles are not those of Table 1: Q_test is 0.939 to 0.970 here for every recording, and the pair lies within 10% of Table 1's for 0 or 1 of 19 recordings. Table 1's Q_test spans 0.545 to 0.960. Each maximum is unique except row 7 with region 2 (5 maxima). Rows 9 and 10, run with each spine's noise estimate, give 0.777 to 0.785 against 0.795.
The share of ground-truth surrogates whose best unit is not their own unit, which Methods 2.8 calls "rare", is 0.064 to 0.439 here: 0.064 to 0.126 for rows 1 to 5, and 0.383 to 0.439 for rows 8, 14, 15, 16, 17 and 20. It is nearly the same in the three noise regions, and recordings with more misidentified surrogates have lower maximum F (Spearman 0.89 between the hit fraction and max F) and a larger gap between the two classifications (up to 0.067, row 14). These are the surrogates the Methods and the code classify differently, so the share decides how far the two readings differ. Whether it is as large in the paper cannot be told from the paper; if it was smaller there, the surrogates here are less distinct between units than the paper's. Rows 16 and 18 keep a unit with no spikes; removing it changes max F by at most 0.003 (0.530 to 0.528 and 0.731 to 0.731 under the Methods' classification, 0.515 to 0.518 and 0.693 to 0.695 under the code's).
Claim mean-max-fscore¶
The simulated ground-truth test gives large F-scores across experiments, a mean maximum F-score of 0.692 ± 0.136 over 19 recordings with 1'000 surrogates per unit.
In [20]:
mf = []
for g in range(3):
for col, name in [("maxF", "Methods"), ("maxF_code", "code")]:
x19 = recs[recs.region == g][col].to_numpy()
x20 = grid[grid.region == g][col].to_numpy()
x19b = grid[(grid.region == g) & (grid.row != 9)][col].to_numpy()
mf.append({"region": g, "reading": name, "mean_19": round(x19.mean(), 3), "sd_19": round(x19.std(), 3),
"mean_19_row10": round(x19b.mean(), 3), "sd_19_row10": round(x19b.std(), 3),
"mean_20": round(x20.mean(), 3), "sd_20": round(x20.std(), 3)})
mf = pd.DataFrame(mf)
mf["match"] = ((mf.mean_19 - 0.692).abs() <= 0.0692) & ((mf.sd_19 - 0.136).abs() <= 0.0136)
print(mf.to_string(index=False))
print("Table 1: mean 0.692, SD 0.136 over 19 recordings (ddof 0)")
region reading mean_19 sd_19 mean_19_row10 sd_19_row10 mean_20 sd_20 match
0 Methods 0.695 0.120 0.696 0.120 0.700 0.119 False
0 code 0.666 0.133 0.666 0.133 0.671 0.131 True
1 Methods 0.695 0.120 0.695 0.120 0.699 0.119 False
1 code 0.665 0.133 0.665 0.133 0.670 0.131 True
2 Methods 0.694 0.120 0.694 0.120 0.699 0.118 False
2 code 0.665 0.132 0.665 0.133 0.670 0.131 True
Table 1: mean 0.692, SD 0.136 over 19 recordings (ddof 0)
Result: partially reproduced. Under the Methods' classification the mean maximum F over 19 recordings is 0.696 with noise region 0 (expected 0.692 ± 0.069) and the SD 0.121 (ddof 0), just below the tolerance of 0.136 ± 0.0136; over the three regions the mean is 0.695 to 0.696 and the SD 0.120 to 0.121. With the code's classification the mean is 0.665 to 0.666 and the SD 0.134, both within tolerance. The F-scores are large across recordings, as the paper says; the spread under the main reading is slightly smaller.
The two classifications differ by the misidentified surrogates of claim table1-optimal-thresholds (6% to 44% of surrogates), so whether the SD falls inside the tolerance depends on that share.
Claim acceptance-9-of-20, reading (b): the optimal quantiles found here¶
In [21]:
acc_b = []
for (r, g), (gt, te) in dist.items():
gr = grid[(grid.row == r) & (grid.region == g)].iloc[0]
decisions = [decide(gt, te, qg, qt, chosen[r]["best_R"]) for qt, qg in gr.all_max]
gthr, tthr, ok = decisions[0]
acc_b.append({"row": r, "region": g, "Q_test": gr.Q_test, "Q_GT": gr.Q_GT, "best_R": round(chosen[r]["best_R"], 3),
"GT_thr": round(gthr, 3), "GT_thr_table": table.GT_thr[r], "test_thr": round(tthr, 3),
"test_thr_table": table.test_thr[r], "accept": ok, "accept_table": table.accept[r],
"same_decision_all_maxima": len({x[2] for x in decisions}) == 1})
acc_b = pd.DataFrame(acc_b)
acc_b["agree"] = acc_b.accept == acc_b.accept_table
print(acc_b[acc_b.region == 0].drop(columns="region").to_string(index=False))
print("\nPer noise region, reading (b): accepted, decisions agreeing with Table 1:")
print(acc_b.groupby("region").agg(accepted=("accept", "sum"), agree=("agree", "sum"),
stable=("same_decision_all_maxima", "all")).to_string())
print("\nDecision per spine and region (1 = accept):")
print(acc_b.pivot_table(index="row", columns="region", values="accept").astype(int).assign(table=table.accept.astype(int)).T.to_string())
row Q_test Q_GT best_R GT_thr GT_thr_table test_thr test_thr_table accept accept_table same_decision_all_maxima agree
1 0.970 0.101 0.522 0.198 0.214 0.311 0.296 True True True True
2 0.970 0.242 0.515 0.368 0.291 0.350 0.271 True True True True
3 0.960 0.242 0.213 0.134 0.112 0.264 0.190 False True True False
4 0.960 0.222 0.643 0.316 0.311 0.489 0.489 True True True True
5 0.949 0.253 0.430 0.385 0.358 0.566 0.520 False False True True
6 0.949 0.323 0.653 0.477 0.448 0.616 0.619 True True True True
7 0.970 0.061 0.344 0.322 0.503 0.692 0.585 False False True True
8 0.970 0.081 0.363 0.290 0.227 0.458 0.475 False False True True
9 0.960 0.343 0.673 0.548 0.507 0.660 0.639 True True True True
10 0.960 0.364 0.762 0.540 0.507 0.644 0.639 True True True True
11 0.960 0.323 0.137 0.304 0.354 0.279 0.351 False False True True
12 0.960 0.263 0.585 0.379 0.353 0.587 0.544 False True True False
13 0.970 0.222 0.143 0.270 0.267 0.109 0.068 False False True True
14 0.939 0.131 0.354 0.430 0.456 0.378 0.435 False False True True
15 0.949 0.212 0.395 0.595 0.436 0.639 0.604 False False True True
16 0.960 0.000 0.537 0.093 0.383 0.567 0.666 False True True False
17 0.960 0.212 0.627 0.346 0.500 0.767 0.707 False False True True
18 0.960 0.313 0.311 0.241 0.214 0.338 0.324 False False True True
19 0.939 0.222 0.222 0.468 0.233 0.206 0.298 False False True True
20 0.960 0.000 0.293 0.030 0.414 0.593 0.364 False False True True
Per noise region, reading (b): accepted, decisions agreeing with Table 1:
accepted agree stable
region
0 6 17 True
1 7 18 True
2 6 17 True
Decision per spine and region (1 = accept):
row 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
region
0 1 1 0 1 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0
1 1 1 0 1 0 1 0 0 1 1 0 1 0 0 0 0 0 0 0 0
2 1 1 0 1 0 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0
table 1 1 1 1 0 1 0 0 1 1 0 1 0 0 0 1 0 0 0 0
Result: partially reproduced. Every reading and noise region gives 17 or 18 decisions that agree with Table 1, the claims file's band for a partial reproduction, so the outcome does not depend on the noise region. With Table 1's quantiles (reading a), 7 spines are accepted with each region, and 18 decisions agree; rows 3 and 16, accepted in Table 1, are rejected here, in each of five seeds and regions. With the optimal quantiles found here (reading b), 7, 7 and 6 are accepted (18, 18 and 17 agree). The stored release probabilities and a Bernoulli draw per spike also give 7 accepted and 18 agreeing. GT thresholds are within 10% of Table 1 for 12 or 13 spines and test thresholds for 10 or 11. Rows 9 and 10, which Table 1 gives identical thresholds (0.507 and 0.639), have different ones here, as each spine has its own noise estimate (0.497 and 0.649 against 0.483 and 0.632 with region 0). The two spines share their best (115) and second-best (44) units, and with row 9's noise estimate row 10's thresholds equal row 9's in every region, so Table 1's identical values fit one noise estimate per recording.
Of the two mismatches, row 3 has best-R 0.221 (Table 1, 0.217) and a GT threshold near Table 1's (0.112 to 0.121 against 0.112), but its test-Rs (of the second-best unit, 23) reach 0.245 to 0.254 at Table 1's Q_test of 0.96, against a threshold of 0.190 there: the upper tail of the test distribution is wider here. Row 16 has GT-Rs with median 0.71, so the GT threshold at Table 1's Q_GT of 0.424 is 0.683 to 0.684 against 0.383, and its best-R is 0.537, not 0.676; its test threshold (0.430 to 0.437 against 0.666) is lower than Table 1's. Both distributions are nearly the same in the three regions, so neither mismatch comes from the noise region.
Claim best-r-accepted-vs-rejected¶
Accepted connections have larger best-R values than rejected ones (0.584 ± 0.149 against 0.331 ± 0.124; Mann-Whitney U = 12.0, n1 = 9, n2 = 11, P < 0.05).
In [22]:
mw = []
for name, acc in [("a", acc_a), ("b", acc_b)]:
for g in range(3):
x = acc[acc.region == g]
A, B = x[x.accept].best_R.to_numpy(), x[~x.accept].best_R.to_numpy()
if len(A) and len(B):
u = stats.mannwhitneyu(A, B, alternative="two-sided")
u1 = stats.mannwhitneyu(A, B, alternative="greater")
mw.append({"reading": name, "region": g, "n_acc": len(A), "n_rej": len(B), "mean_acc": round(A.mean(), 3),
"sd_acc": round(A.std(), 3), "mean_rej": round(B.mean(), 3), "sd_rej": round(B.std(), 3),
"U_acc": u.statistic, "U_min": min(u.statistic, len(A) * len(B) - u.statistic),
"p_two_sided": round(u.pvalue, 4), "p_one_sided": round(u1.pvalue, 4)})
else:
mw.append({"reading": name, "region": g, "n_acc": len(A), "n_rej": len(B)})
mw = pd.DataFrame(mw)
print(mw.to_string(index=False))
print("Paper: 0.584 ± 0.149 against 0.331 ± 0.124, U = 12.0, n1 = 9, n2 = 11, P < 0.05")
reading region n_acc n_rej mean_acc sd_acc mean_rej sd_rej U_acc U_min p_two_sided p_one_sided
a 0 6 14 0.628 0.086 0.354 0.147 78.0 6.0 0.0015 0.0008
a 1 7 13 0.622 0.081 0.336 0.138 86.0 5.0 0.0005 0.0002
a 2 7 13 0.622 0.081 0.336 0.138 86.0 5.0 0.0005 0.0002
b 0 6 14 0.628 0.086 0.354 0.147 78.0 6.0 0.0015 0.0008
b 1 7 13 0.622 0.081 0.336 0.138 86.0 5.0 0.0005 0.0002
b 2 6 14 0.613 0.085 0.360 0.158 74.0 10.0 0.0064 0.0032
Paper: 0.584 ± 0.149 against 0.331 ± 0.124, U = 12.0, n1 = 9, n2 = 11, P < 0.05
Result: partially reproduced. With this attempt's classification (7 accepted, 13 rejected, reading (a), any noise region), accepted spines have larger best-R: 0.622 ± 0.081 against 0.337 ± 0.137, Mann-Whitney U = 5 (smaller U; 86 for the accepted group), two-sided p = 0.0005 (one-sided 0.0002). Both means are within tolerance of 0.584 and 0.331 and the test gives P < 0.05 with the accepted group higher, but U is 5, not 12 ± 1.2; with 7 and 13 spines in place of 9 and 11 the U values are not comparable. Reading (b) with region 2 gives U = 10 (6 and 14 spines, two-sided p = 0.0064).
Claim best-r-not-predictive¶
The best-R alone does not predict the outcome of the acceptance test; connections with R as high as about 0.6 were rejected, while an R of about 0.2 was enough for acceptance in another case.
In [23]:
npred = []
for name, acc in [("a", acc_a), ("b", acc_b)]:
for g in range(3):
x = acc[acc.region == g]
mx_rej = x[~x.accept].best_R.max() if (~x.accept).any() else np.nan
mn_acc = x[x.accept].best_R.min() if x.accept.any() else np.nan
npred.append({"reading": name, "region": g, "max_rejected": mx_rej, "row_max_rejected": x[~x.accept].row.iloc[x[~x.accept].best_R.argmax()] if (~x.accept).any() else None,
"min_accepted": mn_acc, "row_min_accepted": x[x.accept].row.iloc[x[x.accept].best_R.argmin()] if x.accept.any() else None,
"overlap": bool(mx_rej > mn_acc)})
npred = pd.DataFrame(npred)
npred["max_rej_within"] = (npred.max_rejected - 0.6).abs() <= 0.06
npred["min_acc_within"] = (npred.min_accepted - 0.2).abs() <= 0.02
print(npred.to_string(index=False))
reading region max_rejected row_max_rejected min_accepted row_min_accepted overlap max_rej_within min_acc_within
a 0 0.627 17 0.515 2 True True False
a 1 0.627 17 0.515 2 True True False
a 2 0.627 17 0.515 2 True True False
b 0 0.627 17 0.515 2 True True False
b 1 0.627 17 0.515 2 True True False
b 2 0.673 9 0.515 2 True False False
Result: partially reproduced. The largest best-R of a rejected spine is 0.627 (row 17), within tolerance of about 0.6, and it exceeds the smallest accepted best-R, so no cut on best-R separates the groups. The smallest accepted best-R is 0.515 (row 2), not about 0.2: row 3, the accepted spine with R near 0.2 in Table 1, is rejected here. Reading (b) with region 2 gives a largest rejected best-R of 0.673 (row 9).
Claims fig5a-rejection-scenarios and fig5a-threshold-order¶
Connections were rejected in three scenarios, best-R below the test threshold only, below the GT threshold only, or below both, which shows the advantage of having two thresholds.
For some recordings the test threshold is larger than the GT threshold, and for others the reverse.
In [24]:
sc_rows = []
for name, acc in [("a", acc_a), ("b", acc_b)]:
for g in range(3):
x = acc[acc.region == g]
rej = x[~x.accept]
sc_rows.append({"reading": name, "region": g, "rejected": len(rej),
"below_test_only": int(((rej.best_R <= rej.test_thr) & (rej.best_R > rej.GT_thr)).sum()),
"below_GT_only": int(((rej.best_R <= rej.GT_thr) & (rej.best_R > rej.test_thr)).sum()),
"below_both": int(((rej.best_R <= rej.GT_thr) & (rej.best_R <= rej.test_thr)).sum()),
"test_above_GT": int((x.test_thr > x.GT_thr).sum()), "GT_above_test": int((x.GT_thr > x.test_thr).sum())})
print(pd.DataFrame(sc_rows).to_string(index=False))
print("Table 1: 4 below test only, 1 below GT only, 6 below both; test above GT for 15 spines, below for 5")
reading region rejected below_test_only below_GT_only below_both test_above_GT GT_above_test
a 0 14 6 4 4 14 6
a 1 13 5 4 4 16 4
a 2 13 5 4 4 15 4
b 0 14 9 2 3 15 5
b 1 13 8 2 3 18 2
b 2 14 9 2 3 15 5
Table 1: 4 below test only, 1 below GT only, 6 below both; test above GT for 15 spines, below for 5
Result: both claims reproduced. All three rejection scenarios occur among the rejected spines under every reading and region: with Table 1's quantiles and region 0, 5 below the test threshold only, 4 below the GT threshold only and 4 below both (Table 1: 4, 1 and 6); the other readings and regions give 5 to 9, 2 to 4 and 3 to 4. The test threshold is above the GT threshold for 14 spines and below it for 6 with region 0 and Table 1's quantiles (Table 1: 15 and 5); the other readings and regions give 15 to 18 and 2 to 5.
Claim table1-network-measures¶
Each recording's spike contrast (0.17 to 0.69) and ratio of spikes outside to inside network bursts (0.060 to 2.003) have the values given in Table 1.
Every combination of the open details: (a) all sorted units or those kept after the hyperactive exclusion, (b) for OUT/IN, burst periods with or without the extensions, (c) the whole MEA recording or the span of the imaging frames. The bursts are those of the carried-forward reading.
In [25]:
net = []
for f, rs in mapping.items():
s = sessions[f]
c = chosen[rs[0]]
drop = set(sc.hyperactive(s, c["strategy"]).tolist())
for units_set in ["all", "clean"]:
idx = list(range(len(s["trains"]))) if units_set == "all" else [u for u in range(len(s["trains"])) if u not in drop]
for span in ["whole", "imaging"]:
t0_, t1_ = (0.0, s["mea_duration"]) if span == "whole" else (s["frame_times"][0], s["frame_times"][-1])
scv = spike_contrast_of(s, t0_, t1_, idx)
for ext in ["extended", "detected"]:
a, b = (c["starts"] - 0.1, c["ends"] + 0.5) if ext == "extended" else (c["starts"], c["ends"])
oi = sc.out_in_ratio([s["trains"][u] for u in idx], a, b, t0_, t1_)
net.append({"row": rs[0], "units": units_set, "span": span, "bursts": ext, "spike_contrast": scv,
"out_in": oi, "sc_table": table.spike_contrast[rs[0]], "oi_table": table.out_in[rs[0]]})
net = pd.DataFrame(net)
net["sc_ok"] = (net.spike_contrast - net.sc_table).abs() <= 0.1 * net.sc_table
net["oi_ok"] = (net.out_in - net.oi_table).abs() <= 0.1 * net.oi_table
net["both"] = net.sc_ok & net.oi_ok
combo = net.groupby(["units", "span", "bursts"]).agg(sc_match=("sc_ok", "sum"), oi_match=("oi_ok", "sum"),
both_match=("both", "sum"), recordings=("row", "count"))
print(combo.to_string())
best_combo = combo.both_match.idxmax()
print("\nCombination with most recordings matching both:", best_combo)
bc = net[(net.units == best_combo[0]) & (net.span == best_combo[1]) & (net.bursts == best_combo[2])].sort_values("row")
print(bc[["row", "spike_contrast", "sc_table", "sc_ok", "out_in", "oi_table", "oi_ok"]].round(3).to_string(index=False))
print("Spike contrast range:", round(bc.spike_contrast.min(), 3), "to", round(bc.spike_contrast.max(), 3),
"| OUT/IN range:", round(bc.out_in.min(), 3), "to", round(bc.out_in.max(), 3))
# Variant: OUT/IN with the bursts of every reading of p_thresh and strategy (5-point smoothing),
# all units, whole recording, extended bursts
oi_var = []
for f, rs in mapping.items():
s = sessions[f]
for p_, st_ in itertools.product(sc.P_THRESH, sc.STRATEGIES):
a, b, _ = sc.detect_bursts(s, p_, st_, "span5")
oi_var.append({"row": rs[0], "reading": f"{p_}{st_}", "out_in": sc.out_in_ratio(s["trains"], a - 0.1, b + 0.5, 0.0, s["mea_duration"])})
oi_var = pd.DataFrame(oi_var).pivot_table(index="row", columns="reading", values="out_in")
oi_var["table"] = [table.out_in[r] for r in oi_var.index]
oi_var["any_within_10pct"] = [any(abs(oi_var.loc[r, c] - oi_var.table[r]) <= 0.1 * oi_var.table[r]
for c in oi_var.columns[:-1]) for r in oi_var.index]
print("\nVariant, OUT/IN under each burst reading:")
print(oi_var.round(3).to_string())
print("Recordings with OUT/IN within 10% of Table 1 under some burst reading:", int(oi_var.any_within_10pct.sum()), "of", len(oi_var))
sc_match oi_match both_match recordings
units span bursts
all imaging detected 19 0 0 19
extended 19 7 7 19
whole detected 19 0 0 19
extended 19 7 7 19
clean imaging detected 16 0 0 19
extended 16 6 6 19
whole detected 16 0 0 19
extended 16 6 6 19
Combination with most recordings matching both: ('all', 'imaging', 'extended')
row spike_contrast sc_table sc_ok out_in oi_table oi_ok
1 0.377 0.38 True 0.737 0.733 True
2 0.171 0.17 True 2.194 2.003 True
3 0.170 0.17 True 1.755 1.717 True
4 0.227 0.23 True 0.993 0.588 False
5 0.260 0.26 True 0.809 0.318 False
6 0.401 0.41 True 0.762 0.189 False
7 0.483 0.48 True 0.428 0.269 False
8 0.619 0.62 True 0.419 0.060 False
9 0.275 0.28 True 0.443 0.425 True
11 0.205 0.20 True 1.171 0.783 False
12 0.294 0.29 True 0.593 0.719 False
13 0.364 0.36 True 0.773 0.425 False
14 0.624 0.64 True 0.272 0.289 True
15 0.692 0.69 True 0.221 0.222 True
16 0.428 0.43 True 0.564 0.171 False
17 0.303 0.31 True 1.119 0.534 False
18 0.243 0.25 True 0.653 0.464 False
19 0.506 0.50 True 0.297 0.272 True
20 0.541 0.53 True 0.213 0.112 False
Spike contrast range: 0.17 to 0.692 | OUT/IN range: 0.213 to 2.194
Variant, OUT/IN under each burst reading: reading 96hz 96pc 97hz 97pc 98hz 98pc 99hz 99pc table any_within_10pct row 1 0.595 0.595 0.737 0.737 0.932 0.932 1.510 1.510 0.733 True 2 1.993 1.841 2.413 2.195 2.778 2.926 3.293 3.301 2.003 True 3 1.150 1.150 1.755 1.755 2.294 2.294 3.887 3.887 1.717 True 4 0.791 0.747 0.991 0.993 1.171 1.175 1.750 1.771 0.588 False 5 0.683 0.696 0.809 0.816 1.076 1.076 2.027 2.055 0.318 False 6 0.329 0.267 0.373 0.308 0.448 0.445 0.762 0.728 0.189 False 7 0.319 0.314 0.340 0.349 0.428 0.428 0.687 0.715 0.269 False 8 0.139 0.139 0.151 0.150 0.169 0.180 0.419 0.417 0.060 False 9 0.444 0.444 0.460 0.460 0.473 0.473 0.744 0.744 0.425 True 11 0.752 0.928 0.752 0.928 0.895 1.054 1.584 1.171 0.783 True 12 0.593 0.600 0.674 0.717 0.776 0.785 0.951 0.960 0.719 True 13 0.545 0.529 0.556 0.553 0.600 0.562 0.753 0.773 0.425 False 14 0.272 0.272 0.286 0.286 0.324 0.324 0.414 0.414 0.289 True 15 0.207 0.207 0.222 0.222 0.268 0.255 0.427 0.397 0.222 True 16 0.297 0.297 0.326 0.318 0.372 0.391 0.561 0.565 0.171 False 17 0.538 0.538 0.589 0.589 0.823 0.823 1.119 1.119 0.534 True 18 0.600 0.652 0.712 0.757 0.940 0.965 1.724 1.780 0.464 False 19 0.265 0.265 0.297 0.297 0.354 0.354 0.519 0.519 0.272 True 20 0.212 0.213 0.229 0.215 0.275 0.265 0.565 0.624 0.112 False Recordings with OUT/IN within 10% of Table 1 under some burst reading: 10 of 19
In [26]:
# Burst merging (Methods 2.6: bursts closer than t_ibi "were merged into a single larger burst").
# The code's single pass drops every burst after the second of a chain of three or more. For the
# carried-forward reading: bursts lost that way, and OUT/IN (the combination above) and best-R with
# every chain merged; OUT/IN also with the code's lowess smoothing (frac 0.0002), same pair.
def oi_combo(s, c, a, b):
drop = set(sc.hyperactive(s, c["strategy"]).tolist())
idx = [u for u in range(len(s["trains"])) if best_combo[0] == "all" or u not in drop]
t0_, t1_ = (0.0, s["mea_duration"]) if best_combo[1] == "whole" else (s["frame_times"][0], s["frame_times"][-1])
a, b = (a - 0.1, b + 0.5) if best_combo[2] == "extended" else (a, b)
return sc.out_in_ratio([s["trains"][u] for u in idx], a, b, t0_, t1_)
mrg = []
for f, rs in mapping.items():
s, c = sessions[f], chosen[rs[0]]
a, b, info = sc.detect_bursts(s, c["p_thresh"], c["strategy"], "span5")
af, bf, infof = sc.detect_bursts(s, c["p_thresh"], c["strategy"], "span5", merge="full")
al, bl, _ = sc.detect_bursts(s, c["p_thresh"], c["strategy"], "frac")
keep_f = sc.keep_frames(s, af, bf)
for r in rs:
rr = sc.pearson_with_units(chosen[r]["demixed"], s["U"], keep_f)
mrg.append({"row": r, "bursts": info["n_bursts"], "lost_by_single_pass": info["chain_lost"],
"bursts_full_merge": infof["n_bursts"], "out_in": round(oi_combo(s, c, a, b), 3),
"out_in_full_merge": round(oi_combo(s, c, af, bf), 3), "out_in_lowess": round(oi_combo(s, c, al, bl), 3),
"oi_table": table.out_in[r], "best_R": round(chosen[r]["best_R"], 3),
"best_R_full_merge": round(rr[sc.rank_units(rr)[0]], 3), "best_R_table": table.best_R[r]})
mrg = pd.DataFrame(mrg)
print(mrg.to_string(index=False))
one = mrg[mrg.row != 10]
for col in ["out_in", "out_in_full_merge", "out_in_lowess"]:
print(f"{col}: within 10% of Table 1 for {int(((one[col] - one.oi_table).abs() <= 0.1 * one.oi_table).sum())} of 19 recordings")
print("best-R with every chain merged, within 10% of Table 1:",
int(((mrg.best_R_full_merge - mrg.best_R_table).abs() <= 0.1 * mrg.best_R_table).sum()), "of 20 (carried forward:",
int(((mrg.best_R - mrg.best_R_table).abs() <= 0.1 * mrg.best_R_table).sum()), ")")
row bursts lost_by_single_pass bursts_full_merge out_in out_in_full_merge out_in_lowess oi_table best_R best_R_full_merge best_R_table 3 76 4 76 1.755 1.747 1.907 1.717 0.213 0.214 0.217 4 201 0 201 0.993 0.993 0.993 0.588 0.643 0.643 0.644 2 46 5 46 2.194 2.016 2.131 2.003 0.515 0.516 0.516 5 237 0 237 0.809 0.809 0.809 0.318 0.430 0.430 0.429 1 106 6 106 0.737 0.735 0.737 0.733 0.522 0.521 0.523 19 81 1 81 0.297 0.296 0.299 0.272 0.222 0.221 0.219 17 25 0 25 1.119 1.119 1.137 0.534 0.627 0.627 0.602 16 32 0 32 0.564 0.564 0.582 0.171 0.537 0.537 0.676 6 43 0 43 0.762 0.762 0.762 0.189 0.653 0.653 0.647 7 42 0 42 0.428 0.428 0.431 0.269 0.344 0.344 0.346 8 31 0 31 0.419 0.419 0.419 0.060 0.363 0.363 0.340 20 58 18 58 0.213 0.205 0.213 0.112 0.293 0.340 0.294 15 39 0 39 0.221 0.221 0.223 0.222 0.395 0.395 0.400 14 42 0 42 0.272 0.272 0.275 0.289 0.354 0.354 0.394 18 124 5 124 0.653 0.650 0.686 0.464 0.311 0.311 0.311 13 76 0 76 0.773 0.773 0.790 0.425 0.143 0.143 0.160 12 99 1 99 0.593 0.593 0.593 0.719 0.585 0.585 0.597 9 102 0 102 0.443 0.443 0.452 0.425 0.673 0.673 0.671 10 102 0 102 0.443 0.443 0.452 0.425 0.762 0.762 0.765 11 62 0 62 1.171 1.171 1.171 0.783 0.137 0.137 0.148 out_in: within 10% of Table 1 for 7 of 19 recordings out_in_full_merge: within 10% of Table 1 for 7 of 19 recordings out_in_lowess: within 10% of Table 1 for 6 of 19 recordings best-R with every chain merged, within 10% of Table 1: 16 of 20 (carried forward: 17 )
Result: not reproduced. Spike contrast matches Table 1 within 10% for all 19 recordings when all sorted units are used (16 of 19 without the hyperactive units). OUT/IN matches for at most 6 of 19 recordings under any combination (all units, extended bursts: rows 1, 2, 9, 14, 15 and 19), so both measures match for 6, below the 10 needed for a partial reproduction. OUT/IN here is mostly higher than Table 1's (row 8, 0.419 against 0.060; row 6, 0.762 against 0.189). With the bursts of any p_thresh and strategy, 10 of 19 recordings can be matched one by one, but no single reading matches them all, and row 8 is at least 0.139. The spike contrast part holds; the OUT/IN values of Table 1 could not be recovered with the burst detection described.
The code's single merging pass drops every burst after the second in a chain of three or more bursts closer than 75 ms. With the carried-forward reading this loses bursts in 7 of 19 recordings (1 to 6 bursts, 18 in row 20). Merging every chain into one burst, as Methods 2.6 says ("merged into a single larger burst"), changes OUT/IN by at most 0.18 (row 2, 2.194 to 2.016), and the code's lowess smoothing (frac 0.0002) with the same pair changes it by at most 0.25 (row 3, 2.293 to 2.547); each still matches 6 of 19 recordings. OUT/IN for the outcome uses the carried-forward bursts only; the match of 10 of 19 recordings, each under its own burst reading, is shown above but does not count under the criterion. For rows 14, 15, 16 and 20 the spike contrast also decided the row-to-file assignment, so its match there is not independent.
Claims fig5c-spike-contrast-vs-fscore and fig5d-out-in-vs-fscore¶
The network's spike contrast is negatively correlated with the maximum F-score of the recording (Spearman R = -0.737, P = 0.0002), so more synchronous networks give poorer classification.
The ratio of spikes outside to inside network bursts is positively correlated with the maximum F-score of the recording (Spearman R = 0.672, P = 0.0012).
In [27]:
rho = []
for g in range(3):
for col in ["maxF", "maxF_code"]:
for combo_name, sub in net.groupby(["units", "span", "bursts"]):
m20 = grid[grid.region == g].set_index("row")[col]
meas = sub.set_index("row")
rows20 = sorted(spine_rows)
x_sc = [meas.spike_contrast[r if r != 10 else 9] for r in rows20]
x_oi = [meas.out_in[r if r != 10 else 9] for r in rows20]
y = [m20[r] for r in rows20]
rows19 = [r for r in rows20 if r != 10]
i19 = [rows20.index(r) for r in rows19]
for nm, x in [("spike_contrast", x_sc), ("out_in", x_oi)]:
r20 = stats.spearmanr(x, y)
r19 = stats.spearmanr(np.array(x)[i19], np.array(y)[i19])
rho.append({"region": g, "F": col, "units": combo_name[0], "span": combo_name[1], "bursts": combo_name[2],
"measure": nm, "R_20": round(r20.statistic, 3), "P_20": r20.pvalue,
"R_19": round(r19.statistic, 3), "P_19": r19.pvalue})
rho = pd.DataFrame(rho)
main_rho = rho[(rho.units == best_combo[0]) & (rho.span == best_combo[1]) & (rho.bursts == best_combo[2])]
print("Combination of claim table1-network-measures:", best_combo)
print(main_rho.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
print("\nOver every combination, region and F reading: R range")
print(rho.groupby("measure")[["R_20", "R_19"]].agg(["min", "max"]).to_string())
print("Largest P over every combination, region and F reading:")
print(rho.groupby("measure")[["P_20", "P_19"]].max().to_string(float_format=lambda v: f"{v:.4g}"))
print("Paper: spike contrast R = -0.737 (P = 0.0002), OUT/IN R = 0.672 (P = 0.0012)")
Combination of claim table1-network-measures: ('all', 'imaging', 'extended')
region F units span bursts measure R_20 P_20 R_19 P_19
0 maxF all imaging extended spike_contrast -0.719 0.0003532 -0.714 0.000595
0 maxF all imaging extended out_in 0.564 0.009574 0.574 0.01022
0 maxF_code all imaging extended spike_contrast -0.741 0.0001848 -0.738 0.0003092
0 maxF_code all imaging extended out_in 0.6 0.005191 0.609 0.005645
1 maxF all imaging extended spike_contrast -0.719 0.0003532 -0.721 0.0004947
1 maxF all imaging extended out_in 0.594 0.005731 0.598 0.006815
1 maxF_code all imaging extended spike_contrast -0.735 0.0002219 -0.735 0.0003333
1 maxF_code all imaging extended out_in 0.614 0.003979 0.606 0.006006
2 maxF all imaging extended spike_contrast -0.723 0.0003187 -0.718 0.0005384
2 maxF all imaging extended out_in 0.561 0.01008 0.57 0.01076
2 maxF_code all imaging extended spike_contrast -0.739 0.0001952 -0.74 0.0002893
2 maxF_code all imaging extended out_in 0.619 0.003612 0.611 0.005497
Over every combination, region and F reading: R range
R_20 R_19
min max min max
measure
out_in 0.561 0.736 0.570 0.732
spike_contrast -0.764 -0.719 -0.767 -0.714
Largest P over every combination, region and F reading:
P_20 P_19
measure
out_in 0.01008 0.01076
spike_contrast 0.0003532 0.000595
Paper: spike contrast R = -0.737 (P = 0.0002), OUT/IN R = 0.672 (P = 0.0012)
In [28]:
fig, axs = plt.subplots(1, 2, figsize=(9, 4))
g0 = grid[grid.region == 0].set_index("row")
for ax, nm, lab in [(axs[0], "spike_contrast", "spike contrast"), (axs[1], "out_in", "OUT/IN ratio")]:
xx = [bc.set_index("row")[nm][r if r != 10 else 9] for r in sorted(spine_rows)]
ax.plot(xx, [g0.maxF[r] for r in sorted(spine_rows)], "o", label="here (region 0)")
ax.plot([table[nm][r] for r in sorted(spine_rows)], [table.max_F[r] for r in sorted(spine_rows)], "x", color="0.5", label="Table 1")
ax.set_xlabel(lab)
ax.set_ylabel("max F")
ax.legend(fontsize=8)
print("Plotted: all points inside the automatic axes")
Plotted: all points inside the automatic axes
Result: fig5c-spike-contrast-vs-fscore reproduced; fig5d-out-in-vs-fscore partially reproduced (resting on an assumption). Over the 20 rows, spike contrast against maximum F gives Spearman R = -0.731 (P = 0.0003) with noise region 0 and the Methods' classification, and -0.726 (P = 0.0004) over 19 recordings. Every combination, region and classification gives -0.726 to -0.784 with P at most 0.0004, all within -0.737 ± 0.074.
OUT/IN against maximum F gives R = 0.585 (P = 0.007) over the 20 rows and 0.598 (P = 0.007) over 19 recordings with region 0 and the Methods' classification, both below 0.672 ± 0.067 (bound 0.605). Region 2 gives 0.582 and 0.595, also below. Region 1 gives 0.615 and 0.623, and the code's classification 0.625 to 0.645 over 20 rows, within tolerance. So under the claims file's main classification the outcome depends on the noise region: region 0, the background square nearest the spine, is taken as the Methods' "background ROI close to the spine". The 19-recording reading leaves out row 10, whose squares 0 and 2 sit on processes, and still gives 0.598 with region 0. Over every combination R is 0.582 to 0.754, always positive with P at most 0.0072. The OUT/IN values themselves differ from Table 1 (claim table1-network-measures).
Summary¶
| claim | outcome |
|---|---|
| data-scale | reproduced |
| table1-best-r | partially reproduced |
| fig1e-best-r-gap | reproduced |
| table1-optimal-thresholds | partially reproduced |
| mean-max-fscore | partially reproduced |
| acceptance-9-of-20 | partially reproduced |
| best-r-accepted-vs-rejected | partially reproduced |
| best-r-not-predictive | partially reproduced |
| fig5a-rejection-scenarios | reproduced |
| fig5a-threshold-order | reproduced |
| table1-network-measures | not reproduced |
| fig5c-spike-contrast-vs-fscore | reproduced |
| fig5d-out-in-vs-fscore | partially reproduced (assumption: noise region) |
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 convolution is computed at the frame times, not at every 20 kHz sample (the largest difference from elephant is 0.3 Hz on peaks of 19 to 45 Hz). The paper chose p_thresh by eye for each recording and does not report it; the claims file's rule picks the reading nearest to Table 1's best-R. The background region for the noise estimate is not in the dataset; three squares near each spine are placed by a fixed rule, and region 0, the nearest, is used where a choice matters (only the OUT/IN correlation of Figure 5(d) depends on it). A burst ends where the smoothed count falls to the baseline; the smoothed counts are rounded to 9 decimals first, because many equal the integer baseline up to rounding error and the comparison otherwise differed between machines. Units with an undefined R are left out of the ranking. The release probabilities are a fresh draw. Table 1 rows are matched to files by culture and unit count, with ties broken by spike contrast; swapping the tied rows leaves the outcome of table1-best-r unchanged. Spike contrast is computed with elephant 1.2.1; the paper used 0.10.0. The burst merging is the code's single pass, which drops bursts in chains of three or more; merging every chain, as the Methods say, leaves the OUT/IN matches at 6 of 19. Every later claim uses the burst reading whose best-R is nearest to Table 1, so its outcome is conditional on that choice. Between 6% and 44% of the simulated ground-truth surrogates have a best unit other than their own, a share the Methods call rare; these decide how far the two classifications of the F-score differ.
What was found about the dataset. The 20 files hold 19 recordings: sub-3112_ecephys+ophys.nwb
is the same blob as sub-3112_ses-20210227T190233. Nothing maps files to the paper's Table 1 rows
or cultures, and two pairs of recordings cannot be told apart by unit count. The background ROIs
the paper used for the noise estimate are not in the files. The ROIs are not labelled as spine or
shaft. The imaging plane's indicator is "unknown", where the paper used jGCaMP7b. The stored Dff
series is not the dF/F the paper describes (differences up to 0.23). Two files each hold one
unit with no spikes, which Table 1 does not count. The two streams share a clock: the shaft trace
correlates best with the population rate at zero lag in 19 of 20 spines and at one frame in the
other.
In [29]:
print(f"Notebook run time: {(time.time() - started) / 60:.1f} min")
Notebook run time: 22.1 min