Combined heterogeneity in cell size and deformability promotes cancer invasiveness

Tumors are comprised of phenotypically distinct subpopulations, with heterogeneity associated with higher tumor growth and poor therapeutic efficiency (Lawson et al., 2018; Arozarena and Wellbrock, 2019). Several factors contribute to phenotypic heterogeneity, including genetic factors such as genomic instability (Burrell et al., 2013), non-genetic factors such as cellular plasticity (Brock et al., 2009), and by microenvironmental cues such as nutrient availability. In addition, during epithelial-mesenchymal transition (EMT), which plays a central role in cancer metastasis, epithelial cells rarely undergo a complete transition to a mesenchymal state, thereby giving rise to phenotypically distinct subpopulations (Pastushenko et al., 2018; Pastushenko and Blanpain, 2019). EMT also leads to the generation of tumor-initiating cells, also referred to as cancer stem cells (CSCs) (Mani et al., 2008). Latest single-cell technologies have proven to be very useful in assessing and implicating molecular-level heterogeneity in tumor growth and relapse (Gan et al., 2018; Patel et al., 2014). The co-existence of multiple stable subpopulations within the same tumor may be indicative of cooperation between individual subpopulations (Tabassum and Polyak, 2015). Indeed, cooperation between phenotypically distinct subpopulations has been shown to play a beneficial role in tumor growth, invasion and drug resistance (Zhou et al., 2017). In the context of invasion, such cooperation between poorly invasive melanoma cells and proteolytically active melanoma cells drives collective invasion (Chapman et al., 2014). Similarly, in luminal breast cancer cells, lesser invasive CD44^low follower cells invade collectively with CD44^hi leader cells (Yang et al., 2019).

Cell size and cell deformability represent two key biophysical attributes of cells that have been shown to modulate cancer invasion and metastasis (Wolf et al., 2013; Kumar et al., 2018b). Intuitively, smaller cell size is expected to provide an invasion advantage not only during migration through the stroma, but also during transmigration through the endothelium. Studies have established small size as a key feature of stem cells as well as CSCs (De Paiva et al., 2006; Lia et al., 2015). In CSCs, the combination of small size and phenotypic plasticity has been shown to enable fast migration in aligned collagen matrices (Ray et al., 2017). In addition to smaller cell size, increased cell deformability has been shown to promote invasiveness in both cancer cell lines and patient-derived cells (Swaminathan et al., 2011; Xu et al., 2012; Harada et al., 2014), and is one of the key features of tumor-initiating CD44^hi CD24^low breast CSCs (Al-Hajj et al., 2003; Sheridan et al., 2006) and drug-resistant ovarian cancer cells (Kapoor et al., 2018), as well as cells migrating in non-proteolytic fashion (Das et al., 2019). Strikingly, studies have also demonstrated that more heterogeneous cellular populations possess a higher propensity for secondary metastasis (Cheung et al., 2016), suggesting that phenotypic heterogeneity may play a key role in mediating metastasis. Nevertheless, the extent of variation in cell size and cell deformability within these heterogeneous cell clusters, and their relative contribution to cancer invasiveness, remains obscure.

The relative importance of cell size and deformability in cancer invasiveness may also depend on the extent of cell-cell adhesion, which determines the mode of migration, i.e. single-cell migration versus collective cell migration. Whereas epithelial cells possessing high levels of E-cadherin exhibit collective migration, mesenchymal cells with low or no E-cadherin exhibit single-cell migration. EMT has long been considered as the key cellular transformation for cancer metastasis (Dongre and Weinberg, 2019). The process of EMT, originally thought to be a binary process, is increasingly understood to represent a spectrum of states with varying expression of epithelial/mesenchymal genes giving rise to the hybrid epithelial/mesenchymal phenotype (Jolly et al., 2018; Liao and Yang, 2020). This hybrid phenotype has been observed in circulating tumor cells (CTCs), which disseminate as tumor cell clusters with intact cell-cell adhesions (Yu et al., 2013; Fabisiewicz and Grzybowska, 2017). In addition to CTCs, the hybrid epithelial/mesenchymal phenotype has been observed in multiple types of cancers, with the hybrid subpopulation tending to localize at the leading edge of primary tumors (Puram et al., 2017). However, the extent to which the biophysical properties of this subpopulation enables its enrichment at the invasive front has not been adequately studied.

In this paper, using MCF-7 and MDA-MB-231 human breast cancer cell lines, which vary in the extent of cell-cell adhesion, we first document the extent of heterogeneity in cell size, deformability and cell motility. We then probe the importance of this biophysical heterogeneity on cancer invasiveness using a computational model of cancer invasion. Our results suggest that combined heterogeneity in cell size and deformability enhances the random scattering of a tumor-mimicking cell cluster independent of the extent of cell-cell adhesion, with comparable invasiveness of subpopulations of different sizes and deformabilities. Although intermixing of different subpopulations is observed even during directed cell migration, our results predict smaller-sized cells of varying deformability as one prominent subpopulation that gets enriched at the invasive front. The phenotypic characteristics of this subpopulation exhibit considerable overlap with the properties of CD44^hi-expressing CSCs in MDA-MB-231 cells. Finally, using spheroid invasion experiments, we demonstrate the temporal enrichment of CSCs at the invasive front.

