Dynamics and implications of models for intermittent androgen suppression therapy

1.   Introduction

Prostate cancer (PCa) is a major health concern in the US. It is the second most common cancer in men with the chance of life-time diagnosis being 1 in 7. PCa is more common in men over 50, and especially over 65. Recent developments in prostate cancer treatments have increased the rate of treatment and survival rate greatly, evidenced by a 5/10/15-year survival rate of 99%/98%/95%. However, for metastatic prostate tumors, the 5-year survival rate is 28%, which is among the lowest of all cancers ^[1,25]. The only known risk factors are age and heredity, while race and living situation are possible risk factors ^[1,27].

Most prostate cancers are found during screening with a prostate-specific antigen (PSA) blood test or a digital rectal exam (DRE). The fact that PSA level can reflect the size of the tumor and act as a predictor of relapse makes PSA the standard biomarker for PCa. Furthermore, PSA is frequently used as a proxy for the progression of a tumor, especially in mathematical modeling ^[14,20]. Androgens like testosterones and 5$\alpha$-dihydrotestosterones (DHT) are released by prostate glands in males with over 90% produced by the testes and the rest by the adrenal gland. These androgens bind to the androgen receptors and are transported to the nucleus. In the nucleus they bind to androgen responsive elements, which results in the transcription of androgen-regulated genes. This can lead to the proliferation and inhibition of apoptosis in prostate cancer cells. For this reason, the standard hormonal therapy for a patient with advanced PCa is androgen deprivation therapy (ADT): a chemical castration process which is preferred by patients over the surgical process of removing the testes ^[10,11]. This treatment is based on Huggins' and Hodes' Nobel Prize-winning discovery that castration causes rapid regression of PCa. Initially, most of the tumor cells are androgen dependent, but during treatment, they mutate or adapt in response to the low-androgen environment. Once the population of such androgen independent cells reaches a certain threshold, the tumor can continue to grow, and this phenomenon is termed treatment resistance ^[11].

For ADT, the reduction in androgen during treatment can cause undesirable side-effects such as impotence, depression, bone demineralization, and an increased risk of dementia. As a counter measure, intermittent androgen suppression (IAS) was introduced in the mid-1980s. Under IAS, the patient is put on ADT until his PSA level decreases to a predetermined lower threshold, which often corresponds to tumor regression and the PSA thresholds can be varied by the clinicians. Then, ADT is temporarily stopped until the PSA level rises to another predetermined upper threshold ^[1]. IAS has been shown to improve quality of life for patients; however, it remains controversial whether IAS is superior to continuous androgen suppression (CAS) in prolonging patient survival ^[4].

The lack of a solid understanding and a gold standard in treatment of prostate cancer drives the need for mathematical models. Especially in the last 15 years, many mathematical models have emerged to help explain the progression dynamics of prostate cancer in hope of answering some of the aforementioned questions ^[20].

Jackson ^[18] developed the first experimentally driven PCa model under CAS in 2004, which inspired later modeling efforts. In 2008, Ideta et al. ^[17] introduced the first mathematical model for PCa under IAS treatment with androgen dependent (AD) cell to androgen independent (AI) cell mutation. Guo et al. ^[13] formulated a partial differential equation version of the Ideta's model. Shimada and Aihara ^[24] studied the competition between different populations of PCa cells. A different approach to the formulation of a competition model based on Ideta's work was proposed by Yang et al. ^[28]. Eikenberry et al. ^[7] extended Ideta's model with intracellular signal of androgens to study its role in PCa. This was extended by Portz et al. ^[22] to a model that can fit clinical PSA data. And Portz's model was later simplified by Baez and Kuang ^[2] to also fit androgen data. Another extension of Ideta's model and the ideas of Eikenberry et al. ^[7] was presented by Jain et al. in ^[26], whose model captures the detailed biochemical dynamics of PCa. On the other hand, Hirata et al. ^[15] constructed a piecewise linear model based on the phenomenological behaviors of PSA to fit and describe the dynamics of PCa.

