# Large-scale mixing and small-scale demixing in a confluent model for biological tissues

###### 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.

L 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 Safran (2002). 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 Frenkel and Louis (1992); Cahn and Hillard (1958). 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 Rowlinson and Swinton (1982).

In biology, demixing at the subcellular scale can lead to compartmentalization within cells Feric et al. (2016), 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 Steinberg (1963), 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 Harris (1976), and later Brodland Brodland (2002), have highlighted the importance of other contributors to interfacial tension, including regulation of the acto-myosin cortex. There is an emerging consensus Yamada and Nelson (2007); Maître et al. (2012); Manning et al. (2010); Amack and Manning (2012); Mertz et al. (2012) 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 Heasman et al. (1993); Song et al. (2016); Ninomiya et al. (2012); Pawlizak et al. (2015); Dong and Klumpp (2018); Landsberg et al. (2009); Cochet-Escartin et al. (2017). Several experiments have suggested that additional processes such as specialized cell-cell signaling Ninomiya et al. (2012) or cellular jamming Pawlizak et al. (2015) 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 Nagai and Honda (2001); Farhadifar et al. (2007); Staple et al. (2010); Barton et al. (2017), 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 Bi et al. (2015, 2016). This prediction has since been realized in experiments Park and et al. (2015) and is distinct from rigidity transitions in particulate systems Sussman and Merkel (2018); Moshe et al. (2018).

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 Graner (2002). However, disperse-in-area foams under large shear strain will mix Cox (2006). 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 Sussman et al. (2018a). 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 Sussman et al. (2018a) 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 Sussman et al. (2018a). 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 Canty et al. (2017), 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 Marchetti and et al. (2013), 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 th cell of type :

(1) |

where denotes the th cell area of type and the th cell perimeter of type is denoted by . Given the quadratic penalty from deviating for a cell’s preferred area and perimeter, and 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 denoting a competition between cortical tension and cell-cell adhesion. The total energy of the tissue is then defined as . An important parameter in these models is the dimensionless shape index . A regular hexagon has a dimensionless shape index of , for example.

To study binary mixtures, we fix and allow the cell types to have different parameters, and (see Fig. 1a). 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 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 . 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 and 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 Graner and Glazier (1992); Canty et al. (2017) and have also been studied in vertex and Voronoi models Barton et al. (2017); Sussman et al. (2018a). In this case, we add an extra term to the cell energy to arrive at:

(2) |

The sum is over all edges, ,between cells and ; the delta function is equal to unity if the cell-types and (of and cell respectively) are different, and zero if they are the same. For simplicity, we assume that the additional tension, , is the same for all heterotypic edges.

Results

Stability and fluidity of shape bidisperse mixtures.
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, .

Previous work on the vertex model has identified a regime in parameter space dominated by a coarsening instability Farhadifar et al. (2007); Staple et al. (2010), where some cells shrink in size and others grow. We expect that heterogeneous 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 , which is sufficient to reduce fluctuations in area from target area to a standard deviation of less than 1 %, preventing the onset of the coarsening instability. Moreover, Fig. S1 shows that increased area stiffness does not significantly impact the fluidity of homogeneous tissues, as measured by the effective diffusivity (Eqs. 6-8), denoted as , 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 with increasing .

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 is by the average value of the shape index, , and the shape disparity, which is the difference between the two values, with . Figure 1b is a heat map of the effective diffusivity of binary mixtures as a function of and . We see that there is a boundary between fluid-like and solid-like, demarcated by the thick solid line, as determined by the threshold. Interestingly, for , 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 Bi et al. (2016). 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. S2a).

Binary mixtures with two target shapes. After understanding how affects diffusivity in a mixture, let us now understand its role in bulk segregation for a fixed . A snapshot of a typical long-time configuration for such a mixture is shown in Fig. 2a. By eye, it appears that segregation does occur at very small scales, due to some clustering of the cells with larger . 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 , 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 as the number of similar neighboring cells and as the total number of neighboring cells,

(3) |

where the brackets denote averaging over all cells in the tessellation. In a completely mixed state, , whereas in a completely sorted mixture, , in the limit of infinite system size.

