Anomalous cell sorting behavior in mixed monolayers discloses hidden system complexities

In tissue development, wound healing and aberrant cancer progression cell–cell interactions drive mixing and segregation of cellular composites. However, the exact nature of these interactions is unsettled. Here we study the dynamics of packed, heterogeneous cellular systems using wound closure experiments. In contrast to previous cell sorting experiments, we find non-universal sorting behavior. For example, monolayer tissue composites with two distinct cell types that show low and high neighbor exchange rates (i.e., MCF-10A & MDA-MB-231) produce segregated domains of each cell type, contrary to conventional expectation that the construct should stay jammed in its initial configuration. On the other hand, tissue compounds where both cell types exhibit high neighbor exchange rates (i.e., MDA-MB-231 & MDA-MB-436) produce highly mixed arrangements despite their differences in intercellular adhesion strength. The anomalies allude to a complex multi-parameter space underlying these sorting dynamics, which remains elusive in simpler systems and theories merely focusing on bulk properties. Using cell tracking data, velocity profiles, neighborhood volatility, and computational modeling, we classify asymmetric interfacial dynamics. We indicate certain understudied facets, such as the effects of cell death & division, mechanical hindrance, active nematic behavior, and laminar & turbulent flow as their potential drivers. Our findings suggest that further analysis and an update of theoretical models, to capture the diverse range of active boundary dynamics which potentially influence self-organization, is warranted.


Introduction
In systems comprising of two or more distinct components, the extent of mixing or demixing of these components can give rise to unique patterns and structures. Such mixing behavior has been studied in detail in solids, fluids, and gases and some generalizable rules can explain observed mixing behaviors [1][2][3]. However, predicting if two populations will mix or demix with limited prior knowledge of these populations is problematic. Additionally, rules that apply to passive systems of particles might not be suitable for active systems with autonomous force generation that can stabilize non-equilibrium configurations [4,5]. Another complexity arises when you have multi-component systems with packing fractions close to unity [6]. In such systems, long-range force transmission and short-range caging effects add further complexity that needs to be resolved to predict mixing and demixing phenomena. Biological tissues are the perfect example of such active and densely packed systems. Mixing, and sorting behavior in these systems is critical not only during embryogenesis, but also in processes involving interactions of distinct cell types such as wound healing, tissue regeneration, and cancer growth and metastasis [7].
Dynamics of tissues have been extensively studied over the past century and several models have been proposed to explain various dynamic features of heterogeneous cell populations. Foremost among them is the differential adhesion hypothesis (DAH) [8]. This theory describes tissue segregation as a phenomenon based on disparities in adhesion interactions of the involved cell populations analogous to the immiscibility of liquids due to differing surface tensions. Beyond differences in adhesion, the DAH framework was expanded to include differences in interfacial tension due to cortical tension, growth (homeostatic) pressure, and preferred cell shape. Overall, this group of tissue surface tension (TST) based models predicts demixing at the tissue level when populations with different interfacial tensions exist [9]. Our past studies on cell sorting in compact 3D multicellular spheroids showed demixing behavior, which partially conformed to these models [10]. However, these TST models are not generalizable to all cell systems. For example, prior experiments also revealed that demixed populations of cells with similar interfacial tensions can remain demixed [11]. Likewise, this study demonstrates that cells with different interfacial tensions can still mix and do not inevitably form independent domains in 2D confluent cell layers.
Apart from the TST models, there is an entire class of models that center on the structural stability of cell populations based on inter-cellular neighborhood dynamics. These self-propelled Voronoi (SPV) models focus on high-density systems (packing fraction one). The motility of the particles and population rearrangements occur via an edge-to-vertex topology transformation, known as a T1 transition. SPV frameworks uniquely model a transition between jammed and unjammed states, a feature that seems to approximate tissue dynamics of multiple confluent cellular systems [12][13][14]. However, while they accurately capture tissue dynamics for jammed systems, which prevent cell reorganization, unjammed systems seem to generally follow the outcomes predicted by the DAH, when active rearrangements are possible and allow for liquid-like cell dynamics [10,15].
In this work, we find that in 2D, tissue-mimetic, confluent co-culture systems, the mixing/demixing behavior of mechanically different cell types follows an anomalous development that does not conform to these existing hypotheses. We analyze the cell migration and rearrangement dynamics to tease out the principles governing the overserved mixing behavior. Our observations suggest additional characteristics of dynamic tissue systems concerning neighborhood-dependent changes in individual and collective cell mechanics at interfaces have to be recognized to fully explain the mixing/demixing behavior of these 2D systems.