Additionally, numerous attempts have been done to study various clinical aspects such as an optimal schedule using an adaptive threshold and a stochastic hybrid automaton model ^[12], or a discrete hybrid cellular automaton model for a similar purpose for bone metastasized PCa ^[5]. Elishmereni et al. developed an algorithm that takes a mechanistic model in combination with the tumor Gleason score, a classification of the severity of prostate cancer, to predict when castration resistance will occur ^[9].

In this work, we rigorously exam and extend two previous approaches, which were developed by Baez and Kuang ^[2], Hirata and his colleagues ^[15], to further our understanding of the dynamics of prostate cancer and provide better tools for clinical applications.

The data we use to study the models comes from Vancouver Prostate Center ^[3], which contains the information of 109 patients. Similar to ^[2], we require the patient data to have at least 20 data points for both androgen and PSA in the first 1.5 cycles, the number of eligible patient data for comparison purpose is reduced to 62. Additionally, we need sufficient data for an additional treatment cycle for the forecasting comparison, so we select 26 of these 62 patient data.

2.   Model and Method

Our model is built upon previous models which were based on a more complex but mechanistic framework. This is necessary because unlike phenomenological and statistical models, mechanistic models are created with sound hypotheses or well-known laws observed in nature ^[19].

2.1.   The BK model

Baez and Kuang ^[2] recently introduced an improved version of the model in ^[22] with the main emphasis of developing a simple model that can be readily used by physicians. This model is the first one to successfully take into account of clinical androgen data via fitting and forecasting. This is an important feature since the incorporation of androgen data fitting allows for better fitting and forecasting. In this model, $x_1$ and $x_2$ represent the AD and (strong or irreversible) AI cell populations respectively. $Q$ is the androgen level and $P$ is the clinical PSA level. $\lambda(Q)$ represents the mutation (or adaptation) rate from AD to AI cells, which is dependent on the androgen level. The lower the androgen level, the higher this mutation/adaptation rate is.

One important parameter constraint in the BK model is the inequality involving the cell quotas of AI and AD cells, $q_2 < q_1$, which represents the fact that AI cells can proliferate in an environment with a low level of androgen better than AD cells. We also want to clarify the common usage of the terms here. AI cells also depend on androgen for growth, but the reliance on androgen in the body/prostate is smaller in comparison to that of AD cells. Note that in the BK model, the AI population is non-reversible, which will be referred to as the strong AI cells. Moreover, the death term expression is adapted from that of Morken et al. ^[21].

2.2.   The three population model

The validation of a model goes beyond how well it fits the observed data. This is especially true for mechanistic model, which can be used as direct test of well-founded scientific theories. Although the BK model has an impressive fitting and forecasting power ^[2], it may not have the appropriate underlying population assumption. In other word, the population structure in the BK model may not be flexible enough to allow for different cancer dynamics. This is especially crucial when the correct prostate cancer cell population structure and interaction are not yet fully known. Thus, this population structure simplicity may compromise the BK model's ability to appropriately represent the dynamics of the cancer growth and its forecasting potential.

As noted by Hirata et al. ^[15], a system with three distinguished populations is perhaps optimal in the case of metastasized prostate cancer because such a system --- under linear rates of change --- can reproduce the bi-phasic decline that is observed in the PSA data when the treatment starts.

The difference between the two population model in ^[2] and the three population model in ^[15] is an intermediate population between an androgen-dependent (AD) cell population and a strongly androgen-independent (strong AI) cell population. We call this new cell population the weakly androgen-independent (weak AI) cell population. Unlike the strong AI cell population, this weak AI population has the ability to revert back to the AD cell population. Due to permanent nature of the transformation, we consider the transition from $x_1$ to $x_3$ mutation, while the transition from $x_1$ to $x_2$ adaptation. Thus the three population model has a similar population structure as in ^[15]. The underlying population structure of this three population model and the parameter definitions are presented in Figure 1 and Table 1, respectively, and the system of equations takes the form below:

