# Differentiable Voronoi Diagrams for Simulation of Cell-Based Mechanical
Systems

Logan Numerow  ETH ZürichSwitzerland<EMAIL_ADDRESS>,  Yue Li  ETH
ZürichSwitzerland<EMAIL_ADDRESS>,  Stelian Coros  ETH ZürichSwitzerland
<EMAIL_ADDRESS>and  Bernhard Thomaszewski  ETH ZürichSwitzerland
<EMAIL_ADDRESS>
###### Abstract.

Navigating topological transitions in cellular mechanical systems is a
significant challenge for existing simulation methods. While abstract models
lack predictive capabilities at the cellular level, explicit network
representations struggle with topology changes, and per-cell representations
are computationally too demanding for large-scale simulations. To address
these challenges, we propose a novel cell-centered approach based on
differentiable Voronoi diagrams. Representing each cell with a Voronoi site,
our method defines shape and topology of the interface network implicitly. In
this way, we substantially reduce the number of problem variables, eliminate
the need for explicit contact handling, and ensure continuous geometry changes
during topological transitions. Closed-form derivatives of network positions
facilitate simulation with Newton-type methods for a wide range of per-cell
energies. Finally, we extend our differentiable Voronoi diagrams to enable
coupling with arbitrary rigid and deformable boundaries. We apply our approach
to a diverse set of examples, highlighting splitting and merging of cells as
well as neighborhood changes. We illustrate applications to inverse problems
by matching soap foam simulations to real-world images. Comparative analysis
with explicit cell models reveals that our method achieves qualitatively
comparable results at significantly faster computation times.

Voronoi Diagram, Power Diagram, Foam Simulation, Differentiable Simulation

††submissionid: 319††copyright: acmlicensed††journal: TOG††journalyear:
2024††journalvolume: 43††journalnumber: 4††article: 69††publicationmonth:
7††doi: 10.1145/3658152††ccs: Computing methodologies Physical simulation

Figure 1. We use differentiable Voronoi diagrams to simulate embryonic
cleavage in an elastic spherical membrane. Bottom: As cells split according to
a predefined schedule, substantial neighborhood changes and cell deformations
occur. Each cell is modeled using only a single Voronoi site carrying four
degrees of freedom. Top: Images from Jan van Ijken’s ”Becoming” (2018),
depicting the embryonic development of a salamander, are included for
comparison.

##  1\. Introduction

A central characteristic of cellular mechanical systems is that cells change
neighborhoods during dynamic motion and growth-induced splitting (as in
biological tissue), or through merging and collapse (as in soap foams). These
topological transitions cause abrupt force changes, which make computational
modeling very challenging.

There are several existing approaches in the literature for simulation of
cell-based systems. One option is to abstract away the cellular structure and
model the large-scale mechanics of the system in a homogenized sense.
Evidently, such models cannot predict interactions at the cellular level
governed by local topology changes. Another approach is to explicitly model
the network of cell interfaces, enabling efficient computation of cell volumes
and other geometric quantities that influence their mechanics. However,
topology changes are difficult to model with explicit network representations,
requiring discrete mesh adaptations and other discontinuous operations.
Another class of models represents each cell individually using, _e.g_., per-
cell finite element discretizations. While this approach allows for accurate
modeling of cell interactions, it requires a very large number of degrees of
freedom to resolve cell deformations, and contact must be handled explicitly
at every cell interface. These aspects render this strategy unsuitable for
large-scale simulations involving hundreds or thousands of cells.

To overcome the limitations of existing methods, we propose a novel cell-
centered approach based on generalized Voronoi diagrams. Each cell is
represented by a Voronoi site, and the topology and shape of the cells are
defined implicitly. This representation uses far fewer degrees of freedom,
avoids explicit handling of contact between cells, and ensures continuous cell
geometry changes through topological transitions. In order to use this Voronoi
model for simulation, we derive closed-form expressions for first and second
derivatives of cell geometry with respect to site positions, allowing us to
use Newton-type optimization solvers for a wide range of energies. We
furthermore obtain equilibrium-state derivatives of the interface network
using sensitivity analysis, thus opening the door to a wide range of inverse
simulation problems. Finally, we extend our differentiable Voronoi formulation
to enable coupling with arbitrary rigid and deformable boundaries.

We demonstrate the potential of our method on a diverse set of examples,
including tissue growth and foam coarsening. We furthermore show applications
of differentiable Voronoi diagrams to inverse problems by optimizing per-cell
pressures to match input images from real-world foams. Our comparison with
explicit cell models indicates that differentiable Voronoi diagrams can
produce qualitatively similar behavior with much faster computation times.

##  2\. Related Work

#### Simulating Cell & Tissue Dynamics

Various computational models have been developed to simulate cell-based
mechanical systems. One approach is to abstract away the cellular structure
and model the mechanical properties of the system in an averaged sense at
coarser scales (Brodland et al., 2006; Kondo et al., 2018; Piatnitski and
Ptashnyk, 2020). Our focus, however, is on problems that are governed by cell-
level phenomena. One option in this context is to explicitly model the network
of cell interfaces using so-called vertex models (Honda et al., 2004; Fletcher
et al., 2014; Alt et al., 2017). Although extensions have been developed to
allow for topology changes (Okuda et al., 2012; Farhadifar et al., 2007), the
discrete mesh modifications involved in these methods are problematic for
differentiable simulation. Moreover, explicit vertex models require additional
constraints to prevent self-intersections and other invalid geometry (Vedel-
Larsen, 2010). Our differentiable Voronoi diagrams, in contrast, ensure a
valid interface network at all times and guarantee continuous topology changes
for any displacement of the Voronoi sites.

Another line of work has explored so-called deformable cell models, in which
cells are represented through separate surface meshes (Kim et al., 2021;
Conradin et al., 2021). While this approach enables accurate modeling of
contact forces and topology changes, it requires explicit handling of
collisions between cells. Moreover, hundreds to thousands of degrees of
freedom are needed per cell, which rapidly leads to all but intractable
problem size. With each cell encoded by the position of a single Voronoi site,
our approach requires dramatically fewer degrees of freedom. The resulting
model is highly scalable, enabling efficient gradient-based optimization with
thousands of cells.

#### Simulating Soap Films, Foams, & Bubbles

The dynamics of bubbles have attracted substantial interest from the graphics
community. Examples include the formation of bubbles in water (Hong et al.,
2008), soap film dynamics (Ishida et al., 2020; Huang et al., 2020; Deng et
al., 2022), and general foams (Busaryev et al., 2012). Various representations
have been explored, ranging from explicit surface meshes (Da et al., 2015) to
mixed particle-and-grid discretizations (Goldade et al., 2020). We show that,
using only a single site per foam cell, differentiable Voronoi diagrams can
produce realistic 3D foam coarsening simulations. Our method can likewise be
used for inverse problems such as matching soap film simulations to real-world
images.

#### Implicit Voronoi Diagrams

Voronoi diagrams have numerous applications in computer graphics, vision, and
beyond. They have been explored as a basis for, _e.g_., fluid (de Goes et al.,
2015) and foam (Busaryev et al., 2012) simulations, image segmentation
(Williams et al., 2020) as well as for structural design. Zhang _et al_.
(2017) design personalized medical braces represented as centroidal Voronoi
diagrams. Their algorithm uses a variant of Lloyd’s method to iteratively
adjust site positions in order to minimize the design objective. Lumpe _et
al_. (2022) investigate power diagrams for topology optimization. While cell
volumes are adapted in each iteration, all simulations are performed on finite
element meshes generated from the power diagrams. Also targeting topology
optimization, Feng _et al_. (2022) introduce a smooth relaxation of Voronoi
diagrams based on softmax functions. However, their approach ultimately
requires rasterization on a Cartesian grid for optimization.

We refer to the aforementioned Voronoi models as implicit Voronoi diagrams, as
the simulated structures are defined implicitly by Voronoi sites rather than
with explicit degrees of freedom. Lloyd’s algorithm (1982) for generating
centroidal Voronoi tessellations (CVT) is arguably the most widely used method
for optimization of Voronoi diagrams. Nevertheless, substantial accelerations
can be achieved when using gradient-based or quasi-Newton methods (Liu et al.,
2009; Yan et al., 2009, 2011; Wang et al., 2015). Unlike these methods, our
work provides explicit derivatives of Voronoi vertices, opening the door to
general energy functions that can be adapted to a large variety of physical
systems. In addition, we provide an algebraic recipe for unified treatment of
Voronoi cell faces and clipping geometry, enabling simulation of cellular
systems coupled to moving or deformable boundaries.

#### Adaptive Meshes

Adapting meshes to ensure given quality criteria is a common strategy in many
simulation contexts, including large elastoplastic deformations (Wicke et al.,
2010; Bargteil et al., 2007; Wojtan et al., 2009), fracture (O’Brien et al.,
2002; Pfaff et al., 2014; Koschier et al., 2014), fluid flow (Misztal et al.,
2014), and the folding of thin sheet materials (Narain et al., 2012, 2013;
Schreck et al., 2016; Zhang et al., 2022). Remeshing is typically applied as a
post-process after each time step, which can disturb force equilibrium and
lead to abrupt configuration changes. Coupling remeshing with time integration
can avoid such artifacts (Ferguson et al., 2023). While mesh adaptation is
often an optional improvement, topological changes that occur, e.g., during
droplet formation make remeshing compulsory (Brochu and Bridson, 2009). The
applications that we consider in this work likewise involve inherent topology
changes as cells split, collapse, or switch neighbors. However, our
differentiable Voronoi diagrams eliminate the need for explicit remeshing, as
the mesh is given implicitly through the Voronoi sites.

