Kinetics of the cellular Potts model revisited

The relaxation kinetics of cell sorting are studied with the cellular Potts model. In contrast to previous reports, the increase in domain size is found to obey a power law (R(t)∝tn). The growth exponent turns out to be n=1/3 for an even mixture of two cell types, where the domains for each cell type are interconnected and the kinetics are dominated by smoothing of the domain boundary. The exponent is n=1/4 for uneven mixtures where cell sorting proceeds via the diffusion–coalescence of circular cell domains. The exponent is explained by the decrease in motility of a cell cluster as a function of its size according to D(R)∝R−2. Our results provide a theoretical framework for elucidating how cell populations migrate within tissue.


Introduction
How does collective cellular behavior arise from the physical properties of the constituent cells? Different types of animal cells, when dissociated from tissues and intermixed, sort themselves into cell domains of the same type [1]. Such behavior, known as cell sorting, is involved in various morphogenetic events including the regeneration of hydra and the formation of tissue boundaries in vertebrates [2]. The simplicity of this experimental system has resulted in valuable insights into the physical processes underlying collective cell behavior. In the 1960s, Steinberg proposed the differential adhesion hypothesis (DAH) 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. Indeed, as expected from the DAH, cell configurations after sorting reflect adhesive strengths in hydra cells [4]. It is also observed that cells with higher surface tension are enveloped by others with lower surface tension in sorting experiments with chicks [5] and zebrafish [6]. Several theoretical models based on the DAH have succeeded in reproducing the eventual cell configurations in sorting experiments [7]- [10]. Among them, the cellular Potts model (CPM), also called the Glazier-Graner-Hogeweg (GGH) model, is one of the most accepted models of a multicellular system [7,8]. This model not only captures cell sorting phenomena, but also explains patterns of various morphogenetic processes among a wide range of taxa such as: cell configuration in the retina of the fruitfly, Drosophila; vasculogenesis in vertebrates; and the development of slime molds [11]- [13].
In contrast to its stationary patterns, the kinetics of the segregation process have been less studied so far. Given that morphogenetic processes should proceed with well-defined order and timing, the kinetics have a large impact on the temporal sequence of developmental processes. Furthermore, a characterization of the kinetics would give us a deeper understanding of the underlying collective dynamics of cell populations. Experimentally, it has been reported that the characteristic size of cell domains increases linearly with time in a mixture of chick epithelial and retinal cells, much like immiscible fluids [14]. Observations of hydra cells also indicate the fluid-like behavior of cell aggregates [15]. In contrast to these experimental results, it has been suggested that the CPM displays logarithmic growth for domain size [7,16], which is much slower than the observations. It should be noted, however, that the sizes of the simulated systems were too small (the number of cells was of the order of 10 3 ) to determine the growth law. Recently, Belmonte et al [17] adopted a self-propelled particle model and studied the relaxation process of cell sorting. They showed that the growth law is not logarithmic but algebraic with exponent n 0.18, and argued that the transition from the logarithmic to the algebraic law is the result of directional and collective cellular movement introduced in the model. However, the obtained exponent is still much lower than that of the experimental results. Moreover, the growth is even slower than that of model B (the Cahn-Hilliard equation), a canonical model of phase separation of binary substrates [18,19], in which the dynamics are driven by diffusion and the domain size increases algebraically with exponent n = 1/3 [20]. To resolve the discrepancies among these results, in the present paper we have closely re-examined the kinetics of the CPM in two dimensions with a system of sufficiently large size.

