Introduction
Prostate cancer (PCa) is the second most common type of cancer in men, with over 160,000 men diagnosed during 2017 in the US alone, of which over 26,000 succumbed[ 1 ]. The majority of patients that die of the disease does so because of distant metastases, 90% of which are to the bone. For this reason, a better understanding of tumor-bone interactions is key if we are to improve how bone metastatic prostate cancer patients are treated.
Improvements in our understanding of the molecular mechanisms involved in this disease have led to the discovery of new therapeutic targets, but it remains lethal and resistance emerges in virtually every patient. Thus understanding the impact of treatments in a complex heterogeneous tumor requires approaches that can incorporate several scales of biological insights. A successful approach needs to recapitulate the emergent process of metastatic prostate cancer establishment and the emergence of resistance to treatment. In the past we demonstrated how a sophisticated agent-based computational model could help us understand the role of PCa-host interactions, as well as the impact of existing and novel treatments [ 2 , 3 ]. While agent-based models can accurately and quantitatively recapitulate cell dynamics in a specific area of a tumor, capturing the relevant intra-tumor heterogeneity and sensitivity to treatments can be a challenging and time-consuming process [ 4 ]. Alternatively, non-spatial population models that capture intra-tumor heterogeneity can be used to represent the entire metastatic burden of the disease and thus could be used as proxies for the patient in optimization algorithms [ 5 ]. Further, non-spatial models lend themselves to easier measurement and operationalization [ 6 ], and avoid the confusion of population-level intuitions for reductive ground truth [ 7 ].
Simple qualitative models can be useful in providing an understanding of how certain interactions can shape evolution and resistance in cancer. Previously, Ryser et al. used a simple ordinary differential equations (ODEs) model to understand bone remodeling [ 8 , 9 ], and to illustrate how tumor cell interactions with bone remodeling cells can shape cancer progression [ 10 ]. Evolutionary game theory (EGT) is a particularly powerful, yet simple and qualitative, approch to focusing on the role of interactions in cancer. It originated in looking at the effects of heterogeneity [ 11 , 12 ] and we have used the it to look at a wide number of dynamics in cancer like go-vs-grow [ 13 , 14 ], Warburg effect [ 15 ], tumour-stroma interaction [ 16 ], and interaction of multiple public goods [ 17 ]. In the process, we have built EGT models of many cancers, including prostate cancer [ 15 , 16 ] – as have several other groups. Dingli et al. [ 18 ] focused on a type of bone cancer, multiple myeloma, using an EGT approach that incorporated both tumor cells as well as bone stroma cells. More recently, West et al. [ 19 ] studied the impact of novel treatments like adaptive therapies in the context of bone metastatic prostate cancer but without including the cellular species that characterize the bone ecosystem.
In this work, we aim to model not only cancerous growth in the bone but capture key aspects of normal bone homeostasis and its co-option. We built an EGT model where PCa cells can co-opt cells from the Bone Modeling Unit (BMU). In brief, the BMU is an area of trabecular bone that is being remodeled via osteoclasts (OCs) and osteoblasts (OBs), which resorb and deposit bone tissue respectively, and the interplay between these and other BMU cell types via molecules such as receptor activator of nuclear-factor kappa-B ligand (RANK-L), which is generated by osteoblasts and promotes osteoclast differentiation and survival, and TGF- β , which is released when bone tissue is removed by Osteoclasts. These cellular species make the bone a very dynamic organ [ 20 ]. Furthermore, studies have shown that metastatic cancer cells can co-opt this process for their own benefit [ 21 ].
Capturing this aspect of regulation and dis-regulation of bone homeostasis allows us to better understand the selection that drives metastatic cancer evolutionary dynamics and the impact of treatments on the bone microenvironment – a topic of clinical significance. Furthermore, we also incorporate ecological aspects by modeling the role of the bone in regulating the fitness of the different cell types [ 22 ]. Our results suggest that we can optimize the combination of treatment on and off periods to limit tumor growth and control bone growth.
Model
The interactions between osteoclasts (OC) and osteoblasts (OB) orchestrate bone remodeling and homeostasis. As such, bone volume is determined by the balance between the density of OBs ( ρ OB ) and density of OCs ( ρ OC ) on the existing bone, where OBs increase the bone volume and OCs decrease the bone volume according to the following discrete time dynamics: (1) Δ B = 1 K ( ρ O B ( 2 − B ) ︷ Bounded Eect of OB − ρ O C B ︸ Bounded Eect of OC ) where OBs and OCs are in the range B ( t ) ∈ [0, 2] and K is a constant for determining time-scales.
Note that the above dynamics are physical and not evolutionary. The dynamics of densities for OB, OC, and the sensitive and resistant tumours ( ρ T , ρ T R ) however are ecological and given by logistic growth: (2) Δ ρ i = ρ i ( ( 1 − ρ ) σ W i − α ) where i ∈ {OB, OC, T , T R } and ρ = Σ i ρ i is the total density (thus 1 – ρ is the amount of remaining space in the niche), σ is the selection strength, α is the death rate, and W i are the fitness functions.
The fitness function for OB is given by: (3) W O B = ( ρ OC ρ OC + ρ OB ( 1 − B ) + ( ρ T + ρ T R ) δ ) ( 1 − 2 s ) where the first summand of W OB corresponds to the benefit conferred to OCs by OBs that can be attributed to the secretion of RANKL by OBs that results in the recruitment of OCs to the remodeling site. This is proportional to the proportion of OCs among healthy cells ( p OC = ρ OC ρ OC + ρ OB ) and lack of bone 1 − B . In the second summand, δ is the benefit that an OB receives from interacting with any PCa cell as the PCa cells secrete modest amounts of TGF- β .
The fitness function for OC is given by: (4) W O C = ρ OC ρ OC + ρ OB ( B − 1 ) ( 1 − 2 s ) which captures the benefit conferred to OB by OC that can be attributed to the secretion of TGF- β as a result of bone resoption. This is proportional to the proportion of OBs among healthy cells ( p OB = ρ OB ρ OC + ρ OB = 1 − p oc ) and overabundance of bone B −1. Since both W OB and W OC are functions of proportions of OC and OB, we could rewrite their coupling as a replicator dynamic. See Kaznatcheev [ 23 ] for an example of this in a similar model.
Finally, the fitness functions for the chemotherapy sensitive and resistant tumour are: (5) W T = ( ρ O C γ + ( ρ T + ρ T R ) ϵ ) ( 1 − 2 s ) (6) W T R = ρ O C γ + ( ρ T + ρ T R ) ϵ − r where γ is the benefit a PCa cell receives from interacting with OCs. The OC-led resorption of the bone allows the neighboring PCa cells to access nutrients and growth factors previously embedded in the bone; ϵ is the benefit that a PCa cell receives from interacting with other PCa cells derived from the secretion of TGF- β . The cost of resistance to treatment is r . Resistance to chemotherapy is common and results from the treatment providing strong selection for PCa cells that can avoid its cytotoxicity. This resistance often comes through the upregulation of drug exporter pumps on the surface of PCa cells. Producing and maintaining these pumps is energetically costly to the PCa cells but allows those that have a sufficient number of them to deal with cytotoxic drugs.
Chemotherapy.
For all the fitness functions above, chemotherapy is implemented by the introduction of a cost s , which corresponds to the the strength of the chemotherapy regiment. Chemotherapy is widely used in treating hormone-sensitive and insensitive metastatic prostate cancer [ 24 ]. As chemotherapy options like docetaxel become more widely used, so does the importance of understanding how to better administer them so as to minimize the possibility that resistant phenotypes emerge. As a simple caricature, it can be described microdynamically as: when a cell tries to undergo mitosis, it will be killed with probability s . This means that at s = 0.5, there will be no growth in the affected population, since half the time that a cell tries to divide, it succeeds and becomes two cells and half the time it dies and becomes zero cells. For 0 < s < 0.5, the growth rate is slowed; and for 0.5 < s < 1 the growth rate is reversed. This is implemented in all fitness functions (except for chemo-resistant T R ) by multiplying fitness by (1 − 2 s ). Note that chemotherapy only impacts the bone via interfering with the cells that control bone perturbation, i.e. OB and OC.
Results
We will illustrate our results with simulations of a particular parameter settings given by death rate α = 0.005, selection strength σ = 0.05, bone adjustment rate 1 K = 2.5 , and tumour-environment ecological interactions of δ = 1.5, γ = 300, ϵ = 0.03.
With the model described we can recapitulate bone homeostasis when there are no PCa cells. Assuming initial densities of the different population as follows: ρ OB (0) = 0.001, ρ OC (0) = 0.01, ρ T ( 0 ) = ρ T R ( 0 ) = 0.0 , B (0) = 1. The plots shown in figure 1 show how the bone is initially resorbted, and then deposited as OCs and OBs work to restore the balance of bone to homeostatic levels after a simulated microfracture.
Tumor Introduction.
We now assume that a metastatic PCa cell has extravasated into an area of the bone that will undergo remodeling. In this case we assume that the initial conditions are: ρ OB = 0.001, ρ OC = 0.01, ρ T = 0.0005, ρ T R = 0.0001 , B (0) = 1. In this case the results can be seen in figure 2 , which shows how the growth of the PCa population leads to an increase of pathological bone and a dominance of the PCa cells over the OBs and the OCs.
Standard Chemotherapy.
Trying to asses the impact of chemotherapy we assumed the following initial conditions: ρ OB (0) = 0.001, ρ OC (0) = 0.01, ρ T (0) = 0.0005, ρ T R ( 0 ) = 0.0001 , B (0) = 1. Additionally, a simple chemotherapy regiment of s = 0.3 is initiated from timesteps 200 to 350. The results can be seen in figure 3 and show that, while bone continues to grow even under treatment, the tumor population takes a sharp decline highlighting the impact of chemotherapy.
Fractionated Treatment.
To understand the effect of varying treatment regimens, the model was used to simulate a wide number of combinations of treatment and treatment holidays. Each regiment is characterized by dividing the treatment into periods, each period representing the same amount of time. Each slot can either have value of (1) to signify that treatment is being applied or (0) when it is not. Thus, assuming a treatment with n different treatment periods we can consider 2 n different treatment regimens based on combinations of treatments and breaks.
The full treatment space also contains regimens that produce favorable results at the cost of eradicating the stromal population - which would be pressumably cytotoxic to the patient. Thus, we removed all treatments that reduced the stromal population below a threshold. This approach allows us to contrast continuous and a spectrum of fractionated regimens. From the space of 16 windows, where each window is 40 timesteps, Table 1 lists the top 5 treatment schedules that yielded the best results in terms of average bone size, as well as the treatment’s average tumor burden. Table 1 also lists two other treatment regimens: (1) a continuous treatment regiment and (2) a periodic treatment regiment for comparison.
Moreover, figure 4 shows all the potential treatments plotted with respect to average bone size and average tumor burden with three representative treatments shown to the side: a periodic treatment, a continuous treatment and a fractionated treatment.
Discussion
Intra-tumor heterogeneity is increasingly recognized as the key driver of evolution in cancer. This heterogeneity also explains the emergence of resistance to standard of care as well as targeted treatments. We presented an EGT model that captures this heterogeneity and the role of dynamic homeostasis of bone remodeling orchestrateding the interactions between bone-producing cells (osteoblasts) and bone-resorbting cells (osteoclasts). This homeostasis can be disrupted by tumor cells, whose growth, in turn can be prevented by the application of cytotoxic drugs such as docetaxel. Our results recapitulate not only the “vicious cycle” of the disruption and co-option of homeostatic mechanisms but also the pathological bone formation that characterizes bone metastatic prostate cancer [ 25 ]. Knowing that certain tumor phenotypes can be relatively immune to the effects of chemotherapy, our model also captures the emergence of therapy resistant tumors.
The conventional clinical application of chemotherapy in most cancers involves the use of the drug until either the tumor enters remission, resistance renders further application of the drug ineffective or the patient can no longer tolerate it. This approach is called Maximum Tolerable Dosage (MTD) and has been recently contrasted to alternative approaches where the aim is to transform the disease into one that could be managed as a chronic condition. In sipport of this, intermittent docetaxel has been proven to be useful in metastatic prostate cancer [ 26 , 27 ].
Our results show that the efficacy of the treatment depends on the heterogeneity of the tumor. MTD works best if the metastases are homogeneous. Assuming heterogeneity with regards to chemotherapy resistance, a realistic scenario in metastatic prostate cancer, our model allows us to explore how different treatment duration and intervals between drug application, impact tumor heterogeneity, fitness, and bone mass (see figure 4 ). Unsurprisingly the model shows that alternative treatment strategies are likely to yield better results if the metastasis contains both susceptible and resistant phenotypes.
Novel therapies where conventional drugs are used while taking into account the tumor’s evolutionary dynamics have been proposed and mathematically explored by Orlando et al. [ 28 ] under the constraint that the microenvironment of the tumor plays a reduced role and that tumor populations does not interact with each other. Our results extend beyond these constraints and show that neither conventional nor fractionated strategies are always the best solution. While our systematic search has not shown optimal treatments resembling the evolutionary-enlightened therapies proposed by Zhang et al. [ 22 ], there is some evidence that they could work in the context of metastatic prostate cancer. Our main result supports a variant: the application of fractionated therapies where full dosage is first applied for a period of time then followed by the alternation between on and off cycles.
While this finding will not have an immediate impact on how fractionated treatments are delivered, it shows that the model captures the basic biological principles of bone metastatic prostate cancer which allowed us to explore treatment strategies in a meaningful way. Models like the one we have presented here can be used to better understand evolutionary-enlightened treatments. EGT is a well known mathematical tool in which to frame evolutionary questions, and our EGT model allows us to include the microenvironmental selection as well as cell-cell interactions in the game dynamics. We are aware that, by its own nature, mathematical models constitute a simplification of the reality being modeled. For example, our model does not take into consideration the spatial interplay which is known to effect EGT dynamics of cancer [ 14 , 29 ] and has been recently shown to play a key role in the efficacy of some innovative applications of conventional treatments in cancer [ 30 ]. We expect that our results could be quantitatively different as we assume different costs of resistance, treatment impact, or spatial structure. But the principle that a period of continuous application followed by periods where on and off cycles is expected to hold when the parameters change to reflect different tumors and treatment efficacies.
We measured the impact of these different treatment schemes by looking at their effect on the two tumor populations but, in bone metastatic prostate cancer, there are other metrics that need to be considered: total dosage and microenvironmental impact. As shown in figure 2 , tumor growth leads to an increase in bone, growth that can, in some cases, be curtailed by chemotherapy (see figure 3 ). Figure 4 shows that the scheduling of chemotherapy can have a substantial impact on the amount of pathological bone formed, a consideration that is clinically important yet rarely included in mathematical models of bone metastases. The examples offer the interesting possibility that a sustained application of chemotherapy followed by treatment holidays and subsequent application of fractionated therapy could lead to control not only on tumor growth but also bone growth ( figure 4F ; although, see Kaznatcheev [ 31 ] for a perspective on the limits of tumour burden control in models like ours). While the total dosage is slightly different, the purely continuous and the (approximately) fractionated treatments shown in figures 4B and 4A demonstrate bone increases over time. Our EGT model, by combining tumor, stroma and bone microenvironment, constitutes a platform in which to investigate various treatment strategies – include, but are not limited, to bisphosphonates, anti-RANKL and hormonal therapies – to the standard of care for metastatic prostate cancer patients that can minimize dosage and bone disease.
Acknowledgements
We would like to acknowledge Dr. Heiko Enderling from Moffitts Integrated Mathematical Oncology department for the organization of the HIP-IMO program that led to this work. AA, CCL, and DB were partly funded by an NCI U01 (NCI) U01CA202958-01 and a Moffitt Team Science Award. AA was partly funded by a Department of Defense Prostate Cancer Research Program (W81XWH-15-1-0184) fellowship.
References
[1] National Cancer Institute . Cancer stat facts: Prostate cancer , 2017 At: https://seer.cancer.gov/statfacts/html/prost.html .
[2] Araujo Arturo , Cook Leah M , Lynch Conor C , and Basanta David . An integrated computational model of the bone microenvironment in bone-metastatic prostate cancer . Cancer Res , 74 ( 9 ): 2391 – 401 , 5 2014 . 24788098 10.1158/0008-5472.CAN-13-2652 PMC4023121
[3] Cook Leah M , Araujo Arturo , Pow-Sang Julio M , Budzevich Mikalai M , Basanta David , and Lynch Conor C . Predictive computational modeling to define effective treatment strategies for bone metastatic prostate cancer . Scientific Reports , 6 , 2016 . 10.1038/srep29384 PMC4944130 27411810
[4] Heppner Gloria H . Tumor Heterogeneity . Cancer Research , 1984 . 6372991
[5] Gallaher Jill , Cook Leah M , Gupta Shilpa , Araujo Arturo , Dhillon Jasreman , Park Jong Y , Scott Jacob G , Pow-Sang Julio , Basanta David , and Lynch Conor C . Improving treatment strategies for patients with metastatic castrate resistant prostate cancer through personalized computational modeling . Clin Exp Metastasis , 31 ( 8 ): 991 – 9 , 12 2014 . 25173680 10.1007/s10585-014-9674-1 PMC5399888
[6] Kaznatcheev Artem . Two conceptions of evolutionary games: reductive vs effective . bioRxiv , page 231993 , 2017 .
[7] Kaznatcheev Artem . Effective games and the confusion over spatial structure . Proceedings of the National Academy of Sciences , page 201719031 , 2018 . 10.1073/pnas.1719031115 PMC5828619 29444859
[8] Ryser Marc D , Nigam Nilima , and Komarova Svetlana V . Mathematical Modeling of Spatio-Temporal Dynamics of a Single Bone Multicellular Unit . Journal of Bone and Mineral Research , 2009 . 10.1359/jbmr.081229 19063683
[9] Ryser Marc D , Qu Yiding , and Komarova Svetlana V . Osteoprotegerin in Bone Metastases: Mathematical Solution to the Puzzle . PLoS Comput Biol , 2012 . 10.1371/journal.pcbi.1002703 PMC3475686 23093918
[10] Buenzli PR , Pivonka P , and Gardiner BS . Modelling the anabolic response of bone using a cell population model . Journal of theoretical …, 2012 . 10.1016/j.jtbi.2012.04.019 22579551
[11] Tomlinson IPM and Bodmer WF . Modelling the consequences of interactions between tumour cells . Brit. Jour. Cancer , 75 ( 2 ): 157 – 60 , 1997 . 9010019 10.1038/bjc.1997.26 PMC2063276
[12] Tomlinson IP . Game-theory models of interactions between tumour cells . European journal of cancer (Oxford, England : 1990) , 33 ( 9 ): 1495 – 1500 , 8 1997 . 10.1016/s0959-8049(97)00170-6 9337695
[13] Basanta D , Hatzikirou H , and Deutsch A . Studying the emergence of invasiveness in tumours using game theory . The European Physical Journal B , 63 ( 3 ): 393 – 397 , 6 2008 .
[14] Kaznatcheev Artem , Scott Jacob G , and Basanta David . Edge effects in game-theoretic dynamics of spatially structured tumours . Journal of The Royal Society Interface , 12 ( 108 ): 20150154 , 2015 . 10.1098/rsif.2015.0154 PMC4528581 26040596
[15] Basanta D , Simon M , Hatzikirou H , and Deutsch A . Evolutionary game theory elucidates the role of glycolysis in glioma progression and invasion . Cell proliferation , 41 ( 6 ): 980 – 987 , 12 2008 . 19040573 10.1111/j.1365-2184.2008.00563.x PMC6495695
[16] Basanta D , Scott JG , Fishman MN , Ayala G , Hayward SW , and Anderson ARA . Investigating prostate cancer tumour-stroma interactions: clinical and biological insights from an evolutionary game . British Journal of Cancer , 106 ( 1 ): 174 – 181 , 1 2012 . 22134510 10.1038/bjc.2011.517 PMC3251863
[17] Kaznatcheev Artem , Velde Robert Vander , Scott Jacob G. , and Basanta David . Cancer treatment scheduling and dynamic heterogeneity in social dilemmas of tumour acidity and vasculature . British Journal of Cancer , 2017 . 10.1038/bjc.2017.5 PMC5355932 28183139
[18] Dingli D , Chalub FACC , Santos FC , van Segbroeck S , and Pacheco JM . Cancer phenotype as the outcome of an evolutionary game between normal and malignant cells . British Journal of Cancer , 101 ( 7 ): 1130 – 1136 , 10 2009 . 19724279 10.1038/sj.bjc.6605288 PMC2768082
[19] West Jeffrey , Ma Yongqian , and Newton Paul K . Capitalizing on competition: An evolutionary model of competitive release in metastatic castration resistant prostate cancer treatment . bioRxiv , 2018 . 10.1016/j.jtbi.2018.07.028 PMC7519622 30048718
[20] Bilezikian John P , Raisz Lawrence G , and Martin T. John . Principles of bone biology . Academic Press/Elsevier , San Diego, Calif. , 3rd ed edition , 2008 .
[21] Cook Leah M , Shay Gemma , Araujo Arturo , Aruajo Arturo , and Lynch Conor C . Integrating new discoveries into the “vicious cycle” paradigm of prostate to bone metastases . Cancer Metastasis Rev , 33 ( 2–3 ): 511 – 25 , 9 2014 . 24414228 10.1007/s10555-014-9494-4 PMC4096318
[22] Zhang Jingsong , Cunningham Jessica J , Brown Joel S , and Gatenby Robert A . Integrating evolutionary dynamics into treatment of metastatic castrate-resistant prostate cancer . Nature communications , 8 ( 1 ): 1816 , 2017 . 10.1038/s41467-017-01968-5 PMC5703947 29180633
[23] Kaznatcheev Artem . Evolutionary dynamics of cancer in the bone . Theory, Evolution, and Games Group , 12 ( 6 ), 2015 At: https://egtheory.wordpress.com/2015/12/06/bone-cancer/ .
[24] Henríquez I , Rodríguez-Antolín A , Cassinello J , San Segundo C Gonzalez , Unda M , Gallardo E , López-Torrecilla J , Juarez A , and Arranz J . Consensus statement on definition, diagnosis, and management of high-risk prostate cancer patients on behalf of the spanish groups of uro-oncology societies uroncor, guo, and sogug . Clinical and Translational Oncology , pages 1 – 10 , 2017 . 28785912 10.1007/s12094-017-1726-y
[25] Keller Evan T and Brown Julie . Prostate cancer bone metastases promote both osteolytic and osteoblastic activity . J Cell Biochem , 91 ( 4 ): 718 – 29 , 3 2004 . 14991763 10.1002/jcb.10662
[26] Beer TM , Garzotto Mark , Henner WD , Eilers KM , and Wersinger EM . Multiple cycles of intermittent chemotherapy in metastatic androgen-independent prostate cancer . British journal of cancer , 91 ( 8 ): 1425 , 2004 . 15467765 10.1038/sj.bjc.6602198 PMC2409932
[27] Mountzios Ioannis , Bournakis Evangelos , Efstathiou Eleni , Varkaris Andreas , Wen Sijin , Chrisofos Michael , Deliveliotis Charalambos , Alamanis Christos , Anastasiou Ioannis , Constantine Constantinides , Intermittent docetaxel chemotherapy in patients with castrate-resistant prostate cancer . Urology , 77 ( 3 ): 682 – 687 , 2011 . 21256546 10.1016/j.urology.2010.08.044
[28] Orlando Paul A , Gatenby Robert A , and Brown Joel S . Cancer treatment as a game: integrating evolutionary game theory into the optimal control of chemotherapy . Phys Biol , 9 ( 6 ): 065007 , 12 2012 . 23197192 10.1088/1478-3975/9/6/065007 PMC3653600
[29] Nanda Mridu and Durrett Richard . Spatial evolutionary games with weak selection . Proceedings of the National Academy of Sciences , page 201620852 , 2017 . 10.1073/pnas.1620852114 PMC5468649 28533405
[30] Bacevic Katarina , Noble Robert , Soffar Ahmed , Orchid Wael Ammar Benjamin Boszonyik , Prieto Susana , Vincent Charles , Michael E Hochberg Liliana Krasinska , and Fisher Daniel . Spatial competition constrains resistance to targeted cancer therapy . Nature communications , 8 ( 1 ): 1995 , 2017 . 10.1038/s41467-017-01516-1 PMC5722825 29222471
[31] Kaznatcheev Artem . Drug holidays and losing resistance with replicator dynamics . Theory, Evolution, and Games Group , 9 ( 2 ), 2016 At: https://egtheory.wordpress.com/2016/09/02/drug-holidays/ .