# An optimal transport model for dynamical shapes, collective motion and
cellular aggregates

Antoine Diez Institute for the Advanced Study of Human Biology (ASHBi) Kyoto
University Institute for Advanced Study, Kyoto University, Yoshida-Konoe-cho,
Sakyo-ku, Kyoto 606-8501, Japan.
Jean Feydy HeKA team, Inria Paris, Inserm, Université Paris-Cité, PariSanté
Campus, 2-10 rue d’Oradour-sur-Glane, 75015 Paris, France

###### Abstract

Many biological 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: phase fields, cellular automata, vertex models, level-
set, finite element simulations, etc. We introduce a new framework based on
optimal transport theory to model particle systems with arbitrary dynamical
shapes and deformability. 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 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:
www.github.com/antoinediez/ICeShOT

Figure 1: Graphical abstract. (a) We model cells on a domain
$\Omega\subset\mathbb{R}^{d}$ using positive volumes $v_{i}>0$ and cost
functions $c_{i}:y\mapsto c_{i}(y)\in\mathbb{R}$ that penalize the assignment
of a location $y$ to the $i$-th cell. To compute the shapes of a collection of
$N$ cells $(v_{1},c_{1}),\dots,(v_{N},c_{N})$, we rely on a decomposition of
the domain $\Omega$ into pixels $y_{1},\dots,y_{M}$. Then, assigning the
points $y_{j}$ to their preferred cells corresponds to a discrete optimal
transport problem that can be solved efficiently in dimension 2 (b) and 3 (c).
Depending on the cost function $c_{i}$, a cell may either be hard (c, blue) or
soft (c, orange), see also SM Video 24. To model dynamical systems in silico,
we alternate between computing cell shapes with optimal transport and updating
the cell parameters $(v_{i},c_{i})$ according to various interaction models.
(d-e) Our method is versatile: it supports arbitrary cell shapes and
interactions, while automatically taking care of the volume exclusion
constraint. Using plain Python code, we show that collective behaviours emerge
from interactions between incompressible cell shapes.

##  1 Introduction

From the smallest intracellular scale [51] to the macroscopic population scale
[83], shape, deformations and congestion effects play a central role in the
dynamics of biological systems. At the scale of micrometers, cells, bacteria
and other microorganisms appear in a wide variety of morphologies. They form
dense, dynamical, complex multicellular structures such as epithelial tissue
or swarming colonies with a macroscopic behavior [95]. Among other examples:
the peculiar shape of some bacteria conditions their collective motion [34,
130]; the membrane deformations of mesenchymal cells induce their migration
[123, 121]; the complex viscoelastic properties of cell-cell contact
interactions [76, 134, 20] give rise to complex sorting patterns [3, 84] and
tissue structures [32, 56, 91, 92, 132, 71]. Crucially, in such crowded
systems, the shapes of cells and their varied degrees of deformability
strongly constrain global movements and morphological changes [6].
Understanding and measuring these effects is a challenging but important
experimental problem, especially in developmental biology [63, 62].

Over the past decades, in complement to in vitro and in vivo experimental
studies, these questions have been studied using mathematical modeling. In
silico experiments let us assess the impact of conjectured biophysical laws in
minimal experimental settings. Note, however, that there is no unique
mathematical way of representing dynamical shapes and multicellular systems.
This has motivated the development of various classes of models. We briefly
review some of them below for later comparison with our novel framework. We
also refer to [103, 21, 135] and to the introductions of [93, 69] for an in-
depth discussion of the different approaches.

  * •

Enriched point-particle systems. Extending the core idea of statistical
physics, the pioneering work of Vicsek [136] has shown that simple but
potentially high-dimensional dynamical systems can accurately describe
important collective motion phenomena in biology. In this class of models, one
uses a group of active particles to represent migrating cells, microorganism
colonies or larger animal societies. Each particle is described by a set of
xyz coordinates coupled with additional basic information variables such as
polarity vectors [136, 30], 3D orientation matrices [39] or finite-dimensional
shape parameters (size, aspect ratio, etc.) [77, 53]. The formal simplicity of
these models can be a drawback, since particles are not detailed enough to
represent realistic morphologies. On the other hand, this also makes particle
models appealing for both coarse-grained mathematical analysis [28, 26] and
large-scale numerical simulations [40, 109].

  * •

Phase fields, level-sets and boundary flows. Representing arbitrary shapes
requires an infinite-dimensional description. Several techniques can be
borrowed from the study of free-boundary and interface problems. For instance,
instead of using a single (enriched) point, a cell can be described as a
smoothed indicator function of the set that represents its full shape. Using a
dynamical Allen-Cahn type motion, this phase-field approach has been applied
successfully many times in biology [99, 74, 120, 97]. Strongly connected
approaches consider the evolution of the boundary of a shape, defined as a
curve [102, 113, 58, 101] or as the level set of some potential function [105,
104, 122, 140, 93, 125].

  * •

Cellular automata and lattice models. With a completely different approach,
Glazier and Graner [49, 52] introduced the Cellular Potts framework as an
extension of lattice models from statistical physics. In this setting, the
spatial domain is discretized into a uniform grid lattice with each “pixel”
representing either one biological cell or a portion of it. The pixel
allocation and flipping dynamics are defined using custom energy minimization
rules and implemented via a stochastic Metropolis-type algorithm.

  * •

Voronoi tessellation. Historically, cellular arrangement in a tissue were
first described as an off-lattice tiling of the space by polygonal shapes [59,
128]. In the classical Voronoi tessellation, each polygonal cell is generated
from a seed point (or centroid) that can be interpreted as a cell nucleus by
considering all the spatial points which are closer to this seed point than to
any other seed point. This optimization problem is relatively simple to solve
and produces sufficiently realistic configurations. Recent experimental
observations [70, 91] have motivated the development of refined optimization
criteria and detailed interaction mechanisms to model contact forces or
locomotion [10, 9]. Our approach generalizes this line of work to flexible,
non-polygonal shapes.

  * •

Vertex models, meshes and finite elements. Unlike Voronoi tessellations for
which the centroids are the main degrees of freedom, vertex models define a
tiling of the space via a set of vertices connected by straight edges in 2D or
plane surfaces in 3D [60]. The motion of each vertex is driven by a force
balance equation that often comes from an energy minimization hypothesis. A
detailed mathematical formulation and successful applications to e.g. the
epithelial tissue are presented in [2, 119] and references therein. We also
mention several recent improvements of the method that build connections with
point-particle and Voronoi tessellation models [5, 8] or that consider general
edge shapes [64]. Another extension of the vertex model applied to
multicellular systems divides each cell into a finite number of subcellular
elements, thus creating an effective meshing of the cell membrane, the cell
interior and/or the cell environment. This leads to powerful computational
approaches originally developed for fluid and fluid-structure interactions,
such as the immersed boundary method [107, 66, 112] and the finite element
method [31, 133, 141].

Importantly for computational biology, all the methods mentioned above are now
available as (mostly) open source software. Among others, we may cite SiSyPHE
and PhysiCell for particle and spheroid models [40, 48], MorphoSim for phase-
fields [73], Morpheus and CompuCell3D for cellular automata models [126, 129],
Tissue Forge for vertex models [119], Chaste [33] for a wide range of cell-
based modelling frameworks, Omniverse™, Interacting Active Surfaces and Goo
for finite element and physics based models [100, 133, 86]. We also mention
cellPACK which solves 3D packing problems with objects of various sizes and
shapes for applications to intracellular crowding [67].

Our approach. Be it for cells in multicellular aggregates, microorganisms or
pedestrians, the first measurable quantity is the volume. As a starting point,
we thus build our method around an independent volume constraint for each
cell. This does not mean that volumes cannot evolve in time, but that initial
volumes are given as an input of the model and are treated as first-class
variables. This contrasts with most of the approaches described above (phase-
fields, cellular automata and vertex models), which preserve cell volumes
using soft constraints and relaxation forces towards a preferred size. Voronoi
tessellation methods do not usually consider a volume constraint, or only
provide little control on this quantity [10]. As a first description, our
method can be seen as a generalized Voronoi tessellation method which includes
strict volume constraints. To achieve this goal, we rely on the notion of
Laguerre tessellation which has recently appeared in a different context, for
optimal transport simulations of incompressible fluid flows.

Classicaly, the motion of an incompressible fluid is described through
macroscopic quantities such as the local density or local speed, which satisfy
the Euler or Navier-Stokes equations. The pioneering work of Brenier [16, 18,
19] has introduced a Lagrangian description of this motion, i.e. a set of
equations which consider the individual motion of punctual fluid particles.
Importantly, Brenier showed that this motion is intimately linked to the
theory of optimal transport. Originally developed for operations research and
economics by Monge and Kantorovich, the theory of optimal transport describes
the most effective way of allocating resources (canonically, piles of sand or
flour) from one location to another while minimizing a transportation cost
that is a function of the distance. We refer to [42, 108, 114, 137, 138] and
SM Appendix 3 for classical references. In the modern implementation of
Brenier’s ideas [46, 80], the notion of Laguerre tessellation appears as the
natural spatial discretization procedure that preserves the incompressibility
(i.e. the volume) constraint. Related ideas are also used in [83, 75] to model
congestion and crowd motion, where volume exclusion also plays a crucial role.

A fundamental insight from the literature is that since the optimal transport
problem generalizes sorting to spaces of dimension $d>1$, it can be solved
efficiently using a divide-and-conquer “multiscale” strategy [42, 118, 117,
87, 16]. This idea has been successfully implemented on massively parallel
Graphics Processing Units (GPUs), with modern solvers now handling millions of
particles in seconds through a convenient Python interface [29, 45]. These
computational gains have led to the development of several scalable methods
based on optimal transport for e.g. 3D shape registration [43, 124],
transcriptomics [116] and astrophysics [78]. As a direct inspiration for our
work, Laguerre tessellations and their connections with optimal transport have
been applied recently to materials science [1, 14, 11]. Similarly to
biological cell aggregates, some materials like steel are made up of a
collection of tiny “crystals” with various shapes and sizes. The ability to
generate in silico structures for an arbitrary set of crystal grains is key to
understanding the mechanical properties of several alloys.

Our framework relies extensively on these progresses and lies at the
intersection of materials science and computational fluid dynamics. In
particular, we show that Laguerre tessellations let us model general shapes:
appropriate transport cost functions encode arbitrary deformability
properties, in a dynamical setting and with a strict preservation of the
volume constraints. Although our model is formally a tessellation model, we
show that it shares important properties with particle systems, level-set
methods and vertex models. As a consequence, our model is remarkably
versatile: within the same framework, we can represent individual particles
with arbitrary shapes (as in level-set methods), their collective motion (as
in point-particle systems) and much denser tissue-like aggregates whose
dynamics is ruled by surface tension and other contact-based interactions (as
in vertex and finite elements models).

######  Remark 1.

Throughout the article, we use the generic word “particle” that, in most
cases, will refer to a biological cell. However, our framework is quite
general and potentially applicable to other contexts: a particle can also
occasionally refer to a bacteria, a pedestrian in a crowd or any general soft
body.

The article is organized as follows. The first section introduces the notion
of Laguerre tessellation in connection with optimal transport in a static
framework. The next section introduces the notion of incompressibility force
and the basic dynamical system formulation for soft bodies with arbitrary
deformability properties. The two last sections showcase the applicability of
our method in two important contexts: first, the collective motion of
particles induced by shape deformations due to collisions and external
stimuli; and second, the dynamics of multicellular aggregates whose self-
organized sorting patterns are induced by surface tension and surface
interactions. The supplementary material has three independent parts. The
first part presents additional numerical experiments, in particular for soft
active particles and crowd motion. The second part investigates mathematically
the formal theoretical links between our agent-based framework and new
macroscopic continuum PDE models via the notion of mean-field limit. Finally,
although the main text of the article is meant to be sufficiently self-
contained, the last part gathers reminders on optimal transport theory and its
numerical implementation.

##  2 Model: Laguerre Tessellation

###  2.1 Particles with a volume

Let us consider $N$ point particles in a bounded domain
$\Omega\subset\mathbb{R}^{d}$, with positions
$\mathcal{X}_{N}=(\bm{x}_{1},\ldots,\bm{x}_{N})$ and volumes denoted by
$\mathcal{V}_{N}=(v_{1},\ldots,v_{N})$. Each particle, generically labelled by
$i\in\\{1,\ldots,N\\}$, is therefore not reduced to a point but occupies a
certain subset $\mathscr{S}_{i}\subset\Omega$ of the available space. The
position $\bm{x}_{i}$ of a particle $i$ should be understood as a positional
indication but as we shall see later, it does not necessarily correspond to
the barycenter of the set $\mathscr{S}_{i}$. It is however directly comparable
to the seed or generator point of a Voronoi tessellation.

Throughout this article, we assume that the domain denoted by $\Omega$ is
connected and bounded, with total volume $|\Omega|=1$ and typical length
$L=1$. The centroid positions $\bm{x}_{i}$ are assumed to be properly
adimensionalized with respect to this length scale. The notation $|\cdot|$
denotes indifferently the volume of a set, understood as its Lebesgue measure,
or the Euclidean norm of a vector. Then we make the following two assumptions.

  * •

Non-overlapping: the sets $\mathscr{S}_{i}$ are disjoint.

  * •

Incompressibility: the volume of each particle is fixed
$|\mathscr{S}_{i}|=v_{i}$.

Naturally, these constraints imply that the particles need to adapt their
shapes to their current environment, especially in crowded and dynamical
situations. The main contribution of the present work is to introduce a new
way of computing optimal shapes depending on the current configuration and on
the individual deformability properties of the particles.

From a mathematical point of view, we are looking at a partition of the space
$\Omega=\cup\mathscr{S}_{i}$ into disjoint components with given volumes
$|\mathscr{S}_{i}|=v_{i}$ and we denote by $V=v_{1}+\ldots+v_{N}$ the total
volume that is occupied by the cells. When cells form a tissue, it is natural
to expect that $V=1$: the particles fill all the available space. In other
situations, there may be some free space between the particles and $V<1$: for
instance when cells are migrating in a medium. In this case, the extra-space
will be seen as one additional component of the partition of the space, i.e.
we add one additional particle with volume $1-V$. Another approach would be to
discretize this free space into a certain number of extra particles each with
a very small volume compared to the $v_{i}$’s. In our framework, both
approaches are actually equivalent as shown in [80]. Using the first method
with a single big particle is much more efficient numerically; on the other
hand, discretizing the medium with tiny particles lets us model a fluid
behaviour that may be relevant in some settings.

###  2.2 Static framework

####  2.2.1 Semi-discrete optimal transport framework

Two important classes of off-lattice models for multicellular systems use the
same idea of building a partition of the space. First, vertex models define a
partition through an arbitrary number of vertices and edges that create a
meshing of the space. Second, (generalized) Voronoi diagrams construct a
partition by solving an optimization problem to assign each point of the space
to one of several “seed” points. As we explain below, our framework extends
this class of models by casting such an optimization procedure as an optimal
transport problem.

A partition of the space that satisfies the incompressibility constraint can
be seen as a map $T:\Omega\to\\{\bm{x}_{i}\\}_{i}$ (or more generally
$T:\Omega\to\Omega$) that satisfies the volume constraint
$|\mathscr{S}_{i}|=v_{i}$ where $\mathscr{S}_{i}:=T^{-1}(\\{\bm{x}_{i}\\})$ is
the set of points in $\Omega$ that are sent to $x_{i}$ by $T$. We call such a
map a partition map. A simple but remarkable property, which is the basis of
our work, is that a partition map can also be understood as the solution of a
semi-discrete transport problem. First, let us state the following lemma –
reminders about measure theory and optimal transport are gathered in SM
Appendix 3:

######  Lemma 2.

Let $\hat{\mu}=\sum_{i=1}^{N}v_{i}\delta_{\bm{x}_{i}}$ be a discrete measure
with weight $v_{i}$ at location $x_{i}$. Then a map $T:\Omega\to\Omega$ is a
partition map if and only if $T\\#\mathrm{Leb}=\hat{\mu}$, i.e. if it pushes
the Lebesgue measure of $\Omega$ onto $\hat{\mu}$.

###### Proof.

For a given bounded continuous function $\phi\in C_{b}(\Omega)$, we recall the
definition of the push-forward measure $T\\#\mathrm{Leb}$,

$\int_{\Omega}\phi(\bm{x})T\\#\mathrm{Leb}(\mathrm{d}\bm{x})=\int_{\Omega}\phi(T(\bm{x}))\mathrm{d}\bm{x}.$

