# Geometric signatures of tissue surface tension in a three-dimensional model
of confluent tissue

Preeti Sahu1,2,∗, J. M. Schwarz1,3 and M. Lisa Manning1,
<EMAIL_ADDRESS><EMAIL_ADDRESS>1 Department of Physics and
BioInspired Syracuse: Institute for Material and Living Systems, Syracuse
University, Syracuse, NY 13244, USA
2 IST Austria, Am Campus 1, 3400 Klosterneuburg, Austria
3 Indian Creek Farm, Ithaca, NY 14850, USA

###### Abstract

In dense biological tissues, cell types performing different roles remain
segregated by maintaining sharp interfaces. To better understand the
mechanisms for such sharp compartmentalization, we study the effect of an
imposed heterotypic tension at the interface between two distinct cell types
in a fully 3D model for confluent tissues. We find that cells rapidly sort and
self-organize to generate a tissue-scale interface between cell types, and
cells adjacent to this interface exhibit signature geometric features
including nematic-like ordering, bimodal facet areas, and registration, or
alignment, of cell centers on either side of the two-tissue interface. The
magnitude of these features scales directly with the magnitude of imposed
tension, suggesting that biologists can estimate the magnitude of tissue
surface tension between two tissue types simply by segmenting a 3D tissue. To
uncover the underlying physical mechanisms driving these geometric features,
we develop two minimal, ordered models using two different underlying lattices
that identify an energetic competition between bulk cell shapes and tissue
interface area. When the interface area dominates, changes to neighbor
topology are costly and occur less frequently, which generates the observed
geometric features.

## Introduction

An important collective phenomenon observed in groups of biological cells is
the process of cell sorting, where cells of different types spontaneously
spatially segregate into separate compartments. These distinct compartments
not only play an integral role during an organism’s formative stages Trinkaus
and Groves (1955); Waites et al. (2017); Unbekandt and Davies (2010); Montero
et al. (2005); Klopper et al. (2010), but are also crucial for the upkeep of
normal functioning of organs Rübsam et al. (2017); Cochet-Escartin et al.
(2017) and for containment of the spread of diseased/infected tissues Friedl
et al. (2004); Foty and Steinberg (2005); Pawlizak et al. (2015); Song et al.
(2016). Broadly speaking, cell sorting mechanisms can be classified into two
generic categories- (a) biochemical/morphogen gradients Turing (1990);
Streichan et al. (2018) and (b) differences in mechanical properties of
individual cells. These mechanical properties can include cell-cell adhesivity
Steinberg (1963, 2007), acto-myosin contractility Harris (1976); Brodland
(2002); Mertz et al. (2012); Maître et al. (2012), a mechano-chemical coupling
between both cell-cell adhesivity and acto-myosin contractility Amack and
Manning (2012); Manning et al. (2010a); Engl et al. (2014), or an explicit
interfacial tension between unlike cells, often called heterotypic interfacial
tension (HIT) Graner and Glazier (1992); Barton et al. (2017); Canty et al.
(2017); Sussman et al. (2018).

Much of the computational and theoretical work on cell sorting has focused on
particle mixture simulations. In such mixtures, the mechanism for sorting
relies heavily on active fluctuations Flenner et al. (2012); Merks and Glazier
(2005); Sun and Wang (2013); Jiang et al. (1998); Palsson (2008). However, an
essential feature of cell sorting that is observed in experimental co-cultures
is that the interface is much sharper than what is expected from a particulate
mixture Landsberg et al. (2009); Manning et al. (2010b); Dahmann et al.
(2011); Nnetu et al. (2012); Monier et al. (2011); Calzolari et al. (2014).
While such a straight and sharp interface is difficult to obtain merely by
diffusive morphogens, heterotypic interfacial tensions in confluent tissues,
where there are no gaps between cells, can generate a sharp interface easily
Kesavan et al. (2020); Sussman et al. (2018); Sahu et al. (2020a). In such
cases, the fact that a confluent monolayer must tessellate space, which is
captured in vertex or Voronoi models for tissues, results in forces that are
discontinuous functions of cell displacements. It is precisely this non-
analytic behaviour resulting from topological interactions between cells that
drives sharpening in two dimensions. And yet, the tissue remains fluid-like
such that the sharp interface is also easily deformable. Such sharp but
deformable interfaces are not observed in particle-based models with metric
interactions between cells.

While confluent monolayers, with a single layer of cells, are biologically
relevant, fully three-dimensional confluent tissues, such as stratified
epithelia or early vertebrate embryos, are even more ubiquitous. Therefore, it
is important to determine if/how the topological nature of the interactions
between the cells also drive sharp but deformable interfaces in three
dimensions. Prior work on topological models for cell sorting in three
dimensions focused on cells coalescing into small clusters via a Rayleigh-
Plateau instability as well as regions of mixed cell types untangling to
facilitate compartmentalization in the absence of fluctuations Hutson et al.
(2008). In this manuscript we instead focus on quantifying the dynamics and
cell geometries in the presence of fluctuations, and identifying the
topological mechanisms driving cell sorting in three dimensions.

One goal is to study the dynamics of cell sorting in confluent fluid-like
tissues by implementing HIT in the presence of fluctuating forces in 3D.
Indeed, we demonstrate that HIT is an efficient sorting mechanism and that,
unsurprisingly, the magnitude of HIT governs the timescale for segregation in
3D.

Unfortunately, it is rather difficult to test this prediction in experiments,
as it is quite difficult to measure the magnitude of HIT. In laser ablation
experiments on monolayers, cuts are made to ablate a cell-cell junction with
the help of a pulsed laser. The resultant relaxation dynamics can help
determine line tensions Sugimura et al. (2016). However, ablating interfaces
and recording retractions along arbitrary interfaces is very difficult in 3D,
and has not been yet used to estimate HIT in 3D. Indirect measures using few-
cell assays, such as double pipette aspiration experiments, suggest that HIT
can create robust cell sorting in 3D tissues Rübsam et al. (2017); Tsai et al.
(2019); Maître et al. (2012).

However, using isolated cells for measuring effective tension may not provide
a complete picture as a confluent neighborhood can significantly change a
cell’s mechanics Canty et al. (2017); Sussman et al. (2018). Recent work has
shown that geometrical properties of interfacial cells in a confluent
neighborhood can be directly affected by increased interfacial contractility
and tension Okuda et al. (2015); Sussman et al. (2018); Yang and Brodland
(2009). Therefore, a second goal is to explore the idea that cellular geometry
can perhaps be used as a simpler and more direct readout of HIT in 3D
mixtures, as geometric features have recently become more accessible in
experiments due to advances in tissue segmentation techniques.

In addition to developing tools for measuring HIT, a third goal of this work
is to identify the mechanisms driving cell sorting in 3D. For a fluid-like
particulate mixture, the mechanism is simple, as the minimum energy state is a
configuration that minimizes the shared surface area. In the presence of large
enough fluctuations the system can find a simple geometry that achieves this
goal via complete spatial segregation with a strong gradient at the interface
Canongia Lopes (2002); Majumder and Das (2011); Chremos et al. (2014). In
confluent tissues, however, there are cell-scale geometric constraints that
may compete with macroscopic interfacial dynamics. Specifically, in vertex
models for isotropic confluent tissues Nagai and Honda (2001); Farhadifar et
al. (2007); Teleman et al. (2007); Staple et al. (2010); Manning et al.
(2010b); Hilgenfeldt et al. (2008); Chiou et al. (2012); Bi et al. (2015);
Fletcher et al. (2014); Merkel and Manning (2018), cells attempt to attain a
preferred cell shape index, which is dimensionless ratio of perimeter P and
area A i.e. $s_{0}^{2D}=P/\sqrt{A}$ in 2D Bi et al. (2015, 2016), and of
surface area S and volume V i.e. $s_{0}^{3D}=S/V^{2/3}$ in 3D Merkel and
Manning (2018). If the cells are able to attain their preferred cell shape,
which generally happens in isotropic tissues for elongated shapes with
$s_{0}^{2D}>~{}3.81$ and $s_{0}^{3D}>~{}5.41$, then the tissue is fluid-like,
while if they are unable to attain that shape the tissue is solid-like.

Moreover, as discussed above, 2D results suggest specific features of the 2D
cell-scale geometry help pin cells at a heterotypic interface, although it is
not trivial to generalize the arguments to 3D. Therefore, we study whether
these cell-scale constraints, which vary with tissue rheology, impact cell
sorting in 3D. We develop simple toy models to demonstrate that a competition
between the bulk cell shape preference and geometric pinning of cell shapes at
the heterotypic interface drive both cell sorting and the formation of
specialized geometric features at the HIT interface. Finally, we confirm the
predictions made by the toy models with full numerical simulations.

## Model

To understand the interfacial mechanics between two cell types, we use the
recently developed 3D Voronoi model Merkel and Manning (2018). A system with
periodic boundaries is created using a Voronoi tessellation of the cell
centers of $N$ cells. Individual cells have preferred volumes $V_{0}$ and
surface area $S_{0}$. The combination of volume incompressibility and surface
area regulation due to adhesion and contractility generate a preferred cell
shape index $s_{0}=S_{0}/V_{0}^{2/3}$. For example, a regular BCC unit cell
(truncated octahedron), has a dimensionless shape index of $s_{0}\sim 5.31$
Lucarini (2009). Here we set $V_{0}$ to 1. Half of the cells are tagged
differently, creating a mixture of two cell types- $\beta=1$ or $\beta=2$,
which are otherwise identical except there is heterotypic interface between
cells of different type. Of course, in experimental systems there are likely
mechanical differences between different cell types in addition to HIT. Here
we ignore those differences to study effects driven by HIT, as our previous
work in 2D suggests that sorting driven by innate mechanical differences is
significantly weaker than HITSahu et al. (2020b). In addition to the original
monodisperse energy functional, we impose an additional surface tension along
the heterotypic interface. Therefore cells minimize their mechanical energy
using the following dimensionless energy functional:

$e=\sum_{i}{\bigg{[}k_{v}(v_{i}-1)^{2}+(s_{i}-s_{0})^{2}\bigg{]}}+\sum_{\langle
i,j\rangle}(1-\delta_{\alpha\beta})\sigma a_{ij}\text{,}$ (1)

where $v_{i}$ denotes the $i$th cell volume and $s_{i}$ denotes its surface
area, non-dimensionalized by $V_{0}$. The unit of length is defined such that
the average cell volume $\langle V_{i}\rangle$ is 1. Additionally,
$k_{s}=K_{V}/K_{S}$ sets the ratio between volume and area stiffnesses, and is
also set to unity. The second summation imposes an additional surface tension
between heterotypic cells, where the sum is over all facets with area $a_{ij}$
shared between cells $i$ and $j$ of types $\alpha$ and $\beta$ respectively.
The surface tension $\sigma$, for simplicity, is assumed to be the same for
all facets. It is non-dimensionalized by $K_{S}V_{0}^{2}$ which is unity for
our system. Biological cells can establish heterotypic tension by co-
regulating the acto-myosin network and adhesion molecules Amack and Manning
(2012). A biologically relevant estimate of the heterotypic to homotypic
tension ratio, based on examination of the contact angles at cell vertices in
ectoderm-mesoderm co-cultures in Xenopus Canty et al. (2017), indicates
$\sigma\sim 2$ in natural units of the system. For systems with fluctuations,
we analyse the dynamics of over-damped self-propelled particles with a high
angular noise, which effectively leads to Brownian dynamics at the timescales
relevant to us. The timescales are reported in units of the self-diffusivity
timescale $\tau_{s}^{0}$, details of which are provided in Supplemental
section S1. While there are other possible dynamical rules, recent work on 2D
mixtures has shown that the properties of an interface between two cells types
with HIT between them is rather robust to the specifics of the dynamical rules
Sussman et al. (2018). For analyzing behavior in the absence of fluctuations,
we use a conjugate gradient minimizer. See Supplemental section S1 for more
details.

