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

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.



There are no comments yet.


page 20

page 21

page 22

page 23

page 26

page 27


Quantitative in vivo imaging to enable tumor forecasting and treatment optimization

Current clinical decision-making in oncology relies on averages of large...

Optimal control of cytotoxic and antiangiogenic therapies on prostate cancer growth

Prostate cancer can be lethal in advanced stages, for which chemotherapy...

Supermodeling of tumor dynamics with parallel isogeometric analysis solver

We show that it is possible to obtain reliable prognoses about cancer dy...

3D CA model of tumor-induced angiogenesis

Tumor-induced angiogenesis is the formation of new sprouts from preexist...

An image-informed Cahn-Hilliard Keller-Segel multiphase field model for tumor growth with angiogenesis

We develop a new four-phase tumor growth model with angiogenesis, derive...

Complex far-field geometries determine the stability of solid tumor growth with chemotaxis

In this paper, we develop a sharp interface tumor growth model to study ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 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 [29]. 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 [77]


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 [77]. 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 properties

and ; 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 [77]. 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




Figure 1: A high nutrient environment energetically favors tumor growth in our model, whereas low nutrient availability and the action of the cytotoxic drug energetically obstruct it. (A) Plot of with respect to the values of for (top) and (bottom). (B) Plot of with respect to the values of for three values of function (positive, zero, and negative) combined with (top) and (bottom). (C) Plot of with respect to the values of for three values of function (positive, zero, and negative) combined with (top) and (bottom).

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 [51]


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 [69].

3 Well-posedness of the model

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


where is given by (2.7) and (2.3); are given real constants; and , , are positive fixed real constants.

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 .

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


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.

Definition 3.1. Let A solution to the system (3.1)-(3.5) is a triplet , with


which satisfies




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


Then, the system (3.1)-(3.5) has a unique solution in the sense of Definition 3.1, such that

If  then Moreover, if


then we have


Finally, if  the solution has the supplementary regularity

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

for all

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 .

Let us fix in equations (3.1), (3.2):


where in (3.20) is the solution to (3.21), corresponding to , with on and in Thus, we can define the mapping

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.,

Due to a general solvability result for linear parabolic problems (the reader may consult, e.g., Ref.  [48]), for and the initial-boundary value problem for (3.21) has a unique solution

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


We are going to prove that (3.25)-(3.26) has a solution. We apply the Banach fixed point theorem, by fixing and studying the problem


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.  [48]), that (3.29), (3.26) has a unique solution

Now, we aim at proving that the operator associating to the solution to (3.29) and (3.26) is a contraction mapping on For the sake of convenience, we introduce the norm

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

and note that for the mapping is a contraction on . Thus, we deduce that (3.25)-(3.26) has a unique solution.

Next, we recall (3.15), that is a.e. in Testing (3.25) by being the negative part of , by a few calculations 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