---
abstract: 'Surface tension governed by differential adhesion can drive fluid particle mixtures to segregate into distinct regions, i.e., demix. Does the same phenomenon occur in vertex models of confluent epithelial monolayers? Vertex models are different from particle models in that the interactions between the cells are shape-based, as opposed to metric-based. We investigate whether a disparity in cell shape or size alone is sufficient to drive demixing in bidisperse vertex model fluid mixtures. Surprisingly, we observe that both types of bidisperse systems robustly mix on large lengthscales. On the other hand, shape disparity generates slight demixing over a few cell diameters, i.e. micro-demixing. This result, can be understood by examining the differential energy barriers for neighbor exchanges (T1 transitions). The robustness of mixing at large scales suggests that despite some differences in cell shape and size, progenitor cells can readily mix throughout a developing tissue until acquiring means of recognizing cells of different types.'
author:
- 'Preeti Sahu*$^{1,*}$*, Daniel M. Sussman*$^1$*, M. Cristina Marchetti*$^2$*, M. Lisa Manning*$^1$*, J. M. Schwarz*$^{1,*}$*'
bibliography:
- '2dmixtures.bib'
title: 'Large-scale mixing and small-scale demixing in a confluent model for biological tissues'
---

iquid-liquid phase separation, i.e., demixing, drives patterning. In materials science, demixing between two liquids is typically driven by the energetics of interfacial tension overcoming entropy-driven mixing [@Safran2002]. By changing the rate at which a material is cooled, one can tune between a mixed state at high temperature and a demixed state at low temperature. Depending on the material and the parameter regime, this transition can occur continuously via spinodal decomposition or discontinuously via nucleation [@FrenkelLouis; @CahnHillard]. In order to distinguish between mechanisms it is often useful to analyze the lengthscales of emergent patterns: nucleation and spinodal decomposition give rise to characteristic lengthscales that then coarsen, while in the absence of interfacial tension fluids will mix down to the scale of single molecules. These and related demixing phenomena have been studied numerically using multicomponent Lennard-Jones mixtures in which particles have a fixed shape and an interaction potential that depends on the distance between. The potential also energetically distinguishes between particles of different types to model interfacial tension [@RowlinsonSwinton].

In biology, demixing at the subcellular scale can lead to compartmentalization within cells [@Feric2016], while in a developing organism, demixing can lead to compartmentalization among cells of different type, otherwise know as cell sorting. In fact, interfacial tension-driven demixing has long been invoked to explain cellular patterning. The first among such ideas is the Differential Adhesion Hypothesis (DAH), proposed by Steinberg in 1963 [@STEINBERG1963], to explain patterns in the spatial segregation of progenitor cells, such as ectoderm and mesoderm, during embryonic development. The DAH postulates that tissues behave like immiscible liquids composed of motile cells that rearrange in order to minimize their interfacial tension caused by differences in cell-cell adhesion. Building on the DAH, Harris [@Harris1976a], and later Brodland [@Brodland2002], have highlighted the importance of other contributors to interfacial tension, including regulation of the acto-myosin cortex. There is an emerging consensus [@Yamada2007; @Maitre2012; @Manning2010a; @Amack2012; @Mertz2012] that adhesive molecules help to regulate cortical acto-myosin, which can strongly impact cell sorting. However, it remains controversial whether differential adhesion or differential cortical tension alone is sufficient to generate the level of cell sorting and compartmentalization observed in embryos and cell culture systems [@Heasman1993; @Song2016; @Ninomiya2012; @Pawlizak2015; @Dong2018; @Landsberg2009; @Cochet-Escartin2017]. Several experiments have suggested that additional processes such as specialized cell-cell signaling [@Ninomiya2012] or cellular jamming [@Pawlizak2015] enhance or disrupt sorting in living tissues.

One major difference between immiscible liquids composed of cells and immiscible liquids composed of soft spheres is that in the latter case, the particles have a metric-based interaction, while in epithelial layers and even some 3D tissues, the cells are confluent – they can change their shape to completely fill space—and so their interaction is shape-based. To reflect this property, confluent tissues have been studied theoretically and computationally using vertex or Voronoi models [@Nagai2001; @Farhadifar2007b; @Staple2010a; @Barton2017], where cells are constructed from tessellations of space with no gaps between cells. As active fluctuations drive cellular rearrangements, cells must deform so that no gaps open up between them. This suggests cells are subject to strong geometrical and topological constraints. For example, in flat 2D tilings with three-fold coordinated vertices, the average number of neighbors must be precisely six. This constraint leads one to predict that a rigidity transition should occur when neighbor exchange between six-sided cells cost zero energy, i.e. when cells can form regular pentagons at zero cost [@Bi2015g; @Bi2016a]. This prediction has since been realized in experiments [@Park2015a] and is distinct from rigidity transitions in particulate systems [@SM2018; @MBM2018].

Does such an interaction potential with non-trivial geometrical and topological constraints affect the fundamental definition of surface tension? Work on bidisperse foams modeled as ordered vertex models demonstrate that, in equilibrium, demixed cells of two different areas have a lower energy than a mixed system and so demixed states are energetically preferred [@Graner2002]. However, disperse-in-area foams under large shear strain will mix [@Cox2006]. If we think of the shear strain as a temperature-like variable, then these findings are similar to particulate systems.

On the other hand, recent work by some of us demonstrates that so-called heterotypic contacts in vertex models can drastically affect the notion of interfacial tension [@Sussman2018b]. Heterotypic contacts, where cells recognize neighbors of a different cell type, can be modeled in two-dimensional vertex models with a higher or lower line tension along interfaces between cells of different types, or heterotypic line tension. Such a rule results in very sharp, yet deformable, interfaces [@Sussman2018b] where surface tension measured by macroscopic deformation of an overall droplet shape gives a value in line with equilibrium expectations, yet, surface tension measured from interfacial fluctuations is at least an order of magnitude larger. This discrepancy is due to discontinuous pinning forces generated during topological rearrangements between cells of different types. That is, it is a consequence of the shape-based nature of the interactions.

Here, we explore the possibility of interfacial-tension-driven demixing in a two-dimensional vertex model in the absence of explicit heterotypic line tension. Particulate mixtures can demix when a miscibility parameter, the ratio of the strength of the metric-based interaction between dissimilar particles as compared to similar particles, becomes less than one. Since in vertex models the interaction is shape-based, it is very natural to ask if binary vertex model fluids consisting of mixtures of cells with different preferred cell shapes and/or sizes, accounting for differential adhesion, cortical tension or volume, demix even in the absence of specialized heterotypic interactions. In other words, is there an emergent effective interfacial tension between two cell types that is sufficient enough to sort cells? Should the answer be yes, then one can imagine that the sorting of progenitor cells occurs very early on in the development process before robust heterotypic interfacial tensions are established. Should the answer be no, then cells must establish heterotypic interactions before sorting can occur, suggesting a more important mechanical role for cell recognition receptors than previously thought. The topological nature of the discontinuous pinning forces stabilizing interfaces in vertex model fluid mixtures tells us that once such recognition is in place, a finite amount of activity is required to overcome the discontinuity [@Sussman2018b]. Interestingly, a recent study with both [*in vitro*]{} experiments and cellular Potts model simulations suggests that a large heterotypic line tension is required for cell sorting [@Canty2017b], although the mechanism was left unresolved.

