Emergent Behaviors from A Cellular Automaton Model for Invasive Tumor Growth in Heterogeneous Microenvironments
Yang Jiao, Salvatore Torquato
1 Physical Science in Oncology Center, Princeton Institute for the Science and Technology of Materials, Princeton University, Princeton New Jersey 08544, USA
2 Department of Chemistry and Physics, Princeton Center for Theoretical Science, Program in Applied and Computational Mathematics, Princeton University, Princeton New Jersey 08544, USA
Email: yjiao@princeton.edu Email: torquato@electron.princeton.edu
Abstract
Understanding tumor invasion and metastasis is of crucial importance for both fundamental cancer research and clinical practice. In vitro experiments have established that the invasive growth of malignant tumors is characterized by the dendritic invasive branches composed of chains of tumor cells emanating from the primary tumor mass. The preponderance of previous tumor simulations focused on noninvasive (or proliferative) growth. The formation of the invasive cell chains and their interactions with the primary tumor mass and host microenvironment are not well understood. Here, we present a novel cellular automaton (CA) model that enables one to efficiently simulate invasive tumor growth in a heterogeneous host microenvironment. By taking into account a variety of microscopicscale tumorhost interactions, including the shortrange mechanical interactions between tumor cells and tumor stroma, degradation of extracellular matrix by the invasive cells and oxygen/nutrient gradient driven cell motions, our CA model predicts a rich spectrum of growth dynamics and emergent behaviors of invasive tumors. Besides robustly reproducing the salient features of dendritic invasive growth, such as least resistance and intrabranch homotype attraction, we also predict nontrivial coupling of the growth dynamics of the primary tumor mass and the invasive cells. In addition, we show that the properties of the host microenvironment can significantly affect tumor morphology and growth dynamics, emphasizing the importance of understanding the tumorhost interaction. The capability of our CA model suggests that welldeveloped in silico tools could eventually be utilized in clinical situations to predict neoplastic progression and propose individualized optimal treatment strategies.
Author Summary
The goal of the present work is to develop an efficient singlecell based cellular automaton (CA) model that enables one to investigate the growth dynamics and morphology of invasive solid tumors. Recent experiments have shown that highly malignant tumors develop dendritic branches composed of tumor cells that follow each other, which massively invade into the host microenvironment and ultimately lead to cancer metastasis. Previous theoretical/computational cancer modeling neither addressed the question of how such chainlike invasive branches form nor how they interact with the host microenvironment and the primary tumor. Our CA model, which incorporates a variety of microscopicscale tumorhost interactions (e.g., the mechanical interactions between tumor cells and tumor stroma, degradation of extracellular matrix by the tumor cells and oxygen/nutrient gradient driven cell motions), can robustly reproduce experimentally observed invasive tumor evolution and predict a wide spectrum of invasive tumor growth dynamics and emergent behaviors in various different heterogeneous environments. Further refinement of our CA model could eventually lead to the development of a powerful simulation tool that can be utilized in the clinic to predict neoplastic progression and propose individualized optimal treatment strategies.
Introduction
Cancer is not a single disease, but rather a highly complex and heterogeneous set of diseases that can adapt in an opportunistic manner, even under a variety of stresses. It is now well accepted that genome level changes in cells, resulting in the gain of function of oncoproteins or the loss of function of tumor suppressor proteins, initiate the transformation of normal cells into malignant ones and neoplastic progression [1, 2]. In the most aggressive form, malignant cells can leave the primary tumor, invade into surrounding tissues, find their way into the circulatory system (through vascular network) and be deposited at certain organs in the body, leading to the development of secondary tumors (i.e., metastases) [3].
The emergence of invasive behavior in cancer is fatal. For example, the malignant cells that invade into the surrounding host tissues can quickly adapt to various environmental stresses and develop resistance to therapies. The invasive cells that are left behind after resection are responsible for tumor recurrence and thus an ultimately fatal outcome. Therefore, significant effort has been expended to understand the mechanisms evolved in the invasive growth of malignant tumors [2, 4, 5, 6, 7] and their treatment [8, 9]. It is generally accepted that the invasive behavior of cancer is the outcome of many complex interactions occurring between the tumor cells, and between a tumor and the host microenvironment [3]. Tumor invasion itself is a complex multistep process involving homotype detachment, enzymatic matrix degradation, integrinmediated heterotype adhesion, as well as active, directed and random motility [4]. In recent in vitro experiments involving glioblastoma multiforme (GBM), the most malignant brain cancer, it has been observed that dendritic invading branches composed of chains of tumor cells are emanating from the primary tumor mass; see Fig. 1. Such invasive behaviors are characterized by intrabranch homotype attraction and least resistance [4].
Although recently progress has been made in understanding certain aspects of the complex tumorhost interactions that may be responsible for invasive cancer behaviors [4, 10, 11, 12], many mechanisms are either not fully understood or are unknown at the moment. Even if all of the mechanisms for cancer invasion could be identified, it is still not clear that progress in understanding neoplastic progression and proposing individualized optimal treatment strategies could be made without the knowledge of how these different mechanisms couple to one another and to the heterogeneous host microenvironment in which tumor grows [13]. Theoretical/computational cancer modeling that integrates apart mechanisms, when appropriately linked with experimental and clinical data, offers a promising avenue for a better understanding of tumor growth, invasion and metastasis. A successful model would enable one to broaden the conclusions drawn from existing medical data, suggest new experiments, test hypotheses, predict behavior in experimentally unobservable situations, and be employed for early detection and prognosis.
Indeed, cancer modeling has been a very active area of research for the last two decades (see Refs. [14] and [13] for recent reviews). A variety of interactions between tumor cells and between tumor and its host microenvironment have been investigated [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 29, 30, 31, 32], via continuum [25, 26, 27, 28, 32], discrete [16, 20, 33] or hybrid [19, 21, 22, 23] mathematical models. Very recently, multiscale mathematical models [22, 23, 28] have been employed to study the effects of the host microenvironment on the morphology and phenotypic evolution of invasive tumors and it has been shown that microenvironmental heterogeneity can dramatically affect the growth dynamics of invasive tumors. Although the simulated tumors showed certain invasive characteristics (e.g., developing protruding surfaces), no dendritic invasive branches emerged.
In response to the challenge to develop an “Ising” model for cancer growth [13], we generalize a cellbased discrete cellular automaton (CA) model that we have developed [16, 17, 18, 20, 21] to investigate the invasive growth of malignant tumors in heterogeneous host microenvironments. To the best of our knowledge, this generalized CA model is the first to investigate the formation of invasive cell chains and their interactions with the primary tumor mass and the host microenvironment. Our cellular automaton model takes into account a variety of microscopicscale tumorhost interactions, including the shortrange mechanical interactions between tumor cells and tumor stroma, the degradation of extracellular matrix by the invasive cells and oxygen/nutrient gradient driven cell motions and thus, it can predict a wide range of growth dynamics and emergent behaviors of invasive tumors. In particular, our CA model robustly reproduces the salient features of dendritic invasive growth observed in experiments, which is characterized by least resistance and intrabranch homotype attraction. The model also predicts nontrivial coupling of the growth dynamics of the primary tumor mass and the invasive cells, e.g., the invasive cells can facilitate the growth of primary tumor in harsh microenvironment. Moreover, we show that the properties of the host microenvironment can significantly affect tumor growth dynamics and lead to a variety of tumor morphology. These emergent behaviors naturally arise due to various microscopicscale tumorhost interactions, which emphasizes the importance of taking into account microenvironmental heterogeneity in understanding cancer. Further refinement of our model could eventually lead to the development of a powerful in silico tool that can be utilized in the clinic. As a demonstration of the capability and versatility of our CA model, we mainly consider invasive tumor growth in two dimensions, although the model is easily extended to three dimensions. Indeed, the algorithmic details of the model are given for any spatial dimension.
Materials and Methods
Biophysical Background of the CA Model
Voronoi Tessellation: The Underlying Cellular Structure
The underlying cellular structure is modeled using a Voronoi tessellation of the space into polyhedra [34], based on centers of spheres in a packing generated by a random sequential addition (RSA) process [16, 21] (see Fig. 2). In particular, nonoverlapping spheres are randomly and sequentially placed in a prescribed region until there is no void space for additional spheres, i.e., saturation is achieved. Such a saturated RSA packing possesses relative small variations in its Voronoi polyhedra and thus, has served as models for many biological systems [35, 36]. We refer to the polyhedra associated with the Voronoi tessellation as automaton cells. These automaton cells may correspond to real cells or tumor stroma (e.g., clusters of the ECM macromolecules). In previous studies, such automaton cells have represented clusters of real cells of various sizes or have implicitly represented healthy tissues [16]. Thus, the Voronoi tessellation associated with RSA sphere centers provides a highly flexible model for realcell aggregates with a high degree of shape isotropy. For example, one can use a variable automaton cell size to simulate avascular tumor growth from a few malignant cells to its macroscopic size [16]. In addition, such a Voronoi tessellation can reduce the undesired growth bias due to the anisotropy of ordered tessellations based on square and simple cubic lattices.
Since our new CA model explicitly takes into account the interactions between a single cell and its neighbors and microenvironment, each automaton cell here represents either a single tumor cell or a region of tumor stroma. Thus, the linear size of a single automaton cell is approximately m and the linear size of the 2D simulation domain is approximately 5 mm, which contains automaton cells. In the current model, we mainly focus on the effects of ECM macromolecule density and ECM degradation by malignant cells on tumor growth. Henceforth, we will refer to the host microenvironment (or tumor stroma) as “ECM” for simplicity. Each ECM associated automaton cell is assigned a particular density , representing the density of the ECM molecules within the automaton cell. A tumor cell can occupy an ECM associated automaton cell only if the density of this automaton cell , which means that either the ECM is degraded or it is deformed (pushed away) by the proliferating tumor cells.
Microenvironment Heterogeneity
The microenvironment in which tumor grows is usually highly heterogeneous, composed of various types of stromal cells and ECM structures. The ECM is a complex mixture of macromolecules that provide mechanical supports for the tissue (such as collagens) and those that play an important role for cell adhesion and motility (such as laminin and fibronectin) [22, 37, 38]. For different individuals with tumors, the ECM in the host microenvironments generally possess distinct mechanical and transport properties. By explicitly representing the ECM using automaton cells with different macromolecule densities, the effects of microenvironment heterogeneity on tumor growth can be very well explored. For example, various distributions of the ECM densities (i.e., the densities of the ECM macromolecules) can be employed to mimic the actual heterogeneous host microenvironment of the tumor. Certain tumor stroma contains fibroblasts, which actively produce ECM macromolecules leading to a higher ECM density around these cells. The automaton cells representing the ECM with larger densities are considered to be more rigid and more difficult to degrade. Since each automaton cell associated with the ECM has its own density, this allows a variation of ECM characteristics on the length scale comparable to that of a single tumor cell.
In addition, the tumor in our model is only allowed to grow in a compact growthpermitting region. This is to mimic the physical confinement of the host microenvironment, such as the boundary of an organ. In other words, only automaton cells within this region can be occupied by the cells of the tumor as it grows. In general, the growthpermitting region can be of any shape that best models the organ shape. Here we simply choose a spherical region to study the effects of the heterogeneous ECM on tumor growth. More sophisticated growthregion shapes have been employed to investigate the effects of physical confinement on tumor growth [20, 21]. Furthermore, we assume a constant radially symmetric nutrient/oxygen gradient in the growthpermitting region with the highest nutrition concentration at the boundary of this region, i.e., it is a vascular boundary. However, this assumption can also be relaxed.
Tumor Cell Phenotypes and Interactions with the Host Microenvironment
For highly malignant tumors, we consider the cells to be of one of the two classes of phenotypes: either invasive or noninvasive. Following Ref.[16], the noninvasive cells remain in the primary tumor and can be proliferative, quiescent or necrotic, depending on the nutrition supply they get. For avascular tumor growth, as we focus here, the nutrition the tumor cells can get are essentially those diffuse into the tumor through its surface. As the tumor grows, the amount of nutrition supply, which is proportional to the surface area of the tumor, cannot meet the needs of all cells whose number increases with the tumor volume, leading to the development of necrotic and quiescent regions. Following Ref.[16], characteristic diffusion lengths are employed to determine the states of a noninvasive cell. For example, quiescent cells more than away from the tumor surface become necrotic (see details in the next section). The diffusion length (also the characteristic thickness of the rim of living tumor cells) depends on the size of the primary tumor.
As a proliferative cell divides, its daughter cell effectively pushes away/degrades the surrounding ECM and occupies the automaton cell originally associated with the ECM [39, 40, 41]. It is easier for a tumor cell to take up an ECM associated automaton cell with lower density (i.e., less rigid ECM regions) than that with higher density (i.e., more rigid ECM regions) and thus, the tumor growth is affected by the ECM heterogeneity through the local mechanical interaction between tumor cells and the ECM. If there is no space available for the placement of a daughter cell within a distance from the proliferative cell, the proliferative cell turns quiescent.
The invasive cells are considered to be mutant daughters of the proliferative cells [42], which can gain a variety of degrees of ECM degradation ability (i.e., the matrixdegradative enzymes) and motility that allow them to leave the primary tumor and invade into surrounding microenvironment [43]. We consider the invasive cells can move from one automaton cell to another only if the ECM in the target automaton cell is completely degraded (i.e., with ). Each trial move of an invasive cell involves the degradation of the ECM in its neighbor automaton cells, followed by a possible move to one of the automaton cells whose ECM is completely degraded; otherwise the invasive cell does not move. The number of trial moves of an invasive cell and to what extent it degrades the ECM are respectively determined by and (see the following section for details). The oxygen/nutrient gradient also drives the invasive cells to move as far as possible from the primary tumor [44], which takes up the majority of oxygen/nutrients. The motility is the maximum possible number of such trial moves. In addition, we assume that the invasive cells do not divide as they migrate.
Algorithmic Details
We now provide specific details for the CA model to study invasive tumor growth in confined heterogeneous microenvironment. In what follows, we will simply refer to the primary tumor as “the tumor” and explicitly use “invasive” when considering invasive cells. After generating the automaton cells by Voronoi tessellation of RSA sphere centers, an ECM macromolecule density is assigned to each automaton cell within the growthpermitting region, which represents the heterogeneous host microenvironment. Then a tumor is introduced by designating any one or more of the automaton cells as proliferative cancer cells. Time is then discretized into units that represent one real day. At each time step:

Each automaton cell is checked for type: invasive, proliferative, quiescent, necrotic or ECM associated. Invasive cells degrade and migrate into the ECM surrounding the tumor. Proliferative cells are actively dividing cancer cells, quiescent cancer cells are those that are alive, but do not have enough oxygen and nutrients to support cellular division and necrotic cells are dead cancer cells.

All ECM associated automaton cells and tumorous necrotic cells are inert (i.e., they do not change type).

Quiescent cells more than a certain distance from the tumor’s edge are turned necrotic. The tumor’s edge, which is assumed to be the source of oxygen and nutrients, consists of all ECM associated automaton cells that border the neoplasm. The critical distance for quiescent cells to turn necrotic is computed as follows:
where is prescribed parameter (see Table 1), is the spatial dimension and is the distance between the geometric centroid of the tumor and the tumor edge cell that is closest to the quiescent cell under consideration. The position of the tumor centroid is given by
where is the total number of noninvasive cells contained in the tumor, which is updated when a new noninvasive daughter cell is added to the tumor.

Each proliferative cell will attempt to divide with probability into the surrounding ECM (i.e., the automaton cells associated with the ECM) by degrading and pushing away the ECM in that automaton cell. We consider that depends on both the physical confinement imposed by the boundary of the growthpermitting region and the local mechanical interaction between the tumor cells and the ECM, i.e.,
where is the base probability of division (see Table 1), is the distance of the dividing cell from the tumor centroid, is the distance between the closest growthpermitting boundary cell in the direction of tumor growth and the tumor’s geometric centroid and is the ECM density of the automaton cell to be taken by the new tumor cell. When a ECM associated automaton cell is taken by a tumor cell, it density is set to be zero. The predefined growth distance () is described in a following bullet point.

