We have implemented a hybrid cellular automata model based on the structure of human prostate that recapitulates key interactions in nascent tumor foci between tumor cells and adjacent stroma. Model simulations show how stochastic interactions between tumor cells and stroma may lead to a structural suppression of tumor growth, modest proliferation, or unopposed tumor growth. The model incorporates key aspects of prostate tumor progression, including transforming growth factor-β (TGF-β), matrix-degrading enzyme activity, and stromal activation. It also examines the importance of TGF-β during tumor progression and the role of stromal cell density in regulating tumor growth. The validity of one of the key predictions of the model about the effect of epithelial TGF-β production on glandular stability was tested in vivo. These experimental results confirmed the ability of the model to generate testable biological predictions in addition to providing new avenues of experimental interest. This work underscores the need for more pathologically representative models to cooperatively drive computational and biological modeling, which together could eventually lead to more accurate diagnoses and treatments of prostate cancer. [Cancer Res 2009;69(17):7111–20]
- mathematical modeling
- reactive stroma
Cellular automata (CA) models have been used to examine three-dimensional growth as well as the response to radiation treatment of brain tumors, such as glioblastoma multiforme ( 1, 2), with results that agreed well with clinical imaging data. Increasingly, CA models have been applied to more complex tumor growth and treatment problems ( 3, 4). For example, computational models are capable of both selecting effective treatments and screening out ineffective treatments for lethal melanomas and glioblastoma multiforme ( 5).
Computational models of avascular and vascular tumors have been extensively reported and reviewed previously ( 6, 7). Such models have been used to show how harsh growth conditions, such as hypoxia, can promote evolution of invasive tumors, and the results have been related in a preliminary manner to in vitro and in vivo tumor growth experiments ( 8). Tumor growth models reported in the literature also consider the role of nutrients, oxygen, necrosis at low pH, aggregation kinetics in the growth, genetic and phenotypical evolution, and morphology of the growing tumor ( 7, 9). Discrete cell simulations are typically validated by carefully controlled in vitro experiments ( 10) because they can be designed to effectively collect the vast amount of information associated with computational models.
CA models are based on discrete representation of space and time on a lattice where only neighboring lattice elements interact. One variant of the CA method is the hybrid discrete continuous CA (HCA), where continuous equations and discrete elements are solved together on the same grid. Several other tumor types have been accurately modeled using such methods. HCA models have been extensively described in the literature ( 8, 11– 14). Although most of the previously described models neglect intracellular features, they are well suited to study biological phenomena at the cellular level due to the simplicity of the logic, flexibility of the technique, and the ability to integrate multiple interacting variables across a range of spatial scales.
The prostate is a glandular sexual accessory organ composed of acinar ducts lined with luminal secretory epithelium surrounded by a layer of basal epithelial cells. These epithelial acini are encompassed by a stromal compartment composed predominantly of layers of smooth muscle (see Fig. 1A ). The principal cellular component of the prostate is illustrated in Fig. 1B.
Prostatic adenocarcinoma is the second most common cause of male cancer deaths in the Western world ( 15). Early discrimination between relatively benign lesions and highly aggressive prostate adenocarcinomas is critical for identifying those patients that require aggressive treatment while avoiding overtreating patients who would otherwise suffer no ill effects from their tumor.
In the adult prostate, paracine cross-talk between the epithelial and surrounding stromal tissues maintains homeostasis ( 16, 17). The smooth muscle and all other stromal cell types are separated from the glandular acini by a collagen- and laminin-rich basement membrane that provides positional information contributing to the maintenance of tissue architecture and differentiation through cellular signaling and structural constraint ( 18).
The loss of homeostatic interactions between organ tissues in disease has been partially attributed to a breakdown of the positional information established during development, which includes the loss of the basement membrane ( 19, 20) and an alteration of the density and type of extracellular matrix (ECM; ref. 21). This matrix is produced by an ever-expanding population of myofibroblasts ( 22). Although there is evidence that strongly implicates the role of the basement membrane and the stromal microenvironment on glandular integrity in tumor progression, the conflicting data and the vast number of factors involved in the regulation of these components limit our understanding of the multiple steps by which prostate tumors grow and invade surrounding tissues.
The effect of transforming growth factor-β (TGF-β) on cancer progression has been the target of much research and debate ( 23). TGF-β normally inhibits the proliferation of epithelia through induction of the cell cycle inhibitors p15 and p21 ( 24). The determination of whether TGF-β will induce cytostasis or apoptosis in normal epithelia depends on the intensity of their proliferative activity in addition to poorly understood microenvironmental determinants ( 23, 25). Stromal production of TGF-β by prostate carcinoma–associated fibroblasts has been shown to increase the growth and invasiveness of initiated prostate epithelia ( 26); however, the effect of epithelial TGF-β production from organized prostate glands is still unclear. Moreover, systemic inhibitors of TGF-β have yielded conflicting data in therapeutic trials ( 27), underscoring the need for further analysis of the complex roles of TGF-β at different stages of neoplastic progression.
The TGF-β family of cytokines is highly pleiotropic. Some of their functions have been accurately modeled computationally, including roles in vascular remodeling, hyperplasia, and wound repair ( 28). Continuum ordinary differential equation tissue growth models of TGF-β–mediated stromal-epithelial interactions were shown to have good qualitative agreement with experimental biological results ( 29).
Prostate tumors have previously been investigated mathematically ( 30– 34), mainly focusing on the role of prostate-specific antigen and androgen levels in progression. Although insightful, these models were purely continuous and nonspatial and did not consider the role of stroma. In this article, we model prostate tumorigenesis using a HCA model that integrates five different cell species (discrete, individuals) with three different microenvironmental chemical species (continuous, concentrations), all of which are thought to play key roles in prostate cancer. Using this model, we have investigated the importance of TGF-β in driving prostate cancer progression. In particular, we examined how TGF-β regulates tumor-stroma interactions in tissues with tumor cells possessing different degrees of organization (as measured in terms of glandular morphology). In addition, following a prediction generated by the model about the role of TGF-β in organized glandular epithelia, we have provided previously unexamined pathologic and experimental support for the ability of the model to appropriately reflect the physiologic role of TGF-β at a specific stage of tumorigenesis.
Materials and Methods
HCA models are mathematical tools in which discrete entities interact with continuous ones. In the model, cells are the discrete entities represented as points on a two-dimensional grid (2 mm × 2 mm slice of a three-dimensional prostate). This grid structure has been used in other cancer-related CA models ( 7, 35, 36) and has proved to be very efficient. This grid hosts three microenvironmental variables, which are treated as continuous concentrations: TGF-β, matrix-degrading enzyme (MDE) expression, and membrane/ECM. The discrete cells are designated as basal epithelial, luminal epithelial, motile stroma (representing bone marrow–derived cells such as macrophages), static stroma (fibromuscular stroma), and tumor cells. Figure 1B shows how cells are initially arranged: the simulated section of prostate contains three glands (Supplementary Fig. S1), each following the structure shown in Fig. 1B and arranged along an axis going from the upper-left to the lower-right corner of the domain (as shown in Supplementary Fig. S1). The space outside the glands is either empty or occupied by static and motile stromal cells. Static cells can represent muscle or fibroblastic lineages, whereas motile cells are predominantly modeled on bone marrow–derived monocyte/macrophage lineages.
Figure 2 shows flowcharts that describe the behavior of the different cell types. A detailed explanation of these flowcharts can be found in the Supplementary Appendix. The microenvironmental variables of the model are TGF-β, MDE expression, and membrane/ECM concentrations (representing both the ECM, which is present in the mesenchyme and is made of elements such as collagen, fibronectin, laminin, and vitronectin, and the epithelial basement membrane). The dynamics of these microenvironmental variables are defined by three partial differential equations that describe how each of them evolves in space and time. The dynamics of TGF-β (Tβ) are as follows: which shows that TGF-β is produced by basal and cancer cells as well as by motile stromal cells in proportion to the local TGF-β concentration. It also shows that TGF-β is consumed by luminal and motile stroma cells. TGF-β also binds to the ECM at a rate that depends on the local concentration of TGF-β and also that there is some natural decay of the ligand.
MDEs (E) are produced by tumor cells (at rate λ; ref. 37), diffuse (at rate DE), and are depleted as they degrade the ECM and the basement membrane (μ).
Basement membrane/ECM (M) is produced by basal cells (depending on the current local concentration of ECM ensuring the density never exceeds the maximum m0) and motile stroma (depending on rate αF, scaled by the local concentration of TGF-β; ref. 38). Finally, the ECM gets degraded by the MDEs at a rate μ:
These microenvironmental factors (equations above), coupled with the discrete cells (defined by the cellular life cycle flowcharts), represent an evolving dynamic hybrid model of the prostate in which tumor progression emerges from the interactions of all the key elements. Cell behavior results from the interactions with other cells in competition for space and via the microenvironmental factors as shown in Fig. 2F. It can be seen that TGF-β is a key factor that determines the survival and growth of both tumor and normal epithelial cells as well as the behavior of stromal cells. The balance of TGF-β, MDEs, and membrane production modulates the competition of the different cellular species.
Time is discrete and is measured in time steps of 24 h. At each time step, cells follow their life cycle ( Fig. 2). The timescale affecting the microenvironmental variables is assumed to be much faster (×100) than that of the cells to accommodate the different timescales affecting processes at the cellular and molecular levels.
The parameters used to characterize the model are assigned the values shown in Table 1 . Many of the values used are estimates, denoted by an asterisk ( 39– 41), which are biologically plausible but not in every case experimentally measured. Because the parameters need to be nondimensionalized, the results of the nondimensionalization (with respect to time, which is taken to have cell cycles of 24 h, and space of 2 mm × 2 mm) are also given. These parameters have the property of making the system homeostatic before cancer initiation.
Therefore, under normal circumstances (in the absence of tumor cells), there is a natural homeostatic state for all the variables with equilibrium concentrations of TGF-β, ECM, and MDE. These represent the initial conditions for each of these microenvironmental variables (Supplementary Fig. S1).
Because the composition of the stromal compartment is statistically correlated with prostate disease recurrence ( 42), the role of stroma on tumors at different stages of malignancy (determined by TGF-β and MDE production) was studied. Each of the simulations was allowed to run for 20,000 time steps, corresponding to ∼55 years in real time. After 10 time steps, six basal epithelial cells (four in the central acinus and one on each of the two remaining acini) become abnormal cells initiating tumorigenesis. The simulations test different stromal configurations in which the proportion of stroma in the domain ranges from 20% to 40% of the total space. The stromal configurations considered are (motile, nonmotile) as follows: (a) high proportion of motile and nonmotile stroma (40%, 40%), (b) high proportion of motile stroma (40%, 10%), (c) high proportion of nonmotile stroma (10%, 40%), and (d) low proportion of motile and nonmotile stroma (10%, 10%).
Sample simulations. Figure 3 (and Supplementary Movies) shows the three main outcomes resulting from the simulations. Figure 3A shows an example of a simulation in which the tumor breaks out from the acini and grows into the surrounding stroma. Initially, the tumor cells fill the lumen inside the gland ( Fig. 3A). Eventually, the basement membrane starts to degrade and TGF-β begins to leak from the gland and attract motile stroma ( Fig. 3B). At the end of the simulation, the tumor has taken a significant portion of the domain but its growth is constrained by the motile stroma responding to the TGF-β gradient ( Fig. 3C). As a result of periodic boundary conditions, tumor cells growing on the upper boundary appear on the bottom; however, this has no effect on the simulation outcome. Figure 3B shows a simulation in which tumor cells quickly break from the gland and grow unopposed. Finally, Fig. 3C shows tumor cells producing excessive quantities of MDE, which leads to early breakdown of the basement membrane and resulting leakage of the TGF-β, without which the tumor cells die.
Varying the microenvironment. To further investigate the three distinct outcomes the model can produce, we examined the effect of varying each of the three microenvironmental parameters. Due to the limitations of our parameter estimates, varying them over 2 orders of magnitude also allows us to test the robustness of the model outcomes as well as the sensitivity (see Fig. 3D). Figure 3D (I–IV) shows results.
MDE production has a dominant effect on the outcome of the simulations, and results show that too little MDE production leaves the tumor confined to the gland, whereas too much means that the basement membrane breaks down (and consequently, TGF-β spills out of the gland) too early, depriving the tumor cells of the required minimum of TGF-β level to survive ( Fig. 3). This hints that there is an optimal range of MDE production that would allow the tumor to escape from the gland and invade the surrounding tissue. This particular range is parameter dependent, modulated by the rate at which the basal membrane is degraded by tumor cells and repaired by basal and stromal cells as well as the susceptibility of the tumor cells to TGF-β concentration.
Contrary to our expectations, tumor TGF-β production does not modify significantly the simulation outcomes. Stromal configuration, on the other hand, has a more significant effect on tumor progression. Prostates with a high proportion of stromal cells are less likely to be taken over by a tumor than those with a lower density of stromal cells. The role of stroma in tumor promotion becomes critical after the tumor has escaped from the gland, and it is also at this point when TGF-β, which mediates the interactions between stroma and tumor, plays a more dominant role. In fact, tumor cell TGF-β production from within contained glands seems to inhibit tumor progression due to the recruitment of matrix-producing motile stroma. This often leads to an increased deposition of stromal-produced ECM (see Fig. 3A and compare Supplementary Figs. S2 and S3 with different amounts of stromal cells). To ground this in pathologic reality, human prostatic intraepithelial neoplasia (PIN) biopsies were immunostained for TGF-β and counterstained with Masson's trichrome to determine whether an adjacent stromal matrix deposition could be observed. Figure 4 shows very little to no TGF-β immunoreactivity in normal human prostate glandular epithelium ( Fig. 4A) compared with high immunoreactivity in some glandular PIN foci ( Fig. 4B). Furthermore, trichrome staining of these sections revealed that the stroma adjacent TGF-β–expressing PIN foci had increased collagen production, confirming the HCA model outcome.
Given the potential implications, it was important to determine whether the prediction that tumor cell TGF-β production does not significantly alter the simulation outcome was biologically realistic. To test this, constitutively active TGF-β1 was overexpressed in a human prostate epithelial cell line previously shown to be capable of malignant transformation (BPH1; refs. 43, 44). Importantly, this cell line is able to form acini when recombined with inductive rat urogenital mesenchyme (rUGM; ref. 44), allowing us to determine whether the expression of TGF-β might result in enhanced tumorigenicity.
A hemagglutinin (HA)-tagged, constitutively active TGF-β1 previously generated and characterized ( 45) was expressed in BPH1 cells (see Materials and Methods; Fig. 5A ). To confirm its functionality, conditioned media collected from parental BPH1 and BPH1TGFβ1 cells were placed on serum-starved human prostate fibroblasts expressing a red fluorescent protein (RFP) TGF-β reporter. Results show an activation of TGF-β signaling by BPH1TGFβ1 cell conditioned media, but not the parental control ( Fig. 5B). To test our simulation prediction, we recombined each epithelial cell line with rUGM for 6 weeks under the kidney capsule of severe combined immunodeficient mice, a technique previously applied to the study of prostate development and tumor formation ( 44). Kidney capsule grafting of tissue recombination with inductive rUGM results shows that TGF-β1 production decreased the overall size of the tissue recombination graft by 3- to 4-fold ( Fig. 5C). Histologically, trichrome ( Fig. 5D) and smooth muscle ( Fig. 5E) staining of these grafts revealed that BPH1TGFβ1 grafts show a decrease in smooth muscle and an increase in periacinar collagen production over controls, consistent with changes seen around premalignant human lesions ( Fig. 4). Moreover, histochemical identification of basement membrane by periodic acid-Schiff (PAS) staining shows that TGF-β1 production by initiated prostate epithelia does not affect glandular containment ( Fig. 6 ).
This article introduces a mathematical model that recapitulates some key elements of tumor progression, such as epithelial invasion and stromal reorganization. The aim of the model is to provide a basic platform from which to study the role of TGF-β and stroma in prostate tumor progression from nontumorigenic homeostasis to preneoplasia to malignancy. To be useful, mathematical models must generate hypotheses that can be tested in biological systems and compared with clinical observations and, further, must be amenable to testing and generating a deeper understanding of biological phenomena.
PIN was first defined 20 years ago and is considered to be the most common premalignant lesion in prostate cancer ( 46). PIN can be made to regress in humans ( 47) and this process can also be documented in experimental animal models ( 48); however, it is unclear if such regression happens spontaneously in patients. Thus, although PIN can progress to cancer, this is not an inevitable progression and the conditions that promote or suppress such a progression are presently unknown.
Under basal conditions, the model presented here is homeostatic. The disruption of homeostasis as a result of the appearance of tumor cells can result in one of four different scenarios depending on the phenotype of the tumor cells as well as the stromal microenvironment. These scenarios ( Fig. 3) are as follows: (a) nascent tumors that develop slowly and remain at a stage that recapitulates aspects of the premalignant lesion PIN, (b) tumors that break out of their acinar containment and invade locally but do not outgrow the other cell types, (c) tumors that develop rapidly and take over the prostate ( Fig. 4B), and (d) tumors that develop quickly but are self-limiting, dying out before progression beyond PIN ( Fig. 4C). These results show both the potential of the model to capture distinct outcomes and the reproduction of the diversity of timescales at which different prostate tumors are known to develop.
This mathematical model generated several hypotheses. In our simulated prostate environment, before tumor breakout, TGF-β production by tumor cells seems to have little effect on the outcome of the progression (see Fig. 2), which was corroborated experimentally (see Fig. 5). TGF-β in this model does recruit motile stromal cells, which structurally inhibit glandular breakdown through matrix production, which was also corroborated by examination of clinical samples of human PIN foci expressing TGF-β (see Fig. 4). This ineffectiveness persists regardless of either the proportion of stromal subtypes or the rates of production of MDEs by the tumor cells. Similar to biological assumptions ( 23, 49), the effects of TGF-β in this model are more likely to be tumorigenic once tumor cells emerge from a contained PIN-like state as determined by the effects of the TGF-β gradient on tumor cell growth and stromal recruitment ( Fig. 3A). This is consistent with previously described biological observations ( 26).
A second hypothesis generated by the model concerns the role of tumor cell MDE production. The results show that there is a critical range of values of MDE production that are optimal for tumor progression. If tumor cells produce MDEs at a suboptimal rate, the time to break out and even the chances of the tumor progressing from PIN are affected negatively ( Fig. 3). Alternatively, excessive production of MDEs might work against the tumor and simulation results suggest that too much MDE production could, under certain conditions, lead to tumor cell death. Degrading the basement membrane too quickly may result in TGF-β diffusing out of the glandular acinus, leaving the tumor cells in the gland with insufficient growth factors to survive (Supplementary Fig. S4). Therefore, it is a prediction of this model that the appropriate resources (e.g., TGF-β) essential for invasion and survival of tumor cells outside of the gland are available once the basement membrane is breached. The specific values that characterize this critical range of MDE production will depend on the rest of the parameters of the model, such as TGF-β production and diffusion, or the rate at which MDE is used to degrade the basal membrane.
There are several implications coming from these hypotheses. The first suggests that, despite the central role that TGF-β played in the conception of the model ( Fig. 2F), the level of production by tumor cells does not affect the initial premalignant stages of tumorigenesis. This hypothesis does not address the role of TGF-β once the tumor progresses beyond containment, where the role of TGF-β seems to be more protumorigenic. A different hypothesis, which could have far-reaching implications, is the role of excessive MDE production by tumor cells in halting tumor progression, which, at present, is untested.
Finally, an interesting insight provided by the model is the dual role of structural constraints in tumor progression. The glandular structure of the prostate both delays and protects tumor growth. Although this arrangement limits the growth in the early stages, it provides a haven in which tumor cells can grow until they reach a critical mass. The role of the stromal configuration has also been shown to be important before and after the tumor breaks out from the gland. In those cases in which the proportion of stroma is high, degradation of the basement membrane is hindered by the presence of motile stroma. Furthermore, when the membrane is breached, the stroma helps to constrain tumor expansion. Thus, the model hints that the abundance of stromal cells, both motile (that block tumor growth and produce ECM) and static (that amplify the TGF-β signal and provide a clearer gradient for motile stroma), hinders tumor growth.
A key value of designing biologically relevant mathematical models is the clarification of deficiencies in biological modeling. At present, there is no valid biological model of the effect of stromal heterogeneity on human prostate glandular stability or tumor progression. As such, adequate experimental validation and model parameterization of the effects of variable stromal compositions under different MDE conditions remains unexamined. However, attempts at developing such a model are ongoing, which will provide a biological tool for examining both mathematical and pathologic predictions about the effects of stromal heterogeneity on tumor progression.
To obtain a tractable model, simplifications were unavoidable. For example, all tumor cells share the same phenotype. Nonetheless, tumors are microcosms of evolution in which several tumor phenotypes coexist as a result of microenvironmental heterogeneity and a complex evolutionary process. Furthermore, it is known that some types of stromal cells, such as fibroblasts, can alter their phenotype in response to microenvironmental changes produced by a tumor. It is likely that these evolutionary dynamics would affect tumor progression. Therefore, it is our goal to explore the effect of the stromal-tumor coevolution in prostate tumor progression.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Grant support: NIH/National Cancer Institute grants of the Microenvironmental Influences in Cancer (2T32 CA009592-21A1), TMEN (1U54 CA126505 and 1U54 CA 126568), and ICBP (5U54 CA113007) programs.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
We thank the VICBC workshop, in particular Lourdes Estrada, for facilitating and motivating this integrative cancer research project.
Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).
- Received October 23, 2008.
- Revision received May 19, 2009.
- Accepted June 17, 2009.
- ©2009 American Association for Cancer Research.