Research ArticleBIOPHYSICS

Frustration-induced phases in migrating cell clusters

See allHide authors and affiliations

Science Advances  12 Sep 2018:
Vol. 4, no. 9, eaar8483
DOI: 10.1126/sciadv.aar8483


Certain malignant cancer cells form clusters in a chemoattractant gradient, which can spontaneously show three different phases of motion: translational, rotational, and random. Guided by our experiments on the motion of two-dimensional clusters in vitro, we developed an agent-based model in which the cells form a cohesive cluster due to attractive and alignment interactions. We find that when cells at the cluster rim are more motile, all three phases of motion coexist, in agreement with our observations. Using the model, we show that the transitions between different phases are driven by competition between an ordered rim and a disordered core accompanied by the creation and annihilation of topological defects in the velocity field. The model makes specific predictions, which we verify with our experimental data. Our results suggest that heterogeneous behavior of individuals, based on local environment, can lead to novel, experimentally observed phases of collective motion.


Collective motion is an emergent phenomenon in large groups of individuals where the motion can arise from purely local interactions. This phenomenon occurs across scales in systems ranging from bacteria to fish (1). Studies of these systems in the thermodynamic limit of infinite size have revealed a number of interesting features including long-range, scale-free correlations and a discontinuous phase transition (2). These thermodynamic limit studies have spurred interest in hydrodynamic and mean field theories to describe these phenomena (3). Finite groups display collective motion that closely models schools of fish or flock of birds with boundaries. Their kinematics are characterized by unique behaviors including guidance by asymmetric boundaries (4) and the ability to simultaneously exhibit different phases of motion (5). Collective motion of groups was found to exhibit three distinct phases: running, rotating, and random (6, 7). In the running phase, the individuals are all more or less aligned, leading to a large translational velocity of the cluster center of mass. In the random or disordered phase, individual velocities are uncorrelated and there is very little overall motion of the cluster. In the rotating phase, on the other hand, the cluster rotates as a whole around a common center. While the running and random phases have analogs in infinite systems, the mechanisms that can give rise to rotations are less clear.

Through simulations, confinement has been shown to be one mechanism that results in uniform populations of self-propelled particles exhibiting rotational modes (8, 9). Simulations of large groups of unconfined agents can also display rotational phases or milling states, where the group rotates in a donut shape under certain conditions (1012). However, to achieve these rotational milling states, the agents interact over a range up to tens of times the size of an individual agent and form a low-density “hole” at the center, where the defect in the velocity orientation field resides. Groups of cells have also been shown to display these rotations, although it is unlikely that cells can interact much beyond their nearest neighbors (13, 14). This rotational motion has been studied both experimentally and using simulations for small groups of cells confined to different geometries (15, 16). Perhaps even more remarkably, unconfined cell clusters have also been observed experimentally to show transient rotational phases (17, 18), and this behavior has been speculated to promote chemotaxis.

Here, we use a generic agent-based swarming model, which only allows short-range, nearest-neighbor interactions and unconfined space, similar to previous models found in the literature (1921), to address the phenomenon of transient rotations in unconfined clusters. We show that a possible mechanism for driving cluster rotations is density-dependent cell propulsion. This density-dependent propulsion may be caused by contact inhibition of locomotion (CIL), whereby cell protrusions are inhibited by the adhesions between cells (22, 23). This contact inhibition causes cells at the cluster core, surrounded by other cells, to move more slowly than those at the edge of the clusters, which have a lower local cell density (24). The result is an outer rim of cells that move faster than central core ones and display stronger alignment interactions. We also find that decoupling the motion of the rim and central cells suppresses any rotational motion, suggesting that it is the coupling of two systems with different motilities (rim and core) that leads to rotational phases. Specifically, rotations arise in this model when the internal noise is such that the rim cells are in an ordered state with respect to velocity alignment, while the core of the cluster is disordered. The coupling of these two systems (rim around core) results in a frustrated state of the ordered rim being pinned by the disordered core, which is unable to move together in a smooth running phase with the rim. The whole coupled system is then able to relieve this frustration most effectively by existing in a rotational phase, where the more ordered rim is able to pull the disordered core, which is pinning it in place, into a rotational phase. This model successfully captures the dynamics of transitions between the modes of motion and proportions of time spent in each phase observed experimentally (17). In the experiments, when the cell clusters are subjected to a chemical gradient (25, 26), there is an increase in the proportion of running phase and a decrease in the rotational phase. This trend is also captured by our model when a chemical gradient is introduced. Furthermore, our model predicts an increase in the proportion of rotating phase with the size of the cluster, which we confirmed with experimental data. Together, our results suggest a novel form of frustrated interactions between behaviorally different parts of the same cluster that can lead to different collective dynamics—a finding that may have applications beyond the context of cellular clusters.

Experimental measurements of cell cluster phases