where "w.AI" and "s.AI" stands for weak AI and strong AI cells. Due to insufficient data to support differentiating "mutation/adaption" rates in prostate cancer cell, we take the "mutation rates" to be the same, e.g. $\lambda_{ij}(Q) = \lambda(Q)$, or $c_{ij} = c$. It is worth noting that the density dependent death term $\delta_{11} x_1$ is a result of assuming interspecies effect is negligible in cell population growth. This is because cells of the same tend to cluster.

2.3.   Model derivation

We adapt the approach of the BK model formulation to derive our three population model. Recall that the three-population-model contains three different types of prostate cancer cells: the androgen dependent type ($x_1$), the reversible (or weak) androgen independent type ($x_2$) and the irreversible (or strong) androgen independent type ($x_3$). The proliferation rates $\mu_i, \, i = 1, 2, 3,$ of the three populations are assumed to be governed by the Droop cell quota equation, where we assume androgen is the limiting nutrient ^[6]:

where $\mu$ represents the maximum proliferation rate, $Q$ is the androgen cell quota, $q_i$ is the minimum required androgen level that would allow the respective cancer cell type to proliferate. Since we know that androgen dependent cell proliferation requires higher androgen levels, we require that $q_1 > q_2$ and $q_1 > q_3$, while $q_2$ and $q_3$ may have similar values ^[11]. The death rate is assumed to also depend on the androgen level and it takes the following form:

where $R_i$ is the half-saturation level and the apoptosis rates $d_1$ is greater than $d_2, d_3$ under the assumption that androgen deprivation therapy has less effect on androgen independent cells than on androgen dependent cells. The mutation rates between different cell populations are assumed to be the same, dependent on androgen level and take the form of a hill function:

where $K$ is the half saturation level and $c$ is the maximum mutation rate. We further include the interspecific and intraspecific density death rate, $\delta_{ij}$, as described in previous section.

We first derive the rate of change for the androgen cell quota $Q$ (nm) for a tumor with a single type of cells. We assume $Q_x$ to be the total androgen inside the tumor $x$ (mm$^3$) and that the androgen cell quota is uniformly distributed in the tumor. This leads to:

where

Let $\gamma$ be the production rate of androgen in serum. Then $\gamma = u(t)\gamma_1 + \gamma_2$, where $\gamma_1$ and $\gamma_2$ are the androgen production rate by the testes and other organs respectively, e.g. $\gamma_1 >> \gamma_2$. $u(t)$ is the indicator of on and off treatment, where $u(t) = 0$ when the treatment is on and $u(t) = 1$ when the treatment is off. Assume the increase in androgen cell quota is via diffusion, then inflow of androgen to the tumor can be approximated as: $\gamma(Q_m - Q(t))x(t)$, where $Q_m$ is the maximum androgen level. Additionally, assume the outflow of androgen in the tumor is mainly due to the death of tumor cells, then it takes the form $\left(\nu \frac{R}{Q+R} + \delta x(t) \right) Q(t) x(t)$, where $\nu$ is the maximum apoptosis rate for androgen. If we assume the conservation law holds, then the rate of change of androgen inside the tumor takes the form:

Differentiate and solve for $Q'(t)$ to arrive at the desired expression. Note that for the derivation of the three-population-model, we do similarly but letting $x(t) = \sum_{i\in\{1, 2, 3\}} x_i(t)$. Then the death rate also becomes

Finally, we assume the PSA level is produced at a baseline rate of $b Q(t)$ by regular cells and $\sigma \sum_{i\in\{1, 2, 3\}} x_i(t) Q(t)$ by cancer cells, while being cleared from serum at rate $\epsilon$.

2.4.   Parameter estimations

For parameter fitting, we use MATLAB (version 2017a) function fmincon, which utilizes interior point algorithm to minimize an objective error function within the range provided in Table 1. The objective error function is calculated as the sum of the mean squared error (MSE) fitting error for PSA and androgen since particularly large error is undesirable in clinical application.