The demixing parameter, shown in Fig. 2b, vanishes at , 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, , increases with increasing shape disparity , as shown in the inset to Fig. 2b, and the length of time required to reach the steady state also increases with increasing (Fig. S2(b)). For comparison, the dashed black line in Fig. 2b 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 by quantifying the average radius of gyration of the dispersed component. In the case of shape bidisperse mixtures, the more fluid-like (larger ) component tends to be dispersed. The average cluster radius (Fig. 2c) 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 . 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 , which is close to the observed steady state value of .

To further quantify the structure of this micro-demixed state, we study the pair correlation function, , 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 separately, defined by the relative position vectors () between two cells of type :

(4) |

For a completely sorted mixture, should exhibit an envelope that falls off exponentially, with a length scale 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. 2d. We extract a length scale of , which is very similar to the steady state average cluster radius shown in Fig. 2c. For comparison, we measure for all shape bidisperse mixtures and compare it to by computing (see inset to Fig. 2d). 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 , where we take as the unit of length. Both types have the same , or . We have taken care to ensure that our area bidisperse mixture are also in a fluid region of the phase diagram (Fig. S1b) by checking that . For the results shown here, the shape index is fixed at to mimic fluid-like cells. We define the ratio of the preferred areas as .

Visual inspection of a snapshot from a simulation of an area bidisperse mixture with high 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. 3). The DP has been measured and is smaller than those found in shape bidisperse mixtures (Fig. 2b). Since the large- cells occupy more than half of the total area, we perform our cluster analysis on cells with . As shown in Fig. 3c, the final clusters have an average cluster radius that is typically less than two cell diameters and becomes smaller as increases. Similarly, Fig. 3d illustrates that also shows no sign of bulk segregation, with a structural length scale that is always less than , and decreases with decreasing , as seen in the inset to Fig. 3d.

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 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 Bitzek et al. (2006) to identify the nearest local energy minimum for realizations in each of the two scenarios.

Figure 4a shows the ratio between the energy of states with sorted initial conditions () and mixed initial conditions () 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 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. S3, Fig. S4) using a method developed previously by some of us Yang et al. (2017). We find that there is no emergent line tension in any of these mixtures, which is consistent with our energy calculations. In Fig. 4b, 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 increases.

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. 2b. 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. 5, 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 Bi et al. (2015, 2016), 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. 5a. The energy barrier we report in Fig. 5b is the difference between the the final energy at the 4-fold vertex and the initial energy , or , averaged over edges with the same topology in small simulated tissues with cells.

Figure 5a illustrates a particular type of (3:1) T1 energy profile where a single cell with shape parameter invades a 3-cell cluster formed by cells with . We define a signed shape disparity to distinguish it from a T1 with cell types swapped. Negative 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 5b highlights that the energy barriers associated with these (3:1) transitions systematically increase as the magnitude of the shape dispersity 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 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 .

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, to the energy barrier disparity as a function of shape dispersity , as shown in Fig. 5c. 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 mixtures and found to vanish at temperatures higher than the differential energy barrier.

A similar analysis can be performed for area bidisperse mixtures as shown in Fig. S5. 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 ), 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 that is less than unity if a larger cell is invading a cluster of smaller cells and greater than unity otherwise. Fig. S5b 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 is highlighted in Fig. S5c showing a correlation between demixing and , 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 Graner (2002). 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 Cox (2006). Understanding differences between foam and vertex models is therefore an interesting area for future study. Foam models lack the contribution to the energy functional (Eq.1) 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 López-Sánchez et al. (2013). Only as the size dispersity increases does sorting occur due to entropic depletion forces Melby et al. (2007). 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 Rodríguez-Liñán et al. (2016). Interestingly, size bidisperse mixtures of active particles, segregate in the absence of any attractive forces Yang et al. (2014). 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 Canty et al. (2017). 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 Sussman et al. (2018a) 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 (Bollobás and Riordan (2006)). 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 Sussman (2017). Other models with directed cell motility are possible, including an Active Vertex Model Czajkowski et al. (2018) and the SPV model where cells are self-propelled due to an active force with persistence Bi et al. (2016); Barton et al. (2017). Specifically, each vertex located at coordinate experiences a Brownian force , on each vertex with where is white noise with zero mean and , where here and denote spatial components. In epithelial layers, we expect that fluctuations are driven by active cytoskeletal components, and hence the is an effective temperature that represents the magnitude of this activity (Ref. Gunst and Fredberg (2003)). The equation of motion for a single vertex therefore takes the form

