1. Introduction

Prostate cancer is the most frequently diagnosed cancer and second leading cause of death from cancer in men. Growth of this cancer is stimulated by androgens, or male sexual hormones. These androgens circulate in the blood and diffuse into the tissue where they stimulate the prostate tumor to grow. One treatment for advanced or metastatic prostate cancer is androgen deprivation therapy (ADT). Androgen deprivation may be achieved by hormone therapy, which inhibits the production of androgens in the testicles [11]. Leuprolide acetate and goserelin are drugs currently used for this treatment that can be delivered continuously over extended time periods using depot injections [14]. The influence of the remaining androgens produced by other sources, such as the adrenal glands, can be eliminated by anti-androgens such as flutamide, bicalutamide, enzalutamide, and nilutamide [11]. A combination of anti-androgens with chemical castration via ADT is known as the maximal androgen blockade (MAB). Both ADT and MAB facilitate apoptosis, the programmed death of androgen-dependent (AD) cancer cells, and quickly induce temporal regression of tumors.

Currently, many patients receive continuous androgen suppression (CAS) therapy of ADT and MAB. However, many of these patients undergo a relapse with an increase of the PSA level within three years after the initiation of ADT [9]. Androgen-independent (AI) cells are thought to be responsible for this recurrent tumor growth. These cells are not sensitive to androgen suppression, but rather they replicate even in an androgen depleted environment. The AI cells are generated by mutation from AD tumor cells. Once the tumor acquires these AI cells, androgen deprivation is unable to inhibit the cancer growth and a relapse becomes inevitable. In one study, 5 of the 61 patients demonstrated this progression within 4 years [13]. It is important to prevent a relapse or at least delay the time to relapse as long as possible. At the same time, it is also important to reduce economic costs and alleviate adverse side effects of prolonged androgen suppression such as osteoporosis, cardiovascular disease, anemia, and metabolic disorders [1].

A possible strategy to delay the progression from the AD state to the AI state is intermittent androgen suppression (IAS), which is a form of androgen ablative therapy delivered intermittently with off-treatment periods. On-treatment periods last for several months until PSA levels fall below a prescribed threshold and then, to avoid emergence of AI cells, the IAS therapy introduces off-treatment periods that serve to maintain the androgen-deprivation sensitivity of the cancer cells and restore their apoptotic potential, which can be induced by androgen deprivation [10]. Fourteen studies of nineteen models published have confirmed improvement in the quality of life during the off-treatment periods and alleviation of adverse side effects such as sexual dysfunction, hot flushes, and fatigue [1]. However, it is unknown how to optimally plan the IAS therapy.

Ideta, et al. [9], introduce a mathematical model commonly referred to as the ITTA model that describes the growth of a prostate tumor under IAS therapy based on monitoring of the serum prostate-antigen. By treating the tumor growth as a mixed assembly of androgen-dependent and androgen-independent cells, they investigate the difference between CAS and IAS in order to understand the factors that result in AI relapse. The ITTA model is known as a hybrid dynamical system because tumor growth is continuous in time whereas the on-and off-treatment protocol is discrete in time. In [9] and [18] bifurcation analysis of the ITTA model is used to characterize parameter regions that distinguish a prevention of relapse (cancer free state) from the occurrence of relapse.

In contrast, Portz, et al. [12], and Everett, et al. [6], utilize models that include serum PSA as a compartment that evolves dynamically in response to a combination of the AD cells, AI cells, intracellular androgen, and an independent clearance rate. Moreover, they include bidirectional mutation rates between the AD and AI compartments. Although these models include more details regarding the interaction of the compartments within the ITTA model, the ultimate conclusion in [6] is that a simpler model may provide more accurate information for predictive use.

Hirata, et al. [8], extend the ITTA model so that the androgen independent cells may exist in either a reversible or an irreversible state. During an on-treatment cycle, an androgen dependent cell may change to an androgen independent cell of either type and a reversible AI cell may change to an irreversible AI cell. During an off-treatment cycle, an AI cell that is in the reversible state may change to an AD cell. They show that their model provides a better description of the dynamics of prostate cancer under IAS than a model with only reversible or irreversible cell types. Suzuki, et al. [17], simplify the ITTA model and assume that the androgen levels are constant (i.e., at steady-state) during both on-and off-treatment cycles. This yields two linear, autonomous dynamical systems: one for the on-treatment phase and one for the off-treatment phase. They propose a region-dividing method which exploits the phase plane dynamics and saddle-point nature of the disease free equilibrium to control the switching between on-and off-treatment cycles and guarantee a cancer free state. They show that the region-dividing method effectively controls both the ITTA model and the Hirata model.

In this paper, we modify the ITTA model in order to create and analyze mathematical models for continuous, intermittent and periodic androgen suppression (CAS, IAS, and PAS respectively). Our models for IAS and PAS retain the qualitative behavior of the system under the feedback control defined in [9], yet are simple enough so that sensitivity analysis may be conducted to determine both the dynamic impact of parameters during treatment and the relative importance of parameter values on relapse time. Each of the studies [9], [18], and [17] include models that describe the prevention of relapse (i.e., the cancer free state), but we consider only those parameter values that allow for disease relapse as we are not aware of any clinical data that support prevention of relapse by IAS therapy. Although sensitivity to initial conditions can provide valuable information within a parameter estimation process, here we are focused on the relative importance of parameters under three different treatment options (CAS, IAS, and PAS) and we note the importance of collecting data at various times within a clinical trial that includes these options in order to improve estimates for those parameters.

We now summarize the basic structure and parameters of the ITTA model. First, it is assumed the administration is alternatively either present during on-treatment periods $(u=1)$ or absent during off-treatment periods $(u=0)$. The serum androgen concentration is denoted $a(t)$ with units of nmol/l. The normal androgen concentration is denoted by $a_0$ (nmol/l), which takes a value of $8\le a_0\le 35$ for typical adult males [9]. The speed of the recovery and decay of the androgen concentration is controlled by the parameter $\gamma$ $(\textrm{days}^{-1})$. The androgen dynamics are modeled as a linear, non-autonomous ordinary differential equation:

$\dot{a}(t) = -\gamma\left[a(t)-a_0\right]-\gamma a_0u(t). \label{androgen_dynamics}$

(1)

The tumor growth is described with changes in the numbers of AD and non-reversible AI cells, and we assume the proliferation and apoptosis rates of AD and AI cells are dependent on the androgen concentration $a(t)$. The tumor dynamics are represented by a system of nonlinear differential equations, where $x_1(t)$ and $x_2(t)$ represent the populations of AD and AI cells, respectively:

$\dot{x}_1(t) = \left[\alpha_1p_1(a(t))-\beta_1q_1(a(t))-m(a(t))\right]x_1(t)$

(2)

$\dot{x}_2(t) = m(a(t))x_1(t)+\left[\alpha_2p_2(a(t))-\beta_2\right]x_2(t)$

(3)

where

$p_1(a) = a(a+k_2)^{-1},$

(4)

$q_1(a) = k_3+(1-k_3)a(a+k_4)^{-1},$