Phenotypic heterogeneity was characterized by measuring variability in cell-cell adhesion, cell size, cell speed and cell deformability of MCF-7 and MDA-MB-231 breast cancer cells. In comparison to strong E-cadherin localization at cell-cell junctions in MCF-7 cells, E-cadherin levels and peripheral localization were substantially reduced in MDA-MB-231 cells, which are more mesenchymal in nature (Fig. 1A,B). Though the two cell lines exhibited overlapping size distributions, the average size of MDA-MB-231 cells was nearly twice that of MCF-7 cells (Fig. 1Ci). The tracking of the size of individual cells over a 24 h period revealed the presence of cells that exhibited a relatively unchanged size, as well as cells that exhibited a large variation in size (Fig. 1Cii,Ciii; Movie 1). The distribution of maximum change in cell size (with respect to cell size at t=0 h) revealed significantly lesser variation in the size of MDA-MB-231 cells compared to MCF-7 cells (Fig. 1D). MDA-MB-231 cells were nearly twice as motile compared to that of MCF-7 cells (Fig. 1E). However, no correlation between cell size and motility was observed in both the cell types, with the values of the correlation coefficients (r_MCF and r_MDA) being close to zero (Fig. 1F). Further quantification of cell stiffness revealed that MDA-MB-231 cells were marginally softer than MCF7 cells (Fig. 1G). Together, these results suggest that highly invasive MDA-MB-231 cells tend to be larger, softer and more migratory compared to non-invasive MCF-7 cells.

To check whether the observed functional-level heterogeneity is also present at the transcript level, publicly available single-cell RNA-seq datasets from breast cancer samples (GSE75688) (Chung et al., 2017) were analyzed by performing single sample enrichment analysis for stemness, EMT and cell biophysical characteristics (refer to Materials and Methods). Although publicly available genesets were used for calculating stemness and EMT scores, characterization of cell mechanical properties using gene set enrichment analysis (GSEA) was performed using a manually defined BiophysicalGeneset comprising of regulators of cytoskeletal dynamics, including actin, vimentin, actin crosslinking proteins, myosin and Rho GTPases and their downstream effectors (refer to Materials and Methods). Additionally, to gain insights into intratumor and intertumor heterogeneity, enrichment scores were studied for samples from individual patients, as well as for the entire dataset. This analysis indeed revealed a wide heterogeneity in stemness, EMT and BiophysicalGeneset scores from within tumor from the same patient, as well as within the entire dataset (Fig. S1A-D). Next, to study the correlation between biophysical characteristics and stemness/EMT, we compared the enrichment score of BiophysicalGeneset with that of stemness and EMT scores for individual cells. Interestingly, this comparison showed a strong co-enrichment of BiophysicalGeneset and stemness, as well as BiophysicalGeneset and EMT scores (Fig. S1E,F). Though previous studies have shown a strong correlation between stemness and EMT (Chung et al., 2017; Mani et al., 2008), the reported relationship with biophysical characteristics is novel and was not clear from previous studies. In line with the heterogeneity observed in human breast cancer samples, a wide heterogeneity in expression levels of genes regulating cellular biophysical properties was observed in MCF7 and MDA-MB-231 cell lines (Fig. S2) (Jonasson et al., 2019; Birts et al., 2020). Comparatively, less heterogeneity was observed in EMT and stemness markers in cell line datasets. Collectively, these results suggest that breast cancer cells exhibit heterogeneity at the phenotypical level, as well as in gene expression patterns.

To probe the importance of phenotypic heterogeneity in cancer invasion, we used our recently developed cellular Potts model (CPM)-based formalism, which allows us to account for multiple factors, including cell size and deformability, cell-cell adhesion, extracellular matrix (ECM) density and ECM remodeling. The simulation framework consists of three types of pixels representing the cell (pink), ECM fibers (blue lines) and the interstitial fluid (white) (Fig. 2A; see Materials and Methods for details) (Kumar et al., 2016b). At the start of simulations (t=0), the cell aggregate is positioned at the center of the simulation grid. Although ECM fibers sterically hinder cell migration, the lateral invasion of cells is mediated by cell secreted matrix metalloproteinases (MMPs) which degrade the fibers and create space for migration. Simulations were performed for the case of a homogeneous cell aggregate and a heterogeneous cell aggregate. For the heterogeneous aggregate, size (or target cell area) and stiffness (or area constraint) of individual cells within the aggregate were assumed to follow Gaussian distributions (Fig. 2B). Furthermore, cells of varying sizes and stiffness were distributed such that there was no bias in their position within the cluster (Fig. S3A). Cell size and cell stiffness of the homogeneous cell aggregate were set to the mean cell size and mean stiffness of the heterogeneous cell aggregate.