To understand the mechanisms that drive cell sorting in confluent tissues, we study fluid mixtures using a vertex model where the degrees of freedom are the vertices that connect edges between each cell. To avoid complexities introduced by motility [@activematterreview], the dynamics of each vertex is Brownian. We ask whether or not large-scale interfacial tensions and, therefore, cell sorting are emergent/collective properties of such mixtures. We see a small-scale demixing in mixtures with differential adhesion, which is not thermodynamic in origin. We find that this behavior arises from dynamical trapping due to differential T1 barriers.

**Model**\
We simulate multicomponent fluids using a 2D Brownian vertex model, which includes mechanical interactions between cells and dynamics generated by fluctuations at vertices. The mechanical interactions between cells are given by the energy functional for the $j$th cell of type $\beta$: $$\label{eq:energy}
      E_{j,\beta}= K_a (A_{j,\beta} - A_{0,\beta})^2+K_p (P_{j,\beta} -P_{0,\beta})^2,$$ where $A_{j,\beta}$ denotes the $j$th cell area of type $\beta$ and the $j$th cell perimeter of type $\beta$ is denoted by $P_{j,\beta}$. Given the quadratic penalty from deviating for a cell’s preferred area and perimeter, $K_a$ and $K_p$ are area and perimeter stiffnesses, respectively, and both are independent of cell type. Physically, the area term represents the bulk elasticity of the cell, while the perimeter represents the contractility of the acto-myosin cortex with $P_{0,\beta}$ denoting a competition between cortical tension and cell-cell adhesion. The total energy of the tissue is then defined as $E=\sum_{j,\beta} E_{j,\beta}$. An important parameter in these models is the dimensionless shape index $s_{0,\beta} = p_{0,\beta}/ \sqrt{A_{0,\beta}} $. A regular hexagon has a dimensionless shape index of $s_0\approx 3.72$, for example.

To study binary mixtures, we fix $\beta=1,2$ and allow the cell types to have different parameters, $A_{0,\beta}$ and $P_{0,\beta}$ (see Fig. \[fig:MSD\]a). We will focus on cases of 50:50 mixtures where there is an equal number of each cell type, with either preferred shape disparity or preferred area disparity. Unless otherwise specified, the two components are uniformly distributed in the initial state. We set $K_p$ to unity for all systems. To simulate the dynamics of the model, each vertex undergoes over-damped Brownian motion at a fixed temperature of 0.01. One simulation unit time is referred to as $\tau$. See the Methods section for more details.

We are interested in comparing the behavior of these bidisperse systems to ones with an explicit heterotypic line tension (HLT), where cell types $1$ and $2$ recognize their joint interface as a heterotypic interface and, therefore, alter the line tension at that interface. Such interactions are common in cellular Potts models [@Graner1992; @Canty2017b] and have also been studied in vertex and Voronoi models [@Barton2017; @Sussman2018b]. In this case, we add an extra term to the cell energy to arrive at: $$\label{eq:gamma}
        E_{HLT}  =  \sum_{j,\beta} K_a (A_{j,\beta} - A_{0,\beta})^2+K_p (P_{j,\beta} -P_{0,\beta})^2  + \sum_{\langle i,j\rangle} \delta_{\alpha\beta} \gamma l_{ij},$$ The sum is over all edges, $l_{ij}$ ,between cells $i$ and $j$; the delta function is equal to unity if the cell-types $\alpha$ and $\beta$ (of $i^{th}$ and $j^{th}$ cell respectively) are different, and zero if they are the same. For simplicity, we assume that the additional tension, $\gamma$, is the same for all heterotypic edges.

**Results**\
**Stability and fluidity of shape bidisperse mixtures.**\[rheo\] To single out the effect of shape dispersity, we first vary the preferred shapes under the constraint that the preferred/target area is the same across cell types, $A_{0,1}=A_{0,2}=1$.

Previous work on the vertex model has identified a regime in parameter space dominated by a coarsening instability [@Farhadifar2007b; @Staple2010a], where some cells shrink in size and others grow. We expect that heterogeneous $s_0$ values might amplify this instability, as heterogeneity amplifies differences between the cells. To prevent area dispersity from affecting the results in these mixtures, we choose $K_a=100$, which is sufficient to reduce fluctuations in area $A$ from target area $A_0$ to a standard deviation of less than 1 %, preventing the onset of the coarsening instability. Moreover, Fig. \[fig:highKa\] shows that increased area stiffness does not significantly impact the fluidity of homogeneous tissues, as measured by the effective diffusivity (Eqs. \[eq:msdDef\]-\[eq:Def\]), denoted as $D_{eff}$, which is the ratio of the diffusion constant in the presence of interactions to that in the absence of interactions. The onset of a finite effective diffusivity as a function of the shape index remains near $s_0\approx 3.81$ with increasing $K_a$.

![[*Vertex model binary mixtures.*]{} (a) Schematic of vertex-based modeling of a tissue: A typical tessellation with two different types of cells highlighted. The energy depends only on a cell’s perimeter ($P_j$) and area ($A_j$). (b) A heat map of $\log_{10}D_{eff}$ as a function of $\Delta$ on the x-axis and $s_{av}$ on the y-axis. The phase points with: fluid-fluid ($s_{0,1},s_{0,2}>3.81$), solid-fluid ($s_{0,1}<3.81$) and solid-solid ($s_{0,1},s_{0,2}<3.81$) components are denoted by circular, square and star-shaped markers, respectively. Black-filled markers, demarcated by a solid black line, denote mixtures with a $D_{eff}$ less than that of the chosen cutoff of $0.01$. Region above this line denotes fluid-like behaviour on average. []{data-label="fig:MSD"}](PavDeltaplusSchematic.pdf){width="99.00000%"}

Next, we investigate how shape disparity affects the fluidity of the tissue. We find that the most effective way to represent the phase space of the two-component system with two shape indices $s_{0,1},s_{0,2}$ is by the average value of the shape index, $s_{av}= (s_{0,1}+ s_{0,2})/2$, and the shape disparity, which is the difference between the two values, $\Delta=s_{0,2}-s_{0,1}$ with $s_{0,2}>s_{0,1}$. Figure \[fig:MSD\]b is a heat map of the effective diffusivity of binary mixtures as a function of $s_{av}$ and $\Delta$. We see that there is a boundary between fluid-like and solid-like, demarcated by the thick solid line, as determined by the $D_{eff}$ threshold. Interestingly, for $\Delta=0.3$, this boundary does not match up with the fluid-solid boundary for the monodisperse case for the Brownian limit of a self-propelled Voronoi model at a similar temperature, which is near $s_0=3.81$ [@Bi2016a]. Moreover, the solid-fluid mixtures depicted by squares, have a fluid-like diffusivity above the boundary line. This indicates that the fluid-like species in the mixture are sufficient to fluidize the entire tissue, which is additionally confirmed by analyzing the diffusivities of each component (Fig. \[fig:CompDeff\]a).

