Reproducing Pakhomov et al. (2023) from dandiset 001027¶
Paper. Pakhomov AG, Gudvangen E, Mangalanathan U, Kondratiev O, Redondo L, Semenov I (2023). Next generation CANCAN focusing for remote stimulation by nanosecond electric pulses. Bioelectrochemistry 152:108437. doi:10.1016/j.bioelechem.2023.108437. The claims file was written from the accepted author manuscript in PubMed Central (PMC10247516), with its Methods and figure legends; the version of record is not open access, and there is no supplement or linked code.
Dandiset. DANDI:001027, version 0.240605.1859. The paper cites no version of the dandiset. This is the only published version, and the draft has the same ten assets and size (3,667,103,151 bytes) and was last modified when this version was published, so the most recent published version is pinned.
What the dandiset holds. Ten NWB files, one dish each, with one stitched 10x YO-PRO-1 image of a CHO monolayer after a pulse protocol in a quadrupole array: Packet 1 (3 dishes), Packet 2 (4) and Packet 3 (3), each rotated once from every electrode, first pulse 4 kV. This is the experiment of Figure 8. The single-cell time-lapse of Figure 4 and the monolayer experiments of Figures 5, 6, 7 and 9 are not in it.
| Claim | Tested here |
|---|---|
| fig4-unipolar-strongest, fig4-bipolar-weaker, fig4-sham-no-uptake, fig4-twenty-fold-cancellation | Blocked, data missing (single-cell time-lapse not deposited) |
| fig5a-trains-stronger-throughout, fig5a-packet-center-peak-e1-stronger, fig5bc-rotation-focuses-center | Blocked, data missing |
| fig6a-more-packets-more-uptake, fig6b-center-to-periphery-by-number, fig6-peak-narrower | Blocked, data missing |
| fig7a-amplitude-scales-uniformly, fig7bd-targeting-independent-of-amplitude, array-size-independence | Blocked, data missing |
| fig8a-equal-center-uptake | Tested |
| fig8a-more-uptake-near-electrodes | Tested |
| fig8b-packet2-best-targeting | Tested |
| fig9-gap-fades-periphery-faster, fig9b-gap-boosts-targeting | Blocked, data missing |
The helper module cancan.py streams the files, finds the electrodes and measures the ROI
profiles; its docstring says where each parameter comes from.
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import cancan
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
dishes = cancan.load_all()
for d in dishes:
d["electrodes"] = cancan.find_electrodes(d["img"], d["valid"])
d["background"], d["background_each"] = cancan.background(d["img"], d["electrodes"])
print(f"loaded {len(dishes)} files from dandiset {cancan.DANDISET_ID} version {cancan.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
loaded 10 files from dandiset 001027 version 0.240605.1859 in 246 s
Data check¶
What was loaded, and each check of the metadata against the values. The files hold one image each and no time series, so there is no trial table or clock to align; the check that the data streams share a clock does not apply.
info = pd.DataFrame([{
"file": d["path"], "subject_id": d["subject_id"], "packet": d["packet"], "first_pulse_kV": d["kv"],
"session_start": d["session_start"][:16], "shape": d["shape"], "dtype": d["dtype"],
"max_value": d["max_value"], "unimaged_frac": round(d["zero_frac"], 3), "sha256": d["sha256"][:12]}
for d in dishes]).sort_values(["packet", "session_start"])
print(info.to_string(index=False))
print("\nDishes per packet (paper, Fig. 8 legend: n = 3-4):")
print(info.packet.value_counts().sort_index().to_string())
print("\nFirst-pulse amplitudes in the descriptions (kV):", sorted(info.first_pulse_kV.unique()))
print("Distinct SHA-256 digests:", info.sha256.nunique(), "of", len(info))
print("\nPacket wording in the session descriptions:")
for p in (1, 2, 3):
desc = next(d["description"] for d in dishes if d["packet"] == p)
print(f" Protocol{p}:", desc[desc.find(f"Protocol{p}"):desc.find("Detailed")].strip())
print("\nOne of the session descriptions in full:\n ", dishes[0]["description"])
file subject_id packet first_pulse_kV session_start shape dtype max_value unimaged_frac sha256 sub-19/sub-19_image.nwb 19 1 4.0 2021-03-24T10:54 (1, 20480, 20480) uint16 65444 0.259 8c5dbdb1eba7 sub-20/sub-20_image.nwb 20 1 4.0 2021-03-24T10:59 (1, 20480, 20480) uint16 65441 0.259 bb3962b16f3d sub-21/sub-21_image.nwb 21 1 4.0 2021-03-24T11:06 (1, 20480, 20480) uint16 65444 0.259 3e3e718eeb1c sub-2/sub-2_image.nwb 2 2 4.0 2021-03-17T11:36 (1, 20480, 20480) uint16 65107 0.259 00c4676edd72 sub-12/sub-12_image.nwb 12 2 4.0 2021-03-17T12:40 (1, 20480, 20480) uint16 65437 0.259 7ea479a33fdb sub-13/sub-13_image.nwb 13 2 4.0 2021-03-17T12:59 (1, 20480, 20480) uint16 65440 0.259 19bb7eb991aa sub-14/sub-14_image.nwb 14 2 4.0 2021-03-17T13:03 (1, 20480, 20480) uint16 65440 0.259 aa1d196a63f9 sub-03/sub-03_image.nwb 03 3 4.0 2021-03-17T11:40 (1, 20480, 20480) uint16 52086 0.259 fe0edd003fbd sub-11/sub-11_image.nwb 11 3 4.0 2021-03-17T12:35 (1, 20480, 20480) uint16 57698 0.259 a45fb816af0f sub-10/sub-10_image.nwb 10 3 4.0 2021-03-18T18:50 (1, 20480, 20480) uint16 65440 0.259 5cec51c2fe4b Dishes per packet (paper, Fig. 8 legend: n = 3-4): packet 1 3 2 4 3 3 First-pulse amplitudes in the descriptions (kV): [np.float64(4.0)] Distinct SHA-256 digests: 10 of 10 Packet wording in the session descriptions: Protocol1: Protocol1 had 4 phases of pulses, 600ns pulse duration, canceling pulses delivered without delay. Protocol2: Protocol2 had 8 phases of pulses, 600ns pulse duration, canceling pulses delivered without delay. Protocol3: Protocol3 had 8 phases of pulses, 600ns pulse duration, canceling pulses delivered with 5us delays before another set of pulses. One of the session descriptions in full: 2. Experiments were conducted in a 4 electrodes setup. CanCan protocol applied at 4.0kV (first pulse) with the subsequent canceling pulses 12.5% lower than pulses before. Packets were repeated 1 times at 1Hz frequency. Than also at 1Hz frequency protocol was shifted and than electrode 2 was starting electrode with noted amplitude of the first pulse. Protocol finished when all 4 electrodes were active pulse electrodes. Protocol2 had 8 phases of pulses, 600ns pulse duration, canceling pulses delivered without delay. Detailed information on pulsing parameters, electrodes and procedures can be found in 10.1016/j.bioelechem.2023.108437.
Every description names a first pulse of 4.0 kV and one rotation (one packet from each electrode), as the Figure 8 legend says, and the packet counts (3, 4, 3) fit its n = 3-4. The descriptions define Protocol 1 as 4 phases, Protocol 2 as 8 phases without delay and Protocol 3 as 8 phases with 5-us delays, as the paper defines Packets 1, 2 and 3 (Section 2.3).
meta = pd.DataFrame([{
"subject_id": d["subject_id"], "unit": d["unit"], "resolution": d["resolution"], "conversion": d["conversion"],
"species": d["species"].split("/")[-1], "strain": d["strain"], "age": d["age"],
"device": d["device"].split("Additional metadata: ")[1]} for d in dishes])
print(meta.to_string(index=False))
# Are any two images the same recording? Compare the reduced images on a common 64-fold grid.
small = {d["subject_id"]: d["img"].reshape(160, 16, 160, 16).mean(axis=(1, 3)).ravel() for d in dishes}
ids = list(small)
corr = pd.DataFrame(np.corrcoef([small[i] for i in ids]), index=ids, columns=ids)
off = corr.where(~np.eye(len(ids), dtype=bool))
print("\nLargest correlation between two different dishes' images:", round(float(np.nanmax(off.values)), 3))
print("Any two reduced images identical:", any(np.array_equal(dishes[i]["img"], dishes[j]["img"])
for i in range(len(dishes)) for j in range(i + 1, len(dishes))))
subject_id unit resolution conversion species strain age device
2 n.a. 96.0 1.0 NCBITaxon_10029 cellline: CHO P1D Microscope: Unknown microscope, Detector: Hamamatsu Hamamatsu ORCA-Flash4.0, Objective: Olympus IX3 Nosepiece NA=0.40000000000000002 Mag=10, Channels: FITC
03 n.a. 96.0 1.0 NCBITaxon_10029 cellline: CHO P1D Microscope: Unknown microscope, Detector: Hamamatsu Hamamatsu ORCA-Flash4.0, Objective: Olympus IX3 Nosepiece NA=0.40000000000000002 Mag=10, Channels: FITC
10 n.a. 1.0 1.0 NCBITaxon_10029 cellline: CHO P1D Microscope: Unknown microscope, Detector: Unknown detector, Objective: Unknown objective, Channels: Unknown channels
11 n.a. 96.0 1.0 NCBITaxon_10029 cellline: CHO P1D Microscope: Unknown microscope, Detector: Hamamatsu Hamamatsu ORCA-Flash4.0, Objective: Olympus IX3 Nosepiece NA=0.40000000000000002 Mag=10, Channels: FITC
12 n.a. 96.0 1.0 NCBITaxon_10029 cellline: CHO P1D Microscope: Unknown microscope, Detector: Hamamatsu Hamamatsu ORCA-Flash4.0, Objective: Olympus IX3 Nosepiece NA=0.40000000000000002 Mag=10, Channels: FITC
13 n.a. 96.0 1.0 NCBITaxon_10029 cellline: CHO P1D Microscope: Unknown microscope, Detector: Hamamatsu Hamamatsu ORCA-Flash4.0, Objective: Olympus IX3 Nosepiece NA=0.40000000000000002 Mag=10, Channels: FITC
14 n.a. 96.0 1.0 NCBITaxon_10029 cellline: CHO P1D Microscope: Unknown microscope, Detector: Hamamatsu Hamamatsu ORCA-Flash4.0, Objective: Olympus IX3 Nosepiece NA=0.40000000000000002 Mag=10, Channels: FITC
19 n.a. 96.0 1.0 NCBITaxon_10029 cellline: CHO P1D Microscope: Unknown microscope, Detector: Hamamatsu Hamamatsu ORCA-Flash4.0, Objective: Olympus IX3 Nosepiece NA=0.40000000000000002 Mag=10, Channels: FITC
20 n.a. 96.0 1.0 NCBITaxon_10029 cellline: CHO P1D Microscope: Unknown microscope, Detector: Hamamatsu Hamamatsu ORCA-Flash4.0, Objective: Olympus IX3 Nosepiece NA=0.40000000000000002 Mag=10, Channels: FITC
21 n.a. 96.0 1.0 NCBITaxon_10029 cellline: CHO P1D Microscope: Unknown microscope, Detector: Hamamatsu Hamamatsu ORCA-Flash4.0, Objective: Olympus IX3 Nosepiece NA=0.40000000000000002 Mag=10, Channels: FITC
Largest correlation between two different dishes' images: 0.864
Any two reduced images identical: False
The ten files are distinct recordings: no two digests or images are identical (the largest
correlation between two dishes' reduced images is 0.86). The image unit is "n.a." and the
resolution field is 96 in nine files and 1.0 in dish 10, neither a pixel size, so the files do
not give the scale of the image. The subject is a dish of the CHO cell line, and the species identifier
NCBITaxon 10029 is Cricetulus griseus, the Chinese hamster; the dandiset's summary names it as
"Cricetulus griseus - Cricetulus aureus". The age "P1D" has no meaning for a cell line. Nine files
name the ORCA-Flash4.0 camera and a 10x objective (Methods: "10x, 0.38 NA"; the files say NA 0.40);
dish 10 names no device, and its session starts on 18 March at 18:50, the evening after the other
Packet 2 and 3 dishes.
print("Rows (pixels) at which the left edge of the imaged area steps (the stitching tile grid):")
for d in dishes[:3]:
st = cancan.tile_step(d["valid"])
print(f" dish {d['subject_id']}: {st}; differences {np.diff(st).tolist()}")
print("ORCA-Flash4 sensor: 2048 x 2048 pixels of 6.5 um; at 10x, 0.65 um per pixel.")
print("Assumed scale:", cancan.UM_PER_PX, "um per pixel; one block of the reduced image =",
cancan.BLOCK * cancan.UM_PER_PX, "um")
geo = []
for d in dishes:
g = cancan.geometry(d["electrodes"])
geo.append({"dish": d["subject_id"], "packet": d["packet"],
"diag_cc_mm": np.round(g["diag_cc_mm"], 2).tolist(),
"diag_edge_mm": np.round(g["diag_edge_mm"], 2).tolist(),
"side_cc_mm": np.round(g["side_cc_mm"], 2).tolist(),
"footprint_diameter_mm": np.round(g["disk_diameter_mm"], 2).tolist(),
"half_edge_minus_5pct_mm": round(float(np.mean(g["diag_edge_mm"]) * (0.5 - cancan.EXCLUDE_FRAC)), 2)})
geo = pd.DataFrame(geo).sort_values(["packet", "dish"])
print("\nElectrode geometry found in each image (mm at the assumed scale):")
print(geo.to_string(index=False))
allg = [cancan.geometry(d["electrodes"]) for d in dishes]
cc = np.concatenate([g["diag_cc_mm"] for g in allg]); ee = np.concatenate([g["diag_edge_mm"] for g in allg])
fp = np.concatenate([g["disk_diameter_mm"] for g in allg])
print(f"\nCenter-to-center diagonal {cc.min():.2f} to {cc.max():.2f} mm; edge-to-edge {ee.min():.2f} to {ee.max():.2f} mm; "
f"footprints {fp.min():.2f} to {fp.max():.2f} mm across (needles 1.5 mm, Methods 2.3).")
print("Paper: ROIs cover +/-1.65 mm (small array) or +/-3.65 mm (large array, 8.4-mm diagonal), "
"excluding <5% of the diagonal next to the electrodes.")
Rows (pixels) at which the left edge of the imaged area steps (the stitching tile grid): dish 2: [1848, 3688, 7376, 13112, 16800, 18640]; differences [1840, 3688, 5736, 3688, 1840] dish 03: [1848, 3688, 7376, 13112, 16800, 18640]; differences [1840, 3688, 5736, 3688, 1840] dish 10: [1848, 3688, 7376, 13112, 14744, 16792, 18640]; differences [1840, 3688, 5736, 1632, 2048, 1848] ORCA-Flash4 sensor: 2048 x 2048 pixels of 6.5 um; at 10x, 0.65 um per pixel. Assumed scale: 0.65 um per pixel; one block of the reduced image = 5.2 um Electrode geometry found in each image (mm at the assumed scale): dish packet diag_cc_mm diag_edge_mm side_cc_mm footprint_diameter_mm half_edge_minus_5pct_mm 19 1 [5.37, 5.31] [3.57, 3.47] [3.52, 3.93, 3.91, 3.75] [1.98, 1.74, 1.61, 1.93] 1.59 20 1 [5.34, 5.34] [3.56, 3.41] [3.53, 3.93, 3.93, 3.73] [1.97, 1.9, 1.59, 1.98] 1.57 21 1 [5.41, 5.44] [3.62, 3.66] [3.68, 3.87, 4.02, 3.77] [1.98, 1.68, 1.59, 1.88] 1.64 12 2 [5.34, 5.33] [3.58, 3.5] [3.67, 3.86, 3.83, 3.73] [1.9, 1.7, 1.64, 1.97] 1.59 13 2 [5.27, 5.39] [3.45, 3.51] [3.52, 3.87, 4.0, 3.71] [2.04, 1.75, 1.61, 2.01] 1.56 14 2 [5.39, 5.37] [3.62, 3.56] [3.66, 3.85, 3.92, 3.79] [1.94, 1.69, 1.6, 1.92] 1.62 2 2 [5.32, 5.42] [3.5, 3.48] [3.64, 3.92, 3.78, 3.85] [1.93, 1.9, 1.71, 1.99] 1.57 03 3 [5.34, 5.53] [3.49, 3.65] [3.68, 3.87, 3.95, 3.88] [2.02, 1.83, 1.68, 1.93] 1.61 10 3 [5.29, 5.45] [3.47, 3.51] [3.7, 3.74, 3.86, 3.88] [1.97, 1.91, 1.67, 1.96] 1.57 11 3 [5.32, 5.27] [3.49, 3.49] [3.64, 3.86, 3.85, 3.62] [2.03, 1.64, 1.62, 1.9] 1.57 Center-to-center diagonal 5.27 to 5.53 mm; edge-to-edge 3.41 to 3.66 mm; footprints 1.59 to 2.04 mm across (needles 1.5 mm, Methods 2.3). Paper: ROIs cover +/-1.65 mm (small array) or +/-3.65 mm (large array, 8.4-mm diagonal), excluding <5% of the diagonal next to the electrodes.
The image is stitched from tiles. In the three dishes shown, the left edge of the imaged area
steps at rows 1840 to 1848 pixels apart or twice that (3688), as tiles about 1844 pixels tall would
make it; one step in each of them is 5736 pixels, which is not a multiple of that spacing, and
dish 10 has steps of 1632 and 2048 in place of the second step of 3688 in the other two. A spacing of about 1844
pixels is about 90% of the camera's 2048, which is consistent with unbinned tiles overlapping by
10% and not resampled, and so with 0.65 um per pixel. At that scale the electrode footprints are
1.6 to 2.0 mm across, a little wider than the 1.5-mm needles, as a crushed or unexposed rim would
make them. The arrays are the small one: the edge-to-edge diagonal is about 3.5 mm, and 5% short of
each electrode leaves about 1.6 mm from the center, which is the +/-1.65 mm span the Methods give
for the small array. The large array (8.4-mm diagonal) would need a scale about 2.4 times larger,
and its footprints would then be about 4.5 mm across, three times the needle. The scale is
consistent with the data but is not stated in them, so the ROIs are also placed by a rule that
needs no scale (rel, below).
bgt = pd.DataFrame([{"dish": d["subject_id"], "packet": d["packet"], "date": d["session_start"][:10],
"background_mean": round(d["background"], 1),
"background_per_electrode": np.round(d["background_each"], 1).tolist()} for d in dishes])
print("Median intensity inside each electrode footprint (central 60% of the radius), raw counts:")
print(bgt.sort_values(["packet", "dish"]).to_string(index=False))
print("\nBy imaging date:")
print(bgt.groupby("date").background_mean.agg(["count", "mean", "min", "max"]).round(1).to_string())
# Saturation within the ROI strip of each dish (the main ROI geometry), from the per-block fraction.
rows = []
for d in dishes:
e = d["electrodes"]; c = np.array([(a, b) for a, b, _ in e]); cen = c.mean(0)
yy, xx = np.mgrid[:d["img"].shape[0], :d["img"].shape[1]]
half = cancan.HALF_SPAN_MM * 1000 / cancan.UM_PER_PX / cancan.BLOCK
mask = np.zeros(d["img"].shape, bool)
for i, j in [(0, 2), (1, 3)]:
u = (c[j] - c[i]) / np.hypot(*(c[j] - c[i])); nrm = np.array([-u[1], u[0]])
along = (yy - cen[0]) * u[0] + (xx - cen[1]) * u[1]
across = (yy - cen[0]) * nrm[0] + (xx - cen[1]) * nrm[1]
mask |= (np.abs(along) <= half) & (np.abs(across) <= half / 46)
rows.append({"dish": d["subject_id"], "saturated_frac_in_rois": float(d["saturated"][mask].mean()),
"saturated_frac_image": float(d["saturated"].mean())})
sat = pd.DataFrame(rows)
print("\nFraction of pixels >= 65000 (of 65535):")
print(sat.to_string(index=False, float_format=lambda v: f"{v:.2e}"))
Median intensity inside each electrode footprint (central 60% of the radius), raw counts:
dish packet date background_mean background_per_electrode
19 1 2021-03-24 304.8 [350.8, 287.9, 278.4, 301.9]
20 1 2021-03-24 284.9 [315.0, 247.5, 281.3, 295.6]
21 1 2021-03-24 334.5 [362.3, 341.6, 311.3, 322.7]
12 2 2021-03-17 150.4 [143.4, 148.2, 163.1, 146.7]
13 2 2021-03-17 153.6 [149.1, 140.9, 183.1, 141.4]
14 2 2021-03-17 168.9 [195.3, 146.0, 159.3, 174.9]
2 2 2021-03-17 187.4 [181.1, 160.3, 186.6, 221.5]
03 3 2021-03-17 171.9 [166.2, 169.8, 160.5, 191.0]
10 3 2021-03-18 150.3 [150.2, 156.6, 162.2, 132.2]
11 3 2021-03-17 138.3 [143.8, 139.9, 119.1, 150.4]
By imaging date:
count mean min max
date
2021-03-17 6 161.8 138.3 187.4
2021-03-18 1 150.3 150.3 150.3
2021-03-24 3 308.1 284.9 334.5
Fraction of pixels >= 65000 (of 65535): dish saturated_frac_in_rois saturated_frac_image 2 0.00e+00 2.38e-09 03 0.00e+00 0.00e+00 10 0.00e+00 1.80e-06 11 0.00e+00 0.00e+00 12 0.00e+00 1.17e-07 13 0.00e+00 1.40e-06 14 0.00e+00 4.15e-07 19 0.00e+00 6.57e-06 20 0.00e+00 5.46e-07 21 0.00e+00 5.77e-06
The footprints' interior is the only region of each image where cells were not exposed, and the files hold no sham dish or cell-free region, so it serves as the background in the background-subtracted reading. The intensity there is about twice as high in the three Packet 1 dishes (imaged on 24 March 2021) as in the seven others (17 and 18 March), so the a.u. scale differs between imaging days. The files do not record exposure or gain. Saturated pixels are too few to affect an ROI mean.
fig, axes = plt.subplots(2, 5, figsize=(20, 8.5))
order = sorted(dishes, key=lambda d: (d["packet"], int(d["subject_id"])))
for ax, d in zip(axes.flat, order):
img = d["img"]
vmax = np.percentile(img[d["valid"]], 99.5)
ax.imshow(img, cmap="gray", vmin=0, vmax=vmax)
e = d["electrodes"]; c = np.array([(a, b) for a, b, _ in e]); cen = c.mean(0)
for k, (cy, cx, r) in enumerate(e):
ax.add_patch(plt.Circle((cx, cy), r, fill=False, color="tab:red", lw=1))
ax.text(cx, cy, f"E{k}", color="yellow", ha="center", va="center", fontsize=8)
half = cancan.HALF_SPAN_MM * 1000 / cancan.UM_PER_PX / cancan.BLOCK
for i, j, col in [(0, 2, "cyan"), (1, 3, "orange")]:
u = (c[j] - c[i]) / np.hypot(*(c[j] - c[i]))
a, b = cen - u * half, cen + u * half
ax.plot([a[1], b[1]], [a[0], b[0]], color=col, lw=1)
ax.set_xlim(cen[1] - 900, cen[1] + 900); ax.set_ylim(cen[0] + 900, cen[0] - 900)
ax.set_title(f"dish {d['subject_id']}, Protocol {d['packet']} (0 to {vmax:.0f})", fontsize=9)
ax.set_xticks([]); ax.set_yticks([])
fig.suptitle("Electrode footprints found (red) and the ROI span of each diagonal (cyan d1: E0-E2, orange d2: E1-E3); "
"grey scale 0 to the 99.5th percentile, shown in each title")
plt.tight_layout(); plt.show()
print("Each panel spans 1800 x 1800 blocks around the array center; every electrode lies inside it.")
Each panel spans 1800 x 1800 blocks around the array center; every electrode lies inside it.
Profiles under every reading¶
The paper gives the span and number of ROIs but leaves open the details below, which the claims file names or which the Methods leave open (author, step 1.6). Each is run both or all ways:
geom: ROIs cover +/-1.65 mm from the geometrical center at 0.65 um per pixel (abs; Methods 2.4), or run to 5% of the edge-to-edge diagonal short of each electrode footprint (rel; Methods 2.4, "<5% of the diagonal length"), which needs no scale.n_per_diag: "92 ... ROI, placed along one or both diagonals" (Methods 2.4) read as 46 per diagonal (92 over both) or 92 per diagonal.center: the mean of the two ROIs flanking 0, or one ROI centered at 0 (claims file).diag: the mean of both diagonals, or one (d1, E0 to E2; d2, E1 to E3). Which electrode is e1 is not in the files, so d1 and d2 are both run for "one diagonal".bg: intensities as recorded, or less the per-dish background inside the electrode footprints. The Methods say nothing of a background for monolayers (Section 2.2 subtracts one for single cells).
The main reading, used for the measured values, is abs, 46 per diagonal, flank, both,
raw: the stated span and ROI count, both diagonals as Figure 5B averages them, and the a.u. values
as recorded, as the claims file says they are judged.
rows, profs = [], []
for d in sorted(dishes, key=lambda d: (d["packet"], int(d["subject_id"]))):
r, p = cancan.dish_readings(d)
rows += r; profs.append((d["subject_id"], d["packet"], d["background"], p))
R = pd.DataFrame(rows)
R["pct"] = 100 * R.center_au / R.periph_au
keys = ["geom", "n_per_diag", "center", "bg", "diag"]
print(len(R), "rows:", R.dish.nunique(), "dishes x", len(R) // R.dish.nunique(), "readings")
main = (R.geom == "abs") & (R.n_per_diag == 46) & (R.center == "flank") & (R.bg == "raw") & (R.diag == "both")
print("\nMain reading, per dish:")
print(R[main][["dish", "packet", "center_au", "periph_au", "periph_min_end_au", "periph_max_end_au", "pct"]]
.round(1).to_string(index=False))
print("\nMain reading, background subtracted, per dish:")
mainbg = (R.geom == "abs") & (R.n_per_diag == 46) & (R.center == "flank") & (R.bg == "bgsub") & (R.diag == "both")
print(R[mainbg][["dish", "packet", "center_au", "periph_au", "pct"]].round(1).to_string(index=False))
480 rows: 10 dishes x 48 readings Main reading, per dish: dish packet center_au periph_au periph_min_end_au periph_max_end_au pct 19 1 598.700012 891.7 632.7 1285.3 67.1 20 1 674.799988 729.8 629.1 898.2 92.5 21 1 672.500000 799.1 475.3 1014.3 84.2 2 2 385.000000 363.4 213.9 518.1 105.9 12 2 332.899994 525.7 448.0 639.0 63.3 13 2 415.000000 535.0 357.3 691.9 77.6 14 2 421.200012 499.3 435.5 631.5 84.4 03 3 413.200012 301.7 222.3 393.3 137.0 10 3 382.700012 306.2 255.6 345.8 125.0 11 3 354.399994 280.7 238.5 335.3 126.3 Main reading, background subtracted, per dish: dish packet center_au periph_au pct 19 1 294.000000 587.0 50.1 20 1 390.000000 444.9 87.6 21 1 338.000000 464.6 72.7 2 2 197.600006 176.0 112.3 12 2 182.500000 375.4 48.6 13 2 261.399994 381.4 68.5 14 2 252.300003 330.4 76.4 03 3 241.300003 129.8 185.9 10 3 232.399994 155.9 149.1 11 3 216.100006 142.4 151.8
fig, axes = plt.subplots(1, 2, figsize=(15, 5))
colors = {1: "tab:blue", 2: "tab:orange", 3: "tab:green"}
ranges = []
for ax, sub in zip(axes, [False, True]):
for p in (1, 2, 3):
ps = [x for x in profs if x[1] == p]
t, prof_ = cancan.mean_profile([x[3] for x in ps], ("abs", 46, "flank"),
bg=[x[2] for x in ps] if sub else None)
m, se = prof_.mean(0), prof_.std(0, ddof=1) / np.sqrt(len(ps))
ax.errorbar(t, m, se, color=colors[p], label=f"Protocol {p} (n={len(ps)})", capsize=2, lw=1)
ranges.append((sub, p, float((m - se).min()), float((m + se).max())))
ax.set_xlabel("distance from the array center along the diagonal (mm)")
ax.set_ylabel("YP emission (a.u.)" + (", less background" if sub else ", raw"))
ax.legend(); ax.set_title("mean +/- s.e. over dishes, mean of both diagonals")
plt.tight_layout(); plt.show()
for sub, p, lo, hi in ranges:
print(f"{'background subtracted' if sub else 'raw':>22}, Protocol {p}: plotted range {lo:.0f} to {hi:.0f} a.u.; no points outside the axes")
raw, Protocol 1: plotted range 430 to 941 a.u.; no points outside the axes
raw, Protocol 2: plotted range 272 to 534 a.u.; no points outside the axes
raw, Protocol 3: plotted range 215 to 500 a.u.; no points outside the axes
background subtracted, Protocol 1: plotted range 135 to 635 a.u.; no points outside the axes
background subtracted, Protocol 2: plotted range 104 to 377 a.u.; no points outside the axes
background subtracted, Protocol 3: plotted range 71 to 356 a.u.; no points outside the axes
# Each dish's profile, main geometry (mean of both diagonals), raw and less its own background.
fig, axes = plt.subplots(2, 3, figsize=(18, 9), sharex=True)
yr = []
for row, sub in enumerate([False, True]):
for col, p in enumerate((1, 2, 3)):
ax = axes[row, col]
for sid, q, bgv, pr in profs:
if q != p:
continue
t, v = cancan.mean_profile([pr], ("abs", 46, "flank"), bg=[bgv] if sub else None)
ax.plot(t, v[0], lw=1, label=f"dish {sid}")
yr.append((sub, p, sid, float(v[0].min()), float(v[0].max())))
ax.set_title(f"Protocol {p}, " + ("less background" if sub else "raw"), fontsize=10)
ax.set_ylabel("YP emission (a.u.)"); ax.legend(fontsize=8)
if row == 1:
ax.set_xlabel("distance from the array center along the diagonal (mm)")
plt.tight_layout(); plt.show()
yr = pd.DataFrame(yr, columns=["bg_subtracted", "packet", "dish", "min_au", "max_au"]).round(1)
print("Plotted range of each dish's profile (a.u.); no points outside the axes:")
print(yr.to_string(index=False))
Plotted range of each dish's profile (a.u.); no points outside the axes:
bg_subtracted packet dish min_au max_au
False 1 19 452.2 1033.9
False 1 20 381.6 895.5
False 1 21 522.2 884.6
False 2 2 257.7 437.0
False 2 12 250.7 543.5
False 2 13 286.4 545.4
False 2 14 281.4 586.6
False 3 03 256.2 561.3
False 3 10 223.6 473.0
False 3 11 206.4 548.0
True 1 19 147.4 729.1
True 1 20 96.8 610.7
True 1 21 187.7 550.1
True 2 2 70.3 249.6
True 2 12 100.4 393.1
True 2 13 132.8 391.7
True 2 14 112.5 417.7
True 3 03 84.3 389.4
True 3 10 73.3 322.7
True 3 11 68.1 409.7
Most dishes show the shape of their group's mean profile. The Protocol 3 dishes peak broadly at the center (about 400 to 470 a.u. raw) and are lowest toward the electrodes, with narrow spikes above this in dish 03 near one end (its maximum, 561.3 a.u.) and in dish 11 at about -0.2 mm (its maximum, 548.0 a.u.) and at about +0.6 mm. Protocol 2 dishes 12, 13 and 14 are highest at the electrode ends (maxima 543.5 to 586.6 a.u. raw), with a lower central peak. Protocol 2 dish 2 is the exception: it is highest at one end (437.0 a.u.), but its other end (about 290 to 325 a.u.) is below its center (about 400), and its normalized center is 105.9% raw and 112.3% with the background subtracted (cell above), against 63.3 to 84.4% raw in the other three Protocol 2 dishes. Even so, every Protocol 3 dish (125.0 to 137.0% raw) is above every Protocol 2 dish, so the group difference does not rest on one dish. The Protocol 1 dishes are high at both ends and in the center, and dish 20 is lower than the other two on one side. Subtracting each dish's background shifts the curves and keeps these shapes.
fig8a-equal-center-uptake¶
Packets 1, 2 and 3 (4 kV, 1Px4) produce the same YP uptake in the center of the array, 215 to 235 a.u. (Results, 3.4, paragraph 2; Figure 8A)
Criterion: the mean over dishes of the center ROI per packet, in a.u., each within 193.5 to 258.5 a.u.; no pair of packets differs by a two-sided Student's t-test at p < 0.05. Reproduced when both hold; partially reproduced when no pair differs but a value is out of range.
def ttest(a, b):
return stats.ttest_ind(a, b, equal_var=True)
def center_outcome(g, labels):
vals = {p: g[g.packet == q].center_au.values for p, q in labels.items()}
means = {p: v.mean() for p, v in vals.items()}
in_range = {p: 193.5 <= m <= 258.5 for p, m in means.items()}
ps = {f"{a}v{b}": ttest(vals[a], vals[b]).pvalue for a, b in [(1, 2), (1, 3), (2, 3)]}
differ = any(p < 0.05 for p in ps.values())
out = "reproduced" if all(in_range.values()) and not differ else (
"partially_reproduced" if not differ else "not_reproduced")
return {**{f"P{p}_center": round(m, 1) for p, m in means.items()},
**{f"p_{k}": float(f"{v:.3g}") for k, v in ps.items()}, "outcome": out}
AS_GIVEN = {1: 1, 2: 2, 3: 3}
SWAPPED = {1: 1, 2: 3, 3: 2} # variant: Packet 2 read from the Protocol 3 files and the reverse
c_tab = R.groupby(keys).apply(lambda g: pd.Series(center_outcome(g, AS_GIVEN))).reset_index()
main_row = c_tab[(c_tab.geom == "abs") & (c_tab.n_per_diag == 46) & (c_tab.center == "flank") & (c_tab.bg == "raw") & (c_tab.diag == "both")]
print("Main reading:"); print(main_row.to_string(index=False))
print("\nOutcomes over all", len(c_tab), "readings:"); print(c_tab.outcome.value_counts().to_string())
print("\nBy background reading:"); print(c_tab.groupby("bg").outcome.value_counts().to_string())
print("\nRange of each packet's center mean over readings (a.u.):")
for p in (1, 2, 3):
for b in ("raw", "bgsub"):
v = c_tab[c_tab.bg == b][f"P{p}_center"]
print(f" Packet {p}, {b}: {v.min():.1f} to {v.max():.1f}")
print("\nLargest p-value of each pair over readings:", {k: c_tab[k].max() for k in ["p_1v2", "p_1v3", "p_2v3"]})
print("Smallest p-value of each pair over readings:", {k: c_tab[k].min() for k in ["p_1v2", "p_1v3", "p_2v3"]})
Main reading:
geom n_per_diag center bg diag P1_center P2_center P3_center p_1v2 p_1v3 p_2v3 outcome
abs 46 flank raw both 648.700012 388.5 383.399994 0.000438 0.000928 0.862 not_reproduced
Outcomes over all 48 readings:
outcome
not_reproduced 48
By background reading:
bg outcome
bgsub not_reproduced 24
raw not_reproduced 24
Range of each packet's center mean over readings (a.u.):
Packet 1, raw: 640.6 to 689.1
Packet 1, bgsub: 332.5 to 381.1
Packet 2, raw: 378.5 to 390.3
Packet 2, bgsub: 213.4 to 225.2
Packet 3, raw: 325.0 to 386.3
Packet 3, bgsub: 171.5 to 232.8
Largest p-value of each pair over readings: {'p_1v2': np.float64(0.0467), 'p_1v3': np.float64(0.0259), 'p_2v3': np.float64(0.989)}
Smallest p-value of each pair over readings: {'p_1v2': np.float64(0.000174), 'p_1v3': np.float64(3.6e-05), 'p_2v3': np.float64(0.174)}
print("Per-dish center values, main reading, raw and background subtracted (a.u.):")
print(R[(main | mainbg)].pivot_table(index=["packet", "dish"], columns="bg", values="center_au").round(1).to_string())
print("\nThe same, as a ratio to each dish's own background (center / background):")
cb = R[main].copy(); cb["bg_level"] = [next(d["background"] for d in dishes if d["subject_id"] == s) for s in cb.dish]
cb["center_over_bg"] = cb.center_au / cb.bg_level
print(cb[["dish", "packet", "center_over_bg"]].round(2).to_string(index=False))
for p in (1, 2, 3):
v = cb[cb.packet == p].center_over_bg
print(f" Packet {p}: mean {v.mean():.2f} +/- {v.std(ddof=1)/np.sqrt(len(v)):.2f} s.e.")
print("p (center/background), 1 v 2:", round(ttest(cb[cb.packet == 1].center_over_bg, cb[cb.packet == 2].center_over_bg).pvalue, 4),
"; 1 v 3:", round(ttest(cb[cb.packet == 1].center_over_bg, cb[cb.packet == 3].center_over_bg).pvalue, 4),
"; 2 v 3:", round(ttest(cb[cb.packet == 2].center_over_bg, cb[cb.packet == 3].center_over_bg).pvalue, 4))
Per-dish center values, main reading, raw and background subtracted (a.u.):
bg bgsub raw
packet dish
1 19 294.000000 598.700012
20 390.000000 674.799988
21 338.000000 672.500000
2 12 182.500000 332.899994
13 261.399994 415.000000
14 252.300003 421.200012
2 197.600006 385.000000
3 03 241.300003 413.200012
10 232.399994 382.700012
11 216.100006 354.399994
The same, as a ratio to each dish's own background (center / background):
dish packet center_over_bg
19 1 1.96
20 1 2.37
21 1 2.01
2 2 2.05
12 2 2.21
13 2 2.70
14 2 2.49
03 3 2.40
10 3 2.55
11 3 2.56
Packet 1: mean 2.11 +/- 0.13 s.e.
Packet 2: mean 2.37 +/- 0.14 s.e.
Packet 3: mean 2.50 +/- 0.05 s.e.
p (center/background), 1 v 2: 0.2671 ; 1 v 3: 0.0471 ; 2 v 3: 0.4665
Result: not reproduced. Under the main reading the center means are 648.7 a.u. (Packet 1), 388.5 (Packet 2) and 383.4 (Packet 3), none within 193.5 to 258.5, and Packet 1 differs from Packets 2 and 3 (p = 0.00044 and 0.00093). All 48 readings give not reproduced: Packet 1 is 332.5 to 381.1 a.u. even with the background subtracted, and its p-values against the other two stay below 0.05 (largest 0.047 and 0.026). Packets 2 and 3 do not differ in any reading (smallest p 0.174); with the background subtracted they are 213.4 to 225.2 and 171.5 to 232.8 a.u. The Packet 1 dishes were imaged on another day with a background about twice as high (mean 308.1 a.u., against 161.8 for the dishes of 17 March and 150.3 for that of 18 March); as a ratio to each dish's own background the centers are 2.11, 2.37 and 2.50, and only Packet 1 against Packet 3 differs (p = 0.047). That ratio is not the paper's measure and does not decide the outcome.
fig8a-more-uptake-near-electrodes¶
Packets 1 and 3 cause much more YP uptake near the electrodes (250 to 350 a.u.) than Packet 2 (about 120 a.u.), p < 0.01. (Results, 3.4, paragraph 2; Figure 8A)
Criterion: near the electrodes is the most peripheral ROI; readings (a) the mean of the ROIs closest to each electrode, (b) the ROI closest to e1; one diagonal or both. Reproduced when, under one reading, Packets 1 and 3 lie in 225 to 385 a.u., Packet 2 in 108 to 132 a.u., and each of Packets 1 and 3 exceeds Packet 2 by a two-sided Student's t-test at p < 0.01; partially reproduced when the tests hold and the values do not. The files do not say which electrode is e1, so reading (b) is run below in three ways: with each image label as e1 in every dish, with the end most favourable to the claim in each dish, and over every assignment of e1.
def near_outcome(g, labels):
vals = {p: g[g.packet == q].periph_au.values for p, q in labels.items()}
means = {p: v.mean() for p, v in vals.items()}
rng = {1: (225, 385), 2: (108, 132), 3: (225, 385)}
in_range = all(rng[p][0] <= means[p] <= rng[p][1] for p in (1, 2, 3))
t13 = {p: ttest(vals[p], vals[2]) for p in (1, 3)}
tests = all(t13[p].pvalue < 0.01 and t13[p].statistic > 0 for p in (1, 3))
out = "reproduced" if (in_range and tests) else ("partially_reproduced" if tests else "not_reproduced")
return {**{f"P{p}_periph": round(m, 1) for p, m in means.items()},
"t_1v2": round(t13[1].statistic, 2), "p_1v2": float(f"{t13[1].pvalue:.3g}"),
"t_3v2": round(t13[3].statistic, 2), "p_3v2": float(f"{t13[3].pvalue:.3g}"), "outcome": out}
n_tab = R.groupby(keys).apply(lambda g: pd.Series(near_outcome(g, AS_GIVEN))).reset_index()
print("Main reading:")
print(n_tab[(n_tab.geom == "abs") & (n_tab.n_per_diag == 46) & (n_tab.center == "flank") & (n_tab.bg == "raw") & (n_tab.diag == "both")].to_string(index=False))
print("\nOutcomes over all", len(n_tab), "readings:"); print(n_tab.outcome.value_counts().to_string())
print("\nReadings in which Packet 3 is above Packet 2 near the electrodes:", int((n_tab.P3_periph > n_tab.P2_periph).sum()), "of", len(n_tab))
print("Range over readings (a.u.):", {f"P{p}": (n_tab[f'P{p}_periph'].min(), n_tab[f'P{p}_periph'].max()) for p in (1, 2, 3)})
print("t (Packet 3 minus Packet 2) over readings:", n_tab.t_3v2.min(), "to", n_tab.t_3v2.max())
print("t (Packet 1 minus Packet 2) over readings:", n_tab.t_1v2.min(), "to", n_tab.t_1v2.max(), "; p", n_tab.p_1v2.min(), "to", n_tab.p_1v2.max())
Main reading:
geom n_per_diag center bg diag P1_periph P2_periph P3_periph t_1v2 p_1v2 t_3v2 p_3v2 outcome
abs 46 flank raw both 806.9 480.9 296.2 5.31 0.00316 -3.88 0.0117 not_reproduced
Outcomes over all 48 readings:
outcome
not_reproduced 48
Readings in which Packet 3 is above Packet 2 near the electrodes: 0 of 48
Range over readings (a.u.): {'P1': (np.float64(423.0), np.float64(1130.5)), 'P2': (np.float64(216.8), np.float64(575.0)), 'P3': (np.float64(100.1), np.float64(337.9))}
t (Packet 3 minus Packet 2) over readings: -5.3 to -1.84
t (Packet 1 minus Packet 2) over readings: 1.61 to 10.35 ; p 0.000145 to 0.168
ends = []
for sid, p, bgv, pr in profs:
sv = cancan.summary_values(pr[("abs", 46, "flank")])
ends.append({"dish": sid, "packet": p, **{f"E{k}": round(sv[f"end_e{k}"], 1) for k in range(4)}})
ends = pd.DataFrame(ends)
print("Most peripheral ROI next to each electrode (main geometry, raw a.u.); E0 top left, then clockwise in the image:")
print(ends.to_string(index=False))
print("\nPer-packet mean of each electrode's end value:")
print(ends.groupby("packet")[["E0", "E1", "E2", "E3"]].mean().round(1).to_string())
Most peripheral ROI next to each electrode (main geometry, raw a.u.); E0 top left, then clockwise in the image:
dish packet E0 E1 E2 E3
19 1 774.400024 1285.300049 874.500000 632.700012
20 1 683.500000 708.299988 898.200012 629.099976
21 1 475.299988 952.000000 754.900024 1014.299988
2 2 518.099976 355.799988 213.899994 365.700012
12 2 448.000000 639.000000 567.799988 448.200012
13 2 357.299988 691.900024 529.500000 561.200012
14 2 435.500000 631.500000 475.799988 454.399994
03 3 222.300003 336.200012 255.000000 393.299988
10 3 255.600006 334.700012 288.600006 345.799988
11 3 238.500000 281.899994 267.100006 335.299988
Per-packet mean of each electrode's end value:
E0 E1 E2 E3
packet
1 644.4 981.9 842.5 758.7
2 439.7 579.6 446.7 457.4
3 238.8 317.6 270.2 358.1
# Reading (b): the ROI next to e1. The files do not say which electrode is e1, and E0 to E3 are
# assigned by angle in each image, so two versions are run: (i) the same label E0, E1, E2 or E3
# taken as e1 in every dish; (ii) a bound, taking in each dish the end most favourable to the
# claim (the lowest end in each Packet 2 dish, the highest in each Packet 1 and Packet 3 dish).
# The center is the mean of both diagonals, or the center of e1's own diagonal (E0, E2 on d1).
e_rows = []
for sid, p, bgv, pr_ in profs:
for (geom, n, center), pr in pr_.items():
sv = cancan.summary_values(pr, roi_at_zero=center == "at0")
for b in ("raw", "bgsub"):
off = bgv if b == "bgsub" else 0.0
for k in range(4):
for cdiag in ("both", "own"):
c = ((sv["center_d1"] + sv["center_d2"]) / 2 if cdiag == "both"
else sv["center_d1"] if k in (0, 2) else sv["center_d2"])
e_rows.append({"dish": sid, "packet": p, "geom": geom, "n_per_diag": n, "center": center,
"bg": b, "cdiag": cdiag, "e1": f"E{k}", "center_au": c - off,
"periph_au": sv[f"end_e{k}"] - off})
E1 = pd.DataFrame(e_rows); E1["pct"] = 100 * E1.center_au / E1.periph_au
ekeys = ["geom", "n_per_diag", "center", "bg", "cdiag"]
def favourable(df, col, low_for):
"""Per dish and reading, the e1 row with the lowest `col` in packets `low_for`, else the highest."""
idx = [g[col].idxmin() if g.packet.iloc[0] in low_for else g[col].idxmax()
for _, g in df.groupby(ekeys + ["dish"])]
return df.loc[idx]
n_e1 = E1[E1.cdiag == "both"].groupby(ekeys + ["e1"]).apply(lambda g: pd.Series(near_outcome(g, AS_GIVEN))).reset_index()
print("(i) Same label as e1 in every dish:", len(n_e1), "readings")
print(n_e1.outcome.value_counts().to_string())
print("Readings with Packet 3 above Packet 2:", int((n_e1.P3_periph > n_e1.P2_periph).sum()), "of", len(n_e1))
print("Main geometry, raw:")
print(n_e1[(n_e1.geom == "abs") & (n_e1.n_per_diag == 46) & (n_e1.center == "flank") & (n_e1.bg == "raw")].to_string(index=False))
fav = favourable(E1[E1.cdiag == "both"], "periph_au", low_for={2})
n_fav = fav.groupby(ekeys).apply(lambda g: pd.Series(near_outcome(g, AS_GIVEN))).reset_index()
print("\n(ii) Most favourable end in each dish:", len(n_fav), "readings")
print(n_fav.outcome.value_counts().to_string())
print("Readings with Packet 3 above Packet 2:", int((n_fav.P3_periph > n_fav.P2_periph).sum()), "of", len(n_fav))
print("t (Packet 3 minus Packet 2):", n_fav.t_3v2.min(), "to", n_fav.t_3v2.max(), "; Packet 2 range (a.u.):",
n_fav.P2_periph.min(), "to", n_fav.P2_periph.max())
print(n_fav.to_string(index=False))
# (iii) Every assignment of e1 to one of the four electrodes in each dish (4^10 = 1,048,576), with
# the two-sided Student's t-tests computed for each assignment.
import itertools
combos = np.array(list(itertools.product(range(4), repeat=10)), dtype=np.int8)
def t_equal_var(a, b):
na, nb_ = a.shape[1], b.shape[1]
sp = ((na - 1) * a.var(1, ddof=1) + (nb_ - 1) * b.var(1, ddof=1)) / (na + nb_ - 2)
t = (a.mean(1) - b.mean(1)) / np.sqrt(sp * (1 / na + 1 / nb_))
return t, 2 * stats.t.sf(np.abs(t), na + nb_ - 2)
n_enum = []
for key, g in E1[E1.cdiag == "both"].groupby(ekeys):
g = g.sort_values(["packet", "dish", "e1"])
pk = g.drop_duplicates("dish").packet.values
Pe = g.periph_au.values.reshape(10, 4)[np.arange(10), combos].astype(float)
v = {p: Pe[:, pk == p] for p in (1, 2, 3)}
tests = np.ones(len(combos), bool)
for p in (1, 3):
t, pv = t_equal_var(v[p], v[2]); tests &= (t > 0) & (pv < 0.01)
rng = {1: (225, 385), 2: (108, 132), 3: (225, 385)}
inr = np.all([(rng[p][0] <= v[p].mean(1)) & (v[p].mean(1) <= rng[p][1]) for p in (1, 2, 3)], axis=0)
n_enum.append({**dict(zip(ekeys, key)), "frac_reproduced": (tests & inr).mean(), "frac_partial": (tests & ~inr).mean(),
"frac_P3_above_P2": (v[3].mean(1) > v[2].mean(1)).mean()})
n_enum = pd.DataFrame(n_enum)
print("\n(iii) All", len(combos), "assignments of e1, in", len(n_enum), "readings:")
print(n_enum.round(4).to_string(index=False))
(i) Same label as e1 in every dish: 64 readings outcome not_reproduced 64 Readings with Packet 3 above Packet 2: 0 of 64 Main geometry, raw: geom n_per_diag center bg cdiag e1 P1_periph P2_periph P3_periph t_1v2 p_1v2 t_3v2 p_3v2 outcome abs 46 flank raw both E0 644.400024 439.799988 238.800003 2.45 0.0583 -5.05 0.00393 not_reproduced abs 46 flank raw both E1 981.900024 579.599976 317.600006 2.42 0.0601 -2.88 0.03450 not_reproduced abs 46 flank raw both E2 842.500000 446.799988 270.299988 3.90 0.0114 -1.86 0.12200 not_reproduced abs 46 flank raw both E3 758.700012 457.399994 358.100006 2.58 0.0496 -2.00 0.10200 not_reproduced (ii) Most favourable end in each dish: 16 readings outcome not_reproduced 16 Readings with Packet 3 above Packet 2: 10 of 16 t (Packet 3 minus Packet 2): -0.5 to 0.71 ; Packet 2 range (a.u.): 166.6 to 377.5 geom n_per_diag center bg cdiag P1_periph P2_periph P3_periph t_1v2 p_1v2 t_3v2 p_3v2 outcome abs 46 at0 bgsub both 888.099976 173.800003 191.199997 4.39 0.007120 0.26 0.808 not_reproduced abs 46 at0 raw both 1196.199951 338.899994 344.700012 5.35 0.003060 0.10 0.927 not_reproduced abs 46 flank bgsub both 757.900024 198.600006 204.600006 4.71 0.005270 0.08 0.937 not_reproduced abs 46 flank raw both 1065.900024 363.700012 358.100006 6.10 0.001720 -0.09 0.936 not_reproduced abs 92 at0 bgsub both 1021.299988 178.500000 230.199997 6.75 0.001080 0.71 0.507 not_reproduced abs 92 at0 raw both 1329.400024 343.500000 383.700012 8.00 0.000493 0.62 0.564 not_reproduced abs 92 flank bgsub both 853.500000 190.399994 211.899994 4.08 0.009490 0.25 0.814 not_reproduced abs 92 flank raw both 1161.500000 355.500000 365.399994 4.97 0.004210 0.13 0.904 not_reproduced rel 46 at0 bgsub both 806.500000 208.399994 185.399994 5.86 0.002060 -0.30 0.779 not_reproduced rel 46 at0 raw both 1114.500000 373.500000 338.899994 7.60 0.000625 -0.50 0.636 not_reproduced rel 46 flank bgsub both 814.700012 212.399994 202.600006 4.88 0.004540 -0.17 0.872 not_reproduced rel 46 flank raw both 1122.699951 377.500000 356.100006 6.31 0.001470 -0.41 0.699 not_reproduced rel 92 at0 bgsub both 951.500000 208.199997 214.100006 5.46 0.002810 0.06 0.952 not_reproduced rel 92 at0 raw both 1259.500000 373.299988 367.600006 6.59 0.001210 -0.07 0.948 not_reproduced rel 92 flank bgsub both 903.900024 166.600006 215.199997 9.51 0.000217 0.67 0.533 not_reproduced rel 92 flank raw both 1211.900024 331.700012 368.700012 13.63 0.000038 0.60 0.574 not_reproduced
(iii) All 1048576 assignments of e1, in 16 readings: geom n_per_diag center bg cdiag frac_reproduced frac_partial frac_P3_above_P2 abs 46 at0 bgsub both 0.0 0.0 0.0007 abs 46 at0 raw both 0.0 0.0 0.0002 abs 46 flank bgsub both 0.0 0.0 0.0004 abs 46 flank raw both 0.0 0.0 0.0000 abs 92 at0 bgsub both 0.0 0.0 0.0016 abs 92 at0 raw both 0.0 0.0 0.0006 abs 92 flank bgsub both 0.0 0.0 0.0006 abs 92 flank raw both 0.0 0.0 0.0002 rel 46 at0 bgsub both 0.0 0.0 0.0000 rel 46 at0 raw both 0.0 0.0 0.0000 rel 46 flank bgsub both 0.0 0.0 0.0000 rel 46 flank raw both 0.0 0.0 0.0000 rel 92 at0 bgsub both 0.0 0.0 0.0002 rel 92 at0 raw both 0.0 0.0 0.0000 rel 92 flank bgsub both 0.0 0.0 0.0035 rel 92 flank raw both 0.0 0.0 0.0015
Result: not reproduced. Near the electrodes the Protocol 3 files show the least uptake, where the claim has Packet 3 above Packet 2: under the main reading 806.9 a.u. (Packet 1), 480.9 (Packet 2) and 296.2 (Packet 3), with Packet 3 below Packet 2 (t = -3.88, p = 0.012). Packet 3 is below Packet 2 in all 48 readings (t from -5.3 to -1.84). Packet 1 exceeds Packet 2 in every reading (t from 1.61 to 10.35, p from 0.00015 to 0.168).
Reading (b), the ROI next to e1, cannot be decided, as the files do not say which electrode is e1. Taking the same image label as e1 in every dish gives not reproduced in all 64 readings, with Packet 3 below Packet 2 in each. Taking in each dish the end most favourable to the claim (the lowest end of each Packet 2 dish, the highest of each Packet 1 and 3 dish) still gives not reproduced in all 16 readings: Packet 2 is 166.6 to 377.5 a.u., never in 108 to 132, and Packet 3 is not above Packet 2 by the test (t from -0.5 to 0.71). Over all 1,048,576 assignments of e1 to the four electrodes of each dish, none gives reproduced or partially reproduced, and at most 0.35% put the Packet 3 mean above the Packet 2 mean. No choice of e1 reproduces the claim, even in part.
fig8b-packet2-best-targeting¶
Normalized to the most peripheral ROI, YP emission in the center is 183% for Packet 2 and drops to 83% for Packet 3 and 74% for Packet 1. (Results, 3.4, paragraph 2; Figure 8B)
Criterion: normalization per profile to its most peripheral ROI; readings (a) normalize each dish then average, (b) normalize the mean profile; peripheral ROI and diagonals as in the previous claim. Reproduced when all three values are within tolerance (164.7 to 201.3%, 74.7 to 91.3%, 66.6 to 81.4%) under one reading; partially reproduced when Packet 2 is the highest and above 100% while Packets 1 and 3 are below 100%, without all values in tolerance.
TOL = {2: (164.7, 201.3), 3: (74.7, 91.3), 1: (66.6, 81.4)}
def pct_outcome(g, labels, norm):
v = {}
for p, q in labels.items():
s = g[g.packet == q]
v[p] = s.pct.mean() if norm == "per_dish" else 100 * s.center_au.mean() / s.periph_au.mean()
allin = all(TOL[p][0] <= v[p] <= TOL[p][1] for p in (1, 2, 3))
partial = v[2] > max(v[1], v[3]) and v[2] > 100 and v[1] < 100 and v[3] < 100
out = "reproduced" if allin else ("partially_reproduced" if partial else "not_reproduced")
return {"P1_pct": round(v[1], 1), "P2_pct": round(v[2], 1), "P3_pct": round(v[3], 1), "outcome": out}
b_rows = []
for norm in ("per_dish", "mean_profile"):
t = R.groupby(keys).apply(lambda g: pd.Series(pct_outcome(g, AS_GIVEN, norm))).reset_index()
t.insert(0, "norm", norm); b_rows.append(t)
b_tab = pd.concat(b_rows, ignore_index=True)
print("Main reading, both normalizations:")
print(b_tab[(b_tab.geom == "abs") & (b_tab.n_per_diag == 46) & (b_tab.center == "flank") & (b_tab.bg == "raw") & (b_tab.diag == "both")].to_string(index=False))
print("\nOutcomes over all", len(b_tab), "readings:"); print(b_tab.outcome.value_counts().to_string())
print("\nRange over readings (%):", {c: (b_tab[c].min(), b_tab[c].max()) for c in ["P1_pct", "P2_pct", "P3_pct"]})
print("Readings in which Packet 2 has the highest normalized center:", int((b_tab.P2_pct > b_tab[["P1_pct", "P3_pct"]].max(axis=1)).sum()), "of", len(b_tab))
print("Readings in which Packet 3 has the highest normalized center:", int((b_tab.P3_pct > b_tab[["P1_pct", "P2_pct"]].max(axis=1)).sum()), "of", len(b_tab))
print("\nPer-dish normalized center, main reading and background subtracted (%):")
print(R[main | mainbg].pivot_table(index=["packet", "dish"], columns="bg", values="pct").round(1).to_string())
pm = R[main]
print("\nTwo-sided t-test of per-dish % (main reading), Protocol 3 v Protocol 2:",
ttest(pm[pm.packet == 3].pct, pm[pm.packet == 2].pct))
Main reading, both normalizations:
norm geom n_per_diag center bg diag P1_pct P2_pct P3_pct outcome
per_dish abs 46 flank raw both 81.3 82.8 129.4 not_reproduced
mean_profile abs 46 flank raw both 80.4 80.8 129.5 not_reproduced
Outcomes over all 96 readings:
outcome
not_reproduced 96
Range over readings (%): {'P1_pct': (np.float64(44.1), np.float64(98.3)), 'P2_pct': (np.float64(52.5), np.float64(123.7)), 'P3_pct': (np.float64(102.1), np.float64(244.3))}
Readings in which Packet 2 has the highest normalized center: 0 of 96
Readings in which Packet 3 has the highest normalized center: 96 of 96
Per-dish normalized center, main reading and background subtracted (%):
bg bgsub raw
packet dish
1 19 50.1 67.1
20 87.6 92.5
21 72.7 84.2
2 12 48.6 63.3
13 68.5 77.6
14 76.4 84.4
2 112.3 105.9
3 03 185.9 137.0
10 149.1 125.0
11 151.8 126.3
Two-sided t-test of per-dish % (main reading), Protocol 3 v Protocol 2: TtestResult(statistic=np.float64(4.24887305826534), pvalue=np.float64(0.00810053154862612), df=np.float64(5.0))
# Reading (b) of the normalization: each profile normalized to the ROI next to e1, as in the
# previous claim: (i) the same label as e1 in every dish; (ii) in each dish the end most favourable
# to the claim (the lowest end in each Packet 2 dish, the highest in each Packet 1 and 3 dish).
b_e1 = pd.concat([E1.groupby(ekeys + ["e1"]).apply(lambda g: pd.Series(pct_outcome(g, AS_GIVEN, norm))).reset_index().assign(norm=norm)
for norm in ("per_dish", "mean_profile")], ignore_index=True)
print("(i) Same label as e1 in every dish:", len(b_e1), "readings")
print(b_e1.outcome.value_counts().to_string())
print("Readings in which Packet 2 has the highest normalized center:",
int((b_e1.P2_pct > b_e1[["P1_pct", "P3_pct"]].max(axis=1)).sum()), "of", len(b_e1))
print("Range (%):", {c: (b_e1[c].min(), b_e1[c].max()) for c in ["P1_pct", "P2_pct", "P3_pct"]})
fav_b = favourable(E1, "periph_au", low_for={2})
b_fav = pd.concat([fav_b.groupby(ekeys).apply(lambda g: pd.Series(pct_outcome(g, AS_GIVEN, norm))).reset_index().assign(norm=norm)
for norm in ("per_dish", "mean_profile")], ignore_index=True)
print("\n(ii) Most favourable end in each dish:", len(b_fav), "readings")
print(b_fav.outcome.value_counts().to_string())
print("Readings in which Packet 2 has the highest normalized center:",
int((b_fav.P2_pct > b_fav[["P1_pct", "P3_pct"]].max(axis=1)).sum()), "of", len(b_fav))
print("Range (%):", {c: (b_fav[c].min(), b_fav[c].max()) for c in ["P1_pct", "P2_pct", "P3_pct"]})
print(b_fav[(b_fav.geom == "abs") & (b_fav.n_per_diag == 46) & (b_fav.center == "flank")].to_string(index=False))
# (iii) Every assignment of e1 to one of the four electrodes in each dish (4^10 = 1,048,576), with
# the fraction of assignments that would give each outcome.
import itertools
combos = np.array(list(itertools.product(range(4), repeat=10)), dtype=np.int8)
enum = []
for key, g in E1.groupby(ekeys):
g = g.sort_values(["packet", "dish", "e1"])
pk = g.drop_duplicates("dish").packet.values
C = g.center_au.values.reshape(10, 4)[np.arange(10), combos]
Pe = g.periph_au.values.reshape(10, 4)[np.arange(10), combos]
for norm in ("per_dish", "mean_profile"):
v = {p: (100 * C[:, pk == p] / Pe[:, pk == p]).mean(1) if norm == "per_dish"
else 100 * C[:, pk == p].mean(1) / Pe[:, pk == p].mean(1) for p in (1, 2, 3)}
allin = np.all([(TOL[p][0] <= v[p]) & (v[p] <= TOL[p][1]) for p in (1, 2, 3)], axis=0)
partial = (v[2] > np.maximum(v[1], v[3])) & (v[2] > 100) & (v[1] < 100) & (v[3] < 100)
enum.append({**dict(zip(ekeys, key)), "norm": norm, "frac_reproduced": allin.mean(),
"frac_partial": (partial & ~allin).mean()})
enum = pd.DataFrame(enum)
print("\n(iii) All", len(combos), "assignments of e1, in", len(enum), "readings:")
print("Largest fraction giving reproduced:", enum.frac_reproduced.max())
print("Largest fraction giving partially reproduced, by ROIs per diagonal:")
print(enum.groupby("n_per_diag").frac_partial.max().round(4).to_string())
(i) Same label as e1 in every dish:
256 readings
outcome
not_reproduced 256
Readings in which Packet 2 has the highest normalized center: 10 of 256
Range (%): {'P1_pct': (np.float32(35.5), np.float32(281.3)), 'P2_pct': (np.float32(-668.8), np.float32(540.7)), 'P3_pct': (np.float32(101.5), np.float32(309.4))}
(ii) Most favourable end in each dish: 64 readings
outcome
not_reproduced 37
partially_reproduced 27
Readings in which Packet 2 has the highest normalized center: 47 of 64
Range (%): {'P1_pct': (np.float32(35.5), np.float32(62.7)), 'P2_pct': (np.float32(-661.8), np.float32(572.3)), 'P3_pct': (np.float32(75.7), np.float32(118.8))}
geom n_per_diag center bg cdiag P1_pct P2_pct P3_pct outcome norm
abs 46 flank bgsub both 47.799999 257.399994 112.500000 not_reproduced per_dish
abs 46 flank bgsub own 47.599998 245.100006 113.900002 not_reproduced per_dish
abs 46 flank raw both 62.700001 116.800003 107.099998 not_reproduced per_dish
abs 46 flank raw own 62.599998 115.500000 108.000000 not_reproduced per_dish
abs 46 flank bgsub both 44.900002 112.500000 112.400002 not_reproduced mean_profile
abs 46 flank bgsub own 44.900002 111.599998 113.699997 not_reproduced mean_profile
abs 46 flank raw both 60.900002 106.800003 107.099998 not_reproduced mean_profile
abs 46 flank raw own 60.799999 106.300003 107.800003 not_reproduced mean_profile
(iii) All 1048576 assignments of e1, in 64 readings: Largest fraction giving reproduced: 0.0 Largest fraction giving partially reproduced, by ROIs per diagonal: n_per_diag 46 0.0000 92 0.0146
Result: not reproduced. Under the main reading the normalized centers are 81.3% (Packet 1), 82.8% (Packet 2) and 129.4% (Packet 3), where the claim has 74, 183 and 83%. Packet 3, not Packet 2, has the highest normalized center in all 96 readings, and Packet 2 is 52.5 to 123.7%. Per dish, the Protocol 3 files exceed the Protocol 2 files (t = 4.25, p = 0.008).
Reading (b), normalizing to the ROI next to e1, cannot be decided, as the files do not say which electrode is e1. The same image label as e1 in every dish gives not reproduced in all 256 readings. Choosing e1 in each dish to suit the claim (the lowest end in each Packet 2 dish, the highest in each Packet 1 and 3 dish) gives partially reproduced in 27 of 64 readings, all with 92 ROIs per diagonal, where a single end ROI is narrow and noisy; with the background subtracted, an end ROI near the background level makes a per-dish percentage unstable (Packet 2 from -661.8 to 572.3%). Over all 1,048,576 assignments of e1 to the four electrodes of each dish, none gives reproduced, none gives partially reproduced with 46 ROIs per diagonal, and at most 1.5% do with 92. The outcome is decided by reading (a), which can be run; that a small share of the possible e1 assignments meets the partial condition is recorded, not taken as the outcome.
Variant: Protocol 2 and Protocol 3 exchanged¶
The files labelled Protocol 3 (8 phases with 5-us delays) show the profile the paper reports for Packet 2, and the reverse. This variant reruns the three criteria with Packet 2 read from the Protocol 3 files and Packet 3 from the Protocol 2 files. It does not decide any outcome: the outcomes above test the dataset as labelled.
c_sw = R.groupby(keys).apply(lambda g: pd.Series(center_outcome(g, SWAPPED))).reset_index()
n_sw = R.groupby(keys).apply(lambda g: pd.Series(near_outcome(g, SWAPPED))).reset_index()
b_sw = pd.concat([R.groupby(keys).apply(lambda g: pd.Series(pct_outcome(g, SWAPPED, norm))).reset_index().assign(norm=norm)
for norm in ("per_dish", "mean_profile")], ignore_index=True)
for name, t in [("fig8a-equal-center-uptake", c_sw), ("fig8a-more-uptake-near-electrodes", n_sw), ("fig8b-packet2-best-targeting", b_sw)]:
print(f"{name}, labels exchanged: outcomes over {len(t)} readings")
print(t.groupby("bg").outcome.value_counts().to_string()); print()
sel = (n_sw.geom == "abs") & (n_sw.n_per_diag == 46) & (n_sw.center == "flank") & (n_sw.diag == "both")
print("Near electrodes, main geometry, labels exchanged:"); print(n_sw[sel].to_string(index=False))
sel = (b_sw.geom == "abs") & (b_sw.n_per_diag == 46) & (b_sw.center == "flank") & (b_sw.diag == "both")
print("\nNormalized center, main geometry, labels exchanged:"); print(b_sw[sel].to_string(index=False))
print("\nNormalized center over all readings, labels exchanged, by background (%):")
print(b_sw.groupby("bg")[["P1_pct", "P2_pct", "P3_pct"]].agg(["min", "max"]).round(1).to_string())
fig8a-equal-center-uptake, labels exchanged: outcomes over 48 readings
bg outcome
bgsub not_reproduced 24
raw not_reproduced 24
fig8a-more-uptake-near-electrodes, labels exchanged: outcomes over 48 readings
bg outcome
bgsub not_reproduced 24
raw not_reproduced 19
partially_reproduced 5
fig8b-packet2-best-targeting, labels exchanged: outcomes over 96 readings
bg outcome
bgsub partially_reproduced 45
reproduced 2
not_reproduced 1
raw partially_reproduced 47
not_reproduced 1
Near electrodes, main geometry, labels exchanged:
geom n_per_diag center bg diag P1_periph P2_periph P3_periph t_1v2 p_1v2 t_3v2 p_3v2 outcome
abs 46 flank bgsub both 498.8 142.7 315.8 7.90 0.001390 3.03 0.0290 not_reproduced
abs 46 flank raw both 806.9 296.2 480.9 10.74 0.000427 3.88 0.0117 not_reproduced
Normalized center, main geometry, labels exchanged:
geom n_per_diag center bg diag P1_pct P2_pct P3_pct outcome norm
abs 46 flank bgsub both 70.2 162.3 76.4 partially_reproduced per_dish
abs 46 flank raw both 81.3 129.4 82.8 partially_reproduced per_dish
abs 46 flank bgsub both 68.3 161.1 70.8 partially_reproduced mean_profile
abs 46 flank raw both 80.4 129.5 80.8 partially_reproduced mean_profile
Normalized center over all readings, labels exchanged, by background (%):
P1_pct P2_pct P3_pct
min max min max min max
bg
bgsub 44.1 98.3 104.1 244.3 52.5 123.7
raw 59.3 93.3 102.1 150.2 66.1 103.6
Result: with the labels exchanged, the order of Figure 8 appears. In the main geometry the normalized centers are 70.2, 162.3 and 76.4% (Packets 1, 2, 3) with the background subtracted and 81.3, 129.4 and 82.8% raw, so fig8b-packet2-best-targeting would be partially reproduced in 92 of 96 readings and reproduced in 2. fig8a-more-uptake-near-electrodes would still not be reproduced in 43 of 48 readings (with the background subtracted, Packet 1 is 498.8 a.u., above 385, and Packet 3 over Packet 2 gives p = 0.029, not below 0.01), and fig8a-equal-center-uptake in all 48, as Packet 1 is unchanged. The data alone cannot say whether the files' labels or the paper's are exchanged; this variant does not change the outcomes above.
Summary¶
| Claim | Outcome |
|---|---|
| fig4-unipolar-strongest, fig4-bipolar-weaker, fig4-sham-no-uptake, fig4-twenty-fold-cancellation | Blocked, data missing |
| fig5a-trains-stronger-throughout, fig5a-packet-center-peak-e1-stronger, fig5bc-rotation-focuses-center | Blocked, data missing |
| fig6a-more-packets-more-uptake, fig6b-center-to-periphery-by-number, fig6-peak-narrower | Blocked, data missing |
| fig7a-amplitude-scales-uniformly, fig7bd-targeting-independent-of-amplitude, array-size-independence | Blocked, data missing |
| fig8a-equal-center-uptake | Not reproduced |
| fig8a-more-uptake-near-electrodes | Not reproduced |
| fig8b-packet2-best-targeting | Not reproduced |
| fig9-gap-fades-periphery-faster, fig9b-gap-boosts-targeting | Blocked, data missing |
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 paper measured ROIs by hand in MetaMorph on the full-resolution images; here the electrodes are found automatically, ROIs are squares placed along both diagonals of the reduced image (8 x 8 pixel means), and the scale of 0.65 um per pixel is taken from the camera and objective, as the files give none. The ROI size, whether a background was subtracted, which diagonal, and which ROI is the center are not given by the paper; each is run every way, and no outcome depends on them. Reading (b) of the last two claims, the ROI next to e1, cannot be decided, as the files do not identify e1; it was run with each electrode label as e1 and with every assignment of e1, and no assignment reproduces either claim, though at most 1.5% of the assignments give fig8b-packet2-best-targeting the partial condition when 92 ROIs are placed per diagonal.
Two points are unresolved. First, the fig8a-equal-center-uptake failure rests largely on the Packet 1 dishes, which were imaged on another day with about twice the background of the others; a gain or exposure difference would scale all values and would not be removed by subtracting a background, and as a ratio to each dish's background Packet 1 does not differ from Packet 2 (p = 0.27) and differs from Packet 3 at p = 0.047. The files record no exposure or gain, so this cannot be settled. Second, the files labelled Protocol 2 and 3 behave as the paper's Packets 3 and 2. Each file's description defines the protocols as the paper defines the packets, so the data are tested as labelled, but a labelling error in the deposit or in the figure cannot be excluded.
What was found about the dataset. It holds only the Figure 8 experiment (10 dishes), not the
single-cell time-lapse of Figure 4 or the monolayer experiments of Figures 5, 6, 7 and 9. The files
labelled Protocol 3 (8 phases with 5-us delays) show the strongest central focusing and the least
uptake near the electrodes, which the paper reports for Packet 2 (8 phases without delays), and the
files labelled Protocol 2 behave as the paper's Packet 3; either the labels or the paper's
assignment appear exchanged. The files give no pixel size (resolution 96, or 1.0 in dish 10), no
position or identity of the electrodes, and no exposure or gain, and the Packet 1 dishes, imaged on
another day, have about twice the background of the others. Dish 10 names no device. The species
summary names two species for one cell line, and the age (P1D) does not apply to a cell line.