In the experiments, malignant B and T type lymphocytes were placed in a chemical gradient of CCL19, where they assemble into clusters and move toward higher CCL19 concentration (see Materials and Methods) (17) . During this motion, these clusters were observed to exhibit running, rotating, and random phases as well as transitions between them (see Fig. 1A and movie S3). To quantify these phases and analyze the cell cluster dynamics, automated analysis of video recordings of the cell clusters was used to extract velocity vectors of individual cells (see Materials and Methods for details), which were then processed as described below. We first identify the mode of motion of the cell cluster by measuring the polarization (O) and angular momentum (A) with respect to the center of mass of the cluster, given byEmbedded Image(1)Embedded Imagewhere Embedded Image is a unit vector toward the cell position from the center of the cluster. Using the extracted cell velocity vectors, we were able to compute the polarization and angular momentum as functions of time. Figure 1A (bottom) shows a time trace of the polarization and angular momentum of a cluster revealing distinct regions, corresponding to phases, marked by specific combinations of high, low, and intermediate polarization and angular momentum values. Using these values and the criteria described in section S3, we can then label the phase of motion of the cluster for each time point. We see all three phases being represented and the spontaneous transitions between them (Fig. 1A and movie S3). Motivated by these results, we develop a model to explain these observations. We then test the predictions of our model regarding cluster size dependence, dynamics of topological defects, fluidity, and response to the chemical gradient with further analysis of our experimental data.

Fig. 1 Analyzing and modeling cell cluster phases.

(A) Top: Experimental images of a cell cluster in each of the three phases, where the blue cells show positions at a certain time and red shows the positions of the same cells 15 s later. These positions are then used to calculate the cell velocities shown in yellow arrows. Bottom: Time series of the magnitudes of group polarization and angular momentum of the cell cluster. The colors along the bottom axis show the phase of the system with time (red, running; blue, rotating; green, random) for experimental data. (B) Schematic of the model. Green direction indicators show the directions of the neighbors of the gray cell, and the green indicator on the gray cell shows the alignment interaction (Embedded Image). The orange arrows show the Lennard-Jones interaction with each neighboring cell, and the red arrow is the total Lennard-Jones interaction (Embedded Image) on the gray cell. Finally, the blue spring denotes the surface tension interaction (Embedded Image). Note that it only exists between the gray cell and its second nearest neighbors that do not have cells interrupting the path between them. (C) Top: The proportion of time that the cluster spends in each phase [simulations (plain) and experiments (crosshatched)], along with a typical illustration of what each phase looks like in the simulations, with velocity vectors as black arrows. The cluster size for simulations is N = 37 cells, while experimental cluster sizes are distributed with a peak between 35 and 40 and a mean of about 50 (see fig. S7A). Bottom: Time series of the magnitudes of group polarization and angular momentum from simulations of a uniform cluster (dashed) and a cluster with behavioral heterogeneity (solid, corresponding to the point marked in Fig. 2B).


Cell clusters are modeled as groups of particles that move with overdamped dynamics in two-dimensional (2D) continuous space (see section S1). Cells are initially arranged in a circular disc, with velocities pointing in random directions. Cell velocities are determined by their internal self-propulsion (with magnitude pi), as well as physical interactions between cells, such as adhesions and collisions (20, 27). All these interactions assume that cells communicate with each other by contact, so that cells interact with other cells that are within a distance 1.26 r, where r is the average cell diameter, which is small enough to only include nearest neighbors. The cell diameter is selected from a Gaussian distribution, as uniform cell sizes lead to crystal lattice effects that are unlikely to exist in the experimental cell system (see section S2 and movies S1 and S2 for comparison). Finally, the velocities of the cells are subject to some uniform and uncorrelated noise (Embedded Image) due to random traction forces with the substrate and the random nature of the protrusions that cells use for propulsion. Cell positions are then updated according to their individually calculated velocitiesEmbedded Image(2)

To determine the velocity of individual cells, a couple of interactions are taken into account. First, cells propel themselves in a direction (Embedded Image, Eq. 3) determined by the memory of their own previous polarization (Embedded Image) and an alignment interaction with the mean orientation of neighboring cells, Embedded Image, with interaction strength αEmbedded Image(3)

Cell velocities are calculated as arising from the forces described here and illustrated in Fig. 1B. The self-propulsion magnitude is set by p. Additionally, cells experience volume exclusion and adhesion with neighboring cells, which are modeled as arising from a Lennard-Jones force (Embedded Image). We also include a spring-like surface tension interaction (Embedded Image) that is longer-range and only acts between cells that are separated by a distance larger than the first nearest-neighbor interaction separation (1.26r) if there are no interrupting cells directly between the two interacting cells (Fig. 1B). This surface tension force only acts occasionally for cluster rim cells that may be drifting away and is necessary to maintain compact and roughly circular clustersEmbedded Image(4)

