An Integrated Computational Model of the Bone Microenvironment in Bone-Metastatic Prostate Cancer

Prostate cancer frequently metastasizes to bone with approximately 90% of the men displaying evidence of skeletal lesions upon autopsy (1). Despite medical advances, prostate to bone metastases remain incurable with treatments being mainly palliative (2). Advances in our knowledge of the molecular mechanisms underlying the disease should provide therapeutic opportunities to improve overall survival rates but on a more microenvironmental scale, predicting how putative therapies will impact multiple cellular responses remains a challenge. However, integrating key biologic findings with the power of computational modeling offers a unique opportunity to assess the impact of putative therapies on the progression of prostate cancer.

Understanding the normal basic multicellular unit (BMU) bone remodeling process is critical for the generation of a robust computational model (3). The initiation of the BMU by local or systemic signals results in retraction of osteoblasts from the bone surface and the formation of a canopy. Local mesenchymal stromal cells (MSC) generate RANKL-expressing osteoblasts precursors that subsequently facilitate osteoclast recruitment, maturation, and bone resorption. Degradation of the bone results in the release of sequestered growth factors such as TGF-β that in turn serve to control the extent of bone degradation and osteoblast expansion. After osteoclast apoptosis, osteoblasts rebuild the bone with a portion terminally differentiating into osteocytes and the remainder reconstituting the bone lining, ready for the next remodeling cycle.

In the metastatic prostate cancer-bone microenvironment, prostate cancer cells perturb the balance of the BMU to generate a “vicious cycle” via the expression of factors such as RANKL thereby inducing excessive bone resorption (4). The release of sequestered growth factors from the bone matrix such as TGF-β in turn stimulates the survival and growth of the metastatic prostate cancer cells. Of note, prostate to bone metastases are also characterized by areas of extensive bone formation/osteogenesis, a phenomenon that is mediated by factors such as, endothelins and bone morphogenetic proteins, BMP-2 and BMP-4 (5).

To date, the majority of our knowledge of the mechanisms driving prostate cancer progression has been garnered by focusing on the role of individual cancer/host-derived molecules. However, this molecular reductionist approach often does not take into account the multiple cellular effects of molecular mechanisms being investigated. Mathematical and computational models are a powerful means with which to study complex in vivo interactions. Numerous advances using this approach have identified roles for tumor heterogeneity in cancer progression and evolution, accurately predicting glioma progression and response to disease and, describing the cellular dynamics of bone remodeling (6). To understand the simultaneous and multiple interactions occurring over time in the metastatic prostate cancer bone microenvironment, we have generated a hybrid cellular automaton (HCA)-based integrated computational model. Using this approach, we tested the model's response to therapies currently used in the clinic to treat prostate to bone metastases (7–9). We posit that computational models will be a powerful means with which to test the efficacy of available or putative therapies for the treatment of prostate to bone metastases.

Parameters for the HCA model were derived from empirical and published data (Table 1; ref. 9). The model is composed of a grid (200 × 50 points) representing 2 × 0.5 mm² of the bone microenvironment. A major advantage of the HCA is in its intimate interconnection with experimental data, where the model and the experiments inform each other. This increases the accuracy of the model abstractions and connectivity of the basic elements, which yields reliable and biologically relevant emergent behaviors. To model the normal sequence of the BMU program, we have focused on understanding the role and behavior of the key regulators of the BMU dynamics. The principal cellular players, bone, MSCs, precursor osteoblasts (pOB), osteoblasts, precursor osteoclasts (pOC), and osteoclasts, have been explicitly modelled as agents in a grid following specific rule sets in a physical microenvironment (see Quick guide to assumptions and equations). Collectively, these components find a natural homeostatic balance that recapitulate the dynamics of bone remodeling, a homeostasis that can be perturbed via the introduction of metastatic prostate cancer cells.

Human prostate to bone metastases samples (5 μm), provided by Dr. Robert Vessella (University of Washington, Seattle, WA), were rehydrated and blocked before the addition of phospho-specific anti-Smad2 (1:200 dilution; Millipore) and pan anti-cytokeratin (1:200 dilution; Sigma) and appropriate immunoglobulin G (IgG) controls. Tissue sections were incubated overnight at 4°C. Subsequently, species-specific secondary AlexaFluor 568 and AlexaFluor 488-conjugated antibodies (1: 1000 dilution for one hour at room temperature; Invitrogen) were added for imaging by microscopy. For semiquantitative analysis, regional images were segmented on the basis of the intensity of staining using Definions Tissue Studio (TS).

