Multicellular simulations with shape and volume constraints using optimal transport

TL;DR

A new framework based on optimal transport theory to model particle systems with arbitrary dynamical shapes and deformability properties is introduced, builds upon the pioneering work of Brenier on incompressible fluids and its recent applications to materials science.

Abstract

Many living and physical systems such as cell aggregates, tissues or bacterial colonies behave as unconventional systems of particles that are strongly constrained by volume exclusion and shape interactions. Understanding how these constraints lead to macroscopic self-organized structures is a fundamental question in e.g. developmental biology. To this end, various types of computational models have been developed. Here, we introduce a new framework based on optimal transport theory to model particle systems with arbitrary dynamical shapes and deformability properties. Our method builds upon the pioneering work of Brenier on incompressible fluids and its recent applications to materials science. It lets us specify the shapes and volumes of individual cells and supports a wide range of interaction mechanisms, while automatically taking care of the volume exclusion constraint at an affordable numerical cost. We showcase the versatility of this approach by reproducing several classical systems in computational biology. Our Python code is freely available at \url{https://iceshot.readthedocs.io/}.

Figures (11)

• Figure 1: Graphical abstract. (a) Laguerre tesselations generalize Voronoi diagrams and level-set approaches with volume, shape and deformation constraints encoded in a cost function $c$ which can be customized and dynamic. (b) Any active point-particle model can be implemented with additional arbitrary softness and deformation properties. (c) The framework is independent of the dimension and is implemented in 3D. (d) Laguerre tessellations are computed as the solution of a semi-discrete optimal transport problem on a discrete grid, resulting in a map $T$ which assigns each voxel to a cell. (e) In 3D, a meshing of each cell boundary is computed in order to implement surface tension effects and cell sorting mechanisms.
• Figure 2: Static and dynamic shapes in 2D. (a) Three Laguerre tessellations: (Left) Voronoi diagram obtained with the $L^2$ cost and random volumes $v_i$ sampled uniformly with a ratio 1/5. (Center) A bubbly tessellation similar to ishimoto_bubbly_2014 obtained with 66 particles with random volumes sampled uniformly with a ratio 1/20 and power costs \ref{['eq:powercost']} with exponents distributed uniformly between 0.5 and 4. Lighter colors indicate lower values of this exponent and correspond to softer shapes. (Right) Voronoi tessellation with random fluctuations similar to miyazaki_mechanism_2023 obtained with 42 identical particles and a randomly perturbed $L^2$ cost. (b) A swarm of rod-shaped particles obtained with the cost \ref{['eq:spherocylindercost']} showing the emergence of long-range alignment due to deformations and non-overlapping interactions. See also SM Video \ref{['video:rodshape']}. (c) Chemotaxis motion induced by shape deformations with two choices of biased costs \ref{['eq:costpotentialchemotaxis']}-\ref{['eq:bias']} See also SM Videos \ref{['video:chemo_long']},\ref{['video:chemo_fan']}
• Figure 3: 3D benchmarking examples. (a) Falling soft-spheres in a hourglass domain. See also Supplementary Video \ref{['video:fallingstuff']}. (b) Exponential growth of a 3D aggregate via successive cell division and growth phases, zooming out from $N=1$ to $N=50,000$ cells. See also SM Video \ref{['video:tissue_growth_3D']}. (c) Final configuration of a deformation-driven run-and-tumble motion for $N=10,100,1000,10000$ deformable ellipsoids with a space discretization grid of size $M=512$. See also SM Video \ref{['video:benchmark']}.
• Figure 4: Active Brownian Particles with deformations (a) Phase diagram in the $(\alpha,c_0)$ plane. (b) Snapshots of the final configuration for different values of $(\alpha,c_0)$ indicated by the labels (A),(B),(C),(D) on the phase diagram. See also SM Videos \ref{['video:abp_a05']},\ref{['video:abp_a2']},\ref{['video:abp_a10']}. (c) (left) Coefficient of determination of the linear fit in the $(\alpha,c_0)$ plane. The center of the color map is set at $r^2=0.9$. (center) Average shape index $\langle\sigma_i\rangle$ in the $(\alpha,c_0)$ plane. The center of the colormap is set at $\sigma=3.81/(2\sqrt{\pi})$. (right) MSD (plain line) and linear fit (dashed line) over time for $c_0=0.5$ and $\alpha=0.2, 1.7, 7.7$. The black line has slope $\sigma_\text{thr}$. Insert: zoom on the MSD curve $\alpha=7.7$ showing saturation when the system has reached a stable hexagonal pattern.
• Figure 5: Sorting patterns in 3D. (a)-(b) Sorting patterns in a homogeneous mixture of two cell types with resp. $\eta=3>1$ and $\eta=0.3<1$, obtained from an initial mixed aggregate of $N=120$ cells, under different conditions on the surface tensions parameters. The phase domain compares relative compactness $\overline{k}$ with relative softness $\overline{\gamma}$. The line $\{\overline{k}\overline{\gamma}=1\}$ defines the boundary between the regions $\{\eta_{oo} \lessgtr \eta_{bb}\}$ and the line $\{\overline{k}\overline{\gamma}\overline{\eta}=1\}$ defines the boundary between the regions $\{\eta_{ob}\lessgtr\eta_{bb}\}$. Two disconnected regions with the same background color are equivalent with blue and orange cells inverted depending on $\overline{\gamma}\gtrless1$. (c) Engulfment of an aggregate of orange cells in an aggregate of blue cells starting from a configuration where the two cell types are separated (region A: $\overline{\eta}=3$, $\overline{\gamma}=2$, $\overline{k}\overline{\gamma}\overline{\eta}=0.8$), see also SM Video \ref{['video:st_engulfment']}.

Theorems & Definitions (2)

• Theorem 1
• Theorem 2