where N is the total number of data points, $y$ is either the PSA or androgen data generated from the model, and $\hat{y}$ is the respective observed data. Note that, similar to ^[2], we assume that the level of intracellular androgen is similar to the serum androgen. The objective error function that is minimized by fmincon is the sum of PSA and Androgen residuals.

The range for the parameter is adapted from ^[2] under the assumption that the parameter constraints for the weak and strong AI cells are the same. The initial guesses of the parameters are chosen by hand within the established range to obtain the best fit.

3.   Analytical results

In this section, we provide some analytical insights into the three population model. Additionally, we explore a condition on the population structure that ensures the consistency with what is observed in clinical studies.

3.1.   Treatment-resistance under continuous therapy

Hatano et al. ^[14] points out one important flaw with the PKN model (the predecessor of the BK model). Under CAS treatment, PKN model fails to reproduce relapse both analytically and numerically within the given parameter constraints. Note that this does not mean conclusively that the model fails intrinsically. In fact, BK model allows for possible relapse under a similar analysis. If the choice of treatment is CAS, then the constant suppression of androgen production will result in a near minimum level of androgen in the body. Thus we can approximate $Q \approx q_1 > q_2$, which heavily limits the growth of the AD cell population. Therefore the BK model can be approximated by

where $T = \mu\left(1 - \frac{q_2}{q_1}\right) >0$ is the positive growth term of the AI cells that, under the right conditions, can give rise to relapse. The linearized system at extinction is then:

There are many ways to force instability at extinction. One such condition arises from having the determinant to be negative, e.g. $T > D_2(Q)$ or $\mu\left(1 - \frac{q_2}{q_1}\right) > d_2 \frac{R_2}{q_1+R_2}$, which is possible under the parameter constraints given in Table 1. Similar conclusion extends directly to the three population model.

Hatano et al. ^[14] pointed out that the reasons behind PKN model not being able to reproduce relapse under CAS is either because it does not have an absorbing state (with respect to the direction of mutation, see Figure 1 for an example of absorbing state) or the parameter constraints are not optimal. However, the BK model is capable of reproducing the relapse within the provided parameter constraints, this suggests that the appropriate structure of the cell dynamics in the case of prostate tumor should have at least one absorbing state.

3.2.   Properties of the three population model

In this section, we provide some standard analysis for our model to justify its biological validity. First, we establish an invariant positive region under biologically reasonable assumptions. Additionally, we show that the population of a cell is also bounded above, which is coherent with reality. While the PSA level is an indicator of the growth of the tumor, it has negligible effect on the growth of the tumor cells and the production of androgen. This is reflected in the mathematical form, where the dynamics of $P$ is dependent on the dynamics of $x_i$ and $Q$, but not the other way around. This means analytical results for $P$ can be inferred from the dynamics of $x_i$ and $Q$. Thus we choose to omit the analysis of $P$ here.

Proposition 1. Assume $q_3 \leq q_2 \leq q_1 < Q_m$ and $\delta_{11} \geq \delta_{22} \geq \delta_{33}$. Then the solutions that starts in the region $\{(x_1, x_2, x_3, Q) : x_i > 0, i = 1, 2, 3, x_1 + x_2 + x_3 \leq \frac{\mu_3(Q_m) - D_m(q_3)}{\delta_{33}}, q_3 \leq Q \leq Q_m \}$, where $D_m(q_3) = \min\{D_1(q_3), D_2(q_3), D_3(q_3)\}$, stays in that region.

Proof. If $x_i(0) > 0$, it is easy to see that each population remains positive as long as it is bounded. Now since

it is clear that $\frac{dQ}{dt}<0$ when $Q = Q_m$. Since $q_3 \leq q_2 \leq q_1,$ one can see that

when $Q = q_3.$ It follows that $q_3 \leq Q(t) \leq Q_m$ for all $t > 0$ and initial condition $q_2 \leq Q(0) \leq Q_m$.

Since $x = x_1 + x_2 + x_3$, $\delta_{11} \geq \delta_{22} \geq \delta_{33}$, we have