Results
We studied the long-term dynamics of co-cultures of epithelial MCF-10A and mesenchymal MDA-MB-231 cells, mesenchymal MDA-MB-231 and myoepithelial MDA-MB-436 cells, and two differently stained populations of mesenchymal MDA-MB-231 cells [16]. By themselves, the epithelial MCF-10A cells are stiffer and form jammed, rigid monolayers with minimal neighbor exchanges [12,17].To facilitate discussion, we have referred to these cells as solid-like (SL) cell collectives. On the other hand, the two mesenchymal cell lines by themselves are softer and form unjammed cell monolayers with significant rearrangements between neighbors, and are therefore referred to as fluid-like (FL) cells [10].

Domain formation dynamics are governed by cell type combinations
Our central discovery is the disparity of domain structuring in dependence on the system composition; the grouping of SL-FL composites (i.e., MCF-10A & MDA-MB-231) resulted in segregated domains of SL cells, which underwent very few neighborhood exchanges, and zones of FL cells, which readily swapped positions with their neighbors. In FL-FL mixtures (i.e., MDA-MB-231 & MDA-MB-436 or MDA-MB-231 & MDA-MB-231) this separation did not develop and both cell types behaved FL and freely exchanged neighbors and ultimately stayed in a mixed system configuration. Our observations are based on two different setups, which employed a wound closure assay as described in the methods section. The first setup was comprised of initially separated populations of SL & FL or FL & FL cells on either side of a wound enclosure. We observed that as the wound closes and the two populations met, the SL-FL system formed a distinct seam that prevented the FL cells from penetrating the SL cell layer ( figure 1(A)). Staining of the actin cytoskeleton (supplementary figure S1 (https://stacks.iop.org/NJP/23/043034/mmedia)) revealed that this long-lasting ridge is mainly composed of a multitude of interconnected actin fibers, similar to the well-known supracellular actomyosin cables found in other wound closing processes [18]. On the other hand, for both the FL-FL populations, as the wound closed, no distinct seam was formed. Instead, the cell populations mixed into each other in a diffusive fashion, seemingly only marginally hindered by contact with their neighboring cells (figures 1(B) and (C)).  The second setup consisted of pre-mixed populations of SL-FL or FL-FL cells on both sides of the wound assay, as described in the methods section 4.4. As the wound closed and homeostasis was achieved, we found that the SL-FL populations segregated from each other into clusters (groups of cells larger than 60 μm in diameter) (figure 2(A)). These clusters were characterized by asymmetric, elongated shapes (median aspect ratio above 1.8, supplementary figure S2), and did not particularly migrate or coalesce with each other within our experimental observation window, which generally lasted more than 48 h. Groups of SL cells once again formed a prominent network of cortical actin around them (supplementary figure S1). In contrast, the FL-FL pre-mixed populations did not identifiably segregate from each other or formed domains of individual cell types throughout the observation. To quantify the extent of mixing or de-mixing, we measured the segregation parameter defined by the number of unlike neighbors over the total number of neighbors, normalized with the global population ratio. This is plotted for all systems in figure 3. We found that in all initial mixed-population systems, there was some amount of segregation over time. However, the extent of segregation depended on the system makeup: FL-FL mixtures with 231 cells were the most mixed at the end of the experiment with a parameter of 0.78 ± 0.02 (mean ± SD), FL-FL mixtures with 436-231 cells were somewhat mixed with 0.65 ± 0.04, and SL-FL systems with 10A-231 cells were the least mixed with a considerably lower segregation parameter of 0.42 ± 0.10. Both the rate of segregation and the absolute difference between initial and final segregation state were vastly more pronounced in SL-FL systems. These composites showed the highest amount of demixing, whereas both FL-FL systems exhibited both less overall segregation and much less progress in demixing over time and remained fairly mixed throughout. Naturally, the 231-231 mixture was the least segregated composite, and it mainly serves as a control system. The main difference between the two cell types is merely the stable expression of GFP-labeled cytosol in one species.

Tracking individual cells highlights distinct motility phenotypes of mixed and segregated solid-like and fluid-like cells
To further understand the factors driving the segregation and mixing behavior in these two populations we tracked individual cells within these systems as explained in the methods section 4.4. This allowed us to quantify the velocity of individuals cells as well as their neighbor exchanges. We evaluated the system behavior based on the tracking of more than 175 000 cells across all system configurations to depict the ensemble values shown in all figures. Based on our analysis, we found that the SL cells migrated nearly two to three times slower than the FL cells in all systems at homeostasis (before the wound was opened), as seen in figures 4(A) and (B). When the wound was opened, the SL-FL system exhibited a characteristic jump in velocity. The SL cells sped up and migrated nearly as fast as the 231 FL cells. Interestingly, the speed of migration of FL cells was system-specific and dependent on the other cell type they were seeded with ( figure 4(B)). Both FL cell populations in the 231-436 cell mixtures migrated more than 1.5 times faster than the FL cells in the 231-231 mixtures and the 10A-231 mixture. Also, the opening of the wound had no observable effect on the migration velocities of cells in both FL-FL systems, as can be witnessed in the rather constant velocity profile throughout the experiment. While the SL cells migrated slower than FL cells, they were not stationary as would be naively expected for these epithelial cells that normally form rigid, jammed monolayers [12,17].
However, neither cell type migrated altogether randomly within their domains. There was a pronounced spatial correlation (figure 4(C)) between the velocities of SL-FL cells and their neighbors across multiple cell lengths, suggesting that they migrate collectively with their surrounding cells over long distances. The temporal correlation (figure 4(D)) of cells in this composite was much more pronounced across all investigated time frames. On the other hand, the correlation between FL cells and their neighbors in both FL-FL systems was much weaker. It dissipated over short distances (single cell neighbors) and receded for periods longer than four hours, indicating that they generally followed their own random paths. Using the individual cell positions as a basis, we extracted Voronoi shapes of all cells and were able to analyze T1 transitions dynamics [13]. The unique role of SL-FL composites was further consolidated by the observation of fewer neighbor exchanges of SL cells ( figure 5(A)). Though, once again, in contrast to expectations, the SL domains were not perfectly jammed and still exhibited neighbor exchanges, as previously reported [19]. However, compared to the SL domains, there were 1.5 to 2 times more FL cell neighbor exchanges (figures 5(A) and (B)). All cells in FL-FL systems exhibited a relatively steady high neighborhood volatility, which nearly coincided for both cell types in the composite (apparent when comparing graphs of same system type in figures 5(A) and (B)). This consistency did not extend to the SL-FL system, where 231 cells transitioned more frequently than 10A cells and also showed an increase during wound healing, similar to the one observed in the velocity (figures 4(A) and (B)). Interestingly, the interface between SL and FL cells (figure 5(C)) was significantly more active in terms of neighbor exchanges, while both FL-FL composites were slightly more engaged in the bulk of the layer (supplemental table S1). This distinctive behavior at the interface also emerged, when we decomposed the mean velocity into parallel and perpendicular components. In both FL-FL systems (figures 5(D) and (E)) perpendicular and parallel components nearly coincided in their mean velocity. MCF-10A cells of the SL-FL composite, on the other hand, showed a small but significant mismatch, as mean velocity components of perpendicular velocities were slightly slower compared to parallel components. Both cell types of the SL-FL system showed a significant predisposition towards motion parallel to the interface: upwards of 9% of the mean velocity was directed parallel versus perpendicular to the interface, more than three times as much as all cell types of FL-FL composites (supplemental table S1).
A comparison of cage-relative mean-square displacement (MSD) [20] (figures 6(A) and (B)) revealed that SL cells (green tracks in figure 6(A)) were also up to an order of magnitude less motile than FL cells. This illustrates the fairly static nature of 10A neighborhoods in the SL-FL composite and their reliance on collective motion, while also verifying that no saturation of MSD and thus no arrest took place. When we compared the cage-relative MSD with the regular MSD (figure 6(A)-inset)) the almost exclusive collective motion of SL cells became more apparent. All FL cell types across all systems (figure 6(B)) did now show particularly neighborhood-limited MSD behavior. This disparate nature is further exemplified in the long-term motile behavior of the cells, seen via cell traces in figures 6(C) and (D). The neighbor-coordinated movement of SL cells culminated in temporally and spatially stable and therefore easily distinguishable domains (figure 6(C)) in SL-FL systems. Tracks of 231 cells (red) were both longer and much less directed, which resulted in substantial random overlap within their domains. Both cell types in FL-FL mixtures (figure 6(D)) showed similar random and interweaved motion as 231 cells in SL-FL systems, which did not result in any pronounced spatial or temporal domain formation.
We quantified this irregularity of motion by comparing the local (100 μm neighborhood) turbulence intensities across all systems. The SL-FL systems were the least turbulent with an average of 1.38 ± 0.33 (mean ± SD), while both FL-FL configurations had an average turbulence of 1.69 ± 0.36 and 1.91 ± 0.41 for 231-231 and 231-436 respectively.

Discussion
Previously we have shown that mixtures of cell populations universally segregate in 3D spheroids [10]. This observation is partially supported by the DAH, which suggests that cells with distinct intercellular adhesion strengths should self-segregate due to dissimilar interfacial tensions between dissimilar cell types. Though, in spheroids, the DAH had to be extended to incorporate the apparent co-regulation of adhesion and cortical tension to explain the sorting hierarchy [10]. However, this same theory does not cover the mixing behavior of 2D systems of dissimilar FL (MDA-MB-436 and MDA-MB-231) cell types presented here. Modeling of various particle systems has revealed many pathways for phase separation, which seem to be caused by a generic equilibrium-breaking mechanism via local energy input [21,22]. However, our FL-FL systems do not adhere to these pathways, even though they can also be regarded in similar theoretical frameworks. Although 436 and 231 cells have very different intercellular adhesion strengths [10], they freely mix with each other in a confluent 2D culture and do not segregate (figures 7(D)-(F)). However, both the demixing of an initially mixed and the lack of mixing in an originally segregated SL-FL system fit the DAH model.
Alternatively, the SPV models that focus on cell motility can be consulted to explain the mixing of FL-FL cell systems observed here. However, these models cannot explain the segregation of initially mixed SL-FL systems into distinct domains (figures 7(A)-(C)), as they predict that at greater 20% population fraction, the SL cells generate impermeable networks that can rigidify and jam the entire system [23].
Thus, the mixing and demixing behaviors observed here in 2D co-cultures of different cell lines present an anomaly that does not fit existing model frameworks. Instead, based on the results described above, we propose four possible mechanisms that future models need to consider: steric mechanical hindrance, FL cooperative dynamics, active nematic behavior, and layer deformations due to cell proliferation and apoptosis. The steric mechanical hindrance occurring at the interface of epithelial cells. A source for mechanical barriers could be the formation of a continuous supra-cellular actin cortex that envelopes large domains of epithelial cells ( figure 7(B)). These actin structures are observed surrounding the MCF-10A domains shown here (supplemental figure S1), as well as reported in several other studies [18,24]. It appears that these actin structures while presenting only a marginal hindrance to individual FL cells (supplemental figure S3), could prevent large-scale mixing of the SL and FL domains. Additionally, purely physical mechanisms also seem to play a role in sustaining tissue boundaries through complex interactions of traction forces and intercellular stresses in jammed cells. These effects occur irrespective of biochemical interactions at the boundary, even epithelial cells of the same type have shown such repulsive interactions. They manifest in long-lived deformation waves propagating some distance away from the boundary, due to cooperative mechanics of jammed cells on account of deformation-and force-differentials [25]. The generic nature of these propagation waves and cooperative mechanics of cell packs alludes to this being another demixing promoting factor, which coincides with the collective behavior and streamlined motility of our SL epithelial cells. The exact mechanics governing these actin networks spanning multiple cells need to be further evaluated to understand their implications in the overall dynamics of heterogenous tissues of various nature.
Motility phenotypes of SL and FL cells might drive mixing and demixing phenomena akin to those observed in fluid dynamics. The SL cells, while motile, migrate with a high degree of spatial and temporal correlation, giving rise to laminar flow-like streamlines (figure 6(C)). On the other hand, the random motion of individual FL cells mirrors the intertwined and unsteady paths of turbulent flow ( figure 6(D)). These differences in flow profile may lead to demixing between turbulent and laminar flow regions (figure 7(C)), while it could facilitate mixing between regions of turbulent flow (figure 7(F)) [26,27], thus offering one interpretation of the observed mixing and demixing behavior. In fact, when comparing collective migration in cells of different invasive potential a key role has been attributed to cooperative migratory dynamics akin to fluid flow [17]. The influence of cooperative cell pack dynamics thus seems to be a non-negligible factor in the fluidization of densely-packed systems. Once again, the fluid dynamic nature of heterogeneous tissues and its influence on system dynamics needs to be further investigated.
A novel study recently revealed another possible sorting mechanism through coordinative inter-and intracellular stresses, which result in collective extensile and contractile active nematic behavior. These differences are sufficient to prompt demixing in systems with mixed extensile and contractile cells and their colony formation dynamics suggest cell-substrate adhesion overcomes line tension-based effects [30]. The decomposition of interfacial motility revealed that SL cells are preferentially moving parallel to the cluster interface (figure 5(F)), while FL cells do not possess this bias (figures 5(D) and (E)). This collective orientation behavior also manifested in the streamlined long-term cell tracks of SL cells (figure 6(C)). The propensity of 10A cells not to move perpendicular to the interface and invade 231 clusters might be indicative of an underlying line-tension. Since cell polarity and traction forces are coupled to cell motility [31] the interfacial orientation of SL cells could indicate polar nematic properties [32]. Such active nematic behavior in velocity-space due to flow-polarity coupling could signify extensile activity of SL cells, which might lead to stress orientation along the green flow profiles in figure 7(C). MDA-MB-231 cells on the other hand have been shown to increase their contractility to resist forces from steric hindrance [33]. These differences in extensile and contractile cell activity might explain the cluster sorting observed in SL-FL composites. Substrate adhesion and coordinative stresses, therefore, seem to be another piece in the puzzle of understanding active sorting behavior.
The presence of other forces such as those arising from cell death or division may result in rearrangements that transcend those predicted by the DAH and SPV models. Even the SL epithelial domains in our systems are not perfectly jammed and undergo migration and neighbor exchanges, which cannot be inferred from these theories. Cell deaths and divisions, which may not be dismissed in our setup, can generate the forces required to unjam epithelial monolayers of SL cells [28]. Additionally, differential rates of death and division between the cell types, and neighborhood-dependent cell death and division rates for individual cells ( figure 7(A)), as shown in other works, can introduce dynamics that may explain the anomalous mixing and demixing behaviors seen here. They can even lead to the forfeiture of pure phases in the DAH, by preventing full isolation and pure intermixing configurations [29]. To test this hypothesis, we applied cell shape-dependent death and division rates in a standard SPV model. While the SPV model without this enhancement did not explain the observed segregation dynamics (supplemental figure S6(A) and (B)), the inclusion of cell deaths and divisions altered the predicted tissue dynamics to match the experimentally observed segregation behavior in SL-FL and FL-FL co-culture systems (supplemental figure S6(C) and (D)). This suggests a need for a deeper dive into the role of cell proliferation and apoptosis rates on the overall dynamics of heterogeneous 2D tissue systems.
Classical interpretations and modern enhancements of cell segregation utilize imbalanced forces, due to differences in single-cell properties between in-domain and surface/interface configurations, to introduce tension-like bulk effects that drive progressional demixing. As such they regard cell composites as Newtonian fluids governed by bulk parameters and neglect some active dynamics. We have reported mixing and demixing behaviors in co-cultures of different cell lines that do not fit the expected dynamics of confluent tissues. In SL-FL systems domain formation into clusters of predominantly FL and SL cells takes place (figures 7(A)-(C)), while both cell types in FL-FL systems (figures 7(D)-(F)) remain rather mixed. Based on individual cell tracking, cell migration profiles, and neighbor exchanges, we propose possible driving factors that have previously been neglected in established models. Our assays highlight the existence of boundary effects, which manifest through asymmetric neighbor exchanges (figures 5(D)-(F) and figure 7(B)). We suggest steric mechanical hindrance, laminar-and turbulent-like cell motility, active nematic behavior due to coordinative stresses, and cell proliferation and apoptosis as some possible drivers of asymmetric interfacial dynamics (compiled in the illustrations of figure 7).
The compiled interfacial effects illustrate the need to expand prevailing theories beyond the mere incorporation of bulk parameters to holistically model cell sorting. It is conceivable that some of the active dynamics described here generate different activation energies for T1 transitions between bulk and interface domains. These differences could cause emergent collective boundary effects that are capable of self-stabilizing the interface. So, while established sorting theories are correct in focusing on interfaces, they are probably neglecting important factors: by selectively fixating on bulk parameters, they are likely overlooking an underlying multi-parameter phase space, which is externalized in the diverse range of active interfacial dynamics displayed by these composite cell systems. The results presented here and their intrinsic mechanics, though not obvious, are extremely important in understanding critical tissue-level processes such as embryogenesis, wound healing, tissue regeneration, and cancer progression and metastasis.
All cell lines were incubated at 37 • C in 95% air/5% CO 2 atmosphere. The culture medium was changed every 2-3 days and cells were passaged every 4-5 days. For cell detachment a PBS solution containing 0.025% (w/v) trypsin and 0.011% (w/v) EDTA (Cat.No. L11-004, PAA) was applied for several minutes. Cells were only used for experiments at a low passage number.

Fluorescence staining
All cells were DNA labeled with SiR-DNA to allow for long-term tracing of cell nuclei even in confluent culture. SiR-DNA was added in 100 nM concentration to the medium to maintain a stable nuclear signal over longer periods of time. The heterogeneous characteristics of the system were investigated by employing permanently transfected MDA-MB-231 cells with a GFP-labeled cytosol (denoted by MDA-MB-231g when distinction was necessary).

Time lapse microscopy
For phase contrast microscopy migratory measurements, cells were washed with phosphate buffer saline, trypsinized, centrifuged and resuspended in the culture medium. Cells were counted in an EVE automatic cell counter (NanoEntek Inc.) and seeded in ibidi two-well culture-inserts (München, Germany; product number: 80209) on ibidi μ-plate 24-well ibiTreat (München, Germany; product number: 82406) tissue culture treated polymer coverslip wells. The cells were allowed to attach, divide and form a confluent monolayer by incubation at 37 • C with 5% CO 2 , after which the insert was removed. At this point medium was switched to a low-glucose version (1400 mg L −1 glucose) 50%/50% mixture of DMEM/H-F12 supplemented with 10% fetal calf serum (PAA A15-151) and 1% 10 000 U ml −1 penicillin/streptomycin (Sigma P0781) to maintain both cell lines at a reduced proliferation rate. The multi-well was placed on an inverted spinning disc confocal microscope (Axio Observer Z1, Carl Zeiss Microscopy GmbH Jena, Germany) surrounded by a commercial incubation chamber with feedback-controlled temperature set to 37 • C and 5% CO 2 air saturation. Frames were captured for each well every 10 min (using an Orca Flash 4.0 camera (Hamamatsu Photonics), 10× objective with NA 0.3 and a 0.5× C-mount).

Cell tracking and analysis in segregation experiments
Automated nucleus tracking via TrackMate was employed to track cells over time [34] and these traces allowed for a determination of a variety of layer parameters such as trajectories, velocity, persistence, cluster development, and transitions via customized MATLAB routines. In co-cultures cells were segmented by thresholding the GFP-cytosol labeled cell signal and binary masking of nuclei in these areas.
This specialized MATLAB analysis empowered the characterization of every track according to a variety of parameters including cell type, location, velocity, and persistence. Spatial and temporal correlations were computed with a Pearson correlation by using a moving average over distance and time respectively [35]. The local environment was determined by Voronoi tessellation based on these cell positions [13]. This allowed for an analysis of neighbor type, segregation parameter, and ultimately T1 transition dynamics.
Segregation was measured for pre-mixed systems, where cells were homogeneously seeded in both parts of the wound healing chambers and allowed to grow until confluence was reached, the wound healing insert was then removed after 24 h and observation generally lasted for more than 48 h, except for three experiments where visual inspection indicated a near halt of cell rearrangements. The initial segregation measurement for solid-fluid systems showed a higher variance, due to the vast difference in cell division rates between the two cell types. The segregation parameter quantifies the number of unlike neighbors over the total number of neighbors, normalized with the global population ratio, rendering it a measure for mixing of the system. In a 2D system each cell has on average 6 neighbors, so assuming a 50%/50% mixture a value close to 1 indicates a perfectly mixed system in which three neighboring cells are of the same type and three are unlike. A value of 0.666 specifies four like neighbors, 0.333 signifies on average five cells of the same type, whereas a value of 0 would point towards a completely unmixed system where neighbors are always the same type. Values lower than 0.333 are quite extreme in 2D systems, since interfaces are by definition composed of 50%/50% cell types and thus domain areas have to be vast to dominate over junctions.
The velocity direction of cells at interfaces was estimated by defining clusters of cells which contained at least one cell that has only neighbors of the same cell type. Cells within these clusters that also have neighbors of the other cell type are then defined to be the cluster boundary. The normal vector of the boundary for a particular cell and a certain time is estimated as the mean of the vectors from the particular cell and its neighbors of the same cell type to all its neighbors of the other cell type. Thus, all vectors between cells which are not neighbors to each other are discarded. The velocity of the cell in question is transferred to the basis defined by this estimation of the local normal vector of the cluster boundary. For the histograms in figures 5(D)-(F), the velocity components of cells with respect to their cluster boundary were averaged over time. The ratio in supplemental table S1 gives a lower bound for the estimation of the extent of directional motion, as all measurement errors involved are homogeneous in direction and therefore reduce the estimated velocity difference.
The cage-relative mean square displacement was determined through subtraction of the mean velocity of Voronoi-based neighbors from the cells' velocity in each frame and translation into pseudo coordinates for MSD calculation [20].
Turbulence intensity, as a measure for the irregularity of cell paths is a dimensionless parameter obtained via Reynolds decomposition and given by the root mean square of the velocity fluctuations normalized by the mean velocity (of the same cell type locally within 100 μm) [36].

Simulations
For the investigation of the effect of cell death and division, we assumed the cells to be elastic shells encompassing a fluid interior, and simulated them as Voronoi polygons. The model is a 2D adaptation of our previous 3D cell-cell interaction and tumor growth model [37].
The energy of each cell in a packed monolayer was calculated as a function of the elastic stretch in cell cortex (first term in equation (1) below), inter-cellular adhesion (second term in equation (1) below) and the work done against the osmotic pressure inside the cells (third term in the equation below-provides cell size control at homeostasis) where k is the cell cortex stiffness, and P i and P ci are the perimeters of the cell and of a circle with the same area as the cell respectively. γ j is the adhesion strength at the interface of contact between the cell i and cell j and l j is the length of the interface between the cells. λ is the internal pressure parameter and A i and A 0 are area of the cell and the mean area available to each cell (total area available to all cells divided by the number of cells) respectively.
Cellular rearrangements in the model system were driven by migratory forces which were assumed to be same for all cell types. The rearrangements were simulated using a Monte Carlo process where ≈1% of the cells were randomly displaced a short distance (<10% of suspended cell's diameter) during each iteration. Iterations/displacements that lower the total tissue energy ( U i ) were always accepted. On the other hand, iterations/displacements that increased the tissue energy were accepted with a diminishing probability based on the Metropolis algorithm. All cells were assigned characteristic death and division probabilities. These probabilities were defined to be cell morphology-specific: larger, elongated cells were assigned a higher likelihood of dividing and a lower chance of dying, while more rounded, compact cells were assigned higher likelihood of death and a low chance of division based on the equations below p death = p (P i −P ci )/(P m −P ci ) 0 (2) where p 0 is the baseline probability for death and division for the average cell. P m is the mean perimeter of all the cells. Cell perimeter is chosen as the key regulator of cell fate here as it is a good measure for cell size, cell shape and cell-cell interface area, all three of which have been shown to influence cell death and division rates [38][39][40]. Cell death was simulated by drastically lowering λ which slowly decreased the cell area and the space it occupied was consumed by neighboring cells. The dying cell was removed from the simulation once its area dropped below a threshold. Cell division was simulated by replacing a dividing cell by two daughter cells placed equidistantly along the long axis of the mother cell. We assumed that the daughter cells inherited the mechanical phenotype of the mother cells. Cell death and division events were interspersed between cell rearrangement iterations. This introduced dynamic features into the model as cells continuously changed their shape, which modified their likelihood for cell division or death. Using this simulation model, we started the simulations with initially evenly mixed systems and observed the effects of enabling or disabling these characteristic death and division processes on the spatial and temporal evolution of the multicellular system.