Impact of Metabolic Heterogeneity on Tumor Growth, Invasion, and Treatment Outcomes

Tumor heterogeneity at the genetic scale has been known for decades and until recently was largely viewed as a whole tumor metric. Historically, molecular techniques average genomic signals from large numbers of cells from single biopsies, thus smoothing and potentially hiding underlying variations. These average signatures have dominated the molecular era of cancer investigation and have driven biomarker development. However, a potential issue was recently highlighted by Gerlinger and colleagues (9), who showed that multiple biopsies from the same tumor display distinct genetic profiles and yet are phenotypically similar. This genotypic divergence and phenotypic convergence has previously been hinted at theoretically (10) and may be a predictable evolutionary consequence of the tumor ecosystem (11, 12). The intricate dialogue between tumor cells and environment selects for clones that are best adapted phenotypically to survive, regardless of specific mutations that may facilitate tumor progression. Furthermore, this environment is temporally and spatially heterogeneous largely due to variations in blood flow, resulting in local fluctuations of nutrients, growth factors, and other cellular populations (e.g., normal cells, stromal cells, and immune cells). These dynamics occurring within the cancer ecosystem are almost impossible to dissect via experimentation alone.

A major open question in biology is how to connect genotype with phenotype. In many genetic diseases, a defect within a single gene manifests through observable changes that usually allows confident diagnosis of the specific mutation based on the phenotypic manifestations. In contrast, cancer cells frequently harbor thousands of mutations so that the connection to phenotype is less clear. Although we understand that the genotype to phenotype mapping is critical for integration of genetic information with a more functional understanding of the phenotypic behavior that it may facilitate, the practical reality is that it is very difficult to quantify phenotypes and far easier to measure genotypes. Genetic heterogeneity within a tumor continues to be of great interest to the cancer community, with the advent of new tools to measure fewer cells in more detail, a serious effort is being put into quantifying this heterogeneity and understanding how it evolves as the tumor progresses and how it relates to overall outcome (13–15). However, we are far from understanding how the microenvironment modulates this heterogeneity and drives the overall phenotypic behavior of the tumor cell population. We strongly believe that only through the integration of mathematical and computational models with careful experimentation can we hope to bridge the gene-centric and microenvironment-centric views of cancer progression. Here, we focus on the phenotypic scale, as it allows us to circumvent much of the complexity observed at molecular scales and instead examine the functional outcome of mutation (be it genetic or epigenetic). Specifically, we consider the metabolic phenotype, in part, because tumors are known to be metabolically very different from their normal tissue counterparts and, in part, because we believe metabolism contributes significantly to tumor progression.

Mathematical models of cancer progression have examined avascular growth (3, 16, 17), angiogenesis (18, 19), vascular growth (20, 21), invasion (22, 23), and metastasis (24, 25); also see these reviews (17, 25–27). Tumor metabolism, specifically glycolytic metabolism and its development has also been extensively studied (2, 28). However, almost all of these models (with a few exceptions) ignore normal tissue and assume the tumor grows in a field of nutrients or extracellular matrix. Fewer models explicitly consider phenotypic heterogeneity (22, 29) and those that do again ignore normal tissue. Our central focus here is to examine the interplay between cellular metabolism and the microenvironment (including normal cells), and how this leads to phenotypic heterogeneity within the tumor cell population. Understanding how this phenotypic heterogeneity develops through space and time (30) will allow us to better understand progression, treatment failure, and crucially lay the groundwork for novel therapies.

In our evolutionary model, tumors emerge from tightly regulated fields of cells in a largely normal tissue. At its simplest, homeostatic regulation of normal tissue ensures that cell death is balanced by birth. This homeostasis is not static but dynamic and remarkably robust to physical and chemical perturbation, for example, wounding. Cancer cells both overcome and exploit these homeostatic mechanisms as they evolve strategies to maximize their proliferation (31).

