Toward an Ising Model of Cancer and Beyond
Abstract
The holy grail of tumor modeling is to formulate theoretical and computational tools that can be utilized in the clinic to predict neoplastic progression and propose individualized optimal treatment strategies to control cancer growth. In order to develop such a predictive model, one must account for the numerous complex mechanisms involved in tumor growth. Here we review resarch work that we have done toward the development of an“Ising model” of cancer. The Ising model is an idealized statisticalmechanical model of ferromagnetism that is based on simple localinteraction rules, but nonetheless leads to basic insights and features of real magnets, such as phase transitions with a critical point. The review begins with a description of a minimalist fourdimensional (three dimensions in space and one in time) cellular automaton (CA) model of cancer in which healthy cells transition between states (proliferative, hypoxic, and necrotic) according to simple local rules and their present states, which can viewed as a strippeddown Ising model of cancer. This model is applied to model the growth of glioblastoma multiforme, the most malignant of brain cancers. This is followed by a discussion of the extension of the model to study the effect on the tumor dynamics and geometry of a mutated subpopulation. A discussion of how tumor growth is affected by chemotherapeutic treatment, including induced resistance, is then described. How angiogenesis as well as the heterogeneous and confined environment in which a tumor grows is incorporated in the CA model is discussed. The characterization of the level of organization of the invasive network around a solid tumor using spanning trees is subsequently described. Then, we describe open problems and future promising avenues for future research, including the need to develop better molecularbased models that incorporate the true heterogeneous environment over wide range of length and time scales (via imaging data), cell motility, oncogenes, tumor suppressor genes and cellcell communication. A discussion about the need to bring to bear the powerful machinery of the theory of heterogeneous media to better understand the behavior of cancer in its microenvironment is presented. Finally, we propose the possibility of using optimization techniques, which have been used profitably to understand physical phenomena, in order to devise therapeutic (chemotherapy/radiation ) strategies and to understand tumorigenesis itself.
pacs:
87.17.Aa, 87.19.xjSalvatore Torquato
Department of Chemistry, Princeton University, Princeton, NJ 08544, USA
Department of Physics, Princeton University, Princeton, NJ 08544, USA
Princeton Center for Theoretical Science, Princeton University, Princeton, NJ 08544, USA
Program in Applied and Computational Mathematics, Princeton University, Princeton, NJ 08544, USA
Princeton Institute for the Science and Technology of Materials, Princeton University, Princeton, NJ 08544, USA
Corresponding author contact information: Salvatore Torquato Tel.: 6092583341 Fax: 6092586746 Email: torquatoprinceton.edu
Short title: Toward an Ising Model of Cancer
Classification numbers: 87.17.Aa, 87.19.xj
Keywords: tumor growth, glioblastoma multiforme, cellular automaton, heterogeneous media, optimization
1 Introduction
Cancer is not a single disease, but rather a highly complex and heterogeneous set of diseases. Dynamic changes in the genome, epigenome, transcriptome and proteome that result in the gain of function of oncoproteins or the loss of function of tumor suppressor proteins underlie the development of all cancers. While some of the mechanisms that govern the transformation of normal cells into malignant ones are rather well understood [1], many mechanisms are either not fully understood or are unknown at the moment. Even if all of the mechanisms could be identified and comprehended, it is not clear progress in understanding cancer could be made without knowledge of how these different mechanisms couple to one another. It has been observed that many complex interactions occur between tumor cells, and between a cancer and the host environment. Multidirectional feedback loops occur between tumor cells and the stroma, immune cells, extracellular matrix and vasculature [2, 3, 4, 5], which are not well understood synergistically. Clearly, our current state of knowledge is insufficient to deduce clinical outcome, not to mention how to control cancer progression in the most malignant forms of cancer.
This suggests that a more quantitative approach to understanding different cancers is necessary in order to control it and increase the lifetime of patients with these deadly diseases. Theoretical/computational modeling of cancer when appropriately linked with experiments and data offers a promising avenue for such an understanding. Such modeling of tumor growth using a variety of different approaches has been a very active area of research for the last two decades or so [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24] but clearly is in its infancy. A diverse number of mechanisms have been explored via such models, and a multitude of computational/mathematical techniques have been employed; see Ref. [25] for a review. These models have the common aim of predicting certain features of tumor growth in the hope of finding new ways to control neoplastic progression. Given a model which yields reproducible and accurate predictions, the effects of different genetic, epigenetic and environmental changes, as well as the impact of therapeutically targeting different aspects of the tumor, can be probed. However, these models must be linked to data from experimental assays in a comprehensive and systematic fashion in order to develop of a quantitative understanding of cancer.
The holy grail of computational tumor modeling is to develop a simulation tool that can be utilized in the clinic to predict neoplastic progression and response to treatment. Not only must such a model incorporate the many feedback loops involved in neoplastic progression, the model must also account for the fact that cancer progression involves events occurring over a range of spatial and temporal scales [15]. A successful model would enable one to broaden the conclusions drawn from existing medical data, suggest new experiments, test new hypotheses, predict behavior in experimentally unobservable situations, and be employed for early detection.
There is overwhelming evidence that cancer of all types are emerging, opportunistic systems [26]. Success in treating various cancers as a selforganizing complex dynamical systems will require unconventional, innovative approaches and the combined effort of an interdisciplinary team of researchers. A lofty longterm goal of such an endeavor is not only to obtain a quantitative understanding of tumorigenesis but to limit and control the expansion of a solid tumor mass and the infiltration of cells from such masses into healthy tissue.
Figure 1: Picture of an Ising model.
Because a comprehensive review of the vast literature concerning biophysical cancer modeling is beyond the scope of this article, we focus on reviewing the work that we have done toward the development of an “Ising model” of cancer. The Ising model is an idealized statisticalmechanical model of ferromagnetism that is based on simple localinteraction rules (see Figure 1), but nonetheless leads to basic insights and features of real magnets, such as phase transitions with a critical point. Toward the goal of developing an analogous Ising model of cancer, we have formulated a fourdimensional (4D) (three dimensions in space and one in time) cellular automaton (CA) model for brain tumor growth dynamics and its treatment [11, 12, 13, 18, 20, 21]. Like the Ising model for magnets, we will see later that this involves local rules for how healthy cells transition into various types of cancer cells. Before describing our computational models for tumor growth, we first briefly summarize several salient features of solid tumor growth as applied to glioblastoma multiforme (GBM), the most malignant of brain cancers.
The rest of paper is organized as follows: In Section 2, some background concerning GBMs and solid tumors in general is presented. In Section 3, a minimalist 4D CA tumor growth model is described in which healthy cells transition between states (proliferative, hypoxic, and necrotic) according to simple local rules and their present states, and apply it to GBMs. This is followed by a discussion of the extension of the model to study the effect on the tumor dynamics and geometry of a mutated subpopulation. How tumor growth is affected by chemotherapeutic treatment is also discussed, including induced resistance. In Section 4, the modification of the CA model to include explicitly the effects of vasculature evolution and angiogenesis on tumor growth are discussed. In Section 5, the effects of physical confinement and heterogeneous environment are described. In Section 6, a simulation tool for tumor growth that merges and improves individual CA models is presented. In Section 7, a descriptions if given of how one might characterize the invasive network organization around a solid tumor using spanning trees. Section 8 discusses some open problems and promising avenues for future research.
2 GBM and Solid Tumor Background
Figure 2: The picture of a tumor in brain.
Glioblastoma multiforme (GBM) (see Figure 2), the most aggressive of the gliomas, is a collection of tumors arising from the glial cells or their precursors in the central nervous system [27].Unfortunately, despite advances made in cancer biology, the median survival time for a patient diagnosed with GBM is only 1215 months, a fact that has not changed significantly over the past several decades [28]. As suggested by its name (i.e., “multiforme”), GBM is complex at many levels of organization [27]. It exhibits diversity at the macroscopic level, having necrotic, hypoxic and proliferative regions. At the mesoscopic level, tumor cell interactions, microvascular remodeling [29] and pseudopalisading necrosis are observed [30]. Further, the discovery that tumor stem cells may be the sole malignant cell type with the ability to proliferate, selfrenew and differentiate introduces yet another level of mesoscopic complexity to GBM [31, 32]. At the microscopic level, GBM cells exhibit point mutations, chromosomal deletions, and chromosomal amplifications [27].
Figure 3: The picture of a MTS.
A substantial amount of research has been conducted to model macroscopic tumor growth either based on microscopic considerations [33, 34, 35]; or in a more phenomenological fashion [6, 36]. One of the early attempts to model empircally the volume of a solid tumor versus time is the Gompertz model, i.e.,
(1) 
where is the volume at time and and are parameters; see Ref. [37] and references therein. Qualitatively, this equation gives exponential growth at small times which then saturates at large times (decelerating growth). In particular, this model considers the tumor as an oversized idealized multicellular tumor spheroid (see Figure 3), which is stage of early tumor growth. We note that modeling an ideal tumor as an oversized spheroid is especially suited for GBM, since this tumor, like a large MTS, comprises large areas of central necrosis surrounded by a rapidly expanding shell of viable cells (Figure 2). However, we note that real tumors always possess much more complex morphology. More importantly, Gompertziangrowth models are very limited; they only capture gross features of tumor growth and cannot explain their underlying “microscopic” mechanisms.
One of the hallmarks of highgrade malignant neuroepithelial tumors, such as glioblastoma multiforme (GBM), is the regional heterogeneity, i.e., the relatively large number of clonal strains or subpopulations present within an individual tumor of monoclonal origin [38, 39, 40]. Each of these strains is characterized by specific properties, such as the rate of division or the level of susceptibility to treatment [41]. Therefore the growth dynamics of a single tumor are determined by the relative behaviors of each separate subpopulation. For example, the appearance of a rapidly dividing strain can substantially bias tumor growth in that direction. Tumors supposedly harbor cells with an increased mutation rate, which indicates that these tumors are genetically unstable [42, 43, 44]. Genetic and epigenetic events throughout the tumor may occur randomly or be triggered by environmental and intrinsic stresses. The continuing existence of a subpopulation, however, depends primarily on the subpopulation’s ability to compete with the dominant population in its immediate vicinity.
Clonal heterogeneity within a tumor has been shown to have very pronounced effects on treatment efficacy [45, 46]. Treatment resistance is itself a complex phenomenon. There is no single cause of resistance, and many biochemical aspects of resistance are poorly understood. Chemoresistant strains can either be resistant to a single drug or drug family (individual resistance), or they can be resistant to an array of agents (multidrug/pleotropic resistance) [47]. Cellular mechanisms behind multidrug resistance include increased chemical efflux and/or decreased chemical influx, such as with Pglycoproteinmediated (Pgp) drug resistance [48, 49].
Complicating the situation further, resistance can arise at variable times during tumor development. Some tumors are resistant to chemotherapy from the onset. This has been described as inherent resistance, because it exists before chemotherapeutic drugs are ever introduced. In other cases, however, treatment initially proves successful, and only later does the tumor prove resistant. This is an example of acquired resistance, as it develops in response to treatment [47]. There are at least two possible mechanisms for this type of tumor behavior. Acquired resistance may result from a small number of resistant cells that are gradually selected for throughout the course of treatment. At the same time, there is also evidence suggesting that chemotherapeutic agents may induce genetic or epigenetic changes within tumor cells, leading to a resistant phenotype [50]. Other studies indicate that chemotherapy may increase cellular levels of Pgp mRNA and protein in various forms of human cancer [51, 52]. A tumor’s response to radiation therapy can also depend on underlying genetic factors. A cell’s inherent radioresistance may stem from the efficiency of DNA repair mechanisms in sublethally damaged cells [53, 54, 55].
While the properties of GBM cells are very important in understanding growth dynamics, just as important are the properties of the environment in which the tumor grows. For example, GBMs grow in either the brain or spinal cord, and are therefore confined by the shape and size of these organs [20]. Another example of the importance of accounting for the host environment relates to the vascular structure of the brain.
Recent research evidence suggests that tumors arising in vascularized tissue such as the brain do not originate avascularly [29], as originally thought. Instead, it is hypothesized that glioma growth is a process involving vessel cooption, regression and growth. Three key proteins, vascular endothelial growth factor (VEGF) and the angiopoietins, angiopoietin1 (Ang1) and angiopoietin2 (Ang2), are required to mediate these processes [29]. A picture of what likely occurs during the process of glioma vascularization has been summarized by Gevertz and Torquato [18]. As a malignant mass grows, the tumor cells coopt the mature vessels of the surrounding brain that express constant levels of bound Ang1. Vessel cooption leads to the upregulation in Ang2 and this shifts the ratio of bound Ang2 to bound Ang1. In the absence of VEGF, this shift destabilizes the coopted vessels within the tumor center and marks them for regression [56, 57]. Vessel regression in the absence of vessel growth leads to the formation of hypoxic regions in the tumor mass. Hypoxia induces the expression of VEGF, stimulating the growth of new blood vessels [58]. This robust angiogenic response eventually rescues the suffocating tumor. Glioma growth dynamics remain intricately tied to the continuing processes of vessel regression and growth.
Tumor cell invasion is a hallmark of gliomas [59]. Individual glioma cells have been observed to spread diffusely over long distances and can migrate into regions of the brain essential for the survival of the patient [27]. While MRI scans can recognize mass tumor lesions, these scans are not sensitive enough to identify malignant cells that have spread well beyond the tumor region [60]. Typically, when a solid tumor is removed, these invasive cells are left behind and tumor recurrence is almost inevitable [27].
Numerous models have been developed to model certain tumor behavior or characteristics with a great deal of mathematical rigor (e.g., in the form of coupled differential equations). However, with such approaches, the sets of equations that govern tumor behavior often do not correspond to the characteristics of individual tumor cells. An important goal of studying tumor development is to illustrate how their macroscopic traits stem from their microscopic properties. In addition, most of the equations are problemspecific, which limits their utility as general tools for tumor study. Another potential challenge is that tumor models should be appreciated by as diverse an audience as possible. Ideally, the mathematical complexity that allows theoreticians to analyze subtle aspects of it should not be an obstacle for clinicians who treat GBM. A model that accounts for complex tumor behavior with relative mathematical ease could be valuable.
To this end, we have developed what appears to be a powerful cellular automaton (CA) computational tool for tumor modeling. Based on a few salient set of microscopic parameters, this CA model can realistically model the macroscopic tumor behavior, including growth dynamics, emergence of a subpopulation as well as the effects of tumor treatment and resistance [11, 12, 13]. This model has been extended to study the effects of vasculature evolution on early tumor growth and to simulate tumor growth in confined heterogeneous environments [18, 20, 21]. We have also developed mathematical models to characterize the invasive network organization around a solid tumor [61].
3 Toward an Ising Model of Cancer Growth
In this section, we describe a fourdimensional (4D) cellular automaton (CA) model that we have developed that describes tumor growth as a function of time, using the fewest number of microscopic parameters [11, 12, 13]. We refer to this as a minimalist fourdimensional (4D) model because it involves three spatial dimensions and one dimension in time with the goal of capturing the salient features of tumor growth with a minimal number of parameters. The algorithm takes into account that this growth starts out from a few cells, passes through a multicellular tumor spheroid (MTS) stage (Figure 3) and proceeds to the macroscopic stages at clinically designated timepoints for a virtual patient: detectable lesion, diagnosis and death. This 4D CA algorithm models macroscopic GBM tumors as an enormous idealized MTS, mathematically described by a Gompertzfunction given by Eq. (1), since this tumor, like a large MTS, comprises large areas of central necrosis surrounded by a rapidly expanding shell of viable cells (Figure 2). In accordance with experimental data, the algorithm also implicitly takes into account that invasive cells are continually shed from the tumor surface and implicitly assumes that the tumor mass is wellvascularized during the entire process of growth. The effects of vasculature evolution are considered explicitly in Sections 5 and 7.
3.1 A 4D Cellular Automaton Model
A CA model is a spatially and temporally discrete model that consists of a grid of cells, with each cell being in one of a number of predefined states. The state of a cell at a given point in time depends on the state of itself and its neighbors at the previous discrete time point. Transitions between states are determined by a set of local rules. The simulation is designed to predict clinically important criteria such as the fraction of the tumor which is able to divide (GF), the nonproliferative ( arrest) and necrotic fractions, as well as the rate of growth (volumetric doubling time) at given radii. Furthermore, this CA model enables one to study emergence of a subpopulation due to cell mutations as well as the effects of tumor treatment and resistance. The general CA model includes both a proliferation routine which models tumor growth by cell division and a treatment routine which models the cell response to treatment and cell mutations. It also incorporates a novel adaptive automaton cell generation procedure. In particular, the CA model is characterized by several biologically important features:

The model is able to grow the tumor from a very small size of roughly 1000 real cells through to a fully developed tumor with cells. This allows a tumor to be grown from almost any starting point, through to maturity.

The thickness of different tumor layers, i.e. the proliferative rim and the nonproliferative shell, are linked to the overall tumor radius by a 2/3 power relation. This reflects a surface area to volume ratio, which can be biologically interpreted as nutrients diffusing through a surface.

The discrete nature of the model and the variable density lattice allow us to control the inclusion of mutant “hot spots” in the tumor as well as variable cell sensitivity/resistance to treatment. The variable density lattice will allow us to look at such an area at a higher resolution.

Our inclusion of mechanical confinement pressure enables us to simulate the physiological confinement by the skull at different locations within the brain differently.
Our CA algorithm can be broken into three parts: automaton cell generation, the proliferation routine and the treatment routine. In the ensuing discussions, we first present the three parts of our algorithm. Then we show that the our model reflects a test case derived from the medical literature very well, proving the hypothesis that macroscopic tumor growth behavior may be modeled with primarily microscopic data.
3.1.1 Cellular Automaton Cell Generation
The first step of the simulation is to generate the automaton cells. The underlying lattice for the algorithm is the Delaunay triangulation, which is the dual lattice of the Voronoi tessellation [11, 62]. In order to develop the automaton cells, a prescribed number of random points are generated in the unit square using the process of random sequential addition (RSA) of hard circular disks. In the RSA procedure, as a random point is generated, it is checked if the point falls within some prescribed distance from any other point already placed in the system [11, 62]. Points that fall too close to any other point are rejected, and all others are added to the system. Each cell in the final Voronoi lattice will contain exactly one of these accepted sites. The Voronoi cell is defined by the region of space nearer to a particular site than any other site. In twodimensions, this results in a collection of polygons that fill the plane.
Figure 4: Voronoi Cells
Because a real brain tumor grows over several orders of magnitude in volume, the lattice was designed to allowed continuous variation with the radius of the tumor. The density of lattice sites near the center was significantly higher than that at the edge. A higher site density corresponds to less real cells per automaton cell, and so to a higher resolution. The higher density at the center enables us to simulate the flat smalltime behavior of the Gompertz curve. In the current model, the innermost automaton cells represent roughly 100 real cells, while the outermost automaton cells represent roughly real cells. The average distance between lattice sites was described by the following relation:
(2) 
in which is the average distance between lattice sites and is the radial position at which the density is being measured. This relation restricts the increase in the number of proliferating cells as the tumor grows. Note that when modeling the effects of vasculature evolution discussed in the following, a a uniform lattice is used for which each automaton cell includes approximately 10 real cancer cells.
3.1.2 Minimalist 4D Proliferation Algorithm
Figure 5: Idealized tumor
Functions within the model (time dependent)  

