Mathematical analysis and simulation study of a phase-field model of prostate cancer growth with chemotherapy and antiangiogenic therapy effects

TL;DR

This paper develops a mathematical phase-field model to simulate prostate cancer growth and treatment effects, enabling personalized therapy testing through in silico experiments.

Contribution

The study introduces a novel phase-field mathematical model incorporating chemotherapy and antiangiogenic therapy effects for prostate cancer growth simulation.

Findings

Model captures tumor growth morphologies accurately.

Simulations replicate known tumor volume and PSA trends.

Model demonstrates effects of combined therapies.

Abstract

Cytotoxic chemotherapy is a common treatment for advanced prostate cancer. These tumors are also known to rely on angiogenesis, i.e., the growth of local microvasculature via chemical signaling produced by the tumor. Thus, several clinical studies have been investigating antiangiogenic therapy for advanced prostate cancer, either as monotherapy or combined with standard cytotoxic protocols. However, the complex genetic alterations promoting prostate cancer growth complicate the selection of the best chemotherapeutic approach for each patient's tumor. Here, we present a mathematical model of prostate cancer growth and chemotherapy that may enable physicians to test and design personalized chemotherapeutic protocols in silico. We use the phase-field method to describe tumor growth, which we assume to be driven by a generic nutrient following reaction-diffusion dynamics. Tumor proliferation and apoptosis (i.e., programmed cell death) can be parameterized with experimentally-determined values. Cytotoxic chemotherapy is included as a term downregulating tumor net proliferation, while antiangiogenic therapy is modeled as a reduction in intratumoral nutrient supply. Another equation couples the tumor phase field with the production of prostate-specific antigen, which is an extensively used prostate cancer biomarker. We prove the well-posedness of our model and we run a series of representative simulations using an isogeometric method to explore untreated tumor growth as well as the effects of cytotoxic chemotherapy and antiangiogenic therapy, both alone and combined. Our simulations show that our model captures the growth morphologies of prostate cancer as well as common outcomes of cytotoxic and antiangiogenic mono and combined therapy. Our model also reproduces the usual temporal trends in tumor volume and prostate-specific antigen evolution observed in previous studies.

Tables

Table 1. Table 1: Parameter selection.

Parameter	Notation	Value	Simulation case
Tumor dynamics
Diffusivity of the tumor phase field	$\lambda$	640 $\mu m$2/day	All
Tumor mobility	$M$	2.5 1/day	All
Net proliferation scaling factor	$m_{ref}$	7.55$\cdot {10}^{-2}$ 1/day	All
Scaling reference for proliferation rate	${\overline{K}}_{\rho }$	1.50$\cdot {10}^{-2}$ 1/day	All
Proliferation rate	$K_{\rho }$	0.8$\cdot {10}^{-2}$ 1/day	Mild tumor
		1.50$\cdot {10}^{-2}$ 1/day	Aggressive tumor
Scaling reference for apoptosis rate	${\overline{K}}_A$	2.10$\cdot {10}^{-2}$ 1/day	All
Apoptosis rate	$K_A$	0.7$\cdot {10}^{-2}$ 1/day	Mild tumor
		1.37 $\cdot {10}^{-2}$ 1/day	Aggressive tumor
Cytotoxic chemotherapy
Mean lifetime of cytotoxic drug	${\tau }_c$	5 day	Cytotoxic chemotherapy and combined therapy
Cytotoxic drug effect	${\beta }_c$	1.59$\cdot {10}^{-2}$ 1/(mg/m2)	Cytotoxic chemotherapy and combined therapy
Cytotoxic drug dose	$d_c$	75 mg/m2	Cytotoxic chemotherapy and combined therapy
Nutrient dynamics
Nutrient diffusivity	$\eta$	$6.4\cdot {10}^4$ $\mu m$2/day	All
Nutrient supply in healthy tissue	$S_h$	$2$ g/L/day	All
Nutrient supply in tumor tissue	$S_c$	2.75 g/L/day	Reference
		3.125 g/L/day	Rich supply
		2.375 g/L/day	Poor supply
Nutrient uptake by healthy tissue	${\gamma }_h$	$2$ g/L/day	All
Nutrient uptake by tumor tissue	${\gamma }_c$	17 g/L/day	Reference
		18 g/L/day	Higher metabolism
		16 g/L/day	Lower metabolism
Antiangiogenic therapy
Mean lifetime of antiangiogenic drug	${\tau }_a$	30 day	Antiangiogenic and combined therapy
Antiangiogenic drug effect	${\beta }_a$	0.04 g/L/day/(mg/kg)	Antiangiogenic and combined therapy
Antiangiogenic drug dose	$d_a$	15 mg/kg	Antiangiogenic and combined therapy
Tissue PSA dynamics
Tissue PSA diffusivity	$D$	640 $\mu m$2/day	All
Healthy tissue PSA production rate	${\alpha }_h$	1.712$\cdot {10}^{-2}$ ng/mL/cc/day	All
Tumoral tissue PSA production rate	${\alpha }_c$	${\alpha }_c=15{\alpha }_h$	All
Tissue PSA natural decay rate	${\gamma }_p$	0.274 1/day	All

Equations

$${\varphi }_t=\lambda \Delta \varphi -\frac{\partial G}{\partial \varphi }(\varphi ,\sigma ,u)$$

$$G(\varphi ,\sigma ,u)=F(\varphi )-h(\varphi )(m(\sigma )-m_{ref}u).$$

$$m(\sigma )=m_{ref}\left(\frac{\rho +A}{2}+\frac{\rho -A}{\pi }\arctan \left(\frac{\sigma -{\sigma }_l}{{\sigma }_r}\right)\right),$$

$$\rho =\frac{K_{\rho }}{\overline{K_{\rho }}}$$

$$A=-\frac{K_A}{\overline{K_A}},$$

$$\frac{\partial G}{\partial \varphi }=2\varphi (1-\varphi )f(\varphi ,\sigma ,u)$$

$$f(\varphi ,\sigma ,u)=M\left[1-2\varphi -3\left(m(\sigma )-m_{ref}u\right)\right].$$

$$u(t)=\sum _{i=1}^{n_c}{\beta }_c{d}_c{e}^{-\frac{t-T_{c,i}}{{\tau }_c}}\mathcal{H}\left(t-T_{c,i}\right),$$

$${\sigma }_t=\eta \Delta \sigma +S_h(1-\varphi )+(S_c-s)\varphi -\left({\gamma }_h(1-\varphi )+{\gamma }_c\varphi \right)\sigma ,$$

$$s(t)=\sum _{i=1}^{n_a}{\beta }_a{d}_a{e}^{-\frac{t-T_{a,i}}{{\tau }_a}}\mathcal{H}\left(t-T_{a,i}\right)$$

$$p_t=D\Delta p+{\alpha }_h(1-\varphi )+{\alpha }_c\varphi -{\gamma }_{p}p$$

$$P_s={\int }_{\Omega }pdx.$$

$$\frac{d{P}_s}{dt}={\alpha }_h{V}_h+{\alpha }_c{V}_c-{\gamma }_p{P}_s,$$

$$\phi_{t}-\lambda\Delta\phi+2\phi(1-\phi)f(\phi,\sigma,u)=0,\mbox{ in }{\color[rgb]{0,0,0}Q_{T}:={}}(0,T)\times\Omega,$$

$$\sigma_{t}-\eta\Delta\sigma+\gamma_{h}\sigma+(\gamma_{c}-\gamma_{h})\sigma\phi{}={\color[rgb]{0,0,0}S_{h}+(S_{c}-S_{h})\phi-s\phi,{}}\mbox{ in }Q_{T},$$

$$p_{t}-D\Delta p+\gamma_{p}p=\alpha_{h}+(\alpha_{c}-\alpha_{h})\phi,\mbox{ in }{\color[rgb]{0,0,0}Q_{T}},$$

$$\phi=0,\mbox{ }\frac{\partial\sigma}{\partial\nu}=\frac{\partial p}{\partial\nu}=0,\mbox{ on }{\color[rgb]{0,0,0}\Sigma_{T}:={}}(0,T)\times\partial\Omega,$$

$$\varphi (0)={\varphi }_0,\text{mbox}\sigma (0)={\sigma }_0,\text{mbox}p(0)=p_0,\text{mbox}in\Omega ,$$

