Emergent Behaviors from a Cellular Automaton Model for Invasive Tumor Growth in Heterogeneous Microenvironments

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 non-invasive (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 microscopic-scale tumor-host interactions, including the short-range mechanical interactions between tumor cells and tumor stroma, degradation of the 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 paths of cells and intrabranch homotype attraction, we also predict nontrivial coupling between 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 tumor-host interaction. The capability of our CA model suggests that sophisticated in silico tools could eventually be utilized in clinical situations 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, integrin-mediated 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 Figure 1. Such invasive behaviors are characterized by intrabranch homotype attraction and least-resistance paths of cells [4].
Although recent progress has been made in understanding certain aspects of the complex tumor-host 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 distinctly different mechanisms for tumorigenesis, 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 [13].
In response to the challenge to develop an ''Ising'' model for cancer growth [13], we generalize here a cell-based 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 host microenvironment. Our generalized cellular automaton model takes into account a variety of microscopic-scale tumor-host 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 paths of cells and intrabranch homotype attraction. The model also predicts nontrivial coupling between 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

Author Summary
The goal of the present work is to develop an efficient single-cell 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 chain-like invasive branches form nor how they interact with the host microenvironment and the primary tumor. Our CA model, which incorporates a variety of microscopic-scale tumorhost interactions (e.g., the mechanical interactions between tumor cells and tumor stroma, degradation of the 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 for clinical purposes capable of predicting neoplastic progression and suggesting individualized optimal treatment strategies.
dynamics and lead to a variety of tumor morphologies. These emergent behaviors naturally arise due to various microscopicscale tumor-host 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 could 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.

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 Figure 2). In particular, nonoverlapping spheres are randomly and sequentially placed in a prescribed region until there is no void space left 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 biological 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 real-cell aggregates with a relatively 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 in two and three dimensions, respectively. Since our new CA model explicitly takes into account the interactions between a single biological 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 15{20 mm and the linear size of the 2D simulation domain is approximately 5 mm, which contains *100 000 automaton cells. In the current model, we mainly focus on the effects of the ECM macromolecule density and ECM degradation by malignant cells on tumor growth. Henceforth, we will refer to the host microenvironment (or tumor stroma) as the ''ECM'' for simplicity. Each ECM associated automaton cell is assigned a particular density r ECM , representing the density of the ECM macromolecules within the automaton cell. A tumor cell can occupy an ECM associated automaton cell only if the density of this automaton cell r ECM~0 , 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 provides mechanical support for the tissue (such as collagen) 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 contain fibroblasts, which actively produce ECM The associated point configuration for the tessellation, generated by randomly placing nonoverlap circular disks in a prescribed region, i.e., the random sequential addition process [34]. doi:10.1371/journal.pcbi.1002314.g002 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 the 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 growth-permitting region. This is done 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 growth-permitting 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 growth-region 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 growth-permitting 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 non-invasive. Following Ref. [16], the non-invasive cells remain in the primary tumor and can be proliferative, quiescent or necrotic, depending on the nutrient supply they get. For avascular tumor growth, our focus here, the nutrients available to the tumor cells are essentially the nutrient concentrations that diffuse into the tumor through its surface. As the tumor grows, the amount of nutrient 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 d n away from the tumor surface become necrotic (see details in the next section). The diffusion length d n (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 d p 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 the ECM degradation ability x (i.e., the matrix-degradation enzymes) and motility m that allow them to leave the primary tumor and invade into surrounding microenvironment [43]. We consider that 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 r ECM~0 ). 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 m and x (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 m 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 r ECM [(0,1) is assigned to each automaton cell within the growth-permitting 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: N 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.
N All ECM associated automaton cells and tumorous necrotic cells are inert (i.e., they do not change type).
N Quiescent cells more than a certain distance d n 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 d n for quiescent cells to turn necrotic is computed as follows: where a is prescribed parameter (see Table 1), d is the spatial dimension and L t is the distance between the geometric centroid (i.e., the center) x c of the tumor and the tumor edge cell that is closest to the quiescent cell under consideration. The position of the tumor centroid x c is given by where N 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.
N Each proliferative cell will attempt to divide with probability p div 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 p div depends on both the physical confinement imposed by the boundary of the growth-permitting region and the local mechanical interaction between the tumor cells and the ECM, i.e., if any ECM associated automaton cell within the predefined growth distance is in the growthpermitting microenvironment 0 if no ECM associated automaton cell within the predefined growth distance is in the growthpermitting microenvironment where p 0 is the base probability of division (see Table 1), r is the distance of the dividing cell from the tumor centroid, L max is the distance between the closest growth-permitting boundary cell in the direction of tumor growth and the tumor's geometric centroid x c and r ECM is the ECM density of the automaton cell to be taken by the new tumor cell. When an ECM associated automaton cell is taken by a tumor cell, its density is set to be zero. The predefined growth distance (d p ) is described in a following bullet point.
N If a proliferative cell divides, it can produce a mutant daughter cell possessing an invasive phenotype with a prescribed probability c (i.e., the mutation rate). The invasive daughter cell gains ECM degradation ability x and motility m, which enable it to leave the tumor and invade into the 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.
N A proliferative cell turns quiescent if there is no space available for the placement of a daughter cell within a distance d p from the proliferative cell, which is given by where b a nutritional parameter (see Table 1), d is the spatial dimension and L t is the distance between the geometric tumor centroid x c and the tumor edge cell that is closest to the proliferative cell under consideration.
N 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 in that automaton cell is completely degraded. For an invasive cell with motility m and ECM degradation ability x, it will make m attempts to degrade the ECM in the neighboring automaton cells and jump to these automaton cells, where m is an arbitrary integer in ½0, m. For each attempt, the surrounding ECM density r ECM is decreased by dr, where dr is an arbitrary number in ½0, x. 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 m and dr). When the ECM in multiple neighboring automaton cells of the invasive cell are completely degraded (i.e., r ECM~0 ), the invasive cell moves in a direction that maximizes the nutrient 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 Figure 3. We note that non-invasive tumor growth can be studied by imposing a mutation rate c~0. This enables us to compare the growth dynamics of invasive and non-invasive tumors and in turn to investigate the effects of the coupling between the growth dynamics of the primary tumor mass and the invasive cells. Although we only consider spherical growthpermitting 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 a 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 b of the invasive area over the tumor area (defined below), the specific surface s of the invasive pattern, the asphericity a of the primary tumor and the angular anisotropy metric y 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 A T is defined as the area of the circumcircle of the primary tumor (see Figure 4(a)) and the invasive area A I is the area of the region between the effective circumcircle of the invasive pattern and the circumcircle of the primary tumor (see Figure 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 b~A I =A T as a function of time t reflects the degree of coupling between the primary tumor and the invasive cells. If b(t) is linear in t, there is no coupling; otherwise the two are coupled.
The specific surface s [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, s is inversely proportional to the size  Since s 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 s with respect to 2=R T for an arbitrary-shaped tumor with effective radius R T (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 s. The asphericity a of the primary tumor is defined as the ratio of the radius of circumcircle R c of the primary tumor over its incircle radius R in [45], i.e., a~R c =R in (see Figure 4(b)). A large a 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 y. In particular, the entire invasive pattern is evenly divided into n a sectors with lines emanating from the tumor center (see Figure 4(c)). The angular anisotropy metric y is defined as where '(i) is the average length of the invasive branches within the ith sector and is the average length of all invasive branches. For tumors with invasive branches of similar lengths that are uniformly angularly distributed, the metric y is small. Large fluctuations of both invasive branch length and angular distribution can lead to large y values. In the following, we use n a~1 6 to compute y for the simulated invasive tumors.

Model Validation
To verify the robustness and predictive capacity of our CA model, we first employ it to reproduce quantitatively the observed invasive growth of a GBM multicellular tumor spheroid (MTS) in vitro [4]. In particular, the boundary of the growth-permitting region is considered to be vascularized, i.e., a growing tumor can receive oxygen and nutrients from the growth-permitting region. A constant radially symmetric nutrient/oxygen gradient in the growth-permitting 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 growth-permitting region with homogeneous ECM and tumor growth is started. This corresponds to an initial MTS with diameter D MTS &310 mm which is consistent with the in vitro experiment set-up [4]. The following values of the growth and invasiveness parameters are used: p 0~0 :384, a~0:58 mm 1=2 , b~0:30 mm 1=2 , c~0:05, x~0:55, m~3. Note that the value of p 0 corresponds to a cell doubling time of 40 hours, which is consistent with the reported experimental data [4]. A small value of the ECM density r ECM~0 :15 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 growth-permitting 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 Figure 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 in our simulations.
The ratio of the invasion area over the primary tumor area b~A I =A T as a function of time for the simulated tumor is computed and compared to the reported experimental data [4] (see Figure 5(c)). One can clearly see that our simulation results agree with experimental data very well. Moreover, the deviation of b(t) from a linear function of t 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 s of the invasive pattern, the sphericity a of the primary tumor and the angular anisotropy metric y for the invasive branches are computed from our simulation results and from the image of invasive MTS in Figure 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 distributions of the ECM density, i.e., homogeneous, random and sinusoidal-like, to systematically study the effects of microenvironment heterogeneity on invasive tumor growth (see Figure 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 than exhibited here, previous 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 by the nutrient/oxygen gradient, our constantgradient approximation is a very reasonable one.
In the beginning of the simulation, a proliferative tumor cell is introduced at the center of the growth-permitting 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 variables 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 p 0~0 :192, 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 tumor-host 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 r ECM~0 :45. In particular, we consider three invasive cases with the same mutation rate c~0:05 and ECM degradation ability x~0:9, but different cell motility m~1,2,3: A non-invasive growth case (i.e., c~0) in the same microenvironment (r ECM~0 :45) 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 Figure S1, S2, S3, S4). The computed metrics for tumor morphology are given in Table 3. The primary tumors for both invasive (Figure 7(b),(c) and (d)) and non-invasive (Figure 7(a)) cases develop necrotic and quiescent regions. For invasive tumors, when the cell motility is small (i.e., m~1), the invasive cells do not form dendritic invasive branches but rather clump near the outer border of the proliferative rim (see Figure 7(b)),  forming bumpy invasive concentric-like shells with relatively small specific surface (e.g., s~1:09 on day 100). Such invasive shells significantly enhance the growth the primary tumor, i.e., the size of the primary tumor in Figure 7(b) is much larger than in Figure 7(a), (c) and (d). A quantitative comparison of the tumor sizes is shown in Figure S5.
By contrast, for larger cell motility, long dendritic invasive branches are developed as manifested by the large specific surface (e.g., s~7:89 on day 100 and s~9:73 on day 120 for m~3). 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 one 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 pioneering 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 Figure 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 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., r ECM~0 :85. The invasiveness parameters used are c~0:05, m~3 and x~0:9. The snap shots of the growing tumor are shown in Figure 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 Figure 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., a~1:40 and y~1:02 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 between the growth dynamics of the invasive and non-invasive tumor cells.
Effects of the ECM heterogeneity: random distribution of the ECM density. The real host microenvironment for tumors are far from homogeneous in general. To investigate how ECM heterogeneity affects the tumor growth dynamics, we use a random distribution of the ECM density, i.e., for each ECM associated automaton cell, its density r ECM is a random number uniformly chosen from the interval ½0, 1 (see Figure 6(b)). The invasiveness parameters used are c~0:05, m~3 and x~0:9 and the snap shots of the growing tumor are shown in Figure 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., a~1:32 on day 50 and a~1:34 on day 80), which reflects the ECM heterogeneity (Figure 9(a) and (b)).
Since the characteristic heterogeneity length scale is comparable to a single cell, its effects are diminished (e.g., a~1:18 on day 100 and a~1:15 on day 120) as the tumor grows larger and larger ( Figure 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., y~0:64 on day 100) even though the primary tumors almost resumes an isotropic shape.
Effects of the ECM heterogeneity: sinusoidal-like distribution of the ECM density. To represent large-scale heterogeneities in the ECM, we use a sinusoidal-like distribution of the ECM density, i.e., for an automaton cell with centroid (x 1 ,x 2 , . . . ,x d ), the associated ECM density is given by where d is the spatial dimension and L is the edge length of the ddimensional cubic simulation box. A two-dimensional sinusoidallike ECM density distribution is shown in Figure 6(c). The red spots correspond to large r ECM and high ECM rigidity; they can be considered as effective obstacles (e.g., brain ventricles) that hinder tumor growth. show the snap shots of invasive tumors growing in the aforementioned ECM on day 80 and day 120, with invasiveness parameters c~0:05, m~1 and x~0:9 and c~0:05, m~3 and x~0:9, respectively. The plots showing the full growth history is given in Figure S6 and S7, 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., a~1:61, y~1:18 for m~1 on day 80; and a~1:67, y~1:33 for m~3 on day 80). For the case with large cellular motility, anisotropy effects are diminished in later growth stages (a~1:21, y~0:23 for m~3 on day 120). For small cellular motility, anisotropy in both primary tumor shape and the invasive pattern persists (a~1:26, y~0:98 for m~1 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'' concentric-like 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 paths of cells 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 between the growth dynamics 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 concentric-like 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 significantly 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 microscopic-scale tumor-host interactions that are incorporated into our CA model, including the short-range mechanical interaction between the tumor cells and tumor stroma, and the degradation of extracellular matrix by the invasive cells.
It is noteworthy that the growth dynamics of tumors 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 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]. Indeed, we have very recently generalized the CA model reported here to explicitly take into account the pressure exerted on the tumor due to the deformation of its surrounding ECM as well as the local geometry of the tumor-host interface to study mechanical-stress induced tumor morphology instability [49].
Moreover, the spatial-temporal 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 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 32-bit 1.56 Gb Memory 1.44 GHz 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 silico 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 could be uniquely determined and thus, the model could produce robust predictions about 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 could some day be utilized clinically. For example, more complicated and realistic host heterogeneities such as the vascular structure, various stromal cells, the corresponding spatial-temporal evolution of the nutrient/oxygen concentrations as well as environmental stress-induced 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 could potentially be valuable for tumor prognosis and proposing individualized treatment strategies.  Figure S5 The scaled volume V T =V 0 of the primary tumors as function of time t. The volume of the tumor V T is scaled with respect to the volume of the growth-permitting region V 0 . Note that invasive cells with small motility m~1 significantly enhance the growth of the primary tumor. For invasive cells with intermediate (m~2) and large (m~3) motility values, since such invaisve cells leave the primary tumor and do not contribute to the tumor volume, the size of primary tumors are smaller than the non-invasive case in early growth stages. In later stages, the invasive cells have degraded a signficant amount of extracellular matrix, leading to faster growth of the primary tumor than the non-invasive case. The final size of the tumor is determined by the volume of the growth-permitting regions. For invasive cases, at later growth stages, many invasive cells clump at the bounary of the growth-permitting region, which do not contribute to the tumor volume. Therefore, invasive tumors plateau at a smaller size than the non-invasive tumor. (EPS) Figure S6 Evolution of a simulated invasive tumor with cellular motility m~1 in the heterogeneous ECM with a sinusoidal-like density distribution. The mutation rate is c~0:05 and ECM degradation ability is x~0:9.