Breast cancer progresses in a multistep process from primary tumor growth and stroma invasion to metastasis. Nutrient-limiting environments promote chemotaxis with aggressive morphologies characteristic of invasion. It is unknown how coexisting cells differ in their response to nutrient limitations and how this impacts invasion of the metapopulation as a whole. In this study, we integrate mathematical modeling with microenvironmental perturbation data to investigate invasion in nutrient-limiting environments inhabited by one or two cancer cell subpopulations. Subpopulations were defined by their energy efficiency and chemotactic ability. Invasion distance traveled by a homogeneous population was estimated. For heterogeneous populations, results suggest that an imbalance between nutrient efficacy and chemotactic superiority accelerates invasion. Such imbalance will spatially segregate the two populations and only one type will dominate at the invasion front. Only if these two phenotypes are balanced, the two subpopulations compete for the same space, which decelerates invasion. We investigate ploidy as a candidate biomarker of this phenotypic heterogeneity and discuss its potential to inform the dose of mTOR inhibitors (mTOR-I) that can inhibit chemotaxis just enough to facilitate such competition.

Significance:

This study identifies the double-edged sword of high ploidy as a prerequisite to personalize combination therapies with cytotoxic drugs and inhibitors of signal transduction pathways such as mTOR-Is.

Invasion and infiltration are hallmarks of advanced cancers, including breast cancer, and accumulating evidence suggests that invasive subclones arise early during tumor evolution (1).

Infiltrating and invasive phenotypes are often observed among high-ploidy cells. Converging evidence from different cancer types, including colorectal, breast, lung, and brain cancers, suggests a strong enrichment of high-ploidy cells among metastatic lesions as compared with the primary tumor (2, 3). Even in normal development, trophoblast giant cells, the first cell type to terminally differentiate during embryogenesis, are responsible for invading the placenta and these cells often have hundreds of copies of the genome (4). Coexistence of cancer cells at opposite extremes of the ploidy spectrum occurs frequently in cancer and is often caused by whole-genome doubling (WGD). Similar to infiltration, the timing of WGD is early in tumor progression across several cancer types (5, 6), including breast cancer. Tetraploid cells resulting from WGD often lose regions of the genome, giving rise to poly-aneuploid cancer cells (PACC). Multiple studies have described a minority population of PACCs with an unusual resilience to stress (7–9). A very recent investigation of evolutionary selection pressures for WGD suggests that it mitigates the accumulation of deleterious somatic alterations (10). However, it is not clear what cost cells with a duplicated genome pay for this robustness.

To address this question, we developed a mathematical model of high- and low-ploidy clones under various energetic contingencies. We calibrated the model to recapitulate doubling times and spatial growth patterns measured for the HCC1954 ductal breast carcinoma cell line via MEMA profiling (11). This includes exposure of HCC1954 cells to hepatocyte growth factor (HGF) in combination with 48 extracellular matrices (ECM), followed by multi-color imaging (12). Sequencing (13) and karyotyping studies (14, 15) of cancer cell lines have shown that technically, sequences obtained from a cancer cell line encode a metagenome (16), because they represent the aggregate genomes of all clones that coexist within the cell line. To capture this heterogeneity, we modeled how high- and low-ploidy clones coevolve and how that affects the invasiveness of the metapopulation. Our results show that long-term coexistence of low- and high-ploidy clones occurs when sensitivity of the latter to energy scarcity is well-correlated to their chemotactic ability to populate new terrain. Higher energy uniformity throughout population expansion steers selection in favor of the low-ploidy clone, by minimizing the fitness gain the high-ploidy clone gets from its chemotactic superiority. Better understanding of how these two phenotypes coevolve is necessary to develop therapeutic strategies that suppress slowly proliferating, invasive cells before cytotoxic therapy favors them.

We first introduce the conceptual framework that led to the model, formulate the model equations, and then we derive analytic and numerical solutions. Finally, we describe drug sensitivity and RNA sequencing (RNA-seq) data analysis of cell lines with different ploidies.

Overall model design

The use of partial differential equations (PDE) over a stochastic-based approach, such as agent-based modeling, permits us to make predictions based on analytic results derived from the subsequent PDEs and an increase in computational efficiency.

We modeled growth dynamics in polyploid populations of various subpopulation compositions. Appealing to a continuity description and assuming a continuum approximation of the cellular and energy concentration is valid, we derived a system of coupled PDEs. Each compartment in the PDE describes the spatio-temporal dynamics of the quantity of interest (e.g., energy or cellular dynamics). At the core of our model lies the assumption that chemotactic response to an energy gradient is a function of the cell's energetic needs. This trade-off implies that heterogeneous populations will segregate spatially, with higher energy-demanding cells leading the front of tumor growth and invasion. In contrast, for an energy-rich environment, we expect the cells to grow in a similar way as they will have no need to search for places of higher energy density.

We modeled competition for energy in a heterogeneous population, consisting of goer and grower subpopulations, to predict their behavior during plentiful and energy sparse conditions. We assumed both goer and grower have the same random cell motility coefficient Γ, the same chemotactic coefficient χ, and maximal growth rate λ. Sensitivity to the available energy was modeled via a Michaelis–Menten type equation with coefficients |$\Phi _i}$| that determine the amount of energy that population, i, needs for a half maximal growth rate. Chemotactic motion is asymmetrically sensitive to the amount of energy available. This is accounted for by |$Ξ _i}$|⁠. The goers (U) are more motile and require more energy compared with the growers (V). This manifests itself mathematically via the parameter relations |$\Phi _U} \gt{\Phi _V}$| and |$Ξ _U} \lt {Ξ _V}$|⁠.

Quantitative estimates of how a cell's growth rate and motility depend on energy availability have been described previously (17–19). Energetic resources come in various forms, and the identities of the limiting resources that ultimately drive a cell's fate decision vary in space and time. We used MEMA profiling to investigate what likely is only a narrow range of that variability, 48 HGF-exposed ECMs (12). HGF stimulates both growth and migration of epithelial and endothelial cells in a dose-dependent manner, whereby maximal growth-stimulating effects have been reported at concentrations twice as high as concentrations that maximize migration (20). In-line with these reports, our model demonstrates a shift from proliferation to migration as resources get depleted.

Mathematical models of a dichotomy between proliferation and migration are numerous (21–23), but whether the two phenotypes are indeed mutually exclusive remains controversial (24). Our efforts to use mathematical modeling to inform what cost high-ploidy cells (goers) pay for their robustness build upon these prior works. We extend it by accounting for differences in the rate at which cells consume energy (δ) and differences between the rate at which energy diffuses in media (coefficient ΓE). For mid-range energy diffusion coefficients, our model describes directed cell motility in response to a gradient of a soluble attractant, that is, chemotaxis. In contrast, small values of ΓE approximate cell motility toward insoluble attractants, that is, haptotaxis. As such, the chosen value for ΓE determines where along the continuum between haptotaxis- and chemotaxis-directed cell movement resides. A special case applies when ΓE is very large relative to cell movement, in which case neither chemotaxis nor haptotaxis occurs. All these energetic contingencies determine whether phenotypic differences between goers and growers manifest as such, and explain why nonproliferative arrested cells can have the same motility as cycling cells (24).

Quick guide to equations

Our model assumes that the energy diffusion coefficient, ΓE, depends on the type of media or surface upon which the cells grow. We also suppose that energy is consumed in proportion to the number of cells present. For cell motility, we assume it is driven both by random cell motion and chemotaxis. This leads to our general coupled system:

formula
formula
formula

The system (Eqs. A–C) is defined on a dish of radius R, subject to no-flux boundary conditions. We assume the goers (U) and growers (V) are initially concentrated at the center with radius |$\rho _0} \le R$| and initial concentrations |$U_0}\ {\rm{and}}\ {V_0}$|⁠. Furthermore, we assume that the energy density is uniformly distributed on the plate with initial value E0. All parameters, except otherwise stated, are independent of any of the state variables (⁠|$\tilde{E},U,{\rm{and}}\ V$|⁠). The energy is consumed at rate δ. Both cells can divide at the maximal rate λ, but are restricted by the energy density |$\tilde{E}$|⁠. The cells can locally grow to a local maximal density given by K. This parameter is often cell line dependent and is related to contact inhibition and a cell's ability to grow on top of other cells.

