# Force-based three-dimensional model predicts mechanical drivers of cell sorting

## Abstract

Many biological processes, including tissue morphogenesis, are driven by cell sorting. However, the primary mechanical drivers of sorting in multicellular aggregates (MCAs) remain controversial, in part because there is no appropriate computational model to probe mechanical interactions between cells. To address this important issue, we developed a three-dimensional, local force-based simulation based on the subcellular element method. In our method, cells are modelled as collections of locally interacting force-bearing elements. We use the method to investigate the effects of tension and cell–cell adhesion on MCA sorting. We predict a minimum level of adhesion to produce inside-out sorting of two cell types, which is in excellent agreement with observations in several developmental systems. We also predict the level of tension asymmetry needed for robust sorting. The generality and flexibility of the method make it applicable to tissue self-organization in a myriad of other biological processes, such as tumorigenesis and embryogenesis.

### 1. Introduction

Self-organization is a widely studied feature of non-equilibrium systems across a great variety of fields [1–3]. In biology, self-organized processes range in length and time scales from protein folding [4] to the dramatic murmurations of starlings [5]. Generically, these systems exhibit emergence of macroscopic order from disordered states by the action of simple, local, inter-component interactions. Here, we focus on one such process of self-organization: cell sorting. Cell sorting in a multicellular aggregate (MCA) is critical for normal development of an embryo, and is also at the heart of pathologies such as metastatic growth of tumours.

The spontaneous separation of mixtures of embryonic cells has long been thought to be driven by cells’ differing ‘affinity’ for one another [6]. Two primary drivers have been proposed to facilitate differential affinity: differential adhesion and differential interfacial tension.

The differential adhesion hypothesis is a widely studied explanation of cell sorting [7–9]. It asserts that, in an MCA with multiple cell types, those with strong mutual adhesion are drawn together better than cells with weaker mutual adhesion. The latter are thus ‘sorted’ to the outside of the MCA. Although the differential adhesion hypothesis has been applied to a number of systems [10], it has been challenged in the light of observations that sorting occurs even when two cell types adhere to one another with equal strength [11], and that contractile forces are more significant than adhesion for tissue morphogenesis in a number of systems [12,13].

The idea that sorting in many biological processes is driven primarily by contractile forces is at the heart of the differential interfacial tension hypothesis [14], the basic postulates of which are: (i) cell surface membranes exhibit local tension variation [15]; (ii) a cell's surface tension is highest when in contact with the external medium and lowest when in contact with cells of the same type [16]. The variation in surface tension at cell interfaces leads to variation in the cell interface area, and hence to mutual affinity (figure 2*a*). According to the differential interfacial tension hypothesis, variations in mutual affinity, driven by differences in cortical tension, are primarily responsible for self-organization in a cell aggregate. The reduction in tension at cell interfaces is proposed to occur by the exclusion of myosin from the actomyosin cortex at the interface. Such exclusion is caused by the intra-membrane proteins that mediate adhesion between the cells [13]. This passive mechanism allows morphogenesis to proceed without active movements by the cells involved and, moreover, provides the local coordination between cells needed for complex global morphogenesis [17].

The importance of intercellular differential interfacial tension for tissue morphogenesis and maintenance has been demonstrated in a number of systems. For example, it is responsible for boundary maintenance in the dorsoventral compartment boundary of the *Drosophila* wing imaginal disc [18,19]. Much of this body of work has focused on two-dimensional epithelial systems, often maintaining boundaries rather than forming boundaries from a mixed aggregate [20]. However, further evidence of the importance of differential interfacial tension comes from experimental work on three-dimensional aggregates, suggesting that local variations in cortical tension are responsible for internalizing the first set of internal cells in the mouse morula [21]. Moreover, reduction in interfacial tension has also been shown to drive morula compaction [22] and allocation of cells to the inner cell mass of the embryo [23].

In order to investigate in detail the effect of differential interfacial tension on three-dimensional MCAs, we constructed a computational model based on the subcellular element method (SCEM) [24]. To validate the method, we compared its predictions to theoretical models of differential interfacial tension in cell doublets [13] (figure 2*c*). Using our model to simulate sorting in three-dimensional MCAs of two cell types, we were able to predict the minimum adhesion magnitude required for tension-driven sorting to occur. We show that, given this minimum adhesion, interfacial tension asymmetries drive sorting in MCAs. We also show that fairly significant differential adhesion is necessary for cell sorting. Ultimately, we also quantitatively predict the amount of tension asymmetry necessary for sorting.

### 2. Model