We have previously examined metabolism (2, 28, 29) as a driver of tumor progression through the glycolytic pathway and associated acidosis. Extracellular pH (pH_e) is a potential treatment target, and we have recently published a series of experiments that validate this view (4, 32, 33). Our intention here is to examine this work in light of the lens of metabolic and microenvironmental heterogeneity by reproducing observation and suggesting novel interpretation and treatment strategies. The focus of the work is on acid-mediated invasion, wherein tumor cells that have abnormally high glycolytic capacity acidify the local environment, causing the surrounding tissue to degrade and allow for tumor growth.

To investigate the effect of metabolism and the microenvironment on tumor growth and phenotype selection, a hybrid cellular automaton (HCA) in two dimensions is used to model cellular processes such as cell proliferation, death, and vascular remodeling. Metabolic activity of individual cells is calculated by solving a set of algebraic equations. The diffusion of molecules is calculated by solving a system of partial differential equations over the domain of simulation. The components of the model are described in the Quick Guide to Model and Major Assumptions. Additional derivations and supporting algorithms of the model equations are given in the Supplementary Data.

In the absence of tumor cells, the dynamic equilibrium of the model is a homeostatic normal tissue with vasculature. The normal cells turnover, maintaining a cellular density that is determined by the choice of parameters for the proliferation algorithm. This homeostasis plays an important role in the model, competing with tumor cells for available space at the tumor edge. To illustrate the homeostatic properties of the normal tissue, two examples are shown in the Supplementary Data.

There are two modes of tumor growth, depending on the activity of the normal tissue surrounding the tumor. The first is homeostatic-limited growth, where tumor growth is limited by the surrounding normal cells. Because the tumor can only grow into space relinquished by the normal cells at the tumor–normal interface, the growth rate is correlated with the normal tissue homeostatic death rate. Tissue that has a faster turnover rate relinquishes more empty space per unit of time, and therefore allows for faster tumor growth (see the Supplementary Data). The second mode of growth is invasive growth. This occurs when the tumor mass causes surrounding normal cells to die at a rate faster than the homeostatic death rate. This results in rapid availability of free space and subsequent tumor growth. There are two possible causes for induced cell death in the present model: nutrient deprivation and acidosis.

Figure 3A–H shows a sequence of images taken from a representative simulation of tumor growth. In Fig. 3A and B, the average tumor phenotype is still relatively normal. A slight shift toward acquired acid resistance is due to intermittent hypoxia caused by blood vessel degradation followed by angiogenesis within the center of the tumor. Each short period of hypoxia-induced acidosis causes a slight selection pressure toward acid-resistant cells.

In Fig. 3C and D, a sustained hypoxic central area leads to increased reliance on glycolysis by the tumor cells growing there. Because cells must develop resistance to increasing acidosis or perish, the central part of the tumor is acid-resistant (blue). This resistant population is surrounded by a metabolically normal population of tumor cells (green). Because they are near to stable vasculature in the normal tissue, this surrounding buffer of metabolically benign cells exists because there is little selection pressure to acquire acid resistance on the tumor edge.

Figure 3E and F shows the development of necrosis in the center of the tumor, and a few small pockets of glycolytic cells (magenta). The increased glycolysis in these cells causes a further decrease in pH_e, but the cells are not in sufficient number nor adjacent to acid-sensitive cells to allow for acid-mediated invasion to occur. There continues to be a buffer of tumor cells (green) between the aggressive phenotype and the acid-susceptible normal tissue.

Figure 3G and H shows the tumor growing invasively via acid-mediated invasion. The acidosis produced by glycolytic cells affects the surrounding normal tissue. Several lobes of very aggressive glycolytic and acid-resistant cells have punctured through the edges of the tumor, growing into the space left open by normal tissue destroyed by acidosis. A thin zone of black empty space between the tumor and normal tissue is evident at the edge of the largest lobe on the lower right, a phenomenon seen in histologic samples by Gatenby and colleagues (28).