## Results

Demixing parameter: To test if HIT leads to significant segregation in 3D
tissues, in a manner similar to that of 2D mixtures Sahu et al. (2020b), we
first focus on a fluid-like parameter regime Merkel and Manning (2018) where
cells undergo diffusive motion ($s_{0}>5.41$). For a fixed shape index
$s_{0}=5.5$, we start from an initially mixed configuration Fig 1(a) with a
system size of $N=512$ cells. We let the system evolve long enough so each
cell on average explores a distance equivalent to the simulation box length.

Let us now understand the role of HIT on the bulk demixing. For a fixed
$\sigma=1$, a final configuration for such a mixture is shown in Fig 1(b). It
is clearly segregated as compared to the initial snapshot. Some fraction of
ensembles are able to create a planar interface as well. To quantify this
demixing, we study the demixing parameter DP Sahu et al. (2020b), which
measures the average neighborhood composition of every cell. Defining $N_{s}$
as the number of homotypic (similar cell type) neighbors and $N_{t}$ as the
total number of neighbors,

$DP=\big{\langle}DP_{i}\big{\rangle}=\bigg{\langle}2\bigg{(}\frac{N_{s}}{N_{t}}-\frac{1}{2}\bigg{)}\bigg{\rangle},$
(2)

where the brackets denote averaging over all cells in the tissue. In a
completely mixed state, $DP=0$, whereas in a completely sorted mixture,
$DP=1$, in the limit of infinite system size. However, as large system sizes
can be time-consuming, we compute the maximum attainable value of demixing
($DP_{max}$) for a particular system size by looking at minimal surface
configurations as shown in Fig. S1. Hence we plot the demixing parameter
relative to this maximum value in Fig 1. The value of demixing is zero at the
beginning as both the cell types are seeded at random positions, but it soon
attains a high value, very close to $DP_{max}$. In the presence of HIT, the
value of demixing increases quickly, indicating that it can efficiently create
robust segregation —very similar to a liquid-liquid particulate mixture and 2D
confluent mixtures.

Figure 1:  Sorting in fluid-like binary tissue with HIT: Initial(a) and
final(b) snapshots of a $s_{0}=5.5$ mixture with high tension of $\sigma=1$
and system size of $N=512$. Both cell types denoted by red and blue polyhedra.
(c) Quantification of segregation: demixing parameter DP versus simulation
time in units of self-diffusivity timescale $\tau_{s}^{0}$ for increasing
tension $\sigma=0.001,0.003,0.010,0.032$ (pink to blue).

In the presence of fluctuations, we find that HIT efficiently leads to
significant segregation in mixed 3D tissues. We also observe that with higher
values of tension, the initial phase of the sorting process becomes faster as
shown in Fig. 1. This confirms that heterotypic surface tension is very
effective at compartmentalization in 3D, as expected.

Geometric features: While biological cells are capable of upregulating tension
cables along heterotypic interfaces via biochemical pathways, it is very
difficult to directly measure this tension within a 3D tissue. Can features of
individual cells at the interface help us quantify such tensions? In this
section, we focus again on fluid-like systems with $s_{0}=5.5$. A visual
inspection of the segregated mixture shown in Fig 1(b) indicates that the
interfacial cells may be more elongated and nematically ordered as compared to
the cells in the interior. This observation hints at a direct relationship
between the applied tension and the surrounding cellular geometry. To delve
deeper into the ways in which the surrounding cells deform, we set up a
maximally segregated mixture. Here, both compartments are placed side by side,
similar to previous work by Sussman et al Sussman et al. (2018). We then study
the cellular geometry as a function of the applied interfacial tension. Shape
elongation along with the prominent stacking of cells (alignment of the
polyhedral long axes in Fig. 2(a)) can be observed here as well. We next
quantify such geometric effects.

The first quantity is the steady-state cell shape index $s$. This helps us
quantify whether or not the otherwise homogeneous cells remain homogeneous
after HIT is established. We measure the shape of both interfacial cells-
cells that directly touch the boundary ($s_{boundary}$), and the interior
cells ($s_{bulk}$). This further helps isolate the shape changes in the
immediate neighborhood of the interface. Individual cells have a final volume
$V_{i}$ and surface area $S_{i}$. Hence, $s_{boundary}$ is defined as:

$s_{boundary}=\bigg{\langle}\frac{S}{V^{2/3}}\bigg{\rangle}_{boundary},$ (3)

and $s_{bulk}$ is similarly averaged over the interior cells. In the absence
of HIT ($\sigma=0$), both shapes have the same value of $s_{0}=5.5$.

Both the shape indices are shown as a function of increasing tension in Fig.
2(b). For small values of tension, the shapes indices are very similar at the
beginning, but they gradually saturate at a higher value of disparity. In
other words, with a higher interfacial tension, the neighboring cells become
more elongated, whereas the interior cells become more compact/round.

To study the alignment between cells, we next measure the orientation of
interfacial polyhedra. We define orientation vector of a cell as the major
axis of its moment of inertia tensor. We then plot the angular distribution of
the angle made by each vector with respect to the normal to the interface
($\theta$), which in this bilayer arrangement is simply the z axis. For a
homogeneous system with no HIT, the angles are very close to the random
distribution density function in 3D, which is proportional to $\sin\theta$.
But with a slight increase in tension, the cells polarize and orient
themselves perpendicular to the interface, as shown in Fig. 2(c).

Lastly, we study polygonal faces that make up the heterotypic interface by
plotting the facet area distribution with respect to increasing tension. This
is in analogy to the measurement of edge lengths in 2D work Sussman et al.
(2018). With no HIT, the distribution is roughly uniform up to a
characteristic length scale, whereas, with increasing tension it becomes
bimodal, as shown in Fig. 2(d). This means that the facets are either large or
vanishingly small at high tensions. With the help of smaller facets, the
interfacial vertices come very close to having a vertex coordination higher
than the normal tetrahedral coordination, similar to the 4-fold vertices
observed along 2D tension cables. The average area of a facet, $<a>$, also
increased with increasing tension as shown in Fig. S2. This quantitatively
confirms that HIT indeed affects the geometry of the surrounding cells,
inducing shape changes and nematic-like ordering in an otherwise homogeneous
collection of cells.

Since some similar geometric features have been observed in 2D models for
confluent tissues Sussman et al. (2018), we hypothesize that the origin of
these signatures in 3D might be based on topological interactions between the
cells which have been implicated in 2D.

Figure 2: Cellular geometry changes around the high-tension interface: (a)
Snapshot of the bilayer arrangement of sorted compartments for a high value of
HIT, $\sigma=20$. Only the blue cell type is emphasized here, colored by major
axis length. Greenish-blue is for elongated and purple for rounder cells. The
white rods denote the long axis of the polyhedron. (Caption continued on next
page.) Figure 2: (Continued from previous page.) (b) Acquired cell shape index
(s), plotted with respect to tension ($\sigma$), is higher for interfacial
cells (solid green curve) as compared to the cells in bulk (dashed blue
curve). (c) Rose plot for the distribution of orientation angle of interfacial
cells is shown for increasing tension ($\sigma$). The control distribution
($\sigma$=0) is superimposed on a faint black curve that represents the
density for uniform distribution i.e $\sin\theta$. (d) The probability
distribution $P(a)$ of heterotypic facet area $a$ is shown for increasing
tension $\sigma=0.32,1.00,2.00$ (magenta to blue). For these data, N=512 and
$s_{0}=5.5$.

In order to determine the mechanism driving these geometric changes, we first
focus on the specific geometry of an interfacial cell and ask: what does it
take for the system to attain this precise geometry?

A typical interfacial neighborhood is shown in Fig. 3(a). A right prism can be
defined as a polyhedron with flat top and bottom facets, and perpendicularly
aligned lateral facets. The cells here seem to closely resemble the geometry
of a one-sided right prism, with the flat side at the interface. The unique
shape can be attained in a Voronoi tessellation only by fulfilling two
conditions : (a) cell heights are arranged in a plane, and (b) interfacial
pairs align their centers so the distance between centers in the XY plane is
minimized. One way to quantify this alignment is to measure it as the
registration- $R$ between the cell centers. We define it as-

$R=1-\frac{d}{l_{0}},$ (4)

where $l_{0}$ is the typical lengthscale and $d$ is the distance between
centers along the interfacial plane as depicted in Fig. 3(a). The value of
$l_{0}$ is set to unity as $V_{0}^{1/3}=1$ in our systems. When cell centers
are completely registered its value attains a maximum value of unity. We
quantify both of the aforementioned criteria and find them both fulfilled for
higher values of tension, as shown in Fig. 3(b-c). Moreover, the insets to
each figure show that the cell geometries approach the right-prism condition
monotonically as a function of increasing HIT.

Figure 3: Interfacial cells are prism-like: (a) A simulation snapshot
highlighting the cells at both sides of the interface. Cell types are tagged
in blue and red and made translucent to make their centers (white solid
spheres) visible. To characterize the high alignment between the heterotypic
cell centers we use cell-cell registration $R$ defined in the main text in
terms of the distance $d$ between the centers along the interface (highlighted
in red) and the characteristic cell length $l_{0}$ (denoted by the black scale
bar). (b) Distribution $pdf(z)$ of the heights of cell centers $z$ for blue
interface cells is plotted for increasing tension $\sigma=0.03,0.10,0.31,1.00$
(magenta to blue). The dashed black line depicts a system with no HIT. The
inset shows the standard deviation $\Delta z$ with as a function of the HIT
$\sigma$. (c) Distribution of registries pdf($R$) between heterotypic cells is
plotted for increasing tension. The inset shows the averaged registry with
respect to tension. The lowest value of tension has the same average distance
as no tension at all. The circle radii corresponds to the deviation from the
mean values.

From work in 2D we understand that perturbations perpendicular to the surface
of the interface are very costly Sussman et al. (2018), but cell-cell
registration requires perturbations to the cell center that are parallel to
the interface. In the SM we extend the 2D calculation to study the effects of
parallel perturbations by changing the cell-cell registration. In the 2D toy
model, we find that parallel perturbation incurs exactly the same energy
penalty as a perpendicular perturbation up to linear order. Similar to the
perpendicular perturbation, the parallel perturbation along the heterotypic
interface requires the system to form new, high tension edges which are
energetically very costly. Therefore, the cell centers are laterally pinned.
In S5 we explore if de-registration in 3D systems creates a similar effect. We
find that in an idealized hexagonal tissue, such a perturbation creates
several new HIT facets, suggesting that a similar mechanism is operating in
3D.

Mechanisms: To better understand the pinning mechanism in 3D, we next develop
two ordered toy models in 3D and study the response in the limit of zero
fluctuations. Although the data shown so far, and work by some of us in 2D
Sussman et al. (2018), focused exclusively on the fluid-like regime with
$s_{0}^{2D}>3.81$ and $s_{0}^{3D}>5.41$, we are also interested in how tissues
close to the fluid-solid transition balance interfacial and bulk effects.
Therefore, our first toy model is initialized in a BCC lattice, as a ground
state for the shape $s_{0}\sim 5.31$. We fix $s_{0}=5.31$ to this ground state
value, which we expect to be solid-like. This is similar to our recent work on
ordered hexagonal monolayers Sahu et al. (2020a), but here we are
investigating a different form of local perturbation in 3D.

In this work, we study the response of the system to a change in the registry
$R$. In these ordered systems, we now use the lattice spacing as the
lengthscale $l_{0}$ in the definition for registry in Eq. 4. We enable a
string of polyhedra to slide past the string below (shown in Fig. 4(a)
insets), with an extra surface tension along the shared interfacial strip
between both the sets of polyhedra. We compute the shared surface area and
total energy of the system. We also find the global minima for different
values of tension.

