A transport-only null model for apparent heterogeneity in diffusively dosed organoid arrays

Technical Summary: A Transport-Only Null Model for Apparent Heterogeneity in Organoid Arrays

Problem Statement
Patient-derived liver-cancer organoids are increasingly used to study tumor biology, yet nominally identical organoids subjected to the same dosing protocol often exhibit substantial phenotypic variation. While this variation is frequently attributed to biological heterogeneity or measurement error, a third mechanism is often overlooked: unequal intracellular exposure driven by spatial transport. In a shared diffusible pool, an organoid's position relative to the dosing source, chamber boundaries, and neighboring organoids dictates the material reaching and entering each compartment. Consequently, spatially heterogeneous delivery can generate heterogeneous phenotypes even when all biological parameters are identical. The paper addresses the lack of a usable "transport-only null" to quantify how much between-organoid variation can be generated solely by geometry and delivery protocols.

Methodology
The authors construct a transport-to-phenotype null model coupling bulk diffusion, surface kinetics, and intracellular state dynamics.

1. Bulk Transport and Surface Kinetics: The model treats the extracellular domain as a bounded chamber containing $N$ small, well-separated compartments (organoids). It couples bulk diffusion with encounter-dependent adsorption, reversible surface residence, and productive internalization.

   • Renewal Theory: The model employs renewal equations to describe repeated adsorption, residence, desorption, and re-adsorption cycles. It distinguishes between "resetting" (where a desorbed molecule forgets its location and re-enters a well-mixed fluid) and "continued search" (where spatial memory is retained).
   • Boundary Conditions: The surface exchange is modeled as a frequency-dependent reactive filter rather than a static sink, incorporating a boundary local time formalism to handle non-Markovian encounter dynamics.

2. Dimensional Reduction via Matched Asymptotics:

   • For two-dimensional domains, matched asymptotics reduce the perforated-domain problem to an $N \times N$ Green-function system evaluated at the organoid centers.
   • The reduction expresses the productive flux in terms of monopole strengths and shape-dependent capacitances (in 3D) or logarithmic interactions (in 2D).
   • The error of this spatial reduction is proven to be $O(\varepsilon^2)$ for circular targets, where $\varepsilon$ is the ratio of organoid radius to chamber size.

3. Intracellular State Dynamics: The productive internalization flux drives a five-dimensional cell-state model (progenitor, differentiated, mature, apoptotic/stressed, and disease score). The model assumes identical biological parameters across all organoids; thus, any variation in the output phenotype arises strictly from the geometry-dependent productive fluxes.

4. Numerical Verification:

   • Spatial: The reduced model is verified against exact radial benchmarks and full finite-element solutions (P2 elements) of the perforated-domain problem for $N=2, 5, 10$ targets.
   • Temporal: Time-domain reconstruction is validated against numerical Laplace inversion. The authors compare a two-moment Gamma surrogate against exact inversion, noting that while the surrogate is highly accurate for Markovian (exponential) residence laws, it can fail for non-Markovian (heavy-tailed) laws.

Key Contributions

1. Transport-Only Null Formulation: The paper formulates a direct recipe to map measured organoid centers and dosing geometry to productive uptake and phenotype variation under identical biology. This provides a layout-specific reference against which observed heterogeneity can be tested.
2. Quantification of Spatial Spread: Across 2000 random ten-organoid arrays with localized dosing, the model predicts a median transport-only maturation coefficient of variation (CV) of 0.623 and a disease-score CV of 0.835. A one-factor design audit shows this median can range from 0.27 to 0.86 depending on geometry.
3. Mechanistic Clarifications:
   • Conservative vs. Lossy Search: In a conservative chamber (reflecting walls, no bulk loss), the total eventual productive uptake probability is 1. Desorption delays delivery and redistributes uptake but does not reduce total eventual uptake unless a competing loss channel exists.
   • Residence Law Shape: Residence laws with equal means but different shapes (e.g., Gamma vs. Exponential vs. Long-tailed) produce distinct transient phenotypes. The transient depends on the shape of the law, not just its mean.
   • Competition Normalization: Raw uptake conflates competition with source distance. A competition-normalized ratio is derived that cancels source forcing, isolating target-target interaction.
   • Identifiability of Remixing: Partial remixing (interpolation between resetting and continued search) is identifiable from endpoint uptake fractions only when retention is heterogeneous; time-resolved uptake is generically informative even under homogeneous retention.
4. Design Implications:
   • Dosing Geometry: Distributed dosing (e.g., four ports) reduces the baseline spread approximately fivefold compared to a single localized port after matching array-mean exposure.
   • Batch Size Limits: Increasing batch size averages only independent within-batch process variation. Shared line and batch effects yield a generally nondegenerate large-batch acceptance probability, meaning larger batches do not necessarily converge to a deterministic outcome if shared effects are present.

Results

• Baseline Variability: For a localized port, the median maturation CV is 0.623. This drops to 0.124 for a spatially uniform source (after matching mean exposure), demonstrating a roughly fivefold reduction in spread.
• Sensitivity: The variability is highly sensitive to target count, target size, chamber aspect ratio, and dosing location (e.g., corner vs. center). It is largely insensitive to bulk loss rates in the conservative limit.
• Convergence: The reduced model achieves $O(\varepsilon^2)$ accuracy for productive flux. At separations of about four target radii, the error in the phenotype CV is below $10^{-3}$, though individual flux errors remain higher ($0.5\% - 2\%$).
• Batch Acceptance: The probability of batch acceptance does not tend to 0 or 1 as batch size increases if shared effects (line or batch) are present. Instead, it converges to an interior limit determined by the distribution of shared effects.

Significance and Claims
The paper claims to provide a geometry-specific null model rather than a calibrated digital twin. Its primary significance is establishing a baseline: before observed heterogeneity is attributed to biological variation, the measured geometry and delivery protocol must be quantified to see how much variation they can explain.

• A measurement exceeding the model's prediction strictly identifies biological or kinetic variation that transport alone cannot account for.
• Failure to exceed the prediction does not prove biological homogeneity; it merely indicates the data are indistinguishable from transport-induced noise.
• The framework offers specific, testable design recommendations: distributed dosing significantly reduces apparent heterogeneity, and increasing batch size has diminishing returns if shared process effects are not addressed.

The authors explicitly state that the mathematical ingredients (renewal theory, matched asymptotics, Dirichlet-to-Neumann reductions) are inherited from existing literature. The novelty lies in their composition into a verified transport-to-phenotype model for liver-cancer organoid arrays and the derivation of the specific consequences regarding forcing equalization, retention interpretation, and batch acceptance limits.