(5) |

with , where E is the total energy and is a non-local effective mechanical force experienced by the th vertex of the 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 . Unless otherwise specified, . The over-damped Euler integration method is used to update the equations.

In vertex models, one needs to take care of cellular rearrangements explicitly Farhadifar et al. (2007); Bi et al. (2016). 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 , then we rotate the edge by 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 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 and temperature . Specifically, these models exhibit rigidity Bi et al. (2015) or glass Sussman et al. (2018b) transitions where the system transitions from more solid-like to more fluid-like. At , the 2D vertex model exhibits a rigidity transition as a function of cell shape parameterized by . Above a critical value of target shape index , 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 . 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 Bi et al. (2016) and a Brownian Voronoi model Sussman et al. (2018b).

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 Giavazzi et al. (2018), we define the displacement of each cell in a time window , , as the distance the cell traveled in time minus the total displacement of the entire system of cells over that same time interval. Then the MSD is defined as

(6) |

where denotes an average over all cells in the tissue and all times . The self-diffusivity , is defined by assuming the long-time behavior of the system is diffusive,

(7) |

To understand whether cells are being constrained by their neighbors, we compare to the bare, or non-interacting, diffusion constant . For a non-interacting Brownian particle at temperature T with mobility , the Fluctuation-Dissipation theorem states that , where is Boltzmann constant. We set to unity. The effective diffusivity is given by

(8) |

Systems with small are more solid-like, while systems with large are more fluid-like. In practice, we use a threshold of to distinguish between these different behaviors, in line with previous work Bi et al. (2016).

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.

## References