We find that the shared surface area is minimized for increased registration
values (shown in Fig. 4(a)). While this suggests that perhaps complete
registration is an energetically preferred state in some regimes,
surprisingly, that is true only after a threshold value of tension. This can
be seen by plotting the change in total energy with respect to increasing
registration. One can observe two kinds of minima in the system —parabolic and
‘cuspy’ Sussman et al. (2018). While the former is common in particle-based
models, with a spring-like potential locally around the minima, the later has
a discontinuity in its slope due to topological pinning, which is expected
from a 2D calculation as discussed in S4. Physically, this means that the
system experiences a steep linear rise in energy due to the formation of a new
interfacial edge along the high tension cable, resulting in discontinuous
pinning forces that are independent of the magnitude of perturbation, as shown
in Fig. S6(b). The parabolic minima, on the other hand, have a continuous an
linear restoring force, the slope proportional to the stiffness of the
parabola. We plot the registration of energy minima as a function of
interfacial tension in Fig. 4(c).

Figure 4: Global minimum becomes registered for higher tension:(a) Shared
surface area (shaded in dark for all snapshots) is computed as a function of
the registration between the different cell types. The string of blue cells is
allowed to move past the string of stationary pink cells. Snapshots for no,
half and complete registration is shown for encircled points. (b) The change
in the total energy of this system is plotted with respect to registry, for
different values of tension ranging from 1 (pink) to 4 (blue) in increments of
0.25. (c) The ground-state registration is plotted for increasing tension. The
solid curve represents the global minima. The dashed curves represent the
local minima, parabolic in blue and cuspy in black. The inset zooms about the
critical tension at which the transitions occurs between both types of minima.
For all panels $s_{0}$ is set to $5.31$, the ground state for the BCC lattice.

We observe two different branches corresponding to the two types of minima
discussed before. For very low tension, shape preference dominates (as shown
in Fig. S5(a)) and the system stays in BCC configuration. For moderate values
of tension the system stays in the parabolic branch, but continuously
transitions to non-zero registry. Just below the critical value of
$\sigma_{c}\sim 3$, both the shape and interfacial energies become comparable
(Fig. S5(b)) and at $\sigma_{c}$ the system discontinuously transitions to the
tension-dominated- cuspy branch. The registration value also jumps to $R=1$.
This branch originates at $\sigma\sim 2$.

Both kinds of minima become more stable as the tension increases. Stability is
parameterized by the curvature or stiffness for parabolic minima in Fig. S6(a)
and by the linear slope or restoring force for the cuspy minima in Fig. S6(c).

In summary, for the solid-like toy model, we find that the physical mechanism
that drives registration at high tension values is very similar to that of 2D
systems. However, the story changes at lower values of tension where shape
frustration begins to play a dominant role. This leads to minima that are
partially-registered and exhibit a spring-like response to small
perturbations, i.e. states are no longer topologically pinned.

In 2D studies focused exclusively on fluid-like tissues Sussman et al. (2018),
only the cusp-like states were observed. This leads us to hypothesize that
perhaps we observe a transition between normal, parabolic restoring forces to
non-analytic cusp-y restoring forces in tissues that are more solid-like.
Perhaps in more fluid-like systems the interfacial costs always dominate over
the cost of cell shapes in the bulk, whereas in more solid-like systems the
bulk effects dominate when the interfacial tension is low.

To test this hypothesis, we develop a second toy model that is also ordered
but not constrained to a string. Instead it is free to move along the 2D
interface to change its registration. With this flexibility, we can explore
the energetics of more elongated cell shapes like that of a uniform hexagonal
prism ($s_{0}\sim 5.72$). The interfacial cells can therefore be much more
elonagated and fluid-like as compared to the minimal-perimeter-cells in a BCC
lattice that have shape values as small as $s_{0}\sim 5.31$. The interfacial
layers are placed in HCP format as shown in Fig. 5(b). There are buffer cells
placed above and below the interface in a disordered fashion and are allowed
to relax during the course of the perturbation.

Analogous to the the previous analysis, we compute the shared surface area and
change in the energy profile with respect to registration, but this time for a
wide range of cell shapes across the fluid-solid transition. We find that
similar to the BCC toy model, the shared surface area decreases with
registration as shown in Fig. S7(a). For cell shapes near the rigidity
transition, shape frustration plays a dominant role in the system’s energy
(Fig. S7(b)). However, it becomes negligible for more fluid-like cell shapes
as shown in Fig. 5. This suggests that in fully disordered fluid-like systems,
interfacial tension dominates shape preferences in determining interface
geometry and response.

Figure 5: Shape frustration is less dominant for more fluid-like cell shapes:
(a,b) Snapshots for both the models are shown where both cell types are
depicted in blue and red. Cell centers are depicted by solid spheres of
respective cell color. (c) The ground-state registration for HCP toy model
(filled hexagrams) is plotted with respect to tension, for increasing values
of cell shape $s_{0}=5.4,5.5,5.6,5.7$ from orange to green. For comparison,
the registration in BCC model (solid circles) is also shown here in red.

Finally, we verify the toy model predictions by first studying the geometry of
disordered mechanically stable states in planar segregated HIT simulations,
for different values of HIT and preferred cell shape. We let the system come
to a steady state in the presence of fluctuations that can help the system
find lower and lower energy metastable minima in the complex potential energy
landscape.

Fig 6 shows the registration as a function of interfacial tension and target
cell shape index $s_{0}$. The data demonstrate that just as in the toy models,
the average registration rises rapidly to unity – its maximum value – when the
interfacial tension increases above a threshold value. Moreover, for solid-
like shapes the onset occurs at a higher threshold tension of order unity,
while for fluid-like tissues the onset occurs when HIT is more than an order
of magnitude lower than the typical tensions between homotypic cells. In the
supplement, we also confirm that the highly registered states are associated
with cusp-like restoring forces, highlighted in Fig. S8.

Figure 6: Transition to complete registration shifts for fluid-like cell shape
in a disordered 3D simulation: A heat-map for the average steady-state
registration is shown as function of cell shape $s_{0}$ and interfacial
tension $\sigma$. Yellow denotes complete registration and blue denotes
partial registry. The registry between heterotypic cell-pairs is averaged over
200 different initializations.

## Discussion and Conclusions

By studying a computational Voronoi model for confluent tissues in the
presence of fluctuations, we show that three-dimensional binary mixtures, with
heterotypic interfacial tension (HIT), sort robustly. This supports the claim
that HIT is an important mechanism driving sharp compartmentalization during
early embryonic development Kesavan et al. (2020); Canty et al. (2017), and
that it may be important for tissue segregation in other situations. In
addition to these collective dynamics, HIT also drives individual cells
towards a prism-like geometry at the interface. We find that the onset of
these geometric signatures depends on a balance between the the magnitude of
interfacial tension and constraints introduced by the the preferred cell shape
that also govern the bulk tissue rheology.

To understand the onset of these geometric signatures, we use cell-cell
registration at the interface as a probe of the stability of these prism-like
structures. We construct two simple toy models and study their energetics with
respect to registry. In solid-like tissues, we find that for an interfacial
tension $\sigma>\sigma_{c}$, the ground state is completely registered, which
gives rise to a prismatic geometry. These states are topologically pinned, due
to cusp-like pinning forces, previously observed in 2D mixtures with HIT. But
for tensions $\sigma<\sigma_{c}$, interfacial energy is dominated by shape
frustration and hence the ground state is less well-registered. The linear
response of these minima is spring-like and not cuspy. Our data suggests that
$\sigma_{c}$ decreases significantly as the tissue becomes more fluid-like.
This has important implications for development and tissue segregation, as it
suggests that as tissues are tuned to be more solid-like, topological pinning
at heterotypic interfaces is greatly reduced, thereby reducing the sharpness
possible at compartment boundaries. In other words, it suggests the somewhat
counterintuitive design principle for confluent systems that fluid-like
rheologies lead to sharper interfaces.

We have also shown that a change in the magnitude of the interfacial tension
can have a pronounced effect on the neighboring cellular geometry, by
elongating interfacial cells into prism-like polyhedra oriented perpendicular
to the interface. The observable facet areas also become larger. With the
advancement of 3D segmentation techniques Khan et al. (2014); Stegmaier et al.
(2016); Browet et al. (2016); Fernández-de Manúel et al. (2017); Morales-
Navarrete et al. (2019), one can use these signatures as a toolkit to probe
interfacial tensions in a 3D tissue, so that cells themselves can tell us
about relative magnitudes of tissue surface tension nearby.

One example case where this might be useful is in detecting the invasiveness
of a carcinoma tumour. Our simple model would predict that if interfacial cell
centers in the tumor are registered to those of the surrounding healthy
tissue, then the interface has a higher surface tension and therefore it may
be more unlikely for cells to exit the tumor and invade their surroundings. It
would be interesting to see if there are any vestiges of this prediction that
occur in pre-cancerous situations in real-world systems, such as Ductal
carcinoma in situ (DCIS).

Another example is that of a stratified epithelium, where one can also study
the interaction between two nearby tissue types, such as the basal and
suprabasal layers, and look for geometric signatures across the interface. A
prism-like geometry would strongly suggest the presence of an interfacial
tension between these two tissue types. The prism-like geometry of cells can
be visually detected using the En Face imaging techniqueYokouchi et al.
(2016); Rübsam et al. (2017). As some current segmentation algorithms can also
make predictions about 3D shape from random 2D cross-sections, these geometric
signatures could potentially be characterized along the cross-sections Sharp
et al. (2019). In general 3D tissue have complex interfacial geometries.
Hence, one of the future avenues of this work would be to study the dependence
on the curvature of the interface.

Additional work should also focus on teasing apart how topological pinning
affects dynamics in more general scenarios. For example, a natural extension
of our general framework is to study two different cell shapes mixed together.
After all, in realistic situations cells of two different tissue types likely
also differ in preferred cell shape. In 2D, some of us have determined that
unique extrusion behaviour can emerge due to differential pinning of cells
Sahu et al. (2020b), and something similar could occur in 3D. Topological
pinning might be affecting the cell sorting dynamics as well. Presumably, less
pinning can lead to seamless coarsening of nearby droplets, while more pinning
can hinder the coalescence, and alter the sorting process. This would also be
an interesting avenue for future study.

While our work demonstrates that changes to individual cellular geometries are
a necessary consequence of HIT and tissue-scale segregation, these changes to
cellular geometry could also be used as a signal to facilitate downstream
patterning near interfaces. For example, the elongation of cells near a high-
tension interface might trigger oriented cell divisions along the long-axis of
these cells. We speculate that perhaps biology can make use of this subtle
feedback for processes like targeted cellular proliferation during early
phases of embryonic development.

#### Acknowledgement

We thank Paula Sanematsu, Matthias Merkel, Daniel Sussman and Cristina
Marchetti for helpful discussions, and M. Merkel for developing and sharing
the original version of the 3D Voronoi code. This work was primarily funded by
NSF-PoLS-1607416 and NSF-PoLS-2014192. PS and MLM acknowledge additional
support from Simons Grant No. 454947.