This representative simulation offers several insights into the development of acid-mediated invasion. The growth curve for the simulation is shown in Fig. 3I. The initial tumor growth has a shallow slope, representing homeostatic-limited growth. From t = 0 to 1,200 days, the average growth rate is about 0.2 mm per year. Here, the tumor takes advantage of space relinquished by normal cells during homeostatic turnover. The growth rate is correspondingly slow. Eventually, the tumor becomes invasive via self-produced acidosis and from t = 1,800 to t = 1,920 days, the average growth rate is about 6 mm per year, a 30-fold change over the homeostatic-limited growth rate. Interestingly, there is no jump in phenotype that causes this sudden change in behavior. Unlike a genetic mutation model, where one could point to the acquisition of a key mutation shortly before a large behavioral change, this model only permits slow phenotype variation. Indeed, comparing the phenotype flow plots before and after invasive growth (third and fourth rows) shows only a slight shift in average phenotype. To understand this change in behavior, it is important to consider not only the specific phenotypes present, but also where they arise spatially and how they develop dynamically. Besides acquiring the invasive phenotype, three other key factors play a role in the development of acid-mediated invasion.

First, the cells must reach a critical mass of acid-producing cells. A single glycolytic, acid-resistant cell cannot sufficiently change the local pH_e to cause lethal acidosis to surrounding cells. Because the model assumes that phenotypes are passed on to daughter cells with some small random change, reaching critical mass is primarily a function of proliferating enough to create a group of similarly aggressive cells in the local microenvironment that can then have sufficient proton production, as a group, to significantly affect the pH_e of the surrounding tissue.

The second key development is that the invasive cells must be near to cells that are relatively acid-sensitive if invasion is to occur. If a highly glycolytic group of cells develops but is surrounded by acid-resistant, low-glycolytic cells, then the lowered pH_e that is caused by the glycolytic cells will have little effect on the surrounding cells and invasion will not occur. Therefore, the importance of the spatial buffer of relatively normal tumor cells on the edge of the tumor becomes clear. Unless the aggressive cells can invade through this protective barrier en masse, the acid-invasive phenotype will be contained interior to the tumor edge. To state it another way, the diffusion gradient of acidosis produced by a group of glycolytic cells must overcome the spatially structured phenotype gradient that exists toward the edge of the tumor.

Finally, the presence of nearby blood vessels can contribute to decreased acidosis. A group of aggressive tumor cells that grows into a well vascularized region will be slowed by the increased buffering capacity of the vasculature. Invasion will only occur if the low pH_e produced by the mass of cells locally overcomes the ability of the vasculature to remove the acidosis through buffering.

The pathway to becoming a metabolically aggressive cell is not a linear march through phenotype space. The diagonal line on the phenotype plots of Fig. 3B, D, F, and H shows the direction of equal, maximal phenotypic variation. However, simulations show that the cells reaching the most aggressive phenotype tend to follow a nonlinear path similar to the black arrow shown in Fig. 2B. This path shows that acid resistance is selected first, after which high-glycolytic capacity develops. This observation is the opposite of the hypothesis that cells first acquire a glycolytic phenotype followed by acquisition of acid resistance (2), represented by the gray arrow of Fig. 2B. The reason for the order of selection in the present model is that low pH_e naturally develops as a result of decreased or intermittent vascular density in the center of a growing tumor. Increased acid resistance is therefore a valuable trait to acquire early in the development of the tumor if vascular delivery is inconsistent or reduced. Although increased glycolysis speeds up the cell cycle, this phenotype change has little advantage during homeostatic growth, because available space is freed on a time scale that is much longer than the cell-cycle time. Furthermore, unless critical mass is reached, acid-mediated invasion by an individual glycolytic cell is not possible.

Interestingly, because of the spatial structure of hypoxia and cell turnover, the most aggressive cells tend to develop in interior regions of the tumor where nutrient conditions are poor, selection pressures are great, and proliferation is more frequent because of cell death. However, to become invasive, these aggressive cells must come in contact with normal tissue, where the excess proton production can cause cell death. The ability of a mass of aggressive cells to penetrate the surrounding layer of less aggressive tumor cells depends, in part, on how acid-resistant the buffering cells are. If the phenotype gradient is steep, then the aggressive cells are more likely to invade and reach the normal tissue than if the intervening cells are already acid-resistant. Clearly, vascular density is an important aspect of the progression of the tumor phenotypes in both space and time. Further exploration of the effect of vascular density is presented in the Supplementary Data.

