Motility-limited aggregation of mammary epithelial cells into fractal-like clusters

Cell Proliferation and Motility Are Arrested in Reduced EGF Conditions.

Human mammary epithelial cells (MCF-10A) were first tracked with single-cell resolution in “growth” media supplemented with 20 ng/mL EGF relative to “assay” media with reduced EGF (0.075 ng/mL). These cells were stably transfected to undergo a controlled EMT through the Snail transcription factor, using an inducible construct responsive to tamoxifen (OHT) (2, 41). Such EMT induction is associated with increased motility and elongated morphologies over 48 to 72 h, due to increased actin polymerization and membrane protrusion at the leading edge (42). In growth media supplemented with OHT, dispersed cells exhibited relatively fast proliferation and formed continuous monolayers that occupied almost all of the available area by 60 h (Fig. 1A and Movie S1). These cells were also highly motile and constantly rearranged themselves, even at near-confluent cell densities, which remained qualitatively consistent over a range of initial cell densities (SI Appendix, Fig. S1 A and B). In contrast, in assay media with EGF reduced to 0.075 ng/mL and OHT, cells displayed significantly slower proliferation with large unoccupied regions after 60 h (Fig. 1B and Movie S2), even from higher initial densities (SI Appendix, Fig. S1 C and D). Despite these relatively sparse cell densities, initially dispersed cells became associated with a multicellular cluster by 24 h, which was statistically significant based on Ripley’s H function (SI Appendix, Fig. S1E). The overall morphology of these clusters remained consistent through 60 h, suggesting that cell motility was arrested and rearrangements were limited. Interestingly, these clusters exhibited highly branched morphologies that merged into spanning networks over time (Fig. 1B).

For fluorescence imaging, these cells were also stably transfected with red fluorescent protein in the nucleus (H2B mCherry) as well as green fluorescent protein (GFP) in the cytoplasm. Proliferation rates were compared across conditions by detecting the number of fluorescent nuclei present over the duration of the experiment. In growth media, MCF-10A cells exhibited roughly exponential growth with an 11-fold increase in cell density over 60 h (Fig. 1C). In assay media, cell proliferation slowed and approached a plateau at 5-fold higher cell density by 50 h, analogous to resource-limited logistic growth. Cell viability remained above 99% after 60 h (SI Appendix, Fig. S1 B and D), and no changes in the (phenol red) media color were observed, indicating minimal changes in pH. Individual cell speeds were then determined by tracking nuclear motion of all cells over time, with an initial increase in speed over $\sim$10 h as cells adhered to the substrate. In growth media, MCF-10A cells initially exhibited a population averaged speed of $\sim$25 $\mathrm{\mu }$m/h, which gradually decreased to $\sim$8 $\mathrm{\mu }$m/h by 60 h, while cells in assay media started out with average speeds of $\sim$10 $\mathrm{\mu }$m/h and gradually slowed to $\sim$2 $\mathrm{\mu }$m/h by 60 h (Fig. 1D and SI Appendix, Fig. S2 A and B). This decrease in average speed occurred with comparable kinetics across all initial cell densities, suggesting that arrested motion was not dependent on overall cell density (Fig. 1D, Inset). Overall, these combined observations indicate that assay media with reduced EGF resulted in the formation of multicellular clusters and spanning networks despite diminished cell proliferation and motility.

Dispersed Individuals Aggregate into Multicellular Clusters with Fractal-Like Morphology.