$$V_0=H_0^1(\Omega ),\text{mbox}{V}_0^{\prime }=(H_0^1(\Omega ){)}^{\prime }=:H^{-1}(\Omega ),\text{mbox}V=H^1(\Omega ),\text{mbox}{V}^{\prime }=(H^1(\Omega ){)}^{\prime },\text{mbox}$$

$$W_0=H^2(\Omega )\cap H_0^1(\Omega ),\text{mbox}W=\left\{y\in H^2(\Omega );\text{mbox}\frac{\partial y}{\partial \nu }=0\text{mbox}on\partial \Omega \right\},$$

$$u\in L^{\infty}(Q_{T}),\mbox{ }\ {\color[rgb]{0,0,0}s\in L^{\infty}(Q_{T})},$$

$$(r_1,r_2,r_3)\to f(r_1,r_2,r_3)\text{mbox}isgloballyLipschitzcontinuouson{\mathbb{R}}^3.$$

$$\gamma_{ch}=\gamma_{c}-\gamma_{h},\quad\alpha_{ch}=\alpha_{c}-\alpha_{h},\quad{\color[rgb]{0,0,0}S_{ch}=S_{c}-S_{h}},$$

$$\varphi$$

$$\sigma$$

$$p$$

$${\int }_0^T{⟨{\varphi }_t(t),{\psi }_1(t)⟩}_{V_0^{\prime },V_0}dt+{\int }_{Q_T}\left\{\lambda \nabla \varphi \cdot \nabla {\psi }_1+2\varphi (1-\varphi )f(\varphi ,\sigma ,u){\psi }_1\right\}dxdt$$

$$+{\int }_0^T{⟨{\sigma }_t(t),{\psi }_2(t)⟩}_{V^{\prime },V}dt+{\int }_{Q_T}(\eta \nabla \sigma \cdot \nabla {\psi }_2+{\gamma }_h\sigma {\psi }_2+{\gamma }_{ch}\sigma \varphi {\psi }_2)dxdt$$

$$\displaystyle=\int_{Q_{T}}{\color[rgb]{0,0,0}(S_{h}\psi_{2}+(S_{ch}-s)\phi\psi_{2})}dxdt,\mbox{ for all }(\psi_{1},\psi_{2})\in L^{2}(0,T;V_{0}\times V),$$

$${\int }_0^T{⟨p_t(t),\psi (t)⟩}_{V^{\prime },V}dt+{\int }_{Q_T}(D\nabla p\cdot \nabla \psi +{\gamma }_{p}p\psi )dxdt$$

$$={\int }_{Q_T}({\alpha }_h+{\alpha }_{ch}\varphi )\psi dxdt,\text{mbox}forall\psi \in L^2(0,T;V),$$

$$(\phi ,\sigma ,p)(0)=({\phi }_0,{\sigma }_0,p_0).$$

$$X_0$$

Full text

Mathematical analysis and simulation study of a phase-field model of prostate cancer growth with chemotherapy and antiangiogenic therapy effects

Pierluigi Colli

Dipartimento di Matematica, Università degli Studi di Pavia and IMATI-C.N.R.,

Via Ferrata 5, 27100 Pavia, Italy

pierluigi.colli@unipv.it

Hector Gomez

School of Mechanical Engineering, Purdue University,

516 Northwestern Avenue, West Lafayette, IN 47907, USA

and

Weldon School of Biomedical Engineering, Purdue University,

206 S. Martin Jischke Drive, West Lafayette, IN 47907, USA

and

Purdue Center for Cancer Research, Purdue University,

201 S. University Street, West Lafayette, IN 47907, USA

hectorgomez@purdue.edu

Guillermo Lorenzo

Dipartimento di Ingegneria Civile e Architettura,

Università degli Studi di Pavia and IMATI-C.N.R.,

Via Ferrata 3, 27100 Pavia, Italy

guillermo.lorenzo@unipv.it

Gabriela Marinoschi

“Gheorghe Mihoc-Caius Iacob” Institute of Mathematical Statistics

and Applied Mathematics of the Romanian Academy,

Calea 13 Septembrie 13, 050711 Bucharest, Romania

gabriela.marinoschi@acad.ro

Alessandro Reali

Dipartimento di Ingegneria Civile e Architettura,

Università degli Studi di Pavia and IMATI-C.N.R.,

Via Ferrata 3, 27100 Pavia, Italy

alessandro.reali@unipv.it

Elisabetta Rocca

Dipartimento di Matematica, Università degli Studi di Pavia and IMATI-C.N.R.,

Via Ferrata 5, 27100 Pavia, Italy

elisabetta.rocca@unipv.it

Abstract

Chemotherapy is a common treatment for advanced prostate cancer. The standard approach relies on cytotoxic drugs, which aim at inhibiting proliferation and promoting cell death. Advanced prostatic tumors are known to rely on angiogenesis, i.e., the growth of local microvasculature via chemical signaling produced by the tumor. Thus, several clinical studies have been investigating antiangiogenic therapy for advanced prostate cancer, either as monotherapy or in combination with standard cytotoxic protocols. However, the complex genetic alterations that originate and sustain prostate cancer growth complicate the selection of the best chemotherapeutic approach for each patient’s tumor. Here, we present a mathematical model of prostate cancer growth and chemotherapy that may enable physicians to test and design personalized chemotherapeutic protocols in silico. We use the phase-field method to describe tumor growth, which we assume to be driven by a generic nutrient following reaction-diffusion dynamics. Tumor proliferation and apoptosis (i.e., programmed cell death) can be parameterized with experimentally-determined values. Cytotoxic chemotherapy is included as a term downregulating tumor net proliferation, while antiangiogenic therapy is modeled as a reduction in intratumoral nutrient supply. An additional equation couples the tumor phase field with the production of prostate-specific antigen, which is a prostate cancer biomarker that is extensively used in the clinical management of the disease. We prove the well-posedness of our model and we run a series of representative simulations leveraging an isogeometric method to explore untreated tumor growth as well as the effects of cytotoxic chemotherapy and antiangiogenic therapy, both alone and combined. Our simulations show that our model captures the growth morphologies of prostate cancer as well as common outcomes of cytotoxic and antiangiogenic mono and combined therapy. Additionally, our model also reproduces the usual temporal trends in tumor volume and prostate-specific antigen evolution observed in experimental and clinical studies.

Keywords: prostate cancer; computational oncology; phase field; semilinear parabolic equations; well-posedness; isogeometric analysis

**AMS Subject Classification: 35Q92, 92C50, 65M60, 35K51, 35K58 **

1 Introduction

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 [55]. 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 [55]. Advanced PCa patients are usually prescribed androgen-deprivation therapy or chemotherapy [55], 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 [74] and an increase in microvascular density is associated to poorer prognosis [53]. 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 [21]) did not show significant benefit in castration-resistant PCa when administered alone [63], but it produced a 50% PSA decline in 75% of the patients when combined with chemotherapy [61]. 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 [35] 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 [38]. 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 $\phi$ and transitions from the value $\phi\approx 0$ in the host tissue to $\phi\approx 1$ in the tumor. The transition is smooth but steep and takes on a hyperbolic tangent profile in the direction perpendicular to the interface [29]. The model also accounts for the dynamics of a critical nutrient, whose concentration is denoted by $\sigma$ and obeys a reaction-diffusion equation. The concentration of PSA in the prostatic tissue $p$ 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 [77]

[TABLE]

where $\phi_{t}$ denotes the partial derivative of $\phi$ with respect to time, $\Delta$ is the Laplace operator, $\lambda$ is the diffusion coefficient of tumor cells and

[TABLE]

The diffusion coefficient of tumor cells can be computed as $\lambda=M\ell^{2}$, where both $M$ and $\ell$ are positive real constants denoting the tumor mobility and interface width, respectively [77]. Here, $F(\phi)=M\phi^{2}(1-\phi)^{2}$ is a double-well potential, i.e., a nonconvex function, typical in phase-field modeling, which allows the coexistence of the tumoral ($\phi\approx 1$) and healthy ($\phi\approx 0$) tissue. The function $h(\phi)=M\phi^{2}(3-2\phi)$ is also common in non-conserved dynamics of phase fields. It is usually called interpolation function because it verifies the properties $h(0)=0$ and $h(1)=1$; another important property of $h$ is that $h^{\prime}(0)=h^{\prime}(1)=0$, where $h^{\prime}$ denotes the derivative of $h$. The function $m(\sigma)$ is normally called tilting function and, in our model, is defined as