A number of techniques have been developed to model cellular systems [25–28], including explorations of differential interfacial tension using the vertex method in two-dimensional epithelia [19,29]. However, many of these methods rely on a global Hamiltonian, which is a strong, abstract assumption for modelling biological systems, which do not necessarily conserve energy. Instead, we set out to develop a computational model that assumes dynamics driven by local forces. We required a method that allows integration of subcellular-scale mechanics, including stiffness, tension and adhesion, with three-dimensional tissue-scale phenomena. Therefore, we chose to implement a model based on the SCEM.

SCEM is a method for multi-scale three-dimensional simulation that can simultaneously model intra- and inter-cellular dynamics, as well as the development of multicellular tissues. It allows us to study the effects of cell shape changes and complex inter- and intra-cell features on tissue-scale dynamics [30]. The basic SCEM method [24] treats each cell as a collection of elements (figure 1*a*) interacting via nearest-neighbour Morse potentials of the form $V(r)={D}_{e}({\text{e}}^{-2a(r-{r}_{e})}-2{\text{e}}^{-a(r-{r}_{e})}+1)$ [31]. Further details, including routines for cell growth and division, are discussed in electronic supplementary material, appendix A.

Nearest-neighbour elements of different cells also interact, facilitating intercellular interactions by the same mechanism. Different element types, however, can have different potential parameters. Element positions are updated at each time step according to these Morse potential interactions using overdamped Langevin dynamics: $\gamma \dot{\underset{\_}{x}}=-\mathrm{\nabla}U(\underset{\_}{x})+\underset{\_}{\zeta}$.

For the purposes of this work, we modified the original SCEM so that each cell possesses a contractile cortex, the properties of which are different from the main bulk of the cell. With this contractile cortex, we can independently control intercellular adhesion, cortical tension and interfacial tension. This modification involved first identifying the surface elements of the cell. This was done by dividing the cell into 32 sectors of equal steradians, emanating from its centre of mass, and defining as cortex type all elements in a sector with a radius greater than 80% that of the furthest element in that sector (figure 1*b*). The allocation of cortex elements is updated at each time step to accommodate element movement or creation. Cortex elements that move to the centre of the cell are reallocated cytoplasm type, and cytoplasm elements that move to the boundary are allocated cortex type. Once the cortex elements are defined, inter-cellular adhesive interactions are mediated exclusively between cortex elements in different cells. Cytoplasm elements have only repulsive interactions with elements in different cells, which prevents cells from overlapping in space.

To ensure that the interaction between cortex elements is tangential to the cell surface, we performed a Delaunay triangulation [32,33] over each cell’s cortex elements, which makes the cell a polyhedron of triangular facets. By its nature, this triangulation favours roughly equilateral triangles, which produces a relatively smooth cell surface (figure 1*c*; electronic supplementary material, appendix C). Two cortex elements sharing such a triangle edge are considered neighbours and they interact along the edge via a tension force. In cellular systems, this is an active pulling force driven by actomyosin contraction, and thus is treated as a force of constant magnitude that does not vary with element separation. The equilibrium separation of elements is defined by a sum of the bulk interaction forces between elements and the additional tension component. Thus defined, the cortex forces are guaranteed to span the full cell surface, be tangential to it, and enable the modelling of intercellular mechanical interactions. The Delaunay triangulation also prevents elements from detaching from the cell; any cytoplasm element that escapes will immediately be allocated cortex type and experience a returning force.

With our computational method, the cortical tension can vary locally across the cell surface in response to contact with surfaces of other cells. Such variation is achieved by changing the magnitude of intra-cell tension forces between cortex elements that lie at a cell interface. Specifically, we define a cell interface by first identifying the cortex elements that share an adhesive inter-cell interaction with a cortex element of a different cell. These elements are then labelled according to whether the other cell is of the same or different type. An interface between the two cells is defined as a region containing only cortex elements of the same label, and any cortical tension force acting between two cortex elements with the same label can be altered (figure 1*d*). The magnitude of interfacial tension can be specified differently depending on cell type and whether the interface is between like or unlike cell types. Although the factor by which tension is updated for a given interface depends upon the types of cells at the interface, the same factor is applied to interfacial tension in both of the two cells at the interface.

Changing the local interaction between cell cortex elements affects not only the local cortical tension but also the distance between the elements and, therefore, the local element density. Since intercellular adhesion is mediated by these elements, an increase in element density, for example, increases the local adhesion strength between neighbouring cells. To compensate for this, we normalize the adhesion magnitude by the local element density (electronic supplementary material, appendix D). Parameters controlling system behaviour in our simulations are outlined in electronic supplementary material, appendix B.

### 3. Simulation of cell doublets

