{'fraction_non_alphanumeric': 0.04983828981620257, 'fraction_numerical': 0.020967105782125108, 'mean_word_length': 4.256758998175485, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal control of cytotoxic and antiangiogenic therapies on prostate cancer growth

Prostate cancer can be lethal in advanced stages, for which chemotherapy may become the only viable therapeutic option. While there is no clear clinical management strategy fitting all patients, cytotoxic chemotherapy with docetaxel is currently regarded as the gold standard. However, tumors may regain activity after treatment conclusion and become resistant to docetaxel. This situation calls for new delivery strategies and drug compounds enabling an improved therapeutic outcome. Combination of docetaxel with antiangiogenic therapy has been considered a promising strategy. Bevacizumab is the most common antiangiogenic drug, but clinical studies have not revealed a clear benefit from its combination with docetaxel. Here, we capitalize on our prior work on mathematical modeling of prostate cancer growth subjected to combined cytotoxic and antiangiogenic therapies, and propose an optimal control framework to robustly compute the drug-independent cytotoxic and antiangiogenic effects enabling an optimal therapeutic control of tumor dynamics. We describe the formulation of the optimal control problem, for which we prove the existence of at least a solution and determine the necessary first order optimality conditions. We then present numerical algorithms based on isogeometric analysis to run a preliminary simulation study over a single cycle of combined therapy. Our results suggest that only cytotoxic chemotherapy is required to optimize therapeutic performance and we show that our framework can produce superior solutions to combined therapy with docetaxel and bevacizumab. We also illustrate how the optimal drug-na\"{i}ve cytotoxic effects computed in these simulations may be successfully leveraged to guide drug production and delivery strategies by running a nonlinear least-square fit of protocols involving docetaxel and a new design drug.


Introduction
Prostate cancer (PCa) is the second most common cancer among men worldwide [21]. Current estimates indicate that the annual number of new cases globally exceeds 1,000,000. PCa is also responsible for over 300,000 worldwide deaths yearly and a major economic burden for our health care systems.
In most cases, PCa is an adenocarcinoma -a form of cancer that originates and develops in epithelial tissues. PCa is easier to treat in its early stages, especially before it spreads out of the prostate [96]. Although the disease rarely produces any symptoms before the tumor has advanced significantly, current screening protocols have helped identify PCa patients early [68,96]. The most common screening method for PCa is the Prostate Specific Antigen (PSA) test -a blood test that measures the serum level of gamma-seminoprotein which is a biomarker of the prostate activity. PSAbased screening is controversial in the urology community, primarily because it is prostate-specific, but not cancer-specific [96]. Despite the success of screening protocols, a significant number of patients eventually develops advanced PCa [9,96]. These include men who are diagnosed late, and patients who initially receive radical treatment (e.g., surgery or radiation) or follow a conservative management protocol (e.g., active surveillance), but experience cancer progression.
The most common treatments for advanced PCa are hormonal therapy and chemotherapy, but the latter may become the only remaining alternative after PCa becomes insensitive to hormonal strategies [68,96]. Here, we will focus on chemotherapy, which normally relies on cytotoxic drugs. These compounds aim at interrupting cell division and promoting cell death, hence contributing to reduce tumor progression. Cytotoxic chemotherapy for advanced PCa usually relies on the periodic delivery of docetaxel every three weeks [19,68,68,83,96]. However, treated tumors may resume growth months after the conclusion of therapy and develop resistance to docetaxel, hence requiring the use of alternative cytotoxic drugs for a second round of treatment. Angiogenesis has also been regarded as a promising therapeutic target for advanced PCa because it is known to play a pivotal role in the progression of advanced PCa and higher microvascular density correlates with poorer prognosis [65,70,95]. Antiangiogenic drugs block the formation of new blood vessels cutting off the supply of critical nutrients to the tumor. Bevacizumab is an extensively studied antiangiogenic drug that has also being investigated for the treatment of advanced PCa either as monotherapy or in combination with cytotoxic chemotherapy [2,71,83,85]. A clinical study showed that castration-resistant PCa patients did not benefit from treatment with bevacizumab alone [78], but another study showed that 75% of the patients treated with bevacizumab experienced a 50% PSA decline when taking bevacizumab in combination with cytotoxic chemotherapy [76]. However, multiple experimental and clinical studies have shown contradictory evidence about the effectiveness of antiangiogenic therapy for advanced PCa either as monotherapy or combined with a cytotoxic agent [2,48,71,74,83,[85][86][87]89].
Thus, developing optimal dosing strategies for chemotherapeutic drugs is an open problem of major medical interest, but it is particularly difficult when two or more drugs are combined because their joint effect may produce unexpected outcomes and toxicities. Recently, the use of mathematical models to describe the growth and treatment of cancer has been increasing our understanding of these pathologies and providing a means to monitor and forecast tumor growth via patient-specific computer simulations, which can assist physicians in clinical decision-making [1,5,10,18,72,79,100,101]. Several studies have focused on the effects of chemotherapeutic strategies for various types of tumors [8,24,34,39,46,50,77,101]. In this context, optimal control theory provides a robust framework to obtain optimal drug protocols according to a given mathematical model of tumor growth and treatment [6,12,15,17,25,45,52,53,73,94]. We argue that these quantitative approaches can play a significant role in the design of optimal drug regimens, providing the best dosing, timing, and even drug pharmacodynamics to treat each patient's tumor [39,45,63,101]. Hence, optimal drug protocols can help us reduce the total amount of drug needed to achieve a target tumor size reduction, which is very important considering the harsh side effects of cytotoxic and antiangiogenic drugs [2,19,68,71,83,83,85]. Because antiangiogenic and cytotoxic drugs attack the tumor using different mechanisms, mathematical modeling can also help us determine the most effective drug dose combination. Additionally, optimal pharmacodynamics derived from optimal control problems could ultimately be exploited to guide the design of new drug compounds [38,40,63,84].
Here, we propose an optimal control framework that relies on a recently proposed model of PCa growth accounting for the effect of cytotoxic and antiangiogenic drugs [16]. Our approach aims at finding the cytotoxic and antiangiogenic drug effects minimizing an array of common quantities of interest in clinical and experimental settings, such as tumor volume and serum PSA. We demonstrate the existence of unique solutions and derive the optimality conditions for our optimal control problem. Then, we propose a set of numerical algorithms to solve this optimal control problem and we use them to run a simulation study. Our computational methods rely on isogeometric analysis (IGA), a recent generalization of the finite element method with superior approximation properties granted by the use of splines as basis functions [37]. Our simulation study results in optimal therapeutic solutions that only feature cytotoxic effects, and hence suggests that antiangiogenic therapies may not be optimal for treating PCa. Because docetaxel and bevacizumab are reference drugs in therapeutic investigations of advanced PCa, we used their properties to build the initial guesses and admissible space of controls. However, our optimal control problem is drug-independent, i.e., it provides optimal cytotoxic and antiangiogenic drug effects without accounting for any specific drug delivery plan. Thus, in this work, we also capitalize on this advantageous property and propose a two-stage methodology to derive optimal drug compounds and delivery plans matching the optimal drug effects targeting a specific tumor: first, we calculate the optimal drug-naïve optimal effects via computer simulation of our optimal control problem, and, second, we run nonlinear least-square fits of an array of elective drug delivery plans and compounds. Our results show that our approach can render interesting optimal drug protocols matching the optimal drug effects to treat the target tumor.
The paper is organized as follows. Section 2 briefly describes our mathematical model for PCa growth featuring cytotoxic and antiangiogenic drug effects, as well as the formulation of the optimal control problem. Sections 3 and 4 present the existence results and the optimality conditions for the optimal control problem, respectively. Section 5 introduces the numerical algorithms that we propose in order to solve it, and we use these computational methods in a simulation study presented in Section 6. Then, we use the resulting optimal drug effects to illustrate how they can guide drug design and the conception of new drug delivery strategies in Section 7. Finally, in Section 8 we discuss the work presented herein, draw conclusions, and address future avenues of research. We leverage a recently proposed mathematical model of PCa growth [16] that includes the effect of cytotoxic and antiangiogenic drugs. Reference [16] includes a detailed description of the biological mechanisms considered in the model, and a biological interpretation of its parameters. Let Ω ⊂ R d be an open bounded set and d the number of spatial dimensions. Let φ 0 , σ 0 , and p 0 be sufficiently smooth given functions defined in Ω. The set Ω has a sufficiently smooth boundary ∂Ω and T is a finite time. By defining the space/time sets Q T := (0, T ) × Ω and Σ T := (0, T ) × ∂Ω, the model can be formulated as where φ is a phase-field that identifies the spatial location of the tumor; σ is the concentration of a vital nutrient; and p is the tissue PSA concentration. The serum PSA P s commonly used in clinical practice can be obtained as P s = Ω p dx [16,61]. The differential operators in the equations are defined as follows: the subscript t indicates partial differentiation with respect to time; ∆ denotes the Laplace operator with respect to the space variables; F denotes the derivative of F and ∂ ∂ν represents the outward normal derivative to ∂Ω. In Eq. (2.1), where M is a positive constant. We will also use the notation The model parameters λ, η, D, γ p , γ c , γ h , α c , α h , S c and S h are positive constants with biological significance; see Ref. [16]. For the solution of the initial boundary value problem (2.1)-(2.5), the functions U (t, x) and S(t, x) are considered to be given. U (t, x) models the inhibiting effect of a cytotoxic drug on tumor dynamics, while S(t, x) represents the reduction of nutrient supply produced by an antiangiogenic drug. The function m(σ) is defined as where m ref is a positive constant, while ρ and A are constants that determine the rate of cell proliferation and apoptosis. We define ρ = K ρ /K ρ and A = −K A /K A . Here, K ρ and K A are, respectively, the proliferation and apoptosis rates of tumor cells while K ρ and K A are scaling positive constants. The positive constants σ l and σ r are, respectively, a reference and a threshold value for the nutrient concentration; see Ref. [16].

Optimal control problem
Our optimal control problem consists in finding the functions U (t, x) and S(t, x) that provide the optimal cytotoxic and antiangiogenic effects to treat a certain PCa case. Therefore, our optimal control problem can be stated as and The quantities U max and S max are two threshold positive values such that S max ≤ S c . In Eq. (2.8), r + represents the positive part of r; the functions φ Q , φ Ω are prescribed targets for the tumor phase field in Q T and in Ω at the final time, respectively, while p Ω represents an upper target function for the spatial mean value of the tissue PSA (i.e., an upper target for the serum PSA). We assume that Notice that we have not imposed any specific drug delivery plan or pharmacodynamic properties for U (t, x) and S(t, x) in our mathematical model and the optimal control problem herein, hence making our approach drug-independent and solely relying on the type of therapy itself. In this context, U ad aims at defining feasible limits for the therapeutic effects that can be realistically achieved with common drugs. The optimal choice of delivery dates, doses and even pharmacodynamic properties can be made at a second stage, aiming at reproducing the optimal drug-naïve effects obtained from the optimal control problem described in this section (see Section 7).
The coefficients k i , i = 1, ..., 7, in Eq. (2.8) are nonnegative constants, such that there exists at least one i with k i > 0. In the numerical simulations presented in Section 6, we will focus on particular examples of the functional J, where only a few of the k i 's are simultaneously nonzero. For example, the first three terms in J, weighted by k 1 , k 2 , and k 3 , respectively, provide a slightly different control of the tumor extent. In our numerical simulations, not all three constants k 1 , k 2 , and k 3 will be nonzero simultaneously. The same argument can be applied to the constants k 4 and k 5 , which provide control of the PSA. The terms weighted by k 6 and k 7 control the relative participation of cytotoxic and antiangiogenic effects, which may ultimately control the total amount of drugs delivered (see Ref. [16] and Sections 6-7). The control of the total amount of drug may be especially important because some therapeutic compounds have very strong side effects. However, specific dose-dependent toxicities can be accounted for during the design of optimal drug protocols aiming at reproducing the optimal drug-naïve effects obtained from our optimal control problem (see Sections 6, 7, and 8).

Existence results
The functional framework involves the space H = L 2 (Ω), which is identified with its dual space H ≡ H, and the Sobolev spaces H 1 (Ω) and H 1 0 (Ω), the last one containing the elements of H 1 (Ω) with null trace on the boundary ∂Ω. We set as well as with the dense and compact injections W 0 ⊂ V 0 ⊂ H ⊂ V 0 and W ⊂ V ⊂ H ⊂ V . For the norms in these spaces we will use the notation · B , where B is the space we are considering. For r ∈ [1, ∞] and z ∈ L r (Ω) or z ∈ L r (Q T ) we simply denote the norms of z by z r . We will also make use of spaces of functions that depend on time with values in a Banach space B. Namely, for r ∈ [1, ∞] we consider the space L r (0, T ; B) of measurable functions t → z(t) such that t → z(t) r B is integrable on (0, T ) (or essentially bounded if r = ∞) and the space C([0, T ]; B) of continuous functions from [0, T ] to B. Perhaps, it is important to note that L 2 (0, T ; H) is a space completely isomorphic to L 2 (Q T ). Moreover, for r ∈ [1, ∞], W 1,r (0, T ; B) will denote the space of functions t → z(t) such that both z and its (weak) derivative z t belong to L r (0, T ; B). We point out that Hypotheses. According to the explanation already given, we emphasize that Let us introduce the spaces and (ϕ, σ, p)(0) = (ϕ 0 , σ 0 , p 0 ) a.e. in Ω. (3.11) Then, system (2.1)-(2.5) has a unique solution (φ, σ, p), with which satisfies the estimate x ∈ Ω, (3.14) 15) and (σ, p) fulfills the estimate Moreover, the solution is continuous with respect to the data, that is, for two solutions for all t ∈ [0, T ]. Finally, if (φ 0 , σ 0 , p 0 ) ∈ V 0 × V × V, the solution has the supplementary regularity (φ, σ, p) ∈ X 0 × X × X , and satisfies the estimate The proof of this result can be found in Ref. [16].
Proof. It is easy to check that the functional J is nonnegative, so that it has an infimum d ≥ 0.
Using the boundedness and the Lipschitz continuity of F and h in the interval [0, 1], we derive that