#### Equilibrium-Constrained Optimization

Differentiable simulation is a core technology for solving equilibrium-
constrained optimization problems in graphics and beyond. Examples include
material parameter estimation (Miguel et al., 2012; Hahn et al., 2019; Yan et
al., 2018; Bickel et al., 2009), elastic shape optimization (Chen et al.,
2014; Ly et al., 2018; Panetta et al., 2017), inflatables (Skouras et al.,
2014; Panetta et al., 2021), clothing design (Umetani et al., 2011; Wang,
2018; Montes et al., 2020), and trajectory optimization (Geilinger et al.,
2020; Hu et al., 2019). Unlike existing methods that impose equilibrium
constraints on explicit triangle or tetrahedron meshes, our formulation
enables the solution of inverse problems defined on Voronoi diagrams.

##  3\. Differentiable Restricted Voronoi Diagrams

###  3.1. Generalized Voronoi Diagrams

The Voronoi cell $\mathcal{R}_{i}$ for a given site $C_{i}$ can be defined as

(1) $\mathcal{R}_{i}=\\{p:d(p,C_{i})<d(p,C_{j})\forall j\neq i\\}\ ,$

where $d(p,C)=d_{\text{Euclidean}}(p,C)$ is the Euclidean distance between
point $p$ and site $C$. A family of more flexible generalized Voronoi diagrams
can be constructed by modifying the distance metric which defines the Voronoi
cell.

Mouzarkel (1997) and Eppstein (2012) draw theoretical links between foam
structures and generalized Voronoi diagrams, independently showing that a two-
dimensional foam at equilibrium is geometrically equivalent to a sectional
multiplicative Voronoi partition (SMVP),

(2)
$d_{\text{SMVP}}(p,C_{i})^{2}=\frac{1}{k_{i}}\left(d_{\text{Euclidean}}(p,C_{i})^{2}-w_{i}\right)\
,$

for some arrangement of sites and weights $w$, $k$. The multiplicative weight
$k$ gives rise to curved interfaces which, while necessary to represent exact
foam structures, significantly complicate computations of geometric properties
of cells (Sec. 3.3). Analytic integration over 3D curved surfaces such as soap
films is typically impossible, and algebraic line-surface intersection
computations, required for boundary clipping, involve solving high-order
polynomials even for simple 3D interpolating surfaces. We instead use the
power diagram,

(3) $d_{\text{Power}}(p,C_{i})^{2}=d_{\text{Euclidean}}(p,C_{i})^{2}-w_{i}\ ,$

for our simulations, which maintains planar interfaces while allowing for cell
size variability. Sec. 4.4 demonstrates the power diagram’s capacity to
express realistic foam geometry.

###  3.2. Voronoi Vertices

We compute the restricted Voronoi diagram, constrained to within a bounding
domain, as described by Yan _et al_. (2009; 2011). Each restricted Voronoi
vertex in $m$ dimensions is an intersection of $m$ hyperplanes, each of which
is either a Voronoi bisector or a facet of the bounding surface. In three
dimensions, there are four types of vertices; (a) unrestricted Voronoi
vertices, _i.e_., intersections between three

(linearly independent) bisectors, (b) boundary face vertices, _i.e_.,
intersections between two bisectors and one boundary facet, (c) boundary edge
vertices, _i.e_., intersections between one bisector and two boundary facets,
and (d) simple boundary vertices. The generating plane equations can be
expressed as simple functions of the input degrees of freedom, _i.e_., the
Voronoi sites and the vertices of the boundary surface. Given as the
intersection of three planes, _i.e_., the solution of a $3\times 3$ linear
system, each 3D Voronoi vertex can be expressed in closed-form and
differentiated with respect to the inputs.

We derive as an example the closed-form expression for the unrestricted
Voronoi vertex in 3D. Given four sites
$\mathbf{x}_{i}=(x_{i},y_{i},z_{i})^{T}$, $i\in\\{0,1,2,3\\}$ of a Delaunay
simplex, the plane equations of Voronoi bisectors are

(4)
$(x_{j}-x_{i})x+(y_{j}-y_{i})y+(z_{j}-z_{i})z=\frac{1}{2}(x_{j}^{2}-x_{i}^{2}+y_{j}^{2}-y_{i}^{2}+z_{j}^{2}-z_{i}^{2})\
.$

The intersection $\mathbf{x}_{c}=(x_{c},y_{c},z_{c})^{T}$ of three bisectors
solves the linear system

$\begin{bmatrix}x_{1}-x_{0}&y_{1}-y_{0}&z_{1}-z_{0}\\\
x_{2}-x_{0}&y_{2}-y_{0}&z_{2}-z_{0}\\\
x_{3}-x_{0}&y_{3}-y_{0}&z_{3}-z_{0}\end{bmatrix}\begin{bmatrix}x_{c}\\\
y_{c}\\\
z_{c}\end{bmatrix}=\frac{1}{2}\begin{bmatrix}x_{1}^{2}-x_{0}^{2}+y_{1}^{2}-y_{0}^{2}+z_{1}^{2}-z_{0}^{2}\\\
x_{2}^{2}-x_{0}^{2}+y_{2}^{2}-y_{0}^{2}+z_{2}^{2}-z_{0}^{2}\\\
x_{3}^{2}-x_{0}^{2}+y_{3}^{2}-y_{0}^{2}+z_{3}^{2}-z_{0}^{2}\end{bmatrix},$

for which the solution can be written using the known closed-form expression
for the $3\times 3$ matrix inverse. This approach extends to all types of
boundary vertices by substituting bisector plane equations in the linear
system for plane equations of boundary facets. Furthermore, it is easily
adapted for the power diagram by adding $\frac{1}{2}(w_{i}-w_{j})$ to the
right hand side of the bisector plane equation (4). In our implementation,
code for evaluating these expressions and their derivatives is generated using
a computer algebra system.

###  3.3. Integration over Voronoi Cells

The mechanics of cell-based systems are driven by the energies of individual
cells. Cell energies relate to properties such as surface area and volume,
which are computed by integrating over the cell geometry. Volume integrals are
also used to compute cell centroids and moments of inertia. Integration over a
polyhedral Voronoi cell is performed piecewise by dividing the surface and
volume into simplices. A surface integral requires triangulating each facet of
the cell and then summing the per-triangle values. A tetrahedralization of the
volume is, however, unnecessary for computing volume integrals. Rather,
tetrahedra are constructed by joining each surface triangle to the origin as
shown in Fig. 2.

Figure 2. Volume integration over a 3D Voronoi cell. From left to right: (1)
3D Voronoi cell. (2) Face triangulation. (3) Integration tetrahedron formed by
connecting a face triangle to the origin. Tetrahedra can have positive (3) or
negative (4) volume.

In addition to avoiding an explicit tetrahedralization, including the origin
in each tetrahredron simplifies the simplex integral to a function of three
vertices instead of four. Importantly, this reduces the number of first and
second derivatives required per simplex. A generic approach for volume
integration in these tetrahedra is given in App. A.

###  3.4. Energy Formulation

We simulate cell mechanics in the implicit Voronoi model by considering a
potential energy of the form

(5) $E=F(\mathbf{c},\mathbf{x}(\mathbf{c})),$

where $\mathbf{c}$ is a vector containing the degrees of freedom of all
Voronoi sites, and $\mathbf{x}$ is a vector of spatial coordinates of all
Voronoi vertices, defined implicitly in terms of $\mathbf{c}$. In general, an
energy depending only on $\mathbf{x}$ is sufficient to describe physical
interactions between cells, and explicit dependence on $\mathbf{c}$ is used
for regularization. The generic gradient of this energy is

(6) $\frac{\mathrm{d}E}{\mathrm{d}\mathbf{c}}=\frac{\partial
F}{\partial\mathbf{x}}\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\mathbf{c}}+\frac{\partial
F}{\partial\mathbf{c}}\ ,$

and the Hessian follows as

(7) $\displaystyle\frac{\mathrm{d}^{2}E}{\mathrm{d}\mathbf{c}^{2}}$
$\displaystyle=\left(\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\mathbf{c}}\right)^{\top}\frac{\partial^{2}F}{\partial\mathbf{x}^{2}}\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\mathbf{c}}+\left(\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\mathbf{c}}\right)^{\top}\frac{\partial^{2}F}{\partial\mathbf{x}\partial\mathbf{c}}$
$\displaystyle+\frac{\partial^{2}F}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\mathbf{c}}+\sum_{i}\left(\frac{\partial
F}{\partial\mathbf{x}_{i}}\frac{\mathrm{d}^{2}\mathbf{x}_{i}}{\mathrm{d}\mathbf{c}^{2}}\right)+\frac{\partial^{2}F}{\partial\mathbf{c}^{2}}\
.$

This flexible formulation permits any differentiable energy function and is
applicable to a wide range of cell-based systems. The same computational
method can be used to simulate foam as well as biological tissue, as
demonstrated in Sec. 4.

###  3.5. Boundary Coupling