Spatiotemporal evolution of the system is based on random movement of individual pixels subject to transition probabilities using the Monte Carlo method, and is based on the total system energy (Eqn 1). One of the major factors contributing to the total energy of the system is dependent on the nature of adhesive interactions between the three different types of pixels in our simulations (cell-cell, cell-matrix and cell-fluid; first term in Eqn 1; Fig. S3B). Although cell-matrix (J_cm) and cell-fluid (J_cf) adhesion energies were held constant (i.e. J_cm=16, J_cf  =32 in units of K_BT/L), cell-cell adhesion energies (J_cc) were varied to mimic different extents of cell-cell adhesion to model the invasive behavior of both epithelial-like cells, such as MCF-7 cells, as well as mesenchymal-like cells, such as MDA-MB-231 cells. Whereas a choice of J_cc=1 led to collective migration mimicking the behavior of MCF-7 cells, a choice of J_cc=40 led to cell scattering mimicking the behavior of MDA-MB-231 cells (Fig. 2C; Movies 2-4). In comparison, J_cc=16 corresponds to cells with intermediate cell-cell adhesion in which both single cells and cell clusters (a cluster corresponds to a group of cells in contact with each other) were observed. The second term contributing to the total system energy is associated with deviations in the sizes of individual cells (Fig. S3C) from their target sizes (Fig. 2B).

The invasiveness of individual cells was quantified by measuring the total distance (d) traveled by a given cell, as well as the net end-to-end distance (D) traveled by the cell between its initial and final positions. Both d and D were measured based on the position of the cell centroid. Although no difference was observed in the population-averaged total distance traveled by cells (Fig. 2D, inset), independent of J_cc, heterogeneity in cell size and deformability led to an increase in population-averaged net translocation compared to the homogeneous case (Fig. 2D; Fig. S3D). Interestingly, the invasiveness of cells possessing heterogeneity only in size (i.e. stiffness of all cells = ) or heterogeneity only in stiffness (i.e. size of all cells =   ) were less invasive compared to cells possessing combined heterogeneity in cell size and deformability (Fig. S3E).

To test whether higher of heterogeneous population is due to increased motility of one particular subpopulation, the population was divided into five subgroups depending upon their size and stiffness. These included cells of intermediate size and deformability (ICSD), small and soft cells, small and stiff cells, large and soft cells and large and stiff cells. ICSD cells correspond to cells possessing size and stiffness within 1 s.d. of the mean cell size and mean deformability of the population (see Materials and Methods). Whereas ‘small cells’/‘large cells’ correspond to cells smaller/larger than the size of ICSD cells, ‘soft cells’/‘stiff cells’ correspond to cells softer/stiffer than ICSD cells. In spite of considerable variation in the magnitude of translocation within each subpopulation, negligible differences were observed in the average translocation of the five subpopulations (Fig. 2E). As these observations were based on population averages, invasiveness behavior and its correlation with cell/size and deformability at a single-cell level remains obscure. Studying these aspects at the single-cell level becomes more important in this study due to intrapopulation heterogeneity. Therefore, to further probe the association between higher invasiveness and cell size/deformability at the single cell level, t-distributed stochastic neighbor embedding (tSNE)-based two-dimensional-embedding of all cells (from 22 independent simulations) was performed (Fig. 2F, Materials and Methods). Overlaying these two-dimensional projections with cell size, d, D and λ for individual cells revealed that high translocation (D>400 μm) was exhibited by cells of varying sizes and deformabilities (Fig. 2G, red regions). Even though the majority of the cells that showed higher d were smaller in size (Fig. 2G, blue regions), no correlation was observed with D, suggesting that decrease in cell size may increase motility but not invasiveness. Collectively, these in silico studies predict that heterogeneity in cell size and deformability can promote cancer invasiveness independent of cell-cell adhesion.

Though the above simulations revealed a synergistic effect of heterogeneity in cell size and deformability on invasiveness, these results may vary in the presence of a directional cue, which is relevant in vivo. Therefore, simulations were performed wherein directional migration of the heterogeneous cell aggregate was tracked in the presence of a stable linear chemical gradient (Fig. 3A). The efficiency of chemotaxis was modulated by varying the strength of chemotaxis (μ; see Materials and Methods). As effective chemotaxis was observed for μ>1000 (Fig. S4), simulations were performed for μ=2000 and μ=5000, with an increase in μ leading to faster directed migration (Movies 5,6) (Kumar et al., 2018b). Under these conditions, cells were found to invade in a collective manner independent of J_cc, with negligible scattering observed even at J_cc=40 (Fig. 3B). Quantification of ⁠, the time at which the first cell touches the right boundary of the lattice in a given simulation, was found to decrease with J_cc (Fig. 3C). Comparison of and at t_min, i.e., lowest across all 22 simulations, revealed the highest translocation at intermediate cell-cell adhesion (J_cc=16) and highest persistence (<D f d>) for J_cc=16 and J_cc=40 (Fig. 3D). In line with this, cellular level comparison for μ=5000 revealed that cells that moved highest within t_min (i.e. highest ⁠), and also showed highest translocation, exhibited a wide variation in cell size and were of mixed deformabilities (Fig. 3E, red region). Additionally, at moderate and low cell-cell adhesion (i.e. J_CC ∈{16,40}), there was a subpopulation that exhibited very high D but moderate d (blue region); a major fraction of cells in this subpopulation were large in size (blue arrows) and of mixed deformabilities. Simulations performed with lower chemotaxis strength (i.e. μ=2000) showed a similar trend to that of μ=5000. Together, these analyses implicate intermediate cell-cell adhesion in promoting invasion efficiency of cells in the presence of directional cues.