Thus, if $T$ is a partition map, its image set is the discrete set
$\\{\bm{x}_{i}\\}_{i}$ so

$\int_{\Omega}\phi(\bm{x})T\\#\mathrm{Leb}(\mathrm{d}\bm{x})=\sum_{i=1}^{N}\phi(\bm{x}_{i})|T^{-1}(\\{\bm{x}_{i}\\})|=\sum_{i=1}^{N}v_{i}\phi(\bm{x}_{i})=\int_{\Omega}\phi(x)\hat{\mu}(\mathrm{d}\bm{x}).$

Reciprocally, if $T\\#\mathrm{Leb}=\hat{\mu}$ then by choosing a test function
whose support does not contain any $\bm{x}_{i}$, one can see that the image
set of $T$ is contained in $\\{\bm{x}_{i}\\}_{i}$. Then, for any $i$ and using
the same formula as above, by choosing a test function $\phi$ which is equal
to one on $\bm{x}_{i}$ and 0 on all the other $\bm{x}_{j}$, it is easily
checked that $|T^{-1}(\\{\bm{x}_{i}\\})|=v_{i}$. ∎

Without further hypothesis, there exists an infinite number of partitions of
the space. With the previous measure transportation formulation in mind, we
consider the partition map that is solution of the “Monge” optimization
problem from optimal transport theory:

$\min_{\begin{subarray}{c}T:\Omega\to\Omega\\\
T\\#\mathrm{Leb}=\hat{\mu}\end{subarray}}\,\,\sum_{i=1}^{N}\int_{T^{-1}(\\{\bm{x}_{i}\\})}c(\bm{x},\bm{x}_{i})\mathrm{d}\bm{x},$
(1)

where the function $c:\Omega\times\\{\bm{x}_{i}\\}_{i}\to[0,+\infty)$ is
called the cost function. Provided that this optimization problem has a unique
solution for a given cost, this provides an arbitrary but clearly defined way
of computing a partition of the space. As shown in the following theorem, it
is well known that the semi-discrete Monge problem (1) has a unique solution
under mild hypotheses on the cost function.

######  Theorem 3.

Let $\hat{\mu}=\sum_{i=1}^{N}v_{i}\delta_{\bm{x}_{i}}$ be a discrete
probability measure on $\Omega$. Let $c$ be a cost function such that for all
particle $i$, $c_{i}:y\in\Omega\mapsto c(y,\bm{x}_{i})\in[0,+\infty)$ belongs
to $C^{1,1}(\Omega)$. Then the Monge problem (1) has a unique solution $T$
given by

$T:\Omega\to\\{\bm{x}_{i}\\}_{i},\,\,\bm{x}\in\mathscr{L}_{i}\mapsto\bm{x}_{i}$

where the partioning sets $(\mathscr{L}_{i})_{i}$ are the (generalized)
Laguerre cells defined by

$\mathscr{L}_{i}:=\\{x\in\Omega\,:\,c(\bm{x},\bm{x}_{i})-w_{i}\leq
c(\bm{x},\bm{x}_{j})-w_{j}\,\,\,\forall j\\}.$ (2)

The intersection of two Laguerre cells has zero Lebesgue measure and the
weights $w_{i}\in\mathbb{R}$, called Kantorovich potentials, are uniquely
defined such that for all $i\in\\{1,\ldots,N\\}$,

$|\mathscr{L}_{i}|=v_{i}.$ (3)

Optimal transport solvers are computer programs that take as input the weights
$v_{1}$, …, $v_{N}$, the cost functions $c_{1},\dots,c_{N}$, and return as
output the unique vector of optimal values for the Kantorovitch potentials
$(w_{1},\dots,w_{N})$.

This result is well-known and proved for instance in [13, 72] (see also SM
Appendix 3 for more details) where it is applied to computational fluid
dynamics and material science problems. To the best of our knowledge, it is
the first time that this result is used for biological applications. In
particular, although the result is valid for very general cost functions, most
applications consider the $L^{2}$ cost function

$c(\bm{x},\bm{y})=|\bm{y}-\bm{x}|^{2}.$ (4)

One of the main focus of the present work is to investigate how the choice of
an appropriate cost function lets us model various biophysical properties and
behaviors. As a first static example, Fig. 2 shows three examples of a
tessellation of a square periodic domain obtained by Theorem 3 with different
cost functions and different particles shapes.

(a) Voronoi

(b) Bubbly

(c) Random

Figure 2: (a) A Voronoi tessellation obtained with the $L^{2}$ cost of Eq. (4)
and 42 particles with random volumes $v_{i}$ sampled uniformly between 0.2 and
1 (before normalization to sum up to 1). (b) A bubbly tessellation similar to
[64] obtained for 66 particles with random volumes sampled uniformly between
0.05 and 1 (before normalization) and a power cost
$c(\bm{x},\bm{x}_{i})=|\bm{x}-\bm{x}_{i}|^{\alpha_{i}}$ where the $\alpha_{i}$
are distributed uniformly between 0.5 and 4. Lighter colors indicate lower
values of $\alpha_{i}$ and correspond to softer shapes. (c) A Voronoi
tessellation with random fluctuations similar to [91] obtained with 42
identical particles and a perturbed $L^{2}$ cost as explained in Section
4.2.2.

######  Remark 4 (Scaling).

Given a cost $c$, one can define the energy of a partition
$\mathscr{S}=(\mathscr{S}_{i})$ of $\Omega$ by

$\mathcal{E}(\mathscr{S}):=\sum_{i}\int_{\mathscr{S}_{i}}c(\bm{x},\bm{x}_{i})\mathrm{d}\bm{x}.$
(5)

Theorem 3 shows that for a given set of points $\bm{x}_{i}$ and volumes
$v_{i}$, the lowest energy is the transport cost between the Lebesgue measure
and the discrete empirical measure $\hat{\mu}$ :

$\mathcal{T}_{c}(\hat{\mu})=\min_{\begin{subarray}{c}\cup\mathscr{S}_{i}=\Omega\\\
|\mathscr{S}_{i}|=v_{i}\end{subarray}}\mathcal{E}(\mathscr{S}).$

With this remark in mind, a natural scaling for the cost is to impose a finite
energy to a base shape. Unless otherwise specified, the cost
$c(\bm{x},\bm{x}_{i})$ will always be rescaled so that, given a dimensionless
scaling parameter $\lambda_{i}>0$ and a base shape $\mathscr{S}_{0}$ with
volume $V_{0}$,

$\int_{\mathscr{S}_{0}}c(\bm{x},\bm{x}_{i})\mathrm{d}\bm{x}=\lambda_{i}V_{0}.$
(6)

####  2.2.2 Laguerre tessellation as generalized Voronoi diagrams with volume
constraints

The Laguerre tessellation of the space given by the sets (2) is reminiscent of
the notion of generalized Voronoi cells defined by

$\mathscr{V}_{i}=\\{\bm{x}\in\Omega,\,\,d(\bm{x},\bm{x}_{i})<d(\bm{x},\bm{x}_{j})\,\,\forall
j\\},$

where $d$ is an arbitrary generalized distance function. In a classical
Voronoi tessellation, the function $d$ is equal to the Euclidean distance
$d(\bm{x},\bm{x}_{i})=|\bm{x}-\bm{x}_{i}|$ and the resulting Voronoi cells are
convex polytopes. For biological applications, this is not always realistic
[70], as cell aggregates sometimes exhibit curved boundaries, possibly even
with random fluctuations [91]. Following this remark, some motivating examples
for the concept of generalized Voronoi diagram are discussed in particular in
[10, 64, 56] and in the references therein. Replacing the Euclidean distance
by an arbitrary function $d$ indeed naturally allows other boundary shapes,
and in particular circular arcs, a topic which is discussed in details in [10]
using the generalized distance functions

$d(\bm{x},\bm{x}_{i})=|\bm{x}-\bm{x}_{i}|^{2}-s_{i}^{2},\quad
d(\bm{x},\bm{x}_{i})=|\bm{x}-\bm{x}_{i}|^{2}/s_{i}^{2},$ (7)

where $s_{i}\in\mathbb{R}$ is a tunable scaling parameter.

Formally, the Laguerre cells defined by Eq. (2) could be seen as generalized
Voronoi cells with the generalized distance function
$d(\bm{x},\bm{x}_{i})=c(\bm{x},\bm{x}_{i})-w_{i}$. Since, at this point, the
cost function $c$ is as arbitrary as the generalized distance function $d$, in
comparison to the introduction of [10], this discussion is mostly about naming
conventions. More importantly, our optimal transport perspective provides an
original point of view on the three following topics:

  1. 1.

The Laguerre cells in Theorem 3 are generalized Voronoi cells that satisfy
strict volume constraints Eq. (3) thanks to an optimal choice of the
Kantorovich potentials $w_{i}$. Note that in [10], the scaling parameter
$s_{i}$ in the definition of generalized distances is also introduced as a way
of controlling the size of the particles – but without a formal link to strict
volume constraints. Since volume is often the first experimentally measurable
quantity, it is desirable to include it as an easily controllable degree of
freedom of the model.

  2. 2.

While most of the literature in computational geometry focuses on the square
Euclidean cost (as in [10]), the optimal transport interpretation naturally
suggests other choices that we investigate below. As shown in Section 2.2.3,
this lets us encompass both tissue-like aggregates (for which generalized
Voronoi tessellations are commonly accepted as accurate models) and individual
particles with arbitrary deformable shapes (which are often modeled using
completely different methods).

  3. 3.

The numerical computation of (generalized) Voronoi tessellations is a
classical but non-trivial problem in computational geometry. Based on the
optimal transport interpretation, we propose an alternative and novel
numerical method that is applicable to arbitrary cost functions in Section
2.3.

####  2.2.3 Laguerre cells as level-sets

As mentioned in the introduction, when the particles do not fill the entire
space, a simple strategy is to add an additional particle labelled by $i=0$
with volume $v_{0}=1-\sum_{i\geqslant 1}v_{i}$ at an arbitrary position
$\bm{x}_{0}\in\Omega$. The exact position of this particle will not matter as
it will be associated to the cost $c(\bm{x},\bm{x}_{0})=0$ for all
$\bm{x}\in\Omega$. With this choice, we may consider individual isolated
particles that only share a boundary with a surrounding medium modeled by the
particle $\bm{x}_{0}$. A direct corollary of Theorem 3 shows that the shape of
an individual particle is always a level-set of the cost function.

######  Lemma 5.

The shape $\mathscr{S}\subset\Omega$ that minimizes the individual energy

$\mathcal{I}(\mathscr{S})=\int_{\mathscr{S}}c(\bm{x},\bm{x}_{i})\mathrm{d}\bm{x}$

is of the form

$\mathscr{S}(w)=\\{\bm{x}\in\Omega,\,\,c(\bm{x},\bm{x}_{i})\leq w\\},$

for some $w>0$.

Note that with the notations of Theorem 3, one can write more precisely the
value $w$ as the difference between the Kantorovich potentials
$w=w_{0}-w_{i}$. Since this is an important property, we provide a direct
self-contained proof:

###### Proof.

First of all, since the function $U:w\mapsto|\mathscr{S}(w)|$ is non-
decreasing, $U(0)=0$ and $U(w)\to 1$ as $w\to+\infty$ there exists a unique
$w^{*}>0$ such that $U(w^{*})=v_{i}$. Then let $\mathscr{S}\subset\Omega$ be
any shape of volume $v_{i}$. Denoting $\mathscr{S}^{*}=\mathscr{S}(w^{*})$, it
holds that

$\displaystyle\mathcal{I}(\mathscr{S})$
$\displaystyle=\mathcal{I}(\mathscr{S}\cap\mathscr{S}^{*})+\mathcal{I}(\mathscr{S}\cap(\mathscr{S}^{*})^{c})$
$\displaystyle=\mathcal{I}(\mathscr{S}^{*})+\mathcal{I}(\mathscr{S}\cap(\mathscr{S}^{*})^{c})-\mathcal{I}(\mathscr{S}^{*}\cap\mathscr{S}^{c})$

By assumption $|\mathscr{S}^{*}|=|\mathscr{S}|=v_{i}$ so it holds that

$|\mathscr{S}\cap(\mathscr{S}^{*})^{c})|=|\mathscr{S}^{*}\cap\mathscr{S}^{c}|=:\tilde{v}_{i}\geq
0,$

and by definition of $\mathscr{S}^{*}$,

$\mathcal{I}(\mathscr{S}\cap(\mathscr{S}^{*})^{c})\geq\tilde{v}_{i}w^{*},\quad\mathcal{I}(\mathscr{S}^{*}\cap\mathscr{S}^{c})\leq\tilde{v}_{i}w^{*},$

which gives

$\mathcal{I}(\mathscr{S})\geq\mathcal{I}(\mathscr{S}^{*}).$

∎

Consequently, any particle shape can be realized using an appropriate cost
function $c_{i}:\bm{x}\in\Omega\mapsto c(\bm{x},\bm{x}_{i})\in[0,+\infty)$.
These potentials are analogous to the signed distance functions that are
commonly used for volumetric rendering in computer graphics [110] – with the
added benefit that the derivative of our cost encodes the local deformability
of the shape. For instance, in dimension 2, any shape defined by a polar
equation $r=r_{0}(\theta)$ can be encoded by the cost

$c(\bm{x},\bm{x}_{i})=\left(\frac{|\bm{x}-\bm{x}_{i}|}{r_{0}(\theta(\bm{x},\bm{x}_{i}))}\right)^{p},$
(8)

where $\theta(\bm{x},\bm{x}_{i})$ denotes the polar angle of the vector
$\bm{x}-\bm{x}_{i}$ and $p>0$ is a hardness parameter whose influence is
illustrated in Figures 2 and 4. Indeed, all the level-sets of the function (8)
are scaled version of the base shape $r=r_{0}(\theta)$, and the only one with
the correct volume is automatically selected by the Kantorovich potential $w$
in Lemma 5. Note that other cost functions could be considered, as long as the
desired shape is a level-set.

For any exponent $p$, this cost satisfies the scaling (6) with
$\lambda_{i}=2/(p+2)$ and the base shape $\mathscr{S}_{0}$, the boundary of
which is defined by the polar equation $r=r_{0}(\theta)$. This construction
extends to any dimension $d$ when $r_{0}$ defines the equation of a closed
$(d-1)$-dimensional surface (in which case $\lambda_{i}=d/(p+d)$).

Representing individual particles with arbitrary shapes as level-sets of some
custom functions is a well-known method in computational biology, see for
instance [140, 93]. Although this idea appears naturally in our framework, it
should be noted that it formally comes from a partition of the space and not
from an interface tracking method as in traditional level-set methods [140,
93]. This has a clear advantage since when two particles are in contact, there
is no particular surgery to consider when computing their contact surface
which is simply the intersection of two Laguerre cells

$\mathscr{L}_{i}\cap\mathscr{L}_{j}\subset\\{\bm{x}\in\Omega,\,\,c(\bm{x},\bm{x}_{i})-w_{i}=c(\bm{x},\bm{x}_{j})-w_{j}\\}.$

Consequently, our framework handles well individual particles, dense tissue-
like aggregates or a mix of both.

###  2.3 Implementation

The pioneering work of Brenier [16, 18, 19] on the connection between optimal
transport and incompressible fluid flows motivated the development of fast and
efficient numerical methods to compute Laguerre tessellations, with remarkable
achievements in computational fluid dynamics [46, 80, 36]. In this context, a
Laguerre tessellation is comparable to a meshing procedure [81]; individual
mesh elements are Lagrangian particles that represent a microscopic portion of
a fluid.

More recently, Lagrangian tessellations have been introduced in materials
science as a description of polycrystalline materials such as steel [1, 11,
14]. In this context, the Laguerre cells represent crystal grains and Laguerre
tessellations provide a static description of the material. A topic of
interest is the generation of tessellations with an appropriate spatial
distribution of grains of different sizes and shapes, as different grain
arrangements will impact the mechanical properties of the full material. Our
framework lies at the intersection of these two ideas: we generate
macrosocopic particles with arbitrary shapes (analogous to crystal grains) but
in a dynamical framework since biological organisms evolve by responding to
internal and external forces in a fluid-like manner.