Many cellular systems interface with external bodies or free space. Biological
tissues grow freely or within deformable membranes, while engineering designs
that utilize cellular structures must support external loads. To handle these
external interactions, we define the bounding domain of the Voronoi diagram by
a triangular mesh with vertices $\mathbf{v}(\mathbf{p})$, where $\mathbf{p}$
are boundary degrees of freedom. We let $\mathbf{y}=(\mathbf{c},\mathbf{p})$
denote the vector holding all degrees of freedom of the coupled system and
consider a total potential energy

(8)
$E=F\left(\mathbf{c},\mathbf{x}\left(\mathbf{c},\mathbf{v}(\mathbf{p})\right)\right)+F_{B}(\mathbf{p})\
,$

where $F$ is the potential energy of the cells and $F_{B}$ is an additional
boundary energy which may encode, _e.g_., the elastic potential of an
enclosing membrane. As for the static boundary formulation, we compute the
first and second derivatives of the above expression to use in gradient-based
optimization. The gradient with respect to $\mathbf{p}$ is given by

(9) $\frac{\mathrm{d}E}{\mathrm{d}\mathbf{p}}=\frac{\partial
F}{\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathbf{p}}+\frac{\mathrm{d}F_{B}}{\mathrm{d}\mathbf{p}}\
.$

The additional blocks of the Hessian matrix are given in App. B.

Figure 3. Two simulations of a rigid body propelled through a tissue-like
assembly of elastic Voronoi cells. Cells are colored based on starting
position. Turbulence-like effects result in irregular displacement of cells
from their starting positions. The asymmetric shape of the second body results
in a curved trajectory due to forces from neighboring cells.

This boundary coupling formulation lends extreme flexibility to the model,
enabling systems involving elastic membranes, free surfaces, and rigid bodies
in addition to fixed boundaries.

###  3.6. Smoothness

In the general configuration, where each Voronoi vertex in $m$ dimensions is
of degree $m+1$, the topology of the Voronoi diagram is constant within a
neighborhood of the given state, and the vertex coordinates $\mathbf{x}$ are a
locally smooth function of the Voronoi sites $\mathbf{c}$. States with higher-
degree vertices, which occur during topological transitions, have undefined
$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\mathbf{y}}$ which gives rise to non-
smooth energies. Cell volume (area in 2D) and other volume integrals of smooth
functions are $C^{1}$-continuous across these transitional states, though
surface area (perimeter in 2D) is $C^{0}$ as noted in (Bogosel and Oudet,
2022) for the 2D case.

If cell energy is monotonic in surface area and the domain boundary is convex,
the gradient discontinuities are strictly concave, forming local maxima in the
energy landscape which do not hinder minimization. Otherwise, the model may
fail to converge using such energy functions. Our simulations with non-convex
domain boundaries (Secs. 4.1 and 4.3) use $C^{1}$-continuous energies and
therefore avoid these problems.

Degenerate configurations exist with lower-order discontinuities, _e.g_. when
a parallel boundary facet and Voronoi bisector overlap (discussed in (Liu et
al., 2009)) or when multiple Voronoi sites coincide. However, we have not
observed these states during simulation and optimization in any of our
examples.

###  3.7. Dynamics

We have so far considered systems at static equilibrium in which the net
internal force vanishes,

(10) $F(\mathbf{y})=-\frac{\mathrm{d}E}{\mathrm{d}\mathbf{y}}=\mathbf{0}\ .$

We now relax the equilibrium constraint to consider dynamic cell-based systems
governed by Newton’s second law,

(11)
$F(\mathbf{y})=\mathbf{M}\frac{\mathrm{d}^{2}\mathbf{y}}{\mathrm{d}t^{2}}\ ,$

where $\mathbf{M}$ is a (typically diagonal) mass matrix. We additionally
introduce viscous forces
$F_{f}(\mathbf{y})=-\boldsymbol{\eta}\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t}$,
where $\boldsymbol{\eta}$ is a diagonal viscosity matrix. The resulting
equations of motion are

(12)
$\mathbf{M}\ddot{\mathbf{y}}+\boldsymbol{\eta}\dot{\mathbf{y}}+\frac{\mathrm{d}E}{\mathrm{d}\mathbf{y}}=\mathbf{0}\
,$

where $\dot{\mathbf{y}}$ and $\ddot{\mathbf{y}}$ denote first and second
derivatives of $\mathbf{y}$ with respect to time. Unless otherwise noted, we
approximate these derivatives using standard finite-differencing, _i.e_.,

(13)
$\begin{gathered}\dot{\mathbf{y}}\approx\frac{1}{h}(\mathbf{y}_{k+1}-\mathbf{y}_{k})\
,\\\
\ddot{\mathbf{y}}\approx\frac{1}{h^{2}}(\mathbf{y}_{k+1}-2\mathbf{y}_{k}+\mathbf{y}_{k-1})\
.\\\ \end{gathered}$

If the energy function is at least $C^{1}$-continuous, quadratic convergence
can be achieved using a second-order discretization (BDF2) for the
accelerations,

(14)
$\ddot{\mathbf{y}}\approx\frac{1}{2h^{2}}(3\mathbf{y}_{k+1}-7\mathbf{y}_{k}+5\mathbf{y}_{k-1}-\mathbf{y}_{k-2})\
.$

Evidence of quadratic convergence under refinement of the simulation time step
is shown in App. C. The equations of motion can be solved for the next state
$\mathbf{y}_{k+1}$ via the optimization problem

(15)
$\mathbf{y}_{k+1}=\operatorname*{argmin}_{\mathbf{y}}\frac{a_{2}}{2}\ddot{\mathbf{y}}^{\top}\mathbf{M}\ddot{\mathbf{y}}+\frac{a_{1}}{2}\dot{\mathbf{y}}^{\top}\boldsymbol{\eta}\dot{\mathbf{y}}+E(\mathbf{y})\
,$

where $a_{1},a_{2}$ are the $\mathbf{y}_{k+1}$-coefficients of the first and
second finite-difference derivatives. Many cell-based systems operate at low
Reynolds numbers, where inertia is dominated by viscous effects. For such
cases, the momentum term can be neglected, leaving

(16)
$\mathbf{y}_{k+1}=\operatorname*{argmin}_{\mathbf{y}}\frac{a_{1}}{2}\dot{\mathbf{y}}^{\top}\boldsymbol{\eta}\dot{\mathbf{y}}+E(\mathbf{y})\
.$

(a) Frame 0, $n=2000$.

(b) Frame 31, $n=794$.

(c) Frame 50, $n=282$.

(d) Frame 65, $n=121$.

(e) Frame 110, $n=29$.

(f) Final equilibrium, $n=2$.

Figure 4. Coarsening of an initially monodisperse dry foam in a flattened box.
2000 cells collapse to a two-cell equilibrium in 350 frames.

###  3.8. Equilibrium-Constrained Optimization

In addition to forward simulation, our approach can be extended to solve
inverse simulation tasks formulated as equilibrium-constrained optimization
problems. The goal is then to find optimal parameter values for a cell-based
system such that the corresponding equilibrium state best satisfies a given
objective. We demonstrate in Sec. 4.4 that this framework is useful for
characterizing real-world cell-based systems from images.

A general equilibrium-constrained optimization problem reads

(17)
$\min_{\mathbf{y},\mathbf{u}}L(\mathbf{y},\mathbf{u})\quad\text{s.t.}\quad\frac{\partial}{\partial\mathbf{y}}E(\mathbf{y},\mathbf{u})=\mathbf{0}\
,$

where $L$ is an objective function, $\mathbf{y}$ are dynamic degrees of
freedom, and $\mathbf{u}$ are optimization variables. The system is
constrained to be at equilibrium, hence the state $\mathbf{y}$ is a function
of $\mathbf{u}$. The gradient is given by

(18) $\frac{\mathrm{d}L}{\mathrm{d}\mathbf{u}}=\frac{\partial
L}{\partial\mathbf{y}}\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}\mathbf{u}}+\frac{\partial
L}{\partial\mathbf{u}}\ ,$

where the simulation derivatives
$\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}\mathbf{u}}$ can be computed using the
implicit function theorem,

(19)
$\frac{\mathrm{d}}{\mathrm{d}\mathbf{u}}\frac{\partial}{\partial\mathbf{y}}E(\mathbf{y},\mathbf{u})=\mathbf{0}\
\rightarrow\
\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}\mathbf{u}}=-\left(\frac{\partial^{2}E}{\partial\mathbf{y}^{2}}\right)^{-1}\frac{\partial^{2}E}{\partial\mathbf{y}\partial\mathbf{u}}\
.$

Having evaluated the objective gradient with respect to the optimization
variables $\mathbf{u}$, we can minimize the $L$ using first-order or quasi-
Newton methods. Whenever the optimization parameters are updated during this
process, we recompute the equilibrium state $\mathbf{y}(\mathbf{u})$ and
ensure monotonic decrease of the objective using a backtracking line search.

###  3.9. Implementation

The method is implemented in C++ using the Eigen library (Guennebaud et al.,
2010) for matrix and vector operations, and the CHOLMOD solver (Chen et al.,
2008) for linear systems. Geometry processing libraries CGAL (Yvinec, 2023)
and Geogram (Lévy, 2015) are used in Voronoi diagram generation. Finally, we
use Polyscope (Sharp et al., 2019) for visualization and rendering of figures.
Our code can be found at https://github.com/lnumerow-ethz/VoronoiCellSim.

