Cell sorting by active forces in a phase-field model of cell monolayers

Cell sorting, the segregation of cells with different properties into distinct domains, is a key phenomenon in biological processes such as embryogenesis. We use a phase-field model of a confluent cell layer to study the role of activity in cell sorting. We find that a mixture of cells with extensile or contractile dipolar activity, and which are identical apart from their activity, quickly sort into small, elongated patches which then grow slowly in time. We interpret the sorting as driven by the different diffusivity of the extensile and contractile cells, mirroring the ordering of Brownian particles connected to different hot and cold thermostats. We check that the free energy is not changed by either partial or complete sorting, thus confirming that activity can be responsible for the ordering even in the absence of thermodynamic mechanisms.


Introduction
Cell sorting, the separation of mixtures of different cell-types into distinct domains reminiscent of phase ordering in fluids, has long been a topic of interest in the biological literature [1][2][3] .It is a vital component of embryogenesis and morphogenesis, relevant to how cells organise into different cell types prior to differentiation, and is conserved across vertebrates and invertebrates.Model systems, for example plated confluent layers which include different cell types also often show localised de-mixing 2,4,5 .A number of mechanisms have been suggested to explain cell sorting in tissues, and investigated both in vitro and in silico.These can appeal to differences in cell adhesion, line tension, or activity [6][7][8] .Other characteristics of a monolayer such as cell motility 9 or area and shape 10 may also play a role in cell sorting.
The Differential Adhesion Hypothesis introduced by Steinberg 2 proposes that a cell preferentially adheres to cells of the same, rather than another, type.Differential adhesion leads to variation in tissue surface tension, and results in a thermodynamic separation akin to that of oil and water.The earliest mechanism proposed to explain cell sorting, the DAH gained currency and has been investigated experimentally in cell aggregates expressing varying levels of cadherins 6 .The hypothesis has been used to describe the sorting of heterogeneous monolayers in vitro 8 , and in silico through agent-based and cellular-Potts model simulations 11,12 .
Complementary to the DAH is the Differential Interfacial Tension Hypothesis, which considers the tension at cell-cell junctions.This hypothesis holds that line tension can modulate and regulate cell sorting, and even permit one tissue to completely surround another 13 .It has been shown that cell junction tension controls cell sorting in Drosophila 7 .In addition, fluctuations in line tension in a vertex model can result in intercalations 14 and cell sorting 15 .
Experiments by Skamrahl et al. show that differences in the line tension at junctions between cells of different types, which the authors refer to as contractility, has an interplay with differential adhesion during cell sorting 8 .In this work, the authors compare cultures of wild-type Madin-Darby Canine Kidney cells to 'dKD' cells, which have a protein knockdown that enhances contractility.To analyse cell sorting, they tag wild-type cells with green fluorescent protein and mix them with either dKD cells or untagged wild-type cells.The mixture of dKD and wild-type cells separates into patches of each type.The authors also observe that demixing proceeds in a two-stage process: a fast stage followed by a slow stage.They attribute the fast time scale to differential contractility, and the slow time scale to differential adhesion.
Other factors have also been implicated in cell sorting in a monolayer.Heine et al. observe the influence of the solid-or fluid-like nature of a cell culture on sorting within the differential adhesion framework 9 .There are also suggestions that the apical surface area of epithelial cells can affect sorting 8,10 .
Moreover, living systems are active.They continuously use chemical energy to drive motility and growth, and hence remain out of thermodynamic equilibrium.Recent work has successfully interpreted collective processes in mechanobiology in terms of the theories of active matter, and in particular of active nematics.Balasubramaniam et al. 4 observed incomplete cell sorting in a mixture of two MDCK cell strains with different strengths of inter-cellular interactions.They interpreted these results in terms of a mixture of dipolar active extensile and active contractile cells and complemented the experimental results with simulations of a phase-field model.A two-fluid continuum model has also been used to show that fluids of different activities can undergo microphase separation 16 .
Thus many different mechanisms can result in a degree of cell sorting, and the interplay of several may be necessary for complete, well-controlled segregation in a given biological process 5,9 .To obtain a more complete understanding of sorting phenomena in cell monolayers and tissues, it is crucial to investigate each of the mechanisms involved.Isolating individual processes is difficult in experiment but modelling can play a useful role in this regard.Therefore here we implement a phase-field model and concentrate on asking whether active forces acting between cells, which arise from cadherin-based intercellular junctions, can result in sorting in a confluent cell monolayer.
In the next section, we introduce the phase-field model, including both the passive dynamics of each cell and the active forces between cells.In section 3 we present results indicating that cells of different activities undergo localised sorting.Our results are summarised and discussed in section 4.