**[Binary mixtures with two target shapes.]{}** After understanding how $\Delta$ affects diffusivity in a mixture, let us now understand its role in bulk segregation for a fixed $s_{av}=3.85$. A snapshot of a typical long-time configuration for such a mixture is shown in Fig. \[fig:shapeG\]a. By eye, it appears that segregation does occur at very small scales, due to some clustering of the cells with larger $s_0$. The system maintains this small-scale structure at long times. No large-scale segregation is observed. Hence, we shall refer to this process as [*micro-demixing*]{}. To quantify micro-demixing and highlight its long-time steady state, we study three observables.

The first is the demixing parameter $DP$, which measures the average environment of each species, quantifying whether it is more likely to be surrounded by similar (homotypic) or dissimilar (heterotypic) cells. Defining $N_s$ as the number of similar neighboring cells and $N_t$ as the total number of neighboring cells, $$\label{eq:dp}
DP=\big \langle DP_j \big \rangle = \bigg \langle 2\bigg(\frac{N_s}{N_t} - \frac{1}{2} \bigg) \bigg \rangle,$$ where the brackets denote averaging over all cells in the tessellation. In a completely mixed state, $DP=0$, whereas in a completely sorted mixture, $DP=1$, in the limit of infinite system size.

The demixing parameter, shown in Fig. \[fig:shapeG\]b, vanishes at $t=0$, since the two cell types are initially seeded at random, and saturates to a small non-zero value at long times. The final steady state value, $DP_f$, increases with increasing shape disparity $\Delta$, as shown in the inset to Fig. \[fig:shapeG\]b, and the length of time required to reach the steady state also increases with increasing $\Delta$ (Fig. \[fig:CompDeff\](b)). For comparison, the dashed black line in Fig. \[fig:shapeG\]b illustrates the demixing parameter as a function of time for a model with heterotypic line tension. In the HLT case, DP rises very quickly to a value close to unity as one species rapidly forms a circular droplet, in a manner similar to that expected for conventional liquid-liquid binary mixtures.

We then measure the average cluster radius $R$ by quantifying the average radius of gyration of the dispersed component. In the case of shape bidisperse mixtures, the more fluid-like (larger $s_0$) component tends to be dispersed. The average cluster radius (Fig. \[fig:shapeG\]c) shows a small growth in time, which appears to saturate at long times, although the data is noisier given the cluster statistics sampling rate. The steady state radius tends to increase with increasing $\Delta$. In all cases studied, clusters have an average radius of less than 2.5. For comparison, the dashed line shows a system with HLT, which we expect to saturate as a nearly circular droplet of one species embedded inside the other species. For the system size we study, this would correspond to a cluster of radius $~8$, which is close to the observed steady state value of $7.2$.

To further quantify the structure of this micro-demixed state, we study the pair correlation function, $g(r)$, which describes the normalized probability of finding a cell center a given distance from another cell center. In homogeneous fluids and amorphous solids, this function exhibits short range order with peaks occurring at distances that are integer multiples of the typical spacing between two cells. The envelope of these peaks falls off with distance and eventually approaches unity, highlighting that these materials are disordered over larger lengthscales. In bidisperse mixtures, we compute the correlation between each species $\beta$ separately, defined by the relative position vectors ($r^{(\beta)}$) between two cells of type $\beta$: $$\label{eq:gr}
g_{\beta\beta}(r) =\frac{1}{2\pi r N_{\beta} \rho_0} \sum\limits_i^{N_{\beta}}\sum\limits_j^{N_{\beta}}\delta(r -r_{ij}^{(\beta)} ).$$

![[*Shape bidisperse fluid mixtures.*]{} (a) Snapshot of a $P_{av}=3.85$, $\Delta=0.4$, $N=1600$ mixture. Scale bar denotes 10 units. Yellow is used for solid-like cells ($s_0=3.65$) and blue for liquid-like ones ($s_0=4.05$). (b)-(d) Various quantifications of demixing in shape bidisperse mixtures (curves colored from green to red in increasing order of shape disparity i.e $\Delta=0.01,0.1,0.2,0.3,0.4$) are compared to a mixture with an extra heterotypic line tension of value 0.1, $s_0=3.97$ (black dashed curve). (b) Demixing Parameter versus $\log(time)$. The final value ($DP_f$) as a function of $\Delta$ is shown in the inset. (c) Average cluster radius (R) versus time. The inset shows the same plot in log-log. (d) Pair correlation function of high-$s_0$ cells $(g_{22})$ versus radial distance for $t=2\times10^5\,\tau$. The dashed grey line shows an exponential decay. The inset shows the decay lengthscale ($\xi$) in terms of the maximum possible lengthscale ($\xi_{HLT}$) with increasing disparity $\Delta$. Simulation details provided in Table S1. ](shapeGrand.pdf){width="99.00000%"}

[\[fig:shapeG\]]{}

For a completely sorted mixture, $g_{\beta \beta}(r)$ should exhibit an envelope that falls off exponentially, with a length scale $\xi$ that corresponds to the average cluster radius. In the HLT mixtures, where a single droplet forms, we see such a structure, as shown by the dashed black line in Fig. \[fig:shapeG\]d. We extract a length scale of $\xi_{HLT}  \sim 7.5$, which is very similar to the steady state average cluster radius shown in Fig. \[fig:shapeG\]c. For comparison, we measure $\xi$ for all shape bidisperse mixtures and compare it to $\xi_{HLT}$ by computing $\xi/ \xi_{HLT}$ (see inset to Fig. \[fig:shapeG\]d). We find this ratio to be quite small, consistent with previous results.

**Binary mixtures with two target areas.** After studying the impact of shape disparity in cell sorting, we next study the effect of dispersity in area. The mixture is now comprised of equal numbers of cells with $A_{0,1}<A_{0,2}$, where we take $\sqrt{A_{0,1}}$ as the unit of length. Both types have the same $s_0$, or $s_{0,1}=s_{0,2}$. We have taken care to ensure that our area bidisperse mixture are also in a fluid region of the phase diagram (Fig. \[fig:highKa\]b) by checking that $D_{eff}>0.01$. For the results shown here, the shape index is fixed at $s_0=3.85$ to mimic fluid-like cells. We define the ratio of the preferred areas as $\alpha=A_{0,2}/A_{0,1}$.