Cell clusters in the experiments are subject to chemoattractant gradients. To implement a chemical gradient into our model, we introduce an additional term into the calculation for the cell propulsion direction, Embedded Image, by replacing Embedded Image with Embedded Image in Eq. 3, whereEmbedded Image(5)and the sum j is over each distinct pair of adjacent neighbors of cell i. Embedded Image is a vector pointing in the direction bisecting the angle subtended by the centers of the cells of the neighbor pair at the center of cell i, with a magnitude equal to the arc length between the two neighbors (see section S1.2). Here, g reflects the strength of the influence on propulsion direction from the chemokine gradient per unit distance of exposed cell edge arc length, c′ is the change in chemokine concentration per unit distance, and y is the distance (in micrometers) from a concentration point of 0 ng/ml . This results in a gradient force in the direction of the most vacant region around a cell, with a magnitude proportional to the size of the vacancy and the strength of the chemoattractant gradient. Positions and velocities of cells within the simulated clusters are processed in an analogous manner as the corresponding experimental quantities to compute metrics for comparison with experiment including polarization, angular momentum, and the time spent in different phases, among others.


Uniform cell clusters do not rotate

In the case where all cells within the cluster behave identically, specifically the propulsion magnitude p is the same for all cells, the cluster remains in a single phase throughout the simulation. The dashed line in Fig. 1C (bottom) shows a time trace for a cluster in the running phase, where the group polarization remains high and the angular momentum remains low throughout. When compared to a time trace of the same quantities measured in experimental cell clusters (Fig. 1A, bottom), the features of the time traces are very different. In the experiments, the group polarization and angular momentum fluctuate from high to low values corresponding to spontaneous transitions of the cluster between various phases of motion (shown in Fig. 1A, bottom). In the simulations, on the other hand, the cluster undergoes a transition from a running phase to a random phase only with increasing noise or decreasing propulsion. Figure 2A shows the proportion of time spent by clusters in each of the three phases plotted against noise and propulsion. The diagonal line of the transition between running at low noise and random at high noise is the well-known noise-driven transition seen in Vicsek swarming models (2, 28).

Fig. 2 Collective phase proportions of density-dependent propulsion model.

(A) Proportion of time spent by the cluster (N = 37 cells) in each of the three phases plotted against propulsion p and noise Embedded Image for a cluster where all cells behave identically. (B) Phase diagram of the proportion of time spent in each of the three phases for a system with neighbor number–dependent propulsion where the rim cells (those with 3.67 neighbors) have a propulsion of prim = 8. pcore is the propulsion of core cells (those with 6 neighbors), and η is the magnitude of the noise. The black “x” shows the point where the time series and phase proportions shown in Fig. 1 (B, bottom, and C) are taken.

The running and random phases seen here are similar to those seen in experiments; however, the transition between the running and random phases in phase space is very sharp, and there is very little overlap or mixing of the phases. Experimentally, cell clusters are observed to spontaneously switch between running, rotating, and random phases, which means that they coexist within this parameter space. Alternatively, one could say that cells change their internal parameters, such as the propulsion, p, or noise, η, so that they cross over the transition between the phases. However, it is implausible for the entire cell cluster to change internal parameters in a coordinated way. Additionally, the uniform clusters show a very low level of cell rearrangement or fluidity within the cluster, whereas, in the experiments, cells are observed to move between the rim and the core of the cluster regularly. We therefore conclude that additional features of the real system must be incorporated into the model to recapitulate a rotational phase and transitions between phases within a single set of parameter values, as well as large-scale cluster rearrangement.

Heterogeneous behavior: Density-dependent propulsion produces rotations

An aspect of cellular behavior that is missing in this description is the possibility that cells may behave differently in different regions of the cluster, such as the periphery or the interior. Rim cells have increased propulsion compared to inner cluster (core) cells due to reduced CIL, which causes cells adhered to other cells to form fewer protrusions than cells that have more open space around them. We implement this effect by scaling the propulsion with the number of neighbors, increasing for cells with fewer neighborsEmbedded Image(6)Here, ni is the number of neighbors around cell i. prim and pcore are the propulsion of the rim cells (average of 3.67 neighbors) and core cells (average of 6 neighbors), respectively. A similar inverse relation between local density and propulsion (and therefore alignment) was explored for a semi-infinite system in (29). This relationship is also consistent with the behavior of experimental cell velocities (though those are indirect measures of propulsion), and it is to be noted that, as long as propulsion is inhibited by increased numbers of neighbors, the exact functional form of the relation does not affect our results significantly (see section S4.2).