## References

  * Trinkaus and Groves [1955] J. P. Trinkaus and P. W. Groves.  Differentiation in culture of mixed aggregates of dissociated tissue cells.  _Proc. Natl. Acad. Sci._ , 1955.  ISSN 0027-8424.
  * Waites et al. [2017] William Waites, Matteo Cavaliere, Elise Cachat, Vincent Danos, and Jamie A. Davies.  Organoid And Tissue Patterning Through Phase Separation: Use Of A Vertex Model To Relate Dynamics Of Patterning To Underlying Biophysical Parameters.  _bioRxiv_ , 2017.
  * Unbekandt and Davies [2010] Mathieu Unbekandt and Jamie A. Davies.  Dissociation of embryonic kidneys followed by reaggregation allows the formation of renal tissues.  _Kidney Int._ , 2010.  ISSN 00852538.
  * Montero et al. [2005] Juan Antonio Montero, Lara Carvalho, Michaela Wilsch-Bräuninger, Beate Kilian, Chingdem Mustafa, and Carl Philipp Heisenberg.  Shield formation at the onset of zebrafish gastrulation.  _Development_ , 2005.  ISSN 09501991.
  * Klopper et al. [2010] A. V. Klopper, G. Krens, S. W. Grill, and C. P. Heisenberg.  Finite-size corrections to scaling behavior in sorted cell aggregates.  _Eur. Phys. J. E_ , 2010.  ISSN 12928941.
  * Rübsam et al. [2017] Matthias Rübsam, Aaron F. Mertz, Akiharu Kubo, Susanna Marg, Christian Jüngst, Gladiola Goranci-Buzhala, Astrid C. Schauss, Valerie Horsley, Eric R. Dufresne, Markus Moser, Wolfgang Ziegler, Masayuki Amagai, Sara A. Wickström, and Carien M. Niessen.  E-cadherin integrates mechanotransduction and EGFR signaling to control junctional tissue polarization and tight junction positioning.  _Nat. Commun._ , 2017.  ISSN 20411723.
  * Cochet-Escartin et al. [2017] Olivier Cochet-Escartin, Tiffany T. Locke, Winnie H. Shi, Robert E. Steele, and Eva Maria S. Collins.  Physical Mechanisms Driving Cell Sorting in Hydra.  _Biophys. J._ , 2017.  ISSN 15420086.
  * Friedl et al. [2004] Peter Friedl, Yael Hegerfeldt, and Miriam Tusch.  Collective cell migration in morphogenesis and cancer, 2004.  ISSN 02146282.
  * Foty and Steinberg [2005] Ramsey A. Foty and Malcolm S. Steinberg.  The differential adhesion hypothesis: A direct evaluation.  _Dev. Biol._ , 2005.  ISSN 00121606.
  * Pawlizak et al. [2015] Steve Pawlizak, Anatol W. Fritsch, Steffen Grosser, Dave Ahrens, Tobias Thalheim, Stefanie Riedel, Tobias R. Kießling, Linda Oswald, Mareike Zink, M. Lisa Manning, and Josef A. Käs.  Testing the differential adhesion hypothesis across the epithelial-mesenchymal transition.  _New J. Phys._ , 17:091001, 2015.  ISSN 13672630.
  * Song et al. [2016] Wei Song, Chih Kuan Tung, Yen Chun Lu, Yehudah Pardo, Mingming Wu, Moumita Das, Der I. Kao, Shuibing Chen, and Minglin Ma.  Dynamic self-organization of microwell-aggregated cellular mixtures.  _Soft Matt._ , 12:5739, 2016.  ISSN 17446848.
  * Turing [1990] A. M. Turing.  The chemical basis of morphogenesis.  _Bull. Math. Biol._ , 1990.  ISSN 00928240.
  * Streichan et al. [2018] Sebastian J Streichan, Matthew F Lefebvre, Nicholas Noll, Eric F Wieschaus, and Boris I Shraiman.  Global morphogenetic flow is accurately predicted by the spatial distribution of myosin motors.  _Elife_ , 2018.
  * Steinberg [1963] Malcolm S. Steinberg.  Reconstruction of tissues by dissociated cells.  _Sci. Sci._ , 1963.  ISSN 00368075.
  * Steinberg [2007] Malcolm S Steinberg.  Differential adhesion in morphogenesis: a modern view.  _Current Opinion in Genetics and Development_ , 17(4):281 – 286, 2007.  ISSN 0959-437X.  URL http://www.sciencedirect.com/science/article/pii/S0959437X07001062.  Pattern formation and developmental mechanisms.
  * Harris [1976] Albert K. Harris.  Is cell sorting caused by differences in the work of intercellular adhesion? A critique of the steinberg hypothesis.  _J. Theor. Biol._ , 1976.  ISSN 10958541.
  * Brodland [2002] G. Wayne Brodland.  The Differential Interfacial Tension Hypothesis (DITH): A comprehensive theory for the self-rearrangement of embryonic cells and tissues.  _J. Biomech. Eng._ , 2002.  ISSN 01480731.
  * Mertz et al. [2012] Aaron F. Mertz, Shiladitya Banerjee, Yonglu Che, Guy K. German, Ye Xu, Callen Hyland, M. Cristina Marchetti, Valerie Horsley, and Eric R. Dufresne.  Scaling of traction forces with the size of cohesive cell colonies.  _Phys. Rev. Lett._ , 2012.  ISSN 00319007.
  * Maître et al. [2012] Jean Léon Maître, Hélène Berthoumieux, Simon Frederik Gabriel Krens, Guillaume Salbreux, Frank Jülicher, Ewa Paluch, and Carl Philipp Heisenberg.  Adhesion functions in cell sorting by mechanically coupling the cortices of adhering cells.  _Science (80-. )._ , 2012.  ISSN 10959203.
  * Amack and Manning [2012] Jeffrey D. Amack and M. Lisa Manning.  Knowing the boundaries: Extending the differential adhesion hypothesis in embryonic cell sorting, 2012.  ISSN 10959203.
  * Manning et al. [2010a] M. Lisa Manning, Ramsey A. Foty, Malcolm S. Steinberg, and Eva Maria Schoetz.  Coaction of intercellular adhesion and cortical tension specifies tissue surface tension.  _Proc. Natl. Acad. Sci. U. S. A._ , 2010a.  ISSN 00278424.
  * Engl et al. [2014] W. Engl, B. Arasi, L. L. Yap, J. P. Thiery, and V. Viasnoff.  Actin dynamics modulate mechanosensitive immobilization of E-cadherin at adherens junctions.  _Nat. Cell Biol._ , 2014.  ISSN 14764679.
  * Graner and Glazier [1992] François Graner and James A. Glazier.  Simulation of biological cell sorting using a two-dimensional extended Potts model.  _Phys. Rev. Lett._ , 69:213–216, 1992.  ISSN 00319007.
  * Barton et al. [2017] D. L. Barton, S. Henkes, C. J. Weijer, and Sknepnek R.  Active vertex model for cell-resolution description of epithelial tissue mechanics.  _PLoS Comput. Biol._ , 13:e1005569, 2017.
  * Canty et al. [2017] Laura Canty, Eleyine Zarour, Leily Kashkooli, Paul François, and François Fagotto.  Sorting at embryonic boundaries requires high heterotypic interfacial tension.  _Nat. Comm._ , 8:157, 2017.
  * Sussman et al. [2018] Daniel M. Sussman, J. M. Schwarz, M. Cristina Marchetti, and M. Lisa Manning.  Soft yet Sharp Interfaces in a Vertex Model of Confluent Tissue.  _Phys. Rev. Lett._ , 120:058001, 2018.
  * Flenner et al. [2012] Elijah Flenner, Lorant Janosi, Bogdan Barz, Adrian Neagu, Gabor Forgacs, and Ioan Kosztin.  Kinetic monte carlo and cellular particle dynamics simulations of multicellular systems.  _Phys. Rev. E_ , 85:031907, Mar 2012.
  * Merks and Glazier [2005] Roeland M.H. Merks and James A. Glazier.  A cell-centered approach to developmental biology.  _Phys. A Stat. Mech. its Appl._ , 2005.  ISSN 03784371.
  * Sun and Wang [2013] Yi Sun and Qi Wang.  Modeling and simulations of multicellular aggregate self-assembly in biofabrication using kinetic Monte Carlo methods.  _Soft Matter_ , 2013.  ISSN 17446848.
  * Jiang et al. [1998] Yi Jiang, Herbert Levine, and James Glazier.  Possible cooperation of differential adhesion and chemotaxis in mound formation of Dictyostelium.  _Biophys. J._ , 1998.  ISSN 00063495.
  * Palsson [2008] Eirikur Palsson.  A 3-D model used to explore how cell adhesion and stiffness affect cell sorting and movement in multicellular systems.  _J. Theor. Biol._ , 2008.  ISSN 00225193.
  * Landsberg et al. [2009] Katharina P. Landsberg, Reza Farhadifar, Jonas Ranft, Daiki Umetsu, Thomas J. Widmann, Thomas Bittig, Amani Said, Frank Jülicher, and Christian Dahmann.  Increased Cell Bond Tension Governs Cell Sorting at the Drosophila Anteroposterior Compartment Boundary.  _Curr. Biol._ , 2009.  ISSN 09609822.
  * Manning et al. [2010b] M. Lisa Manning, Ramsey A. Foty, Malcolm S. Steinberg, and Eva-Maria Schoetz.  Coaction of intercellular adhesion and cortical tension specifies tissue surface tension.  _Proc. Natl. Acad. Sci._ , 2010b.  ISSN 0027-8424.
  * Dahmann et al. [2011] Christian Dahmann, Andrew C. Oates, and Michael Brand.  Boundary formation and maintenance in tissue development, 2011.  ISSN 14710056.
  * Nnetu et al. [2012] Kenechukwu David Nnetu, Melanie Knorr, Josef Käs, and Mareike Zink.  The impact of jamming on boundaries of collectively moving weak-interacting cells.  _New J. Phys._ , 2012.  ISSN 13672630.
  * Monier et al. [2011] Bruno Monier, Anne Pélissier-Monier, and Bénédicte Sanson.  Establishment and maintenance of compartmental boundaries: Role of contractile actomyosin barriers, 2011.  ISSN 1420682X.
  * Calzolari et al. [2014] Simone Calzolari, Javier Terriente, and Cristina Pujades.  Cell segregation in the vertebrate hindbrain relies on actomyosin cables located at the interhombomeric boundaries.  _EMBO J._ , 2014.  ISSN 14602075.
  * Kesavan et al. [2020] Gokul Kesavan, Anja Machate, Stefan Hans, and Michael Brand.  Cell-fate plasticity, adhesion and cell sorting complementarily establish a sharp midbrain-hindbrain boundary.  _Development_ , 2020.  ISSN 14779129.
  * Sahu et al. [2020a] Preeti Sahu, Janice Kang, Gonca Erdemci-Tandogan, and M. Lisa Manning.  Linear and nonlinear mechanical responses can be quite different in models for biological tissues.  _Soft Matter_ , 2020a.  ISSN 17446848.
  * Hutson et al. [2008] M. Shane Hutson, G. Wayne Brodland, Justina Yang, and Denis Viens.  Cell sorting in three dimensions: Topology, fluctuations, and fluidlike instabilities.  _Phys. Rev. Lett._ , 101:148105, Oct 2008.
  * Sugimura et al. [2016] Kaoru Sugimura, Pierre François Lenne, and François Graner.  Measuring forces and stresses in situ in living tissues.  _Dev._ , 2016.  ISSN 14779129.
  * Tsai et al. [2019] Tony Y.-C. Tsai, Mateusz Sikora, Peng Xia, Tugba Colak-Champollion, Holger Knaut, Carl-Philipp Heisenberg, and Sean G. Megason.  An adhesion code ensures robust pattern formation during tissue morphogenesis.  _bioRxiv_ , 2019.
  * Okuda et al. [2015] Satoru Okuda, Yasuhiro Inoue, and Taiji Adachi.  Three-dimensional vertex model for simulating multicellular morphogenesis.  _Biophys. Physicobiology_ , 2015.  ISSN 2189-4779.
  * Yang and Brodland [2009] Justina Yang and G. Wayne Brodland.  Estimating interfacial tension from the shape histories of cells in compressed aggregates: A computational study.  _Ann. Biomed. Eng._ , 2009.  ISSN 00906964.
  * Canongia Lopes [2002] José N. Canongia Lopes.  Microphase separation in mixtures of Lennard-Jones particles.  In _Phys. Chem. Chem. Phys._ , 2002.
  * Majumder and Das [2011] Suman Majumder and Subir K. Das.  Universality in fluid domain coarsening: The case of vapor-liquid transition.  _EPL_ , 2011.  ISSN 02955075.
  * Chremos et al. [2014] Alexandros Chremos, Arash Nikoubashman, and Athanassios Z. Panagiotopoulos.  Flory-Huggins parameter $\chi$, from binary mixtures of Lennard-Jones particles to block copolymer melts.  _J. Chem. Phys._ , 2014.  ISSN 00219606.
  * Nagai and Honda [2001] Tatsuzo Nagai and Hisao Honda.  A dynamic cell model for the formation of epithelial tissues.  _Phil. Mag. B: Phys. Cond. Matt., Stat. Mech., Elec. Opt. Mag. Prop._ , 81:699–719, 2001.
  * Farhadifar et al. [2007] 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.  _Curr. Biol._ , 17(24):2095–2104, 2007.  ISSN 09609822.
  * Teleman et al. [2007] A. A. Teleman, L. Hufnagel, H. Rouault, B. I. Shraiman, and S. M. Cohen.  On the mechanism of wing size determination in fly development.  _Proc. Natl. Acad. Sci._ , 2007.  ISSN 0027-8424.
  * Staple et al. [2010] D. B. Staple, R. Farhadifar, J. C. Röper, B. Aigouy, S. Eaton, and F. Jülicher.  Mechanics and remodelling of cell packings in epithelia.  _Eur. Phys. J. E_ , 33:117–127, 2010.
  * Hilgenfeldt et al. [2008] S. Hilgenfeldt, S. Erisken, and R. W. Carthew.  Physical modeling of cell geometric order in an epithelial tissue.  _Proc. Natl. Acad. Sci._ , 2008.  ISSN 0027-8424.
  * Chiou et al. [2012] Kevin K. Chiou, Lars Hufnagel, and Boris I. Shraiman.  Mechanical stress inference for two dimensional cell arrays.  _PLoS Comput. Biol._ , 2012.  ISSN 1553734X.
  * Bi et al. [2015] Dapeng Bi, J. H. Lopez, J. M. Schwarz, and M. Lisa Manning.  A density-independent rigidity transition in biological tissues.  _Nat. Phys._ , 11(12):1074–1079, 2015.  ISSN 17452481.
  * Fletcher et al. [2014] Alexander G. Fletcher, Miriam Osterfield, Ruth E. Baker, and Stanislav Y. Shvartsman.  Vertex models of epithelial morphogenesis, 2014.  ISSN 15420086.
  * Merkel and Manning [2018] Matthias Merkel and M. Lisa Manning.  A geometrically controlled rigidity transition in a model for confluent 3D tissues.  _New J. Phys._ , 2018.  ISSN 13672630.
  * Bi et al. [2016] Dapeng Bi, Xingbo Yang, M. Cristina Marchetti, and M. Lisa Manning.  Motility-driven glass and jamming transitions in biological tissues.  _Phys. Rev. X_ , 2016.  ISSN 21603308.
  * Lucarini [2009] Valerio Lucarini.  Three-dimensional random voronoi tessellations: From cubic crystal lattices to poisson point processes.  _J. Stat. Phys._ , 2009.  ISSN 00224715.
  * Sahu et al. [2020b] Preeti Sahu, Daniel M. Sussman, Matthias Rübsam, Aaron F. Mertz, Valerie Horsley, Eric R. Dufresne, Carien M. Niessen, M. Cristina Marchetti, M. Lisa Manning, and J. M. Schwarz.  Small-scale demixing in confluent biological tissues.  _Soft Matter_ , 2020b.  ISSN 17446848.
  * Khan et al. [2014] Zia Khan, Yu Chiun Wang, Eric F. Wieschaus, and Matthias Kaschube.  Quantitative 4D analyses of epithelial folding during Drosophila gastrulation.  _Dev._ , 2014.  ISSN 14779129.
  * Stegmaier et al. [2016] Johannes Stegmaier, Fernando Amat, William C. Lemon, Katie McDole, Yinan Wan, George Teodoro, Ralf Mikut, and Philipp J. Keller.  Real-Time Three-Dimensional Cell Segmentation in Large-Scale Microscopy Data of Developing Embryos.  _Dev. Cell_ , 2016.  ISSN 18781551.
  * Browet et al. [2016] A. Browet, C. De Vleeschouwer, L. Jacques, N. Mathiah, B. Saykali, and I. Migeotte.  Cell segmentation with random ferns and graph-cuts.  In _Proc. - Int. Conf. Image Process. ICIP_ , 2016.  ISBN 9781467399616.
  * Fernández-de Manúel et al. [2017] Laura Fernández-de Manúel, Covadonga Díaz-Díaz, Daniel Jiménez-Carretero, Miguel Torres, and María C. Montoya.  ESC-track: A computer workflow for 4-D segmentation, tracking, lineage tracing and dynamic context analysis of ESCs.  _Biotechniques_ , 2017.  ISSN 19409818.
  * Morales-Navarrete et al. [2019] Hernán Morales-Navarrete, Hidenori Nonaka, André Scholich, Fabián Segovia-Miranda, Walter de Back, Kirstin Meyer, Roman L. Bogorad, Victor Koteliansky, Lutz Brusch, Yannis Kalaidzidis, Frank Jülicher, Benjamin M. Friedrich, and Marino Zerial.  Liquid-crystal organization of liver tissue.  _Elife_ , 2019.  ISSN 2050084X.
  * Yokouchi et al. [2016] Mariko Yokouchi, Toru Atsugi, Mark Van Logtestijn, Reiko J. Tanaka, Mayumi Kajimura, Makoto Suematsu, Mikio Furuse, Masayuki Amagai, and Akiharu Kubo.  Epidermal cell turnover across tight junctions based on Kelvin’s tetrakaidecahedron cell shape.  _Elife_ , 2016.  ISSN 2050084X.
  * Sharp et al. [2019] Tristan A. Sharp, Matthias Merkel, M. Lisa Manning, and Andrea J. Liu.  Inferring statistical properties of 3D cell geometry from 2D slices.  _PLoS One_ , 2019.  ISSN 19326203.
  * Marchetti et al. [2013] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, Madan Rao, and R. Aditi Simha.  Hydrodynamics of soft active matter.  _Rev. Mod. Phys._ , 85:1143–1189, Jul 2013.  URL https://link.aps.org/doi/10.1103/RevModPhys.85.1143.