This implies the limit supremum of $x(t)$ is bounded by $\frac{\mu_3(Q_m) - D_m(q_3)}{\delta_{33}}$. Since $x(0) = x_1(0)+x_2(0) + x_3(0) \leq \frac{\mu_3(Q_m) - D_m(q_3)}{\delta_{33}}$, it implies $x(t)$, or $x_1(t)+x_2(t)+x_3(t)$, is bounded above by $\frac{\mu_3(Q_m) - D_m(q_3)}{\delta_{33}}$.

The number and locations of nonnegative steady states are natural but challenging mathematical questions that maybe of interest to readers. Here our focus is to exam the existence and local stability of an absorbing steady state where both the AD and reversible AI cells are extinct.

Proposition 2. Assume $q_3 \leq q_2 \leq q_1 < Q_m$ and $\delta_{11} \geq \delta_{22} \geq \delta_{33}$, then $E_1 = (0, 0, \frac{\mu_3(Q_*) - D_3(Q_*)}{\delta_{33}}, Q_*)$ where $Q_* = \frac{\gamma Q_m + \mu q_3}{\gamma + \mu}$, is the only nonnegative boundary steady state for the system with zero as the first or second component.

Proof. Let $E_1 = (x_{1*}, x_{2*}, x_{3*}, Q_*)$ be a steady state of the system with $x_{1*}x_{2*} = 0.$ Assume first that $x_{1*} = 0$ and $x_{2*} > 0$. Then at this equilibrium, $\frac{dx_1}{dt} = \lambda_{21} (Q_*) x_{2*} > 0$, which is a contradiction. Similarly if $x_{1*} > 0$ and $x_{2*} = 0$, then $\frac{dx_2}{dt} = \lambda_{12} (Q_*) x_{1*} > 0$, also a contradiction. Hence we have $x_{1*} = x_{2*} = 0.$ Solving the $\frac{dQ}{dt} = 0$ yields $Q_* = \frac{\gamma Q_m + \mu q_3}{\gamma + \mu}$. Solving $\frac{dx_3}{dt}(E) = 0$ for $x_3^*$ yields $x_3^* = 0$ or $x_3^* = \frac{\mu_3(Q_*) - D_3(Q_*)}{\delta_{33}}$.

Thus $E_1 = (0, 0, \frac{\mu_3(Q_*) - D_3(Q_*)}{\delta_{33}}, Q_*)$ is the only nonnegative boundary steady state for the system with zero as the first or second component. It is also the only steady state where only $x_3$ survives.

The next proposition presents a set of intuitive but unrealistic sufficient conditions for the extinction of some or all cancer populations.

Proposition 3. Define $T_1(Q) = \mu_1(Q) - D_1(Q) - \lambda_{13}(Q)$, $T_2(Q) = \mu_2(Q) - D_2(Q) - \lambda_{23}(Q)$, $S_3(Q) = \mu_3(Q) - D_3(Q)$. If $T_i(Q_m) < 0, i = 1, 2$, then $x_1, x_2$ eventually go extinct. If in addition $S_3(Q_m) < 0$, then $x_3$ also goes extinct.

Proof. Define

Then the rate of change of each population can be expressed as

By our definition of $T_1(Q)$ and $T_2(Q)$, a derivative test reveal that they are both strictly increasing with respect to the variable $Q$. Furthermore,

Since $Q(t)$ is bounded above by $Q_m$, in the case that $T_1(Q_m) < 0$ and $T_2(Q_m) < 0$, we see that

This leads to the estimate

This implies the extinction of the AD and weak-AI populations.

Additionally if $S_3(Q_m) < 0$, then asymptotically $x_3'(t) \leq S_3(Q_m)x_3(t)$, so the equation for strong-AI population is bounded above by $x_3(0) e^{-|S_3(Q_m)|t}$, which approaches to $0$.

For stability analysis, the following proposition extends Proposition 2 by establishing sufficient conditions of asymptotic stability of the equilibrium state with only the strong AI population.

