The Impact of Intermittent Androgen Suppression Therapy in Prostate Cancer Modeling

Previous studies on prostate cancer modeling under hormonal therapy successfully fit clinical serum androgen data, under the assumption that the levels of intracellular and serum androgen are similar. However, such an assumption may not hold throughout the course of treatment. In this paper, we propose a model that directly accounts for serum androgen and its interaction with intracellular androgen. We establish biological links between the model and clinical data, and discuss in detail parameter ranges and the initialization of model variables. We further investigate parameter sensitivity over time, which gauges the maximum effect of varying each parameter and allows us to fix some parameters, to increase the robustness of the parameter fitting process. By relying on the characteristics of intermittent androgen suppression therapy (IAS), we employ a two-part weighted error function for fitting. We also carry out mathematical analyses to study the dynamic aspects of the system with different androgen thresholds. We find that the proposed model shows superior forecasting ability, compared to its predecessor. Furthermore, we demonstrate the impact of androgen on the dynamics of the androgen-dependent and -independent cancer cells, which suggests the discrete description of androgen dependency may not give a realistic characterization of the cancer population. We show that IAS has certain characteristics that need to be considered for parameter estimation. Our results demonstrate that the model and the fitting scheme are viable for similar applications of prostate cancer modeling under hormonal therapy.


Introduction
Prostate cancer (PCa) is the second most common cancer found in American men.The estimated lifetime chance of being diagnosed with prostate cancer is 1 in 7.While prostate cancer is commonly found in men over 50 years old, it can start to develop in men in their twenties and takes decades to progress.Survival rates for localized prostate cancer have been improved significantly, due to advances in treatment; however, once cancer metastasizes beyond the prostate, the survival rate drops to about 28% [1,2].It has been noted that the growth of prostate cancer is highly dependent on androgen, a male characteristic hormone like testosterone.Thus, for metastasized cancer, the standard hormonal treatment is androgen deprivation therapy (ADT), which dates back to Huggins' Nobel Prize-winning discovery that castration causes rapid regression of PCa.While the treatment is successful in making the tumor regress initially, the cancer cells eventually mutate and adapt to be able to grow in an environment with low-androgen level, leading to the ineffectiveness of hormonal therapy or hormonal resistance [3].Additionally, ADT comes with many side effects, for instance, impotence, depression, bone demineralization, and an increased risk of dementia.Thus, intermittent androgen suppression therapy (IAS) was introduced.During IAS, the patients are put on treatment until their prostate-specific antigen (PSA-a prostate cancer biomarker) level reaches a desired low threshold, then they are taken off the treatment until the PSA level rises above a predetermined threshold.The process is then repeated until the patients reach hormonal resistance [1].
Various mathematical models have been developed, in the past two decades, to study many aspects of prostate cancer, especially under androgen suppression therapy.The first experimentallydriven model of PCa under continuous androgen suppression therapy (CAS) is credited to Jackson in 2004 [4].This inspired a series of works that study PCa under IAS and compare the effect of CAS versus IAS [5][6][7][8].Furthermore, theoretical perspectives on population dynamics and spatial structure have also been explored in subsequent works [9][10][11].In addition, a biologically detailed model with a complex biochemical network has also been developed and verified with experimental data [12].On a more applicable note, optimal scheduling, complex algorithms to predict hormonal resistance, and stochastic models have been used to provide insights into clinical aspects of prostate cancer treatment [13][14][15][16].For a detailed review of recent work, the authors refer to [17].
In this paper, we focus on a previous approach by Baez and Kuang, which was the first successful attempt at using both PSA and androgen data for fitting and forecasting [8].A recent work by Phan et al. [18] shows that, with the available types of data (e.g., PSA and androgen data), the model in Baez and Kuang (BK model) has the optimal population structure to describe and to forecast the dynamics of prostate cancer under IAS.On the other hand, the model fails to completely capture the growth phase of the PSA level when the patient is taken off treatment.Additionally, fitting of the androgen data is done by assuming similar levels of body androgen and intracellular androgen, which may not hold throughout the course of treatment.Thus, we extend the BK model to directly take into account serum androgen separately from intracellular androgen.We carry out numerical comparisons between the two models, in terms of fitting and forecasting.Finally, we study the mathematical properties of the new model, to show the effect that androgen has on the survival of cancer cells.
The data is taken from a clinical study from the Vancouver Prostate Center [19].Similar to [18], we select 26 patients with sufficient data for fitting and forecasting comparisons, over 2.5 cycles of treatment.