## Electronic Supplemental Material

###  S1 Details about self-propulsion dynamics

For dynamical simulations, we use the static model developed for 3D confluent
tissues Merkel and Manning [2018], but with added cellular activity. The model
essentially allows active fluctuations during time evolution of cells. Cells
are polyhedra derived from a space-filling Voronoi tessellation of the
periodic simulation box. The energy functional given by Eq. 1, provides the
mechanical force $\mathbf{f}_{i}=\partial e/\partial\mathbf{r}_{i}$ on cells
due to changes in cell shape and/or shared interface area. For introducing
activity in this dynamics, one needs to choose the frame of reference – for
example a static extra-cellular fluid. In a self-propelled system, the
$i_{th}$ cell has a polarization vector $\hat{\mathbf{n}}_{i}$ and exerts an
active force of $v_{0}/\mu$ on the static media, where $\mu$ is the mobility
of the cell and $v_{0}$ is the self-propulsion speed. Setting $\mu=1$, the
dynamical equation is :

$\dfrac{d\mathbf{r}_{i}}{dt}=v_{0}\hat{\mathbf{n}}_{i}+\mathbf{f}_{i}.$ (5)

The polarization vector evolves via white Gaussian noise on a unit sphere with
diffusion coefficient $D_{r}$, or

$\dfrac{d\hat{\mathbf{n}}_{i}}{dt}=\sqrt{2D_{r}}\big{(}\mathbf{E}-\hat{\mathbf{n}}_{i}\hat{\mathbf{n}}_{i}\big{)}\dot{\xi}_{i},$
(6)

where $\mathbf{E}$ is the $3\times 3$ identity tensor, the dyadic product of
the polarization vector with itself is given by
$\hat{\mathbf{n}}_{i}\hat{\mathbf{n}}_{i}$ and $\xi_{i}$ is the white Gaussian
noise with $\langle\xi_{i}\rangle=0$ and
$\langle\xi_{i}\rangle(t)\langle\xi_{j}\rangle(t^{\prime})=\delta_{ij}\delta(t-t^{\prime})\mathbf{E}$.

Timescales in our dynamic equations are nondimensionalized by the natural time
unit $\tilde{t}=1/K_{V}V_{0}^{4/3}$, which is unity for our choice of
parameter values. We use a self-propulsion speed of $v_{0}=0.1$ and a
rotational diffusion coefficient $D_{r}=1.0$. For such high values of
rotational diffusion, the transition to a Brownian regime happens rather
quickly, i.e $t/\tilde{t}>1/D_{r}$. In the absence of any cell-cell
interactions, a self-propelled particle with that choice of parameters would
diffuse over a characteristic self-diffusivity timescale $100\tilde{t}$
Marchetti et al. [2013]. However, cell-cell interactions slow down diffusion,
and we find that in our model the characteristic diffusion timescale for a
fixed shape index of 5.5 is $\tau_{s}^{0}\equiv 1000\tilde{t}$. Therefore, for
simulations in the presence of fluctuations, we use an integration step size
$\Delta t=0.01\tilde{t}$ and a time-span $T_{SS}$ of $1\times 10^{5}$ steps
which is sufficient to allow a system of this size to achieve a steady state.
While this holds true for the entire range of the explored tension, values in
the lower-moderate regime ($\sigma<1$), achieve a steady state sooner, i.e. in
just $1\times 10^{4}$ 10000 steps. In the sorting simulations, we extend this
timescale even further to allow cells to travel across the length of the
simulation box more than 60 times, ensuring that even long-timescale
patterning processes are allowed.

###  S2 Maximum demixing value with respect to system-size

For a small system-size, the demixing cannot attain the maximum possible value
of unity as the number of heterotypic facets is not negligible as compared to
the number of homotypic facets. Therefore we look at demixing in a segregated
arrangement as a function of system-size.

The simple assumption that a given interfacial cell shares only one-third of
its facets with heterotypic neighbors, i.e. $DP_{final}\propto N^{1/3}$ is in
fairly good agreement with the observed final demixing values seen in our
simulations as shown in Fig S1.

Figure S1: Plot for the observed final demixing parameter in simulations
($DP_{final}$) as a function of the cubic root of the system size. Data is
consistent with a straight line.

The small deviations from this prediction we see in simulations are not
surprising; even though the bilayer geometry is the least-surface-area
configuration for a cubic simulation space like ours, for a mixed initial
conditions with periodic boundary conditions it is actually difficult to
achieve this configuration. In general, compartments with different topologies
can be observed – for example with one or two holes. Many of these dynamically
transform into the most-stable planar compartment that has no holes.

###  S3 Increase in the average observable facet area

We observe that with HIT, heterotypic facet areas become bimodal i.e. there
are a majority of vanishingly small faces and the rest are larger faces. From
an experimental point of view, vanishingly small facets are difficult to
detect and hence, we plot the average of larger facets in Fig. S2. We see that
it increases directly with tension.

Figure S2: (a) Plot for average observable interfacial facet area $\langle
a\rangle$ as a function of tension $\sigma$.

###  S4 Change in 2D registry leads to cusp-like response