By (3.2) and (3.3) we obtain
Then, we can pass to the limit in the variational formulations (3.6)-(3.8) for the solution (φ n , σ n , p n ) and deduce that (φ * , σ * , p * ) is the solution to (2.1)-(2.5) corresponding to (U * , S * ). From (3.13) we see that φ n (T ), p n (T ) are uniformly bounded in H, so that φ n (T ) → φ * (T ) and p n (T ) → p * (T ) weakly in H by (3.21) and (3.23), and consequently (3.24) Also, with the help of (3.20) and the Lipschitz continuity of the positive part, it turns out that We aim at showing that J(U * , S * ) ≤ lim inf n→∞ J(U n , S n ) = lim n→∞ J(U n , S n ) = d by (3.19), whence J(U * , S * ) = d and so (U * , S * ) is an optimal control in (P ). Indeed, the second and the last two terms in J(U n , S n ) are weakly lower semicontinuous and On the other hand, we point out that the other four therms in J(U n , S n ) converge, as previously specified. This concludes the proof.
Proof. We note that equations (4.1)-(4.3) are linear parabolic and have the coefficients of Y, Z and P in L ∞ (Q T ) and the source terms in L 2 (0, T ; H), according to (3.1)-(3.3) and to the property 0 ≤ φ * ≤ 1 a.e. in Q T . The initial data are zero, then smooth. Therefore, by the general results concerning the existence and uniqueness of the solutions to linear parabolic systems (cf. Ref. [57]), we deduce that there is a unique triplet (Y, Z, P ) ∈ X 0 × X × X solving (4.1)-(4.3) and satisfying estimate (4.6).
Theorem 4.2. The following convergence properties hold and so they show that (4.1)-(4.5) is precisely the system of first order variations related to (2.1)-(2.5).
Proof. We write the system satisfied by (Y µ , Z µ , P µ ) : where Here, φ µ int , σ µ int and φ int,µ are measurable functions (see Appendix A) occurring in the Taylor ex- with similar expressions for σ µ int and φ int,µ . Since U µ → U * and S µ → S * strongly in L 2 (Q T ) as µ → 0, then by estimate (3.17) we see that We test (4.10) by Y µ , (4.11) by Z µ , (4.12) by P µ and sum up all of them, obtaining We see that a 1 , a 3 , b 1 , b 3 are uniformly bounded in L ∞ (Q T ) and, by (4.8), we have that |u| ≤ U max , |s| ≤ S max a.e in Q T . Then, by the Young inequality, it follows that Next, by applying the Hölder inequality and the continuous embedding V 0 ⊂ L 4 (Ω), we easily infer that where we used (4.6). The following remaining terms on the right-hand side of (4.22) can be treated in the same way, in order to find out that Going back to (4.22), we obtain by using the Lipschitz continuity of the function h . An application of the Gronwall lemma leads to Recalling (4.20), (4.21) and the boundedness and Lipschitz continuity of F , h , h in [0, 1] and of m in R, it is not difficult to check that a 2 , a 4 , b 2 and the last term in (4.24) tend to zero strongly in L 2 (Q T ), as µ → 0. Thanks to (4.24) this implies that Therefore, we see that where the limit is understood in the sense of the strong convergence in C([0, T ]; H) ∩ L 2 (0, T ; V ). This concludes the proof.
At this point, we can introduce the adjoint system in terms of the adjoint variables w, z, q : where k i , i = 1, ..., 5, and φ Q , φ Ω , p Ω are the coefficients and the functions involved in the cost functional (2.8).
Proof. Using the transformation t → T −t, we can rewrite (4.25)-(4.29) as a linear parabolic system with initial conditions in H and with the coefficients of w, z and q in L ∞ (Q T ). Moreover, the source terms are in L 2 (0, T ; H), according to (3.1)-(3.3) and 0 ≤ φ * ≤ 1 a.e. in Q T . Therefore, by the general results concerning the existence and uniqueness of the solutions to linear parabolic systems (cf. Ref. [57]), we deduce that there is a unique triplet (w, z, q) ∈ X 0 × X × X solving (4.25)-(4.29) and satisfying estimate (4.30).
We introduce the convex closed subsets of L 2 (Q T ): We also specify that the normal cone at U * to K 1 is given by Similarly, we denote by N K 2 (S * ) the normal cone at S * to K 2 .
Theorem 4.4. Let (U * , S * ) ∈ U ad be an optimal control for (P ) with the corresponding state (φ * , σ * , p * ). Then, the first order optimality conditions hold true, w and z being the first and third components of the solution (w, z, q) ∈ X 0 × X × X to the system (4.25)-(4.29).
Proof. Since (U * , S * ) ∈ U ad is an optimal control for (P ), then we have that where U µ and S µ were defined in (4.7), namely, Replacing the expressions of U µ and S µ in (4.33), by a straightforward calculation we obtain the optimality conditions in terms of the solution (Y, Z, P ) of the system in variations (4.1)-(4.5) corresponding to (u, s). We easily deduce that Now we test (4.1) by w, (4.2) by z, (4.3) by q, sum up and integrate by parts taking into account the boundary conditions (4.4) and (4.28). Using also the initial conditions (4.5) of the system in variations and the final conditions (4.29) of the adjoint system, we obtain the equation Using again the adjoint system (4.25)-(4.27) in (4.35) we obtain By comparison between (4.34) and (4.36) we easily deduce for all (ū,s) ∈ U ad . Then, takings = S * andū arbitrary in (4.37) we obtain (4.31). Finally, settinḡ u = U * ands arbitrary we derive (4.32), as claimed.
Remark 4.5. We note that the optimality conditions do not depend on the component q of the solution to the adjoint system. This is essentially due to the biological meaning of the model, because p, the third component of the state system, does not influence the evolution of the other components φ and σ, but is a result of their evolution.
We also recall that if Hence, if the coefficients k 6 and k 7 are positive, then (4.31) and (4.32) entail that The relations (4.40) can be written as and On the other hand, if a coefficient is zero, let us say, if k 6 = 0, (4.31) becomes A similar result can be deduced for S * if k 7 = 0.
Remark 4.7. We discuss here how the previous results change in the case that the controllers U and S are considered time dependent only. This formulation might be more realistic for the therapy and, hence we will use it to explore the behavior of the optimal control problem in the simulations presented in Section 6. First, the functional J becomes and the minimization problem is and by remaking the calculation from (4.36) we get Remark 4.8. For a later use we also recall that from a condition of optimality written as we can formally read that where ∇ U J(U * , S * ) and ∇ S J(U * , S * ) are the derivatives of J with respect to the first and second variables, respectively, calculated at (U * , S * ). Consequently, we have 5 Numerical method