Clustering dynamics were then analyzed from single-cell tracks over the duration of the experiment. Initially dispersed cells migrated randomly, but their motion was arrested upon encountering other cells (Fig. 2A). Individual cells were defined as cells whose nuclei were located at least 75 $\mathrm{\mu }$m from any other nuclei. Based on this definition, roughly 10% to 30% of individuals were typically observed at the start of the experiment, which decreased to 5% or less over $\sim$20 h (Fig. 2B). These individuals exhibited random-walk–like migration, which was quantitatively analyzed from the mean-squared displacement (MSD) per cell $⟨\mathrm{\Delta }{r}^2\left(\tau \right)⟩$. These MSDs could be fitted to a power law of the form $K_{\alpha }{\tau }^{\alpha }$, where $\alpha =1$ corresponds to random diffusion-like motion and $\alpha =2$ corresponds to persistent ballistic-like motion (43). Note that $K_{\alpha }$ is a generalized diffusion coefficient with units of $\mathrm{\mu }{\mathrm{m}}^2\cdot {\mathrm{h}}^{-\alpha }$. The statistical distribution of these fitted parameters was $\overline{\alpha }=1.23\pm 0.40$ for n = 44 individuals (from 3 independent experiments) (Fig. 2B, Inset). Nevertheless, it should also be noted that $K_{\alpha }$ also varied considerably with $\alpha$, resulting in large variations in cell motility (SI Appendix, Fig. S2C). Although there was substantial cell-to-cell heterogeneity, individual cell migration was more random than directed.

Multicellular clusters were then defined by computationally linking nuclei located less than 75 $\mathrm{\mu }$m apart, which was verified by segmenting cluster morphologies in the cytoplasmic channel (SI Appendix, Fig. S2 D and E). This relatively high size cutoff was chosen to correctly associate elongated cells with a cluster, but more compact cells within a cluster were typically located less than 50 $\mathrm{\mu }$m apart (37). For each field of view, $\sim$10 clusters were typically observed, a number which remained relatively stable over $\sim$24 h, indicating that cells aggregated irreversibly and that clusters did not dissociate (Fig. 2C). Nevertheless, the number-averaged cluster size increased steadily over time, reaching $\sim$100 cells by 60 h (Fig. 2D). This increase could be attributed to migratory individuals being “captured” at the cluster periphery, some proliferation events, and finally the merging of local clusters, resulting in a sharp decrease in the number of clusters after 20 h (Fig. 2C).

At later times, cells organized into large clusters ($>$10 cells) that exhibited a dendritic, noncompact morphology with geometrically “rough” features at the periphery (Fig. 1B). In particular, cells at the interior were tightly connected with many neighbors, while cells at the periphery were less connected with fewer neighbors, indicating a decrease in local cell density with increasing radial distance. Both the cluster morphology and this radially decreasing density were strikingly reminiscent of fractal-like structures associated with the diffusion-limited aggregation of colloids (23–25). A quantitative signature of fractals is the fractal dimension $D_f$, which scales the radius of gyration $R_g$ with the cluster size $M$ as a power law: $R_g=M^{1/D_f}$. For these multicellular clusters, the number of nuclei per cluster scaled with the radius of gyration as $D_f=1.74\pm 0.03$ (95% confidence interval for power-law fit, ordinary nonparametric bootstrap analysis) for all clusters of at least 4 nuclei over varying initial densities and early times (Fig. 2E). This fractal dimension was also evaluated independently using a box-counting analysis of the cell morphology over varying length scales, which yielded a similar $D_f=1.65\pm 0.04$ (mean $\pm$ standard deviation [SD]) (SI Appendix, Fig. S2D). It should be noted that the first analysis of fractal dimension is based on the discrete spatial distribution of cell nuclei, whereas the second analysis is based on the continuous spatial distribution of cell cytoplasm. The slight differences may be explained by the elongated cell morphologies at the peripheral branches relative to nuclei positions. Remarkably, $D_f$ = 1.7 is consistent with previous experimental and computational work on diffusion-limited aggregation of nonliving colloids, which exhibit random diffusive motion but adhere irreversibly on contact (23–25). Thus, this motility-driven aggregation of living cells into stationary fractal-like clusters exhibits striking similarities to the diffusion-limited aggregation of colloids.

Transient Collective Migration Connects Clusters into Spanning Networks.