[TABLE]

where $\rho$ and $A$ represent constant proliferation and apoptosis indices, respectively. These nondimensional parameters are related to the proliferation and apoptotis rates in tumoral tissue as follows:

[TABLE]

and

[TABLE]

where $K_{\rho}$ is the proliferation rate of tumor cells, $K_{A}$ is the apoptosis rate of tumor cells, $\bar{K_{\rho}}$ is a scaling reference value for the proliferation rate, and $\bar{K_{A}}$ is a scaling reference value for the apoptosis rate. Therefore, $m(\sigma)$ can be interpreted as a function describing the tumor net proliferation rate. The positive constant $m_{ref}$ scales the strength of the tilting function within our phase-field framework. The constants $\sigma_{r}$ and $\sigma_{l}$ are, respectively, a reference and a threshold value for the nutrient concentration [77]. For nutrient concentrations lower than $\sigma_{l}$ healthy tissue is energetically more favorable than tumor tissue and viceversa. The function $u$ in (2.2) represents the tumor-inhibiting effect of a cytotoxic drug and is described in Section 2.2.2. When $|m(\sigma)-m_{ref}u|<1/3$, the function $G$ is a double-well potential with local minima at $\phi=0$ and $\phi=1$. Within this range, low values of the nutrient concentration (or large values of $u$) produce a lower energy level (value of $G$) in the healthy tissue ($\phi=0$) than in the tumoral tissue ($\phi=1$). The opposite is true for high values of the nutrient concentration (or low values of $u$). The described behavior of $G$ is further illustrated in Figure 1. From (2.2), we can obtain

[TABLE]

where

[TABLE]

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

[TABLE]

where $n_{c}$ is the number of chemotherapy cycles, $\beta_{c}$ measures the cytotoxic effect of the treatment on tumor dynamics per unit of drug dose delivered, $d_{c}$ is the drug dose, $T_{c,i}$ with $i=1,\dots,n_{c}$ are the times of drug delivery, $\tau_{c}$ is the mean lifetime of the chemotherapeutic drug, and $\mathcal{H}$ denotes the Heaviside function. These parameters are required to be strictly positive. The experimental and clinical values observed for $\beta_{c}$ ensure that the condition $|m(\sigma)-m_{ref}u|<1/3$ is verified; see Section 2.2.1. According to (2.2), $u(t)$ 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 $\tau_{c}$, $\beta_{c}$ and $d_{c}$ accordingly.

2.2.3 Nutrient dynamics

Nutrient dynamics is controlled by the equation

[TABLE]

where $\eta$ is the diffusion coefficient of the nutrient, $S_{h}$ is the nutrient supply rate in the healthy tissue, $S_{c}$ stands for the nutrient supply rate in the cancerous tissue, $s$ is a given function yielding the reduction in nutrient supply caused by antiangiogenic therapy (see Section 2.2.4), and $\gamma_{h}$, $\gamma_{c}$ are positive constants that represent the nutrient uptake rate in the healthy and cancerous tissue, respectively. We require that $S_{h}$, $S_{c}$, and $s$ are non-negative and $s$ should satisfy the constraint $s\leq S_{c}$.

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):

[TABLE]

where $n_{a}$ is the number of antiangiogenic treatment cycles, $\beta_{a}$ measures the antiangiogenic effect of the treatment on the nutrient supply per unit of drug dose delivered, $d_{a}$ is the drug dose, $T_{a,i}$ are the times of drug delivery ($i=1,\dots,n_{a}$), $\tau_{a}$ stands for the mean lifetime of the antiangiogenic drug, and $\mathcal{H}$ 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 $\tau_{a}$, $\beta_{a}$ and $d_{a}$.

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 $\gamma_{p}$. Hence, we propose the equation [51]

[TABLE]

where $D$ is the diffusion constant. The constants $\alpha_{h}$ and $\alpha_{c}$ represent, respectively, the tissue PSA production rate of healthy and malignant cells. The serum PSA $P_{s}$ can be defined as the integral of the tissue PSA over the prostate $\Omega$, that is,

[TABLE]

By using (2.11), one can show that for free-flux boundary conditions

[TABLE]

where $V_{h}$ is the volume of prostatic healthy tissue and $V_{c}$ 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 [69].

3 Well-posedness of the model

In this section we shall provide an existence and uniqueness result for this system

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $f$ is given by (2.7) and (2.3); $M,m_{ref},\rho,A,\sigma_{l},\sigma_{r}$ are given real constants; and $\lambda,$ $\eta,$ $D,$ $\gamma_{c},$ $\gamma_{h},$ $\gamma_{p},$ $\alpha_{c},$ $\alpha_{h}$, $S_{c}$, $S_{h}$ are positive fixed real constants.

Here, $\Omega$ is an open bounded subset of $\mathbb{R}^{N}$, $N\leq 3$, with a sufficiently smooth boundary $\partial\Omega,$ $\frac{\partial}{\partial\nu}$ represents the outward normal derivative to $\partial\Omega$ and $T$ denotes a final time.

3.1 Functional framework

The problem will be treated in a functional framework involving the space $H=L^{2}(\Omega),$ which is identified with its dual space $H^{\prime}\cong H,$ and the Sobolev spaces

[TABLE]

[TABLE]

and with the dense and compact injections $W_{0}\subset V_{0}\subset H\subset V_{0}^{\prime}$ and $W\subset V\subset H\subset V^{\prime}.$ $H^{1}_{0}(\Omega)$ contains the elements of $H^{1}(\Omega)$ with null trace on the boundary $\partial\Omega$. For the norms in these spaces we will use the notation $\|\,\cdot\,\|_{B}$, where $B$ is the space we are considering.

For $q\in[1,\infty]$ and $z\in L^{q}(\Omega)$ or $z\in L^{q}(Q_{T})$, we simply denote the norm of $z$ in these $L^{q}$ spaces by $\left\|z\right\|_{q}.$

We will also make use of spaces of functions that depend on time with values in a Banach space $B$. Namely, for $q\in[1,\infty]$ we consider the space $L^{q}(0,T;B)$ of measurable functions $t\mapsto z(t)$ such that $t\mapsto\|z(t)\|_{B}^{q}$ is integrable on $(0,T)$ (or essentially bounded if $q=\infty$) 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 isomorph to $L^{2}(Q_{T})$. Moreover, for $q\in[1,\infty]$, $W^{1,q}(0,T;B)$ will denote the space of functions $t\mapsto z(t)$ such that both $z$ and its (weak) derivative $z_{t}$ belong to $L^{q}(0,T;B)$. We point out that $W^{1,q}(0,T;B)\subset C([0,T];B)$ for all $q\in[1,\infty]$.

**Hypotheses. **Following the considerations presented in the Introduction and extending the particular choices of $u$ and $s$ in (2.8) and (2.10), we let

[TABLE]

and it can be easily checked that

[TABLE]

We also set

[TABLE]

and point out that the problem parameters $\lambda,$ $\eta,$ $D,$ $\gamma_{c},$ $\gamma_{h},$ $\gamma_{p},$ $\alpha_{c},$ $\alpha_{h}$, $S_{c}$, $S_{h}$ are positive.

In the sequel, by $C$ we shall denote a constant, that may change from line to line, depending on the problem parameters, the domain $\Omega,$ the final time $T,$ the norms of the initial data and possibly on the norms of $u$ and $s$. Moreover, we assume that (3.6)-(3.8) hold.

Definition 3.1. Let $(\varphi_{0},\sigma_{0},p_{0})\in H\times H\times H.$ A solution to the system (3.1)-(3.5) is a triplet $(\phi,\sigma,p)$, with

[TABLE]

which satisfies

[TABLE]

[TABLE]

and

[TABLE]

Please note that taking in (3.10) $\psi_{2}$ first, and $\psi_{1}$ then, equal to the null function (and this is of course a suitable choice in both cases), we obtain two separate variational equalities for $\phi$ and $\sigma$, as it is (3.11) for $p$.

Let us denote

[TABLE]

**Theorem 3.2. ***Let * **

[TABLE]

[TABLE]