Whether it is in computational fluid dynamics or in material science, the
semi-discrete optimal transport problem is considered mostly for the $L^{2}$
cost $c(\bm{x},\bm{y})=|\bm{x}-\bm{y}|^{2}$. Extremely efficient numerical
methods, based on a damped Newton algorithm, have been developed for this case
and are now available as open source software [79, 12, 88, 90]. Unfortunately,
these implementations leverage the structure of the squared Euclidean distance
to compute Laguerre cells efficiently [111] and cannot be applied with general
cost functions.

This prevents us from using them for biological applications: the freedom in
the choice of arbitrary, individualized cost functions $c_{i}$ is of upmost
importance to us. To unlock the use of Laguerre cells with general shapes and
deformability properties, we work in a fully discrete setting with a discrete
approximation of the Lebesgue measure

$\hat{\mu}_{0}=\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{y}_{j}}\approx\mathrm{Leb},$

that is supported by “pixel centers” $\bm{y}_{j}$ on a uniform grid
discretization of the space $\Omega=[0,1]^{d}$. We choose $M\gg N$ so that
each particle is effectively discretized with a sufficiently high number of
pixels $\bm{y}_{j}$ with volume $1/M$. Solving discrete optimal transport
problems is by now classical owing to the numerous applications in data
science [108]. Our optimal solver is based on a numerically stable,
symmetrized and annealed variant of the renowned Sinkhorn algorithm that is
described in SM Appendix 3. In order to handle arbitrary cost functions and
fine discretizations with $M\sim 10^{6}$ pixels, we rely on the semi-symbolic
lazy tensors introduced by the KeOps and GeomLoss libraries [29, 45, 44]. This
lets our plain Python code perform scalable geometric computations on the GPU,
with automatic differentiation and without memory overflows.

A drawback of our discretization procedure is that there is generally no Monge
map between two discrete measures. As a replacement for the Monge map defined
in Theorem 3, we consider the following discrete allocation map

$T:\\{\bm{y}_{j}\\}_{j\in\\{1,\ldots,M\\}}\to\\{\bm{x}_{i}\\}_{i\in\\{1,\ldots,N\\}},\quad\bm{y}_{j}\mapsto\bm{x}_{j^{*}},$
(9)

with

$j^{*}:=\operatorname*{arg\,min}_{i\in\\{1,\ldots,N\\}}\\{c(\bm{y}_{j},\bm{x}_{i})-w_{i}\\}.$
(10)

Other choices such as the popular barycentric mapping [108, Remark 4.11] could
also be considered, but the definition (9-10) directly extends that of Theorem
3 and provides the best results in practice.

We use NumPy, PyTorch, KeOps and GeomLoss for our simulations [55, 106, 29,
44, 42] while relying on Matplotlib, PyVista and Paraview for visualizations
[61, 127, 4]. Our source code, containing scripts to reproduce all of our
numerical experiments, is freely available at

https://github.com/antoinediez/ICeShOT

###  2.4 Dynamical setting

So far, the centroid positions $\bm{x}_{i}$ and volumes $v_{i}$ are given and
Theorem 3 shows that for any such configuration, there is always a unique
partition of the space that satisfies the volume constraint. In a dynamical
setting (discretized in time), at each discrete time step, forces are applied
to the centroid positions $\bm{x}_{i}$ and a new Laguerre tessellation of the
space is recomputed. We can also update the volumes and/or the cost functions
to model morphological changes, or to include other events such as the
creation or deletion of particles. Similarly to (generalized) Voronoi models
[9], the computation of the space tessellation at each time step is the most
expensive part. However, thanks to the GPU implementation discussed above,
each step typically takes from a few seconds to less than a minute, which
remains manageable.

Importantly, the forces applied to the centroids or the evolution of the cost
function can depend on the geometrical properties of the Laguerre
tessellation. Extracting the relevant information from the Laguerre cells and
translating it into biophysical properties is the main modeling task. It is
the object of the remaining of this article, which is organized along the
following three main ideas.

  1. 1.

A repulsion force ensures a realistic spatial distribution of the particles.
The competition between the repulsion force and the other forces depend on how
easy particles can deform or be displaced, which is encoded in the cost
function.

  2. 2.

Cost functions can be enriched by considering for instance a polarity, a
deformation bias in some directions or protrusions locations, in order to
model the collective motion of particles with a complex oriented structure.

  3. 3.

Interaction forces can depend on the contact surface between the Laguerre
cells in order to model surface tension phenomena in cellular aggregates.

##  3 Incompressibility force and deformation cost

###  3.1 Damped Lloyd algorithm

In a dynamical framework, a reasonable assumption is that, in the absence of
other forces, each particle should relax towards its preferred shape (computed
by Lemma 5). In the formula Eq. (5) of the energy of a partition, each term
has an absolute minimum when $\mathscr{S}_{i}$ is the shape given by Lemma 5
(that is when the particle is only in contact with the medium associated to
the zero cost). Consequently, for the whole particle system, the most stable
centroid configuration $\bm{x}_{1},\ldots,\bm{x}_{N}$ is the one which
minimizes the total cost

$\min_{\bm{x}_{1},\ldots,\bm{x}_{N}}\mathcal{T}_{c}(\hat{\mu}).$

This suggests to introduce the following gradient descent motion for each
particle:

$\dot{\bm{x}}_{i}=-\tau_{i}\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu}).$
(11)

The gradient descent step $\tau_{i}>0$ may depend on each particle and can be
interpreted as an inertia parameter controlling how easy it is to move a
particle from its current location.

######  Remark 6.

As it is customary in computational biology, particles are assumed to evolve
in an over-damped regime so that only first order equations are considered.

The force on the right-hand side of Eq. (11) can also be understood as a
repulsion force between the particles (see Fig. 3 and SM Video 1). However,
this force strongly differs from the classical repulsion forces in point
particle systems, which are typically modeled as a sum of binary interactions

$\dot{\bm{x}}_{i}=-\tau_{i}\nabla_{\bm{x}_{i}}\sum_{j=1}^{N}V(\bm{x}_{j}-\bm{x}_{i}),$

where $V$ is a repulsive potential. Contrary to this repulsion force, the
gradient descent Eq (11) takes into account the shapes and the global
configuration of the system. In computational fluid dynamics [46], the
gradient of the transport cost can be directly related to the pressure
gradient in an incompressible flow. We will therefore rather refer to
$-\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu})$ as an incompressibility
force.

(a) $t=0$

(b) $t=99$

Figure 3: (a) A four-cell configuration showing the incompressibility forces
acting on the centroids (red arrows) for the $L^{2}$ cost. (b) Equilibrium
configuration after 99 Euler discretization steps ($\Delta t=0.1$) of Eq.
(11). See also SM Video 1

On the practical side, it has been shown in [139] that the gradient of the
total transport cost has a simple computable expression, namely

$\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu})=\int_{\mathscr{L}_{i}}\nabla_{\bm{x}_{i}}c(\bm{x}_{i},\bm{x})\mathrm{d}\bm{x}.$

In particular, when $c(\bm{x},\bm{y})=\ell(|\bm{y}-\bm{x}|)$ is an isotropic
cost which depends only on the distance through a function
$\ell:\mathbb{R}\to\mathbb{R}$,

$\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu})=\int_{\mathscr{L}_{i}}\frac{\ell^{\prime}(|\bm{x}-\bm{x}_{i}|)}{|\bm{x}-\bm{x}_{i}|}(\bm{x}_{i}-\bm{x})\mathrm{d}\bm{x}.$

It is a force in the direction of the weighted barycenter $\bm{\bar{x}}_{i}$
of the Laguerre cell

$\bm{\bar{x}}_{i}:=\frac{\int_{\mathscr{L}_{i}(w)}\frac{\ell^{\prime}(|\bm{x}-\bm{x}_{i}|)}{|\bm{x}-\bm{x}_{i}|}\bm{x}\,\mathrm{d}\bm{x}}{\int_{\mathscr{L}_{i}(w)}\frac{\ell^{\prime}(|\bm{x}-\bm{x}_{i}|)}{|\bm{x}-\bm{x}_{i}|}\,\mathrm{d}\bm{x}}.$
(12)

In computational geometry, the gradient descent Eq. (11) is interpreted as a
damped Lloyd algorithm [139, 89] with damping parameter $\tau_{i}$. In
comparison, the classical Lloyd algorithm is an iterative algorithm commonly
used to compute well-shaped Voronoi tessellations by repositioning at each
iteration the centroids $\bm{x}_{i}$ at the barycenter $\bm{\bar{x}}_{i}$. In
other words, Lloyd algorithm is the explicit Euler discretization of Eq. (11)
with
$\tau_{i}^{-1}=\int_{\mathscr{L}_{i}(w)}\frac{\ell^{\prime}(|\bm{x}-\bm{x}_{i}|)}{|\bm{x}-\bm{x}_{i}|}\,\mathrm{d}\bm{x}$
and a unit time step $\Delta t=1$.

Although we will typically use Lloyd algorithm to compute realistic initial
configurations, the damped version Eq. (11) is preferable as it can be
advantageously coupled with additional force or noise terms, leading to
general mean-field like stochastic differential equation systems, for
instance:

$\mathrm{d}\bm{x}_{i}=\bm{b}(\bm{x}_{i},\hat{\mu})\mathrm{d}t-\tau_{i}\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu})\mathrm{d}t+\sqrt{2\sigma_{i}}\mathrm{d}\bm{B}^{i}_{t},$
(13)

where $\bm{b}$ is a force term which may depend on the other particles,
$\bm{B}^{i}_{t}$ is a Brownian motion and $\sigma_{i}>0$ is a diffusion
coefficient. In other words, the incompressibility force which takes into
account the shape of the particles and the global configuration can be added
to any point-particle dynamical system. Note that the cost function $c$ in
this equation can also depend on time and evolve depending on the geometrical
properties of the system, as explained in the following. In the rest of this
article, the dynamics of the system will always be described by a (stochastic)
differential system of the form Eq. 13 or a variation of it. The time
discretization is always a forward Euler method with a step size denoted by
$\Delta t$.

The natural question of the mean-field limit when $N\to+\infty$ will be
shortly discussed in SM Appendix 2.

######  Remark 7 (Boundary conditions).

By definition, the Laguerre cells are always contained in $\Omega$ and no
particular procedure needs to be considered to handle boundary conditions,
which allows arbitrary domain shapes (see for instance SM Video 2). It is also
possible to consider periodic boundary conditions by simply replacing the
Euclidean distance $|\cdot|$ everywhere it appears by the periodized version
$d_{P}$ along each dimension. For $\Omega=[0,1]^{d}$ with periodic boundaries,

$d_{P}(\bm{x},\bm{y})^{2}:=\sum_{k=1}^{d}\delta_{P}(x_{k},y_{k})^{2},$

where $\bm{x}=(x_{1},\ldots,x_{d})^{\mathrm{T}}$ and for two scalar values
$x,y\in[0,1]$,

$\delta_{P}(x,y):=\min(|y-x|,1-|y-x|).$

No particular constraint is imposed on the centroids $\bm{x}_{i}$. Actually,
Theorem 3 remains valid even if $\bm{x}_{i}\notin\Omega$, and it does not even
imply that $\bm{x}_{i}\in\mathscr{L}_{i}$ either. It will nevertheless always
be the case in merely all the examples presented in this article thanks to the
incompressibility force which tends to reposition the centroids at the
(weighted) barycenter of their Laguerre cell. Although it may seem
questionable at first sight, this fact may actually also be seen as a
(necessary) feature of the model to handle arbitrary deformations up to very
degenerate non-convex shapes where the centroids $\bm{x}_{i}$ should be
considered more as numerical parameters than as physical quantities.

###  3.2 Deformation and cost comparison

In order to assess the influence of the cost function in the gradient descent
Eq. (11), let us consider the family of isotropic power costs of the form

$c(\bm{y},\bm{x})=\lambda^{-1}_{\alpha}|\bm{y}-\bm{x}|^{\alpha},$ (14)

depending on the parameters $\alpha,\lambda_{\alpha}>0$. For any value of the
parameter $\alpha$, the base shape of the particles is always a ball. We
assume that all the particles have equal volumes $v_{i}=\pi R_{0}^{2}$ for a
given $R_{0}>0$ (in dimension 2) and we set the scaling parameter
$\lambda_{\alpha}=\frac{2}{\alpha+2}R_{0}^{\alpha}$ so that for any $\alpha$,

$\int_{B(\bm{x}_{i},R_{0})}c(\bm{x},\bm{x}_{i})\mathrm{d}\bm{x}=\pi
R_{0}^{2}=v_{i}.$

Moreover, the gradient descent step in Eq. (11) is normalized as
$\tau\equiv\tau/R_{0}$ so that the parameters are independent of the size
$R_{0}$ of the particles. Unless otherwise specified, we always use this
normalization throughout the article.

Although the individual cost of an isolated spherical particle does not depend
on $\alpha$, it will not be the case in crowded areas where particles need to
adopt deformed shapes. Indeed, when $\alpha>1$, the cost is a convex function
of the distance between the point $\bm{x}$ and the centroid position
$\bm{x}_{i}$. Large deformations thus require to allocate points $\bm{x}$ far
from the centroid position $\bm{x}_{i}$ and are thus strongly penalised. On
the contrary, when $\alpha<1$, the cost is a concave function of the distance,
which allows larger deformations and potentially shapes far from the base ball
shape. From this observation, the parameter $\alpha$ appears as a hardness
parameter. In the three numerical experiments below, we confirm this
interpretation by applying a constant force load to particles and measuring
the deformation response depending on $\alpha$. Other quantitative and
qualitative examples are presented in SM Appendix 1 (in particular for a
system of active Brownian particles and a model of crowd motion).

######  Remark 8 (Quotient method).

Another way to tune the deformation properties of the particles is to consider
the same cost for all the particles but different scaling factors
$\lambda_{i}$ in Eq. (6). This is similar to the so-called quotient method in
[10]. An example is shown in the SM Video 2.

####  3.2.1 Falling soft spheres

We first apply a constant downward force to the particles, that is Eq. (13)
with $\sigma_{i}=0$ and $\bm{b}=-F_{0}\bm{\vec{e}}_{z}$. In a crowded
situation, illustrated in Fig. 4 and SM Videos 4,5,6, there is a competition
between the repulsive incompressibility force and this gravity-like force.
When $\alpha$ is small, the particles manage to squeeze their way by adopting
elongated columnar shapes. When $\alpha$ increases, particles behave like
hard-spheres and end up in a typical hexagonal packing situation.

Figure 4: Initial configuration (top) and equilibrium configuration (bottom)
of a system of falling soft spheres for three values of the deformation
parameter $\alpha$. For $\alpha=0.75$, large deformations lead to more
elongated and non-convex shapes (left). For $\alpha=2$, the particles are
convex and the equilibrium configuration is close to an hexagonal lattice
(middle). For $\alpha=10$, deformations from the base spherical shape are
strongly penalized and the particles thus behave as hard-spheres (right).
Parameters: grid size $512\times 512$, $\tau=1.5$, $F_{0}=0.4$, $\Delta
t=0.001$. See also SM Videos 4,5,6.

####  3.2.2 Mixed spheres

Contrary to the cost function which encodes the deformability properties, the
gradient step $\tau_{i}$ in Eq. (13) represents the magnitude of the
incompressibility force (relatively to the other forces) and may be
interpreted as an inertia parameter, typically the inverse of a mass. To
illustrate this, we consider the same system of falling spheres but with half
of the spheres being softer and various values of the parameter $\tau$, which
impact the final configuration (see Fig. 5 and SM Videos 7,8,9).

(a) Light ($\tau_{o}=8$)

(b) Same ($\tau_{o}=3$)

(c) Heavy ($\tau_{o}=1$)

Figure 5: A mixed system of 30 falling spheres with 15 harder blue spheres
($\alpha=2$) and 15 softer orange spheres ($\alpha=1$). The gradient descent
step of the blue spheres is fixed to $\tau_{b}=3$ and both the orange and blue
spheres are subject to the same downward force with magnitude $F_{0}=0.4$. The
gradient descent step of the orange particles $\tau_{o}$ is analogous to the
inverse of a mass. Consequently heavier orange particles tend to be lower.
Parameters: grid size $512\times 512$, $\Delta t=0.001$. See also SM Videos
7,8,9.

####  3.2.3 Micropipette aspiration