This variation of cell propulsion causes the rotating phase to emerge and to coexist with the other phases, as seen in experiments. There is now a region in parameter space where there is a peak in the rotational phase at low values of pcore and intermediate noise (Fig. 2B). In this region, there are proportions of all three phases that are close to those seen in experiments. The point highlighted in Fig. 2B is an example of a location in parameter space where the simulations closely match the experiments. The time series for the group polarization and angular momentum at this point in parameter space is shown in Fig. 1C (bottom) and is very similar to the experimental time series (Fig. 1A, bottom). Furthermore, when the proportion of time spent in each of the three phases is compared to the experimentally measured values, they match very closely (Fig. 1C, top). We next investigate, more closely, the mechanism that drives the rotating phase.

Stable rotations emerge in a simple rim-core model

The values of propulsion for individual cells in our simulation with heterogeneity are mostly close to pcore, except for those near the periphery where it rapidly climbs to an average of prim (see section S4). This prompts us to consider whether the behavior of the system can be understood as arising from the coupling of two different systems—a ring-like rim with a higher propulsion and a uniform core with a lower propulsion. We first examine the rim cell system by confining the cell positions to lie only on the circumference of a circle without any core cells and then assigning them a fixed value of propulsion (= prim) independent of neighbor number. This system is described by the phase diagram shown in Fig. 3, where the contours are boundaries of regions where the proportion of time spent in the corresponding phase exceeds 30 and 50%. The diagonal dashed contour lines belong to the isolated rim system. This phase space shares some characteristics with the uniform cluster phase diagram (Fig. 2A), such as the transition from a running to a random phase with increasing noise or decreasing propulsion, and the lack of a rotational phase. Because of the lower number of neighbors in the rim case, however, the slope of the diagonal running-random transition line is smaller compared to the case with a uniform cluster.

Fig. 3 Collective phase proportions of the coupled rim-core model.

Proportion of time spent by the system in each of the three phases as a function of propulsion p and noise Embedded Image for a ring of 18 cells confined to a circle with propulsion prim = p. Dashed contour lines indicate regions (shaded) where the proportion of time spent in the corresponding phase exceeds 30% (blue) and 50% (red). Solid contour lines show the same contours but for the rim confined to a circle with prim = 8, coupled with a core of cells with pcore = p, and a full cluster size of N = 37. Note that the rotational phase only has nonzero values for the coupled system. The horizontal dashed line marks the noise value below which the rim alone would be ordered (greater than 30% running phase) with prim = 8, and the diagonal solid line marks the region above which a core with an average propulsion set by Eq. 7 (with pcore = p) is disordered (greater than 30% random phase).

We next couple the rim cells to the core, resulting in a ring of cells confined to a circle with propulsion prim = 8, positioned around a core of cells with propulsion pcore. We now consider the phase behavior of two uniform systems to compare with our coupled system—the isolated rim system and a uniform cluster of the same size as the coupled system. The black dashed line in Fig. 3 shows the noise value below which a ring of cells with the same prim ( = 8) would be ordered (or in the running phase more than 30% of the time). Our choice for the uniform cluster for comparison purposes is one with the same number of cells but with a uniform propulsion across all cells that is equal to the average propulsion for the heterogeneous cluster given byEmbedded Image(7)Below the black solid line, we expect this cluster, which we consider the uniform equivalent of our heterogeneous coupled system, to be ordered (running phase greater than 30%), and above it, disordered. This expectation is based on the transition shown in Fig. 2A where the propulsion is equal to the average propulsion in Eq. 7 with prim = 8, Nrim = 18, pcore = p, and Ncore = 19 for the full cluster of N = 37 cells. We notice that there is a triangular region between the dashed and solid black lines at low pcore and intermediate noise where the rim should be in its ordered state, while the uniform equivalent of our full cluster would be in the disordered state. This observation also suggests that the interior core (with a lower p = pcore) of the full system on its own would also be disordered here. In this region of phase space, the overall heterogeneous cluster shows rotations. This in turn suggests the possibility that when the overall system cannot exist in a running phase but the rim alone would be in its ordered phase, the rim is effectively pinned in place by the disordered core. Such a cluster could relieve the frustration and maximize order by existing in a rotational phase, where the rim is moving around and pulling the core with it.

When we couple the ordered rim to the random phase core in this parameter regime, we find a peak in the rotating phase (Fig. 3, solid contours). The fact that the peak in the rotational phase does not exist for the rim or core alone but emerges when the two systems are coupled suggests that the rotational phase is driven by the coupling of the ordered rim to the disordered core.

There are some differences between the ring-disc confined system and the full unconfined model shown in Fig. 2B, but these differences are mainly due to effects arising from the confinement of the rim cells to a circle. Relaxing the confinement of the rim cells, while still treating the cluster as two different coupled systems of the rim cells with higher propulsion around core cells with lower propulsion, recovers the original, fully heterogeneous model phase space (see section S5). Together, these results suggest that it is the coupling between the disordered core and ordered rim that is the mechanism behind the cell cluster rotations seen experimentally.