Within stationary clusters, small groups of cells often exhibited transient but highly correlated motions (Fig. 3A). This “dynamic heterogeneity” was elucidated using the 4-point susceptibility ${\chi }_4\left(\tau \right)$, which resolves correlated density fluctuations in space and time (44). Typically, ${\chi }_4\left(\tau \right)$ exhibits a peak that corresponds approximately to the characteristic size and timescale of collective motion. For this analysis, cells were classified as “motile” with 0% overlap if they traveled more than 1 nuclear diameter (10 $\mathrm{\mu }$m) over some time interval $\tau$, analogous to “mobile” nonliving particles (45). This cutoff distance is comparable to the 15% of a cell diameter used elsewhere for jamming in epithelial monolayers (36). The self-overlap function $⟨Q\left(\tau \right)⟩$ was calculated by averaging these overlap values for all cells within a time frame and then ensemble averaging over a range of start times (45). Essentially, $⟨Q\left(\tau \right)⟩$ represents the fraction of cells that have displaced less than 1 nuclear diameter after some time duration $\tau$ has elapsed. For instance, $⟨Q\left(\tau \right)⟩=1$ at $\tau =0$ h but decreased to $⟨Q\left(\tau \right)⟩=0.5$ at $\tau =2$ h at early times ($<12$ h), indicating that roughly half of the cells were motile (nonoverlapping) after 2 h (Fig. 3B). At later times ($>12$ h), $⟨Q\left(\tau \right)⟩$ increased, corresponding to a decreasing fraction of motile cells and an overall arrest of cell migration over time. Next, to account for temporal correlations between motile cells, the 4-point susceptibility ${\chi }_4\left(\tau \right)=N\left[⟨Q{\left(\tau \right)}^2⟩-{⟨Q\left(\tau \right)⟩}^2\right]$ was calculated from the moments of $⟨Q\left(\tau \right)⟩$ and the total number of cells $N$. At early times ($<12$ h), the peak ${\chi }_4^{*}\left(\tau \right)$ occurred at $\tau \sim$ 1.5 h, which corresponds to a characteristic lifetime of correlated motion (Fig. 3B). At later times, ${\chi }_4^{*}\left(\tau \right)$ shifted slightly to $\tau \sim$ 2 h. The peak height ${\chi }_4^{*}\left(\tau \right)$ decreased from 6 cells to 2 cells, which can be attributed to the increase in cell numbers over time.

Some highly motile cells ($>$7 $\mathrm{\mu }$m/h) often exhibited an elongated morphology with directionally persistent motion, reminiscent of a leader cell phenotype (Fig. 3C) (8, 9, 14). Indeed, these leader-like cells exhibited highly concerted motion with followers (Fig. 3A) and retained some cell–cell contacts. Immunofluorescence staining at the completion of the experiment revealed E-cadherin localization at cell–cell junctions, which was absent in growth media (Fig. 3C). Cells at the cluster periphery also exhibited increased vimentin expression, a classical mesenchymal biomarker (46). These highly motile cells therefore represent a leader cell phenotype that retains both epithelial and mesenchymal biomarker expression, consistent with a partial EMT.

Over time, leader cells from adjacent clusters often guided collectively migrating strands to merge together into a continuous multicellular structure. Once cell–cell contact occurred, leader cells lost their motile phenotype with front–back polarization and reverted to an adherent epithelial morphology with cell–cell junctions, consistent with previous reports (9). As a consequence, isolated clusters became connected together as spanning, space-filling networks, analogous to the gelation of nonliving particles (26). The motion of these multicellular strands away from the cluster was highly persistent, which could be directly visualized using a kymograph analysis (Fig. 3D and SI Appendix, Fig. S3 A and B). The typical persistence time of these strands based on lamellipodial extension was $\sim$4 h, which greatly exceeded the persistence time of 2.5 h for the largely random motion of other cells within the cluster (Fig. 3E). Indeed, these transient bouts of leader-driven collective migration often reached lamellipodial speeds of 10 $\mathrm{\mu }$m/h or more, which greatly exceeded the population-averaged speed of $<$5 $\mathrm{\mu }$m/h at later times, when the migration of most cells was largely arrested (Fig. 3F). Interestingly, the number of leader cells at the cluster periphery increased with increasing cluster size. For instance, smaller clusters of less than 10 cells were typically associated with only 1 leader cell (Fig. 3G). Moreover, larger clusters with 10 to 30 cells were associated with 2 or more leader cells, etc. The number of leader cells scales approximately linearly with the radius of gyration, with an average separation of $\sim$245 $\mathrm{\mu }$m between leader cells along the periphery (SI Appendix, Fig. S3 C and D), in agreement with previous theoretical predictions (10, 15). Overall, the emergence of leader cells at the cluster periphery resulted in transient collective motion that connected isolated clusters into a continuous, space-filling architecture.