To further understand how heterogeneity in cell size and deformability imparts high invasiveness to invading populations, the distance from the right boundary at ⁠, i.e. δ, was compared across cells of varying sizes and deformabilities (Fig. 4A). For μ=5000, population averaged δ or was comparable for all the three J_cc values. However, for μ=2000, was lowest for J_cc=1, suggesting that at moderate chemotactic strengths high cell-cell adhesion enables collective invasion of the whole population (Fig. 4B). As before, invasiveness of the population was highest, i.e., was lowest, for the case of combined heterogeneity in cell size and deformability instead of heterogeneity in either cell size or in deformability (Fig. S5A).

To probe whether invading fronts were enriched in cells of particular sizes/deformabilities, the distribution of cell size and deformability was determined for cells that reached within 50 μm from the right edge, i.e., δ≤50 μm (Fig. 4C). Although 40–50% of cells at the invasive front were of intermediate size and deformability, enrichment analysis, i.e. comparison of the proportion of cells of a given subpopulation at the invasive front in comparison to that in the entire population, provided interesting insight. Specifically, small cells of varying deformabilities, small and intermediate-sized soft cells, and intermediately stiff cells of varying cell sizes represent three heterogeneous subpopulations that were enriched at one or more values of J_cc (Fig. 4D).

To probe the possible mechanism(s) mediating enrichment, for all cells within the invasive front, the positions of cells at were plotted as a function of their position at t=50 Monte Carlo Step (MCS) with different colors and shapes depicting different subpopulations (Fig. 4E,F). As the cells were positioned in a rectangular grid-like fashion at t=0 MCS (Fig. S3A), a time point of t=50 MCS was chosen to allow for some initial intermixing of the cells. For J_cc=1, cells that ended up at the invasive front at ⁠, were situated at distances δ_i (distance of i^th cell from right boundary) ranging from 400 to 600 µm at t=50 MCS. The distance 400 µm corresponds to positions at the front of the cell aggregate closest to the right boundary, whereas 600 µm corresponds to positions at the back of the cell aggregate farthest from the right boundary. In comparison, the initial distance ranged from 400 to 500 µm for J_cc=16 and J_cc=40. Consistent with these observations, the average number of neighbor changes of individual cells was highest in the case of J_cc=1 and lowest in the case of J_cc=40 (Fig. S5B). In combination, these results suggest that the initial position of cells and intermixing represent two mechanisms mediating enrichment, with the initial position playing a more dominant role at intermediate and low cell-cell adhesion, and intermixing playing a prominent role when cell-cell adhesions are intact. Interestingly, for J_cc=1, where enrichment of cells from the rear of the cell aggregate was observed, the common feature of these cells was that they were either soft or of intermediate stiffness (Fig. 4F).

The above analysis has helped us identify small cells of varying deformabilities, soft cells of varying cell sizes, and intermediately stiff cells of varying cell sizes, as three heterogeneous subpopulations that are enriched at the invasive front depending on J_cc. To further probe which of these subpopulations is maximally invasive, simulations were performed using heterogeneous populations (A, B and C), wherein cell size/deformability of these populations were drawn from uniform distributions covering the same range of cell size/deformability as that of the above three subpopulations. In addition, a population D was generated by combining populations A, B and C (Fig. 5A). Although values were comparable for the four populations, for all values of J_cc, lowest was observed for population A, i.e. small cells of varying deformabilities (Fig. 5B; Movie 7). However, within this population, no dependence of δ on deformability was observed (Fig. 5C).

To compare and contrast the biophysical properties of population A with that of CSCs, a biophysical characterization of CD44^hi CD24⁻ mesenchymal CSCs (herein referred to as mCSCs) and CD44^hi CD24⁺  hybrid CSCs (herein referred to as hCSCs) was undertaken (Fig. 5D,E; Fig. S6A). MDA-MB-231 cells were found to contain ∼60% mCSCs and ∼6% hCSCs (Fig. S6A). Although mCSCs were smaller and softer compared to the parental population, hCSCs exhibited size and stiffness values comparable to the parental population (Fig. 5F,G). In addition to highlighting differences between the biophysical properties of mCSCs and hCSCs, our results are indicative of an overlap in the biophysical properties of the maximally invasive population A predicted by our simulations with that of CSCs.