Phase-Field Model
Phase-field models are simplified approaches to cellular dynamics [17][18][19][20][21][22] which have been used to model the motility of single and multiple cells and of confluent cell layers.A particular advantage of phase-field approaches is that it is possible to distinguish passive and active contributions to cell dynamics, and to implement both intracellular forces and forces between cells mediated through cell junctions in a transparent way.
Each cell i is represented by a phase field, ϕ (i) (x).The phase field takes a value 0 outside the cell and 1 in the interior and varies smoothly across the cell boundary.The phase field ϕ (i) (x) moves with a velocity v (i) (x), and evolves according to where F is a free energy.The right-hand side of Equation (1) describes the dissipative relaxation dynamics of the cells to a free energy minimum at a rate J 0 .Assuming over-damped dynamics, the velocity v (i) (x) of a cell is determined by the local force density acting on the cell, where we distinguish active and passive contributions to the force density and ξ is a friction coefficient.The force balance holds everywhere in the domain: any force imbalance is countered by the dynamic friction which arises from substrate interactions, representing focal adhesions 23 .The passive force density is modelled by three free energy contributions: The first of these is a Cahn-Hilliard term, which favours phase separation into regions of ϕ (i) = 0, 1 with an interface width λ .The parameter γ controls the line tension for the interface of a single cell, which governs how easily the cell can deform.The second free energy term imposes a soft area constraint with energy scale µ, where πR 2 is the equilibrium area of a cell and R its target radius.The third term penalises cell-cell overlap with an energy scale κ.The length unit in the system is equal to the lattice unit It is next necessary to define the active forces acting in the system.The polar motility of individual epithelial cells arises from lamellipodia and other protrusions 24,25 that generate polar forces.However, it has been argued that the formation of lamellipodia is suppressed in confluent epithelia, a phenomenon known as contact inhibition of locomotion 26 , suggesting a reduction in the strength of any persistent polar force in the monolayer.There is evidence that individual epithelial cells can create dipolar stresses 4 , and many cell monolayers exhibit active flows very similar to those predicted by active nematic models 19 .Therefore we choose to concentrate on dipolar (balanced) active forces.
To calculate the active dipolar contribution to the force density, f (i) active (x), we first calculate the deformation tensor that quantifies the shape of each cell 19,27 , where ⊥ are the orthonormal eigenvectors of the D (i) , along and perpendicular to the elongation axis of the cell respectively, normalised so that d We next define a director, n (i) , associated with each cell i and assume that n (i) relaxes towards d where J n controls the time scale of relaxation.We assume no noise in the relaxation of the director.
In the absence of any unbalanced active forces, the dipolar contribution to the active stress acting on cell i is split into a selfinduced stress and a stress due to the cell's neighbours, related to the director field by 28,29 where each cell is assigned a nematic tensor The force density arising from the active stress is then ∥ at rate J n = 0.1.The timescale for the shape of a cell to relax according to the Cahn-Hilliard phase separation term is (γ∆ℓJ 0 ) −1 ∼ 10 2 , where ∆ℓ = 1 is the lattice unit.
The phase-field model was implemented with periodic boundary conditions.The system was initiated by placing cells randomly within relevant areas of the simulation domain to give a confluent layer.The system's initial configuration was obtained by relaxing the randomly-placed cells as though they were passive for 5 × 10 3 time steps, at which point activity was turned on (denoted by t = 0) and data were recorded to t = 1×10 5 or t = 5 × 10 5 timesteps.The dynamical equations were solved using a forward Euler scheme with one predictor-corrector step.The unit of time is 1, with each step sub-divided into five substeps, and the unit of length is the lattice spacing.We note that the model differs from the first implementation of this software 19 by the omission of an adhesive contribution to the free energy and by the coupling of Q (i) not directly to D (i) but to its principal eigenvector d (i) ∥ .Intercellular interactions, treated separately from intracellular dipolar stresses, have been investigated more recently 29 , and the treatment here follows this recent work.