Transitions between phases are characterized by defect dynamics

While we have shown that the running, rotating, and random phases can coexist within our model with heterogeneous neighbor-dependent propulsion, we have not yet examined the dynamics of the transitions between these phases. To investigate these dynamics, we quantify these transitions by monitoring changes in the overall topological properties of the phases, which are easier to track. In particular, note that in condensed matter systems, including active matter (30, 31), phase transitions may be driven by the interactions and dynamics of topological defects (32). We first take a coupled rim-core cluster with rim cells confined to a circle and project the rim cell velocities onto the confining circle. We then identify a defect in this effectively 1D velocity field as a point where the velocity projections switch directions. Note that the defects, as we have defined them, exist when there is no defect in the full velocity field of the entire cluster in its running phase and vanish when there is a vortex in the cluster in its rotating phase (±1 defect in the director field).

Figure 4A shows these defects for a cluster in the running phase (top) and the rotating phase (bottom). In the running phase, the cluster has two defects of opposite signs at roughly opposite sides of the cluster, while in the rotating phase, there are no defects present. The formation and spreading apart of a defect pair coincide with the transition from rotating to running phase, while the annihilation of the pair results in the running to rotating phase transition (see movie S4). By measuring the frequency of occurrence of any given number of defect pairs in each cluster phase, we investigate the correlation between the phase and number of defect pairs (Fig. 4B, solid bars). We see a large peak in the rotating phase for zero defect pairs and a peak in the running phase for one defect pair. The random phase has a much broader peak around two or three defect pairs, suggesting that the random phase could occur when multiple defect pairs spontaneously form. We then investigate whether the correlation between the number of defect pairs and the cluster phase appears in the experimental cell clusters as well (Fig. 4B, crosshatched bars). We find that the peaks in defect number seen in the simulations also exist in the experimental clusters, suggesting that the tracking of topological defects in the velocity field of the rim cells can be used to characterize phases and transitions in active matter clusters.

Fig. 4 Defect dynamics and the transitions between phases.

(A) Velocities of the rim cells of a 37-cell cluster that are confined to a circular shape projected onto the circle for a simulated cell cluster (left), as well as an experimental cell cluster with rim cell velocities projected onto a circle relative to the center of mass of the cluster (right). In the running phase (red panels), there are two defects of opposite signs in the velocity field, denoted by the orange and blue points. There are no defects in the rotating phase (blue panels). (B) Proportion of the number of defect pairs for each phase, with a peak at zero defect pairs for the rotating phase (blue) and one defect for the running phase (red). Simulation defect pairs are shown in solid bars, while the experimental defect pair counts are shown in the crosshatched bars. (C) Pair distribution function plotted against the separation between two defects when only one defect pair exists for parameters where the cluster primarily displays a running phase [note that g(r) is calculated over the whole simulation, independent of specific phases at any given point in time]. Inset: Pair distribution function for points in parameter space dominated by rotating (blue) and random (green) phases. Experimentally measured values of g(r) are shown as black points, with a linear fit shown by the dashed black line. A similar comparison of the results of the full-model simulation to the experiments shows no major difference (see section S11.)

To observe the effective interactions of defects with each other, we calculate the pair distribution function (g(r)) for the spacing between two individual defects when a single pair exists. Figure 4C shows the pair distribution function calculated over all time throughout a simulation, independent of what phase the cluster is in at any particular point in time, with the separation (r) normalized by the maximum separation possible (half of the cluster circumference). We see that for parameter values when the cluster is predominantly in the running phase, the pair distribution indicates that the two defects will repel and largely exist at maximum separation. When the cluster is mostly in the rotating or random phase (Fig. 4C, inset, blue and green, respectively), the interaction is also repulsive, although the slope is much smaller for both of these cases than in the running case, so the defects only repel weakly, increasing their chance of annihilating and transitioning out of the running phase. The black points in the inset indicate the value of g(r) measured experimentally with a linear fit shown as a dashed black line. We see that the experimental cluster defect interactions are close to those seen in the simulations when the rotating phase emerges. The stochastic transitions between phases are driven by noisy processes that can be characterized by the dynamics of topological defects. Thus, system-wide parameters have a significant influence on the interactions between topological defects, which, in turn, controls the dynamics of the defects, the formation and annihilation of which are correlated with cluster phase transitions.

Cluster size dependence: Larger clusters exhibit a higher proportion of rotations

We next examine the effect of cluster size on the phase diagram. Figure 5A shows regions of the parameter space with a proportion of rotating phase greater than 30% for different cluster sizes. Compared to the predictions of the simulations regarding the proportion of the rotating phase, the predictions regarding the running and random phases have a very large spread across the parameters tested. Thus, we focus on the more significant rotational phase shown here (see section S6 for other phases). Again, the dashed black line shows the noise value below which the rim alone would be in an ordered phase, while the colored dashed lines show the transition between the running and random phases for a uniform system of average p (Eq. 7) for each different system size. Our results indicate that larger clusters have a higher proportion of rotational phase, while smaller clusters are less likely to rotate.