Viscoelastic properties of individual biological cells are often quantified
using micropipette aspiration techniques [54, 50]: a single cell is first
placed at the tip of a thin micropipette tube and a controlled pressure
difference then creates an aspiration force which sucks the cell inside the
micropipette. The biomechanical properties are quantitatively evaluated by
measuring the portion of the cell which effectively travels through the tube,
this aspiration length ranging from zero for solid-like cell to the full tube
for liquid-like cells. This experiment can be easily mimicked in silico by
considering a micropipette-shaped domain $\Omega$. Then, for a given set of
fixed parameters (cell size, micropipette width, force magnitude $\tau_{i}$),
increasing the value of the deformability parameter $\alpha$ indeed allows to
interpolate between particles with a liquid behaviour to particles which acts
as solid spheres. Tuning appropriately the cost function thus allows to
reproduce any visco-elastic cell (see Fig. 6 and SM Video 10).

(a) $t=0$

(b) Equilibrium

(c)

Figure 6: Micropipette experiment. (a)-(b) A single cell is sucked inside a
thin tube by a constant force until it is blocked by the incompressibilty
force ($\alpha=1$). (c) Aspiration length over time between 0 (outside the
tube) and 1 (completely inside the tube) for different values of the
deformability parameter $\alpha$. See also SM Video 10.

##  4 Dynamical costs

In this section, we present two ways of making shapes evolve with time by
dynamically changing the properties of the cost function.

###  4.1 Orientation and alignement

The position $\bm{x}_{i}$ of a particle is often insufficient to describe the
motion of active particles which also need to be endowed with an orientation.
Typically, the orientation is a unit vector $\bm{v}_{i}$ which defines the
direction of motion of the particle, as in flocking models. In cell biology,
under the generic name of polarity, this vector may encode various structural
differences which specify e.g. the migration direction, the division axis or
some shape anisotropy. In this latter case, the anisotropy can be further
quantified by introducing an aspect ratio $a_{i}:b_{i}$ in dimension 2 or more
generally a covariance matrix $\Sigma_{i}$ in any dimension, which we always
consider normalized to have determinant one in the following. We recall that
in dimension 2, the matrix $\Sigma_{i}$ is defined in terms of the polarity
vector $\bm{v}_{i}$ and the aspect ratio $s_{i}=a_{i}/b_{i}>1$ by the
eigenvector decomposition

$\Sigma_{i}\bm{v}_{i}=\frac{a_{i}}{b_{i}}\bm{v}_{i},\quad\Sigma_{i}\bm{v}^{\perp}_{i}=\frac{b_{i}}{a_{i}}\bm{v}^{\perp}_{i},$

where $\bm{v}_{i}^{\perp}$ is the direct orthogonal vector to $\bm{v}_{i}$.

Oriented particles with anisotropic shapes can be realized in our optimal
transport framework by considering appropriate cost functions which depend on
$s_{i},\bm{v}_{i}$ or $\Sigma_{i}$. For instance, the following cost is a
direct generalization of the isotropic $L^{2}$ cost Eq. (14) with scaling
factor $\lambda>0$,

$c(\bm{x},\bm{y})=\lambda^{-1}(\bm{y}-\bm{x})^{\mathrm{T}}\Sigma^{-1}(\bm{y}-\bm{x}).$

It corresponds to an ellipsoid shape with covariance matrix $\Sigma$. More
refined shapes can be realized using cost functions of the form Eq. (8) with
an appropriate function $r_{0}$. A classical example in biology is the
spherocylinder shape of 2D bacilli. Given a long axis vector
$(\cos\theta_{0},\sin\theta_{0})^{\mathrm{T}}$ and an aspect ratio $s=a/b\geq
1$, this shape can be realized using the cost function associated to the
function $r_{0}(\theta)$ defined by

$r_{0}(\theta)=\left\\{\begin{array}[]{l}\frac{b}{|\sin(\theta-\theta_{0})|},\quad\text{if}\,\,|\tan(\theta-\theta_{0})|\geq\frac{1}{s-1},\\\
b\left(|\cos(\theta-\theta_{0})|(s-1)+\sqrt{1-((s-1)\sin(\theta-\theta_{0}))^{2}}\right),\quad\text{otherwise.}\end{array}\right.$

Coming back to a particle system, each particle is now endowed with its own
cost function which may depend on the particle orientation vector
$\bm{v}_{i}$, aspect ratio $s_{i}\geq 1$ or covariance matrix $\Sigma_{i}$.
Just as the positions $\bm{x}_{i}$ which evolve in time, there are several
possible ways to model the evolution of these shape and orientation
parameters. Since the pioneering work of Vicsek [136, 30], in the literature
on point particle systems, a popular class of models consider local alignment
interactions by introducing a force typically of the form

$\dot{\bm{v}}_{i}=\left(\mathrm{Id}-\bm{v}_{i}\bm{v}_{i}^{\mathrm{T}}\right)\sum_{j\sim
i}\bm{v}_{j},$

where the sum is over the indices $j$ of the particles which are close to
$\bm{x}_{i}$. This kind of interaction could of course be directly implemented
in our model but it does not fully take into account the shapes of the
particles. Actually, in polymer physics or in the modelling of swarms of rod-
like bacteria, the local alignment observed experimentally is often understood
as a consequence of collisions events and volume exclusion effect between
anisotropic shapes. In the Vicsek model and its numerous variants which focus
on the consequences of alignment, it is then assumed as a model feature. Since
our modelling framework can fully take into account the shape of the
particles, we rather want to focus on the cause of alignment and confirm the
idea that it emerges from the sole volume exclusion interactions and a natural
assumption of local re-orientation due to collisions. A similar idea has
recently been developed in [77, 71] in a computational model of deformable
ellipsoidal particles representing fibroblasts.

In order to illustrate how to model changes of orientation due to volume
exclusion in our framework, we simulate a simplistic model of bacterial swarm
[34, 130] made of active spherocylindrical particles (Fig. 7 and SM Video 16).
The equations of motion only contain a self-propulsion term and the
incompressibility force which models volume exclusion:

$\dot{\bm{x}}_{i}=c_{0}\bm{v}_{i}-\tau\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu}),$

where $c_{0}>0$ is the motion speed and
$\bm{v}_{i}=(\cos\theta_{i},\sin\theta_{i})^{\mathrm{T}}$ is a unit vector in
the direction given by the angle $\theta_{i}$ which defines the cost function
$c(\bm{y},\bm{x}_{i})$ as above. The “collision” between two particles, which
formally means that the associated Laguerre cells share a boundary, induces a
shape change from the base shape which can be quantified by computing the
empirical covariance matrix

$\hat{\Sigma}_{i}=\frac{1}{N_{i}-1}\sum_{\bm{y}\in\mathscr{L}_{i}}(\bm{y}-\bm{x}_{i})(\bm{y}-\bm{x}_{i})^{\mathrm{T}},\quad
N_{i}:=\\#\\{\bm{y}\in\mathscr{L}_{i}\\}.$

We can therefore consider instantaneous plastic deformations, by defining the
covariance of particle $i$ at all time as

$\Sigma_{i}=\frac{1}{\det(\hat{\Sigma}_{i})^{1/d}}\hat{\Sigma}_{i}.$

In dimension 2, this dynamics encodes at the same time the evolution of the
polarity vector $\bm{v}_{i}$ and the aspect ratio. In the numerical experiment
below, particles represent bacteria with a fixed aspect ratio. Thus, only the
direction of motion $\bm{v}_{i}$ evolves, due to the bending effects of the
collisions. Assuming that this directional change is much faster than the
self-propelled motion, $\bm{v}_{i}$ is defined at each time as the normalized
eigenvector associated to the largest eigenvalue of the matrix
$\hat{\Sigma}_{i}$ (with the sign correctly chosen to avoid sign flip).
Alternatively, one could also consider a nematic alignment process, for
instance the relaxation equation

$\dot{\bm{v}}_{i}=\frac{\nu}{2}\nabla_{\bm{v}_{i}}\big{(}(\bm{v}_{i}\cdot\bm{\hat{v}}_{i})^{2}\big{)},$

where $\bm{\hat{v}}_{i}$ is the normalized eigenvector associated to the
largest eigenvlue of $\hat{\Sigma}_{i}$, $\nabla$ denotes the spherical
gradient and $\nu>0$ is the relaxation rate. Both choices lead to long-range
alignment and ordered motion.

(a) $t=0$

(b) $t=3.2$

(c) $t=24$

Figure 7: A swarm of rod-shaped particles with fixed aspect ratio 3 filling
65% of a periodic square domain. The direction of motion is defined at each
time as the largest eigenvector of the empirical covariance matrix of the
Laguerre cells. This collision process produces long-range alignment.
Parameters: grid size $512\times 512$, $c_{0}=0.3$, $\tau=0.25$, $\Delta
t=0.002$. See also SM Video 16.

###  4.2 Surface fluctuations

For some biological cells, the movement is induced by the deformation of their
membrane, namely elongated protrusions, created randomly or in the direction
of a stimulus, that the cell uses to pull itself forward. A fully realistic
model of this complex process would require a detailed description of the cell
internal structure (in particular its actin cortex) and its medium (the extra-
cellular matrix as a fibrous viscoelastic material), like in [57, 140]. Taking
into account all this complexity would certainly be challenging in our model
and we do not claim that our model can achieve a comparable level of details.
Aiming at a minimal phenomenological model, and recalling that the cell
membrane can be modeled as the level-set of the cost function, we model the
local membrane deformations of a base spherical shape using a biased cost
function. Since the cost is a nonnegative multiplicative quantity, we rather
introduce the logarithmic potential

$\varphi(\bm{x},\bm{x}_{i})=-\log c(\bm{x},\bm{x}_{i}).$

This potential $\varphi$ will play a role similar to its homonym in the level-
set method formalism [140]. In this classical setting [105, 104, 122], the
dynamic on $\varphi$ is given by PDE, for instance, like in [140], by

$\partial_{t}\varphi(\bm{x},\bm{x}_{i})+\bm{v}(\bm{x},\bm{x}_{i})\cdot\nabla_{\bm{x}}\varphi(\bm{x},\bm{x}_{i})=0,$

where $\bm{v}$ is a velocity field constructed by incorporating all the
internal and external forces. Although theoretically possible in our setting,
this approach poses technical difficulties linked to the implementation
framework which requires an analytical formula for the cost function (see SM
Appendix 3). Such formula typically cannot be obtained by solving the above
PDE. However, with a somehow simpler approach, we consider cost potentials of
the form

$\varphi=\varphi_{0}+f,$

where $f$ is a (small) perturbation of the potential $\varphi_{0}=-\log c_{0}$
associated to the base $L^{2}$ cost
$c_{0}(\bm{x},\bm{x}_{i})=|\bm{x}-\bm{x}_{i}|^{2}/R_{0}^{2}$ (for simplicity).
This perturbation approach is similar to the model of deformable particles
introduced in [102, 58, 113, 101] where $f$ is constructed from a Fourier
expansion with time-evolving coefficients. In our setting, the volume and non-
overlapping constraints are automatically taken into account by the optimal
transport solver and the incompressibility force, which gives more freedom on
the choice of $f$. We will consider a perturbation $f$ which models the shape
deformations in the two following classical situations.

####  4.2.1 Chemotaxis movement

Many cells, like Dictyostelium, move in response to changes in the
concentration of some chemical quantities in their medium. As explained in
[140] and modeled using a level-set method, the movement is induced by
morphological changes that we aim at reproduce. Similarly to [140], the level-
set deformations are defined to match two known biological shapes.

We consider a chemo-attractant density $u(\bm{x})$ in $\Omega$ and we
classically assume that a particle located at $\bm{x}_{i}$ can sense the local
gradient along the directions $\bm{x}-\bm{x}_{i}$, defined through the finite
difference formula

$\delta
u(\bm{x},\bm{x}_{i})=\frac{u(\bm{x})-u(\bm{x}_{i})}{|\bm{x}-\bm{x}_{i}|}.$

Then we introduce the biased cost potential

$\varphi(\bm{x},\bm{x}_{i})=\varphi_{0}(\bm{x},\bm{x}_{i})+\beta
f\big{(}\delta u(\bm{x},\bm{x}_{i})\big{)},$

for some constant $\beta>0$ and some function $f$ which models how the
gradient affects the deformation. Denoting by $\delta_{+}=\max(0,\delta)$ the
positive part of a number $\delta$, a natural choice for $f$ is

$f(\delta)=-\delta_{+},$ (15)

which corresponds to lowering the potential and thus allowing deformations in
the directions of large nonnegative gradients. This produces elongated shapes
(see Fig 8 and SM Video 17). Another choice is to consider

$f(\delta)=\delta_{-}^{2}+\delta_{+},$ (16)

where $\delta_{-}=\max(0,-\delta)$ is the negative part of the number
$\delta$. This corresponds to a larger cost in the direction of larger
gradient but even larger for negative gradient. This produces fan-like shapes
(see Fig 8 and SM Video 18).

In both situations, we consider the simplest possible dynamics of the
centroids

$\dot{\bm{x}}_{i}=-\tau\nabla_{\bm{x}_{i}}\mathcal{T}_{c_{0}}(\hat{\mu}),$

which, corresponds to movement in the direction of the barycenter (which is
the gradient of the base $L^{2}$ cost, see Eq. (12)). Since the morphological
changes are biased in the direction of increasing gradient, this choice
induces chemotactic migration. An important point to note is that not only the
volume is preserved regardless of the deformation, but this model also
includes non-overlapping and cell-cell repulsion, thanks to the
incompressibility force. As illustrated many times, crowded areas also induce
other morphological changes which would typically be difficult to model in the
level-set formalism (where mostly a single cell or non-interacting cells are
considered).

(a) Elongated shape

(b) Fan-like shape

Figure 8: Chemotaxis motion induced by shape changes with two choices of
perturbation functions. The chemo-attractant density is given by a Gaussian
law with variance $\sigma=0.2^{2}$ and is indicated by the underlying color
field (a) Elongated shape obtained with Eq. (15). (b) Fan-like shape obtained
with Eq. (16). Parameters: $\beta=0.02$, $\tau=1$, grid size $512\times 512$,
$\Delta t=0.01$. See also SM Videos 17,18.

####  4.2.2 Exploration via random protrusions

Mesenchymal cell migration is caused by a complex coordination between
protrusion, adhesion and retraction [123]. Here we assume that cells evolve
randomly without any external stimuli by following randomly created
protrusions. We thus consider a cost perturbation of the form

$f=-\sum_{k=1}^{K}\xi_{k}G_{\theta_{k}}(\theta),$

where $K$ is the number of protrusions, $\theta_{k}$ are the polar angles
defining their location in the membrane and $G_{\theta_{k}}(\theta)$ is the
Gaussian density with center $\theta_{k}$ and a given variance which sets the
width a protrusion. We assume that the elongation range (amplitude) of the
protrusion $\xi_{k}$ is constrained by $\xi_{k}\geq 0$ and follows the
stochastic differential equation

$\mathrm{d}\xi_{k}=-\xi_{k}\mathrm{d}t+\mathrm{d}B^{k}_{t},$

that is, protrusions are created randomly and naturally retract. To model the
cytoskeleton of the cell, we also introduce an elastic force in the direction
of the protrusions. The final equations of motion are thus given by:

$\dot{\bm{x}}_{i}=\int_{\mathscr{L}_{i}}\frac{\bm{x}-\bm{x}_{i}}{|\bm{x}-\bm{x}_{i}|}\,\mathbbm{1}_{|\bm{x}-\bm{x}_{i}|\geq
R_{0}}\,\mathrm{d}\bm{x}-\tau\nabla_{\bm{x}_{i}}\mathcal{T}_{c_{0}}(\hat{\mu}).$

The result is presented in SM Video 19.

##  5 Surface interactions

###  5.1 Goals and previous approaches

When several biological cells are in contact, or in contact with the medium,
various interfacial forces lead to potential complex deformations and
meaningful sorting phenomena. Such phenomena are of course of primary
importance in developmental biology and various computational models have thus
been proposed. The most archetypal results are obtained using cellular Potts,
vertex or finite elements models which are particularly well-suited to
describe the dynamics of tissues and dense aggregates. Due to their
resemblance to soap bubbles, droplets and multi-phase systems, biological cell
aggregates have often been modeled using the physics of capillarity and
wetting phenomena in mind [35, 82, 94]. Over the last decades, the precise
biological mechanisms and their implementation in mathematical models have
been the subject of various debates, some of them are reviewed in [21, 103].
In a series of famous articles [31, 23, 22, 20], Chen and Brodland have
proposed, and verified, the so-called Differential Interfacial Tension
Hypothesis (DITH) which, in short, assumes that the various forces created by
internal pressure, Cell Adhesion Molecules (CAM) etc, can all be resolved into
a difference between equivalent surface tension parameters. This difference
creates tangential forces at the interface between the cells (or the cells and
the medium) which drive their local deformation and their spatial re-
organisation.