Results
We concentrate primarily on dipolar intercellular interactions, choosing ζ self = 0 and returning to consider intracellular forces at the end of this section.To investigate whether cell sorting can be driven by differences in active inter-cellular forces, we prepared an equal mixture of cells with extensile and contractile activities, ζ inter = 0.4 and ζ inter = −0.4,respectively.All other simulation parameters were identical for both cell types.The simulation was started by placing 672 cells of each type randomly in the simulation domain which measured 560 by 560 lattice units.This number density corresponds to a cell density where the cell layer is confluent, but cell motion is still possible.Although a 'packing fraction' is not uniquely defined for soft particles, we describe the target packing fraction as the total area fraction of cells divided by the area of the domain, which is ∼ 0.86.  1 (c), (d).The cells show no tendency to mix on the time scale of the simulation (here to t = 1 × 10 5 ), therefore both microphase and macrophase separation persist for long times.Note that contractile cells are elongated slightly by the dipolar forces and demonstrate nematic alignment.We seek to quantify the phase separation and describe its mechanism.There are a number of ways to measure the segregation in a system 4,8,11,30 .Following 8 , we define a segregation index by counting the neighbours of each cell.Two phase-field cells are identified as neighbours when the corresponding phase fields take a value greater than a certain threshold, set here to be 0.1, on the same lattice site.The segregation index is then defined by where n is the number of neighbours of a cell with the same activity, n is the number of neighbours of a cell with opposite activity and the average is taken over all cells.With our choice of parameters, the phase-field cells are stable and have sharp interfaces so no spurious neighbour pairs are identified.
The segregation index for a 1:1 mixture is plotted in Figure 2  (a).The microphase separation proceeds quickly after the onset of dipolar activity at t = 0 but then slows.It is unclear from the data whether microphase separation is arrested or continues over very long time scales, since the SI appears to saturate at late times.
An alternative measure of phase separation is the density autocorrelation function where the density field ρ ζ for the subset of cells in the system with ζ inter = ζ is defined by Figure 2 shows the moving average of the lengthscale of the density autocorrelation ρ ζ , taken over 10 4 timesteps.Recalling that each cell has a nominal diameter ∼16 lattice units, the length scale of the density autocorrelation is initially roughly one cell, consistent with a well-mixed monolayer, while the length scale at long times increases to slightly more than 4 cells, consistent with microphase separation into clusters roughly 4 cells wide.These results demonstrate quantitatively that the mixture of extensile and contractile cells sorts according to activity.The data suggest that the cluster size saturates, but we cannot rule out a continued extremely slow growth of the clusters.The segregation index and density autocorrelation indicate microphase separation but yield no information about the dynamics of the system.We next examine the mean-square displacement of the phase-field cells, comparing the fully phase separated state (Fig. 1d) to the microphase-separated state (Fig. 1b).
Results for the phase-separated configuration shown in Fig. 3a indicate that, when surrounded by cells of the same type, the motion of the contractile cells is almost entirely arrested, while extensile cells move more freely.This is true also in pure monolayers: Figures SI3 and SI4 show, respectively, snapshots of layers and plots of mean-square displacement for homogeneous monolayers of extensile and contractile cells.This is a result of the different intercellular interactions due to the active dipolar forces.Contractile cells elongate and align to form a solid-like, nematic configuration.Extensile cells prefer to lie at right angles which leads to frustration and gives unstable configurations.Pairs of cells which lie perpendicular to each other have a polarity which results in net migration 29 .
In the microphase separated system the diffusion of contractile cells is enhanced and that of extensile cells reduced (Fig. 3b).This is because less motile clusters of contractile cells constrain the paths of extensile cells, while the extensile cells tend to push the contractile cells around the simulation domain.As part of the diffusion, individual extensile cells are able to squeeze between contractile clusters from one extensile cluster to another.Given the apparent saturation of the SI and the density autocorrelation lengthscale at long times, the system may be entering a dynamic steady state.It is unclear whether the layer can sort completely, as in Figure 1 (d), at much longer times.
Taking the diffusion constant in the system to be that of the extensile cells, ⟨x 2 ⟩ ζ =+0.4 = 4Dt, the timescale associated with diffusion is τ D = R 2 /D ∼ 10 4 .Madin-Darby Canine Kidney cells migrate on the order of 10 microns per hour 31 , which is on the order of a cell diameter.The diffusion timescale in the 1 : 1 mixed system suggests that 1 × 10 5 timesteps in silico corresponds to order a day of real time, which is the duration of the experiments of Balasubramaniam et al. 4 .
We note that differential diffusivity has been shown to sort otherwise passive soft particles 32 , even at high packing fractions 33 .
Here, extensile cells are akin to hot particles and contractile cells to cold ones.Weber, Weber and Frey 32 attribute the phase separation in a mixture of particles at two different temperatures to an effective attraction between cold particles; here the active forces act to attract the contractile cells into more coherent clusters.
Fig. 3 Mean-square displacement for contractile (blue) and extensile (red) cells in (a) a 1:1 phase separated monolayer (Fig. 1d) (b) a 1:1 microphase-separated state (Fig. 1b).When the cells are completely phase separated, the extensile cells behave as a fluid while the contractile cells behave as a solid; extensile cells diffuse three orders of magnitude as quickly as contractile cells.In the microphase-separated state, however, the presence of contractile cells slows the diffusion of extensile cells, which in turn push the contractile cells around the system.The MSD has not saturated by t = 5 × 10 5 , which indicates that the clusters are still evolving and rearranging.
In systems that sort, there is often a thermodynamic basis for the phase separation.Consequently, we investigate whether the Table 1 Mean Ftot and standard deviation of free energy per cell for the microphase-separated state (Fig. 1a,b) and the sorted state (Fig. 1c,d) over time, starting at t = 1 × 10 3 , shortly after activity is turned on.The clustering of the phase fields over time is not associated with a reduction in free energy at long times, in contrast to phase separation according to mechanisms such as differential adhesion.clustering is associated with any changes in the total free energy of our model, Equations 4. For the microphase separated system (Fig. 1a,b), the mean and standard deviation of the total free energy per cell from times t = 1 × 10 3 to t = 5 × 10 5 are listed in Table 1.As a comparison, data for the system with phaseseparated initial condition (Fig. 1c,d) from times t = 1 × 10 3 to t = 1×10 5 are also listed.In both cases, the times run from shortly after the turning on of dipolar active stresses to the end of the simulation.The free energy shows no dependence on time or, indeed, cluster size, giving additional support to an active origin to the cell sorting.