Proposition 4. The steady state $E_1 = (0, 0, \frac{\mu_3(Q_*) - D_3(Q_*)}{\delta_{33}}, Q_*)$ is locally asymptotically stable if $S_3(Q_*) > 0$ and

Furthermore, a set of sufficient conditions for the locally asymptotically stability of $E_1$ is

Proof. The Jacobian matrix evaluated at $E_1$ is given by

where

and $M_{12}$ is the zero matrix. This format of the $J(E_1)$ immediately implies its eigenvalues are $-S_3(Q_*)$ and $-(\gamma+\mu)$ and the eigenvalues of $M_{11}$. Denote $M_{11} = \begin{pmatrix} a&b\\ c&d \end{pmatrix}$, then its eigenvalues are given by $\lambda_{1, 2} = \frac{1}{2} \left(a+d \pm \sqrt{(a-d)^2 + 4bc} \right)$. $E_1$ is locally asymptotically stable if $\lambda_{1, 2}$ are negative. Hence, we require $-(a+d) > \sqrt{(a-d)^2 +4bc}$. This holds if $a + d < 0$ and $a d > b c$, which are the proposed conditions.

4.   Numerical result

In this section, we present key numerical results comparing the two and three population models against clinical data. Additionally, we provide a sensitivity analysis to further enlighten the important parameters of the three population model. The sensitivity analysis for the two population model can be found in ^[2].

4.1.   Numerical comparison

We calculate the error of the fitting and forecasting parts for a quantitative comparison, where we used the predetermined on-treatment time given in the data to make the forecast. Table 2 shows comparable errors in the fitting of the first 1.5 cycles of both models. Additional criterion for comparison, such as the Akaike Information Criterion (AIC) that was used in similar contexts ^[14], could be used to strengthen the comparison. It is worth pointing out that information theoretic comparison like the AIC does not directly take into account of how well a model describes the dynamics of the natural phenomenon, for instance how well a model can follow the trend of the data. On the other hand, Table 3 shows the forecasting ability of both models are in comparable range. Instances where the models give large forecasts errors are likely due to the uncertainty associated with the number of parameters. Examples of the fitting and forecasting are presented in Figure 2. In addition, we also studied the case where $\lambda_{ij}(Q)$ varies between different transition pathways and the cases where either the mutation from $x_1$ to $x_3$ or from $x_2$ to $x_3$ is negligible. While each case shows good fitting, the forecasting is not improved likely due to the increase in uncertainty associated with the higher number of parameters.

4.2.   Fitting and forecasting of androgen and cell population

We present selected figures of androgen levels and cell populations corresponding to the forecasting in Figure 2.

4.3.   Sensitivity analysis

Sensitivity serves an important part in the analysis of any mechanistic model. Sometimes, by performing sensitivity analysis, we can pinpoint parameters that do not serve important parts in the dynamics of the models. While for nonlinear models, global sensitivity analysis techniques are often preferred, local sensitivity analysis serves as a good first step to disclose the relationship of small variation in parameter values and the variable of interest. Since sensitivity analysis is often affected by the magnitude of the parameter or the variable, we will use the normalized sensitivity as presented in section 5.5 in ^[23], which takes the form:

Here, $p$ is the parameter and $x$ is the variable of interest, for example PSA level. To carry out this process, we vary each parameter by 1% separately while fixing all other parameters. The normalized sensitivity is evaluated at 120 days. The result is shown in Figure 5. A positive peak shows a positive correlation between the parameter and the variable of interest, whereas a negative peak shows a negative correlation.

Figure 5 shows a dominating influence of the maximum androgen $Q_m$ on all variables. The cell populations are sensitive to the growth rate $\mu$ and the respective death rate half-saturation level $R$, while small variation in initial population has negligible effect on the overall dynamics. Other intrinsic parameters to the cell population such as the minimum cell quota $q$ and maximum cell death rate $d$ have minimal effects. The quantity of most interest, PSA level, is most sensitive to $Q_m$ and the clearance rate $\epsilon$. Curiously, $R_1$ also appears to have a considerable effect on the PSA level.