Spatial discretization
We use Isogeometric Analysis [37], a recent generalization of Finite Element Analysis, to discretize in space the forward and adjoint problems. Our spatial discretization of the forward problem is based on the following weak form of Eqns.
We also introduce a weak form of the adjoint problem defined by Eqns. (4.25)-(4.27), which is stated as: find w ∈ V 0 , z ∈ V , and q ∈ V such that where, The above-defined weak forms are discretized by defining finite-dimensional spaces V h ⊂ V and V h 0 ⊂ V 0 . We construct the discrete spaces using C 1 -continuous quadratic B-splines [37]. The space V h is defined as V h = span{N A } A=1,...,n b , where n b = dim(V h ) and the N A 's are multivariate splines. We use the superscript h to denote finite-dimensional approximations to the exact solution of the forward and adjoint problems.
where the φ A 's, called control variables, are the unknowns of the problem. The functions σ h , p h , w h , z h , and q h are defined analogously. The functions that belong to V h 0 will have some control variables constrained to ensure that Dirichlet boundary conditions are satisfied. The Neumann boundary conditions relevant to this problem are naturally enforced within the weak form.

Time integration
Our time integration scheme is based on the generalized-α algorithm [13,44]. Let us call Φ the global vector of degrees of freedom associated with the unknown φ h , i.e., Φ = {φ A } A=1,...,n b . Likewise, we will also make use of the vectors Σ = {σ A } A=1,...,n b and P = {p A } A=1,...,n b . We introduce the residual vector for the forward problem We call U n = {Φ n , Σ n , P n } the time-discrete approximation to the control variables of the forward problem at time t n . In the forward problem, we calculate U n+1 from U n by enforcing the equation where ) Here, t n+1 − t n = ∆t n > 0 is the time step, and α m , α f , and γ are real-valued parameters that define the accuracy and stability of the algorithm. For the adjoint problem, we define the global vectors of degrees of freedom W , Z, and Q, corresponding, respectively, to the discrete functions w h , z h , and q h . The residual vector for the adjoint problem is Let us call Y n = {W n , Z n , Q n } the global vector of degrees of freedom for the adjoint problem at t = t n . The adjoint problem is solved backwards in time starting from data at t = T , so for this problem ∆t n = t n − t n+1 < 0. The equation that allows us to determine Y n from Y n+1 is As shown in [44], the generalized-α algorithm can be made A-stable and second-order accurate by taking ρ ∞ ∈ [0, 1] and All the calculations presented in this paper were performed taking ρ ∞ = 1/2 and using Eq. (5.29).
We linearized the nonlinear algebraic equations defined by (5.17)-(5.20) and (5.25)-(5.28) by using the Newton-Raphson algorithm. The convergence criterion to advance from one time step to the next one was that the individual residuals (e.g., R φ , R σ , and R p for the forward problem) are reduced to N L of its initial value. The linear systems that result after linearization are solved using GMRES [80] with diagonal preconditioner up to a predefined tolerance L or a maximum number of iterations.