(5)

$m(a) = m_1(1-a/a_0), \textrm{ and}$

(6)

$p_2(a) = 1-\kappa a/a_0. \label{Ideta_functions}$

(7)

Here we use the parameter $\kappa$ to distinguish between two different models investigated by Ideta, et al. [9]:

$\mbox{(i)} \, \kappa=0, \mbox{(ii)} \, \kappa=1-\beta_2/\alpha_2. \label{Ideta_kappa}$

(8)

The coefficients $\alpha_1, \beta_1$, and $\alpha_2,$ are parameters which depend on the bone metastatic site, i.e., the rate at which the tumor cells spread to the bones. The functions $\alpha_1 p_1$ and $\alpha_2 p_2$ represent the proliferation rates for the AD and AI cells, respectively. Similarly, $\beta_1 q_1$ and $\beta_2$ represent apoptosis rates for each compartment, where the apoptosis rate of the AI cells is assumed to be constant in time since these cells are unaffected by treatment. Finally, the function $m$ defines the rate at which AD cells mutate to the AI state, and the parameter $m_1$ denotes the maximum mutation rate.

The functional forms and parameters in (1)-(8) were defined in [9] for both qualitative and quantitative accuracy and, where possible, based on experimental data. For example, since AD cells do not multiply in the absence of androgen and the proliferation rate is known to be approximated by $\alpha_1=0.0204$ days$^{-1}$ for large values of $a$, the function $p_1(a)$ was chosen so that $p_1(0)=0$ and $p_1$ approaches 1 as $a$ increases. A more detailed discussion for all of these functional forms and parameters may be found in [9]. Parameter values, which assume a bone metastasis site, are summarized in Table 1, and the initial conditions used are

$a(0)=30, x_1(0)=15, x_2(0)=0.01 . \label{ICs}$

(9)

The suggested range of values for each parameter was created to stay within values available in the literature and in the cases where no such information was available, a standard deviation was used to ensure that the sample coefficient of variation remained under 30%. In Section 3 we further discuss the use of these parameter ranges to simulate data for sensitivity analysis.

Recall that the serum PSA is an effective biomarker for the prostate tumor growth. Its concentration is the only observable output of the system and is used as a basis for the intermittent administration in IAS therapy. Since a large amount of serum PSA is secreted by cancer cells, the PSA concentration $y(t)$ is assumed to be represented by the total sum of the subpopulations of cancer cells, i.e.,

$y(t) = x_1(t)+x_2(t). \label{PSA}$

(10)

The administration is suspended when the serum PSA concentration falls below $r_0$ (ng/ml) during on-treatment periods and it is re-instituted when the concentration exceeds $r_1$ (ng/ml) during off-treatment periods. Therefore, how to plan the IAS therapy or how to appropriately set the adjustable parameters $r_0$ and $r_1$ under the condition of $r_1> r_0>0$ is a significant issue for clinical practice. The model for IAS therapy in [9] consists of equations (1)-(10) coupled with an algorithmic feedback control determined by