##  4\. Results

###  4.1. Intercellular Navigation with Rigid Body

In this first 2D experiment, we simulate the propulsion of a rigid body
through a tissue-like assembly of elastic Voronoi cells. The cell energy is
given by

(20)
$E=a_{0}(A-\bar{A})^{2}+a_{1}\int_{x\in\mathcal{R}}\|x-\bar{x}\|^{2}dA+a_{2}\|c-\bar{x}\|^{2}\
,$

where $A$ and $\bar{A}$ are the current and target area of the cell, $\bar{x}$
is its centroid, and $c$ is the Voronoi site position. In this and the
following examples, the final term is a regularizer on the distance between
the site and its cell centroid. The rigid body, modelled as an internal
clipping geometry with translation and rotation degrees of freedom (see Sec.
3.5), is subject to a constant force $F_{x}$ to the right. We use a dynamics
model with momentum (Eq. 15) to resolve turbulence-like effects in the system,
the results of which can be observed in Fig. 3. The curved trajectory of the
rigid body in the second example is a result of its asymmetric shape and
forces applied by the surrounding cells. This demonstrates the dynamic
coupling between the cells and boundary degrees of freedom.

For this and the following experiments, refer to App. D for parameter values
and additional details for reproducibility.

###  4.2. Foam Coarsening

Due to pressure differences between bubbles in an equilibriated liquid foam,
gas diffuses slowly through the interfaces from smaller, higher-pressure
bubbles into larger, lower-pressure ones. As a result, smaller foam cells
collapse and the total number of foam cells in the system decreases over time
in a process known as coarsening. This phenomenon has seen extensive research
(Vedel-Larsen, 2010; Thomas et al., 2015) due to its detrimental effect in
applications such as enhanced oil recovery. Foam coarsening presents a
challenge for vertex models due to the topological changes which result from
the collapse and removal of a cell, particularly in 3D. Existing 2D vertex
models (Vedel-Larsen, 2010) rely on heuristics to remove collapsing cells and
reconstruct topology, resulting in discontinuous changes to neighboring cell
volumes and to the system energy.

We simulate dry foam coarsening using a momentumless dynamic model as
described in Eq. 16. Each site has five degrees of freedom: three spatial
coordinates $(x,y,z)$, a power diagram weight $w$ and a volume target
$\bar{V}$ representing the mass of air contained in the cell. The energy of a
cell is

(21)
$E=a_{0}\left(\frac{V}{\bar{V}}-1\right)^{2}+a_{1}A+a_{2}\|c-\bar{x}\|^{2}\ ,$

where $V$ and $A$ are the volume and surface area, $\bar{x}$ is the centroid
and $c$ is the Voronoi site position. The viscosity coefficients (diagonal
entries of $\boldsymbol{\eta}$ in Eq. (16)) of the site degrees of freedom are
very small, such that the system evolves quasi-statically. The larger
viscosity coefficient, corresponding to the volume target $\bar{V}$,
represents the resistance to diffusion of air through the foam interfaces.
Fig. 4 shows the result of a large-scale 3D coarsening simulation beginning
with $n=2000$ monodisperse foam cells. Variation of power diagram weights in
our differentiable Voronoi model allows cells to collapse smoothly to zero
volume. Furthermore, our model implicitly performs the complete topological
restructuring of the local neighborhood that results from each cell collapse.
It should be noted that the geometric accuracy of this simulation, and other
foam simulations involving wide variation of cell sizes and pressures, is
limited due to the lack of curved interfaces.

###  4.3. Biological Tissue Growth

Embryonic cleavage is the process in which a complex organism grows from a
single cell via repeated cell division. It involves an exponential increase in
the number of cells and therefore a rapid increase in the complexity of an
organism. Cell proliferation, a similar process that involves the growth and
subsequent division of cells, is another primary driver of tissue development.
Simulation of such processes is a common use case for computational models
(Conradin et al., 2021; Ghaffarizadeh et al., 2018). Our first experiment
simulates embryonic cleavage in a spherical membrane, beginning with a single
Voronoi cell. By the dynamic coupling formulation described in Sec. 3.5, the
domain boundary is deformable and subject to its own elastic forces, with a
boundary mesh edge of length $\ell$ having energy $E_{\text{B}}=k\ell^{2}$.
The cell energy is

(22)
$E=a_{0}\left(\frac{V}{\bar{V}}-1\right)^{2}+\frac{a_{1}}{\bar{V}^{4/3}}\int_{x\in\mathcal{R}}\|x-\bar{x}\|^{2}dV+\frac{a_{2}}{\bar{V}^{2/3}}\|c-\bar{x}\|^{2},$

where $V$ and $\bar{V}$ are the current and rest volumes, $\bar{x}$ is the
centroid and $c$ is the Voronoi site position. We simulate the evolution of
the system using a momentumless dynamic model as described in Eq. (16). After
every $k$ simulation time steps, each Voronoi site is replaced by two daughter
sites at positions $c\pm\beta\bar{V}^{1/3}\mathbf{n}$, where $\mathbf{n}$ is
the unit normal vector to the cleavage plane. The first three cleavages are
orthogonal, after which the cleavage planes are chosen randomly. The daughter
sites are initialized with the same power diagram weight as the parent site,
and rest volume $\frac{1}{2}\bar{V}_{\text{parent}}$. The system is simulated
over 700 time steps with a division every $T=60$ frames, resulting in a final
cell count of 4096. In Fig. 1, the simulated embryo is compared to frames from
(van Ijken, 2018) depicting embryonic development of a salamander. While this
comparison is not intended to be quantitative, it shows that our method is
able to produce simulations qualitatively similar to real-world cellular
systems undergoing many topological changes.

A second experiment simulates cell proliferation in a cylindrical container of
radius 1, with an added gravity term and free upper surface (Fig. 5). After
each simulation step of size $h$, each cell divides with probability
$p=\alpha\bar{V}h$, where $\alpha$ is the proliferation rate of the tissue.
The daughter sites are initialized with rest volume
$\frac{1}{2}\bar{V}_{\text{parent}}$ increasing linearly to
$\gamma\bar{V}_{\text{parent}}$ over time period $\tau$.

(a) Frame 0, $n=1$.

(b) Frame 100, $n=5$.

(c) Frame 250, $n=24$.

(d) Frame 450, $n=87$.

(e) Frame 736, $n=250$.

(f) Frame 1100, $n=620$.

Figure 5. Cell proliferation in a cylindrical container. The free boundary is
coupled to the cells using a weak elastic membrane model. Random cell
divisions induce frequent and irregular topology changes throughout the
simulation.

###  4.4. Characterization of Foam from Image

We apply our equilibrium-constrained optimization formulation (Sec. 3.8) to
characterize a 2D foam from an image and create a matching computational
model. The image (Fig. 6a) was taken in our lab and depicts a soap foam
between two glass plates. We manually annotate the image, defining the vertex
coordinates and the edges in the desired tessellation. We assume a cell energy
of the form

(23) $E=a_{0}(A-\bar{A})^{2}+a_{1}P,$

where $A$ and $P$ are the cell area and perimeter, and $\bar{A}$ is a target
area. This is analogous to the 3D dry foam energy used in our coarsening
simulation (Sec. 4.2). The model is initialized by placing sites at the
centroids of the cells in the annotated image, and assigning $\bar{A}$ equal
to their area. The system is allowed to converge to equilibrium, resulting in
the pre-optimization model shown in Fig. 6b.

(a) Input image.

(b) Before optimization. Objective value 0.28.

(c) After optimization. Objective value 0.034.

Figure 6. Construction of a Voronoi foam model from an input image using
equilibrium-constrained optimization.

We then perform equilibrium-constrained optimization on the area targets
$\bar{A}$ to find cell pressures matching the observed configuration. The
optimization objective is to minimize the sum of squared distances between the
annotation vertices and the corresponding Voronoi vertices. The corresponding
Voronoi vertex does not need to exist in the tessellation; the generating
expression (the circumcenter of the three neighboring sites, or similar for
boundary vertices) is always computable. The optimization converged in 77
L-BFGS iterations and 13 seconds of computation time. As can be seen from Fig.
6c, the optimization successfully reduces the discrepancy between simulated
and real-world foams. It should be noted that this procedure converges
reliably to a good solution only when the topology of the initial equilibrium
state matches that in the image.

###  4.5. Comparison to Deformable Cell Model

To compare the accuracy and performance of our model with other simulation
methods, we construct a deformable cell model (DCM) using Incremental
Potential Contact (IPC) (Li et al., 2020) to resolve contact between the
cells. The DCM cells have 30 vertices and hence 60 degrees of freedom, while
ours are power diagram cells with 3 degrees of freedom. We use a per-cell
energy of the form

(24) $E=a_{0}(A-\bar{A})^{2}+a_{1}P^{2},$

where $A$ and $P$ are the cell area and perimeter, and $\bar{A}$ is a target
area. We use a combined target of 1.2 times the domain area such that the DCM
cells are in compression and fill the entire space. The DCM uses an additional
quadratic penalty per edge length for regularization as well as the IPC
energy.

(a) Initial state, 1x1 domain.

(b) Final state, 1.50x0.67 domain.