We now convert our system to dimensionless form that was used for all subsequent simulations and analysis via appropriate rescaling (Supplementary Materials and Methods),

formula
formula
formula

Hereby, rescaling simplified the system from 13 (10 parameters and three initial conditions) to nine (seven parameters and two initial conditions), with dimensionless variables: |$\gamma _E} = \cr {\Gamma _E}/\chi ,\;a = \delta K/\lambda ,\;b = \Gamma /\chi$|⁠, |$\phi = \Phi /{E_0}$|⁠, and |$\xi = Ξ /{E_0}$|⁠.

Analytic estimates of infiltration

The desire for cells to move is inherently tied to the availability of nutrients and space. To this end, we define |$Ψ ( \tau ): = [ {\rho ( \tau ) - {\rho _0}} ]/{\rho _0}$|⁠, where |$\rho ( 0 ) = {\rho _0}$|⁠, the initial radius of cell seeding density. Ψ can be thought of as a nondimensional measure of infiltration attained after time τ. This dimensionless measure has the added benefit of being scale independent. An inherent difficulty with random cell motility and calculating infiltration is that the system always reaches the boundary of the dish in finite time. Instead, we defined the maximum degree of infiltration to be given by the time needed for the total energy to be below a threshold, |$\varepsilon \ll 1$|⁠.

In general, the maximum degree of infiltration is difficult to predict analytically, so we only considered the single subpopulation case when obtaining our analytic estimates. We also made use of the simplification that most energy type molecules (e.g., glucose) have a diffusion coefficient that is very large, relative to cell movement. This allowed us to write a reduced model which has energy homogeneous in space (see Supplementary Eqs. S3a–S3b in Supplementary Materials and Methods). We estimated infiltration analytically through two different approaches. First, we extracted an ordinary differential equation (ODE) system that couples nutrient consumption to the wave front location. Second, we derived estimates of the infiltration achieved by using the cell concentration after it has become uniform.

No chemotaxis infiltration estimate

We can derive an estimate for infiltration in the absence of chemotaxis by appealing to Supplementary Eqs. S3a–S3b. These logistic growth reaction diffusion models often exhibit complex dynamics. One such example is that of a traveling wave solution, where one state, typically the stable state, travels (infiltrates) through the domain. The canonical example of this phenomenon is the Fisher–KPP. In contrast to the Fisher–KPP and other traveling wave problems typically studied, our model has a decaying growth rate and so the magnitude of the nonlinearity that caused the traveling wave tends to zero. Therefore, in the classic sense, our system does not admit a traveling wave. We here extend the theory by assuming a separation of time scales between consumption of energy (e.g., decay of energy-dependent growth rate) and the speed of the traveling wave.

To begin, we made the assumptions that the wave speed is a slow function of r and τ. The solution obtained will verify that these assumptions are valid for our system. Our ansatz takes the form |$u( {r,\tau } ) = U( {r - \eta \tau } ) = U( z )$|⁠. Note that in spatial equilibrium, u = 1 is stable and u = 0 is an unstable steady state. If the unstable state is what governs the wave speed, then the wave is said to be “pulled,” otherwise it is “pushed” (25, 26). Following Witelski and colleagues (27), the resulting analysis yielded a coupled system of ODEs that govern the speed of the front (Supplementary Materials and Methods):

formula
formula

where we have the dimension n of the system (e.g., |$n = 2,\ 3$|⁠) entering in both the consumption of energy and influencing the second term of the wave front location ρ. We saw that our assumptions on the behavior of the traveling wave were verified, because ρ is assumed much larger than 0 and |$0 \ll \rho \ll R$|⁠. This shows that, E is a slowly varying function of time and |$\eta = d\rho /d\tau$| is a slowly varying function of time and its current distance from the center. In other words, a traveling wave will only form when the initial seeding radius is relatively large (⁠|$\rho _0} \sim 1$| was sufficient in most simulations).

Estimating the degree of infiltration from equilibration

The previous section yields a system that can be integrated to track the evolution of the cell front over time. However, we may be more interested in how far it will ultimately travel before energy exhaustion and not the speed at which it gets there. An interesting alternative to tracking the wave over time is to only assume it travels as a wave, but only record the density after the system has reached uniformity. This is possible if the death rate (which has been neglected) is much smaller than the time it would take the cells to spread uniformly. If this is the case, we can bound the degree of infiltration from only knowing the uniform value |$\bar{u}$| at the end of the experiment (see Supplementary Materials and Methods for details):

formula

where Λ is the transition width, that is, the length scale on which cell concentration goes from u = 1 to u = 0, and |$\bar{u}$| is the concentration at the end of the experiment (at equilibration). Solving Eq. I gives us the estimated infiltration as:

formula

Because we assumed that the transition width is unknown, we can bound ψ [or ρ(T)] by considering the lower and upper bounds |$\Lambda = 0,\ 1$|⁠, respectively.

Numerical estimates of infiltration

When directed cell movement is not negligible, analytic approximations are more difficult to obtain, and numerical simulations are preferred to estimate the degree of infiltration. For simulations where we measured infiltration, we took the 1-norm |$\parallel\!\! E{\parallel _1} = \mathop \sum_i | {{E_i}} | \,\lt\,\varepsilon = 10^{- 4}$| as the threshold defining minimum energy requirements. We calibrated the model to recapitulate doubling times and spatial growth patterns measured for the HCC1954 ductal breast carcinoma cell line via MEMA profiling. The dataset included exposure of HCC1954 cells to HGF in combination with 48 ECMs in a DMSO environment (i.e., no drug was added to the media). Between 13 and 30 replicates of each ECM were printed on a rectangular MEMA array as circular spots (Supplementary Fig. S1A), adding up to a total of 692 spots (12).

MEMA data analysis

An average of 31 cells (90% confidence interval, 23–41) were seeded on each 350-μm spot and grown for 3 days. The median number of cells imaged at day 3 was 121, which falls within the range expected from the estimated seeding density and a doubling time reported for this cell line (72–128 cells). Confluence at seeding was calculated from the ratio between the cumulative area of cells and the area of the spot (see Supplementary Fig. S1B and S1C).

Quantification of segmented, multi-color imaging data obtained for each spot 3 days after seeding was downloaded from Synapse (synID: syn9612057; plateID: LI8C00243; and well: B03). Cells tended to grow in a symmetric, toroidal shape (Supplementary Fig. S1D), albeit considerable variability was observed across the ECMs. We binned cells detected within a given spot according to their distance to the center of the spot and calculated the confluence of each bin (Supplementary Materials and Methods). This was then compared with the confluence obtained from the simulations as described below (Supplementary Data S1–S3).

Simulation environment

Simulations were run to recapitulate logistics and initial growth conditions on the MEMA array (Supplementary Fig. S1A), the spatial domain was circular with radius |$R = 1\!,\!750$| μm and the temporal domain was 3 days. Cells were seeded uniformly at the center of this domain along a radius |$\rho _0} = 175$| μm at 36% confluence. To recapitulate the configuration of the MEMA profiling experiment, cells leaving the ρ0 domain can no longer adhere to the ECM and die. This was implemented by introducing the carrying capacity, K(x), rapidly approaching zero when x > ρ0. This setup can lead to energy attracting cells to the periphery of a MEMA spot and beyond. We ran 161,000 simulations (Supplementary Data S1) at variable energy consumption rates, chemotactic/haptotactic coefficients, energetic sensitivities, and diffusion rates of the growth-limiting resource (i.e., ECM bound HGF; Table 1; Supplementary Table S1; Supplementary Data S2). Ninety-two percent of these simulations were completed successfully. For each simulation/ECM pair, we compared spatial distributions of in silico and in vitro confluence (Supplementary Data S3) using the Wasserstein metric (28).

Table 1.

Measured and inferred simulation parameters for one-population model (second column); and two-population model (third column).

