Intermittent Hormone Therapy Models Analysis and Bayesian Model Comparison for Prostate Cancer

The prostate is an exocrine gland of the male reproductive system dependent on androgens (testosterone and dihydrotestosterone) for development and maintenance. First-line therapy for prostate cancer includes androgen deprivation therapy (ADT), depriving both the normal and malignant prostate cells of androgens required for proliferation and survival. A significant problem with continuous ADT at the maximum tolerable dose is the insurgence of cancer cell resistance. In recent years, intermittent ADT has been proposed as an alternative to continuous ADT, limiting toxicities and delaying time-to-progression. Several mathematical models with different biological resistance mechanisms have been considered to simulate intermittent ADT response dynamics. We present a comparison between 13 of these intermittent dynamical models and assess their ability to describe prostate-specific antigen (PSA) dynamics. The models are calibrated to longitudinal PSA data from the Canadian Prospective Phase II Trial of intermittent ADT for locally advanced prostate cancer. We perform Bayesian inference and model analysis over the models’ space of parameters on- and off-treatment to determine each model’s strength and weakness in describing the patient-specific PSA dynamics. Additionally, we carry out a classical Bayesian model comparison on the models’ evidence to determine the models with the highest likelihood to simulate the clinically observed dynamics. Our analysis identifies several models with critical abilities to disentangle between relapsing and not relapsing patients, together with parameter intervals where the critical points’ basin of attraction might be exploited for clinical purposes. Finally, within the Bayesian model comparison framework, we identify the most compelling models in the description of the clinical data.

shown that patients are responsive to multiple hormone therapy cycles, eventually delaying the androgen independence insurgence (Klotz et al. 1986;Larry Goldenberg et al. 1995;Bruchovsky et al. 2006).
We consider models of intermittent therapy due to clinical interest and solve the inference problem using longitudinal PSA data from the Canadian Prospective Phase II Trial of IAD for locally advanced prostate cancer. This work aims to present the first systematic comparative study of IAD models, emphasizing their ability to disentangle relapsing and not relapsing patients and compare the models in the Bayesian framework. The goal is to detect the single model (or the group of models) that best represent the information in the considered dataset and, therefore, if possible, the most promising biological frame representing them. A general and historical review of the available prostate cancer models can be found elsewhere (Phan et al. 2020).
In Sect. 2, we present the data included in our analysis. In Sect. 3, we introduce the statistical framework used to analyze the data. Section 4 presents an analysis of the models and their performance over the dataset utilizing the framework considered. Section 5 compares the performance of all the models, and Sect. 6 concludes and discusses the paper's findings and future developments.