The BK Model
Baez and Kuang [8] introduced a two-subpopulation model for prostate cancer undergoing androgen suppression therapy.This model is a simplification from the model in Portz et al. [7].The model takes the form: tumor's PSA production In this model, x 1 and x 2 represent the androgen-dependent (AD) and androgen-independent (AI) cells, respectively.P is the PSA level, and Q is the intracellular androgen level.In [8], the authors assumed that the level of intracellular androgen is the same as that of the serum androgen.Since only serum androgen data is available (and not intracellular androgen data), this assumption allows data fitting using clinical androgen data with the variable Q, along with the data fitting of tumor biomarker PSA.The model assumes q 1 > q 2 , which biologically means that the AI cell subpopulation can grow in an environment with lower androgen levels than AD cells can.A more thorough derivation of this model can be found in [8,18].As pointed out in [18], the mutation from AD to AI cells is sufficient to describe the transformation between these two cells, due to the lack of clinical data types.
As mentioned in [18], the equation governing the dynamics of Q has a significant biological flaw.Since the diffusion is scaled with the difference between Q m and the current value of Q, it means that immediately after the patient is taken off treatment, the serum androgen level reaches its max (e.g., Q m ).Scaling the growth of intracellular androgen in this way is perhaps one of the causes for the sudden jump in PSA level between the on-and off-treatment intervals.

The New Model
In the BK model, the assumption of similarity in the levels of serum androgen and intracellular androgen may not hold consistently throughout the treatment.This motivates the following model, which directly takes into account body androgen production and its diffusion to intracellular androgen.Additionally, previous work by Portz et al. [7] extrapolated serum androgen data with an exponential growth/decay function.The serum androgen values are then fed into a five-compartmental model, known as the PKN model.Extending upon this idea, we explicitly include a compartment that governs the dynamics of serum androgen.The following model is presented as an extension of the BK model: suppression of production (10) tumor's PSA production In contrast to the BK model, the variable A in this model represents serum androgen.While the treatment is on, the rate of change of A is strongly influenced by a negative exponential term, indicating the effect of the total androgen blockage drug.While the treatment is off, the production of A takes the form of a logistic-type growth, where A 0 is the maximum serum androgen, which is akin to the carrying capacity in the logistic equation.The uptake rate of intracellular androgen, Q, depends on the current level of serum androgen.This is a biological improvement over the BK model as, after the patient is taken off treatment, it takes time for the serum androgen to recover to its maximum.As a consequence, the intracellular androgen level will take more time to recover its normal level.

Sensitivity Analysis for the New Model
Sensitivity analysis of the parameters plays an important part in evaluating and improving the robustness of a model.Previous studies focus on the local sensitivity at a fixed time point [8,18].While this can provide valuable information, the dynamics of IAS therapy are broken down into onand off-treatment intervals, with distinct characteristics in both parts.Additionally, understanding the effect each parameter has on the dynamics of the model can lead to the fixing of certain parameters, during the parameter estimation process.If some of the less sensitive parameters can be fixed, then we essentially reduce the uncertainty in the remaining parameters, which increases their identifiability.Here, we employ the normalized sensitivity equation from section 5.5 in [20], to carry out the sensitivity analysis of the new model: where p is the parameter, and x is the variable of interest; for example, µ and PSA level, respectively.Note that, since we will be varying each parameter by 1% with respect to their value, we can use the following approximation of the sensitivity equation: where x varied is the value of x after increasing a parameter by 1%.The result of the sensitivity analysis is summarized in Table 1.
While the sensitivity process is perhaps more relevant when carried out on patient data, there is a large variation in the data between patients, and so finding a consensus among all patients can be difficult.Instead, we: (1) Synthesize a set of parameters within the range in Table 2, and a treatment schedule that produces the expected dynamics observed in clinical data.(2) Using these parameters, we carry out the sensitivity process by increasing one parameter at a time by 1% of its value.(3) We approximate S p at each time point for 1 treatment cycle using Equation (13).The process is repeated for all parameters.Finally, we record the maximum order of sensitivity in Table 1.
Table 1 shows the maximum effect that changing a parameter by 1% has on the variables.For instance, by looking at the intersection between c and A, we see that changing c by 1% can cause A to change by a maximum order of −8.To compare the sensitivity of these parameters, we add up their sensitivity coefficients.Doing so, we find that the parameters with combined order less than −10 (least sensitive) are c, K, δ 1 , δ 2 , and d 2 .Since d 1 and d 2 serve similar purposes and d 1 is a sensitive parameter, we choose to not to fix d 2 .In summary, for the fitting process of the new model, we fix c, K, δ 1 , δ 2 , and γ 2 .Note that we do not carry out the sensitivity analysis for γ 2 , due to its small value; in fact, we fix γ 2 to 0 for the new model (although a small value would suffice).The values that we fix for c, K, δ 1 , and δ 2 are 0.00015, 1, 5, and 5, respectively.Table 1.Order of sensitivity.The number represents the maximum order of sensitivity of a variable with respect to a parameter.For instance, at the intersection of µ 1 and x 1 , 0 means changing µ 1 by 1% will affect at most a change of order 0 in x 1 .In this table, 0 corresponds to most sensitive and −8 means least sensitive.The most insensitive parameters with combined order less than −10 are c, K, δ 1 , and δ 2 .Thus, we fix them in the new model.Note that d 2 also has a combined order of sensitivity less than −10; however, since d 1 is sensitive, we do not fix d 2 .Also, we do not carry out the sensitivity for γ 2 , since we fix it as 0 for the new model, due to its small range.