The simplest experimental system to define the mechanical phenotype of cells and their interfaces is a cell doublet, for which *in vitro* experiments and theoretical models exist. This makes this system a suitable test case to validate our method. We expect sorting to be driven by changes in relative affinity, reflected by changes in equilibrium interfacial contact area (or, analogously, contact angle, which is more tractable to measure experimentally) between cells. This interfacial contact area depends upon the adhesion magnitude between cells (*A*_{M}) the cortical tension of the cells (*γ*_{m}) and the interfacial tension (*γ*_{c}).

Assuming a cell doublet is in mechanical equilibrium and has cylindrical symmetry, the dependence of the contact angle between two cells on their surface tensions can be deduced using linear force balance at the vertex [13] (figure 2*c*). Further assuming, based on experimental observations [13,34,35], that the effect of adhesion is negligible in the high-tension limit, one can derive a relationship between *γ*_{m}, *γ*_{c} and the doublet contact angle *θ*: $\beta =\mathrm{cos}(\theta /2)$ where $\beta ={\gamma}_{c}/{\gamma}_{m}$. From here, we call this model, which neglects adhesion and assumes that cell–cell affinity reflects a balance of cortical and interfacial tension, the linear force balance model.

To test the effect of system parameters on cell doublets in our simulations, we allow an initial cell to divide into two, with both daughter cells having the same properties, and specify a corresponding value of *β* for the interface formed between them. We then allow the system to reach mechanical equilibrium without growth or division, producing a doublet of identical cells, adhered at a shared interface (figure 2*b*). An alternative protocol—bringing two cells into contact and allowing an interface to form from nothing—was assessed to be less representative of behaviour in real MCAs.

Experimentally, a straightforward indicator of cell affinity in doublets is the contact angle between the two cells in a doublet, *θ*. While this measurement is intractable in our simulations, we could straightforwardly measure the cell–cell interface area as the sum of the areas of all the Delaunay triangles whose vertex elements all interact with elements in the neighbouring cell, and express this area as the ratio, *I*_{P} of the interface area to the total cell surface area. Using simple trigonometry to relate interfacial area to contact angle, *θ*, we found the corresponding prediction of the linear force balance model for *I*_{P}: *I*_{P} = (1 − *β*)/(3 − *β*). This allows direct comparison of our simulations to the linear force balance model at cell interfaces and provides a validation test of our method’s predictions (figure 2*d*). The value of *θ* can also be measured in experiments. The validation consisted of simulating cell doublets, from which we obtained measurements of *I*_{P} for values of *β* between 0.25 and 1. We define an adhesion magnitude *A*_{M} for our simulations. The value of *A*_{M} is the maximum force applied between elements across the interaction range of the adhesive potential. For these doublet simulations, we used a representative set of *A*_{M} and *γ*_{m} values, corresponding to low-tension and high-tension regimes. The resulting *I*_{P} values were then compared to the theoretical predictions of the linear force balance model (figure 2*d*).

We found very good agreement between the simulation results and the theoretical prediction of the force balance model in the high-tension regime (solid lines in figure 2*d*), with adhesion producing a small constant offset across *β*. The magnitude of this offset increases gradually with adhesion. By contrast, in the low-tension regime, adhesion cannot be ignored when considering the linear force balance in a cell doublet. Moreover, two cells with no adhesion cannot form an interface, and thus the linear force balance model is unsurprisingly not appropriate to use for low-adhesion conditions. Indeed, the simulation predicts a sharp drop in the interface at low adhesion. Most importantly, we conclude that due to the excellent agreement with the linear force balance model in the high-tension regime, our method can be used to predict the mechanical interactions within MCAs.

By simulating a large number of cell doublets at different parameters, we constructed phase diagrams for the behaviour of the cell doublet interface area as a function of *A*_{M} and *γ*_{m} for *β* = 0.5, 0.75 and 1.00 (figure 2*e*–*g*). The phase diagrams indicate that the highest *I*_{P} value achieved for any parameter set is approximately 0.32. This value is in good agreement with the theoretical limit for the interface between two hemispheres, which is exactly 1/3. Our doublet simulations also show that, for each value of *β*, the measurement of *I*_{P} drops sharply with increasing *γ*_{m}, tending to a minimum (figure 2*j*). The minimum is determined by *A*_{M}, with very low adhesions producing little to no interface. Increasing adhesion from zero produces a sharp increase in the interface, but beyond a threshold adhesion the increase of interface with adhesion magnitude has a much smaller gradient (figure 2*h*).