Figure 7. Comparison between implicit Voronoi model (background) and
deformable cell model (white overlay) in a 100-frame quasi-static simulation.
Identical topology changes occur in both models. Cells are colored to show
correspondence between initial and final states.

We perform a quasi-static simulation with 30 cells starting in a square
domain. The initial state for the DCM is constructed by offsetting each
Voronoi cell inward and distributing vertices evenly around the perimeter,
followed by energy minimization to reach equilibrium. In each frame, the
domain is reshaped and both models are allowed to converge to a new
equilibrium. The two models behave similarly, achieving the same final
topology as shown in Fig. 7. However, the DCM required 1000 times more
computation time and nearly 40 times as many average Newton iterations to
converge (see Table 1), demonstrating the efficiency of our implicit Voronoi
model.

###  4.6. Runtime

Runtime statistics for all simulation experiments are collected in Table 1.
Experiments are performed using a workstation with an AMD Ryzen Threadripper
PRO 5995WX CPU.

Table 1. Runtime statistics for experiments. Each frame represents one simulation time step and consists of solving a single optimization problem. Experiment | # Cells | # DOF | # Frames | # Iter / Frame (avg / max) | Time / Iter (avg) [ms] | Total time [s]
---|---|---|---|---|---|---
Comparison (Ours) | 30 | 90 | 100 | 5/25 | 2.32 | 1.23
Comparison (DCM) | 30 | 1800 | 100 | 195/2621 | 65 | 1258
Tissue Growth (membrane) | 4096 | 24070 | 700 | 3/12 | 3830 | 8521
Tissue Growth (cylinder) | 620 | 6067 | 1101 | 10/47 | 1383 | 14970
Coarsening | 2000 | 10000 | 350 | 51/522 | 1275 | 22996
Rigid Body 1 | 2000 | 4003 | 396 | 3/4 | 212 | 251
Rigid Body 2 | 2000 | 4003 | 304 | 3/5 | 213 | 194

The number of iterations per frame depends highly on the number and complexity
of topological changes that occur, resulting in large discrepancies among the
3D experiments. In the embryonic cleavage example, few topological changes
occur except in the frames immediately following each simultaneous cell
division. The coarsening example is particularly expensive because many cells
collapse per frame, and cell collapse results in more complex topological
changes than intercalation. Furthermore, the energy gradient in the coarsening
example is discontinuous across topology changes, slowing the convergence of
Newton’s method, while the other examples use $C^{1}$-continuous energies.

We perform an additional experiment to measure the runtime scaling of our
method with number of cells. In a square (2D) or cube (3D) with $n$ randomly
placed Voronoi sites, we perform a single Newton iteration and measure the
total runtime, along with the runtime of significant subroutines including the
generation of the Voronoi diagram, computation of the energy Hessian and
solution of the linear system (Fig. 8).

Figure 8. Runtime of a Newton iteration vs. number of cells for 2D (left) and
3D (right) models. In the range of problem sizes considered, the Voronoi
diagram generation (including evaluation of derivatives e.g.
$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\mathbf{c}}$) and Hessian computation
exhibit roughly linear scaling, while the cost of solving the linear system
dominates for larger problems.

As can be observed from the results, 3D cases are significantly more expensive
than 2D cases for similar numbers of cells. This is primarily due to the
increase in average number of neighbors for each cell, from 6 in 2D to
$\sim$15.5 (Meijering, 1953) in 3D, leading to a much denser Hessian and a
larger number of Voronoi vertices.

##  5\. Conclusions

We have proposed a novel simulation approach for mechanical cellular systems
based on differentiable Voronoi diagrams. The method successfully addresses
challenges faced by existing models in handling topological transitions, while
using an extremely compact state representation that implicitly defines the
shape and topology of the interface network. We demonstrated this using a
diverse set of examples, including simulations of complex biological and
physical processes as well as comparisons to alternative models and real-world
systems.

###  5.1. Limitations & Future Work

The development of informative biological simulations will require more
complex energy models and the incorporation of real-world data to determine
accurate parameters. Deriving a smoother approximation to the surface area of
Voronoi cells would enable further applications, including adhesion effects at
biological cell interfaces. Our model currently uses an isotropic distance
metric that restricts the space of possible cell shapes, and the method
assumes planar or piecewise linear cell interfaces. However, epithelial
tissues, for example, exhibit non-convex scutoid-shaped cells (Gómez-Gálvez et
al., 2018) which permit different cell packing patterns. Extending our method
towards more complex cell shapes while retaining a compact site-based
representation is an interesting challenge for future work. Voronoi diagrams
are also a natural design space for foam-like 3D-printed metamaterials. This
direction has been explored, _e.g_., by Martinez _et al_. (2018), with
extensions to orthotropic foams based on Voronoi diagrams with spatially
varying distance metrics (Martínez et al., 2019). An extension of our approach
may permit gradient-based optimization of these materials. Finally, we have
investigated a simple instance of intercellular navigation, where a rigid
object is driven by a constant force. Designing control algorithms to guide an
object through cellular systems towards a desired target is another exciting
future direction.

###### Acknowledgements.

We thank the anonymous reviewers for their valuable feedback. This work was
supported by the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program (grant agreement No. 866480), and
the Swiss National Science Foundation through SNF project grant 200021_200644.