To experimentally determine the spatial position of CSCs during invasion, MDA-MB-231 spheroids were implanted in three-dimensional collagen gels and the spatiotemporal positioning of CSCs was assessed. Although the proportion of mCSCs within the spheroids was nearly identical to that in MDA-MB-231 cells cultured in two-dimensional collagen-coated substrates, the proportion of hCSCs was higher (Fig. S6B,C). The temporal pattern of cell scattering was quantified 1, 3 and 5 days after implantation (day 0) (Fig. 6A). Spheroids were implanted in 1.5 mg/ml collagen gels with an average pore size of ∼2 μm (Fig. 6B). To probe the initial spatial distribution of CSCs and non-CSCs (nCSCs) within the spheroid, cells were stained for CD44, F-actin and DAPI (Fig. 6C). Quantification of integrated CD44 and F-actin intensity as a function of radial position from the spheroid center revealed a uniform distribution, suggesting that CSCs (i.e. CD44^hi cells) were randomly positioned within the spheroid (Fig. 6D). Temporal tracking of spheroids revealed an increase in spheroid size, as well as an outward scattering of cells from the spheroid core (Fig. 6E).

To probe for the possible enrichment of CSCs at the invasive front, the spheroid core of day 5 spheroids were divided into seven zones, and the average integrated CD44 intensity was computed for each zone (Fig. 6F). Zones 0-2 were chosen such that they spanned the area occupied by day 0 spheroids. Plots of integrated CD44, intensity normalized with respect to the average CD44 intensity in the innermost zone (i.e. zone 0) at day 1, day 3 and day 5, revealed time-dependent enrichment of CD44^hi cells at the outermost zone (Fig. 6G). A threefold enrichment in the outermost zone (i.e. zone 4) at day 1 increased to fivefold enrichment at day 3 in zone 5, and remained unchanged thereafter at day 5 in the outermost zone 6.

To probe the mechanism underlying the enrichment of CD44^hi cells at the invasive front, interconversion experiments were performed wherein mCSCs, hCSCs and nCSCs (CD44^−/low cells) were isolated from parental MDA-MB-231 cells by fluorescence-activated cell sorting (FACS) and individually expanded up to a duration of 48 h (Fig. 7A). FACS analysis of the individually grown cultures allowed us to determine the rates of interconversion in these three subpopulations. Although both nCSCs and hCSCs converted to mCSCs, mCSCs mostly exhibited symmetric division (Fig. 7B). Collectively, these results suggest that enrichment of CD44^hi cells at the invasive front may be partly associated with transdifferentiation of nCSCs into mCSCs.

In this paper, we began by mapping biomechanical heterogeneity in breast cancer cell lines by measuring cell size and cell deformability of two cell types, which vary in the extent of cell-cell adhesion. Although both cell types exhibit overlapping size and deformability distributions, MDA-MB-231 cells are on an average larger in size and more motile compared to MCF-7 cells. Our simulations predict that cancer invasiveness is maximum when heterogeneity in both cell size and cell deformability is present. Although a directional cue may drive collective invasion of the entire cell aggregate independent of the extent of cell-cell adhesion, invasiveness is most efficient at intermediate cell-cell adhesion. Analysis of subpopulations present at the invasive front suggests that small cells of varying deformabilities represent the most invasive subpopulation and may correspond to a subpopulation of CD44^hi CSCs. In addition to highlighting the potential role of combined heterogeneity in cell size and deformability in hastening cancer invasion, our results suggest biophysical heterogeneity drives invasiveness (Fig. 7C). Although our simulations were performed keeping the extent of cell-cell adhesion fixed for a given simulation, variable extent of cell-cell adhesion, which is likely to be more relevant in vivo, is expected to generate a combination of single cells and cell clusters at the invasive front.

Our results predict that heterogeneity in cell size alone or deformability alone is not sufficient to positively influence cancer invasiveness; instead, heterogeneity in both cell size and deformability is required. Combined heterogeneity may promote invasiveness by allowing cells to continuously reposition themselves within the cell cluster. Such internal reorganization is particularly evident during collective migration with intact cell-cell adhesions (i.e. J_cc=1), wherein cells cover the largest distances to reach the boundary and also undergo maximum neighbor changes. This type of re-arrangement enabled the enrichment of cells from the tumor interior, with these cells being softer than the bulk population. In contrast, at moderate and low cell-cell adhesion, as cells are more prone to scattering, the scope for re-arrangement may be lesser. Under these conditions, initial position of cells is likely to play a more important role in determining their presence at the invasive front. Highest population invasiveness observed for J_cc=16 may be attributed to the synergistic effects of intermixing and initial positioning. Our observation of fastest chemotaxis observed at intermediate cell-cell adhesion is consistent with the notion that groups of cells exhibit better gradient-sensing via collective guidance through the regulation of contact inhibition of locomotion (CIL) (Camley et al., 2016).