Then, the system (3.1)-(3.5) *has a unique solution $(\phi,\sigma,p)$ in the sense of *Definition 3.1, such that

[TABLE]

If $(\sigma_{0},p_{0})\in L^{\infty}(\Omega)\times L^{\infty}(\Omega),$ then $(\sigma,p)\in L^{\infty}(Q_{T})\times L^{\infty}(Q_{T}).$ Moreover, if

[TABLE]

then we have

[TABLE]

Finally, if $(\phi_{0},\sigma_{0},p_{0})\in V_{0}\times V\times V,$ the solution has the supplementary regularity $(\phi,\sigma,p)\in X_{0}\times X\times X,$ and satisfies the estimate

[TABLE]

In addition, the solution is continuous with respect to the data, that is, for two solutions $(\phi_{i},\sigma_{i},p_{i})$ corresponding to the data (\phi_{0}^{i},\sigma_{0}^{i},p_{0}^{i},u_{i},{\color[rgb]{0,0,0}s}_{i}),* $i=1,2,$ we have*

[TABLE]

for all $t\in[0,T].$

**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

[TABLE]

which is a complete metric space provided we take the distance induced by some norm in $C([0,T];H)$.

Let us fix $z\in\mathcal{M}$ in equations (3.1), (3.2):

[TABLE]

[TABLE]

where $\sigma^{z}$ in (3.20) is the solution to (3.21), corresponding to $z$, with $\frac{\partial\sigma}{\partial\nu}=0$ on $\Sigma_{T}$ and $\sigma(0)=\sigma_{0}$ in $\Omega.$ Thus, we can define the mapping

[TABLE]

where $\phi^{z}$ is the solution to (3.20), corresponding to $z$ and $\sigma^{z}$, with the homogeneous Dirichlet boundary condition and the initial datum $\phi_{0}.$ We shall prove that $\Psi(\mathcal{M})\subset\mathcal{M}$ and that $\Psi$ is a contraction mapping on $\mathcal{M}$.

First, we treat the initial-boundary value problem for equation (3.21). For all $t\in[0,T]$ we introduce the operator $A(t):V\rightarrow V^{\prime}$ defined by

[TABLE]

and observe that it is continuous

[TABLE]

and quasi-monotone from $V$ to $V^{\prime}$, i.e.,

[TABLE]

Due to a general solvability result for linear parabolic problems (the reader may consult, e.g., Ref.  [48]), for $\sigma_{0}\in H$ and {\color[rgb]{0,0,0}s}\in L^{\infty}(Q_{T}) the initial-boundary value problem for (3.21) has a unique solution

[TABLE]

A first estimate performed by testing (3.21) by $\sigma^{z}$ and integrating over $(0,t)$ yields, after a standard calculation,

[TABLE]

for all $t\in[0,T],$ where we used the fact that $z\in\mathcal{M}.$ As specified before, $C$ is a constant, that may vary from line to line, depending on the problem parameters and $T$: indeed, the Gronwall lemma has been applied to derive (3.22).

Now, if we take two solutions $\sigma_{1}$ and $\sigma_{2}$ corresponding to (z_{1},{\color[rgb]{0,0,0}s_{1}},\sigma_{0}^{1}) and (z_{2},{\color[rgb]{0,0,0}s_{2}},\sigma_{0}^{2}), thanks to the Hölder and Young inequalities we obtain for the differences $\sigma:=\sigma_{1}-\sigma_{2}$ and $z:=z_{1}-z_{2}:$

[TABLE]

By the Sobolev embedding theorems, it turns out that $V\subset L^{4}(\Omega)$ with continuous embedding. Then, by exploiting the Young inequality and (3.22) for $\sigma_{1}$ we infer that

[TABLE]

Going back to (3.23) we deduce that

[TABLE]

where C_{01}:=C\left({\color[rgb]{0,0,0}1+{}}\left\|\sigma_{0}^{1}\right\|_{H}^{2}+\left\|{\color[rgb]{0,0,0}s}_{1}\right\|^{2}_{\infty}\right)\!. Finally, by using Gronwall lemma it is not difficult to conclude that

[TABLE]

where $C$ stands for another constant depending additionally on the initial data, $T$, and the source terms.

Next, we treat the initial-boundary value problem for the equation (3.20). To this end, for $z\in\mathcal{M}$ and $\sigma^{z}$ defined as above we consider the intermediate equation

[TABLE]

with

[TABLE]

where

[TABLE]

We observe that $g$ is bounded and Lipschitz continuous on $\mathbb{R}$. Recalling (2.7) we note that $f(z,\sigma^{z},u)\in L^{\infty}(Q_{T})$ and

[TABLE]

We are going to prove that (3.25)-(3.26) has a solution. We apply the Banach fixed point theorem, by fixing $\zeta\in C([0,T];H)$ and studying the problem

[TABLE]

together with (3.26). We note that the right-hand side of (3.29) is in $L^{\infty}(Q_{T})$ and that $\phi_{0}\in H.$ It is clear, again by the solvability result for linear parabolic problems (see Ref.  [48]), that (3.29), (3.26) has a unique solution

[TABLE]

Now, we aim at proving that the operator associating $\zeta$ to the solution $\phi$ to (3.29) and (3.26) is a contraction mapping on $C([0,T];H).$ For the sake of convenience, we introduce the norm

[TABLE]

which is equivalent to the standard norm in $C([0,T];H).$

Let us consider two solutions $\phi_{1}$ and $\phi_{2}$ to (3.29), (3.26) corresponding to $\zeta_{1}$ and $\zeta_{2},$ respectively. We multiply the difference of the equations (3.29) by $\phi_{1}-\phi_{2}$ and integrate over $\Omega\times(0,t).$ With the help of (3.26) and (3.28) we have that

[TABLE]

Hence, the Young inequality and the Gronwall lemma allow us to get

[TABLE]

We multiply this inequality by $e^{-2\gamma t}$ and have successively

[TABLE]

Taking the supremum with respect to $t\in[0,T]$, we obtain

[TABLE]

and note that for $\gamma>\frac{C}{2}$ the mapping $\zeta\mapsto\phi$ is a contraction on $C([0,T];H)$. Thus, we deduce that (3.25)-(3.26) has a unique solution.

Next, we recall (3.15), that is $0\leq\phi_{0}\leq 1$ a.e. in $\Omega.$ Testing (3.25) by $-\phi^{-}$ $(\phi^{-}$ being the negative part of $\phi)$, by a few calculations we obtain

[TABLE]

since $\phi_{0}^{-}=0$ a.e. in $\Omega$, $\phi^{-}=0$ a.e. in the set where $\phi\in[0,1]$ and $g(\phi)=0$ a.e. in the set where $\phi\notin[0,1].$ This implies that $\phi^{-}(t)=0,$ that is $\phi(t)\geq 0$ a.e. in $\Omega$, for all $t\in[0,T].$

Next, we multiply (3.25) by $(\phi-1)^{+}$ and by similar calculations we deduce that

[TABLE]

Here, we used the fact that $\phi=0$ a.e. on $\Sigma_{T},$ so that $(\phi-1)^{+}=0$ a.e. on $\Sigma_{T}.$ Thus, we have shown that $\phi(t)\leq 1$ a.e. in $\Omega$, for all $t\in[0,T].$

By these two last results we actually proved that $0\leq$ $\phi\leq 1$ a.e. on $Q_{T},$ implying that $g(\phi)=-2\phi(1-\phi).$ It turns out that this solution $\phi$ actually solves equation (3.20). Moreover, we have that $\phi$ satisfies (cf. (3.9))

[TABLE]

Since the problem (3.20), (3.26) may also have other solutions fulfilling (3.30), we have to prove a uniqueness result. For that, let us consider two solutions $\phi_{1}$ and $\phi_{2}$ corresponding to the same data, and compute the difference of the respective equations (3.20) tested by $\phi_{1}-\phi_{2}.$ Owing to a straightforward calculation and (3.30), we deduce that

[TABLE]

which by Gronwall lemma leads to the desired uniqueness property.

The solution $\phi$ to (3.20), (3.26) we have found is actually what we have denoted by $\phi^{z}=\Psi(z).$ We have already shown that $\phi^{z}\in C([0,T];H)$ and that $0\leq\phi^{z}\leq 1$ a.e. in $Q_{T}$, whence it follows that $\Psi(\mathcal{M})\subset\mathcal{M}.$ Thus, for the last step it suffices to check that $\Psi$ is a contraction mapping on $C([0,T];H)$.

