Review: Mathematical Modeling of Prostate Cancer and Clinical Application

We review and synthesize key findings and limitations of mathematical models for prostate cancer, both from theoretical work and data-validated approaches, especially concerning clinical applications. Our focus is on models of prostate cancer dynamics under treatment, particularly with a view toward optimizing hormone-based treatment schedules and estimating the onset of treatment resistance under various assumptions. Population models suggest that intermittent or adaptive therapy is more beneficial to delay cancer relapse as compared to the standard continuous therapy if treatment resistance comes at a competitive cost for cancer cells. Another consensus among existing work is that the standard biomarker for cancer growth, prostate-specific antigen, may not always correlate well with cancer progression. Instead, its doubling rate appears to be a better indicator of tumor growth. Much of the existing work utilizes simple ordinary differential equations due to difficulty in collecting spatial data and due to the early success of using prostate-specific antigen in mathematical modeling. However, a shift toward more complex and realistic models is taking place, which leaves many of the theoretical and mathematical questions unexplored. Furthermore, as adaptive therapy displays better potential than existing treatment protocols, an increasing number of studies incorporate this treatment into modeling efforts. Although existing modeling work has explored and yielded useful insights on the treatment of prostate cancer, the road to clinical application is still elusive. Among the pertinent issues needed to be addressed to bridge the gap from modeling work to clinical application are (1) real-time data validation and model identification, (2) sensitivity analysis and uncertainty quantification for model prediction, and (3) optimal treatment/schedule while considering drug properties, interactions, and toxicity. To address these issues, we suggest in-depth studies on various aspects of the parameters in dynamical models such as the evolution of parameters over time. We hope this review will assist future attempts at studying prostate cancer.


Introduction
The past 20 years have seen an accelerating development of mathematical models for prostate cancer. One reason is the success of earlier work to explain and replicate experimental and clinical data. Another reason is the availability of experimental and clinical data due to collaborations with physicians and advancements in technology. Thus far, modeling efforts have utilized a broad range of model types and methods to uncover some fundamentals about prostate cancer progression and the efficacy of its treatments. However, there has not been a systematic and comprehensive synthesis It is difficult to state exactly what contributes to the development of prostate cancer; however, age, race, and hereditary factors are the most firmly supported risk factors for prostate cancer [5]. While prostate cancer is most commonly diagnosed in men over 50, a benign tumor can start its development decades earlier. The slow progression rate of prostate cancer means that often, patients diagnosed with prostate cancer will not die because of it. Furthermore, the side effects of prostate cancer treatment can be significant, so active surveillance without treatment is the recommended management for some cases of low-risk disease [6].

Prostate Cancer: Physiology and Treatments
In this section, we aim to give a summary of the biological and clinical details of prostate cancer and its treatments.
The prostate is a walnut-sized gland that is part of the endocrine system and the male reproductive system. Its main job is to secrete seminal fluids that are used to nourish and transport sperm. The function, proliferation, and atrophy of healthy and cancerous prostate cells depend highly on the circulating level of androgen, particularly testosterone [7]. Androgens are steroid hormones that regulate the development and maintenance of male characteristics by binding to androgen receptors (AR), which is a cytoplasmic receptor that is activated by binding to androgens, which induce the AR to migrate to the nucleus and bind to the Androgen Response Elements (ARE) within the DNA [8]. The testicular-hypothalamic-pituitary axis, also known as the hypothalamus-pituitary-gonadal axis-where gonads indicate the reproductive system-produces androgen in a negative feedback loop. This means that the higher the androgen level in the body, the lower its production rate becomes [9]. Approximately 90%-95% of the body's androgen is produced by the testes and the remainder by the adrenal gland [10]. Once androgens, or testosterones, enter the prostate, about 90% is converted to the more potent form 5α-Dihydrotestosterone (DHT) by the 5α-reductase enzyme [11]. DHT has been observed to be a better agonist to the androgen receptors, with about 2-3 and 15-30 times higher affinity than testosterone and adrenal androgens, respectively [12,13]. Additionally, the unbinding rate of DHT and AR is roughly five times slower than that of testosterone and AR [14]. Once androgens bind to the AR, they are transported to the nucleus to bind to the ARE. This process results in the transcription of androgen-regulated genes, which signals the proliferation and inhibition of apoptosis in both healthy and cancerous prostate cells. The proliferation of prostate cells produces a prostate-specific antigen (PSA). Most of this byproduct of prostate cell growth is normally contained within the prostate wall; however, in the presence of prostate cancer cells, the cell wall is disrupted, which allows for an increased leakage of PSA into the bloodstream. For this reason, a simple blood test can be used to measure the level of PSA, which acts as a biomarker for prostate cancer and is often used as a proxy for tumor volume in mathematical models [15]. A simplified schematic of prostate cancer dynamics is described in Figure 2. A simplified schematic of the dynamics of androgen and prostate-specific antigen in prostate cancer-adapted from [16]. (1) Serum androgen flows into the tissue of prostatic epithelial cells. (2) It is then converted to Dihydrotestosterone (DHT) by 5-al pha reductase. (3) Both testosterones and DHT bind to androgen receptors. (4) This further activates the proliferative and survival signals. (5) Cancer activity results in the production of prostate-specific antigen (PSA), which is leaked out of the cells (6). Note that the gray compartments indicate measurable concentrations that are relatively easier to obtain.
The diagnosis of prostate cancer usually starts with a blood test to measure the PSA level. To verify the presence of a prostate tumor, a digital rectal exam, a transrectal ultrasound, and/or a prostate biopsy may be carried out. If cancer is confirmed, clinicians also assign the patient a Gleason score, which measures the resemblance of biopsy cancer tissues to healthy tissues. The Gleason score is often used as part of the clinical process to determine which treatment would be best for the patient [17], but it is not often utilized in mathematical models. For localized prostate cancer, the main options are: watchful waiting and active surveillance, surgical removal of the prostate, and radiation. If the cancer has metastasized, hormonal therapy is the standard treatment. Initially, hormonal therapy was administered as a continuous treatment, also known as continuous androgen suppression therapy (CAS). However, hormonal therapy cuts off the production of male hormones, so it often associates with a high level of toxicity. Some of the common side effects are impotence, depression, bone demineralization, and dementia [1]. Furthermore, cancer cells eventually develop resistance to hormonal therapy. Therefore, instead of CAS, intermittent androgen suppression therapy (IAS) was introduced as an alternative to reduce the side effects of CAS and potentially delay the time to treatment failure [18]. In contrast with CAS where the patient is put on treatment until it becomes ineffective, during IAS, patients go on and off of therapy according to either a fixed interval of time or threshold PSA level. The on-off cycles are repeated until the treatment is ineffective. While IAS gives patients a better quality of life, it is controversial whether it is better than CAS for delaying the onset of treatment resistance [19]. Recently, the concept of adaptive therapy was introduced into clinical studies of prostate cancer, which has generated some exciting results and opportunities [20]. On the other hand, while hormonal therapy is the focus of many modeling efforts, immunotherapy and chemical therapies can also be used to treat prostate cancer-usually after the cancer has developed treatment resistance to hormonal therapy [21,22].