Parameter Range
Previous work employed various parameter ranges without a detailed discussion of how the parameter range is obtained.This caused an unnecessary inconsistency across the literature.Here, we present our method of approximation of the parameter range used for fitting in Table 2.This is possible, partially because our model has clear biological links to clinical data and previous work.In doing so, we hope to build upon this, in future work, for coherency across models.First, we note that the units of x i , Q, A, and P are respectively L, nmol/L, nmol/L, and µg/L.From this, the following ranges and unit of parameters are established: Berges et al. [21] studied the proliferation and death rates of prostate cancer cells.We interpret their result for metastatic prostate cancer cells from hormonally untreated patients to be suitable for AD cells because, before hormonal treatment, most cancer cells should be androgen-dependent.This gives us a range of (0.004-0.081) for µ 1 .Similarly, for µ 2 (the maximum proliferation rates for AI cancer subpopulation), its range (0.001-0.046) is derived from the result for hormonally failing patients.This is because we expect AI cells to take over, once the patients become hormonal refractory.Taking into account experimental and survey errors (especially in this case, where the cancer population size is small), we use (0.001-0.09) as the range for both µ 1 and µ 2 in the parameter fitting process.• d 1 is the maximum death rate for AD cells.Its unit is [day] −1 .Its range is also derived from Berges et al. [21], under similar assumptions in the case of the proliferation rates.The ranges for d 1 and d 2 are (0.001, 0.0525) and (0.015, 0.0775), respectively.For parameter fitting, we use the ranges (0.001, 0.09) and (0.01, 0.09) for d 1 and d 2 , respectively.• q 1 is the minimum androgen cell quota for AD cells.It has the same unit as Q, [nmol/L].Tsutomu Nishiyama [22] showed in his review paper that hormonal suppression therapy aims to reduce serum androgen to below 0.69 nmol/L, which is considered the castration level (this threshold varies in practice, but no more than 1.73 nmol/L).The lowest threshold recorded in Nishiyama [22] is 0.41 nmol/L.This gives us an approximate bound for q 1 , (0.41-1.73) because, once serum androgen drops below this range, the growth of AD cells is repressed.q 2 , the minimum androgen cell quota for AI cells, should be smaller than 0.41 nmol/L, to allow AI cells some growth during castration.Since no information is available on the lower bound of q 2 , we approximate the range to be (0.01-0.41).Since the value of q i is patient-specific, a fixed range will not suffice for every patient.Thus, we find the minimum serum androgen (min A ) in the first 1.5 cycles (from the data) and use the ranges (min A + 0.2, min A + 0.5) and (0.01, min A + 0.2), respectively, for q 1 and q 2 to reflect the condition q 1 > q 2 and to represent the effect of cell quota on growth.

•
The δ i 's are the density death rates.Their unit is No experimental information is available for the density death rate of cancer cells.Additionally, our sensitivity results in Table 1 indicate that their importance is negligible, thus we fix them at 10 −5 [mm 3 ] −1 [day] −1 , or 10 [L] −1 [day] −1 using the value taken from Baez and Kuang [8].We keep δ i in our model to avoid unrealistic overgrowth of the cancer cell population.• K is the half-saturation level of the androgen-dependent mutation rate from AD cells to AI cells.Its unit is [nmol][L] −1 .Baez and Kuang used the range (0-1), while Portz et al. [7] used (0.08-1.7).Since no experimental validation is available, and as K is one of the least sensitive parameters, we fix K to be 1 [nmol][L] −1 .• R i are the half-saturation level of androgen-dependent death rates for AD and AI cells.Its unit is nmol/L.While no experimental information is available, computation results in [8,23] show that (0-3) and (1-6) are reasonable ranges for R 1 and R 2 , respectively.Thus, we employ these ranges.• c is the maximum mutation rate from AD cells to AI cells.Its unit is [day] −1 .While there is no experimental ranges for c, numerical experimentation using a mathematical prostate model from Ideta et al. [5] established a reasonable range for c to be (10 −5 -10 −4 ).This is often fixed at about 0.00015 in the literature [7,23].Since c is one of the most insensitive parameters, we choose to fix it at 0.00015 Baez and Portz [7,8] used two different ranges for b.Taking the lower and upper limits of both, we will use (0.0001-0.1) as the range for b.We note that, intuitively, the baseline production of PSA is produced by normal prostate cells; however, to simplify the model, we assume that this production is proportional to the intracellular androgen level, because the lack or abundance of androgen directly affects the product from prostate cells.• σ is the androgen-dependent PSA production rate of the cancer cells.[23] used several ranges within (0.001-1), so we employ (0.001-1) as the range for σ.