Microphase-Separated
Thus far, we have considered mixtures of extensile and contractile cells, with ζ inter = +0.4,−0.4.It is possible, in addition, to examine active-passive systems.An extensile-passive system, with ζ inter = +0.8,0.0, fails to phase separate because the passive cells have no active forces which can hold them together (Figs.SI5 and  SI6) -the contractile cells form a solid-like phase in addition to being less diffusive than extensile cells.A contractile-passive system, with ζ inter = −0.8,0.0, also fails to phase separate (Figs.SI7  and SI8).In this system, contractile cells are not able substantially to deform each other and the geometry remains frustrated.Any randomly-initialised patches of contractile cells adopt a nematic configuration, while any passive cells not subject to contractile stresses remain isotropic.The failure of a mixture of active and passive cells to phase separate is consistent with simulations of active Brownian particles mixed with passive particles 34 .
Finally we mention intracellular dipolar activity ζ self ̸ = 0 (Figs.SI9, SI10 and SI11).As observed in 29 , values of ζ self and ζ inter with the same sign tend to cancel.This reduction in effective activity reduces the diffusivity of the active phase-field cells and slows sorting.

Discussion
There is considerable experimental evidence of sorting of different cell types both in vitro, in confluent cell layers 1,2,4,8,30 or cellular spheriods 5 , and in vivo, for example during morphogenesis 7,26 .Many different physical factors are candidates for driving the sorting.Without doubt equilibrium effects, such as the dependence of cell-cell adhesion or line tension of cell-cell junctions on different cell neighbours, can lead to ordering akin to thermodynamically-driven phase ordering in liquid-liquid mixtures 6,13,15 .However, biological systems are naturally out of equilibrium, and it has also been suggested that different forms of activity can result in the sorting of different cell types 4,35 .Balasubramaniam et al. in particular show that cells with a mixture of extensile and contractile force dipoles are able to phase separate on short length scales 4 .
Isolating the different contributions to cell sorting is difficult in experiments, but easier in the context of computational models of cell motility.Therefore in this paper we use a multi-phase field model of a confluent cell layer to study the influence of dipolar active interactions between neighbouring cells on cell sorting.We demonstrate that a mixture of extensile and contractile active dipolar cells, which are otherwise identical, can undergo partial sorting.We interpret this as an out-of-equilibrium effect resulting from the different dynamics of the two cell populations.Extensile cells are smaller and more motile whereas contractile cells tend to elongate and form static, solid-like nematic patches.We caution, however, that both the microphase separated state and a fully sorted state are stable on the timescale of the simulations, and it is unclear whether further coarsening will occur on times we cannot access.
Any model of cell mechanics must still be viewed with caution as there are still many questions about the model details needed to faithfully reproduce the cells' interactions and dynamics.Here we have focussed on forces, mediated by adherens junctions between cell cortices, that act across cell-cell boundaries, which we have modelled as balanced dipolar forces.However, fluctuating polar forces 36 or active forces which act along cellcell junctions 14 may be relevant, as may apical-basal asymmetry if a monolayer is modelled in three dimensions 37 .