Average overall tumor radius (see Appendix)  
Proliferative rim thickness (determines growth fraction)  
Nonproliferative thickness (determines necrotic fraction)  
Probability of division (varies with time and position)  
Growth parameters  
Base probability of division, linked to celldoubling time (0.192)  
Base necrotic thickness, controlled by nutritional needs ()  
Base proliferative thickness, controlled by nutritional needs ()  
Maximum tumor extent, controlled by pressure response () 
The proliferation algorithm is designed to allow a tumor consisting of a few automaton cells, representing roughly 1000 real cells, to grow to a full macroscopic size. An idealized model of a macroscopic tumor is an almost spherical body consisting of concentric shells of necrotic, nonproliferative and proliferative regions (see Figure 5). The four microscopic growth parameters of the algorithm are , , , and reflecting, respectively, the rate at which the proliferative cells divide, the nutritional needs of the nonproliferative and proliferative cells, and the response of the tumor to mechanical pressure within the skull. In addition, there are four key timedependent quantities that determine the dynamics of the tumor, i.e., , , , giving, respectively, the average overall tumor radius, proliferative rim thickness, nonproliferative thickness and probability of division. These quantities are based on the four parameters (, , , ) and are calculated according to the following algorithm.

Initial setup: The cells within a fixed initial radius of the center of the grid are designated proliferative. All other cells are designated as nontumorous.

Time is discretized and incremented, so that at each time step:

Each cell is checked for type: nontumorous or (apoptotic and) necrotic, nonproliferative or proliferative tumorous cells.

Nontumorous cells and tumorous necrotic cells are inert.

Nonproliferative (growtharrested) cells more than a certain distance, , from the tumor’s edge are turned necrotic. This is designed to model the effects of a nutritional gradient. The edge of the tumor is taken to be the nearest nontumorous cell, i.e.,
(3) 
Proliferative cells are checked to see if they will attempt to divide according to the probability of division, , which is influenced by the location of the dividing cell, reflecting the effects of mechanical confinement pressure. This effect requires the use of an additional parameter, the maximum tumor extent, . is given by
(4) 
If a cell attempts to divide, it will search for sufficient space for the new cell beginning with its nearest neighbors and expanding outwards until either an empty (nontumorous) space is found or nothing is found within the proliferation radius, . The radius searched is calculated as:
(5) 
If a cell attempts to divide but cannot find space it is turned into a nonproliferative cell.


After a predetermined amount of time has been stepped through, the volume and radius of the tumor are plotted as a function of time.

The type of cell contained in each grid are saved at given times.
Figure 6: An illustration of the proliferation algorithm.
The above simulation procedure is also illustrated in Figure 6. We note that the redefinition of the proliferative to nonproliferative transition implemented in the algorithm is one of the most important new features of the model. They allow a larger number of cells to divide, since cells no longer need to be on the outermost surface of the tumor to divide. In addition, it ensured that cells that cannot divide are correctly labeled as such. Table 1 summarizes the important timedependent functions calculated by the proliferation algorithm and the constant growth parameters used. The readers are referred to Ref. [11] for the detailed description of the algorithm and parameters.
3.1.3 Extending the 4D CA Model to Study Emergence of a Subpopulation
Malignant brain tumors such as GBM generally consist of a number of distinct subclonal populations. Each of these subpopulations, arising from the constant genetic and epigenetic alteration of existing cells in the rapidly growing tumor, may be characterized by its own behaviors and properties. However, since each single cell mutation only leads to a small number of offspring initially, very few newly arisen subpopulations survive more than a short time. Kansal et al [12] have extended the CA to quantify “emergence,” i.e. the likelihood of an isolated subpopulation surviving for an extended period of time. Only mutations affecting the rate of cellular division were considered in this rendition of the model. In addition, only competition between clones was taken into account; there were no cooperative effects included, although such effects can easily be incorporated.
The simulation procedure is as follows: an initial tumor composed entirely of cells of the primary clonal population is introduced, which is allowed to grow using the proliferation algorithm until it reaches a predetermined average overall radius. Then, a single (or a small number of) automaton cell is changed from the primary strain to a secondary strain with an altered probability of division, which represents very small fractions of the total population of proliferative tumor cells and the tumor is allowed to continue to grow using the proliferation algorithm. It is important to note that this does not represent a single mutation event but rather a mutation event that results in a subpopulation reaching a size dictated by the limits of the lattice resolution employed (i.e., a specified number of cells).
The behavior of the secondary strain was characterized in terms of two properties: the degree and the relative size of the initial population of mutated cells, i.e.,
(6) 
which represents the ratio between the base probability of division of the new clone, , and that of the original clone, ; and
(7) 
Positive, negative and no competitive advantages are respectively conferred for , , and . The initial value , i.e., , is a parameter of the model reflecting the size of the mutated region introduced.
3.1.4 Extending the 4D CA Model to Study Treatment
Besides the four growth parameters in the minimalist 4D CA model, three additional parameters for treatment were subsequently introduced: , , and , the values of which reflect, respectively, the proliferative cells’ treatment sensitivity, the nonproliferative cells’ treatment sensitivity, and the mutational response of the tumor cells to treatment [13]. Furthermore, there are three additional timedependent quantities , and , giving respectively fraction of proliferative cells that die upon treatment (equivalent to ), fraction of nonproliferative cells that die upon treatment (equivalent to ) and volume fraction of mutated living cells. These parameters are summarized in Table 2 and a detailed discussion is given in Ref. [13].
Treatment parameters  

Governs the proliferative cells’ response at each  
instance of treatment ()  
Allows for different treatment responses between  
proliferative and nonproliferative cells ()  
Fraction of surviving proliferative cells  
that mutate in response to treatment ()  
Other Terms  
Fraction of proliferative cells that die  
upon treatment; equivalent to  
Fraction of nonproliferative cells that die  
upon treatment; equivalent to  
Volumetric fraction of living cells (proliferative and  
nonproliferative) belonging to the secondary strain 
In the simulation, treatment was introduced as “periodic impulse”, i.e., a small tumor mass is introduced which is intended to represent a GBM after successful surgical resection and allowed to grow using the proliferation algorithm; then treatment is applied and considered effective at discrete time points. In particular, the simulation proceeds through the proliferative steps until every week timepoint, at which time the treatment routine is introduced:

After the last round of cellular division, each proliferative cell is checked to see if it is killed by the treatment. The probability of death for a given proliferative cell is given by
(8) where is the proliferative treatment parameter. Dead proliferative automaton cells are converted to healthy cells.

Each nonproliferative cell is checked to see if it is killed. The probability of death for a given nonproliferative cell is given by
(9) where is the nonproliferative treatment parameter and is a fraction of . A nonproliferative cell is converted to a necrotic cell upon death.

Each surviving nonproliferative cell is checked to see if it is within the proliferative thickness of a healthy cell (i.e. the tumor surface). If so, the nonproliferative cell is converted back to a proliferative cell.