•
is the PSA clearance rate.Its unit is [day] −1 .While no clinical experiment supports a range for , Portz et al. [7] used 0.08, and Baez and Kuang [8] extended this to (0.001-0.01).We further extend the range to (0.0001-0.1) for our parameter estimation.
• γ 1 is the androgen production rate associated with the testes, the primary source of androgen production.Its unit is [day] −1 .Previously, the value used in Baez et al. [8] seemed large, which resulted in an instantaneous increase in androgen level after the on-treatment interval stopped [18].On the contrary, Ideta et al. [5] set γ 1 to be 0.08.Thus, for our purpose, we use (0.008-0.8) as a possible range for γ 1 .• γ 2 is the secondary androgen production rate, produced from sources other than the testes, such as the adrenal gland.The range established from [8] is (0.001-0.1) for γ 2 , where their γ 1 = 20.However, the value of γ 2 should only be a small fraction of γ 1 .Thus, in our case, the value of γ 2 is negligible.Hence, we fix it at 0 for the fitting process.• m is the diffusion rate that links serum androgen level with intracellular androgen level.It is unitless.This rate is taken at 1 implicitly in [7,8].Here, We allow m to vary from 0.01 to 0.9, to incorporate the level of synergy between A and Q.

•
A 0 is the maximum serum androgen level.Clinical data shows that typical serum androgen range is below 27 [nmol][L] −1 [19].Furthermore, computational work often assumed [5,8].Since this range is patient-specific, for parameter fitting we find the maximum serum androgen (max A ) in the first 1.5 cycles (from data) and use the range (max A − 5, max A + 5) for A 0 .
Aside from the ranges for parameters, we also discuss the initial conditions for the model.A(0) and P(0) are taken as the first recorded serum androgen and PSA level (pre-treatment), respectively, from patient data.We assume Q(0) to be 90% of A(0).This assumption means, prior to treatment, the level of serum androgen and intracellular androgen are approximately the same.For x 1 (0) and x 2 (0), we use experimental data on the effect of androgen suppression therapy on prostate volume (not tumor) from Bruchovsky et al. [19].To derive this value, we compare the prostate volume pre-treatment (mean 24.7 cm 3 ) and after the first treatment (mean 14.7 cm 3 ).We assume that the AD cells are mostly responsible for this initial reduction, or an approximate initial size (for metastasized tumor) of 10 cm 3 , or 0.01 L. This indicates a biologically reasonable range for the initial volume of x 1 .Since other factors can affect this value, we use a range of (0.009-0.02) as our possible range for x 1 (0).Additionally, if we assume that the initial subpopulation of x 2 (0) is about 1% of x 1 (0), then we obtain the range (0.00001-0.001)for x 2 (0).The parameter ranges and initialization are summarized in Table 2.
Table 2. Parameter definition and range, established via clinical data and previous computational work.Note that the parameters with an asterisk next to their fitting range are fixed in the new model, e.g., c, K, δ 1 , δ 2 , and γ 2 are fixed at 0.00015, 1, 5, 5, and 0, respectively.The ranges of parameters with a double asterisk are dependent on the data from the first 1.5 cycles for each patient.