In order to validate a computational model, the archetypal system considers an
aggregate of two populations of cells immersed in a medium. Depending on the
surface tension differences, Chen and Broadland have proposed a set of
conditions which should lead to various sorting phenomena observed in biology,
such as homophilic sorting, dissociation, engulfment, checkerboard patterning
[3, 132]. The DITH can be implemented in various computational frameworks. In
their work, Chen and Brodland have advocated the finite element method.
Similarly to vertex models, in a finite element discretization the interfaces
(edges) and junction points (vertices or nodes) are indeed handled directly as
degrees of freedom. This represents a clear advantage to compute local force
balances and interface deformation. However, it may come at the price of
technical difficulties, in particular to deal with topological and
connectivity changes, which are usually handled by-hand with ad hoc rules.
Other powerful approaches which directly model the interface include the
recent [93] which is based on a level-set method and a global energy
minimization scheme of the contact surface. With a different perspective,
point particle systems are intrinsically simpler, both theoretically and
computationally, since mostly only binary attraction-repulsion forces between
points are modeled. Unfortunately, with this somehow too simplified
representation, it is typically impossible or difficult to reproduce the
important sorting phenomena mentioned above [103]. More or less recent
approaches, and in particular the so-called Self-Propelled Voronoi (SPV) and
Active Vertex Models (AVM) [5, 9], are hybrid in the sense that they consider
forces between the centers of a tessellation (thus treated as point particles)
but the forces themselves are computed depending on the properties of a
tessellation which is recomputed at each time step.

Our framework enters this latter hybrid approach where the positions
$\bm{x}_{i}$ are the main degrees of freedom. In [5, 9], the authors consider
Voronoi tessellation and they compute the forces applied to the cell centers
as the negative gradient of some energy, typically built as a function of the
cell perimeters. In comparison, we use a different tessellating procedure,
which ensures a strict volume constraint but which makes the computation of
such gradient much more difficult. Indeed, although theoretically
implementable, computing the gradient of the interface lengths, or other
related quantities, would require to compute the gradient of the Kantorovich
potentials which appear in the definition of the Laguerre cells. Similar
questions also appear in another context in [81] which introduces a meshing
algorithm based on Laguerre tessellations. The differentiability properties of
the solution of an optimal transport problem is a well-known but potentially
involved question [37, 108] that we leave for future work. In the coming
section, we follow a somehow simpler path, which is directly inspired by the
physical approach of [20]. Instead of minimizing an energy, we compute
pressure, surface tension and other interfacial forces that are applied
directly to the centroids $\bm{x}_{i}$. Naturally, the level of precision of
our method cannot be compared to the one of finite element models which can
describe with an arbitrary precision the local deformations of a single cell
[133] thanks to an arbitrary fine mesh. In comparison, in our framework, a
particle is described by a single “macroscopic” mesh element. Nevertheless,
for cell aggregates, where the cells are sufficiently small compared to the
size of the system, we show that our method can accurately reproduce the
complex sorting phenomena obtained with a finite element method in [20, 21],
with a level-set method in [93], with a cellular Potts model in [52, 49, 7]
and with a Voronoi tessellation method [128].

###  5.2 Forces

In our model, a particle is represented by a single Laguerre cell and should
thus rather be considered as a macroscopic immersed body, where the
surrounding environment is made of other cells or the medium. Following basic
physical laws, we thus introduce a pressure-like force acting orthogonally to
the boundary of each Laguerre cell and we assume that the resulting force acts
on the centroid. More precisely, at the interface between two Laguerre cells,
denoted by $\Gamma_{ij}=\mathscr{L}_{i}\cap\mathscr{L}_{j}$, (with possibly
$j=0$ in which case it is an interface between a particle and the medium), we
define the following force which applies on $\bm{x}_{i}$

$\bm{\vec{F}}_{j\rightarrow
i}=-\frac{\gamma_{ij}}{|\bm{x}_{j}-\bm{x}_{i}|}\int_{\Gamma_{ij}}\bm{\vec{n}}\,\mathrm{d}\sigma,\quad\bm{\vec{F}}_{0\rightarrow
i}=-\gamma_{i0}\int_{\Gamma_{i0}}\frac{\bm{\vec{n}}}{|\sigma-\bm{x}_{i}|}\,\mathrm{d}\sigma,$
(17)

where $\vec{\bm{n}}$ is the outward normal to $\Gamma_{ij}$ and
$\gamma_{ij}>0$ is an arbitrary coefficient. In dimension 3, $\gamma_{ij}$ is
a force per unit of length and can thus be directly understood as the surface
tension between the particles $i$ and $j$. In the 2D examples below, we keep
this analogy and we also note that with this modelling choice, the equilibrium
of a so-called cell doublet system is given by the Young-Dupré equation as it
should be [84, 82, 94].

######  Lemma 9.

Let us consider two identical particles at positions $\bm{x}_{i}$ and
$\bm{x}_{j}$ and the $L^{2}$ cost $c(\bm{x},\bm{y})=|\bm{y}-\bm{x}|^{2}$. We
denote respectively by $\gamma_{ij}$ and $\gamma_{0}$ the surface tensions
between the particles $i$ and $j$ and between any of the particles and the
medium. The equilibrium configuration, shown in Fig. 9, is made of two
spherical caps with contact angle $\theta$ which satisfies

$\cos\theta=\frac{\gamma_{ij}}{2\gamma_{0}}.$

###### Proof.

The fact that the particles are spherical caps is a direct consequence of
Theorem 3 (as Lemma 5). We denote by $R$ their radius and $\theta$ the contact
angle. The total force applied on $\bm{x}_{i}$ is thus equal to
$\bm{\vec{F}}_{j\to i}+\bm{\vec{F}}_{0\to i}$ with

$\displaystyle\bm{\vec{F}}_{j\to i}$
$\displaystyle=-\frac{\gamma_{ij}}{|\bm{x}_{j}-\bm{x}_{i}|}\int_{\Gamma_{ij}}\vec{\bm{e}}_{ij}\mathrm{d}\sigma=-\gamma_{ij}\tan\theta\,\bm{\vec{e}}_{ij},$
$\displaystyle\vec{\bm{F}}_{0\to i}$
$\displaystyle=-\frac{\gamma_{0}}{R}\int_{\theta}^{2\pi-\theta}R(\cos\alpha\,\bm{\vec{e}}_{ij}+\sin\alpha\,\bm{\vec{e}}_{ij}^{\perp})\mathrm{d}\alpha=2\gamma_{0}\sin\theta\,\bm{\vec{e}}_{ij}.$

where $\bm{\vec{e}}_{ij}=(\bm{x}_{j}-\bm{x}_{i})/|\bm{x}_{j}-\bm{x}_{i}|$ and
$\bm{\vec{e}}_{ij}^{\perp}$ is the direct orthogonal vector to
$\bm{\vec{e}}_{ij}$. The conclusion follows. ∎

In [20, 22], the authors consider a force acting on each triple junction point
$\bm{z}_{ijk}$ between three distinct particles $i,j,k$ (with possibly $k=0$),
which is given by

$\bm{\vec{\gamma}}_{ijk}=\bm{\vec{\gamma}}_{ij}+\bm{\vec{\gamma}}_{jk}+\bm{\vec{\gamma}}_{ki},$

where $\bm{\vec{\gamma}}_{ij}$ is the vector of magnitude $\gamma_{ij}>0$
oriented outward along the boundary $\Gamma_{ij}$. In the cell-doublet system
of Lemma 9, there are two triple junction points and
$\bm{\vec{\gamma}}_{ijk}=0$ when the contact angle satisfies Young-Dupré
equation. As explained above, the triple junction points are computed
automatically so we need to report these forces on the centroids. We simply
assume a linear spring force between the centroids and the junction points,
which can be interpreted as a very basic cytoskeleton model. This leads to the
force

$\bm{\vec{F}}_{\bm{z}_{ijk}\to
i}=\left(\bm{\vec{\gamma}}_{ij}\cdot\frac{\bm{z}-\bm{x}_{i}}{|\bm{z}-\bm{x}_{i}|}\right)\frac{\bm{z}-\bm{x}_{i}}{|\bm{z}-\bm{x}_{i}|},$
(18)

where $\cdot$ denotes the Euclidean dot product.

Figure 9: Geometry of a cell doublet

###  5.3 Cell sorting experiments

We consider the dynamical system given by

$\dot{\bm{x}}_{i}=\alpha\sum_{j=0}^{N}\bm{\vec{F}}_{j\to
i}+\beta\sum_{j=1}^{N}\sum_{k=0}^{N}\bm{\vec{F}}_{\bm{z}_{ijk}\to
i}-\tau\nabla_{\bm{x}_{i}}\mathcal{T}(\hat{\mu}),$ (19)

where the force in the first sum is equal to zero if the Laguerre cells $i$
and $j$ do not share a boundary and similarly in the second sum if the three
Laguerre cells do not share a triple junction point. Moreover, although
theoretically not necessary in this situation, the incompressibility force is
added in the experiments below to ensure that the configuration remains
sufficiently realistic and in particular that $\bm{x}_{i}\in\mathscr{L}_{i}$.

Following a classical experimental set up introduced in [23] and reviewed in
[103, 21, 93], we consider an initial mixed aggregate with two types of cells
(blue and orange) and we vary the five surface tension parameters, namely
$\gamma_{b}$, $\gamma_{o}$, $\gamma_{bo}$, $\gamma_{bw}$ and $\gamma_{ow}$,
respectively between blue particles, orange particles, blue and orange
particles, blue particles and the medium and orange particles and the medium.
The simulation results are presented in Fig. 10 and SM Videos 20,21,22, where
our model correctly captures the three following important configurations.

  * •

Homophilic sorting with internalization of one of the population.

  * •

Homophilic sorting with separation by the medium.

  * •

Checkerboard patterning.

Figure 10: Sorting patterns in a mixture of two cell types, obtained from the
same initial condition (top) under different conditions on the surface
tensions. (Left) Internalization when $\gamma_{bo}<\gamma_{bw}-\gamma_{ow}$
($\gamma_{b}=\gamma_{o}=1$, $\gamma_{bo}=5$, $\gamma_{ow}=30$,
$\gamma_{bw}=7$). (Middle) Checkerboard when
$\gamma_{bo}<\frac{1}{2}\gamma_{b}$ ($\gamma_{b}=\gamma_{o}=7$,
$\gamma_{bo}=1$, $\gamma_{bw}=\gamma_{ow}=30$). (Right) Separation when
$\gamma_{bo}>\gamma_{bw}+\gamma_{ow}$ ($\gamma_{b}=\gamma_{o}=1$,
$\gamma_{bo}=12$, $\gamma_{bw}=\gamma_{ow}=5$). Parameters: $\tau=0.5$,
$\alpha=\beta=0.14$, grid size $600\times 600$, $\Delta t=0.0003$. See also SM
Videos 20,21,22.

In addition, Fig. 11 and SM Video 23 show an engulfment phenomenon where the
two populations are initially separated with a small contact interface before
one population becomes internalized.

(a) $t=0$

(b) $t=0.033$

(c) $t=0.23$

(d) $t=1.23$

Figure 11: Engulfment of the orange cells inside the mass of blue cells.
Parameters: $\tau=0.5$, $\alpha=\beta=0.14$, $\gamma_{b}=\gamma_{o}=1$,
$\gamma_{bo}=5$, $\gamma_{ow}=30$, $\gamma_{bw}=7$, grid size $512\times 512$,
$\Delta t=0.0005$. See also SM Video 23.

##  6 Discussion and conclusion

In this article we have introduced a novel approach to model, within a same
framework, diverse biological systems ranging from individual deformable cells
to complex tissue aggregates. The core idea is based on semi-discrete optimal
transport theory and the notion of Laguerre tessellation, which together
provide a natural way to represent arbitrary shapes and physical deformations
with a strict volume constraint. We have verified the validity of our method
by reproducing a wide range of classical models in computational biology.
Although our method is formally an extension of the idea of Voronoi
tessellation, we have highlighted formal links with other independent methods,
in particular point-particle systems, level-set methods and vertex models. We
also provide an efficient GPU implementation and a detailed documentation.

Since, to the best of our knowledge, our approach is new in the literature,
the main focus of the present article was to validate it by comparison to
other classical computational biology models, thus providing an independent
validation of their biological conclusions. The next important step is to seek
a direct biological validation of our method by the confrontation to
experimental data. In particular, our method seems well suited to the in
silico modeling of the development of blastocysts. Recent experimental
findings [63, 92] have identified several complex phenomena with morphological
changes, sorting phenomena and dynamical topological properties via the
formation of lumen. In order to understand the regulation processes between
these phenomena, we need a mathematical model that can handle both individual
cell deformations and their behaviours among dense aggregates. This seem
particularly challenging using traditional methods, or at the price of an
unreasonable computational cost, especially in a three-dimensional setting
(note that although the present paper mostly focuses on 2D situations for
simplicity, the theory readily extends to 3D). This biological question was
the initial motivation for the development of our method and is still an on-
going work.

On the theoretical and numerical sides, several open questions and future
improvements could be considered:

  * •

Since our model has an intrinsic point-particle description, a natural
question is to extend the active vertex theory developed in [5, 9, 128] to the
case of Laguerre tessellations with volume constraint. This lead to
theoretical questions on the differentiability properties of the solution of
an optimal transport problem.

  * •

Another clear advantage of a point-particle description is the potential to
derive coarse-grained PDE models that are better suited to mathematical
analysis. This seems to be a largely open question, with the notable exception
of [114, 98] on the notion of interacting Laguerre cells as a particle method
to discretize particular cases of gradient flow PDEs. Our formal analysis
(gathered in the SM Appendix 2) rather suggests a methodology based on mean-
field theory and should lead to the study of novel Fokker-Planck-Monge-Ampère
coupled systems of PDEs.

  * •

The theory of semi-discrete optimal transport theory is historically linked to
the design of numerical schemes for fluid flows. Thus, it provides a direct
way to include fluid-structure interactions between macroscopic cells immersed
in a fluid environment. This corresponds to a more complex setting than the
one considered in the present article, but which is typically considered in
the so-called immersed boundary methods [107, 66, 112].

  * •

Although our GPU implementation is general and reasonably efficient, several
performance improvements are still required for real-time large-scale
simulations, especially in 3D. First of all, as detailed in SM Appendix 3
(Section C), we must adapt annealing schedules that were designed for the
squared Euclidean cost to our much more general setting. Second, implementing
the block-sparse truncation method that corresponds to Algorithm 3.6 of [42]
should result in a significant performance boost. Finally, complex dynamical
costs pose a bigger challenge, especially if they are defined implicitly as
the solution of a PDE as in level-set methods [140].

##  7 Supplementary Information

######  SM Appendix 1 (Section A).

This section gathers additional numerical experiments which complement the
conclusions in the main text and showcase other applications. We consider a
simple but classical model of tissue growth, show how to handle motion in
arbitrary domains, confirm an important glass transition result for active
Brownian particles and reproduce a crowd simulation leading to congestion
effects.

######  SM Appendix 2 (Section B).

This rather mathematical section formally explores the links between our
agent-based model and coarse-grained PDE descriptions. We draw a parallel with
related known results on crowd motion modeling and introduce a new general
class of coupled Fokker-Planck-Monge-Ampère systems.

######  SM Appendix 3 (Section C).

This section gathers reminders and classical results in measure theory and
optimal transport.

The following supplementary videos can be found on Figshare at

https://doi.org/10.6084/m9.figshare.25240669.v1

######  SM Video 1 (Incompressibility).

Related to Fig. 3: simulation of a simple gradient descent dynamics in a four-
cell system, where the negative gradient of the cost function (called
incompressibility force) is displayed in red.

######  SM Video 2 (Run-and-tumble in a disk).

Related to SM Appendix 1: simulation of a system of self-propelled particles
with a constant speed motion, random direction changes and random energies in
a disk domain.

######  SM Video 3 (Tissue growth).

Related to SM Appendix 1: simulation of a basic tissue growth model based on
division and cell growth starting from a single cell.

######  SM Video 4 (Falling spheres (1/6)).