Experimentally measured parametersSymbolHomogeneous cells on any ECMHeterogeneous cellsUnit
Maximal growth rate λ 0.56 0.56 Per day 
Goer's seeding density u0 36.0 3.8 % confluence 
Grower's seeding density v0  1.3 % confluence 
Simulation time t 127 Days 
Radius of domain |$\hat{\bf R}$| 1,750 1,750 μm 
Seeding radius ρ0 175 175 μm 
Growth beyond ρ0  No Yes  
Experimentally measured parametersSymbolHomogeneous cells on any ECMHeterogeneous cellsUnit
Maximal growth rate λ 0.56 0.56 Per day 
Goer's seeding density u0 36.0 3.8 % confluence 
Grower's seeding density v0  1.3 % confluence 
Simulation time t 127 Days 
Radius of domain |$\hat{\bf R}$| 1,750 1,750 μm 
Seeding radius ρ0 175 175 μm 
Growth beyond ρ0  No Yes  
Inferred parametersSymbolHomogeneous cells on GAP43Heterogeneous cellsUnit
Energy diffusion coefficient on ECM ΓE 9,500 8,982 μm2/day 
Energy consumption rate δ 57 23 Per cells/day 
Energy needed for half-maximal growth Φu (0.01, 3) 0.310  
 Φv  (0.05, 0.16)  
Chemotactic coefficient χ 5,500 1,919 μm2/day 
Energy deficit inducing chemotaxis ξu (2, 300) (1, 8)  
 ξv  140  
Inferred parametersSymbolHomogeneous cells on GAP43Heterogeneous cellsUnit
Energy diffusion coefficient on ECM ΓE 9,500 8,982 μm2/day 
Energy consumption rate δ 57 23 Per cells/day 
Energy needed for half-maximal growth Φu (0.01, 3) 0.310  
 Φv  (0.05, 0.16)  
Chemotactic coefficient χ 5,500 1,919 μm2/day 
Energy deficit inducing chemotaxis ξu (2, 300) (1, 8)  
 ξv  140  

Note: Parentheses indicate ranges with equally good fits. Corresponding nondimensional parameters are shown in Supplementary Table S1. Net chemotaxis of cells is given by the interaction of the chemotactic coefficient (χ) and energy deficit inducing chemotaxis (ξu), that is, effective cell chemotaxis is approximately χ/ξu. The larger ξu, the later the cells sense an energy deficit, that is, the larger the energy gradient must be in order for them to accelerate their movement in response to it.

