Large-scale simulations of biological cell sorting driven by differential adhesion follow diffusion-limited domain coalescence regime

Cell sorting, whereby a heterogeneous cell mixture segregates and forms distinct homogeneous tissues, is one of the main collective cell behaviors at work during development. Although differences in interfacial energies are recognized to be a possible driving source for cell sorting, no clear consensus has emerged on the kinetic law of cell sorting driven by differential adhesion. Using a modified Cellular Potts Model algorithm that allows for efficient simulations while preserving the connectivity of cells, we numerically explore cell-sorting dynamics over very large scales in space and time. For a binary mixture of cells surrounded by a medium, increase of domain size follows a power-law with exponent n = 1/4 independently of the mixture ratio, revealing that the kinetics is dominated by the diffusion and coalescence of rounded domains. We compare these results with recent numerical studies on cell sorting, and discuss the importance of algorithmic differences as well as boundary conditions on the observed scaling.


Introduction
Collective cell behaviors are involved in many morphogenetic events, such as embryo development, organ regeneration, wound healing, and the progression of metastatic cancer [1]. One of the simplest and best studied examples of such behavior is the spontaneous separation of two randomly mixed cell populations, in a process called cell sorting [2]. In the 1960s, Steinberg proposed the differential adhesion hypothesis (DAH) as a mechanism to explain cell sorting [3]. The DAH postulates that the surface tension of tissue arises from cell adhesion, and that cell configurations are determined by minimization of the surface energy of cell aggregates in a manner similar to phase separation processes, although surface energy minimization is an active process for cell sorting: cellular activity is necessary to untrap from the many local minima of the rough energy landscape of a cellular system. Since Steinberg's hypothesis statement, several numerical studies have proven that the DAH could reproduce cell sorting phenomena similar to those observed experimentally [4][5][6][7][8]. However, no clear consensus has emerged from these studies on the kinetics law of cell sorting driven by differential adhesion: in their seminal papers [4,5], Graner and Glazier found-on a somewhat smaller system-that the boundary length, defined as the number of mismatched cell contacts, decays logarithmically. More recently, Nakajima et al. [9] reported a power-law decay, with an exponent which depends on the proportion of the two cell types: from n = 1/3 for even mixture to n = 1/4 for uneven mixture. Cochet-Escartin et al. [8] also reported power-law decays, but with significantly higher exponents: n 2 [0.49, 0.74] from experiments, while n 2 [0.55, 0.59] from simulations.
Meanwhile, alternative mechanisms have been proposed in recent years to explain cell sorting phenomena-such as differences in cell motility that are either intrinsic [10] or dependent on the cells' local environment [11,12], or a combination of locally coherent motion with differential adhesion [13]-leading to different scaling behaviours. Although differential adhesion is sufficient to explain cell segregation, several recent studies show that the tendency of cells to actively follow their neighbors reduces segregation time scales [11,13,14].
The kinetics of cell sorting is crucial for morphogenetic processes, which must proceed with well-defined order and timing. Most importantly, characterizing the kinetics can help us to discriminate between the wide variety of proposed models, and thus would give us a deeper understanding of the underlying mechanisms that induce cell sorting. In this context, more extensive simulations of DAH-based cell sorting are needed to find out which type of kinetics law actually governs cell sorting. In this paper, we present extensive DAH-based cell sorting simulations for two-dimensional binary mixtures of cells surrounded by medium, using a recently modified Cellular Potts Model (CPM) algorithm that allows for simulations on very large scales in space and time. On long times, reported kinetics obeys algebraic scaling law with exponent n = 1/4 independently of the mixture ratio. We discuss the importance of algorithmic differences, boundary conditions, initial cluster geometry, space dimension, and finite size effects on the observed scaling.