Parameter Estimations
For parameter fitting , we use function fmincon from MATLAB (version 2018b), which utilizes an interior point algorithm to minimize an objective error function within the range provided in Table 2. Previously, the objective error function was calculated as the sum of the mean squared error (MSE) fitting error for PSA and androgen, because MSE punishes particularly large errors (which is undesirable in clinical applications).This emphasis on good fit can be problematic for forecasting, as the peak PSA level (during the off-treatment interval) is much larger than the PSA level during the on-treatment interval.Specifically, MSE will prioritize large errors, thus fmincon will estimate parameter values to achieve high PSA level to match the data at peak level, which leads to overestimated forecasts.This problem is escalated, as the time for on-and off-treatment intervals are taken from data.Furthermore, the problem is not resolved using other error measures.For instance, using the relative error metric causes the fitting to skew toward smaller data points due to their number, and so forecasting is often significantly underestimated [18].Instead, taking into account the characteristic of the treatment schedule and its effect on PSA level, we employ a weighted least square error that will allow emphasis on the recent data, and a controlled weight on the error contribution from the peak level.First, we introduce the weighted error function: where N is the total number of data points, y is either the PSA or androgen data generated from the model, and ŷ is the respective observed data.t f is the final time of the second treatment period, and t i is the time of data point ŷi .Here α is the exponential weight, which is usually taken to be between 0 and 1.The larger α is, the more memory-less the process becomes, e.g., earlier data points become weightless.Since we want to control the error contribution from before and after the peak PSA level, and as the data are collected somewhat irregularly, we employ the following scheme to give sufficient weights to the data before and after the peak: Here, W 1 and W 2 are the desired contributions of error from prior and after the peak PSA level, respectively.Similarly, N 1 and N 2 are the total number of data points, for the respective part.We find that setting α 1 = 0.01 = α 2 and W 1 = W 2 = 0.5, meaning equal contribution of error from both parts, provides consistent results.Note that, in practice, the values of α i and W i should be specialized to each patient.Here, we only wish to show their effectiveness in fitting and forecasting.The objective error function that is minimized by fmincon is the sum of PSA and androgen residuals, which takes the form: Objective MSE = Error MSE,PSA + Error MSE,androgen .

Numerical Results
For comparison purposes, we calculate the errors from the fitting and forecasting of BK and the new model.Table 3 shows a comparable good fit in the first one and a half cycles, for both models.On the other hand, Table 4 shows that, while both models produce good forecasting, the new model outperforms the BK model in term of forecasting.Note that, for the BK model, we fit 17 parameters and 2 initial conditions, while the new model fitting accounts for 14 parameters and 2 initial conditions.Thus, information criteria (like the AIC) can be used to further the comparison between the two models, as done in [24].However, such comparison does not evaluate the biological mechanism of the model, thus we omit it here.We select representatives of fitting and forecasting results, and present it in Figures 1 and 2.
Table 3. Fitting mean squared error (MSE) comparison (median, 1st quartile, 3rd quartile).Two models are compared, using 26 patient's data.Both models show comparable good fitting results.Note that the fitting ability of model 2 has been shown to be superior, among existing work [8,18].Additionally, the incorporation of serum androgen in the new model does not significantly improve the fitting of androgen.

PSA Androgen
Median  [8,18].Additionally, the forecasting of androgen is comparable between the two models.

PSA Androgen
Median  The blue line divides the fitting and forecasting period.Both models are fitted using data from the first 1.5 cycles of treatment (left side of blue line).Then, we forecast out one additional cycle.The data are represented by the black circles.Note that, while we do not experience the sudden jump in PSA level mentioned in [18], the forecast still fails to capture a delay in PSA relapse immediately after entering the off-treatment interval, which occurs in some patients; for example, patient 77.

Fitting and Forecasting of Androgen and Cancer Cell Dynamics
We provide an example of fitting and forecasting comparisons for androgen on the left figure in Figure 2. Note that while both models fit and forecast the dynamics of androgen comparatively well, the new model exhibits a higher flexibility, in contrast to the square-like shape of androgen in the BK model.dWe also provide an example of cancer cell dynamics, in the right figure in Figure 2. The result from cancer dynamics suggests that the onset of hormonal resistance can develop early into the treatment.Additionally, since cancer subpopulations can go into near extinction (during the on-treatment interval) before bouncing back, this suggests that the mechanism used to incorporate androgen into growth and death rates needs improvement, or perhaps the problem asks for an optimization of parameter values using analytical means.Finally, since the cancer population can stay under control after AI cells overtake AD cells in number, it suggests against the intuitive definition that hormonal resistance occurs when AI cells overtake AD cells.

Note on Numerical Fitting Using Two-Part Weighted Error
Different patients react differently to the treatment, especially during the first cycle of treatment-the time for PSA to fully relapse varies between patients.Figure 3 shows an example where PSA level takes over two years to reach the upper PSA threshold.Recall that our fitting scheme utilizes a fixed exponential weight that works well for most patients, but fails in this case due to the abnormally long time it takes PSA to relapse (left figure).This indicates that our choices of α i = 0.01 and W i = 0.5 are not suitable for this particular case.As a counter-example, by reducing the exponential weight and readjusting the contribution of errors from each section of data (e.g., adjusting α 1 to 0.001, W 1 = 0.7, and W 2 = 0.3), the fitting and forecasting results significantly improve (right figure).