The acid-mediated invasion seen in the model has been observed experimentally in a murine window-chamber model (4). In that experiment, a droplet of tumor cells was cultured and implanted into a dorsal wound chamber. The pH_e and tumor growth rates were measured radially, and the direction of maximal growth was correlated with the direction of low pH_e, as seen in Fig. 3J. A similar correspondence between pH_e at the tumor edge and radial growth was observed in the model, as seen in Fig. 3K.

This representative simulation is produced consistently by the model with the same parameters. Additional repeats of the simulation are shown in the Supplementary Data, including further exploration of the impact of normal cell turnover rate and normal cell density. Simulations were also performed with larger CA domains (up to 1 cm) to explore slower rates of phenotypic drift (see the Supplementary Data). Although the tumor takes more time to become invasive with a slow phenotypic drift rate, the qualitative nature of the evolution through phenotype space is the same as for simulations run in a smaller domain with faster evolution as described above. Because large domains become computationally expensive, we limit our present results and analysis to the CA domains of approximately 3 to 5 mm per side.

Recent experimental results from our group have shown that application of systemic pH buffers to mice can prevent or slow down tumor development and metastatic growth (32). In a murine TRAMP model of spontaneous tumor growth, sodium bicarbonate has a preventative effect on tumor formation in these mice, but only if it is administered before a certain age (33).

A generic pH buffering therapy has been simulated using the model (see the Supplementary Data for discussion of the assumptions for implementation). To compare the simulation with the murine buffer study (33), an untreated control simulation is compared with same simulation treated at different time points (t_rel days). Figure 4 compares the untreated control with two buffer-treated simulations. The tumor that is treated earlier (Fig. 4B, E, H, and K, middle) shows a significant delay in growth and acquisition of an invasive tumor phenotype compared with the untreated control (Fig. 4A, D, G, and L, left). On the other hand, the tumor treated at a later time point (Figure 4C, F, I, and L, right) is only partly affected by the buffer. There is some slowing of the growth, but the invasive phenotype persists and the tumor growth rate is comparable with the untreated case. Supplementary Fig. S11 shows the growth of the untreated tumor (black plot) and three different starting times for buffer therapy. The earlier the treatment is administered, the more delay of growth is observed.

These results suggest a possible mechanism for the preventative effect of early buffer treatment seen in the TRAMP model. If the therapy is administered before the protective layer of metabolically benign cells being breached (i.e., acid-mediated invasion has commenced), then the therapy can prevent the aggressive cells from reaching this invasive state for long periods. If the invasive state is already reached, however, the therapy has only a minimal effect. To test this hypothesis, we suggest that tumors from early and late treated mice be examined for metabolic signatures and phenotypic structure that could indicate the presence of acid-mediated invasion. The model also suggests that buffer therapy may be of limited use to patients with a primary tumor that is already glycolytic and exploiting acid-mediated invasion. However, buffer therapy may still be beneficial to these patients because micrometastases that have yet to develop an invasive phenotype may be prevented by this systemic therapy. Even if the metastatic seed was highly glycolytic, if the micr-metastasis has not reached the critical mass needed to induce acid-mediated invasion, then the buffer therapy will further delay the acquisition of that invasive state by effectively raising the necessary threshold of critical mass. These results are also consistent with the TRAMP mouse experiments where late-treated mice, in contrast to untreated mice, did not develop metastasis (4).

Antiangiogenic therapy was developed on the premise that if a tumor is starved of nutrients, the tumor will grow slower, become quiescent, and perhaps die away. Although this seems like a straightforward response to decreased nutrient availability, this nutrient-deprivation theory does not take into account the effect of waste products produced as a result of worsened metabolic conditions within a tumor. Specifically, reduction of perfusion with antiangiogenic therapy will cause pockets of hypoxia, which, in turn, will lead to increased glycolysis in all cells and increased acidosis. The result of this sequence of events is that cells that may have previously developed acid resistance due to an avascular tumor growth phase may suddenly be at an advantage when vascular density is decreased. The selection pressure due to therapy could have the effect of speeding up the development of the most aggressive subset of tumor cells at the time of treatment.