Fig. 5 Cluster size and chemical gradient dependence.

(A) Proportion of time spent in the rotating phase by a system with neighbor number–dependent propulsion as a function of pcore and noise Embedded Image for prim = 8. The black horizontal dashed line marks the noise value below which the rim alone would be ordered (corresponding to the dashed transition line in Fig. 3 with prim = 8). The diagonal lines show the noise value above which a uniform system with average propulsion p would be disordered (red, N = 19; blue, N = 37; green, N = 61). The shaded regions are where the clusters spend at least 30% time in the rotational phase, with the same color scheme (blue, N = 37; green, N = 61; note that there is no red shaded region). Inset: Dependence of the proportion of rotating phase on system size. The shaded red region shows the range of dependence for a spread of parameter values marked with black crosses in the main figure. The experimental measurements are shown as the blue points. (B) Experimental image of a cell cluster where cells are colored red or green for visualization. The red cell labeled by the white arrow moves from the rim of a cluster into the core between the top and bottom images over a 2-min time period. (C) Fluidity of the cluster measured as the rate of exchange between the core and rim cells of the cluster, for several system sizes, for both simulations (plain bars) and experimental data (crosshatched bars). Inset: Contours for the fluidity of the cluster over the pcore-η parameter space. (D) Simulated proportion of each phase (see legend) plotted with increasing chemical gradient (gcr in Eq. 5, where r is the cell diameter), along with experimental data in the inset. The concentration gradient of chemokine in the experiments is measured in (ng/ml)/mm and shown on the x axis for the inset.

This is consistent with the idea that the coupling-induced rotational phase only exists in the area of phase space where the rim propulsion would result in an ordered state (running phase greater than 30% of the time; below the dashed black line), while the average p of the whole cluster would lead to a disordered state (random phase greater than 30% of the time; above the colored dashed lines for each size). Larger systems have a higher proportion of core cells and therefore will have a larger proportion of the phase space where the average p results in a disordered cluster (because pcore is always lower than prim), while the rim remains ordered, leading to a larger overlap of the two and a more stable rotational phase. Figure 5A (inset) shows the comparison between experimental and simulation values, where the blue points are the experimentally measured proportions of rotating phase against cluster size, and the red shaded region shows the size dependence of the simulations, with the parameter values marked by black crosses in Fig. 5A. Although these parameter values were chosen based only on the proportions of all phases exhibited by a cluster of size N = 37, the dependence on system size seen in the simulations is very similar to that of the experiments. The fact that the cluster size dependence, which is a predicted consequence of our model, agrees with experiment lends further support to the validity of the model.

Cluster fluidity decreases with increasing cluster size

Exchanges between the periphery and the interior of the cell clusters were proposed to have an important functional role in exposing “fresh” cells with unsaturated receptors to the chemical gradient (17). Figure 5B shows a cell in an experimental cluster labeled by the white arrow moving from the rim of the cluster (top image) to the core (bottom image) over a short time. To examine this feature of our clusters, we look at the fluidity of a cluster as measured by the rate of exchange between core and rim cells. A rim cell is any cell whose exposed edge is larger than a cell diameter, whereas the exposed edge is defined as the arc length of the largest uninterrupted cell-free segment of a circle around the cell with the radius of a cell diameter (see section S7). We then average the number of cells that switch between the rim and the core in each time step to obtain a measurement for the cluster fluidity.

We measure this exchange rate for the original, fully heterogeneous neighbor number–dependent propulsion model (phase space in Fig. 2B). Figure 5C (inset) shows contours for the fluidity across parameter space. It turns out that the fluidity of the cluster is significantly higher when the cluster shows a larger proportion of the running or rotating phase compared to the random phase. This is presumably due to the fact that the large-scale rearrangements or rim-core cell exchanges happen when the rim cells move past the slow-moving core cells and then mix back into the cluster. This is why the fluidity drops for high pcore values, approaching prim. Figure 5C shows the dependence of the rim-core exchange on cluster size for the simulations, as well as the experiments (crosshatched bars). The trend of decreasing rim-core exchange with increasing cluster size is more marked in the simulations but is maintained between both experiments and simulation, further supporting a rim-core coupling as the mechanism for rotational phases in cell clusters.