Dynamics of the New Model
In this section, we carry out mathematical analyses for the newly proposed model, which help to justify its coherence with biological expectations.Previous studies [8,18] have already examined some standard properties of the BK model and its variations, so we only focus on the new model here.We remark that, while PSA serves as an important indicator of how good the model is against clinical data, the dynamics of the system are independent of the dynamics of PSA.In addition, if the dynamics of x i , Q, and A are known, one can immediately infer the dynamics of P. Thus we will only carry out our analysis with x i , Q, and A. The first proposition establishes a positively invariant region, which is also bounded.This supports the biological validity of the model since we expect x i , Q, and A to stay non-negative and never become infinite.
and assume q 2 ≤ q 1 < A and γ 2 > γ 1 q 2 .Furthermore, we assume that when androgen is abundant, the growth rate is larger than the death rate, meaning µ i (1 , we see that A (q 2 ) ≥ 0 and A (A) ≤ 0. It follows that q 2 ≤ A(t) ≤ A for all t > 0 with initial condition in the same region.Similarly, since , we find that Q (q 2 ) ≥ 0 and Q (A) ≤ 0. Thus, a similar conclusion follows.Consider the rate of change of x 1 , This implies that the limit supremum of x 1 (t) is bounded by 1 In the limit t → ∞, x 2 is bounded above by The function h(x 2 ) is a concave-down parabola.Thus lim sup t→∞ x 2 (t) is bounded above by the larger solution of h(x 2 ) = 0. Solving for h(x) = 0, we obtain the upper bound of the limit supremum of x 2 , Since x 1 (0) < M x and x 2 (0) < M y , x 1 and x 2 are bounded above by M x and M y , respectively.Additionally, if x i (0) > 0, then the positivity of x i (t) follows immediately.
The second proposition shows that the system has a boundary equilibrium where AD cells become extinct, during the off-treatment interval.Additionally, it establishes sufficient conditions for the existence of a unique positive steady-state.Proposition 2. In the case of u = 1, assume q 1 > q 2 , then where Q * = mA * +µ 1 q 1 m+µ 1 and A * = (γ 2 + γ 1 A 0 )/γ 1 , is the only equilibrium state where the AD subpopulation becomes extinct.Furthermore, suppose conditions (1) and (2) hold: where Then the system has a unique positive steady-state (x * 1 , x * 2 , Q * , A * ).
Proof.For u = 1 (off-treatment interval), we first consider the boundary equilibria of the system.Since (x 1 , x 2 ) = (0, 0) is undefined for the system and x 2 * = 0 implies x 1 * = 0, the only equilibrium involving an extinction of a cancer subpopulation takes place when x 1 * = 0.In this case, solving . Additionally, solving Q = 0 and A = 0 yields For the positive steady state, suppose that x * 1 , x * 2 , Q * , A * > 0. Then solving for x 1 = 0 and x 2 = 0 gives: where A i are the coefficients of x * i , e.g., Dividing the previous two equations to obtain then by Descartes's rule of sign, the system has a unique solution y, regardless of the sign for A 2 .The uniqueness of y implies that the ratio We proceed to establish the positivity of A 1 .Explicitly, A 1 takes the form Note that, since the treatment is off, biologically Q(t) >> q 1 , which allows strong growth of the AD cell subpopulation.We will establish later that there is indeed a lower bound of Q(t) satisfying lim inf Q(t) > q 1 , under reasonable biological constraints.For the moment, we assume such lower bound exists.Then where Q m >> q 1 is a lower bound of Q(t) for the off-treatment interval.We would like for A 1 > 0, thus rearranging terms give the condition: Next we examine the expression for Q (t) to find an appropriate lower bound.Note that, since as t → ∞, A(t) → γ 2 γ 1 + A 0 =: A * , we have the following: Then as t → ∞, the exponential term approaches 0, so This means that, eventually, Q(t) will be bounded below by which is true because of our biologically reasonable condition mA 0 m+max{µ 1 ,µ 2 } > q 1 (e.g., A 0 >> q 1 ).Thus, we establish that indeed there is a lower bound of Q(t) that satisfies lim inf Q(t) > q 1 .Now, the condition in Equation ( 17) is satisfied in the limit of t → ∞, if mA 0 +min{µ 1 q 1 ,µ 2 q 2 } m+max{µ 1 ,µ 2 } is greater than the larger solution for Q m (or smaller than the smaller solution) in Equation 17.Therefore, A 1 > 0 requires either or, where µ 1 max{R 1 , K}.Since the condition in Equation ( 19) is not possible (e.g., the right hand side is negative, but the left hand side is positive), the only possible condition is Equation (18), which is the condition in the proposition.This establishes the existence of positive x * 1 and x * 2 , which have a fixed ratio y * .Solving Q = 0 for Q * gives the positive steady state Finally, to show uniqueness, recall that the ratio of x * 1 and x * 2 is a fixed constant, y * .Thus, substitute x * 2 = x * 1 /y * in Equation ( 20) to obtain . This shows that Q * is unique.Finally, from Equation ( 14), . It follows that x * 1 and x * 2 are also unique.Since lim t→∞ A(t) = γ 2 γ 1 + A 0 =: A * , the existence of a unique positive equilibrium follows.
The third proposition shows that when the patient is on treatment, under biologically reasonable conditions, the boundary equilibrium with the extinction of AD cells is the only equilibrium of the system.In this case, it establishes a sufficient and necessary condition for the boundary equilibrium to be globally stable.Additionally, it provides sufficient conditions for both cancer subpopulations to eventually become extinct.Proposition 3. In the case of u = 0, if then the system does not have a positive steady-state and x 1 becomes extinct as t → ∞.Furthermore, Q(t) approaches Q * = mA * +µ 1 q 1 m+µ 1 , where A * = γ 2 γ 1 .In addition, if , Q * , A * ) is globally stable.On the other hand, if either then both x 1 and x 2 become extinct in the limit of t → ∞.
Proof.Similar to the previous proof, we establish that in the limit of t → ∞, A → γ 2 γ 1 and Thus there is no positive steady-state for x 1 (t) and lim t→∞ x 1 (t) = 0. Then in the limit of x 1 → 0 and A → γ 2 γ 1 , and E 1 is globally stable.

