Prostate cancer (PCa) is a major health problem striking over one million men annually and is responsible for over 300,000 deaths [20, 55]. PCa is usually an adenocarcinoma, a form of cancer that originates in the epithelial tissue of the prostate. The evolution of a tumor depends on the genetic alterations that originated it and on its microenvironmental conditions. A key process controlling the conditions of the tumor microenvironment is tumor-induced angiogenesis [22, 34], i.e., the growth of new blood vessels from pre-existing ones via chemical signals produced by the tumor.
Arguably, the most effective way to combat PCa is a combination of prevention and regular screening for early detection. PCa screening is usually accomplished by way of periodic digital rectal exams (DREs) and prostate-specific antigen (PSA) tests . The DRE is a physical test in which a doctor palpates the rectal wall next to the prostate to search for hard, lumpy, or abnormal areas typically indicative of cancer. The serum level of PSA is a biomarker of the prostate activity that rises during PCa. The PSA test is a blood test that measures the amount of PSA in the bloodstream. If either the PSA or the DRE test are positive, the patient will be recommended to undergo a biopsy. Together with DRE and serum PSA, the results of the biopsy establish the basis for a diagnosis of PCa. Patients diagnosed with localized PCa normally receive a radical treatment with curative intent, such as surgery or radiation . Advanced PCa patients are usually prescribed androgen-deprivation therapy or chemotherapy , which are administered systemically. In addition to its use in advanced PCa, chemotherapy has also been proposed as an efficient neoadjuvant therapy before radical prostatectomy [12, 60] that aims at reducing the severity of the disease before surgery by shrinking the tumor and eliminating micrometastases that may have developed.
Chemotherapy for PCa is mostly based on cytotoxic drugs (e.g., docetaxel, cabazitaxel), which obstruct tumor growth by inhibiting cell proliferation and promoting tumor cell death [55, 18, 67]. While antiangiogenic therapies for PCa are being actively investigated [2, 67, 68], they are outside of the standard of care for PCa. This is somehow counterintuitive because angiogenesis is known to play a central role in the progression of castration-resistant PCa  and an increase in microvascular density is associated to poorer prognosis . There are multiple studies with contradicting evidence about the effectiveness of antiangiogenic therapy [2, 67, 68, 57]. However, we would like to highlight that bevacizumab (the most common antiangiogenic drug ) did not show significant benefit in castration-resistant PCa when administered alone , but it produced a 50% PSA decline in 75% of the patients when combined with chemotherapy . Thus, bevacizumab and other antiangiogenic drugs are currently being regarded as promising agents to use in combination with cytotoxic chemotherapy. However, prostatic tumors have continuously evolving heterogeneous genetic profiles, which may show varying rates of resistance to the prescribed chemotherapeutic treatment and these may even increase during conventional protocols [67, 44, 34, 26]. This is known to be behind the widely varying and sometimes contradictory results of chemotherapeutic clinical studies, also making extremely complex to know whether a patient will benefit from chemotherapy and which drugs are best to treat his tumor.
Recently, the computational modeling and simulation of cancer has shown promise to extend our understanding of these pathologies as well as in forecasting tumor growth and treatment outcomes [1, 78, 17, 51, 49]. In this context, several studies have been focusing on studying the effects of chemotherapy on tumors through mathematical models and computational simulations [5, 9, 26, 37, 45, 62]. This approach would enable to test in silico alternative drug protocols and combinations, hence assisting physicians finding optimal chemotherapeutic plans for each patient. Additionally, computational models of cancer growth and treatment may contribute to better comprehend the intricate mechanisms of drug resistance and find early predictors of chemotherapy failure.
This paper proposes a model of PCa growth and chemotherapy, where the advancement of the tumor is controlled by a critical nutrient. The effect of cytotoxic and antiangiogenic drugs is incorporated by downregulating tumor net proliferation and reducing nutrient supply, respectively. The model also produces as an output the time evolution of serum PSA, which is information commonly used in clinical practice to monitor tumor evolution or the disease’s response to a particular treatment. Compared to our previous work on PCa [51, 50, 49], for the model presented herein we perform a complete mathematical analysis, while retaining critical features such as the ability of the model to predict the morphological shift of the tumor under certain circumstances, which is compatible with in vitro experiments of PCa cell lines  and clinical observations [19, 58]. We prove that the model is well posed and develop an algorithm to solve the equations numerically. Our computational method is based on isogeometric analysis (IGA), a recently-proposed generalization of the finite element method that uses splines as basis functions . Our computational results show complex tumor dynamics, matching previous observations in both computational and clinical studies.
The paper is organized as follows: Sections 2 and 3.1 describe, respectively, the model equations and the functional framework that we utilize. We prove the well-posedness of the model in Section 3. Section 4 presents our computational method and Section 5 shows representative simulations. We draw conclusions in Section 6.
2 Model of prostate cancer growth with chemotherapy
2.1 Modeling approach
The model describes the tumor dynamics using a phase field, i.e., a continuous field that defines the time evolution of the tumor’s location and geometry. The phase-field method has been extensively used to describe tumor growth in the computational literature [10, 14, 16, 15, 25, 28, 27, 51, 54, 24, 76, 77]. In this work, the phase field is denoted by and transitions from the value in the host tissue to in the tumor. The transition is smooth but steep and takes on a hyperbolic tangent profile in the direction perpendicular to the interface . The model also accounts for the dynamics of a critical nutrient, whose concentration is denoted by and obeys a reaction-diffusion equation. The concentration of PSA in the prostatic tissue is governed by a linear reaction-diffusion equation.
2.2 Governing equations
2.2.1 Tumor dynamics
The tumor dynamics is described by the equation 
where denotes the partial derivative of with respect to time, is the Laplace operator, is the diffusion coefficient of tumor cells and
The diffusion coefficient of tumor cells can be computed as , where both and are positive real constants denoting the tumor mobility and interface width, respectively . Here, is a double-well potential, i.e., a nonconvex function, typical in phase-field modeling, which allows the coexistence of the tumoral () and healthy () tissue. The function
is also common in non-conserved dynamics of phase fields. It is usually called interpolation function because it verifies the propertiesand ; another important property of is that , where denotes the derivative of . The function is normally called tilting function and, in our model, is defined as
where and represent constant proliferation and apoptosis indices, respectively. These nondimensional parameters are related to the proliferation and apoptotis rates in tumoral tissue as follows:
where is the proliferation rate of tumor cells, is the apoptosis rate of tumor cells, is a scaling reference value for the proliferation rate, and is a scaling reference value for the apoptosis rate. Therefore, can be interpreted as a function describing the tumor net proliferation rate. The positive constant scales the strength of the tilting function within our phase-field framework. The constants and are, respectively, a reference and a threshold value for the nutrient concentration . For nutrient concentrations lower than healthy tissue is energetically more favorable than tumor tissue and viceversa. The function in (2.2) represents the tumor-inhibiting effect of a cytotoxic drug and is described in Section 2.2.2. When , the function is a double-well potential with local minima at and . Within this range, low values of the nutrient concentration (or large values of ) produce a lower energy level (value of ) in the healthy tissue () than in the tumoral tissue (). The opposite is true for high values of the nutrient concentration (or low values of ). The described behavior of is further illustrated in Figure 1. From (2.2), we can obtain
2.2.2 Cytotoxic chemotherapy
Docetaxel is a usual cytotoxic drug in the clinical management of advanced PCa [55, 67, 18]. The standard cytotoxic chemotherapy based on docetaxel usually consists of up to 10 drug doses delivered every three weeks. We consider that the action of the cytotoxic drug on tumor dynamics depends linearly on the drug concentration [32, 15, 27, 9, 8, 45, 37, 62]. The pharmacodynamics of this cytotoxic drug shows an approximately exponential decrease in drug concentration after delivery of the dose at the systemic level [3, 70]. Hence, we propose the formulation
where is the number of chemotherapy cycles, measures the cytotoxic effect of the treatment on tumor dynamics per unit of drug dose delivered, is the drug dose, with are the times of drug delivery, is the mean lifetime of the chemotherapeutic drug, and denotes the Heaviside function. These parameters are required to be strictly positive. The experimental and clinical values observed for ensure that the condition is verified; see Section 2.2.1. According to (2.2), can be interpreted as a cytotoxic drug-induced decrease in tumor net proliferation that decays as the drug concentration in time. Thus, the action of the cytotoxic drug energetically disfavors tumor tissue (see Figure 1). Notice that the formulation paradigm in (2.8) is very similar to the modeling of radiation effects on tumor dynamics [47, 59, 17]. In principle, the formulation proposed in (2.8) is extensible to other cytotoxic drugs by adapting , and accordingly.
2.2.3 Nutrient dynamics
Nutrient dynamics is controlled by the equation
where is the diffusion coefficient of the nutrient, is the nutrient supply rate in the healthy tissue, stands for the nutrient supply rate in the cancerous tissue, is a given function yielding the reduction in nutrient supply caused by antiangiogenic therapy (see Section 2.2.4), and , are positive constants that represent the nutrient uptake rate in the healthy and cancerous tissue, respectively. We require that , , and are non-negative and should satisfy the constraint .
2.2.4 Antiangiogenic therapy
Bevacizumab is a common antiangiogenic drug that has been actively investigated for the treatment of PCa [2, 67, 68, 21]. This antiangiogenic therapy is usually delivered simultaneously with the cytotoxic chemotherapy or following a similar schedule in case it is the only treatment. Again, we consider that the action of this antiangiogenic drug depends linearly on its concentration [45, 33, 5, 6]. Pharmacodynamic studies of bevacizumab also revealed an exponential decay in drug concentration following systemic delivery of the prescribed dose [30, 52]. Therefore, we choose an analogous formulation to (2.8):
where is the number of antiangiogenic treatment cycles, measures the antiangiogenic effect of the treatment on the nutrient supply per unit of drug dose delivered, is the drug dose, are the times of drug delivery (), stands for the mean lifetime of the antiangiogenic drug, and denotes the Heaviside function. These parameters are required to be strictly positive. Again, the formulation in (2.10) may be extended to other antiangiogenic drugs by appropriately calibrating , and .
2.2.5 Tissue PSA dynamics
Both healthy and cancerous prostatic cells release PSA, although tumor cells do so at a much larger rate. The PSA is assumed to diffuse through the prostatic tissue and decay naturally at rate . Hence, we propose the equation 
where is the diffusion constant. The constants and represent, respectively, the tissue PSA production rate of healthy and malignant cells. The serum PSA can be defined as the integral of the tissue PSA over the prostate , that is,
By using (2.11), one can show that for free-flux boundary conditions
where is the volume of prostatic healthy tissue and is the tumor volume. Equation (2.13), which we derived from the tissue level equation (2.11), was proposed as a phenomenological model of serum PSA dynamics that fits clinical data .
3 Well-posedness of the model
In this section we shall provide an existence and uniqueness result for this system
Here, is an open bounded subset of , , with a sufficiently smooth boundary represents the outward normal derivative to and denotes a final time.
3.1 Functional framework
The problem will be treated in a functional framework involving the space which is identified with its dual space and the Sobolev spaces
and with the dense and compact injections and contains the elements of with null trace on the boundary . For the norms in these spaces we will use the notation , where is the space we are considering.
For and or , we simply denote the norm of in these spaces by
We will also make use of spaces of functions that depend on time with values in a Banach space . Namely, for we consider the space of measurable functions such that is integrable on (or essentially bounded if ) and the space of continuous functions from to . Perhaps, it is important to note that is a space completely isomorph to . Moreover, for , will denote the space of functions such that both and its (weak) derivative belong to . We point out that for all .
and it can be easily checked that
We also set
and point out that the problem parameters , , are positive.
In the sequel, by we shall denote a constant, that may change from line to line, depending on the problem parameters, the domain the final time the norms of the initial data and possibly on the norms of and . Moreover, we assume that (3.6)-(3.8) hold.
Please note that taking in (3.10) first, and then, equal to the null function (and this is of course a suitable choice in both cases), we obtain two separate variational equalities for and , as it is (3.11) for .
Let us denote
Theorem 3.2. Let
If then Moreover, if
then we have
Finally, if the solution has the supplementary regularity and satisfies the estimate
and satisfies the estimate
In addition, the solution is continuous with respect to the data, that is, for two solutions corresponding to the data we have
Proof. We shall prove the existence of the solution by using the Banach fixed point theorem combined with a variational approach. In this respect, we will be formal in the following sense: when referring to the weak solutions of equations like (3.2)-(3.4) coupled with the related boundary conditions, let us write directly the equations instead of their variational formulations. This would allow us to come quickly to the point in our argumentations.
Let us set
which is a complete metric space provided we take the distance induced by some norm in .
where is the solution to (3.20), corresponding to and , with the homogeneous Dirichlet boundary condition and the initial datum We shall prove that and that is a contraction mapping on .
First, we treat the initial-boundary value problem for equation (3.21). For all we introduce the operator defined by
and observe that it is continuous
and quasi-monotone from to , i.e.,
A first estimate performed by testing (3.21) by and integrating over yields, after a standard calculation,
for all where we used the fact that As specified before, is a constant, that may vary from line to line, depending on the problem parameters and : indeed, the Gronwall lemma has been applied to derive (3.22).
Now, if we take two solutions and corresponding to and thanks to the Hölder and Young inequalities we obtain for the differences and
By the Sobolev embedding theorems, it turns out that with continuous embedding. Then, by exploiting the Young inequality and (3.22) for we infer that
Going back to (3.23) we deduce that
where Finally, by using Gronwall lemma it is not difficult to conclude that
where stands for another constant depending additionally on the initial data, , and the source terms
Next, we treat the initial-boundary value problem for the equation (3.20). To this end, for and defined as above we consider the intermediate equation
We observe that is bounded and Lipschitz continuous on . Recalling (2.7) we note that and
together with (3.26). We note that the right-hand side of (3.29) is in and that It is clear, again by the solvability result for linear parabolic problems (see Ref. ), that (3.29), (3.26) has a unique solution
which is equivalent to the standard norm in
Let us consider two solutions and to (3.29), (3.26) corresponding to and respectively. We multiply the difference of the equations (3.29) by and integrate over With the help of (3.26) and (3.28) we have that
Hence, the Young inequality and the Gronwall lemma allow us to get
We multiply this inequality by and have successively
Taking the supremum with respect to , we obtain
since a.e. in , a.e. in the set where and a.e. in the set where This implies that that is a.e. in , for all
Next, we multiply (3.25) by and by similar calculations we deduce that
Here, we used the fact that a.e. on so that a.e. on Thus, we have shown that a.e. in , for all
By these two last results we actually proved that a.e. on implying that It turns out that this solution actually solves equation (3.20). Moreover, we have that