Multicellular simulations with shape and volume constraints using optimal transport

Technical Summary: Multicellular Simulations with Shape and Volume Constraints Using Optimal Transport

Problem Statement
Living systems, ranging from intracellular components to macroscopic tissues and bacterial colonies, exhibit complex dynamics driven by volume exclusion, shape interactions, and deformability. Understanding how these constraints lead to macroscopic self-organization is a fundamental challenge in developmental biology and biophysics. Existing computational models face a trade-off between physical accuracy, numerical cost, and mathematical tractability. Common approaches include:

• Point-particle systems: Mathematically simple but often lack realistic shape representation.
• Phase-field and level-set models: Capable of representing arbitrary shapes but can be computationally expensive and challenging for contact interactions.
• Cellular automata and lattice models: Discretize space into voxels, potentially limiting geometric flexibility.
• Voronoi and vertex models: Efficient for tissue-like structures but often rely on soft constraints for volume or lack control over specific shape deformations.

A significant gap exists in modeling systems with strict volume constraints, arbitrary dynamical shapes, and high particle counts (tens to tens of thousands) within a unified framework that supports both individual deformability and collective interactions.

Methodology
The authors introduce a mesoscale simulation framework grounded in optimal transport theory, specifically utilizing Laguerre tessellations.

1. Static Model (Laguerre Tessellation):
   The space occupied by each particle $i$ is modeled as a Laguerre cell $L_i$, defined as the set of points $x$ satisfying:
   $$L_i = \{x \in \Omega, c(x, x_i) - w_i \leq c(x, x_j) - w_j \text{ for all } j\}$$
   Here, $x_i$ is the particle position, $v_i$ is its target volume, and $c(x, x_i)$ is a cost function encoding the particle's shape and deformability properties. The weights $w_i$ (Kantorovich potentials) are uniquely determined to satisfy the strict volume constraint $|L_i| = v_i$.

   • This formulation generalizes Voronoi diagrams (where $w_i=0$ and $c$ is Euclidean distance) to allow for curved boundaries and non-spherical shapes by choosing arbitrary cost functions (e.g., power costs, anisotropic costs).
   • Empty space is modeled as an additional particle with zero cost.

2. Dynamical Model:
   The evolution of particle positions $x_i$ follows a first-order gradient descent equation:
   $$\dot{x}_i = -\tau_i \nabla_{x_i} T_c + F_i^{\text{point}} + F_i^{\text{surf}} + \dot{\xi}_i$$

   • Incompressibility Force: The term $-\nabla_{x_i} T_c$ acts as a repulsive force derived from the total transport cost $T_c$. It ensures particles do not overlap and naturally enforces volume constraints. Unlike pairwise forces, this force depends on the global geometry of all contact surfaces.
   • External/Interaction Forces ($F^{\text{point}}$): Includes self-propulsion, gravity, or attraction, similar to point-particle models.
   • Surface Forces ($F^{\text{surf}}$): Models pressure and surface tension effects along interfaces, allowing for phenomena like cell sorting and adhesion.
   • Shape Dynamics: Shapes can evolve via parametrized costs (e.g., ellipsoids defined by covariance matrices) or level-set potential functions that respond to external fields (e.g., chemotaxis).

3. Numerical Implementation:

   • The problem is solved as a semi-discrete optimal transport problem on a discrete grid (voxels).
   • The authors employ a divide-and-conquer multiscale strategy and a quasi-Newton solver (L-BFGS-B) to compute the Kantorovich potentials $w_i$ efficiently.
   • The implementation is massively parallelized on GPUs using the KEOPS and GEOMLOSS libraries, enabling simulations of up to 50,000 particles in 3D with state-of-the-art performance.

Key Contributions

• Unified Framework: The method unifies concepts from Voronoi tessellations, level-set methods, and vertex models under the umbrella of optimal transport. It supports arbitrary shapes, strict volume constraints, and complex interaction mechanisms within a single framework.
• Strict Volume Enforcement: Unlike many existing models that use soft constraints or relaxation forces, this approach mathematically guarantees volume constraints at all times through the optimization of Kantorovich potentials.
• Scalability: The GPU-accelerated solver allows for the simulation of large-scale systems (up to 50,000 particles) in 3D, which is computationally prohibitive for many traditional finite-element or phase-field methods.
• Versatility in Modeling: The framework successfully models diverse phenomena including:
  • Fluid-solid phase transitions in active Brownian particles driven by deformability.
  • Emergent orientational order in swarming rod-shaped particles solely due to volume exclusion.
  • Chemotactic migration driven by shape deformations.
  • Cell sorting and engulfment phenomena in 3D heterogeneous aggregates, reproducing results consistent with the Differential Interfacial Tension Hypothesis (DITH).

Results
The authors validate the framework through several benchmarks:

• Growth and Division: A 3D aggregate grows from a single cell to 50,000 cells via division, filling a domain in approximately one day of computation.
• Phase Transitions: In a system of active Brownian particles, the deformability parameter $\alpha$ drives a transition from a fluid phase (diffusive behavior) to a solid phase (jammed, hexagonal packing). The model captures the vanishing effective diffusion and the saturation of mean-square displacement in the solid phase.
• Rod-Shape Alignment: Simulations of self-propelled rod-shaped particles demonstrate the spontaneous emergence of long-range orientational order due to collision-induced deformations, consistent with polymer physics theories.
• 3D Cell Sorting: The model reproduces various sorting patterns (engulfment, separation, checkerboard) in 3D mixtures of two cell types with different surface tensions and softness, matching theoretical predictions and experimental observations.
• Crowd Dynamics: The model simulates crowd evacuation, showing that deformable particles (fluid-like) evacuate safely, while hard-sphere particles (solid-like) form stable arches and jam at exits.

Significance and Claims
The paper claims to introduce a novel approach for modeling multicellular systems that bridges the gap between detailed physical descriptions and computational efficiency.

• Validation: The primary contribution is the validation of this new framework against classical computational models, demonstrating that it can reproduce established results (e.g., cell sorting, phase transitions) while offering superior handling of volume constraints and 3D complexity.
• Biological Relevance: The authors highlight the method's potential for in silico modeling of complex developmental processes, such as blastocyst development, where both individual cell deformations and dense aggregate behaviors are critical.
• Theoretical Extension: The work suggests a new optimal transport perspective for existing methods (point-particle, level-set, vertex) and provides a pathway toward deriving coarse-grained continuum models (mean-field approximations) for systems with volume interactions, a largely open question in mathematical biology.
• Accessibility: The authors provide a freely available Python implementation, making the method accessible for further research in computational physics and biology.

The authors remain modest, noting that while the framework is new, the immediate focus was on validation against existing models. They identify future directions, such as extending active vertex theory to Laguerre tessellations, refining the derivation of continuum limits, and improving the handling of complex dynamical costs and surface interactions for real-time simulation.