5.   Discussion

There is a need for mathematical modeling effort for understanding and treating (metastasized) prostate cancer. The aim of this work is to compare two different cell population structures and explore key aspects that could help future prostate cancer modeling effort. We summarize key results and discuss our findings here:

(A1) The data fitting part for both models are similarly accurate, especially when there is sufficient and clear regularity in the data as can be seen in the fitting and forecasting of patient 28 data in Figure 2. On the other hand, the two population presents better forecasting ability, perhaps due to having a smaller number of parameters. Since the incorporation of an additional cell population does not significantly improve fitting and worsen forecasting, it suggests that to model clinical PSA dynamics in the absence of extensive data, the cancer cells can be assumed to fall into two types, androgen dependent and androgen independent.

(A2) A delay of PSA relapse from the onset of the off-treatment period can be seen for many patients. The sudden switch in the androgen available between the off and on treatment periods does not account for this, leading to situations where the forecast could have been stronger had such measures have been taken, for instance patients 77 and 100 in Figure 2. The delay, which can be up to several months if present, can be seen as a sign of drug residual in the patient's body interfering with the production of androgen, or perhaps a more complex mechanism that delays the availability of androgen in the tumor. While incorporating a delay term to account for the observed delay is straightforward, a more biologically motivated alternative would be to use the interpolation of androgen introduced in ref. ^[22]; however, the latter approach may require additional assumption on the shape of androgen data for forecasting.

(A3) It is possible to obtain highly accurate fitting and forecasting with a small number of PSA data points (around 10), example can be found in simulation of patient 100 in Figure 2. Aside from the obvious assistance of concurrent fitting of androgen data, having a balanced ratio of high and low data points and regularity in patient data may play the key role. Although using a small number of data points may lead to over fitting, most parameters are relatively not sensitive. Thus a thorough sensitivity analysis follow by fixing the non-sensitive parameters is a possible method to deal with the problem of problem of over fitting. On the other hand, the disproportional number of data points will possibly effect the fitting and forecasting (see the three population model fitting and forecasting for patient 29). This suggests that a resampling technique should be applied and recent approaches in this direction can be found in ^[16].

Preliminary dynamical and sensitivity analysis of the model are also provided, which shows that the model structure is biologically sound. While some conditions for local stability and existence of a coexistence of all cancer cells are established, a global study is still needed to explore this further. The simulations in Figure 3 and Figure 4 suggest the global stability for the on and off-treatment connected by a switch. According to the simulation results, periodic switching of treatment may allow the model to generate a stable limit cycle solution. Furthermore, although a local technique has been used to explore the relationship between variables and parameters, global techniques is highly desirable for the sensitivity analysis due to the nonlinearity in the model ^[23]. In addition to data, it will be helpful to perform parameter identification analysis to exam the credibility of the parameter values identified via model data fitting prcess ^[8].

From the perspective of physicians and patients, the following questions could be useful for planning the course of treatment to optimize the effect. What are the ideal thresholds for PSA level in IAS? Under optimal treatment schedule and dosage, can IAS prolong a patient's life in comparison to CAS? What is a good metric to compare different treatments? When will the tumor become resistant to hormonal therapy? What triggers the development of hormonal refractory and how to prevent it? Each question by itself poses difficult practical and theoretical problems, for example the information on the types of cancer cells and how they interact is difficult to obtain. But perhaps, the most serious problem is the lack of a good practical method to track the development of a metastasized prostate tumor. In practice, PSA is primarily used as a direct correspondence to the total tumor population. This has led to the almost exclusive use of PSA-data-fitting as the only justification for mathematical models ^[20].

Acknowledgments

This work is partially supported by NSF grants DMS-1148771, 1518529 and 1615879. We are profoundly grateful to Dr. Nicholas Bruchovsky for the clinical data set and his encouragement. In addition, the authors would like to thank the reviewers for many helpful comments that improved the presentation of this manuscript.

Conflict of interest

All authors declare no conflicts of interest in this paper.