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]

Major Findings

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.

Quick Guide to Equations and Assumptions
HDC model
Model Assumptions
Cell life cycle assumptions

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.

Figure 1.

HDC model. Cell life cycle (A) and ECM domain (B). See Quick Guide for details. * in Dn indicates (× 10−4).

Figure 1.

HDC model. Cell life cycle (A) and ECM domain (B). See Quick Guide for details. * in Dn indicates (× 10−4).

Close modal
ECM and cell movement assumptions

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.

Data Input Equations
Proliferation (M)

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).

Migration (Dn)

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.

Game Theory Model
Model assumptions

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.

Equations

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

$p=2c+h−1c.$
Equation (1):

The magnitude of this coexistence region is dependent on both mEI fitness and mE parameters and defined as (see Supplementary Appendix)

$1−h2≤c≤1−h$
Equation (2):

The cost for relying on the mE is traded against the cost for mE independence (related to fitness, h).

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.

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).

Table 1.

Parameter values for phenotypic traits

Cell line
pBabeneuNneuTAT1CA1aCA1d
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−4Dn 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 line
pBabeneuNneuTAT1CA1aCA1d
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−4Dn 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

*Doubling times in brackets were obtained using serum-free, growth factor–free cell culture conditions.

Value normalized to MCF10A-p120-knockdown.

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).

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.

See Quick Guide for diffusion coefficient calculations.

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 (1114). 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.

Figure 2.

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.

Figure 2.

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.

Close modal

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).

Figure 3.

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.

Figure 3.

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.

Close modal

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).

Figure 4.

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.

Figure 4.

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.

Close modal

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.

Figure 5.

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.

Figure 5.

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.

Close modal

### 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.

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 (3032). 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, 3340). 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.

No potential conflicts of interest were disclosed.

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.