Related to Fig. 4: falling spheres with the deformation parameter
$\alpha=0.75$. Initially round shape cells adopt a columnar shape due to the
external force and a high tolerance to deformation.

######  SM Video 5 (Falling spheres (2/6)).

Related to Fig. 4: falling spheres with the deformation parameter $\alpha=2$.
The initially round particles are not so prone to deformation: they thus keep
a convex shape and organize into a Voronoi-like configuration.

######  SM Video 6 (Falling spheres (3/6)).

Related to Fig. 4: falling spheres with the deformation parameter $\alpha=10$.
The particles behave as hard-spheres.

######  SM Video 7 (Falling spheres (4/6)).

Related to Fig. 5: falling spheres where the magnitude of the
incompressibility force of the orange and blue spheres are respectively
$\tau_{o}=8$ and $\tau_{b}=3$.

######  SM Video 8 (Falling spheres (5/6)).

Related to Fig. 5: falling spheres where the magnitude of the
incompressibility force of the orange and blue spheres are respectively
$\tau_{o}=3$ and $\tau_{b}=3$.

######  SM Video 9 (Falling spheres (6/6)).

Related to Fig. 5: falling spheres where the magnitude of the
incompressibility force of the orange and blue spheres are respectively
$\tau_{o}=1$ and $\tau_{b}=3$.

######  SM Video 10 (Micropipette aspiration).

Related to Fig. 6: micropipette experiment for the deformability parameters
$\alpha=0.5$, $\alpha=0.75$, $\alpha=1$, $\alpha=2$. The aspiration length
decreases with the deformability.

######  SM Video 11 (Active Brownian Particles (1/3)).

Related to SM Appendix 1: simulation of a system of active Brownian particles
with deformability parameter $\alpha=0.5$. The mean-square displacement of the
particles is smaller but of the same order of magnitude than the theoretical
value for independent point particles.

######  SM Video 12 (Active Brownian Particles (2/3)).

Related to SM Appendix 1: simulation of a system of active Brownian particles
with deformability parameter $\alpha=2$. As the deformability decreases, the
mean-square displacement is drastically reduced.

######  SM Video 13 (Active Brownian Particles (3/3)).

elated to SM Appendix 1: simulation of a system of active Brownian particles
with deformability parameter $\alpha=10$. The particles behave as hard-spheres
and self-organize into an optimal sphere packing configuration with almost no
movement.

######  SM Video 14 (Crowd motion (1/2)).

Related to SM Appendix 1: simulation of crowd motion towards a single exit
point with deformability parameter $\alpha=0.5$. All the particles manage to
escape at the price of large deformations.

######  SM Video 15 (Crowd motion (2/2)).

Related to SM Appendix 1: simulation of crowd motion towards a single exit
point with deformability parameter $\alpha=8$. The particles behave as hard-
spheres and end up in a stable congested arch state.

######  SM Video 16 (Rod-shape particles).

Related to Fig. 7: simulation of a system of rod-like particles with
directional changes induced by deformation and bending effects. Long-range
order spontaneously emerges after some time.

######  SM Video 17 (Chemotaxis motion (1/2)).

Related to Fig. 8: simulation of a system of 10 particles where the movement
is induced by the shape deformation in response to a chemo-attractant
gradient. With a particular choice of cost perturbation, particle migrate by
adopting elongated shapes in the direction of increasing chemo-attractant
gradient.

######  SM Video 18 (Chemotaxis motion (2/2)).

Related to Fig. 8: simulation of a system of 10 particles where the movement
is induced by the shape deformation in response to a chemo-attractant
gradient. With a particular choice of cost perturbation, particle migrate by
adopting fan-like shapes orthogonal to the direction of increasing chemo-
attractant gradient.

######  SM Video 19 (Random protrusions).

Simulation of a random motion induced by random protrusions.

######  SM Video 20 (Cell sorting (1/3)).

Related to Fig. 10: simulation of a two-population cell aggregate under the
differential interfacial tension hypothesis. The orange cells are internalized
inside the mass of blue cells.

######  SM Video 21 (Cell sorting (2/3)).

Related to Fig. 10: simulation of a two-population cell aggregate under the
differential interfacial tension hypothesis. The orange and blue cells self-
organize into a checkerboard pattern state.

######  SM Video 22 (Cell sorting (3/3)).

Related to Fig. 10: simulation of a two-population cell aggregate under the
differential interfacial tension hypothesis. The orange and blue cells cluster
with their own type and are separated by the medium.

######  SM Video 23 (Engulfment).

Related to Fig. 11: simulation of a two-population cell aggregate under the
differential interfacial tension hypothesis. the orange cells are initially
completely outside the aggregate of blue cells but a small initial interface
lead to the total engulfment of the orange cells inside the mass of blue
cells.

######  SM Video 24 (Corner in 3D).

3D example where all the particles are attracted to a corner with a constant
force. The blue particles are less deformable than the red particles. The
video on the left-hand side shows the raw data computed on a grid of size
$100\times 100\times 100$. In the video on the right-hand side the surface of
each cell is smoothened: the support of each cell is first binarized and
smoothened by a Gaussian of standard deviation equal to 1.5 pixel. The smooth
surface is then defined as the level-set of value 0.5.

######  SM Video 25 (Run-and-tumble in 3D).

This is the same run-and-tumble system as in SM Video 2 but in a 3D cube. The
darker particles are more deformable than the lighter ones. The smoothening
procedure is the same as in SM Video 24.

## Acknowledgments

This paper was inspired by discussions with Takafumi Ichikawa on the modelling
of embryonic development and with David Bourne, Maciej Buze and Steven Roper
on the modelling of polycrystalline metals. The authors wish to thank Bruno
Lévy, Roya Mohayaee, Steffen Plunder, Tsubasa Sukekawa, Tetsuya Hiraiwa and
Sungrim Seirin-Lee for their insightful advice and for the useful discussions.
The work of AD is supported by the following grant: KAKENHI Grant-in-Aid for
Early-Career Scientists (Grant number 23K13015).

Supplementary material

##  Appendix A Additional experiments

###  A.1 Tissue growth

We consider a simplistic model of tissue growth due to cell division and
volume expansion. This is a classical test case for computational models of
cell aggregates [99, 41, 5, 69].

We assume that the volume of the particles increases linearly

$\dot{v}_{i}=\lambda_{i}\gamma$

where $\lambda_{i}$ is a uniform random variable between 0.5 and 2,
independent for all the particles. The growth continues until the volume
reaches the maximum $v_{\mathrm{max}}$. All the particles with a volume above
80% of $v_{\mathrm{max}}$ divide in two along a random axis at a rate $\nu$.
In addition to the incompressibility force, all the particles are subject to
an external cohesive force of magnitude $F_{0}$ oriented towards the center.
The result is shown in Fig 12 and SM Video 3.

(a) $t=0$

(b) $t=1.3$

(c) $t=3.3$

(d) $t=4.42$

Figure 12: Growing tissue from 1 to 400 cells. Parameters: $\gamma=1.8\cdot
10^{-3}$, $\nu=5$, $\tau=1$, $F_{0}=0.2$, grid size $512\times 512$, $\Delta
t=0.01$. See also SM Video 3.

###  A.2 Run-and-tumble in a disk with energy differences

In this system, the particles are constrained to move inside a disk. They all
have the same volume $v_{0}$ but unlike most of the other numerical
experiments so far, the deformations are not handled by changing the cost
function but rather by changing its scaling. We consider the $L^{2}$ cost

$c(\bm{x},\bm{x}_{i})=\frac{2\lambda_{i}}{R_{0}^{2}}|\bm{x}-\bm{x}_{i}|^{2},$

where the scaling factor $\lambda_{i}>0$ is sampled independently for each
particles, so that

$\int_{B(\bm{x}_{i},R_{0})}c(\bm{x},\bm{x}_{i})\mathrm{d}\bm{x}=\lambda_{i}v_{0}.$

Consequently, the boundary between two particles is always a circular arc as
advocated in [10] under the name “quotient method”. However, particles
associated to a lower $\lambda_{i}$ have a smaller base cost and are thus
prone to deformations.

The particle motion is given by a simple self-propelled motion

$\dot{\bm{x}_{i}}=c_{0}\bm{v}_{i}-\tau\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu}),$

where $c_{0}>0$ and the direction $\bm{v}_{i}$ is a unit vector which is re-
sampled uniformly at rate $\nu$ (run-and-tumble). The result is shown in Fig.
13 and SM Video 2.

Figure 13: Particles filling 87.5% of a disk domain of radius 0.5 and
different energies given by independent random variables $\lambda_{i}$ sampled
uniformly in the interval $(0.4,3.5)$. Lighter color indicates lower energy
(higher deformability). The red arrow indicates the direction of motion of a
sample of five particles. Parameters: $c_{0}=0.4$, $\tau=0.2$, $\nu=1$, grid
size $512\times 512$, $\Delta t=0.0025$. See also SM Video 2.

The same system in a three-dimensional setting is shown in SM Video 25.

###  A.3 Deformable active Brownian particles and glass transition

In [113, 9], a system of deformable active Brownian particles is considered
and it is shown that when only little deformation is allowed (or when the
self-propulsion force is too small), then the system end up in a solid state
characterised by the optimal hexagonal packing of spherical particles. In our
framework, the particle motion is given by the following system of equations:

$\displaystyle\mathrm{d}\bm{x}_{i}$
$\displaystyle=c_{0}\bm{v}_{i}\mathrm{d}t-\tau\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu})\mathrm{d}t$
$\displaystyle\mathrm{d}\bm{v}_{i}$
$\displaystyle=(\mathrm{Id}-\bm{v}_{i}\bm{v}_{i}^{\mathrm{T}})\circ\sqrt{2\sigma}\mathrm{d}B^{i}_{t},$

where $\bm{v}_{i}$ is a unit vector which defines the direction of motion and
$B^{i}_{t}$ is an independent Brownian motion in $\mathbb{R}^{2}$. The
projection operator $\mathrm{Id}-\bm{v}_{i}\bm{v}_{i}^{\mathrm{T}}$ ensures
that $|\bm{v}_{i}|=1$ at all time (this corresponds to Brownian motion on the
sphere) and this stochastic differential equation should be understood in the
Stratonovich sense. As in the main text, in order to model deformability, the
cost function is of the form

$c(\bm{y},\bm{x})=\lambda^{-1}_{\alpha}|\bm{y}-\bm{x}|^{\alpha},$

with the same normalization constant $\lambda_{\alpha}$ as in the main text.
The motion of the particles is measured by the Mean Square Displacement (MSD).
For point particles, it is well known that the MSD grows linearly in time with
a known theoretical rate. For soft spheres particles, the MSD seems to still
grow linearly in time but with reduced rate (which seems to decreases
geometrically with $\alpha$). For $\alpha=8$, the motion is almost impossible
and we do observe a stable optimal sphere packing state as in [113]. See Fig
14 and SM Videos 11,12,13.

(a) $\alpha=0.5$

(b) $\alpha=2$

(c) $\alpha=8$

(d) MSD

Figure 14: (a)-(c) Configuration at $t=5$ of a system of 250 soft spheres
filling 94% of a square periodic domain, for three values of the deformation
parameter $\alpha$. (d) Mean-square displacement as a function of $\alpha$ and
linear fitting. Parameters: $c_{0}=0.3$, $\sigma=20$, $\tau=3$, grid size
$512\times 512$, $\Delta t=0.002$. See also SM Videos 11,12,13.

A natural question in this framework is to prove theoretically that the MSD
grows linearly in time regardless of $\alpha$ (or not) and if it is the case,
to compute the effective diffusion rate. For point particles, such computation
can be carried out easily, for instance by studying the Fokker-Planck equation
associated to the stochastic differential equation defining the motion of an
active Brownian particle. In the present system, the non-overlapping
interactions make this analysis much more difficult because there is no single
Fokker-Planck equation describing the motion of a typical particle. Such
equation could formally be obtained in the limit $N\to+\infty$ as discussed in
SM Appendix 2. However, due to the nonlinear structure of this equation, the
computation of the MSD is still an open question for us.

###  A.4 Crowd motion and congestion

In this numerical experiment, particles represent pedestrians in a crowd who
all wish to escape a room via a single exit point. The motion is given by the
same system of equations as in the case of falling spheres except that the
force is directed towards the unique exit point. Some noise is also added to
the direction of motion.

$\displaystyle\mathrm{d}\bm{x}_{i}$
$\displaystyle=c_{0}\bm{v}_{i}\mathrm{d}t-\tau\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu})\mathrm{d}t$
$\displaystyle\mathrm{d}\bm{v}_{i}$
$\displaystyle=(\mathrm{Id}-\bm{v}_{i}\bm{v}_{i}^{\mathrm{T}})\circ(\bm{\vec{b}}_{0}(\bm{x}_{i})\mathrm{d}t+\sqrt{2\sigma}\mathrm{d}B^{i}_{t}),$

where $\bm{\vec{b}}_{0}(\bm{x}_{i})$ is a unit vector pointing towards the
exit of the room and the cost belongs to the family of power costs as in the
previous system of active Brownian particles. When the deformation parameter
$\alpha$ is small, all the particles manage to exit the room at the price of
strong deformations (SM Video 14). When $\alpha$ is too big, particles behave
like hard-spheres and end up in a congested state, forming a stable arch
around the exit point (Fig 15 and SM Videos 15).

(a) $t=0$

(b) $t=0.5$

(c) $t=5$

Figure 15: Simulation of a an evacuation situation where each particle moves
towards the exit (red dot at the center right part of the domain) up to some
noise. There are initially 111 particles filling 42% of the rooom. The last
arch configuration at $t=5$ is stable and no particle manage to exit.
Parameters: $\alpha=8$, $c_{0}=0.4$ $\sigma=0.2$, $\tau=3$, grid size
$512\times 512$, $\Delta t=0.001$. See also SM Videos 15,14.

The mathematical modelling of crowd motion and congestion effect is
notoriously difficult [83]. Usually, particles, representing pedestrians, are
modelled by spheres with a tunable overlapping parameter. Our framework is
slightly different as we rather consider physical deformations with a volume
constraint (which may appear more natural). Both approaches seem to agree
asymptotically in the sense that when no overlapping or no deformation is
allowed, then the system is composed of hard-spheres and congestion effects
such as stable arches appears. An important but difficult mathematical
question is to model these phenomena at the continuum PDE level, obtained as
the scaling limit of an infinite number of particles. This question is
discussed in SM Appendix 2 in relation to the previous works [83, 75, 115].

##  Appendix B Formal mean-field limit: the Fokker-Planck-Monge-Ampère system

There is a fundamental duality between agent-based and continuum PDE models in
mathematical biology. Each description has its own strengths and weaknesses.
Typically, agent-based models provide the finest level of details and modeling
freedom and are naturally fit for in silico computations with potentially
stochastic components. However, they provide little room for theoretical
mathematical analysis due to their complexity and usually high-dimensional
nature. On the contrary, continuum PDE models provide a synoptic point of view
at a statistical scale, benefit from all the analytical mathematical
machinery, are often closer to physical laws and can eventually lead to much
stronger “proved” conclusions. Bridging the gap between the two descriptions
in a rigorous mathematical manner is a fundamental issue in mathematical
biology but its feasibility strongly depends on the modeling framework [24].
Point-particle systems are classically the easier to coarse-grain, as they can
often be written in a statistical physics framework for which many theoretical
results are now available [27, 28]. Particles with a shape, regardless of the
description (vertex models, cellular automata, phase-fields etc), seem much
more difficult to handle, since no universal or natural continuum limit and
scaling can be easily deduced. As such, the framework introduced in the
present article can be seen as an intrinsically agent-based model. However
formal computations (at this stage) and earlier works in computational optimal
transport suggest that the underlying point-particle description can actually
naturally lead to coarse-grained continuum models.

Let us recall the general form of the equations of motion of a particle system
under the incompressibility force, external, possibly interacting, forces and
independent Brownian noise:

$\mathrm{d}\bm{x}_{i}=\bm{b}(\bm{x}_{i},\hat{\mu})\mathrm{d}t-\tau\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu})\mathrm{d}t+\sqrt{2\sigma_{i}}\mathrm{d}\bm{B}^{i}_{t}.$
(20)