If a proliferative cell divides, it can produce a mutant daughter cell possessing an invasive phenotype with a prescribed probability (i.e., the mutation rate). The invasive daughter cell gains ECM degradation ability and motility , which enables it to leave the tumor and invade into surrounding ECM. The rules for updating invasive cells are given in a following bullet point. If the daughter cell is noninvasive, it is designated as a new proliferative cell.

A proliferative cell turns quiescent if there is no space available for the placement of a daughter cell within a distance from the proliferative cell, which is given by
where a nutritional parameter (see Table 1), is the spatial dimension and is the distance between the geometric tumor centroid and the tumor edge cell that is closest to the proliferative cell under consideration.

An invasive cell degrades the surrounding ECM (i.e., those in the neighboring automaton cells of the invasive cell) and can move from one automaton cell to another if the associated ECM is completely degraded. For an invasive cell with motility and ECM degradation ability , it will make attempts to degrade the ECM in the neighboring automaton cells and jump to these automaton cells, where is an arbitrary integer in . For each attempt, the surrounding ECM density is decreased by , where is an arbitrary number in . Using random numbers for ECM degradation ability and cellular motility is to take into account tumor genome heterogeneity, which is manifested as heterogeneous phenotypes (such as different and ). When the ECM in multiple neighboring automaton cells of the invasive cell are completely degraded (i.e., ), the invasive cell moves in a direction that maximizes the nutrients and oxygen supply. Here we assume that the migrating invasive cells do not divide. The degraded ECM shows the invasive path of the tumor.
The aforementioned automaton rules are briefly illustrated in Fig. 3. We note that noninvasive tumor growth can be studied by imposing a mutation rate . This enables us to compare the growth dynamics of invasive and noninvasive tumors and in turn to investigate the effects of coupling growth of the primary tumor mass and the invasive cells. Although we only consider sphericalgrowthpermitting regions here, the CA rules given above allow growthpermitting regions with arbitrary shapes. The important parameters mentioned in the bullet points above are summarized in Table 1. In the following, we will employ our CA model to investigate the growth dynamics of malignant tumors with different degrees of invasiveness in variety of different heterogeneous microenvironments.
Quantitative Metrics for Tumor Morphology
To characterize quantitatively the morphology of simulated tumors, we present several scalar metrics that capture the salient geometric features of the primary tumor, dendritic invasive branches or the entire invasive pattern. These metrics include the ratio of invasive area over tumor area (defined below), the specific surface of the invasive pattern, the asphericity of the primary tumor and the angular anisotropy metric for the invasive branches. The metrics are computed for all simulated tumors and compared to available experimental data. We note that the invasive pattern associated with a neoplasm includes both the primary tumor and the invasive branches.
Following Ref. [4], the tumor area is defined as the area of the circumcircle of the primary tumor (see Fig. 4(a)) and the invasive area is the area of the region between the effective circumcircle of the invasive pattern and the circumcircle of the primary tumor (see Fig. 4(a)). The radius of the effective circumcircle of the invasive pattern is defined to be the average distance from the invasive branch tip to the tumor center. The ratio as a function of time reflects the degree of coupling between the primary tumor and the invasive cells. If is linear in , there is no coupling; otherwise the two are coupled.
The specific surface [34] for the invasive pattern is defined as the ratio of the total length of the perimeter of the invasive pattern over its total area. In general, is inversely proportional to the size of the tumor and thus, large tumors have small values. Moreover, given the tumor size, tumors with a large number of long dendritic invasive branches possess a large value of . And is minimized for perfectly circular tumors with , where is the radius. Since depends on the size of the tumor, which makes it difficult to compare tumors with different sizes, in the calculations that follow we employ a normalized with respect to for an arbitraryshaped tumor with effective radius (i.e., the average distance from tumor edge to tumor center). For simplicity, we will still refer to the normalized specific surface as “specific surface” and designate it with symbol .
The asphericity of the primary tumor is defined as the ratio of the radius of circumcircle of the primary tumor over its incircle radius [45], i.e., (see Fig. 4(b)). A large value indicates a large deviation of the shape of primary tumor from that of a perfect circle, i.e., the tumor is more anisotropic.
To quantify the degree of anisotropy of the invasive branches, we introduce the angular anisotropy metric . In particular, the entire invasive pattern is evenly divided into sectors with lines emanating from the tumor center (see Fig. 4(c)). The angular anisotropy metric is defined as
(1) 
where is the average length of the invasive branches within the th sector and
(2) 
is the average length of all invasive branches. For tumors with invasive branches of similar lengths that are uniformly angularly distributed, the metric is small. Large fluctuations of both invasive branch length and angular distribution can lead to large values. In the following, we use to compute for the simulated invasive tumors.
Results
Model Validation
To verify the robustness and predictive capacity of our CA model, we first employ it to quantitatively reproduce the observed invasive growth of a GBM multicellular tumor spheroid (MTS) in vitro [4]. In particular, the boundary of the growthpermitting region is considered to be vascularized, i.e., a growing tumor can receive oxygen and nutrients from the growthpermitting region. A constant radially symmetric nutrient/oxygen gradient in the growthpermitting region with the highest nutrient/oxygen concentration at the vascular boundary is used. Initially, approximately 250 proliferative tumor cells are introduced at the center of the growthpermitting region with homogeneous ECM and tumor growth is started. This corresponds to an initial MTS with diameter m which is consistent with the in vitro experiment setup [4]. The following values of the growth and invasiveness parameters are used: , mm, mm, , , . Note that the value of corresponds to a cell doubling time of 40 hours, which is consistent with the reported experimental data [4]. A small value of ECM density is used, which corresponds to the soft DMEM medium used in the experiment [4]. In the visualizations of the tumor that follow, we use the following convention: the ECM in the growthpermitting region is white, and gray outside this region. The ECM degraded by the tumor cells is blue. In the primary tumor, necrotic cells are black, quiescent cells are yellow and proliferative cells are red. The invasive tumor cells are green.
Figure 5(a) and (b) respectively show the morphology of simulated MTS and a magnification of its invasive branches with increasing branch width towards the proliferative core. Specifically, one can clearly see that within the branches, chains of cells are formed as observed in experiments [4] (see Fig. 1). The invasive cells tend to follow one another (which is termed “homotype attraction”) since paths of degraded ECM are formed by pioneering invasive cells and it is easier for other cells to follow and enhance such paths than degrading ECM to create new paths by themselves. In other words, invasive cells tend to take paths with “least resistance”. We note that no CA rules are imposed to force such cellular behaviors. Instead, they are emergent properties that arise from our simulations.
The ratio of invasion area over primary tumor area as a function of time for the simulated tumor is computed and compared to reported experimental data [4] (see Fig. 5(c)). One can clearly see that our simulation results agree with experimental data very well. Moreover, the deviation of from a linear function of indicates that the growth of primary tumor and the invasive branches are strongly coupled [4] Other metrics for tumor morphology such as the specific surface of the invasive pattern, the sphericity of the primary tumor and the angular anisotropy metric for the invasive branches are computed from our simulation results and from the image of invasive MTS in Fig. 1(a) at 24 hours after initialization. The values are given in Table 2, from which one can see again a good agreement. Thus, we have shown that our CA model is both robust and quantitatively accurate with properly selected parameters.
Simulated Invasive Growth in Heterogeneous Miroenvironments
Having verified the robustness and predictive capacity of our CA model, we now consider three types of ECM density distributions, i.e., homogeneous, random and sinelike, to systematically study the effects of microenvironment heterogeneity on invasive tumor growth (see Fig. 6). These ECM density distributions represent real host microenvironments in which a tumor grows. (Details about these ECM distributions are given in the following sections.) Again, the boundary of the growthpermitting region is considered to be vascularized with a constant radially symmetric nutrient/oxygen gradient in the growthpermitting region pointing to the tumor center. We note that although generally the nutrient/oxygen concentration field in vivo is more complicated, previously numerical studies that considered the exact evolution of nutrient/oxygen concentrations have shown a decay of the concentrations toward the tumor center [22, 23]. Since the directions of cell motions are determined only by the nutrient/oxygen gradient, our constantgradient approximation is a very reasonable one.
In the beginning, a proliferative tumor cell is introduced at the center of the growthpermitting region and tumor growth is initiated. The growth parameters for the primary tumors in all cases studied here are the same and are given in Table 1. The invasiveness parameters and ECM densities are variable and specified in each case separately. The values of the growth parameters for the CA model were chosen to be consistent with GBM data from the medical literature [16, 20]. Specifically, the value of the base probability of division is , which corresponds to a cell doubling time of 4 days [46, 47]. This value is used for all of the cases of invasive growth that follow. Since our CA model takes into account general microscopic tumorhost interactions, we expect that the general growth dynamics and emergent behaviors predicted by the model will qualitatively apply to other solid tumors. We note that all of the reported growth dynamics and emergent properties of the simulated tumors for any specific set of growth and invasiveness parameters are repeatedly observed in 25 independent simulations.
Effects of Cellular Motility
We first simulate the growth of malignant tumors with different degrees of invasiveness in a homogeneous ECM with . In particular, we consider three invasive cases with the same mutation rate and ECM degradation ability , but different cell motility . A noninvasive growth case (i.e., ) in the same microenvironment () is also studied for comparison purposes.
Figure 7 shows the simulated growing tumors 100 days after initiation (plots showing the full growth history of the tumors are given in the Supplementary Information). The computed metrics for tumor morphology are given in Table 3. The primary tumors for both invasive [Figs. 7(b),(c) and (d)] and noninvasive [Fig. 7(a)] cases develop necrotic and quiescent regions. For invasive tumors, when the cell motility is small (i.e., ), the invasive cells do not form dendritic invasive branches but rather clump near the outer border of the proliferative rim [see Fig. 7(b)], forming bumpy invasive concentriclike shells with relatively small specific surface (e.g., on day 100). Such invasive shells significantly enhance the growth the primary tumor, i.e., the size of the primary tumor in Fig. 7(b) is much larger than in Figs. 7(a), (c) and (d). A quantitative comparison of the tumor sizes is given in the Supplementary Information.
By contrast, for larger cell motility, long dendritic invasive branches are developed as manifested by the large specific surface (e.g., on day 100 and on day 120 for ). In particular, one can clearly see that within the branches the cells tend to follow one another to form chains, as observed in experiments [4]. We emphasize that no rules are imposed to force the cells to follow on another in our CA model. This homotype attraction is purely due to the mechanical interaction between the invasive cells and the ECM, i.e., once a path of invasion is established by a leading invasive cell (by degrading the ECM), other invasive cells nearby turn to follow and enhance this path since the resistance for migration is minimized on a existing path. Furthermore, we can see that larger cell motility (i.e., high malignancy) leads to more invasive branches [see Figs. 7(c) and (d)] and thus, a larger specific area of the invasive pattern.
Effects of the ECM Rigidity
It is not very surprising that isotropic tumor shapes and and invasive patterns are developed in a homogeneous ECM with relative low density (i.e.,the ECM is soft) compared to the ECM degradation ability of the invasive tumor cells. However, real tumors are rarely isotropic, primarily due to the host microenvironment in which they grow, which we now explore.
Consider the invasive growth of a tumor in a much more rigid ECM than that in the previous section, i.e., . The invasiveness parameters used are , and . The snap shots of the growing tumor are shown in Fig. 8 and the tumor morphology metrics are given in Table 3. It can be clearly seen that both the size of the primary tumor and the extent of its invasive branches are much smaller than those of the tumors growing in a softer ECM (see Fig. 7). Importantly, although the ECM is still homogeneous, due to its high rigidity, the primary tumor develops an anisotropic shape with protrusions in the proliferation rim caused by the invasive branches (e.g., and on day 100). Since the invasive cells have degraded the ECM either completely or partially along the invasive branches, it is easier for the proliferative cells in the primary tumor to take these “weak spots” than to push against the rigid ECM themselves. Again, we emphasize that we do not force the cells to behave this way by imposing special CA rules; this behavior results purely from the mechanical interaction between the tumor and its host and the coupling dynamics of invasive and noninvasive tumor cells.
Effects of the ECM Heterogeneity: Random Distribution of ECM Density
The real host microenvironment for tumors are far from homogeneous in general. To investigate how heterogeneity of ECM affects the tumor growth dynamics, we use a random distribution of ECM density, i.e., for each ECM associated automaton cell, it density is a random number uniformly chosen from the interval (see Fig. 6(b)). The invasiveness parameters used are , and and snap shots of the growing tumor are shown in Fig. 9. The tumor morphology metrics are given in Table 4. Note that the primary tumor develops a rough surface and slightly anisotropic shape in the early growth stages (e.g., on day 50 and on day 80), which reflects the ECM heterogeneity [Figs. 9(a) and (b)]. Since the characteristic heterogeneity length scale is comparable to a single cell, its effects are diminished (e.g., on day 100 and on day 120) as the tumor grows larger and larger [Fig. 9(c) and (d)]. (In other words, on large length scales, the ECM is still effectively homogeneous.) However, the anisotropy in the invasive pattern (i.e., the extents of invasive branches in different directions) still persists (e.g., on day 100) even though the primary tumors almost resumes an isotropic shape.
Effects of the ECM Heterogeneity: Sinelike Distribution of ECM Density
To represent largescale heterogeneities in the ECM, we use a sinelike distribution of ECM density, i.e., for an automaton cell with centroid , the associated ECM density is given by
(3) 
where is the spatial dimension and is the edge length of the dimensional cubic simulation box. A twodimensional sinelike ECM distribution is shown in Fig. 6(c). The red spots correspond to large and high ECM rigidity; they can be considered as effective obstacles (e.g., brain ventricles) that hinder tumor growth.
Figures 10(a),(b) and (c),(d) show the snap shots of invasive tumors growing in the aforementioned ECM on day 80 and day 120, with invasiveness parameters , and and , and , respectively. The plots showing the full growth history is given in Supplementary Information and the tumor morphology metrics are given in Table 4. We can see that in the early growth stage, both the primary tumor and invasive pattern in the two cases are significantly affected by the ECM heterogeneity. In particular, the tumors are highly anisotropic in shape and the invasive branches clearly favor two orthogonal directions associated with low ECM densities (e.g., , for on day 80; and , for on day 80). For the case with large cellular motility, anisotropy effects are diminished in later growth stages (, for on day 120). For small cellular motility, anisotropy in both primary tumor shape and the invasive pattern persists (, for on day 120). Furthermore, one can see that again invasive cells with low motility significantly facilitate the growth of the primary tumor. However, instead of forming “bumpy” concentriclike shells as in homogeneous ECM, the invasive cells form large invasive cones, protruding into the ECM. These invasive cones are followed by weak protrusion of the proliferative rim, leading to bumpy surface of the primary tumor. The fact that such complex growth dynamics are only observed for tumors growing heterogeneous ECM emphasizes the crucial importance of understanding the effects of physical heterogeneity in cancer research.
Discussion
We have developed a novel cellular automaton (CA) model which, with just a few parameters, can produce a rich spectrum of growth dynamics for invasive tumors in heterogeneous host microenvironment. Besides robustly reproducing the salient features of branched invasive growth, such as least resistance and intrabranch homotype attraction observed in in vitro experiments, our model also enables us to systematically investigate the effects of microenvironment heterogeneity on tumor growth as well as the coupling growth of the primary tumor and the invasive cells. In particular, we have shown that in homogeneous ECM with low densities (i.e., soft microenvironment), both the shape of the primary tumor and invasive pattern are isotropic. For high cellular motility cases, the invasive cells form extended dendritic invasive branches; while for low cellular motility cases, the invasive cells clump near the primary tumor surface and form a bumpy concentriclike shell that facilitates the growth of the primary tumor. Tumors growing in a highly rigid homogeneous ECM can develop anisotropic shapes, facilitated by the invasive cells that degrade the ECM; both the tumor size and the extent of invasive branches are much smaller. In heterogeneous ECM, both the primary tumor and invasive pattern are significantly affected during the early growth stages, i.e., anisotropic shapes and patterns are developed to avoid high density/rigid regions of the ECM. If the characteristic length scale of the heterogeneities is comparable to the macroscopic tumor size, such effects can persist in later growth stages. In addition, invasive cells with large motility can significant diminish the anisotropy effects by their ECM degradation activities. We emphasize that we did not manipulate the behavior of cells by imposing artificial CA rules to give rise to these complex and rich growth dynamics. Instead, these are emergent behaviors that naturally arise due to various microscopicscale tumorhost interactions that are incorporated into our CA model, including the shortrange mechanical interaction between tumor cells and tumor stroma, and the degradation of extracellular matrix by the invasive cells.
Note that the growth dynamics of tumors growing in a heterogeneous microenvironment is distinctly different than those in a homogeneous microenvironment. This emphasizes the importance of understanding the effects of physical heterogeneity of the host microenvironment in modeling tumor growth. Here we just make a first attempt to take into account a simple level of host heterogeneity, i.e., by considering the ECM with variable density/rigidity. Currently, the invasion of the malignant cells into the host microenivronment is considered to be a consequence of invasive cell phenotype gained by mutation, and is not triggered by environmental stresses. However, the effects of environmental stresses can be easily taken into account. For example, a CA rule can be imposed that if the division probability of a malignant cell is significantly reduced by ECM rigidity, i.e., it is extremely difficult to push away/degrade ECM to make room for daughter cells, the malignant cell leaves the primary tumor and invades into soft regions of surrounding ECM. This would lead to reduced tumor invasion (i.e., development of the dendritic invasive branches) in soft microenvironments but enhanced invasion in rigid microenvironments [48].
Moreover, the spatialtemporal evolution of more complicated and realistic nutrient/ oxygen fields can be incorporated into our CA model. This can be achieved by solving the coupled nonlinear partial differential equations governing the evolution of the nutrient/oxygen concentrations as was done in Refs. [19] and [21]. Since the CA rules are given for any spatial dimension, our model is readily generalized to three dimensions. In addition, the model can be easily modified to incorporate other host heterogeneities, such as stromal cells, blood vessels and the shape anisotropy of the host organ [20, 21]. As currently implemented, a single 2D simulation takes less than 0.5 hours on a 32bit 1.56Gb Memory 1.44GHz dual core Dell Workstation. We expect that a 3D simulation will take no longer than 24 hours on a supercomputer when a proper parallel implementation is used.
Such an in silica tool not only enables one to investigate tumor growth in complex heterogeneous microenvironment that closely represents the real host microenvironments but also allows one to infer and even reconstruct individual host microenvironment given limited growth data of tumors (such as shape and size at various times). Such microstructural information of the individual host would be extremely valuable for developing individualized treatment strategies. For example, based on the host microstructure one can design special encapsulation and transport agents that maximize drug delivery efficiency [13].
In our current CA model, the microscopic parameters governing tumor invasion are variable and can be arbitrarily chosen within a feasible range as given in Table 1. Given sufficient and reliable experimental data of invasive tumor growth, the parameters in our CA model can be uniquely determined and thus, the model can produce robust predictions on the neoplastic progression. Although the current CA model is specifically implemented to reproduce and predict the growth dynamics of invasive solid tumors in vitro, further refinement of the model could eventually lead to the development of a powerful simulation tool that can be utilized clinically. For example, more complicated and realistic host heterogeneities such as the vascular structure, various stromal cells, the corresponding spatialtemporal evolution of the nutrient/oxygen concentrations as well as environmental stressinduced mutations should be incorporated as we described earlier. If the robustness of the refined model could be validated clinically, we would expect it to produce quantitative predictions for in vivo tumor growth, which are valuable for tumor prognosis and proposing individualized treatment strategies.
Acknowledgments
The authors are grateful to Robert Gatenby and Bob Austin for valuable comments on our manuscript. The authors are also grateful to the anonymous reviewers for their valuable comments.
Figure Legends
Time dependent terms  