For a perfect prism-like geometry along 3D interfaces, we need cells to be at
the same distance from the interface as well as the cell pairs across the
interface must have their centers aligned, i.e. registered. While the former
condition is already understood using simple square lattice geometry
calculation in a recent work Sussman et al. [2018], the latter remains to be
studied. Hence we use the exact same setup with 9 neighboring cells, but with
a different local perturbation. As shown in Fig. S3, we consider a small
perturbation of the central cell along the interface by an amount $\epsilon$,
redraw the Voronoi tessellation and calculate final energy.

Figure S3: A sketch of the 9-cell set up where the different cell types are
colored with pink and blue and the heterotypic boundary highlighted in red.
The Voronoi centers are shown in grey-filled circles . The black edges show
the initial cell-cell boundaries. The central cell’s boundary after the
perturbation is highlighted in blue.

The initial energy of the 9-cell system is:

$E_{i}=9k_{A}(a-a_{0})^{2}+9(p-p_{0})^{2}+3\gamma_{0}l_{ij},$ (7)

where $\gamma_{0}$ is the interfacial tension, $l_{ij}$ is the length of edges
shared between unlike cells. By making the assumptions of $s_{0}=1$ and
$p_{0}=4$ to start from the ground state, we reduce the expression to:

$E_{i}=3\gamma_{0},$ (8)

After displacing the central cell to the right, the new total energy to
leading order becomes-

$\displaystyle E_{f}=3\gamma_{0}+\gamma_{0}(\sqrt{2}-1)\epsilon_{x}$ (9)
$\displaystyle+\bigg{(}\frac{3k_{A}}{2}+(22-14\sqrt{2}+\gamma_{0}\big{(}\frac{20\sqrt{2}-8}{32\sqrt{2}}\big{)})\bigg{)}\epsilon_{x}^{2},$

We see that to linear order, the energy for a perturbation parallel to the
interface has exactly the same cuspiness (e.g. linear scaling with distance of
perturbation $\epsilon_{x}$) that was observed for a perturbation
perpendicular to the interface. Moreover, the coefficient is
$\gamma_{0}\epsilon^{2}$, which is twice that of coefficient for a
perpendicular perturbation Sussman et al. [2018], indicating the non-linear
stiffness is higher for a parallel (or sideways) perturbation.

###  S5 Change in 3D registry creates several new HIT facets

The calculation in the previous section establishes that a cell adjacent to a
high-tension interface has a cusp-like response to both kinds of perturbations
in 2D: perpendicular as well as parallel to the interface.

Figure S4: A hexagonal right prism is perturbed both (a) perpendicular and (c)
parallel to the shared horizontal interface at the top, as depicted by the
blue arrows. The Voronoi cell centers are shown in blue circles. This gives
rise to several new facets that can be seen from the top view on the right sub
figures in (b) and (d). The original facet is labelled as ‘old neighbor’ and
the new facets are labelled as ‘new’ and highlighted by dark edges. The
perpendicular perturbation creates 6 new facets as shown in (b). The parallel
perturbation creates 8 new facets as shown in (d). As the system is no longer
symmetric about z-axis, three of the neighboring cells are colored differently
to emphasize which cell the new facet belongs to. The perturbed cell is shown
in pink.

A qualitative similarity is only expected because the same number of new edges
is created along the high tension cable for both kinds of perturbations. In
order to extend this observation to 3D tissues one needs to check if new
facets are being created along the high-tension interface, the same needs to
be verified for a cell adjacent to a high-tension interface. To explore this
we create a Voronoi tessellation of a group of hexagonal right prisms, stacked
about their hexagonal faces. The isotropic hexagonal face of unit area has a
2D shape index of 3.72 and a the height of the prism is set to unity. This
results in a 3D shape index of 5.72. In a manner very analogous to the 2D
calculation performed in the previous section, we push a central prism both
perpendicularly and parallely to the interface as shown in Fig. S4(a,c). The
parallel perturbation shown here refers to perturbing the cell along the
lateral direction that connects two neighboring centers.

As one would expect, the perpendicular perturbation leads to formation of 6
new facets as highlighted in the top view shown in Fig. S4(b). The parallel
perturbation, on the other hand, leads to 8 new facets. Therefore, cusp-like
behaviour is expected to accompany changes in registry in 3D systems in both
cases.

###  S6 BCC toy model

We further explore the properties of both kinds of minima in Fig S6. We also
look at the contribution from shape frustration and surface energy to the
total energy change of the system in Fig. S5.

Figure S5:  Shape versus surface energies at extreme tension values : Plot of
change in total energy (in dashed orange), blue cell’s mechanical energy (in
solid blue) and interfacial energy (in solid black) with respect to registry
is shown for (a) $\sigma=0.001$ and (b) $\sigma=1$.  Figure S6: Properties of
parabolic and cuspy minima : (a) Plot of stiffness (magnitude of curvature) of
the parabolic minima with respect to increasing tension. (b) The restoring
force ($F_{R}$) for perturbations away from the cuspy minima is independent of
the displacement on y axis for different values of tension ($\sigma$), except
for the very largest displacements. (c) The value of restoring force at a
fixed value of displacement is plotted as a function of tension, exhibiting a
linear relationship.

###  S7 HCP toy model

For the HCP toy model, we explore the shared surface area and the system’s
energy profiles as a function of increased sigma. For shapes near transition,
the behaviour is qualitatively very similar to the BCC toy model. The shared
surface area is a minimum for complete registry as shown in Fig S7(a).
However, the ground state registry is unity only for very high tension values.
For lower values, it depends on the cell shape. For a more fluid-like shapes
the effect of shape frustration becomes negligible and the state prefers being
registered for almost all values of tension.

Figure S7:  Global minima becomes registered for higher tension: (a) Shared
surface area between both interfacial layers is computed as a function of the
registration between the layers. The layer of blue cells is allowed to move
past the string of layer of pink cells. Snapshots for no, half and complete
registration is shown for encircled points. For better visualization, one cell
from each layer is highlighted using solid-edge triangulation. (b) The change
in the total energy of this system is plotted with respect to registry, for
different values of tension $\sigma=0.001,0.003,0.010,0.031,0.10,0.31$ shown
in pink to blue, for a fixed shape value of 5.5.

###  S8 Cusp-like restoring forces in bulk systems

For the second part of this verification process i.e.,to check if the
restoring forces are Hookean or non-Hookean in nature, one needs to analyse
the qualitative scaling of the restoring forces with respect to perturbations
about the ground state. This computation requires the fluid-like systems to be
at mechanical equilibrium, though this is difficult in our simulations as we
describe below. We let the steady state configurations found above to relax at
zero fluctuations ($v_{0}=0$) using standard steepest descent dynamics. The
integration time step is reduced to $dt=0.001$ and we integrate over a longer
time window of $1\times 10^{6}$ steps for a system-size $N=216$. With these
choices, the noise floor, i.e. the magnitude of fluctuations in force that
quantify how far the configuration is from the local minimum in the energy at
the end of this many relaxation steps, reduces to a value of approximately
$10^{-4}$. Unfortunately, we can not go below this value without major
alterations to our code, because inhomogeneity in tensions in confluent models
gives rise to higher fold coordinated vertices. In the current version of the
3D code, such configurations are not allowed and are thus computationally
inaccessible. In order to decrease the number of such instances, we decrease
the distance tolerance of the Voronoi algorithm by an order of magnitude to
$10^{-12}$. Nevertheless, we still see a significant fraction of instances
(e.g. for $\sigma=0.16$, $\sim 65\%$ of the minimization procedures survive
and for $\sigma=0.32$ $\sim 20\%$ survive), and the number of such instances
increases if we allow the system to relax for more time steps. Therefore, the
choices for this minimization procedure have been chosen to balance between
minimizing the number of many-fold-coordinated vertices while still having a
small enough noise floor to be able to study truly minimized states. At the
end of every successful minimization procedure, we then record the response of
an interfacial cell to discrete changes in its cell-cell registry denoted by
$\epsilon$.

For the cases that survive, a randomly chosen interfacial cell is perturbed
away from being registered by a discrete amount $\epsilon$ and we record the
restoring force $f_{R}$ on the cell while keeping every other cell intact. The
distribution of these forces is shown as a function of increasing tension in
Fig. S8 for the lowest value of perturbation magnitude $\epsilon=10^{-6}$.
Similar plots exist for higher values of $\epsilon$. For all cases, we observe
a bimodal distribution of forces. The lower peak in the distribution
corresponds to forces below the noise floor; in other words, these forces are
consistent with systems that have not equilibrated even after $10^{6}$ time
steps, and so our probe displacement of a single cell does not probe the
energy minimum. The higher peak is above the noise floor and represents real
measurements of forces in response to perturbations away from a state that is
truly in mechanical equilibrium. Therefore, we focus on the statistics of this
higher-force peak in what follows. We expect that larger tensions are able to
help equilibrate the system more effectively, shifting more configurations to
the higher force bump. That is what is observed in our simulations: for
$\sigma=0.16$, only $\sim 15\%$ of the systems are on the right-hand-side
Gaussian, while for $\sigma=0.32$, it rises to $\sim 50\%$. The scaling of the
average force in this higher peak (according to the best-fit Gaussian) is then
shown as a function of perturbation magnitude $epsilon$. It is non-Hookean as
shown in the inset of Fig. S8. Therefore, the computationally accessible
ground states have a cusp-like behaviour.

Figure S8:  Cusp-like restoring forces observed at high tension : The
probability distribution of restoring forces $f_{R}$ for the lowest
$\epsilon=10^{-6}$ is shown as a function of tension values increasing from
0.16 (light blue) to 0.32 (deep blue) in steps of 0.04. In the inset, the mean
value of $f_{R}$ for the second Gaussian (denoted by *) on the right is
plotted for increasing $\epsilon$ for $\sigma=0.16$ in open circles and
$\sigma=0.32$ in asterisk. The dashed line a Hookean fit to the large
displacement data.


# DENT-DDSP: Data-efficient noisy speech generator using differentiable
digital signal processors for explicit distortion modelling and noise-robust
speech recognition

###### Abstract

The performances of automatic speech recognition (ASR) systems degrade
drastically under noisy conditions. Explicit distortion modelling (EDM), as a
feature compensation step, is able to enhance ASR systems under such
conditions by simulating the in-domain noisy speeches from the clean
counterparts. Yet, existing distortion models are either non-trainable or
unexplainable and often lack controllability and generalization ability. In
this paper, we propose a fully explainable and controllable model: DENT-DDSP
to achieve EDM. DENT-DDSP utilizes novel differentiable digital signal
processing (DDSP) components and requires only 10 seconds of training data to
achieve high fidelity. The experiment shows that the simulated noisy data from
DENT-DDSP achieves the highest simulation fidelity compared to other baseline
models in terms of multi-scale spectral loss (MSSL). Moreover, to validate
whether the data simulated by DENT-DDSP are able to replace the scarce in-
domain noisy data in the noise-robust ASR tasks, several downstream ASR models
with the same architecture are trained using the simulated data and the real
data. The experiment shows that the model trained with the simulated noisy
data from DENT-DDSP achieves similar performances to the benchmark with a 2.7%
difference in terms of word error rate (WER). The code of the model is
released online111https://github.com/guozixunnicolas/DENT-ddsp.

Index Terms: explicit distortion modelling, DDSP, noise-robust ASR

##  1 Introduction

The end-to-end automatic speech recognition (ASR) systems have achieved a
remarkable success, yet, the systems are prone to noisy conditions. Several
datasets [1, 2, 3, 4] containing clean speeches, noisy speeches and the
corresponding transcribed texts have been collected and utilized to boost ASR
performances under noisy conditions. Yet, such datasets are usually scarce.
Explicit distortion modelling (EDM) [5], as an alternative, is able to
simulate noisy speeches from their clean counterparts. However, distortion
models with high simulation fidelity are hard to be obtained.