Data Cohort
We consider data from the Canadian Prospective Phase II Trial of intermittent ADT for biochemically recurrent prostate cancer (Bruchovsky et al. 2006(Bruchovsky et al. , 2008. The total patient number is N pat 101. Their median pretreatment serum testosterone is 13.0 µg L −1 , ranging between 0.4 and 23.0 µg L −1 . Over a maximum of n 5 intermittent ADT cycles, a median of 35.1-36.0 weeks is spent on-treatment (depending on n), and 25.6-53.7 weeks (e.g., n 5 and n 1 respectively) are off-treatment during the 6-year study. An example of a PSA profile for an individual patient is shown in Fig. 1a. This patient responded to treatment during the first two treatment cycles (τ 1 and τ 2 ) and progressed in his third cycle of treatment (τ 3 ). The oscillatory dynamics demonstrate the effect of the intermittent treatment, with a decrease in PSA during treatment and an increase once treatment is turned off. Each data point is assigned with an error of 1 day in time (i.e., the time resolution of the dataset) and a maximal PSA error value e max assigned of e max 0.1 µg L −1 assumed from the literature (Borros 2009).
We set the minimal PSA detection threshold equal to 1 0.1 µg L −1 , i.e., any patient data below this threshold is set to 0.1 µg L −1 . Patients with a minimal perday fluctuation below 2.0 µg L −1 , i.e., a minimal per-day fluctuation of the KLK3 glycoprotein enzyme of PSA of a typical man (Morgentaler and Conners 2015), are excluded because such small fluctuations are considered natural and not pathological. To only consider PSA concentrations above Poisson noise, patients with less than 2 N pat (i.e., the sample shot/Poisson noise) data points are also excluded. These exclusions result in our analysis considering data from 89 (N pat 89) rather than 101 patients. The patients' distribution per number of data points used after the selection process is shown in Fig. 1b compared to the original distribution.

Data Interpretation
The PSA trend shown in Fig. 1a is based on the interplay between cellular populations, i.e., with a compartment modeling approach. An androgen-dependent set of N D ≥ 1 cell populations (each with a concentration n D,k n D,k (t), k 1, . . . , N D representing the compartment concentration [μgL −1 ], t ∈ R time [day]) is assumed to contribute to the oscillatory behavior of PSA. Additionally, a set of time-dependent androgen-independent cellular populations of N I ≥ 0 cell populations n I ,l n I ,l (t), l 0, . . . , N I , also contributes to the PSA profile, such that PSA concentration c PSA is given by c PSA f n D,k , n I ,l , where f ∈ C 0 is a function belonging to the class of continuous solution C 0 (not necessarily smooth) of a suitably designed ODEs system. Any further dependence on space, temperature, and pressure is generally neglected in the IAD models' compartment approach. Furthermore, f is often assumed to be a linear combination of the N D and N I compartments, e.g., c PSA k w k n D,k + l w l n I ,l for some weights w l and w k .
By assuming ADT to be highly effective in the first treatment interval τ 1 , we can set n D (t ∈ τ 1 ) ∼ 0 for some t as an initial condition (hereafter i.c.). This approach does not necessarily hold for τ i with i > 1: Generally, ∀i where c PSA ∼ 0, we can equally well assume this setting for the i.c. of the n D,k n D,k (t) equations. Equivalently we can assume that a non-holonomic (i.e., with inequalities) condition for the fitting procedure holds at the beginning of the patient time series n D (t) ≤ n I (t) for some t ∈ τ 1 . Furthermore, in most of the models that we accounted for, these considerations are articulated with the addition of a few extra equations that interpret, at a local or global level in the parameter space, the contribution to c PSA (t) by the androgen quota, cellular plasticity, staminal cells populations, or other model specificities.
Finally, we know from biological arguments that, under treatment, the models' equations are designed to permit, at least for c P S A , to tend asymptotically to the value c PSA 0. Any model that does not permit the phase state c PSA c PSA (t) to reach approximatively null values for any t, i.e., t R: c PSA (t) ∼ 0, would fail to reproduce the patients whose first treatment is always successful (see Fig. 1a). We elaborate more on this in Sect. 4. Therefore, it is worth investigating if the models allow for stationary equilibria outside the treatment intervals and then if any of the patient best fit values have fallen close to those equilibria (when they exist). This behavior would imply a stationary or recurrent solution for the dynamics and, therefore, a constrained PSA's evolution if this "basin of attraction" is achievable in a biological time of interest. We stress that this mathematical behavior does not imply that the patient can effectively reach the point of equilibrium on the biological timescale of interest or a plausible point regarding toxicity levels.

Bayesian Inference
The Bayesian regression approach stems from the concept of probability as a measure of the plausibility of a model given the truth of the information in the data presented above. First, we encode the prior state of knowledge about the parameters considered p { p 1 , p 2 , . . .} into a prior distribution function Pr(p|I ), where I represents any available information. Typically, this can be achieved with a flat, uniform, and not informative prior at the beginning or with a sharper prior when the model is better trained. We return to this point in Sect. 3.1. Secondly, we consider the dataset, D, through the likelihood L(p, D) Pr(D|p, I ). Finally, the inference problem is solved, studying the probability distribution function encoding the knowledge of the prior and the information encoded in the likelihood of the data Pr(p|D, I ) ∝ Pr(p|I )L(p, D).
Standard techniques to achieve this result are fully analytical (e.g., for some linear regression), approximated (e.g., asymptotic approximation, Laplacian approximation, Gaussian approximation, etc.), iterative (e.g., Levenberg-Marquardt), or fully numerical (e.g., simulated annealing genetic algorithms). The choice between these techniques depends on the nature of the problem. Here, we start using Laplacian approximation with hyperparameters (Hutter et al. 2011;Murphy 2012;Theodoridis 2015), as a few of the mathematical models that we consider herein are nested, to solve the inference problem (i.e., to search for the optimal set of parameters p that best represent the data). In order to confirm the inference results and to perform the Bayesian model comparison numerically, we test the results both against the nested sampling approach to the global likelihood (hereafter evidence) (Skilling 2004;Mukherjee et al. 2006;Feroz and Hobson 2008) and with the differential evolution search (Feoktistov 2006;Goode and Annin 2015) with up to aggressive scaling factors (≤ 0.9) and cross probabilities (≥ 0.1). For the Bayesian model comparison part of our work, see Sect. 5, the nested sampling-based approach will embed the results in a natural framework.
Finally, we note that substantial limitations in the fitting procedure came from the sparse and irregular temporal sampling in the clinical data. This irregularity impacts the parameter space exploration due to the lack of condition on the PSA trend's derivative.
The partial derivative ∂ t c P S A (p; t) is not smooth, thus inhibiting using some straightforward optimization techniques based on the PSA curves' gradients or convexity (Theodoridis 2015).

The Priors Pr p|I
While robust approximations or numerical tools have been adopted for the Bayesian framework, special attention is paid to the use of priors. As mentioned, Bayesian inference requires the use of the priors, Pr(p|I ), for parameter estimation. With initially unknown priors, we implement uniform priors over the parameters' full ranges (Fig. 2a). By requiring all model parameters to be positive, we can assume the Heaviside step function θ θ (p) as (unnormalized) prior, this approach is generally referred to as "improper prior" as it is unbounded above, it cannot be normalized and therefore it does not have a mean, standard deviation, median, or quantiles. We set an upper bound for each parameter to be p < p max with a max value p max < +∞ ∀p strictly. An alternative functional tested is the non-informative Jeffreys prior, Pr(p|I ) ∝ √ det(F(p)) with F symbol referring to the Fisher Information matrix (Jeffreys 1946) and "det" to the matrix determinant.
Extra than testing with flat/Jeffreys priors, in the numerical nested sampling approach, we explore the parameter space logarithmically to avoid divergences, and once we reach a statistically significative sample, i.e., above the Poisson noise fluctuation (∼ N pat ) shaping the posterior PDF, we proceed to implement the posterior as a prior for the patients analyzed in the dataset; finally, we reiterate by implementing a recursive determination of the prior (Fig. 2b, c). Further details can be found in Pasetto et al. (2021), where we discuss Bayesian analysis of retrospective data to guide clinical decisions.

Analysis of IADT Mathematical Models
Here we consider only IAD models due to current clinical interest. Each model is presented and justified in a biological and mathematical sense in the original paper where the model was first presented, and we refer the reader to that papers for detailed model derivations. Similarly, the sensitivity analysis of the model parameters is presented in each paper individually, and we elaborate on it here only where necessary. We will refer to the relapsing patient set as -set and not relapsing to as relapse ¬ -set.
We parametrize the individual IAD data with a patient-specific control function T ps defined as follows: T ps (t) n i 1 1 τ i (t), 0 < t ∈ [t min , t max ] 1 with t min and t max minimum and maximum patient-specific treatment under consideration (e.g., Fig. 1a) and t min generally after the first treatment drop; n ≥ 1 is the number of intervals τ i considered. τ i ⊆ ]t min , t max [∀i is referred to as the "i th treatment cycle," 1 We use compact set notation here, e.g., 0 < t ∈ [t min , t max ] means all the possible values of t, positive, between t min and t max , i.e., 0 < t min ≤ t ≤ t max . Open brackets will exclude the borders, e.g., soon after τ i ⊆ ]t min , t max [ means that the interval τ i , e.g., τ [a, b] is properly included between t min and t max , but the limits of t min and t max are excluded: t min < τ i < t max . This allows us to work with the domain of existence of the indicator functions but arbitrarily truncate it at t min or t max . and 1 τ i is the indicator function 2 for the interval τ i (defined as 1 τ i 1 for t ∈ τ i , 0 otherwise). For modeling purposes, the weights/errors, e i for each data i, have been assigned either uniformly e i cnst.∀i or with a linear decreasing relevance from the last PSA concentration c P S A peak, sayĉ P S A (e.g., in τ 3 of Fig. 1) with e i c PSAi −ĉ P S A t t , i.e. at t t, and e i ĉ PSA − t − t + |c PSAi | for t t. Finally, we performed sensitivity analysis on all the models included here. Comments on the technique adopted are technical and left to Supplement A.

Ideta et al. (2008)
The model by Jackson (2004) can be considered the continuous ADT model prototype. Its extension to IAD therapy of interest was presented by Ideta et al. (2008). In this model (hereafter, I08), the authors drop the dependence of Jackson's model on the spatial distribution, which is only of theoretical interest but not resolved in clinical PSA data. Model simulations predict that intermittent ADT can only prevent progression if normal androgen levels decrease the growth rate of AI cells, which may be biologically unlikely since AI cells have androgen receptors with increased sensitivity (Grossmann et al. 2001). We consider the I08 model in the following form: with initial conditions n D (t 0D ) n D0 , n I (t 0I ) n I 0 . 3 As previously mentioned, n D and n I are the androgen-dependent and -independent population number of cells (or concentration). γ i and δ i , i ∈ {D, I } are growth and apoptosis rates for AD and AI cells, given, respectively, by: (2) In Eq.
(2) γ Dmax and δ Dmax are the maximal AD proliferation and apoptosis rates, δ D A is a control parameter on the effect of low androgen levels on the AD apoptosis rate, k D Aγ 0 is the AD half-saturation rate, k D Aδ 0 is the AD apoptosis rate dependence on androgen. Finally, δ I A and γ I A 0 modulate hormonally patient failing death and growth. Mutation from AD to AI cells are allowed at a mutation rate: thus, the mutation rate decreases as the androgen (here normalized at its homeostatic level c A0 0) approaches its max value μ DImax . A decoupled ODE model of the serum androgen concentration under treatment c A is given by: with initial condition c A (t 0A ) c A0 0, where δ c A is the androgen clearance rate. Here T ps T ps (t) is the patient treatment-specific function as defined in Sect. 2.1. Finally, the PSA density concentration of interest to us, c P S A , is a linear combination with weight w i of the population densities

I08A in the Context of the Data
We noted that the system of equations (hereafter SoE) composed by Eqs. (1)-(4) decouples in the androgen concentration c A . The analysis of the system results in a line of infinite equilibria on the intersection of the plane n D 0 with the plane c A c A0 − c A0 T ps in the space of phase-state variables (n D , n I , c A ). Thus, c A c A0 off-treatment and c A 0 on-treatment. Standard linear stability analysis (Wiggins 2003) shows that the Jacobian of the system produces a null generalized eigenvalue λ 1 0, a negative one λ 2 −δ A , and a more complicate third generalized eigenvalue that takes, off-treatment, the elegant form:λ off . The sign of λ off 3 can be evaluated for the best-fit parameter values that result from the inference works of Sect. 3 in the patients' cohort considered here (Sect. 2), resulting in being always positive for all the patients. Therefore, the abovefound equilibria lines represent a 1D nonstable manifold, and further investigations (e.g., in the context of the central manifold theory) are not of additional interest to us.
We are indeed more interested to further exploit the characteristics of the present dataset in the context of this model by using the decoupled nature of the serum androgen concentration c A . All the patients are considered from their first cycle of treatment, starting with T ps (t) 1 for t ∈ τ 1 . Hence, we can emulate with a Heaviside step function T ps θ (−t) a cycle of treatment followed by the off-treatment period for a suitable cyclic interval (on-off, on-off, on-off, and so forth) around the off-treatment start, set at t 0. Within this approach, the general solution of Eq. (4) is algebraic and reads: This equation is monotonic on the two phases on/off-treatment because the deriva- , on-treatment nor for t ≥ 0, off-treatment. By splitting the treatment in on/off-time, we can always reverse the bilinear map off-treatment 4 for c A c A0 , and c A 0 and δ A 0. We exploit Eq. (6) to obtain the probability distribution function (PDF) of the orbits over all the sets of patients remapping each cycle over the phase space section (O, n D , n I ). We take advantage by the sharp c A passage from its homeostasis value c A0 to null and vice versa in conjunction with the bijection map just found. Figure 3a shows the c A profile for a representative patient. While time is a monotonic increasing function, the map we considering is one-to-one only over the treatment cycle T ps 1 and the off-cycle T ps 0, respectively, and in these two tracks we can write the SoE as n D n D (t(c A )) n D (c A ) and n I n I (c A ). As c A sharply switches from c A c A0 and c A 0, we can limit ourselves to a first-order solution of the SoE. After simple algebra, we arrive at the approximate solution of the SoE in the form: to the first order in c A (and where means asymptotic-to). As evident, the second equation remains close to its initial value n I 0 , while the first is perturbed away, suggesting that we can first sample the PDF of the dataset for fixed values in n I around n I 0 and then investigate the PDF as sampled from the best fit obtained by the patient in the trial with Eq. (7). The results are shown in Fig. 3b. The trend of the two distributions 4 Note how we could consider the resulting SoE as function of the variable c A to reach a fully algebraic solution of the system by taking the ration of . Nevertheless, it is more fruitful to look at the trend of Eq. (6) as obtained by the best fit procedure introduced in the next section.
for the development of resistance and continuing response patients is comparable as above the starting value n D n D0 , while the trend diverges for smaller values n D . Because we assume n D is a proxy for c P S A at small values of n I , as evicted from Eqs. (5) and (7), if the model correctly interprets the data, then a patient with an initial PSA-drop below 10% of its initial value is highly likely to be a continuous responder. The risk of resistance development grows to about 50% when the initial drop in PSA is around 30%.

I08B in the Context of the Data
For I08B, where δ I A γ I A , the ratio presented in Eq.
(3) evidences the structural non-identifiability of the SoE. The treatment of the equilibria and their stability is more straightforward in this model form than in I08A. The only equilibrium point is given by {n D , n I , c A } eq 0, 0, c A0 − c A0 T ps with generalized eigenvalues λ i of the Jacobian at the equilibrium given by λ 1 . Following the I08A assumptions, we investigate the model under the conditions δ Dmax > γ Dmax , δ D A > 1, and by requiring that μ max < γ I − δ I to avoid the annihilation of the populations. Under these conditions, we can prove that λ 1 ≤ 0 and λ 2 ≤ 0, and that for λ 3 it holds the same consideration as for I08A due to the non-stable nature of the resulting equilibrium manifold.
Analogous consideration on the non/relapsing treatment holds for I08B as for I08A, but with more straightforward treatment for I08B than for I08A: the two equilibria at homeostasis c A c A0 and at null androgen concentration, c A 0, attract the dynamics as for I08A and self-explain the orbit profiles. Therefore, the identical results from the inference of I08B on patients' trials can be obtained for the PDF but are not depicted again.

Eikenberry et al. (2010)
The model developed by Eikenberry et al. (2010, hereafter E10) was an attempt to describe the interaction between testosterone (T, the primary androgen in the serum), its enzyme 5a-reductase to dihydrotestosterone (DHT), and their binding (T:AR and DHT:AR) with the androgen receptors (AR) in the prostate. Because of model E10's versatility, we have included it in the IAD treatment model comparison. Of note, the authors have not proposed the model to fit data, and here we reinterpret E10 beyond the scope of the original paper. The modulation due to intermittent IAD is assumed in testosterone time modulation. While a linear relation might not be readily available from the literature between testosterone and PSA level (Elzanaty et al. 2017), we recode the testosterone concentration n T , in E10 as follows: which we couple with the original system of equations: with five nominals initial conditions: Here, the treatment function T ps modulates testosterone influx into the prostate-function ϒ(n S ) original in E10 and that we are going to adopt here, where n S is the testosterone serum concentration. Furthermore, we consider the androgen receptor concentration n R and the dihydrotestosterone concentration n D H T together with two quota concentrations q T :R and q D H T :R (Droop 1968), here, taken to be the T:AR complex and the DHT:AR complex concentration, respectively. γ R is the AR production rate, δ R is the AR degradation rate, δ T is the testosterone-specific degradation rate, and δ D H T is the dihydrotestosterone degradation rate. The mass-action constants for the androgendependent component (testosterone) and dihydrotestosterone binding the AR are , and the 5α reductase converts T to DHT by Michaelis-Menten enzyme kinetics with concentration n 5α , turnover number μ cat and constant k M 0.

The Model in the Context of the Data
If we set a ≡ μ cat n 5α − δ T k M , b ≡ 1 − T ps ϒ(n s ), and ≡ (a + b) 2 − 4bδ T k M , then two critical points can be isolated at the intersection of the nullclines hyperplanes of the phase space. On-treatment, the first point While only the second of these equilibria is of biological interest, it is not a stable equilibrium. Obtaining the complete set of generalized eigenvalues requires a cumbersome solution of three cubic equations, yet the check for the stability requires much less effort once we realize that one of the generalized eigenvalues from the characteristic equations and that it proves to be always positive for all the inference results in the trial patients.
Finally, we note how the model could represent an essential instrument for investigating the relapsing mechanism evidenced in some patients, which remains one of the goals of this work for its potential clinical implications. We identify three over five state variables by inspecting the phase-state space with a striking separation between and ¬ . Figure 4a shows the 3D probability distribution function of n T , n R , and q T :R . The density map of the temporal evolution of and ¬ sets clusters (over the orbital evolution spanned by the patients analyzed) on a well distinct area of the phase , n T and n R for a representative patient. The optimal fit c PSA dynamics (i.e., for optimal parameters p) and corresponding data are shown by the red curve and black dots with error bars, respectively; dashed green curves and the corresponding shadows show the sensitivity of c PSA when the parameters are increased/decreased by (Color figure online) space, splitting in the n T vs. n R space and at least partially in the orthogonal q T :R space.
In Fig. 4, panels b, c, and d, we exploited the Direct Differential Method (DDM) for sensitivity analysis to track the time dependence of the sensitivity S c PSA j ≡ computed at the best fit parameter valuesp, where p {n T0 , n R0 , q T:R0 }, respectively. As shown in Fig. 4b-d, a slight variation of the parameters does not dramatically affect the trend of c P S A . Thus, there is minimal sensitivity of c P S A to the parameters. This result shows that the PDF of the combination of parameters investigated might be an excellent tool to explore the origin of the resistance with the E10 model. The sensitivities were computed using the Direct Differential Method (DDM), as mentioned at the beginning of Sect. 4 and it is reported more in details in Supplement A. As evident from Fig. 4b-d, different parameters have different sensitivity on a different phase orbit with n T 0 more sensitive under treatment and n R or q T :R more sensitive out of treatment. DDM not only demonstrates the stability of the results obtained but also adds extra information on when a model is sensitive to a parameter change. This result is significant when dealing with models with varying behavior on-and off-treatment.

Hirata et al. (2010)
A series of studies (Tanaka et al. 2010;Hirata et al. 2012;Hirata and Aihara 2015) motivated the model by Hirata et al. 2010 (hereafter model H10) to capture intermittent ADT dynamics. The model is based on the coupled AD-AI population cells, supplemented with a population of irreversible AI cells, AI-Irr representing the first three-compartment model in the literature (Fig. 5a). Here we report the mathematical formulation in the proposed framework's formalism and refer to the original paper for a detailed model description. The SoE reads with our generalized notation: n D (t 0D ) n D0 , n I (t 0I ) n I 0 , n Irr (t 0Irr ) n Irr0 , where terms retain the identical biological meaning as previously described and the two irreversible, and reversible changes in the AI cell population are considered with the relative growth rate γ on/off i on-and off-treatment with i ∈ {D, I , Irr}. The serum concentration is computed as in Eq. (5) for i ∈ {D, I , Irr}.

The Model in the Context of the Data
During both on-and off-treatment cycles, nullclines analysis leads to {n D , n I , n Irr } eq {0, 0, 0} as the only equilibrium point. By setting a ≡ γ off D +γ off we can write the generalized eigenvalues of the Jacobian at the equilibrium as λ 1 γ off Irr + T ps γ on Irr − γ off Irr and the other two λ 2,3 in a compact form as λ 2,3 1 2 T ps (b − a) + a ± . This result implies that the equilibrium is stable on-treatment and unstable off-treatment.
The phase space shows that responsive and resistant patients cluster differently on the phase-state variables. Figure 5b shows that the probability density function for the best-fit patient groups around the initial value for n I ∼ n I 0 and n Irr ∼ 2.1n Irr0 . Thus, the irreversible component of the model offers a potential tool to disentangle patient responses from the model fitting. As the resistant patients are expected to increase their irreversible cell component (i.e., asymptotically n Irr n Irr0 with " " meaning asymptotic greater), we note that n I n I 0 in responsive patients. The model structure allows for the simulation of various PSA profiles thanks to the introduction of a new degree of freedom carried out with the third-compartment equations. Figure 5c shows the phase space plane for an example taken from the set of patients (Patient #33), while Fig. 5d shows the quality of the captured PSA concentration c P S A profile achieved by this model.

Portz et al. (2012)
The Portz et al. (2012) model is based on the cell quota concept (Droop 1968), which is modeled as: with q(t 0i ) q 0i for i ∈ {D, I }. The cell quota can grow to the maximum cell quota rate γ max and degrades at a constant rate δ q , with q max representing the shared max cell quota, v max the maximum cell quota uptake rate, q imin < q max the minimum cell quota for androgen, and 1 k q/2 > 0 the uptake rate half-saturation level (Packer et al. 2011). The authors allow mutation between both cell populations, from AD to AI and vice versa, at rates μ DI and μ ID given, respectively, by the Hill's equations of index m 2: where μ DImax is the maximum AD to AI mutation rate, μ IDmax is the maximum AI to AD mutation rate, and k m DI/2 and k m ID/2 are the cells mutation rate half-saturation level. The model follows the evolution of AD/AI cell populations, n D and n I , respectively, with the following equations: for q i (t) 0∀t and i.c. n D (t 0D ) n D0 and n I (t 0I ) n I 0 . The cell apoptosis and proliferation rates are, respectively, given by δ i and γ i for i {D, I }. The authors model the quota for both AD and AI cell populations independently. In general, we assume q I min < q Dmin to ensure that AI cells have a greater proliferation capacity in low androgen environments and n D (t 0D ) ∼ 0 with t 0D soon after treatment, as well as n I (t 0I ) ∼ 0 at t 0I at the beginning of the first treatment. Furthermore, a communal maximum proliferation rate γ max between the two populations is assumed. Both AD and AI cells produce PSA at a baseline rate γ PSA0 under the androgen dependence specified by: with c P S A (t 0P S A ) c P S A0 , and where k P S A,i/2 are the half-saturation rates and γ P S A,i the growth rates, for i {D, I }. Several variants of this quota model can be found in the literature. In the present work, we consider only a couple of them. A detailed comparison between  and (Portz et al. 2012) can be found elsewhere . The model's complexity is demonstrated with a tube plot (Fig. 6a, b).

P12A in the Context of the Data
The model is an extension of the models by Ideta et al.,shown in Sect. 4.1, where the equation of the quota decouples from the two cell populations behavior. Nevertheless, the quota evolution q q(t), common to n D and n I , is generally The quota profiles are represented as a cross section. The orbital evolution corresponds to the gray box time interval in panel a. The blue and green arrows represent the independent and dependent quota levels of the model, respectively. On the orbit section, the tube cross section has been computed considering q D along the normal N and q I along the binormal B and by solving the Frenet-Serret formulas (assuming the vector field of the system under consideration being along the tangent T ). c, d The distribution of the eigenvalues λ off 1 and λ off 2 off-treatment for models P12B and P12A, respectively (Color figure online) smoother than c A c A (t) in I08A or I08B hence not justifying the approximations worked out in those models. In P12A, the only equilibrium point is at {n D , n I } on/off eq {0, 0} and for the decoupled quota equation at q on eq γ Dmax q min γ Dmax +δ q and q off eq γ Dmax( k q/2 +1)q min (q max −q min )+q max v max
Nevertheless, we note that from the plot in Fig. 6c, how the best-fit solutions obtained from our inference work for all patients with this model falls in the area where λ i > 0, for both i 1 and i 2, i.e., we are never in the presence of an attractor (off-treatment). Therefore, patient dynamics never intercept an area of the parameters' space defined by the hyperplane that would (eventually asymptotically) lead to the annihilation of the n D and n I cell population, i.e., a steady state or a reduction of the disease present under the detection threshold. This plot is compared with the companion model in the next section, which simplifies P12A.

P12B in the Context of the Data
In this model, the authors extend the use of the quota concept to both n D and n I individually, i.e., fully exploiting Eq. (11), but retaining the same proliferation rate γ max . The large number of parameters required by the model makes the posterior maximization time-consuming and computationally expensive in the Bayesian framework, especially in a fully numerical nested sample approach (Skilling 2004) or using differential evolution optimization tool (Feoktistov 2006;Goode and Annin 2015). For this reason, a first inference approach has been performed within Laplace approximation and followed up at the patient-specific level where judged necessary. 5 As in P12A, the P12B critical points are {n D , n I , c P S A } eq {0, 0, 0} both onand off-treatment, while for the decoupled quota equation-stability points are found at q on i,eq γ max δ q +γ max q imin and q off i,eq a −1 μ max q imin k q/2 + 1 (q max − q imin ) + q max v max with a ≡ v max − k q/2 + 1 δ q + γ max (q imin − q max ) 0, δ q 0 and γ max 0 and for i ∈ I , D. As in P12A, the three generalized eigenvalues of the Jacobian at the equilibrium are always negative. Equations for the generalized eigenvalues λ on/off i for i 1, 2 along the quota directions are analytically available but slightly cumbersome; vice versa, more interesting is the plot of λ off i for i 1, 2 shown in Fig. 6d. The P12B solutions distribute a small number of patients in the P12A inaccessible area of double negative generalized eigenvalues (orange square in Fig. 6c, d). In this zone of the P12B parameter space off-treatment, the model predicts a constrained (or asymptotically constrainable) tumor cell population. Finally, we note how P12A is nested in P12B. Thus, P12B always obtains a better score in the same data representation but suffers from overfitting. We investigate this problem and offer a solution in Sect. 5 in the context of the Bayesian model comparison.

Morken et al. (2014)
In Morken et al. (2014), the authors extend model P12B by adding ADT-induced apoptosis of prostate cancer cells in addition to the inhibition of their growth and proliferation. Therefore, the model (hereafter M14) implements the per capita mortality of androgen-dependent and independent populations introduced in the previous section with the equation: where k i/2 for i ∈ {D, I } are the apoptosis and half-saturation levels for the dependent and independent populations, respectively. We consider the SoE in the form of: for q i (t) 0∀t and i.c. n D (t 0D ) n D0 and n I (t 0I ) n I 0 , together with the equivalent of Eq. (14): with i.c. c PSA (t 0PSA ) c PSA0 . Furthermore, the same notation as in the models by Portz et al. is followed and not repeated here.

The Model in the Context of the Data
The analytical treatment is analogous to P12B but enriched in the dynamics variety for the extra parameters introduced in Eq. (15), although without changing equilibrium points. 6 Our model analysis did not report other notable features.

Baez and Kuang (2016)
The model by Baez and Kuang (2016) presents a variant of the P12A model that is able to fit PSA and androgen dynamics, thus improving PSA trend forecasting. Two models are presented in the authors' work and considered here. The first (hereafter B16A) is a single population model of cellular concentration n, and two equations are coupled with it, for δ max the time-dependent (over a timescale τ δ max ) maximum baseline cell death rate and c P S A the PSA concentration that are modeled as: and a decoupled equation for androgen level: with n(t 0n ) n 0 , c P S A (t 0P S A ) c P S A0, δ max (t 0δmax ) δ 0max , and q t 0q q 0 > 0 strictly. The quota q 0∀t is produced at a rate γ γ 1 T ps + γ 2 .
In the same work, the authors also presented a two-populations model tracking both sensitive n D and independent n I cell evolution (hereafter B16B). By implementing their SoE within the approximation that all the cells have, on average, the same mass and density, we can recast their SoE in the form: for n i , i ∈ {D, I } never contemporaneously null, with initial conditions n D (t 0D ) n D0 , n I (t 0I ) n I 0 , q t 0q q 0 , and c PSA (t 0PSA ) c PSA0 . The maximum AD to AI mutation rate is given by μ DI max . Furthermore, because AI cells, n I , proliferate at lower androgen level it is assumed that q I min < q Dmin , and δ Dmax > δ I max because independent cells are less susceptible to apoptosis by androgen deprivation than sensitive cells.

B16A in the Context of the Data
The decoupled quota equation presents an equilibrium at q eq γ max (q min −q max ) γ max +(γ 2 +γ 1 T ps ) + q max when γ max + (γ 2 + γ 1 T ps ) 0, belonging to the positive hyper-quadrant of the phase space (i.e., it is of biological interest). The remaining set in Eq. (18) shows two equilibria at {n, δ max , c P S A } (1) eq 0, 0, γ PSA,0 q eq δ PSA , which are always in the positive quadrant of the phase space of interest and {n, δ max , c P S A } (2) eq γ max( q eq −q min) q eq δ , 0, δγ PSA,0 q eq +γ max γ PSA,1( q eq −q min) δ PSA δ with q eq 0, δ PSA 0 and δ 0, which is also biologically meaningful. By studying the generalized eigenvalues, we see that the first of the equilibrium presents three negative generalized eigenvalues, one of which is always positive (i.e., it is a saddle point); the second equilibrium point produces the eigenvalues λ −γ max − γ 2 + γ 1 T ps which are all always negative, thus representing a stable point of attraction.
Due to the stability of the second equilibrium (on-and off-treatment), it is worth investigating the proximity of the patients' orbits to the equilibria on the Poincare sections involving the PSA concentration c PSA we obtained from Eq. (18). Nevertheless, the low quality of the likelihood, L(p) Pr(D|p, I ), see Sect. 3) in the set of patients, demotivates further analysis. A single population n seems to not adequately capture disease progression, which remain the primary focus of our work, making the model less attractive for clinical implications and therefore not pushed forward here.

B16B in the Context of the Data
The model presents cubic dependence on q and quadratic on n D . 7 We select to investigate only the null-equilibrium point of independent and dependent cells. It is evident that n D 0 is an equilibrium for the first of Eq. (20). Therefore, by assuming n D 0 (and n I > 0 strictly), we can confirm the existence of two equilibria, the first located at {n I , q, c PSA } (1) eq 0, q max + γ max (q I min −q max ) γ max +(γ 2 +γ 1 T ps ) , γ PSA0 δ PSA q eq , for γ max + (γ 2 + γ 1 T ps ) 0 and δ PSA 0, which is of biological interest. The second, algebraically more cumbersome, reduces its nonnegativity condition to the simple one δ I max + γ γ max (q I min −q max ) γ max q I min +γ q max + γ γ max (q I min −q max ) ≤ 0, that is verified over all studied patients.
Again, as explored in previous models, we are interested in the existence of negative generalized eigenvalues of the Jacobian at the equilibria off-treatment, i.e., a point of equilibrium with an asymptotic constrained expansion of the tumoral cell population. Despite the model complexity, it is easy to prove numerically that the Jacobian for both equilibrium points has at least one positive generalized eigenvalue, making these points saddle points that are not of interest to us.

Elishmereni et al. (2016)
The Elishmereni et al. (Elishmereni et al. 2016) model accounts for two dynamics: disease dynamics represented by PSA used as a proxy for tumor volume and the pharmacology dynamics combined with the emergence of resistant cells from androgen receptor-independent n I and testosterone androgen receptor-dependent n I AR mechanism. The PSA concentration c P S A of interest to us, is governed by the following numerically highly complex SoE 8 :

The Model in the Context of the Data
The system has no equilibria influencing its dynamics, as evident from the 6th of Eq. (21). Further analysis is done in Sect. 5 to determine how well the model performs in the Bayesian model comparison.

Zhang et al. (2017)
Zhang et al. (2017) present a three-population competition model, based on Lotka-Volterra (LV) dynamics, where androgen-dependent n D , androgen producing n P , and androgen-independent cells n I , are considered. Basing the approach on game theory, the authors derive a competition matrix α α i j i, j ∈ {D, I , P} based on the parametrization of growth rates γ i and carrying capacities K i with i ∈ {D, I , P} resulting in this set of algebraic-differential equations: dn D dt γ D n D 1 − α 11 n D +α 12 n p +α 13 n I n p (β−Tps+1) , dn P dt γ P n p 1 − α 21 n D +α 22 n p +α 23 n I K P , dn I dt γ i n I 1 − α 31 n D +α 32 n p +α 33 n I K I , where ADT is modeled by the decreasing carrying capacity with β < 1 or supporting androgen-dependent cells with β > 1. The authors considered several constraints derived from the literature and researchers' experience to shape the model parameter influence: α ii 1∀i, α 31 > α 21 , α 32 > α 12 , α 13 > α 23 , α 13 > α 21 , α 32 > α 31 , and α i j ∈ ]0, 1[∀i j. Finally, the PSA dynamics is governed by: with δ the PSA clearance rate.

The Model in the Context of the Data
With the coupling of Eq. (24), the system presents four equilibria, but only two are of biological interest: {n D , n P , n I , c PSA } (1) where these ratios exist. For the first equilibrium, the eigenvalues of the Jacobian are positive in the n D , n P and c PSA phase space and therefore of marginal interest. Vice versa, by setting a ≡ 1 + β, b ≡ β − α 12 + 1, d ≡ α 21 (β − α 12 + 1) + 1 and e ≡ β + α 21 (β − α 12 + 1) 2 − α 12 + 1 together with the discriminant squared 2 (eγ D + βγ P + γ P ) 2 − 4adeγ D γ P , we can write the four eigenvalues of the Jacobian for the second equilibrium offtreatment as: λ , where the ratios exist, which are always negative for the fitted parameters, hence representing a stable equilibrium and opening the possibility to achieve an equilibrium off-treatment.

Phan et al. (2019)
The model (hereafter P19) presented by Phan et al. (Phan et al. 2019) is a variant of the work of Sect. 4.6 (Baez and Kuang 2016) in which the third population of weakly dependent cells, n wD , is added to investigate the influence of extra degrees of freedom added by the new population. The death term is also adapted from Eq. (16). Retaining the notation used in Sect. 4.6, we can recast the model in the following form: q + k I /2 − δ I n 2 I dq dt q − γ 2 + γ 1 T ps − γ max + γ max (q Dmin n D + q I min n I + q w Dmin n w D ) c PSA0 together with the required biological inequalities q Dmin > q w Dmin and q Dmin > q I min .

P19 in the Context of the Data
The idea of a third population is not new and already advanced with success in the model by Hirata et al. (2010). Nevertheless, the structure of the equations in Eq. (10) is very different from the Hirata et al. model in Eq. (10), with significantly more parameters not readily justifiable within the present dataset quality. Similar considerations were already worked out by Phan et al. We remark only that the complexity of the analysis, already evident in Sect. 4.6.2, is pushed further in this context, where only numerical investigation is available for equilibria and stability. The only off-treatment equilibrium accessible by the orbits is the one for n D , n w D, n I , q, c PSA off eq 0, 0, 0, γ max q I min +γ 2 q max γ 2 +γ max , γ PSA0 (γ max q I min +γ 2 q max ) δ PSA (γ 2 +γ max ) and δ PSA 0 which is always positive with always negative eigenvalues λ off 1 −γ 2 −γ max and γ off 2 −δ PSA . This is of limited biological interest as it is not compatible with the irreversibility nature of n I , if not by surgical castration.

Brady-Nicholls et al. (2020)
The Brady-Nicholls et al. (2020) model (hereafter B20) is based on the hypothesis that prostate cancer stem cells' enrichment induces resistance. The model correlates stem cell proliferation with serum PSA through SoE for the prostate cancer stem cells n S , the non-stem (differentiated) cells n D , and for PSA serum concentration c P S A . We report the system in the following way: with initial conditions n S (t 0S ) n S0 , n D (t 0D ) n D0 and c P S A (t 0P S A ) c P S A0 . It is assumed that stem cells divide at rate log(2), and the division is either symmetric yielding two stem cells (Enderling 2015) or asymmetric, where the stem cell produces one stem and one differentiated cell. The parameter that governs this effect is p s . The PSA differentiated cell production rate and PSA clearance rate are given by γ P S A and δ P S A , respectively, and T ps is the patient-specific treatment function (see Sect. 2.1).

The Model in the Context of the Data
The SoE presents an infinite set of equilibrium points when off-treatment T ps (t) 0 in the intersection of the plane n S (t) 0 with the plane given by c PSA (t) γ PSA n D (t) δ PSA conditional to n D 0 and δ P S A 0 and the generalized eigenvalues of the Jacobian results in a double-zero generalized eigenvalue λ 1 0, λ 2 0 and a third negative eigenvalue λ 3 −δ P S A . Standard center manifold computation (Wiggins 2003) shows slow-2D-manifold dynamics that can be integrated to prove that the equilibria are unstable, and therefore not of interest.

Bayesian Model Comparison
Maybe the most vital point of the Bayesian framework, and the reason for its increasing popularity, is its innate model comparison ability, based on logic as an instrument for selection. We exploit this feature here using the Bayesian factor to compare the different models in their ability to simulate the data. It should be noted that this framework innate penalizes models based on the number of parameters required. This phenomenon is sometimes referred to as the Occam's razor factor (Jefferys and Berger 1992).
Starting from the classical Bayesian theorem, the Bayes factor β i j for PSA model M i over the PSA model M j is computed as a ratio of the probabilities of the two models (the odd-ratio, O i j ) such that, because N m i 1 Pr(M i |D, I ) 1 (with N m number of models to compare) if we are interested in how a model, say M 1 , compares to the other models M j , we arrive at We implement Eq. (28) to compare one patient at a time in one model against all the other models individually. For example, we are going to implement the comparison between M 1 , and every other M 2 as Pr(M 2 |D, I ) 1 1+O −1 21 , and we proceed iteratively.
We first explore the Laplace approximation framework under the assumption of equally-prioritized models, i.e., assuming that no previous preference can be accorded to any of the PSA models considered. We can exploit the asymptotic approximation (Murphy 2012;Theodoridis 2015) to the global likelihood, i.e., the evidence of the i th model, Pr(D|M i ), writing Pr(D|M i ) ∫ dp Pr(p|M, I )L(p|I ) with F being the information matrix introduced in Sect. 3.1. A classical result of Bayesian analysis is to consider the limit of the previous expression, but for an increased number of data points (N p → ∞) and flat priors, i.e., to compute the popular BIC index against AIC (Akaike 1974;Schwarz 1978). As the number of patient data points is often limited N p ∞ and we make explicit use of priors, BIC or AIC indices are not justifiable for model comparison. Instead, we build up a model-ofmodel function (Pasetto et al. 2021) to encode prior information as soon as available. Furthermore, as introduced above, we verified the Laplace approximation with fully numerical integration based on nested sampling algorithms (Skilling 2004;Mukherjee et al. 2006;Feroz and Hobson 2008), i.e., a numerical technique designed explicitly to compute the global likelihood of models with different degrees of freedom.  problem to all the models. The simulated disease dynamics vary significantly between the different models, and discrepancies between different models and patient data may indicate likely or unlikely biological mechanisms driving individual patients' resistance.

Single Patient Comparison Results
Model evidence (Fig. 7b) demonstrates that no single model represents all patient data accurately, suggesting that different biology drive individual patients' responses or that no model correctly faces the PSA problem. It may also imply that the PSA dynamics alone may be insufficient to discriminate between the different biological models. For some patients, model selection identifies models with a higher probability than others, but selection varies on a per-patient basis. As a classical proof-of-concept of the Bayesian technology employed, we report for the best performing model, E16, for patient #60 the unnormalized posterior marginalized PDF for each parameter in Fig. 7c. The PDFs are almost unimodal (but not for all parameters, see next Sect. 6), suggesting that this model represents fairly the patient and that the Laplace approximation could be justified. The credible intervals for the log parameters are also plotted and superimposed to the x-axis.

Overall Model Selection
We calculate the Bayesian maximum a posteriori performance for all the patients for each model (Fig. 7d), resulting in the Elishmereni et al. (2016) model marginally performing better on most patients. This result does not surprise us, as it is a model designed on clinical necessities, i.e., it was crafted with careful handling of the medical treatment. Nevertheless, as mentioned before, in the case of model comparison on a patient-to-patient basis, we could not identify a model that performed statistically better than the others thus eventually indicating the correct biological mechanics governing PSA dynamics. Figure 7d shows that E16 is preferred only on 10% of the patients, and eight of the 13 models have scores above the 8%.

Conclusion and Discussion
This work considers several mathematical models (Table 1) to simulate PSA dynamics of prostate cancer response to IADT in a prospective clinical trial. We exploit Bayesian continuous and discrete inference to interpret the data and identify the model with the highest likelihood of simulating the clinically observed dynamics. Using the PSA biomarker and the comparison between the different models, we1) identify several models that can separate responding patients and patients that develop resistance to intermittent ADT through the model fitting, 2) performed the Bayesian model comparison and demonstrated that the model by Elishmereni et al. (2016) performed slightly better than the others, i.e., as a better representative of most patients in the trial. Nevertheless, as evidenced in the example of Fig. 7c, the marginalized posterior PDF is often not all optimally single-peaked, casting shadows in an attempt to use this model to solve forecast problems. While we have focused on the models' inference to evaluate the possible connection with their underpinned biology, we will explore the potentiality and limitation of the models' forecasting ability to predict clinical PSA trends in a follow-up paper (Pasetto et al. 2021, in preparation).
The models analyzed herein synonymously use longitudinal PSA data to infer biological mechanisms underlying the observed PSA dynamics. PSA alone limited the potentiality of the presented approach and did not identify a single dominant model. Further information is necessary to simulate accurately and ultimately predict patientspecific PSA trajectories and the corresponding biological drivers of resistance. PSA alone might not be a helpful biomarker due to several dominant environmental factors outside the models' scopes that influence its evolution under treatment. The use of PSA as a surrogate marker for prostate cancer burden is indeed controversial. Overexpression of the PCA3 gene obtained from the mRNA in urine samples is proposed to be more suited to monitoring the cancer evolution (Bussemakers et al. 1999, p. 3;Laxman et al. 2008;Neves et al. 2008;Hessels and Schalken 2009, p. 3;Borros 2009).
Two alternative directions might to improve our understanding of the PSA as a prostate cancer monitor biomarker. From one side, a deeper understanding of the connection between PSA and tumor burden throughout model investigation might present the opportunity for a new class of models. Recently, the role of immature blood vessels formed under angiogenesis cues has been investigated to decrypt the Table 1 Model compartment sketches with phase variables and parameters modeled in the inference process relation between an increased tumor burden contemporaneous to decreasing PSA concentration (Barnaby et al. 2021). Additionally, models that include both PSA and androgen concentrations might present some advantages in the future. The modest but significant evidence of the E16 model over the other models might indicate a more critical relevance of the dormancy whose biology and mathematics are worth undoubtedly deeper understanding.
Exploring PSA model probability distributions to disentangle responsive and resistant patient cohorts in a clinical setting could be investigated through cross-correlations with PCA3 biomarkers. Such cross-correlation would provide independent verification of the analytical findings herein that remain, for the moment, data-driven and, therefore, entirely dependent on the one dataset utilized in all discussed models.
Alternatively, PSA could be a perfect biomarker, but inter-patient heterogeneity in resistance mechanisms may disallow identifying a single model for all patients. Additionally, different resistance mechanisms may evolve in an individual patient, with their respective contribution to the observed response dynamics changing during therapy. More complex models and dynamic adaptive weighting of different variables, terms, and parameters may be necessary. Such models, however, would be non-identifiable with the presently available data. A close dialogue between biologists, statisticians, and mathematical and genitourinary oncologists may help identify which data should be collected in future clinical studies to help detangle the complex prostate cancer response dynamics to intermittent ADT. Initial condition (i.c.) are not reported, but phase variables are. α {k DI/2 , k I D/2 , k q/2 , q max , v max , δ D , δ q , μ DImax , μ I D max }. Simbols used for the compartments are self-explanatory While the Bayesian framework is an invaluable tool to estimate model parameters and fit model dynamics to clinical measurements, the goodness of a fit informs neither the reliability of the estimated parameters nor the likelihood of a model representing the data chosen for the valid biological reason. Relatively invariant PSA profiles can be obtained for a significant range in each parameter, as it is the case of a weakly sensitive-highly non-identifiable parameter. This fact is often omitted in the modeling literature, where the results are often presented without structural or practical identifiability analysis. Many of the herein discussed models have not demonstrated structural identifiability, hence jeopardizing the attempt to claim the inference's practical identifiability herein. Nevertheless, we stress that a model's value may also be found in its interpretative role (Enderling and Wolkenhauer 2021). The complexity of the mechanism involved in the biological responses to intermittent ADT can be captured correctly for a single patient but fail for others. Therefore, the model comparison is not intended to provide an absolute ranking; instead, it provides an instrument to explore the different biological mechanisms implemented in mathematical models in clinically observed treatment response and progression dynamics.

Supplement
We have performed a sensitivity analysis for all the models included in the paper. However, as this analysis overlaps with the original papers' work, we do not include those results here. Our sensitivity is motivated by: (1) to understand the dependence of our results on the parameters. For example, if we can claim the possibility to split between relapsing ( ) and not to relapse (¬ ) patients by exploiting some specific model parameter combination, then the robustness of our result worth be investigated on the same parameter sensitivity to assign it the correct relevance and to evaluate its possibility to be applied to clinical tumor forecasting.
(2) The technique implemented for the sensitivity analysis investigates the parameter's sensitivity and the best-fit orbital integration, i.e., over the available longitudinal data. This approach enhances our understanding of when a particular /¬ segregating technique is more useful during or off-treatment, with consequent indications on the role that a model splitting potentiality might or might not have (and when) on a per-patient base. (3) Continuous but not differentiable functions might need particular attention in the computation of the sensitivities because of their definition as the Jacobian matrix's function. This approach represents a current research field often omitted in the mathematical oncology literature and is worth being brought to light.
Therefore, in what follows, we exploit the direct differential method (DDM) for sensitivity analysis (Gu and Wang 2013) to track the time dependence of the sensitivity S i j ∂ x i (t, p) ∂ p j , where in general it is x i c PSA and p j the generic parameter dependent on the particular model in the exam. For a generic vector field ∂x(p;t) ∂t f (x, p) with x(t 0 , p) x 0 we couple the integration of the SoE defining the model with: Generalized sensitivity (Stechlinski et al. 2018), based on the concept of generalized derivative for non-smooth c PSA profiles (Clarke 1990) and used because of the loss of differentiability at the bifurcation points T ps {0, 1} on the treatment parameter, has also been considered. We will not report the DDM analysis if not relevant to strength our specific results and we refer to the original model paper for general sensitivity analysis of the presented models.