Clustering and Jamming Are Governed by Local Cell Density and EGF.

Cells were observed to accumulate greater numbers of neighbors over time during the transitions from dispersed individuals to aggregated clusters. As a measure of local cell density, the number of neighbors (“bonds”) $⟨B\left(t\right)⟩$ was determined by counting all nuclei located within 75 $\mathrm{\mu }$m of a given cell and then ensemble averaging over all cells within a field of view. For the lowest initial density condition ($\sim$30 cells/${\mathrm{m}\mathrm{m}}^2$), the transition from dispersed individuals to aggregated clusters occurred at $⟨B⟩\hspace{0.17em}\sim \hspace{0.17em}4$ when $t\sim$ 30 h (Fig. 4A). Similarly, for a higher-density condition ($\sim$ 60 cells/${\mathrm{m}\mathrm{m}}^2$), this transition to aggregated clusters occurred at $⟨B⟩\hspace{0.17em}\sim \hspace{0.17em}4$ when $t\hspace{0.17em}\sim \hspace{0.17em}20$ h, while the subsequent formation of spanning networks occurred at $⟨B⟩\sim$ 7 when $t\sim$ 40 h. This transition has analogies to a percolation threshold associated with gelation (47), although it exceeds the critical bond number associated with an idealized square lattice ($\approx$4), due to increased connectivity since clusters are linked by “strands” that are several cells wide. In general, higher-density conditions formed clusters even more rapidly ($<$20 h), which connected together soon afterward (by 30 h). It should be noted that cells typically packed together more closely than 75 $\mathrm{\mu }$m within epithelial monolayers, so that this analysis likely overestimates the number of cells in direct proximity at higher densities. Overall, $⟨B\left(t\right)⟩$ captures the transitions to clusters and spanning networks across conditions.

The experimental regime where clustering and arrested motility occurred was characterized by systematically varying the initial EGF concentration from 0 ng/mL up to 20 ng/mL in growth media. For low [EGF] $<0.5$ ng/mL, cells were observed to form clusters ($⟨B⟩<$ 4) and spanning networks ($⟨B⟩<$ 7) (Fig. 4B), with proliferation arrested at subconfluent cell densities. Moreover, for [EGF] $\le$ 0.075 ng/mL, cell motility was largely arrested for both clusters and spanning networks ($<$7 $\mathrm{\mu }$m/h), while at [EGF] = 0.1 ng/mL, cell motility was arrested only in spanning networks. At high [EGF] $\ge$ 0.5 ng/mL and above, cells typically migrated much more rapidly and did not form stable clusters (Fig. 4C), but proliferated exponentially (SI Appendix, Fig. S4 A and B). Motility was consistently arrested at $⟨B⟩\hspace{0.17em}\sim \hspace{0.17em}13$, which correlated with the formation of a confluent, space-filling monolayer. It should be noted that these EMT-induced cells showed weak E-cadherin expression at cell–cell junctions for high concentrations of EGF, suggesting that this slowdown in motility may occur through transient cell–cell contacts, consistent with contact inhibition of locomotion through active repulsion or random deflection (48). These results indicate that increasing EGF enhances both motility and proliferation and that cluster formation and arrested motility occur for [EGF] $\le$ 0.1 ng/mL and $⟨B⟩\hspace{0.17em}<\hspace{0.17em}13$.