To achieve EDM, traditional methods such as parallel model combination (PMC)
[6] and vector taylor series (VTS) [7, 8, 9] are widely used to estimate the
distortion model before the deep learning era. An intuitive method is proposed
in [10] where static noise within the limited collected noisy data are
aggregated and added to the clean speeches. Another method [11] adds pre-
recorded noise to the clean speeches and passes them through pre-defined
codecs to obtain the simulated noisy data. These two static methods [10, 11],
however, are untrainable thus cannot be adapted to other circumstances.
Recently, GAN-based methods [12, 13] are utilized to achieve EDM and have
obtained promising results. We find SimuGAN [12] the closest match to our
work. It operates directly on magnitude spectrograms and uses GAN and
contrastive learning methods to distort the clean spectrograms. Yet, GAN-based
methods contain large amounts of parameters making the distortion models
unexplainable and uncontrollable.

Recently with the advent of the differentiable digital signal processing
(DDSP) [14], traditional digital signal processors (DSP) become differentiable
and trainable and have seen success in the field of speech synthesis [15, 16].
Inspired by the success of DDSP in the synthesis domain, we find DDSP a viable
solution to achieve EDM. In this paper, we focus on simulating the VHF/UHF
transmitted data (e.g., air traffic control). By comparing the characteristics
of the VHF/UHF transmitted data with their clean counterparts, we observe that
the VHF/UHF transmitted data are distorted, compressed to a fixed dynamic
range, equalized and also contain colored noise. We hence find DDSP capable of
achieving such exact conversion.

In this work, we propose a data-efficent noisy speech generator using novel
DDSP components (DENT-DDSP) to simulate the VHF/UHF transmitted data. To our
best knowledge, DENT-DDSP is the very first model that utilizes DDSP to
achieve explicit distortion modelling. It only requires 10 seconds of parallel
training data, yet, the trained distortion model achieves high fidelity
simulation. As a result, the VHF/UHF transmitted data can easily be simulated
from clean speeches. The simulated data can be further used for other
downstream tasks (e.g., noise-robust ASR, speech enhancement).

More specifically, DENT-DDSP consists of two signal chains: an audio signal
chain and a noise signal chain. Each signal chain consists of trainable DDSP
components connected in series. The audio signal chain is able to distort,
compress and equalize the clean speeches. The noise signal chain is able to
transform the input white noise to the desired spectral characteristics. The
outputs from the two signal chains are added to form the final simulated noisy
speech. In the audio signal chain, two novel DDSP components are proposed:
waveshaper and computionally-efficient dynamic range compressor (DRC). The
experiment has shown that the proposed computationally-efficient DRC
drastically improves the computational efficiency while maintaining the
simulation fidelity.

Extensive experiments have shown that DENT-DDSP achieves the highest
simulation fidelity among all distortion models in terms of MSSL. Moreover, it
has shown strong generalization abilities over unseen training data and
outperforms other GAN-based models. To further validate whether the labor-
intensive way of collecting the VHF/UHF transmitted data can be replaced by
simulating these data from DENT-DDSP, a downstream noise-robust ASR model
trained using simulated data from DENT-DDSP is benchmarked against a same
model trained with real noisy data. The experiment has shown that the model
trained with simulated data from DENT-DDSP has achieved similar WER to the
benchmark model with 2.7% difference and has outperformed the best baseline
model by 7.3% in terms of WER.

##  2 Model Architecture

Figure 1: DENT-DDSP architecture. DRC stands for dynamic range compressor and
EQ stands for equalizer. The parameters of the DDSP components are trainable
and controllable.

The model architecture of DENT-DDSP is shown in Figure 1. DENT-DDSP contains 2
parallel signal chains. The audio signal chain contains a waveshaper, a
computationally-efficient DRC and an equalizer and receives clean speech:
$s_{in}(t)$ as input. The noise signal chain contains an equalizer and
receives white noise: $n_{in}(t)$ as input. In each signal chain, the input
signal will pass through each DDSP component in series. The simulated noisy
speech: $s_{simulated}(t)$ will be formed by adding the audio signal chain
output: $s_{out}(t)$ and the weighted noise signal chain output:$n_{out}(t)$
following Equation 1. $\lambda$ is a weighting and non-trainable parameter and
is set to $1$ during training. However, it can be adjusted during the
generation phase to simulate noisy data with different SNRs in order to
achieve data augmentation. A spectral loss function: MSSL (multi-scale
spectral loss) [14] will be calculated between the simulated noisy speech:
$s_{simulated}(t)$ and the real noisy data to update the model’s parameters
via backpropagation. We calculate MSSL using the following FFT sizes: 2048,
1024, 512, 256, 128, 64, to reflect the spectral distance in different
spectral resolutions.

$s_{simulated}(t)=s_{out}(t)+\lambda\cdot n_{out}(t)$ (1)

We then formulate the DDSP components in the signal chains and introduce their
functionalities. For simplicity, we represent each component’s input as $x(t)$
and output as $y(t)$. $x_{dB}(t)$ and $y_{dB}(t)$ represents the input and
output in the $dB$ scale.

###  2.1 Waveshaper

The waveshaper is able to introduce a distortion effect to the input signal as
shown in Equation 2. The distortion gain $g_{distort}$ represents the amount
of distortion added to the input signal. With $g_{distort}$ close to $0$, the
transfer function of a waveshaper is almost linear hence little distortion is
applied. However, with an increasing $g_{distort}$, the transfer function
becomes non-linear and saturates quickly. Such non-linearity introduces the
distortion to the input signal. The distortion gain $g_{distort}$ is set to be
a trainable parameter.

$y(t)=\frac{2}{\pi}\arctan(g_{distort}\cdot\frac{\pi}{2}\cdot
x(t))\;\;,\;\;g_{distort}>0$ (2)

###  2.2 Computationally-efficient DRC

Dynamic range compressors(DRC) [17] are widely used in music production to
limit the dynamic ranges of different audio tracks. By observation, the
VHF/UHF transmitted data are compressed to a fixed dynamic range in the time
domain. We hence find DRC a suitable DDSP component to achieve such limiting
behaviour. A detailed survey [17] has introduced different kinds of DRCs and
we have chosen a hard-knee DRC as our DDSP backbone. We then formulate our
novel computationally-efficient DRC in Equation 3-7.

Equation 3 calculates the amplitude reduction gain $g(t)$ based on the
compression threshold $T$ and the compression ratio $R$. Equation 4-6 shows
the proposed efficient implementation of gain smoothing using a novel
companding operation. The amplitude reduction gain $g(t)$ is firstly
downsampled via a linear interpolation by a downsampling factor $ds\\_factor$
in Equation 4. The downsampled gain function $g_{d}(t)$ is then smoothed by
$\alpha_{A}$ and $\alpha_{R}$ which represents attack and release time
constant respectively in Equation 5. The smoothed gain function $g_{ds}(t)$ is
then upsampled with overlapping hanning windows by an upsampling factor
$us\\_factor$ in Equation 6. $ds\\_factor$ and $us\\_factor$ are set to be
equal such that $g_{es}(t)$ will have the same shape as $g(t)$. $g_{ds}(t)$ is
the main performance bottleneck in gain calculation since the operation is
auto-agressive and recursive with time complexity of $O(L/ds\\_factor)$ where
$L$ is the total number of audio samples. Hence, the time complexity will
increase linearly with a decreasing $ds\\_factor$. We will later prove the
companding operation in Equation 4 and 6 drastically improve the efficiency
while maintaining the simulation quality in Section 3.3. Finally, the output
$y_{dB}(t)$ is obtained in Equation 7 with an additional makeup gain
$g_{makeup}$. The following parameters are trainable: $T$, $R$, $\alpha_{A}$,
$\alpha_{R}$, $g_{makeup}$.

$g(t)=\begin{cases}0&\text{$x_{dB}(t)\leq T$}\\\
\displaystyle\frac{x_{dB}(t)-T}{R}&\text{$x_{dB}(t)>T$}\\\ \end{cases}$ (3)
$g_{d}(t)=downsample(g(t),ds\\_factor)$ (4)
$g_{ds}(t)=\begin{cases}\alpha_{A}g_{ds}(t-1)+(1-\alpha_{A})g_{d}(t)&\text{$g_{d}(t)>g_{ds}(t-1)$}\\\
\alpha_{R}g_{ds}(t-1)+(1-\alpha_{R})g_{d}(t)&\text{$g_{d}(t)\leq
g_{ds}(t-1)$}\\\ \end{cases}$ (5) $g_{es}(t)=upsample(g_{ds}(t),us\\_factor)$
(6) $y_{dB}(t)=x_{dB}(t)+g_{es}(t)+g_{makeup}$ (7)

###  2.3 Equalizer(EQ)

By observation, the VHF/UHF transmitted data are band-limited and equalized in
the frequency domain. To achieve this, we adopt the equalizer (EQ)
implementation from [14] where the EQ is implemented using a linear time
invariant FIR filter. For each EQ, the trainable parameters are set to be the
magnitude of the filter frequency response: $FR_{mag}$ and the number of the
frequency bins are set to be 1000. Additionally, for the EQ in the noise
signal chain, the noise amplitude is made trainable in order to adjust the
volume of the filtered noise.

##  3 Experiments and results

###  3.1 Dataset

The Robust automatic transcription of speech (RATS) project [1] has collected
a parallel dataset which contains clean speeches, the corresponding VHF/UHF
transmitted noisy speeches and the transcribed texts. To obtain the dataset,
pre-recorded conversational speeches are broadcasted over 8 radio channels and
the transmitted noisy audio are captured concurrently. We select Channel A
from the RATs dataset which contains 57.4 hours of data for training and
testing. To address the data scarcity problem mentioned in Section 1, only
less than 60 seconds of parallel audio are selected as training data. During
training, the data are batched into 1-second chunks.

###  3.2 Experiment setup

####  3.2.1 Training data selection and model comparison

To obtain the best performing and most efficient distortion model, several
DENT-DDSP models are trained and compared using different amounts of parallel
training data from 10 to 60 seconds with various speech to total ratio: $s2t$.
$s2t$ is defined as the duration of active speech over the total duration in
each 1-second chunk. Active speech duration is calculated based on the clean
audio energy in a sliding window with a threshold of $50dB$. By analyzing the
dataset, with $s2t<0.4$, the data contains mostly non-speech audio or silence.
Hence, we only select data with $0.4\leq s2t<1.0$. To acquire the desired
training data, the 1-second data chunks with the desired $s2t$ are randomly
chosen and aggregated to the desired total duration. To evaluate the
simulation fidelity, testing clean speeches are input to the trained
distortion model and the MSSL[14] is calculated between the real and simulated
noisy data. We compare the proposed DENT-DDSP with the following baseline
distortion models in terms of simulation fidelity:

  * •

Clean augment model[10]: simulated noisy speeches are obtained by adding the
aggregated stationary noise from RATs Channel A to the clean speeches.

  * •

G.726 augment model[11]: The clean speeches are first passed through the G.726
codec. The aggregated stationary noise from RATs Channel A is added
subsequently to the codec outputs.

  * •

Codec2 augment model[11]: Codec choice in the G.726 augment model is replaced
to Codec2 with 700 bits/s.

  * •

SimuGAN[12]: simulated noisy speeches are obtained by feeding the clean
speeches to a trained SimuGAN.

####  3.2.2 Noise-robust automatic speech recognition (ASR)