- Safran (2002) S. A. Safran, “Statistical thermodynamics of soft surfaces,” Surf. Sci. 500, 127–146 (2002).
- Frenkel and Louis (1992) D. Frenkel and A. A. Louis, “Phase separation in binary hard-core mixtures: An exact result,” Phys. Rev. Lett. 68, 3363 (1992).
- Cahn and Hillard (1958) J. W. Cahn and J. E. Hillard, “Free energy of a nonuniform system. i. interfacial free energy,” J. Chem. Phys. 28, 258 (1958).
- Rowlinson and Swinton (1982) J. S. Rowlinson and F. L. Swinton, Liquids and Liquids Mixtures, 3rd ed. (Elsevier, 1982).
- Feric et al. (2016) M. Feric, N. Vaidya, T. S. Harmon, D. M. Mitrea, L. Zhu, T. M. Richardson, R. W. Kriwacki, R. V. Pappu, and C. P. Brangwynne, “Coexisting Liquid Phases Underlie Nucleolar Subcompartments,” Cell 165, 1686–1697 (2016).
- Steinberg (1963) M. S. Steinberg, “Reconstruction of tissues by dissociated cells. Some morphogenetic tissue movements and the sorting out of embryonic cells may have a common explanation.” Science 141, 401–408 (1963).
- Harris (1976) Albert K. Harris, “Is cell sorting caused by differences in the work of intercellular adhesion? a critique of the steinberg hypothesis,” J. Theor. Biol. 61, 267–285 (1976).
- Brodland (2002) G. W. Brodland, “The Differential Interfacial Tension Hypothesis (DITH): A Comprehensive Theory for the Self-Rearrangement of Embryonic Cells and Tissues,” AMSE J. Biomech. Eng. 124, 189–197 (2002).
- Yamada and Nelson (2007) Soichiro Yamada and W. James Nelson, “Localized zones of Rho and Rac activities drive initiation and expansion of epithelial cell-cell adhesion,” J. Cell Biol. 178, 517 (2007).
- Maître et al. (2012) Jean Léon Maître, Hélène Berthoumieux, Simon Frederik Gabriel Krens, Guillaume Salbreux, Frank Jülicher, Ewa Paluch, and Carl Philipp Heisenberg, “Adhesion functions in cell sorting by mechanically coupling the cortices of adhering cells,” Science 338, 253–256 (2012).
- Manning et al. (2010) M. Lisa Manning, Ramsey A. Foty, Malcolm S. Steinberg, and Eva-Maria Schoetz, “Coaction of intercellular adhesion and cortical tension specifies tissue surface tension,” Proc. Nat, Acad. Sci. USA 107, 12517–12522 (2010).
- Amack and Manning (2012) J. D. Amack and M. L. Manning, “Knowing the boundaries: Extending the differential adhesion hypothesis in embryonic cell sorting,” Science 338, 212 (2012).
- Mertz et al. (2012) Aaron F. Mertz, Shiladitya Banerjee, Yonglu Che, Guy K. German, Ye Xu, Callen Hyland, M. Cristina Marchetti, Valerie Horsley, and Eric R. Dufresne, “Scaling of traction forces with the size of cohesive cell colonies,” Phys. Rev. Lett. 108, 198101 (2012).
- Heasman et al. (1993) J Heasman, D Ginsberg, B Geiger, K Goldstone, T Pratt, C Yoshida-Noro, and C Wylie, “A functional test for maternally inherited cadherin in Xenopus shows its importance in cell adhesion at the blastula stage,” Develop. 120, 49 (1993).
- Song et al. (2016) Wei Song, Chih Kuan Tung, Yen Chun Lu, Yehudah Pardo, Mingming Wu, Moumita Das, Der I. Kao, Shuibing Chen, and Minglin Ma, “Dynamic self-organization of microwell-aggregated cellular mixtures,” Soft Matt. 12, 5739 (2016).
- Ninomiya et al. (2012) H. Ninomiya, R. David, E. W. Damm, F. Fagotto, C. M. Niessen, and R. Winklbauer, “Cadherin-dependent differential cell adhesion in Xenopus causes cell sorting in vitro but not in the embryo,” J. Cell Sci. 125, 1877 (2012).
- Pawlizak et al. (2015) Steve Pawlizak, Anatol W. Fritsch, Steffen Grosser, Dave Ahrens, Tobias Thalheim, Stefanie Riedel, Tobias R. Kießling, Linda Oswald, Mareike Zink, M. Lisa Manning, and Josef A. Käs, ‘‘Testing the differential adhesion hypothesis across the epithelial-mesenchymal transition,” New J. Phys. 17, 091001 (2015).
- Dong and Klumpp (2018) J. J. Dong and S. Klumpp, “Simulation of colony pattern formation under differential adhesion and cell proliferation,” Soft Matter 14, 1908 (2018).
- Landsberg et al. (2009) Katharina P. Landsberg, Reza Farhadifar, Jonas Ranft, Daiki Umetsu, Thomas J. Widmann, Thomas Bittig, Amani Said, Frank Jülicher, and Christian Dahmann, “Increased Cell Bond Tension Governs Cell Sorting at the Drosophila Anteroposterior Compartment Boundary,” Curr. Biol. 19, 1950 (2009).
- Cochet-Escartin et al. (2017) Olivier Cochet-Escartin, Tiffany T. Locke, Winnie H. Shi, Robert E. Steele, and Eva Maria S. Collins, “Physical Mechanisms Driving Cell Sorting in Hydra,” Biophys. J. 113, 2827 (2017).
- Nagai and Honda (2001) Tatsuzo Nagai and Hisao Honda, “A dynamic cell model for the formation of epithelial tissues,” Phil. Mag. B: Phys. Cond. Matt., Stat. Mech., Elec. Opt. Mag. Prop. 81, 699–719 (2001).
- Farhadifar et al. (2007) R. Farhadifar, B. Röper, J. C.and Aigouy, S. Eaton, and F. Jülicher, “The Influence of Cell Mechanics, Cell-Cell Interactions, and Proliferation on Epithelial Packing,” Curr. Biol. 24, 2095–2104 (2007).
- Staple et al. (2010) D. B. Staple, R. Farhadifar, J. C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, “Mechanics and remodelling of cell packings in epithelia,” Eur. Phys. J. E 33, 117–127 (2010).
- Barton et al. (2017) D. L. Barton, S. Henkes, C. J. Weijer, and Sknepnek R., “Active vertex model for cell-resolution description of epithelial tissue mechanics,” PLoS Comput. Biol. 13, e1005569 (2017).
- Bi et al. (2015) D. Bi, J. H. Lopez, J. M. Schwarz, and M. L. Manning, “A density-independent rigidity transition in biological tissues,” Nat. Phys. 11, 1074–1079 (2015).
- Bi et al. (2016) D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, “Motility-driven glass and jamming transitions in biological tissues,” Phys. Rev. X 6, 021011 (2016).
- Park and et al. (2015) Jin Ah Park and et al., “Unjamming and cell shape in the asthmatic airway epithelium,” Nat. Mat. 14, 1040–1048 (2015).
- Sussman and Merkel (2018) D. M. Sussman and M. Merkel, “No unjamming transition in a voronoi model of biological tissue,” Soft Matt. 14, 3397–3403 (2018).
- Moshe et al. (2018) M. Moshe, M. Bowick, and M. C. Marchett, “Geometric frustration and solid-solid transitions in model 2d tissue,” Phys. Rev. Lett. 120, 268105 (2018).
- Graner (2002) Francois Graner, “Two-Dimensional Fluid Foams at Equilibrium,” in Morphology of Condensed Matter - Physics and Geometry of Spatially Complex Systems (Springer-Verlag, Berlin, 2002).
- Cox (2006) S. J. Cox, “The mixing of bubbles in two-dimensional bidisperse foams under extensional shear,” J. Non-Newton. Fluid Mech. 137, 39–45 (2006).
- Sussman et al. (2018a) Daniel M. Sussman, J. M. Schwarz, M. Cristina Marchetti, and M. Lisa Manning, “Soft yet Sharp Interfaces in a Vertex Model of Confluent Tissue,” Phys. Rev. Lett. 120, 058001 (2018a).
- Canty et al. (2017) L. Canty, E. Zarour, L. Kashkooli, P. François, and F. Fagotto, “Sorting at embryonic boundaries requires high heterotypic interfacial tension,” Nat. Comm. 8, 157 (2017).
- Marchetti and et al. (2013) M. C. Marchetti and et al., “Hydrodynamics of soft active matter,” Rev. Mod. Phys. 85, 1143–1189 (2013).
- Graner and Glazier (1992) François Graner and James A. Glazier, “Simulation of biological cell sorting using a two-dimensional extended Potts model,” Phys. Rev. Lett. 69, 213–216 (1992).
- Bitzek et al. (2006) E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, “Structural relaxation made simple,” Phys. Rev. Lett. 97, 170201 (2006).
- Yang et al. (2017) Xingbo Yang, Dapeng Bi, Michael Czajkowski, Matthias Merkel, M. Lisa Manning, and M. Cristina Marchetti, “Correlating cell shape and cellular stress in motile confluent tissues,” Proc. Nat. Acad. Sci. USA 114, 12663–12668 (2017).
- López-Sánchez et al. (2013) Erik López-Sánchez, César D. Estrada-Álvarez, Gabriel Pérez-Ángel, José Miguel Méndez-Alcaraz, Pedro González-Mozuelos, and Ramón Castañeda-Priego, “Demixing transition, structure, and depletion forces in binary mixtures of hard-spheres: The role of bridge functions,” J. Chem. Phys. 139, 104908 (2013).
- Melby et al. (2007) Paul Melby, Alexis Prevost, David A. Egolf, and Jeffrey S. Urbach, “Depletion force in a bidisperse granular layer,” Phys. Rev. E 76, 051307 (2007).
- Rodríguez-Liñán et al. (2016) Gustavo M. Rodríguez-Liñán, Yuri Nahmad-Molinari, and Gabriel Pérez-Ángel, “Clustering-induced attraction in granular mixtures of rods and spheres,” PLoS ONE 11, e0156153 (2016).
- Yang et al. (2014) Xingbo Yang, M. Lisa Manning, and M. Cristina Marchetti, “Aggregation and segregation of confined active particles,” Soft Matt. 10, 6477–6484 (2014).
- Bollobás and Riordan (2006) Béla Bollobás and Oliver Riordan, “The critical probability for random voronoi percolation in the plane is 1/2,” Prob. Th. Rel. Fields 136, 417–468 (2006).
- Sussman (2017) D. M. Sussman, “cellGPU: Massively parallel simulations of dynamic vertex models,” Comp. Phys. Commun. 219, 400–406 (2017).
- Czajkowski et al. (2018) M. Czajkowski, D. Bi, M. L. Manning, and M. C. Marchetti, “Hydrodynamics of shape-driven rigidity transitions in motile tissues,” Soft Matt. 14, 5628–5642 (2018).
- Gunst and Fredberg (2003) Susan J. Gunst and Jeffrey J. Fredberg, “The first three minutes: smooth muscle contraction, cytoskeletal events, and soft glasses,” J. Appl. Physiol. 95, 413–425 (2003).
- Sussman et al. (2018b) Daniel M. Sussman, M. Paoluzzi, M. Cristina Marchetti, and M. Lisa Manning, “Anomalous glassy dynamics in simple models of dense biological tissue,” Europhys. Lett. 121, 36001 (2018b).
- Giavazzi et al. (2018) F. Giavazzi, C. Malinverno, G. Scita, and R. Cerbino, ‘‘Tracking-Free Determination of Single-Cell Displacements and Division Rates in Confluent Monolayers,” Frontiers in Physics 6, 120 (2018).

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