![ [*Area bidisperse fluid mixtures.*]{} (a) Snapshot of a $s_0=3.85$, $\alpha=2.5$, $N=1600$ mixture. Scale bar denotes 10 units. Yellow is used for larger cells ($A_0=1.43$) and blue for smaller ones ($A_0=0.57$). (b)-(d) Various quantifications of demixing in area bidisperse mixtures (curves colored from green to red in increasing order of size disparity i.e $\alpha=1.0,1.5,2.0,2.5$) are compared to a mixture with an extra heterotypic line tension of value 0.1, $s_0=3.97$ and $A_0=1.0$ (black dashed curve). (b) Demixing Parameter versus $\log(time)$. The final value ($DP_f$) as a function of $\alpha$ is shown in the inset. (c) Average cluster radius (R) versus time. (d) Pair correlation function of small-$A_0$ cells $(g_{ss})$ versus radial distance in units of the smallest lengthscale for $t=2\times10^5\tau$. The dashed grey line shows an exponential decay. The inset shows the decay lengthscale ($\xi$) in terms of the maximum possible lengthscale ($\xi_{HLT}$) with increasing disparity $\Delta$. Simulation details provided in Table S2. []{data-label="fig:areaR"}](areaGrand.pdf){width="99.00000%"}

Visual inspection of a snapshot from a simulation of an area bidisperse mixture with high $\alpha=2.5$ at long times demonstrates that observing cluster formation by eye is difficult, particularly given the disparity in area fraction between the two cell types (see Fig. \[fig:areaR\]). The DP has been measured and is smaller than those found in shape bidisperse mixtures (Fig. \[fig:shapeG\]b). Since the large-$A_0$ cells occupy more than half of the total area, we perform our cluster analysis on cells with $A_{0,1}$. As shown in Fig. \[fig:areaR\]c, the final clusters have an average cluster radius that is typically less than two cell diameters and becomes smaller as $\alpha$ increases. Similarly, Fig. \[fig:areaR\]d illustrates that $g_{ss}(r)$ also shows no sign of bulk segregation, with a structural length scale that is always less than $0.2*\xi_{HLT}$, and decreases with decreasing $\alpha$, as seen in the inset to Fig. \[fig:areaR\]d.

**Zero-temperature energy configurations.** Our finite-temperature simulations suggest that large-scale sorting is not preferred in these mixtures. To understand this, we study an ensemble of energy minimized states. [*If the mixed state has a lower energy at zero temperature, then we expect that energetics cannot drive demixing at finite temperature.*]{} Therefore, we compare the energy of two initial states of $N=400$ cells: a sorted system where all of the cells with cell centers in the left half of the box are labeled type 1, and the remainder are labeled type 2, and a mixed system where cell types are randomly assigned. We use FIRE minimization [@FIRE] to identify the nearest local energy minimum for $250$ realizations in each of the two scenarios.

Figure \[fig:cost\]a shows the ratio between the energy of states with sorted initial conditions ($E_s$) and mixed initial conditions ($E_m$) in the case where type 1 and 2 cells have different shape parameters. At larger system sizes, there is a clear trend that the sorted states typically possess higher energies than the mixed states, so that the ratio rises above unity as the shape disparity $\Delta$ increases. This indicates that there are no energetic forces driving the segregation in larger systems. We have also quantified the effective interfacial line tension (Fig. \[fig:shapeTension\], Fig. \[fig:areaTension\]) using a method developed previously by some of us  [@Yang12663]. We find that there is no emergent line tension in any of these mixtures, which is consistent with our energy calculations. In Fig. \[fig:cost\]b, which shows the ratio of energies between sorted and mixed states for cells with area dispersity, the trend is even clearer. Again, sorted states have significantly higher energy compared to mixed states as the area dispersity $\alpha$ increases.

![[*Minimal energy configurations.*]{} Systems with $N=100,400,900$ cells (green to black) are energy minimized using the FIRE algorithm to get the total energy of the configurations- mixed($E_m$) and sorted($E_s$), for a $<s_{av}>=3.85$, $K_a=100$ with increasing disparity. The ratio $E_s/E_m$ versus $\Delta$. (b) The ratio $E_s/E_m$ is plotted versus $\alpha$ for $s_0=3.85$ and $K_a=1$. Simulation details provided in Table S3.[]{data-label="fig:cost"}](ratio.png){width="99.00000%"}

**Zero-temperature T1 energy barriers.** Although the zero temperature energy calculations above help us understand the lack of macroscopic segregation in mixtures with no heterotypic interfacial tension, they do not explain the small-scale demixing/segregation seen, for example, in Fig. \[fig:shapeG\]b. Since both cell types are subject to the same geometrical and topological constraints and rearrange via T1 transitions, we now turn to an energetic analysis of T1 transitions for the bidisperse system. Specifically, we study the statistics of energy barriers in bidisperse systems, where there are nine types of T1 transitions possible. While we present data in the supplemental information for symmetric cases where two of the cells are of type 1 and two of are type 2, we focus here on asymmetric systems where 3 of the cells are of one type and one is of another type. As illustrated by the 4-cell cluster diagrams in Fig. \[fig:corrSp\], such 3:1 arrangements naturally represent the cost of one cell type invading an interface composed of cells of a distinct type, which determines the dynamic stability of such an interface.

Similarly to previous work [@Bi2015g; @Bi2016a], we compute the T1 energy barrier height by measuring the global tissue energy as we force a single edge in our bidisperse simulation to shrink to zero length while minimizing the energy and allowing the other degrees of freedom to relax, as shown in Fig. \[fig:corrSp\]a. The energy barrier $E_b$ we report in Fig. \[fig:corrSp\]b is the difference between the the final energy $E(l=0)$ at the 4-fold vertex and the initial energy $E_0$, or $E_b=E(l-0)-E(l)$, averaged over $250$ edges with the same topology in small simulated tissues with $N=80$ cells.

![[*Differential energy barriers in shape bidisperse fluid mixtures.*]{} (a) Energy $E(l)$ relative to $E_0$ is plotted against T1 edgelength $l$ for a typical shape bi-disperse T1 pair ($\Delta=0.4, s_{av}=3.85$). (b) Energy Barrier $E_b$ is plotted against signed disparity in shape $\Delta_s$. Positive and negative $\Delta_s$ values imply stiffer cluster in yellow and floppier cluster in blue respectively. Each solid curve represents the barrier for a heterotypic cell to get out of the cluster for $s_{av}= 3.79,3.85,3.88)$ (from solid-like (orange) to liquid-like (green) (c) Correlation plot for $s_{av}=3.85$ between Differential Energy Barriers on the right y-axis $\delta(E_b)$ (in maroon triangles) and demixing relative to mixed scenario $DP_f$ on the left y-axis (in black discs). Shape difference $\Delta$ is plotted on x-axis. Simulation details provided in Table S4.[]{data-label="fig:corrSp"}](corrln_shape_Feb_Fromppt.png){width="99.00000%"}