These transitions between motile individuals and arrested clusters were recapitulated by a minimal model of self-propelled particles that incorporated 3 mechanisms. First, each particle moved randomly with a polarity force of constant magnitude, which changed to a randomly chosen direction at constant intervals (offset to different starting times) (SI Appendix, Fig. S4E). Second, pairs of particles adhered via a short-ranged attraction, but could not approach each other closer than a fixed particle size (SI Appendix, Fig. S4F). This simplified interaction was inspired by theoretical models of diffusion-limited aggregation in nonliving colloidal particles with strong adhesion (47) and neglects more complex friction-like biomolecular interactions that could occur at cell–cell or cell–matrix adhesions (35). Third, cells proliferated with a cell cycle of constant duration (offset to different starting times). However, cells also obey contact inhibition of proliferation and could not divide if they had 4 or more neighbors. For a representative simulation with fixed adhesion ($A$ = 0.03) and varying polarization (P = 0.0025 to 0.0043), the initially dispersed individuals transitioned over time to multicellular clusters and then to a spanning network (Fig. 4C and Movie S3). These values further corresponded to slower proliferation, decreasing fraction of individuals, decreasing number of clusters, increasing cluster size, and a fractal dimension of $D_f=1.74\pm 0.003$ (SI Appendix, Fig. S4 G–K), in good agreement with experiments (Figs. 1C and 2 B–E). The corresponding increase in the ensemble-averaged number of neighbors $⟨B⟩$ associated with clusters ($>$2) and spanning networks ($>$4) was qualitatively consistent with experimental observations (Fig. 4A). This minimal model was also modified to address leader cell formation based on a local curvature-dependent mechanism that sustained cell polarization (SI Appendix, Fig. S3 E–H), analogous to previous work (15). For simplicity, proliferation was controlled only by contact inhibition, so that particles continued to proliferate until a spanning network was formed. In comparison, some experiments with lower initial cell densities showed minimal proliferation of clusters, indicating that proliferation may exhibit a more complicated time dependence (Fig. 4A). Future work will systematically explore how time-varying proliferation rates affect clustering in this minimal model.

More generally, particles in this regime exhibited a gradual arrest of cell motility with increasing ensemble-averaged number of neighbors $⟨B⟩$ (Fig. 4D and SI Appendix, Fig. S4 L and M), in qualitative agreement with the low-EGF regime (Fig. 4B). A second regime with slightly decreased adhesion ($A$ = 0.02) resulted in a transition from subconfluent individuals to a confluent monolayer with corresponding arrest of motility (SI Appendix, Fig. S4 N and O and Movie S4), consistent with the observed behaviors in the higher-EGF regime (Fig. 4B). It should be noted that the crossover between these 2 regimes could not be captured by a single adhesion value in our model (SI Appendix, Fig. S4F) and likely occurred because the short-ranged adhesion was very sensitive to changing polarization value ($A$ = 0.03). This theoretical result suggests that experimentally increasing EGF may both increase cell motility and decrease cell–cell adhesion, maintaining cells as individuals rather than clusters. Finally, it should be noted that this model is based on monodisperse particles, which cluster together with uniform separations. However, living cells have deformable shape and varying sizes, which results in varied nuclear separations. We have shown that a 75-$\mathrm{\mu }$m threshold is optimal for linking cells into clusters, but this likely overcounts the number of nearest neighbors relative to the minimal model.

Clustering Occurs with Decreased EGF Signaling for Growth-Factor–Dependent Mammary Epithelial Cells.

These clustering behaviors were most pronounced with OHT treatment, but qualitatively similar behaviors were observed under a number of experimental conditions where EGF signaling was reduced. For instance, clustering also occurred with MCF-10A cells without EMT induction, when cultured in assay media with a matched concentration of dimethyl sulfoxide (DMSO) (corresponding to the concentration needed to dissolve OHT) (SI Appendix, Fig. S5 and Movies S5 and S6). Indeed, considerably fewer leader cells were observed without EMT induction (DMSO), resulting in clusters that were more compact with reduced connectivity at longer times (SI Appendix, Figs. S6 and S7) relative to conditions with EMT induction (OHT). Moreover, clustering also occurred after drug treatment to inhibit the EGF receptor (i.e., gefitinib, in growth media) (SI Appendix, Fig. S8 and Movies S7 and S8). Cell viability was minimally affected by culture in assay media and gefitinib relative to growth media (SI Appendix, Figs. S5 C and D and S8 B and G). It should be noted that the proliferation and migration of MCF-10A cells are highly dependent on growth factors (e.g., EGF) (49), unlike transformed breast cancer cell lines, which are less sensitive to growth factor concentration. For a meaningful comparison, a second growth factor-dependent mammary epithelial cell line (hTERT-HME1) also exhibited clustering in reduced growth factor conditions across a range of starting densities and after gefitinib treatment (SI Appendix, Fig. S9). In contrast, clustering was not observed with the highly metastatic breast adenocarcinoma cell line, MDA-MB-231, which is transformed and exhibits reduced sensitivity to EGF (SI Appendix, Fig. S10). These experiments suggest that clustering occurs only for growth factor-dependent mammary epithelial cells and is driven by decreased EGFR signaling, through either reduced EGF or chemical inhibition.