It should be noted that the simulations slightly overestimate the proportion of the running phase in smaller clusters at the expense of the random phase. We speculate that this overestimate might be a signature of cellular clusters maintaining a roughly constant effective prim across cluster sizes, perhaps by sensing curvature. In this case, our simulations are essentially overestimating the effective prim value because the rim cells in smaller clusters have fewer neighbors, and p increases linearly with decreasing number of neighbors. This would result in a slight overestimation of the running phase with decreasing cluster size at the expense of the random phase seen in the experiments (section S8). Similarly, for large cluster sizes, the simulations underestimate the effective prim, leading to a higher proportion of the random phase at the expense of the running phase, as well as an underestimation of the rim-core exchange in simulations due to the fact that about 60% of the exchanges that take place occur during the running phase in both experiments and simulations. To summarize, our model predicts that the exchange of cells between the rim and bulk of the cluster occurs when the cluster is in the running phase, as opposed to the rotating phase. Because the rotating phase increases in proportion for larger clusters, this offers a mechanism that explains the observed decrease in cellular exchanges for larger clusters (Fig. 5C).

Chemical gradients favor the running phase

Cell clusters can chemotax robustly up a chemical gradient (17), and it has been shown that such a collective chemotactic motion can be obtained by cells at the cluster rim having a propulsive force normal to the surface of the cluster with a magnitude that depends on the local concentration of the chemokine (17, 33, 34). We are interested in how such a gradient would affect the proportion of time that the clusters spend in each of the different phases. We implement this in the model (see the “Model” section) by a gradient force in the direction of the most vacant region around a cell, with a magnitude proportional to the size of the vacancy and the strength of the chemoattractant gradient, which causes a large outward force on rim cells and negligible force on core cells, resulting in an overall upward drift due to the unbalanced forces (17).

We introduce such a chemical gradient to a system with running, rotating, and random proportions close to those measured experimentally in the absence of a chemical gradient. We find that the gradient leads to cluster motion up the gradient, as anticipated and observed experimentally. We then measure the changes in the proportion of phases as a function of the gradient (shown in Fig. 5D). With increasing gradient, we find an increase in the proportion of the running phase and a decrease in the amount of the rotating phase, similar to what is seen experimentally (Fig. 5D, inset), while the random phase proportion stays essentially unchanged.


Cell clusters exhibit running, rotating, and random phases in experiments. We have identified, using a theoretical model, the possible cause for the coexistence of these phases in a cluster that has only short-range cohesive and alignment interactions. Similar models have been shown to have a milling rotational phase with a low-density core (10); however, these models rely heavily on long-range interactions that are inaccessible to cells. Our model accounts for the tendency of cells at the cluster rim to have increased propulsion due to less contact inhibition, and therefore stronger alignment interactions, compared to cells at the core of the cluster. We have identified a likely mechanism by which this increased propulsion can lead to rotational phases in cell clusters.

This effect involves the effective existence of two different systems within the cell cluster—a high-propulsion, ordered rim system and a low-propulsion, disordered core. When these two systems are disconnected, there is no significant rotational phase present in either of them. However, when they are coupled together, the rotational phase appears robustly, indicating that it is the coupling of these two systems that leads to the observed rotations. We find that the ordered rim is capable of dragging the disordered core with it, resulting in a solid body–like rotation of the entire cluster (see section S9). This behavior, whereby the ordered phase induces large-scale coherence in the adjacent disordered phase, is somewhat reminiscent of the coupling between superconducting and normal metals (proximity effect) (35).

Our simulations exhibit multiple features that are seen experimentally, including spontaneous transitions between the different phases along with correlated topological defect dynamics, the cluster size dependence of the proportions of phases and the fluidity, as well as the response to chemical gradients. Our results indicate that larger clusters show an increased proportion of the rotational phase. An increase in rotational motion with cluster size has also been observed in simulations of flocks and experiments with fish schools (6), suggesting that this mechanism for rotations may extend to systems other than only cell clusters.

With increasing chemokine concentration gradient, the cell cluster spends an increased proportion of time in the running phase. An interesting consequence of this is that because most of the rim-core exchanges take place while the cluster is in the running phase, we would expect the exchange to also increase with increasing chemical gradient. This is what we see in the simulations (see section S10) with a 50% increase in exchange over a fourfold increase in the gradient. These results are supportive of the conjectured functional benefit of exchanges of rim and core cells in maintaining robust chemotaxis (17). An increase in rim-core exchange allows the cells on the rim in high-concentration gradients to shuffle back into the center of the cluster to replenish their chemical receptors, which become saturated while they are on the exposed rim of the cluster. At the same time, this brings core cells with unsaturated receptors to the rim, allowing more chemokine sensitivity. Thus, the cluster can use its collective dynamics to ensure a more robust response to gradients compared to individual cells. The emergence of exchanges between the periphery and interior, especially with increasing directional input at the periphery, might be of importance to flocks and swarms, where sharing the inherent advantages/disadvantages of being at the core/rim (like temperature extremes in penguin colonies or the threat of predation in fish schools) is beneficial to the group as a whole.