To illustrate this, we compare a simulation without any therapy with one in which antiangiogenesis was started at a given time step. For simplicity, we assume that antiangiogenic therapy simply turns off the angiogenic aspect of the model: after the initiation of therapy, no new blood vessels are created. Existing blood vessels remain, and vessel degradation by the growing tumor continues. As expected, this leads to decreased vessel density within the tumor, so that nutrients will primarily diffuse from the surrounding vasculature within the normal tissue.

Figure 5 shows the comparison between treated and untreated cases for a representative simulation. The effect of the therapy is to greatly increase the selection pressure for acid-resistant cells and speed evolution by increased cellular turnover. This pressure drives the phenotype of the tumor cells toward a more aggressive state, and therefore acid-mediated invasion occurs much sooner than in the untreated simulation. The addition of therapy in this case is clearly a worse outcome than allowing the tumor to progress without therapeutic intervention. This result suggests a possible explanation for the historical failure of many antiangiogenic therapies (34), including the recent withdrawals of antiangiogenic therapies for both glioblastoma (35) and breast cancer (36). It is of interest to note that shortly after therapy is applied, the number of tumor cells decreases, as cells that are far from the remaining blood vessels will be starved of nutrients. However, as the remaining cells adapt to the new environment, the tumor will grow back and become invasive sooner than the control.

The diffusible nature of therapies, such as chemotherapy, has implications for the spatial distribution of the drug, and therefore selection pressures on the tumor will vary depending on the position of a cell relative to nearby blood vessels. In principle, the most aggressive tumor cells tend to develop far from blood vessels. At the same time, these distal positions will be the hardest locations for the drug to reach by diffusion. These observations suggest that a cytotoxic agent may cause spatially variable selection, an effect that can cause enhanced development of aggressive cells in certain circumstances.

To test this prediction, we administered a diffusible cytotoxic drug into the model (see the Supplementary Data for additional discussion of assumptions). Figure 6G shows the tumor size over time for three representative simulations. The solid black line is an untreated tumor. The dashed line is the result of starting the cytotoxic drug regimen at a relative time of t = 0, and the dotted line is the result of starting the regimen at t = 105. For the earlier treated case, the result is a significant decrease in the growth rate of the tumor, while the opposite is true for the later treated case: the tumor has grown faster than our untreated control.

To understand the reason for this difference, key time points from the simulation are shown in Fig. 6A–F. For the case of the early treatment, the therapy has decreased the size of the tumor and essentially set back the time course of growth. Furthermore, a subtle restructuring of the spatial heterogeneity of cell phenotypes has limited the ability of aggressive cells near the center to invade. However, in the later treated case, the therapy has allowed preexisting aggressive cells in the center of the tumor to emerge and become active much sooner than in the other simulations. In the analysis of the untreated growth in “Tumors acquire acid resistance followed by glycolytic capacity”, it was noted that the metabolically benign cells (green) near the tumor boundary acted as a containment mechanism to prevent the more aggressive cells in the interior of the tumor from initiating acid-mediated invasion. In the case of the late administration of chemotherapy, the spatial diffusion of the cytotoxic agent has stripped away this relatively benign protective layer of cells. The aggressive cells situated near the center of the tumor benefit from several advantages under this type of therapy. First, they are further from active blood vessels, and therefore the concentration of drug is decreased because of diffusion distances and consumption by cells closer to the vasculature. Second, limited space and nutrients near the center limits the division rate of these cells, and this also allows them to escape the effects of a drug that targets proliferation. Finally, once the outer layers of cells are stripped away by the drug, these aggressive cells find themselves in a perfect environment to initiate acid-mediated invasion, because they can quickly reach critical mass through rapid proliferation posttherapy and then directly interface with normal tissue cells. This combination causes the tumor phenotype to quickly transform into the acid-mediated invasion type.

The behavior of a cancer as a whole is the relevant clinical metric that informs prognosis and therapy decisions. The actions of a single malignant cell are only relevant if they lead to an aggressive cancer, and this can only happen in the context of the collective behavior of tumor cells, normal tissue, and surrounding microenvironment. Genetic, epigenetic, signaling, phenotypic, temporal, and spatial heterogeneity all play a role in shaping the dynamics of the complex system that drives tumor growth. Therefore, tumor prognosis and response to therapy can only truly be understood if we consider their impact of these heterogeneity mechanisms as a spatially and temporally distributed system.