Mathematical Models for Prostate Cancer
Due to a lack of detailed knowledge of the complex mechanisms of PCa and the effects of numerous therapies, some researchers have begun to use mathematical models to better explain the observations from experimental and clinical studies [15]. Whereas statistical models can help to establish the relationships between variables in a study, mathematical models pursue the dynamics of the interaction of those variables. One could argue that the studies using mathematical models should come after the observed relationships of the variables. Additionally, mathematical models can be developed using known physical, biological, or empirical laws and incorporate control effects (such as drugs), making them suitable for a variety of purposes, such as testing a hypothesis or finding an optimal combination of drugs. This means that mathematical models often contain assumptions on the nature of the phenomenon being studied. Furthermore, mathematical models can range wildly in their complexity, which may result in either highly unrealistic or mathematically and computationally untractable models. There have been numerous reviews of the growing literature on mathematical oncology in the past 15 years [15,[23][24][25][26]. However, none of these reviews exclusively and extensively cover mathematical models for prostate cancer. This review aims to address the current lack of an extensive review of the rapidly-developing area of mathematical modeling for prostate cancer. To do so, we exhaustively collect existing mathematical models developed specifically for prostate cancer and synthesize their results. First, we present conclusions that are independent of the details of the models. This is carried out by categorizing the models based on some important characteristics and then finding the common conclusion of models that belong in the same category. Secondly, we provide a summary of key findings and their implications for each model. Finally, we point out relevant but lesser explored concepts related to prostate cancer and its treatment that may be addressed by future research. Before carrying out our analysis, we first introduce the terminology and key characteristics of existing models.

Elements of Mathematical Models for Clinical Applications
In this section, we introduce and discuss the elements of mathematical models that are important for transitioning from theoretical studies to clinical applications. While we do not consider all aspects of mathematical models, the chosen categories should serve as a solid foundation for a fruitful discussion.

Population Structure and Dynamics
To have a complete understanding of the structure of prostate cancer and how the subpopulations of cells interact, with or without treatment, would be similar to having solved the problem of cancer from a theoretical perspective. A system-biology approach treats cancer as a unified population of cells that is part of a larger ecosystem. On the other hand, cancer can also be looked at like a society with many interacting sub-types of cancer cells. While it would perhaps be ideal to understand the details of each aspect of prostate cancer, it is impossible to do so. Instead, each model has to make certain simplifying assumptions that aim to answer a specific set of questions. Thus, we will examine the underlying assumptions of the population structure and dynamics on which each model was built and analyzed. Specifically, we aim to explore the following features of mathematical models in more detail.
Tumor heterogeneity and evolution: treatment failure due to cancer cells developing resistance is an unsolved problem in treating advanced prostate cancer. The most likely explanation is that over the course of treatment, selective pressure allows selection for resistant cell types. Thus, all existing population-type prostate cancer models incorporate some mechanisms that describe this heterogeneity and the evolution of the cancer population during treatment. For this feature, we study the existing literature to obtain conclusions that explore the heterogeneity and evolution of a tumor and its effects.
Tracking the progression of tumor size: having accurate and consistent measurements for the growth of the tumor that can be incorporated into mathematical models is crucial for clinical applicability and model assessment. While imaging data is perhaps best suited when it comes to measuring the changes in the tumor, the challenges of collecting high-quality spatial data hinder this possibility. Instead, time-serial PSA data has often been used as a proxy for the progression of prostate cancer. In this regard, we synthesize suggestions from existing work to provide an assessment of PSA as a biomarker for prostate cancer and explore other alternatives.
Model types and dynamics: some of the most successful prostate cancer models make use of the characteristics of the cancer progression or the observed data, while others construct more mechanistic models that make use of either fundamental or physical laws, like conservation laws, and/or the biochemical networks of cancers. In essence, mathematical modeling of prostate cancer is at the interface between simple classical mathematical models and complex system-biology type models, which offers rich opportunities to study interesting properties of these models from a mathematical, biological, or computational perspective. Additionally, while modeling with ordinary differential equations still dominates existing studies, the potential for using other types of dynamics has also been explored. For this reason, we compare and contrast different types of models and how they may offer insights for treating prostate cancer.

Applicability in Clinical Settings
While theoretical work offers insights into the fundamentals of prostate cancer, the complex interactions of various components in a real situation blur the implications from theories. Moreover, in the case of mathematical models, the predictions of the progression of the tumor must be validated against real-world clinical data, which often means fitting and forecasting using incomplete sets of data. Aside from these theoretical issues, clinical applications come with another set of questions that are often case-specific and not explored in depth in theoretical settings. For example, how does one predict whether a certain treatment will be effective for a specific patient? If so, how long will it take for the treatment to be effective? If not completely, then for how long would such treatment be effective? What is the dose level of certain drugs that should be used to achieve optimal results? What is the best combination of drugs to optimize effectiveness? Cancer treatments can be invasive and expensive, which can harm patients physically and financially [1,6]. Although treating cancer is a nonlinear process, this is not the case in practice; different adjustments to the treatment should be made based on the specific situation. Existing modeling studies attempt to give some hints to possible treatments. However, other features of treatment, such as dosage, combinations of drugs and therapies, can also be optimized to achieve better results. For these reasons, we discuss several important gaps between theoretical work and the clinical application of mathematical models.
Real-time estimability: in theoretical studies, researchers often have the luxury of knowing the complete set of a patient's data. However, in the clinical setting, researchers must make do with a smaller but increasing set of a patient's data. This means the parameters for the model cannot be identified with high certainty at the earlier stage, rendering implications from mathematical models unreliable.
Uncertainty, identifiability, and sensitivity: another concern that arises with models of high complexity is the sensitivity of the model outputs to variations in the parameters. If the variation is small, it is usually referred to as sensitivity analysis, whereas uncertainty quantification refers to larger variation in the parameters, perhaps due to the large uncertainty in the data. Since measurements and models are imperfect, uncertainty quantification allows physicians to put confidence bounds on a model's predictions. On the other hand, sensitivity analysis can help pinpoint the important parameters in a model, which can aid with decision making during treatment.
In close connection with sensitivity and uncertainty, the term identifiability refers to the ability to determine a unique set of model parameters with respect to the given data. This concept is further broken down into structural identifiability and practical identifiability. Structural identifiability refers to the relationship between the inherent structure of the model and perfect data. In essence, structural identifiability asks whether all of the parameters in a model can be uniquely identified given some perfect set of data on the variables in the model. On the contrary, when the data contains errors, the same question refers to the practical identifiability of the model. While these are important topics, there have been few attempts to study these closely related aspects in the case of prostate cancer.
Optimal schedule, optimal treatment, and patient classification: while the previous issues deal with the validation of mathematical models in clinical settings, optimal schedule, optimal treatment, and patient classification are important contributions of mathematical models. For instance, patients can be categorized in advance based on whether they would react positively toward a certain type of treatment using mathematical models. Furthermore, to fully draw out the effect of a treatment, a mathematical model can be used to study how to best administer the treatment as pertaining to treatment schedule or in combination with other treatments. Additionally, other factors, such as cost, could also be taken into account for an optimal treatment study.

Population-Based Models of Tumor Relapse
In this section, we discuss the implications from mathematical models. A synthesis of our findings is presented in Table 1.

Tracking the
The most commonly used biomarker for prostate cancer growth (PSA) is useful but can be unreliable. Instead, other measurements can be used in place or in concurrence with PSA to track tumor progression.
The accuracy of tracking the progression of tumors using multiple biomarkers needs to be examined. While accuracy is key, the availability of such biomarkers should also be taken into account.
progression of tumor Section 3.5 Model types Ordinary differential equations build the foundation for studying prostate cancer. However, the lack of various modes of modeling implies many aspects of prostate cancer are left unexplored.
As more data becomes available, especially imaging data, spatial, memory-based, and stochastic models will become useful in capturing spatial patterns in cancer progression and interaction, specifically, the metastatic processes. and dynamics Section 3.6 The highly heterogeneous nature of prostate cancer, especially after it has metastasized, is one important reason why population-type models take up a large percentage of existing work. This includes simple models that describe only the interactions among sub-populations of cells, or more complex models that also account for external factors, such as environmental changes due to treatment. The population-type models, in general, focus on the interactions of the cancer sub-populations with or without external effects. On the other hand, kinetic-type models tend to view cancer as a whole and how cancer interacts within the human physiological network, often via biochemical pathways.
Yorke et al. (1993). The earliest approach to model prostate cancer progression mathematically is credited to Yorke et al. [27]. In their work, they developed a simple kinetic model based on Gompertzian growth that aims to explain the progression of prostate cancer, its metastasis process, and how the metastasis process affects local treatments (surgical and radiation therapy). By comparing their model simulations with the patient's survival data in Fuk et al. [28], Yorke and colleagues suggested that local relapse is associated with more aggressive cancer phenotypes. Furthermore, due to the aggressive cancer phenotypes observed in local relapse, they suggest that local control (e.g., the complete eradication of the primary tumor and lymph node metastasis) can significantly improve the disease-free survival and metastasis-free survival time. However, after the cancer has metastasized, standard treatments show early positive results but fail to contain tumor relapse.
Jackson (2004). The relapse following positive responses to hormonal castration has long been hypothesized to be a result of prostate cancer being made up of different subpopulations of cells, each with its own level of dependency on the extracellular androgen. This concept is first demonstrated with mathematical modeling in 2001 by Jackson [29]. Jackson balances proliferation, apoptosis, and fluxes of two types of cells under radial symmetry to construct a system of partial differential equations for androgen deprivation therapy in prostate cancer, as shown in Equations (1)- (2).
In this model, u(R, t) is the net rate of collective cellular motion, R(t) is the tumor radius and p(R, t), q(R, t) are the volume fractions of androgen dependent (AD) and androgen independent (AI) prostate cancer cells, respectively. The model assumes that a high level of androgen increases the proliferation rate of androgen dependent cells but does not affect the proliferation rate of androgen independent cells. On the other hand, a low level of androgen decreases the apoptosis rate of androgen dependent cells and increases the apoptosis rate of androgen independent cells. With these assumptions, the models were able to capture the qualitative dynamics observed from experimental data of LuCaP 23 in mice [30]. Thus, Jackson concludes that the model supports the notion that cancer relapse occurs because AI cells have a lower apoptosis rate, rather than a higher proliferation rate, a conclusion that agrees with some experimental data [31]. This conclusion supports the hypothesis that prostate cancer cells pay a cost to evolve resistance to treatment-thus having a disadvantage against AD cells when no treatment is applied [20].
Additionally, Jackson shows that androgen deprivation therapy is only effective for a small region of parameters-which agrees with clinical data showing that androgen deprivation therapy is not curative, even in early stage disease (see Figure 3). As with all models, Jackson's modeling effort has certain limitations. First, the assumption of spherical symmetry in the model formulation limits the application of the model to the early stage, where the cancer can be approximated with a sphere. Secondly, the assumption that increasing the level of androgen increases the apoptosis rate of androgen independent cells is arguable because androgen independent cells may still be dependent on AR, which is known as the outlaw pathway to resistance [11,15,32]. Regardless, Jackson's work has paved the way for future modeling studies.

Ideta et al. (2008).
Inspired by Jackson's work, Ideta et al. [33] studied a system of ordinary differential equations describing the interaction of androgen dependent and independent cells, while explicitly treating the androgen level as a dynamical variable.
Here, a represents the androgen level, the index i runs between 1 and 2, which represents the androgen dependent cells and androgen independent cells, respectively. p i (a), q i (a), and m(a) are functions that represent the effects of androgen on the proliferation, death, and mutation rates of the respective cells. The authors consider three different assumptions of how androgen affects the growth and death rates of the androgen independent population (p 2 (a) is different for each assumption). Furthermore, the model incorporates the idea that androgen dependent cancer cells can transform into androgen independent cancer cells due to survival pressure.
Ideta et al. used this model to compare the effects of intermittent androgen deprivation therapy verses continuous androgen deprivation therapy, making it the first mathematical modeling work to carry out such a comparison in the decade-old debate of CAS and IAS. They observed that relapse can be averted only under the assumption that the AI cell population decreases when the androgen level is normal. Their comparison concludes that if androgen dependent cancer cells have a competitive advantage over androgen independent cancer cells without treatment, then IAS is superior in prolonging the onset of treatment resistance than CAS. It is noteworthy to point out that both Jackson and Ideta make the assumption that AI cells do not have a competitive advantage over AD cells. This assumption is partially justified due to the rare occurrence of AI cells in early PCa, yet AI cells are quite common after the cancer has metastasized.
While the mechanisms leading to the occurrence of AI relapse are still somewhat of a mystery, the study by Zhang and colleagues [20] shows strong evidence to support the assumption that AI cells are out-competed by AD cells in an androgen rich environment. In addition, Zhang et al. hypothesizes that in order for cancer cells to obtain treatment resistance, they must expend part of their system to develop some mechanisms to do so, which makes them less fit than the cells that do not have such mechanisms. However, Jackson's and Ideta et al. models do not directly incorporate the competitive effect of the two populations.

Shimada and Aihara (2008), Yang et al. (2016).
To address this limitation, Shimada and Aihara [34] and Yang et al. [35] proposed models that directly account for competition between the two cancer phenotypes. The competition model of Shimada and Aihara introduces a competition term with a constant rate, whereas the competition terms in Yang et al. depend on the androgen level and are asymmetrical. Furthermore, Yang et al. also introduce intraspecific competition, which addresses the theoretical issue of having an environment that can support an infinite amount of cancer cells in Ideta's work. Numerical simulations for both models show the potential successes of IAS, which increases with higher competitive advantages for AD cells.

Guo et al. (2008), Tao et al. (2010).
Aside from approaches using ordinary differential equations, Jackson's model inspired spatial models by Guo et al. [36] and Tao et al. [37]. The work by Guo et al. is a spatial extension of the Ideta et al. model using some assumptions from Jackson's formulation. Their results give similar conclusions to Ideta's. Tao et al., on the other hand, analyzed the case where a mutation inhibitor is used as part of the regime for androgen deprivation therapy. They conclude that CAS can neither control nor cure PCa, even with a mutation inhibitor. However, their simulations suggest that mutation inhibition may delay the relapse of cancer.

Friedman and Jain (2013), Lorenzo et al. (2016)
In addition to these studies, Friedman and Jain [38] construct a spatial model of prostate cancer based on a similar framework and prove the existence and uniqueness of solutions for the model. Taking on a more computational approach, the phase-field method, which accounts for the dynamics and co-existence of healthy and cancerous cells, has also been applied in Lorenzo et al. [39], where the authors simulated several observed patterns of growth in prostate tumor (in vitro) with their model.

Hirata et al. (2010).
These early approaches suffer from a common problem: they are observation driven but not data validated, meaning the model is not tested against real-world data quantitatively. This brings us to the work by Hirata, Bruchovsky, and Aihara (HBA) [40]. In 2010, Hirata et al. took a different approach to develop a model for prostate cancer. Hirata and colleagues observed that when treatment starts, the PSA level first plummets to a certain level, then decreases more slowly to an undetectable level. Hirata et al. proposed that the decay trend of PSA is most naturally modeled by a three-population model with interacting populations via linear rates, which results in a matrix representation. To account for the effect of treatment, Hirata et al. use a discrete switching mechanism. When the patient is on treatment, the HBA model takes the form: and during the off-treatment period, The parameters w ij are the transformation rates from x i to x j and for simplicity, the PSA level is taken as the sum of x 1 , x 2 , x 3 . The model is tested against different sets of clinical data [40,41]. While the model is phenomenologically based, it has successfully reproduced the highly nonlinear trends in PSA data in many instances [23,24]; see Figure 4. Hirata and colleagues also use the model to explore the classification of cancer patients to determine which treatment (IAS or CAS) may benefit a patient the most. This is done by looking at the numerically estimated parameters for each patient based on their PSA data. Furthermore, the optimal treatment schedule has also been studied using their model [42].

Portz et al. (2012).
While purely phenomenological models, such as the HBA model, are useful in establishing the core structure of cancer dynamics, they are perhaps better viewed as a building block to a more comprehensive and biologically tractable model. In this regard, and Portz, Kuang, and Nagy (PKN) [43] formulated a mechanistic model of prostate cancer in 2012. The model utilizes the idea of cell quota as previously studied by Droop in a marine ecology context [44] and has gained traction since. The PKN model assumes only two subpopulations of cancer cells that can mutually transform from one to the other. The model assumes that androgen is the limiting nutrient (or driving force) in the proliferation and death of cancer cells and incorporates it as a cell quota. Specifically, the model presumes that there is a minimum quantity of androgen required for the growth of both AD and AI cells.
In this model, x 1 and x 2 are the AD and AI cancer cells, respectively. Q i is the intracellular androgen level, in contrast to the serum androgen level A. PSA is represented by the variable P. While the most commonly presented form of the PKN model is in Equation (7), the cell-quota assumption in the model means that the growth term takes a more accurate form of µ m max{0, 1 − q i Q i }x i . In this form, the biological meaning of the model is preserved when Q i dips below q i . This should be noted in both mathematical analysis and data fitting to avoid negative population simulations [45]. The model is first validated using the data from Akakura [46]. The PKN model requires the input of serum androgen data, so an exponential model is used to interpolate serum androgen levels. Similar to the HBA model, the PKN model also reproduces the time series observed in clinical data (see Figure 5). Additionally, a comparison of the HBA and the PKN models shows that they are equally capable of fitting and forecasting clinical PSA data [23,24]. shows the level of intracellular androgen for AD and AI cells (Q 1 and Q 2 , respectively) relative to the minimum cell quotas for AD and AI cells (q 1 and q 2 , respectively).

Baez and Kuang (2016).
The work by Portz et al. has inspired a series of modifications to address various limitations of the model. First, the PKN model requires interpolation of androgen data, so an assumption of the future androgen level is required, which makes accurate forecasts difficult to obtain. Thus, Baez and Kuang (BK) [47] introduced an improved version of the model where the main difference is the direct incorporation of androgen as a variable in the model, similar to the Ideta et al. model. This allows for the simultaneous fitting of androgen and PSA data. Furthermore, by having androgen as a variable, forecasts can be compared directly to clinical data.
Aside from addressing the limitation of forecasting in the PKN model, the BK model also implicitly addresses another limitation, which is pointed out by Hatano et al. [23]. Using approximation, Hatano et al. showed that the PKN model is not able to reproduce PSA relapse under continuous androgen deprivation therapy. This was missed by the original effort, most likely because their attention was focused on intermittent therapy. Phan et al. [48] show that the main reason for this biological limitation is due to the PKN model having a reversible transformation between the two subpopulations. On the other hand, by having only a transformation from AD to AI cells, which inversely depends on the serum androgen level, the BK model avoids this biological limitation and can produce relapse for CAS. This leads to the suggestion that the simplest way to ensure this biological aspect is by formulating the model with a unidirectional transformation from AD to AI cells. Furthermore, their conclusion shows that the success of cancer treatment is often heavily influenced by the structure of the model. Specifically, if a model assumes reversible transformation between the cancer sub-populations, then the control of the tumor tends to be simpler. Whereas a unidirectional transformation to a resistant population almost guarantees resistance will take place.

Phan et al. (2019).
Phan et al. [48] also compared a two subpopulation version of the BK model and a three subpopulation version with a similar structure to the HBA model. They observed that the limitation of data types and data points hinders significant improvements in using a more complex model. This concept is further analyzed using a Fisher information matrix and profile likelihood in Wu et al. [49], where a similar conclusion is reached. Additionally, in the work by Baez and Kuang, upon off treatment, the level of androgen grows very rapidly-causing the model to overshoot in its forecasts of androgen and PSA levels. Phan et at. [50] pointed out that this is perhaps due to two main reasons. First, the parameter associated with androgen production seems to be much larger than previous estimates. Second, the BK model assumes instantaneous diffusion of serum androgen to intracellular androgen. Thus, a new model is constructed by adding a compartment for serum androgen. By addressing these two features, Phan and colleagues show an improvement in the fitting and forecasting of androgen and PSA using clinical data. This suggests that separating intracellular and serum androgen is an effective and natural method to improve model fitting. As hormonal therapy is the standard of care for metastatic prostate cancer, it is crucial to obtain a good understanding of how the drugs used in hormonal therapy affects the androgen level and consequently the prostate cancer cells. To this end, the work by Barton and Andersen in 1998 [9] has paved the way for an initial framework of androgen regulation of prostate growth. Its extension and validation using rat data was later carried out in Potter et al. [51]. Recently, Reckell et al. [52] also formulated a pharmacokinetics model that incorporates the specific properties of a drug on androgen production in a prostate cancer model and tests it with clinical patient data. While Reckell et al. focus on the effects of combination hormonal drugs (Cyproterone acetate and Leuprolide acetate) on a macroscopic level using phenomenological functional responses, the stochastic differential equation phamacokinetics model by Cerasuolo et al. [53] attempts to capture the molecular interactions (oxygen consumption and protons production of cancer cells) between the cells and a newer hormonal drug (enzalutamide) using mice data. In essence, these works are successful in their attempts to replicate the drug effects on the dynamics of androgen and prostate growth; however, the complexity of these models often make them mathematically untractable and highly difficult to fit with high certainty. To address this issue, in both pharmacokinetics attempts, the parameter estimations are carried out in a multi-level process. a multi-scale model to study the effect of androgen in the evolutionary process from the benign to treatment stage of prostate cancer. The model couples the AR kinetics model with a state-transition model, where 100 states included with each state representing a different strain of prostate cells with varying AR expression. This reflects the observation that AR up-regulation is one of the most important pathways in which cancer cells becomes AI. Their results demonstrate how a heterogeneous population of prostate cancer cells can be skewed to select for androgen independent phenotypes in a low androgen environment. This suggests that while low levels of androgen may delay the appearance of malignant cancer cells, it may increase the chance of more aggressive cancer phenotypes. This result offers a reasonable explanation to the observation that finasteride, a 5-α reductase inhibitor, can reduce the overall rate of prostate cancer but may increase the rate of high-grade prostate cancer [55].  [16] applied the framework into a two cell types model under androgen deprivation therapy (both intermittent and continuous). Since their model is rather complex, they carry out multi-level fitting to obtain the estimates for the parameter values separately. By fitting average patient data reported in Goldenberg et al. [56], they show that their model not only contains patient specific parameters but is also capable of reproducing a variety of clinically observed dynamics of cancer progression. The model predicts similar conclusions as previous work: intermittent androgen suppression therapy is superior compared to continuous therapy when the AD cells have a competitive advantage over the AI cells and vice versa. Due to the complexity of the model, not much analytical work can be done. Furthermore, the large number of parameters presented in the model means clinical parameter estimation for each patient is difficult. Therefore, in subsequent work, Jain and Friedman [57] simplify their model and carry out mathematical analysis on the simplified version. By defining rigorously the definition of treatment viability and failure, they were able to compare the efficacy of continuous vs. intermittent therapy, which yields a similar conclusion to previous work. Additionally, they discover that even if it is possible to control a tumor with an optimal schedule of intermittent therapy, a sub-optimal one may still lead to the emergence of treatment resistance.

Zhang et al. (2017).
In the same spirit of using a simple phenomenological model, but in contrast to Hirata et al. where the focus is to examine the mutation process that governs the dynamics of the cancer progression, Zhang et al. [20] utilized a three population Lotka-Volterra type model to investigate solely the competition aspect of cancer sub-populations. The model takes the form: and Here, a ij is the competition matrix, where x 1 , x 2 , and x 3 represent AD cells, treatment-sensitive AI cells, and treatment-resistant AI cells, respectively. The PSA level is given by the variable P. The model is used in concurrence with their clinical study. Their work comes at a time where there is no conclusive evidence of whether IAS is better than CAS in terms of prolonging the time to the relapse of cancer. They argue that the result of previous clinical comparisons of intermittent and continuous treatments that study whether intermittent therapy gives any benefit in terms of delaying the onset of treatment resistance [58,59] are not well supported. Zhang and colleagues were able to offer clinical evidence that the previous protocol for intermittent therapy, in essence, has the same effect as continuous therapy. Furthermore, they show that by changing the standard protocol, they were able to obtain significant improvement in the delay of treatment resistance for intermittent therapy with abiraterone over previous results. Their results show that AD cells are likely to have a significant competitive advantage over AI cells when treatment is not applied. In the same work, they support an alternative form of therapy, namely adaptive therapy. The main idea behind adaptive therapy is not to cure cancer but to obtain a stable tumor burden and maintain it indefinitely [60]. Compared to standard treatment where a stable tolerable dose is used, adaptive therapy uses varying drug dosages that change, depending on the response of the tumor to the drug, resulting in a lower overall drug usage than the standard treatment over the same treatment period. By using a smaller dosage, the treatment does not wipe out the drug-sensitive population. Hence, the treatment can be said to rely on the sensitive cancer population to control the resistance population due to their competitive advantage. West et al. (2018, 2019). Subsequently, the work by Zhang et al. was extended further in West et al. [61], where an additional compartment was added in order to implement multi-drug therapy. By noting the limitation of parameter identification, West et al. put forward an assumption that aids the parameter fitting process: parameters corresponding to cancer cells are similar across different patients. What distinguishes the cancer progression between different patients is the initial relative size of the cancer subpopulations. Using this assumption, West and colleagues parametrized the parameters and fixed them across patients. They then carry out the data fitting process only on the initial size of each cancer subpopulation. They argue that the identifiability of their parameter estimation process is supported by the theorem in Sontag's paper [62]. Their work shows that the incorporation of an additional drug (adding of docetaxel to abiraterone) may further delay the onset of treatment failure. A similar concept is demonstrated in the work by West et al. [63] but for chemotherapy in prostate cancer.

Peng et al. (2016).
While many of the mathematical modeling efforts for metastasized prostate cancer focus on hormonal therapy, several researchers have also taken note of the recent development of immunotherapy (vaccination) for prostate cancer. One prominent example is Provenge (or sipuleucel-T), which has been approved by the Food and Drug Administration for treating prostate cancer. A notable example is the work by Peng et al. [64] in 2016. In this work, Peng and colleagues construct a system of differential equations consisting of castrate resistant and sensitive tumor cells and the immune micro-environment. After parametrizing the models using mouse data, they study the treatment efficacy of combination therapy between four different treatments (castration, vaccination, cytokine interleukin-2/IL2 neutralization and regulatory T cells/Treg depletion). The concurrent use of castration and vaccination is motivated by the review of Ching et al. [65], providing evidence that ADT can increase the efficacy of immunotherapy. On the other hand, the incorporation of IL2 neutralization and Treg depletion is motivated by a study [66], showing that ADT is followed by the activation of cytotoxic T lymphocytes (killer T-cells/CTLs); however, the production of IL-2 and Treg may inhibit the activity of CTLs in the prostate lymph nodes. Their work shows the potential of using system biology type models to address complex multi-drug approaches for prostate cancer. While their study shows the potential of a system biology type approach in the modeling of prostate cancer, their limitation lies in the lack of a sufficient source of data for validation of their complex model. [67] propose an alternative approach to modeling immunotherapy motivated by the work of Kirschner and Panetta [68]. Aside from the AD and AI cells, the model also considers the number of activated T-cells, the concentration of cytokines, the concentration of androgen, and the number of dendritic cells. In computational studies, they show a small advantage of combining IAS with immunotherapy over CAS with immunotherapy. Subsequently, Rutter and Kuang [69] extend their model to study the effects that drug dosage amounts and frequencies of administration on the time to relapse. Their computational and mathematical analysis of the model show the possibility of tailoring the vaccine dosage to the patient's effective immune system to maximize the effectiveness of treatment. In another work, Kronik et al. [70] formulated a simple mathematical model that describes the immune response in prostate cancer patients receiving immunotherapy. Kronik and colleagues fit the model to individual patient data and carry out the forecasting using the estimated parameters. Their model produces robust predictions with regards to the data set of 15 patients used for validation.

Elishmereni et al. (2016), Stura et al. (2016).
The idea of predicting treatment failure time is also presented in the work of Elishmereni et al. [71], where a simple model of cancer growth that includes three types of cancer cells with testosterone is used for hormonal therapy. The model is fitted to all patient data. Then the parameters for an individual patient are obtained using Markov Chain Monte Carlo (MCMC) from the distribution obtained from the parameter estimates of the full cohort. Using the individual parameters, 1000 simulations are carried out to predict the biochemical failure, or the onset of castrate resistance, timing for each patient. The model shows high accuracy in predicting the timing of biochemical failure for ADT. In preference of a simpler model for predicting time to cancer relapse in prostatectomized patients, Stura et al. [72] used a form of the generalized Von Bertalanffy growth law to model prostate cancer growth. Using statistical analysis, they highlight the importance of the growth parameter for PSA as a means to predict cancer relapse. Interestingly, they also note that this growth parameter is larger in the case without androgen deprivation therapy as opposed to with androgen deprivation therapy-a conclusion that is supported by other research.

Swanson et al. (2001).
Mathematical models of prostate cancer rely on clinical and experimental data for their validation. For clinical purposes, a model is often validated using the byproduct of prostate cell activity, PSA. However, at the patient level, the correlation between PSA level and cancer volume is questionable. Motivated by this observation, Swanson et al. [73] examined how well PSA represents the prostate (both healthy and cancerous) volume.
This simple differential equation assumes linear production rates of PSA based on the proportion of cancerous/healthy cells coupled with a linear degradation rate. y(t) represents the serum PSA with linear production rates β h and β c from the healthy and cancerous cancer cells, respectively. The volume of healthy cells is assumed to be constant (V h is constant), while the volume of cancerous cells is governed by exponential growth (V c = V 0 e ρt ). The parameter values are estimated using human-derived mouse xenograft LuCaP 23 published in Ellis et al. [30]. Although there are limitations in the biological basis of the model, it offers a valuable explanation for the abnormality in the PSA representation of prostate tumor volume. They conclude this abnormality can be explained by the ratio of the PSA degradation rate and tumor growth rate.

Vollmer et al. (2002), Vollmer and Humphrey (2003).
On the other hand, by adjusting the parameter values in Swanson et al. to account for differences in mouse and human, Kuang et al. [15] (section 5.3.1) argue that while variation in cancer growth rates can explain some of the poor correlation, the key parameters are the PSA production rates by cancer cells and PSA degradation rates. This conclusion is further supported by the work of Vollmer and Humphrey in 2003 [74]. Due to the limitation of PSA as a biomarker of prostate cancer, further investigations have been carried out to explore alternative measurements to PSA in predicting cancer progression. For example, the work by Vollmer et al. [75] in 2001 on the dynamics of PSA during watchful waiting show that PSA amplitude and relative velocity are better predictors of cancer progression-a conclusion supported by clinical studies [76].

Dimonte (2010), Dimonte et al. (2012).
Existing modeling studies sometimes use additional clinical measurement besides PSA as either input or validation for the model. One such example is the work by Dimonte [77], where the author constructed a cell kinetics model to track prostate cancer progression from diagnosis to final outcome. In this case, the patient's Gleason score is used to obtain the transition rate of cells. The model was used to study the author's and another patient's data. In a subsequent work, Dimonte and colleagues [78] simplify the model and use it to explain the variability in the recurrence time of prostate cancer patients. They reach a similar conclusion to use PSA doubling time to improve predictive power. There are various potential biomarkers for prostate cancer that are being studied clinically. For a list of potential biomarkers for prostate cancer, the readers are referred to [79]. Modeling work aiming to incorporate multiple biomarkers can potentially result in increasing the accuracy of the model's predictive power. Lorenzo et al. (2019), Farhat et al. (2017), Liu et al. (2015). Aside from the aforementioned work that focuses on hormonal therapy, there are other works that utilize mathematical models to study other aspects of prostate cancer. For instance, Lorenzo et al. [80] use mathematical models to study personalized treatment with radiation therapy, where they suggest several potential prognostic measurements that can be obtained from the model. A recent study by Farhat et al. [81] formulates a mathematical model of metastatic prostate cancer while taking into account the bone micro-environment to investigate several possible therapeutic strategies. Both modeling approaches are novel; however, their work lacks substantial validation to support their findings. Taking on a more computational approach, Liu et al. [82] use a hybrid automaton model to find a personalized therapeutic strategy. . Furthermore, while existing modeling work for prostate cancer relies strongly on deterministic models, stochasticity in the model has also been considered numerically by Tanaka et al. [83] and Cerasuolo et al. [53]. Furthermore, analytical consideration of stochastic modeling has been carried out by Zazoua and Wang [84]. Their results show that the stochastic models can capture the statistical components of the dynamics of PSA time serial data. Additionally, extensions using delay differential equations and fractional differential equations for prostate cancer modeling have been carried out by Baez [85] and Mizrak et al. [45], which show some improvement in data fitting. Although these approaches show early promise, follow-up studies are needed to verify and expand their implications.

Mathematical Models in Clinical Settings
Many modeling efforts for prostate cancer diverge from one another, where previous results are often not utilized. In this section, we discuss the limitations of current work in a clinical setting and how to potentially address them. A summary of our discussion is presented in Table 2.

Real-Time Estimability
Mathematical models need to be validated against data prior to any application. In theoretical studies, the complete set of patient data is often known to the researchers and is used for model validation purposes. This makes sense because existing mathematical models often contain a large number of parameters, so having more data gives a better chance of estimating the parameters. A study of model identification methods by Hirata et al. [86] shows that 1.5 cycles, where a cycle refers to an onand off-treatment period in IAS, of data is the minimum requirement for model identification in most cases. Yet in clinical settings, a model should be useful even when a limited amount of information is available; however, this is often not the case. If just a few data points are used to parametrize a model, its predictions would be highly unreliable [49]. On the other hand, if researchers wait for more data as the treatment progresses, it may be too late, perhaps because the cancer becomes resistant and the patient must switch to a different treatment. Furthermore, this issue affects even the classification systems that are introduced along with some mathematical models [40,87].
One might assume that obtaining data more quickly, perhaps by way of a self-measuring device or enticing the patient to visit the hospital more regularly, would resolve this issue; however, this is deceptive. The shape of the patient's data trend (PSA or androgen) is unique. Thus, forcing a collection of data quickly does not expose completely the shape of this trend; instead, this would only oversample certain sections leading to biases.

Real-time estimability
The estimation of parameters in mathematical models often require a large quantity of data. However, the nature of data collection in real-time means that reliable estimation of parameters for patients may not be possible at the early stages of treatment.
Some parameters share similar values across patients, while others are more patient-specific. This distinction should be studied in detail. Utilizing multiple data sets is another possibility to allow early estimates of parameters. Finally, parameter evolution can be accounted for to address the limitation of data availability. Local sensitivity analysis and uncertainty quantification should be studied for each patient. Clear links between each parameter and its physical interpretation should be established, which potentially allows for laboratory estimates/bounds to resolve identifiability.
and identifiability Section 4.2 Optimal schedule, Studies on optimal schedule and treatment yield useful information on how intermittent, adaptive, and combination therapies should be carried out. Patient classification based on treatment effectiveness can be done using model parameters. However, both aspects are affected heavily by the estimability of the parameters and the uncertainty in the model's forecasts.
The usefulness of optimal studies and classification hinges on how well the uncertainty in the model can be quantified, which relates to previous issues. In addition, the objective of optimal studies may be extended to include drugs' properties, cost, and important features of each treatment.
optimal treatment, and patient classification Section 4.3 There are two proposed ways to address the issue of data limitation in practice. West et al. [61] noticed the issue of the high ratio of the number of parameters to data points, so they assumed that the parameters in the model are the same across patients and the only difference is the initial population of cancer cells. If this assumption holds, the number of parameters needed to be estimated for existing models would decrease dramatically-allowing mathematical models to be useful in clinical studies. However, this is generally not the case. Various studies on personalized medicine, or the idea that each patient should have his own unique set of parameters, show that the parameters vary significantly across patients [88]. This variability between patients makes intuitive sense. The characteristics of a tumor and how the patient reacts to certain drugs should depend on the specificity of the patient's physiology and the tumor's composition, metastatic site, etc. However, it may be possible to show that their assumption holds for patients within a certain group, for example, patients who share certain physical traits. Additionally, some parameters may in fact be relatively constant among patients. Of course, if such categorization is possible, it could potentially be used to resolve this issue.
The second possibility is not as clear cut as the first; however, it is inspired by classical physical studies. In this approach, we could consider models that contain mostly (if not only) parameters that can be measured or estimated from laboratory testing. Thus, most of the parameters would not come from data fitting but from actual physical testing of the patients. However, this method has its own set of problems. Perhaps the biggest problem is the construction of the model. One cannot mindlessly add all known mechanisms/factors into the model, There are many limitations by doing so, such as increasing the complexity of the model, existence of unknown biological details, not well tested biological details that may turn out to be incorrect in the future, etc. In order to create such a model, the model formulation rests on the expertise in both mathematical modeling and biological knowledge. One compromise using this approach is to carry out studies that focus on multi-level fitting as done in previous work. This would allow for more accurate estimations of certain parameter values. However, one must carefully isolate the parameters with respect to the data to avoid ambiguity in the biological interpretations of the parameters.
An alternative approach to resolve the limitation of data in clinical settings is to account for the evolution of parameters as the treatment progresses, instead of relying extensively on data fitting. The rationale is that the parameters in dynamical models are often stand-ins for more complex processes, which means their values represent the average values over a certain time interval of the underlying processes. Hence, as time goes on, there would be changes accumulated from the underlying processes of the parameters. The effect of the underlying process may be minimal in many cases, but it could also be substantial. This is a complex subject that is akin to a time-scale analysis. Sometimes, fast changing processes can be accounted for using an average as long as the interval of time that the average is taken over is sufficiently long enough for the effect of the fast changing process to be negligible. On the other hand, other processes take a long time for their effects to be noticeable. When it comes to parameter estimations, there should be consideration for the time interval used for data fitting. For example, the extended version of the BK model (Phan et al. [50]) contains a parameter that stands in for the maximal level of serum androgen. While they take this as constant, the observed maximal level of androgen as treatment goes on tends to decrease with each cycle, see Figure 6. The reason for this phenomenon is not known; however, it is suggested to be due to an accumulated damage from drugs or the patient's behavior changes over time after learning of the cancer. Furthermore, Phan et al. [50] demonstrate that by focusing on the more recent data using a time-weighted objective, significant improvements on the fitting and forecasts can be obtained. This observation supports the implementation of parameter evolution within dynamical models. The first study to directly incorporate this approach into the prostate cancer model treats the parameter associated with the resistance of cancer as a dynamic variable [47]. This not only simplifies the model but also results in better fitting compared to more complex models. From a computational perspective, the evolution of parameters can be thought of as updating the parameter values as new data becomes available by using Kalman filters as shown in Wu et al. [49]. The concept of parameter evolution over time recently was incorporated in Brady et al. [89] and is explored conceptually in a more general ecological context by Loladze [90]. The data is taken from Bruchovsky et al. [91]. Note that the trend of the maximal level of androgen goes down over the course of treatment.

Uncertainty, Identifiability, And Sensitivity
Bootstrapping and the ensemble Kalman filter are two standard methods that have been used to study the uncertainty in prostate cancer models [41,49]. As shown in those works, existing models often contain large uncertainty in their prediction, especially when longer forecasts are needed to make a decision on optimal schedule, see Figure 7a. Other variance-based uncertainty quantification methods can also be applied for similar purposes, such as the work by Elishmereni et al. [71]. In the case of sensitivity analysis, a standard approach is to vary the parameters one by one by some percentage and evaluate the corresponding changes at a fixed time point [47]. This method can also be extended to study the sensitivity over the entire treatment [50,92]. The advantage of studying the sensitivity of a parameter over time is evident in intermittent and adaptive therapy. Since there are orders of magnitudes difference for certain variables (such as PSA level) for different phases of the treatment, the sensitivity of the parameter is also affected, see Figure 8. Thus, having a better understanding of how the sensitivity of each parameter changes as the treatment progresses can provide a better tool to optimize treatment. Additionally, sampling-based methods to account for simultaneous effects of all parameters is also possible; however, this should be done in a relatively small range after the parameters have been identified for a specific patient.
There are a variety of methods to estimate the parameters of a model. However, recently, when the identifiability of some prostate cancer models were examined, Wu et al. [49] show that these models are not practically identifiable. They further show that for an unidentifiable model, different sets of parameters may yield an indistinguishable fitting but result in vastly different forecasting, see Figure 7b. Such a result is troublesome as it could undermine the applicability of mathematical models in a clinical setting. To take into account these issues, some researchers rely on the 2n + 1 law by Sontag [62]. However, this law requires models to be structurally identifiable-a condition that is often not examined when applying the law. Alternatively, if one carries out the sensitivity analysis of a model before hand, then information revealing the least sensitive parameters could be used to enhance the identifiability of the model. The issue of identifiability ties back to the complexity of the model and the limited data in clinical settings. It essentially requires a sufficient number of data points and data types, to uniquely determine the parameters of a model. In Wu et al. [49], the authors suggest the use of an observer experiment to quantify the minimum required data that would allow the model to be identifiable. The suggested technique is to use the Fisher information matrix because of its ability to test different combinations of data sets easily. While the process is tedious, such information can provide insight into the types of data needed for a model to be useful. Additionally, if the functional forms of some parameters are available (either through study of parameter evolution or study of the nonlinear relationships between parameters), it may directly address the issue of model identifiability [93].

Optimal Schedule and Patient Classification
Due to the lack of a gold standard when it comes to determining a treatment schedule for IAS, its full benefits may not be realized [20]. Furthermore, in some instances, patients may benefit more from CAS. To tackle this issue, Hirata and colleagues studied patient classification based on whether a patient would benefit more or less from IAS as compared to CAS [40]. Another classification system was introduced by Morken et al. [87], where the type of treatment resistance is studied based on cell death rate analysis. In treating prostate cancer, mathematical models have been used to study potential optimal treatment schedules that may give the best chance for patients. While various studies have been done on this topic, the goal of an optimal schedule has varied among them. For instance, Suzuki and Aihara [94] choose the objective to be the minimization of the amount of time a patient is on treatment while still keeping IAS effective in controlling the cancer progression. Hirata et al. [42] instead focus on delaying the relapse of cancer as long as possible, and Cunningham et al. [95] study three potential objectives: minimizing average tumor volume, tumor mass variance, or average density of androgen independent cells. While these are obvious objectives to be minimized for prostate cancer (or cancer in general), other creative alternatives such as minimizing cancer activity at any given time, or PSA doubling rate, may provide better control of cancer.  The study of drug combinations can be considered as a subset of the optimal schedule. Various studies have been done on finding the best combinations of drugs for certain cases [52,61]. While these studies hold promises, their implications hinge on how well the mathematical models represent the underlying biological system, which means drawing definitive conclusions is difficult due to the aforementioned problems with model validation. Furthermore, finding an optimal time (free terminal time) for a treatment may aid the design of a treatment, especially in the case of adaptive therapy. A comprehensive collection of tutorials for optimal control application in bio-sciences are presented by Lenhart and Workman [96].

Data
Various relevant data sets exist, ranging from experimental studies of cancer cells to clinical studies of patients undergoing hormonal therapy. In this section, we summarize some of these data.
Perhaps the most used data set for prostate cancer model validation at the patient level comes from a clinical trial at Vancouver Prostate Center. The study admitted patients who showed a rising serum PSA level after undergoing radiotherapy without evidence of distant metastasis or being previously subjected to hormonal therapy, with the exception of less than three months of neoadjuvant hormonal therapy [91]. Additionally, all patients exhibit high serum PSA levels (≥ 6 µg/L) prior to the study. Another set of data comes from a clinical study [46], where three stage C and four stage D prostate cancer patients were treated with intermittent hormonal therapy for 21 to 47 months (two to four cycles). All patients show decreasing PSA levels during the course of the study. Moreover, there are various clinical studies that have been used in mathematical modeling. For instance, patients' data sets were used in the study by Draghi et al. [88], which can also be extracted directly from their paper. Another study [41] also utilized additional data sets from different clinical studies in the United States and Japan [97][98][99][100]. Additionally, we would like to point out a clinical trial [20,61] that was the first for prostate cancer that utilizes adaptive therapy to treat metastatic castrate-resistant prostate cancer. However, as this is ongoing, the data may not be readily available. For a review of clinical studies that focus more on the statistical correlation of different variables for prostate cancer, we refer the readers to the work by Dimonte [77].

Parameter Ranges
One of the main issues that is common to existing work is determining the values of model parameters. Here, we also provide references with regard to some important parameters.
Cancer maximum proliferation rates: by making the assumption that measurements of cell growth and death for hormonally untreated patients to be that of AD cells, the study of tumor doubling time at various stages in Berges et al. [101] can be used to suggest that the range for the proliferation rates of AD cells is (0.004-0.081)[day] −1 . Similarly, we can find a range of (0.001-0.046)[day] −1 for AI cells. Note that these ranges fall within the expectation that AD cells out-compete AI cells in an androgen-rich environment. Additionally, growth rates for specific cell types can be obtained in vitro [20,61]. However, appropriate scaling or fitting is necessary to account for differences in tumor environment. Since the environment in a laboratory is ideal for the growth of cancer cells, their growth rate should be lower in practice. Thus, we consider these ranges to be for maximum proliferation rates.
Cancer cell death rates: similar to cancer growth rates, we estimate the ranges of cancer death rates for AD and AI cells to be (0.001-0.0525)[day] −1 and (0.015-0.0775)[day] −1 , respectively.
Cancer cell maximal transformation rate: Robust estimate of cancer transformation rates from experimental data is lacking.
However, numerical experiments and fitting show that (10 −5 -10 −4 )[day] −1 is an appropriate range for the maximal transformation rate [24,33]. The transformation rate encompasses the mutation rate. While there are several estimates of the PCa cell mutation rate, mathematical models often use transformation rates instead of the more specific mutation rate.
PSA clearance rate and production rates by healthy and cancerous cells: initial estimates [73] were for human-derived mouse xenograft sublines LuCaP 23.1, 23.8, and 23.12. These values were later extrapolated for human prostate cancer in Section 5.3.1 [15] using the study of Berges et al. [101] and Vesely et al. [102]. The PSA production rate from healthy cells is ( [73]. Additionally, the PSA clearance rate is within the range of (0.1754-0.4030)[day] −1 [103]. We summarize these in Table 3. Table 3. Ranges for some commonly used parameter values in mathematical models for prostate cancers. Note that for more accurate ranges, the specificity of the situation needs to be taken into account, for example, sublines of cancer cells.  [73] Some parameters are difficult to estimate robustly, for instance, the competition rates of cancer cells. In the case of difficulty in estimating certain parameters from laboratory data, researchers also rely on experts' opinions or fix the parameters in an ad hoc way [20,43]. Alternatively, the model can also be derived from physical laws (e.g., conservation laws), which introduce competition rates without adding additional explicit parameters for competition between cells [47,48,50]. For readers interested in how some parameters are estimated in existing models, we refer to studies [16,51], where binding rates are estimated from multi-level processes.

Conclusions
Mathematical models play a vital role in the study of prostate cancer dynamics. Over the past two decades, various models have been developed to examine different aspects of prostate cancer, including treatment options and schedules. In this review, we collect and synthesize the results of existing studies to share conclusions that are well agreed upon and raise questions that are in need of further investigation.
The study of prostate cancer at the multi-cell level gives insights on the treatment efficacy with respect to the competition among the cells. Specifically, the consensus among existing work is that IAS is superior to CAS at delaying the relapse of metastasized cancer, if AD cells have some competitive advantages over AI cells in an androgen-rich environment. While some studies have shown that AD cells can indeed out-compete AI cells, there is no extensive evidence that this is true in general. Instead of suggesting that the competition advantage is an intrinsic property of AD cells, this property needs to be examined on a case-by-case basis because many factors such as the physiology of the cancer or the specific cancer phenotype can affect the relative competition between cancer cells. Additionally, we find that lowering the androgen level as a preventive measure of prostate cancer may result in selection for aggressive androgen independent phenotypes, if the tumor forms. On the other hand, studies that focus on tools used to track prostate cancer progression, such as PSA, show that while PSA is an overall good biomarker for cancer growth, it can display poor correlation in many cases due to the patient's specific PSA production rates by cancer cells and PSA degradation rates. Instead, the relative velocity of PSA appears to be a better alternative for monitoring tumor progression. Furthermore, to tackle the problem of identifiability of models, multiple biomarkers should be used to enhance the accuracy of models.
Thus far, there exists mathematical models for almost all clinical stages of prostate cancer. Some models are even flexible enough to be used across many stages. However, the application of each model is still limited to certain stages. Having a comprehensive model that can describe tumor development as a whole may be useful for clinical practice. Treatment effectiveness is one of the most examined issues using mathematical models with a focus on hormonal therapy. This is often coupled with the study of optimal schedule, especially for IAS and adaptive therapies. While the studies so far show potential, they lack key components for use in clinical settings. First, the drugs being considered in a mathematical model are often generalized-meaning the specific properties of the drug are not accounted for. Secondly, for combinations of drugs, the actual dosages and toxicity are not considered. This is perhaps due to the fact that previous studies are dominated by the traditional way of administering hormonal treatment (fixed dosage). However, as adaptive therapy is becoming more prominent, the actual dosages and toxicity of multi-therapy should be accounted for in modeling efforts. Furthermore, in the case of optimal schedule, the best objective is often not clear, which is perhaps due to the lack of a credible set of biomarkers that can track cancer growth.
Perhaps the most urgent issue with existing efforts in modeling prostate cancer is the validation of mathematical models. In the early stage of development, models are validated using observations of clinical trends, which can only give qualitative descriptions of the process. As the field progresses, there have been more attempts to validate models against clinical data in a quantitative way. While this is an improvement, efforts are lacking and parameter identification becomes a big problem for the validity of many mathematical models. This issue is further evident when mathematical models attempt to reach clinical application because of the lack of retrospective data sets. The identification of model parameters is coupled with issues of uncertainty, sensitivity, and stochasticity in making any conclusions or forecasts from mathematical models, while useful in providing an indication of likely scenarios, problematic. As many of the existing models are phenomenological, the issue of model identification is further emphasized. Instead, we suggest the use of mechanistic models, which tend to be more complex, but parameters can often be linked to experimental data. The formulation of mechanistic models should be based on well-tested structures from phenomenological models, while acknowledging the biological details. In this way, previous results can aid future studies.
As optimal treatments require model predictions far into the future, it severely limits the usage of mathematical models due to increasing uncertainty, especially for long term treatments such as hormonal therapy. To address this issue, we observe that parameters in dynamical systems are often an average of an underlying process over a period of time (e.g., the period of time used in parameter estimation). Thus, the parameters in dynamical models should be updated regularly. There are several means to incorporate parameter evolution in dynamical models, such as applying the Kalman filter [49], incorporating the dynamics of the parameters into the model [47], using weighted fitting to emphasize more recent data [50], or studying the trends of parameter evolution [89]. If the trends of certain parameters can be established for a specific treatment, it can allow for better prediction without higher certainty. An in-depth exploration of parameter evolution could significantly bridge the gap between mathematical models and clinical practices. Furthermore, it could provide a new set of interesting dynamical questions [104].
In summary, mathematical modeling efforts for prostate cancer so far have shown that (1) IAS or adaptive therapy is better than CAS at delaying cancer relapse, if AD cells hold competitive advantages over AI cells, and (2) cancer progression is better associated with the relative velocity of PSA due to variation in PSA production by the tumor and its degradation rate. Moving forward, to close the gap between modeling work and clinical application, the following questions need to be addressed: (1) how to obtain the individualized model parameters without extensive data, (2) what are the most important parameters that affect cancer growth and how to use this knowledge in clinical applications, (3) what is the most appropriate objective for an optimal treatment/schedule study, and (4) how to best quantify the uncertainty in model prediction? Many physicians have found value in mathematical models for improving our understanding of prostate cancer progression and creating better treatment for patients. However, many issues stand in the way of a complete theory and a clinically applicable model. In this paper, we have highlighted the general consensus among existing studies and limitations that should be addressed in the future. We hope that this review can aid future attempts at studying prostate cancer. Moreover, we note that while the findings and issues emphasized in this paper focus on various aspects of mathematical modeling of prostate cancer, they are also relevant to other cancers and their respective treatments. For instance, breast cancer shares many similarities with prostate cancer. The issue of modeling treatments in the case of co-occurring cancers is not well explored and we are not aware of any attempt to model treatment of prostate cancer in the presence of another cancer. However, since prostate cancer is more prevalent in older men (55+), when cancer is more likely and detectable, modeling attempts in the case of co-occurrence of cancers can be beneficial. Thus, we hope this review can provide insights into future modeling work in other areas of mathematical oncology.