Results
The CPM was introduced as a modified Potts model for reproducing the biophysical properties of cells such as excluded volumes and deformations. In the model, M cells are indexed by the 3 Potts spin σ ∈ {1, . . . , M}, and the ith site is labeled by the index of the occupying cell as σ i . Thus, each cell is represented as a set of sites on a lattice. The type of the σ th cell is denoted by τ (σ ). The dynamics of the CPM are based on the process of minimizing the effective energy of the model, which is given as follows [13,21]: Here σ , σ , τ and τ are abbreviations for σ i , σ i , τ (σ i ) and τ (σ i ), respectively. V σ and L σ are the volume and surface area (peripheral length in two dimensions) of the σ th cell, respectively. The first term represents the adhesive interactions among neighboring cells (the Moore neighborhood is considered here). J τ τ indicates the surface tension per unit length between the ith and i th sites, which depends on the interacting cell types. The second and third terms represent the energy costs for a cell to deviate from the preferred volume V 0 and the length L 0 , respectively. The total number of sites in the system is V 0 M. The time evolution of the model is obtained by Monte Carlo simulations with the Metropolis algorithm. First, the cell index of a randomly chosen site is substituted with that of a neighboring site as a trial. Subsequently, the energy difference before and after the trial, H ≡ H after − H before , is calculated. The new configuration is accepted with probability p = min {1, exp(−β H)}. β is the inverse of the effective temperature of the system, which is related to the random ruffling activity of cell membranes in the CPM where a high (low) temperature corresponds to active (less active) cell motion. In the present study, we considered the sorting process of two cell types (τ = ±1) in a two-dimensional (2D) square lattice. The parameters were set as follows: Subsequently, cells of the same type began to aggregate and form clusters (figures 1(b) and (e)), and the sizes of the clusters grew with time (figures 1(c) and (f)). The domains for each cell type continued to increase in size until they grew to a size comparable to the system size.
To quantify the kinetics of this growth process, we calculated two physical quantities during the sorting process. One is the domain boundary length γ (t), defined as [7] γ (t) = 1 √ mixture ratio number of heterogeneous bonds total number of bonds where M − is the number of minority cells and M is the total number of cells. The prefactor was added for normalization and the summation was taken over pairs of neighboring sites. Another index is the spatial correlation function of cell types, Here, · denotes the average of all sites. The typical domain size R(t) was defined as the first zero-crossing point of C(r, t) [20].  [7,16]. However, on a longer time scale, the temporal behavior did not show logarithmic decay. Instead, γ (t) decreased algebraically according to the relation γ (t) ∼ t −n with n = 0.313 ± 0.007 (figure 2(a)). The time dependence of R(t) for the even mixture also showed power-law growth with n = 0.334 ± 0.006 instead of logarithmic growth ( figure 2(b)). We also checked the dynamical scaling of the CPM, given that domain growth in phase separation with power-law kinetics is supposed to obey dynamical scaling [19]. Dynamical scaling implies that the correlation function satisfies C(r, t) ≈ F(r/R(t)) with a unique function F(x); thus, the system rescaled by R(t) exhibits the same statistical property. In figure 2(c), C(r, t) is plotted as a function of r/R(t). The data for different time points turned out to lie on the same curve. Therefore, scaling was found to hold for the sorting process of the CPM. Taken together, cell sorting for the even mixture obeyed power-law growth and proceeded via smoothing of domain boundaries as in the case of model B dynamics, where the diffusive transport of molecules (or cells) is driven by surface tension [22]. The obtained exponents suggest that the CPM belongs to the same class of models as model B whose exponent is n = 1/3.
Uneven mixture ratio. We next examined cases in which the mixture ratio of the two cell types is not 50 : 50. Here again, it was found that γ (t) and R(t) exhibited power-law (rather than logarithmic) dynamics in the late stages (figures 2(a) and (b)), although their growth  exponents were different from those of the even mixture. We evaluated the growth exponents of γ (t) and R(t) from least-square fits to the data, as summarized in table 1. Whereas the exponents for the even mixture were close to 1/3, those for uneven mixtures were in the range of 0.24 < n < 0.26. To ascertain the correct growth exponents from the asymptotic behavior of R(t), we calculated the time-dependent effective exponent n eff (t), which was defined as  with c = 5. The broken lines indicate n eff = 1/4 and 1/3. The data are within the time interval 0 < MCS < 4 × 10 6 . For M = 128 2 and 256 2 , the data were averaged over four independent runs. For M = 512 2 , the data were averaged over six runs in 0 < MCS < 1 × 10 6 and over two runs in 1 × 10 6 < MCS < 4 × 10 6 . Error bars indicate SD of independent trials. [20]. Figure 3 shows n eff as a function of 1/R(t) for the mixture ratios 50 : 50 and 20 : 80. For the even mixture, the effective exponent increased towards 1/3 and stayed along near 1/3 for large values of R(t) (1/R(t) → 0). For the 20 : 80 mixture, the effective exponent was about 1/4 for large values of R(t). For both even and uneven mixture ratios, n eff (t) increasingly coincided with the expected values as the size of the system increased. As discussed below, the difference between the exponent of the uneven mixture and that of the even mixture came from the distinct morphologies of cell domains; for the even mixture, the domains of both cell types showed interconnected structures (figures 1(a)-(c)), whereas for the uneven mixture, those of minority cells exhibited circular shapes (gray cells in figures 1(d)-(f)). Note that for a nearly even mixture, the system showed cross-over of the exponent from 1/3 (10 4 < MCS < 10 5 for 40 : 60 in figure 2(a)) to 1/4 (10 5 < MCS). This crossover coincided with morphological changes in the cell domains (data not shown).
The potential mechanisms of domain growth for uneven mixtures include evaporationcondensation [23] and direct domain coalescence (diffusion-and-coalescence) processes [24]. The evaporation-condensation process is a continuous transport of molecules (or cells) from small to large domains, driven by the gradient of chemical potential among domains. This is the dominant mechanism for the domain growth of uneven mixtures in model B [22]. In the present case, however, this mechanism apparently was not engaged because cells were barely separated from a domain to which they belonged. To understand this, consider a situation in which two neighboring cells of the same type are surrounded by cells of another type. To separate the neighboring cells, it is necessary to overcome an energy barrier to break the adhesion between cells. Given that the energy cost to break a single site homotypic bond is , the cost to separate two cells of the same type is H cell ∼ l H site , where l is the number of adhesion contacts between the cells. With the parameters we used ( H site = 0.2, β = 2.5 and l ∼ L 0 ), the probability that the two neighboring cells will separate  was estimated as p cell = exp(−β H cell ) ∼ exp(−10) = 4.5 × 10 −5 . This is an event that rarely occurs because the average wait time (τ ∼ p −1 cell ) is comparable to the time required for the entire simulation. Thus, the evaporation-condensation process is not dominant in the sorting process of the CPM. Instead, domain growth proceeds via diffusion and coalescence of the cell domains. In this mechanism, the rate-limiting process is the diffusion of the domains, where the time dependence should be written as | r| 2 ∼ Dt. Suppose that the diffusion constant of a cell cluster, D, depends on its size according to the relation, D(R) ∼ R −γ ; then the time dependence of the domain size is written as R(t) ∼ t 1/(2+γ ) , as indicated by dimensional analysis.
To test this scenario, we calculated D(R) of the CPM by numerical simulations. We evaluated the diffusion constant by measuring the displacement of a cell cluster over time. Figure 4 shows the diffusion constant D as a function of cluster size M 0 (the number of cells in the cluster). Here, the diffusion constants were measured under two conditions to test the effects of the surrounding cells: a cluster of cells surrounded either by the medium or by cells of the other type. For the former case, the cluster can diffuse freely in space because the energy constraint of the second and third terms of equation (1) where d is a spatial dimension and D 0 is the noise intensity of ξ .