All proliferative cells (including newlydesignated ones) are checked for mutations for treatment resistance with probability . A new is randomly generated for mutated cells while remains constant.
Clinically, GBM treatment consists of both radiation therapy and chemotherapy. However, in our model we do not distinguish between the separate effects of these two methods. The tumors’ response to all treatment is captured by the treatment algorithm. Moreover, this response is assumed to be instantaneous at each fourweek time point.
4 Putting the 4D CA Model Through Its Paces
4.1 A Test Case for Proliferation Algorithm
The tumor growth data generated via the minimalist 4D CA proliferation algorithm was compared with available experimental data for an untreated GBM tumor from the medical literature [11]. The parameters compared were cell number, growth fraction, necrotic fraction and volumetric doubling time, which are used to determine a tumor’s level of malignancy and the prognosis for its future growth. Because it is impossible to determine the exact time a tumor began growing, the medical data are listed at fixed radii. The different cell fractions used were extrapolated from the spheroid level and compared to data published for cell fractions at macroscopic stages.
Summarized in Table 3 is the comparison between simulation results and data (experimental, as well as clinical) taken from the medical literature (see Ref. [11] for detailed references). The simulation data were created using a tumor which was grown from an initial radius of 0.1 mm. The following parameter set (see Table 1) was used:
, , ,
This value of corresponds to a celldoubling time of 4.0 days. The parameters and have been chosen to give a growth history that quantitatively fits the test case. The specification of these parameters corresponds to the specification of a clonal strain. The parameter was similarly chosen to match the test case history. In this case, however, the fit is relatively insensitive to the value of , as long as the parameter is somewhat larger than the fatal radius in the test case. On the whole, the simulation data reproduces the test case very well. The virtual patient would die roughly 11 months after the tumor radius reached 5 mm and 3.5 months after the expected time of diagnosis. The fatal tumor volume is about .
Figure 7: Crosssections of a Growing MonoClonal Tumor.
Figure 8: A cutaway view of the simulated tumor.
Figure 9: The Volume and Radius of the Developing Tumor.
Central crosssections of the tumors are shown in Figure 7, in which the growth of the tumor can be followed graphically over time. Here necrotic cells are labeled with black, nonproliferative tumorous cells with light gray and proliferative tumor cells with dark gray. A cutaway view of the simulated tumor is shown in Figure 8. As expected in this idealized case, the tumor is essentially spherical, within a small degree of randomness. The high degree of spherical symmetry ensures that the central crosssection is a representative view. The volume and radius of the developing tumor are shown versus time in Figure 9. Note that the virtual patient dies while the untreated tumor is in the rapid growth phase.
4.2 Modeling the Emergence of A Subpopulation
Recall that the parameter reflects the degree of advantage of the mutated subpopulation over the primary clone (positive, negative and no competitive advantages are respectively conferred for , , and ) and the initial value , i.e., , is a parameter of the model reflecting the size of the mutated region introduced. A subpopulation is considered to have emerged once it comprises of the actively dividing cell population or if it remains in the actively dividing state once the tumor has reached a fully developed size. Numerous simulations (at least 100) were run at each parameter set by Kansal et al [12] in order to calculate the expected probabilities of emergence, i.e.,
(10) 
along with confidence intervals, , defined as
(11) 
where represents the observed probability of emergence in trials. We note that the probability of emergence is actually a conditional probability: it is the probability that a subpopulation with a mutation of degree emerges given that a region of relative size has mutated.
Figure 10: P vs. alpha and cutaway view of simulated tumor with a subpopulation.
The results represented were run with a parameter set in which
, , ,
,
for the primary strain, in a simulation in which each time step represents one day [12]. Figure 10 depicts the observed probability of emergence, , for a subpopulation of initial size as a function of , which gives an approximation of the true, asymptotic, probability of emergence. Also shown in Figure 10 is a cutaway view of the simulated tumor with a subpopulation. Not surprisingly, is a monotonically increasing function that tends to 0 for and to 1 as become significantly greater than 1.
Perhaps the most striking feature of these results is that there is a nonzero probability of emergence for a very small population with no growth rate advantage, or even with a small disadvantage (i.e. ). This suggests that a mutated subpopulation may arise even without any growth advantage. These populations could represent “dormant” clones which confer an advantage not being selected for at the time. An example would be the appearance of hypoxia tolerant or even treatment resistant clones. It should be stressed that populations with less competitive advantages over other tumor strains can have a nonzero probability of emergence especially if they are localized in space, which leads to a minimum surface area between the two populations per unit tumor volume. In this way, the population with smaller competitive advantage can compete more effectively. We will see in the next subsection that this same principle is at work when resistance is induced due to treatment. It was also found that the emergence probability is a monotonically increasing function in and has a logarithmic dependence on [12].
Figure 11: Effects of the subpopulation on the tumor geometry.
Figure 12: Effects of the subpopulation on growth history.
Figure 11 shows the effects of growth of the subpopulation on the tumor geometry. It can be seen clearly that the center of mass of the tumor is significantly shifted by the emergence of the subpopulation. Another example of the importance of subpopulations is depicted in Figure 12 [12]. In this example, a diagnosis was made (on day ) giving information about the macroscopic size and growth rate of the tumor. From this information three possible growth histories of the tumor are plotted. One is the time history of the tumor with an emergent subpopulation. The others represent limiting cases, each with a monoclonal tumor of either the primary (“base ”) or secondary (“high ”) clonal strain. Note that at the time of diagnosis all three scenarios have very similar dynamics. So any of the three histories is a reasonable prediction given only size and growth rate information. However, estimating a fatal tumor volume to be 65 and defining the survival time to be the time required to reach this volume, the base case mispredicts survival times to be 90 days, which is 30 days more than the 60 days of the “true” course.
It is noteworthy that from this perspective the overall future growth dynamics of the entire tumor closely follows that of the most aggressive case, indicating that the more aggressive clone dominates overall outcome and should therefore also define the appropriate treatment. This finding supports the current practice in pathology of grading tumors according to the most malignant area (i.e. population) found in any biopsy material. Although of less clinical significance, the high case similarly mispredicts the past history of the tumor. If the diagnosis had been made earlier, the base case would yield still worse future predictions. Similarly the “high” case would yield worse past predictions for a diagnosis made at a later time. The predictive errors arising from the assumption of a monoclonal tumor indicate how important an accurate estimate of the clonal composition of a tumor is in establishing a complete history and prognosis. Note that the numbers given here are intended to show the scale of the inaccuracy possible, not to reflect any data extracted from actual patients.
4.3 Modeling the Effects of Tumor Treatment and Resistance
Combining the proliferation algorithm and the treatment algorithm, the behavior of tumors that are able to develop resistance throughout the course of treatment were investigated [13]. Recall that additional parameters were introduced in the treatment routine: , , and , the values of which reflect, respectively, the proliferative cells’ treatment sensitivity, the nonproliferative cells’ treatment sensitivity, and the mutational response of the tumor cells to treatment (see Table 2).
These investigation consisted of three individual case studies. In Case 1, the growth dynamics of monoclonal tumors are studied to determine how tumor behavior is affected by the treatment parameters and . Case 2 builds upon this information, analyzing the behavior of twostrain tumors. Here, a secondary treatmentresistant strain exists alongside a primary treatmentsensitive strain. A secondary subpopulation was introduced at the onset of each simulation, initializing it in different spatial arrangements and at several (small) relative volumes. In both Cases 1 and 2, no additional subpopulations arise in the tumors once the simulation has begun (i.e. ). In Case 3, however, tumors were studied that were capable of undergoing resistance mutations in response to each round of treatment (). In these simulations, the growth and morphology of the tumors were analyzed in relation to the fraction of mutating cells.
Here we only report on the results of Cases 2 and 3. In Case 2, the smaller subpopulation of a secondary treatmentresistant strain was initially spatially distributed in two different ways on the tumor surface that primarily consist of the primary treatmentsensitive strain: a localized and scattered scenario, reflecting possible effects of the result of surgery, for example (see Figure 13). In the simulation, the tumors were initialized as a single strain, i.e., monoclonally with and and treatment was introduced every four weeks while the tumor is growing from a small mass with a radius of 4mm, corresponding to approximately of surgical volume resection. For the scattered resistance scenario, the resistant strain was found to compete more effectively with the sensitive strain and it was shown that the initial number of resistant cells were not a significant indicator or prognosis.
Figure 13: Spatial distributions of the resistance strain.
These conclusions may at first glance seem to contradict the those reported by by Kansal et al [12]. Recall that in this work the selection pressure was different (growthrate competition versus treatment effects). Moreover, the roles of the primary and secondary strains are reversed in the Case 2 example: the primary strain possessed a competitive advantage over the secondary strain. Nevertheless, the conclusions of both papers [12, 13] follow precisely the same principle. The proliferative ability of a strain with a competitive advantage varies directly with its contact area with the less comptetive strain per cell.
Unlike Case 2, the tumors in Case 3 begin the simulations as a single strain. Here, however, treatment can induce the appearance of mutant strains (). In these simulations, the growth and morphology of the tumors were analyzed in relation to the fraction of mutating cells. The tumors in Case 3 are all initialized monoclonally with and . With this initial value, nearly every mutant strain that arises from the initial population will posses a lower value. This is not to suggest that all induced mutations must possess increased resistance. This fact here merely stems from the initial sensitive tumor under consideration.
At first, the tumors in Case 3 will develop like treatmentsensitive, monoclonal tumors; growth will then accelerate as resistant cells begin to dominate. This corresponds to a case of acquired resistance via induced (genetic and epigenetic) mutations. Overall, the tumor dynamics here are more variable than in Cases 1 and 2. When a new strain appears, it begins as a single automaton cell. Unlike Case 2, not all new strains will be able to proliferate to an appreciable extent. Some are overwhelmed by the parent strain from which they arise.
The mean survival time of the tumors were determined as a function of and these data are summarized in Figure 14. plots this data; From to , the survival times vary nearly logarithmically with . When , the mean time is near 27 months, as most tumors remain monoclonal (or nearly monoclonal) with , . As increases, resistant strains appear more commonly and survival times fall.
Figure 14: Survival times associated with continuously mutating tumors.
One of the more intriguing observations in this case involves the gross morphology of the mutating tumors. Their threedimensional geometries exhibit an interesting dependence on the value of . Figure 15 presents representative images of the fullydeveloped tumors for small, intermediate and large fractions of mutated proliferative cells after treatment. For small (left panel of Figure 15), some tumors develop a secondary strain while others do not. The tumors that remain monoclonal maintain their spherical geometry. When a resistant subpopulation does develop, it appears as a lobe on the parental tumor. For intermediate , resistant subpopulations consistently arise from the parental strain. The middle panel of Figure 15 depicts a typical tumor, whose geometries consistently deviate from an ideal sphere. These tumors are multilobed in appearance, and the original strain is commonly overwhelmed. However, when is large, the geometric trend reverses, i.e., the tumors (right panel of Figure 15) again appear more spherical, despite the fact that they experience the greatest fraction of mutations per treatment event. These images suggest that extreme mutational responses can lead to similar macroscopic geometries. Nonspherical geometries result from intermediate values.
Figure 15: Images of continuously mutating tumors.
5 Modeling the Effects of Vasculature Evolution
As pointed out in the Introduction, there are complex interactions occurring between between a tumor and the host environment, which makes it very difficult in predicting clinical outcome, even if mutations responsible for oncogenesis that determine tumor growth are beginning to be understood. These interactions include the effects of vasculature evolution on tumor growth, the organimposed physical confinement as well as the host heterogeneity. While the three studies described in the previous section were successful at analyzing and characterizing GBM growth both with and without treatment, in each case, the CA model made the simplifying assumption that the tumor mass was wellvascularized (the vascular network and angiogenesis were implicitly accounted for) and the effects of mechanical confinement were limited to one parameter (), which allowed for growth of spherically symmetric tumors with a maximum radius. Sphericallike growth is realistic provided that the environment is effectively homogeneous, but heterogeneous environments will cause apshericallyshaped tumors.
In order to incorporate a greater level of microscopic detail, a 3D (two dimensions in space and one in time) hybrid variant of the original CA model that allows one to study how changes in the tumor vasculature due to vessel cooption, regression and sprouting influence GBM was developed by Gevertz and Torquato [18]. This computational algorithm is based on the cooptionregressiongrowth experimental model of tumor vasculature evolution [29, 56]. In this model, as a malignant mass grows, the tumor cells coopt the mature vessels of the surrounding tissue that express constant levels of bound angiopoietin1 (Ang1). Vessel cooption leads to the upregulation of the antagonist of Ang1, angiopoietin2 (Ang2). In the absence of the antiapoptotic signal triggered by vascular endothelial growth factor (VEGF), this shift destabilizes the coopted vessels within the tumor center and marks them for regression [29, 56]. Vessel regression in the absence of vessel growth leads to the formation of hypoxic regions in the tumor mass. Hypoxia induces the expression of VEGF, stimulating the growth of new blood vessels.
A system of reactiondiffusion equations was developed to track the spatial and temporal evolution of the aforementioned key factors involved in blood vessel growth and regression [18] (see Section 6 for a detailed description). Based on a set of algorithmic rules, the concentration of each protein and bound receptor at a blood vessel determines if a vessel will divide, regress or remain stagnant. The structure of the blood vessel network, in turn, is used to estimate the oxygen concentration at each cell site. Oxygen levels determine the proliferative capacity of each automaton cell. The reader is referred to [18] for the full details of this algorithm. The model proved to quantitatively agree with experimental observations on the growth of tumors when angiogenesis is successfully initiated and when angiogenesis is inhibited. Further, due to the biological details incorporated into the model, the algorithm was used to explore tumor response to a variety of single and multimodal treatment strategies [18].
6 Modeling the Effects of Physical Confinement and Heterogeneous Environment
An assumption made in both the original CA algorithm and the one that explicitly incorporates vascular evolution is that the tumor is growing in a spherically symmetric fashion. In a study performed by Helmlinger et al [63], it was shown that neoplastic growth is spherically symmetric only when the environment in which the tumor is developing imposes no physical boundaries on growth. In particular, it was demonstrated that human adenocarcinoma cells grown in a 0.7% gel that is placed in a cylindrical glass tube develop to take on an ellipsoidal shape, driven by the geometry of the capillary tube. However, when the same cells are grown in the same gel outside the capillary tube, a spherical mass develops [63]. This experiment clearly highlights that the assumption of radially symmetric growth is only valid when a tumor grows in an unconfined or spherically symmetric environment.
Since many organs, including the brain and spinal cord, impose nonradially symmetric physical confinement on tumor growth, the original CA algorithm was modified to incorporate boundary and heterogeneity effects on neoplastic progression [20]. The first modification that was made to the original algorithm was simply to specify and account for the boundary that is confining tumor growth. Several modifications were made to the original automaton rules to account for the impact of this boundary on neoplastic progression. The original CA algorithm imposed radial symmetry in order to determine whether a cancer cell is proliferative, hypoxic, or necrotic. The assumption of radially symmetric growth was also utilized in determining the probability a proliferative cell divides. In order to allow tumor growth in any confining environment, all assumptions of radial symmetry from the automaton evolution rules were removed. It was demonstrated that models that do not account for the geometry of the confining boundary and the heterogeneity in tissue structure lead to inaccurate predictions on tumor size, shape and spread (the distribution of cells throughout the growthpermitting region). The readers are referred to Ref. [20] for the details of this investigation, but an illustration of confinement effects are given in the next section.
7 A Merged Tool for Growing Heterogeneous Tumors In Silico
7.1 Algorithmic Details
Each of the previously discussed algorithms were designed to answer a particular set of questions and successfully served their purpose. Hence, Gevertz and Torquato [21] merged each algorithm into a single cancer simulation tool that would not only accomplish what each individual algorithm had accomplished, but had the capacity to have emergent properties not identifiable prior to model integration. In developing the merged algorithm, some modifications were made to the original automaton rules to more realistically mimic tumor progression. The merged simulation tool is summarized as follows:

Automaton cell generation: A Voronoi tessellation of random points generated using the nonequilibrium procedure of random sequential addition of hard disks determines the underlying lattice for our algorithm [11, 62]. Here a uniform density lattice is used instead of the lattice with variable density. Each automaton cell created via this procedure represents a cluster of a very small number of biological cells ().

Define confining boundary: Each automaton cell is divided into one of two regimes: nonmalignant cells within the confining boundary and nonmalignant cells outside of the boundary.

Healthy microvascular network: The blood vessel network which supplies the cells in the tissue region of interest with oxygen and nutrients is generated using the random analog of the Krogh cylinder model detailed in Ref. [18]. One aspect of the merger involved limiting blood vessel development to the subset of space in which tumor growth occurs.

Initialize tumor: Designate a chosen nonmalignant cell inside the growthpermitting environment as a proliferative cancer cell.

Tumor growth algorithm: Time is then discretized into units that represent one real day. At each time step:

Solve PDEs: A previouslydeveloped system of partial differential equations [18] is numerically solved one day forward in time. The quantities that govern vasculature evolution, and hence are included in the equations, are concentrations of VEGF (), unoccupied VEGFR2 receptors (), the VEGFR2 receptor occupied with VEGF (), Ang1 (), Ang2 (), the unoccupied angiopoietin receptor Tie2 (), the Tie2 receptor occupied with Ang1 () and the Tie2 receptor occupied with Ang2 (). The parameters in these equations include diffusion coefficients of protein (), production rates and , carrying capacities , association and dissociation rates ( and ) and decay rates . Any term with a subscript denotes an indicator function; for example, is a proliferative cell indicator function. It equals 1 if a proliferative cell is present in a particular region of space, and it equals 0 otherwise. Likewise, is the hypoxic cell indicator function, is necrotic cell indicator function and is the endothelial cell indicator function. The equations solved at each step of the algorithm are:
(12) (13) (14) (15) (16) (17) (18) (19) In these equations, represents the concentration of hypoxic cells and represents the endothelial cell concentration per blood vessel. The system of differential equations contains 21 parameters, 13 of which were taken from experimental data. Parameters were unable to be found in the literature were estimated. For more details on the parameter values, as well as information on the initial and boundary conditions and the numerical solver, the reader is referred to Ref [18].

Vessel Evolution: Check whether each vessel meets the requirements for regression or growth. Vessels with a concentration of bound Ang2 six times greater than that of bound Ang1 regress [57], provided that the concentration of bound VEGF is below its critical value. Vessel tips with a sufficient amount of bound VEGF sprout along the VEGF gradient.

Nonmalignant Cells: Healthy cells undergo apoptosis if vessel regression causes its oxygen concentration to drop below a critical threshold (more particularly, if the distance of a healthy cell from a blood vessel exceeds the assumed diffusion length of oxygen, 250 m ). Further, nonmalignant cells do not divide in the model. While nonmalignant cell division occurs in some organs, a hallmark of neoplastic growth is that tumor cells replicate significantly faster than the corresponding normal cells. Hence, we work under the simplifying assumption that nonmalignant division rates are so small compared to neoplastic division rates that they become relatively unimportant in the time scales being considered. In the cases where this assumption does not hold, nonmalignant cellular division would have to be incorporated into the model.

Inert Cells: Tumorous necrotic cells are inert. This assumption is certainly valid for the tumor type that motivated this modeling work, glioblastoma multiforme. In the case of glioblastoma, the presence of necrosis is an important diagnostic feature and, in fact, negatively correlates with patient prognosis.

Hypoxic Cells: A hypoxic cell turns proliferative if its oxygen level exceeds a specified threshold [18] and turns necrotic if the cell has survived under sustained hypoxia for a specified number of days. In the original algorithms, the transition from hypoxia to necrosis was based on an oxygen concentration threshold. However, given that cells (both tumorous and nonmalignant alike) have been shown to have a limited lifespan under sustained hypoxic conditions, a temporal switch more accurately describes the hypoxic to necrotic transition. Thus, a novel aspect of the merged algorithm is a temporal hypoxic to necrotic transition. It has been measured that human tumor cells remain viable in hypoxic regions of a variety of xenografts for 410 days [18]. In our simulations, we will use the upperend of this measurement and assume that tumor cells can survive under sustained hypoxia for 10 days.

Proliferative Cells: A proliferative cell turns hypoxic if its oxygen level drops below a specified threshold. However, if the oxygen level is sufficiently high, the cell attempts to divide into the space of a viable nonmalignant cell in the growthpermitting region. The probability of division is given by:
(20) where is the base probability of division, is the distance of the dividing cell from the geometric center of the tumor and is the distance between the closest boundary cell in the direction of tumor growth and the tumor’s center. In the original implementations of the algorithm, was fixed to be 0.192, giving a celldoubling time of 4 days. In the merged algorithm proposed here, we wanted to account for fact that tumor cells with a higher oxygen concentration likely have a larger probability of dividing than those with a lower oxygen concentration. For this reason, we have modified the algorithm so that depends on the distance to the closest blood vessel (which is proportional to the oxygen concentration at a given cell site). The average value of was fixed to be 0.192, and we have specified that takes on a minimum value of 0.1 and a maximum value of 0.284. This means that a proliferative cell in the model can have a cell doubling time anywhere in the range of three to seven days. The formula used to determine is
(21) where is the diffusion length of oxygen, taken to be 250 m [20, 18]. Both and depend on the average probability of division. If this average probability changes, so does and .

Tumor Center and Area: After each cell has evolved, recalculate the geometric center and area of the tumor.