Local tumor radius (varies with cell positions)  
Local maximum tumor extent (varies with cell positions)  
Characteristic proliferative rim thickness  
Characteristic livingcell rim thickness (determines necrotic fraction)  
Probability of division (varies with cell positions)  
Growth parameters  
Base probability of division, linked to celldoubling time (0.192 and 0.384)  
Base necrotic thickness, controlled by nutritional needs ( mm)  
Base proliferative thickness, controlled by nutritional needs ( mm)  
Invasiveness parameters  
Mutation rate (determines the number of invasive cells, 0.05)  
ECM degradation ability ()  
Cell motility (the number of “jumps” from one automaton cell to another, )  
Other terms  
ECM density (determines the ECM rigidity and varies with positions, ) 
Summarized here are definitions of the parameters for tumor growth and invasion, and all other (timedependent) quantities used in the simulations. For each parameter, the number(s) listed in parentheses indicates the value or range of values assigned to the parameters during the simulations. The values of the parameters are chosen such that the CA model can reproduce reported growth dynamics of GBM from the medical literature [4, 16, 20].
Metrics  Simulated MTS  Experimental data 

Specific surface  9.24  9.78 
Asphericity  1.09  1.12 
Angular anisotropy metric  0.17  0.19 
Noninvasive tumor in ECM with
Metrics  Day 50  Day 80  Day 100  Day 120 