Figure \[fig:corrSp\]a illustrates a particular type of (3:1) T1 energy profile where a single cell with shape parameter $s_0^{1}$ invades a 3-cell cluster formed by cells with $s_0^{cluster}$. We define a [*s*igned]{} shape disparity $\Delta_s=s_0^{1}-s_0^{cluster}$ to distinguish it from a T1 with cell types swapped. Negative $\Delta_s$ indicates that a more stiff cell is invading a cluster of floppy cells. We have checked that energy barriers are statistically identical for cells entering or exiting a cluster. Because cells are as likely to leave a cluster as to enter it, this suggests that clusters of a given cell type will not grow or shrink over long-time or length-scales.

Figure \[fig:corrSp\]b highlights that the energy barriers associated with these (3:1) transitions systematically increase as the magnitude of the shape dispersity $\Delta_s$ increases. In other words, it becomes energetically more difficult for a single cell to invade or leave a cluster of a different cell type as the shape dispersity between the two types increases. Perhaps more importantly, it also shows that these energy barriers are not symmetric around zero; there is a systematic difference between a stiffer cell invading a floppier cluster and vice-versa, especially for lower values of $s_{av}$ as the system approaches the jamming transition. Stiffer clusters tend to be more difficult to break up than floppier clusters. To characterize this effect, we define the energy barrier disparity between invading stiffer and floppier clusters as $\delta E_b (\Delta)= E_b (\Delta)- E_b (-\Delta)$.

To test whether this mechanism might be relevant for the micro-segregation we observed in our finite-temperature simulations, we directly compare the demixing parameter associated with the final, steady state in each simulation, $DP_f$ to the energy barrier disparity $\delta E_b$ as a function of shape dispersity $\Delta$, as shown in Fig. \[fig:corrSp\]c. This plot shows a quite strong correlation between the two quantities, suggesting that this mechanism is a very likely driver of micro-demixing. To further test this idea, we have increased the temperature for the $\Delta=0.1$ mixtures and found $DP$ to vanish at temperatures higher than the differential energy barrier.

A similar analysis can be performed for area bidisperse mixtures as shown in Fig. \[fig:corr\]. An important difference from the shape bidisperse case is that while there is a clear connection between cell shape and tissue rheology (stiffer cells have a smaller $s_0$), there is no such connection between area and rheology. Moreover, there is very little evidence for micro-segregation, and so we expect the signal to be much weaker. Nevertheless, we can define a quantity $\alpha_s=A_0^{cluster}/A_0^{1}$ that is less than unity if a larger cell is invading a cluster of smaller cells and greater than unity otherwise. Fig. \[fig:corr\]b suggests that large-cell clusters are more difficult to invade than small-cell clusters, although the differential energy barrier is quite a bit smaller than for the case of shape bidispersity. In particular, the case where $s_0=3.85$ is highlighted in Fig. \[fig:corr\]c showing a correlation between demixing and $\delta E_b$, although the amplitude of both effects is quite small.

**Discussion** Using Brownian vertex model simulations, we show that two-dimensional mixtures bidisperse in preferred shape and in preferred area both robustly mix at large scales in the absence of an explicit heterotypic line tension distinguishing between the two cell types. Energy minimization at zero temperature supports our finding for both kinds of dispersity. In particular, mixed systems have lower energy than sorted ones, so that bidispersity is not sufficient to energetically stabilize an interface between the two fluids. For shape bidisperse mixtures, we find that, in spite of having solid-like cells making up half the mixture, the mixtures are still able to fluidize. We also discover micro-demixing in shape bidisperse mixtures, where clustering of the same cell type over several cell diameters is observed.

To understand micro-demixing in shape bidisperse mixtures, we establish a correlation between micro-demixing and zero-temperature differential energy barriers for neighbor exchanges (T1 transitions) between four cells at the heterotypic boundaries. Specifically, we find that the energy barriers for a fluid cell type to invade a cluster of three stiff cells is typically higher than for a stiff cell to invade a cluster of three fluid cells. Hence stiff cells tend to cluster. This stiff cell clustering allows them to percolate throughout the system and thereby shift the location of the rigidity transition for temperatures smaller than the height of the energy barriers. For area bidisperse systems, the differential energy barriers for neighbor exchanges are smaller than for the shape bidisperse case, and we find a rather negligible amount of micro-demixing. Our differential energy barrier calculations at zero temperature also yields a prediction for the temperature above which the micro-demixing does not occur—a prediction that has indeed been verified in our simulations.

The observation of robust mixing on large scales for both types of mixtures is rather surprising, given that the shape-based interaction distinguishes between the two cell types just as changing the strength of the metric-based interaction between two particles of different types in Lennard-Jones mixtures. In the particulate case, there can be large-scale demixing (depending on the miscibility), while in the cellular case, there is not. The robust mixing observed in area bidisperse systems at zero temperature also seems counter-intuitive when compared with area bidisperse foams in ordered hexagonal states [@Graner2002]. In this case, the system demixes at zero temperature given an additional perturbative energetic cost to an interface between cells of slightly different areas. Only for large applied shear strains, do area bidisperse foams mix [@Cox2006]. Understanding differences between foam and vertex models is therefore an interesting area for future study. Foam models lack the $P^2$ contribution to the energy functional (Eq.\[eq:energy\]) and this restricts the fluid-like phase space accessible to such models, perhaps contributing to differences between them.

On the other hand, athermal two-dimensional bidisperse particulate mixtures of different size discs with purely repulsive forces, such as models for granular particles with no (or little) friction, are not expected to sort at small size disparities [@Lopez-Sanchez2013a]. Only as the size dispersity increases does sorting occur due to entropic depletion forces [@Melby2007]. Entropic depletion forces do not apply to a confluent tessellation in which the packing fraction is fixed at unity, though may to some extent apply to Voronoi models. Depletion forces also drive segregation in vertical vibrated shape bidisperse mixtures of rods and spheres [@Rodriguez-Linan2016]. Interestingly, size bidisperse mixtures of active particles, segregate in the absence of any attractive forces [@Yang2014a]. The segregation here is due to an asymmetry in the energy barrier between one smaller particle passing through two larger particles as compared to one larger particle passing through two smaller ones. Given the above analogy, a vertex model fluid mixture perhaps has more in common with an athermal disordered binary packing than with a fluid mixture with differential adhesion. The robust demixing we observe may be more naturally understood if we assume that target shape is controlled by the cell cortical tension, while adhesion simply enforces the confluency condition.

Our results bring an understanding to earlier work demonstrating that sorting at embryonic boundaries requires high heterotypic interfacial tension [@Canty2017b]. Given our T1 energy barrier analysis encoding both the topological and geometrical constraints of confluent packings, we now understand why these mixtures robustly mix. This robustness suggests that despite some difference in shape and size, progenitor cells can readily mix throughout the embyro. To demix (or sort), progenitor cells have learned to work around such constraints by eventually developing means of recognizing whether neighboring cells are of the same type or a different types. And while a small amount of heterotypic line tension can generate stable interfaces [@Sussman2018b] in the absence of fluctuations, fluctuations may be able to overcome such barriers. Our analysis gives a new way to understand bulk behavior based on cellular rearrangements in such confluent mixtures. In other words, based on the analysis of T1 energy barriers between different cell types, experimentalists can predict whether or not different cell types will mix or not mix. T1 transition rates are presumably a reasonable proxy for energy barriers.