The results presented here attempt to elucidate how metabolic heterogeneity affects tumor progression and treatment outcome. First, development of localized hypoxia due to vascular instability leads to the selection of acid-resistant tumor cells, which is then followed by acquisition of the glycolytic phenotype. This nonlinear trajectory through phenotype space is most clear in situations of prolonged hypoxia, in which the metabolic response of normal tissue is the Pasteur effect with subsequent acidification of the microenvironment. Extrapolating to in vivo situations, the case of prolonged hypoxia might be representative of ductal carcinomas, where angiogenesis is not possible in the initial stages of tumor growth due to the barrier of the basement membrane. On the other hand, a metastasis that immediately coopts and disrupts the normal vasculature of a tissue may follow the profile of intermittent hypoxia. In both cases, however, the model predicts the same pathway through phenotype space.

Although the present work do not impose any costs on the different phenotypes, we have run simulations where we imposed a simple linear cost to the acquisition of acid resistance (results not shown). In this version of the model, tumors cells that increased their acid resistance had proportionally slower cell cycles. Various settings of the slope of this cost function had little impact on the overall evolutionary dynamic. The trajectory through phenotype space was conserved, and the cost only served to slow down the evolutionary pace. We hypothesize that this is because the competition between cells in hypoxic areas depends much more on the acquisition of acid resistance rather than the ability to proliferate quickly, because the selection pressure is high while the turnover rate is relatively low. A full exploration of evolutionary costs is beyond the scope of the present work.

The simulations of treatments in the model offer some insights as to why these therapies may fail in some cases and be more successful in others. In aggregate, all of the results presented here argue for treatments that are aimed at maintaining the protective layer of less aggressive cells rather than eradicating them. Therapies that buffer the pH_e appear to be preventative rather than regressive. Indeed, there is no clear mechanism by which increasing extracellular pH should lead to cytotoxicity of tumor cells. Rather, the buffer changes the microenvironment so that the acid-meditated invasion phenotype cannot use its acid production against susceptible neighbors. However, the model suggests that there are limits to this preventative action that depend on both buffer dose and the mix of phenotypes in the tumor. A more detailed exploration of the interactions between the phenotypes of the tumor cells and the amount of buffering could perhaps define a quantifiable relationship between dose and growth delay, but this is beyond the scope of this work.

The results suggest that the effectiveness of cytotoxic therapies is strongly influenced by spatial heterogeneity of both the drug diffusion and the tumor structure. This may result in selective intratumoral lethality in which the drug kills less aggressive cells positioned adjacent to blood vessels but not the more malignant populations distant from blood vessels, therefore selecting for more resistant populations and removing a potential barrier to their proliferation by eliminating the perivascular cells. Alternative strategies to optimize the delivery of both dose and schedule to prevent the outgrowth of more aggressive clones posttherapy will be the subject of future work with the model, including drug combinations.

The model examined tumors grown from benign single cells. In the case of metastatic disease, the seeding cell may be of any phenotype present in the primary tumor. Bernards and Weinberg (37) have proposed that metastases arise early in the development of the primary and continue through the life of the tumor. This suggests that early disseminated cells will potentially be less aggressive, but will have more time to establish themselves in the metastatic niche. More aggressive cells, while arising later in the evolution of the primary, may reach invasiveness faster due to their phenotypic advantage. Future work will explore the trade-off between early but phenotypically benign metastases versus late, aggressive metastases.

All models incorporate limiting assumptions. We have developed the model in two dimensions, representing a slice through a tissue. Although three-dimensional modeling has been useful for exploring the specifics of tumor morphology and blood vessel remodeling, our interest here is in the evolution of tumor phenotypic heterogeneity as a function of environmental selection. We expect that there would be little advantage to implementing the model in 3D with respect to the insights provided on acquisition of heterogeneity. However, the principles of the model can be directly extended into three dimensions in the future.

