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

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.

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 [27][28][29][30][31]. However, an essential feature of cell sorting that is observed in experimental cocultures is that the interface is much sharper than what is expected from a particulate mixture [32][33][34][35][36][37]. 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 [26,38,39]. * preeti.sahu@ist.ac.at; mmanning@syr.edu 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 particlebased 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 [40]. 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 [41]. 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 [6,19,42].
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 [25,26]. Recent work has shown that geometrical properties of interfacial cells in a confluent neighborhood can be directly affected by increased interfacial contractility and tension [26,43,44]. 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 [45][46][47]. 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 [33,[48][49][50][51][52][53][54][55][56], cells attempt to attain a preferred cell shape index, which is dimensionless ratio of perimeter P and area A i.e. s 2D 0 = P/ √ A in 2D [54,57], and of surface area S and volume V i.e. s 3D 0 = S/V 2/3 in 3D [56]. If the cells are able to attain their preferred cell shape, which generally happens in isotropic tissues for elongated shapes with s 2D 0 > 3.81 and s 3D 0 > 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 [56]. 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 For example, a regular BCC unit cell (truncated octahedron), has a dimensionless shape index of s 0 ∼ 5.31 [58].
Here we set V 0 to 1. Half of the cells are tagged differently, creating a mixture of two cell types-β = 1 or β = 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 HIT [59]. 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: where v i denotes the ith 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 V i 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 α and β respectively. The surface tension σ, for simplicity, is assumed to be the same for all facets. It is non-dimensionalized by K S V 2 0 which is unity for our system. Biological cells can establish heterotypic tension by co-regulating the acto-myosin network and adhesion molecules [20]. A biologically relevant estimate of the heterotypic to homotypic tension ratio, based on examination of the contact angles at cell vertices in ectodermmesoderm co-cultures in Xenopus [25], indicates σ ∼ 2 in natural units of the system. For systems with fluctuations, we analyse the dynamics of over-damped selfpropelled 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 τ 0 s , 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 [26]. 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 [59], we first focus on a fluid-like parameter regime [56] 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 σ = 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 [59], 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, 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.
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 fluidlike 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 [26]. 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 cellscells 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 sur-face area S i . Hence, s boundary is defined as: and s bulk is similarly averaged over the interior cells. In the absence of HIT (σ = 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 (θ), 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 θ. 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 [26]. 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 [26], 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.
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 (Continued from previous page.) (b) Acquired cell shape index (s), plotted with respect to tension (σ), 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 (σ). The control distribution (σ=0) is superimposed on a faint black curve that represents the density for uniform distribution i.e sin θ. (d) The probability distribution P (a) of heterotypic facet area a is shown for increasing tension σ = 0.32, 1.00, 2.00 (magenta to blue). For these data, N=512 and s0 = 5.5.
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- 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 1/3 0 = 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(bc). Moreover, the insets to each figure show that the cell geometries approach the right-prism condition monotonically as a function of increasing HIT.
From work in 2D we understand that perturbations perpendicular to the surface of the interface are very costly [26], 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 [26], focused exclusively on the fluidlike regime with s 2D 0 > 3.81 and s 3D 0 > 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 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 l0 (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 σ = 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 ∆z with as a function of the HIT σ. (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. state for the shape s 0 ∼ 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 [39], 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' [26]. While the former is common in particlebased 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).
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 σ c ∼ 3, both the shape and interfacial energies become comparable (Fig. S5(b)) and at σ c the system discontinuously transitions to the tension-dominated-cuspy branch. The registration value also jumps to R = 1. This branch originates at σ ∼ 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 [26], only the cusp-like states were observed. This leads us to hypothesize that perhaps we observe a transition between normal, parabolic restoring forces to nonanalytic 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 ∼ 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 ∼ 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.
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 valuewhen 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 fluidlike 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. FIG. 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 s0 and interfacial tension σ. 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 [25,38], 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 σ > σ c , the ground state is completely registered, which gives rise to a prismatic geometry. These states are topologically pinned, due to cusplike pinning forces, previously observed in 2D mixtures with HIT. But for tensions σ < σ 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 σ 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 [60][61][62][63][64], 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 technique [6,65]. 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 [66]. 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 [59], 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.

S1. Details about self-propulsion dynamics
For dynamical simulations, we use the static model developed for 3D confluent tissues [56], 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 f i = ∂e/∂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 extracellular fluid. In a self-propelled system, the i th cell has a polarization vectorn i and exerts an active force of v 0 /µ on the static media, where µ is the mobility of the cell and v 0 is the self-propulsion speed. Setting µ = 1, the dynamical equation is : The polarization vector evolves via white Gaussian noise on a unit sphere with diffusion coefficient D r , or where E is the 3 × 3 identity tensor, the dyadic product of the polarization vector with itself is given byn ini and ξ i is the white Gaussian noise with ξ i = 0 and ξ i (t) ξ j (t ) = δ ij δ(t − t )E. Timescales in our dynamic equations are nondimensionalized by the natural time unitt = 1/K V V 4/3 0 , which is unity for our choice of parameter values. We use a selfpropulsion 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/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 100t [67]. 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 τ 0 s ≡ 1000t. Therefore, for simulations in the presence of fluctuations, we use an integration step size ∆t = 0.01t and a time-span T SS of 1 × 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 (σ < 1), achieve a steady state sooner, i.e. in just 1 × 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 f inal ∝ N 1/3 is in fairly good agreement with the observed final demixing values seen in our simulations as shown in Fig S1. 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.

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 [26], 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 , redraw the Voronoi tessellation and calculate final energy.
The initial energy of the 9-cell system is: where γ 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: After displacing the central cell to the right, the new total energy to leading order becomes- 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 x ) that was observed for a perturbation perpendicular to the interface. Moreover, the coefficient is γ 0 2 , which is twice that of coefficient for a perpendicular perturbation [26], 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.
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.

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.

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 fluidlike systems to be at mechanical equilibrium, though this FIG. 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 σ = 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.
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 × 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 σ = 0.16, ∼ 65% of the minimization procedures survive and for σ = 0.32 ∼ 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 .
For the cases that survive, a randomly chosen interfacial cell is perturbed away from being registered by a discrete amount 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 = 10 −6 . Similar plots exist for higher values of . For all cases, we observe a bimodal distribution of forces. The lower peak in the distribu-tion 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 σ = 0.16, only ∼ 15% of the systems are on the right-hand-side Gaussian, while for σ = 0.32, it rises to ∼ 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.