Neighbor-enhanced diffusivity in dense, cohesive cell populations

The dispersal or mixing of cells within cellular tissue is a crucial property for diverse biological processes, ranging from morphogenesis, immune action, to tumor metastasis. With the phenomenon of ‘contact inhibition of locomotion,’ it is puzzling how cells achieve such processes within a densely packed cohesive population. Here we demonstrate that a proper degree of cell-cell adhesiveness can, intriguingly, enhance the super-diffusive nature of individual cells. We systematically characterize the migration trajectories of crawling MDA-MB-231 cell lines, while they are in several different clustering modes, including freely crawling singles, cohesive doublets of two cells, quadruplets, and confluent population on two-dimensional substrate. Following data analysis and computer simulation of a simple cellular Potts model, which faithfully recapitulated all key experimental observations such as enhanced diffusivity as well as periodic rotation of cell-doublets and cell-quadruplets with mixing events, we found that proper combination of active self-propelling force and cell-cell adhesion is sufficient for generating the observed phenomena. Additionally, we found that tuning parameters for these two factors covers a variety of different collective dynamic states.


Introduction
Biological cells are the fundamental building blocks of all lifeforms, ranging from single-cell animal-like amoeba to more complex multi-cellular organisms like us human beings, in which they form various organs and complex physical structures. Indeed, cells are a remarkable material with an amazingly wide range of versatility and flexibility. Some densely packed populations of cells are rigid enough to be viewed as a solid (for example, imagine a piece of hardwood or animal bone), while in some other cases cells prefer to stay alone and move erratically, like gas particles as in populations of Dictyostelium discoidium amobae in their vegetative state [1] or microglia (immune cells of the brain) which tend to avoid each other upon a physical contact in culture [2]. Cells within a solid-like population are caged by their immediate neighbors, maintaining their relative positions and orientations [3][4][5][6][7][8]. By far more common situation is that cells form soft tissues or tissue-like structures which are analogous to amorphous materials like glass, or perhaps, more accurately, something in between glass and liquid in which cells can diffuse, flow, and flock [3][4][5][6][7][8][9][10].
Ice melts to water, and water evaporates to vapor, as the temperature rises: In other words, the state of matter transforms as some system parameter(s) change. Likewise, the biophysical state of the cellular population can also undergo transitions [11], and the topic has drawn a great deal of attention in biology in association with morphogenesis and tissue-remodeling [12][13][14], wound repairing [15,16], and tumor growth [17,18]. Biological tissues are a collection of interacting 'active' cells in nonequilibrium states, therefore, in principle, cell-population can support numerous different non-epithelial dynamic states, to which an initially sedentary tissue state can switch. One of such states is flocking and a good example is the recent experimental study of Mitchel et al. [19], where air pressure was used as a control parameter for transforming an airway epithelial tissue to cooperatively migratory flocking cells. The authors referred to the transformation as an unjamming transition (UJT). The same study could also induce a different type of transition, which the authors termed as partial epithelial-to-mesenchymal transition (pEMT): Here, a chemical agent called TGF-β1 was used to transform an initial epithelial state into a hybrid state having a mixture of epithelial and mesenchymal characteristics. The authors have concluded that UJT and pEMT, respectively, yielded two rather different liquid-like states having very divergent dynamic, structural and molecular marker characteristics. The study of Mitchel et al. is an excellent example suggesting that for fully unfolding the dynamics of the dense cell population the relevant phase space needs to be at least two-dimensional (i.e., requires two independent parameters).
Until now, a great deal of attention has been paid mainly to the (phase space) parameter regime in which interesting large-scale collective waves and flocks could form [16,[20][21][22]. In the meanwhile, a more prevalent condition, especially for tumors, could involve the pEMT states near the border to sedentary epithelial states. Unlike a large-scale flock, whose spatial correlation length spans multiple cell-diameters, pEMT states may have a very small range of cell-cell interaction [3,23] with a spatial correlation length just about the size of a single-cell. So far, the spatiotemporal population dynamics within this regime is not well characterized at all except for the fact that cellular rearrangement can occur via cage-breaking with junctional change and rearrangements [5,24].
In this paper, we showed that cells in a densely packed two-dimensional cell layer in a hybrid (pEMT) regime could support some remarkable motile properties: They exhibited a higher diffusivity and longer persistent time in the moving direction than those of freely crawling single cells. Moreover, we demonstrated that the fascinating phenomena observed in small cell-clusters, namely 'cell-doublet rotation' and 'position swapping in cell-quadruplet', account for this dynamic property. The set of unusual observations was initially made with two-dimensional cell cultures of MDA-MB-231 breast cancer cell lines. Then, they were reproduced faithfully with a cellular Potts Model (CPM), in which we tuned two key parameters associated with cell-to-cell adhesiveness (or stickiness) and self-propulsion, in a systematic fashion. Importantly, we identified the relevant parameter regime, where the CPM simulations reproduced the core experimental observations in a self-consistent way. Interestingly, a significant amount of cell-cell adhesion was mandatory for the observed enhancement of persistence and diffusivity of cells. The phenomenon of neighbor-enhanced diffusivity (NED) of confluent population of MDA-MB-231 has a strong similarity to 293-MOCA expressing human kidney 293T cells, which become super-diffusive only in the presence of cell-cell adhesive interactions [25].

MDA-MB-231 cell-population vs. freely crawling cells
MDA-MB-231 cell-lines are a famous triple-negative breast cancer cell known for their aggressive, metastatic potential [26][27][28]. We chose MDA-MB-231 cell lines for our investigation since they have been popularly used for many different cancer studies. More importantly, recent studies [29,30] clearly suggested that they pertain to the aforementioned hybrid epithelial-mesenchymal properties. Monoclonal MDA-MB-231 cells were cultured following a typical laboratory culture protocol, harvested, and uniformly plated on a culture dish for imaging (see Methods for more details). A small central square area of the densely packed population is shown as a black and white image on the x-y plane of Fig 1A. The same figure also illustrates exemplary trajectories of five moving cells that were initially located adjacent to each other in the middle of the population. As they diffused out, they intermittently made sharp turns and fluctuations in their instantaneous speeds (0~7 μm/min), which are represented by the thickness of the colored 'tubes' in Fig 1A. Similarly, Fig 1B illustrates a superposition of trajectories of five randomly chosen, freely moving cells that experienced little or no cell-cell interaction. In both cases, the cells were super-diffusive, with the diffusivity of the cells in densely packed population being even stronger than that of the freely moving cells. This property is apparent in Fig 1C and 1D. Fig 1C shows a log-log plot of mean squared displacement (MSD) hδ 2 i versus time t, where solid red and blue dots trace the average profiles for freely crawling cells and cells in confluent population, respectively, while the background shadows represent the associated standard deviations. The upper panel of Fig 1D plots the diffusion exponent α for the marked time window in Fig 1C, which is blown up in lower Fig 1D. If we assume hδ 2 i = Dt α with D diffusion coefficient, α = 1 stands for a normal diffusion, α<1 (α>1) represents a subdiffusion (super-diffusion), and α = 2 stands for a ballistically moving object. As the plot of α in Fig 1D clearly illustrates, the value of α for the cells in confluent population was significantly larger than that of freely crawling cells for the entire range of time. In addition, the time range of super-diffusivity (α>1) for freely migrating cells terminated around t~10 2.4 (~251) min, quite earlier than the time, approximately 10 2.7 (~501) min, beyond which cells in confluent population lose their super-diffusivity.
The difference in the degree of diffusivity discussed in Fig 1C and 1D was also consistent with the analysis of auto-correlation function of unit tangent vector, c auto ¼ hêðt o Þ �êðt o þ tÞi (see Fig 1E), which measures the average degree of moving directional correlation between two unit-tangent vectors obtained at two different time instances which are separated by time t. The function c auto decayed significantly faster for freely migrating cells than for the cells in confluent populations (average decaying time constants of the ensemble: hτ 2,free i = 66 ± 104 min vs. hτ 2,pop i = 116 ± 150 min. In other words, cells in a confluent population had a stronger tendency of moving along (or a longer memory of) the direction to which they were heading than freely migrating cells. Subsequently, considering the mean migration speed of hv free i = 0.8 ± 0.3 μm/min and hv pop i = 0.6 ± 0.2 μm/min, the matching persistence lengths are, respectively, hl free i = 52 ± 81 μm and hl pop i = 72 ± 99 μm. According to [31], one may require α > 1 for the time greater than the directional persistence time τ 2 as a necessary condition for superdiffusivity. Since the transition time points to the normal diffusion were measured to be around 251 min (for freely crawling cells, indicated by the 1st arrow in Fig 1D upper panel) and around 501 min (for cells in confluency, indicated by the 2nd arrow), which are much larger than hτ 2,free i = 66 min and hτ 2,pop i = 116 min, respectively, it seems appropriate to use the term super-diffusivity.
The analyses given in Fig 1A-1D univocally pointed to the fact that the enhanced diffusivity was not at all due to a large-scale cooperative flocking movement. As a matter of fact, the range of the cell-to-cell interaction turned out to be rather short, spanning only about one celldiameter (~30 μm). For this estimation, we identified all individual cells and their centroids for a stack of a time-lapse movie (700 frames, 2 mins interval) and computed mean spatial velocity-velocity correlation map of c cross ¼ hêð r  Fig 1F, we may interpret that the cell motility in confluency is a random mixture of a few cells co-moving in line like a "two-cart vehicle" or rotating (i.e., moving past each other) together in a time-shared manner. Finally, Fig 1G depicts the overall shape of MDA-MB-231 cell in a confluent population. Simply, it is the average point density map of the centroids of all cells in a given population with respect to the reference cell, whose position was again brought to the origin (0,0) and whose instantaneous moving direction was aligned toward the positive y-axis as it was done for Fig 1F: The oval shape with its long axis along the y-axis is rather clear; this is because the migrating cells tended to elongate along the direction of movement. This oval shape hides the heterogeneities and time-fluctuations of individual cells' shapes within the population. In fact, despite the monoclonal nature, the shapes of MDA-MB-231 cell-lines in confluent culture were quite heterogeneous as clearly shown by the broad distribution function of shape index p (defined as perimeter= ffi ffi ffi ffi ffi ffi ffi ffi ffi area p ) in Fig 1H (edges of the cells in confluency were obtained by thresholding the phase-contrast images and identifying contours in the binary images).

Description of cell-doublets
A useful insight into the nature of NED within dense population could be obtained by analyzing the motile behavior of dispersed 'doublets.' Low-density cultures of MDA-MB-231 cells generally exhibited many dispersed pairs, as the cells divided following a cell-cycle (period3 2 hr) and preferred to stay together during migration (see S2 Fig for details). Doublets could also be formed by random encounters of individual cells. Accordingly, we tracked multiple doublets and analyzed the moving trajectories of the cells involved. Remarkably, the MDA-MB-231 doublets in our culture almost invariably exhibited a robust rotational movement as shown by the example in Fig 2A. This doublet rotated in a highly periodic fashion about 7 times in 900 minutes, during which one abrupt reverse-turn occurred between t = 240 min and 300 (marked by an arrow in Fig 2A; also see the top frame of Fig 2B).  cells maximally stretched. Thus, the doublet rotation must not be viewed as a simple rigid body rotation but more as a continuous position swapping between two contacting cells.
The observed MDA-MB-231 pair-rotation was a robust phenomenon (see Fig 2C and S3 Movie): 57 pairs (pooled from three different cultures) maintained their integrities for more than 212 min, which was the minimum time duration of tracked doublets. For some doublets, the rotation lasted as long as 1324 min beyond which they replicated or separated. Note that the slope of unwrapped angle θ vs. time t plot in Fig 2C varies significantly from one to the others. In fact, the heterogeneity in the pool of analyzed doublets could be quantified with a number of different measurements. Mean rotation period T for each doublet was defined to be T ¼ 2p=joj and its broad probability distribution function (PDF) is shown in the upper left frame of Fig 2D. Also, we could characterize the oscillatory behavior in d by calculating the peak position time τ d of the Fourier transform of d. The PDF of τ d is shown in the lower-left position of Fig 2D, which is also broadly distributed. We noted that the ensemble mean of T (= 305 min) was very close to twice as much as the mean of τ d (= 153 min) confirming that two stretches (and contractions) had occurred during one complete rotation. corr pair , which we defined as the temporal mean of ;ðtÞ ¼ê 1 ðtÞ �ê 2 ðtÞ, as well as the Pearson correlation coefficient r exhibited a broad PDF as shown in the lower (upper) right frame of Fig 2D, respectively.

Description of cell-quadruplets
The rotational behavior of MDA-MB-231 doublets was conserved even as they grew to a cellquadruplet, yet with additional complexity and intrigue that were incurred by having multiple neighbors (see Fig 3). Four cells forming a quadruplet cluster could sustain a steady 'unimpaired' rotation (which we refer to as a rotation with no change in cyclic positional order) but typically only for some short duration of time: See, for example, the dynamics during the time window t = 260~440 min which is highlighted as cyan shading in Fig 3A and 3B; and its matching sequence of snapshots is given in Fig 3C, top row. This counterclockwise unimpaired rotation of four (pseudo-color coded) cells during t = 260~440 min was robust with their phase ordering (red, yellow, blue, violet) maintained, but was not a rigid body rotation as their relative distances to the centroid of the quadruplet significantly varied (see green lines in Fig 3B). Typically, the unimpaired rotation was interdigitated by frequent position swapping events: A good exemplary time window is highlighted as red shading in Fig 3A and 3B, and its matching sequence of snapshots is given in Fig 3C, bottom row. During the time window of t = 60~240 min, two events of pairwise position swapping took place (first, between the purple and the yellow, and second, between the blue and the yellow), which are marked by two dashed black circles in Fig 3B (top). The swapping events were rather quick in the sense that the angular speed |ω| (black lines) peaked very sharply as shown in Fig 3B and the peaks' bandwidths were a mere fraction of the typical period (~640 min) of quadruplet rotation. In addition, there were concurrent changes in d but at a much slower time scale, suggesting that there were some significant cell-stretching and contraction dynamics involved during and beyond the abrupt change in θ associated with a swapping action. Thus, we may as well view this swapping action as a 'mixing' process within a rotating cell-quadruplet. Note that in the aftermath of the two swapping events, the initial clockwise unimpaired rotation (for example, during t = 0~70 min) of the cluster became counterclockwise (for t ≳ 240 min).
The discussed MDA-MB-231 quadruplet rotation phenomenon was universal (see S4 Movie), and once they formed a cluster, they rarely disintegrated. However, as it was the case for the rotating doublets, the analyzed quadruplets exhibited quite heterogeneous physical properties. Fig 3D shows temporal evolution of mean rotation angle (phase) hθi, which was  Fig 2C, all global reverse turns were unfolded to straighten out the function hθ(t)i, more or less to a line, so that the overall slope of each line represents the average angular speed of its matching quadruplet. Position swapping events created jitters in the line of hθ(t)i and caused phase lags; thus, they somewhat slowed down the overall rotation speed. Fig 3D well illustrates the wide range of slopes (i.e., angular speeds) that MDA-MB-231 quadruplets could support. Subsequently, the PDFs of mean rotation period T given in the upper-left frame of Fig 3E show a broad spectrum. Note that its ensemble mean value of 638 min was much larger than the ensemble mean of rotation period of the doublets, which was 305 min: This was of course natural since cells in a quadruplet must travel a longer distance than those in a doublet to complete a turn as the cluster (area) size had approximately doubled; furthermore, once again there were many phase-delaying position swapping events for the quadruplets. The other PDFs in Fig 3C all consistently show a broad spectrum. Interestingly, the mean value of τ d (183 ± 75 min) and that of r (-0.3 ± 0.2) both matched to those of doublets discussed earlier (153 ± 64 min and -0.3 ± 0.1, respectively) very closely. This could be an indication that the observed swapping events might not be a stochastic process but related to the innate nature of contacting doublets making rotational movement. Incidentally, we find the observed position-swapping events within MDA-MB-231 cell quadruplets almost exclusively took place not along 'diagonal direction' but along 'normal directions' [see Fig 3E (lower right frame) and the schematic illustrations shown in Fig 3F].

Simulated cell-population vs. freely crawling cells
Most of the features from our experimental observations could be faithfully recapitulated by a CPM of active cells which can generate self-propulsive movement [32][33][34]. As we will demonstrate, self-propulsion strength factor S and interfacial energy E were indeed key parameters governing the collective behavior of confluent population as well as the coordinated movement of cell-doublets (see Methods for more details about other parameter values used for the simulation). Before we present a complete two-dimensional phase-diagram spanned by S and E showing various modes of cell motility, we first discuss our specific simulation results that closely matched those that we observed in experiments. The particular choice of S and E used for the simulation was validated by a "self-consistency check" guided by our experimental data. In other words, almost all physical properties of rotating doublets, mixing cell-quadruplets, and non-flocking, super-diffusive migration of cell-population, were recapitulated successfully with a single set of S and E.    Fig 5B); robustness of the rotation (Fig 5C); and the fact that mean modulation period τ d of the cell-center displacement was approximately half the mean rotation period T (Fig 5D). The spread of the unwrapped θ lines shown in Fig 5C reflects   This discrepancy was mainly due to the fact that MDA-MB-231 cells experienced position swapping events more frequently than the CPM cells. Another factor for the discrepancy could be that the actual cells were softer and more ramified than the model cells; we hypothesize that such properties could result in more severe shape deformation that in turn slows down the coherently directed rotation. Finally, the bar graph in the lower right of Fig 6E confirms that the swapping between diagonal neighbors was almost non-existent as observed in the experiments. failing to make a turn (S5 Movie shows the process at S = 0, E = 0). As S became large (≳ 5) with E~0, doublets separated (S6 Movie at S = 6, E = 0). Doublet rotation was visible as E became more negative while S being non-negligible (S7 Movie at S = 5, E = -100). The threshold value of S for rotation tended to decrease slightly as E became more negative for the region S ≲ 5 (see S4B Fig). In the upper right corner of the phase diagram, more stretched rotation (S8 Movie at S = 10, E = -100) were visible.  measure of the diffusivity of the cells, and p was chosen since it has emerged in recent studies as an important structural signature for the migratory transition of a confluent tissue [5,7]. Clearly, S was critical for bringing a jammed population (S5 Movie) to an unjammed state for the entire range of E, while decreasing E facilitated the unjamming for the case of a small (~2) S (Figs 7C and S5A). For the migratory state around S~5 and E~0 (i.e., no adhesiveness), diffusion exponent α calculated at a long time-scale showed super-diffusivity with greater than 1 (S5A Fig), and concurrent rise in the average size of cohesively moving cell-clusters was visible ( S5B Fig and S6 Movie). As E decreased, cells migrated more diffusively and individually (S7 Movie), while interfaces became more ragged (Fig 7D). Finally, all four level curves shown in Fig 7A-7D were put together in Fig 7E, and they almost meet together at a common point (red dot, S = 2.8, E = -65), whose values of S and E were used for all CPM simulations given in Figs 4-6.

Confluent population phase diagram
Increasing the self-propulsion parameter S, with a fixed value of E = -65 corresponding to the red dot in Fig 7E, intermittent co-movement (flocking) of the doublets became more frequent (S6 Fig). As it happened, the instantaneous angular velocity decreased while linear velocity increased. The time-fraction of doublets' flocking as a function of S indicates that there is a continuous transition between the rotation and flocking. The time-shared mixed behavior is somewhat similar to the mixed motility (running, rotation, and random individual motion) mode of larger cell-clusters under the gradient of chemoattractant as discussed in [35,36]. Subsequently, we postulated that the cells in a confluent population might possess the two different modes of cell motion which were supported by doublets: namely, the flocking and rotation. Indeed, the velocity correlation maps (Figs 1F and 4F) strongly supported that idea: positive (negative) correlation along the migration (perpendicular) axis is evident. Besides, the correlation length along the migration axis (y-axis) of the cells within confluent population also revealed a similar transition as a function of S, which eventually led to a correlation length comparable to a single-cell's diameter (S7A Fig). The negative correlation along the x-axis, reminiscent of rotation with adjacent neighbor, existed over a wide range of S (see S7B Fig).

Discussion
Super-diffusivity of freely crawling animal cells (on 2D substrate) is not an uncommon phenomenon, being supported by many different types of cells, and its underlying mechanism has been explored in numerous previous studies [2,[37][38][39][40]. One of the essential properties for the super-diffusivity is the self-propelling force, generated by polymerization of actins at the moving front, and depending on the strength (and time duration) of this force (as well as other factors influencing cortical actin dynamics) the phenotypic cell motility would be rather different. What is surprising in this work is that the diffusion exponent α is larger for the cells in confluent populations compared to that of freely crawling cells.
The observed phenomena can be realized with cells that have a proper degree of directional persistence and are just sticky enough to hold on to each other. For the chosen set of (S, E), the cells forming a doublet rotate about each other. But in other regime the cells can be co-moving (flocking). The transition between the two different modes is continuous. The enhanced directional persistence and super-diffusivity in confluency is a consequence of acquiring a small but finite velocity-velocity correlation length along the direction of instantaneous moving direction. The small correlation length along the y-axis is a reminiscence of the co-moving cells in a doublet forming a "two-cart vehicle". On average, the cell motility in confluency can be viewed as a random mixture of cell doublets co-moving or rotating together in a time-shared manner.
Cell-doublets had a dramatically different motility from that of a freely moving single cell: The doublets rotated periodically without much drift. The doublet-rotation must be a cooperative work of the constituent cells, since isolated MDA-MB-231 cells do not exhibit rotation by themselves although there are some cells exhibiting intrinsic chiral rotation [41]. We note that similar coherent angular motion was reported earlier with human breast epithelial cells from mammoplasty and the phenomenon was discussed in connection to acini formation [42]: The authors found that blocking E-cadherin function naturally disrupts cell-to-cell adhesion to stop the rotation and that the doublet rotation involves non-trivial cortical actin dynamics. Incidentally, MDA-MB-231 cell lines are also derived from human breast epithelial tumors. So, the doublet rotation that we observed could be a salient feature of breast epithelial cells.
Rotating epithelial cell-clusters within geometric confinements were reported earlier [43,44], in which the coherence of the rotation was controlled by the size of the confinement. The authors found coherent cell rotation became disrupted in large confinements as the cells in the middle initiated an instability throughout the whole cluster [44]. In another study, freely migrating lymphocytic cell-cluster under the influence of chemical gradient exhibited intermittent rigid-body rotation amid persistent movement and directionless individual cell movements [35]. Incidentally, the frequent position-swapping between the rim-cells (or leader cells) and core-cells during the rotation is rather analogous to the cell position-swapping discussed in this paper, although in our case there existed neither extrinsic chemotactic influence nor cell-cell heterogeneity (i.e., rim-vs. core-cells). A self-propelled agent-based model [36], implementing each cell's density-dependent self-propulsion, elucidated the heterogeneous propulsion between the rim and core. The coupling between the cells of these subdivisions was found to be critical as for generating rotational behavior of a cluster.
The 'position swapping' event in the quadruplet motility that we discussed in this paper was somewhat analogues to the differential rim-core motility of a cell cluster discussed in [36]. That is, for the case of cell quadruplet rotation, the "rim-bulk graded motility" might apply transiently when there comes a position-swapping event, in which 4 cells forming a rim-only state produces a 3-cell rim and 1-cell bulk state, which is a reminiscence of frustrated state where the cell in the core drags and slows the 3 cells in the rim. However, in our case there was no pre-assigned, space-dependent graded motility.
Investigations on cell doublet rotation were conducted earlier via phase-field model simulations of doublets under confinements [45,46]. Studying the effect of different polarity mechanisms on rotation demonstrated that it can be generated either by front-front inhibition (mechanism avoiding front-front adhesion) or by polarity-velocity alignment (tendency of cell aligning its polarity to instantaneous migration velocity) [45]. As an extension of this model, implementing intercellular diffusion of inhibitory component at the contact site, confinementfree rotation of adhesive cell-doublet was shown [46]: Here, the rotation was viewed as a consequence of continuous 'walk-past' movement, crawling in the direction away from the contact while maintaining cell-cell adhesion. Fundamentally, these models employed mechanisms in which cells adhered to their partner differentially depending on the latter's cellular polarity. On the other hand, our current CPM model had no pre-built-in elements reflecting how single cell (internal) properties change upon cell-cell contact or depend on local cell density. In fact, we had no experimental data to address such issues for the MDA-MB-231 cells and incorporate them realistically into our CPM. The same was true for cell polarization: Simply, we do not know if there exists any direct polarization-polarization interaction among neighbors in our MDA-MB-231 populations.
In our CPM model, neighboring polarization vectors interacted only indirectly by the cell population dynamics, as the neighboring cells in contact could alter the cell shapes each other. For example, the interaction between two polarity vectors of a cell doublet is mediated by the interface between two neighboring cells (see Fig 8): Initially, there is a pair of initially unpolarized, mirror-symmetric, adherent cells; then, imagine their interface starts to deform due to random fluctuation (a basic harmonic mode is assumed); with that, their shapes (centroids) also change (move); consequently, the cells acquire polarity vectors (Fig 8A); the polarity vectors then enhance the interfacial deformation to initiate a rotation via the CPM energetics part incorporating a polarity-velocity alignment tendency (Fig 8B); the same energetics keeps amplifying the harmonic deformation ( Fig 8C) to a steady profile and maintains The result of the pixel updates will be (probabilistically) highly like to be the image as depicted on the right [the yellow (green) pixel replaced by the green (yellow) within the red box on top (bottom)]. Consequently, a new Δx is created (for each cell) towards which the nextp will be attracted. (c) As the CPM update repeats over the time along the whole interface, the initial small-amplitude harmonic fluctuation will grow and rotate in the clockwise direction. Of course, this polarity driven (pure) rotation will compete with 'curvature-smoothing" energetic components in the CPM. In fact, the competition seems to make the actual rotation not like a rigid-body rotation but more like a periodic sequence of position swapping events as discussed in Fig 5. https://doi.org/10.1371/journal.pcbi.1009447.g008 the rotation. Of course, the above sequence of events happens only for an appropriate parameter regime of (S, E).
One limitation of current CPM was that it did not have a built-in contact-inhibition-locomotion (CIL) mechanism, where CIL refers to the phenomenon in which the physical contact between (colliding) two (more more) cells initiates some chemical reactions that affect either the mechanics or the adhesive property of the cells involved. As discussed in the two earlier papers [45,46], CIL could be a critical element as important as cell mechanics (e.g., shapes), propulsion strength, cell polarity, or adhesion, as for governing the behavior of cells in contact. Unfortunately, we do not know if MDA-MB-231 cells had any significant CIL and if it they did how they work. The same was true for the cell polarization interaction. In our model, the polarization (as well as velocity) vectors interact only indirectly by the cell population dynamics (i.e., cell contacts): Each polarization vector changes as the cell's centroid position, which is affected by the shape of the cell that is influenced by its adjacent neighbors, moves in space. Simply, we do not know if there exists any "direct" polarization-polarization interaction among neighbors in our MDA-MB-231 populations.
Another limitation of our CPM model was that the polarization vector was not emergent from the CPM-calculated cell shape but was set externally. In other words, the polarization vector was not directly born out from the cell shape. We should indicate that there are variants of the CPM where polarization and persistence emerge from the underling actin dynamics [47,48], and therefore motility and persistence are highly shape dependent. Currently, it is not clear how critical to have the polarization vector evaluated from the cell shape directly. In the future, it would be interesting to simulate these variant CPMs that are based on realistic actin dynamics. Likewise, it will be also worthwhile to add on additional elements such as CIL and/ or direct polarity-polarity interactions in the CPM to see how our current observations change.
The enhanced diffusivity of the confluent MDA-MB-231 cells in this report draws important similarities to the enhanced motility of human kidney cells expressing MOCA protein [25]. MOCA-expressing cells exhibited cluster formations due to enhanced cell-cell adhesion in contrast to the control. Higher diffusivity of the interacting MOCA-expressing cells was attributed to their correlated movement evidenced by the longer correlation length (amounts to roughly a cell-diameter) than that of the control cells. In the meantime, recently reported contact enhancement of locomotion (CEL) [49] considered active Brownian particles which switched to a ballistic mode of motility upon a collision, recapitulating the experimental observation of Dictyostelium discoideum cells exhibiting greater persistence in higher cell-density. The CEL study employed cell cultures in which the cells can freely explore the free space before contacting another cell, contrary to the confluent situation of our densely packed samples. In other words, CEL does arise not by cell-cell adhesion but by intracellular changes brought by physical contacts. On the other hand, the NED was attributed to the local adhesion-mediated interaction. In any case, all of these studies clearly underscore the importance of understanding population dynamics as for the cell motility: Clearly, "More is different" [50].
If the NED effect were demonstrated here on MDA-MB-231 cells, we would expect it to be relevant to many different amoeba-like immune as well as tumor cells in a much broader context. Generally, NED property would endow efficient exploration or invasion inside a closepacked tissue environment, which may be critical for immune or metastatic cells, as well as in morphogenesis. As mentioned, the observed NED effect has a strong similarity to the CEL [49]. We suspect that CEL and NED could share common biochemical origins at their microscopic scale. In fact, our earlier observation on the contact interaction between a regular MDA-MB-231 cell and an enormously expanded senescent MDA-MB-231 cell clearly indicates that MDA-MB-231 cells prefer to move along the other cells' boundaries [51].
Although MDA-MB-231 cell lines are monoclonal, the cells are known to be quite heterogeneous in many different aspects [52][53][54]. Indeed, the tumor cell populations in our cell cultures exhibited a wide spectrum of diffusion exponent α as shown in Fig 1C. In principle, this heterogeneity in motility can be caused by extrinsic as well as intrinsic stochastic gene expressions. The simulation results in Figs 4-6 were based on a fixed set of parameters for the sake of simplicity. We have also carried out more realistic simulations by incorporating distributed values of parameter S across the population of cells. Resulting MSD profiles of confluent as well as freely crawling cells became more similar to the experimental results ( S8 Fig). Notably, the variation of MSD profiles of confluent population was significantly smaller than that of freely crawling cells (S8A Fig). Also, having a wide dispersal of E (with a fixed S) for the population created no significant variation in MSD profiles; thus, the self-propulsion force S (over E) seems to be more responsible for the broad distributions seen in the experimental results.

Conclusion
Our bottom-up approach successfully delineated the physical behaviors of freely crawling cells, cell-doublets, cell-quadruplets to the cells within confluent population in a systematic way. More specifically, we presented the counter-intuitive experimental observation of enhanced diffusivity of MDA-MB-231 cells in confluent 2D cell-population. Also, we provided unequivocal evidence of cell-doublet rotation, which could be viewed as a continuous sequence of cell position swapping events, rather than a steady rigid-body rotation. The doublet rotation phenomenon was conserved even in cell-quadruplets, yet with an additional complexity: Unimpaired, coherent rotation of four cohesive cells was often hampered by intermittent events of pairwise cell-position swapping, each of which caused a reordering in the cyclic rotation sequence of the four cells involved. This swapping event within cell-quadruplets might be a reminiscence of sudden doublet rotation reversal. In the same light, we might view the uncaged, individual migration of MDA-MB-231 cells in confluency as an intermittent sequence of position swapping with immediate neighbors based on the observation of short spatial correlation length, and negatively correlated movement with immediate neighbors within population.
Essentially, all experimental observations made with MDA-MB-231 cells were successfully recapitulated by a CPM with a suitable choice of S and E, which were chosen self-consistently so that cell-doublet behavior, as well as cell-population dynamics, could match experimentally measured characteristics. With the same set of parameter values, the remarkable dynamics of cell-quadruplet was also reproduced successfully. Again, what seems to be the most significant of our current CPM study is that even without having any complex CIL or direct polaritypolarity interactions, the toy model can still well reproduce all n = 1, n = 2, n = 4, and n!1 phenomena rather well based on a single set of parameters (S, E). Finally, we should point out that our current work did not address any detailed intracellular programs or molecular pathways regulating cell-cell adhesion or directional persistence of crawling cells, on which there have been numerous studies by others [55][56][57][58].

Cellular potts model description
Briefly, CPM is a lattice-based cell (or cell population) model, in which cells are defined as a simply-connected group of lattice sites, evolving in space and time based on Monte Carlo simulation with Metropolis algorithm. It has been widely used for describing various phenomena in cellular biology [59]. The update process of the model is the following: 1) we randomly pick a (uniformly distributed) lattice site and an 'neighbor' site from a local area of given radius centered around this site. 2) The probability of the initially chosen site accepting the neighbor site's cell-type is dependent on Metropolis probability p accept which is a function of the change in the total Hamiltonian generated when the chosen site actually accepted the new cell-type. The p accept is given as 1 if the change in total Hamiltonian ΔH was negative (always accept), and decays exponentially (e −ΔH/T ) if ΔH was positive with the temperature T determining the decay rate.
In this study, the Hamiltonian includes three energy components associated with the cells' target surface areas and volumes, interfacial (e.g., cell-to-cell) adhesion, and the degree of a cell's activeness. Mathematically, where A target is the target area and s target is the target roughness defined as the fraction of circumference of a circle which has the same area as the given cell. Greater s target produces more tortuosity and ambiguity in interfaces. λs are strengths of the energies associated with area and surface constraints which penalize deviations from the respective target values. E t i t j defines the interfacial adhesiveness between two cell-types τ i and τ j . There is an additional term in the change of Hamiltonian ΔH motion which gives bias to the migration in the direction of the 'polarity' p ! defined for each cell: represents displacement of the cell's centroid due to the hypothetical cell-type acceptance described. The polarity updates at each Monte Carlo step (MCS) according to D p i represents a displacement of the centroid from the previous MCS and t memory sets the decaying rate of the polarity. This update rule represents alignment of the polarity (which can also be viewed as propulsion 'force') to the accumulated moving velocity, in which the constant 1 À 1 t memory controls the effective number of accumulated velocities: If t memory !1, then the polarity update rule dictates p ! i ðtÞ ¼ p ! i ðt À 1Þ þ D r ! i meaning entire history of a cell's trajectory equals the polarity. If t memory = 1, which is its minimum value, then the polarity is identical to the instantaneous displacement ( p ! i ðtÞ ¼ D r ! i ) which is highly stochastic. In the simulation, 1 MCS corresponds to updating sites the number of times equal to the number of total sites in the system.
The set of parameter values we fixed for all CPM simulations are: t memory = 4, λ area = 1, λ surface = 1, A target = 100, s target = 0.9, E medium−medium = 0, E cell−medium = 0, and the temperature was set to 10. Throughout this paper we varied only two parameters: E = E cell−cell , the interfacial energy factor representative of the strength of adhesiveness between cells of the same type, and S which controls the strength of the energy associated with the alignment between polarity and velocity as described in [34]. S determines not only the directional persistence but also the migration speed of a cell. Physical units for the simulation were set as the following. First of all, we set the grid mesh size to be 3 μm: It was determined based on an average diameter of MDA-MB-231 cells (~30 μm) and the model cell's target area A target which was fixed to 100. Assigning a physical time unit to one MCS before we run CPM simulation was in principle not possible, and it came only after we established the 2 phase diagrams (Fig 7A and 7B) for which we did not need to know the physical time scale beforehand. We identified the set of S and E where two relevant level curves crossed, used the corresponding values of S and E to compute the model cells' instantaneous speed v, which then was compared with the experimentally measured one, and finally set 1 MSC = 2 min.