Osteoblast (MC3T3) and MSC migration was assessed using a modified Boyden Chamber assay. Cells (5 × 10⁵) were seeded in the upper chamber and their migration to TGF-β rich PAIII conditioned media in the presence of a TGF-β inhibitor (1D11; Genzyme) or isotype control (13C4) at a concentration of 5 μg/mL was measured over a 5-hour period at 37°C. Migrated cells were stained with hematoxylin and air dried. The number of migrated cells was counted in three random ×20 fields for each condition. All experiments were performed in triplicate.

All animal experiments were done with University of South Florida (Tampa, FL) Institutional Animal Care and Use Committee approval (CCL; #R3886). PAIII luciferase expressing cells were injected (1 × 10⁵ in 10 μL volume) into the tibia of anesthetized immunocompromised mice (recombinase activating gene-2 null; RAG-2^−/−; ref. 20). Tumor growth was imaged using bioluminescence imaging and quantitated with IVIS Living Image software. After 2 weeks, tumor-bearing tibias were excised, soft tissue removed and processed for histology and histomorphometry as we have previously described (19, 21).

For cell culture and statistical methods, please see Supplementary Materials and Methods.

Key for the generation of the BMU was the incorporation of major cellular players namely the prostate cancer metastases, MSCs, pOBs, osteoblasts, pOCs, and osteoclasts. Furthermore, based on the literature, RANKL, TGFβ, and BDFs were considered main mediators driving the cellular interactions (14). Although RANKL has been well characterized in the context of prostate to bone metastases, few studies have explored the role of TGF-β signaling in this setting. Therefore, we initially examined the activity of TGF-β signaling in prostate to bone metastases in human specimens. Using phosphorylated SMAD2 (pSMAD2) as a readout for TGF-β receptor activity, our results indicated that TGF-β signaling is active in human prostate to bone metastases (Fig. 2A). Quantitative analysis revealed that in the human specimens, TGF-β signaling was highest in the prostate cancer cells but we also observed strong staining in stromal cells, including osteoblasts (Fig. 2B). We also examined the effects of TGF-β on the proliferation and migration of the cellular components of the vicious cycle including prostate cancer cells, osteoclast precursors, and osteoblast precursors, including MSCs. Of note, we observed that TGF-β significantly impacted the migration of MSCs and the osteoblast precursor cell line, MC3T3-E1, suggesting a role of the growth factor in the recruitment of cells that could impact prostate cancer-induced osteogenesis (Fig. 2C and D).

On the basis of these empirical data and the literature, we parameterized a HCA computational model of the homeostatic BMU (Fig. 3A and Supplementary Fig. S1 and Movie 1). On the basis of multiple simulations (N = 25) of the computational model, we observed little variation in each population cell number. In some instances (N = 2), the BMU failed to initiate, in part, due to spatial and cytokine gradient differences between the different simulations (data not shown) but, we expect that persistent remodeling initiation stimuli would eventually lead to the formation of the BMU in vivo. We also observed that in a subset of BMU simulations (N = 2), simultaneous osteoclast fusion resulted in two sites of resorption. However, the generation of BDFs by the osteoclasts sufficiently increased osteoblast numbers and returned the BMU to baseline (data not shown). Importantly, the typical interactions between the different elements of the computational model result in a homeostatic BMU. It is important to note that each of the cells behave as autonomous agents and possess the inherent ability to respond to the surrounding environment independently.

Next, we introduced a single metastatic prostate cancer cell expressing TGF-β ligand and receptors into the BMU. Interestingly, we observed that in many of the simulations (N = 18 of 25), the metastases failed to generate a lesion. We anticipate that the introduction of prostate cancer emboli would significantly enhance the “take rate.” This “take rate” result emphasizes the stochastic nature of the model and reflects the in vivo reality where not every metastatic cancer cell that successfully invades a BMU would result in a viable lesion. In the simulations where lesions were initiated (N = 7 of 25), we observed that after 10 days, the presence of the prostate cancer cells resulted in the integrity of the canopy being compromised and a resultant increase in osteoclast recruitment and maturation (Fig. 3B). Over time, the prostate cancer cells expanded, resulting in MSC infiltration and osteoblast-mediated bone formation ultimately recapitulating the “vicious cycle” paradigm (Fig. 3B and Supplementary Fig. S2 and Movie 2).

Analysis of the computational prostate cancer-bone microenvironment at day 240 revealed striking histologic similarities to an in vivo model of the osteogenic/osteolytic prostate to bone metastasis environment (Fig. 4A and B). We noted that the number of prostate cancer cells varied amongst simulations (8,625 ± 4,580; N = 5) but in general, the growth rate predicted by the computational model was comparable with the growth rate of the prostate cancer cells in vivo (Supplementary Fig. S3). We also noted similar proportions of stromal cell populations at the computational and biologic study endpoints (Fig. 4C and D). In the computational model, distinct phases of cell activity were observed. For example, the numbers of adult osteoblasts increased over time but notable plateaus before increases in cell number existed (Fig. 4C). These plateaus in adult osteoblasts corresponded with dips in the pOB and pOC precursor population. MSC numbers, however, gradually increased over time and in general, paralleled increases in cancer cell number. In fact, the model predicts that MSCs are crucial for the progression of the metastases because relaxing the probability of recruitment greatly impacts the growth of the metastases (Supplementary Fig. S4). Furthermore, we observed that osteoclasts were critical for cancer progression in the model with numbers changing over time from zero to 12 (Fig. 4E). Although the in vivo model output has a similar number of osteoclasts per field (Fig. 4D and F) at the study endpoint, the phasic nature of osteoclast involvement is not apparent.

To test the applicability of the model in treating prostate to bone metastases, we applied two standard-of-care treatments, namely, a bisphosphonate and an anti-RANKL inhibitor, that induce osteoclast apoptosis during resorption and inhibit osteoclastogenesis, respectively. To mimic the clinical scenario, we applied the bisphosphonate at a time where metastases had established (day 80), although it should be noted that therapies could be applied at any juncture to the model (Fig. 5A and B and D; N = 5 simulations per group). During bisphosphonate treatment, osteoclasts still formed but typically died within 24 hours of initiating bone resorption (Fig. 5B, bottom). However, residual bone resorption during bisphosphonate therapy was sufficient to sustain the metastases with the number of cancer cells on average at the final time point being 7,138 ± 2,343 (n = 5) compared with nontreated control of 8,624 ± 4,580 (n = 5; P > 0.05). Interestingly, the application of an anti-RANKL inhibitor halted cancer growth with no cancer cells detectable after 20 days of administration (Fig. 5C and D; P < 0.05). These data suggest that treatment with anti-RANKL inhibitors should be curative in the clinical setting. However, in clinical trials with anti-RANKL inhibitors such as Xgeva, there is a slight but not significant increase in survival of patients on the therapy compared with bisphosphonates (2). In our simulations, we found that reducing the efficacy of the anti-RANKL inhibitor from 100% to 40% resulted in an output of cancer cells (3,157 ± 3,037; N = 5) that was comparable with that of the bisphosphonate-treated group (7,138 ± 2,343 cancer cells; N = 5; Fig. 5B and Supplementary Fig. S5). These data suggest that improved efficacy of anti-RANKL delivery into the prostate cancer bone microenvironment could be curative.

In the current study, we have generated a faithful computational model of the BMU. It is important to note that the homeostatic behavior is not hardcoded but emerges from the interactions between the different primary cell types of the bone in response to TGF-β, RANKL, and BDFs. Furthermore, informed by experimental evidence, the introduction of a simulated TGF-β ligand and receptor expressing prostate cancer cell into the BMU resulted in a vicious cycle that yielded mixed osteogenic/osteolytic lesions over clinically relevant periods of time. Key findings arising from the computational model include: (i) the ability to assess temporal changes in cellular populations and dissect complex dynamics that are difficult to determine in vivo, (ii) the phasic osteolytic/osteogenic nature of the metastases, (iii) the application of clinically used therapies such as bisphosphonates and anti-RANKL therapies illustrate the usefulness of the model in predicting the efficacy of targeted inhibitors, and (iv) the impact of inhibitors at varying doses on the progression of prostate to bone metastases.

In our study, we selected the BMU as the primary target for metastatic prostate cancer cells in which to establish and grow. This is logical because reports have shown that high rates of bone turnover correlate with poorer prognoses for patients with prostate cancer and that the metastases can utilize many of the bioavailable factors present in the remodeling environment (22–24). On the basis of this rationale, we introduced a single metastatic prostate cancer cell into the BMU that over time recapitulated the vicious cycle paradigm. The “histologic” output of the computational model at the endpoint was remarkably similar both in morphology and in cellular population proportion to an in vivo model used by our group to study mixed prostate to bone metastases. In contrast, however, the computational model illustrates that the cellular population in the prostate cancer-bone microenvironment is dynamic and changes occur both exponentially (tumor growth) and in phases (osteoclast, osteoblasts). For example, our model illustrates that osteoclast numbers rarely exceed a total of 20 out of approximately 2 × 10⁴ cells per computational field of view (Fig. 4C and E). Reports have shown that high levels of TGF-β can hinder osteoclastogenesis therefore, limiting the number of osteoclasts that can form in the tumor-bone microenvironment even in the presence of RANKL-producing osteoblast precursors (25). In fact, our results show that osteoclast-mediated bone resorption is critical for the induction of the osteogenic metastases, an observation that supports the use of anti-RANKL therapies in men with osteogenic prostate to bone metastases.

In a number of model iterations, we observed that the recruitment of MSCs to the prostate cancer-bone microenvironment was essential for the generation of osteoblast precursors and the development of osteogenic lesions. Our in vitro data also show that TGF-β contributes to MSCs and osteoblast precursor migration therefore, providing a means through which these cell types can be recruited to areas of prostate to bone metastases and contribute to their progression (Fig. 2). Although the role of MSCs in the prostate cancer-bone microenvironment has not been explicitly explored thus far, our model predicts an important role for this cell type in tumor growth. It is important to note that we did not consider the trans-differentiation of prostate cancer cells into osteoblasts, but “osteomimicry” is a distinct possibility that could also be integrated into future iterations of the model (26).

The ability to dissect changes in cellular composition in the computational model provides key insights into how the cancer cells, MSCs, osteoblast precursors, osteoblasts, osteoclast precursors, and osteoclasts are interacting with each other over a clinically relevant period of more than 200 days. In vivo models are typically analyzed at endpoint or at predetermined time steps to assess cancer growth and changes in the microenvironment in control and test groups. For example, the MDA-MB-231 progression in bone is often measured at weekly time points (27). The computational model generated in this study clearly illustrates that a number of the host microenvironment components, notably the mature osteoclasts and osteoblasts, undergo phases of activity and rest (Figs. 3 and 4). This level of resolution is not available in biologic models but demonstrates that the time points or endpoints chosen for in vivo models are a “snapshot” that may not be truly reflective of what has, or is about to happen in the tumor-bone microenvironment. Knowledge of the dynamic changes occurring over time in the computational cancer bone-microenvironment could lead to a better understanding of when to apply inhibitors or what happens to the cell populations over time once inhibitors have been applied.

Bisphosphonates and more recently anti-RANKL therapies are used as treatment strategies to protect patients with prostate to bone metastases from skeletal-related events (SRE). Studies have shown that bisphosphonates can extend the average time to SRE for patients, and anti-RANKL–based therapies are significantly better than bisphosphonates in extending that time to first SRE (2, 28). However, neither treatment increases overall survival. We applied these inhibitors to our computational model. Assuming an efficacy of 100%, application of a bisphosphonate after a period in which the metastasis is actively growing (day 80) demonstrated an impact on osteoclast activity and on tumor growth. This time point was chosen based on the prostate metastases being established and actively growing but therapeutics could be applied at any juncture, a useful feature to study the response of multiple bone metastases at various stages of progression. Subsequent to the application of bisphosphonates to the model at day 80, we observed that osteoclasts still formed and that the residual production of TGF-β and BDFs was sufficient to sustain tumor growth, albeit to a lesser extent compared with the nontreatment control arm (Fig. 5). Significantly, we also assumed that the dosing of bisphosphonates was at 100% efficacy over the course of the therapy, but, in reality, there are most likely changes in concentrations. Gradients in therapy concentration can be accounted for and pharmacokinetics can be explicitly simulated in the computational model and it will be interesting to determine the impact of dosing gradients on the behavior of the cell populations over time in the computational microenvironment.

Application of an anti-RANKL inhibitor, again at a presumed efficacy of 100%, significantly impacted the prostate tumor-bone microenvironment by preventing osteoclast formation, and subsequently, tumor growth. The model therefore, predicts that anti-RANKL inhibitors should be curative in the clinical setting. However, the clinical reality is that anti-RANKL inhibitors do not significantly extend overall survival in men with metastatic prostate cancer (2). Interestingly, reducing the efficiency of the anti-RANKL therapy to 40% closely mimics that of the bisphosphonate treatment and suggests that increasing doses or better targeting of the anti-RANKL inhibitors to the bone could enhance the efficacy of the drug provided that there are no or minimal increases in noted side effects such as osteonecrosis of the jaw or cataract formation (Supplementary Fig. S3). Clinically, the dosage and frequency of administration of anti-RANKL therapy such as Xgeva are based on trials that demonstrated the optimal balance of efficacy, as determined by a more than 70% decrease in urine N-terminal collagen fragments (NTX), and tolerability was 120 mg subcutaneously delivered every 4 weeks (29). It is possible that higher doses may prove more efficient in significantly enhancing overall survival in patients with prostate to bone metastases as suggested by the computational model but this increase in dose may be outweighed by increased risks of side effects. A major advantage of the computational model is the application of combination or putative therapies to study cellular behavior in the prostate bone-microenvironment over time. For example, our results highlight the role of active TGF-β signaling in the cancer and host cells of human prostate to bone metastases and in the migration of MSCs and osteoblasts (Fig. 2). This observation is in keeping with other studies and underscores the key role TGF-β signaling plays in the bone microenvironment in regards to promoting the progression of prostate to bone metastases. The computational model could therefore easily test the efficacy of TGF-β inhibitors applied to the prostate cancer-bone microenvironment and predictions used to inform preclinical and ultimately clinical trials.

There are a number of caveats to the computational model described herein. Quantitative predictions from computational models are typically dependent on the information used to parameterize it. The key values used to parameterize the computational model presented in this paper are based on TGF-β and RANKL. The flexibility of the model ensures that ranges in the concentration and balance of other factors and cells that can impact the BMU in the normal and disease setting can be easily integrated. Enhancing these qualities will improve the accuracy of the generated predictions, but our existing model is already quite robust to changes in the parameterization (Supplementary Fig. S4). Although our model is relatively sophisticated, simpler less computationally intensive mathematical models have a number of advantages in that they are easier to understand and analyze. Furthermore, having fewer parameters, they are amenable to the fitting of existing experimental data using techniques such as approximate Bayesian computation (ABC). Fitting extant experimental data to mathematical models has been used successfully to drive a model to closely represent a set of observations in cases such as imatinib response in patients with leukemia (30). These approaches can also provide exact results when sufficient summaries are used (31). Simpler mathematical models are usually preferable, especially when values for biologic parameters are largely unknown (32, 33). In this case, a more complex model was necessary to capture the bone homeostasis emerging from the interactions between MSCs, osteoclasts, osteoblasts, pOCs, and pOBs. The existence of reliable empirical data made it possible to properly parameterize such model. This has allowed us to elucidate the mechanisms involved in the vicious cycle of prostate to bone metastases. One drawback of more complex computational models is that they are computationally intensive and therefore, limiting the ability to perform large numbers of simulations that statistically can explore the robustness of each of the chosen parameters. However, in multiple simulations (N = 25), we have ensured that the plausible plasticity of these parameters fall within the currently accepted biologic empirical values. As such, we have focused on investigating changes on the parameters that may vary experimentally; especially molecules such as RANKL and TGF-β that are difficult to measure in vivo that have yielded interesting insights (Fig. 5 and Supplementary Materials and Methods). Follow-up experiments will use increased numbers of simulations to enhance statistical analysis and use new biologic parameters being empirically determined to improve the accuracy of the model predictions.

In conclusion, using empirical and published data, we have generated a hybrid discreet model of the BMU and shown that the introduction of single active metastatic prostate cancer cell into the BMU is sufficient to generate osteogenic lesions that are similar in pathophysiology to those in an animal model of the disease. Furthermore, the application of existing clinical therapies to the computational model underscores the value of this new approach for testing the impact of combining available therapies or putative targeted therapies for the treatment of prostate to bone metastases. Clinically, the versatility of the equations used to build the computational model ensures that it can be quickly individualized and be a powerful tool for the delivery of precision medicine to better treat and cure men with prostate to bone metastases.

No potential conflicts of interest were disclosed.

Conception and design: A. Araujo, L. Cook, C.C. Lynch, D. Basanta

Development of methodology: A. Araujo, L. Cook, C.C. Lynch, D. Basanta

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): L. Cook

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): A. Araujo, L. Cook, C.C. Lynch, D. Basanta

Writing, review, and/or revision of the manuscript: A. Araujo, L. Cook, C.C. Lynch, D. Basanta

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A. Araujo, D. Basanta

Study supervision: A. Araujo, L. Cook, C.C. Lynch, D. Basanta

The authors thank Scott Lonning at Genzyme for the 1D11 inhibitor and Lizzie Atomi Pamen for her technical support.

This work was supported by an NCI U01 (U01CA151924-01A1; D. Basanta), the Department of Defense Hypothesis Development Award (W81XWH-12-1-0016; C.C. Lynch and D. Basanta), the National Cancer Institute (RO1CA143094; C.C. Lynch), and the American Cancer Society Post-Doctoral Fellowship (PF-13-175-01-CSM; L.M. Cook). This work was also supported, in part, by the Moffitt Analytical Microscopy Core facility at the Moffitt Cancer Center.

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.