In skin and mammary primary tumors, screening of a large panel of cell surface markers identified different tumour subpopulations exhibiting different extents of EMT, with complete EMT requiring cells to transit through multiple intermediate hybrid states (Pastushenko et al., 2018). The hybrid epithelial/mesenchymal state is associated with higher tumor initiating potential, increased metastasis and chemoresistance (Grosse-Wilde et al., 2015; Biddle et al., 2016; Navas et al., 2020). In breast cancer cells, although CD44^hi CD24⁻ CSCs are mesenchymal in nature, ALDH⁺ CSCs, as well as CD44^hi CD24⁺  CSCs, exhibit a hybrid epithelial/mesenchymal phenotype (Liu et al., 2014; Grosse-Wilde et al., 2015). In our simulations, J_cc=16 may correspond to this hybrid epithelial/mesenchymal state in which both single cells and cell clusters (i.e. groups of cells with intact cell-cell contacts) were observed. Cells of intermediate deformability, predicted by our simulations to be enriched at the invasive front for J_cc=16, may correspond to a subset of the CD44^hi CD24⁺  hCSCs. Although nCSCs gave rise to predominantly mCSCs in our experiments, conversion of nCSCs to a CD44^hi CD24⁺  drug tolerant state in the presence of the chemotherapeutic drug docetaxel highlights the importance of stress-induced phenotypic plasticity (Goldman et al., 2015). The extent to which phenotypic plasticity is influenced by physical stresses felt by the cell nucleus during three-dimensional migration remains to be explored (Mukherjee et al., 2020).

A correlation between smaller cell size and stemness has been speculated, both in normal and diseased contexts (Lia et al., 2015). Consistent with this notion, compared to larger PC3 pancreatic cancer cells, smaller-sized cells were found to possess more CSC-like characteristics, including clonogenic and tumorigenic properties (Li et al., 2008). Similarly, spheres generated by culturing A431 epidermoid cancer cells under non-adherent conditions were enriched with podoplanin⁺ small cells with increased tumorigenic potential (Bortolomai et al., 2010). In SUM149 inflammatory breast cancer cells, ALDH⁺ CSCs have been shown to be softer compared to the parental population (Chen et al., 2019). Here, we show that in MDA-MB-231 cells, the two CSC subpopulations differ considerably in their size and deformability. Although CD44^hi CD24⁻ mCSCs were significantly smaller and softer than parental cells, similar to the above studies, CD44^hi CD24⁺  hCSCs exhibited a wide variation in their size and stiffness. Interconversion experiments suggest that although both nCSCs (CD44^−/low cells) and hCSCs convert primarily to mCSCs, mCSCs exhibit symmetric division. This may partly explain the enrichment of CD44^hi cells at the invasive front. As mCSCs represent the majority fraction of CD44^hi cells, it is likely that the majority of the CD44^hi cells at the spheroid invasive front are mCSCs. In line with this, a computational study predicted that EMT-inducing factors, such as TGFβ, secreted by stromal cells can give rise to the spatial segregation of different CSC subsets with mCSCs localized at the migration front (Bocci et al., 2019b). However, whether or not such spatial segregation of CSCs also occurs in cell types with moderate or strong cell-cell adhesion remains to be explored.

Using computational modeling, we have previously shown that, in addition to the proportion of CSCs, the extent of their spatial scattering also influences tumor growth (Kumar et al., 2016a). In our spheroid invasion assays, CSCs and nCSCs were randomly localized at day 0. However, within 24 h, enrichment of CD44^hi cells in the outermost zone of the spheroids was clearly observed, suggesting that CSCs serve as leader cells. This is in line with previous reports in primary human breast cancers wherein CD44^hi CD24⁻ mesenchymal cells localized at the tumor invasive front, and ALDH⁺ cells were found in the tumor interior (Liu et al., 2014). Although similar localization of CD44^hi cells at the invasive front was observed during the collective invasion of luminal breast cancer cells, these cells were found to exhibit a hybrid epithelial/mesenchymal gene signature (Yang et al., 2019). Although Notch signaling is known to induce EMT via its ligands Delta and Jagged, EMT induction through Jagged gives rise to clusters of hybrid epithelial/mesenchymal cells (Boareto et al., 2016). The proportion of these hybrid epithelial/mesenchymal cells and the size of the cell clusters has been predicted to depend on the rate at which cells undergo EMT and cell-cell adhesion (Bocci et al., 2019a). This model can successfully reproduce the typical cluster size distributions observed in CTCs isolated from patient blood.

In conclusion, in this study we have established the advantage of phenotypic heterogeneity in breast cancer invasion. Future experiments taking cell cycle stage into account can help us assess the contribution of cell cycle heterogeneity to our measurements of biophysical heterogeneity. Our multiscale stochastic computational model recapitulating experimental observations predicts a potential role of combined heterogeneity in cell size and deformability in promoting cancer invasion, with gradual enrichment of small cells at the invasive front, which may experimentally correspond to CSCs. Future studies focused on probing interactions between the different subpopulations may provide us with strategies to inhibit heterogeneity-driven cancer invasion.