$u(t) = \left\{ \begin{array}{lllll} 0\to1 & \textrm{when} & y(t)=r_1 & \textrm{and} & \dot{y} > 0, \\ 1\to0 & \textrm{when} & y(t)=r_0 & \textrm{and} & \dot{y} < 0. \end{array}\right. \label{u_dynamics}$

(11)

Based on this model, Ideta, et al., use numerical simulations and bifurcation analysis to show how tumor growth and relapse time are influenced by the net growth rate of the androgen-independent cells and the mutation rate from androgen-dependent cells to androgen-independent cells. In [12] and [6], a direct sensitivity to androgens for both the AD and AI compartments is assumed and their models include a response to intracellular androgen concentrations for both types of cells. However, the coupling of the differential equations (1)-(3) allows for indirect effects of androgen on even the AI compartment through changes in the parameter space. It is important to quantify these effects using sensitivity analysis in order to make a comparison with the outcomes of other models.

The remainder of this paper is organized as follows. In Section 2 we discuss the model under an assumption of continuous androgen suppression that will be used as a baseline for sensitivity analysis that will be extended to the intermittent and periodic models. In Section 3 we construct and analyze a dynamic model for IAS that treats the administration of anti-androgen therapy in (11) as a continuous state variable on the same time scale as the other compartments. In Section 4, we investigate a model for periodic androgen suppression (PAS) where $u$ is modeled as a periodic function of time. We compute both relapse times and on-treatment times as $r_0$ varies in the IAS model and the period varies in the PAS model. Floquet theory is used to determine feasible parameter regions for relapse to occur. We derive an approximate relapse time and use this expression to rank the relative sensitivity of the relapse time with respect to the model parameters. Finally, in Section 5 we summarize our main results and discuss model behavior, treatment protocol, and the medical significance of sensitivity analysis with regard to performing clinical trials and data collection.

---

2. Continuous androgen suppression

To effectively analyze the effects of Intermittent Androgen Suppression (IAS), we first look at a Continuous Androgen Suppression (CAS) therapy model in order to establish a baseline for which sensitivity analysis can be performed in a straightforward manner. Moreover, this treatment regimen is currently the only recommended option as a standard of care outside of Europe [4]. CAS is modelled by letting $u=1$ in equation (1). In Figure 1 we present the solutions under this assumption ($\kappa=0$) using the baseline parameter values presented in Table 1. Here the dashed line represents the androgen concentration, $a(t)$, which decays exponentially. The time constant for this decay is $\gamma^{-1}=12.5$ days so that the decay rate is relatively fast on a yearly time scale. The solid line is a plot of the concentration of AD cells, $x_1(t)$, which initially grows but then quickly decays as androgen is removed from the system. Finally, the dash-dot line illustrates the initial slow growth of concentration in AI cells, $x_2(t)$, followed by an eventual steep rise after the AD cells have been removed and all cells are in the AI state.

Given the coupled system of differential equations (1)-(3), an understanding of the influence that parameters within the model have on solutions is important since changes in the values of those parameters (e.g., the androgen concentration rate, $\gamma$) may be under investigation in a medical trial (e.g., how different pharmaceutical options affect the suppression of androgen). If a model is sensitive to a particular parameter value, then altering the value of that parameter may have a dramatic influence on the output of the model. One way to quantify this influence are the time-dependent sensitivity functions defined for a generic parameter $\theta$ and a generic differential equation, $\dot{x}=f(t, x)$, by a sensitivity equation that models the change in model output relative to changes in the parameter value [2]:

$\dot{s} = s(t)\frac{\partial f}{\partial x} + \frac{\partial f}{\partial\theta}, \label{sen_func}$

(12)

$s(0) = 0.$

(13)

Each equation in a dynamical system yields a corresponding sensitivity equation for each parameter within the model. Thus, for a model that has $n$ nonlinear equations with $m$ total parameters, the sensitivity functions are determined by solving $2nm$ equations. Each sensitivity equation is subject to trivial initial conditions. In the CAS model, an illustrative example is provided by the following sensitivity equations for each compartment with respect to the parameter $\gamma$:

$\dot{s}_{\gamma}^a = -\gamma s_{\gamma}^a-a, \\ \dot{s}_{\gamma}^{x_1} = \left(\alpha_1\frac{dp_1}{da}-\beta_1\frac{dq_1}{da}-\frac{dm}{da}\right)x_1 s_{\gamma}^a+[\alpha_1p_1(a)-\beta_1q_1(a)-m(a)] s_{\gamma}^{x_1}, \\ \dot{s}_{\gamma}^{x_2} = \left(\frac{dm}{da}x_1+\alpha_2\frac{d p_2}{da}x_2\right)s_{\gamma}^a+m(a)s_{\gamma}^{x_1}+[\alpha_2p_2(a)-\beta_2]s_{\gamma}^{x_2}, \\ $

where

$\frac{dp_1}{da}=\frac{k_2}{(a+k_2)^2}, \frac{dq_1}{da}=\frac{k_4(1-k_3)}{(a+k_4)^2}, \frac{dm}{da}=-\frac{m_1}{a_0}, \frac{d p_2}{da}=-\frac{\kappa}{a_0}.$

We note here that the derivative of $p_2$ with respect to $\alpha_2$ and $\beta_2$ will depend on the form of $\kappa$. Solutions to these equations provide important information about the original system when the chosen parameter value is perturbed about some a priori fixed value of that parameter. Here we use the baseline values in Table 1 for those fixed values. A large magnitude in the sensitivity function indicates a point in time when data collection will be most informative during a clinical trial for estimating a particular parameter [3]. Furthermore, a spike in the sensitivity function provides information that is helpful in producing high quality numerical solutions for the model. Whenever parameter sensitivity is high for any parameter in the model, increasing the length of subdivisions in the discretization in time may lead to inaccurate solutions at that point in time and at future times [16].

The sensitivity solution for androgen concentration with respect to the parameter $\gamma$ can be seen in Figure 2 There is a spike in the graph between $t=0$ and $t=50$ days where the magnitude of sensitivity increases significantly, indicating that the androgen concentration is sensitive with respect to $\gamma$ during this time period. No significant difference between the cases $\kappa=0$ and $\kappa=1-\beta_2/\alpha_2$ was observed, so in this section we only present the results for the first case. The negative values for the sensitivity function in Figure 2 indicate that as the value of $\gamma$ is increased by a small amount, the response of the androgen concentration compartment is decreasing. Data should be collected at times during the 0-50 day interval to optimally determine the recovery/decay rate of the androgen concentration ($\gamma$). The sensitivity of androgen concentration with respect to other parameters remains constant as in this case the androgen dynamics depending only on $\gamma$.

Next, we consider the sensitivity of the androgen dependent cells to changes in the parameter space. In Figure 3 we present the sensitivity functions for the AD population with respect to the parameters $\alpha_1$ and $\beta_1$, respectively. Here we see a very large magnitude of sensitivity to $\alpha_1$ and to $\beta_1$ during the first four months of treatment. Further analyzing the extrema of these graphs, an initial surge of sensitivity to the AD proliferation rate occurs first, followed by a slightly longer period of higher sensitivity to the apoptosis rate. The sensitivity of the AD population to the parameters $\gamma$ and $m_1$ mimics that of the sensitivity of this compartment to the apoptosis rate, but at a lower magnitude.

In contrast, sensitivity of the AI population to changes in the parameter space remains relatively low until many months of treatment have elapsed. In Figure 4 we see that sensitivity of this compartment to the androgen concentration rate $\gamma$ has magnitude similar to that of the effects of this parameter on the androgen concentration, but now the effects are near the end of our simulation time frame. Likewise, the sensitivity of the AI population to the proliferation and apoptosis rates for the AD population ($\alpha_1$ and $\beta_1$, respectively) is comparable in magnitude and direction with the respective sensitivities of the AD population to these parameters, but again the steep increase in sensitivity is delayed until over a year of continuous treatment has passed. This time frame surpasses the typical induction period of ADT during an administration of IAD [15]. The sensitivity of the AI population to the proliferation and apoptosis rates for the AI population ($\alpha_2$ and $\beta_2$, respectively), behave in much the same way but at a higher magnitude due to their direct impact on the AI population itself. Moreover, the sensitivity of the AI population to the maximum mutation rate follows this pattern and achieves the highest magnitude of all the sensitivity functions, as shown in Figure 5. As a result, data collection efforts aimed at estimating the parameters in this model must remain consistent in frequency for roughly three full years of CAS treatment, even though a relapse is likely to have occurred.

In Section 5, we further explore the CAS model in order to examine the effect of our choice for the value of $\kappa$ on the length of time before disease relapse occurs, along with a similar analysis for models that include intermittent treatment. Having here established a baseline for sensitivity of the model compartments to the parameter space, we now consider how that sensitivity changes under a new model for IAS.

---

3. Modeling intermittent androgen suppression

In order to compare the parameter sensitivity in this scenario with what we computed for CAS treatment, we first reformulate the algorithmic definition as an ordinary differential equation. First, we consider the switch from on-treatment to off-treatment (from $u=1$ to $u=0$) represented in the second half of (11). We model the requirement $y(t)=r_0$ as a product of step functions:

$H(y-(r_0-\epsilon))H(r_0+\epsilon - y), \label{y_r0}$

(14)

where $\epsilon>0$ and

$H(y-K) = \left\{ \begin{array}{cc} 0 & \textrm{if } y \le K, \\ 1 & \textrm{if } y > K. \\ \end{array}\right. \label{step function}$

(15)

Because of numerical inaccuracies, $y(t)=r_0$ may not be achieved on a particular interval of time. Therefore, a small value of $\epsilon$ is used to account for small differences in $|y(t)-r_0|$ during a simulation. In order to meet the condition that $\dot{y} <0$, We multiply (14) by $H(-\dot{y})$.

To create a smooth transition from $u=1$ to $u=0$, we use a logistic function with stable equilibrium $u=0$ and unstable equilibrium $u=1+\hat{\epsilon}$,

$f_{10}(u)=u(u-(1+\hat{\epsilon})), \label{u_1}$

(16)

where $\hat{\epsilon}$ represents a different small positive number. Thus, if $0 <u(t) <1+\hat{\epsilon}$ at some time $t$, then $u$ will approach zero.

Finally, the dynamics of $u$ for the switch from on-treatment to off-treatment ($u=1$ to $u=0$) are modeled by combining these terms:

$\dot{u}_{10}=H(y-(r_0-\epsilon))H(r_0+\epsilon - y)H(-\dot{y})u(u-(1+\hat{\epsilon})). \label{u_10}$

(17)

Similarly, we model the switch from off-treatment to on-treatment (i.e., from $u=0$ to $u=1$) shown by the first half of (11) as

$\dot{u}_{01}(t)=H(y-(r_1-\epsilon))H(r_1+\epsilon - y)H(\dot{y})(u+\hat{\epsilon})(1-u). \label{u_01}$

(18)

We model the back and forth switching in (11) by adding (17) and (18). Hence, we have the following ordinary differential equation:

$\dot{u} = \lambda[H(y-(r_0-\epsilon))H(r_0+\epsilon-y)H(-\dot{y})u(u-(1+\hat{\epsilon})) \\ +\: H(y-(r_1-\epsilon))H(r_1+\epsilon-y)H(\dot{y})(u+\hat{\epsilon})(1-u)]. \\ $

(19)

Here we use a parameter $\lambda>0$ to control the rate at which $u$ transitions from 1 to 0 and vice versa. The value of the parameter $\lambda$ is chosen using numerical experimentation. For example, we found that for $\lambda=1000,$ $u$ oscillates but becomes negative. We simulated (19) using many different values and found $\lambda=1100$ to accurately represent the on and off switching in the model of Ideta, et al., as shown in Figure 6. In Figure 7, we present the behavior of the serum PSA level $y(t)$ using this model. The parameters $\epsilon$ and $\hat{\epsilon}$ are chosen small enough so that our simulations demonstrate the same qualitative behavior as that produced by the algorithmic model (1)-(11) used in [9] (see Figures 8(a) and 9(a)). They both produce a PSA that oscillates between $y=5$ ng/ml and $y=20$ ng/ml and then increases without bound, representing a relapse.

To accurately conduct parameter sensitivity computations, we need to ensure that the right side of (19) is differentiable. Thus, we will use the arctangent function to approximate the step functions in (19). In particular, if $z$ depends on a parameter $\theta$, we define

$G[z(\theta)] = \frac{1}{\pi}\tan^{-1}[10000 \, z(\theta)]+\frac{1}{2},$

(20)

so that

${d \over d\theta}G[z(\theta)] = \frac{10000}{\pi}(1+\left[10000 \, z(\theta)\right]^2)^{-1} {dz \over d\theta},$

(21)

and hence continuity of the derivative of $z$ will ensure that $G$ is differentiable. We replace each instance of $H$ in (19) with an appropriate form of $G$. Our model for IAS then becomes

$\dot{a} = f^{(1)}(a)+g^{(1)}(u), \label{IASandrodyn}$

(22)

$\dot{x_1} = f^{(2)}(a),$

(23)

$\dot{x_2} = f^{(3)}(a),$

(24)

$\dot{u} = g^{(2)}(u).$

(25)

where

$f^{(1)}(a) = -\gamma(a(t)-a_0)$

(26)

$g^{(1)}(u) = -\gamma a_0u(t),$

(27)

$f^{(2)}(a) = [\alpha_1p_1(a(t))-\beta_1q_1(a(t))-m(a(t))]x_1(t),$

(28)

$f^{(3)}(a) = m(a(t))x_1(t)+[\alpha_2p_2(a(t))-\beta_2q_2(a(t))]x_2(t),$

(29)

$g^{(2)}(u) = \lambda[G(y-(r_0-\epsilon))G(r_0+\epsilon-y)G(-\dot{y})u(u-(1+\hat{\epsilon})) \\ +\: G(y-(r_1-\epsilon))G(r_1+\epsilon-y)G(\dot{y})(u+\hat{\epsilon})(1-u)].$

(30)

We use the following system for the sensitivity to each parameter $\theta$, where $\dot{s}_{\theta}^a$, $\dot{s}_{\theta}^{x_1}$, $\dot{s}_{\theta}^{x_2}$, and $\dot{s}_{\theta}^u$ determine the sensitivity functions for $a(t)$, $x_1(t)$, $x_2(t)$, and $u(t)$, respectively.

$\dot{s}_{\theta}^a = f_a^{(1)}(a) s_\theta^a(t)+ f_\theta^{(1)}(a)+g_u^{(1)}(u) s_\theta^u(t) +g_\theta^{(1)}(u),$

(31)

$\dot{s}_{\theta}^{x_1} = f_{x_1}^{(2)}(a) s_\theta^{x_1}(t) + f_\theta^{(2)}(a),$

(32)

$\dot{s}_{\theta}^{x_2} = f_{x_2}^{(3)}(a)s_\theta^{x_2}(t) + f_\theta^{(3)}(a),$

(33)

$\dot{s}_{\theta}^{u} = g_u^{(2)}(u) s_\theta^{u}(t) + g_\theta^{(2)}(u).$

(34)

with initial conditions

$s_{\theta}^a(0)= s_{\theta}^{x_1}(0)= s_{\theta}^{x_2}(0)=s_{\theta}^u(0)=0. \label{sensitivity system initial}$

(35)

The sensitivity solutions for the parameter $\gamma$ can be seen in Figure 8. Sensitivity of $a(t)$ to $\gamma$ spikes as we shift from on-treatment to off-treatment and vice-versa, with negative sensitivity during the on-treatment times (as in the CAS model) and positive sensitivity during the off-treatment times. Note that the magnitude of these spikes increases over time and that this behavior extends further in time for the case $\kappa=1-\beta_2/\alpha_2$ since $a(t)$ continues to oscillate well past 700 days, the approximate time at which this sensitivity decays to zero for $\kappa=0$. Similar behavior is observed for sensitivity of $x_1(t)$ to $\gamma$. However, note that the negative sensitivity during on-treatment periods of time has a much larger magnitude and the rate of change in sensitivity is much larger than positive sensitivity during off-treatment periods. Finally, $x_2(t)$ remains mostly insensitive to $\gamma$ until after relapse has occurred, similar to the behavior observed under CAS but with small oscillations occurring between on-treatment and off-treatment times.

In contrast, the sensitivities for the mutation parameter $m_1$, shown in Figure 9 remain relatively low until after a year of treatment has elapsed, when mutation has a chance to have a significant effect on the AD compartment. Afterwards, sensitivity of $a(t)$ spikes in the positive direction when switching to on-treatment and in the negative direction when switching to off-treatment. These spikes continue to grow in magnitude over time until relapse occurs. For the AD compartment, $x_1(t)$, the expected negative sensitivity to $m_1$ during periods of off-treatment steadily decreases over time, interrupted by increasingly large spikes of positive sensitivity when switching to on-treatment periods. Sensitivity for the AI compartment, $x_2$, remains relatively low for the first year of treatment and then eventually grows exponentially. Similar behavior is observed for many other parameters in the model, but the magnitude of spikes in sensitivity to $m_1$ were the most significant.

In Figure 10 we see the sensitivities to the parameter $r_0$, which determines when off-treatment periods may begin. Although the behavior here is similar to that for $\gamma$, we note the asymmetry in the magnitude for sensitivities of $a(t)$ and $x_1(t)$ during periods of on-treatment, which are negative, compared with periods of off-treatment when these sensitivities are positive. Moreover, note that sensitivity spikes are possible even after treatment has been permanently turned on once relapse has occurred. Somewhat surprisingly, the sensitivities to the parameter $r_1$, which determines when on-treatment periods must begin, were all very low in magnitude at all times.

Thus, we generally expect the same behavior of our sensitivity functions during periods of on-treatment as are observed in the CAS model, but an opposite behavior is generally observed during periods of off-treatment. These swings in sensitivity generally spike at points in time where treatment is turned either on or off. For $a(t)$ and $x_1(t)$, sensitivities tend to decay after relapse has occurred, while for $x_2(t)$ sensitivities tend to remain low until after relapse. This indicates the importance of continuing to collect data in a clinical trial of IAS even after the decision has been made to remain on-treatment permanently. Additionally, the magnitude of the spikes in sensitivities for $a(t)$ and $x_1(t)$ frequently double within a year, as seen in Figures 8 and 9, indicating an increase in the importance of collecting data at these switching points. However, as seen in Figure 10, we may reach a point in time of maximum sensitivity for $a(t)$ and/or $x_1(t)$ many months prior to relapse.

We also examine the relative sensitivity of our model to each parameter in order to determine which parameters have the biggest impact on the long-term behavior of our solution components. Using the range of parameter values presented in Table 1 and an assumption that each parameter was normally distributed in its range, independent of the other parameter values, 100 randomized parameter sets were generated and used to generate simulated data sets $\hat{\mathbf{z}} = (\hat{a}, \hat{x}_1, \hat{x}_2, \hat{u})^T$ at the same points in time as our computed solutions using the baseline values of our parameters, $\mathbf{z}(t) = (a(t), x_1(t), x_2(t), u(t))^T$. We define a cost function

$J(\mathbf{q})=\sum\limits_k||\mathbf{z}(t_k;\mathbf{q}_0)-\hat{\mathbf z}(t_k;\mathbf{q})||_{{2}}^{{2}},$

representing a squared error in the Euclidean norm of the difference between our solutions generated using the baseline parameters (represented by $\mathbf{q}_0$) and solutions generated with one of the randomized parameter sets (represented by $\mathbf{q}$). Hence,

$\frac{\partial J}{\partial \mathbf{q}}(\mathbf{q}_0)=2\sum\limits_k\left(\mathbf{z}(t_k;\mathbf{q}_0)-\hat{\mathbf z}(t_k;\mathbf{q})\right)^T\frac{\partial \mathbf{z}}{\partial \mathbf{q}}(t_k;\mathbf{q}_0).$

For each randomized parameter set, we calculate the relative sensitivity to a given parameter $\theta$ as

$\Upsilon^J_{\theta} := \frac{\partial J}{\partial\theta}\times\frac{\theta}{J}.$

In Figure 11, we present a summary of these values for the case $\kappa=0$. The parameters $r_0$, $\alpha_2$, and $\beta_2$ had the highest median relative sensitivity, in magnitude, and the widest range of relative sensitivities over the simulated data sets. Similar results are presented for the case $\kappa=1-\beta_2/\alpha_2$ in Figure 12, and in Table 2 we present the median values of our computed relative sensitivities for both cases. The model is also somewhat sensitive to $\alpha_1$ and $\beta_1$ for both cases, and is sensitive to $m_1$ when $\kappa=0$. We note that other than $r_0$, the list of parameters for which our model is most sensitive depend a great deal on the site of metastasis. We expect that data from clinical trials would need to be carefully collected in attempts to better estimate values of these parameters.

---

4. Periodic androgen suppression and relapse time

In Ideta, et al., Figures 8 and 9 show simulations of the IAS model for $\kappa$ as in (8). With all other parameters in the model fixed, the time to relapse is observed to be dependent upon $r_0$. In this section, we formulate a periodic androgen suppression (PAS) model in order to investigate the relapse time and its sensitivity to parameters in the model. We use a simple definition of relapse time:

$T_R=\left\{t>0 \left| \right. y(t)=\overline{y}\right\} \label{relapse_time},$

(36)

where $\overline{y}$ is a pre-determined threshold value of total PSA, $y(t)$. We note that a clinical definition of relapse, for example as described in [9], is different from this mathematical definition.

The IAS model (1)-(10) controls the treatment indicator function $u$ through the algorithmic feedback control (11). When $r_1=15$, a wide range of values of $r_0$ results in an initial periodic-like behavior in $y$. Disease relapse occurs as these oscillations cease, the AD cells are extinguished, and the AI cells increase without bound (see Figure 7 above and Figures 8 and 9 in [9]). Figure 13 shows the relapse time as $r_0$ varies using the IAS model with on-and off-treatment $u$ controlled by algorithmic feedback (11) and our model (19). This figure also shows the total time the patient is on-treatment which we denote $T_{on}$. We can compute $T_{on}$ using the treatment indicator function $u(t)$ as follows

$T_{on}=\int_0^{T_R} u(\tau) \, d\tau. \label{on_treatment_time}$

(37)

It can be seen in Figure 13 that in general $r_0$ should be small (when all other parameters are fixed as in Table 1) in order to maximize relapse time and minimize on-treatment time. For $\kappa=0$, the maximum difference between relapse time and on-treatment time occurs at $r_0 \approx 1.5$, when relapse time is approximately $818$ days and on-treatment time is approximately $362$ days. For $\kappa=1-\beta_2/\alpha_2$, the maximum difference between relapse time and on-treatment time occurs at $r_0 \approx 2.6$, when relapse time is approximately $1709$ days and on-treatment time is approximately $755$ days.

The periodic nature of $y$ observed prior to relapse motivates our consideration of periodic cycling of androgen suppression. Specifically, we model $u$ as a T-periodic function such that $u(t)=\hat{u}(t+T)$ for all $t \ge 0$, where

$\hat{u}(t) = \left\{ \begin{array}{lllll} 1 & \textrm{when} & 0 \le t <T/2, \\ 0 & \textrm{when} & T/2 \le t <T. \end{array}\right.$

(38)

The PAS model consists of equations (1)-(8) with treatment control (38). After substituting (38) into (1), the androgen dynamics are $a(t)=\hat{a}(t+T)$ for all $t \ge 0$ where

$\hat{a}(t) = \left\{ \begin{array}{lllll} a_0e^{-\gamma t}, & \textrm{when} & 0 \le t <T/2, \\ a_0e^{-\gamma t}+a_0\left[1-e^{-\gamma (t-T/2)}\right], & \textrm{when} & T/2 \le t <T, \end{array}\right.$

(39)

with $a(0)=a_0$ as in Ideta, et al. Next we substitute (39) into (2) and (3) to get a linear system of two equations with periodic coefficients. This system has the matrix form ${dX/dt}=A(t)X$, with $X=(x_1, x_2)^T$ and

$A(t)=\left[\begin{array}{cc} \alpha_1p_1(a)-\beta_1q_1(a)-m(a) & 0 \\ m(a) & \alpha_2p_2(a)-\beta_2 \end{array} \right]. \label{A_define}$

(40)

Denoting initial conditions as $X(0)=X_0$, the solution is determined by standard methods [7] to be $X(t)=\Phi(t)X_0$ where $\Phi(t)$ is the fundamental matrix which solves $d\Phi/dt=A(t)\Phi$, subject to $\Phi(0)=I$, with $I$ the $2 \times 2$ identity matrix. That is, the columns of $\Phi$ are determined by solving ${dX/dt}=A(t)X$ subject to $X(0)=(1, 0)^T$ and $X(0)=(0, 1)^T$. Straightforward calculations yield

$\Phi(t)=\left[\begin{array}{cc} \exp\left(\int_0^t \left[f_1(\tau)-m(\tau)\right] \, d\tau\right) & 0 \\ \psi(t) & \exp\left(\int_0^t f_2(\tau) \, d\tau\right) \end{array} \right], \label{Phi_define}$

(41)

where

$f_1(\tau)= \alpha_1p_1(a(\tau))-\beta_1q_1(a(\tau)), \label{f1_def}$

(42)

$f_2(\tau)= \alpha_2p_2(a(\tau))-\beta_2, \label{f2_def}$

(43)

$\psi(t) = e^{\int_0^t f_2(\tau) \, d\tau} \int_0^t m(a) e^{\int_0^t \left(f_1(\tau)-m(\tau)-f_2(\tau)\right) \, d\tau} \, d\tau.$

(44)

Figure (14) depicts the serum PSA $y(t)=x_1(t)+x_2(t)$ computed from the fundamental matrix solution (41) when $T=200$ and using the values of $\kappa$ in (8). The serum PSA level typically rises as $a(t)$ rapidly decays towards zero and increases (in the case of $\kappa=0$) or remains relatively constant (in the case of $\kappa=1-\beta_2/\alpha_2$) when $a(t)$ is near its threshold value of $a_0$.

Simulation of the PAS model reveals that for $T$ smaller than a certain threshold value, the AD cell population may oscillate with an increasing average value. Typically, in the IAS model disease relapse occurs as AD cells mutate to AI cells and the AD cell population becomes negligible so that eventually all cells are in the AI state. In this case, $x_1(t)$ and $x_2(t)$ obey

$\lim\limits_{t \rightarrow \infty} x_1(t)=0, \lim\limits_{t \rightarrow \infty} x_2(t)=\infty. \label{x1_x2_conditions}$

(45)

The following theorem details conditions under which (45) holds for the PAS model.

Theorem 4.1. If the parameters are chosen so that

$\int_0^T \left[f_1(\tau)-m(\tau)\right] \, d\tau <0 \; and \; \kappa <\rho^{-1}(1-\beta_2/\alpha_2), \label{lambda_conditions}$

(46)

where $\rho=1/2+(\gamma T)^{-1}$ then (45) holds.

Proof. The equation for $x_1(t)$ is given by model equation (2), a first order, linear, homogenous differential equation with a periodic coefficient $f_1(t)$ defined in (42). Solving this problem with initial condition $x_1(0)$ and using Floquet theory [7], [19] yields that

$x_1(t)= x_1(0)\exp\left(\int_0^t \left[f_1(\tau)-m(\tau)\right] \, d\tau \right)=p_1(t) \exp\left({\zeta_1 t}\right) \label{x1_Floquet}$

(47)

where

$\zeta_1=T^{-1}\int_0^T \left[f_1(\tau)-m(\tau)\right] \, d\tau, \label{zeta1_def}$

(48)

and $p_1(t)$ is a $T$ periodic function. Thus assuming the first inequality in (46) holds, then the limit for $x_1(t)$ in (45) must hold.

We next determine $x_2(t)$ from the model equation (3) which is a linear, first order, non-homogeneous differential equation. The solution of (3) with initial condition $x_2(0)$ (which may be determined by the integrating factor method), is that

$x_2(t)=x_2(0)\exp\left(\int_0^t f_2(\tau) \, d\tau\right)+\psi(t) \label{x2_sol}$

(49)

with $f_2(t)$ defined in (43) and $\psi(t)$ as in (44). The first term on the right side of (49) solves the homogeneous problem $dx_2/dt=f_2(t)x_2(t)$ with periodic coefficient function $f_2(t)$, which again after applying Floquet theory, may be expressed as

$x^h_2(t)= x_2(0)\exp\left(\int_0^t f_2(\tau) \, d\tau\right)=p_2(t) \exp\left({b_2 t}\right) \label{x2_Floquet}$

(50)

where $b_2=T^{-1}\int_0^T f_2(\tau) \, d\tau$ and $p_2(t)$ is a $T$ periodic function. After substituting (43) and integrating, we get that

$b_2=T^{-1}\int_0^T f_2(\tau) \, d\tau=\alpha_2(1-\kappa\bar{a}/a_0)-\beta_2, \label{lambda_conditions_a}$

(51)

where

$\bar{a}=T^{-1}\int_0^T a(\tau) \, d\tau={a_0 \over \gamma T}\left(e^{-\gamma T/2}-e^{-\gamma T}\right)+{a_0 \over 2}. \label{average_a}$

(52)

Using the bound $\bar{a}/a_0 <1/2+(\gamma T)^{-1}$ it follows from (51) that $b_2>0$ when the second inequality in (46) holds. Thus we see from (50) that $x_2^h(t) \rightarrow \infty$ as $t \rightarrow \infty$. Inspection of (39) and (50) shows that $m(a(t))x_1(t)>0$ and it follows that $x_2(t)>x_2^h(t)$ for $t>0$ (see Corollary 6.3 of Theorem 6.1 in [7]). This implies that the second limit in (45) must hold.

In order to determine a parameter set such that the inequalities (46) hold, we first compute $\zeta_1$ from (48) to get

$\begin{array}{cc} \zeta_1=\gamma^{-1}T^{-1} [-\alpha_1v_2-\alpha_1 w_{12}+a_0\alpha_1 c_2^{-1}s_{21}-k_3 \beta_1 \gamma T+(1-k_3)\beta_1v_4&\\ +\beta_1(1-k_3) w_{13}-a_0 \beta_1(1-k_3) c_3^{-1}s_{31}-m_1 \gamma T/2&\\-m_1(e^{-\gamma T/2}-1) + m_1(e^{\gamma T/2}-1)(e^{-\gamma T}-e^{-\gamma T/2})], \label{zeta1_compute} & \end{array}$

(53)

where

$v_i=\ln\left({a_0e^{-\gamma T/2}+k_i \over a_0+k_i}\right), \\ w_{ij}=\ln\left({ c_i e^{-\gamma T}+c_j \over c_i e^{-\gamma T/2}+c_j}\right), \\ s_{ij}=\ln\left({ c_i e^{\gamma T}+c_j \over c_i e^{\gamma T/2}+c_j}\right),$

and

$c_1=a_0\left(1-e^{\gamma T/2}\right), c_2=a_0+k_2, c_3=a_0+k_4.$

Numerically, we find a critical value $T_c\approx 123$ at which $\zeta_1$ changes sign. The value of $\zeta_1$ is positive for $T <T_c$ and negative for $T>T_c$ assuming all other parameters in (53) are fixed at the baseline values in Table 1. To satisfy the inequalities (46) we let $T>T_c$ so that $\zeta_1 <0$, ensuring the first inequality holds, and choose $\kappa$ small enough so that the second inequality holds. Note that the second case for $\kappa$ in (8) is small enough to satisfy this condition.

We next use the results in Theorem 4.1 to derive a closed form expression that approximates the relapse time $T_R$ given by (36) for the PAS model and then compute the parameter sensitivity of $T_R$ using this approximation. First we note that for $t < O(\gamma^{-1})$ it can be seen that $a(t) \approx a_0$ and $x_1(t)$ initially increases at a rate $\alpha_1 p_1(a_0)-\beta_1 q_1(a_0)$. The results in Theorem 4.1 show that $x_1(t)$ may be expressed as in (47) where $\zeta_1$ is defined in (48) and computed in (53). Since $p_1(t)$ is $T$-periodic, it follows that $\left|\zeta_1\right|$ controls the decay rate of $x_1(t)$, assuming (46) hold so that $\zeta_1<0$. To approximate the relapse time for the periodic model when $t > O(\gamma^{-1})$ we then assume $x_1(t) \propto e^{-\left|\zeta_1\right|t}$ on this interval. With these approximations the solution of (2) is

$\hat{x}_1^{(1)}(t)=x_1(0)e^{t/\tau_1}, 0 \le t < \gamma^{-1}, \textrm{ and} \label{x1cont_1}$

(54)

$\hat{x}^{(2)}_1(t)=\hat{x}^{(1)}_1(\gamma^{-1}) \, e^{-\left|\zeta_1\right|(t-\gamma^{-1})}, \; t \ge \gamma^{-1},$

(55)

where

$\tau_1^{-1}=\alpha_1 p_1(a_0)-\beta_1 q_1(a_0) \\ =\alpha_1 a_0(a_0+k_2)^{-1}-\beta_1\left[k_3+(1-k_3)a_0(a_0+k_4)^{-1}\right]. \label{tau1_def}$

(56)

We next approximate $x_2(t)$ on $0 \le t < T/2$ assuming $\gamma^{-1} < T/2$. In this case, we let $a=a_0$ on $0 \le t < \gamma^{-1}$ and $a=0$ on $\gamma^{-1} \le t < T/2$. We solve (3) with $x_1(t)$ as in (54) and (55) to get that

$\hat{x}^{(1)}_2(t)=\left(x_2(0)-\theta_1\right)e^{b_1 t}+\theta_1e^{t/\tau_1}, \; 0 \le t < \gamma^{-1}, \textrm{ and}$

(57)

$\hat{x}^{(2)}_2(t)=\left(\hat{x}^{(1)}_2(\gamma^{-1})-\theta_2\right)e^{b_2(t-\gamma^{-1})} \\ + \theta_2 e^{-\left|\zeta_1\right|(t-\gamma^{-1})}, \; \gamma^{-1} \le t < T/2, \label{x2per_2}$

(58)

where the constants are defined as

$b_1=\alpha_2 p_2(a_0)-\beta_2=\alpha_2(1-\kappa)-\beta_2, b_2=\alpha_2 p_2(0)-\beta_2=\alpha_2-\beta_2,$

(59)

$\theta_1=m_1 x_1(0)\left[\tau_1^{-1}-b_1\right]^{-1}_, \theta_2=m_1 \hat{x}^{(1)}_1(\gamma^{-1})\left[\left|\zeta_1\right|-b_2\right]^{-1}_.$

(60)

For $\kappa=0$ the constants $b_1$ and $b_2$ are equal and $b_1=b_2=\alpha_2-\beta_2$. In this case we may extend the interval in (58) to $t \ge \gamma^{-1}$. If $T_R-\gamma^{-1}>>\left|\zeta_1\right|^{-1}$ and $m_1<<1$ we approximate $T_R$ by ignoring the second term in (58) and solving $x_2(t) \approx y(t)=\overline{y}$ to get that

$T_R \approx \, {1 \over \alpha_2-\beta_2}\ln \left( {\overline{y} \over \hat{x}^{(1)}_2(\gamma^{-1})-\theta_2}\right)+\gamma^{-1}. \label{periodic_relapse_kappazero}$

(61)

For $\kappa=1-\beta_2/\alpha_2$ it follows that $b_1=0$ when $a=a_0$ so that $x_2(t)$ is constant on $\left(n+1\right)T/2 \le t < (n+1)T$, $n=0, 1, 2, \ldots$, and the solution of (3) with $x_1(t)$ as in (54) and (55) is

$\tilde{x}^{(1)}_2(t)=\hat{x}^{(2)}_2(T)\exp\left\{(\alpha_2-\beta_2) \, \left(t-{n+1 \over 2}T\right)\right\}, \\ nT \le t < \left(n+{1 \over 2}\right)T, \textrm{ and}$

(62)

$\tilde{x}^{(2)}_2(t)=\hat{x}^{(2)}_2(T)\exp\left\{ (\alpha_2-\beta_2) \, {1 \over 2} \, n T\right\}, \\ \left(n+{1 \over 2}\right)T \le t < (n+1)T,$

(63)

where $n$ is any natural number. To approximate the relapse time, we solve $x_2(t)=\overline{y}$ and using (62), (63) for $x_2(t)$ yields that

$T_R \approx {1 \over \alpha_2-\beta_2}\ln\left({\overline{y} \over \hat{x}^{(2)}_2(T)}\right)+{N_c+1 \over 2}T,$

(64)

where

$N_c= \left \lceil {2 \over T (\alpha_2-\beta_2)}\ln\left({\overline{y} \over \hat{x}^{(2)}_2(T)}\right) \right \rceil.$

(65)

Figure 15 shows comparisons of the actual and approximated relapse times from (61) and (64) as $T$ varies and $\overline{y}=20$. The re-setting behavior that occurs in the relapse time for $\kappa=1-\beta_2/\alpha_2$ is a result of the plateaus in $y(t)$ during the on-treatment cycles when $a(t) \approx a_0$ as observed in Figure 14. If the period is such that $y(t)$ plateaus at $\overline{y}$, then the relapse time re-sets to a smaller value. In addition, Figure 15 shows the on-treatment time from (37). Here it can be seen that $T$ should be small in order to maximize relapse time and minimize on-treatment time, which mimics the response of relapse time to changes in $r_0$ as observed in Figure 13. For $\kappa=0$, the maximum difference between relapse time and on-treatment time occurs at $T \approx 206$ days, when relapse time is approximately $824$ days and on-treatment time is approximately $413$ days. For $\kappa=1-\beta_2/\alpha_2$, the maximum difference between relapse time and on-treatment time occurs at $T \approx 409$ days, when relapse time is approximately $1643$ days and on-treatment time is approximately $825$ days.

To approximate the relapse time for the continuous model, where we can solve (1) analytically, we have $a(t)=a_0e^{-\gamma t}$ for $t \ge 0$ and we note that this is simply the limiting case $T\rightarrow \infty$ for the periodic model. The approximations of $x_2(t)$ in (57) and (58) depend on $T$ through the parameter $\zeta_1$, and using (53) we get that

$\lim\limits_{T\rightarrow \infty} \, \left|\zeta_1\right|=\beta_1\left(k_3+{1 \over 2}c_3^{-1}a_0(1-k_3)\right)-{1 \over 2}c_2^{-1}\alpha_1-{1 \over 2}m_1. \label{zeta1_limit}$

(66)

Figure 16 shows the relapse time for $\kappa \in (0, 1)$ estimated from (61) with $\left|\zeta_1\right|$ replaced by (66) and also shows the relapse time determined by solving the continuous model directly (using an integrating factor method with $a(t)=a_0e^{-\gamma t}$ and Matlab's adaptive Lobatto quadrature routine to perform the integration). Note that the relapse time increases with $\kappa$ as $b_1$ depends on $\kappa$.

For a generic parameter $\theta$, we define the relative sensitivity index [5] in the relapse time as

$\Upsilon_{\theta}^{T_R} := {\partial T_R \over \partial \theta} \times {\theta \over T_R}.$

(67)

We use the sensitivity index to investigate the sensitivity of the relapse time to the parameters in the model and summarize these results in Figures 17 and 18 and Table 3. For the periodic model we compute $\Upsilon_{\theta}^{T_R}$ using 100 random draws of the period $T$ distributed normally. Here we use a mean value of 425 days and a standard deviation of 25 days so that we expect most of the data to be in the range $350 \le T \le 500$. The lower bound of this range is well above the critical value $T_c = 123$ for model feasibility that is guaranteed by the conditions in Theorem 1. Inspection of Figure 15 shows that this range of $T$ is in the optimal treatment region that ensures relatively large relapse times and relatively short on-treatment cycles.

For the continuous model we compute $\Upsilon_{\theta}^{T_R}$ using 100 random draws of $\kappa$ distributed normally. Here we use a mean value of 0.5 and a standard deviation of 0.165 so that we expect most of the data to be in the range $0.005 \le \kappa \le 0.995$, corresponding (approximately) to the range of $\kappa$ depicted in Figure 16.

Table 3 displays the mean values of the sensitivity indices. In all cases, the mean and median values were very close. These values inform us as to how changes in each parameter influence relapse time. For example, if $\kappa=0$ (with all other parameters fixed at their baseline values), increasing $\alpha_2$ by 1 % decreases the relapse time by 3.12 %. Figures 17 and 18 and Table 3 show that the magnitude of the sensitivity index of $T_R$ to both $\alpha_2$ and $\beta_2$ are significantly larger than the other sensitivity indices in each case. In the periodic model for $\kappa=1-\beta_2/\alpha_2$, the sensitivity index to the period is also relatively large.

---

5. Conclusions

In this paper we have investigated mathematical models that describe the treatment of advanced prostate cancer by androgen suppression. Our dynamical system models, derived from the hybrid dynamical system presented in [9], allow us to simulate continuous androgen suppression (CAS), intermittent androgen suppression (IAS), and periodic androgen suppression (PAS). These models describe the behavior of the androgen dynamics as well as the dynamics of both androgen dependent (AD) and non-reversible androgen independent (AI) cancer cells. We alter feedback mechanism for the treatment control function $u$ from [9] (shown in (11)), using a differential equation model in equation (19) for IAS and a periodic model of $u$ in equation (38) for PAS. With these models we analyze parameter sensitivity and relapse time for the two main values of $\kappa$ in (8), where $\kappa$ controls the proliferation rate of the AI cells.

Our dynamic sensitivity analysis for the CAS model produced an expected result -compartments in the model were highly sensitive to parameters influencing their growth and decay, but only during times when the solution variables for those compartments had a large enough magnitude that we could expect their presence to be measurable in a clinical experiment. When considering the dynamic sensitivities for the IAS model, this pattern continues in an expected way as we switch between on-treatment and off-treatment periods of time. However, we learned that the spikes in the sensitivity of androgen ($a$) and the AD cell population ($x_1$) at these switching points actually increases in magnitude over time for most of our parameters until $a$ and $x_1$ permanently decay in value, indicating a relapse. Simultaneously, the sensitivity of the AI cell population ($x_2$) remains relatively low until after relapse, as with the CAS model. This underscores the necessity of maintaining data collection efforts each time a switch is made between on-treatment and off-treatment and then continuing well after the decision has been made to cease drug holidays. We also observed, through relative sensitivities based on simulated data, that the IAS model is much more sensitive to the parameter $r_0$, which controls the PSA level at which off-treatment times begin, than any other parameter in the model. Surprisingly, the model is least sensitive to the parameter $r_1$, which controls the PSA level at which on-treatment times begin. Thus, our model indicates that calibrating how we decide to begin a so-called drug holiday may be much more important than how we choose to end those drug holidays. We expect that data from a clinical trial can therefore be useful in estimating not only $r_0$, but also parameters influencing cancer cell proliferation, apoptosis, and mutation.

Unlike intermittent androgen suppression where the on and off treatment is governed by certain threshold levels of total PSA, treatment with periodic androgen suppression requires that the patient switch between on and off treatment at certain intervals defined by the period $T$. Specifically, the patient remains on-treatment for the first half of the period and off-treatment for the second half of the period. Here we solved the periodic equations and established conditions on $\kappa$ and $T$ that guarantee relapse. We derived expressions for the relapse time, $T_R$, for the two main values of $\kappa$ and computed the relative sensitivity of $T_R$ over a parameter region that is optimal for patient care (that is, large relapse time and small on-treatment time). Finally, we compared $T_R$ and on-treatment time $T_{on}$ (see Figures 13 and 15) as the main control parameters $r_0$ (IAS) and $T$ (PAS) vary. When we compare IAS and PAS at parameter values for which the difference between relapse time and on-treatment time is maximum, for $\kappa=0$ the PAS model produced a slightly better relapse time of 824 days versus 818 days (0.7% larger) while the IAS model produced a better on-treatment time of 362 days versus 413 days (12% smaller). However, if $\kappa=1-\beta_2/\alpha_2$, the IAS model produced better results in both instances: a relapse time of 1709 days versus 1643 days for PAS (4% larger) and an on-treatment time of 755 days versus 825 days for PAS (8% smaller). This indicates that IAS is generally better (although not by a large amount) when considering the optimal treatment protocol. Finally, one advantage of PAS when compared with IAS is that periodic androgen suppression does not require regular blood samples from the patient in order to monitor PSA and determine when treatment should be switched on to off and vice versa. To our knowledge, PAS is not currently used in clinical practice. Our observations on the modest differences in relapse times suggest that a scheduled periodic treatment regimen may warrant further investigation.

In this paper we derived our models using only the foundational model of Ideta, et al., from [9]. In this case, the AI cells are not reversible and the AD cells are independent of the AI cells. This allowed us to approximate relapse time in a relatively straightforward manner as the model equations are only forward coupled. Future work may include investigations of IAS and PAS models where the AI cell types are both reversible and irreversible. In this case, the AD cell proliferation will also depend on the AI cells during the off-treatment cycle as described in [8], and a more detailed analysis of the relapse time is necessary as the model equations are backward coupled.

---