### Simulation Parameters

Here we provide tables for the parameters used for each aspect of the simulations. For our dynamical simulations, the systems are equilibrated for time, , 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 .

Parameters | values |
---|---|

1. Ensembles | 50 |

2. | 3.79- 3.91 |

3. | 0.0-0.4 |

4. | 1000 |

5. | 0.001 |

6. | (100,1) |

7. | 0.01 |

8. | 400 |

9. Total time | |

10. | 0.04 |

11. HLT | 0.1 |

Parameters | values |
---|---|

1. Ensembles | 50 |

2. | 3.85 |

3. | 1.0-2.5 |

4. | 1000 |

5. | 0.01 |

6. | (1,1) |

7. | 0.01 |

8. | 400 |

9. Total time | |

10. | 0.04 |

11. | 1 |

Parameters | values |
---|---|

1. Ensembles | 250 |

2. | 3.85 |

3. | 1.0-2.5 |

4. | 0-0.12 |

5. | 0.01 |

6. | 1 & 100 |

7. | 0.01 |

8. | 100,400,900 |

9. Maximum FIRE steps | |

10. | 0.04 |

Parameters | values |
---|---|

1. Ensembles | 250 |

2. | 3.79-3.88 |

3. | 1.0-2.5 |

4. | 3.82-3.88 |