MCF-7 and MDA-MB-231 breast cancer cell lines were obtained from National Center for Cell Science (Pune, India) and cultured in high glucose Dulbecco's modified Eagle medium (DMEM, Invitrogen, AL007A) containing 10% fetal bovine serum (FBS, HiMedia, RM9952), maintained at 37°C at 5% CO₂ humidified atmosphere, and passaged using 0.25% trypsin-EDTA (HiMedia, TCL099). For experiments, cells were cultured for 24 h on glass coverslips coated with rat tail collagen I (Sigma-Aldrich, C3867) at a concentration of 25 μg/ml. For E-cadherin staining, cells were plated sparsely on collagen-coated glass coverslips. After 72 h, cells were fixed, permeabilized with 0.5% Triton-X 100 (Sisco Research Laboratories, 64518) and stained with anti-mouse E-cadherin antibody (Cell Signal Technology, 14472, 1:400) overnight at 4°C. Next day, after washing, cells were incubated with goat-anti-mouse IgG Alexa Fluor 555 (Invitrogen, A21422, 1:1000) for 2 h at room temperature. After staining nuclei with Hoechst 33342 (Thermo Fisher Scientific, H3570), images were captured using an Olympus IX81 inverted microscope at 40× magnification. Quantification of E-cadherin intensity was performed using Fiji Image J software by tracking the line intensities across the cell edges and quantifying the ratio of membrane to cytoplasmic E-cadherin intensity.

For cell size measurement, after 24 h of culture on collagen-coated glass coverslips, cells were imaged at 10× magnification using an inverted microscope (Olympus IX71). Cell area was quantified using ImageJ software by manually outlining cells. For tracking temporal evolution of cell size after 24 h culture, cells were imaged every 10 mins over a duration of 24 h (Zeiss Spinning Disk Microscope). The same movies were analyzed for measuring random cell motility based on the positions of cell centroids determined by manually outlining the periphery of individual cells in ImageJ. Cell cortical stiffness measurements were obtained by indenting cells with soft pyramidal probes of nominal spring constant 160 pN/nm (22 kHz, Olympus, TR800PB), and fitting the first 500 nm of indentation data with the Hertz model. For quantification of cell size, cell motility and cell cortical stiffness, experiments were repeated three times and a comparable number of measurements from independent experiments were pooled.

For flow cytometric analysis, MDA-MB-231 cells were stained with fluorescence-labeled antibodies CD44-FITC (BD biosciences, 555478) and CD24-PE (BD biosciences, 555428) as per the manufacturer's protocol. In brief, after trypsinization, MDA-MB-231cells were incubated with CD44-FITC and CD24-PE for 45 min. After labeling, cells were washed with FACS buffer (PBS with 1% FBS) and cell sorting was performed in FACS Aria Fusion (Becton Dickinson), using FACS Diva software (Becton Dickinson). Suitable negative isotype controls were used to rule out the background fluorescence. Sorted cells were collected in complete medium (DMEM high glucose plus 10% FBS plus Antibiotic-Antimycotic solution (HiMedia, A002A) to perform further experiments. The percentage of each positive population was determined using quadrant or region statistics. Data were analyzed using FACS Diva Software (Beckton Dickinson).

For spheroid invasion experiments, MDA-MB-231 cells were trypsinized and resuspended in medium. Spheroids were generated using the hanging drop method by suspending 4000 cells in 15 µl medium supplemented with 6.25 µg/ml rat tail collagen I (Corning, 354249). The drops were then incubated at 37°C in 5% CO₂ for 48 h to generate spheroids. After isolation, spheroids were seeded on a bed of 1.5 mg/ml collagen gels, and then incubated with three-dimensional collagen solution at 37°C for gelation to occur. After 30–45 mins, wells were flooded with medium. Cells were fixed and stained for nuclei, F-actin (Thermo Fisher Scientific, A22282) and CD44 (Novus, NBP1-47386) at day 0 (2 h after spheroids were embedded in collagen gels), day 1, day 3 and day 5. For quantification of CD44 intensity, images were captured using an LSM microscope (Zeiss) at 10× magnification under identical gain and exposure settings for all the conditions. Captured images were processed and analyzed using Fiji ImageJ. Raw integrated intensities of single cells were determined by manually outlining individual cells using the selection tool in Fiji ImageJ. For obtaining the normalized zone-wise distribution of CD44 intensity, the average integrated CD44 intensity of each zone, determined using the raw intensities of cells within that zone, was divided by the average integrated CD44 intensity of the innermost zone. This work was approved by the biosafety committee of the institute.