Discussion
Mathematical models for prostate cancer have made significant progress since the work of Jackson [25].Consequently, recent approaches have been focused on the application aspects in clinical settings.However, there are still obvious questions without an easy answer; for example, how one would define hormonal resistance in a rigorous mathematical sense and yet still meaningful in applications, or the identifiability of the model parameters.These difficult questions, that link theoretical work and application, ask for mechanistic models that can test theories against data.To this end, we aim to build upon a series of previous works to create a model with direct links to clinical data.We study the model and explore various aspects relating to using it.In this section, we summarize our findings and point out the limitations of our work.
(A1) From Table 3, both the BK model and the new model achieved similar accuracy for fitting.Since we used a two-part fitting scheme with exponential weight, the fitting was good when PSA quickly relapsed after the patient was put off treatment; see patients 39, 77, and 79 in Figure 1.On the other hand, when PSA took a longer time to reach the peak, the accuracy of fitting was not as good as in the previous cases; for instance, in patient 28 in Figure 1.In other words, because we chose to use a fixed exponential weight for fitting, the earlier the data points, the less weight that those data points contributed in the parameter estimation process.Thus, for patients with an abnormally long first cycle, most of the earlier data were essentially not accounted for during fitting; see patient 66 in Figure 3.To resolve this problem, a possible scheme is to use a linear weight.Specifically, we would first normalize the weight contribution from the peak PSA level to be 1.Then, the weight contribution of the points before the peak is the relative height with respect to the peak.This would ensure that earlier points have a more significant contribution to the fitting.While it is well-suggested that individual patients should have their unique parameter set, we further suggest that the process of studying the patient data should be unique, in order to capture the patient's characteristics; for example, a patient-specific weight for fitting in this case.
(A2) From Table 4, the new model slightly outperformed the BK model, in terms of prediction, which is noteworthy as the new model used fewer parameters for fitting.An interesting observation is that the new model did not overshoot forecasting, perhaps due to the parameter ranges established in Table 2.This can become a problem, if there is a stable equilibrium below the upper threshold PSA level required for putting the patient on treatment.In such a case, the model would be unrealistic for the parameter values which influence the equilibrium.While this problem can be mediated by forcing hormonal resistance within five years, as done in [6], it is perhaps better done by optimization of parameter ranges via mathematical analysis.We further note that when the data had good regularity, in terms of the similarity in the shape of data for different cycles (see patients 28, 77, and 79), the forecasting was quite accurate.However, when unforeseen changes occur in the data (patient 39), the model may fail to capture this change.Hence, the use of a stochastic model or uncertainty quantification is necessary for such applications.
(A3) Phan et al. [18] noted the problem of the sudden jump in PSA between the transition of the on-and off-treatment intervals in the fitting and forecasting of the BK model.Our result suggests that the main cause was the abnormally large value of the rate of primary androgen production (γ 1 ).
By adjusting this value in accordance to [6], we resolved this problem, as seen for patients 28, 39, 77, and 79 in Figure 1.Phan et al. [18] also mentioned a delay that sometimes occurs during the relapsing period of PSA; for example, in patient 77.While this can potentially be a problem with the parameter, it is more likely caused by the drug residue or the mechanism of how the patient's body reacts to the drug.Thus, we suggest future attempts take the drug dynamics into account.
(A4) We carried out a mathematical analysis of the model, to validate its biological properties and to show that it is possible to drive the cancer cells to extinction, as seen in the limit of time approaches infinity (in an environment of sufficiently low androgen level).
The right figure in Figure 2 shows that the cancer cell subpopulations can indeed reach near-extinction, prior to the on-treatment phase.We note that, while the cause of this is likely due to setting the secondary production term to 0 (e.g., γ 2 = 0) in our numerical study of the new model, it suggests that if the androgen level drops below the cell quota of AI cell, the AI subpopulation can become near-extinction for the duration of on-treatment.However, this is only possible in a theoretical sense.In practice, it is difficult to find clinical studies that support this theoretical result.This is because the androgen level required for the extinction of cancer cells may be much lower than the limit of hormonal therapy.Furthermore, the result will not hold if there is a variation of the cancer cells that gives some cells near-complete independence from androgen, which would make the condition in proposition 3 impossible to achieve.To resolve this in a theoretical sense, we could replace the two cell subpopulations with N subpopulations.Taking the limit as N → ∞ represents the continuous variation in cancer cells.However, such complexity resulting from adding N compartments of cancer cells is unlikely to be either insightful or practical.Instead, we suggest future work to incorporate a mechanism that takes into account the degree of androgen dependency of the tumor in a continuous spectrum, without adding unnecessary compartments.This may allow for a more useful characterization of the condition for cancer extinction under hormonal treatment.
In addition, we also studied the sensitivity of the parameters over time.While this method was one at a time, this helped to gauge the maximum effect that changing a parameter by a small percentage could have on the dynamics of the system.By doing so, we highlighted the most insensitive parameters and fixed them, which improved the robustness of the fitting and forecast.Since the model contains many parameters, the question of identifiability is of significance, especially in the context of precision medicine.An example of analysis involving parameter identifiability for cancer modeling can be found in [26].Furthermore, the uncertainty in the data and model means that uncertainty quantification is a necessary component in clinical setting [27].Finally, the parameter range, provided with the method used to obtain them, is not meant to be extensive or absolute; however, it serves as a unified framework for future improvement.In conclusion, future work should consider the identifiability and improve the accuracy of the model, in terms of both computational robustness and biological relevance.In the future, we will carry out comparisons of this new model across several different models by other researchers, to further advance the field.Additionally, mathematical analyses will be significant to help rigorously define concepts and propose hypotheses that can be tested to improve our understanding of prostate cancer and its treatment.If a model is identifiable (or if we know the relationship between the parameters), then we can use the fitted parameters and mathematical analysis to classify cancer types.Finally, uncertainty quantification for the prediction of the model is necessary to quantify how likely an additional cycle of treatment will be useful or not.opportunity.In addition, the authors would like to thank the reviewers for many helpful comments that improved the presentation of this manuscript.