Culture of MDA-MB-231 and sample preparation
Once MDA-MB-231 cells growing on a culture dish in an incubator which maintains 37˚C and 5% of CO2 reached~70% confluency, the cells were trypsinized with Trypsin-EDTA solution 10X (Cat. No. 59418C, Sigma-Aldrich) diluted to 5X with DPBS (Cat. No. D8537, Sigma-Aldrich) for 3 min in the incubator, after removing the DMEM culture media (10% FBS). After being centrifugated, supernatant was discarded, then DPBS was added and well-mixed. This washing procedure was done twice. Subsequently, culture media was added, mixed and counted using hemocytometer. The suspended cells were added on a 35 mm culture-dish and kept in the incubator overnight. For imaging, the sample was prepared by fully adding culture media to the dish and covered it with a thin PDMS to prevent evaporation.

Time-sequence imaging
Live-cell imaging of the MDA-MB-231 cells was carried out by placing the sample inside a custom-made small, thin cylindrical incubator which was placed on the x-y stage of a microscope (inverted, IX71, Olympus). The incubator was heated by an insulated heating pad which was connected with a temperature controller (SDM9000, Sanup, Korea), while the temperature inside was monitored by a (PT100) thermometer to maintain 37.5 ± 0.1˚C. To prevent water condensation, two indium tin oxide (ITO)-coated windows along the imaging axis were electrically heated. Mixture of CO2 (5%) and air (95%) were continuously perfused to the incubator. Acquisition of the time-sequence images were done using PCO CCD Camera (1600s, Germany) with 4X (NA 0.13) objective lens and Micromanager software. The time-interval between each acquisition was fixed to 2 min throughout the experiments in this study.

Statistical analysis
Statistical analysis and mathematical fitting was performed using Matlab and home-built softwares. All statistical data were presented as mean ± standard deviation and assessed for statistical significance using a one-way ANOVA. where N is the number of cells, and A is the total area. The unit length is μm. This result illustrates that cells tend to adhere to each other forming small colonies; and this tendency becomes less pronounced as the cell population gets confluent.