5. | 0.01 |

6. | 1 & 100 |

7. | 0.01 |

8. | 80 |

9. Maximum FIRE steps | |

10. | 0.04 |

11. | 0-0.12 |

### Effect of area stiffness on fluidity

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 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 values. We find that barely affects the diffisivities and that the large changes in curvature of versus remain close to such that larger values of do not significantly affect the fluidity of the cells. See Fig. S1.

### Diffusivity of area bidisperse mixtures

Monodisperse systems with have a fluid-like diffusivity. Here we check diffusivity for mixtures having the same for all cell types but bidisperse in size. We see that the average fluid-like diffusivity remains unchanged. See Fig. S1b.

### Component-wise diffusivity and timescales in shape bidisperse mixtures

We study the diffusivities of individual components for mixtures with fixed . 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. S2a). 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 , the DP saturates to a final value. We check if the timescale associated with this saturation increases with dispersity since Fig. S2a demonstrates that the solid components (of high dispersity mixtures) do not diffuse as much. We define 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 , as shown in the inset of Fig. S2b.

### Cortical tension for sorted vs. mixed configurations

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. 4), 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. Yang et al. (2017).

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. S3 and S4.

### Differential T1 energy barriers in area bidisperse mixtures

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. S5. We also present data for other types of topologies for both shape and area bidisperse mixtures for completeness (see Fig. S8).

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 using MATLAB. The area stiffness is kept very high and the initial condition is recursively fed from a longer energy minimized configuration to the subsequent shorter . We can accommodate different sizes and shapes as long as cells of different types are positioned symmetrically about both and 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 , where is the edge length of a uniform hexagon with unit area.

To study the effect of shape bidispersity (Fig. S6), 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. S7). 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.