Without the incompressibility force, this system is known as a McKean-Vlasov
system and has been studied in great details since the seminal work of McKean
[85] following the program introduced by Kac [68]. Review articles with a
discussion of the numerous applications include [131, 96, 65, 27, 28]. The
main idea, which originates from statistical physics, is to consider the limit
of the sequence of random empirical measures
$\hat{\mu}\equiv\hat{\mu}_{N}:=\frac{1}{N}\sum_{i=1}^{N}\delta_{\bm{x}_{i}}$
when $N\to+\infty$. When $\tau=0$ and under some assumption on $\bm{b}$, it
can be shown that this sequence has a deterministic limit $\hat{\mu}_{N}\to f$
which is the solution of the (nonlinear) Fokker-Planck equation

$\partial_{t}f=-\nabla_{\bm{x}}\cdot(\bm{b}(x,f)f)+\Delta_{\bm{x}}f.$

This type of limit, called mean-field limit, is deeply linked to the notion of
molecular chaos which originates from the work of Boltzmann and which has been
formalised by Kac [68, 27, 28]. Within this framework, the solution $f$ is the
probability distribution of the state of a single typical particle and thus
models the system at a statistical scale. Many important mathematical results,
which are relevant in biology, have been obtained with this method, which
could not all be cited here, see for instance [26] as an example for cell
sorting and [28] for a review of other applications, in particular in
collective dynamics.

Coming back to the Eq. (20), a first difficulty comes from the fact that the
incompressibility force $\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu})$ does
not appear as a simple function of $\hat{\mu}$ due to the gradient in
$\bm{x}_{i}$. To overcome this issue, we recall the notion of first variation
of a functional $\mathcal{F}:\mathcal{P}(\Omega)\to\mathbb{R}$ on the space of
probability measure $\mathcal{P}(\Omega)$ : this is the function
$\frac{\delta\mathcal{F}}{\delta\mu}(\mu):\Omega\to\mathbb{R}$ defined such
that

$\frac{\mathrm{d}}{\mathrm{d}h}\mathcal{F}(\mu+h\chi)\Big{|}_{h=0}=\int_{\Omega}\frac{\delta\mathcal{F}}{\delta\mu}(\mu)\mathrm{d}\chi,$

for every measure $\chi$ such that $\mu+h\chi\in\mathcal{P}(\Omega)$ for small
enough $h$. With this notion, it is possible to apply the chain rule to
compute the gradient

$\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu})=N^{-1}\nabla\frac{\delta\mathcal{T}_{c}}{\delta\mu}(\hat{\mu})(\bm{x}_{i}),$

where the gradient on the right-hand side is the usual gradient of a real-
valued function on $\mathbb{R}^{d}$. Thus taking the scale
$\tau\equiv\tau_{N}=N\tau_{0}$, Eq. (20) can be rewritten,

$\mathrm{d}\bm{x}_{i}=\bm{b}(\bm{x}_{i},\hat{\mu})\mathrm{d}t-\tau_{0}\nabla\frac{\delta\mathcal{T}_{c}}{\delta\mu}(\hat{\mu})(\bm{x}_{i})\mathrm{d}t+\sqrt{2\sigma_{i}}\mathrm{d}\bm{B}^{i}_{t}.$
(21)

This latter equation enters the classical framework and lead to the formal
mean-field limit

$\partial_{t}f=-\nabla\cdot(\bm{b}(\bm{x},f)f)+\tau_{0}\nabla\cdot\left(\nabla\frac{\delta\mathcal{T}_{c}}{\delta\mu}(f)f\right)+\Delta
f.$ (22)

At this stage, it remains to identify the first variation of the transport
cost $\mathcal{T}_{c}$. We first note that if $f$ satisfies Eq. (22), it
should have a density with respect to the Lebesgue measure. We also assume
that the cost is of the form $c(\bm{x},\bm{y})=\ell(\bm{y}-\bm{x})$ for a
convex function $\ell$. Since the work of Brenier, McCann, Gbangbo, Cafferelli
and others [17, 47, 25, 38], it is well-known that the Monge problem
$T\\#f=\mathrm{Leb}$ has a unique solution which is of the form

$T(\bm{x})=\bm{x}-(\nabla\ell)^{-1}(\nabla\Phi(\bm{x})),$

where $\Phi$ is a Kantorovich potential and it satisfies the Monge-Ampère
equation [38]

$|\det(\nabla T)|=f.$

Then it can be shown [114, Proposition 7.17] that

$\frac{\delta\mathcal{T}_{c}}{\delta\mu}(f)=\Phi.$

Consequently, and after a few more computations [47], Eq. (22) can be
rewritten as a coupled system of Fokker-Planck-Monge-Ampère equations

$\displaystyle\partial_{t}f=-\nabla\cdot(\bm{b}(\bm{x},f)f)+\tau_{0}\nabla\cdot\left(\nabla\Phi
f\right)+\Delta f.$ (23a)
$\displaystyle\det\big{(}\mathrm{I}-\nabla^{2}\ell^{*}(\nabla\Phi)\nabla^{2}\Phi\big{)}=\pm
f,$ (23b)

where $\ell^{*}$ denotes the Legendre transform of $\ell$ defined by

$\ell^{*}(\bm{y})=\sup_{\bm{x}\in\mathbb{R}^{d}}\\{\bm{x}\cdot\bm{y}-\ell(\bm{x})\\}.$

For the $L^{2}$ cost $c(\bm{x},\bm{y})=\frac{1}{2}|\bm{x}-\bm{y}|^{2}$, it
holds that $\ell^{*}=\ell$ and thus $\nabla^{2}\ell^{*}$ is the identity and
Eq. (23b) reduces to

$\det(\mathrm{I}-\nabla^{2}\Phi)=\pm f.$

In this case, this system can also be derived less abstractly by using the
direct expression, shown in [139],

$\dot{\bm{x}}_{i}=-\tau\nabla_{\bm{x}_{i}}\mathcal{T}_{c}(\hat{\mu})=\tau\int_{L_{i}}(\bm{x}-\bm{x}_{i})\mathrm{d}\bm{x}=\frac{\tau}{N}(\bm{\bar{x}}_{i}-\bm{x}_{i})\,$
(24)

where $\bm{\bar{x}}_{i}=N\int_{L_{i}}\bm{x}\,\mathrm{d}\bm{x}$ is the
barycenter of the Laguerre cell. In this semi-discrete setting, there is no
Monge map sending $\hat{\mu}$ to the Lebesgue measure. However, it is natural
to think that as $N\to+\infty$, the Laguerre shrinks to a single point and if
$\hat{\mu}\to f$, then the barycentric map $\bm{x}_{i}\mapsto\bm{\bar{x}}_{i}$
should behave like the map $T$ which sends the continuous measure $f$ to the
Lebesgue measure.

Actually, the same particle system Eq. (24) has been studied in [75] with the
important difference that the Laguerre cell is associated to the more involved
semi-discrete optimal transport problem where the Lebesgue measure is replaced
by the projection of $\hat{\mu}$ on the set of probability density functions
$\rho$ which satisfies $\rho\leq 1$. In this case, the limit
$\hat{\mu}\to\rho$ is proved (for a particular scaling) where $\rho$ is a weak
solution of the system

$\displaystyle\partial_{t}\rho+\nabla\cdot(\rho v)=0,$ $\displaystyle
v=-\nabla p,$ $\displaystyle p\geq 0,\,p(1-\rho)=0.$

This latter system is introduced as a model of crowd motion with the
congestion constraint $\rho\leq 1$. We also refer to [115] and to [83] for an
in-depth discussion of this subject. See also [98] for a similar approach
leading to porous media equations. Note also that in these works, the particle
system is not introduced per se but rather as a particle discretization scheme
of a PDE model, following the fruitful gradient flow theory developed by Otto,
Jordan and Kinderlehrer [114, Chapter 8]

Our approach is somehow simpler, as we only consider the Lebesgue measure as
the source measure of the semi-discrete optimal transport model. However, our
starting point is the particle system and for arbitrary costs, the more
abstract formulation Eqs. 21-(23) seems more appropriate. In the numerical
experiments on crowd motion in SM Appendix 1, the typical congestion effect,
materialized by the formation of stable arches, seems to also appear by taking
the hard-sphere limit ($\alpha\to+\infty$) starting from deformable soft
spheres. It would be interesting to study if such limit can be derived in the
Fokker-Planck-Monge-Ampère system Eq. (23). This system is not so common in
literature, with the exception of [15] where a special case is studied in
connection with the Euler equations. This idea thus seems to open many
interesting, rather theoretical questions, in particular regarding the
following topics.

  * •

The rigorous mean-field limit $\hat{\mu}\to f$, in the presence of stochastic
terms: to the best of our knowledge, it seems to be an open question, although
it may be linked to optimal control systems (see [115]).

  * •

Are there other scaling limits, in particular when the cost depends on $N$, as
in the finite energy scaling discussed in the main text? Related questions for
the solutions of Eq. (23) are also interesting: in particular when
$\nabla\Phi=o(1)$ in which case the Monge-Ampère equation reduces to a simpler
elliptic equation and the full system Eq. (23) to a parabolic-elliptic system
which can be seen as a Keller-Segel system with dispersion instead of
concentration.

  * •

What happens for more complex transport costs? In particular, in the numerical
experiments, the transport cost is rather of the form

$\hat{\mu},\hat{\mu}_{0}\mapsto\mathcal{T}_{c}(V\hat{\mu}+(1-V)\hat{\mu}_{0}),$

where $V<1$ is the volume occupied by the particle and $\hat{\mu}_{0}$ is a
discretization of the free space (either with one particle with zero cost or
many small particles).

##  Appendix C Optimal transport reminders

We gather here some important definitions and results that may be useful for a
better understanding of the main text. Optimal transport theory is by now
well-established and we refer to the classical references [114, 137, 138] for
a deeper introduction to the subject from a mathematical perspective and to
[108, 42] for a focus on numerical methods and applications in data sciences
and computational geometry.

Historically, Monge considered the problem of transporting sand piles located
at some positions $\bm{y}_{1},\ldots,\bm{y}_{M}\in\Omega\subset\mathbb{R}^{d}$
and with respective masses $\beta_{1},\ldots,\beta_{M}$, to some target
positions $\bm{x}_{1},\ldots,\bm{x}_{N}\in\Omega$ that should respectively
receive a volume $\alpha_{1},\ldots,\alpha_{N}$ of sand. The total volume
$\sum_{j}\beta_{j}=\sum_{i}\alpha_{i}=1$ is classically normalized to 1. The
transport of a sand grain between the locations $\bm{y}$ and $\bm{x}$ comes
with a cost $c(\bm{y},\bm{x})$; the optimal transport problem is to find an
optimal way of moving all of the sand at a minimal cost. When $N=M$ and
$\alpha_{i}=\beta_{j}=1/N$, this problem can be seen as a matching problem,
which consists in finding the one-to-one map (or equivalently a permutation of
the set $\\{1,\ldots,N\\}$)

$T:\\{\bm{y}_{1},\ldots,\bm{y}_{N}\\}\to\\{\bm{x}_{1},\ldots,\bm{x}_{N}\\},$

that minimizes the total transport cost

$\mathcal{T}_{c}:=\frac{1}{N}\sum_{i=1}^{N}c(\bm{y}_{i},T(\bm{y}_{i})).$

In the general case with $N\neq M$ and arbitrary volumes, such a Monge map may
not exist since we may have to split the sand piles and send their different
parts to different locations. The optimal transport problem then only makes
sense for the more general notion of transport plan introduced by Kantorovich.
In this setting, a transport plan is a $N\times M$ matrix $\pi=(\pi_{ij})$
whose non-negative coefficients $\pi_{ij}\geq 0$ indicate how much of the sand
pile $j$ is send to the location $i$. Consequently, the transport plan $\pi$
should satisfy the constraints

$\forall j,\,\,\sum_{i=1}^{N}\pi_{ij}=\beta_{j},\quad\forall
i,\,\,\sum_{j=1}^{M}\pi_{ij}=\alpha_{i}.$

The optimal transport problem then consists in finding the matrix $\pi$ that
satisfies this constraints and minimizes the total transport cost

$\mathcal{T}_{c}:=\sum_{i=1}^{N}\sum_{j=1}^{M}\pi_{ij}c(\bm{y}_{j},\bm{x}_{i}).$

Going further, this problem can be generalized further to both continuous and
discrete probability measures. The sand piles and target locations are
respectively modeled by the discrete probability measures

$\nu=\sum_{j=1}^{M}\beta_{j}\delta_{\bm{y}_{j}},\quad\mu=\sum_{i=1}^{N}\alpha_{i}\delta_{\bm{x}_{i}}.$
(25)

The Monge problem corresponds to the minimization problem