## References

  * (1)
  * Alt et al. (2017) Silvanus Alt, Poulami Ganguly, and Guillaume Salbreux. 2017.  Vertex models: From cell mechanics to tissue morphogenesis.  _Philosophical Transactions of the Royal Society B: Biological Sciences_ 372 (05 2017), 20150520.  https://doi.org/10.1098/rstb.2015.0520
  * Bargteil et al. (2007) Adam W. Bargteil, Chris Wojtan, Jessica K. Hodgins, and Greg Turk. 2007.  A finite element method for animating large viscoplastic flow. In _ACM SIGGRAPH 2007 Papers_ (San Diego, California) _(SIGGRAPH ’07)_. Association for Computing Machinery, New York, NY, USA, 16–es.  https://doi.org/10.1145/1275808.1276397
  * Bickel et al. (2009) Bernd Bickel, Moritz Bächer, Miguel A. Otaduy, Wojciech Matusik, Hanspeter Pfister, and Markus Gross. 2009.  Capture and modeling of non-linear heterogeneous soft tissue.  _ACM Trans. Graph._ 28, 3, Article 89 (jul 2009), 9 pages.  https://doi.org/10.1145/1531326.1531395
  * Bogosel and Oudet (2022) Beniamin Bogosel and Edouard Oudet. 2022.  Longest minimal length partitions.  _Interfaces and Free Boundaries_ 24, 1 (Jan. 2022), 95–135.  https://doi.org/10.4171/IFB/468
  * Brochu and Bridson (2009) Tyson Brochu and Robert Bridson. 2009.  Robust Topological Operations for Dynamic Explicit Surfaces.  _SIAM J. Sci. Comput._ 31, 4 (jun 2009), 2472–2493.  https://doi.org/10.1137/080737617
  * Brodland et al. (2006) G. Wayne Brodland, Daniel I-Li Chen, and Jim H. Veldhuis. 2006.  A cell-based constitutive model for embryonic epithelia and other planar aggregates of biological cells.  _International Journal of Plasticity_ 22, 6 (2006), 965–995.  https://doi.org/10.1016/j.ijplas.2005.05.002
  * Busaryev et al. (2012) Oleksiy Busaryev, Tamal K. Dey, Huamin Wang, and Zhong Ren. 2012.  Animating Bubble Interactions in a Liquid Foam.  _ACM Trans. Graph._ 31, 4, Article 63 (jul 2012), 8 pages.  https://doi.org/10.1145/2185520.2185559
  * Chen et al. (2014) Xiang Chen, Changxi Zheng, Weiwei Xu, and Kun Zhou. 2014.  An asymptotic numerical method for inverse elastic shape design.  _ACM Trans. Graph._ 33, 4, Article 95 (jul 2014), 11 pages.  https://doi.org/10.1145/2601097.2601189
  * Chen et al. (2008) Yanqing Chen, Timothy A Davis, William W Hager, and Sivasankaran Rajamanickam. 2008.  Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate.  _ACM Transactions on Mathematical Software (TOMS)_ 35, 3 (2008), 1–14.
  * Conradin et al. (2021) Raphaël Conradin, Christophe Coreixas, Jonas Latt, and Bastien Chopard. 2021.  PalaCell2D: A framework for detailed tissue morphogenesis.  _Journal of Computational Science_ 53 (2021), 101353.  https://doi.org/10.1016/j.jocs.2021.101353
  * Da et al. (2015) Fang Da, Christopher Batty, Chris Wojtan, and Eitan Grinspun. 2015.  Double Bubbles Sans Toil and Trouble: Discrete Circulation-Preserving Vortex Sheets for Soap Films and Foams.  _ACM Trans. on Graphics (SIGGRAPH North America 2015)_ (2015).
  * de Goes et al. (2015) Fernando de Goes, Corentin Wallez, Jin Huang, Dmitry Pavlov, and Mathieu Desbrun. 2015.  Power Particles: An Incompressible Fluid Solver Based on Power Diagrams.  _ACM Trans. Graph._ 34, 4, Article 50 (jul 2015), 11 pages.  https://doi.org/10.1145/2766901
  * Deng et al. (2022) Yitong Deng, Mengdi Wang, Xiangxin Kong, Shiying Xiong, Zangyueyang Xian, and Bo Zhu. 2022.  A moving eulerian-lagrangian particle method for thin film and foam simulation.  _ACM Trans. Graph._ 41, 4, Article 154 (jul 2022), 17 pages.  https://doi.org/10.1145/3528223.3530174
  * Eppstein (2012) David Eppstein. 2012.  The Graphs of Planar Soap Bubbles.  https://doi.org/10.48550/ARXIV.1207.3761
  * Farhadifar et al. (2007) Reza Farhadifar, Jens-Christian Röper, Benoit Aigouy, Suzanne Eaton, and Frank Jülicher. 2007.  The Influence of Cell Mechanics, Cell-Cell Interactions, and Proliferation on Epithelial Packing.  _Current Biology_ 17, 24 (2007), 2095–2104.  https://doi.org/10.1016/j.cub.2007.11.049
  * Feng et al. (2022) Fan Feng, Shiying Xiong, Ziyue Liu, Zangyueyang Xian, Yuqing Zhou, Hiroki Kobayashi, Atsushi Kawamoto, Tsuyoshi Nomura, and Bo Zhu. 2022.  Cellular topology optimization on differentiable Voronoi diagrams.  _Internat. J. Numer. Methods Engrg._ 124, 1 (sep 2022), 282–304.  https://doi.org/10.1002/nme.7121
  * Ferguson et al. (2023) Zachary Ferguson, Teseo Schneider, Danny Kaufman, and Daniele Panozzo. 2023.  In-Timestep Remeshing for Contacting Elastodynamics.  _ACM Trans. Graph._ 42, 4, Article 145 (jul 2023), 15 pages.  https://doi.org/10.1145/3592428
  * Fletcher et al. (2014) Alexander G. Fletcher, Miriam Osterfield, Ruth E. Baker, and Stanislav Y. Shvartsman. 2014.  Vertex Models of Epithelial Morphogenesis.  _Biophysical Journal_ 106, 11 (2014), 2291–2304.  https://doi.org/10.1016/j.bpj.2013.11.4498
  * Geilinger et al. (2020) Moritz Geilinger, David Hahn, Jonas Zehnder, Moritz Bächer, Bernhard Thomaszewski, and Stelian Coros. 2020.  Add: Analytically differentiable dynamics for multi-body systems with frictional contact.  _ACM Transactions on Graphics (TOG)_ 39, 6 (2020), 1–15.
  * Ghaffarizadeh et al. (2018) Ahmadreza Ghaffarizadeh, Randy Heiland, Samuel H. Friedman, Shannon M. Mumenthaler, and Paul Macklin. 2018.  PhysiCell: An open source physics-based cell simulator for 3-D multicellular systems.  _PLOS Computational Biology_ 14, 2 (02 2018), 1–31.  https://doi.org/10.1371/journal.pcbi.1005991
  * Goldade et al. (2020) Ryan Goldade, Mridul Aanjaneya, and Christopher Batty. 2020.  Constraint bubbles and affine regions: reduced fluid models for efficient immersed bubbles and flexible spatial coarsening.  _ACM Trans. Graph._ 39, 4, Article 43 (aug 2020), 15 pages.  https://doi.org/10.1145/3386569.3392455
  * Gómez-Gálvez et al. (2018) Pedro Gómez-Gálvez, Pablo Vicente-Munuera, Antonio Tagua, Cristina Forja, Ana M Castro, Marta Letrán, Andrea Valencia-Expósito, Clara Grima, Marina Bermúdez-Gallardo, Óscar Serrano-Pérez-Higueras, et al. 2018\.  Scutoids are a geometrical solution to three-dimensional packing of epithelia.  _Nature communications_ 9, 1 (2018), 2960.
  * Guennebaud et al. (2010) Gaël Guennebaud, Benoît Jacob, et al. 2010\.  Eigen v3.  http://eigen.tuxfamily.org.
  * Hahn et al. (2019) David Hahn, Pol Banzet, James M Bern, and Stelian Coros. 2019.  Real2sim: Visco-elastic parameter estimation from dynamic motion.  _ACM Transactions on Graphics (TOG)_ 38, 6 (2019), 1–13.
  * Honda et al. (2004) Hisao Honda, Masaharu Tanemura, and Tatsuzo Nagai. 2004.  A three-dimensional vertex dynamics cell model of space-filling polyhedra simulating cell behavior in a cell aggregate.  _Journal of Theoretical Biology_ 226, 4 (2004), 439–453.  https://doi.org/10.1016/j.jtbi.2003.10.001
  * Hong et al. (2008) Jeong-Mo Hong, Ho-Young Lee, Jong-Chul Yoon, and Chang-Hun Kim. 2008.  Bubbles alive.  _ACM Trans. Graph._ 27, 3 (aug 2008), 1–4.  https://doi.org/10.1145/1360612.1360647
  * Hu et al. (2019) Yuanming Hu, Jiancheng Liu, Andrew Spielberg, Joshua B Tenenbaum, William T Freeman, Jiajun Wu, Daniela Rus, and Wojciech Matusik. 2019.  Chainqueen: A real-time differentiable physical simulator for soft robotics. In _2019 International conference on robotics and automation (ICRA)_. IEEE, 6265–6271.
  * Huang et al. (2020) Weizhen Huang, Julian Iseringhausen, Tom Kneiphof, Ziyin Qu, Chenfanfu Jiang, and Matthias B. Hullin. 2020.  Chemomechanical simulation of soap film flow on spherical bubbles.  _ACM Trans. Graph._ 39, 4, Article 41 (aug 2020), 14 pages.  https://doi.org/10.1145/3386569.3392094
  * Ishida et al. (2020) Sadashige Ishida, Peter Synak, Fumiya Narita, Toshiya Hachisuka, and Chris Wojtan. 2020.  A model for soap film dynamics with evolving thickness.  _ACM Trans. Graph._ 39, 4, Article 31 (aug 2020), 11 pages.  https://doi.org/10.1145/3386569.3392405
  * Kim et al. (2021) Sangwoo Kim, Marie Pochitaloff, Georgina Stooke-Vaughan, and Otger Campàs. 2021.  Embryonic Tissues as Active Foams.  _Nature Physics_ 17 (07 2021), 1–8.  https://doi.org/10.1038/s41567-021-01215-1
  * Kondo et al. (2018) Yohei Kondo, Kazuhiro Aoki, and Shin Ishii. 2018.  Inverse tissue mechanics of cell monolayer expansion.  _PLOS Computational Biology_ 14 (03 2018), e1006029.  https://doi.org/10.1371/journal.pcbi.1006029
  * Koschier et al. (2014) Dan Koschier, Sebastian Lipponer, and Jan Bender. 2014.  Adaptive Tetrahedral Meshes for Brittle Fracture Simulation.. In _Symposium on Computer Animation_ , Vol. 3.
  * Lévy (2015) Bruno Lévy. 2015.  Geogram.  https://github.com/BrunoLevy/geogram.
  * Li et al. (2020) Minchen Li, Zachary Ferguson, Teseo Schneider, Timothy Langlois, Denis Zorin, Daniele Panozzo, Chenfanfu Jiang, and Danny M. Kaufman. 2020.  Incremental Potential Contact: Intersection-and Inversion-Free, Large-Deformation Dynamics.  _ACM Trans. Graph._ 39, 4, Article 49 (aug 2020), 20 pages.  https://doi.org/10.1145/3386569.3392425
  * Liu et al. (2009) Yang Liu, Wenping Wang, Bruno Lévy, Feng Sun, Dong-Ming Yan, Lin Lu, and Chenglei Yang. 2009.  On Centroidal Voronoi Tessellation—Energy Smoothness and Fast Computation.  _ACM Trans. Graph._ 28, 4, Article 101 (sep 2009), 17 pages.  https://doi.org/10.1145/1559755.1559758
  * Lloyd (1982) S. Lloyd. 1982.  Least squares quantization in PCM.  _IEEE Transactions on Information Theory_ 28, 2 (1982), 129–137.  https://doi.org/10.1109/TIT.1982.1056489
  * Lumpe et al. (2022) Thomas S Lumpe, Michael Tao, Kristina Shea, and David I W Levin. 2022.  Computational design and fabrication of active 3D-printed multi-state structures for shape morphing.  _Smart Materials and Structures_ 32, 1 (dec 2022), 015008.  https://doi.org/10.1088/1361-665X/aca5d6
  * Ly et al. (2018) Mickaël Ly, Romain Casati, Florence Bertails-Descoubes, Mélina Skouras, and Laurence Boissieux. 2018.  Inverse elastic shell design with contact and friction.  _ACM Trans. Graph._ 37, 6, Article 201 (dec 2018), 16 pages.  https://doi.org/10.1145/3272127.3275036
  * Martínez et al. (2018) Jonàs Martínez, Samuel Hornus, Haichuan Song, and Sylvain Lefebvre. 2018.  Polyhedral Voronoi Diagrams for Additive Manufacturing.  _ACM Trans. Graph._ 37, 4, Article 129 (jul 2018), 15 pages.  https://doi.org/10.1145/3197517.3201343
  * Martínez et al. (2019) Jonàs Martínez, Mélina Skouras, Christian Schumacher, Samuel Hornus, Sylvain Lefebvre, and Bernhard Thomaszewski. 2019.  Star-Shaped Metrics for Mechanical Metamaterial Design.  _ACM Trans. Graph._ 38, 4, Article 82 (jul 2019), 13 pages.  https://doi.org/10.1145/3306346.3322989
  * Meijering (1953) J. L. Meijering. 1953.  Interface area, edge length, and number of vertices in crystal aggregates with random nucleation.  _Philips Research Reports_ 8 (1953), 270–290.
  * Miguel et al. (2012) E. Miguel, D. Bradley, B. Thomaszewski, B. Bickel, W. Matusik, M. A. Otaduy, and S. Marschner. 2012.  Data-Driven Estimation of Cloth Simulation Models.  _Comput. Graph. Forum_ 31, 2pt2 (may 2012), 519–528.
  * Misztal et al. (2014) Marek K. Misztal, Kenny Erleben, Adam Bargteil, Jens Fursund, Brian Bunch Christensen, J. Andreas Bærentzen, and Robert Bridson. 2014.  Multiphase Flow of Immiscible Fluids on Unstructured Moving Meshes.  _IEEE Transactions on Visualization and Computer Graphics_ 20, 1 (jan 2014), 4–16.  https://doi.org/10.1109/TVCG.2013.97
  * Montes et al. (2020) Juan Montes, Bernhard Thomaszewski, Sudhir Mudur, and Tiberiu Popa. 2020.  Computational design of skintight clothing.  _ACM Trans. Graph._ 39, 4, Article 105 (aug 2020), 12 pages.  https://doi.org/10.1145/3386569.3392477
  * Moukarzel (1997) C. Moukarzel. 1997.  Geometrical consequences of foam equilibrium.  _Phys. Rev. E_ 55 (Jun 1997), 6866–6880.  Issue 6. https://doi.org/10.1103/PhysRevE.55.6866
  * Narain et al. (2013) Rahul Narain, Tobias Pfaff, and James F. O’Brien. 2013.  Folding and crumpling adaptive sheets.  _ACM Trans. Graph._ 32, 4, Article 51 (jul 2013), 8 pages.  https://doi.org/10.1145/2461912.2462010
  * Narain et al. (2012) Rahul Narain, Armin Samii, and James F. O’Brien. 2012.  Adaptive anisotropic remeshing for cloth simulation.  _ACM Trans. Graph._ 31, 6, Article 152 (nov 2012), 10 pages.  https://doi.org/10.1145/2366145.2366171
  * O’Brien et al. (2002) James F. O’Brien, Adam W. Bargteil, and Jessica K. Hodgins. 2002.  Graphical modeling and animation of ductile fracture.  _ACM Trans. Graph._ 21, 3 (jul 2002), 291–294.  https://doi.org/10.1145/566654.566579
  * Okuda et al. (2012) Satoru Okuda, Yasuhiro Inoue, Mototsugu Eiraku, Yoshiki Sasai, and Taiji Adachi. 2012.  Reversible network reconnection model for simulating large deformation in dynamic tissue morphogenesis.  _Biomechanics and modeling in mechanobiology_ 12 (09 2012).  https://doi.org/10.1007/s10237-012-0430-7
  * Panetta et al. (2021) Julian Panetta, Florin Isvoranu, Tian Chen, Emmanuel Siéfert, Benoît Roman, and Mark Pauly. 2021.  Computational Inverse Design of Surface-Based Inflatables.  _ACM Trans. Graph._ 40, 4, Article 40 (July 2021), 14 pages.  https://doi.org/10.1145/3450626.3459789
  * Panetta et al. (2017) Julian Panetta, Abtin Rahimian, and Denis Zorin. 2017.  Worst-case stress relief for microstructures.  _ACM Trans. Graph._ 36, 4, Article 122 (jul 2017), 16 pages.  https://doi.org/10.1145/3072959.3073649
  * Pfaff et al. (2014) Tobias Pfaff, Rahul Narain, Juan Miguel de Joya, and James F. O’Brien. 2014.  Adaptive tearing and cracking of thin sheets.  _ACM Trans. Graph._ 33, 4, Article 110 (jul 2014), 9 pages.  https://doi.org/10.1145/2601097.2601132
  * Piatnitski and Ptashnyk (2020) Andrey Piatnitski and Mariya Ptashnyk. 2020.  Homogenization of biomechanical models of plant tissues with randomly distributed cells.  _Nonlinearity_ 33 (09 2020), 5510–5542.  https://doi.org/10.1088/1361-6544/ab95ab
  * Schreck et al. (2016) Camille Schreck, Damien Rohmer, Stefanie Hahmann, Marie-Paule Cani, Shuo Jin, Charlie C. L. Wang, and Jean-Francis Bloch. 2016.  Nonsmooth Developable Geometry for Interactively Animating Paper Crumpling.  _ACM Trans. Graph._ 35, 1, Article 10 (dec 2016), 18 pages.  https://doi.org/10.1145/2829948
  * Sharp et al. (2019) Nicholas Sharp et al. 2019\.  Polyscope.  www.polyscope.run.
  * Skouras et al. (2014) Mélina Skouras, Bernhard Thomaszewski, Peter Kaufmann, Akash Garg, Bernd Bickel, Eitan Grinspun, and Markus Gross. 2014.  Designing inflatable structures.  _ACM Trans. Graph._ 33, 4, Article 63 (jul 2014), 10 pages.  https://doi.org/10.1145/2601097.2601166
  * Thomas et al. (2015) Gilberto Thomas, Julio Belmonte, François Graner, James Glazier, and Rita de Almeida. 2015.  3D simulations of wet foam coarsening evidence a self similar growth regime.  _Colloids and Surfaces A: Physicochemical and Engineering Aspects_ 473 (02 2015).  https://doi.org/10.1016/j.colsurfa.2015.02.015
  * Umetani et al. (2011) Nobuyuki Umetani, Danny M. Kaufman, Takeo Igarashi, and Eitan Grinspun. 2011.  Sensitive couture for interactive garment modeling and editing. In _ACM SIGGRAPH 2011 Papers_ (Vancouver, British Columbia, Canada) _(SIGGRAPH ’11)_. Association for Computing Machinery, New York, NY, USA, Article 90, 12 pages.  https://doi.org/10.1145/1964921.1964985
  * van Ijken (2018) Jan van Ijken. 2018.  Becoming.  https://www.janvanijken.com/films/becoming/
  * Vedel-Larsen (2010) Bue Krogh Vedel-Larsen. 2010.  _Quasi-static Simulation of Foam in the Vertex Model_.  Master’s thesis. University of Copenhagen, Denmark Department of Computer Science.
  * Wang (2018) Huamin Wang. 2018.  Rule-free sewing pattern adjustment with precision and efficiency.  _ACM Trans. Graph._ 37, 4, Article 53 (jul 2018), 13 pages.  https://doi.org/10.1145/3197517.3201320
  * Wang et al. (2015) Xiaoning Wang, Xiang Ying, Yong-Jin Liu, Shi-Qing Xin, Wenping Wang, Xianfeng Gu, Wolfgang Mueller-Wittig, and Ying He. 2015.  Intrinsic computation of centroidal Voronoi tessellation (CVT) on meshes.  _Computer-Aided Design_ 58 (2015), 51–61.
  * Wicke et al. (2010) Martin Wicke, Daniel Ritchie, Bryan M. Klingner, Sebastian Burke, Jonathan R. Shewchuk, and James F. O’Brien. 2010.  Dynamic local remeshing for elastoplastic simulation.  _ACM Trans. Graph._ 29, 4, Article 49 (jul 2010), 11 pages.  https://doi.org/10.1145/1778765.1778786
  * Williams et al. (2020) Francis Williams, Jerome Parent-Levesque, Derek Nowrouzezahrai, Daniele Panozzo, Kwang Moo Yi, and Andrea Tagliasacchi. 2020.  Voronoinet: General functional approximators with local support. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops_. 264–265.
  * Wojtan et al. (2009) Chris Wojtan, Nils Thürey, Markus Gross, and Greg Turk. 2009.  Deforming meshes that split and merge.  In _ACM SIGGRAPH 2009 papers_. 1–10.
  * Yan et al. (2009) Dong‐Ming Yan, Bruno Levy, Yang Liu, Feng Sun, and Wenping Wang. 2009.  Isotropic Remeshing with Fast and Exact Computation of Restricted Voronoi Diagram.  _Comput. Graph. Forum_ 28 (07 2009), 1445–1454.  https://doi.org/10.1111/j.1467-8659.2009.01521.x
  * Yan et al. (2011) Dong‐Ming Yan, Wenping Wang, Bruno Levy, and Yang Liu. 2011.  Efficient Computation of Clipped Voronoi Diagram for Mesh Generation.  _Computer-Aided Design_ 45 (09 2011).  https://doi.org/10.1016/j.cad.2011.09.004
  * Yan et al. (2018) Guowei Yan, Wei Li, Ruigang Yang, and Huamin Wang. 2018.  Inexact descent methods for elastic parameter optimization.  _ACM Transactions on Graphics (TOG)_ 37, 6 (2018), 1–14.
  * Yvinec (2023) Mariette Yvinec. 2023.  2D Triangulations.  In _CGAL User and Reference Manual_. CGAL Editorial Board.  https://doc.cgal.org/5.5.2/Manual/packages.html#PkgTriangulation2
  * Zhang et al. (2022) Jiayi Eris Zhang, Jérémie Dumas, Yun (Raymond) Fei, Alec Jacobson, Doug L. James, and Danny M. Kaufman. 2022.  Progressive Simulation for Cloth Quasistatics.  _ACM Trans. Graph._ 41, 6, Article 218 (nov 2022), 16 pages.  https://doi.org/10.1145/3550454.3555510
  * Zhang et al. (2017) Xiaoting Zhang, Guoxin Fang, Chengkai Dai, Jouke Verlinden, Jun Wu, Emily Whiting, and Charlie C.L. Wang. 2017.  Thermal-Comfort Design of Personalized Casts. In _Proceedings of the 30th Annual ACM Symposium on User Interface Software and Technology_ (Québec City, QC, Canada) _(UIST ’17)_. Association for Computing Machinery, New York, NY, USA, 243–254.  https://doi.org/10.1145/3126594.3126600