The model was implemented in C++ (standard C++11). The armadillo package (ARMA version: 9.860.1; ref. 29) was used for simulation of the PDEs. Simulations were run on a Intel Core i7 MacBook Pro, 2.6 GHz, 32 GB RAM. The source code is available at the GitHub repository for the Integrated Mathematical Oncology department: GoOrGrow (https://github.com/MathOnco/GoOrGrow_PloidyEnergy). An R Markdown script was included, enabling replication of comparisons between experimental data and model simulations (Supplementary Data S1–S3).

Ploidy as biomarker of phenotypic divergence

We identified 44 breast cancer cell lines of known ploidy (30) and with available RNA-seq data in Cancer Cell Line Encyclopedia (CCLE; ref. 31) and analyzed their drug sensitivity and expression profiles as follows.

Drug sensitivity analysis

We used growth rate inhibition (GR) metrics as proxies of differences in drug sensitivities between cell lines. Unlike traditional drug sensitivity metrics, like the IC50, GR curves account for unequal division rates, arising from biological variation or variable culture conditions, a major confounding factor of drug response (32). Previously calculated GR curves and metrics were available for 41 of 44 breast cancer cell lines. A total of 46 drugs had been profiled on at least 80% of these cell lines and their area over the curve (GRAOC) drug–response metric (33) was downloaded from GRbrowser (http://www.grcalculator.org/grtutorial/Home.html). For each drug, we calculated the z-score of GRAOC across cell lines to compare drugs administered at different dose ranges. Of these 46 drugs, 39 could be broadly classified into two categories as either cytotoxic (25 drugs) or inhibitors of signaling pathways (14 drugs; Supplementary Table S2). We then evaluated a cell line's ploidy as a predictor of its GRAOC value using a linear regression model. Because molecular subtype of breast cancer cell lines is known to influence drug sensitivity, we performed a multivariate analysis, including the molecular subtype as well as an interaction term between ploidy and drug category into the model.

RNA-seq analysis

The molecular subtype classification of all cell lines was available from prior studies (34–37). Of these 44 cases, four were suspension cell lines and were excluded from further analysis. Of the remaining 40 cell lines, 20 originated from primary breast cancer tumors and were the focus of our analysis. Gene expression data were downloaded from CCLE. We used gene set variation analysis to model variation in pathway activity across cell lines (38). Pathways for which less than 10 gene members were expressed in a given cell lines were not quantified. The gene membership of 1,417 pathways was downloaded from the REACTOME database (ref. 39; v63) for this purpose.

High-ploidy breast cancer cell lines have increased metabolic activity and cell motility

To better understand the phenotypic profile of high-ploidy cells, we compared the ploidy of 41 breast cancer cell lines with their response to 46 drugs. For a drug–response metric, we used the integrated effect of the drug across a range of concentrations estimated from the GRAOC (32, 33). We observed that cytotoxic drugs and drugs inhibiting signal transduction pathways were at opposite ends of the spectrum (Fig. 1A). Namely, ploidy was negatively correlated with the GRAOC for several cytotoxic drugs and positively correlated with the GRAOC of various mTOR inhibitors, suggesting high-ploidy breast cancer cell lines tend to be resistant to DNA-damaging agents, while sensitive to drugs targeting nutrient sensing and motility.

Figure 1.

Ploidy, pathway activity, and drug sensitivity across breast cancer cell lines from CCLE. A, High-ploidy breast cancer cell lines are resistant to cytotoxic drugs, but tend to be more sensitive to inhibitors of mTOR, EGFR, and MAPK signaling pathways. Hereby ploidy is defined as the number of chromosomes in the cell line's consensus karyotype, weighted by chromosome size. Only drugs with a Pearson correlation coefficient at or above 0.2 are shown here. B, Distribution of ploidy within and across three molecular breast cancer subtypes. C, Regression coefficient of ploidy as predictor of GRAOC has opposite signs depending on drug category across all subtypes. D, Distribution of ploidy across 20 primary, adherent breast cancer cell lines from CCLE. Ploidy is correlated with the activity of pathways involved in metabolism of vitamins and cofactors (E) and hyaluronan metabolism (F). One cell line with available MEMA profiling data, HCC1954, is highlighted (red arrow).

Figure 1.

Ploidy, pathway activity, and drug sensitivity across breast cancer cell lines from CCLE. A, High-ploidy breast cancer cell lines are resistant to cytotoxic drugs, but tend to be more sensitive to inhibitors of mTOR, EGFR, and MAPK signaling pathways. Hereby ploidy is defined as the number of chromosomes in the cell line's consensus karyotype, weighted by chromosome size. Only drugs with a Pearson correlation coefficient at or above 0.2 are shown here. B, Distribution of ploidy within and across three molecular breast cancer subtypes. C, Regression coefficient of ploidy as predictor of GRAOC has opposite signs depending on drug category across all subtypes. D, Distribution of ploidy across 20 primary, adherent breast cancer cell lines from CCLE. Ploidy is correlated with the activity of pathways involved in metabolism of vitamins and cofactors (E) and hyaluronan metabolism (F). One cell line with available MEMA profiling data, HCC1954, is highlighted (red arrow).

Close modal

We built a multivariate regression model of drug sensitivities to test the hypothesis that the relationship between ploidy and GRAOC was different for cytotoxic drugs than for inhibitors of cell signaling pathways. Molecular subtype alone (Fig. 1B), could explain 0.4% of the variability in GRAOC z-scores across cell lines (adjusted R2 = 0.0044; P = 0.026). Including ploidy into the model did not improve its predictive accuracy (adjusted R2 = 0.0037; P = 0.058). However, an interaction term between ploidy and drug category (cytotoxic: 27 drugs vs. signaling: 16 drugs) increased accuracy to explain 2.6% of variability in drug sensitivity across cell lines (adjusted R2 = 0.026; P < 1e-5; Fig. 1C). The same improvement from an interaction term between ploidy and drug category was observed in an independent dataset of half-maximal inhibitory concentration (IC50) values of 34 cytotoxic drugs and 51 signaling inhibitors obtained from the Genomics of Drug Sensitivity in Cancer database (Supplementary Fig. S2; ref. 40).

We then focused on a subset of the aforementioned 41 cell lines, namely those that had been established from primary breast cancer tumors as adherent cells (20 cell lines; Fig. 1D), and we quantified their pathway activity (see Materials and Methods). A total of 27 pathways were correlated to ploidy at a significant P value (⁠|$| {{\rm{Pearson}}\,r} | \ge 0.44$|⁠; P ≤ 0.05; Supplementary Table S3). The strongest correlations were observed for metabolic pathways such as hyaluronan metabolism and metabolism of vitamins (Fig. 1E and F). Hyaluronic acid is a main component of the ECM and its synthesis has been shown to associate with cell migration (41, 42).

These results support a model that connects high ploidy with both the chemotactic ability and metabolic energy deficit of a cell.

Infiltration of homogeneous populations

Model design was guided by the goal to describe growth dynamics along two axes: from random to directed cell motility and from homogeneous to heterogeneous cell compositions (Materials and Methods, Eqs. D–F). While the last section will step into the second axis, the following two subsections distinguish scenarios along the first axis: (i) “homogeneous nutrient environments” are environments in which random cell motility dominates throughout population growth and (ii) “heterogeneous nutrient environments” imply the formation of gradients, which cause cells to move in a directed fashion, as is the case during cellular growth on an ECM.

Homogeneous nutrient environments

When the diffusion of the nutrient occurs at a much faster time scale than the actions taken by cells, we can assume that at the time scale of cells, the nutrient is essentially uniform in space. This simplification allows us to neglect chemotactic/haptotactic motion and consider only random cell motility as the driving force that spreads the cell density throughout the dish.

We obtained analytic estimates for the degree of infiltration in a homogeneous environment that lay the groundwork for new predictions. To arrive at Eq. G–Eq. H, we employed the approximations that diffusion of energy molecules (e.g., glucose) is fast relative to cell movement and that the cells' movement through the dish can be approximated by a traveling wave (43). These assumptions were verified by comparing the front estimates with results from the full numerical model (Eqs. D–F; Fig. 2A and B).

Figure 2.

Comparing analytic approximations of the degree of infiltration with those obtained from simulations. A and B, Traveling wave solutions at energy consumption rates a = 1.5 (A) and a = 3.5 (B). C, Top and bottom boundaries of traveling wave solutions estimated from equilibration are shown as a function of consumption rate. Approximation is found by bottom and top bound Λ = 0,1 from equation D. D, Phase diagram of energy consumption and front location using the derived coupled system (Eqs. G and H). A–D, All approximations and simulations assume energy is uniformly distributed at all times, that is, chemotaxis does not take place. Parameter values for initial seeding radius (ρ0), dish radius (R), and sensitivity to low energy (ϕ) were set to 3, 10, and 0.05 respectively. Red, leading edge of wave (estimated by finding value of cell concentration closest to 0.01); blue, midpoint of wave (estimated by finding value of cell concentration closest to 0.5); purple, average of red and blue; and black lines are approximations based on analytic solutions. Time is in units of the maximal growth rate of the given cell line. Front location is given in units of the characteristic length, |$\sqrt {χ / λ}$|⁠.

Figure 2.

Comparing analytic approximations of the degree of infiltration with those obtained from simulations. A and B, Traveling wave solutions at energy consumption rates a = 1.5 (A) and a = 3.5 (B). C, Top and bottom boundaries of traveling wave solutions estimated from equilibration are shown as a function of consumption rate. Approximation is found by bottom and top bound Λ = 0,1 from equation D. D, Phase diagram of energy consumption and front location using the derived coupled system (Eqs. G and H). A–D, All approximations and simulations assume energy is uniformly distributed at all times, that is, chemotaxis does not take place. Parameter values for initial seeding radius (ρ0), dish radius (R), and sensitivity to low energy (ϕ) were set to 3, 10, and 0.05 respectively. Red, leading edge of wave (estimated by finding value of cell concentration closest to 0.01); blue, midpoint of wave (estimated by finding value of cell concentration closest to 0.5); purple, average of red and blue; and black lines are approximations based on analytic solutions. Time is in units of the maximal growth rate of the given cell line. Front location is given in units of the characteristic length, |$\sqrt {χ / λ}$|⁠.

Close modal

An interesting alternative to tracking the wave over time is to assume it travels as a wave, but only record the density after the system has reached uniformity. If the death rate is much smaller than the time needed for the cells to spread uniformly, we can bound the degree of infiltration that occurred by only knowing the uniform density of cells (equation; ref. 4; Fig. 2C).

These analytic solutions point to scaling relationships for the speed of the moving front. For highly efficient energy-using cell lines (⁠|$\phi \ll 1$|⁠; Fig. 2D), the front will evolve at a speed nearly independent of energy. In contrast, for large |$\phi \gg 1$|⁠, the speed of the front falls off as |$1/\sqrt \phi$|⁠. These predictions of the behavior of infiltration on parameters can be investigated experimentally.

Heterogeneous nutrient environments

Assumptions made in the prior section apply to standard cell cultures of adhesive cells in a typical cell culture dish, where energetic resources diffuse so fast that gradients do not form. These assumptions break down during cellular growth on an ECM. Binding to the ECM can cause soluble factors (like HGF) to act and signal as solid phase ligands (44, 45). Proteolytic degradation of these ECMs then creates haptotactic gradients. Figure 1 includes HCC1954, a near-tetraploid breast cancer cell line, whose growth on various ECMs has been measured via MEMA profiling. We analyzed HCC1954 and looked to determine whether our mathematical model can explain its spatial growth patterns. MEMA profiling resulted in considerable variability of growth patterns across different ECMs (Fig. 3A and B).

Figure 3.

Model calibration using MEMA profiling of HCC1954 cells. A–C, Experimentally measured data. Variability in cell growth patterns across ECMs is demonstrated via two example ECMs, CDH1 (A) and GAP43 (B). The average local cell densities are displayed for both (color legend). C, Projecting their 2D spatial distributions onto 1D reveals enrichment of cells at the edge of the ECM spot for GAP43, but not for CDH1. D–F, Simulated data. Comparing prior distribution of chemotactic coefficients (D) and sensitivity to low energy (E) to ECM-specific posterior distributions reveals clear differences between CDH1 and GAP43 for both parameters. F, Maximum-likelihood parameter choices for CDH1 and GAP43 result in distinct spatial distributions between the two ECMs, each of which resemble the measured distributions (C). G–I, ECM-specific model parameters. Simulations were compared with each measured ECM-specific growth pattern and ranked by their maximum similarity. G and H, The five model parameters from the top 2.3% simulations were projected onto UMAP space, revealing three clusters. Color-coding simulations by the ECM responsible for their presence in the top simulations suggest enrichment of most ECMs to only a single cluster (G). H, This was confirmed when comparing cluster membership across the 12 represented ECMs. I, All five model parameters (x-axis) show significant differences between the three clusters. ****, P ≤ 0.0001. J and K, Parameter sensitivity analysis. J, The % variance explained per parameter shows significant contribution of all five model parameters to at least one of the first three principal components (PC). K, The Sobol index (x-axis) tells us which parameter best explains which aspect of the cells' spatial distribution (color code), its skewness, confluence, or gradient near the edge.

Figure 3.

Model calibration using MEMA profiling of HCC1954 cells. A–C, Experimentally measured data. Variability in cell growth patterns across ECMs is demonstrated via two example ECMs, CDH1 (A) and GAP43 (B). The average local cell densities are displayed for both (color legend). C, Projecting their 2D spatial distributions onto 1D reveals enrichment of cells at the edge of the ECM spot for GAP43, but not for CDH1. D–F, Simulated data. Comparing prior distribution of chemotactic coefficients (D) and sensitivity to low energy (E) to ECM-specific posterior distributions reveals clear differences between CDH1 and GAP43 for both parameters. F, Maximum-likelihood parameter choices for CDH1 and GAP43 result in distinct spatial distributions between the two ECMs, each of which resemble the measured distributions (C). G–I, ECM-specific model parameters. Simulations were compared with each measured ECM-specific growth pattern and ranked by their maximum similarity. G and H, The five model parameters from the top 2.3% simulations were projected onto UMAP space, revealing three clusters. Color-coding simulations by the ECM responsible for their presence in the top simulations suggest enrichment of most ECMs to only a single cluster (G). H, This was confirmed when comparing cluster membership across the 12 represented ECMs. I, All five model parameters (x-axis) show significant differences between the three clusters. ****, P ≤ 0.0001. J and K, Parameter sensitivity analysis. J, The % variance explained per parameter shows significant contribution of all five model parameters to at least one of the first three principal components (PC). K, The Sobol index (x-axis) tells us which parameter best explains which aspect of the cells' spatial distribution (color code), its skewness, confluence, or gradient near the edge.

Close modal

We projected two-dimensional (2D) spatial distributions measured on the MEMA array onto one dimension (1D; Fig. 3C), rendering them comparable with those obtained from our simulations (Fig. 3D and E). For each simulation/ECM pair, we calculated the distance between in silico and in vitro spatial cell distributions using the Wasserstein metric (Fig. 3C and F) and ranked simulations by their minimum distance across ECMs. The top 2.3% simulation parameters were then stratified by the ECM whose spatial pattern they best resembled and compared with uniform prior parameter distributions (Fig. 3D and E). Seventy-five percent of the ECMs accounted for only 1.55% of the top simulations. The tendency of cells on these ECMs to grow in a toroidal shape was strong, suggesting it may be the consequence of nonuniform printing of ECMs onto the array. We concluded that, our model cannot explain the growth patterns on those ECMs well and focused our attention on the remaining 12 ECMs, represented by 3,429 simulations. We refer to the parameters of these simulations as inferred parameter space. Principal component analysis, uniform manifold approximation and projection (UMAP; ref. 46), and density clustering (47) of the inferred parameter space revealed three clusters (Fig. 3G), with different ECMs segregating mainly into different clusters (Fig. 3H).

The two largest clusters differed mostly in their chemotactic/haptotactic and energy diffusion coefficients; while the small cluster stood out by a high sensitivity to low energy and fast chemotactic/haptotactic response (Fig. 3I). Overall, all five model parameters showed significant differences between the three clusters, suggesting they all contribute to distinction between ECM growth patterns (Fig. 3I). This was further affirmed when looking at the % variance explained per principal component per parameter (Fig. 3J). To formalize parameter sensitivity analysis independent of ECMs, we also calculated the Sobol index (48) of each parameter. The Sobol index quantifies how much of the variability in spatial cell concentration is explained by each parameter, while accounting for all its interaction effects. Each parameter contributed to significant variance (Sobol index > 0.02; ref. 48) in at least one of three spatial statistics (Fig. 3K): skewness, confluence, and gradient near the edge of the ECM.

We observed substantial differences in chemotactic/haptotactic coefficients and energy consumption rates between ECMs (Supplementary Fig. S3). To query the biological significance of this variability, we quantified the expression of the 12 ECMs in the HCC1954 cell line (Materials and Methods). Two of the five inferred ECM-specific model parameters were correlated with RNA-seq–derived expression of the corresponding ECM: energy consumption rate (Pearson r = −0.657; P = 0.028) and sensitivity to low energy (Pearson r = 0.562; P = 0.071), although latter fell short above significance (Supplementary Fig. S4).

In summary, the posterior distributions of model parameters represent a substantial departure from the uniform priors and could explain a significant proportion of growth conditions on the HGF-exposed MEMA array. This approach identified regions of interest in the parameter search space, allowing us to focus further simulations on biologically relevant chemotactic/haptotactic coefficients and energy diffusion rates.

Infiltration of heterogeneous, chemotactic populations

Growth of cells in a given ECM environment was measured across 13–30 replicates on the MEMA platform. While our model, when calibrated to the corresponding ECM environment, could explain the observed growth pattern in the majority of these replicates, a substantial fraction could not be explained by fixed choices of sensitivity to low energy and directed cell motility (Supplementary Fig. S3). One possibility that may explain this is HCC1954 is a heterogeneous cell line, with clones of variable phenotypes coevolving. Representation of these clones among the 31 cells that were on an average sampled for each replicate may vary (Supplementary Fig. S5). This hypothesis is supported by a bimodal distribution of DNA content observed among replicating HCC1954 cells on individual ECM spots (Fig. 4A and B; Supplementary Fig. S6A and S6B). If the HCC1954 population was homogeneous, we would expect a unimodal distribution of DAPI intensity among S-phase cells of this cell line. The observation of a bimodal distribution among S-phase cells suggests that HCC1954 is likely a polyploid cell line, that is, clones of variable ploidies coexist in this cell line.

Figure 4.

Internal competition of coexisting subpopulations for same space slows down invasion of the metapopulation. A, DNA content and cell-cycle state of 162 cells growing on HGF-exposed ICAM1. B, DAPI intensity of 58 replicating (EdU+) cells shows a bimodal distribution, indicating the presence of two subpopulations, a low-ploidy population (grower) comprising approximately 55% cells and a high-ploidy population (goer) comprising approximately 45% cells. C, Arms race between the grower's energetic sensitivity (y-axis) and the goer's chemotactic ability (x-axis) reduces infiltration distance (color bar). Red circles outline parameter combinations of interest explored in D–F. D–F, Spatial distribution of goer and grower for parameter values outlined in C. Dotted lines outline extreme trajectories of expected cell concentrations due to incertitude in initial goer/grower proportions, as estimated from the silhouette coefficient of cells in B (see also Supplementary Fig. S6B). D, High chemotactic motility will cause the goer to leave the center of the dish too soon, leaving room for the grower to expand there. E, With an intermediate motility, the goer succeeds in maintaining high representation both at the center and edge of the dish. F, Low motility will prevent the goer from gaining a sufficient spatial lead from the grower while energy is still abundant, and it will lose dominance at the edge of the dish once energy becomes sparse. Red arrows indicate maximum infiltration distance achieved by either of the two populations.

Figure 4.

Internal competition of coexisting subpopulations for same space slows down invasion of the metapopulation. A, DNA content and cell-cycle state of 162 cells growing on HGF-exposed ICAM1. B, DAPI intensity of 58 replicating (EdU+) cells shows a bimodal distribution, indicating the presence of two subpopulations, a low-ploidy population (grower) comprising approximately 55% cells and a high-ploidy population (goer) comprising approximately 45% cells. C, Arms race between the grower's energetic sensitivity (y-axis) and the goer's chemotactic ability (x-axis) reduces infiltration distance (color bar). Red circles outline parameter combinations of interest explored in D–F. D–F, Spatial distribution of goer and grower for parameter values outlined in C. Dotted lines outline extreme trajectories of expected cell concentrations due to incertitude in initial goer/grower proportions, as estimated from the silhouette coefficient of cells in B (see also Supplementary Fig. S6B). D, High chemotactic motility will cause the goer to leave the center of the dish too soon, leaving room for the grower to expand there. E, With an intermediate motility, the goer succeeds in maintaining high representation both at the center and edge of the dish. F, Low motility will prevent the goer from gaining a sufficient spatial lead from the grower while energy is still abundant, and it will lose dominance at the edge of the dish once energy becomes sparse. Red arrows indicate maximum infiltration distance achieved by either of the two populations.

Close modal

To better understand the growth dynamics in a polyploid population, we used the two subpopulation version of our model, whereby variable chemotactic abilities and energetic sensitivities of goer and grower subpopulations compete with one another (Eqs. D–F). We used fixed values for energy diffusion and consumption rates as informed by model calibration (Fig. 3B) and varied sensitivity to low energy and chemotactic ability of both goer and grower, subject to Eqs. D–F (Table 1; Supplementary Table S1). We initially used the same spatial and temporal domains as during model calibration, but concluded that the implied duration of the experiment (3 days) was too short for dynamics between the two populations to manifest. Each MEMA spot has a low capacity, whereby confluence is reached at no more than a few hundreds of cells. Such a small number of cells will not exhibit wave-like behavior and, therefore, will not suffice for spatial structure to emerge. We, therefore, extended temporal and spatial domains of our simulations, seeding cells at a lower confluence and letting them grow onto the entire energy domain until they consumed all available energy (average of 127 days; Table 1).

We observed a nonmonotonic relation between the goer's chemotactic ability and the speed with which the metapopulation invades the dish, with intermediate values being the least beneficial to its growth and spread (Fig. 4C). Temporal analysis of the simulations (Supplementary Data S4–S6) revealed that if the goer's chemotactic motility is too high, it will leave the center of the dish too soon, leaving room for the grower to expand locally (Fig. 4D). In contrast, if the goer's motility is too low, it will miss the window of opportunity to ensure its dominance further away from the center of the dish while energy is still abundant. As a consequence, it will be outgrown by the grower at the edge of the dish once energy becomes sparse (Fig. 4E). Only when the goer has an intermediate motility, does the grower persistently coexist with it, both at the center and edge of the dish (Fig. 4F).

Models of infiltration are typically formulated under two critical assumptions. First, energy production and consumption are nonuniform, leading to the formation of an energy gradient (49–51); or second, energy consumption is very slow compared with production, leading to an essentially infinite energetic resource (52). Here, we formulated a generalized model of infiltration when energy is finite and investigated its behavior along a spectrum of scenarios, from permanent energy uniformity to scenarios where this uniformity is gradually lost. The model derivation does not assume a particular dimension (e.g., 2D in vitro experiments vs. in vivo or 3D spheroid experiments). Many parameters that were valid in 2D will also extend naturally to 3D. For example, we would not expect a difference in consumption rates or half maximal growth rates of the cells. However, energy diffusion (ΓE) or random cell motion (Γ) will be higher because of the increased degree of freedom (53). When energy is uniformly distributed at all times and the time scale for cell death is substantially longer than that of cell motility and birth, our results suggest that the degree of infiltration can be approximated using the cells' density at equilibration of movement and growth (Fig. 2C).

With an energy gradient that becomes steeper over time, our analytic approximations no longer hold, as directed cell movement becomes nonnegligible. For this scenario, we leveraged MEMA profiling to inform regions of interest in the parameter search space. These regions of interest are relevant for cellular growth on a variety of HGF-exposed ECM proteins. We observed correlations between inferred model parameters and RNA-seq–derived signatures, even though the latter were not used during parameter inference. A potential explanation for the negative correlation between ECM-specific energy consumption and expression is that our model does not account for the possibility that cells can replace the ECM they degrade. The slower the rate of this replacement is, the higher the consumption rate appears to be. On the other hand, the more dependent cells are on an ECM for growth, the faster they must replace it, potentially explaining a positive correlation between ECM-specific expression and sensitivity to low energy.

When calibrating our model to a given ECM environment, growth patterns of a substantial fraction of replicates of that ECM could not be explained by fixed choices of sensitivity to low energy (Supplementary Fig. S3). A potential explanation for this is variable cell compositions across experimental replicates. An alternative explanation is that this variability stems from artifacts that arise during nonuniform printing of ECMs onto the array—the so-called ring effect. However, a bimodal distribution was also observed in the DNA content of replicating cells, which is not affected by potential printing artifacts. The second peak of this bimodal distribution was wider, consistent with the fact that high-ploidy cells with more DNA need longer to replicate.

The cell line, HCC1954, is described as a hypertetraploid cell line with an average DNA content of 4.2 (30). However, this average value may be misleading, as suggested by stark variability in nuclei sizes (Fig. 4A). Despite a wealth of genomic information generated for this cell line (30), to the best of our knowledge, no prior reports indicate whether or not the cell line is polyploid. We and others have found that high ploidy is an aneuploidy-tolerating state that accompanies intratumor heterogeneity in vivo and in vitro (5, 54, 55). Our results suggest that HCC1954 is likely polyploid.

One event that could have led to this polyploid state is WGD. In contrast to cell lines, WGD events in primary tumors are mostly clonal, not subclonal (5, 6, 10), clones carrying a doubled genome often sweep over the population, such that by the time the tumor is detected, the diploid ancestor no longer exists. A related scenario is advanced, therapy-exposed tumors shown to revert to genomic stability, potentially bringing a WGD population back to a genomic state that more closely resembles its diploid ancestral state (56). The model presented here can investigate how dynamics between the two subpopulations unfold in both of these scenarios—early, shortly after the WGD or late, after therapy exposure. This would characterize what circumstances prevent the WGD carrying clone from becoming dominant or from retaining its dominance and could help explain WGD incidence in primary and recurrent tumors.

If spatial and temporal domains were to be extended beyond the configuration of MEMA spots, our simulations predict that spatial segregation of two coexisting subpopulations according to their ploidy is a likely scenario and depends on the energy consumption rate. Our model can easily be extended to more than two subpopulations, for example, to include a subpopulation of normal cells. For each additional cell type, a new compartment can be added to the model, with growth- and motion-related parameters that are specific to the corresponding cell type. However, the model currently does not incorporate mutations, that is, the process of generating new clonal lines. A next step will be to extend our model to include mutation events, specifically chromosome mis-segregations that contribute extensively to diversify ploidy of a population (57, 58). The additional DNA content of high-ploidy cells, although energetically costly, brings a masking effect against the deleterious consequences of chromosome losses (10). This duality may explain the higher sensitivity of high-ploidy cells to glycolysis inhibitors and their lower sensitivity to cytotoxic drugs reported previously in glioblastoma (59).

In-line with prior reports, we found that increased resistance of breast cancer cell lines to cytotoxic drugs was associated with high ploidy. In contrast, high-ploidy breast cancer cell lines were sensitive to inhibitors of signal transduction pathways, including EGFR and especially mTOR signaling. A commonality among those pathways is their contribution to a cell's chemotactic response (60–62), suggesting opportunities to tune chemotaxis. mTOR inhibitors (mTOR-I), such as rapamycin, significantly decrease migration of breast cancer cells in a dose-dependent manner (63, 64). Rapamycin inhibits cell motility by a mechanism similar to that by which it inhibits cell proliferation (65), suggesting that the mTOR pathway lies at the intersection of a cell's decision between proliferation and migration. If high ploidy is indeed a characteristic specific to goer-like cells, then mTOR-Is are likely affecting this cell type (Fig. 1A, E, and F), and could be used to inhibit its chemotactic response, thereby moving the population up the x-axis of Fig. 4C. Delaying chemotactic response of highly chemotactic cells could slow down invasion by maximizing competition within a polyploid population. If, on the other hand, chemotactic response of high-ploidy cells is already at an intermediate level, our simulation suggests that, further reduction may accelerate invasion of low-ploidy cells. For such scenarios, therapeutic strategies that include an mTOR-I may not be successful. Experiments will be needed to verify these in silico results in vitro. Knowing how coexisting clones with differential drug sensitivities segregate spatially can offer opportunities to administer these drug combinations more effectively.

No potential conflicts of interest were disclosed.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH or the H. Lee Moffitt Cancer Center and Research Institute

G.J. Kimmel: Conceptualization, software, formal analysis, validation, visualization, methodology, writing-original draft, writing-review and editing. M. Dane: Resources, validation, visualization, writing-review and editing. L.M. Heiser: Resources, data curation, validation, methodology. P.M. Altrock: Supervision, funding acquisition, visualization, writing-original draft, writing-review and editing. N. Andor: Conceptualization, data curation, supervision, funding acquisition, visualization, methodology, writing-original draft, writing-review and editing.

We thank the researchers of the Department of Integrated Mathematical Oncology at Moffitt Cancer Center for fruitful discussions and early feedback on parts of the article. This work was supported by the NCI grant R00CA215256 awarded to N. Andor. P.M. Altrock and N. Andor acknowledge support through the NCI, part of the NIH, under grant number P30-CA076292. L.M. Heiser acknowledges support through NIH research grants 1U54CA209988 and U54-HG008100, Jayne Koskinas Ted Giovanis Foundation for Health and Policy, and Breast Cancer Research Foundation.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Spiteri
I
,
Caravagna
G
,
Cresswell
GD
,
Vatsiou
A
,
Nichol
D
,
Acar
A
, et al
Evolutionary dynamics of residual disease in human glioblastoma
.
Ann Oncol
2019
;
30
:
456
63
.
2.
Brastianos
PK
,
Carter
SL
,
Santagata
S
,
Cahill
DP
,
Taylor-Weiner
A
,
Jones
RT
, et al
Genomic characterization of brain metastases reveals branched evolution and potential therapeutic targets
.
Cancer Discov
2015
;
5
:
1164
77
.
3.
Angelova
M
,
Mlecnik
B
,
Vasaturo
A
,
Bindea
G
,
Fredriksen
T
,
Lafontaine
L
, et al
Evolution of metastases in space and time under immune selection
.
Cell
2018
;
175
:
751
65
.
4.
Hannibal
RL
,
Chuong
EB
,
Rivera-Mulia
JC
,
Gilbert
DM
,
Valouev
A
,
Baker
JC
. 
Copy number variation is a fundamental aspect of the placental genome
.
PLoS Genet
2014
;
10
:
e1004290
.
5.
Dewhurst
SM
,
McGranahan
N
,
Burrell
RA
,
Rowan
AJ
,
Grönroos
E
,
Endesfelder
D
, et al
Tolerance of whole-genome doubling propagates chromosomal instability and accelerates cancer genome evolution
.
Cancer Discov
2014
;
4
:
175
85
.
6.
Bielski
CM
,
Zehir
A
,
Penson
AV
,
Donoghue
MTA
,
Chatila
W
,
Armenia
J
, et al
Genome doubling shapes the evolution and prognosis of advanced cancers
.
Nat Genet
2018
;
50
:
1189
95
.
7.
Amend
SR
,
Torga
G
,
Lin
KC
,
Kostecka
LG
,
de Marzo
A
,
Austin
RH
, et al
Polyploid giant cancer cells: unrecognized actuators of tumorigenesis, metastasis, and resistance
.
Prostate
2019
;
79
:
1489
97
.
8.
Pienta
KJ
,
Hammarlund
EU
,
Axelrod
R
,
Brown
JS
,
Amend
SR
. 
Poly-aneuploid cancer cells promote evolvability, generating lethal cancer
.
Evol Appl
2020
;
13
:
1626
34
.
9.
Zhang
S
,
Mercado-Uribe
I
,
Xing
Z
,
Sun
B
,
Kuang
J
,
Liu
J
. 
Generation of cancer stem-like cells through the formation of polyploid giant cancer cells
.
Oncogene
2014
;
33
:
116
28
.
10.
López
S
,
Lim
EL
,
Horswell
S
,
Haase
K
,
Huebner
A
,
Dietzen
M
, et al
Interplay between whole-genome doubling and the accumulation of deleterious alterations in cancer evolution
.
Nat Genet
2020
;
52
:
283
93
.
11.
Watson
SS
,
Dane
M
,
Chin
K
,
Tatarova
Z
,
Liu
M
,
Liby
T
, et al
Microenvironment-mediated mechanisms of resistance to HER2 inhibitors differ between HER2+ breast cancer subtypes
.
Cell Syst
2018
;
6
:
329
42
.
12.
Keenan
AB
,
Jenkins
SL
,
Jagodnik
KM
,
Koplev
S
,
He
E
,
Torre
D
, et al
The library of integrated network-based cellular signatures NIH program: system-level cataloging of human cells response to perturbations
.
Cell Syst
2018
;
6
:
13
24
.
13.
Andor
N
,
Lau
BT
,
Catalanotti
C
,
Sathe
A
,
Kubit
MA
,
Chen
J
, et al
Joint single cell DNA-seq and RNA-seq of gastric cancer cell lines reveals rules of in vitro evolution
.
NAR Genom Bioinform
2020
;
2
:
lqaa016
.
14.
Chun
YH
,
Kil
JI
,
Suh
YS
,
Kim
SH
,
Kim
H
,
Park
S-H
. 
Characterization of chromosomal aberrations in human gastric carcinoma cell lines using chromosome painting
.
Cancer Genet Cytogenet
2000
;
119
:
18
25
.
15.
Rege-Cambrin
G
,
Scaravaglio
P
,
Carozzi
F
,
Giordano
S
,
Ponzetto
C
,
Comoglio
PM
, et al
Karyotypic analysis of gastric carcinoma cell lines carrying an amplified c-met oncogene
.
Cancer Genet Cytogenet
1992
;
64
:
170
3
.
16.
Andor
N
,
Graham
TA
,
Jansen
M
,
Xia
LC
,
Aktipis
CA
,
Petritsch
C
, et al
Pan-cancer analysis of the extent and consequences of intratumor heterogeneity
.
Nat Med
2015
;
22
:
105
13
.
17.
Pirkmajer
S
,
Chibalin
AV
. 
Serum starvation: Caveat emptor
.
Am J Physiol Cell Physiol
2011
;
301
:
C272
9
.
18.
Lambert
K
,
Pirt
SJ
. 
Growth of human diploid cells (strain MRC-5) in defined medium; replacement of serum by a fraction of serum ultrafiltrate
.
J Cell Sci
1979
;
35
:
381
92
.
19.
Bartholomew
JC
,
Yokota
H
,
Ross
P
. 
Effect of serum on the growth of Balb oT3 A31 mouse fibroblasts and an SV40‐transformed derivative
.
J Cell Physiol
1976
;
88
:
277
86
.
20.
Cai
W
,
Rook
SL
,
Jiang
ZY
,
Takahara
N
,
Aiello
LP
. 
Mechanisms of hepatocyte growth factor-induced retinal endothelial cell migration and growth
.
Invest Ophthalmol Vis Sci
2000
;
41
:
1885
93
.
21.
Hatzikirou
H
,
Basanta
D
,
Simon
M
,
Schaller
K
,
Deutsch
A
. 
‘Go or grow’: the key to the emergence of invasion in tumour progression
?
Math Med Biol
2012
;
29
:
49
65
.
22.
Kathagen-Buhmann
A
,
Schulte
A
,
Weller
J
,
Holz
M
,
Herold-Mende
C
,
Glass
R
, et al
Glycolysis and the pentose phosphate pathway are differentially associated with the dichotomous regulation of glioblastoma cell migration versus proliferation
.
Neuro Oncol
2016
;
18
:
1219
29
.
23.
Dhruv
HD
,
McDonough Winslow
WS
,
Armstrong
B
,
Tuncali
S
,
Eschbacher
J
,
Kislin
K
, et al
Reciprocal activation of transcription factors underlies the dichotomy between proliferation and invasion of glioma cells
.
PLoS One
2013
;
8
:
e72134
.
24.
Vittadello
ST
,
McCue
SW
,
Gunasingh
G
,
Haass
NK
,
Simpson
MJ
. 
Examining go-or-grow using fluorescent cell-cycle indicators and cell-cycle-inhibiting drugs
.
Biophys J
2020
;
118
:
1243
7
.
25.
Wang
CH
,
Matin
S
,
George
AB
,
Korolev
KS
. 
Pinned, locked, pushed, and pulled traveling waves in structured environments
.
Theor Popul Biol
2019
;
127
:
102
19
.
26.
Bayliss
A
,
Volpert
VA
. 
Complex predator invasion waves in a Holling–Tanner model with nonlocal prey interaction
.
Physica D
2017
;
346
:
37
58
.
27.
Witelski
TP
,
Ono
K
,
Kaper
TJ
. 
On axisymmetric traveling waves and radial solutions of semi-linear elliptic equations
.
Nat Resour Model
2008
;
13
:
339
88
.
28.
Wasserstein
LN
. 
Markov processes over denumerable products of spaces describing large systems of automata
.
Problems Inform. Transmission
1969
;
5
:
47
52
.
29.
Sanderson
C
,
Curtin
R
. 
Armadillo: a template-based C++ library for linear algebra
.
J Open Source Softw
2016
;
1
:
26
.
30.
van der Meer
D
,
Barthorpe
S
,
Yang
W
,
Lightfoot
H
,
Hall
C
,
Gilbert
J
, et al
Cell model passports—a hub for clinical, genetic and functional datasets of preclinical cancer models
.
Nucleic Acids Res
2019
;
47
:
D923
9
.
31.
Barretina
J
,
Caponigro
G
,
Stransky
N
,
Venkatesan
K
,
Margolin
AA
,
Kim
S
, et al
The cancer cell line encyclopedia enables predictive modelling of anticancer drug sensitivity
.
Nature
2012
;
483
:
603
7
.
32.
Hafner
M
,
Niepel
M
,
Chung
M
,
Sorger
PK
. 
Growth rate inhibition metrics correct for confounders in measuring sensitivity to cancer drugs
.
Nat Methods
2016
;
13
:
521
7
.
33.
Hafner
M
,
Heiser
LM
,
Williams
EH
,
Niepel
M
,
Wang
NJ
,
Korkola
JE
, et al
Quantification of sensitivity and resistance of breast cancer cell lines to anti-cancer drugs using GR metrics
.
Sci Data
2017
;
4
:
170166
.
34.
Heiser
LM
,
Sadanandam
A
,
Kuo
W-L
,
Benz
SC
,
Goldstein
TC
,
Ng
S
, et al
Subtype and pathway specific responses to anticancer compounds in breast cancer
.
Proc Natl Acad Sci U S A
2012
;
109
:
2724
9
.
35.
Daemen
A
,
Griffith
OL
,
Heiser
LM
,
Wang
NJ
,
Enache
OM
,
Sanborn
Z
, et al
Modeling precision treatment of breast cancer
.
Genome Biol
2013
;
14
:
R110
.
36.
Bräutigam
K
,
Mitzlaff
K
,
Uebel
L
,
Köster
F
,
Polack
S
,
Pervan
M
, et al
Subtypes of triple-negative breast cancer cell lines react differently to eribulin mesylate
.
Anticancer Res
2016
;
36
:
2759
66
.
37.
Dai
X
,
Cheng
H
,
Bai
Z
,
Li
J
. 
Breast cancer cell line classification and its relevance with breast tumor subtyping
.
J Cancer
2017
;
8
:
3131
41
.
38.
Hänzelmann
S
,
Castelo
R
,
Guinney
J
. 
GSVA: gene set variation analysis for microarray and RNA-seq data
.
BMC Bioinformatics
2013
;
14
:
7
.
39.
Croft
D
,
Mundo
AF
,
Haw
R
,
Milacic
M
,
Weiser
J
,
Wu
G
, et al
The reactome pathway knowledgebase
.
Nucleic Acids Res
2014
;
42
:
D472
7
.
40.
Yang
W
,
Soares
J
,
Greninger
P
,
Edelman
EJ
,
Lightfoot
H
,
Forbes
S
, et al
Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells
.
Nucleic Acids Res
2013
;
41
:
D955
61
.
41.
Chen
WYJ
,
Abatangelo
G
. 
Functions of hyaluronan in wound repair
.
Wound Repair Regen
1999
;
7
:
79
89
.
42.
Ellis
IR
,
Schor
SL
. 
Differential effects of TGF-beta1 on hyaluronan synthesis by fetal and adult skin fibroblasts: Implications for cell migration and wound healing
.
Exp Cell Res
1996
;
228
:
326
33
.
43.
Fisher
RA
. 
The wave of advance of advantageous genes
.
Annals of Eugenics
1937
;
7
:
355
69
.
44.
Liotta
LA
. 
Tumor invasion and metastases–role of the extracellular matrix: Rhoads Memorial Award Lecture
.
Cancer Res
1986
;
46
:
1
7
.
45.
Hynes
RO
. 
Extracellular matrix: not just pretty fibrils
.
Science
2009
;
326
:
1216
9
.
46.
Becht
E
,
McInnes
L
,
Healy
J
,
Dutertre
C-A
,
Kwok
IWH
,
Ng
LG
, et al
Dimensionality reduction for visualizing single-cell data using UMAP
.
Nat Biotechnol
2018
Dec 3 [Epub ahead of print]
.
47.
Xu
C
,
Su
Z
. 
Identification of cell types from single-cell transcriptomes using a novel clustering method
.
Bioinformatics
2015
;
31
:
1974
80
.
48.
Sobol
IM
. 
Global sensitivity indices for nonlinear mathematical models
.
Mathematics and Computers in Simulation
. 
2001
;
55
:
271
80
.
49.
Lim
S
,
Nam
H
,
Jeon
JS
. 
Chemotaxis model for breast cancer cells based on signal/noise ratio
.
Biophys J
2018
;
115
:
2034
43
.
50.
Keller
EF
,
Segel
LA
. 
Model for chemotaxis
.
J Theor Biol
1971
;
30
:
225
34
.
51.
Anderson
ARA
,
Rejniak
KA
,
Gerlee
P
,
Quaranta
V
. 
Microenvironment driven invasion: a multiscale multimodel investigation
.
J Math Biol
2009
;
58
:
579
624
.
52.
Kimmel
GJ
,
Gerlee
P
,
Altrock
PM
. 
Time scales and wave formation in non-linear spatial public goods games
.
PLoS Comput Biol
2019
;
15
:
e1007361
.
53.
Einstein
A
. 
On the motion of small particles suspended in liquids at rest required by the molecular-kinetic theory of heat
. 
1905
;
17
:
208
.
54.
Andor
N
,
Simonds
EF
,
Czerwinski
DK
,
Chen
J
,
Grimes
SM
,
Wood-Bouwens
C
, et al
Single-cell RNA-seq of follicular lymphoma cancers reveals malignant B cell types and co-expression of T cell immune checkpoints
.
Blood
2019
;
133
:
1119
29
.
55.
Birkbak
NJ
,
Eklund
AC
,
Li
Q
,
McClelland
SE
,
Endesfelder
D
,
Tan
P
, et al
Paradoxical relationship between chromosomal instability and survival outcome in cancer
.
Cancer Res
2011
;
71
:
3447
52
.
56.
Morrissy
AS
,
Garzia
L
,
Shih
DJH
,
Zuyderduyn
S
,
Huang
X
,
Skowron
P
, et al
Divergent clonal selection dominates medulloblastoma at recurrence
.
Nature
2016
;
529
:
351
7
.
57.
Laughney
AM
,
Elizalde
S
,
Genovese
G
,
Bakhoum
SF
. 
Dynamics of tumor heterogeneity derived from clonal karyotypic evolution
.
Cell Rep
2015
;
12
:
809
20
.
58.
Elizalde
S
,
Laughney
AM
,
Bakhoum
SF
. 
A Markov chain for numerical chromosomal instability in clonally expanding populations
.
PLoS Comput Biol
2018
;
14
:
e1006447
.
59.
Donovan
P
,
Cato
K
,
Legaie
R
,
Jayalath
R
,
Olsson
G
,
Hall
B
, et al
Hyperdiploid tumor cells increase phenotypic heterogeneity within glioblastoma tumors
.
Mol Biosyst
2014
;
10
:
741
58
.
60.
Bailly
M
,
Wyckoff
J
,
Bouzahzah
B
,
Hammerman
R
,
Sylvestre
V
,
Cammer
M
, et al
Epidermal growth factor receptor distribution during chemotactic responses
.
Mol Biol Cell
2000
;
11
:
3873
83
.
61.
Gulhati
P
,
Bowen
KA
,
Liu
J
,
Stevens
PD
,
Rychahou
PG
,
Chen
M
, et al
mTORC1 and mTORC2 regulate EMT, motility, and metastasis of colorectal cancer via RhoA and Rac1 signaling pathways
.
Cancer Res
2011
;
71
:
3246
56
.
62.
James
RG
,
Davidson
KC
,
Bosch
KA
,
Biechele
TL
,
Robin
NC
,
Taylor
RJ
, et al
WIKI4, a novel inhibitor of tankyrase and Wnt/β-catenin signaling
.
PLoS One
2012
;
7
:
e50457
.
63.
Du
L
,
Li
X
,
Zhen
L
,
Chen
W
,
Mu
L
,
Zhang
Y
, et al
Everolimus inhibits breast cancer cell growth through PI3K/AKT/mTOR signaling pathway
.
Mol Med Rep
2018
;
17
:
7163
9
.
64.
Alalem
M
,
Ray
A
,
Ray
BK
. 
Metformin induces degradation of mTOR protein in breast cancer cells
.
Cancer Med
2016
;
5
:
3194
204
.
65.
Liu
L
,
Li
F
,
Cardelli
JA
,
Martin
KA
,
Blenis
J
,
Huang
S
. 
Rapamycin inhibits cell motility by suppression of mTOR-mediated S6K1 and 4E-BP1 pathways
.
Oncogene
2006
;
25
:
7029
40
.