Cell Culture.

MCF-10A mammary epithelial cells stably transfected with an inducible Snail expression construct fused to an estrogen receptor response element were a gift from D. A. Haber (Massachusetts General Hospital, Boston, MA) (41). One variant overexpressed fluorescent proteins in the nucleus (mCherry-H2B) and cytoplasm (GFP) for live cell tracking. Cells were routinely cultured as described in SI Appendix. For these experiments, cells were cultured on collagen-coated 96-well plates at varying EGF concentrations (0 to 20 ng/mL) with DMSO and OHT, with a total media volume of 125 $\mathrm{\mu }$L. No changes in media color were observed after 60 h. Cells tested negative for Mycoplasma contamination (PlasmoTest; Invivogen) and an absence of additional cell culture contaminants was routinely verified by culturing cells in antibiotic-free media. At least 3 independent experiments were conducted for each condition.

Time-Lapse Microscopy.

Cell proliferation, clustering, and migration were measured with an inverted epifluorescence microscope (Nikon TiE) with a light-guide coupled white light illumination system (Lumencore Spectra-X3) under environmentally controlled conditions ($\mathrm{37}{}^{○}$C, 5% ${\mathrm{C}\mathrm{O}}_2$, humidified). Using Nikon Elements software, images were acquired every 15 min with 12-bit resolution using an sCMOS camera (Andor Neo), a 10$\times$ Plan Apo objective (NA 0.45, 4 mm working distance), a GFP/FITC Filter Set (Chroma 49002), and a TRITC/DSRed Filter Set (Chroma 49004). Images were recorded under consistent acquisition parameters (e.g., exposure time, camera gain/gamma control, and microscope aperture).

Image Analysis and Cell Tracking.

Automated cell tracking based on fluorescent nuclei was performed using Bitplane Imaris 8.2. First, nuclei were detected based on an estimated $x-y$ diameter of 11.5 $\mathrm{\mu }$m and then linked together across time points using an autoregressive motion algorithm, with a maximum displacement of 30.0 $\mathrm{\mu }$m, a max gap size of 3 frames, and fill gaps enabled. Consistency in cell phenotype (e.g., confluency, clustered/nonclustered, morphology, network architecture) was first confirmed qualitatively across all experiments. Subsequently, 3 representative experiments were chosen for the bulk of semiautomated analyses (single-cell tracking and manual verification) and 2 representative experiments were chosen for labor-intensive manual analyses (cluster threshold determination and leader cell identification).

Proliferation and Velocity Analyses.

Cell counts were determined from detected nuclei in each time-lapse image. The normalized cell count was calculated by dividing cell counts at all times by the starting cell count.

Single-cell velocities were calculated based on the displacement of nuclei over some time interval, which was averaged over 1 h (4 time frames) to reduce noise from nuclear shape changes. Thus, the velocity ${\mathbf{v}}_i\left(t_j\right)=\left[{\mathbf{r}}_i\left(t_j\right)-{\mathbf{r}}_i\left(t_{j-4}\right)\right]/\left(t_j-t_{j-4}\right)$, where ${\mathbf{r}}_i$ was the position of cell $i$ at time frame $j$ (for $j>4$); cells tracked for fewer than 5 frames were necessarily discarded. The population-averaged velocity was calculated by averaging over all single-cell velocities at a given time frame. Significance was computed at the $P<0.05$ level using a paired-sample $t$ test between conditions with $n=3$ and was reported hourly for simplicity. The time interval used to calculate velocity was also varied (15 min, 30 min, and 2 h compared to 1 h) to show that trends between conditions were robust.

Mean-Squared Displacement, Cluster Mass, and Fractal Dimension.