Specific surface  1.23  1.13  1.09  1.04 
Asphericity  1.21  1.18  1.08  1.12 
Invasive tumor with in ECM with
Metrics  Day 50  Day 80  Day 100  Day 120 

0.66  0.38  0.21  0.19  
Specific surface  1.76  1.48  1.26  1.18 
Asphericity  1.23  1.12  1.08  1.06 
Angular anisotropy metric  0.13  0.29  0.33  0.17 
Invasive tumor with in ECM with
Metrics  Day 50  Day 80  Day 100  Day 120 

1.28  2.12  2.67  2.08  
Specific surface  1.94  3.92  3.67  3.28 
Asphericity  1.42  1.38  1.16  1.23 
Angular anisotropy metric  0.86  0.67  0.64  0.45 
Invasive tumor with in ECM with
Metrics  Day 50  Day 80  Day 100  Day 120 

2.14  2.43  2.64  2.89  
Specific surface  1.71  4.28  7.89  9.73 
Asphericity  1.38  1.27  1.13  1.8 
Angular anisotropy metric  1.25  0.68  0.41  0.18 
Invasive tumor with in ECM with
Metrics  Day 50  Day 80  Day 100  Day 120 

  5.17  3.89  2.63  
Specific surface  1.21  3.40  3.91  5.79 
Asphericity  1.33  1.36  1.40  1.56 
Angular anisotropy metric    1.32  1.02  0.65 
Invasive tumor with in random ECM
Metrics  Day 50  Day 80  Day 100  Day 120 