First, we shall deduce a general estimate.

Let us consider two pairs of data (\phi_{0}^{i},\sigma_{0}^{i},u_{i},{\color[rgb]{0,0,0}s_{i}};z_{i}), $i=1,2,$ and the corresponding solutions $(\phi_{i}=\phi^{z_{i}},\sigma_{i}=\sigma^{z_{i}}),$ $i=1,2.$ Then, we take the difference of the equations (3.20), add $\lambda(\phi_{1}-\phi_{2})$ to both sides, then multiply by $\phi_{1}-\phi_{2}$ and integrate over $\Omega\times(0,t).$ Here, let us denote $\phi:=\phi_{1}-\phi_{2},$ $\phi_{0}:=\phi_{0}^{1}-\phi_{0}^{2},$ and take into account that $f$ is Lipschitz continuous with a Lipschitz constant $L_{f}$ depending on $M,$ $m_{ref},$ $\rho,$ $A,$ $\sigma_{r}$. We have

[TABLE]

Therefore, recalling (3.28) and using the Young inequality, we infer that

[TABLE]

where $C$ depends on $\left\|u_{1}\right\|_{\infty}$ as well. Next, we use (3.1) and obtain

[TABLE]

Using the Gronwall lemma and estimating the right-hand side, we deduce that

[TABLE]

for all $0\leq s<t\leq T.$ Taking now the supremum with respect to $0\leq s\leq t$, we have that

[TABLE]

for all $t\in[0,T].$ At this point, if we choose the same initial data and the same $u_{i}$ and $s_{i}$, $i=1,2$, we obtain

[TABLE]

Hence, in the light of the argument used before, by considering here the norm

[TABLE]

and playing with (3.32), we plainly obtain

[TABLE]

Then, if $\gamma>\frac{C}{2},$ the mapping $z\mapsto\phi^{z}=\Psi(z)$ is a contraction on $\mathcal{M}$ and $\phi^{z}=z$ is the unique fixed point. Thus, problem (3.1)-(3.2) with the initial and boundary conditions has a unique solution, with all properties inherited from the intermediate solutions occurred in the proof.

Now, we can replace $z_{i}$ by $\phi_{i}$, $i=1,2$, in (3.31) and apply the Gronwall lemma to infer that

[TABLE]

Combining this estimate with (3.1) finally leads to

[TABLE]

Let us now pass to (3.3) in which the right-hand side is in $L^{\infty}(Q_{T}).$ It is easily seen, relying on the arguments before for $\sigma$, that the probem based on (3.3), with the initial datum $p_{0}\in H$ and the homogeneous Neumann boundary condition, has a unique solution

[TABLE]

Moreover, for the difference of two solutions corresponding to different data we have the estimate

[TABLE]

By combining (3.1) and (3.34), we obtain the continuity of the solution with respect to the data, as claimed by (3.1).

Next, we prove the boundedness of $\sigma$ and $p$ and the property (3.17). We rewrite eq. (3.2) as

[TABLE]

and consider $\sigma_{0}\in L^{\infty}(\Omega).$ By Theorem 7.1, p. 181, in Ref.  [46], we find that

[TABLE]

Let $\sigma_{0}(x)\geq 0,$ $s(t,x)\leq S_{c}$ a.e. and multiply (3.2) by $-\sigma^{-}.$ We have

[TABLE]

since $0\leq\phi\leq 1$ and {\color[rgb]{0,0,0}(S_{c}-s)}\geq 0 a.e. in $Q_{T}.$ By Gronwall lemma we obtain that $\sigma^{-}(t)=0,$ hence $\sigma(t)\geq 0$ a.e. in $\Omega,$ for all $t\in[0,T].$

By a similar argument we deduce that $p\in L^{\infty}(Q_{T})$ and $p(t)\geq 0$ a.e. in $\Omega$, for all $t\in[0,T].$

Finally, we prove the further regularity of the solution and the estimate (3.1). The properties $\phi\in X_{0}$, $\sigma\in X$, $p\in X$ can be compared with (3.9) and we can observe that, e.g., for $\sigma$ and $p$ the triplet of spaces $(H,V,W)$ now replaces $(V^{\prime},H,V)$. This fits into the frame of the structural regularity of solutions to parabolic systems (see, e.g., Ref.  [48]).

Let $\phi_{0}\in V_{0}$ and write (3.1) as

[TABLE]

Since the operator $B_{0}:W_{0}\subset H\rightarrow H,$ $B_{0}=-\lambda\Delta$ is potential (that is, $B_{0}\phi=\partial j(\phi),$ $j:V_{0}\rightarrow[0,+\infty],$ $j(\phi)=\frac{\lambda}{2}\int_{\Omega}\left|\nabla\phi\right|^{2}dx)$, according e.g. to Theorem 4.18, p. 179, in Ref.  [4], it follows that

[TABLE]

Then, we can test (3.1) by $\phi_{t}$; taking into account that (cf. (3.1))

[TABLE]

we have

[TABLE]

since $\left\|\phi\right\|_{\infty}\leq 1.$ This yields

[TABLE]

By comparison in the equation for $\phi$ we realize that $\Delta\phi\in L^{2}(0,T;H)$. Then, using the elliptic regularity theory, due to the homogeneous Dirichlet boundary condition for $\phi$, we can conclude that

[TABLE]

Thus, we have obtained $\phi\in X_{0}$.

Assume that $\sigma_{0}\in V$ and consider again the equation (3.2) written as

[TABLE]

The operator $-\eta\Delta+\gamma_{h}I:W\subset H\to H$ ($I$ standing for the identity in $H$) is potential as well, and then $\sigma\in W^{1,2}(0,T;H).$ Multiplying (3.2) by $\sigma_{t}$ and integrating over $\Omega\times(0,t)$, we easily arrive at

[TABLE]

for all $t\in[0,T]$. Also, from a comparison in (3.38) it follows that $\Delta\sigma\in L^{2}(0,T;H)$. This property, together with the homogeneous Neumann boundary condition for $\sigma$, entails

[TABLE]

A similar argumentation for $p$ yields

[TABLE]

By collecting (3.36), (3.37), (3.39)-(3.41) it is easy to get (3.1), as claimed. Let us add the comment that, under the additional regularity specified by the estimate (3.1), the triplet $(\phi,\sigma,p)$ turns out to be a strong solution to the system (3.2)-(3.6), that is, it actually satisfies equations and boundary conditions almost everywhere in $Q_{T}$ and on $\Sigma_{T}$, respectively. By this remark, we conclude our proof. $\square$

4 Numerical method

4.1 Space discretization

We begin by deriving a weak form of the governing equations, in preparation for an isogeometric discretization based on the Galerkin method. The weak form is derived multiplying the governing equations with suitable smooth functions, integrating over the domain, and performing integration by parts. We recall that the spaces $V=H^{1}(\Omega)$ and $V_{0}=H^{1}_{0}(\Omega)$, defined in Section 3. The weak form can be stated as follows: find $\phi\in{V}_{0}$, $\sigma\in{V}$ and $p\in{V}$ such that for all $\{w_{1},w_{2},w_{3}\}\in V_{0}\times V\times V$,

[TABLE]

Once the equations have been written in weak form, they can be discretized by restricting the weak form to finite dimensional spaces ${V}^{h}\subset{V}$ and ${V}_{0}^{h}\subset{V}_{0}$. In this work, we construct these spaces using quadratic B-spline basis functions [38]. As an example, we take ${V}^{h}=\hbox{\rm span}\{N_{j}\}_{j=1,\dots,n_{b}}$ where the $N_{j}$’s are multivariate B-splines. The function $\phi^{h}$, which is a finite dimensional approximation to $\phi$, is defined as $\phi^{h}(t,{x})=\sum_{j=1}^{n_{b}}\phi_{j}(t)N_{j}({x})$ where the $\phi_{j}$’s are referred to as control variables. The functions $\sigma^{h}$, $p^{h}$, $w_{i}^{h}$ for $i=1,2,3$ are defined analogously. Because $\phi^{h}\in V^{h}_{0}$, the control variables $\phi_{j}$ corresponding to the control points $j$ on the boundary must be zero to satisfy the homogeneous Dirichlet boundary conditions.