Finally, the micro-demixing effect could be utilized in biology to create more subtle patterning. For instance, when randomly tagging a tessellation half with one cell type (and half with another cell type), both cell types percolate individually through the system ([@Bollobas2006]). However, if there is now some spatial correlation in the tagging introduced even at the small scale, such as the tagging of one cell type is positively correlated with tagging a neighboring cell of the same type, then the percolation transition point can be finely tuned to readily transition from a tenuous spanning structure to one that is not tenuous in order to efficiently respond to changes in the environment.

**Methods**\
**Simulation details:** We simulate a vertex model where the degrees of freedom evolve according to over-damped Brownian dynamics [@Sussman2017a]. Other models with directed cell motility are possible, including an Active Vertex Model [@Czajkowski2018] and the SPV model where cells are self-propelled due to an active force with persistence [@Bi2016a; @Barton2017]. Specifically, each vertex $i$ located at coordinate $\bf{r}$ experiences a Brownian force $\mathbf{F}^B$, on each vertex $i$ with $\mathbf{F}_{i}^{B}=$[$\xi$]{}$_i$\[eq:dynbm\] where [$\xi$]{}$_i$ is white noise with zero mean and $\langle\xi_{\gamma i}(t)\xi_{\lambda k}(t')\rangle=2T\delta_{\gamma\lambda}\delta_{ik}\delta(t-t')$, where here $\gamma$ and $\lambda$ denote spatial components. In epithelial layers, we expect that fluctuations are driven by active cytoskeletal components, and hence the $T$ is an effective temperature that represents the magnitude of this activity (Ref. [@Gunst2015]). The equation of motion for a single vertex therefore takes the form $$\dot{\mathbf{r}}_{i}=\mu \mathbf{F}_{i} +\mu \mathbf{F}_{i}^B,\label{eq:dyn}$$ with $\mathbf{F}_{i} =-$[$\nabla$]{}$_i E$, where E is the total energy and $\mathbf{F}_i$ is a non-local effective mechanical force experienced by the $i$th vertex of the $j^{th}$ cell and hence represents the cell-cell interactions. In the absence of mechanical interactions, an isolated cell performs a random walk with a long time effective diffusion rate of $T/\mu$. Unless otherwise specified, $\mu=1$. The over-damped Euler integration method is used to update the equations.

In vertex models, one needs to take care of cellular rearrangements explicitly [@Farhadifar2007b; @Bi2016a]. In the absence of cell division or death, such rearrangements correspond to T1 transitions in which one edge shrinks to zero length and two new cells are connected via a new growing edge. In simulations, if an edge length falls below some threshold length $l_c$, then we rotate the edge by $\pi/2$ and reconnect the topology of the surrounding cells to generate a local neighbor-exchange. Such neighbor-exchange moves are only accepted if they reduce the global energy. Unless otherwise specified $l_c$ is set to 0.04. The noise is controlled by temperature (T) which is set to 0.01.

Past work has demonstrated that the mechanical properties of vertex models depend sensitively on the shape parameter $s_0$ and temperature $T$. Specifically, these models exhibit rigidity [@Bi2015g] or glass [@Sussman2018c] transitions where the system transitions from more solid-like to more fluid-like. At $T=0$, the 2D vertex model exhibits a rigidity transition as a function of cell shape parameterized by $s_0$. Above a critical value of target shape index $s_0^* \sim 3.81$, cells are able to move past each other with very small energy cost and below which they cannot. To understand this transition, one analyzes the energetics of how cells move past each other via T1 transitions. A minimal four cell calculation with fixed unit area hexagonal cells revealed that if the two cells that would no longer share an edge after the edge swap formed regular pentagons, then the energy barrier for the formation of a four-vertex vanishes, suggesting that pentagon shape formation is a geometrically compatible transition pathway for three-fold coordinated lattices. Interestingly, the shape parameter for a regular pentagon is $s_0=2 \sqrt{5} (5-2\sqrt{5})^{1/4}\approx 3.812\approx s_0^*$. In the presence of activity or temperature, vestiges of this zero-temperature rigidity transition have been found in a glassy transition between fluid-like and more solid-like behavior in an active Self-Propelled Voronoi (SPV) model [@Bi2016a] and a Brownian Voronoi model [@Sussman2018c].

Given the complex phase behavior of such vertex models, we want to ensure that the mixtures are fluid-like. To do so, we first measure the Mean-Squared Displacement (MSD). To account for global tissue motion possible in these types of models [@Giavazzi2018], we define the displacement of each cell in a time window $t$, $\mathbf{x}(t)$, as the distance the cell traveled in time $t$ minus the total displacement of the entire system of cells over that same time interval. Then the MSD is defined as $$\label{eq:msdDef}
MSD(\Delta t)  \equiv  \langle (\mathbf{x}(t + \Delta t)-\mathbf{x}(t))^2 \rangle,$$ where $\langle \cdot \rangle$ denotes an average over all cells in the tissue and all times $t$. The self-diffusivity $D_{s}$, is defined by assuming the long-time behavior of the system is diffusive, $$\label{eq:msdDs}
 D_s=\lim_{t\rightarrow\infty} \frac{MSD(t)}{4t}.$$ To understand whether cells are being constrained by their neighbors, we compare $D_{s}$ to the *bare*, or non-interacting, diffusion constant $D_0$. For a non-interacting Brownian particle at temperature T with mobility $\mu$, the Fluctuation-Dissipation theorem states that $D_0=\mu k_B T$, where $k_B$ is Boltzmann constant. We set $\mu k_B$ to unity. The effective diffusivity is given by $$\label{eq:Def}
D_{eff}= \frac{D_s}{D_0 }.$$ Systems with small $D_{eff}$ are more solid-like, while systems with large $D_{eff}$ are more fluid-like. In practice, we use a threshold of $0.01$ to distinguish between these different behaviors, in line with previous work [@Bi2016a].\
**Code availability:** The codes are programmed using open source cellGPU code available at <https://gitlab.com/dmsussman/cellGPU>.\
**Data availability:** The data that support the findings of this study are available from the corresponding authors upon reasonable request.\

**Acknowledgments** We would like to thank Matthias Merkel for useful discussions. This work was supported by NSF-POLS-1607416 (MCM, MLM, JMS) and Simons Foundation-454947 (MLM) and Simons Foundation-342345 (MCM). All authors also acknowledge support of the Syracuse University Soft and Living Matter Program.\
**Author contributions**: All authors conceived the study; P.S. conducted the simulations; all authors interpreted the results and wrote the paper.\
**Additional information**:\
**Competing interests:** The authors declare no competing financial or non-financial interests.

