Abstract
Tumor-microenvironment interactions are increasingly recognized to influence tumor progression. To understand the competitive dynamics of tumor cells in diverse microenvironments, we experimentally parameterized a hybrid discrete-continuum mathematical model with phenotypic trait data from a set of related mammary cell lines with normal, transformed, or tumorigenic properties. Surprisingly, in a resource-rich microenvironment, with few limitations on proliferation or migration, transformed (but not tumorigenic) cells were most successful and outcompeted other cell types in heterogeneous tumor simulations. Conversely, constrained microenvironments with limitations on space and/or growth factors gave a selective advantage to phenotypes derived from tumorigenic cell lines. Analysis of the relative performance of each phenotype in constrained versus unconstrained microenvironments revealed that, although all cell types grew more slowly in resource-constrained microenvironments, the most aggressive cells were least affected by microenvironmental constraints. A game theory model testing the relationship between microenvironment resource availability and competitive cellular dynamics supports the concept that microenvironmental independence is an advantageous cellular trait in resource-limited microenvironments. [Cancer Res 2009;69(22):8797–806]
These results suggest that a critical feature of the process of tumor progression is selection of cells that can escape from resource limitations by achieving a relative microenvironmental independence.
Each cell can move, proliferate, degrade ECM, and adhere to other cells on the two-dimensional lattice. The first step in the cell life cycle (Fig. 1) is examination of the ECM concentration and the number of immediate neighbors. Next, depending on the phenotype, cells proliferate, become quiescent, or migrate. The six cell phenotypes, along with nondimensionalized parameters for phenotypic traits (M, doubling time; Dn, migration; Ai, cell-cell adhesion; δ, ECM degradation), are shown in colored boxes.
HDC model. Cell life cycle (A) and ECM domain (B). See Quick Guide for details. * in Dn indicates (× 10−4).
HDC model. Cell life cycle (A) and ECM domain (B). See Quick Guide for details. * in Dn indicates (× 10−4).
As shown in Fig. 1, ECM is modeled as a heterogeneous mixture with pixels modeled at the cell scale (diameter h) and with concentrations from 0% to 100% (heat map). When space constraints are modeled, the “ECM threshold” represents the concentration below which cells are allowed to move or proliferate. Cells (Ni) can modify the ECM concentration (f) through localized degradation (δ). Cells can only move one step to their four orthogonal neighboring lattice sites with probability P1 (left), P2 (right), P3 (down), or P4 (up) or remain stationary (P0) each time step (k). Further detailed description of all HDC assumptions is given in Supplementary Table S1.
We scaled calculated doubling times by the timescale (τ = 24 h) multiplied by the simulation time step (k = 0.00005) to obtain the final value. Thus, the proliferation parameter M = doubling time/(τ × k) [i.e., pBabe in full serum, M = 22.36/(24 × 0.00005) ∼ 18,600 time steps]. Time units are 1/10th of the timescale (e.g., 2.4 h).
To obtain diffusion coefficients, we used the relationship Dn = (ts2)/2, where s is the experimentally measured speed of motion for each cell line for a given time t(9). Dn has units of μm2/h. To nondimensionalize this parameter, we scaled it by the time (24 h) and length scales (1 cm2) of the model to obtain Dn = s2 × 12/108. To implement unbiased migration from a diffusion coefficient, we used the relationship between a random walk and the diffusion equation (Supplementary Table S1). This relationship allows definition of movement probabilities for cells modulated by diffusion coefficients (Fig. 1) and has been formally defined in the HDC model (10) and for reinforced random walks (11).
To scale experimental dissociation values, normalized measurements were divided into four equally spaced bins derived from setting the highest dissociation value (p120-knockdown = 100.0%) to 0, and the lowest percent dissociation value (AT1 = 10.73%) to 3, to define outer limits of the data range. The four established bins were set to 0 (77.69–100%), 1 (55.37–77.68%), 2 (33.05–55.36%), and 3 (10.73–33.04%); all phenotypes were then assigned an integer (between 0 and 3) according to the appropriate bin.
We assume that the tumor is made of cells that do (mED) or do not (mEI) depend on their microenvironment for survival. The interactions between mEI and mED cells can be summarized as follows: Although the maximum fitness for any cell is 1, mEI cells are assumed to have a constant fitness (h < 1), regardless of the presence of other cells. Conversely, mED cells change their fitness if they are surrounded by mEI cells (which do not take resources from the mE) versus other mED cells (which do compete for those resources). Thus, the change in fitness of a mED cell interacting with mEI cells is 1 − c, whereas this value is 1 − nc if an mED cell interacts with other mED cells. c represents the cost of relying on the mE (c ∼ 0 in a microenvironment with easy access to growth factors, space, or other resources); n represents the increase in microenvironmental dependence cost due to available mE resources being shared between multiple mED cells.
The fitness of microenvironment-independent (mEI) and microenvironment-dependent (mED) cells is as follows: h, fitness of mEI cells and is independent of cellular interactions; c, cost of extracting resources from the microenvironment; n, a constant that accounts for greater extraction of resources by various (more than one) mED competing cells.
To calculate the proportion p of mEI cells in a population with both types of cells, we assume that in equilibrium, the fitness of the two populations is equal (see Supplementary Appendix for details). Assuming that the cost of relying on the mE is twice as much in environments in which more than one mED cell has to share available resources (n = 2), then the proportion p of mEI cells is
The magnitude of this coexistence region is dependent on both mEI fitness and mE parameters and defined as (see Supplementary Appendix)
The cost for relying on the mE is traded against the cost for mE independence (related to fitness, h).
Introduction
The microenvironment can either suppress or promote tumor progression, depending on the context (1). Due to the complex, multivariate nature of tumor-microenvironment interactions, traditional reductionist approaches may not be sufficient to achieve a complete mechanistic understanding of the process. A complementary method to probe tumor progression mechanisms is mathematical and computational modeling, to integrate and test in silico interactions between tumor cells and various aspects of the microenvironment.
We previously described a hybrid discrete-continuum (HDC) mathematical model that drives computer simulations of tumor growth by integrating cellular and microenvironmental variables (2, 3). In the HDC model, individual cells with different phenotypic traits interact spatially with each other and with microenvironmental variables, such as extracellular matrix (ECM) or oxygen. Analysis of evolving simulated tumors can provide mechanistic understanding by allowing in silico testing of factors that affect tumor progression. For example, the interaction of cell phenotypes with “harsh” microenvironments (limiting oxygen, heterogeneous ECM) was found to produce tumors with invasive morphology, dominated by few cancer cell phenotypes with aggressive traits (high proliferation, low adhesion, high motility; ref. 2). In contrast, under “mild” microenvironmental conditions (spatially homogeneous, abundant oxygen), there was little clonal selection and noninvasive morphologies.
In our previous work, theoretical phenotypes were parameterized with diverse values from the literature (2, 3). However, an important step is to incorporate experimental data gathered under comparable conditions to connect simulation predictions with biological context. For this purpose, we collected a homogeneous experimental phenotypic trait data set from a panel of related breast epithelial cell lines with varying aggressiveness. Surprisingly, differences between measured phenotypic traits of normal, transformed, and tumorigenic cell lines were much smaller than those previously used for the theoretical cell phenotypes. Nonetheless, in simulated tumors run with growth factor–deprived and space-constrained microenvironments, tumorigenic phenotypes became dominant. Conversely, simulations run with a nonlimiting microenvironment instead favored transformed cell types. Analysis of the performance of single-phenotype tumors in diverse microenvironments revealed that all cell types grew more slowly in resource-limited environments; however, the most aggressive cells were least affected. A nonzero sum game theory analysis confirms the general advantage conferred by microenvironmental independence for cells competing in resource-constrained microenvironments. These results suggest that tumorigenesis is an evolving process for which the increasing acquisition of independence from microenvironmental constraints is a critical feature.
Materials and Methods
MCF10A parent, pBabe, neuN, and neuT were from Dr. Joan Brugge (Harvard Medical School, Boston, MA), which were authenticated by three-dimensional culture and cultured as described (4); MCF10AT1, CA1a and CA1d were from the original source Dr. Fred Miller (Karmanos Institute, Detroit, MI; ref. 5). Experimental raw data are shown in Supplementary Figs. S1 to S3.
Extraction of model parameters from experimental data, scaling, and fitting were as follows.
Proliferation
Data were obtained from in vitro experiments with cells grown under sparse, exponential growth conditions. In vivo growth rates would likely be difficult to use in simulations because they represent not only intrinsic cell doubling times but also responses to contact inhibition and unknown factors in the tumor microenvironment (e.g., growth factors and hypoxia). Experimental growth data (Supplementary Fig. S1) were assessed for a best fit model using the “curve estimation tool” in GraphPad Prism. Exponential growth fits were significantly better than other tested models (R2 values between 0.904 and 0.995). From experimental data, we fitted exponential growth curves of the form y = aegt using SPSS, where y is cell number, a is a scaling constant, g is growth rate, and t is time (in hours). Doubling times (Table 1) were calculated from ln2/g (i.e., pBabe in full serum, g = 0.031, ln2/0.031 = 22.36 h). See Quick Guide for model input of doubling times (M).
Parameter values for phenotypic traits
. | . | Cell line . | |||||
---|---|---|---|---|---|---|---|
pBabe . | neuN . | neuT . | AT1 . | CA1a . | CA1d . | ||
Tumorigenic | No | No | No | Rarely | Yes | Yes | |
Experimental parameters | |||||||
Proliferation* (h) | M | 22.36 [—] | 17.34 [30.14] | 17.33 [33.01] | 18.73 [24.76] | 15.75 [24.76] | 16.12 [21.66] |
Unbiased migration (μm/h) ± SE | s | 46.96 ± 1.41 | 58.56 ± 1.62 | 51.06 ± 1.70 | 47.12 ± 1.99 | 44.59 ± 1.49 | 57.26 ± 2.56 |
ECM degradation/cell area (10−3) ± SE (10−3) | δ | 0.32 ± 0.09 | 4.90 ± 1.32 | 4.59 ± 3.17 | 2.75 ± 1.57 | 5.82 ± 1.45 | 12.07 ± 3.98 |
Cell-cell adhesion (% dissociation)† ± SE | Ai | 34.70 ± 0.78 | 33.37 ± 2.63 | 49.10 ± 2.86 | 10.73 ± 1.08 | 13.46 ± 1.41 | 14.29 ± 1.77 |
Nondimensional parameters | |||||||
Proliferation* (time steps) | M | 18,600 [—] | 14,400 [25,117] | 14,400 [27,508] | 15,600 [20,633] | 13,200 [20,633] | 13,400 [18,050] |
Unbiased migration (× 10−4) | Dn | 2.64 | 4.11 | 3.99 | 2.66 | 2.38 | 3.93 |
ECM degradation | δ | 0.0016 | 0.0245 | 0.023 | 0.0138 | 0.0291 | 0.0604 |
Cell-cell adhesion index | Ai | 2 | 2 | 1 | 3 | 3 | 3 |
. | . | Cell line . | |||||
---|---|---|---|---|---|---|---|
pBabe . | neuN . | neuT . | AT1 . | CA1a . | CA1d . | ||
Tumorigenic | No | No | No | Rarely | Yes | Yes | |
Experimental parameters | |||||||
Proliferation* (h) | M | 22.36 [—] | 17.34 [30.14] | 17.33 [33.01] | 18.73 [24.76] | 15.75 [24.76] | 16.12 [21.66] |
Unbiased migration (μm/h) ± SE | s | 46.96 ± 1.41 | 58.56 ± 1.62 | 51.06 ± 1.70 | 47.12 ± 1.99 | 44.59 ± 1.49 | 57.26 ± 2.56 |
ECM degradation/cell area (10−3) ± SE (10−3) | δ | 0.32 ± 0.09 | 4.90 ± 1.32 | 4.59 ± 3.17 | 2.75 ± 1.57 | 5.82 ± 1.45 | 12.07 ± 3.98 |
Cell-cell adhesion (% dissociation)† ± SE | Ai | 34.70 ± 0.78 | 33.37 ± 2.63 | 49.10 ± 2.86 | 10.73 ± 1.08 | 13.46 ± 1.41 | 14.29 ± 1.77 |
Nondimensional parameters | |||||||
Proliferation* (time steps) | M | 18,600 [—] | 14,400 [25,117] | 14,400 [27,508] | 15,600 [20,633] | 13,200 [20,633] | 13,400 [18,050] |
Unbiased migration (× 10−4) | Dn | 2.64 | 4.11 | 3.99 | 2.66 | 2.38 | 3.93 |
ECM degradation | δ | 0.0016 | 0.0245 | 0.023 | 0.0138 | 0.0291 | 0.0604 |
Cell-cell adhesion index | Ai | 2 | 2 | 1 | 3 | 3 | 3 |
*Doubling times in brackets were obtained using serum-free, growth factor–free cell culture conditions.
†Value normalized to MCF10A-p120-knockdown.
Cell-cell adhesion
Experimental cell-cell adhesion measurements were obtained using in vitro dispase assays (Supplementary Fig. S3; ref. 6) and are presented as average percentage of cell-cell dissociation, where 100% is equivalent to complete dissociation after pipetting. Because MCF10A-p120-knockdown cells exhibit very little cell-cell adhesion (∼96% dissociation), they were used as a positive control. Experimental dissociation values were normalized against p120-knockdown values to achieve recalculated percent dissociations (i.e., out of 100%).
Cell-cell adhesion (Ai) is modeled by designating a value for each cell that represents its likelihood to adhere to other cells. Movement is restricted based on these values and on the number of neighboring cells (see Supplementary Table S1 and refs. 2, 3 for details). These internal cell-cell adhesion values are represented by four integers, ranging from 0 (no adhesion) to 3 (maximal adhesion), which must be derived from experimental values before input (see Quick Guide).
ECM degradation
Experimental ECM degradation measurements were obtained using an invadopodia assay (Supplementary Fig. S2; ref. 7) and presented as average degradation area per cell area. Because a cell is represented as a single lattice site, the cell area should match that of the model and the timescale is reasonably close. However, in simulations using these raw values, ECM degradation occurred slowly. Thus, to minimize computation time, we used an arbitrary scaling of 5 to enhance the degradation ability (δ) of each cell type. Control simulations run under various ECM constraints without this scaling and for a 5-fold longer time period gave the same clonal dominance ranking, but a slightly more ruffled edge due to greater effect of ECM heterogeneity (results not shown).
Unbiased migration
Single-cell velocities (μm/h; Table 1) were measured as described (8) and converted into diffusion rates (Dn) by using the relationship between cell velocity and diffusion. We find that cells, under these conditions, follow a random walk, as defined by a linear relationship between mean squared displacement and time (9), and by 4 h enter a fully diffusive regimen.4
4M. Harris, P.L. Frick, D. Tyson, S. Garbett, B. Weidow, E. Koeler, Y. Shyr, and V. Quaranta. Unpublished data.
Results
To experimentally parameterize the HDC model, we adopted a panel of cultured cell lines (Table 1) based on the nontumorigenic, immortalized breast epithelial cell line MCF10A (4). These cell lines were produced by introduction of genes known to be involved in breast cancer and, in some cases, by passaging the cells through mice: pBabe (vector control); neuN (overexpressed normal c-neu oncogene); neuT (oncogenic mutant of neu; ref. 10); AT1 (activated ras); CA1a and CA1d (ras plus multiple in vivo passages). pBabe is immortalized but nontransformed, neuN and neuT are transformed but not tumorigenic, AT1 is poorly tumorigenic, and CA1a and CA1d are consistently tumorigenic and metastatic in tail vein injection assays (5). Thus, these cell lines can be taken to represent different tumor progression stages.
Experimental measurements and data integration
Key parameters in the HDC model include rates of proliferation (M), ECM degradation (δ), motility (Dn), and cell-cell adhesion values (Ai; Table 1; Fig. 1; ref. 2). These parameters were measured in each cell line under conventional tissue culture conditions. Raw experimental values and detailed graphs and images are shown in Table 1 and Supplementary Figs. S1 to S3, respectively. Interestingly, the measured trait values do not immediately distinguish cell lines with respect to known tumorigenicity (Table 1). Only ECM degradation was higher in the consistently tumorigenic CA1a and CA1d cell lines compared with other cell lines. Proliferation in growth factor–deprived media was the highest in the CA1d cells but was equivalent between the strongly and weakly tumorigenic CA1a and AT1 cells, respectively.
To introduce these values into the model, the data were analyzed and nondimensionalized (Materials and Methods). Each simulated “cell” is endowed with a collection of experimentally derived phenotypic traits and can interact with other cells (through cell-cell adhesion) and the microenvironment (Fig. 1; Supplementary Table S1). Cells can move and proliferate if there is adjacent lattice space and if their internal cell-cell adhesion value is less than the number of cellular neighbors. This latter rule gives cells with low cell-cell adhesion values a higher probability of moving away from neighboring cells than cells with high cell-cell adhesion values (ref. 3; Supplementary Table S1). Movement at each time step is also modulated by the migration speed of each phenotype. Migration and proliferation are considered as continuous processes with one following the other at each time step (e.g., first a cell moves then it proliferates) as constrained by cell phenotype (Table 1) and availability of adjacent lattice sites.
Incorporation of experimental data inevitably led to model changes. We originally used two partial differential equations to model ECM degradation with cellular production and diffusion of proteases coupled to ECM disappearance. The direct measurement of cell-associated ECM degradation obviated the need for separate equations and made the process more realistic because transmembrane proteases such as membrane type 1 matrix metalloproteinase are necessary and sufficient for cell invasion and ECM degradation (11–14). In addition, haptotaxis and oxygen consumption are not modeled in this study. It is unclear whether haptotaxis is a critical component of tumor cell motility because it has only been measured in vitro, and in tissues, proteolytic activity is more likely to result in rearrangement of ECM fibers than generation of smooth haptotactic gradients (15, 16). Oxygen consumption is not modeled because it became apparent that it should be treated as a separate study with measurement of metabolic parameters and incorporation of glycolytic switch equations (17, 18). Cell death is not explicitly modeled, but rather implicitly modeled using net cell doubling times, experimentally gathered in the presence or absence of growth factors.
Transformed and tumorigenic phenotypes coexist in a resource-rich microenvironment
To examine tumor population dynamics in a resource-rich microenvironment, simulations of heterogeneous tumor growth were run in which equal numbers of the six experimentally parameterized cell types were mixed together. We define a microenvironment as “resource-rich” when the microenvironment imposes no limitations; that is, cells compete only for space on the lattice. Analysis of tumor cell dynamics shows that transformed neuT and neuN phenotypes become dominant over time, with some competition from the tumorigenic CA1d phenotype. By contrast, the normal pBabe and respectively weakly and strongly tumorigenic AT1 and CA1a cells are quickly suppressed by other phenotypes (Fig. 2A). Similar trends are evident in single-phenotype simulations (Fig. 2B), with the most “fit” phenotypes being neuT, neuN, and CA1d. However, small differences in fitness are critical in the competition for space during heterogeneous tumor development (Fig. 2A). Normalization of cell-cell adhesion values to either 1 (neuT value; Fig. 2C) or 3 (CA1d value; not shown) leads to a shift in the order of phenotype dominance, with CA1d > neuN > neuT. In all simulations, AT1, CA1a, and pBabe were least fit regardless of cell-cell adhesion, indicating lower fitness with regard to other traits.
Tumor evolution in a resource-rich microenvironment. A, heterogeneous tumor simulations. Left, screen shot from representative simulation depicts the spatial phenotype distribution within the tumor at t = 120 (time units). Coloration of tumor cells represents phenotype, as indicated. Right, fraction of each phenotype in the population as a function of time. B, single-phenotype tumor simulations. Left, sample screen shots (t = 120) from representative simulations. Cell coloration reflects proliferating (green) or quiescent (blue) cells. Right, individual growth curves from single-phenotype tumors. C, heterogeneous tumor simulations with equalized cell-cell adhesion values (Ai = 1). Left, sample screen shot depicting spatial distribution of phenotypes at t = 120. Right, tumor fraction plot.
Tumor evolution in a resource-rich microenvironment. A, heterogeneous tumor simulations. Left, screen shot from representative simulation depicts the spatial phenotype distribution within the tumor at t = 120 (time units). Coloration of tumor cells represents phenotype, as indicated. Right, fraction of each phenotype in the population as a function of time. B, single-phenotype tumor simulations. Left, sample screen shots (t = 120) from representative simulations. Cell coloration reflects proliferating (green) or quiescent (blue) cells. Right, individual growth curves from single-phenotype tumors. C, heterogeneous tumor simulations with equalized cell-cell adhesion values (Ai = 1). Left, sample screen shot depicting spatial distribution of phenotypes at t = 120. Right, tumor fraction plot.
In the resource-rich microenvironment, tumor growth in HDC simulations depends heavily on migration and proliferation. The classic work of Fisher (19) considers a gene spreading through a population as an invading wave travelling at constant speed and similarly driven by diffusion and proliferation. In the case where only migration and proliferation differentially affect cellular behavior (Fig. 2C, equalized cell-cell adhesion), it is possible to calculate Fisher wave speeds and obtain the same ranking of phenotype dominance (results not shown). However, this relationship breaks down when cell-cell adhesion and degradation are modifying factors and phenotypes are placed under microenvironmental constraints.
Comparison of simulation results with known biological behavior leads to an obvious discrepancy: neuT and neuN cells are not tumorigenic yet occupy the largest fraction of simulated tumors (together ∼80%; Fig. 2A). Because the simulations were run under environmental conditions without resource limitations, and our previous work implicated microenvironmental selection in tumor invasiveness (2), we hypothesized that a more realistic microenvironment might alter competitive cellular dynamics during tumor evolution. In general, tumor microenvironments are considered to be stressful, with limited availability of growth factors, nutrients, or oxygen (20, 21), as well as spatial constraints likely to limit tumor growth (11). To test the effects of resource limitation on tumor dynamics and composition, we modeled growth factor deprivation and ECM constraints, as detailed below.
Implementation of resource limitation promotes dominance of the tumorigenic CA1d cells
To model the effect of growth factor deprivation on heterogeneous tumor evolution, growth factor–deprived cell doubling times (bracketed parameters for M, Table 1) were used in simulations along with the other phenotypic parameters. The pBabe phenotype was unable to grow due to lack of proliferation. AT1 and CA1a phenotypes initially grew well; however, their levels declined after time unit ∼40 due to increasing dominance of the CA1d phenotype and “trapping” inside the tumor. neuN and neuT each maintained a steady fraction of 16% and 12% in the growing tumor, respectively (Fig. 3A).
Growth factor– or space-constrained environments promote dominance of tumorigenic CA1d phenotype. A, tumor simulations in a growth factor–deprived microenvironment. Left, sample screen shot (t = 120). Right, fraction of each phenotype in tumor population as a function of time. B, tumor simulations in a space-constrained microenvironment. ECM constraints are modeled by designating a maximum ECM concentration that allows cell movement or proliferation. Left, screen shots (t = 120) under growth factor–rich conditions and increasing ECM constraints (98%, 90%, 80%, and 60%). Note transition in tumor composition with increasing space constraints. Right, line graphs show cell number of each phenotype in the population at t = 140 as a function of ECM constraint. C, tumor simulations with equalized ECM degradation. Left, sample screen shots (t = 120) from simulated tumors in which ECM degradation values of all phenotypes were set equal to that of CA1d. Right, line graph of cell number at t = 140 for each phenotype in the population under various ECM constraints.
Growth factor– or space-constrained environments promote dominance of tumorigenic CA1d phenotype. A, tumor simulations in a growth factor–deprived microenvironment. Left, sample screen shot (t = 120). Right, fraction of each phenotype in tumor population as a function of time. B, tumor simulations in a space-constrained microenvironment. ECM constraints are modeled by designating a maximum ECM concentration that allows cell movement or proliferation. Left, screen shots (t = 120) under growth factor–rich conditions and increasing ECM constraints (98%, 90%, 80%, and 60%). Note transition in tumor composition with increasing space constraints. Right, line graphs show cell number of each phenotype in the population at t = 140 as a function of ECM constraint. C, tumor simulations with equalized ECM degradation. Left, sample screen shots (t = 120) from simulated tumors in which ECM degradation values of all phenotypes were set equal to that of CA1d. Right, line graph of cell number at t = 140 for each phenotype in the population under various ECM constraints.
In comparison with the resource-rich microenvironment, growth factor deprivation leads to an earlier and larger degree of clonal dominance. Thus, the tumorigenic CA1d phenotype occupies 63% of the tumor at the simulation end, with no other phenotype occupying more than 16% (Fig. 3A). By contrast, under unconstrained conditions, neuN, neuT, and CA1d coexist to a greater extent, with respective tumor fractions of 50%, 35%, and 12% (Fig. 2A).
To model tissue space constraints, the movement and proliferation of cells were limited by ECM concentrations above a certain threshold. ECM is represented in a generalized manner, without specification of molecular details or structure (Fig. 1; refs. 2, 3). Cellular ECM degradation lowers the ECM level at a given lattice site, allowing movement or proliferation if the new ECM level falls below the threshold. Although the ECM threshold is the same for all cells within a given simulation, each cell type has a different ability to degrade ECM based on the measured trait parameter δ (Table 1; Fig. 1B). At stringent thresholds (e.g., ECM concentrations ≤80% allow cell movement or proliferation), tumors become smaller and more compact (Fig. 3B). Furthermore, neuT no longer dominates, giving way to neuN at intermediate (90–86% ECM threshold) and to CA1d under high (≤84% ECM threshold) spatial constraints. Thus, similar to the growth factor–deprived environment, imposition of high space constraints leads to dominance of CA1d (Fig. 3B), with suppression of other phenotypes. Similarly, in recently published in vitro colony expansion assays, we find that CA1d colonies expand further than MCF10A colonies in the presence, but not in the absence, of space constraints, and CA1d cells dominate the periphery of mixed MCF10A:CA1d space-constrained cultures (22). Similar to cell-cell adhesion in the resource-rich environment, ECM degradation accounts for the dominance switch between neuN/neuT/CA1d under high ECM constraints, as equalizing cellular degradation rates leads to coexistence of neuN/neuT/CA1d under high ECM constraints (Fig. 3C).
Tumor progression in a dually constrained microenvironment
In vivo, tumors grow and progress under multiple microenvironmental constraints. In an environment that is both growth factor and space constrained, CA1d once again becomes dominant with the highest degree of suppression of other phenotypes (Fig. 4A; ∼90% of the tumor population in growth factor–deprived + highest ECM constraint, compared with ∼60% of the tumor population in the growth factor–deprived + lowest ECM constraint). This finding validates the aggressiveness of the tumorigenic CA1d phenotype and indicates that microenvironment constraints can combine to further drive clonal dominance and tumor progression. However, this is not an entirely surprising result because the CA1d phenotype is more successful in each of the separate growth factor–deprived and space-constrained environments (Fig. 3).
Dually constrained microenvironment drives tumor progression. A, tumor simulations under ECM constraints (threshold: 98%, 90%, 80%, 60%) and growth factor deprivation. Population dynamics for each simulation are plotted in a line graph under the appropriate spatial image (t = 120). B, tumor simulations in the absence of CA1d with ECM constraints and growth factor deprivation. Note that although CA1a had no obvious advantage in the presence of CA1d (A), it now dominates simulations in the 60% ECM constraint.
Dually constrained microenvironment drives tumor progression. A, tumor simulations under ECM constraints (threshold: 98%, 90%, 80%, 60%) and growth factor deprivation. Population dynamics for each simulation are plotted in a line graph under the appropriate spatial image (t = 120). B, tumor simulations in the absence of CA1d with ECM constraints and growth factor deprivation. Note that although CA1a had no obvious advantage in the presence of CA1d (A), it now dominates simulations in the 60% ECM constraint.
To test whether tumor evolution would occur differently in the absence of CA1d cells, simulations were run with a mixture of the remaining five phenotypes under growth factor deprivation and in the presence or absence of space constraints (Fig. 4B). Interestingly, under minimal space constraints (90–98% ECM threshold), AT1 and CA1a initially dominate equally but are overtaken quickly by the better-migrating neuN and neuT phenotypes. Under intermediate space constraints (80% threshold), AT1 and CA1a gain an increasing share of the tumor, to the disadvantage of neuT. However, as stringent space constraints are imposed (≤60% threshold), CA1a becomes dominant, holding a >50% fraction of the tumor. This clonal dominance is likely due to the slight advantage that CA1a holds in ECM degradation capacity (Table 1) over the other phenotypes. Thus, in the absence of CA1d, the next phenotype that is predicted to dominate under highly selective microenvironmental conditions is the only other consistently tumorigenic cell line, CA1a. These findings suggest that a key feature of tumorigenicity is the ability to survive under adverse circumstances, and that the process of tumor progression selects for cells that can do so.
Remarkably, CA1a was never competitive for dominance in simulations when all six phenotypes were mixed (compare Fig. 4B to Figs. 2, 3, and 4A). Instead, the main competition for dominance appeared to take place between CA1d and the neuN/neuT phenotypes, which constituted the largest fraction of the tumor. Thus, when both CA1a and CA1d are present, the similarity between phenotypes with a slight advantage to CA1d may lead to intense competition within the same tumor niche and suppression of CA1a by its next of kin, CA1d.
Independence from microenvironmental constraints as an aggressive cancer cell trait
To understand why resource constraints give an advantage to tumorigenic phenotypes, we analyzed the relative fitness of each phenotype in resource-constrained environments. To do so, we divided the final cell number at the end of resource-constrained simulations by the final cell number at the end of resource-rich simulations (microenvironment independence values, shown underneath appropriate simulation images in Fig. 5A). To isolate microenvironment effects, we analyzed single-phenotype simulation results. Interestingly, although all cells grow slowly in resource-constrained microenvironments, normal cells are severely affected and tumorigenic cells are least affected by resource restriction (Fig. 5A). Thus, in the dually constrained microenvironment, the aggressive CA1d phenotype achieves 38% of its performance in the resource-rich environment. By contrast, normal pBabe cells perform at 1% in this comparison and the other cells are intermediate, with a slight advantage to CA1a. In singly constrained environments, the best separation between tumorigenic and transformed cell lines occurs with ECM constraint, consistent with invasion as a late stage in tumor progression and the association of growth factor independence with early stages of cellular transformation. Overall, these results suggest that microenvironmental independence is a critical cellular feature acquired during the process of tumor progression.
Tumorigenic cells exhibit relative independence from microenvironment constraints. Single-phenotype simulations were run in each of the four microenvironments, as indicated. A, final spatial images from each simulation. Green, proliferating cells; blue, quiescent cells. Analysis of the microenvironmental independence of each phenotype was done by dividing the final number of cells in simulations in resource-constrained environments (as indicated) by the number of cells in the resource-rich environment (e.g., growth factor–rich, no space constraint). These ratios of cell numbers in constrained environment/rich environment are shown for each appropriate environment underneath the spatial images. A fully microenvironment-independent cell should have a value of 1 for its relative performance, whereas a very microenvironment-dependent cell will have a value close to 0. Although in the dually constrained microenvironment (bottom row) no cell type achieves a value above 0.5, there is a correlation of cellular aggressiveness with relative microenvironmental independence (e.g., compare 0.38 for CA1d to 0.01 for pBabe). B and C, game theoretical model results. B, two-dimensional plots show analysis of the long-term tumor composition: The proportion of mEI or mED cells is plotted as a function of cost (c) when very unfit (h = 0.1; left) and fit (h = 0.9; right) mEI phenotypes are considered and n = 2. C, three-dimensional plot shows variance of mED cell fraction with mEI cell fitness and environment cost in the long term.
Tumorigenic cells exhibit relative independence from microenvironment constraints. Single-phenotype simulations were run in each of the four microenvironments, as indicated. A, final spatial images from each simulation. Green, proliferating cells; blue, quiescent cells. Analysis of the microenvironmental independence of each phenotype was done by dividing the final number of cells in simulations in resource-constrained environments (as indicated) by the number of cells in the resource-rich environment (e.g., growth factor–rich, no space constraint). These ratios of cell numbers in constrained environment/rich environment are shown for each appropriate environment underneath the spatial images. A fully microenvironment-independent cell should have a value of 1 for its relative performance, whereas a very microenvironment-dependent cell will have a value close to 0. Although in the dually constrained microenvironment (bottom row) no cell type achieves a value above 0.5, there is a correlation of cellular aggressiveness with relative microenvironmental independence (e.g., compare 0.38 for CA1d to 0.01 for pBabe). B and C, game theoretical model results. B, two-dimensional plots show analysis of the long-term tumor composition: The proportion of mEI or mED cells is plotted as a function of cost (c) when very unfit (h = 0.1; left) and fit (h = 0.9; right) mEI phenotypes are considered and n = 2. C, three-dimensional plot shows variance of mED cell fraction with mEI cell fitness and environment cost in the long term.
Game theory model
To investigate further the role of microenvironmental independence, we used a game theoretical model that studies how different microenvironments (mEs, defined by cost of obtaining resources from them) explain tumor population composition. Game theoretical models have been used to study evolutionary dynamics in biological populations (23), including cancer (24). Model assumptions and equations are described in Quick Guide and Supplementary Appendix.
The game theoretical analysis shows that mEs with a low cost of reliance favor mED cells (Fig. 5B and C). However, resource-limited mEs favor mEI cells. The two-dimensional plots in Fig. 5B show the expected proportion of mED or mEI cells as a function of varying cost (c) for extracting resources. Note the region of coexistence of mEI and mED populations [which is broader for low than for high fitness (h) values] as the cost of obtaining microenvironment resources increases, similar to the crossover that occurs in the plot of tumor cell number versus ECM threshold in Fig. 3B. In both cases, mEI cells are not as fit as mED cells in a low-cost microenvironment; however, within a high-cost microenvironment, they are comparatively fitter. The three-dimensional plot in Fig. 5C shows how the proportion of mED cells in the tumor population varies based on the cost (c) of relying on external mE factors and the fitness (h) of mEI cells.
Discussion
To understand the role of the microenvironment in tumor progression, we gathered experimental parameter data from breast epithelial cells with normal, transformed, and tumorigenic properties for integration into the HDC model. We find that tumorigenic cells are predicted to dominate simulated tumors only under resource-constrained (growth factor–deprived or space-constrained) conditions. Analysis of tumor growth of each cell type in diverse simulated environments revealed that all phenotypes grew less well with limited resources; however, biological aggressiveness correlated with a relative resistance to microenvironmental constraints. A game theory analysis suggests that, in general, the ability to achieve microenvironment independence without a large fitness loss is advantageous to cells and likely to be selected for in high-cost environments. In aggregate, our results suggest that independence from microenvironmental constraints is a critical feature of tumor progression.
Microenvironmental independence as a key feature of carcinogenesis
Microenvironmental independence as a defining feature of cancer is consistent with cancer as a disease of deregulation, in which proliferation is unchecked by normal controls. Thus, it is not rampant proliferation that characterizes cancer cells (as is frequently assumed), but rather the ability to continue proliferating despite microenvironmental context that would inhibit normal cells. In considering the “hallmarks of cancer” (20), almost all of them can be connected to microenvironmental independence. Cell-intrinsic changes that allow cells to ignore their microenvironment include evading apoptosis, self-sufficiency in growth signals, and insensitivity to antigrowth signals. Active modification of the microenvironment to make it more favorable includes angiogenesis and tissue invasion (20). In light of our study, we speculate that the acquisition of microenvironment-independent traits should correlate with each step of tumor progression. Indeed, aggressive metastatic cells are extremely microenvironment independent, colonizing diverse tissue sites.
In colon cancer, a specific series of genetic changes (the “vogelgram”) was shown to correlate with the progression from a normal cell to advanced carcinoma and includes loss of APC, activation of k-ras, loss of DCC, and loss of p53 genes (25). As opposed to a model in which genetic alterations remove intrinsic cell barriers to tumorigenesis that are unrelated to local microenvironment, our study is instead consistent with a model in which microenvironment-dependent selection accounts for the frequency of tissue-specific genetic changes. In light of that concept, how might one consider the well-defined “vogelgram”? APC mutations usually represent the first hit in colon tumorigenesis and decouple transit-amplifying and crypt base cell responses from extracellular Wnt ligand concentrations by upregulating β-catenin levels (26). Activation of the signaling protein k-ras occurs most frequently at the transition from early to intermediate adenomas and allows growth factor independence (20, 25). DCC is frequently lost at the intermediate to large adenoma transition and is a “dependence receptor” that promotes apoptosis when not bound to its ligand netrin (25, 27). Thus, loss of DCC confers microenvironment independence to evolving adenomas by allowing survival in areas of low netrin (e.g., at the tip of intestinal villi; ref. 27). Loss of p53 function is a late event in many cancers, occurring at the late adenoma to carcinoma transition, and regulates cell cycle arrest, apoptosis, and cellular senescence. Although many p53 functions appear to be cell intrinsic (e.g., cell cycle arrest in response to DNA damage), p53 is important in regulating cellular responses to many microenvironmental stresses including hypoxia and osmotic stress (25, 28, 29).
Recent studies have shown that normal host cells can promote tumor progression (30–32). Our findings are not inconsistent with that view. We would argue that cancer cells in fact do not evolve independently from their microenvironment, but rather that a normal tissue environment provides a variety of microenvironmental barriers that must be overcome during the process of tumorigenesis and progression. Whether those barriers are overcome entirely by cancer cells themselves or through the recruitment and activation of professional cells, such as inflammatory cells or fibroblasts, is a matter of strategy. The creation of a permissive environment may be more easily achieved via helper cells such as activated fibroblasts, which can lay down matrix, secrete angiogenic and growth factors, and/or remodel matrix (31), than through multiple changes in the carcinoma cell itself. Over time, such a permissive environment may allow fixation of stable genetic changes in cancer cells (e.g., progression).
Modeling approach
Many individual-based models of tumor growth and progression exist (18, 33–40). We focused on the current model because it fits easily into current biological and experimental frameworks. The theoretical focus of the HDC model on cell-microenvironment interactions is relevant to current debates on mechanisms of tumor progression (32), clonal expansion (25, 41), rates of mutation acquisition (42), and metastasis (43). Model limitations are primarily associated with representation of cells as lattice points, such that ECM is not given elastic or force properties, molecular specificity, or tissue architecture, and cell-cell adhesion is represented in an abstract manner, without cell geometry and associated forces. Modeling approaches that more naturally represent these elements include the lattice-based cellular Potts model (38, 40), the lattice-free, center-based model (34, 44), and the immersed boundary model (45). Nonetheless, the model has other advantages, including easy incorporation of individual phenotypes and experimental data, which allow analysis of competition between cellular phenotypes. Future inclusion of more accurate representations of ECM structure and mechanics and cell shape could be achieved by linking with models that already have these in place, albeit with fewer cells (46).
Implications and future directions
If microenvironmental independence is a critical, defining feature of cancer, therapies that reverse this process might be effective. One possibility is differentiation therapy to limit the cellular response range by inducing maturation (e.g., all-trans retinoic acid for promyelocytic leukemia or potentially histone deacetylase inhibitors +/− peroxisome proliferator–activated receptor-γ ligands for lung adenocarcinomas; refs. 47, 48). Another approach is inhibition of genes or processes that promote adaptability to the microenvironment.
Future testing of microenvironmental independence as a driving force for tumor progression will require characterizing cell performance in additional microenvironments, including acidic, hypoxic, inflammatory, or drug treated, along with in vivo validation of simulation results. Another future direction could include investigation of the role of tissue architecture and other tissue-specific microenvironment barriers in driving tumor progression. We anticipate that future integration of experimentation with modeling will continue to provide insight into the process of tumor progression.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Acknowledgments
Grant support: U54CA113007 (V. Quaranta) and K22CA109590 (A.M. Weaver).
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.