Abstract
A quantitative model of carcinogenesis based on methods from population biology and game theory demonstrates normal cells in vivo occupy a ridgeshaped maximum in a welldefined tissue fitness landscape, a novel configuration that allows cooperative coexistence of multiple cellular populations. This state, although necessary for development of functioning multicellular organisms, is subject to invasion by fitter, mutant phenotypes permitting somatic evolution of cancer. The model demonstrates carcinogenesis is an emergent phenomenon requiring a sequence of evolutionary steps as cellular proliferation follows successful adaptation to varying environmental constraints. In the initial development of preneoplastic lesions, cellular proliferation is controlled exclusively by interactions with other cells, the extracellular matrix, and soluble or insoluble growth factors so that gain of function mutations in oncogenes, loss of function mutations in tumor suppressor genes, and disruption of normal senescence pathways will permit clonal expansion. This provides explicit selection mechanisms for the mutations depicted in the classical FearonVogelstein model of colorectal carcinogenesis. The model demonstrates neoplastic cellular proliferation can also be promoted by alterations in the somatic landscape that reduce inhibitory signals produced by the normal cells and extracellular matrix. This is consistent with experimental evidence for a strong microenvironmental influence in tumorigenesis independent of genomic changes in the neoplastic populations. However, we find that these changes alone produce only selflimited neoplastic growth because cellular crowding alters system dynamics so that proliferation is dependent on substrate availability. Consequent cellular competition for critical nutrients results in Darwinian selection pressures favoring phenotypes that increase substrate delivery (e.g., angiogenesis) or uptake (e.g., amplify membrane transporters). These previously unknown substrate dynamics in the later stages of carcinogenesis provide a mechanism for evolution of cellular properties typically found in invasive cancers including the angiogenic and glycolytic phenotypes.
INTRODUCTION
Quantitative methods are essential for defining system dynamics and critical parameters in complex, fundamentally nonlinear processes such as carcinogenesis (1 , 2) . Mathematical models provide a framework of understanding for interpretation, organization, and clinical application of extant data as well as integration of rapidly accumulating new information. Although many of the genetic events in carcinogenesis are known (3 , 4) , their precise interactions with environmental factors that control clonal expansion and malignant progression are unclear in the absence of a generalized model. Similarly, without understanding the mechanisms by which cellular properties confer growth advantages on mutant clones, biologically and clinically important genetic changes cannot be distinguished from the large background of evolutionarily neutral mutations.
Fearon and Vogelstein (5) have proposed a phenomenological model of colorectal carcinogenesis correlating specific genetic events with evolving tissue morphology. This conceptual approach describes a system that proceeds linearly from normal mucosa to a small polyp to a large polyp to an invasive cancer, with each step driven by welldefined alterations in the genome. This word model, although limited by the linear logic typical of that approach, may form the schematic framework for mechanistic, nonlinear quantitative models that more fully describe the complex biological system. Because the outcomes of such nonlinear interactions cannot be determined by intuitive reasoning alone, they must be computed from analysis of appropriate mathematical models.
A variety of mathematical approaches to carcinogenesis have been proposed during the past several decades, as nicely summarized by Bogen (6) . Much of the more recent work has focused on controversial role of the mutator phenotype and the dynamics of accumulating gene mutations in the transition from normal cells to invasive cancer (7, 8, 9, 10) . Other quantitative research has used both population ecology approaches (11) and game theory (12 , 13) . These modeling studies have generally focused on intracellular processes such as genomic mutations or the interactions of the mutant cells with other cellular populations. However, little work has focused on the interaction of the mutant phenotypes with environmental selection parameters. That is, favorable mutations are typically defined simply as those that confer “a selective growth advantage.” Extant models have not generally attempted to identify the explicit environmental selection forces that control cellular proliferation or the dynamics of the interactions of these selection factors with mutant cellular properties. In addition, very little attention has been directed to experimental evidence that alterations in the microenvironment, in the absence of cellular mutations, strongly affect the evolution of the malignant phenotype (14 , 15) . Furthermore, these theoretical models generally ignore the effects tumor cells have on their own microenvironment and how the presence of mutant cells may alter environmental selection parameters for subsequent phenotypes generated by random accumulating genomic mutations. It is likely that these short and longterm interactions of evolving mutant cellular populations with a changing microenvironment must be defined to allow full understanding of the dynamics of carcinogenesis. Finally, virtually no theoretical work in carcinogenesis has integrated dynamics that produce two critical properties of invasion: angiogenesis and the near universal adoption of the glycolytic phenotype.
Here we propose a comprehensive, mechanistic, quantitative model of carcinogenesis based on methods developed in population biology and game theory that mathematically frames the FearonVogelstein model as a sequence of competing populations subject to random mutations resulting in evolution of optimal proliferative strategies in a changing adaptive landscape. The purpose of this work is to clarify and extend the FearonVogelstein conceptual model by (a) identifying environmental selection parameters that control cell proliferation and clonal expansion, (b) defining the interactive dynamics by which these environmental selection parameter produce the typical characteristics of the malignant phenotype, and (c) understanding the environmental effects induced by successive generations of progressively more transformed cellular populations.
Through this model, we wish to establish the beginning of a quantitative, theoretical framework for carcinogenesis that is biologically realistic and consistent with extant data. However, it is important to emphasize that development of mathematical models is a dynamic, interactive process, and many aspects of the current work may require revisions as new information accumulates, and as the model predictions are tested against experimental data. Nevertheless, as in the physical sciences, these mathematical models are essential tools for understanding the complex, nonlinear systems in carcinogenesis and provide the necessary next step beyond current theories based on verbal reasoning and linear intuition. Our main objective is to illustrate the potential usefulness of this approach in the study of carcinogenesis.
MATERIALS AND METHODS
General Cellular Population and Substrate Dynamics in Multicellular Organisms.
Population and substrate dynamics equations (16) assume cellular populations in a multicellular organism are subject to two broad control mechanisms: tissue developmental factors and substrate availability. The former consists of the complex cellular and environmental interactions that form functioning tissue from multiple individual cells. The latter simply requires that a cell proliferates only if substrate uptake exceeds basal metabolic demands.
Here we use a modified form of the general state equations from a previous study (16) . Each volume of tissue contains p different populations with N_{p} individuals. We assume two broad mechanisms controlling cellular proliferation: (a) tissue factors communicated through cellular interactions with other cells, the extracellular matrix, and soluble growth factors; and (b) availability of substrate sufficient to generate the energy and macromolecule necessary for formation of a new cell. We include only the substrate dynamics of glucose because it has been extensively studied, and critical parameter values can be readily obtained from this literature. Assume, for simplicity, that all neoplastic cells retain the proliferative capacity. We use a MichaelisMenten substrate uptake and a fixed basal consumption rate m. We introduce a constant B that converts substrate quantity to proliferation (i.e., the quantity of glucose required to successfully proceed through mitosis producing one more individual), a maximum proliferation term b, and a death (or apoptosis) rate term d. The resultant normal cellular population dynamics (designated with the subscript n) is given by Accumulating mutations produce a variety of genotypes and thus multiple potential tumor populations (designated with the subscript t) modeled by where N_{i} is the number of individual in each population (e.g., N_{1} is normal tissue and N_{2}, N_{3}, ⋮ is tumor tissue), p is the total number of different distinct cellular populations, and the dot is derivative with respect to time.
The substrate is assumed to be delivered by the local vascular network and consumed by the cells with the following dynamics: where r is glucose delivery rate that is dependent on vascular density and flow rate according to where r_{e} is a delivery factor (>1). This delivery factor must be greater than one; otherwise, glucose could not be delivered at a rate to match the demand requirements. However, there is a maximum rate at which the body can deliver glucose, so we include the following:
Normal tissue proliferative controls are encompassed in the N_{max} term, which is a lumped phenomenological term that sums the positive and negative intra and extracellular growth factors within the tissue environment.
The potential constraint from substrate availability is encompassed in the second term so that cellular proliferation occurs only if substrate uptake quantitatively exceeds the amount necessary for maintenance of normal cellular functions. We anticipate substrate dynamics will ordinarily be regulated through physiologic control mechanisms so that the amount delivered, r, equals the amount metabolized, so that the interstitial concentration, R, is constant.
The process of carcinogenesis represents progressive drift of the cellular genotype and phenotype from the population of origin due to accumulating genomic mutations that produce multiple novel genotypes. Each new genotype represents an evolutionary “experiment” and will produce a new population only if the resulting phenotype yields sufficient growth advantage to allow clonal expansion. That is, a mutated cell N_{2} arising from a preexisting population N_{1} is clinically relevant only if Ṅ_{2} > 0 under the current stable conditions in which Ṅ_{1} = 0.
Cellular Evolution.
To convert the above model to an evolutionary model, we identify two adaptive parameters that interact with the dominant proliferative constraints in the system. The first, u_{i1}, is the number of substrate transporter on the cell surface. It is normalized so that u_{i1} = 0.5 corresponds to the average number of transporters/cm^{2} in a healthy cell. It directly determines the substrate uptake terms E_{n} and E_{t} as given in the Appendix. The second parameter, u_{i2}, is related to cellular response to normal growth constraints. It is normalized so that u_{i2} = 0.5 corresponds to the sensitivity of normal cells in responding to environmental signals for proliferation. It directly determines the distribution and density of each cell population, terms N_{nmax} and N_{tmax} as given in the Appendix. All other parameters are assumed constant. This approach allows us to explicitly examine, both separately and together, mutations that affect the two primary cellular growth constraints as outlined above. Details of the evolutionary model and parameter values used are included in the Appendix.
We first examine the evolutionary stability of normal cells through the use of the fitness generating function (Gfunction). A brief discussion of the Gfunction method is given in the Appendix. Basically, the Gfunction is used to determine the fitness of all extant strategies or phenotypes as well as the fitness of any phenotype the current populations could adopt through an evolutionary process. In this way we can determine the optimal strategies or phenotypes available within the specific fitness landscape and thus predict the outcome and dynamics of the somatic evolution of the cellular population, assuming the environment remains stable (possible confounding effects of environmental instability are discussed below). The Gfunction is a function of all population numbers, strategies, and a virtual strategy vector (17) (the “virtual” strategy allows determination of the fitness of strategies both present and not present in the extant populations but available through evolution). A strategy vector used by a cell type i (e.g. [u_{i1} u_{i2}] is evolutionarily stable if the Gfunction is a proper global maximum with respect to the virtual strategy vector. If one plots the Gfunction versus all possible values for the virtual strategy, then this plot provides a graphic way of determining evolutionary stability, because the global requirement is easily verified. Such a plot is called an adaptive landscape.
All of the figures are generated using model details described in the Appendix using the parameter values listed there.
Normal Tissue Dynamics.
We propose that normal cells in the in vivo environment of a multicellular organism interact with definable adaptive landscapes with characteristics dependent on time and spatial location within the organism. To briefly illustrate the potential dynamics, we examine the differentiation pathway of normal pluripotent cells located at a specific tissue site during, for example, embryogenesis or tissue repair. This process may be defined as a warping of the current local tissue adaptive landscape as shown in Fig. 1 ⇓ , which shows four different snapshots of the adaptive landscape. These snapshots illustrate how the landscape changes as the cells grow from an initial value of 200 at time t = 0 to a final equilibrium value of 1000. Note that the fitness of any cell using a specific strategy is determined by the point on the landscape corresponding to the two components of the strategy vector. A possible initial state for a pluripotent normal cell is indicated by the dot. This strategy corresponds to evolutionary equilibrium for normal cells and remains a maximum point for every snapshot as the landscape changes shape underneath it. Because of the interactive, nonlinear dynamics, the characteristics of the growing cell population alter the adaptive landscape. In the case of normal tissue, the differentiated cell populations ultimately find themselves at a point on a ridgeshaped maximum.
Interestingly, this novel configuration allows stable coexistence of multiple cellular “species” defined by the same Gfunction but with different “strategies” along the edge of the ridge. These results provide a simple, general control mechanism for simultaneous direction of pluripotent stem cells into a population with a specific differentiated phenotype and maintenance of a society of stable, noncompeting cellular populations in functioning tissue, a necessary condition for formation of multicellular organisms. This general mechanism is supported by experimental data demonstrating that differentiation in embryogenesis and tissue repair is controlled by signals from the microenvironment (18, 19, 20, 21) .
Note the orientation of the ridge in nonequilibrium tissue (at t = 0) along the v_{2} axis, which defines the cellular interaction with normal environmental proliferation constraints. This indicates that coexistent, noncompeting populations must possess different strategies defined by these normal tissue controls (i.e., proliferation constraints by different combinations of oncogenes and tumor suppressor genes). That is, organization of function in developing or remodeling normal tissue (i.e., tissue not at a stable steady state) requires neighboring cells to possess different sets of growth control parameters (i.e., different “strategies”). The requirement for distinctly different growth control strategies in coexisting cellular populations provides a potential explanation for the large number of oncogenes and tumor suppressor genes found in the human genome and the diversity of their mutations found among cancers of different cell types.
It appears that invasive cancer is the “price” of the fitness landscape that allows formation of multicellular organisms. Because normal cells do not achieve a proper maximum at equilibrium, there is an opportunity for tissue invasion by mutant cellular phenotypes. As long as the mutants are derived by the same Gfunction, they do not change the shape of the adaptive landscape, and they will simply coexist at low numbers. However, if a mutant should have a different Gfunction, the situation is quite different (22) . In fact, because the system dynamics are highly nonlinear, the introduction of a mutant population initially at a higher fitness, as discussed below, can deform the adaptive landscape such that normal cells now occupy a local minimum. This situation is illustrated in Fig. 2 ⇓ , where it is shown how the presence of a tumor cell deforms the normal cell adaptive landscape so that, at equilibrium, the normal cells are at a local minimum on the adaptive landscape. The tumor cells are introduced in small numbers (x_{2} = 10 at t = 0) and grow to x_{2} = 218 at equilibrium. Because neither the normal or tumor cells are allowed to evolve, the tumor remains at a this equilibrium value and does not affect the number of normal cells. However, because the normal cells are now at a minimum, the stability of normal tissue requires proliferative controls that are both redundant (i.e., several oncogenes and tumor suppressor genes) and robust (i.e., maintenance of a sufficiently low basal rate of cellular mutations to prevent significant cellular evolution during the reproductive life of the organism).
We have seen that normal tissue may be invaded if the mutant population can change the adaptive landscape in a detrimental way (to normal cells). We will now show that if the mutant cells can evolve (toward a peak on the mutant adaptive landscape), the combined effect allows for the mutant type to replace the normal cells. We propose the transition of cellular populations from their normal configuration on the adaptive landscape to one that promotes emergence of a mutant phenotype constitutes the somatic evolution of an invasive cancer.
Evolution of Invasive Cancer.
The model demonstrates that the evolutionary steps of invasive cancer take place at bifurcations in the system dynamics that result in a change of the equilibrium state. We now examine this process using the evolution of colorectal cancer as a specific example. Because we find that each normal cell population must be controlled by different sets of oncogenes and tumor suppressor genes, favored mutations in specific genetic loci will vary from one cell type to another. In colorectal carcinogenesis these targets appear to include, among others, the APC and KRAS genes (3 , 23 , 24) . The former is a multifunctional gene that acts as a tumor suppressor through downregulation of βcatenin, which is a transcription activator (25) found at intercellular junctions (from which it transits to the nucleus), suggesting a role in transmitting information regarding contact with neighboring cells and other aspects of the microenvironment. The majority of KRAS mutations are gainoffunction missense leading to increasing growth stimulatory signals, although they also affect cell adhesion, cell cycle, and cell metabolism (26) .
Thus, KRAS or APC mutations will diminish environmental control of proliferation in the tumor cells through increased translation of growth promoter signals and decreased inhibition from cellcell contact (27) . This is incorporated into the model by increasing the mean tissue carrying capacity of tumor cells.
With a normal background mutation rate of about 2 × 10^{−7} mutations/gene/cell division (28) , ς_{n} = ς_{t} = 0.0001. Mutations in u will be very close to the mean value, resulting in a very small rate of evolution. This is accounted for in the figures by setting SD of mutations among the normal and tumor cells to a small value. In generating Fig. 2 ⇓ , all parameter values for the normal cells are the same as those used in Fig. 1 ⇓ , and all parameter values for the tumor cells are the same as the normal cells, except we increase the mean tissue carrying capacity of tumor cells. This change in carrying capacity may be due to two general reasons. First, the neoplastic cell may experience a mutation in loci such as APC or KRAS that interferes with its ability to detect, process, or respond to normal tissue growthinhibitory signals. Second, perturbations in the surrounding tissue could result in a decrease in the production of these signals by adjacent normal cells or extracellular matrix. The former encompasses the cellcentric approach epitomized by the FearonVogelstein model. The latter is consistent with work by Bissell, BarcellosHoff, and others demonstrating a strong microenvironmental influence on carcinogenesis (14 , 15) .
Because of this change in tissue carrying capacity, the mutants’s fitness can not be determined from the normal cells’ Gfunction, and a second Gfunction is introduced for the mutant. Under these conditions, numerical simulations demonstrate only selflimited growth of the mutant population, as illustrated in Fig. 3 ⇓ . Note that this coexistence result is possible only because r is limited, and thus, cellular proliferation is constrained by substrate availability. In other words, we find that even though normal cellular proliferation constraints are lost due to intracellular mutations or environmental changes, proliferation of the mutant populations is constrained by diminished substrate availability due to cellular crowding and vascular compromise.
In the next phase, the evolutionary constraints on the tumor cells are removed by assuming an increased mutation rate due to microsatellite instability, chromosomal instability, or mutagenic environmental factors such as chronic inflammation or infection. We now examine the mechanisms by which the tumor population can evolve to an ESS, 2 driving the normal cells to extinction, resulting in an invasive cancer.
There are two possible evolutionary pathways in the model: further reduce social constraints through increases in mean tissue carrying capacity or loss of normal senescence pathways or increasing substrate acquisition and delivery. The former might result from mutations in the DCC, telomerase, and p53 genes that are typically found in colon cancers and advanced preneoplastic lesions (5) . We simulate this stage by increasing the rate of evolution (28) in tumor cells but not normal cells, increasing mean substrate uptake for tumor cells, and increasing the value of the tumor adaptive parameters for maximum mean tissue carrying capacity and mean substrate uptake. This reflects the fact that the greater potential offered by increasing N_{tmax} and E_{tmax} is obtainable through evolution. All other constants held the same. Fig. 4 ⇓ illustrates the results. The crowding of the cells is now more apparent during the early “plateau” stage. Thus, significant tumor progression from a small adenoma to a large adenoma has occurred with the addition of cellular dynamics permissive for evolution. Fig. 5 ⇓ shows how the adaptive parameters evolve through time.
An important insight from these simulations is obtained by noting the isolated effects of accumulating mutations in tumor suppressor genes and oncogenes in the later stages of carcinogenesis. We find these genetic changes contribute to tumorigenesis in colorectal cancer by increasing N_{tmax} and therefore increasing the size of the neoplastic population (i.e., small polyp or large polyp). However, as noted above, changes only in N_{tmax} produce a quasistable steady state that is not an ESS. We find the critical final steps allowing mutant populations to form an ESS and, therefore, invasive cancer are mediated by an environment dominated by substrate competition. This is illustrated below.
From Fig. 4 ⇓ it is evident that the transformed cells in the preneoplastic lesion (i.e., the first plateau) consume more resources than can be delivered from the altered environment reducing the local concentrations of substrate (glucose). Thus, there is a transition such that cellular proliferation is no longer constrained by tissue social factors as encompassed in the N_{tmax} term, and system dynamics are dominated by the second term (i.e., the availability of adequate substrate, equivalent to the bifurcation of solutions E and N, see Appendix).
The mutant populations now need to evolve new strategies to maximize fitness within an adaptive landscape dominated by substrate limitation. Available strategies include increasing substrate delivery through angiogenesis, increased efficiency of substrate uptake, or both. As is evident from Fig. 5 ⇓ , these new selection parameters drive cellular evolution to increasingly efficient acquisition of glucose in an environment of low glucose concentration by increasing the value of E_{t}. An equilibrium is finally achieved when a mutant population is able to maximize its substrate uptake and in the process further decrease substrate concentrations such that the other populations are driven to extinction. Thus, only after evolving strategies that allow for increased substrate uptake will the system come to equilibrium with complete destruction of all normal cells, i.e., an invasive cancer. Because its strategy is at a maximum on its adaptive landscape, the cancer population cannot be dislodged through the introduction of other mutant strategies (derived from same Gfunction). This prediction of a decline in somatic evolution and, therefore, decreased cellular heterogeneity in an established invasive cancer has been observed (28) and is further discussed below.
DISCUSSION
A variety of quantitative methods, including matrix evolutionary models, have been applied to carcinogenesis (10 , 12 , 13 , 29) . In general, this work has focused on the controversial role of the mutator phenotype. Our approach attempts to develop a reasonably comprehensive model of the entire process and focuses primarily on the specific cellular phenotypic changes and environmental selection parameters that control the dynamics during somatic evolution of cancer. This work represents an initial effort, and the model will undoubtedly undergo revision and expansion to accommodate the full range of cellular and environmental parameters (e.g., immune response) that were not included. Nevertheless, we find that these elementary models develop features that mimic many observed properties of carcinogenesis. We do not consider this coincidental for two reasons: first we are working from wellknown mechanistic models, and second, we are using actual data for estimating parameter values. Indeed, the model is highly sensitive to these parameters, and we find that our simulation results are strongly dependent on utilization of physiologically realistic values. In addition, these simulations allow several explicit predictions regarding carcinogenesis that can be compared with extant data.
First, a small nonmalignant neoplasm will result from an isolated mutation in an oncogene or tumor suppressor gene. Examples of these tumors likely include skin nevi or small adenomas in the colon. Because most of these lesions will be due to the background mutation rate, they will increase monotonically with the age of the host. Each of these benign lesions is predicted to consist of a monoclonal population possessing a single stable mutation in an oncogene or tumor suppressor gene. This prediction is consistent with recent observations from microdissection studies (30) . A germline mutation will produce a large number of small lesions in those tissues in the host which are controlled primarily by the mutated gene [e.g., familial adenomatous polyposis (26)] .
Second, small lesions resulting from isolated mutations will progress to invasive cancer (i.e., to an ESS) only when allowed to evolve. This is consistent with the hypothesis that the mutator phenotype is frequently an early and (in the absence of environmental mutagenic factors) necessary development in carcinogenesis (28 , 31) . The critical role of increased mutation rate and more rapid evolution in mammalian cells is supported by in vitro studies in which exposure to a mutagen decreased the time necessary for emergence of neoplastic transformants from cultured normal mammalian cells (32) . We note, however, that environmental mutagenesis may initially substitute for the mutator phenotype so that a small adenoma may progress if subjected to a mutagenic environment such as chronic infection or inflammation (33 , 34) . Furthermore, we note that induction of tumor growth can also be accomplished by alterations in the local microenvironmental tissue growth controls. This is consistent with experimental observations that the invasive phenotype can be induced or suppressed by changes in the environment both in vitro and in vivo in the absence of genetic changes in the target population (14 , 15 , 35 , 36) .
Third, a second stage of carcinogenesis, dominated by competition for substrate, is predicted. The resulting dynamics select tumor populations with increased glucose uptake, consistent with decades of experimental and clinical observations (37 , 38) . The predicted transition to the glucose avid phenotype in the evolution from advanced premalignant lesions to invasive cancer is consistent with observations that increased expression of glucose transporters is a late event during malignant progression in both colon and esophageal cancers (39 , 40) . In addition, 18Ffluorodeoxyglucose positron emission tomography (FDGPET) imaging studies have demonstrated that the transition from large polyp to invasive cancer coincides with adoption of increased glucose uptake (41) . Similar increases in GLUT expression have also been noted in the transition from premalignant lesions to cervical cancer (42) and in cholangiocarcinoma. Interestingly, a decrease in glycolysis and glucose uptake has been shown to correlate extremely well with response to therapy in both breast and prostate cancer (43 , 44) .
Finally, the model makes the novel and somewhat counterintuitive prediction that premalignant lesions should demonstrate more phenotypic heterogeneity than the subsequent invasive cancer. This is because each new genotype in the premalignant state is not an ESS. Therefore, its growth is selflimited, and it will typically achieve a quasistable steady state with other populations. The net result will be a mosaic of different genotypes and phenotypes within the evolving premalignant lesion. However, a genotype allowed to evolve to an ESS will invade and destroy adjacent nonESS populations. Furthermore, a cellular strategy at an ESS it will have no further pressure to evolve at least within the current microenvironment. Thus, the transition from a premalignant state to an invasive cancer should be accompanied by a transition from heterogeneous cellular population to one that is comparatively homogeneous. This is consistent with analysis of prostate invasive cancers and premalignant lesions (prostatic intraepithelial neoplasia) that demonstrated significantly greater phenotypic heterogeneity in prostatic intraepithelial neoplasia than invasive cancer (45 , 46) .
We present a quantitative model of carcinogenesis that can provide a framework of understanding for organization of extant data, integration and interpretation of new information, and guidance for future experiments. The model provides a clear mechanisms for the sequence of genetic mutations and corresponding tumor growth patterns described in the phenomenological “Vogelgram.” The additional or alternative role of alteration in the normal microenvironment as described by Bissell and others can be readily incorporated into this evolutionary sequence. Finally, the possible role of cellular competition for substrate during the later stages of carcinogenesis is demonstrated, providing a mechanism for selection of the angiogenic and glycolytic phenotypes, characteristic properties of invasive cancers that are not explained by any current model of tumorigenesis.
It is important to emphasize that a major motivation for development of this comprehensive carcinogenesis model is stimulation of further theoretical and experimental studies aimed at forming progressively more complete and accurate integrative models of this complex process. Thus, this work should be regarded as an initial effort rather than a final product and will undoubtedly “evolve” as inconsistencies with empirical studies require revisions. This interactive approach of quantitative models and experimental results will likely be necessary for a complete understanding of the complex, fundamentally nonlinear cellular and microenvironmental dynamics that result in carcinogenesis.
Introducing Adaptive Parameters into the Model
The model represented by Eqs. 1 2 3 4 can not evolve without specifying how it depends on the adaptive parameters. We have identified u_{i1} and u_{i2} as adaptive parameters. Because adaptive parameters will generally have a normal distribution of values within a population, we use the following distribution functions where

N_{nmean} = Mean tissue carrying capacity of normal cells

N_{tmean} = Mean tissue carrying capacity of tumor cells

E_{nmean} = Mean substrate uptake for normal cells

E_{tmean} = Mean substrate uptake for tumor cells

u_{Nnormal} = Value of u_{i2} for largest N_{nmax} i = 1

u_{Ntumor} = Value of u_{i2} for largest N_{tmax} i = 2, …, n

u_{Enormal} = Value of u_{i1} for maximum E_{n} i = 1

u_{Etumor} = Value of u_{i1} for maximum E_{t} i = 2, …, n

ς_{N} = Variance in N distributions

ς_{E} = Variance in E distributions
The maximum values for substrate uptake for normal cells are at a maximum when u_{11} = u_{Enormal}, and the substrate uptake for tumor cells is at a maximum when u_{i1} = u_{Etumor}. Likewise, the tissue carrying capacity for normal cells is at a maximum when u_{12} = u_{Nnormal}, and the tissue carrying capacity for tumor cells is at a maximum when u_{i2} = u_{Ntumor}.
Parameter Values
To make the model as realistic as possible, we assume that the critical substrate is glucose and use parameter values available in the literature (38) . The following values are converted to consistent units for use with the model. Other values not specified here but required for simulation will be noted and supplied as required. Interestingly, we find the models highly sensitive to use of realistic parameter values. Substitution of values outside of the normal physiologic range typically produced physically unrealistic results.

B_{n} = B_{t} = 0.06 × 10^{6} cells/μmol glucose

(b_{n} − d_{n}) = (b_{t} − d_{t}) = 0.8333

N_{nmean} = 10^{9} cells/cm^{3}

N_{tmean} > 10^{9} cells/cm^{3} (range: 10^{9} − 2 × 10^{9})

E_{nmean} = 14.668 μmol glucose/10^{6} cells/day

E_{tmean} > 14.668 μmol glucose/10^{6} cells/day (range: 15–120)

m_{n} = 2.544 μmol glucose/10^{6} cells/day

m_{t} > 2.544 μmol glucose/10^{6} cells/day (range: 7.92–28.8)

R_{n0} = 3330 μmol glucose/cm^{3}

R_{t0} < 3330 μmol glucose/cm^{3} (range: 3330–150)
Fitness Generating Functions
The Gfunction method (17 , 47 , 48) is used to examine the evolutionary stability characteristics of this model. We briefly review this method for completeness.
The full characterization of a cell will generally require a vectorvalued strategy 3 whose dimension equals the number of independent, heritable traits available to the cell. In our model, we have restricted the number of traits to two. Each pair of traits represents a strategy, so that in a population of n strategies, there will be a total of 2n heritable traits or parameters represented by the vector where the strategies are partitioned using  to clarify which components belong to each strategy type. The normal cell components are contained in the first partition, and the tumor cell components are in the remaining partitions. Different strategies in the population define differences among individual cells. Thus, for each strategy u_{i}, there is an associated number of cells N_{i} that possess it. The total number of cells in the population, N, is simply the sum of all the different types of cells as defined by their strategies. We may also think of as a row vector where N_{i} make up the n components.
The individual fitness H_{i} for the cell type N_{i} using strategy u_{i}, will be a function of u, N, and R, which we write as H_{i}(u, N, R). Once the fitness functions have been defined, as per above, the population dynamics for each of the cell types are given by equations of the form We assume that given a constant strategy u*, there exists at least one nonzero equilibrium solution N* (that is, not every component of N* is zero) and R*. There are two ways we can find an equilibrium solution for N*. One way is to solve for N* from the system of equations Note that, in general, a solution to this system of equations will require a numerical procedure. If more than one equilibrium solution exists, then the particular solution obtained will depend on the initial guess made for the solution. A second method for determining an equilibrium solution (when such a solution is asymptotically stable), is to choose u* and an initial condition N(0) and simply let the solution to Eqs. 1 and 2 determine the equilibrium point by integrating the differential equations. Numerical solutions to these systems of simultaneous equations maybe obtained using programs such as Matlab.
At equilibrium, one or more components of N* must be positive and nonzero for there to be a viable solution. The corresponding strategies are those that can coexist in the population of cells. But, does this solution predict the outcome of evolution by natural selection, and is this an efficient way to model evolution within a population? The answer to both questions is ‘No’. The equilibrium solution to the above system of equations only considers the outcome for those strategies already resident in the population and does not consider, nor can it consider, the potentially infinite number of feasible strategies that may occur in the future via selection and/or mutation. Furthermore, advantage is not taken of the fact that each individual cell shares an evolutionary history and context with others in the population via common ancestry. Whereas individuals within a population of cells are distinctive in the current strategy they are using, they are distinctive in other ways as well. We may group them not only by strategy, but by evolutionary potential. That is to say, individuals within a population, while using different strategies, may be evolutionarily identical. These individuals draw their strategies from the same set of evolutionarily feasible strategies, 𝒰 (in general, this set would depend on biological constraints; however, in this paper, the strategies are assumed to be unconstrained).
The fact that there will be one or more groups of evolutionarily identical individuals in the population creates a type of symmetry such that evolutionarily identical individuals are completely interchangeable. Only their strategies are relevant for determining an individual’s expected fitness. When a population of individuals are evolutionarily identical, just a single fitness generating function is needed to describe fitness within the population. Given the set of all evolutionarily feasible strategies, the fitness generating function (Gfunction) (49 , 50) is defined as: where v is a virtual strategy giving the Gfunction the property of being identically equal to the individual fitness functions, H_{i}, when the virtual strategy is set equal to the strategy u_{i}. In other words, once the virtual strategy is set equal to any population strategy u_{i}, the Gfunction becomes identical to the fitness function for those individuals using the strategy u_{i}. In this context, the Gfunction describes the evolutionary potential for all evolutionarily identical individuals.
Since the parameters for normal cells and tumor cells are generally different, we need two Gfunctions for modeling carcinogenesis, one for the normal cells, G_{n}, and one for the tumor cells, G_{t}. Using the Gfunction definition, we obtain where v = [v_{1} v_{2}]. The first component of v determines E according to and the second component of v determines N according to
Conditions Promoting Carcinogenesis
The population vector N is composed of both normal and tumor cells. Assuming there is only one type of normal cells, N_{1}, and any number of tumor cell types, N_{2}, N_{3}, …, N_{n}, the population dynamics for these cells, in terms of the Gfunctions G_{n} and G_{t}, become The resource dynamic equations remain unchanged with and in the calculation of r. Finally, because of the distribution of strategies available to each cell type, there exists strategy dynamic equations (17 , 47) that, to first order, are given by where ς_{n} and ς_{t} are related to the standard deviation of mutations among the normal and tumor cells.
We may now examine the conditions in normal tissue that prevent or support carcinogenesis by solving the cell dynamic equations. In particular, we seek equilibrium solutions to the above equations that satisfy an ESS maximum principle (17) . This principle requires that the equilibrium solutions be ecologically stable (return to equilibrium for a fixed u) and evolutionarily stable (return to equilibrium under changes in u). Equilibrium solutions that do not satisfy the ESS maximum principle are, in general, subject to invasion. The ESS maximum principle may be visualized geometrically by plotting G (v, u*, N*, R*) as a function of v. That is, plotting G as a function of v at equilibrium conditions. For example, u*_{1} will satisfy the ESS maximum principle only if the plot of G(v, u*, N*, R*) takes on a zero maximum when v = u_{1}. If the plot shows that G(v, u*, N*, R*) takes on a zero minimum or a saddle when v = u_{1}, the resultant solution is evolutionarily unstable.
The nature of the equilibrium solutions depend on n. Consider first the case of n = 1. In this case, assuming N*_{1} ≠ 0, equilibrium requires where r must satisfy Eqs. 11 and 12 . Generally, all four equations must be solved simultaneously for the equilibrium values of N*_{1}, R*, and u*_{1}. However, we may think of using Eqs. 17 and 18 as determining u*, and we examine the equilibrium possibilities as determined from Eqs. 15 and 16 . For a given u*, there are two equilibrium solutions possible, depending on how G_{n}(v, u*_{1}, N*_{1}, R*)_{v=u*1} = 0. We will refer to these as equilibrium N and equilibrium E.
Equilibrium N is possible if the second term in Eq. 6 is zero to yield R* is obtained from Eq. 16 . However, two equilibriums are possible. If 1 < r < r_{max}, then r in Eq. 16 is replaced by Eq. 11 to solve for R*. On the other hand, if r = r_{max}, then r in Eq. 16 is replaced by r_{max} to solve for R*. Equilibrium E is possible only if the third term in Eq. 6 is zero and r = r_{max}. N*_{1} is obtained from Eqs. 15 and 16 , and R* is obtained by setting the third term in Eq. 6 equal to zero. The expressions obtained are given in Table 1 ⇓ . Note that the equilibrium number of cells is independent of u*_{1} under equilibrium E.
We can now examine both the ecological and evolutionary stability of these equilibrium solutions. Using an eigenvalue analysis to examine the ecological stability and the ESS maximum principle (17) to examine evolutionary stability, we find that both equilibrium solutions of equilibrium N are ecologically stable but only marginally evolutionarily stable. Equilibrium E is both ecologically and evolutionarily unstable.
Parameter Values Used in Figures
We need to specify parameters used to define the distribution functions as well as initial conditions (for time plots) or final conditions for adaptive landscape plots.
Fig. 1 ⇓ .
Values used for normal cell case are given in the Table 2 ⇓ . Units have been defined previously and are not repeated here for brevity. However, N has been normalized to 10^{6} cells (i.e., N = 1000 ⇒ 10^{9} cells).
When n = 1, the critical value for r_{e} is given by r_{e} = 1. It can be shown through simulation (when n = 1) that the normal cells evolve to u*_{1} = [0.5 0.5] corresponding to equilibrium N. The system of Eqs. 8 , 13 , and 10 (when n = 1) drives the system to equilibrium N with values using any feasible initial conditions in the neighborhood of this solution.
Figs. 2 ⇓ and 3 ⇓ .
Tumor cells are introduced by setting n = 2 with the initial conditions , where 1 refers to normal cells, and 2 refers to tumor cells. A background mutation rate is simulated by setting All parameter values for the tumor cells are the same as the normal cells, except that N_{tmean} = 1500.
Figs. 4 ⇓ and 5 ⇓ .
The same initial conditions are used as in Fig. 3 ⇓ with the following changes in parameter values
Footnotes

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 To whom requests for reprints should be addressed, at Department of Radiology, University Medical Center, 1501 North Campbell Avenue, Tucson, AZ 85724. Phone: (520) 6265725; Fax: (529) 6269981; Email: rgatenby{at}radiology.arizona.edu

↵2 The abbreviation used is: ESS, evolutionarily stable strategy.

↵3 J. S. Brown, Y. Cohen, and T. L. Vincent, Adaptive dynamics with vectorvalued strategies, manuscript in preparation.
 Received February 19, 2002.
 Revision received July 7, 2003.
 Accepted July 23, 2003.
 ©2003 American Association for Cancer Research.