The doublet simulation results are even clearer if we take slices through the phase diagram (figure 2*h*–*j*). As expected, we find that, for both low and high tensions, little to no interface is produced at low adhesion values. However, for both low and high tension, the interface increases approximately linearly with increasing adhesion. In the low-tension regime (figure 2*h*), interfacial contact area increases sharply with adhesion until a value of approximately *A*_{M} = 0.1, after which the interfacial contact area plateaus. The value at which the interfacial contact area plateaus is approximately equal to the theoretical limit for two hemispheres (0.33) at low *β*. In the high-tension regime (figure 2*i*), we observe two adhesion regimes. Below a value of *A*_{M} = 0.2, the interfacial contact area increases slowly with adhesion; beyond that value, the interfacial contact area depends very weakly on adhesion.

Our doublet simulations demonstrate that our method is an accurate predictor of the mechanical interactions between cells. Furthermore, differences between simulation results and theoretical predictions further highlight the constraints of using the linear force balance model in situations for which adhesion cannot be neglected. Importantly, our method does not require the low adhesion assumption and can be used in any limit of adhesion or tension.

### 4. Quantitative predictions of sorting in multicellular aggregates

Having validated our method and used it to study cell doublet mechanics, we now demonstrate the influence of adhesion and interfacial tension on the sorting of two cell types in MCAs.

In order to assess quantitatively the effect of mechanics on the extent and speed of sorting in a cell aggregate, we developed three measures of sorting. These measures of sorting, as described below, are continuously calculated during the evolution of an MCA, and compared to a randomized reference distribution to produce a sorting index.

(1) | Radius measure. We can consider sorting in spherical MCAs, as studied here, to be the aggregation of one cell type within the MCA. This can be quantified by taking the mean radius of all cells of that type from the centre of mass of that cell type. This quantity is our first sorting measure, referred to as the radius measure, | ||||

(2) | Neighbour measure. In a sorted system, we expect a higher number of same-type cell neighbours. We consider two cells to be neighbours if at least one element in one cell interacts with an element of the other. It is via such interactions that neighbour cells apply forces to one another, but the cells are counted as neighbours even if the magnitude of this force is zero. The neighbour measure, | ||||

(3) | Surface measure. For spherical inside–outside sorting, we expect one cell type to occupy a greater proportion of the external surface of the MCA than the other. With our model, we are able to calculate the total cell surface area across all cells that is not part of a cell–cell interface. Therefore, for the surface measure, we calculate the fraction, |

To contextualize the measures described above, we introduced uniformly randomized control systems, against which the measures could be compared. To generate such controls, the type of all cells in a system was randomly reassigned at each data output interval, while maintaining the same number of each cell type and the spatial configuration of cells, at which point the sorting measures were calculated again for these randomized systems (figure 3*a*). Repeating this procedure 10^{4} times (electronic supplementary material, appendix E) produces a normal distribution of the sorting measures (figure 3*b*) from which we can calculate the mean, *E*(*X*_{i}), maximum, *X*_{i,max}, minimum, *X*_{i,min} and standard deviation, *σ*(*X*_{i}), of each sorting measure *X*_{i}. After these repeats, the MCA is returned to its original state and the simulation continues. The sorting measure as calculated from the simulation, *X*_{i}, is compared to the distribution of randomized measure values by calculating a sorting index, *S*_{i}: *S*_{i} = (*X*_{i} − *E*(*X*_{i}))/3*σ*(*X*_{i}). Using *σ* in the divisor rather than the full range of randomized system values ensures that highly deviant results in the randomized distribution do not overly affect the sorting index. Defined in this way, we expect the sorting indices to run roughly from 0 to 1, with values near 0 indicating a randomly mixed system, and values near 1 indicating a sorted system.