Discussion
In summary, we investigated the kinetics of cell sorting by adopting a commonly used multicellular model, the CPM, in which the cell movement is diffusive and driven by mechanical interactions between cells. We found that there exist two kinds of kinetics depending on the mixture ratio. The domain growth for the even mixture proceeded by smoothing the interface of interconnected domains driven by surface tension. This process obeyed power-law growth with exponent n = 1/3. For uneven mixtures, the growth proceeded via the diffusionand-coalescence process of circular cell domains, which obeyed power-law growth with exponent n = 1/4 in two dimensions. These results indicate that cell sorting in the CPM proceeds much faster than was previously believed [7,16].
We should note that the sorting process in chick cells is faster than that of the CPM results obtained here [14]. This difference implies that the migration of chick cells is driven, for instance, by active dynamics such as chemotaxis and/or coordinated movement among neighboring cells. Directed and collective movements of cells were addressed in the sorting processes of hydra [15,25] and zebrafish cells [6]. Thus, we might need to assume active effects to model these systems. In addition to the systems above, the kinetics of cell migration are crucial ingredients of metastasis and many other developmental processes. Recent imaging techniques allow researchers to track cells in a population and measure their migration dynamics quantitatively [26]. The present results on the sorting kinetics and diffusion of cell clusters will provide a reference framework for interpreting such experimental observations in the future.