For studying the impact of phenotypic heterogeneity on cancer invasion, we adapted our recently developed CPM-based two-dimensional model of cancer invasion, taking cell size and cell deformability into account (Kumar et al., 2016b). In our model, a cell aggregate (comprised of 69 cells) is situated at the center of a two-dimensional ECM lattice (1 mm×1 mm size) comprised of randomly positioned ECM fibers, with each ECM fiber being 2 µm in thickness and 30-40 µm in length. In simulations taking phenotypic heterogeneity into account, size and area constraint (λ_a) of individual cells were assigned values from normal distributions as depicted in Fig. 2B. Although cell sizes were drawn from a normal distribution with mean size (⁠ and s.d. σ_A = 86.1 μm², cell deformabilities (i.e. λ_a’s) were drawn from a normal distribution with and s.d. ⁠. For homogeneous cell clusters, cell size and cell deformability of all the cells were set to the means of the above distributions, i.e. A=517.5 μm² and λ=1E/L⁴. Based on cell size and deformability, for initial analysis, the heterogeneous population was divided into five groups, namely small and soft cells (⁠⁠), small and stiff cells (⁠⁠), cells with intermediate size and deformability (⁠⁠), large and soft cells (⁠⁠, ⁠), and large and stiff cells (⁠⁠, ⁠). Apart from ECM fibers and cells, the remaining pixels in the lattice were considered as ‘fluid’ pixels mimicking the interstitial fluid. The different extent of adhesion between the different types of pixels was accounted for by using distinct adhesion energies (J_xy) between entities x andy. Although cell-matrix (J_cm), cell-fluid (J_cf), matrix-fluid (J_mf) and fluid-fluid (J_ff) adhesion energies were held constant (J_cm=16 k_BT/L, J_cf  =32 k_BT/L, J_mf  =  J_ff  =35 k_BT/L), cell-cell adhesion energies (J_cc) were varied to mimic different extents of cell-cell adhesion, with J_cc=1 k_BT/L mimicking collective cell migration and J_cc = 40 k_BT/L mimicking single cell migration.

Cancer cells are known to invade by MMP-mediated degradation of the surrounding ECM depending on the physicochemical properties of the ECM, including ligand density and stiffness (Das et al., 2013, 2017). Soluble MMPs activated by membrane-bound MT1-MMP then diffuse into the extracellular space and degrade the surrounding ECM (Kumar et al., 2018a). To incorporate this effect, we have integrated our CPM model with reaction-diffusion MMP dynamics. To mimic the ECM density-dependent MMP secretion profile observed in breast cancer cells, in our model, the number of MMP molecules secreted per unit time by a cell was assumed to be 0.05 sec⁻¹ at the site of cell-fiber contact. Diffusion and degradation of soluble MMPs were incorporated in our model using the reaction diffusion formalism given by the equation ⁠, where represents the concentration of MMP molecules at point and time t, D represents the diffusion coefficient, and δ the degradation rate of soluble MMPs. D and δ were chosen to be 1.0×10⁻⁹cm²sec⁻¹ and 0.002 sec⁻¹, respectively.

The complete simulation framework was implemented using the open source package CompuCell3D (CC3D) (Swat et al., 2012). For visualization, *.vtk files were generated from the CC3D simulations and visualized in CC3D itself. For quantifying migration trajectories, the cell centroid was tracked and logged to CSV files. These files were then processed in R using custom written scripts to extract different invasiveness metrics. The CC3D code used to implement the complete model is available at https://github.com/sandeep13712/senlab-phenotypicheterogeneity.

The breast cancer dataset was downloaded from the GEO repository (GSE75688) (Chung et al., 2017). Non-tumor and bulk samples were removed from the dataset. Then, control RNA readings were removed, and multiple readings for the same gene were added up to get a dataset with unique gene names. An obtained preprocessed dataset was then used for GSEA using the single sample GSEA method (Barbie et al., 2009). Stemness and EMT enrichment analyses were performed using MALTA_CURATED_STEMNESS_MARKERS and HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION gene sets downloaded from the Molecular Signatures Database (Subramanian et al., 2005; Liberzon et al., 2011; Malta et al., 2018). For calculating the enrichment of biophysical characteristics, a custom gene set BiophysicalGeneset comprising MYH9, MYH10, MYO10, MYO5B, MYO5C, RAC1, RAC2, ACTB, ACTG1, ACTN1, ACTN4, VIM, ROCK1, ROCK2 and CDC42 was manually curated and used. For tSNE visualization, data were first log transformed and then the first ten principal components were used for computing the two-dimensional-tSNE projection. To quantify biophysical heterogeneity in the MCF7 and MDA-MB-231 cell lines from the GSE115869 and GSE124989 datasets (Birts et al., 2020; Jonasson et al., 2019), downloaded datasets were log-transformed after adding 1 to all transcript counts and then the distributions of genes relevant for different aspects of cellular mechanics were quantified. Analysis was performed using Python scripts [available in github repository (https://github.com/sandeep13712/senlab-phenotypicheterogeneity)].

We thank IIT Bombay-Industrial Research and Consultancy Centre for providing Bio-AFM, FACS and confocal microscopy facilities.

The peer review history is available at https://journals.biologists.com/jcs/article-lookup/134/7/jcs250225/