A downstream noise-robust ASR task is performed to validate the effectiveness
of DENT-DDSP. The 57.4-hour parallel data from Channel A are split into 3
folds [18]: 44.3-hour data for training; 4.9-hour data for validation and
8.2-hour data for testing. We adopted the conformer-based [19] dual-path ASR
system from [12] with 12 conformer layers in the encoder and 6 transformer
layers in the decoder. Such dual-path ASR architecture is specifically
designed for noise-robust ASR [20, 21]. We pre-train a language model [22]
with 2 RNN layers using the existing transcribed texts and utilize it during
text decoding. To benchmark the performance, a benchmark model is firstly
trained with real noisy data. To address the data scarcity problem
aforementioned in the real world, the benchmark model only incorporates the
single-path setting in [12] where only real noisy data and the transcribed
text are input to the network . Another set of model is trained using the
simulated noisy data from various distortion models using the dual-path
setting in [12]. These sets of models will be tested on a common test set
consists of real noisy data and the performances in terms of WER will be
compared.

###  3.3 Results

####  3.3.1 Effect of different amounts of training data with various $s2t$
and model comparison

Table 1 shows the MSSL and between the real and simulated audio using
different amounts of training data with various $s2t$. In general, MSSL
decreases with fewer amounts of training data and increasing $s2t$. The best
DENT-DDSP model is hence obtained using 10 seconds of training data with
$0.8\leq s2t<1.0$. Comparing DENT-DDSP with other distortion models, we
observe that DENT-DDSP is trainable, contains much fewer trainable parameters
and requires much fewer amount of training data. Moreover, noisy audio
simulated by DENT-DDSP has the lowest MSSL which indicates DENT-DDSP has the
highest simulation fidelity among all baseline models. Readers are strongly
encouraged to listen to the simulated noisy speeches and compare them with the
real noisy data from222https://guozixunnicolas.github.io/DENT-DDSP-demo/. To
visualize the simulation fidelity, example spectrograms of the clean audio,
audio simulated from various distortion models and the real noisy data are
plotted in Fig 2.

Table 1: Simulation quality of DENT-DDSP using different amounts of training
data with various s2t and model comparison

model no. of trainable param. s2t amount of training data MSSL DENT-DDSP
(ours) 2k $[0.8,1,0)$ 60-sec parallel 0.184 $[0.6,0.8)$ 60-sec parallel 0.187
$[0.4,0.6)$ 60-sec parallel 0.186 $[0.8,1,0)$ 40-sec parallel 0.177
$[0.6,0.8)$ 40-sec parallel 0.187 $[0.4,0.6)$ 40-sec parallel 0.188
$[0.8,1,0)$ 20-sec parallel 0.170 $[0.6,0.8)$ 20-sec parallel 0.189
$[0.4,0.6)$ 20-sec parallel 0.189 $[0.8,1,0)$ 10-sec parallel 0.170
$[0.6,0.8)$ 10-sec parallel 0.181 $[0.4,0.6)$ 10-sec parallel 0.189 SimuGAN
14M - 10-min unparallel 0.173 Clean augment non-trainable - - 0.192 G.726
augment non-trainable - - 0.197 Codec2 augment non-trainable - - 0.197

Figure 2: Magnitude spectrogram comparison between clean audio, audio
simulated from various distortion models and the real noisy data.

####  3.3.2 Effectiveness of the companding operation in the proposed
computationally-efficient DRC

Table 2 shows the effectiveness of the proposed companding operation in the
novel computationally-efficient DRC. With decreasing $ds\\_factor$s using the
best DENT-DDSP obtained, the training time increases linearly, however,
without obvious improvements on the MSSL. This justifies the effectiveness of
the companding operation which is able to increase computational efficiency
while maintaining the generation quality.

Table 2: Training time and MSSL comparison among DENT-DDSPs with different
$ds\\_factor$s

ds factor training time MSSL $16$ 7min 0.170 $8$ 13.5min 0.169 $4$ 25.5min
0.170 $2$ 60.2min 0.169

####  3.3.3 Generalization ability

To reflect the generalization ability of the trainable distortion models:
DENT-DDSP and SimuGAN over different types of audio, non-speech audio are used
to test the distortion models. The non-speech audio (available onlinefoot2)
include ambient noise, guitar and piano music and synthesized audio. The
spectrograms of the testing audio and the simulated noisy audio from the two
models are shown in Figure 3. By observation, data simulated by DENT-DDSP have
consistent noise spectrums and the audio are distorted to the desired spectral
behaviour: the spectrogram becomes blurred with certain frequencies being
dampened. Yet, data simulated by SimuGAN fail to produce such consistent
distortion behaviour. Moreover, data simulated by SimuGAN either contain
artefacts or become unintelligible. Hence, DENT-DDSP shows stronger
generalization ability compared to SimuGAN due to the use of explainable DDSP
components.

Figure 3: Noisy data simulated by DENT-DDSP and SimuGAN using unseen non-
speech data

####  3.3.4 Benchmarking the noise-robust ASR task

To compare against the benchmark model, several additional models are obtained
as follows: Data simulated by DENT augment model is obtained by setting
$\lambda=0.79,1.26$ in Equation 1 in a trained DENT-DDSP model so that the
noise gain of the $n_{out}(t)$ become $\pm 2dB$; Additionally, a clean model
is trained only using clean speeches. Table 3 shows the WER of these models
tested on a common test set consists of real noisy data. By observation, the
DENT augment model achieves the closest WER compared to the benchmark model
with a difference of 2.7% and outperforms the best baseline model: SimuGAN by
7.3%. This proves that DENT-DDSP is able to simulate noisy data with similar
characteristics to the real noisy data and offers an alternative to simulate
the VHF/UHF transmitted data thus greatly reducing the effort to collect these
data manually on a large scale.

Table 3: WER comparison using same ASR model with simulated data from
different distortion models and real noisy data

model noisy speech source WER Clean - 93.4% Clean augment clean speech +
aggregated noise 73.8% G.726 augment g726 codec speech + aggregated noise
75.6% Codec2 augment codec2 speech + aggregated noise 83.2% SimuGAN SimuGAN
simulated 65.9% DENT-DDSP (ours) DENT simulated 66.4% DENT augment (ours) DENT
simulated with noise gain $=\pm 2dB$ 58.6% Benchmark real noisy data 55.9%

##  4 Conclusion

We propose a fully explainable and controllable model DENT-DDSP based on novel
DDSP components which only requires 10 seconds of training data to achieve
explicit distortion modelling. Besides the existing DDSP components, we
propose 2 novel DDSP components: waveshaper and computationally-efficient DRC
which can be easily integrated to other DDSP models. Experiments have shown
that data simulated by DENT-DDSP achieve the highest simulation fidelity among
all distortion models. Moreover, the noise-robust ASR system trained with the
data simulated by DENT-DDSP achieves similar WER compared to the benchmark
trained with real noisy data with a 2.7% difference in terms of WER. It has
also achieved 7.3% WER improvement over the best baseline distortion model.

## References

  * [1] D. Graff, K. Walker, S. Strassel, X. Ma, K. Jones, and A. Sawyer, “The RATS collection: Supporting HLT research with degraded audio data,” in _Proc. of the 9th Int. Conf. on Lang. Resour. and Eval. (LREC)_ , 2014, pp. 1970–1977.
  * [2] S. Badrinath and H. Balakrishnan, “Automatic speech recognition for air traffic control communications,” _Transportation Research Record_ , vol. 2676, no. 1, pp. 798–810, 2022.
  * [3] T. Pellegrini, J. Farinas, E. Delpech, and F. Lancelot, “The airbus air traffic control speech recognition 2018 challenge: towards atc automatic transcription and call sign detection,” in _Proc. Int. Speech Commun. Assoc.(Interspeech)_ , 2018, pp. 2993–2997.
  * [4] J. Godfrey, E. Holliman, and J. McDaniel, “Switchboard: telephone speech corpus for research and development,” in _Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_ , vol. 1, 1992, pp. 517–520 vol.1.
  * [5] J. Li, L. Deng, Y. Gong, and R. Haeb-Umbach, “An overview of noise-robust automatic speech recognition,” _IEEE ACM Trans. on Audio, Speech, and Lang. Process._ , vol. 22, no. 4, pp. 745–777, 2014.
  * [6] M. J. F. Gales, “Model-based techniques for noise robust speech recognition,” in _Ph.D. thesis, University of Cambridge_ , 1995.
  * [7] Y. Minami and S. Furui, “A maximum likelihood procedure for a universal adaptation method based on hmm composition,” in _Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_ , vol. 1, 1995, pp. 129–132.
  * [8] A. Sankar and C.-H. Lee, “Robust speech recognition based on stochastic matching,” in _Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_ , vol. 1, 1995, pp. 121–124.
  * [9] A. Acero, l. Deng, T. Kristjansson, and J. Zhang, “Hmm adaptation using vector taylor series for noisy speech recognition,” 2000, pp. 869–872.
  * [10] D. Ma, G. Li, H. Xu, and E. S. Chng, “Improving code-switching speech recognition with data augmentation and system combination,” in _Proc. Asia-Pacific Signal and Inf. Process. Assoc. Annu. Summit and Conf. (APSIPA ASC)_ , 2019, pp. 1308–1312.
  * [11] M. Ferràs, S. Madikeri, P. Motlicek, S. Dey, and H. Bourlard, “A large-scale open-source acoustic simulator for speaker recognition,” _IEEE Signal Processing Letters_ , vol. 23, no. 4, pp. 527–531, 2016.
  * [12] C. Chen, N. Hou, Y. Hu, S. Shirol, and E. S. Chng, “Noise-robust speech recognition with 10 minutes unparalleled in-domain data,” in _Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_ , 2022, pp. 4298–4302.
  * [13] H. Hu, T. Tan, and Y. Qian, “Generative adversarial networks based data augmentation for noise robust speech recognition,” in _Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_ , 2018, pp. 5044–5048.
  * [14] J. H. Engel, L. Hantrakul, C. Gu, and A. Roberts, “DDSP: differentiable digital signal processing,” in _in Proc. Int. Conf. on Learn. Representations(ICLR)_ , 2020.
  * [15] G. Fabbro, V. Golkov, T. Kemp, and D. Cremers, “Speech synthesis and control using differentiable dsp,” 2020. [Online]. Available: https://arxiv.org/abs/2010.15084
  * [16] J. Alonso and C. Erkut, “Latent space explorations of singing voice synthesis using ddsp,” in _Proc. of 18th Sound and Music Computing Conf._ , 2021.
  * [17] D. Giannoulis, M. Massberg, and J. Reiss, “Digital dynamic range compressor design—a tutorial and analysis,” _J. of the Audio Eng. Soc.(AES)_ , vol. 60, pp. 399–408, 2012.
  * [18] C. Chen, Y. Hu, N. Hou, X. Qi, H. Zou, and E. S. Chng, “Self-critical sequence training for automatic speech recognition,” in _Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_ , 2022, pp. 3688–3692.
  * [19] A. Gulati, J. Qin, C. Chiu, N. Parmar, Y. Zhang, J. Yu, W. Han, S. Wang, Z. Zhang, Y. Wu, and R. Pang, “Conformer: Convolution-augmented transformer for speech recognition,” in _Proc. Int. Speech Commun. Assoc.(Interspeech)_ , 2020, pp. 5036–5040.
  * [20] Y. Hu, N. Hou, C. Chen, and E. S. Chng, “Interactive feature fusion for end-to-end noise-robust speech recognition,” in _Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_ , 2022, pp. 6292–6296.
  * [21] Y. Hu, N. Hou, C. Chen, and E. Chng, “Dual-path style learning for end-to-end noise-robust speech recognition,” in _preprint arXiv2203.14838_ , 2022.
  * [22] S. Karita, N. Chen, T. Hayashi, T. Hori, H. Inaguma, Z. Jiang, M. Someki, N. E. Y. Soplin, R. Yamamoto, X. Wang, S. Watanabe, T. Yoshimura, and W. Zhang, “A comparative study on transformer vs rnn in speech applications,” _IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_ , Dec 2019.