( 9 )
Biologically, ζ self is a measure of the magnitude of stresses linked to intracellular myosin motors, while ζ inter indicates the strength of intercellular stresses mediated by adherens junctions acting between neighbouring cell cortices.The parameters assigned to the free energy terms (4) were γ = 1.4,λ = 2.0, µ = 120 and κ = 1.5.The coefficient of friction in Equation (2) was ξ = 3.0 and the target cell radius was R = 8.0.Relaxation according to Equation (1) was controlled by the rate J 0 = 5 × 10 −3 .Each cell's nematic director relaxed towards its long axis d (i)

Figure 1 (
Figure 1(a) and (b) illustrate the configurations at t = 1000 and after t = 5 × 10 5 time steps.The snapshots show evidence of par-

Fig. 1
Fig. 1 Snapshots of segregation in a 1:1 mixture of extensile (red) and contractile (blue) cells.(a) t = 1000, (b) t = 5×10 5 for a fully mixed initial condition showing microphase separation.(c) t = 1000, (d) t = 1×10 5 for a state that is fully sorted at the beginning of the simulation and which does not mix.Note that the contractile cells develop nematic order.

Fig. 2
Fig. 2 Segregation in a 1:1 mixture of extensile and contractile cells.(a) Segregation index SI versus time.The onset of segregation is on a timescale t ∼ 10 3 , and .The SI increases past ∼ 0.7 on the timescale of the simulation; the inset showing SI on log-log axes indicates the system coarsens steadily.(b) Time-averaged correlation length of ρ ζ for extensile and contractile cells in the simulation illustrated in Fig. 1 (a) and (b).The correlation length grows up to t ≈ 3 × 10 5 , then appears to saturate on a scale commensurate with the diameter of ∼ 4 cells.

Figure 2 (
Figure2shows the moving average of the lengthscale of the density autocorrelation ρ ζ , taken over 10 4 timesteps.Recalling that each cell has a nominal diameter ∼16 lattice units, the length scale of the density autocorrelation is initially roughly one cell, consistent with a well-mixed monolayer, while the length scale at long times increases to slightly more than 4 cells, consistent with microphase separation into clusters roughly 4 cells wide.These results demonstrate quantitatively that the mixture of extensile and contractile cells sorts according to activity.The data suggest that the cluster size saturates, but we cannot rule out a continued extremely slow growth of the clusters.Figure2(a) is reproduced, alongside plots of C ζ (r) at t = 1×10 3 and t = 1×10 5 , in Figure SI1, while Figure SI2 provides a comparison of SI and C(r) for a fully sorted system.The segregation index and density autocorrelation indicate microphase separation but yield no information about the dynamics of the system.We next examine the mean-square displacement of the phase-field cells, comparing the fully phase separated state (Fig.1d) to the microphase-separated state (Fig.1b).Results for the phase-separated configuration shown in Fig.3aindicate that, when surrounded by cells of the same type, the motion of the contractile cells is almost entirely arrested, while extensile cells move more freely.This is true also in pure monolayers: FiguresSI3 and SI4show, respectively, snapshots of layers