2.54  4.14  2.13  2.78  
Specific surface  2.47  3.97  4.65  8.98 
Asphericity  1.32  1.34  1.18  1.15 
Angular anisotropy metric  0.63  0.87  0.64  0.19 
Invasive tumor with in sinelike ECM
Metrics  Day 50  Day 80  Day 100  Day 120 

2.78  2.46  1.28  0.86  
Specific surface  1.89  2.95  2.73  1.92 
Asphericity  1.42  1.61  1.49  1.26 
Angular anisotropy metric  1.23  1.18  1.09  0.98 
Invasive tumor with in sinelike ECM
Metrics  Day 50  Day 80  Day 100  Day 120 

5.24  3.86  3.13  2.96  
Specific surface  2.51  4.12  6.13  8.76 
Asphericity  1.19  1.67  1.48  1.21 
Angular anisotropy metric  1.36  1.33  0.88  0.23 
References
 Coffey DS (1998) Self organization, complexity and chaos: The new biology for medicine. Nat Med 4: 882885.
 Hanahan D, Weinberg RA (2000) The hallmarks of cancer. Cell 100:5770.
 Fearon ER, Vogelstein B (1990) A genetic model for colorectal tumorigenesis. Cell 61:759767.
 Deisboeck TS, Berens ME, Kansal AR, Torquato S, Rachamimov A, Louis DN, Chiocca EA (2001) Patterns of selforganization in tumor systems: Complex growth dynamics in a novel brain tumor spheroid model. Cell Proliferat. 34:115134.
 Fidler IJ (2003) The pathogenesis of cancer metastasis: the “seed and soil” hypothesis revisited. Nat Rev Cancer 3:453458.
 Kerbel RS (1990) Growth dominance of the metastatic cancer cell: cellular and molecular aspects. Adv Cancer Res 55:87132.
 Liotta LA, Kohn EC (2003) Cancer’s deadly signature. Nat Genet 33:1011.
 Chen ZJ, Gillies GT, Broaddus WC, Prabhu SS, Fillmore H, Mitchell RM, Corwin FD, Fatouros PP (2004) A realistic brain tissue phantom for intraparenchymal infusion studies. J. Neurosurgery 101:314322.
 Frieboes HB, Edgerton ME, Fruehauf JP, Rose FR, Worrall LK, Gatenby RA, Ferrari M, Cristini V (2009) Prediction of drug response in breast cancer using integrative experimental/computational modeling. Cancer Res. 69:44844492.
 Crossa FR, Tinkelenberga AH (1991) A potential positive feedback loop controlling cln1 and cln2 gene expression at the start of the yeast cell cycle. Cell 65:875883.
 Brand U, Fletcher JC, Hobe M, Meyerowitz EM, Simon1dagger R (2000) Dependence of stem cell fate in arabidopsis on a feedback loop regulated by clv3 activity. Science 289:617619.
 Kitano H (2004) Cancer as a robust system: implications for anticancer therapy. Nat Rev Cancer 4:227235.
 Torquato S (2011) Toward an Ising model of cancer and beyond. Phys Biol 8:015017.
 Byrne HM (2010) Dissecting cancer through mathematics: from the cell to the animal model. Nat Rev Cancer 10:221230.
 Anderson ARA, Chaplain MAJ (1998) Continuous and discrete mathematical models of tumorinduced angiogenesis. Bull Math Biol 60:857â900.
 Kansal AR, Torquato S, Harsh GR, Chiocca EA, Deisboeck TS (2000) Simulated brain tumor growth using a threedimensional cellular automaton. J Theor Biol 203:367382.
 Kansal AR, Torquato S, Chiocca EA, Deisboeck TS (2000) Emergence of a subpopulation in a computational model of tumor growth. J Theor Biol 207:431441.
 Schmitz JE, Kansal AR, Torquato S (2002) A Cellular Automaton Model of Brain Tumor Treatment and Resistance. J Theor Med 4:223239.
 Gevertz JL, Torquato S (2006) Modeling the effects of vasculature evolution on early brain tumor growth. J Theor Biol 243:517531.
 Gevertz JL, Gillies GT, Torquato S (2008) Simulating tumor growth in confined heterogeneous environments. Phys Biol 5:036010.
 Gevertz JL, Torquato S (2009) Growing heterogeneous tumors in silico. Phys Rev E 80:051910.
 Anderson ARA (2005) A hybrid mathematical model of solid tumor invasion: the important of cell adhesion. Math Med Biol 22:163186.
 Anderson ARA, Weaver AM, Cummings PT, Quaranta V (2005) Tumor morphology and phenotypic evolution driven by selective pressure from microenvironment. Cell 127:905915.
 Bankhead III A, Heckendorn RB (2007) Using evolvable genetic cellular automata to model breast cancer. Genet Program Evolvable Mach 8:381393.
 Gatenby RA (1996) Application of competition theory to tumour growth: implications for tumour biology and treatment. Eur J Cancer 32:722726.
 Gatenby RA, Gawlinski ET (1996) A reactiondiffusion model of cancer invasion. Cancer Res. 56:57455753.
 Gatenby RA, Gawlinski ET, Gmitro AF, Kaylor B, Gillies RJ (2006) Acidmediated tumor invasion: a multidisciplinary study. Cancer Res. 66:52165223.
 Frieboes HB, Zheng XM, Sun CH, Tromberg B, Gatenby R, Gristini V (2006) An integrated computational/experimental model for tumor invasion. Cancer Res. 66:15971604.
 Bellomo N, Preziosi L (2000) Modelling and mathematical problems related to tumor evolution and its interaction with the immune system. Math Comput Model 32:413452.
 Scalerandi M, Sansone BC, Condat CA (2001) Diffusion with evolving sources and competing sinks: Development of angiogenesis. Phys Rev E 65:011902.
 Scalerandi M, Sansone BC (2002) Inhibition of vascularization in tumor growth. Phys Rev Lett 89:218101.
 Kim Y, Friedman A (2010) Interaction of tumor with its microenvironment: A mathematical model. Bull Math Biol 72:10291068.
 Stein AM, Demuth T, Mobley D, Berens M (2007) A mathematical model of glioblastoma tumor spheroid invasion in a threedimensional in vitro experiment. Biophys J 92:356365.
 Torquato S (2002) Random Heterogeneous Materials: Microstructure and Macroscopic Properties. (SpringerVerlag, New York).
 Torquato S, Stillinger FH (2010) Jammed hardparticle packings: From Keplerto Bernal and Beyond. Rev Mod Phys 82:26332672.
 Patel AB, Gibson WT, Gibson MC, Nagpal R (2009) Modeling and inferring cleavage patterns in proliferating epithelia. PLoS Compt Biol 5:e1000412.
 Gevertz JL, Torquato S (2008) A novel threephase model of brain tissue microstructure. PLoS Comput. Biol. 4:e1000152.
 Burridge K, ChrzanowskaWodnicka M (1996) Focal adhesions, contractability, and signalling. Annu. Rev. Cell Dev. Biol. 12:463518.
 Sarntinoranont M, Rooney F, Ferrari M (2003) Interstitial stress and fluid pressure within a growing tumor. Annals of Biomedical Engineering 31:327335.
 Gordon VD, Valentine MT, Gardel ML, AndorArdo D, Dennison S, Bogdanov AA, Weitz DA, Deisboeck TS (2003) Measuring the mechanical stress induced by expanding multicellular tumor system: a case stduy. Experimental Cell Res. 289:5866.
 Liotta LA, Rao CN, Barsky SH (1983) Tumor invasion and the extracellular matrix. Lab. Invest. 49:636649.
 Boyle JO, Hakim J, Koch W (1993) The incidence of p53 mutations increases with progression of head and neck cancer. Cancer Res. 53:44774480.
 StetlerStevenson WG, Aznavoorian S, Liotta LA (1993) Tumor cell interactions with the extracellular matrix during invasion and metastasis. Annu. Rev. Cell Biol. 9:541573.
 Lawrence JA, Steeg PS (1996) Mechanisms of tumor invasion and metastasis. World J. Urol. 14:124130.
 Torquato S, Jiao Y (2009) Dense Packings of the Platonic and Archimedean Solids. Nature 460:876879.
 Hoshino T, Wilson CB (1979) Cell kinetic analyses of human malignant brain tumors (gliomas). Cancer 44:956962.
 Pertuiset B, Dougherty D, Cromeryer C, Hoshino T, Berger M, Rosenblum ML (1985) Stem cell studies of human malignant brain tumors. Part 2: Proliferation kinetics of braintumor cells in vitro in earlypassage cultures. J. Neurosurg. 63:426432.
 Guiot C, Pugno N, Delsanto PP, Deisboeck TS (2007) Physical aspects of cancer invasion. Phys. Biol. 4:P1P6.