For all following sorting simulations, we used our method to simulate MCAs growing from 10 to 30 cells with two different cell types. Once the system reached 30 cells, simulations were stopped and the final state of the system at that point was analysed. We define cell type 1 to be that expected to sort to the inside of the MCA, and cell type 2 to be that expected to sort to the outside (figure 3*a*). We calculated the radius, neighbour and surface sorting measures, and hence the respective sorting indices, at regular intervals. We can define three interfacial tension factors: like–like interfaces between both cell types (*β*_{1,1}, *β*_{2,2}) and unlike interfaces (*β*_{1,2}). For simplicity in all simulations, we set *β*_{1,2} = *β*_{2,2} = 1.0 and vary *β*_{1,1}, henceforth simply referred to as *β*. Therefore, cell type 1 is the reference cell type for calculating the sorting indices. Both cell types were given the same cortical tension. Furthermore, we restricted ourselves to the high cortical tension limit, and define the dimensionless adhesion parameter *α* = *A*_{M}/*γ*_{m} to simulate the dynamics of MCAs for a wide range of values of *α* and *β*. We also assume, unless otherwise stated, the same adhesion magnitude for both cell types, for both like and unlike adhesion. Experimentally, *α* is difficult to measure, but it is thought to be in the range of 0.2–0.25 at the most [12,13,36], justifying our use of the high-tension regime as the biologically relevant regime.

To test also the effect of a division bias, we ran two types of simulations for each pair of *α* and *β* values: (i) a symmetric division, in which in which each division produced cells of the same type as the parent; (ii) an asymmetric division, each dividing cell had a 50% change of producing two daughter cells of the same type as the parent, and 50% chance of producing two daughter cells of different types. Each parameter set was repeated four times in order to estimate variability in the simulations.

We first used our method to test the hypothesis that doublet mechanics is a good predictor of sorting in an MCA, and then went on to use it to establish a better understanding of the mechanical drivers of cell sorting.

#### (a) Doublet mechanics predicts sorting in multicellular aggregates

One of the main issues impeding the experimental validation of the mechanical drivers of sorting in MCAs is the difficulty of measuring forces within an MCA. Therefore, cell doublets are often used to quantify the forces that can be used to predict sorting in an MCA. We thus sought to test the hypothesis that cell doublets are a good and sufficient model for MCAs.

Based on ranges of published results, we ran the simulations for several random values *β* less than or equal to 1.0 and *α* less than or equal to 0.4 for both symmetric and asymmetric division. We quantified the sorting index for the final state of 30 different systems. We compared these values with the results of doublet testing (figure 2) with the same parameters. In this way, we can assess how measurements taken from cell doublets predict the extent of sorting in MCAs. Figure 4 shows final sorting indices plotted against the difference in doublet interface proportion (Δ*I*_{P}). This difference is measured for like–like doublets of two cell types with the same parameters as the corresponding MCA. So, if the parameters used in the MCA simulation correspond to like–like doublet interface proportions of *I*_{P,11} for cell type 1 and *I*_{P,22} for cell type 2, where cell type 1 is expected to sort to the inside, Δ*I*_{P} = *I*_{P,11} − *I*_{P,22} for this MCA. The resulting scatter plots clearly show a nearly perfect positive correlation between the difference in interface proportion and the final sorting index across all sorting measures. This correlation is especially strong for systems with symmetric division. Importantly, this strong correlation establishes that experimentally tractable cell–cell affinity in a doublet, regardless of the underlying forces giving rise to that affinity, can be used to accurately predict the degree of sorting of MCAs. It also acts as further validation of the precision of our method in sorting MCAs.

#### (b) Mechanical drivers of cell sorting in multicellular aggregates

In the following, we will use our simulations to make quantitative predictions about what adhesion and tension asymmetries are necessary for MCAs to sort. We performed simulations over a wide range of *α* and *β* values, and constructed a phase diagram in the *α* − *β* parameter space from the final value of sorting indices calculated at the end of each simulation (figure 5*a*–*c*). From these phase diagrams, we can deduce the limits on the values of *α* and *β* that drive sorting.

Figure 5*a*–*c* demonstrates that there is good agreement between sorting indices. Moreover, the extent of sorting increases as *β* is reduced. This is predictable, given that a reduction of *β* is indicative of a relaxation of interfacial tension at the expense of surface cortical tension for cell type 1. Significantly, with our method, we can now quantitatively predict the maximum of *β* (≈0.6) necessary to drive sorting at experimentally realistic values of *α*. We also see in figure 5*d* that the kinetics of sorting are faster for smaller values of *β*, with very little sorting occurring for *β* = 0.7. It is also clear that a minimum adhesion is required in order for the *β* parameter to have any effect: for values of *α* < 0.2, little to no sorting is observed regardless of the value of *β*. This is interesting because *α* = 0.2 is in agreement with experimental observations showing adhesion magnitude being around 0.2 times cortical tension magnitude in systems where cortical tension has been proposed as a significant factor in morphogenesis [12,13].

We also use our method to explore the effects of differential adhesion. As seen in (figure 5*f*) differential adhesion can drive sorting, but requires a large difference in adhesion magnitude between the cell types—approximately a factor of 10 for full sorting. We also see in figure 5*e* that sorting by differential adhesion with *A*_{M,2}/*A*_{M,1} = *R* = 0.25 is slower than that driven by differential interfacial tension with *β* = 0.5. This suggests, in conjunction with the fact that adhesion forces are much smaller than tension forces, that tension asymmetry is a more significant factor in MCA sorting than adhesion asymmetry.

### 5. Discussion

We introduced a simulation method, an extension of SCEM, to investigate the mechanics of cell sorting. The model is based on local intra- and inter-cellular forces, and includes cell growth and division, and therefore evolution of different size MCAs. It also allows us to define different regions of a cell and the mechanical properties of those regions. Therefore, our method is highly flexible and can model intercellular and intracellular mechanical interactions in a broad range of multicellular systems.

We first employed the experimentally tractable cell doublet system both to validate our numerical model and to probe the effect of adhesion and tension on the interfacial contact area between two cells. In a doublet, this contact area is determined by the balance of the adhesion force, cortical tension, and interfacial tension. We studied the effects of these parameters on doublet interfacial contact area and found that it decreases sharply with interfacial tension before reaching a stationary phase in the high-tension regime. We also found that no interface can form when the adhesion drops below a threshold, above which the interfacial contact area depends only weakly on adhesion. Informed by experimental observations that tension forces are at least four times stronger than adhesion [12,13], we studied the high-tension regime and found that the variation of interfacial contact area with the tension parameter, *β*, agrees remarkably well with the predictions of using local force balance at cell interfaces, neglecting adhesion. This validates our method, showing that it captures important aspects of the underlying mechanics.

Next, we simulated the evolution of MCAs, growing from 10 to 30 cells, to study the effects of adhesion and tension on the sorting process. Using the novel, unbiased sorting index that we developed, we observed a very strong correlation between the sorting indices in the MCAs and the interfacial contact area between cells in doublets. This observation is significant because, while experimental determination of mechanics within MCAs is difficult, measurements of cortical tension and contact angle in cell doublets are relatively straightforward. Perhaps just as importantly, this observation confirms that even though interfacial tension, cortical tension and adhesion can all affect cell–cell affinity, it is only the combination of their contributions in the form of cell–cell affinity that affects cell sorting. Furthermore, these plots (figure 4*a*–*c*) suggest that a relatively large difference in doublet interface area between the two cell types is required to produce complete sorting reliably.

Finally, in order to disentangle the relative effects of differential interfacial tension and differential adhesion, we constructed phase diagrams for the sorting indices as a function of *α* and *β*. These phase diagrams show that DIT alone is sufficient to produce full sorting of MCAs with a sufficiently small *β*. Importantly, measurements of this experimentally accessible parameter in cell doublets can be used as a test of this prediction. Furthermore, we found that there is a lower bound for adhesion necessary for sorting of an MCA which is in agreement with experimentally observed values. In contrast to differential interfacial tension, the amount of differential adhesion necessary to drive cell sorting alone is comparatively large (approx. a factor of 10). This suggests that tension asymmetries are the primary driver of cell sorting.

Our results lead us to conclude that, for relatively small MCAs, similar to the sizes of those seen in early development, sorting can occur only if there are relatively large differences in interfacial tension, very large differences in adhesion or both between the different cell types. This could have vast implications in, for example, the early developing embryo where only a few cells are present in an MCA.

Our method can be extended in several ways. For example, one can endow intra-cellular regions with different properties, making it possible to model the nucleus and the plasma membrane. Even more generally, this method can be extended straightforwardly to test aggregates that are not confined to roughly spherical geometries, many of which occur in biological systems. For these reasons, our method can be useful well beyond studying tissue morphogenesis and could be applied to investigate a range of biological processes and pathologies, such as cancer cell metastasis.

### Data accessibility

Computational code and accompanying comments can be found at https://github.com/chris-revell/SEM.

### Competing interests

We declare we have no competing interests.

### Funding

This work was financially supported by the Medical Research Council and the Wellcome Trust. K.J.C. is a University Research Fellow of the Royal Society, which also provided financial support for this work. This work was performed using the Darwin Supercomputer of the University of Cambridge High Performance Computing Service (http://www.hpc.cam.ac.uk/), provided by Dell Inc. using Strategic Research Infrastructure Funding from the Higher Education Funding Council for England and funding from the Science and Technology Facilities Council.

## Acknowledgements

The authors thank Tim Newman for providing the basic SCEM code for updating, and Ewa Paluch for valuable advice and input.

### References

- 1.
Ashby WR . 1947 Principles of the self-organizing dynamic system.**J. Gen. Psychol.**, 125-128. (doi:10.1080/00221309.1947.9918144) Crossref, PubMed, ISI, Google Scholar**37** - 2.
Bak P . 1996**How nature works: the science of self-organized criticality**. New York, NY: Copernicus. Crossref, Google Scholar - 3.
Camazine S . 2003**Self-organization in biological systems**. Princeton, NJ: Princeton University Press. Google Scholar - 4.
Dobson CM . 2003 Protein folding and misfolding.**Nature**, 884. (doi:10.1038/nature02261) Crossref, PubMed, ISI, Google Scholar**426** - 5.
Couzin ID, Krause J . 2003 Self-organization and collective behavior in vertebrates.**Adv. Study Behav.**, 1-75. (doi:10.1016/S0065-3454(03)01001-5) Crossref, ISI, Google Scholar**32** - 6.
Townes PL, Holtfreter J . 1955 Directed movements and selective adhesion of embryonic amphibian cells.**J. Exp. Zool.**, 53-120. (doi:10.1002/(ISSN)1097-010X) Crossref, Google Scholar**128** - 7.
Steinberg MS . 1962 On the mechanism of tissue reconstruction by dissociated cells, I. Population kinetics, differential adhesiveness, and the absence of directed migration.**Proc. Natl Acad. Sci. USA**, 1577-1582. (doi:10.1073/pnas.48.9.1577) Crossref, PubMed, ISI, Google Scholar**48** - 8.
Steinberg MS . 1962 Mechanism of tissue reconstruction by dissociated cells. II. Time-course of events.**Science**, 762-763. (doi:10.1126/science.137.3532.762) Crossref, PubMed, ISI, Google Scholar**137** - 9.
Steinberg MS . 1962 On the mechanism of tissue reconstruction by dissociated cells, III. Free energy relations and the reorganization of fused, heteronomic tissue fragments.**Proc. Natl Acad. Sci. USA**, 1769-1776. (doi:10.1073/pnas.48.10.1769) Crossref, PubMed, ISI, Google Scholar**48** - 10.
Cochet-Escartin O, Locke TT, Shi WH, Steele RE, Collins ES . 2017 Physical mechanisms driving cell sorting in*Hydra*.**Biophys. J.**, 2827-2841. (doi:10.1016/j.bpj.2017.10.045) Crossref, PubMed, ISI, Google Scholar**113** - 11.
Moore R, Cai KQ, Escudero DO, Xu X . 2009 Cell adhesive affinity does not dictate primitive endoderm segregation and positioning during murine embryoid body formation.**Genesis**, 579-589. (doi:10.1002/dvg.v47:9) Crossref, PubMed, ISI, Google Scholar**47** - 12.
Chan EHY, Shivakumar PC, Clément R, Laugier E, Lenne P . 2017 Patterned cortical tension mediated by N-cadherin controls cell geometric order in the Drosophila eye.**eLife**, e22796. (doi:10.7554/eLife.22796) Crossref, PubMed, ISI, Google Scholar**6** - 13.
Maître JL, Berthoumieux H, Krens SFG, Salbreux G, Jülicher F, Paluch E, Heisenberg CP . 2012 Adhesion functions in cell sorting by mechanically coupling the cortices of adhering cells.**Science**, 253-256. (doi:10.1126/science.1225399) Crossref, PubMed, ISI, Google Scholar**338** - 14.
Harris AK . 1976 Is cell sorting caused by differences in the work of intercellular adhesion? A critique of the Steinberg hypothesis.**J. Theor. Biol.**, 267-285. (doi:10.1016/0022-5193(76)90019-9) Crossref, PubMed, ISI, Google Scholar**61** - 15.
Chalut KJ, Paluch EK . 2016 The actin cortex: a bridge between cell shape and function.**Dev. Cell**, 571-573. (doi:10.1016/j.devcel.2016.09.011) Crossref, PubMed, ISI, Google Scholar**38** - 16.
Brodland GW . 2002 The differential interfacial tension hypothesis (DITH): a comprehensive theory for the self-rearrangement of embryonic cells and tissues.**J. Biomech. Eng.**, 188-197. (doi:10.1115/1.1449491) Crossref, PubMed, ISI, Google Scholar**124** - 17.
Lecuit T, Lenne P, Munro E . 2011 Force generation, transmission, and integration during cell and tissue morphogenesis.**Annu. Rev. Cell Dev. Biol.**, 157-184. (doi:10.1146/annurev-cellbio-100109-104027) Crossref, PubMed, ISI, Google Scholar**27** - 18.
Landsberg KP, Farhadifar R, Ranft J, Umetsu D, Widmann TJ, Bittig T, Said A, Jülicher F, Dahmann C . 2009 Increased cell bond tension governs cell sorting at the*Drosophila*anteroposterior compartment boundary.**Curr. Biol.**, 1950-1955. (doi:10.1016/j.cub.2009.10.021) Crossref, PubMed, ISI, Google Scholar**19** - 19.
Aliee M, Röper J, Landsberg KP, Pentzold C, Widmann TJ, Jülicher F, Dahmann C . 2012 Physical mechanisms shaping the*Drosophila*dorsoventral compartment boundary.**Curr. Biol.**, 967-976. (doi:10.1016/j.cub.2012.03.070) Crossref, PubMed, ISI, Google Scholar**22** - 20.
Umetsu D, Dahmann C . 2010 Compartment boundaries: sorting cells with tension.**Fly**, 241-245. (doi:10.4161/fly.4.3.12173) Crossref, PubMed, ISI, Google Scholar**4** - 21.
Samarage CR, White MD, Álvarez YD, Fierro-González JC, Henon Y, Jesudason EC, Bissiere S, Fouras A, Plachta N . 2015 Cortical tension allocates the first inner cells of the mammalian embryo.**Dev. Cell**, 435-447. (doi:10.1016/j.devcel.2015.07.004) Crossref, PubMed, ISI, Google Scholar**34** - 22.
Maître J, Niwayama R, Turlier H, Nédélec F, Hiiragi T . 2015 Pulsatile cell-autonomous contractility drives compaction in the mouse embryo.**Nat. Cell Biol.**, 849-855. (doi:10.1038/ncb3185) Crossref, PubMed, ISI, Google Scholar**17** - 23.
Maître J, Turlier H, Illukkumbura R, Eismann B, Niwayama R, Nédélec F, Hiiragi T . 2016 Asymmetric division of contractile domains couples cell positioning and fate specification.**Nature**, 344. (doi:10.1038/nature18958) Crossref, PubMed, ISI, Google Scholar**536** - 24.
Newman TJ . 2007 Modeling multicellular structures using the subcellular element model. In**Single-cell-based models in biology and medicine**(edsAnderson ARA, Chaplain MAJ, Rejniak KA ), pp. 221-239. Berlin, Germany: Springer. Crossref, Google Scholar - 25.
Paluch EK . 2015 After the greeting: realizing the potential of physical models in cell biology.**Trends Cell Biol.**, 711-713. (doi:10.1016/j.tcb.2015.08.001) Crossref, PubMed, ISI, Google Scholar**25** - 26.
Tanaka S . 2015 Simulation frameworks for morphogenetic problems.**Computation**, 197-221. (doi:10.3390/computation3020197) Crossref, ISI, Google Scholar**3** - 27.
Anderson A, Rejniak K 2007**Single cell-based models in biology and medicine**. Basel, Switzerland: Birkhäuser. Crossref, Google Scholar - 28.
Osborne JM, Fletcher AG, Pitt-Francis JM, Maini PK, Gavaghan DJ . 2017 Comparing individual-based approaches to modelling the self-organization of multicellular tissues.**PLoS Comput. Biol.**, e1005387. (doi:10.1371/journal.pcbi.1005387) Crossref, PubMed, ISI, Google Scholar**13** - 29.
Schwarz US, Dunlop CM . 2012 Developmental biology: a growing role for computer simulations.**Curr. Biol.**, R441-R443. (doi:10.1016/j.cub.2012.04.038) Crossref, PubMed, ISI, Google Scholar**22** - 30.
Sandersius SA, Newman TJ . 2008 Modeling cell rheology with the subcellular element model.**Phys. Biol.**, 15002. (doi:10.1088/1478-3975/5/1/015002) Crossref, PubMed, ISI, Google Scholar**5** - 31.
Morse PM . 1929 Diatomic molecules according to the wave mechanics. II. Vibrational levels.**Phys. Rev.**, 57-64. (doi:10.1103/PhysRev.34.57) Crossref, Google Scholar**34** - 32.
Renka RJ . 1997 Algorithm 772: Stripack: delaunay triangulation and voronoi diagram on the surface of a sphere.**ACM Trans. Math. Softw. (TOMS)**, 416-434. (doi:10.1145/275323.275329) Crossref, ISI, Google Scholar**23** - 33.
Delaunay B . 1934 Sur la sphere vide.**Izv. Akad. Nauk SSSR, Otdelenie Matematicheskii i Estestvennyka Nauk**, 1-2. Google Scholar**7** - 34.
Susienka MJ, Wilks BT, Morgan JR . 2016 Quantifying the kinetics and morphological changes of the fusion of spheroid building blocks.**Biofabrication**, 045003. (doi:10.1088/1758-5090/8/4/045003) Crossref, PubMed, ISI, Google Scholar**8** - 35.
Youssef J, Nurse AK, Freund LB, Morgan JR . 2011 Quantification of the forces driving self-assembly of three-dimensional microtissues.**Proc. Natl Acad. Sci. USA**, 6993-6998. (doi:10.1073/pnas.1102559108) Crossref, PubMed, ISI, Google Scholar**108** - 36.
Amack JD, Manning ML . 2012 Knowing the boundaries: extending the differential adhesion hypothesis in embryonic cell sorting.**Science**, 212-215. (doi:10.1126/science.1223953) Crossref, PubMed, ISI, Google Scholar**338**