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

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 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.

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 tumor-host 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.

Citation: Jiao Y, Torquato S (2011) Emergent Behaviors from a Cellular Automaton Model for Invasive Tumor Growth in Heterogeneous Microenvironments. PLoS Comput Biol 7(12): e1002314. https://doi.org/10.1371/journal.pcbi.1002314

Editor: Feilim Mac Gabhann, Johns Hopkins University, United States of America

Received: May 13, 2011; Accepted: November 2, 2011; Published: December 22, 2011

Copyright: © 2011 Jiao, Torquato. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Funding: The research described was supported by the National Cancer Institute under Award NO. U54CA143803. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Cancer Institute or the National Institutes of Health. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

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 , . 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) .

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 , – and their treatment , . 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 . 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 . 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 . Such invasive behaviors are characterized by intrabranch homotype attraction and least-resistance paths of cells .

Download:

Figure 1. GBM multicelluar tumor spheroid (MTS) gel assay showing dendritic invasive branches.

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 , based on centers of spheres in a packing generated by a random sequential addition (RSA) process , (see ). 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 , . 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 . 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 . 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.

Download:

Figure 2. A 2D Voronoi tessellation and the associated point configuration.

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) , , . 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 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 , . 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. , 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. , characteristic diffusion lengths are employed to determine the states of a non-invasive 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 –. 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 , which can gain a variety of degrees of the ECM degradation ability (i.e., the matrix-degradation enzymes) and motility that allow them to leave the primary tumor and invade into surrounding microenvironment . 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 ). 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 , 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 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:

• 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 ), is the spatial dimension and is the distance between the geometric centroid (i.e., the center) 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 bywhere 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 growth-permitting region and the local mechanical interaction between the tumor cells and the ECM, i.e.,where is the base probability of division (see ), is the distance of the dividing cell from the tumor centroid, is the distance between the closest growth-permitting 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 an ECM associated automaton cell is taken by a tumor cell, its 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 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.
• 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 bywhere a nutritional parameter (see ), 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 in that automaton cell 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 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.

Download:

Table 1. Parameters and terms in the CA model.

Download:

Figure 3. Illustration of cellular automaton rules.

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 the invasive area over the 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. , the tumor area is defined as the area of the circumcircle of the primary tumor (see ) 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 ). 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.

Download:

Figure 4. Schematic illustration of the quantities in the definitions of tumor morphology metrics.

Results

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 . 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 which is consistent with the in vitro experiment set-up . The following values of the growth and invasiveness parameters are used: , , , , , . Note that the value of corresponds to a cell doubling time of 40 hours, which is consistent with the reported experimental data . A small value of the ECM density is used, which corresponds to the soft DMEM medium used in the experiment . 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.

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 (see ). 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.

Download:

Figure 5. Simulated invasive growth of MTS in vitro.

Download:

Table 2. Comparison of tumor morphology metrics associated with simulated MTS and experimental data at 24 hours after tumor initialization.

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 ). 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 growth-permitting region is considered to be vascularized with a constant radially symmetric nutrient/oxygen gradient in the growth-permitting 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 , . Since the directions of cell motions are determined by the nutrient/oxygen gradient, our constant-gradient approximation is a very reasonable one.

Download:

Figure 6. Different distributions of the ECM densities.

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 non-invasive growth case (i.e., ) in the same microenvironment () is also studied for comparison purposes.

shows the simulated growing tumors 100 days after initiation (plots showing the full growth history of the tumors are given in , , , ). The computed metrics for tumor morphology are given in . The primary tumors for both invasive () and non-invasive () 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 ), forming bumpy invasive concentric-like 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 is much larger than in . A quantitative comparison of the tumor sizes is shown in .

Download:

Figure 7. Simulated growing tumors in homogeneous ECM with on day 100 after initiation.

Download:

Table 3. Morphology metrics for simulated tumors growing in homogeneous ECM.

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., . The invasiveness parameters used are , and . The snap shots of the growing tumor are shown in and the tumor morphology metrics are given in . 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 ). 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 between the growth dynamics of the invasive and non-invasive tumor cells.

Download:

Figure 8. Evolution of simulated tumor in homogeneous ECM with .

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 is a random number uniformly chosen from the interval (see ). The invasiveness parameters used are , and and the snap shots of the growing tumor are shown in . The tumor morphology metrics are given in . 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 (). 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 (). (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.

Download:

Figure 9. Evolution of simulated tumor in random ECM.

Download:

Table 4. Morphology metrics for simulated tumors growing in heterogeneous ECM.

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 , the associated ECM density is given bywhere is the spatial dimension and is the edge length of the -dimensional cubic simulation box. A two-dimensional sinusoidal-like ECM density distribution is shown in . 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.

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 and , and the tumor morphology metrics are given in . 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” 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.

Download:

Figure 10. Simulated growing tumors in the ECM having a sinusoidal-like density distribution with different cellular motilities ().

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 . 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 .

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. and . 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 , . 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 .

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 . 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.

Supporting Information

 

The scaled volume of the primary tumors as function of time . The volume of the tumor is scaled with respect to the volume of the growth-permitting region . Note that invasive cells with small motility significantly enhance the growth of the primary tumor. For invasive cells with intermediate () and large () 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.

https://doi.org/10.1371/journal.pcbi.1002314.s005

(EPS)

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.

Author Contributions

Conceived and designed the experiments: ST. Performed the experiments: YJ ST. Analyzed the data: YJ ST. Contributed reagents/materials/analysis tools: YJ ST. Wrote the paper: YJ ST.

References

1. Coffey DS (1998) Self organization, complexity and chaos: The new biology for medicine. Nat Med 4: 882–885.
2. Hanahan D, Weinberg RA (2000) The hallmarks of cancer. Cell 100: 57–70.
3. Fearon ER, Vogelstein B (1990) A genetic model for colorectal tumorigenesis. Cell 61: 759–767.
4. Deisboeck TS, Berens ME, Kansal AR, Torquato S, Rachamimov A, et al. (2001) Patterns of self-organization in tumor systems: Complex growth dynamics in a novel brain tumor spheroid model. Cell Prolif 34: 115–134.
5. Fidler IJ (2003) The pathogenesis of cancer metastasis: The “;seed and soil” hypothesis revisited. Nat Rev Cancer 3: 453–458.
6. Kerbel RS (1990) Growth dominance of the metastatic cancer cell: Cellular and molecular aspects. Adv Cancer Res 55: 87–132.
7. Liotta LA, Kohn EC (2003) Cancer's deadly signature. Nat Genet 33: 10–11.
8. Chen ZJ, Gillies GT, Broaddus WC, Prabhu SS, Fillmore H, et al. (2004) A realistic brain tissue phantom for intraparenchymal infusion studies. J Neurosurg 101: 314–322.
9. Frieboes HB, Edgerton ME, Fruehauf JP, Rose FR, Worrall LK, et al. (2009) Prediction of drug response in breast cancer using integrative experimental/computational modeling. Cancer Res 69: 4484–4492.
10. 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: 875–883.
11. 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: 617–619.
12. Kitano H (2004) Cancer as a robust system: Implications for anticancer therapy. Nat Rev Cancer 4: 227–235.
13. Torquato S (2011) Toward an Ising model of cancer and beyond. Phys Biol 8: 015017.
14. Byrne HM (2010) Dissecting cancer through mathematics: From the cell to the animal model. Nat Rev Cancer 10: 221–230.
15. Anderson ARA, Chaplain MAJ (1998) Continuous and discrete mathematical models of tumor- induced angiogenesis. Bull Math Biol 60: 857900.
16. Kansal AR, Torquato S, Harsh GR, Chiocca EA, Deisboeck TS (2000) Simulated brain tumor growth using a three-dimensional cellular automaton. J Theor Biol 203: 367–382.
17. Kansal AR, Torquato S, Chiocca EA, Deisboeck TS (2000) Emergence of a subpopulation in a computational model of tumor growth. J Theor Biol 207: 431–441.
18. Schmitz JE, Kansal AR, Torquato S (2002) A Cellular Automaton Model of Brain Tumor Treat- ment and Resistance. J Theor Med 4: 223–239.
19. Gevertz JL, Torquato S (2006) Modeling the e®ects of vasculature evolution on early brain tumor growth. J Theor Biol 243: 517–531.
20. Gevertz JL, Gillies GT, Torquato S (2008) Simulating tumor growth in confined heterogeneous environments. Phys Biol 5: 036010.
21. Gevertz JL, Torquato S (2009) Growing heterogeneous tumors in silico. Phys Rev E 80: 051910.
22. Anderson ARA (2005) A hybrid mathematical model of solid tumor invasion: The important of cell adhesion. Math Med Biol 22: 163–186.
23. Anderson ARA, Weaver AM, Cummings PT, Quaranta V (2005) Tumor morphology and pheno- typic evolution driven by selective pressure from microenvironment. Cell 127: 905–915.
24. Bankhead A III, Heckendorn RB (2007) Using evolvable genetic cellular automata to model breast cancer. Genet Program Evolvable Mach 8: 381–393.
25. Gatenby RA (1996) Application of competition theory to tumour growth: Implications for tumour biology and treatment. Eur J Cancer 32: 722–726.
26. Gatenby RA, Gawlinski ET (1996) A reaction-di®usion model of cancer invasion. Cancer Res 56: 5745–5753.
27. Gatenby RA, Gawlinski ET, Gmitro AF, Kaylor B, Gillies RJ (2006) Acid-mediated tumor invasion: A multidisciplinary study. Cancer Res 66: 5216–5223.
28. Frieboes HB, Zheng XM, Sun CH, Tromberg B, Gatenby R, et al. (2006) An integrated computational/experimental model for tumor invasion. Cancer Res 66: 1597–1604.
29. Bellomo N, Preziosi L (2000) Modelling and mathematical problems related to tumor evolution and its interaction with the immune system. Math Comput Model 32: 413–452.
30. Scalerandi M, Sansone BC, Condat CA (2001) Diffusion with evolving sources and competing sinks: Development of angiogenesis. Phys Rev E 65: 011902.
31. Scalerandi M, Sansone BC (2002) Inhibition of vascularization in tumor growth. Phys Rev Lett 89: 218101.
32. Kim Y, Friedman A (2010) Interaction of tumor with its micro-environment: A mathematical model. Bull Math Biol 72: 1029–1068.
33. Stein AM, Demuth T, Mobley D, Berens M (2007) A mathematical model of glioblastoma tumor spheroid invasion in a three-dimensional in vitro experiment. Biophys J 92: 356–365.
34. Torquato S (2002) Random Heterogeneous Materials: Microstructure and Macroscopic Properties. New York: Springer-Verlag.
35. Torquato S, Stillinger FH (2010) Jammed hard-particle packings: From Kepler to Bernal and Beyond. Rev Mod Phys 82: 2633–2672.
36. Patel AB, Gibson WT, Gibson MC, Nagpal R (2009) Modeling and inferring cleavage patterns in proliferating epithelia. PLoS Compt Biol 5: e1000412.
37. Gevertz JL, Torquato S (2008) A novel three-phase model of brain tissue microstructure. PLoS Comput Biol 4: e1000152.
38. Burridge K, Chrzanowska-Wodnicka M (1996) Focal adhesions, contractability, and signalling. Annu Rev Cell Dev Biol 12: 463–518.
39. Sarntinoranont M, Rooney F, Ferrari M (2003) Interstitial stress and fluid pressure within a growing tumor. Ann Biomed Eng 31: 327–335.
40. Gordon VD, Valentine MT, Gardel ML, Andor-Ardo D, Dennison S, et al. (2003) Measuring the mechanical stress induced by expanding multicellular tumor system: A case stduy. Exp Cell Res 289: 58–66.
41. Liotta LA, Rao CN, Barsky SH (1983) Tumor invasion and the extracellular matrix. Lab Invest 49: 636–649.
42. Boyle JO, Hakim J, Koch W (1993) The incidence of p53 mutations increases with progression of head and neck cancer. Cancer Res 53: 4477–4480.
43. Stetler-Stevenson WG, Aznavoorian S, Liotta LA (1993) Tumor cell interactions with the extra- cellular matrix during invasion and metastasis. Annu Rev Cell Biol 9: 541–573.
44. Lawrence JA, Steeg PS (1996) Mechanisms of tumor invasion and metastasis. World J Urol 14: 124–130.
45. Torquato S, Jiao Y (2009) Dense Packings of the Platonic and Archimedean Solids. Nature 460: 876–879.
46. Hoshino T, Wilson CB (1979) Cell kinetic analyses of human malignant brain tumors (gliomas). Cancer 44: 956–962.
47. Pertuiset B, Dougherty D, Cromeryer C, Hoshino T, Berger M, et al. (1985) Stem cell studies of human malignant brain tumors. Part 2: Proliferation kinetics of brain-tumor cells in vitro in early-passage cultures. J Neurosurg 63: 426–432.
48. Guiot C, Pugno N, Delsanto PP, Deisboeck TS (2007) Physical aspects of cancer invasion. Phys Biol 4: P1–P6.
49. Jiao Y, Torquato S (2011) Diversity of dynamics and morphologies of invasive solid tumors. AIP Advances. in press.