##  Appendix A Volume integration in tetrahedra

Here we provide a typical approach for volume integration over a tetrahedron,
as used for analytically evaluating certain cell energy functions in this
work. We consider a tetrahedron with vertices ($x_{i}$, $y_{i}$, $z_{i}$),
$i=0,1,2$ and the origin. The tetrahedron can be parameterized with
coordinates ($u$, $v$, $w$) $\in$ [0,1] such that

(25) $\displaystyle x=u(vx_{0}+(1-v)(wx_{1}+(1-w)x_{2})),$ $\displaystyle
y=u(vy_{0}+(1-v)(wy_{1}+(1-w)y_{2})),$ $\displaystyle
z=u(vz_{0}+(1-v)(wz_{1}+(1-w)z_{2})).$

The Jacobian matrix of this transformation is

(26) $J=\begin{bmatrix}\frac{\partial x}{\partial u}&\frac{\partial
x}{\partial v}&\frac{\partial x}{\partial w}\\\ \frac{\partial y}{\partial
u}&\frac{\partial y}{\partial v}&\frac{\partial y}{\partial w}\\\
\frac{\partial z}{\partial u}&\frac{\partial z}{\partial v}&\frac{\partial
z}{\partial w}\end{bmatrix}$

where the partial derivatives are

(27) $\displaystyle\frac{\partial x}{\partial u}=$ $\displaystyle
vx_{0}+(1-v)(wx_{1}+(1-w)x_{2}),$ $\displaystyle\frac{\partial x}{\partial
v}=$ $\displaystyle u(x_{0}-wx_{1}-(1-w)x_{2}),$ $\displaystyle\frac{\partial
x}{\partial w}=$ $\displaystyle u(1-v)(x_{1}-x_{2}).$