Optimal control algorithm
We solve numerically the optimal control problem using the steepest-descent gradient method; see, e.g., Ref. [3] (Algorithm 2.2). As noted in Remark 4.7, we will consider the controls U and S to be exclusively time-dependent in our simulations, and so we build the steepest-descent gradient algorithm accordingly here. We determine the optimal functions U and S constructing a sequence of approximations {U k } k≥1 and {S k } k≥1 with U 0 and S 0 given. In what follows, we describe our algorithm to calculate (U k+1 , S k+1 ) from (U k , S k ), which involves 7 main steps: Step 1: Compute the functions (φ k , σ k , p k ) solving the forward problem (2.1)-(2.5) with U = U k and S = S k .
Step 3: Evaluate the gradient of J at (U k , S k ) using Eqns. (4.49)-(4.50), i.e., Step 4: Check if any of the following convergence criteria is satisfied: Criterion 2: where ε SD1 and ε SD2 are predefined tolerances. If any of the convergence criteria is satisfied, the iterative process to compute the optimal functions ends. If none of the convergence criteria are satisfied, we proceed to Step 5.
Step 5: Compute a pool of potential values for (U k+1 , S k+1 ) using the updates where µ j = j/N µ with j = 1, . . . , N µ . Here, N µ > 1 is an integer parameter of the algorithm.
6 Simulation study of the optimal control problem 6.1 Description In this section, we explore the behavior of the optimal control problem stated in Section 2.2 for time-dependent controls U (t) and S(t). To this end, we carry out a simulation study featuring a prostatic tumor treated with combined cytotoxic and antiangiogenic therapy. Thus, the aim of these simulations is to compute the optimal drug-naïve cytotoxic and antiangiogenic effects, U (t) and S(t), respectively, to effectively treat the tumor according to our model of PCa growth. We run this study over the time of a single cycle of combined therapy with docetaxel and bevacizumab, i.e., T = 21 days [16,48,68]. In Section 7 we will use the results of this study to explore alternative therapeutic strategies for specific drugs.
In the following, we present values of the model parameters used in the simulations, the choice of terms and corresponding weighting constants in the objective functional J, the construction of initial guesses and the admissible space for U (t) and S(t), and the computational details of the specific implementation of the algorithms presented in Section 5 for this simulation study.

Parameters of the PCa growth model
We consider an aggressive case of PCa. This scenario corresponds to a tumor having a high Gleason score, which is a routine clinical variable associated to cancer aggressiveness [68]. To model this instance of PCa growth, we select ρ and A in Eq. (2.7) based on the average values of tumor cell proliferation and apoptosis, tumor doubling times, and serum PSA doubling times previously reported in the literature for tumors with high Gleason score [7,16,82].
The values of the other parameters in our PCa growth model have been adopted from previous studies [16,[59][60][61]99]. Table 1

Objective functional
As explained in Section 2.2, we selected only some k i , i = 1, . . . , 7 to be nonzero in our simulations.
In particular, we consider three versions of the functional J defined in Eq. (2.8) corresponding to three different sets of nonzero constants k i , i = 1, . . . , 7, as follows: and We choose φ Q = 0 and φ Ω = 0, such that the optimization problem aims at killing the tumor. Consequently, we set p Ω = α h vol(Ω)/γ p , which is the value of serum PSA that would be obtained if the complete computational domain Ω were comprised of healthy tissue. We did not consider the term multiplied by k 5 in the original definition of J in Eq. (2.8) because, in general, serum PSA can be obtained with reasonable time resolution in experimental and clinical settings [68].
In this simulation study, we will explore the performance of our optimal control formulation under an array of values of the weighting constants in J 1 , J 2 , and J 3 . In simulations with J 1 , we impose k 1 = k 2 due to the similarity of the terms with k 1 and k 2 . Indeed, notice that the term with k 2 introduces in the transversality condition for the adjoint variable w (Eq. (4.29)) a similar term to the forcing induced by the k 1 term in the equation for w (Eq. (4.25)). We further choose k 4 = k 1 such that we use the same weighting for the terms accounting for tumor volume and serum PSA in the objective functional. Likewise, we impose k 4 = k 2 and k 4 = k 3 for J 2 and J 3 , respectively. We also fix k 6 = 1 and k 7 = 1 in all simulations in Section 6.3, such that we do not favor one therapy over the other within our modeling framework. These choices make J 1 , J 2 , and J 3 solely dependent on k 1 , k 2 , and k 3 , respectively; which facilitates the analysis of this first simulation study of the optimal control problem presented herein. The specific values of the weighting constants employed in each simulation are reported along with the corresponding results in Section 6.3.