1
Witz
IP
.
Yin-yang activities and vicious cycles in the tumor microenvironment
.
Cancer Res
2008
;
68
:
9
13
.
2
Anderson
AR
,
Weaver
AM
,
Cummings
PT
,
Quaranta
V
.
Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment
.
Cell
2006
;
127
:
905
15
.
3
Anderson
ARA
.
A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion
.
Math Med Biol
2005
;
22
:
163
86
.
4
Debnath
J
,
Muthuswamy
SK
,
Brugge
JS
.
Morphogenesis and oncogenesis of MCF-10A mammary epithelial acini grown in three-dimensional basement membrane cultures
.
Methods
2003
;
30
:
256
68
.
5
Santner
SJ
,
Dawson
PJ
,
Tait
L
, et al
.
Malignant MCF10CA1 cell lines derived from premalignant human breast epithelial MCF10AT cells
.
Breast Cancer Res Treat
2001
;
65
:
101
10
.
6
Huen
AC
,
Park
JK
,
Godsel
LM
, et al
.
Intermediate filament-membrane attachments function synergistically with actin-dependent contacts to regulate intercellular adhesive strength
.
J Cell Biol
2002
;
159
:
1005
17
.
7
Clark
ES
,
Whigham
AS
,
Yarbrough
WG
,
Weaver
AM
.
Cortactin is an essential regulator of matrix metalloproteinase secretion and extracellular matrix degradation in invadopodia
.
Cancer Res
2007
;
67
:
4227
35
.
8
Bryce
NS
,
Clark
ES
,
Leysath
JL
,
Currie
JD
,
Webb
DJ
,
Weaver
AM
.
Cortactin promotes cell motility by enhancing lamellipodial persistence
.
Curr Biol
2005
;
15
:
1276
85
.
9
Potdar
AA
,
Lu
J
,
Jeon
J
,
Weaver
AM
,
Cummings
PT
.
Bimodal analysis of mammary epithelial cell migration in two dimensions
.
Ann Biomed Eng
2009
;
37
:
230
45
.
10
Bargmann
CI
,
Hung
MC
,
Weinberg
RA
.
Multiple independent activations of the neu oncogene by a point mutation altering the transmembrane domain of p185
.
Cell
1986
;
45
:
649
57
.
11
Hotary
KB
,
Allen
ED
,
Brooks
PC
,
Datta
NS
,
Long
MW
,
Weiss
SJ
.
Membrane type I matrix metalloproteinase usurps tumor growth control imposed by the three-dimensional extracellular matrix
.
Cell
2003
;
114
:
33
45
.
12
Sabeh
F
,
Ota
I
,
Holmbeck
K
, et al
.
Tumor cell traffic through the extracellular matrix is controlled by the membrane-anchored collagenase MT1-MMP
.
J Cell Biol
2004
;
167
:
769
81
.
13
Artym
VV
,
Zhang
Y
,
Seillier-Moiseiwitsch
F
,
KM
,
Mueller
SC
.
Dynamic interactions of cortactin and membrane type 1 matrix metalloproteinase at invadopodia: defining the stages of invadopodia formation and function
.
Cancer Res
2006
;
66
:
3034
43
.
14
Weaver
AM
.
Invadopodia: specialized cell structures for cancer invasion
.
Clin Exp Metastasis
2006
;
23
:
97
105
.
15
Friedl
P
,
Wolf
K
.
Tube travel: the role of proteases in individual and collective cancer cell invasion
.
Cancer Res
2008
;
68
:
7247
9
.
16
Provenzano
PP
,
Eliceiri
KW
,
Campbell
JM
,
Inman
DR
,
White
JG
,
Keely
PJ
.
Collagen reorganization at the tumor-stromal interface facilitates local invasion
.
BMC Med
2006
;
4
:
38
.
17
Gatenby
RA
,
Smallbone
K
,
Maini
PK
, et al
.
Cellular adaptations to hypoxia and acidosis during somatic evolution of breast cancer
.
Br J Cancer
2007
;
97
:
646
53
.
18
Gerlee
P
,
Anderson
ARA
.
An evolutionary hybrid cellular automaton model of solid tumour growth
.
J Theor Biol
2007
;
246
:
583
603
.
19
Fisher
RA
.
.
Ann Eugenics
1937
;
7
:
355
69
.
20
Hanahan
D
,
Weinberg
RA
.
The hallmarks of cancer
.
Cell
2000
;
100
:
57
70
.
21
Smalley
KS
,
Brafford
PA
,
Herlyn
M
.
Selective evolutionary pressure from the tissue microenvironment drives tumor progression
.
Semin Cancer Biol
2005
;
15
:
451
9
.
22
Kam
Y
,
Karperien
A
,
Weidow
B
,
L
,
Anderson
AR
,
Quaranta
V
.
Nest expansion assay:a cancer systems biology approach to in vitro invasion measurements
.
BMC Res Notes
2009
;
2
:
130
.
23
Maynard Smith
J
.
Evolution and the theory of games
.
Cambridge; New York
:
Cambridge University Press
;
1982
.
24
Basanta
D
,
Deutsch
A
.
A game theoretical perspective on the somatic evolution of cancer
. In:
Bellomo
N
, editor.
Selected topics on cancer modeling: genesis, evolution, immune competition, therapy
.
Boston
:
Birkhauser
;
2008
.
25
Vogelstein
B
,
Kinzler
KW
.
The multistep nature of cancer
.
Trends Genet
1993
;
9
:
138
41
.
26
Reya
T
,
Clevers
H
.
Wnt signalling in stem cells and cancer
.
Nature
2005
;
434
:
843
50
.
27
Arakawa
H
.
Netrin-1 and its receptors in tumorigenesis
.
Nat Rev Cancer
2004
;
4
:
978
87
.
28
Murray-Zmijewski
F
,
Slee
EA
,
Lu
X
.
A complex barcode underlies the heterogeneous response of p53 to stress
.
Nat Rev Mol Cell Biol
2008
;
9
:
702
12
.
29
Staib
F
,
Robles
AI
,
Varticovski
L
, et al
.
The p53 tumor suppressor network is a key responder to microenvironmental components of chronic inflammatory stress
.
Cancer Res
2005
;
65
:
10255
64
.
30
Chung
LW
,
Baseman
A
,
Assikis
V
,
Zhau
HE
.
Molecular insights into prostate cancer progression: the missing link of tumor microenvironment
.
J Urol
2005
;
173
:
10
20
.
31
Kalluri
R
,
Zeisberg
M
.
Fibroblasts in cancer
.
Nat Rev Cancer
2006
;
6
:
392
401
.
32
Rubin
H
.
Microenvironmental regulation of the initiated cell
.
2003
;
90
:
1
62
.
33
Dormann
S
,
Deutsch
A
.
Modeling of self-organized avascular tumor growth with a hybrid cellular automaton
.
In Silico Biol
2002
;
2
:
393
406
.
34
Drasdo
D
,
Hohme
S
.
A single-cell-based model of tumor growth in vitro: monolayers and spheroids
.
Phys Biol
2005
;
2
:
133
47
.
35
Jiang
Y
,
Pjesivac-Grbovic
J
,
Cantrell
C
,
Freyer
JP
.
A multiscale model for avascular tumor growth
.
Biophys J
2005
;
89
:
3884
94
.
36
Rejniak
KA
.
A single cell approach in modeling the dynamics of tumor microregions
.
Math Biosci Eng
2005
;
2
:
643
55
.
37
Stott
E
,
Britton
NF
,
Glazier
JA
,
Zajac
M
.
Stochastic simulation of benign avascular tumor growth using the Potts model
.
J Theor Biol
1999
;
30
:
183
98
.
38
Turner
S
,
Sherratt
JA
.
Intercellular adhesion and cancer invasion: a discrete simulation using the extended Potts model
.
J Theor Biol
2002
;
216
:
85
100
.
39
Zhang
L
,
Athale
CA
,
Deisboeck
TS
.
Development of a three-dimensional multiscale agent-based tumor model: simulating gene-protein interaction profiles, cell phenotypes and multicellular patterns in brain cancer
.
J Theor Biol
2007
;
244
:
96
107
.
40
Poplawski
NJ
,
Agero
U
,
Gens
JS
,
Swat
M
,
Glazier
JA
,
Anderson
AR
.
Front instabilities and invasiveness of simulated avascular tumors
.
Bull Math Biol
2009
;
71
:
1189
227
.
41
Jones
S
,
Chen
WD
,
Parmigiani
G
, et al
.
Comparative lesion sequencing provides insights into tumor evolution
.
Proc Natl Acad Sci U S A
2008
;
105
:
4283
8
.
42
Beerenwinkel
N
,
Antal
T
,
Dingli
D
, et al
.
Genetic progression and the waiting time to cancer
.
PLoS Comput Biol
2007
;
3
:
e225
.
43
Bernards
R
,
Weinberg
RA
.
A progression puzzle
.
Nature
2002
;
418
:
823
.
44
Ramis-Conde
I
,
Drasdo
D
,
Anderson
AR
,
Chaplain
MA
.
Modeling the influence of the E-cadherin-β-catenin pathway in cancer cell invasion: a multiscale approach
.
Biophys J
2008
;
95
:
155
65
.
45
Rejniak
KA
.
An immersed boundary framework for modelling the growth of individual cells: an application to the early tumour development
.
J Theor Biol
2007
;
247
:
186
204
.
46
Rejniak
KA
,
Anderson
ARA
.
A computational study of the development of epithelial acini: I. Sufficient conditions for the formation of a hollow structure
.
Bull Math Biol
2008
;
70
:
677
712
.
47
Chang
TH
,
Szabo
E
.
Enhanced growth inhibition by combination differentiation therapy with ligands of peroxisome proliferator-activated receptor-γ and inhibitors of histone deacetylase in adenocarcinoma of the lung
.
Clin Cancer Res
2002
;
8
:
1206
12
.
48
Nowak
D
,
Stewart
D
,
Koeffler
HP
.
Differentiation therapy of leukemia—three decades of development
.
Blood
2009
;
113
:
3655
65
.