The mean-squared displacements were determined as a function of time interval $\tau$ as $MSD\left(\tau \right)=⟨{|{\mathbf{r}}_i\left(t+\tau \right)-{\mathbf{r}}_i\left(t\right)|}^2⟩$, averaged over all early starting times $t$ and computed separately for each individual cell $i$, using the open-source MATLAB code MSDANALYZER (52). To be considered in the analysis, a cell must remain separated from all neighboring cells by greater than 75 $\mathrm{\mu }$m for at least 5 h over early times (0 h $<$ t $<$ 20 h) and subsequently $MSD\left(\tau \right)$ was computed for $\tau >$ 1 h to eliminate noise from nuclear motion and $\tau <$ 10 h since individual cells quickly associated into clusters. The resulting curve was fitted according to a power law of the form $K_{\alpha }{\tau }^{\alpha }$ using the MATLAB function FIT with option POWER1. The data were plotted alongside a line of $\alpha$ = 1 and $\alpha$ = 1.5 to guide the eye. The mean and SD of fitted parameter $\alpha$ were quantified for the resulting curves of individuals with high goodness of fit (R² $>$ 0.90) across 3 independent experiments (n = 44). The generalized diffusion coefficient $K_{\alpha }$ was plotted against exponent $\alpha$ for all individuals meeting these criteria to show the spread of the data.

Multicellular clusters were defined at each time frame as groups of 2 or more cells where neighboring pairs of nuclei were separated by 75 $\mathrm{\mu }$m or less; anything not meeting these criteria was defined as a single cell. The value 75 $\mathrm{\mu }$m was chosen as the optimal nuclear linking distance based on visual inspection of linked clusters across a range of threshold distances relative to the actual cytoplasmic connectivity (detailed in SI Appendix). The fraction of individuals was computed across the time lapse by taking the quotient of the number of single cells and the total number of cells. Cluster size was calculated by dividing the total number of cells in clusters by the total number of clusters.

For the purpose of computing the fractal dimension, multicellular clusters were defined at each time frame as groups of 4 or more cells where successive pairs of nuclei were separated by 75 $\mathrm{\mu }$m or less, to reduce the noise associated with especially small clusters. For each cluster with a “mass” of $M$ cells, the centroid was calculated as ${\mathbf{r}}_{cm}=\left(1/M\right){\sum }_i^M{\mathbf{r}}_i$. The radius of gyration $R_g$ was then calculated from $R_g^2=\left(1/M\right){\sum }_i^M{\left({\mathbf{r}}_i-{\mathbf{r}}_{cm}\right)}^2$. The scaling of the radius of gyration $R_g$ vs. $M$ as a power-law fit with exponent $1/D_f$ was first calculated from this dataset using the MATLAB function FIT with option POWER1. Next, an ordinary, nonparametric bootstrap analysis was conducted to determine confidence intervals for this power-law scaling using the R function BOOT. Briefly, the dataset was randomly resampled 1,000 times (with replacement runs) and each sample was reanalyzed to empirically determine a 95% confidence interval for the power-law fit.

Clustering was independently verified by direct segmentation of cell cytoplasm and nuclei using custom MATLAB code. First, cell nuclei were segmented (by adaptive thresholding) and connected using mathematical morphology to obtain a rough foreground binary mask. Next, MATLAB’s implementation of the Chan–Vese algorithm was used to obtain segmented cell boundaries. The number of nuclei inside each cluster was estimated by dividing the total area of the nuclear segmentation by typical nucleus size. Finally morphological features, including fractal dimension, were calculated for each segmented cluster. The segmentation quality was manually evaluated and only data obtained from correctly segmented images were used in subsequent analysis.

Dynamic 4-Point Susceptibility Function.