Together, our results show that the rotations induced by rim-core coupling hold across a range of system sizes, propulsion strengths, noise values, and even in the presence of directional forcing. They may even extend to 3D rotations (3639), suggesting that the coupling between two swarming systems that are in different ordered phases can lead to interesting behaviors not seen in either system alone. Heterogeneous behavior within a single group is a robust mechanism that cells or other types of swarming organisms may use to enhance rotating phases or other phases that would be unlikely or impossible to achieve otherwise.


The human chronic lymphocytic leukemia (CLL)–derived cell line JVM3 was obtained from the DSMZ (Deutsche Sammlung von Mikroorganismen und Zellkulturen) collection. They were authenticated by analysis of B cell surface markers and tested regularly for CCR7 expression and to assess their mycoplasma-free state. Recombinant human CCL19 chemokines were from PeproTech, and Alexa Fluor 647–labeled CCL19 chemokines were from Almac. Chemokines were aliquoted and stored according to the manufacturer’s instructions. We analyzed the motility of cells exposed to chemokine gradients by video microscopy using collagen IV–coated 2D chemotaxis slides from ibidi, as described (17). Briefly, cells (0.5 × 106 in 10 μl of culture medium) were loaded into the central transversal channel and incubated at 37°C for 30 min to allow attachment. Gradients of CCL19 were generated following the manufacturer’s instructions. We verified the linearity of the gradients using a 10% dextran–fluorescein isothiocyanate solution. Cell migration was recorded during that time interval at a rate of one picture every 15 to 120 s using wide-field microscopes equipped with an incubation chamber for temperature and CO2 control. Images were acquired using the MetaMorph software (Molecular Devices). Migration tracks of individual cells and clusters were obtained using the Manual Tracking plugin (F. Cordelires, Institut Curie, Orsay, France) of the ImageJ software (W. S. Rasband, National Institutes of Health, Bethesda, MD). Chemotaxis plots and migration parameters were obtained with the Chemotaxis and Migration plugin from ibidi, while nuclei velocity and its magnitude were obtained as described in (17).


Supplementary material for this article is available at

Section S1. Model details

Section S2. Lattice-induced rotations

Section S3. Phase characterization

Section S4. Density-dependent propulsion

Section S5. Rim-core cluster model

Section S6. Size dependence

Section S7. Rim cell exposed edge

Section S8. Varying rim propulsion

Section S9. Solid body–like rotations

Section S10. Fluidity versus gradient

Section S11. Comparison of defects in experiment to full model

Fig. S1. Schematic illustration of the chemical gradient force on the gray cell.

Fig. S2. Collective phase proportions of a cluster with monodisperse cell sizes.

Fig. S3. The propulsion of cells plotted against distance from the cluster center of mass.

Fig. S4. Effects of density-propulsion relationship on phases.

Fig. S5. Rim-core model phase proportions, with the rim cells confined to a circle.

Fig. S6. Rim-core model phase proportions, with the rim cells unconfined.

Fig. S7. Cluster size dependence of all phases.

Fig. S8. Schematic for rim cell definition.

Fig. S9. Collective phase proportions with varying rim propulsion.

Fig. S10. Rotational slip of outer rim around the inner core.

Fig. S11. Cluster fluidity as a function of chemical gradient.

Fig. S12. Defect dynamics and the transitions between phases for the full model.

Movie S1. Lattice-induced rotations for a crystalline cell cluster, which only occurs when the cells are of identical sizes and noise is sufficiently low.

Movie S2. A system with the same parameters as movie S1 but with polydisperse cell sizes with a spread of 10% of the average cell size.

Movie S3. Experimental cell cluster transitioning between the three phases of motion: running, rotating, and random.

Movie S4. Defect dynamics as a cluster transitions from the rotating phase to the running phase and back again.

Reference (40)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: A.G. and K.C. were partially supported by NSF grants EF-1038697 and DMS-1616926, a James S. McDonnell Foundation Award, and, in part, by the NSF-CREST: Center for Cellular and Bio-molecular Machines at University of California Merced (NSF-HRD-1547848). N.G. gratefully acknowledges funding from the Israel Science Foundation (grant no. 580/12). Work in G.S. laboratory was partially supported by the Associazione Italiana per la Ricerca sul Cancro (AIRC 10168), the Worldwide Cancer Research (AICR-14-0335), and the European Research Council (Advanced-ERC-268836). Author contributions: N.G. and A.G. designed the research. G.M.-E. and G.S. performed and interpreted experiments. K.C., G.M.-E., and W.Y. analyzed the experimental data and identified various phases of motion. K.C., N.G., and A.G. developed theory and performed the modeling. K.C. performed the simulations. K.C., N.G., and A.G. analyzed the simulation results. K.C., N.G., and A.G. wrote the paper. All authors discussed the results and commented on and edited the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article