$\inf_{T:\Omega\to\Omega}\left\\{\int_{\Omega}c(\bm{y},T(\bm{y}))\nu(\mathrm{d}\bm{y})\,\Big{|}\,\,\,T\\#\nu=\mu\right\\},$
(26)

where we recall that the push-forward measure $T\\#\nu$ is defined as the
measure on $\Omega$ such that

$T\\#\nu(\mathscr{B})=\nu(T^{-1}(\mathscr{B})),$

for all open sets $\mathscr{B}$.

The Kantorovich problem corresponds to the minimization problem

$\inf_{\Pi\in\mathcal{P}(\Omega\times\Omega)}\left\\{\int_{\Omega}\int_{\Omega}c(\bm{y},\bm{x})\Pi(\mathrm{d}\bm{y},\mathrm{d}\bm{x})\,\Big{|}\,\,\Pi_{1}=\nu,\,\,\Pi_{2}=\mu\right\\},$
(27)

where $\mathcal{P}(\Omega\times\Omega)$ denotes the set of probability
measures on $\Omega\times\Omega$ and $\Pi_{1}$ and $\Pi_{2}$ are the first and
second marginals of $\Pi\in\mathcal{P}(\Omega\times\Omega)$.

Although it can be checked that with the choice Eq. (25), the Monge and
Kantorovich problems reduce to the optimization problems introduced above, the
formulations Eqs. 26-(27) are much more general and are not restricted to the
case of discrete measures $\mu,\nu$. In particular, $\mu$ and/or $\nu$ can be
continuous measures (i.e. defined by a probability density function), in which
case the following fundamental theorem due to Brenier, Gangbo, McCann,
Caffarelli and others [17, 47, 25] can be proved.

######  Theorem 10.

If the source measure $\nu$ is a continuous measure and under some assumptions
on the cost function $c$, then the Monge problem Eq. (26) has a unique
solution $T$.

In our approach, we consider the semi-discrete case where $\nu$ is a
continuous measure and $\mu=\sum_{i}\alpha_{i}\delta_{\bm{x}_{i}}$ is a
discrete measure. In this case, a direct application of the previous theorem
shows that the optimal transport problem is equivalent to the computation of a
partition of the space $\Omega=\cup_{i}\mathscr{S}_{i}$ into the disjoint
sets:

$\mathscr{S}_{i}:=T^{-1}(\\{\bm{x}_{i}\\}).$

Moreover, a simple computation, presented in the main text in the special case
where $\nu$ is the Lebesgue measure, shows that
$\nu(\mathscr{S}_{i})=\alpha_{i}$. A graphical illustration of the discrete,
continuous and semi-discrete problems is shown in Fig. 16.

Discrete solver. To solve our discrete optimal transport problem efficiently,
we rely on the symmetric Sinkhorn algorithm with simulated annealing that is
presented in Algorithm 3.5 of [42]. As detailed in [44] and its appendices,
the Sinkhorn iterations of Lines 6 and 7 can be implemented in parallel on
both CPUs and GPUs, with a negligible memory footprint and in a numerically
stable way using the KeOps library.

We note that although the initializations of Lines 1 to 3 and the annealing
schedule of Line 4 have been tested extensively for the squared Euclidean cost
$c(x_{i},y_{j})=\|x_{i}-y_{h}\|^{2}$ on the unit domain $[0,1]^{d}$, they are
certainly sub-optimal for other cost functions. For the sake of safety, we
thus run all of our experiments with a very large value of
$n_{\text{iter}}=3000$. Going forward, designing an optimization procedure
that works well with mixtures of arbitrary cost functions for each cell is one
of our main priorities.

(a) Discrete OT

(b) Continuous OT

(c) Semi-discrete OT

Figure 16: Discrete, continuous and semi-discrete optimal transport. (a) A
discrete optimal transport problem is a matching point problem but which may
require mass splitting in which case there is no Monge map. (b) If the source
measure has a density, then there is a Monge map. (c) In the semi-discrete
case, the target measure is discrete and the Monge map defines a partition of
the space.

1:Algorithm 1:  The symmetric Sinkhorn algorithm, with $\varepsilon$-scaling

2:Input: number of iterations $n_{\text{iter}}$

3:Input: real-valued cost function $c(x_{i},y_{j})$ on a domain
$\Omega\subset\mathbb{R}^{d}$

4:Input: seed locations $x_{1},\dots,x_{N}$ and positive weights
$v_{1},\dots,v_{N}$

5:Input: pixel locations $y_{1},\dots,y_{M}$ associated to a uniform weight
$1/M$

6:Output: dual potentials $(f_{1},\dots,f_{N})$ and $(g_{1},\dots,g_{M})$

7:$R~{}\leftarrow~{}\tfrac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{M}c(x_{i},y_{j})\,v_{i}/M$

8:$f_{i}~{}\leftarrow~{}\sum_{j=1}^{M}c(x_{i},y_{j})/M~{}-~{}R$

9:$g_{j}~{}\leftarrow~{}\sum_{i=1}^{N}c(x_{i},y_{j})\,v_{i}~{}-~{}R$

10:for $k=0$ to  $n_{\text{iter}}$ do

11:
$\varepsilon~{}\leftarrow~{}4^{2}\cdot(1/4^{2}M^{2/d})^{k/n_{\text{iter}}}$
$\triangleright$ Exponential cool down from $4^{2}$ to $1/M^{2/d}$.

12:
$\widetilde{f}_{i}~{}\leftarrow~{}-\varepsilon\log\sum_{j=1}^{M}\exp\big{[}\,\tfrac{1}{\varepsilon}\big{(}g_{j}-c(x_{i},y_{j})\big{)}+\log(1/M)\,\big{]}~{}$

13:
$\widetilde{g}_{j}~{}\leftarrow~{}-\varepsilon\log\sum_{i=1}^{N}\exp\big{[}\,\tfrac{1}{\varepsilon}\big{(}f_{i}-c(x_{i},y_{j})\big{)}+\log(v_{i})\,\big{]}~{}$

14:
$f_{i},~{}g_{j}~{}\,\leftarrow~{}\tfrac{1}{2}(f_{i}+\widetilde{f}_{i}),~{}\tfrac{1}{2}(g_{j}+\widetilde{g}_{j})$

15:return $f_{i},~{}g_{j}$ $\triangleright$ The $f_{i}$’s correspond to the
potentials $w_{i}$ of Theorem 3.

## References

  * [1] Andreas Alpers, Andreas Brieden, Peter Gritzmann, Allan Lyckegaard and Henning Friis Poulsen “Generalized Balanced Power Diagrams for 3D Representations of Polycrystals” In _Philosophical Magazine_ 95.9, 2015, pp. 1016–1028 DOI: 10.1080/14786435.2015.1015469
  * [2] Silvanus Alt, Poulami Ganguly and Guillaume Salbreux “Vertex Models: From Cell Mechanics to Tissue Morphogenesis” In _Phil. Trans. R. Soc. B_ 372.1720, 2017, pp. 20150520 DOI: 10.1098/rstb.2015.0520
  * [3] Peter B. Armstrong “Cell Sorting Out: The Self-Assembly of Tissues _In Vitro_ ” In _Critical Reviews in Biochemistry and Molecular Biology_ 24.2, 1989, pp. 119–149 DOI: 10.3109/10409238909086396
  * [4] Utkarsh Ayachit “The paraview guide: a parallel visualization application” Kitware, Inc., 2015
  * [5] Daniel L. Barton, Silke Henkes, Cornelis J. Weijer and Rastko Sknepnek “Active Vertex Model for Cell-Resolution Description of Epithelial Tissue Mechanics” In _PLoS Comput. Biol._ 13.6, 2017, pp. e1005569 DOI: 10.1371/journal.pcbi.1005569
  * [6] Clemens Bechinger, Roberto Di Leonardo, Hartmut Löwen, Charles Reichhardt, Giorgio Volpe and Giovanni Volpe “Active Particles in Complex and Crowded Environments” In _Rev. Mod. Phys._ 88.4, 2016, pp. 045006 DOI: 10.1103/RevModPhys.88.045006
  * [7] Roman Belousov, Sabrina Savino, Prachiti Moghe, Takashi Hiiragi, Lamberto Rondoni and Anna Erzberger “When Time Matters: Poissonian Cellular Potts Models Reveal Nonequilibrium Kinetics of Cell Sorting”, 2023 arXiv:2306.04443 [cond-mat, physics:physics]
  * [8] Dapeng Bi, J. H. Lopez, J. M. Schwarz and M. Lisa Manning “A Density-Independent Rigidity Transition in Biological Tissues” In _Nature Physics_ 11.12, 2015, pp. 1074–1079 DOI: 10.1038/nphys3471
  * [9] Dapeng Bi, Xingbo Yang, M. Cristina Marchetti and M. Lisa Manning “Motility-Driven Glass and Jamming Transitions in Biological Tissues” In _Physical Review X_ 6.2, 2016, pp. 021011 DOI: 10.1103/PhysRevX.6.021011
  * [10] Martin Bock, Amit Kumar Tyagi, Jan-Ulrich Kreft and Wolfgang Alt “Generalized Voronoi Tessellation as a Model of Two-dimensional Cell Tissue Dynamics” In _Bull. Math. Biol._ 72.7, 2010, pp. 1696–1731 DOI: 10.1007/s11538-009-9498-3
  * [11] D. P. Bourne, P. J. J. Kok, S. M. Roper and W. D. T. Spanjer “Laguerre Tessellations and Polycrystalline Microstructures: A Fast Algorithm for Generating Grains of given Volumes” In _Philosophical Magazine_ 100.21, 2020, pp. 2677–2707 DOI: 10.1080/14786435.2020.1790053
  * [12] David P. Bourne “MATLAB-SDOT” URL: https://github.com/DPBourne/MATLAB-SDOT
  * [13] David P. Bourne, Bernhard Schmitzer and Benedikt Wirth “Semi-Discrete Unbalanced Optimal Transport and Quantization”, 2018 arXiv:1808.01962
  * [14] D.P. Bourne, M. Pearce and S.M. Roper “Geometric Modelling of Polycrystalline Materials: Laguerre Tessellations and Periodic Semi-Discrete Optimal Transport” In _Mechanics Research Communications_ 127, 2023, pp. 104023 DOI: 10.1016/j.mechrescom.2022.104023
  * [15] Y. Brenier and G. Loeper “A Geometric Approximation to the Euler Equations: The Vlasov-Monge-Ampère System” In _GAFA, Geom. funct. anal._ 14.6, 2004, pp. 1182–1218 DOI: 10.1007/s00039-004-0488-1
  * [16] Yann Brenier “A Combinatorial Algorithm for the Euler Equations of Incompressible Flows” In _Computer Methods in Applied Mechanics and Engineering_ 75.1-3, 1989, pp. 325–332 DOI: 10.1016/0045-7825(89)90033-9
  * [17] Yann Brenier “Polar Factorization and Monotone Rearrangement of Vector‐valued Functions” In _Comm. Pure Appl. Math._ 44.4, 1991, pp. 375–417 DOI: 10.1002/cpa.3160440402
  * [18] Yann Brenier “Derivation of the Euler Equations from a Caricature of Coulomb Interaction” In _Comm. Math. Phys._ 212.1, 2000, pp. 93–104 DOI: 10.1007/s002200000204
  * [19] Yann Brenier “Optimal Transportation of Particles, Fluids and Currents”, 2015, pp. 59–85 DOI: 10.2969/aspm/06710059
  * [20] G. Wayne Brodland “The Differential Interfacial Tension Hypothesis (DITH): A Comprehensive Theory for the Self-Rearrangement of Embryonic Cells and Tissues” In _Journal of Biomechanical Engineering_ 124.2, 2002, pp. 188–197 DOI: 10.1115/1.1449491
  * [21] G Wayne Brodland “Computational Modeling of Cell Sorting, Tissue Engulfment, and Related Phenomena: A Review” In _Applied Mechanics Reviews_ 57.1, 2004, pp. 47–76 DOI: 10.1115/1.1583758
  * [22] G Wayne Brodland and Helen H Chen “The Mechanics of Cell Sorting and Envelopment” In _Journal of Biomechanics_ 33.7, 2000, pp. 845–851 DOI: 10.1016/S0021-9290(00)00011-7
  * [23] G. Wayne Brodland and Helen H. Chen “The Mechanics of Heterotypic Cell Aggregates: Insights From Computer Simulations” In _Journal of Biomechanical Engineering_ 122.4, 2000, pp. 402–407 DOI: 10.1115/1.1288205
  * [24] Andreas Buttenschön and Leah Edelstein-Keshet “Bridging from Single to Collective Cell Migration: A Review of Models and Links to Experiments” In _PLoS Comput. Biol._ 16.12, 2020, pp. e1008411 DOI: 10.1371/journal.pcbi.1008411
  * [25] Luis Caffarelli “Allocation Maps with General Cost Functions” In _Partial Differential Equations and Applications. Collected Papers in Honor of Carlo Pucci_ CRC Press, 1996
  * [26] Jose A. Carrillo, Hideki Murakawa, Makoto Sato, Hideru Togashi and Olena Trush “A Population Dynamics Model of Cell-Cell Adhesion Incorporating Population Pressure and Density Saturation” In _J. Theoret. Biol._ 474, 2019, pp. 14–24 DOI: 10.1016/j.jtbi.2019.04.023
  * [27] Louis-Pierre Chaintron and Antoine Diez “Propagation of Chaos: A Review of Models, Methods and Applications. I. Models and Methods” In _Kinet. Relat. Models_ 15.6, 2022, pp. 895–1015 DOI: 10.3934/krm.2022017
  * [28] Louis-Pierre Chaintron and Antoine Diez “Propagation of Chaos: A Review of Models, Methods and Applications. II. Applications” In _Kinet. Relat. Models_ 15.6, 2022, pp. 1017–1173 DOI: 10.3934/krm.2022018
  * [29] B. Charlier, J. Feydy, J. A. Glaunès, F.-D. Collin and G. Durif “Kernel Operations on the GPU, with Autodiff, without Memory Overflows” In _J. Mach. Learn. Res._ 22.74, 2021, pp. 1–6 URL: https://www.kernel-operations.io/
  * [30] Hugues Chaté, Francesco Ginelli, Guillaume Grégoire and Franck Raynaud “Collective Motion of Self-Propelled Particles Interacting without Cohesion” In _Phys. Rev. E_ 77.4, 2008, pp. 046113 DOI: 10.1103/PhysRevE.77.046113
  * [31] Helen H. Chen and G. Wayne Brodland “Cell-Level Finite Element Studies of Viscous Cells in Planar Aggregates” In _Journal of Biomechanical Engineering_ 122.4, 2000, pp. 394–401 DOI: 10.1115/1.1286563
  * [32] Roie Cohen, Liat Amir-Zilberstein, Micha Hersch, Shiran Woland, Olga Loza, Shahar Taiber, Fumio Matsuzaki, Sven Bergmann, Karen B. Avraham and David Sprinzak “Mechanical Forces Drive Ordered Patterning of Hair Cells in the Mammalian Inner Ear” In _Nat. Commun._ 11.1, 2020, pp. 5137 DOI: 10.1038/s41467-020-18894-8
  * [33] Fergus Cooper, Ruth Baker, Miguel Bernabeu, Rafel Bordas, Louise Bowler, Alfonso Bueno-Orovio, Helen Byrne, Valentina Carapella, Louie Cardone-Noott, Jonathan Cooper, Sara Dutta, Benjamin Evans, Alexander Fletcher, James Grogan, Wenxian Guo, Daniel Harvey, Maurice Hendrix, David Kay, Jochen Kursawe, Philip Maini, Beth McMillan, Gary Mirams, James Osborne, Pras Pathmanathan, Joe Pitt-Francis, Martin Robinson, Blanca Rodriguez, Raymond Spiteri and David Gavaghan “Chaste: Cancer, Heart and Soft Tissue Environment” In _Journal of Open Source Software_ 5.47, 2020, pp. 1848 DOI: 10.21105/joss.01848
  * [34] Nicholas C. Darnton, Linda Turner, Svetlana Rojevsky and Howard C. Berg “Dynamics of Bacterial Swarming” In _Biophysical Journal_ 98.10, 2010, pp. 2082–2090 DOI: 10.1016/j.bpj.2010.01.053
  * [35] Pierre-Gilles De Gennes, Françoise Brochard-Wyart and David Quéré “Capillarity and Wetting Phenomena” New York, NY: Springer New York, 2004
  * [36] Fernando De Goes, Corentin Wallez, Jin Huang, Dmitry Pavlov and Mathieu Desbrun “Power Particles: An Incompressible Fluid Solver Based on Power Diagrams” In _ACM Trans. Graph._ 34.4, 2015, pp. 1–11 DOI: 10.1145/2766901
  * [37] Frédéric De Gournay, Jonas Kahn and Léo Lebrat “Differentiation and Regularity of Semi-Discrete Optimal Transport with Respect to the Parameters of the Discrete Measure” In _Numer. Math._ 141.2, 2019, pp. 429–453 DOI: 10.1007/s00211-018-1000-4
  * [38] Guido De Philippis and Alessio Figalli “The Monge–Ampère Equation and Its Link to Optimal Transportation” In _Bull. Amer. Math. Soc._ 51.4, 2014, pp. 527–580 DOI: 10.1090/S0273-0979-2014-01459-4
  * [39] Pierre Degond, Amic Frouvelle and Sara Merino-Aceituno “A New Flocking Model through Body Attitude Coordination” In _Math. Models Methods Appl. Sci._ 27.06, 2017, pp. 1005–1049 DOI: 10.1142/S0218202517400085
  * [40] Antoine Diez “SiSyPHE: A Python Package for the Simulation of Systems of Interacting Mean-Field Particles with High Efficiency” In _Journal of Open Source Software_ 6.65, 2021, pp. 3653 DOI: 10.21105/joss.03653
  * [41] Reza Farhadifar, Jens-Christian Röper, Benoit Aigouy, Suzanne Eaton and Frank Jülicher “The Influence of Cell Mechanics, Cell-Cell Interactions, and Proliferation on Epithelial Packing” In _Current Biology_ 17.24, 2007, pp. 2095–2104 DOI: 10.1016/j.cub.2007.11.049
  * [42] Jean Feydy “Geometric Data Analysis, beyond Convolutions”, 2020
  * [43] Jean Feydy, Benjamin Charlier, François-Xavier Vialard and Gabriel Peyré “Optimal transport for diffeomorphic registration” In _Medical Image Computing and Computer Assisted Intervention- MICCAI 2017: 20th International Conference, Quebec City, QC, Canada, September 11-13, 2017, Proceedings, Part I 20_ , 2017, pp. 291–299 Springer
  * [44] Jean Feydy, Alexis Glaunès, Benjamin Charlier and Michael Bronstein “Fast geometric learning with symbolic matrices” In _Advances in Neural Information Processing Systems_ 33, 2020, pp. 14448–14462
  * [45] Jean Feydy, Thibault Séjourné, François-Xavier Vialard, Shun-ichi Amari, Alain Trouve and Gabriel Peyré “Interpolating between Optimal Transport and MMD Using Sinkhorn Divergences” In _Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics_ PMLR, 2019, pp. 2681–2690 URL: https://proceedings.mlr.press/v89/feydy19a.html
  * [46] Thomas O. Gallouët and Quentin Mérigot “A Lagrangian Scheme à La Brenier for the Incompressible Euler Equations” In _Foundations of Computational Mathematics_ 18.4, 2018, pp. 835–865 DOI: 10.1007/s10208-017-9355-y
  * [47] Wilfrid Gangbo and Robert J. McCann “The Geometry of Optimal Transportation” In _Acta Math._ 177.2, 1996, pp. 113–161 DOI: 10.1007/BF02392620
  * [48] Ahmadreza Ghaffarizadeh, Randy Heiland, Samuel H. Friedman, Shannon M. Mumenthaler and Paul Macklin “PhysiCell: An Open Source Physics-Based Cell Simulator for 3-D Multicellular Systems” In _PLoS Comput. Biol._ 14.2, 2018, pp. e1005991 DOI: 10.1371/journal.pcbi.1005991
  * [49] James A. Glazier and François Graner “Simulation of the Differential Adhesion Driven Rearrangement of Biological Cells” In _Phys. Rev. E_ 47.3, 1993, pp. 2128–2154 DOI: 10.1103/PhysRevE.47.2128