Spatially heterogeneous dynamics were quantified based on the self-overlap of cells, which defined motile or nonmotile subpopulations. A cutoff distance of $L=10\hspace{0.17em}\mathrm{\mu }$m was selected, which corresponds to 1 nuclear diameter and roughly 20% of a cell diameter, comparable to previous analyses on epithelial monolayers (36). The cutoff function $w_i$ is 1 if the displacement of cell $i$ is less than $L$ over some time interval $\tau$ (nonmotile) and 0 otherwise (motile). An instantaneous self-overlap parameter is given by $Q\left(L,\tau \right)=\frac{1}{N}{\sum }_{i=1}^N{w}_i$. A 4-point susceptibility function can be calculated from the moments of $Q\left(L,t\right)$ averaged over all starting times and all cells: ${\chi }_4\left(\tau \right)=⟨N⟩\left[⟨Q{\left(L,\tau \right)}^2⟩-{⟨Q\left(L,\tau \right)⟩}^2\right]$ (45), where $⟨N⟩$ is the average number of cells present over the time range indicated.

Leader Cell Analysis.

Leader cells were manually identified by 3 individuals (S.E.L., Z.J.N., J.Y.S.) based on the following criteria: 1) Leader cells were located at the periphery of clusters, 2) leader cells drove collective migration of follower cells (while maintaining physical contact) over at least several hours, 3) leader cells exhibited elongation in the direction of cell migration, and 4) leader cells tended to present ruffled lamellipodia and an enlarged cell size relative to follower cells. The combined number of leaders identified across all individuals, or the total number of unique leader cells, was then used to determine the threshold at which multiple leaders emerge and the dependency on cluster size. Kymographs were constructed using Icy (Institut Pasteur, 2011) and used to determine net leader speed (start to end position of leading edge over time range of interest) and persistence time (the duration over which cells traveled along the same direction).

Immunofluorescence Staining.

Immunofluorescence staining protocols are described in SI Appendix. Briefly, cells were fixed and labeled with E-cadherin primary antibody (mouse mAB; BD 610181) and vimentin primary antibody (rabbit mAB; CST 5741) and then incubated with Hoechst (ThermoFisher H3569), goat anti-mouse, and goat anti-rabbit secondary antibodies (Alexa Fluor 488 and 555; ThermoFisher). Cells were then washed several times with 1$\times$ PBS prior to imaging.

Local Neighbor Density and Jamming Phase Diagram.

The ensemble-averaged local neighbor density was calculated for each cell at each time point by finding the number of neighboring cells within a 75-$\mathrm{\mu }$m radius and then averaging over all cells detected for that time point. The jamming phase diagram was then mapped by comparing this local neighbor density with the ensemble averaged velocity $⟨v_{rms}⟩$ (as previously calculated), over a range of experimental EGF concentrations.

Computational Self-Propelled Particle Model.

A minimal physical model treated cells as self-propelled particles which randomly polarize in new directions and can interact with other cells through a tunable attractive interaction and a short-range repulsion corresponding to the cell radius. Particle dynamics were simulated in a square domain with periodic boundary conditions. Further details are available in SI Appendix.

Statistical Analysis.

Experiments were repeated at least 3 times, with the number of independent experiments used for analysis indicated in the figure legends. For major findings of population-averaged measurements over time (e.g., normalized cell count [Fig. 1C and SI Appendix, Figs. S5E, S8C, S8H ], cell speed [Fig. 1D and SI Appendix, Figs. S5F, S8D, S8I ], fraction of individuals [Fig. 2B and SI Appendix, Fig. S6A ], number of clusters [Fig. 2C and SI Appendix, Fig. S6B ], and cluster size [Fig. 2D and SI Appendix, Fig. S6C ]), the mean and SD were computed across 3 biological replicates from independent experiments. Statistical significance between groups was determined and plotted at hourly intervals, where applicable, at the 5% significance level (P values $<$0.05) using the paired-sample t test. For significance testing between distributions of individual cells from manual analyses (e.g., leader/follower cells [Fig. 3 E and F and SI Appendix, Fig. S7 D–F ]), data from 2 representative experiments were pooled and results from the 2-sample Kolmogorov–Smirnov test were reported (5% significance level). For a subset of findings, representative data from a single experiment were shown for ease of visualization (4-point susceptibility function [Fig. 3B and SI Appendix, Fig. S7A ], Ripley’s H function [SI Appendix, Figs. S1E, S5G, S8E, S8J ]). In these cases, data were rigorously analyzed over several replicates and across at least 2 experiments, to ensure consistency in the observed trends.