A volume integral can be transformed as follows, where $f^{*}$ is obtained by
substituting equations 25 for ($x,y,z$) in $f$:

(28) $\iiint
f(x,y,z)dV=\int\limits_{0}^{1}\\!\\!\\!\\!\int\limits_{0}^{1}\\!\\!\\!\\!\int\limits_{0}^{1}\det(J)\hskip
2.84526ptf^{*}(u,v,w)\hskip 2.84526ptdudvdw$

##  Appendix B Hessian for Boundary Coupling

Taking the derivative of (9), we obtain the second order terms

(29) $\begin{gathered}\begin{aligned}
\frac{\mathrm{d}^{2}E}{\mathrm{d}\mathbf{c}\mathrm{d}\mathbf{p}}&=\left(\frac{\partial\mathbf{x}}{\partial\mathbf{c}}\right)^{\top}\frac{\partial^{2}F}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathbf{p}}+\left(\frac{\partial\mathbf{x}}{\partial\mathbf{c}}\right)^{\top}\frac{\partial^{2}F}{\partial\mathbf{x}\partial\mathbf{p}}\\\
&+\frac{\partial^{2}F}{\partial\mathbf{c}\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathbf{p}}+\left(\sum_{i}\frac{\partial
F}{\partial\mathbf{x}_{i}}\frac{\partial^{2}\mathbf{x}_{i}}{\partial\mathbf{c}\partial\mathbf{v}}\right)\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathbf{p}}+\frac{\partial^{2}F}{\partial\mathbf{c}\partial\mathbf{p}},\end{aligned}\\\
\begin{aligned}
\frac{\mathrm{d}^{2}E}{\mathrm{d}\mathbf{p}^{2}}&=\left(\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathbf{p}}\right)^{\top}\frac{\partial^{2}F}{\partial\mathbf{x}^{2}}\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathbf{p}}+\left(\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathbf{p}}\right)^{\top}\frac{\partial^{2}F}{\partial\mathbf{x}\partial\mathbf{p}}\\\
&+\frac{\partial^{2}F}{\partial\mathbf{p}\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial\mathbf{v}}\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathbf{p}}+\left(\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathbf{p}}\right)^{\top}\left(\sum_{i}\frac{\partial
F}{\partial\mathbf{x}_{i}}\frac{\partial^{2}\mathbf{x}_{i}}{\partial\mathbf{v}^{2}}\right)\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathbf{p}}\\\
&+\left(\sum_{i}\left(\frac{\partial
F}{\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial\mathbf{v}_{i}}\right)\frac{\mathrm{d}^{2}\mathbf{v}_{i}}{\mathrm{d}\mathbf{p}^{2}}\right)+\frac{\mathrm{d}^{2}F_{B}}{\mathrm{d}\mathbf{p}^{2}}\
,\end{aligned}\\\ \end{gathered}$

which must be added to the Hessian (7) of the fixed-boundary energy.

##  Appendix C Dynamic simulation convergence analysis

Convergence of simulations under refinement of the time step is important for
assessing the accuracy of numerical solutions. To evaluate the convergence of
dynamic simulations using our method, we perform a simple 2D simulation with
four cells and one topological transition, illustrated in Fig. 9.

Figure 9. Approximate initial (left) and final (right) states of a dynamic
simulation used to assess the convergence of the implicit Voronoi model.

The simulation is repeated for the same total time using a range of time
steps. Error measurements in Fig. 10 show that given a sufficiently smooth
energy function and using second-order finite difference approximations
(BDF2), the dynamic simulation method presented in section 3.7 converges
quadratically. Linear convergence is achieved for $C^{0}$-continuous energy
functions (such as perimeter and surface area) or using only a first-order
accurate acceleration.

Figure 10. Final state error vs. time step for a simple 2D case including a
topological transition. Final state from highest resolution simulation is
taken as ground truth $\bar{\mathbf{y}}$. Solid lines represent exact linear
and quadratic convergence, with steeper slope indicating higher-order
convergence.

##  Appendix D Experiment Details

In this section we provide additional details and parameters required to
reproduce the examples in Sec. 4.

#### Intercellular Navigation with Rigid Body

The domain is rectangular with width 4m and height 2m, containing $n=2000$
cells along with the rigid body. Other parameters are $\bar{A}=0.004$ m^2,
$a_{0}=0.05$ N/m^3, $a_{1}=100.0$ N/m^3, $a_{2}=1.0$ N/m, $F_{x}=0.035$ N,
$h=0.01$ s, $m_{\mathbf{c}}=0.0003$ kg, $m_{\mathbf{p}}=0.0003$ kg,
$\eta_{\mathbf{c}}=0.0003$ N.s/m, $\eta_{\mathbf{p}}=0.003$ N.s/m. Here and in
the following examples, $m$ and $\eta$ are diagonal entries of the mass matrix
$\mathbf{M}$ and viscosity matrix $\boldsymbol{\eta}$ corresponding to the
site ($\mathbf{c}$) and boundary ($\mathbf{p}$) degrees of freedom.

#### Foam Coarsening

The domain is a 2m$\times$2m$\times$20cm box starting with $n=2000$ cells.
Other parameters are $\bar{V}_{\text{initial}}=0.0004$ m^3, $a_{0}=0.002$ N.m,
$a_{1}=0.01$ N/m, $a_{2}=0.1$ N/m, $\eta_{\mathbf{c}}=0.0001$ N.s/m,
$\eta_{\bar{V}}=1.0$ N.s/m^3. The time step $h=\frac{1.0}{n}$ s is adapted
each frame to the number of remaining cells $n$, increasing as the rate of
cell collapse slows.

#### Embryonic Cleavage

The initial domain is an icosphere with 4 subdivisions, approximating a sphere
with radius 1.34m. Other parameters are $\bar{V}_{\text{initial}}=10.25$ m^3,
$a_{0}=50.0$ N.m, $a_{1}=5.0$ N, $a_{2}=1.0$ N.m, $h=0.01$ s,
$\eta_{\mathbf{c}}=5.0$ N.s/m, $\eta_{\mathbf{p}}=0.5$ N.s/m, $\beta=0.1$.

#### Tissue Proliferation in Cylinder

The initial domain is a cylinder of radius 1m whose upper deformable face
comprises 2432 triangles. Other parameters are
$\bar{V}_{\text{initial}}=0.113$ m^3, $a_{0}=10.0$ N.m, $a_{1}=10.0$ N,
$a_{2}=1.0$ N.m, $h=0.01$ s, $\eta_{\mathbf{c}}=0.5$ N.s/m,
$\eta_{\mathbf{p}}=0.003$ N.s/m, $\alpha=30$ m^-3.s^-1, $\beta=0.1$,
$\gamma=0.75$, $\tau=0.11$ s.