4.2 Time discretization

We integrate in time using the generalized-$\alpha$ method [13, 42]. To define our time integration scheme, we introduce the global vector of degrees of freedom $\boldsymbol{\Phi}=\{\phi_{j}\}_{j=1,\dots,n_{b}}$. Analogously, we define $\boldsymbol{\Sigma}$ and $\boldsymbol{P}$ as the global vectors of control variables for $\sigma^{h}$ and $p^{h}$, respectively. Let us introduce the following residual vectors

[TABLE]

where

[TABLE]

Let us call $\boldsymbol{U}_{n}=\{\boldsymbol{\Phi}_{n},\boldsymbol{\Sigma}_{n},\boldsymbol{P}_{n}\}$ the time discrete approximation to the control variables at time $t_{n}$. Our time-integration algorithm can be defined as follows: given $\boldsymbol{U}_{n}$, $\dot{\boldsymbol{U}}_{n}$ and the current time step $\Delta t_{n}=t_{n+1}-t_{n}$, find $\boldsymbol{U}_{n+1}$, $\dot{\boldsymbol{U}}_{n+1}$ such that

[TABLE]

The parameters $\alpha_{m}$, $\alpha_{f}$ and $\gamma$ define the stability and accuracy of the algorithm. The generalized-$\alpha$ method can be proven to be second-order accurate and $A$-stable by taking $\rho_{\infty}\in[0,1]$ and

[TABLE]

In our calculations, we take $\rho_{\infty}=1/2$. The resulting nonlinear system of equations is linearized using the Newton-Raphson algorithm. We converged all the individual residuals $\boldsymbol{R}_{\phi}$, $\boldsymbol{R}_{\sigma}$, and $\boldsymbol{R}_{p}$ up to a certain tolerance $\varepsilon_{NR}$. This strategy ensured that the fulfilment of the convergence criterion by the global residual $\boldsymbol{R}$ was not due to potentially different scaling of the equations of the system (4.5)-(4.7). Within each iteration of the Newton-Raphson algorithm, we solved the corresponding linearized system by means of the GMRES [65] algorithm with a diagonal preconditioner. We chose a certain tolerance $\varepsilon_{GMRES}$ and a maximum number of iterations as convergence criteria for the GMRES algorithm.

5 Simulation study

5.1 Description and parameter selection

We consider two different cases of PCa in this simulation study: a mild tumor and an aggressive tumor. These different tumor behaviors can be implemented by appropriately choosing $\rho$ and $A$ within the function $m(\sigma)$ in (2.3). In general, tumor growth is usually characterized by a high proliferation index and a low apoptotis index [34], which would translate into choosing $|\rho|>|A|$ within the context of our model. Then, we may represent mild tumors with lower values of these parameters for which the difference $|\rho|-|A|$ is also low, while more aggressive tumors can be modeled with higher values and such that the difference $|\rho|-|A|$ is also larger. For the simulations in this study, we calibrated $\rho$ and $A$ using the proliferation and apoptosis rates previously reported in the literature [7, 66] in equations (2.4)-(2.5). We further computationally calibrated $m_{ref}$ in (2.2)-(2.3) to ensure $|m(\sigma)-m_{ref}u|<1/3$ (see Section 2.2.1).

We begin our simulation study by analyzing mild and aggressive PCa growth in an untreated scenario. The metabolic profile of PCa is known to vary as it progresses towards a more malignant disease [34, 71]. Thus, we further consider four additional scenarios of untreated mild and aggressive PCa growth. We explore the effects of an effective and a poor nutrient supply within the tumor, i.e., high and low $S_{c}$ in (2.9), respectively. Moreover, we analyze the effect of tumor metabolism on its dynamics by choosing a larger or a smaller value of $\gamma_{c}$ in (2.9).

After studying untreated tumor growth, we proceed to analyze the effects of cytotoxic and antiangiogenic therapy alone and combined for both the mild and the aggressive tumor. As introduced in Section 2.2.2, we adopt the docetaxel-based cytotoxic chemotherapeutic protocol that is commonly used in the clinical management of advanced PCa. Hence, we consider 10 equal doses $d_{c}=$ 75 mg/m2 of docetaxel delivered every three weeks [55, 43, 18]. We set $\tau_{c}=5$ days, according to pharmacodynamic studies of this drug [3, 70]. We choose bevacizumab as antiangiogenic drug and assume that the antiangiogenic therapy also consists of 10 equal doses $d_{a}=15$ mg/kg every three weeks [2, 43, 61]. We set $\tau_{a}=30$ days [30, 52, 21]. We computationally calibrated the value of the constants $\beta_{c}$ and $\beta_{a}$ to produce a meaningful chemotherapeutic response as observed in clinical literature.

Table 1 summarizes the parameter selection discussed so far in this Section. The values of the other parameters in the model have been previously justified in the literature [51, 50, 49, 77].

5.2 Computational setup

All simulations were run on a 2D square domain with side length of $L_{d}=3000$ $\mathrm{\SIUnitSymbolMicro m}$ and 256 isogeometric elements per side. We chose a constant time step of $\Delta t=0.1$ days in all simulations. The convergence of the Newton-Raphson method was set at tolerance $\varepsilon_{NR}=10^{-3}$, while for the GMRES algorithm was set at $\varepsilon_{GMRES}=10^{-3}$ or a maximum of 500 iterations.

The initial tumor phase field is approximated as an ellipsoidal tumor placed in the center of the domain and semiaxes $a=150$ $\mathrm{\SIUnitSymbolMicro m}$ and $b=200$ $\mathrm{\SIUnitSymbolMicro m}$ parallel to the sides of domain. To implement this initial condition, we $L^{2}$-projected the hyperbolic tangent function

[TABLE]

over the quadratic B-spline space supporting our spatial discretization. This operation provides the control variables $\phi_{0,A}$ of the spline representation of the phase-field initial condition, i.e., $\phi^{h}_{0}(\boldsymbol{x},t)=\sum_{A=1}^{n_{b}}\phi_{0,A}(t)N_{A}(\boldsymbol{x})$ (see Section 4.1).

The initial conditions for the nutrient and the tissue PSA are approximations based on $\phi_{0}$ given by

[TABLE]

and

[TABLE]

The constants $c^{0}_{\sigma}$, $c^{1}_{\sigma}$, $c^{0}_{p}$, and $c^{1}_{p}$ are computationally calibrated [50], such that $\sigma_{0}$ and $p_{0}$ reproduce a constant value of the nutrient and tissue PSA within the tumor and the host tissue. Hence, we choose $c^{0}_{\sigma}=1$ g/L, $c^{1}_{\sigma}=-0.8$ g/L, $c^{0}_{p}=0.0625$ ng/mL/cc, and $c^{1}_{p}=0.7975$ ng/mL/cc.

In the simulations with the chemotherapeutic treatments, we initially let the tumor grow until $t=60$ days, when we deliver the first drug dose. Then, we apply the remaining 9 cycles of chemotherapy included in all treatment strategies (i.e., cytotoxic, antiangiogenic, and combined) until $t=249$ days, when the last drug dose is delivered. Finally, we let the simulation proceed until completing 1 year, i.e., $t=365$ days.

5.3 Simulation of untreated prostate cancer growth

Figures 2 and 3 respectively depict the evolution of the mild tumor and the corresponding nutrient distribution in all untreated simulation scenarios considered in this study. We observe that the mild tumor grows at a slow pace and maintaining a round, spheroidal morphology that matches previous studies in the literature [35, 19, 51, 77]. An increase in intratumoral nutrient supply or a metabolic adjustment decreasing the tumor nutrient uptake accelerates growth but does not alter the spheroidal morphology (see Figure 2B,E). Conversely, a more scarce nutrient supply or a higher metabolic dependence on the nutrient arrests tumor growth and progressively lead to its extinction (see Figure 2C,D). These results also align with the metabolic profile of mild prostatic tumors [34, 71], which have not acquired a more competitive metabolism to support a highly proliferative behavior and still show a considerable dependence on the baseline nutrient in the host tissue.