7.2 Simulating Heterogeneous Tumor Growth
The 3D cancer simulation tool described here was employed to study tumor growth in a confined environment: a twodimensional representation of the cranium in space as a function of time [21]. The cranium is idealized as an elliptical growthpermitting environment with two growthprohibiting circular obstacles representing the ventricular cavities. Tumor growth is initiated in between a ventricular cavity and the cranium wall. In this setting, we find that the earlytime characteristics of the tumor and the vasculature are not significantly different than those observed when radial symmetry is imposed on tumor growth. In particular, after 45 days of growth (Figure 16(a)), vessels associated with the radially symmetric tumor begin to regress and hypoxia results in the tumor center. Twenty days later (Figure 16(b)), a strong, disordered angiogenic response has occurred in the still radially symmetric tumor. Over the next 50 days of growth (Figure 16(c) and (d)), the disorganized angiogenic blood vessel network continues to vascularized the growing tumor, but the tumor’s shape begins to deviate from that of a circle due to the presence of the confining boundary. The patterns of vascularization observed are consistent with the patterns observed in the original vascular model [18], suggesting that the merged algorithm maintains the functionality of the original vascular algorithm.
Figure 16: Tumor growing in a 2D representation of the cranium.
However, if the results of this simulation are compared with those of the environmentallyconstrained algorithm without the explicit incorporation of the vasculature [20], we find that the merged model responds to the environmental constraints in a way that is more physically intuitive. In the original environmentallyconstrained algorithm [20], the tumor responds quickly and drastically to the confining boundary and ventricular cavities. This occurs because the original evolution rules not only determine the probability of division based on the distance to the boundary, but also determine the state of a cell based on a measure of its distance to the boundary. In the merged model which explicitly incorporates the vasculature, the state of each cell depends on the blood vessel network, and only the probability of division directly depends on the boundary. For this reason, the merged algorithm exhibits an emergent property in that it grows tumors that respond more gradually and naturally to environmental constraints than does the algorithm without the vasculature.
The tumor growth in a twodimensional irregular region of space that truly allows the neoplasm to adapt its shape as it grows in time (i.e., a 3D model) was also investigated by Gevertz et al [20]; see also Gevertz and Torquato [21]. As with the twodimensional representation of the cranium in space, an emergent property of the merged algorithm in which that a more subtle and natural response to the effects of physical confinement is found occurs. The studies taking into account mutations responsible for phenotypic heterogeneity have been carried out by Gevertz and Torquato [21], to which the readers are referred for more details. We note that all the results presented in this section need to be validated experimentally.
8 Analysis of the Invasion Network: Minimal Spanning Trees
It is well known that cancer cells can break off the main tumor mass and invade healthy tissue. For many cancers, this process can eventually result in metastases to other organs. Tumorcell invasion is a hallmark of glioblastomas, as individual tumor cells have been observed to spread diffusely over long distances and can migrate into regions of the brain essential for the survival of the patient [27]. In certain cases, the invading tumor cells form branched chains (see Figure 3), i.e., tree structures [61]. The brain offers these invading cells a variety of pathways they can invade along (such as blood vessel and white fiber tracts), which may be interpreted as the edges of an underlying graph with the various “resistances” values along these pathways playing the role of edge weights. The underlying physics behind the formation of the observed patterns are only beginning to be understood.
Figure 17: Examples of weighted graph and the resulting minimal spanning tree.
The competition between local and global driving forces is a crucial factor is determining the structural organization in a wide variety of naturally occurring branched networks [64, 65, 66]. As an attempt toward a model of the invasive network emanating from a solid tumor, Kansal and Torquato [61] investigated the impact of a global minimization criterion versus a local one on the structure of spanning trees. Spanning trees are defined as a loopless, connected set of edges that connect all of the nodes in the underlying graph (see Figure 17). In particular, these authors considered the generalized minimal spanning tree (GMST) and generalized invasive spanning tree (GIST), because they generally offer extremes of global (GMST) and local (GIST) criteria. Both GMST and GIST are defined on graphs in which the nodes are partitioned into groups and each edge has an assigned weight. GMST is refined (relative to that of a spanning tree) such that the requirement that every node of the graph is included in the tree is replaced by the inclusion of at least one node from each group with the additional requirement that the total weight of tree is minimized [67]. GIST can be constructed by growing a connected cluster of edges by “invading” the remaining the edge with the minimal weight at its boundary with the requirement of the inclusion of at least one node from each group in the final tree [61].
Kansal and Torquato [61] have developed efficient algorithms to generate both GMST and GIST structures, as well as a method to convert GIST structure incrementally into a more globally optimized GMSTlike structure (see Figure 18). The readers are referred to the original paper for more algorithmic details. These methods allow various structural features to be observed as a function of the degree to which either criterion is imposed and the intermediate structures can then serve as benchmarks for comparison when a real image is analyzed.
Figure 18: Examples of GMST and GIST.
We note that a general procedure by which information extracted from a single, fixed network structure can be utilized to understand the physical processes which guided the formation of that structure is highly desirable in understanding the invasion network of tumor cells, since the temporal development of such a network is extremely difficult to observe. To this end, Kansal and Torquato [61] examined a variety of structural characterizations and found that the occupied edge density (i.e., the fraction of edges in the graph that are included in the tree) and the tortuosity of the arcs in the trees (i.e., the average of the ratio of the path length between two arbitrary nodes in the tree and the Euclidean distance between them) correlate well with the degree to which an intermediate structure resembles the GMST or GIST. Since both characterizations are straightforward to determine from an image (e.g., only the information of the tree is required for tortuosity and additional information of underlying graph is needed for occupied edge density), they are potentially useful tools in the analysis of the formation of invasion network structures. Once the distribution of the invasive cells in the brain is understood, a cellular automaton simulation tool for glioblastoma that is useful in a clinical setting could be developed. This of course would apply more generally to invasion networks around other solid tumors.
9 Conclusions and Future Work
In this paper, we have reviewed the work that we have performed to attempt to develop an Ising model of cancer. We began by describing a minimalist 4D cellular automaton model of cancer in which healthy cells transition between states (proliferative, hypoxic, and necrotic) according to simple local rules and their present states, which can viewed as a strippeddown Ising model of cancer [11, 12, 13]. Using four proliferation parameters, this model was shown to reflect the growth dynamics of a clinically untreated tumor very well [11]. This was followed by discussion of an extension of the model to study the effect on the tumor dynamics and geometry of a mutated subpopulation [12] and how tumor growth is affected by chemotherapeutic treatment, including induced resistance, with additional three treatment parameters [13]. An improved CA model that explicitly accounts for angiogenesis [18] as well as the heterogeneous and confined environment in which a tumor grows [20] were discussed. A general cancer simulation tool that merges, adapts and improves all of the aforementioned mechanism into a single CA model was also presented and applied to simulate the growth the GBM in a vascularized confined cranium [21]. Finally, we touched on how one might characterize the invasive network organization (local versus global organization) around a solid tumor using spanning trees [61]. However, we must move well beyond the improved CA model as well as other computational models of cancer in order to make real progress on controlling this dreaded set of diseases.
9.1 The Obvious but Necessary
Formulating theoretical and computational tools that can be utilized clinically to predict neoplastic progression and propose individualized optimal treatment strategies to control cancer growth is the holy grail of tumor modeling. Although the development of our most comprehensive cellular automaton model is potentially a useful step towards the longterm goal of an Ising model for cancer, numerous complex mechanisms involved in tumor growth and their interactions needs to be identified and understood in order to truly achieve this goal.
For example, an effective Ising model of cancer must incorporate molecularlevel information via a better understanding of the cellular origin of the tumor. Such information might become available if imaging techniques for spatial statistics of cell/molecular heterogeneity can be developed. This would enable an improved understanding of invading cancer cells: cell motility, cellcell communication and phenotypes of invading cells. Such knowledge is crucial in order to predict the effects of treatment and tumor recurrence. The incorporation of stem cells, oncogenes and tumor suppressor genes in computational models would aid in our understanding of tumor progression.
In addition, we must quantitatively characterize the biological (host) environment (i.e., a heterogeneous material/medium) in which cancer evolves, including both the microstructure and the associated physical properties. For example, a better knowledge of diffusion and transport of nutrients, drugs, etc. would significantly improve the accuracy of the model simulating the effects of vasculature evolution and treatment. Similarly, cell mechanics and mechanical stresses must be understood. In such cases, imaging of the biological environment over a wide spectrum of length and time scales will be crucial.
Figure 19: A cartoon of a twophase medium.
It is important to emphasize that the theory heterogeneous media is a huge field within the physical sciences that can be brought to bear to better understand the host heterogeneous microenvironment of cancer and metastases (see Figure 19). For example, there exist powerful and sophisticated theoretical/computational techniques to characterize the microstructure of heterogeneous materials and predict their physical properties [62]. Specifically, the details of complex microstructures are described in terms of various statistical descriptors (different types of correlation functions), which in turn determine the physical properties of the heterogeneous materials [62]. In particular, the effective properties that have been predicted include the diffusion coefficient [68], reaction rates [69, 70], elastic/viscoelastic moduli [71, 72], thermal conductivity [73], thermal expansion coefficient [74], fluid permeability [75], and electrical conductivity [76, 77]. Accurate characterizations of these properties of the host environment and tumor mass are essential in order to significantly improve models for tumor growth and invasion. For example, a knowledge of the elastic properties enables one to better model the effects of physical confinement and the mechanical response of solid tumor; while the diffusion coefficient and fluid permeability are crucial to model transport of nutrients and proteins, delivery of drugs and even the migration of cancer cells. These techniques have been used to propose a novel biologically constrained threephase model of the brain microstructure [78].
Given such information, the CA model can be modified accordingly to take into account the available cell/molecular details of the tumor mass, its invasion network and the host heterogeneity (e.g., the capillary vasculature and adaptive physical confinement). Realtime tumor growth and treatment simulations can be carried out to generate data of clinical utility. For example, instead of only producing data which qualitatively reflects the general effects of tumor treatment and resistance, one could use the model to make reliable prognosis and to optimize individual treatment strategy.
It would be fascinating to see if a more refined Ising model for cancer predicted a “phase transition” phenomenon, which would be in keeping with the behavior of the standard Ising model for spin systems. For example, it is not hard to imagine that part of the tumorigenesis process involves a “phase transition” between premalignant cells and malignant cells.
9.2 Not So Obvious: Optimization and Cancer
We also note that variational principles [62, 79] and optimization techniques [80, 81, 82, 83] have been fruitfully applied to design structures with optimal properties. Can optimization techniques be applied to understand and control cancer? Although optimization methods have begun to be employed for such purposes, there full potential has yet to be realized. For tumor treatment, for example, optimization techniques could be employed to design chemotherapy/radiation strategies depending on tumor size, genomic information and the heterogeneous environment as well as the optimal durations of treatment and rest periods. Given sufficient patientspecific information, optimized treatment strategies can be designed for individual patients. A variety of optimization techniques could be brought to bear here, including simulated annealing methods, and linear and nonlinear programming techniques.
Figure 20: Minimal surface structure
We have developed an optimization methodology that provides a means of optimally designing multifunctional composite microstructures [81, 83]. We have shown how the competition between two different performance demands (thermal versus electrical behaviors or electrical versus mechanical behaviors) results in unexpected microstructures, namely, minimal surfaces [81, 82] (see Figure 20), which also appear to optimal for fluid transport [84] as well as diffusioncontrolled reactions [85]. This work suggests that it may be fruitful to explore the development of cancer, which not only involves competition but cooperation, from a rigorous multifunctional optimization viewpoint. Cancer processes involve a competition between the primary clone, subclones, healthy tissue, immune system, etc. as well as a cooperation between different cells types (e.g., stroma cells and cancer cells) in a heterogeneous environment. This competition/cooperation can be translated into an optimization problem in space and time. Adaptation of this multifunctional optimization approach to cancer modeling could provide an alternative to gametheory approaches to understanding cancer [86].
9.3 The Far Out
Figure 21: Diamond and disordered ground state.
Even more challenging and intriguing questions can be asked: Can we exploit the unique properties of normal stem cells [87] to control cancer (e.g., to deliver therapy to tumors or to have them compete with the tumor)? Can we use inverse optimization methods to design “hypothetical” cancers or stem cells with particular the cellcell interactions to yield targeted behaviors and then make them in the lab? These “inverse” problems are motivated by their analog in statistical mechanics [88, 89, 90, 91]. In statistical mechanics, the “forward problem” is one in which a Hamiltonian (interaction potential) for a manybody system is specified and then structure of the system and its thermodynamics are predicted. By contrast, the “inverse” problem of statistical mechanics seeks to find the “optimal” interaction potential that leads spontaneously to a novel “targeted” structure (or behavior). We have discovered optimal interaction potentials that yield unusual or counterintuitive targeted ground (zerotemperature) states, e.g., lowcoordinated diamond crystal [89] and disordered states [90] with only isotropic pair potentials (see Figure 21). Ground states are those manyparticle configurations that arise as a result of slowly cooling a liquid to absolute zero temperature. The aforementioned obtained targeted ground states are so unusual because much of our experience involves ground states that are highlycoordinated crystal structures [91]. An extremely challenging and fascinating question is whether we can devise inverse optimization techniques to control cancer?
It is clear that theoretical methods based in the physical and mathematical sciences offer many different fruitful ways to contribute to tumor research. However, for this approach to be successful, intensive interactions with cell biologists, oncologists, radiologists, clinicians, physicists, chemists, engineers, and applied mathematicians are essential. Such an interdisciplinary approach appears to be necessary in order to control this deadly disease. This could be achieved most effectively if we could have an analog of the “Manhattan Project” in which there was a single facility with such an interdisciplinary team of scientists dedicated to this supreme achievement.
Acknowledgments
The author thanks Yang Jiao and Jana Gevertz for very useful discussions and their critical reading of this manuscript. The research described was supported by Award Number U54CA143803 from the National Cancer Institute. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Cancer Institute or the National Institutes of Health.
References
References
 [1] Hanahan D and Weinberg R A. The hallmarks of cancer. Cell, 100:57–70, 2000.
 [2] Crossa F R and Tinkelenberga A H. A potential positive feedback loop controlling cln1 and cln2 gene expression at the start of the yeast cell cycle. Cell, 65:875–883, 1991.
 [3] Deisboeck T S, Berens M E, Kansal A R, Torquato S, Rachamimov A, Louis D N, and Chiocca E A. Patterns of selforganization in tumor systems: Complex growth dynamics in a novel brain tumor spheroid model. Cell Proliferat., 34:115–134, 2001.
 [4] Brand U, Fletcher J C, Hobe M, Meyerowitz E M, and Simon1dagger R. Dependence of stem cell fate in arabidopsis on a feedback loop regulated by clv3 activity. Science, 289:617–619, 2000.
 [5] Kitano H. Cancer as a robust system: implications for anticancer therapy. Nature Reviews Cancer, 4:227–235, 2004.
 [6] Chaplain M A J and Sleeman B D. Modelling the growth of solid tumours and incorporating a method for their classification using nonlinear elasticity theory. J. Math. Biol., 31:431–473, 1993.
 [7] Tracqui P, Cruywagen G C, Woodward D E, Bartoo G T, Murray J D, and Alvord E C. A mathematical model of glioma growth: The effect of chemotherapy on spatiotemporal growth. Cell Proliferation, 28:17–31, 1995.
 [8] Gatenby R A. Applications of competition theory to tumour growth: Implications for tumour biology and treatment. Eur. J. Cancer, 32:722–726, 1996.
 [9] Anderson A R A and Chaplain M A J. Continuous and discrete mathematical models of tumorinduced angiogenesis. Bull. Math. Biol., 60:857â900, 1998.
 [10] Panetta J C. A mathematical model of drug resistance: Heterogeneous tumors. Mathematical Biosciences, 147:41–61, 1998.
 [11] Kansal A R, Torquato S, Harsh G R, Chiocca E A, and Deisboeck T S. Simulated brain tumor growth using a threedimensional cellular automaton. J. Theo. Bio., 203:367–382, 2000.
 [12] Kansal A R, Torquato S, Chiocca E A, and Deisboeck T S. Emergence of a subpopulation in a computational model of tumor growth. J. Theo. Bio., 207:431–441, 2000.
 [13] Schmitz J E amd Kansal A R and Torquato S. A cellular automaton model of brain tumor treatment and resistance. J. Theo. Medicine, 4:223–239, 2002.
 [14] McDougall S R, Anderson A R A, Chaplain, M A J, and Sherrat J A. Mathematical modelling of flow through vascular networks: implications for tumorinduced angiogenesis and chemotherapy strategies. Bull. Math. Biol., 64:673–702, 2002.
 [15] Alarcón T, Byrne H M, and Maini P K. A multiple scale model for tumor growth. Multiscale Model. Simul., 3:440–475, 2005.
 [16] Cristini V, Frieboes H B, Gatenby R, Caerta S, Ferrari M, and Sinek J. Morphological instability and cancer invasion. Clin. Cancer Res., 11:6772–6779, 2005.
 [17] Frieboes H B, Zheng X, Sun ChungHo, Tromberg B, Gatenby R, and Cristini V. An integrated computational/experimental model of tumor invasion. Cancer Res., 66:1597–1604, 2006.
 [18] Gevertz J L and Torquato S. Modeling the effects of vasculature evolution on early brain tumor growth. J. Theo. Bio., 243:517, 2006.
 [19] Swanson K R, Rockne R, Rockhill J K, and Alvord E C Jr. Combining mathematical modeling with serial mr imaging to quantify and predict response to radiation therapy in individual glioma patient. NeuroOncology, 9:575, 2007.
 [20] Gevertz J L, Gillies G, and Torquato S. Simulating tumor growth in confined heterogeneous environments. Physical Biology, 5:036010, 2008.
 [21] Gevertz J and Torquato S. Growing heterogeneous tumors in silico. Phys. Rev, E, 80:051910, 2009.
 [22] Gerlee P and Anderson A R. Evolution of cell motility in an individualbased model of tumour growth. J. Theor. Biol., 259:67–83, 2009.
 [23] RamisConde I, Chaplain M A, Anderson A R, and Drasdo D. Multiscale modelling of cancer cell intravasation: the role of cadherins in metastasis. Phys. Bio., 6:16008, 2009.
 [24] HB Frieboes, F Jin, YL Chuang, SM Wise, JS Lowengrub, and V Cristini. Threedimensional multispecies nonlinear tumor growth II: Tumor invasion and angiogenesis. J. Theor. Biol., 264:1254–78, 2010.
 [25] Byrne H. Dissecting cancer through mathematics: from the cell to the animal model. Nature Review, 10:221–231, 2010.
 [26] Coffey D S. Selforganization, complexity and chaos: the new biology for medicine. Nat. Med., 4:882–885, 1998.
 [27] Holland E C. Glioblastoma multiforme: the terminator. PNAS, 97:6242–6244, 2000.
 [28] C J Wheeler, K L Black, G Liu, M Mazer, XX Zhang, S Pepkowitz, D Goldfinger, H Ng, D Irvin, and J S Yu. Vaccination elicits correlated immune and clinical responses in i glioblastoma multiforme patients. Cancer Res., 68:5955–5964, 2010.
 [29] Holash J, Maisonpierre P C, Compton D, Boland P, Alexander C R, Zagzag D, Yancopoulos G D, and Weigand S J. Vessel cooption, regression, and growth in tumors mediated by angiopoietins and vegf. Science, 284:1994–1998, 1999.
 [30] Brat D J, CastellanoSanchez A A, Hunter S B, Pecot M, Cohen C, Hammond E H, Devi S N, Kaur B, and Van Meir E G. Pseudopalisades in glioblastoma are hypoxic, express extracellular matrix proteases, and are formed by an actively migrating cell popultion. Cancer Res., 64:920–927, 2004.
 [31] Singh S K, Clarke I D, Hide T, and Dirks P B. Cancer stem cells in nervous system tissues. Oncogene, 23:7267–7273, 2004.
 [32] Singh S K, Hawkins C, Clarke I D, Squire J A, Bayani J, Hide T amd Mekelman R M, Cusimano M D, and Dirks P B. Identification of human brain tumour initiating cells. Nature, 432:396–401, 2004.
 [33] Düchting W and Vogelsaenger T. Recent progress in modelling and simulation of threedimensional tumor growth and treatment. Biosystems, 18:79–91, 1985.
 [34] Qi A S, Zheng X, Du C Y, and An B S. A cellular automaton model of cancerous growth. J. Theor. Biol., 161:1–12, 1993.
 [35] Smolle J and Stettner H. Computer simulation of tumour cell invasion by a stochastic growth model. J. Theor. Biol., 160:63–72, 1993.
 [36] Wasserman R, Acharya R, Sibata C, and Shin K H. A patientspecific in vivo tumor model. Math. Biosci., 136:111–140, 1996.
 [37] Steel G G and Peckham M J. Exploitable mechanisms in combined radiotherapychemotherapy: The concept of additivity. Int. J. Rad. Oncol., 5:85–91, 1997.
 [38] Berkman R A, Clark W C, Saxena A, Robertson J T, Oldfield E H, and Ali I U. Clonal composition of glioblastoma multiforme. J. Neurosurg., 77:432–437, 1992.
 [39] Coons S W and Johnson P C. Regional heterogeneity in the DNA content of human gliomas. Cancer, 72:3052–3060, 1993.
 [40] Paulus W and Peiffer J. Intratumoral histologic heterogeneity of gliomas. a quantitative study. Cancer, 64:442–447, 1989.
 [41] Yung W A, Shapiro J R, and Shapiro W R. Heterogeneous chemosensitivities of subpopulations of human glioma cells in culture. Cancer Res., 42:992–998, 1982.
 [42] Nowell P C. The clonal evolution of tumor cell populations. Science, 194:23–28, 1976.
 [43] Loeb L A. Microsatellite instability: Marker of a mutator phenotype in cancer. Cancer Res., 54:5059–5063, 1994.
 [44] Lengauer C, Kinzler K W, and Vogelstein B. Genetic instabilities in human cancers. Nature, 396:643–649, 1998.
 [45] Heppner G H and Miller B E. Therapeutic implications of tumor heterogeneity. Bull. Mathemat. Biol., 16:91–105, 1989.
 [46] Schnipper L E. Clinical implications of tumor cell heterogeneity. New Eng. J. Med., 314:1423–1431, 1986.
 [47] Bredel M. Anticancer drug resistance in primary human brain tumors. Brain Research Reviews, 35:161–204, 2001.
 [48] Endicott L C and Ling V. The biochemistry of the pglycoproteinmediated multidrug resistance. Annual Reviews in Biochemistry, 58:137–171, 1989.
 [49] German U A. Pglycoprotein–a mediator of multidrug resistance in tumour cells. European Journal of Cancer, 32:927–944, 1996.
 [50] Poppenborg H, Munstermann G, Knopfer M M, Hotfilder M, and Wolff J E A. C6 crossresistant to cisplatin and radiation. Anticancer Research, 17:2073–2077, 1997.
 [51] Chadhary P M and Roninson I B. Induction of multidrug resistance in human cells by transient exposure to different chemotherapeutic drugs. Journal of the National Cancer Institute, 85:632–639, 1993.
 [52] Gekeler V, Beck J, Noller A, Wilisch A, Frese G, Neumann M, Handgretinger R, Ehninger G, Probst H, and Niethammer D. Druginduced changes in the expression of mdrassociated genes  investigations on cultured cell lines and chemotherapeutically treated leukemias. Annals of Hematology, 69:S19–S24, 1994.
 [53] Gerweck L E, Kornblinth P L, Burlett P, Wang J, and Sweigert S. Radiation sensitivity of cultured glioblastoma cells. Radiology, 125:231–241, 1977.
 [54] Kayama T, Yoshimoto T, Fujimoto S, and Sakurai Y. Intratumoral oxygen pressure in malignant brain tumors. Journal of Neurosurgury, 74:55–59, 1991.
 [55] Zhang W, Hara A, Sakai N, Andoh T, Yamada H, anf P H. Gutin Y N, and Kornblinth P L. Radiosensitization and inhibition of deoxyribonucleic acid repair in rat glioma cells by long term treatment with 12otetradecanoylphorbol 13acatate. Neurosurgury, 32:432–437, 1993.
 [56] Holash J, Wiegand S J, and Yancopoulos G D. New model of tumor angiogenesis: dynamic balance between vessel regression and growth mediated by angiopoietins and vegf. Oncogene, 18:5356–5362, 1999.
 [57] Maisonpierre C S, Suri C, Jones P F, Bartunkova S, Wiegand S J, Radziejewski C, Compton D, McClain J, Aldrich T H, Papadlopoulous N, Daly T J, Davis S, Sato T N, and Yancopoulos G D. Angiopoietin2, a natural antagonist for tie2 that disrupts in vivo angiogenesis. Science, 277:55–60, 1997.
 [58] Secomb T W, Hsu R, Beamer N B, and Coull B M. Theoretical simulation of oxygen transport to brain by networks of microvessels: effects on oxygen supply and demand on tissue hypoxia. Microcirculation, 7:237–247, 2000.
 [59] Giese A and Manfred W. Glioma invasion in the central nervous system. Neurosurgery, 39:235–252, 1996.
 [60] Visted T, Enger P O, LundJohansen M, and Bjerkvig R. Mechanisms of tumor cell invasion and angiogenesis in the central nervous system. Front. in Bioscience, 8:289–304, 2003.
 [61] Kansal A R and Torquato S. Globally and locally minimal weight branched tree networks. Physica A, 301:601–619, 2001.
 [62] Torquato S. Random Heterogeneous Materials: Microstructure and Macroscopic Properties. SpringerVerlag, New York, 2002.
 [63] Helmlinger G, Netti P A, Lichtenbeld H C, Melder R J, and Jain R K. Solid stress inhibits the growth of multicellular tumor spheroids. Nat. Biotechnol., 15:778–83, 1997.
 [64] T. Sun, P. Meakin, and T. Jossang. Minimum energy dissipation model for river basin geometry. Phys. Rev. E, 49:4865â4872, 1994.
 [65] I. RodriguezâIturbe and A. Rinaldo. Fractal River Basins. Cambridge University Press, Cambridge, 1997.
 [66] G. B. West, J.H. Brown, and B. J. Enquist. A general model for the origin of allometric scaling laws in biology. Science, 276:122â126, 1997.
 [67] Haouari M Dror M and Chaouachi J. Generalized spanning trees. Eur. J. Oper. Res., 120:583–592, 2000.
 [68] Torquato S, Kim I C, and Cule D. Effective conductivity, dielectric constant, and diffusion coefficient of digitized composite media via firstpassagetimeequations. J. Appl. Phys., 85:1560–1572, 1999.
 [69] S B Lee, I C Kim, C A Miller, and S Torquato. Randomwalk simulation of diffusioncontrolled processes among static traps. Phys. Rev. B, 39:11833–11839, 1989.
 [70] Torquato S. Diffusion and reaction among traps: some theoretical and simulation results. J Stat. Phys., 65:1173–1207, 1991.
 [71] Quintanilla J and Torquato S. New bounds on the elastic moduli of suspensions of spheres. J. Appl. Phys., 77:4361–4372, 1995.
 [72] Torquato S. Exact epression for the effective elastic tensor of disordered composites. Phys. Rev. Lett., 79:681–685, 1997.
 [73] Kim I C and S Torquato. Effective conductivity of suspensions of hard spheres by Brownian motion simulation. J. Appl. Phys., 69:2280–2289, 1991.
 [74] Gibiansky L V and Torquato S. Thermal expansion of isotropic multiphase composites and polycrystals. J. Mech. Phys. of Solids, 45:1223–1252, 1997.
 [75] Torquato S and Lu B. Rigorous bounds on the fluid permeability: effect of polydispersivity in grain size. Phys. Fluids A, 2:487–491, 1990.
 [76] Torquato S. Effective electrical conductivity of twophase disordered composite media. J. Appl. Phys., 58:3790–3797, 1985.
 [77] Sen A K and Torquato S. Effective electrical conductivity of twophase disordered anisotropic composite media. Phys. Rev. B, 39:4504–4516, 1989.
 [78] J. L. Gevertz and S. Torquato. A novel threephase model of brain tissue microstructure. Plos. Comput. Bio., 4:e100052, 2008.
 [79] Torquato S and Pham D C. Optimal bounds on the trapping constant and permeability of porous media. Phys. Rev. Lett., 92:255505, 2004.
 [80] Sigmund O and Torquato S. Design of materials with extreme thermal expansion using a threephase topology optimization method. J. Mech. Phys. Solids, 45:1037–1067, 1997.
 [81] Torquato S, Hyun S, and Donev A. Multifunctional composites: optimizing microstructures for simultaneous transport of heat and electricity. Phys. Rev. Lett., 89:266601, 2002.
 [82] S. Torquato and A. Donev. Minimal surfaces and multifunctionality. Proc. R. Soc. Lond. A, 460:1849–1856, 2004.
 [83] Torquato S. Optimal design of heterogeneous materials. Ann. Rev. Mater. Res., 40:101–129, 2010.
 [84] Y. Jung and S. Torquato. Fluid permeabilities of triply periodic minimal surfaces. Phys. Rev. E, 92:255505, 2005.
 [85] J. L. Gevertz and S. Torquato. Mean survivial times of absorbing triply periodic minimal surfaces. Phys. Rev. E, 80:011102, 2009.
 [86] D Dingli, F A C C Chalub, F C Santos, S Van Segbroeck, and J M Pacheco. Cancer phenotype as the outcome of an evolutionary game between normal and malignant cells. Brit. J. Cancer, 101:1130–1136, 2009.
 [87] Reya T, Morrison S J, Clarke M F, and Weissman1 I L. Stem cells, cancer, and cancer stem cells. Nature, 414:105–111, 2001.
 [88] Rechtsman M, Stillinger F H, and Torquato S. Optimized interactions for targeted selfassembly: application to honeycomb lattice. Phys. Rev. Lett., 95:228301, 2005.
 [89] Rechtsman M, Stillinger F H, and Torquato S. Synthetic diamond and wurtzite structures selfassemble with isotropic pair interactions. Phys. Rev. E, 75:031403, 2007.
 [90] Batten R D, Stillinger F H, and Torquato S. Novel lowtemperature behavior in classical manyparticle systems. Phys. Rev. Lett., 103:050602, 2009.
 [91] Torquato S. Inverse optimization techniques for targeted selfassembly. Soft Matter, 5:1157–1174, 2009.
Glossary

Neoplasm: A neoplasm is a synonym for a tumor.

Glioma: A collection of tumors arising from the glial cells or their precursors in the central nervous system.

Cellular automaton: A spatially and temporally discrete model that consists of a grid of cells, with each cell being in one of a number of predefined states. The state of a cell at a given point in time depends on the state of itself and its neighbors at the previous discrete time point. Transitions between states are determined by a set of local rules.

Ising model: The Ising model is an idealized statisticalmechanical model of ferromagnetism that is based on simple localinteraction rules, but nonetheless lead to basic insights and features of real magnets, such as phase transitions with a critical point.

Voronoi cell: Given a set of points, the Voronoi cell is the cell that is formed about an arbitrary point in the set by finding the region of space closer to that point than any other point in the system (Torquato, 2002).

Delaunay triangulation: Given a Voronoi graph (a set of Voronoi cells), the Delaunay graph is its dual that results from joining all pairs of sites that share a Voronoi face. If this graph consists of only simplices, the graph is called a Delaunay triangulation (Torquato, 2002).

Quiescent: A cell is considered quiescent if it is in the G0 phase of the cell cycle and is not actively dividing.

Necrotic: A cell is considered necrotic if it has died due to injury or disease, such as abnormally low oxygen levels.