Initial guess and maximum values of the controls
Pharmacodynamic studies of common cytotoxic and antiangiogenic drugs for cancer treatment usually show an exponential decay in drug concentration following the systemic delivery of the prescribed dose [4,27,62,90]. Additionally, previous efforts to model cytotoxic and antiangiogenic drug effects usually rely on a linear dependence on the drug concentration [6,8,14,16,25,34,46,50,77]. We adopted this approach for the formulation of the initial guess of cytotoxic drug effects on tumor dynamics U 0 (t), which is hence given by where β c measures the effect of the cytotoxic drug on tumor dynamics per unit of drug dose delivered, d c is the prescribed dose of the cytotoxic drug, and τ c stands for the mean lifetime of the chemotherapeutic drug. Likewise, we define the initial guess of the antiangiogenic effect on nutrient supply S 0 (t) as S 0 (t) = β a d a e − t τa , (6.5) where β a measures the effect of the antiangiogenic drug on the nutrient supply per unit of drug dose delivered, d a is the prescribed dose of antiangiogenic drug, and τ a denotes the mean lifetime of the antiangiogenic drug. Notice that this choice of U 0 (t) and S 0 (t) matches our previous modeling approach [16]. Additionally, observe that under the assumption of linear dependence of U (t) and S(t) on the cytotoxic and antiangiogenic drug concentrations, respectively, the last two terms in our optimal control problem functional J in Eq. (2.8) are penalizing large concentrations of cytotoxic and antiangiogenic drugs [6,8,14,16,25,34,46,50,77]. Docetaxel is considered the gold standard drug in cytotoxic chemotherapy of advanced PCa [19,48,68]. Additionally, bevacizumab has been extensively investigated in single-drug antiangiogenic or combined cytotoxic-antiangiogenic drug protocols for advanced PCa [2,48,76]. Therefore, we use the standard dosage and pharmacodynamic properties of these two drugs to calculate U 0 (t) and S 0 (t) in our simulations, respectively. The values of the parameters involved in the dynamics of drug effects in Eqns. (6.4)-(6.5) are provided in Table 2.
Additionally, the admissible values of the controls U (t) and S(t) are respectively bounded by maximum values U max and S max in our formulation of the optimal control problem (see Section 2.2). We set U max = 0.012 1/day and S max = 0.80 g/L/day. These values respectively correspond to While we have used docetaxel and bevacizumab to define U 0 (t), S 0 (t), U max , and S max , we want to remark that the optimal controls U (t) and S(t) describe the optimal cytotoxic and antiangiogenic effects according to our PCa growth model, which may also be obtained with different dosages of other drugs showing different pharmacodynamics. In Section 7, we provide and example of how to estimate some drug protocols that would approximately yield the optimal U (t) and S(t) obtained in the simulations of our optimal control problem.

Computational setup
We implemented the numerical algorithms introduced in Section 5 to resolve the forward problem, the adjoint problem, and the steepest descent gradient method by extending our in-house isogeometric codes to simulate PCa growth [59][60][61]. These codes were built following the general directions in [37].

Space and time discretization
The computational domain in all simulations is a square with side length of L d = 3000 µm and 256 isogeometric elements per side. The time step was set to a constant value of ∆t n = 0.1 days for the forward problem and ∆t n = −0.1 days for the adjoint problem.

Convergence of numerical algorithms
The convergence of the Newton-Raphson method was set to tolerance ε N L = 10 −3 , while for the GMRES algorithm was set to ε L = 10 −3 or a maximum of 500 iterations. The convergence for the steepest-descent gradient algorithm was set to ε SD1 = ε SD2 = 10 −6 or a maximum of 100 iterations. To find the update for the optimal U (t) and S(t) in each step of the steepest-descent gradient algorithm, we evaluated N µ = 10 evenly-spaced values of µ in (0, 1], i.e., µ j = j/10 for j = 1, . . . , 10 (see Section 5.3).

Initial conditions of the forward problem
We approximate the initial tumor phase field as an ellipsoidal tumor placed in the center of the domain with semiaxes a = 150 µm and b = 200 µm parallel to the domain sides. We implement this initial condition by L 2 -projecting the hyperbolic tangent function over the quadratic B-spline space supporting our spatial discretization. This operation yields the control variables φ 0,A = φ A (0), A = 1, . . . , n b , for the spline representation of the phase-field initial condition, i.e., φ h . The initial conditions for the nutrient and the tissue PSA are estimated from φ 0 as and The constants c 0 σ , c 1 σ , c 0 p , and c 1 p are computationally estimated [60], such that σ 0 and p 0 represent a constant value of the nutrient and tissue PSA within the tumor and the host tissue. Hence, we choose c 0 σ = 1 g/L, c 1 σ = −0.8 g/L, c 0 p = 0.0625 ng/mL/cc, and c 1 p = 0.7975 ng/mL/cc. In all simulations, we initially let the tumor grow untreated for 60 days, i.e., U (t) = 0 and S(t) = 0. This enables us to have a reference for the dynamics of the untreated tumor and obtain a good estimate of fields φ, σ, and p according to our PCa growth model away from the estimated initial conditions [16]. Then, we proceed to add cytotoxic and antiangiogenic drug effects. To facilitate the ensuing implementation of the steepest-descent algorithm, we reset t = 0 at the date of drug delivery and we use the values of φ, σ, and p obtained at the end of the untreated growth phase as initial conditions for the forward problem within our optimal control framework. Hence, the forward and dual problems are run between t = 0 and t = T = 21 days, corresponding to the common duration of a combined therapy cycle for PCa [16,48,68]. Figures 1-3 respectively show the optimal U (t) and S(t) distributions that we obtained in simulations of our optimal control problem using the objective functionals J 1 , J 2 , and J 3 for an array of values of k 1 , k 2 , and k 3 . These figures also depict the evolution of tumor volume and serum PSA obtained for each optimal solution, along with the tumor contours corresponding to the isosurface φ = 0.5 at t = T . Larger values of the driving weighting constant k i , i = 1, 2, 3, in each corresponding objective functional J i , i = 1, 2, 3, increasingly penalize tumor volume and serum PSA (see Section 6.1.2). Consequently, the minimization of J i , i = 1, 2, 3, will require at least one of the controls to also increase for larger choices of k i , i = 1, 2, 3, because this provides a more intense inhibitory effect on tumor dynamics that ultimately leads to smaller tumor volume and lower serum PSA, as shown in Figs. 1-3.