The growth of the aggressive tumor and the corresponding nutrient distribution in all untreated simulation scenarios are depicted in Figures 4 and 5, respectively. The aggressive tumor initially grows with a round morphology, but it eventually undergoes a morphological transformation by which it develops multiple branches and hence becomes more invasive [35, 19, 58]. This phenomenon has also been observed in previous studies [51, 28, 76]. Our simulations suggest that this shape instability is a tumor response to escape starvation because branching facilitates the access to the nutrient and leads to a nutrient redistribution that limits minimal concentrations in the inner regions of the tumor, which would potentially hamper PCa growth (see Figure 5). Hence, branching may be regarded as a mechanism of aggressive tumors to adapt to their harsh local environment. The finger-like structures are thicker when we increase the intratumoral nutrient supply or reduce the metabolic dependence of the tumor on the considered nutrient (see Figure 4B,E), whereas they become thinner in the opposite scenarios (see Figure 4C,D). We have also observed that the morphological shift takes place sooner when nutrient availability is scarce or the tumor’s metabolic dependence on the nutrient is larger. Conversely, the simulations with higher nutrient supply within the tumor and lower tumor nutrient consumption rate show a later appearance of branching.

Figure 6 displays the plots of tumor volume and serum PSA during each simulation of untreated mild PCa growth, while Figure 7 depicts the corresponding results for the aggressive cancer case. In general, aggressive tumors grew faster and reached larger volumes in comparison to the corresponding simulation cases for the mild tumor, as shown in Figure 7. In the reference case of untreated aggressive PCa, the tumor volume increase slightly slowed down by the end (see Figure 7A). This trend was also observed in the first half of the simulations with decreased intratumoral nutrient supply and higher nutrient consumption by the tumor (see Figure 7C,D). However, in these simulations we observe that the evolution of tumor volume exhibits a transient deceleration during branching. This marks an inflexion point that is later followed by a more rapid increase in tumor volume by the end of the simulation. This faster dynamics corresponds to the thriving growth of the newly-formed branches with a more favorable nutrient distribution. Hence, these observations suggest that branching was not completed in the reference simulation for the aggressive tumor, and the faster growing trend in tumor volume would ensue during the second year of simulation.

Figures 6 and 7 show that serum PSA followed a similar time evolution to tumor volume in all simulations of untreated PCa growth. This result is consistent with the extensive clinical use of this biomarker as a surrogate for the patient’s tumor burden [55, 51, 69, 73]. However, we observe that there were two cases in which tumor volume dynamics were noticeably faster than serum PSA dynamics, which correspond to the simulations featuring the aggressive tumor with increased intratumoral nutrient supply and lower nutrient uptake by cancerous tissue (see Figure 7B,E). Theses cases exhibited the fastest tumor growth dynamics as well as the largest tumor volumes registered in our simulation study and rapidly growing tumors have been previously reported to render comparatively lower serum PSA values [69].

5.4 Simulation of treatment plans

All treatment plans successfully removed mild PCa. Figure 8 depicts the evolution of tumor volume and serum PSA under cytotoxic chemotherapy, antiangiogenic therapy, and combined therapy for the reference mild tumor. Cytotoxic chemotherapy was very effective, as it completely eliminated the tumor right after the second drug dose (see Figure 8A). Antiangiogenic therapy took comparatively longer to eliminate the tumor, as this required eight doses instead (see Figure 8B). We observe that the cytotoxic drug induces an immediate and steep decrease in tumor volume, which gradually smoothens between consecutive cytotoxic cycles as the concentration of the drug decays. Conversely, the antiangiogenic drug first slowly decelerates tumor growth and then drives its dynamics toward its extinction. Consequently, combined therapy was mostly driven by cytotoxic effects and, again, the tumor disappeared after the second cycle (see Figure 8C). Hence, our simulations show no advantage of combined therapy over cytotoxic chemotherapy alone for the mild tumor. In all simulations including drug effects, the reference mild tumor maintained a spheroidal, round morphology until vanishing and serum PSA followed similar dynamics to tumor volume evolution.

Tumor volume and serum PSA dynamics during all simulated treatment plans for the reference aggressive tumor are plotted in Figure 9. Figures 10 and 11 respectively show the morphological evolution of the aggressive tumor and the corresponding nutrient distribution during each treatment. Cytotoxic chemotherapy was able to temporarily control tumor growth and reduce the volume of the aggressive cancer, but did not completely eliminate it (see Figure 9A). Initially, each round of cytotoxic treatment produced an instantaneous and sheer decrease in tumor volume, as in the simulations with the mild tumor. However, the aggressive tumor managed to counteract this trend as the intensity of cytotoxic action decayed with the drug concentration and eventually resumed growth between consecutive cytotoxic cycles. This rendered a saw-like tumor volume evolution during the length of chemotherapy [45, 62, 32, 36, 44, 11]. Serum PSA followed parallel dynamics to tumor volume. After the conclusion of the prescribed chemotherapeutic plan, the aggressive tumor resumed the growth dynamics exhibited before chemotherapy.

Figure 10B shows that the tumor underwent branching during cytotoxic chemotherapy, leading to the development of two vertically connected cancerous masses by the end of the simulation. This result supports our interpretation of branching as a landmark feature of aggressive tumors. Figure 11B suggests a possible mechanism that partially explains the failure of the simulated cytotoxic chemotherapy. The small tumor volumes achieved during this treatment do not globally consume as much nutrient as larger untreated tumors. Hence, the nutrient concentration around and within the tumor is comparatively higher, which remarkably favors cancer growth once the cytotoxic drug concentration has decayed enough. By the end of the chemotherapy plan, the cytotoxic drug is only able to reduce the tumor to a minimum limit volume, at which the tumor-promoting effect of the higher nutrient concentration overcomes the decaying cytotoxic action of the drug. This translates into the steady periodic pattern in tumor volume evolution and serum PSA dynamics during the last four cytotoxic cycles (see Figure 9A).

Antiangiogenic therapy was also unable to kill the aggressive cancer. While it provided a transient control on tumor volume too, it was less effective than cytotoxic chemotherapy (see Figure 9B). The first antiangiogenic drug dose slowed down tumor growth. Then, the ensuing ones induced a gradual and global decreasing trend in tumor volume, similar to that observed in the mild PCa case. The tumor-growth inhibitory effect of the antiangiogenic drug showed a little peak at the time of delivery of each dose, but it rapidly decreased afterwards and the tumor even grew momentarily between doses. Serum PSA exhibited similar dynamics to tumor volume. Figure 10C shows that the tumor developed two branches that grew vertically during the antiangiogenic treatment, without exhibiting further branching except until the end of the simulation. This suggests that branching helped the tumor counteract the poorer intratumoral nutrient supply induced by antiangiogenic therapy, hence arising again as an adaptive feature of aggressive tumors. Indeed, the morphology shift enabled the tumor to escape the inhibitory effect of antiangiogenic therapy before the end of treatment because tumor volume reaches a minimum between the eighth and ninth cycles and it exhibits a global increasing trend afterwards (see Figure 9B). Hence, antiangiogenic therapy had already failed before the delivery of the last two drug doses. Similar trends in tumor volume evolution during antiangiogenic therapy have been reported in the literature [6, 33].

Finally, combined therapy succeeded in eliminating the aggressive tumor right after the tenth cycle. Figure 9C shows that tumor volume evolution and serum PSA dynamics exhibited a global consistent decrease, even though both quantities briefly grew between consecutive doses when the aggressive tumor dynamics overcame the decaying combined tumor-inhibiting effect of the cytotoxic and antiangiogenic drugs. The tumor initially evolved as an ellipsoid. During the third and fourth cycles, it quickly developed an ellongated vertical shape showing an incipient morphology shift. However, the tumor did not fully undergo branching and had reverted its morphology back to an ellipsoid by the fifth dose. During the second half of treatment and thereafter the tumor adopted a round, spheroidal morphology (see Figure 10D). We believe that the absence of branching in this simulation may be interpreted as an early sign of successful treatment outcome. Additionally, Figure 11D shows that the nutrient concentration within small tumor volumes during combined therapy was lower than with cytotoxic chemotherapy. This is a consequence of the action of the antiangiogenic drug. Our simulations suggest that this reduction of intratumoral nutrient availability was a key enabling feature to kill the aggressive tumor.

6 Conclusions