Figure 1 .
Figure 1.Fitting and forecasting comparison of the BK model (magenta) and the new model (green).The blue line divides the fitting and forecasting period.Both models are fitted using data from the first 1.5 cycles of treatment (left side of blue line).Then, we forecast out one additional cycle.The data are represented by the black circles.Note that, while we do not experience the sudden jump in PSA level mentioned in[18], the forecast still fails to capture a delay in PSA relapse immediately after entering the off-treatment interval, which occurs in some patients; for example, patient 77.

Figure 2 .
Figure 2. (1) The left figure shows the fitting and forecasting of serum androgen for the BK model (magenta) and the new model (green).The blue line divides the fitting the forecasting periods.(2) The right figure shows the cancer cell dynamics, over the course of 2.5 cycles.The solid green curve represents the AD cell subpopulation, while the dashed green curve represents the AI cell subpopulation.

Figure 3 .
Figure 3.The figures show a case for a patient-specific fitting scheme.The figures represent the fitting and forecasting of the BK model (magenta) and the new model (green).This particular case sees PSA taking an abnormally long time to reach the upper PSA threshold for treatment.(1) Left figure: The weight used across all patients do not work in this case.(2) Right figure: Using a more appropriate weight produces significant improvement.

Table 4 .
Forecasting mean squared error (MSE) comparison (median, 1st quartile, 3rd quartile).Two models are compared, using 26 patient's data.The new model demonstrates a slightly better forecasting ability, compared to the BK model.Furthermore, the forecasting ability of model 2 has been shown to be superior, among existing work under different metrics