Simulation results of the optimal control problem
In all cases run with J 1 , the optimal therapy obtained in our simulations involves only a monotonically decreasing cytotoxic effect and virtually no antiangiogenic effect, i.e., S(t) ≈ 0. The optimal cytotoxic effect even saturates to U max at early times for k 1 = 5, whereas S(t) still takes negligible values. Thus, these results suggest that strong cytotoxic effects suffice to optimally control the volume and serum PSA of the simulated tumor according to our PCa model. A similar result is obtained for the simulations involving J 2 and J 3 depicted in Figs. 2-3. However, the optimal cytotoxic effects U (t) obtained for functionals J 2 and J 3 are approximately constant for early times and increase in time towards t = T , especially for large values of k 2 and k 3 . As the tumor phase field φ drives serum PSA dynamics in our PCa model (see Section 2.1), the tumor volume becomes the most powerful clinical quantity of interest participating in the objective functional J. Therefore, the increase in optimal U (t) by t = T for larger k 2 and k 3 is a consequence of J 2 and J 3 controlling for tumor volume only at t = T .  (C) Evolution of tumor volume v φ (top) and serum PSA Ps (bottom) using the optimal U (t) and S(t) obtained via simulation for each value of k1 compared to standard combined therapy, i.e., v φ,0 and Ps,0 respectively (gray dotted lines). (D) Tumor contours obtained at t = T using the optimal U (t) and S(t) obtained via simulation for each value of k1. The gray area denoted with φ0 is the tumor region at the onset of the optimal control problem (t = 0). The black line is the tumor contour at t = T for the standard combined therapy (U0, S0). Bottom: Detail of optimal S(t) for k2 = 1, 20, 100 calculated numerically (solid lines) and with Eq. (4.45) (dashed lines). (C) Evolution of tumor volume v φ (top) and serum PSA Ps (bottom) using the optimal U (t) and S(t) obtained via simulation for each value of k2 compared to standard combined therapy, i.e., v φ,0 and Ps,0 respectively (gray dotted lines). (D) Tumor contours obtained at t = T using the optimal U (t) and S(t) obtained via simulation for each value of k2. The gray area denoted with φ0 is the tumor region at the onset of the optimal control problem (t = 0). The black line is the tumor contour at t = T for the standard combined therapy (U0, S0).

Tumor volume
Optimal cytotoxic effects U(t) Optimal antiangiogenic effects S(t) S (g/L/day) 10 Bottom: Detail of optimal S(t) for k3 = 1, 10, 50 calculated numerically (solid lines) and with Eq. (4.45) (dashed lines). (C) Evolution of tumor volume v φ (top) and serum PSA Ps (bottom) using the optimal U (t) and S(t) obtained via simulation for each value of k3 compared to standard combined therapy, i.e., v φ,0 and Ps,0 respectively (gray dotted lines). (D) Tumor contours obtained at t = T using the optimal U (t) and S(t) obtained via simulation for each value of k3. The gray area denoted with φ0 is the tumor region at the onset of the optimal control problem (t = 0). The black line is the tumor contour at t = T for the standard combined therapy (U0, S0). Figs. 1-3 show that our optimal control problem enables the calculation of an optimal U (t) providing a similar control of tumor volume and serum PSA as the standard combined therapy with docetaxel and bevacizumab, i.e., U 0 (t) and S 0 (t) respectively. Figures 1-3 also show that we can obtain optimal U (t) and S(t) exhibiting an improved therapeutic performance over the standard combined protocol as we increase the driving weighting constant k i , i = 1, 2, 3. However, we notice that the optimal U (t) obtained with J 2 and J 3 would require a prolonged exposition to the cytotoxic drug. This would entail a high toxicity that may eventually provoke serious side-effects and may potentially lead to abort the treatment prematurely. Because J 1 accounts for tumor volume for all times, we believe that this objective functional enables the computation of superior optimal drug effects precisely tuning treatment to tumor dynamics. We note that the curves of optimal cytotoxic effects U (t) shown in Fig. 1 do not resemble the exponential trends usually caused by the pharmacodynamics of usual cytotoxic drugs, e.g., the docetaxel curve corresponding to U 0 (t) in Figs. 1-3. In Section 7, we illustrate how these optimal U (t) time distributions can be leveraged to design novel drug protocols.

The simulations presented in
Additionally, Figures 1-3 show that the optimal U (t) found in our simulations matches the corresponding theoretical formulation in Eq. (4.45) for each value of the driving weighting constant k i , i = 1, 2, 3, in J i , i = 1, 2, 3, respectively. Figures 1-3 also show a detail of optimal S(t) distributions computed numerically and using the theoretical formulation from Eq. (4.45) for some simulations. In general, these two optimal S(t) distributions differ in a value that would not be captured by the tolerance of our implementation of the steepest-descent gradient algorithm (see Section 5). However, the numerical and theoretical solutions for k 1 = 5 virtually coincide. We also observe increasingly better agreement between both distributions for larger values of the driving weighting constant k i , i = 1, 2, 3, respectively in J i , i = 1, 2, 3. Hence, further reducing the tolerance in the steepest-descent gradient algorithm would make both distributions of optimal S(t) converge. Nevertheless, we want to remark that both the theoretical and numerical optimal coincide for therapeutic purposes here, as both expressions of S(t) ≈ 0 regardless of the objective functional considered.
Finally, Appendix B shows ancillary simulations in which we decrease k 7 by several orders of magnitude using J i , i = 1, 2, 3. For each of the objective functionals, simulations with k 7 = 0.01 and k 7 = 0.001 render a non-negligible optimal S(t) that matches the corresponding theoretical estimates provided by Eq. (4.45). For the simulations with J 1 these optimal S(t) solutions show a decreasing profile in time, while for J 2 and J 3 the corresponding optimal S(t) solutions are mostly constant or slightly increasing towards t = T . The optimal cytotoxic effects U (t) remain practically constant in all simulations for each choice of the objective functional as we vary k 7 and the time distribution is virtually the same as the corresponding simulations shown in Figs. 1-3. However, the increase of optimal S(t) has limited effect on tumor volume and serum PSA compared to the parallel increase observed in optimal U (t) in Figs. 1-3. This suggests that the choices of k 7 in Appendix B might artificially bias our optimal control problem and produce unrealistic solutions of limited therapeutic use. As we discuss in Section 8, a rigorous analysis of the values of the weighting constants k i , i = 1, . . . , 7, is required to identify choices representing clinically relevant scenarios for our optimal control framework.
7 Design of drug protocols inspired by the computed optimal effects 7.1 Methodology To illustrate the applicability of the results of our optimal control problem, we calculate an array of model-inspired drug protocols aiming at producing the optimal cytotoxic effect U (t) obtained in one of the simulations of our optimal control problem reported in Section 6. In this section, we assume S(t) = 0, in accordance with the results showed in the previous one. We consider four different drug protocols, which we model leveraging the same paradigm that we employed to define U 0 (t) in Eq. (6.4). We consider a single dose of docetaxel, for which we simply calculate d c in Eq. (6.4). Then, we consider a new design drug for which we compute τ c and the docetaxel equivalent dose d c , such that we can use the same β c in Eq. (6.4). Finally, we will consider a three-dose protocol of either docetaxel or the new design drug. In this case, we extend the formulation of Eq. (6.4) to where H(t) is the Heaviside function. We compute the three doses d c,i , i = 1, 2, 3 for either drug, together with τ c for the new design drug. Additionally, while we fix the delivery of the first dose at t c,1 = 0, we calculate the delivery times of the second and third doses for either drug as well (i.e., t c,2 and t c,3 , respectively). We compute each of these drug protocols by running a nonlinear least-square fit of the optimal U (t) obtained using J 1 (U, S) and k 1 = k 2 = k 4 = 2, which is shown in Fig. 1. We use the trust-region method as provided by the Curve Fitting Toolbox in MATLAB (Release R2020a, The Mathworks, Inc., Natick, Massachusetts, US). The starting dose is set at d c = 75 mg/m 2 in singledose protocols. In 3-dose protocols the starting point is three equal doses d c = 25 mg/m 2 delivered weekly, i.e., with t c,2 =7 days and t c,3 =14 days. The starting value for τ c for the new design drug is 5 days, corresponding to docetaxel (see Table 2). The admissible value ranges are [0, 100] mg/m 2 for drug doses, [0, 21] for the delivery times of the second and third doses in 3-dose protocols, and [1,20] for the characteristic decay time of the new design drug. Additionally, we run a simulation of the forward problem with the fitted drug protocols using the computational setup described in Section 6.2 and compare the dynamics of tumor volume and serum PSA with those provided by the optimal U (t). We assess the goodness of fit by means of R 2 and the root mean squared error (RMSE). Table 3 shows the result of the nonlinear fitting of the parameters participating in the formulation of each protocol. Figure 4 shows the cytotoxic effects generated by each of the investigated protocols, as well as the tumor volume and serum PSA evolution that they would produce according to the corresponding forward simulations of our PCa growth model. Table 4 shows the fit quality of the cytotoxic effects, tumor volume, and serum PSA for each protocol with respect to those provided by the optimal U (t). The single-dose docetaxel protocol uses a slightly larger dose than the standard protocol (75 mg/m 2 , see Table 2), but, in absence of the antiangiogenic drug, the corresponding control of tumor volume and serum PSA is poorer than in the standard protocol. The single-dose docetaxel protocol also produced the worse fitting of the optimal cytotoxic effects. Conversely, the 3-dose docetaxel protocol produces an enhanced fitting of the optimal U (t) and a good control of  Figure 4: Design of drug protocols with docetaxel (Dtx) or a new design drug (New) aiming at reproducing the optimal cytotoxic effect U * (t) as obtained from the simulation of our optimal control problem using J1(U, S) and k1 = k2 = k4 = 2. (A) Cytotoxic drug-induced effects for each of the calculated protocols compared to those of the standard docetaxel protocol U0(t) and the target optimal cytotoxic effects obtained in our simulations of the optimal control problem U * (t). (B) Evolution of tumor volume v φ and serum PSA Ps for each calculated drug protocol, the standard docetaxel protocol (v φ,0 , Ps,0), and the optimal control problem solution (v * φ , P * s ). tumor volume and serum PSA. The best fitting of optimal cytotoxic effects and tumor control was provided by the protocols involving a new design drug. Our calculations suggest that this design drug should have a larger τ c , approximately twice the corresponding value for docetaxel. This results in slower decay and hence a prolonged cytotoxic effect in time in comparison to docetaxel. We also observe that this new drug would require lower dosage for both the single-dose and 3-dose cases compared to the corresponding protocols using docetaxel. Since the control of tumor volume and serum PSA is virtually the same with the single-dose and 3-dose protocols using the new design drug, we believe that the single-drug protocol would be more suitable for clinical implementation as it would require less visits of the patient to the clinic.