The present work uses a lattice-based approach without mechanical deformation. Tumor cells are limited to grow only when there is space available in the lattice. Cells are not displaced volumetrically, as may be the case in a real tumor. This assumption imposes a strict limit on the growth of a tumor, whereas a deformable tissue may permit tumor growth through normal tissue even if available adjacent space is absent, thereby escaping the encapsulation and corresponding homeostatic growth phase. To investigate these effects, however, would require a significantly more complex model with mechanical forces and off-lattice modeling, and therefore is beyond the scope of the present work.

We have only considered acid-mediated invasion, and other forms of invasion such as MMP secretion and EMT transition are also important. The model can be extended to examine the effects of phenotypic changes on these additional axes in the future. Furthermore, cell migration without proliferation was not considered, and the effect of migration on spatial heterogeneity is likely to play a role in how the tumor develops. There have been previous models in relation to metabolism that incorporate these effects (22, 38, 39). However, the effects of acid-mediated invasion rely on critical mass, and dispersion of individual cells through the normal tissue due to migratory elements may not necessarily lead to faster development of the invasive phenotype. Clearly, one logical next step is to examine the influence of cell migration on the properties that lead to acid-mediated invasion. Escaping cells might be more likely to penetrate the protective shell of relatively benign tumor cells seen in the present work, but migration may also mean more cellular dispersion, which would spatially dilute the acid-production rate of a group of aggressive cells, limiting the effectiveness of acid-mediated invasion. Indeed, migration may be less important when trying to establish a niche, but more important in reaching that niche in the first place, such as with micrometastases.

In this work, we chose a point-source vasculature model. This is representative of a field of parallel blood vessels passing perpendicular to the two-dimensional tissue slice and is a simplification that has been considered before (40). There is a significant amount of literature devoted to modeling angiogenesis and nutrient delivery (18, 19, 41). However, including these complex mechanisms are beyond the scope of this work. We opted for a simply defined vascular system that would deliver temporally and spatially heterogeneous nutrient and waste gradients with a minimal number of parameters. A dynamic point-source vasculature serves as a good surrogate for a more complicated model that would accomplish the same goal. Indeed, much of the behavior seen in this model is sensitive to localized nutrient gradients, and it is unlikely that a complex vasculature would produce gradient structures on a local level that are significantly different from a point-source system.

For the present work, we confined the phenotypic variation to small unbiased changes occurring only at the time of division. Other methods of phenotypic variation can easily be implemented in the model in the future. For example, major mutations to key metabolic genes could be modeled by adding rare but large jumps in these parameters. Cells could also be allowed to alter their phenotype independently of proliferation.

In summary, we have developed a hybrid multiscale model of tumor growth within a normal, homeostatic tissue, and demonstrated the mechanisms by which phenotypic, temporal, and spatial heterogeneity affect growth of a tumor and the outcomes of treatments. The model predictions are consistent with clinical and experimental observations and provide insight into the mechanisms that result in treatment failure. This work highlights the importance of the phenotypic heterogeneity of tumor cells and how this heterogeneity varies in time and space. A key prediction of our model is that early stages of tumor development (either primary or metastatic) maintain a phenotypically spatially structured population, where less aggressive clones spatially suppress more aggressive counterparts. Importantly, we demonstrate that standard treatment modalities may selectively destroy this structured population and facilitate subsequent progression. Controlling tumor progression by maintaining rather than destroying this suppressive tumor layer appears to be more effective that conventional high-dose density therapy that aims to kill the maximum possible number of tumor cells.

No potential conflicts of interest were disclosed.

Conception and design: M. Robertson-Tessi, R.J. Gillies, R.A. Gatenby, A.R.A. Anderson

Development of methodology: M. Robertson-Tessi, A.R.A. Anderson

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M. Robertson-Tessi, R.J. Gillies

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M. Robertson-Tessi, A.R.A. Anderson

Writing, review, and/or revision of the manuscript: M. Robertson-Tessi, R.J. Gillies, R.A. Gatenby, A.R.A. Anderson

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): R.J. Gillies, A.R.A. Anderson

Study supervision: A.R.A. Anderson

The authors thank Veronica Estrella for providing the mouse data.

This work was sponsored in part by the Moffitt Cancer Center PSOC, NIH/NCI U54CA143970. This work was supported in part by the Cancer Center Support Grant (P30-CA076292) from the National Cancer Institute.

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.