#!/usr/bin/env python3 """ nta_sim.py - How faithfully does nanoparticle tracking analysis (NTA) recover an extracellular-vesicle (EV) size distribution? Monte Carlo of the whole NTA measurement chain: ground truth -> optical detection -> 2-D Brownian tracks -> minimum-track-length filter -> per-track D -> Stokes-Einstein diameter Ground truth : mixture of lognormals (small EVs + larger vesicles). Detection : logistic P(detect | d). Scattering scales ~ d^6, so small, low-refractive-index EVs fall below the camera threshold. Physics : Stokes-Einstein D = kT / (3 pi eta d). Each frame the particle steps N(0, 2 D dt) per axis; every recorded position carries Gaussian localization error. Tracks : geometric track lengths; tracks shorter than MIN_TRACK_FRAMES are discarded, as NTA software does. Sizing : D from an ordinary-least-squares fit of the time-averaged MSD at lags 1-4 with a free intercept (the intercept absorbs the 4*sigma^2 localization offset). The covariance-based estimator (CVE; Vestergaard, Blainey & Flyvbjerg, Phys Rev E 89:022726, 2014) is computed alongside for comparison. D <= 0 cannot be converted to a size: those tracks are rejected and counted, as are sizes outside the instrument range SIZE_RANGE_NM. Not modelled: motion blur, particles leaving the focal depth (which shortens tracks of small, fast EVs), refractive-index spread, drift, camera-level and threshold settings, and FTLA-style post-processing done by NanoSight software. Outputs: nta_simulation.png next to this file + a summary table on stdout. Runs in about a second. Needs numpy, scipy, matplotlib. The helpers diameter_to_D, D_to_diameter, simulate_tracks, estimate_D_msd, estimate_D_cve and summarize are self-contained and unit-explicit (nm, s, nm^2/s) so they can be lifted straight into VesiScope. """ import time from pathlib import Path import numpy as np from scipy.stats import gaussian_kde import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.lines import Line2D from matplotlib.patches import Circle, Patch from matplotlib.ticker import FixedLocator, FuncFormatter, MultipleLocator, NullLocator # ============================================================================= # PARAMETERS (edit here) # ============================================================================= SEED = 2026 # Ground truth: lognormal mixture, (median diameter nm, geometric SD, weight) N_PARTICLES = 4000 COMPONENTS = ( (90.0, 1.35, 0.75), # small EVs (190.0, 1.25, 0.25), # larger vesicles ) # Optical detection: logistic in diameter, P = 1 / (1 + exp(-(d - D50) / W)) DET_D50_NM = 65.0 # 50 % detected DET_WIDTH_NM = 15.0 / np.log(9.0) # 6.8 nm -> 10 % at 50 nm, 90 % at 80 nm # Physics: water at 25 C KB_J_PER_K = 1.380649e-23 TEMP_K = 298.15 VISCOSITY_PA_S = 0.890e-3 # Camera and tracking FPS = 30.0 # frame rate; dt = 1 / FPS SIGMA_LOC_NM = 25.0 # localization error, per axis, per frame MEAN_TRACK_FRAMES = 25.0 # geometric track-length distribution (positions) MIN_TRACK_FRAMES = 10 # shorter tracks are discarded MSD_LAGS = (1, 2, 3, 4) # lags used in the MSD fit SIZE_RANGE_NM = (10.0, 1000.0) # instrument sizing range; outside -> rejected # Reporting and figure SMALL_CUTOFF_NM = 70.0 BINS_NM = np.arange(0.0, 401.0, 10.0) EXAMPLE_SIZES_NM = (70.0, 120.0, 250.0) EXAMPLE_FRAMES = 60 # 2 s window for panel A OUT_PNG = Path(__file__).absolute().with_name("nta_simulation.png") C_TRUE, C_DET, C_REC = "#2a78d6", "#1baf7a", "#eb6834" C_EXAMPLES = ("#a48be0", "#6f4bc4", "#3a1f7a") # one hue, light = small INK, INK_2, INK_3, GRID = "#1f2328", "#4b5563", "#8a919c", "#e6e8eb" # ============================================================================= # PHYSICS # ============================================================================= DT_S = 1.0 / FPS # kT / (3 pi eta), converted from m^3/s to nm^3/s, so that D[nm^2/s] = SE / d[nm] SE_NM3_PER_S = KB_J_PER_K * TEMP_K / (3.0 * np.pi * VISCOSITY_PA_S) * 1e27 def diameter_to_D(d_nm): """Stokes-Einstein: hydrodynamic diameter (nm) -> D (nm^2/s).""" return SE_NM3_PER_S / np.asarray(d_nm, float) def D_to_diameter(D_nm2_s): """Inverse Stokes-Einstein: D (nm^2/s) -> diameter (nm). D <= 0 -> NaN.""" D = np.asarray(D_nm2_s, float) with np.errstate(divide="ignore", invalid="ignore"): return np.where(D > 0, SE_NM3_PER_S / D, np.nan) def detection_probability(d_nm): """Probability that a particle of diameter d (nm) is seen by the camera.""" return 1.0 / (1.0 + np.exp(-(np.asarray(d_nm, float) - DET_D50_NM) / DET_WIDTH_NM)) def sample_population(rng, n, components=COMPONENTS): """Draw n diameters (nm) from the lognormal mixture (exact component counts).""" w = np.array([c[2] for c in components], float) counts = np.floor(w / w.sum() * n).astype(int) counts[0] += n - counts.sum() parts = [np.exp(rng.normal(np.log(med), np.log(gsd), k)) for (med, gsd, _), k in zip(components, counts)] return rng.permutation(np.concatenate(parts)) def simulate_tracks(rng, d_nm, n_frames, sigma_loc_nm=SIGMA_LOC_NM, dt=DT_S): """2-D Brownian tracks with localization error. d_nm: diameters (nm); n_frames: positions per track. Returns observed positions, shape (n_tracks, max_frames, 2), in nm. Frames past each track's end are NaN. """ d_nm = np.atleast_1d(np.asarray(d_nm, float)) n_frames = np.atleast_1d(np.asarray(n_frames, int)) n, fmax = d_nm.size, int(n_frames.max()) step_sd = np.sqrt(2.0 * diameter_to_D(d_nm) * dt) # per axis, nm steps = rng.standard_normal((n, fmax - 1, 2)) * step_sd[:, None, None] true_xy = np.concatenate([np.zeros((n, 1, 2)), np.cumsum(steps, axis=1)], axis=1) obs = true_xy + rng.normal(0.0, sigma_loc_nm, true_xy.shape) obs[np.arange(fmax)[None, :] >= n_frames[:, None]] = np.nan return obs # ============================================================================= # ESTIMATORS # ============================================================================= def time_averaged_msd(xy, lags=MSD_LAGS): """Per-track time-averaged 2-D MSD (nm^2) at each lag (overlapping windows).""" out = np.empty((xy.shape[0], len(lags))) for j, k in enumerate(lags): sq = np.sum((xy[:, k:] - xy[:, :-k]) ** 2, axis=2) # NaN past track end out[:, j] = np.nanmean(sq, axis=1) return out def estimate_D_msd(xy, lags=MSD_LAGS, dt=DT_S): """OLS fit of MSD(t) = 4 D t + b over the given lags, free intercept b. Returns D (nm^2/s) per track. Short, noisy tracks can give D <= 0. """ t = np.asarray(lags, float) * dt tc = t - t.mean() msd = time_averaged_msd(xy, lags) slope = msd @ tc / (tc @ tc) # sum(tc) = 0, so no need to center msd return slope / 4.0 # 2-D: MSD = 4 D t def estimate_D_cve(xy, dt=DT_S): """Covariance-based estimator (Vestergaard et al. 2014), averaged over x, y. Per axis: D = / (2 dt) + / dt. Unbiased under localization error and motion blur; can be <= 0. """ dx = np.diff(xy, axis=1) var = np.nanmean(dx ** 2, axis=1) cov = np.nanmean(dx[:, 1:] * dx[:, :-1], axis=1) return np.mean(var / (2.0 * dt) + cov / dt, axis=1) # ============================================================================= # SUMMARY STATISTICS # ============================================================================= def kde_mode(d_nm, grid=np.arange(0.0, 1000.5, 0.5)): """Mode (nm) = peak of a Gaussian KDE with a robust Silverman bandwidth.""" d = np.asarray(d_nm, float) sd = d.std(ddof=1) iqr = np.subtract(*np.percentile(d, [75, 25])) h = 0.9 * min(sd, iqr / 1.34) * d.size ** -0.2 return grid[np.argmax(gaussian_kde(d, bw_method=h / sd)(grid))] def summarize(d_nm): """count, mean, mode, D10, D50, D90 (nm) and fraction below SMALL_CUTOFF_NM.""" d = np.asarray(d_nm, float) d10, d50, d90 = np.percentile(d, [10, 50, 90]) return dict(count=d.size, mean=d.mean(), mode=kde_mode(d), d10=d10, d50=d50, d90=d90, frac_small=np.mean(d < SMALL_CUTOFF_NM)) def size_ratio_bands(n_frames, D_hat, D_true, n_groups=20): """Running 10/50/90 % of d_hat/d_true in equal-count track-length groups. Quantiles are taken on D_hat/D_true, which exists for every track (even D_hat <= 0), then mapped to the size ratio D_true/D_hat. A rejected track therefore counts at the 'infinitely large' end instead of silently vanishing. """ order = np.argsort(n_frames, kind="stable") x, q10, q50, q90 = [], [], [], [] for g in np.array_split(order, n_groups): r90, r50, r10 = np.percentile(D_hat[g] / D_true[g], [90, 50, 10]) x.append(np.median(n_frames[g])) q10.append(1 / r90 if r90 > 0 else np.inf) q50.append(1 / r50 if r50 > 0 else np.inf) q90.append(1 / r10 if r10 > 0 else np.inf) return [np.array(v) for v in (x, q10, q50, q90)] # ============================================================================= # FIGURE # ============================================================================= def set_style(): plt.rcParams.update({ "font.family": "DejaVu Sans", "font.size": 11, "axes.titlesize": 13, "axes.titleweight": "bold", "axes.titlelocation": "left", "axes.titlepad": 10, "axes.labelsize": 12, "axes.labelcolor": INK, "xtick.labelsize": 11, "ytick.labelsize": 11, "xtick.color": INK_2, "ytick.color": INK_2, "text.color": INK, "axes.spines.top": False, "axes.spines.right": False, "axes.edgecolor": "#9aa1ab", "axes.linewidth": 0.9, "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.8, "axes.axisbelow": True, "legend.fontsize": 11, "legend.frameon": False, "figure.facecolor": "white", "axes.facecolor": "white", }) def make_figure(r): set_style() fig = plt.figure(figsize=(12, 20 / 3), dpi=150) # 1800 x 1000 px gs = fig.add_gridspec(2, 2, width_ratios=[1, 1.62], left=0.07, right=0.985, top=0.835, bottom=0.085, wspace=0.2, hspace=0.52) axA, axB, axC = fig.add_subplot(gs[:, 0]), fig.add_subplot(gs[0, 1]), fig.add_subplot(gs[1, 1]) # ---- A: example trajectories over the same time window ---------------- rng_ex = np.random.default_rng(SEED + 1) t_win = EXAMPLE_FRAMES / FPS extent, handles = 0.0, [] for d, c in zip(EXAMPLE_SIZES_NM, C_EXAMPLES): xy = simulate_tracks(rng_ex, d, EXAMPLE_FRAMES + 1)[0] / 1000.0 # um D_um = diameter_to_D(d) * 1e-6 # um^2/s rms = np.sqrt(4 * D_um * t_win) axA.add_patch(Circle((0, 0), rms, fill=False, lw=1.0, ls=(0, (3, 2.5)), ec=c, zorder=2)) axA.plot(xy[:, 0], xy[:, 1], color=c, lw=1.7, zorder=3, solid_joinstyle="round") axA.scatter(*xy[-1], s=44, color=c, ec="white", lw=1.2, zorder=5) extent = max(extent, np.abs(xy).max(), rms) handles.append(Line2D([], [], color=c, lw=2.2, label=f"{d:.0f} nm D = {D_um:.1f} µm²/s")) axA.scatter(0, 0, s=30, color=INK, zorder=6) lim = np.ceil(extent * 1.08) axA.set(xlim=(-lim, lim), ylim=(-lim, lim), xlabel="x (µm)", ylabel="y (µm)") axA.set_aspect("equal", adjustable="box") axA.xaxis.set_major_locator(MultipleLocator(5)) axA.yaxis.set_major_locator(MultipleLocator(5)) axA.set_anchor("N") axA.set_title(f"A {t_win:.0f}-s window: small EVs move farther") handles.append(Line2D([], [], color=INK_3, lw=1.0, ls=(0, (3, 2.5)), label=f"RMS displacement √(4Dt) at {t_win:.0f} s")) axA.legend(handles=handles, loc="upper center", bbox_to_anchor=(0.5, -0.12), ncol=1, handlelength=2.2, borderaxespad=0) # ---- B: size distributions on shared bins ----------------------------- s_true, s_det, s_rec = r["s_true"], r["s_det"], r["s_rec"] axB.hist(r["d_true"], BINS_NM, histtype="stepfilled", color=C_TRUE, alpha=0.13, lw=0, zorder=1) for d, c, z in ((r["d_true"], C_TRUE, 2), (r["d_det"], C_DET, 3), (r["d_rec"], C_REC, 4)): axB.hist(d, BINS_NM, histtype="step", color=c, lw=2, zorder=z) for s, c in ((s_true, C_TRUE), (s_det, C_DET), (s_rec, C_REC)): axB.axvline(s["mode"], color=c, lw=1.4, ls=(0, (4, 3)), zorder=5) axB.set(xlim=(0, 400), xlabel="Diameter (nm)", ylabel="Particles per 10-nm bin") axB.set_ylim(0, axB.get_ylim()[1] * 1.04) axB.set_title("B Size distributions: true, detected, NTA-recovered") hB = [Patch(fc=(0.165, 0.47, 0.84, 0.13), ec=C_TRUE, lw=2, label=f"True population (n = {s_true['count']:,}) · mode {s_true['mode']:.0f} nm"), Line2D([], [], color=C_DET, lw=2, label=f"Detected subset (n = {s_det['count']:,}) · mode {s_det['mode']:.0f} nm"), Line2D([], [], color=C_REC, lw=2, label=f"NTA-recovered (n = {s_rec['count']:,}) · mode {s_rec['mode']:.0f} nm"), Line2D([], [], color=INK_3, lw=1.4, ls=(0, (4, 3)), label="Mode of each (dashed)")] axB.legend(handles=hB, loc="upper right", borderaxespad=0.2) n_off = int(np.sum(r["d_rec"] > BINS_NM[-1])) if n_off: axB.annotate(f"+{n_off} recovered tracks > {BINS_NM[-1]:.0f} nm (off-scale)", xy=(0.995, 0.43), xycoords="axes fraction", ha="right", va="top", color=INK_2, fontsize=11) # ---- C: per-track error vs track length -------------------------------- rng_plot = np.random.default_rng(SEED + 2) ratio, f = r["ratio"], r["n_frames"] ok = np.isfinite(ratio) y_lo, y_hi = 0.25, 10.0 axC.axvspan(0, MIN_TRACK_FRAMES - 0.5, color="#9aa1ab", alpha=0.14, lw=0, zorder=0) axC.text((MIN_TRACK_FRAMES - 0.5) / 2, np.sqrt(y_lo * y_hi), f"discarded (< {MIN_TRACK_FRAMES} frames)", ha="center", va="center", rotation=90, fontsize=11, color=INK_2) jit = f + rng_plot.uniform(-0.35, 0.35, f.size) inside = ok & (ratio >= y_lo) & (ratio <= y_hi) axC.scatter(jit[inside], ratio[inside], s=7, color=C_REC, alpha=0.3, lw=0, zorder=2) over, under = ok & (ratio > y_hi), ok & (ratio < y_lo) # off-scale: edge triangles axC.scatter(jit[over], np.full(over.sum(), y_hi / 1.06), s=16, marker="^", color=C_REC, alpha=0.7, lw=0, zorder=2, clip_on=False) axC.scatter(jit[under], np.full(under.sum(), y_lo * 1.06), s=16, marker="v", color=C_REC, alpha=0.7, lw=0, zorder=2, clip_on=False) x, q10, q50, q90 = r["bands"] axC.fill_between(x, np.clip(q10, y_lo, y_hi), np.clip(q90, y_lo, y_hi), color=C_REC, alpha=0.22, lw=0, zorder=3) axC.plot(x, np.clip(q50, y_lo, y_hi), color=INK, lw=2.2, zorder=4) axC.axhline(1.0, color=INK_2, lw=1.1, ls=(0, (1.5, 2)), zorder=1) axC.set_yscale("log") axC.set_ylim(y_lo, y_hi) axC.yaxis.set_major_locator(FixedLocator([0.25, 0.5, 1, 2, 4, 8])) axC.yaxis.set_minor_locator(NullLocator()) axC.yaxis.set_major_formatter(FuncFormatter(lambda v, _: f"{v:g}×")) axC.set_xlim(0, np.percentile(f, 99.5) + 5) axC.set(xlabel=f"Track length (frames at {FPS:.0f} fps)", ylabel="Recovered / true diameter") axC.set_title("C Per-track sizing error: short tracks scatter widely") hC = [Line2D([], [], ls="", marker="o", ms=5, color=C_REC, alpha=0.6, label=f"Track (n = {ok.sum():,})"), Line2D([], [], color=INK, lw=2.2, label="Running median"), Patch(fc=C_REC, alpha=0.22, lw=0, label="10 to 90% of tracks"), Line2D([], [], color=INK_2, lw=1.1, ls=(0, (1.5, 2)), label="Exact sizing (1×)")] axC.legend(handles=hC, loc="upper right", ncol=2, borderaxespad=0.2, columnspacing=1.4) axC.annotate(f"▲▼ off-scale · {r['n_nonpos']} tracks with D ≤ 0 cannot be sized", xy=(0.995, 0.03), xycoords="axes fraction", ha="right", va="bottom", color=INK_2, fontsize=11) # ---- title + subtitle -------------------------------------------------- fig.text(0.07, 0.965, "How faithfully does NTA recover an EV size distribution?", fontsize=17, weight="bold", ha="left", va="top") (m1, g1, w1), (m2, g2, w2) = COMPONENTS fig.text(0.07, 0.912, f"{N_PARTICLES:,} EVs (lognormal: {w1:.0%} {m1:.0f} nm + {w2:.0%} {m2:.0f} nm) · " f"50% detected at {DET_D50_NM:.0f} nm · {FPS:.0f} fps · σloc {SIGMA_LOC_NM:.0f} nm · " f"tracks ≥ {MIN_TRACK_FRAMES} frames · MSD lags {MSD_LAGS[0]} to {MSD_LAGS[-1]} fit", fontsize=11.5, color=INK_2, ha="left", va="top") fig.savefig(OUT_PNG, dpi=150, facecolor="white") plt.close(fig) # ============================================================================= # MAIN # ============================================================================= def main(): t0 = time.perf_counter() rng = np.random.default_rng(SEED) # 1. Ground truth d_true = sample_population(rng, N_PARTICLES) # 2. Optical detection detected = rng.random(d_true.size) < detection_probability(d_true) d_det = d_true[detected] # 3. Track lengths, then the minimum-track-length filter n_frames = rng.geometric(1.0 / MEAN_TRACK_FRAMES, d_det.size) kept = n_frames >= MIN_TRACK_FRAMES d_trk, f_trk = d_det[kept], n_frames[kept] # 4. Simulate tracks, estimate D per track xy = simulate_tracks(rng, d_trk, f_trk) D_true = diameter_to_D(d_trk) D_msd, D_cve = estimate_D_msd(xy), estimate_D_cve(xy) # 5. Back to diameters; reject D <= 0 and sizes outside the instrument range lo, hi = SIZE_RANGE_NM d_msd, d_cve = D_to_diameter(D_msd), D_to_diameter(D_cve) ok_msd = np.isfinite(d_msd) & (d_msd >= lo) & (d_msd <= hi) ok_cve = np.isfinite(d_cve) & (d_cve >= lo) & (d_cve <= hi) d_rec, d_rec_cve = d_msd[ok_msd], d_cve[ok_cve] # ---- summary table ------------------------------------------------------ stats = {"True population": summarize(d_true), "Detected (true sizes)": summarize(d_det), "NTA-recovered (MSD fit)": summarize(d_rec), " alt: recovered (CVE)": summarize(d_rec_cve)} head = (f"{'':25s}{'count':>7s}{'mean':>8s}{'mode':>8s}{'D10':>8s}{'D50':>8s}" f"{'D90':>8s}{'<' + format(SMALL_CUTOFF_NM, '.0f') + ' nm':>9s}") print("\nNTA recovery of an EV size distribution (sizes in nm)") print(head) print("-" * len(head)) for name, s in stats.items(): print(f"{name:25s}{s['count']:7d}{s['mean']:8.1f}{s['mode']:8.1f}{s['d10']:8.1f}" f"{s['d50']:8.1f}{s['d90']:8.1f}{100 * s['frac_small']:8.1f}%") s_t, s_d, s_r, s_c = stats.values() shift = lambda s: 100 * (s["d50"] / s_t["d50"] - 1) small = d_true < SMALL_CUTOFF_NM n_nonpos_msd, n_nonpos_cve = int(np.sum(D_msd <= 0)), int(np.sum(D_cve <= 0)) n_out_msd = int(np.sum(np.isfinite(d_msd) & ~ok_msd)) n_out_cve = int(np.sum(np.isfinite(d_cve) & ~ok_cve)) print(f"\nD50 shift vs true: detection only {shift(s_d):+.1f}% | " f"detection + tracking (MSD fit) {shift(s_r):+.1f}% | with CVE {shift(s_c):+.1f}%") print(f"Never detected: {100 * (1 - detected.mean()):.1f}% of all true particles " f"({100 * (1 - detected[small].mean()):.1f}% of those < {SMALL_CUTOFF_NM:.0f} nm)") print(f"Funnel: {d_true.size} true -> {d_det.size} detected -> {kept.sum()} tracks >= " f"{MIN_TRACK_FRAMES} frames ({(~kept).sum()} discarded) -> {ok_msd.sum()} sized") print(f"Non-physical D <= 0: {n_nonpos_msd} tracks (MSD fit), {n_nonpos_cve} (CVE); " f"outside {lo:.0f}-{hi:.0f} nm: {n_out_msd} (MSD fit), {n_out_cve} (CVE)") if n_nonpos_msd: print(f" longest track with D <= 0 (MSD fit): {f_trk[D_msd <= 0].max()} frames") sub = ok_msd & (d_msd < SMALL_CUTOFF_NM) print(f"Recovered < {SMALL_CUTOFF_NM:.0f} nm: {sub.sum()} tracks, of which " f"{100 * np.mean(d_trk[sub] >= SMALL_CUTOFF_NM):.0f}% are truly >= {SMALL_CUTOFF_NM:.0f} nm " f"(sized small because D was overestimated)") ratio = d_msd / d_trk # NaN where D <= 0 bands = size_ratio_bands(f_trk, D_msd, D_true) for lab, sel in ((f"{MIN_TRACK_FRAMES}-14 frames", f_trk < 15), (">= 50 frames", f_trk >= 50)): rr = D_msd[sel] / D_true[sel] r90, r50, r10 = np.percentile(rr, [90, 50, 10]) inv = lambda v: 1 / v if v > 0 else np.inf print(f"Per-track d_rec/d_true, tracks {lab:>12s} (n={sel.sum():4d}): median {inv(r50):.2f}, " f"10-90% band {inv(r90):.2f}-{inv(r10):.2f}") make_figure(dict(d_true=d_true, d_det=d_det, d_rec=d_rec, s_true=s_t, s_det=s_d, s_rec=s_r, ratio=ratio, n_frames=f_trk, bands=bands, n_nonpos=n_nonpos_msd)) print(f"\nFigure: {OUT_PNG}") print(f"Runtime: {time.perf_counter() - t0:.2f} s (seed {SEED})") if __name__ == "__main__": main()