Supplemental Material {#supplemental-material .unnumbered}
=====================

Simulation Parameters {#simulation-parameters .unnumbered}
---------------------

Here we provide tables for the parameters used for each aspect of the simulations. For our dynamical simulations, the systems are equilibrated for time, $t_{eqb}$, and subsequently run for a longer time. For our FIRE simulations, simulations typically run until the maximum force experienced by a vertex reduces below a threshold value of $10^{-13}$.

  Parameters                                               values
  ---------------------------------------------- ---------------------------
  1\. Ensembles                                              50
  2\. $s_{av}$                                           3.79- 3.91
  3\. $\Delta$                                             0.0-0.4
  4\. $t_{eqb}$                                          1000$\tau$
  5\. $dt=\tau$                                             0.001
  6\. $(K_a,K_p)$                                          (100,1)
  7\. $T$                                                   0.01
  8\. $N$                                                    400
  9\. Total time                                  $2\times\,10^5 + t_{eqb}$
  10\. $l_c$                                                0.04
  11\. HLT $(\gamma)\,\,\, for\,\,\, p_0=3.97$               0.1

  : Area Bi-disperse Dynamical Simulations

\[tab:shape\_dynamical\]

  Parameters                             values
  ---------------------------- ---------------------------
  1\. Ensembles                            50
  2\. $s_{0}$                             3.85
  3\. $\alpha$                           1.0-2.5
  4\. $t_{eqb}$                        1000$\tau$
  5\. $dt$                                0.01
  6\. $(K_a,K_p)$                         (1,1)
  7\. $T$                                 0.01
  8\. $N$                                  400
  9\. Total time                $2\times\,10^5 + t_{eqb}$
  10\. $l_c$                              0.04
  11\. $\langle A_0 \rangle$                1

  : Area Bi-disperse Dynamical Simulations

\[tab:area\_dynamical\]

  Parameters                           values
  ------------------------ ------------------------------
  1\. Ensembles                         250
  2\. $s_{av},s_0$                      3.85
  3\. $\alpha$                        1.0-2.5
  4\. $\Delta$                         0-0.12
  5\. $dt$                              0.01
  6\. $K_a$                 1$(\alpha)$ & 100 $(\Delta)$
  7\. $T$                               0.01
  8\. $N$                           100,400,900
  9\. Maximum FIRE steps              $10^5 $
  10\. $l_c$                            0.04

  : T1 energy barriers

\[tab:InterfaceCost\]

  Parameters                           values
  ------------------------ ------------------------------
  1\. Ensembles                         250
  2\. $s_{av}$                       3.79-3.88
  3\. $\alpha$                        1.0-2.5
  4\. $s_0$                          3.82-3.88
  5\. $dt$                              0.01
  6\. $K_a$                 1$(\alpha)$ & 100 $(\Delta)$
  7\. $T$                               0.01
  8\. $N$                                80
  9\. Maximum FIRE steps              $10^5 $
  10\. $l_c$                            0.04
  11\. $\Delta$                        0-0.12

  : T1 energy barriers

\[tab:T1FIRE\]

Effect of area stiffness on fluidity  {#effect-of-area-stiffness-on-fluidity .unnumbered}
-------------------------------------

High shape-disparity can amplify coarsening in mixtures, resulting in further enhanced disparity in cell areas. To prevent this coarsening from occuring, we increase the area stiffness $K_a$ to 100. To make sure this does not affect the fluid phase seen in monodisperse mixtures, we study the effective diffusivity as a function of the target shape parameter for several $K_a$ values. We find that $K_a$ barely affects the diffisivities and that the large changes in curvature of $D_{eff}$ versus $s_0$ remain close to $3.81$ such that larger values of $K_a$ do not significantly affect the fluidity of the cells. See Fig. \[fig:highKa\].

![(a)[*Effective diffusivity ($D_{eff}$) with respect to target shape parameter $s_0$.*]{} Different curves are for monodisperse systems with $K_a$ varying from 1 to 100. The solid horizontal line represents the cutoff at 0.01 used previously. The vertical dashed line denotes $s_0*\approx 3.81$. (b)[*Effective diffusivity in area bidisperse mixtures.*]{} Plot of the effective diffusivity ($D_{eff}$) with respect to increasing area dispersity $\alpha$. Parameter details provided in Table S2.](Deffs_monoKa_area_v2.png){width="\textwidth"}

[ \[fig:highKa\]]{}

Diffusivity of area bidisperse mixtures {#diffusivity-of-area-bidisperse-mixtures .unnumbered}
---------------------------------------

Monodisperse systems with $s_0>3.81$ have a fluid-like diffusivity. Here we check diffusivity for mixtures having the same $s_0=3.85$ for all cell types but bidisperse in size. We see that the average fluid-like diffusivity remains unchanged. See Fig. \[fig:highKa\]b.

Component-wise diffusivity and timescales in shape bidisperse mixtures {#component-wise-diffusivity-and-timescales-in-shape-bidisperse-mixtures .unnumbered}
----------------------------------------------------------------------

We study the diffusivities of individual components for mixtures with fixed $s_{av}=3.85$. Although increased dispersity signals a solid-fluid mixture, we see that the average behavior remains fluid-like up to high dispersities. Hence, we measure the diffusivity of each component to determine if the solid-like cells diffuse (Fig. \[fig:CompDeff\]a). We find that a fluid-like component is indeed able to help the solid-like cells diffuse.

For the demixing observed in shape bidisperse mixtures, as mentioned in the main text, we observe that for most of the $\Delta$, the DP saturates to a final value. We check if the timescale associated with this saturation increases with dispersity since Fig. \[fig:CompDeff\]a demonstrates that the solid components (of high dispersity mixtures) do not diffuse as much. We define $t_{1/2}$ as the average time taken by the system to get to half of its final DP. We observe that this half-time increases exponentially with $\Delta$, as shown in the inset of Fig. \[fig:CompDeff\]b.

![[*Component-wise diffusivities and timescale to approach steady state.*]{} (a) Plot for effective diffusivity ($D_{eff}$) with respect to increasing shape dispersity $\Delta$. The solid lines are for the two different components with triangles and circles representing higher $s_0$ (type 2) and lower $s_0$ (type 1) respectively. The dashed curve represents the averaged $D_{eff}$. (b) The average time it takes for the system to achieve half of its steady state DP value, or $t_{1/2}$, is plotted against $\Delta$. The solid curves from \[fig:shapeG\](b) are used to compute $t_{1/2}$. The inset shows log-log plot of the same, with a linear fit in solid yellow line is $y=3x+3$.](CompDeff_Halftime.png){width="\textwidth"}

[ \[fig:CompDeff\]]{}

Cortical tension for sorted vs. mixed configurations {#cortical-tension-for-sorted-vs.-mixed-configurations .unnumbered}
----------------------------------------------------

An emergent line tension between two different kinds of cells must show a high line tension along heterotypic edges and lower line tension along the homotypic edges. Hence, for both the sorted and mixed scenarios (Fig. \[fig:cost\]), we study a *line tension map* where the thickness of the edge is linearly proportional to its line tension. A positive value is colored in red and a negative value is colored in blue. The cortical tension for each edge can be computed using the method suggested in Ref. [@Yang12663].

![[*Cortical tension in shape bidisperse mixtures of $\Delta=0.4$.*]{} Left and right panels shows line tension maps for sorted and mixed scenarios for a $N=900$ system respectively. Heterotypic edges are shown in dash-dot lines. Yellow and blue cells have $p_0=3.65$ and $4.05$ respectively. They are followed by histograms for heterotypic(in black) and homotypic(colored) edges. Vertical lines show the mean values for each curve in their respective colors.](shape_Tension_v2.png){width="99.00000%"}

[ \[fig:shapeTension\]]{}

![[*Cortical tension in area bidisperse mixtures of $\alpha=2.5$.*]{} Left and right panels shows line tension maps for sorted and mixed scenarios for a $N=900$ system respectively. Heterotypic edges are shown in dash-dot lines. Yellow and blue cells have $A_0=0.57$ and $1.43$ respectively. They are followed by histograms for heterotypic(in black) and homotypic(colored) edges. Vertical lines show the mean values for each curve in their respective colors.](area_Tension_v2.png){width="99.00000%"}

[ \[fig:areaTension\]]{}

The cortical tension analysis conveys the fact that there is no emergent line tension due to bidispersity in the mixtures we study. The mean heterotypic line tension (black vertical line) is less than or equal to the mean homotypic line tension (colored vertical line) for all the scenarios. See Figs. \[fig:shapeTension\] and  \[fig:areaTension\].

Differential T1 energy barriers in area bidisperse mixtures {#differential-t1-energy-barriers-in-area-bidisperse-mixtures .unnumbered}
-----------------------------------------------------------

We present data supporting the notion that the differential energy barriers are smaller for the area bidisperse mixtures as compared to the shape bidisperse ones. We focus on larger cells trying to invade a cluster of smaller cells and vice versa to determine the stability (or lack thereof) an interface. See Fig. \[fig:corr\]. We also present data for other types of topologies for both shape and area bidisperse mixtures for completeness (see Fig. \[fig:otherT1s\]).

![[*Differential energy barriers in area bidisperse mixtures.*]{} (a) Energy $E(l)$ relative to $E_0$ versus T1 edgelength $l$ for a typical size bi-disperse T1 pair ($\alpha=2.5, s_{av}=3.85$). (b) Energy Barrier $E_b$ is plotted against area disparity $\alpha_s$ where large values on right and small values on left imply large-cell cluster in yellow and small-cell cluster in blue respectively. Each solid curve represents the barrier for a heterotypic cell to get out of the cluster for a fixed $s_{0}$ (varied from solid-like (orange) to liquid-like (green) - 3.82,3.85,3.88) (c) Correlation plot for $s_{0}=3.85$ between Differential Energy Barriers on the right y-axis $\delta(E_b)$ (in maroon triangles) and demixing relative to mixed scenario $DP_f$ on the left y-axis (in black discs). Size ratio $\alpha$ is plotted on x-axis. Simulation details provided in Table \[tab:T1FIRE\]. []{data-label="fig:corr"}](corrln_area_Feb.pdf){width="99.00000%"}

Finally, to study differential T1 energy barriers in a simplified setting, we consider four cells connected to each other symmetrically. The energy is minimized with respect to a diminishing T1 edge length $l$ using MATLAB. The area stiffness is kept very high and the initial condition is recursively fed from a longer $l$ energy minimized configuration to the subsequent shorter $l$. We can accommodate different sizes and shapes as long as cells of different types are positioned symmetrically about both $x$ and $y$ axis and make the cells sharing the T1 edge (T1 pair) have different properties from the non-T1 pair. The formula used to compute energy barrier is $E(l=l_{H})-E(l=0)$, where $l_{H}$ is the edge length of a uniform hexagon with unit area.

![[*Symmetric 4-cell T1 energy barriers for shape bidispersity.*]{} On the left is the color plot of energy barrier as a function of independently tunable shapes of T1 pair and non-T1 pair is plotted along x-axis and y-axis respectively. The dashed line represents monodisperse calculation ie for $\Delta=0$. As expected it is red till it reaches the monodisperse transition point $s*_0=3.813$, after which it becomes blue. Off-diagonal phase points depict bidisperse mixtures ie $\Delta\neq0$. We see that it is necessary for the T1-pair to be fluid like, for vanishing barrier. On the right is a cross-section of the phase diagram on left. Energy barrier is plotted against area disparity for increasing values of $s_0$ 3.79 to 3.85.](matlab_shape.png){width="80.00000%"}

[ \[fig:4cellShape\]]{}

![[*Symmetric 4-cell T1 energy barriers for area bidispersity.*]{} On the left is the color plot of energy barrier as a function of independently tunable sizes of T1 pair (blue polygons) and non-T1 pair(yellow polygons). On the left is when blue polygons are bigger than yellow. On the right is smaller blue cells sandwiched between yellow (BssB). $\alpha$ and $P_0$ are the area ratios and preferred shape index respectively. Dashed black line represents the monodisperse transition point $s*_0=3.813$. This graph predicts the energy barriers to vanish at a shape index higher than $s*_0$ in highly bidisperse systems. On the right is a cross-section of the cumulative phase diagram on left. Energy barrier is plotted against area disparity for increasing values of $s_0$ 3.79 to 3.85.](matlab_area.png){width="99.00000%"}

[ \[fig:4cellarea\]]{}

To study the effect of shape bidispersity (Fig. \[fig:4cellShape\]), the energy barrier (red when non-zero and blue when vanishingly small) is plotted with respect to the shape of T1 pair (x-axis) and the shape of the non-T1 pair (y-axis), which can be independently varied. A similar analysis is done for mixtures with bidisperse areas (Fig. \[fig:4cellarea\]). We observe differential energy barriers in both cases with, again, the size of the barrier generically larger in the shape bidisperse case as compared to the area bidisperse case, even in this simplified calculation. Re-phrasing this in terms of invading a cluster of the opposite type, one can think of these as invading doublets of opposite kind.

![[*T1 transitions in shape and area bidisperse mixtures.*]{} (a)-(c) T1 topologies (shown as cartoons on axis extremities) and their barrier statistics for shape bidisperse mixtures. (d)-(f) T1 topologies (shown as cartoons on axis extremities) and their barrier statistics for size bidisperse mixtures. Parameters used are in Table \[tab:T1FIRE\]. ](RestT1s_v2.png){width="\textwidth"}

[ \[fig:otherT1s\]]{}