Cellular Potts model
Cellular Potts model (CPM), also called the Glazier-Graner-Hogeweg (GGH) model, is one of the most accepted models of a multicellular system. It was initially created to simulate the behavior of single cells during cell sorting driven by differential adhesion [4,5]. Since then, this model has been extended to reproduce many other collective behaviors of cells [15][16][17][18]. Cellular activity is modeled with an effective temperature T, and the update algorithm is based on Metropolis-like dynamics [6]. The main justification for it is that average time evolution of a cell configuration then obeys the overdamped force-velocity relation F = μv, where v is the cell velocity, μ is an effective mobility, and F ¼ À rH is the acting force [19]. The CPM is a lattice based modeling technique: each cell is a subset of lattice sites sharing the same cell ID σ. A cell type τ(σ) is also defined for each cell. In our simulations cells are of two types, noted τ = B and τ = Y, where B and Y stand for Blue and Yellow, the respective colors of the cells in our graphical outputs (see Fig 1). To simulate a fluid medium that surrounds the cell aggregate, a third cell type τ = M is introduced (where M stands for Medium). Following [4,5], the CPM Hamiltonian H we use reads: The first sum in Eq (1) is carried over neighbouring sites hk, li and represents the boundary energy: each pair of neighbours having unmatching indices determines a boundary and contributes to the boundary energy. Here, σ k and σ l are the site values of site k and l, respectively. τ and τ 0 are abbreviations for τ(σ k ) and τ(σ l ). J τ,τ 0 (= J τ 0 , τ ) is the energy per unit contact length between cell types τ and τ 0 . The second sum in Eq (1) represents the compressive energy of the cells. B is the effective bulk modulus of a cell, which captures its out-of-plane elongation [20,21], A i is the area of cell i, and A 0 the nominal area.
Eq 1 is the minimal Hamiltonian whose minimization induces cell sorting phenomena. However, various additional terms can be added to account for real cell behaviors, or application of external fields [15,[22][23][24]. Especially, it is common to add a perimeter penalty term to Eq 1 to include the effects of cortical tension [9,19,23]. As this term does not depend on how the two cell types are disposed within the aggregate, it has no perceptible effect on the kinetics of cell sorting, which only requires that the adhesive energy of two heterotypic cells is less than that of two homotypic cells [4,5]. Let us emphasize that, as we investigate kinetics of cell sorting driven by differential adhesion exclusively, we just use Eq 1 in our simulations; we do not add extra terms that would account for spatial or temporal correlations in cell motion [10][11][12][13].
CPM presents a number of advantages over competing numerical models for multicellular systems, such as vertex-or voronoi-based models [25]: interfaces can have arbitrary shapes, T1 topological rearrangements or neighbor exchanges, which are crucial for dynamics in such systems, are naturally included, and free boundaries with the medium are handled as simply as boundaries between adjacent cells. However, a notorious weakness of CPM-when used with the standard updating rule-is that it does not guarantee connectivity of cells: as temperature increases, small fragments detach from cells, biasing the cell perimeters and cell centers of mass. Ideally, CPM cell sorting simulations should be run in a temperature range such that interface fluctuations are high enough for generating T1 neighbors swapping, but low enough for avoiding cell fragmentation. However, such a temperature range is almost nonexistent: the thermal energy required to generate a T1 in a hexagonal pattern (honeycomb) is estimated to be a fraction of the typical energy of one side: DE T 1 � 0:1gL [26], where γ is the boundary energy per unit length. In comparison, the minimal energy to create a fragment with size l is *γl [27]. Then, fragments of size up to l ' 0.1L show up before first T1 rearrangements occur. To prevent the visual appearance of fragments, quenching process is usually performed before snapshots. However, the creation of fragments, inefficient for the cell sorting process, artificially slows down simulations [27].

Modified algorithm that preserves cell connectivity
In a recent work [27], we proposed a modified CPM algorithm that forbids fragmentation and improves substantially computational efficiency for a same running temperature: the time wasted in checking cell connectivity is more than offset by the time spared in fragments moves, evidencing that fragmentation occurs even at low temperature. Moreover, efficiency can be further enhanced by working at temperature range that are unattainable with the standard algorithm, due to the proliferation of fragments. Like the standard CPM algorithm, this modified version follows a modified Metropolis updating rule: the value of a lattice site, chosen randomly (and called candidate site), is changed to a value chosen randomly in its neighborhood (target value)) with probability p = min(1, e −ΔE/T ) (where ΔE is the change of energy associated with the modification), and under the extra-requirement that this change preserves the local connectivity of the target and candidate cells. That is, all the sites with target (resp. candidate) value within the Moore neighborhood of the central site are adjacent to each other. Starting with an initial configuration where cells form simply connected domains, this test (which is performed at every copy attempt) ensures that connectivity of cells is preserved, i.e. it prevents the formation of fragments or multiply connected cells. Because our algorithm also belongs to the Metropolis-like class, we expect to observe the same kinetics than with the standard algorithm. Indeed, assuming that the frequency of fragmentation events that occur with the standard algorithm is homogeneous in time, the real time-correspondence of one MCS for both algorithms must be proportional to each other, for a given set of the different parameter values.
In [27], we performed a few DAH-based cell sorting simulations on small systems to benchmark this modified algorithm and compare its efficiency with that of the standard CPM algorithm. Here we use this efficient algorithm to fully characterize DAH-based cell sorting kinetics on unprecedentedly large systems. We performed simulations for binary mixtures of N cells, with N ranging from 200 to 320,000. All data were averaged over five independent runs. Moreover, in order to smooth out completely the anisotropy of the underlying square lattice, we used a 4th order neighborhood (composed of 20 lattice sites) in the evaluation of the Hamiltonian, whereas previous studies used a 2nd order neighborhood, composed of 8 lattice sites only [4,5,9]. Consistently, we also used a larger cell size (the target area was set to

Results
Fig 1 shows different snapshots of the simulation for an even (50:50) mixture of cells. Initially, the two types of cells were distributed randomly in a rounded cluster (Fig 1A). Subsequently, cells of the same type began to aggregate and form clusters (Fig 1B), and the sizes of the clusters grew with time ( Fig 1C).
Choice of running temperature has a strong influence on the cell sorting timescale: fluctuations must be large enough to jump from a local minimum of energy to a lower one, but not too large to avoid jumping to a higher one, so there is likely an optimal temperature for which cell sorting timescale is minimal. Most importantly, there is a full temperature range for which the kinetics of cell sorting is preserved. Indeed, the net rate of transition from configuration i to configuration j obtained with Metropolis-like algorithm is r ¼ 1 À e DH=T ' À DH=T for small relative change of energy DH=T. In this regime, a change in temperature uniformly scales the rate of transition between configurations, leaving the kinetics of cell sorting unchanged. Note that with the parameter values we used, DH=T ' ðJ=TÞ ' 0:2.
Following previous studies [4,5,28], we quantify the kinetics of cell sorting by tracking the evolution of the total boundary length Γ(t), defined as the number of mismatched cell contacts. Fig 2 shows the evolution of Γ(t) for even (50:50) mixtures of various sizes. We see that before complete cell sorting, curves for N � 2 × 10 4 cells follow the same algebraic variation: Γ(t) * t −n with n ' 1/4, attesting that this exponent is not affected by finite-size effects. Unexpectedly, our exponent value is different from the one reported by Nakajima and Ishihara [9] for even mixture ratio (n = 1/3), but is identical to the value they obtain for uneven mixture ratios.
We then study how the mixture ratio affects the power-law exponent in our simulations. Fig 3 shows the typical time evolution of cell sorting for a 20:80 mixture. As for the 50:50 mixture, cells of the same type aggregate into clusters whose size grows with time. We note however that at same simulation time, clusters are much smaller and rounder than for the 50:50 mixture. Fig 4 compares the time evolution of boundary length for a 50:50 and a 20:80 mixture ratio. Both curves rapidly follow the same power law with exponent n = 1/4. We also compare the evolution of the typical cluster size R(t) for the two mixture ratios. Following [29], we define R(t) as the position of the first zero of the autocorrelation function Cðr; tÞ ¼ hsðr 0 ; tÞsðr 0 þ r; tÞi r 0 À hsðr 0 ; tÞi where σ(r 0 , t) is the value at time t of the lattice site located at r 0 , and h. . .ir 0 denotes the average over all lattice site locations r 0 . Typical autocorrelation functions are shown in the inset of Fig 5. In agreement with [9], when they are plotted with the rescaled distance r/R(t), they all superimpose on a master curve, which shows that the system rescaled by R(t) exhibits the same statistical property. Note that boundary length and cluster size are related for circular clusters: is the number of clusters, whereas mass conservation implies N * R −2 , yielding Γ * R −1 . The evolution of R(t) for the two mixture ratios is shown in the inset of   this scaling regime is reached much more rapidly. For the even mixture, the curve follows a quite long transitory regime R(t)*t 1/3 , indicating that clusters are not rounded yet.
Our results then suggest a universal scaling n = m = 1/4 on long timescales, regardless of the mixture ratio, in contradiction with results of Nakajima and Ishihara [9]. Our exponent values also differ from those reported more recently by Cochet-Escartin et al. [8]. We discuss in the next section the origins of these discrepancies.

Discussion
There are two major mechanisms for phase demixing and domain growth. The first one consists in the "evaporation" then "condensation" of the units (atoms, or here, cells) from the smaller clusters to the larger ones. In this evaporation-and-condensation mechanism, centers of mass of clusters are fixed. The second mechanism is the diffusion then coalescence of clusters. In this diffusion-and-coalescence mechanism, centers of mass of clusters are free to move. In real systems, combination of these two idealized mechanisms can occur. For each of these two mechanisms, kinetics result from two competing processes: the growth rate of the clusters on one side, and the rounding of their interfaces on the other side. The well known Ostwald ripening corresponds to evaporation-and-condensation mechanism when kinetics is limited by the transfer of material from smaller to larger clusters, which then have a rounded shape. Lifshitz, Sliozov and Wagner (LSW) theory then predicts that the domain growth follows a power law *t m with exponent m = 1/3 independently of space dimension [30]. Although LSW theory assumes that the volume fraction of one of the components of the mixture is vanishingly small, it has been reported to hold for larger volume fraction as well [30,31].
As for the diffusion-and-coalescence mechanism, it is well-established [14,32,33] that in the diffusion-limited regime-in which clusters have rounded shape-domain growth also follows a power-law with exponent m = 1/(2 − dα), where d is the space dimension and α the scaling factor between diffusing velocity and mass of aggregates. For cells with aligning interactions, the value of α is known to depend on the alignment strength [14]. For cells with no aligning interactions, as it is the case in our DAH simulations, α = −1 [9,14], yielding m = 1/4 at two dimensions.
These results suggest that in our simulations cell sorting is achieved through diffusion-andcoalescence mechanism and that on long times the kinetics is limited by the diffusion of (rounded) clusters, for both even and uneven mixture ratios. Actually, this should come as no surprise, since cluster diffusion has a slower kinetics than cluster rounding, which generally is characterized by an exponential rather than algebraic law [7,8].
At this point one question is left: how to explain the n = 1/3 power-law reported by Nakajima and Ishihara for even mixture, while they did observe the n = 1/4 regime for uneven mixture? Although a n = 1/3 power-law is characteristic of the LSW theory, this mechanism has been ruled out by Nakajima and Ishihara to explain the observed scaling, as the probability of detachment of two homotypic cells is vanishingly small with the parameter values they used. Note that with the parameter values we used, the probability of detachment is p cell ¼ e À DE cell =T with ΔE cell * (2J YB − J YY − J BB )zL * 7.5, yielding p cell * 5.5 × 10 −4 . Hence, such events are unlikely in our simulations as well. Besides, we also checked that detached cells are never observed in the simulation snapshots.
Rather than the detachment of homotypic cells, evaporation-and-condensation could proceed through the detachment of fragments of cells. Obviously, such events do not occur with our non-fragmenting algorithm, but they are unavoidable with the standard CPM algorithm [27]. Such a mechanism is particularly favored in even mixtures, in which the distance between homotypic clusters is minimized. That may explain why Nakajima and Ishihara observe the n = 1/3 power-law for even mixtures, but not for uneven ones.
Apart from the use of a non-fragmenting algorithm, the other difference between both studies resides in the boundary conditions which are used: Nakajima and Ishihara used periodic boundary conditions, while in the present study, free boundary conditions are used, meaning that cells are surrounded by the medium. This difference in boundary conditions, combined with finite size effects, could also explain why Nakajima and Ishihara do not observe the n = 1/4 power-law for even mixtures. Indeed, as this is illustrated in Fig 1B, there is a long transient regime in our simulations during which clusters form entangled structures in even mixtures (a point which is also discussed in [34]): because of the high concentration of clusters, cluster growth process is faster than cluster rounding up process. The m = 1/3 power-law variation of the cluster size R(t) we observe (see inset of Fig 4) corresponds to this transient regime where clusters grow in size but are not circular. As these clusters keep growing in size, their diffusion slows down. Eventually, on long times, rounding is faster and the kinetics follows the n = 1/4 power-law, typical of the regime of diffusion and coalescence of circular clusters. When simulations are performed with periodic boundary conditions, these entangled structures rapidly span the simulation lattice and close on themselves via the periodic boundary conditions, as illustrated in Fig 6. Such spanning clusters then stop diffusing immediately (as their displacements require a global move of the whole system), and slow down interface smoothing, because of a competition between spanning clusters that tend to form stripes with flat boundaries, and smaller clusters that tend to form circular boundaries. Nakajima and Ishihara attribute the observed n = 1/3 power law to this interface smoothing process [9]. However, on long times, the n = 1/4 diffusion process, which is slower than the smoothing process, should still prevail, unless there are no small diffusive clusters left because of the finite size of the simulated system. But then the mean cluster radius R(t) should not evolve much as clusters stop growing, which is not what is reported in [9].
To get more insight on the exact cause of the discrepancy in power-law exponent, and discriminate between difference in algorithm and difference in boundary conditions, we performed cell sorting simulations of even mixtures with our modified CPM algorithm and using periodic boundary conditions. Contrary to Nakajima and Ishihara, we obtain a n = 1/4 power-law for the evolution of boundary length (see Fig 7), exactly as we obtained with the free boundary conditions. This result suggests that evaporation-and-condensation of cell fragments is at the origin of the n = 1/3 power-law reported in [9]. Accordingly, the long-term scaling law is independent of boundary conditions and mixture ratio when using the nonfragmenting algorithm.
Initial geometry of the aggregate and dimension of space also affect the kinetics of cell sorting, as pointed out by Cochet-Escartin et al [8]. In their study, cells were initially placed on a sheet in a 3D space. The exponents they observed in this configuration, either experimentally or numerically, were significantly higher than 1/4. However, the initial flat configuration cannot affect the long time scaling of domain growth, and the evolution in a 3D space should rather decrease the exponent (passing, in the diffusion-limited regime, from n = 1/(2 − dα) = 1/4 in 2D to n = 1/5 in 3D). Moreover, the exponents for Γ(t) and R(t) were quite different in their experimental data, suggesting that clusters did not have time to round up. Finite size effects, but also hydrodynamic couplings (for experiments) [34,35], could explain the larger exponent values they obtained.
The good agreement between the scaling law reported in the present study and in [9] for uneven binary mixtures confirms that the real time-correspondence of 1 Monte Carlo Step (MCS) for the standard CPM algorithm and ours, which prevents fragmentation, are proportional to each other t 0 1MCS ¼ at 1MCS , where the prefactor α depends only on the parameters that characterize the system (here J ij , B and A 0 ). This linear relationship can be easily explained: as the frequency of fragmentation events that occur with the standard algorithm is homogeneous in time, the amount of time spent in the creation of fragments is homogeneous all along the sorting process.

Conclusion
Using a modified CPM algorithm that forbids cell fragmentation and increases computational efficiency, we were able to assess the kinetics of cell sorting driven by differential adhesion on very large systems (up to 320,000 cells). Our results show that on long times, kinetics is controlled by diffusion of rounded clusters binary mixture for cell aggregate surrounded by medium. The growth exponent is n = 1/4 for two-dimensional systems, regardless of the mixture ratios of the two cell types. Two explanations have been proposed to rationalize the apparent contradiction with exponent value reported previously: the use of a modified CPM algorithm which forbids cell fragmentation in one hand, and the use of different boundary conditions on the other hand. However, when using our non-fragmenting algorithm with periodic boundary conditions, we obtain the same n = 1/4 power-law as with free boundary conditions, suggesting that the change in algorithm is at the origin of the difference in results. In order to confirm this result, further investigations must be conducted, with varying mixture ratio and simulation temperature. If this is verified, previous studies using the standard algorithm for simulating cell kinetics might have to be reconsidered, as they could be affected by the same nonrealistic cell fragmentation issue.
As a concluding remark, we note that both 1/3 and 1/4 power laws are slower kinetics than those reported in recent experimental observations [8,36,37]. This, however, does not necessarily reject the differential adhesion hypothesis. Indeed, these exponents are obtained for ideal planar geometries, and in a steady regime, requiring systems with a large number of cells. These two assumptions are not easily fulfilled in experiments, and change of cluster geometry or finite size effects can alter the observed kinetics. These limitations will be taken into account in future research work. Other measures at the cell level, such as velocity distribution in the two cell populations or correlations in cell velocity can help to discriminate the different models. Using such complementary measures [8] to support the DAH over the differential motility model as the mechanism that drives regeneration of Hydra. However, other experimental studies support collective cell migration as the main mechanism of cell segregation [37]. It is therefore likely that cell sorting processes, which occur in a large variety of systems, cannot be explained by a single universal mechanism. Given that DAH-based cell sorting has slower kinetics, it must be involved in small multicellular systems primarily, while other mechanisms such as collective motility or chemotaxis might occur in larger systems.