Here, we present a mathematical model to describe the growth of prostatic tumors and their treatment using cytotoxic and antiangiogenic drugs. We describe the coupled dynamics of healthy and tumor tissue by leveraging the phase-field method and we assume that PCa growth is driven by a generic nutrient that follows reaction-diffusion dynamics [29, 51, 50, 49, 77]. Our modeling approach enables the direct use of experimentally-determined tumor proliferation and apoptotis rates within the tumor phase-field equation. Cytotoxic chemotherapy is included as a tumor-inhibiting term in the phase-field equation. Antiangiogenic therapy is modeled as a term reducing the intratumoral nutrient supply in the nutrient dynamics equation. The adopted formulation for both treatment strategies aligns with previous mathematical models of chemotherapy and the pharmacodynamic behavior of common cytotoxic and antiangiogenic drugs used in PCa treatment. To make the problem prostate-specific, we further couple the dynamics of the tumor to that of PSA [51, 50, 69], which is an extensively used biomarker of PCa in clinical management of the disease [55].

We prove that this model is well posed, that is, its solution exists and is unique in an appropriate functional space. In order to prove the existence of the solution, we apply a non-trivial argument based on the determination of fixed points in the systems of equations and initial and boundary conditions. Moreover, we derive further properties of the solution which are physically relevant. In particular, the solution component $\phi$, representing the phase variable, turns out to lie between [math] (healthy phase) and $1$ (tumor phase), and the other components $\sigma$ and $p$ are non-negative and uniformly bounded, as expected from the model.

Our simulations show that the model reproduces the common morphologies of untreated mild and aggressive tumors observed in previous computational, experimental, and clinical studies [35, 19, 58, 51, 49]. We have also identified a branching instability as a feature characterizing a malignant behavior, which enables the tumor to adapt its morphology to facilitate the access to nutrients and avoid limit concentrations that would hamper its growth [51, 28, 76]. By varying the intratumoral nutrient supply and the tumor metabolic dependence on the nutrient, we were able to reproduce the starvation of a mild tumor and the adaptive branching morphologies of aggressive tumors driven by nutrient availability. Both the evolution of tumor volume and serum PSA obtained in our simulations match the corresponding trends observed in clinical practice [55, 73]. Our model was also able to reproduce the poorer correlation between serum PSA and tumor volume in aggressive cancers exhibiting fast tumor dynamics [69].

Additionally, our results align with the use of cytotoxic chemotherapy as the reference approach to treat advanced PCa [55, 67, 18], which may show further benefit from combination with an antiangiogenic drug in order to treat aggressive tumors [2, 67, 68]. Our simulations show that a mild tumor could be effectively treated with cytotoxic chemotherapy, antiangiogenic therapy, and combined therapy. However, our results suggest that cytotoxic drugs may suffice to effectively treat mild tumors. Conversely, only combined therapy succeeded in eliminating the aggressive tumor in our simulations. We observed that the combination of cytotoxic action with the reduction of intratumoral nutrient availability provided by antiangiogenic therapy was a pivotal mechanism to kill the aggressive tumor. Indeed, the tumor could not complete the morphology shift under combined therapy and reverted to a round, spheroidal tumor, which is more common in mild tumors. Cytotoxic and antiangiogenic therapy alone provided some tumor volume control, especially the cytotoxic approach. This reduction of tumor volume could be interesting in a neoadjuvant scenario, i.e., to optimize the outcome of a subsequent radical treatment (e.g., surgery, radiotherapy) [55]. Interestingly, our simulations suggest that the branching morphology shift may enable aggressive PCa to resist cytotoxic and antiangiogenic monotherapies. Chemoresistance is commonly attributed to a reduction in the drug effect on the tumor, but in our simulations we do not change the drug effect on PCa (i.e. $\beta_{c}$ and $\beta_{a}$). Hence, the study of tumor morphology may provide early information on the clinical outcome of the treatment with cytotoxic and antiangiogenic drugs.

The mathematical model presented herein could be extended to include the evolving local angiogenic microvasculature [22, 34]. This feature would improve the model by providing (1) a direct way to characterize intratumoral nutrient supply, (2) a direct target for antiangiogenic therapy, and (3) the possibility to refine the modeling action of cytotoxic and antiangiogenic therapy because these drugs would reach the tumor through this microvasculature network [45, 62, 33, 32, 5, 6, 37]. Hence, the introduction of tumor microvasculature would enable the study of drug supply to the tumor and would permit to explore the interactions between both forms of treatment considered herein, for instance, whether the reduction of microvasculature due to antiangiogenic therapy may significantly obstruct the effective delivery of cytotoxic drugs to the tumor and whether it would be more effective to normalize the tumor microvasculature to boost the supply of cytotoxic drugs [40, 56]. Microvasculature may be modeled following alternative approaches [72]. Hybrid models enable a precise description of the morphology and evolution of microvasculature, but they are more computationally expensive [24, 72, 77]. Instead, a continuous reaction-diffusion equation describing the density of microvasculature would facilitate the computational coupling with our model and ensuing mathematical analysis at the cost of losing geometrical precision on the microvasculature [45, 33, 5, 6, 72].

Additionally, we plan to explore our mathematical model of PCa growth and cytotoxic and antiangiogenic drug effects within the context of optimal control problems [10, 16, 27, 5] to study the drug distributions that render an optimal treatment outcome. We also plan to explore in silico the effects of cytotoxic and antiangiogenic drugs in 3D organ-scale, patient-specific scenarios using the same numerical approach described herein [51, 50, 49]. Computational models of tumor growth including the effect of mechanical deformation on cancer development have been found to provide superior predictions of pathological outcome [75]. We have recently observed in a computational study that the mechanical stresses created by cancer growth and benign prostatic hyperplasia may obstruct prostatic tumor growth, which would explain the more favorable features of tumors arising in larger prostates [49]. Thus, the extension of our model of PCa growth and drug effects to a poroelastic framework would not only improve the description of prostatic tumor dynamics, but it would also enable to include the effect of mechanical stress and fluid pressure on the delivery of nutrient and drugs to the tumor [23, 41, 64]. This would permit to increase our understanding of the drug delivery and action in the complex tumor local environment as well as to refine current treatment strategies accordingly. Our mathematical model could also be extended to a multiphase approach featuring various tumor species with varying sensitivities to the prescribed drugs [32, 39] and the modeling of the tumor-inhibiting effects of the drugs can be further refined [31, 37]. Radiation effects could be included in the tumor phase-field equation as an additional cytotoxic term [62, 17, 47] within our modeling framework.

Finally, we assumed that proliferation, apoptosis, and the drug effect rates $\beta_{c}$ and $\beta_{a}$ remained constant during the simulation. However, these parameters may evolve due to treatment action and the phenotypic evolution of the tumor [44, 67, 34]. Longitudinal series of PSA and medical images would enable the periodic update of these parameters during treatment, which could refine model predictions and provide early information on treatment outcome to guide further clinical decision-making (e.g., continuity of the treatment plan or shift to other treatment options). Additionally, this approach could contribute to extend our understanding of the mechanisms of chemoresistance in prostatic tumors and provide unique guidance in the design of optimal treatment protocols accordingly. In this context, our model of PCa growth with cytotoxic and antiangiogenic drug effects offers a mathematically robust framework with a vast modeling potential in order to explore personalized drug-based treatment strategies in silico, which may assist physicians to successfully treat advanced PCa in the future.

Acknowledgements

This research activity has been performed in the framework of the Italian-Romanian collaboration agreement “Control and stabilization problems for phase field and biological systems” between the Italian CNR and the Romanian Academy. The financial support of the project Fondazione Cariplo-Regione Lombardia MEGAs-TAR “Matematica d’Eccellenza in biologia ed ingegneria come acceleratore di una nuova strateGia per l’ATtRattività dell’ateneo pavese” is gratefully acknowledged. The paper also benefits from the support of the GNAMPA (Gruppo Nazionale per l’Analisi Matematica, la Probabilità e le loro Applicazioni) of INdAM (Istituto Nazionale di Alta Matematica) for PC and ER. GL and AR have been also partially supported by the MIUR-PRIN project XFAST-SIMS (no. 20173C478N). The present paper benefits from the support of the of Ministry of Research and Innovation, CNCS –UEFISCDI, project number PN-III-P4-ID-PCE-2016-0011, for GM. The authors acknowledge the Rosen Center for Advanced Computing at Purdue University (USA) for providing HPC resources that contributed to the results presented in this paper.