Conclusions
We present an optimal control theory enabling the calculation of the optimal drug-naïve cytotoxic and antiangiogenic effects aiming at successfully treating advanced PCa, which can then be used to design optimal drug compounds and delivery plans. We build this framework departing from our model of PCa growth with cytotoxic and antiangiogenic drug therapies [16]. We define an objective functional featuring diverse measurements of tumor morphology, tumor volume, and serum PSA, which are common quantities of interest in experimental and clinical studies of advanced PCa. The functional contains terms accounting for the differences between the computed solution for the prostate tumor or PSA and prescribed values of them, which should be reached during the simulated time or at a final endpoint by an optimal treatment procedure. We also include penalty terms for cytotoxic and antingiogenic effects, which could be ultimately linked to the corresponding drug concentrations following the usual linear paradigm in the literature [6,8,14,16,25,34,46,50,77]. From the mathematical viewpoint we prove that the functional reaches its minimum at the optimal values of the cytotoxic and antiangiogenic effects, that is, it provides the optimal design of the treatment which may possibly lead to a final desired result regarding the decrease of the tumor and/or of the PSA. The necessary conditions that should be satisfied by the optimal values are determined by means of the adjoint problem.
In this work, we also propose an algorithm implementing the steepest-descent gradient method to solve our optimal control problem [3], which relies on IGA [37] and the generalized-α method [13,44] to discretize the involved forward and adjoint problems in space and time, respectively. We then use this algorithm in a simulation study to analyze the performance of our optimal control problem. In general, we observed a remarkable agreement between our numerical approximations of the optimal drug effects and the theoretical formulations found in Section 4. Our results show that accounting for the spatio-temporal evolution of the tumor morphology in the objective functional (i.e., J 1 in Section 6) led to superior optimal drug effects, which precisely adapt the optimal therapy to tumor dynamics and would lead to drug protocols with expected lower toxicity compared to the optimal drug effects obtained by only accounting for the final tumor morphology or volume at t = T (i.e., J 2 and J 3 in Section 6 respectively). A previous simulation study on our PCa model [16] revealed that drug-induced changes in tumor morphology during chemotherapy may contribute to chemoresistance, which supports the inclusion of spatio-temporal dynamics of tumor morphology in the objective functional to refine the design of optimal drug therapies.
The simulated cases of our optimal control problem in Section 6 feature S(t) ≈ 0. We think that these results might align with the studies suggesting that antiangiogenic therapy might not be optimal, and hence that a cytotoxic monotherapy would suffice to optimally treat advanced PCa [2,48,71,74,83,[85][86][87]89]. However, several methodological shortcomings have been detected in clinical studies of antiangiogenic therapies, requiring a closer monitoring of pathological events to account for a more direct control of tumor evolution than survival endpoints and serum PSA [2,71,83,85]. Some of these monitoring variables are already considered or could be added to the objective functional in our optimal control framework. Furthermore, we plan to study the selection of the values of the weighting constants k i , i = 1, . . . , 7, in the objective functional J to define a constant set enabling the consistent design of personalized optimal drug protocols. We believe that such a selection requires to account for the intrinsic scales in our PCa growth model. Consequently, this update may result in optimal solutions with relevant antiangiogenic effects (e.g., Appendix B shows some simulations with non-zero optimal S(t) for 0 < k 7 1). Additionally, our PCa model only includes antiangiogenic effects as a decrease in nutrient supply. We plan to overcome this simplifying limitation by extending our model to account for the dynamics of the evolving tumor-supporting vascular network [32]. This feature would enable a superior description of nutrient and drug supply to the tumor and provide a clear target for antiangiogenic therapy [6,29,30,34,50,51,77]. Accounting for the vascular network would also enable the study of the synergistic or antagonistic effect of combining cytotoxic and antiangiogenic therapies, e.g., whether the reduction of tumor microvasculature caused by antiangiogenic therapy may hamper the effective delivery of cytotoxic drugs to the tumor or whether normalizing the tumor microvasculature would effectively improve the supply of cytotoxic drugs [42,69]. The dynamics of the tumor-induced vascular network can be modeled following different approaches [93]. Hybrid models can capture the evolving morphology of tumor-induced microvasculature with great detail, but usually require intensive computational resources [23,51,54,75,91,99]. Conversely, the local density of the vascular network can be modeled using a continuous formulation [6,30,36,50,88,93], which would dramatically reduce the computational cost of simulations, especially at tissue and organ scale. This advantage comes at the cost of losing geometrical precision, but it would still enable to account for the tumor-induced vasculature observable in current magnetic resonance imaging [36,98].
Our model of PCa growth could also be refined by incorporating the effect of tumor-induced mechanical stresses on tumor dynamics, which has been shown to improve tumor forecasting [55,97]. Indeed, a recent computational study suggests that the mechanical stresses created by PCa and coexisting benign prostatic hyperplasia may obstruct tumor growth [59]. Accounting for tumorinduced mechanical deformation can also improve the modeling of the tumor-supporting vascular network [36,81]. Furthermore, a poroelastic formulation would enrich the description of nutrient and drug dynamics in the tumor region by taking into account how these phenomena are affected by local changes in mechanical stress and fluid pressure [22,43,51]. This approach could provide new insights on drug delivery and action in the complex tumor environment, enabling us to refine therapeutic strategies accordingly. Our PCa growth model may also include multiphase formulations accounting for several tumor species with varying responses to the prescribed therapy [24,29,41,54,79,94]. To address the computational challenges of these model extensions and hence rationalize the computational resource demand, our numerical methods could be accelerated by implementing adaptive time stepping [26,92], dynamic local adaptivity in the spatial discretization [11,60], and superior algorithms to solve the optimal control problem [3,35].
Finally, we have also presented a strategy that decouples the problem of finding optimal treatment solutions in two phases: first, we compute the optimal therapeutic effect in a drug-independent optimal control framework, and, second, we calibrate alternative drug protocols for specific drug dosage, effects, and pharmacodynamics targeting the calculated optimal therapeutic effect. The application of this approach revealed interesting drug protocols that remarkably reproduced the desired optimal effects simulated herein. We illustrate the design of a generic new drug in Section 7 by employing the equivalent dose and effects of docetaxel because this is the most common compound used for the chemotherapy of PCa, hence providing a reference for comparison during drug design. However, our framework can also accommodate specific drug doses and effects, e.g., by respectively choosing d c and calibrating β c with data in Eq. (7.1). The design of therapeutic solutions could be further refined in multiple directions, e.g., by exploring alternative paradigms to model the dynamic effects of specific cytotoxic and antiangiogenic drugs [28,34,39,53,73,101], by explicitly including drug toxicity (e.g., throughout the time integral of drug concentration [6,53]), by accounting for the synergistic action of drug combinations [39,66], or by considering other forms of cancer treatment (e.g., the cytotoxic action of radiation therapy [18,33,56,77]). Ultimately, the optimal pharmacodynamic properties calculated with our approach could assist in the development of new drug compounds that are more efficacious, show lower toxicities, and may even target specific cancer subtypes or patient-specific tumors [38,40,63,84]. We also plan to explore the design of drug protocols for 3D PCa growth scenarios and considering longer simulation times, which are closer to the reality of experimental and clinical settings [16,48,[59][60][61]68]. Additionally, the acquisition of longitudinal PSA and imaging data during treatment would enable to recalibrate key parameters in our framework, for example: tumor proliferation and apoptosis in our PCa model, or drug effect rates β c and β a during drug protocol design. Although we have assumed these parameters to be constant, they are known to vary due to treatment action and the phenotypic evolution of the tumor [24,32,39,49,83]. Data-driven reparameterization would also enable us to update the optimal controls, adapt the design of personalized drug protocols accordingly, and early detect the emergence of chemoresistance patterns. In sum, we believe that our optimal control framework can provide a versatile and powerful approach to investigate patient-specific therapeutic strategies in silico to succesfully treat advanced PCa.
Appendix A. Supplementary mathematical results.
The notation we use in this appendix is somehow independent from the one in the paper. We are going to check that if u, v are two measurable functions defined on Ω and F : R → R is a C 1 function, then there exists a measurable function w : Ω → R attaining intermediate values between the ones of u and v and such that F (u(x)) − F (v(x)) = (u(x) − v(x))F (w(x)) a.e. x ∈ Ω. (A.1) The whole argument is due to Vittorino Pata, with many thanks from the authors.
Theorem A.1. Let f : R → R be a continuous function. Then f has a Borel measurable right inverse k defined on (the interval) Im(f ).
For the proof of Theorem A.1 we need a lemma. Proof. The latter equality is obvious, since f is continuous, and so the infimum is actually a minimum. Let us prove the first assertion. Let t be fixed, such that k(t) = a. We will show that k is either left or right continuous at t. This is enough to ensure that k is Borel measurable. To this end, call y = k(t) ∈ (a, b]. By definition, f (τ ) = y whenever τ < t. Assume then that f (τ ) < y for all τ < t. We prove that k is left continuous at y (if instead f (τ ) > y for all τ < t then k is right continuous at y, the proof being the same). By contradiction, if k is not left continuous, there is y n ↑ y such that t n = k(y n ) → k(y) = t. By compactness, there ist < t such that t n →t, up to a subsequence. But since f is continuous, y n = f (t n ) → f (t), so that f (t) = y. This contradicts the fact that t is the smallest element of [a, b] for which the equality f (t) = y holds.
We now modify the sets A n in order to have them disjoint. Hence we put and so on. Possibly, some A n can be empty, and we simply ignore it. Clearly, we have n A n = Im(f ).
Let now f n be the restriction of f on [n, n + 1], and let k n : f ([n, n + 1]) → [n, n + 1] be given by k n (t) = inf{f −1 n (t)}, From the Lemma, k n is Borel measurable and f n (k n (t)) = t for all t ∈ f ([n, n + 1]). Finally, define the (Borel measurable) function k(t) = k n (t), t ∈ A n .
Then, for t ∈ A n , we have that f (k(t)) = f n (k n (t)) = t.
Remark A.3. With inessential changes in the proof, Theorem A.1 still holds if f is defined on a generic interval I.
Corollary A.4. Given a domain Ω ⊂ R n , let µ : Ω → R be a measurable function (finite everywhere), and let f : R → R be a continuous function such that Im(f ) ⊃ Im(µ). Then there exists a measurable function q : Ω → R such that f (q(x)) = µ(x), for all x ∈ Ω.
Proof. By Theorem A.1, let k be the Borel measurable right inverse of f , and define q(x) = k(µ(x)), x ∈ Ω.
q is Lebesgue measurable, being the composition of a Borel measurable function and a Lebesgue measurable one.
We finally conclude with the desired application involving the Lagrange mean value theorem.
Corollary A.5. Given a domain Ω ⊂ R n , n ∈ N, let u, v : Ω → R be measurable functions (finite everywhere), and let F : R → R be a C 1 function. Then there exists a measurable function w : Ω → R such that (A.1) holds.
Proof. We define Then the problem amounts to finding a measurable function w such that F (w(x)) = µ(x), x ∈ Ω.
The Lagrange theorem implies that Im(F ) ⊇ Im(µ). Hence, the conclusion follows from Corollary A.4 by taking w = q.  and serum PSA Ps (right) using the optimal U (t) and S(t) obtained via simulation for each value of k7 compared to standard combined therapy, i.e., v φ,0 and Ps,0 respectively (gray dotted lines).