Mathematical analysis and simulation study of a phase-field model of prostate cancer growth with chemotherapy and antiangiogenic therapy effects
Pierluigi Colli Hector Gomez Guillermo Lorenzo Gabriela Marinoschi Alessandro Reali Elisabetta Rocca
math.AP cs.CE q-bio.TO 35Q92, 92C50, 65M60, 35K51, 35K58




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




§ INTRODUCTION


Prostate cancer (PCa) is a major health problem striking over one million men annually and is responsible for over 300,000 deaths <cit.>. 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 <cit.>, 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 <cit.>. 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 <cit.>. 
Advanced PCa patients are usually prescribed androgen-deprivation therapy or chemotherapy <cit.>, 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 <cit.> 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 <cit.>.
While antiangiogenic therapies for PCa are being actively investigated <cit.>, 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 <cit.> and an increase in microvascular density is associated to poorer prognosis <cit.>. There are multiple studies with contradicting evidence about the effectiveness of antiangiogenic therapy <cit.>. However, we would like to highlight that bevacizumab (the most common antiangiogenic drug <cit.>) did not show significant benefit in castration-resistant PCa when administered alone <cit.>, but it produced a 50% PSA decline in 75% of the patients when combined with chemotherapy <cit.>. 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 <cit.>. 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 <cit.>.
In this context, several studies have been focusing on studying the effects of chemotherapy on tumors through mathematical models and computational simulations <cit.>. 
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 <cit.>, 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 <cit.> and clinical observations <cit.>. 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 <cit.>. 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.



§ MODEL OF PROSTATE CANCER GROWTH WITH CHEMOTHERAPY



 §.§ 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 <cit.>. In this work, the phase field is denoted by ϕ and transitions from the value ϕ≈0 in the host tissue to ϕ≈1 in the tumor. The transition is smooth but steep and takes on a hyperbolic tangent profile in the direction perpendicular to the interface <cit.>. 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 p is governed by a linear reaction-diffusion equation.



 §.§ Governing equations



  §.§.§ Tumor dynamics

The tumor dynamics is described by the equation <cit.>

    ϕ _t=λΔϕ -∂ G/∂ϕ(ϕ
    ,σ ,u)

where ϕ_t denotes the partial derivative of ϕ with respect to time, Δ is the Laplace operator, λ is the diffusion coefficient of tumor cells and

    G(ϕ,σ,u)=F(ϕ)-h(ϕ)(m(σ)-m_refu).

The diffusion coefficient of tumor cells can be computed as λ=Mℓ^2, where both M and ℓ are positive real constants denoting the tumor mobility and interface width, respectively <cit.>.
Here, F(ϕ)=Mϕ^2(1-ϕ)^2 is a double-well potential, i.e., a nonconvex function, typical in phase-field modeling, which allows the coexistence of the tumoral (ϕ≈1) and healthy (ϕ≈0) tissue. 

The function h(ϕ)=Mϕ^2(3-2ϕ) 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^'(0)=h^'(1)=0, where h^' denotes the derivative of h. The function m(σ) is normally called tilting function and, in our model, is defined as

    m(σ )=m_ref( ρ +A/2+ρ -A/πarctan( σ -σ _l/σ _r) ),

where ρ 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:

    ρ=K_ρ/K̅_̅ρ̅

and

    A=-K_A/K̅_̅A̅,

 where K_ρ is the proliferation rate of tumor cells, K_A is the apoptosis rate of tumor cells, K̅_̅ρ̅ is a scaling reference value for the proliferation rate, and K̅_̅A̅  is a scaling reference value for the apoptosis rate. Therefore, m(σ) 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 σ_r and σ_l are, respectively, a reference and a threshold value for the nutrient concentration <cit.>. For nutrient concentrations lower than σ_l healthy tissue is energetically more favorable than tumor tissue and viceversa. The function u in (<ref>) represents the tumor-inhibiting effect of a cytotoxic drug and is described in Section <ref>. When |m(σ)-m_refu|<1/3, the function G is a double-well potential with local minima at ϕ=0 and ϕ=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 (ϕ=0) than in the tumoral tissue (ϕ=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 <ref>. From (<ref>), we can obtain

    ∂ G/∂ϕ=2ϕ (1-ϕ )f(ϕ,σ ,u)

where

    f(ϕ ,σ ,u)=M[ 1-2ϕ -3( m(σ)-m_refu) ].






  §.§.§ Cytotoxic chemotherapy
 

Docetaxel is a usual cytotoxic drug in the clinical management of advanced PCa <cit.>. 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
 <cit.>. The pharmacodynamics of this cytotoxic drug shows an approximately exponential decrease in drug concentration after delivery of the dose at the systemic level <cit.>. Hence, we propose the formulation

    u(t)=∑_i=1^n_cβ_c d_c e^-t-T_c,i/τ_cℋ(t-T_c,i),

where n_c is the number of chemotherapy cycles, β_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,…,n_c are the times of drug delivery, τ_c 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 β_c ensure that the condition |m(σ)-m_refu|<1/3 is verified; see Section <ref>. According to (<ref>), 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 <ref>). Notice that the formulation paradigm in (<ref>) is very similar to the modeling of radiation effects on tumor dynamics <cit.>. In principle, the formulation proposed in (<ref>) is extensible to other cytotoxic drugs by adapting τ_c, β_c and d_c accordingly.



  §.§.§ Nutrient dynamics


Nutrient dynamics is controlled by the equation

    σ_t=ηΔσ+S_h(1-ϕ)+(S_c-s)ϕ-( γ_h(1-ϕ) + γ_cϕ) σ,

where η 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 <ref>), and γ_h, γ_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≤ S_c.










  §.§.§ Antiangiogenic therapy


Bevacizumab is a common antiangiogenic drug that has been actively investigated for the treatment of PCa <cit.>. 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 <cit.>. Pharmacodynamic studies of bevacizumab also revealed an exponential decay in drug concentration following systemic delivery of the prescribed dose <cit.>. Therefore, we choose an analogous formulation to (<ref>):

    s(t)=∑_i=1^n_aβ_a d_a e^-t-T_a,i/τ_aℋ(t-T_a,i)

where n_a is the number of antiangiogenic treatment cycles, β_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,…,n_a), τ_a 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 (<ref>) may be extended to other antiangiogenic drugs by appropriately calibrating τ_a, β_a and d_a.



  §.§.§ 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 γ_p. Hence, we propose the equation <cit.> 

    p_t = DΔ p + α_h(1-ϕ) + α_c ϕ - γ_p p

where D is the diffusion constant. The constants α_h and α_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 Ω, that is,

    P_s=∫_Ω p   dx.

By using (<ref>), one can show that for free-flux boundary conditions

    d P_s/d t=α_h V_h + α_c V_c - γ_p P_s,

where V_h is the volume of prostatic healthy tissue and V_c is the tumor volume. Equation (<ref>), which we derived from the tissue level equation (<ref>), was proposed as a phenomenological model of serum PSA dynamics that fits clinical data <cit.>.



§ WELL-POSEDNESS OF THE MODEL



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

    ϕ _t-λΔϕ +2ϕ (1-ϕ )f(ϕ ,σ ,u)=0,Q_T := (0,T)×Ω ,


    σ _t-ηΔσ +γ _hσ +(γ _c-γ
    _h)σϕ= S_h + (S_c -S_h)ϕ -s ϕ ,Q_T ,


    p_t-DΔ p+γ _pp=α _h+(α _c-α _h)ϕ ,
    Q_T ,


    ϕ =0,∂σ/∂ν=∂ p/
    ∂ν=0,Σ_T:= (0,T)×∂Ω ,


    ϕ (0)=ϕ _0,σ (0)=σ _0,p(0)=p_0,Ω ,

where f is given by (<ref>) and (<ref>); M, m_ref, ρ, A, σ_l, σ_r are given real constants; and  λ , η , D, γ
_c, γ _h, γ _p, α _c, α _h, S_c, S_h are positive fixed real constants.

Here, Ω is an open bounded subset of ℝ^N, N≤3, with a
sufficiently smooth boundary ∂Ω , ∂/∂ν represents the outward normal derivative to ∂Ω and T
denotes a final time. 




 §.§ Functional framework


The problem will be treated in a functional framework involving the space H=L^2(Ω ), which is identified with its dual space H^'≅ H, and the Sobolev spaces

    V_0=H_0^1(Ω ),V_0^'=(H_0^1(Ω
    ))^'=:H^-1(Ω ),V=H^1(Ω ),V^'=(H^1(Ω ))^',


    W_0=H^2(Ω )∩ H_0^1(Ω ),W={ y∈
    H^2(Ω );∂ y/∂ν=0∂Ω} ,

and with the dense and compact injections W_0⊂ V_0⊂ H⊂
V_0^' and W⊂ V⊂ H⊂ V^'. 
H^1_0(Ω) contains the elements of H^1(Ω) with null trace on the boundary ∂Ω.
For the norms in these spaces we will use the notation · _B, where 
B is the space we are considering.

For q∈ 1,∞ ] and z∈ L^q(Ω ) or z∈
L^q(Q_T), we simply denote the norm of z in these L^q spaces by ‖ z‖
_q.

We will also 
make use of spaces of functions that depend on time with values in a Banach space B. Namely, for 
 q∈ 1,∞ ] 
we consider the space L^q (0,T;B) of measurable functions t↦ z(t) such that 
t ↦ z(t)_B^q is integrable on (0,T) (or essentially bounded  if q=∞) 
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∈ 1,∞ ],  W^1,q (0,T;B)
will denote the space of functions t↦ z(t) such that both z and its (weak) derivative z_t
belong to L^q (0,T;B). We point out that  W^1,q (0,T;B) ⊂  C([0,T]; B)
for all  q∈ 1,∞ ]. 


Hypotheses. Following the considerations presented in the
Introduction and extending the particular choices of u and 
s in (<ref>) and (<ref>), we let


    u∈ L^∞(Q_T), s∈ L^∞(Q_T),

and it can be easily checked that 

    (r_1,r_2,r_3)→ f(r_1,r_2,r_3)ℝ^3.

We also set 

    γ _ch=γ _c-γ _h,  α _ch=α _c-α_h, 
      S_ch= S_c - S_h,

and point out that the problem parameters λ , η , D, γ
_c, γ _h, γ_p, α _c, α _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 Ω , the final time T, the norms of the initial data and possibly on
the norms of u and s. Moreover, we assume that (<ref>)-(<ref>) hold.





Definition 3.1. Let (φ _0,σ _0,p_0)∈
H× H× H. A solution to the system (<ref>)-(<ref>) is a triplet (ϕ,σ,p), with

    ϕ   ∈   W^1,2(0,T;V_0^')∩ C([0,T];H)∩
    L^2(0,T;V_0)∩ L^∞(Q_T), 
    σ   ∈   W^1,2(0,T;V^')∩ C([0,T];H)∩ L^2(0,T;V),
     
    
    p    ∈   W^1,2(0,T;V^')∩ C([0,T];H)∩ L^2(0,T;V),

which satisfies 

    ∫_0^T⟨ϕ _t(t),ψ _1(t)⟩
    _V_0^',V_0dt+∫_Q_T{λ∇ϕ·∇ψ
    _1+2ϕ (1-ϕ )f(ϕ ,σ ,u)ψ _1} dxdt
    
       +∫_0^T⟨σ _t(t),ψ _2(t)⟩
    _V^',Vdt+∫_Q_T(η∇σ·∇ψ
    _2+γ _hσψ _2+γ _chσϕψ _2)dxdt 
    
       =∫_Q_T(S_hψ _2 + (S_ch - s) ϕψ_2)dxdt,(ψ _1,ψ _2)∈
    L^2(0,T;V_0× V),


    ∫_0^T⟨ p_t(t),ψ (t)⟩ _V^',Vdt+∫_Q_T(D∇ p ·∇ψ +γ _ppψ )dxdt
    
          =∫_Q_T(α _h+α _chϕ )ψ dxdt,ψ∈ L^2(0,T;V),

and 

    (φ ,σ ,p)(0)=(φ _0,σ _0,p_0).

Please note that taking  in (<ref>) ψ_2 first, and ψ_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 
ϕ and σ, as it is (<ref>) for p.


Let us denote

    X_0   =   W^1,2(0,T;H)∩C([0,T];V_0)∩ L^2(0,T;W_0),
    
    
    X    =   W^1,2(0,T;H)∩C([0,T];V)∩ L^2(0,T;W).




Theorem 3.2. Let   

    (ϕ _0,σ _0,p_0)∈ H× H× H,


    0≤ϕ _0(x)≤ 1x∈Ω .

Then, the system (<ref>)-(<ref>) has a unique solution (ϕ,σ,p) in
the sense of Definition 3.1, such that 

    0≤ϕ (t,x)≤ 1(t,x)∈ Q_T.

If (σ _0,p_0)∈ L^∞(Ω )× L^∞(Ω ), then (σ ,p)∈ L^∞(Q_T)×
L^∞(Q_T). Moreover, if 

    σ _0(x)≥ 0,p_0(x)≥ 0x∈Ω
    ,s(t,x)≤ S_c(t,x)∈ Q_T.

then we have 

    σ (t,x)≥ 0,p(t,x)≥ 0(t,x)∈ Q_T.

Finally, if (ϕ _0,σ _0,p_0)∈ V_0× V×
V, the solution has the supplementary regularity (ϕ ,σ
,p)∈ X_0× X× X, and satisfies the estimate 

    ‖ϕ‖ _X_0+‖σ‖
    _X+‖ p‖ _X
       ≤   C(‖ϕ _0‖ _V_0^2+‖σ
    _0‖ _V^2+‖ p_0‖ _V^2+‖
    u‖ _L^2(0,T;H)^2+‖s‖ _L^2
    (0,T;H)^2+1).

In addition, the solution is continuous with respect to the data,
that is, for two solutions (ϕ _i,σ _i,p_i)corresponding to the data (ϕ _0^i,σ _0^i,p_0^i,u_i, s_i),
 i=1,2, we have 

    ‖ (ϕ _1-ϕ _2)(t)‖ _H^2+‖ (σ
    _1-σ _2)(t)‖ _H^2+‖
    (p_1-p_2)(t)‖ _H^2
          +‖ϕ _1-ϕ _2‖
    _L^2(0,T;V_0)^2+‖σ _1-σ _2‖
    _L^2(0,T;V)^2+‖ p_1-p_2‖ _L^2(0,T;V)^2
       ≤   C( ‖ϕ _0^1-ϕ _0^2‖
    _H^2+‖σ _0^1-σ _0^2‖
    _H^2+‖ p_0^1-p_0^2‖ _H^2.  
          . +‖ u_1-u_2‖
    _L^2(0,T;H)^2+‖s_1- s_2‖
    _L^2(0,T;H)^2)

for all t∈ 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  (<ref>)-(<ref>) 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 

    ℳ={z∈ C([0,T];H):0≤ z(t,x)≤ 1
    (t,x)∈ Q_T},

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 equations (<ref>), (<ref>):

    ϕ _t-λΔϕ =-2ϕ (1-ϕ )f(z,σ ^z,u),
    Q_T,


    σ _t-ηΔσ +γ _hσ +γ _chσ z= S_h+ (S_ch - s ) z,
    Q_T,

where σ ^z in (<ref>) is the solution to (<ref>),
corresponding to z, with ∂σ/∂ν=0 on Σ _T and σ (0)=σ _0 in Ω . Thus, we can define
the mapping 

    z↦ϕ ^z=:Ψ (z)

where ϕ ^z is the solution to (<ref>), corresponding to z and σ^z,
with the homogeneous Dirichlet boundary condition and the initial datum ϕ _0. We shall prove that Ψ (ℳ)⊂ℳ and
that Ψ is a contraction mapping on ℳ.

First, we treat the initial-boundary value problem for equation (<ref>). For all t∈ 0,T] we
introduce the operator A(t):V→ V^' defined by 

    ⟨ A(t)σ ,ψ⟩ _V^',V=∫_Ω(η∇σ·∇ψ +γ _hσψ +γ
    _chzσψ )dx,σ, ψ∈ V,

and observe that it is continuous 

    ‖ A(t)σ‖ _V^'≤max{η ,γ
    _h,|γ _ch|}‖σ‖ _V

and quasi-monotone from V to V^', i.e.,

    ⟨ A(t)σ ,σ⟩ _V^',V≥min{η ,γ _h}‖σ‖ _V^2-|γ _ch|‖σ‖ _H^2.

Due to a general solvability result  for linear parabolic problems (the reader 
may consult, e.g., Ref.  <cit.>), 
for σ _0∈ H and s∈ L^∞(Q_T) the
initial-boundary value problem for (<ref>) has a unique solution 

    σ ^z∈ W^1,2(0,T;V^')∩ C([0,T];H)∩
    L^2(0,T;V).

A first estimate performed by testing (<ref>) by σ ^z and
integrating over (0,t) yields, after a standard calculation,

    ‖σ ^z(t)‖ _H^2+∫_0^t‖σ
    ^z(τ )‖ _V^2dτ≤ C(1 +‖σ
    _0‖ _H^2+‖s‖ _L^2(0,T;H)^2),

for all t∈ 0,T], where we used the fact that z∈ℳ. 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 
(<ref>).

Now, if we take two solutions σ _1 and σ _2
corresponding to (z_1,s_1 ,σ_0^1) and (z_2,s_2,σ_0^2), 
thanks to the Hölder and Young inequalities
we obtain for the differences σ :=σ _1-σ _2
and z:=z_1-z_2: 

    1/2‖σ (t)‖ _H^2+min{η ,γ
    _h}∫_0^t‖σ (τ )‖ _V^2dτ
       ≤1/2‖σ _0^1-σ_0^2‖
    _H^2+|γ _ch|∫_0^t∫_Ω{| zσ _1| +|σ
    z_2|}|σ| dxdτ
         +∫_0^t∫_Ω| S_chzσ| dxdτ 
    +∫_0^t∫_Ω(| s_1 z| + | (s_1 - s_2) z_2| )|σ|  dxdτ
       ≤1/2‖σ _0^1-σ_0^2‖
    _H^2+|γ _ch|∫_0^t‖ z(τ
    )‖ _H‖σ _1(τ )‖ _4‖σ (τ )‖ _4dτ
         +|γ _ch|∫_0^t‖σ (τ )‖ _H^2dτ+ |S_ch|^2 ∫_0^t‖ z (τ )‖ _H^2dτ
    + 2 ∫_0^t s_1_∞^2 ‖ z (τ )‖ _H^2dτ
         +2∫_0^t‖ (s_1-s_2)(τ )‖ _H^2dτ
    +1/2∫_0^t‖σ (τ )‖ _H^2dτ .

By the Sobolev embedding theorems, it turns out that V⊂ L^4(Ω )
with continuous embedding. Then, by exploiting the Young inequality and (<ref>) for σ_1 we infer that 

    |γ _ch|∫_0^t‖ z(τ
    )‖ _H‖σ _1(τ )‖ _4‖σ (τ )‖ _4dτ≤ C∫_0^t‖ z(τ
    )‖ _H‖σ _1(τ )‖ _V‖σ (τ )‖ _Vdτ
       ≤   1/2min{η ,γ _h}∫_0^t‖σ
    (τ )‖ _V^2dτ +C∫_0^t‖ z(τ
    )‖ _H^2‖σ _1(τ )‖ _V^2dτ
       ≤   1/2min{η ,γ _h}∫_0^t‖σ
    (τ )‖ _V^2dτ +C‖ z‖
    _C([0,t];H)^2∫_0^t‖σ _1(τ )‖
    _V^2dτ
       ≤   1/2min{η ,γ _h}∫_0^t‖σ
    (τ )‖ _V^2dτ +C‖ z‖
    _C([0,t];H)^2( 1+‖σ _0^1‖ _H^2+‖
    s_1‖ _L^2(0,T;H)^2).

Going back to (<ref>) we deduce that

    ‖σ (t)‖ _H^2+min{η ,γ
    _h}∫_0^t‖σ (τ )‖ _V^2dτ
       ≤   ‖σ _0^1-σ_0^2‖
    _H^2+4∫_0^t‖ (s_1-s_2)(τ )‖
    _H^2dτ +C_01‖ z‖
    _C([0,t];H)^2+C∫_0^t‖σ (τ )‖
    _H^2dτ,

where C_01:=C( 1+ ‖σ _0^1‖
_H^2+‖s_1‖^2_∞). Finally,
by using Gronwall lemma it is not difficult to conclude that

    ‖σ _1-σ _2‖ _C([0,t];H)∩
    L^2(0,t;V)^2
       ≤   C( ‖σ _0^1-σ _0^2‖
    _H^2+‖s_1-s_2‖ _L^2(0,t;H)^2+‖
    z_1-z_2‖ _C([0,t];H)^2) ,t∈ 0,T],

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 (<ref>). To this end, for z∈ℳ and σ^z defined as above 
we consider the intermediate equation

    ϕ _t-λΔϕ =g(ϕ )f(z,σ ^z,u),Q_T,

with 

    ϕ =0, Σ _T;  ϕ (0)=ϕ _0, Ω ,

where 

    g(r)={[ -2r(1-r),r∈ [ 0,1],;                  0, ]
    .

We observe that g is bounded and Lipschitz continuous on ℝ.
Recalling (<ref>) we note that f(z,σ ^z,u)∈ L^∞(Q_T)
and

    ‖ f(z,σ ^z,u)‖ _∞≤ C(1+‖
    u‖ _∞).

We are going to prove that (<ref>)-(<ref>) has a solution. We
apply the Banach fixed point theorem, by fixing ζ∈ C([0,T];H) and
studying the problem 

    ϕ _t-λΔϕ =g(ζ )f(z,σ ^z,u),Q_T,

together with (<ref>). We note that the right-hand side of (<ref>) is in L^∞(Q_T) and that ϕ _0∈ H. It is clear, again by the 
solvability result for linear parabolic problems (see Ref.  <cit.>), that 
(<ref>), (<ref>) has a unique solution 

    ϕ∈ W^1,2(0,T;V_0^')∩ C([0,T];H)∩ L^2(0,T;V_0).


Now, we aim at proving that the operator associating ζ to the solution ϕ
to (<ref>) and (<ref>) is a contraction mapping on C([0,T];H). For 
the sake of convenience, we introduce the norm 

    ‖ζ‖ _B:=sup_t∈ 0,T]e^-γ
    t‖ζ (t)‖ _H,γ >0,ζ∈
    C([0,T];H),

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

Let us consider two solutions ϕ _1 and ϕ _2 to (<ref>), (<ref>) corresponding to ζ _1 and ζ _2, respectively. We
multiply the difference of the equations  (<ref>) by ϕ _1-ϕ _2 and
integrate over Ω× (0,t). With the help of (<ref>) and (<ref>) we have that 

    1/2‖ (ϕ _1-ϕ _2)(t)‖ _H^2+λ∫_0^t‖∇ (ϕ _1-ϕ _2)(τ )‖
    _H^2dτ
       ≤   ∫_0^t∫_Ω| g(ζ _1)-g(ζ
    _2)|| f(z,σ ^z,u)||ϕ
    _1-ϕ _2| dxdτ
       ≤   C∫_0^t‖ (ζ _1-ζ _2)(τ )‖
    _H‖ (ϕ _1-ϕ _2)(τ )‖ _Hdτ .

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

    ‖ (ϕ _1-ϕ _2)(t)‖ _H^2≤
    C∫_0^t‖ (ζ _1-ζ _2)(τ )‖
    _H^2dτ .

We multiply this inequality by e^-2γ t and have successively

    e^-2γ t‖ (ϕ _1-ϕ _2)(t)‖ _H^2≤
    Ce^-2γ t∫_0^te^2γτe^-2γτ‖
    (ζ _1-ζ _2)(τ )‖ _H^2dτ
       ≤   Ce^-2γ t∫_0^te^2γτ‖ζ
    _1-ζ _2‖ _B^2dτ≤C/2γ
    (1-e^-2γ t)‖ζ _1-ζ _2‖ _B^2,t∈ 0,T].

Taking the supremum with respect to t∈ 0,T], we obtain

    ‖ϕ _1-ϕ _2‖ _B≤√(C/2γ)
    ‖ζ _1-ζ _2‖ _B

and note that for γ >C/2 the mapping ζ↦ϕ
is a contraction on C([0,T];H). Thus, we deduce that (<ref>)-(<ref>) has a unique solution.

Next, we recall (<ref>), that is 0≤ϕ _0≤ 1 a.e. in Ω . Testing (<ref>) by -ϕ ^- (ϕ ^- being the
negative part of ϕ ), by a few calculations we obtain

    1/2‖ϕ ^-(t)‖ _H^2+λ∫_0^t‖∇ϕ ^-(τ )‖ _H^2dτ
       ≤1/2‖ϕ _0^-‖
    _H^2-∫_0^t∫_Ωϕ ^-g(ϕ )f(z,σ
    ^z,u)dxdτ =0,

since ϕ_0 ^-=0 a.e. in Ω,
ϕ ^-=0 a.e. in the set where ϕ∈ 0,1] and g(ϕ )=0 a.e. in the set where ϕ∉ 0,1]. This implies
that ϕ ^-(t)=0, that is ϕ (t)≥ 0 a.e. in Ω, 
for all t∈ 0,T].

Next, we multiply (<ref>) by (ϕ -1)^+ and by similar calculations
we deduce that

    1/2‖ (ϕ -1)^+(t)‖ _H^2+λ∫_0^t‖∇ (ϕ -1)^+(τ )‖ _H^2dτ
    =0.

Here, we used the fact that ϕ =0 a.e. on Σ _T, so that (ϕ
-1)^+=0 a.e. on Σ _T. Thus, we have shown that 
ϕ (t)≤ 1 a.e. in Ω, for all t∈ 0,T].

By these two last results we actually proved that 0≤ ϕ≤ 1
a.e. on Q_T, implying that g(ϕ )=-2ϕ (1-ϕ ). It turns out
that this solution ϕ actually solves equation (<ref>). Moreover,
we have that ϕ satisfies (cf. (<ref>))

    ϕ∈ W^1,2(0,T;V_0^')∩ C([0,T];H)∩
    L^2(0,T;V_0)∩ L^∞(Q_T).


Since the problem (<ref>), (<ref>) may also have other solutions fulfilling (<ref>),
we have to prove a uniqueness result. For that, let us consider two
solutions ϕ _1 and ϕ _2 corresponding to the same data, and compute the
difference of the respective equations  (<ref>) tested by ϕ _1-ϕ
_2. Owing to a straightforward calculation and (<ref>), we deduce that 

    1/2‖ (ϕ _1-ϕ _2)(t)‖ _H^2+λ∫_0^t‖∇ (ϕ _1-ϕ _2)(τ )‖ _H^2dτ
       =   -2∫_0^t∫_Ω(ϕ _1(1-ϕ _1)-ϕ _2(1-ϕ
    _2))f(z,σ ^z,u)(ϕ _1-ϕ _2)dxdτ
       ≤   2‖ f(z,σ ^z,u)‖ _∞∫_0^t∫_Ω(ϕ _1-ϕ _2)^2| 1-(ϕ
    _1+ϕ _2)| dxdτ
       ≤   C∫_0^t‖ (ϕ _1-ϕ _2)(τ )‖
    _H^2dτ ,t∈ 0,T],

which by Gronwall lemma leads to the desired uniqueness property.

The solution ϕ to (<ref>), (<ref>) we have found is
actually what we have denoted by ϕ ^z=Ψ (z). We have already shown
that ϕ ^z∈ C([0,T];H) and that 0 ≤ϕ ^z≤1 a.e. in Q_T, whence it follows that Ψ (ℳ)⊂ℳ. 
Thus, for the last step it suffices to check that Ψ is a contraction mapping on C([0,T];H).

First, we shall deduce a general estimate.

Let us consider two pairs of data 
(ϕ _0^i,σ _0^i,u_i,s_i;z_i), i=1,2, and the corresponding
solutions (ϕ _i=ϕ ^z_i,σ _i=σ ^z_i), i=1,2.
Then, we take the difference of the equations (<ref>), 
add λ(ϕ _1-ϕ _2) to both sides, then
multiply by ϕ _1-ϕ _2 and integrate over Ω× (0,t) . Here, let us denote 
ϕ :=ϕ _1-ϕ _2, ϕ _0:=ϕ _0^1-ϕ _0^2, and take into account that f is Lipschitz continuous with
a Lipschitz constant L_f depending on M, m_ref, ρ , A, σ _r. We have 

    1/2‖ϕ (t)‖ _H^2+λ∫_0^t‖ϕ (τ )‖ _V^2dτ
       ≤   1/2
    ‖ϕ _0‖ _H^2+λ∫_0^t‖ϕ
    (τ )‖ _H^2dτ+2∫_0^t∫_Ωϕ^2 | 1-ϕ _1-ϕ
    _2|| f(z_1,σ
    _1,u_1)| dxdτ
          +2∫_0^t∫_Ω| f(z_1,σ
    _1,u_1)-f(z_2,σ _2,u_2)||ϕ
    _2(1-ϕ _2)||ϕ|
    dxdτ
       ≤   1/2‖ϕ _0‖ _H^2+( λ
     +4‖
    f(z_1,σ _1,u_1)‖ _∞) ∫_0^t‖ϕ (τ )‖ _H^2dτ
          +2L_f∫_0^t(‖ (z_1-z_2)(τ )‖
    _H+‖ (σ _1-σ _2)(τ )‖ _H+‖
    (u_1-u_2)(τ )‖ _H)‖ϕ (τ )‖
    _Hdτ .

Therefore, recalling (<ref>) and using the Young inequality, we infer
that 

    1/2‖ϕ (t)‖ _H^2+λ∫_0^t‖ϕ (τ )‖ _V^2dτ
       ≤1/2‖ϕ _0‖
    _H^2+C∫_0^t‖ϕ (τ )‖ _H^2dτ 
    +2L_f∫_0^t‖ (z_1-z_2)(τ )‖
    _H^2dτ
         +2L_f( ∫_0^t‖ (σ _1-σ _2)(τ
    )‖ _H^2dτ +∫_0^t‖ (u_1-u_2)(τ
    )‖ _H^2dτ) ,

where C depends on ‖ u_1‖ _∞ as well. Next,
we use (<ref>) and obtain

    1/2‖ϕ (t)‖ _H^2+λ∫_0^t‖ϕ (τ )‖ _V^2dτ
          ≤1/2‖ϕ _0‖
    _H^2+C∫_0^t‖ϕ (τ )‖ _H^2dτ
    +C∫_0^t‖ z_1-z_2‖ _C([0,τ ];H)^2dτ
            +C( ‖σ _0^1-σ _0^2‖
    _H^2+‖s_1-s_2‖ _L^2(0,T;H)^2+‖
    u_1-u_2‖ _L^2(0,T;H)^2) .

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


    ‖ϕ (s)‖ _H^2+∫_0^s‖ϕ (τ
    )‖ _V^2dτ
    
    
    
    
    
          ≤ C∫_0^t‖ z_1-z_2‖ _C([0,τ ];H)^2dτ
    +C ( ‖ϕ _0‖ _H^2+‖σ
    _0^1-σ _0^2‖ _H^2) 
            +C(‖s_1-s_2‖ _L^2(0,T;H)^2+‖
    u_1-u_2‖ _L^2(0,T;H)^2),


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



    ‖ϕ _1-ϕ _2‖
    _C([0,t];H)^2+∫_0^t‖ (ϕ _1-ϕ _2)(τ
    )‖ _V^2dτ
          ≤ C∫_0^t‖
    z_1-z_2‖ _C([0,τ ];H)^2dτ +C(‖ϕ _0^1-ϕ _0^2‖ _H^2+‖σ _0^1-σ _0^2‖ _H^2) 
            + C(‖s_1-s_2‖ _L^2(0,T;H)^2+‖ u_1-u_2‖
    _L^2(0,T;H)^2),

for all t∈ 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

    ‖ϕ _1-ϕ _2‖ _C([0,t];H)^2≤
    C∫_0^t‖ z_1-z_2‖ _C([0,τ ];H)^2dτ .

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

    ||| z |||
    =sup_t∈ 0,T]e^-γ t‖ z
    ‖ _C([0,t];H),γ >0,z ∈ C([0,T];H),

and playing with (<ref>), we plainly obtain 

    |||ϕ _1-ϕ _2 ||| ≤√(C/2γ) |||z_1-z_2 ||| .

Then, if γ >C/2, the mapping z↦ϕ ^z=Ψ
(z) is a contraction on ℳ and ϕ ^z=z is the unique fixed
point. Thus, problem (<ref>)-(<ref>) 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 ϕ _i, i=1,2, in (<ref>) and apply the
Gronwall lemma to infer that 

    ‖ϕ _1-ϕ _2‖ _C([0,T];H)∩
    L^2(0,T;V)^2
       ≤    C(‖ϕ _0^1-ϕ _0^2‖
    _H^2+‖σ _0^1-σ _0^2‖
    _H^2+‖ u_1-u_2‖ _L^2(0,T;H)^2
    +‖s_1-s_2‖ _L^2(0,T;H)^2) .

Combining this estimate with (<ref>) finally leads to 

    ‖ϕ _1-ϕ _2‖ _C([0,T];H)∩
    L^2(0,T;V)^2+‖σ _1-σ _2‖
    _C([0,T];H)∩ L^2(0,T;V)^2
       ≤   C(‖ϕ _0^1-ϕ _0^2‖
    _H^2+‖σ _0^1-σ _0^2‖
    _H^2+‖ u_1-u_2‖ _L^2(0,T;H)^2
    +‖s_1-s_2‖ _L^2(0,T;H)^2).


Let us now pass to (<ref>) in which the right-hand side is in L^∞(Q_T). It is easily seen, relying on the arguments before for σ,
that the probem based on (<ref>), with the initial datum p_0∈ H and 
the homogeneous Neumann boundary condition, has a unique solution 

    p∈ W^1,2(0,T;V^')∩ C([0,T];H)∩ L^2(0,T;V).


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

    ‖ p_1-p_2‖ _C([0,T];H)∩ L^2(0,T;V)^2
       ≤‖ p_0^1-p_0^2‖ _H^2+C‖ϕ
    _1-ϕ _2‖ _L^2(0,T;H)^2
       ≤ C (‖ p _0^1-p _0^2‖
    _H^2+‖σ _0^1-σ _0^2‖
    _H^2 +‖ϕ _0^1-ϕ _0^2‖_H^2 .
    
             .+‖ u_1-u_2‖ _L^2(0,T;H)^2
    +‖s_1-s_2‖ _L^2(0,T;H)^2)

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

Next, we prove the boundedness of σ and p and the property (<ref>). We rewrite eq. (<ref>) as 

    σ _t-ηΔσ =  S_h(1-ϕ) +(S_c -s) ϕ 
    -γ _hσ -γ _chσϕ∈ L^∞(0,T;H),

and consider σ _0∈ L^∞(Ω ). By Theorem 7.1, p. 181, in 
Ref.  <cit.>, we find that 

    σ∈ L^∞(Q_T).

Let σ _0(x)≥ 0, s(t,x)≤ S_c a.e. and multiply (<ref>) by -σ ^-. We have 

    1/2‖σ ^-(t)‖ _H^2+η∫_0^t‖∇σ ^-(τ )‖ _H^2dτ
    +γ _h∫_0^t‖σ ^-(τ )‖
    _H^2dτ
       ≤   1/2‖σ _0^-‖ _H^2+|γ _ch|∫_0^t∫_Ω(σ ^-)^2ϕ
    dxdτ -∫_0^t∫_Ω( S_h (1-ϕ) + (S_c - s) ϕ)σ ^-dxdτ    
       ≤   |γ _ch|∫_0^t‖σ
    ^-(τ )‖ _H^2dτ ,

since 0≤ϕ≤ 1 and (S_c - s)≥ 0 a.e. in Q_T. By Gronwall lemma
we obtain that σ ^-(t)=0, hence σ (t)≥ 0 a.e. in Ω, for all t∈ 0,T].

By a similar argument we deduce that p∈ L^∞(Q_T) and p(t)≥
0 a.e. in Ω, for all t∈ 0,T].

Finally, we prove the further regularity of the solution and 
the estimate (<ref>). The properties ϕ∈ X_0, σ∈ X, p ∈ X
can be compared with  (<ref>) and we can observe that, e.g., for σ and p 
the triplet of spaces (H,V,W) now replaces (V', H,V). This fits into the frame 
of the structural regularity of solutions to parabolic systems (see, e.g., Ref.  <cit.>).

Let ϕ _0∈ V_0 and write (<ref>) as 

    ϕ _t-λΔϕ =-2ϕ (1-ϕ )f(ϕ ,σ ,u)∈
    L^2(0,T;H).

Since the operator B_0:W_0⊂ H→ H, B_0=-λΔ is
potential (that is, B_0ϕ =∂ j(ϕ ), j:V_0→[0,+∞], j(ϕ )=λ/2∫_Ω|∇ϕ| ^2dx),
according e.g. to Theorem 4.18, p. 179, in Ref.  <cit.>,  it
follows that 

    ϕ∈ W^1,2(0,T;H).

Then, we can test (<ref>) by ϕ _t; taking into account that (cf. (<ref>)) 

    ‖ f(ϕ ,σ ,u)‖ _∞≤ C(1+‖ϕ‖ _∞+‖ u‖ _∞),

we have

    ∫_0^t‖ϕ _t (τ )‖ _H^2dτ +
    λ/2‖ϕ (τ )‖ _V_0^2≤
    λ/2‖ϕ _0‖
    _V_0^2+C∫_0^t(1+‖ 
    u‖ _∞)‖ϕ _t(τ)‖_Hdτ ,

since ‖ϕ‖ _∞≤ 1. This yields 

    ‖ϕ _t‖ _L^2(0,T;H)+‖ϕ
    (t)‖ _V_0^2≤ C(‖ϕ _0‖
    _V^2+‖ u‖ _L^∞(Q_T)^2+1),t∈ 0,T].

By comparison in the equation for ϕ we realize that Δϕ∈
L^2(0,T;H). Then, using the elliptic regularity theory, due to the homogeneous Dirichlet
boundary condition for ϕ, we can conclude that 

    ‖ϕ‖ _L^2(0,T;W_0)^2≤ C(‖ϕ
    _0‖ _V^2+‖ u‖ _L^2(0,T;H)^2+1).

Thus, we have obtained ϕ∈ X_0.

Assume that σ _0∈ V and consider again the equation (<ref>)
written as 

    σ _t-ηΔσ +γ _hσ =S_h + ( S_ch - s) ϕ -γ _chσϕ∈ L^2(0,T;H).

The operator -ηΔ +γ _hI : W⊂ H → H (I standing for the identity in H) is potential as well, and then σ∈
W^1,2(0,T;H). Multiplying (<ref>) by σ _t and integrating over Ω× (0,t), we easily arrive at

    ∫_0^t‖σ _t(τ )‖ _H^2dτ
    +‖∇σ (t)‖ _H^2≤ C(‖σ
    _0‖ _V^2+‖s‖ _L^2(0,T;H)^2+1),

for all t∈ 0,T]. Also, from a comparison in (<ref>) it follows that Δσ∈ L^2(0,T;H). 
This property, together with the homogeneous Neumann boundary condition for σ, entails

    ‖σ‖ _L^2(0,T;W)^2≤ C(‖σ
    _0‖ _V^2+‖s‖ _L^2(0,T;H)^2+1).

A similar argumentation for p yields 

    ‖ p‖ _W^1,2(0,T;H) ∩ C([0,T];V)  ∩ L^2(0,T;W)^2≤ C(‖
    p_0‖ _V^2+‖ϕ‖ _L^2(0,T;H)^2+1).

By collecting  (<ref>), (<ref>), (<ref>)-(<ref>) it is easy to get (<ref>), as claimed. Let us add the comment that, under the additional regularity 
specified by the estimate (<ref>), the triplet (ϕ,σ, p) turns out to be a strong solution
to the system (<ref>)-(<ref>), that is, it actually satisfies equations and boundary conditions almost everywhere in Q_T and on Σ_T, respectively. By this remark, we conclude our proof.
□




§ NUMERICAL METHOD




 §.§ 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(Ω) and V_0=H^1_0(Ω), defined in Section <ref>.
The weak form can be stated as follows: find ϕ∈V_0, σ∈V and p∈V such that for all {w_1,w_2,w_3}∈ V_0 × V × V,

    ∫_Ω w_1 [ ϕ_t + 2ϕ(1-ϕ)f(ϕ,σ,u) ]dx + ∫_Ωλ∇ w_1·∇ϕ dx  = 0, 
       ∫_Ω w_2 [ σ_t + ( γ_h(1-ϕ) + γ_cϕ ) σ - S_h(1-ϕ) - (S_c-s)ϕ] dx
            +  ∫_Ωη∇ w_2·∇σ dx  = 0, 
       ∫_Ω w_3 [ p_t + γ_p p -α_h (1-ϕ) - α_cϕ] dx + ∫_Ω D∇ w_3·∇ p dx  = 0.

Once the equations have been written in weak form, they can be discretized by restricting the weak form to finite dimensional spaces V^h⊂V and V_0^h⊂V_0. In this work, we construct these spaces using quadratic B-spline basis functions <cit.>. As an example, we take V^h= span{N_j}_j=1,…,n_b where the N_j's are multivariate B-splines. The function ϕ^h, which is a finite dimensional approximation to ϕ, is defined as ϕ^h(t,x)=∑_j=1^n_bϕ_j(t)N_j(x) where the ϕ_j's are referred to as control variables. The functions σ^h, p^h, w_i^h for i=1,2,3 are defined analogously.
Because ϕ^h∈ V^h_0, the control variables ϕ_j corresponding to the control points j on the boundary must be zero to satisfy the homogeneous Dirichlet boundary conditions.



 §.§ Time discretization


We integrate in time using the generalized-α method <cit.>. To define our time integration scheme, we introduce the global vector of degrees of freedom Φ={ϕ_j}_j=1,…,n_b. Analogously, we define Σ and P as the global vectors of control variables for σ^h and p^h, respectively. Let us introduce the following residual vectors

    R={R_j };   R_j = {R_ϕ,j, R_σ,j, R_p,j},

where

    R_ϕ,j = ∫_Ω N_j [ ϕ^h_t + 2ϕ^h(1-ϕ^h)f(ϕ^h,σ^h,u) ]dx + ∫_Ωλ∇ N_j·∇ϕ^h dx,   
        R_σ,j =∫_Ω N_j [ σ^h_t + ( γ_h(1-ϕ^h) + γ_cϕ^h ) σ^h  - S_h(1-ϕ^h) - (S_c-s)ϕ^h ] dx 
             +∫_Ωη∇ N_j·∇σ^h dx,   
        R_p,j =∫_Ω N_j [ p^h_t + γ_p p^h -α_h (1-ϕ^h) - α_cϕ^h ] dx + ∫_ΩD ∇ N_j·∇ p^h dx.

Let us call U_n={Φ_n, Σ_n, P_n } the time discrete approximation to the control variables at time t_n. Our time-integration algorithm can be defined as follows: given U⃗_n, U̇⃗̇_n and the current time step Δ t_n=t_n+1-t_n, find U⃗_n+1, U̇⃗̇_n+1 such that

    R⃗(U̇⃗̇_n+α_m,U⃗_n+α_f)=0⃗,  
          U⃗_n+1=U⃗_n+Δ t_nU̇⃗̇_n+γΔ t_n(U̇⃗̇_n+1-U̇⃗̇_n),  
          U̇⃗̇_n+α_m=U̇⃗̇_n+α_m(U̇⃗̇_n+1-U̇⃗̇_n),  
          U⃗_n+α_f=U⃗_n+α_f(U⃗_n+1-U⃗_n).

The parameters α_m, α_f and γ define the stability and accuracy of the algorithm. The generalized-α method can be proven to be second-order accurate and A-stable by taking ρ_∞∈[0,1] and

    α_m=1/2(3-ρ_∞/1+ρ_∞),  α_f=1/1+ρ_∞,  γ=1/2+α_m-α_f.

In our calculations, we take ρ_∞=1/2. The resulting nonlinear system of equations is linearized using the Newton-Raphson algorithm. We converged all the individual residuals R_ϕ, R_σ, and R_p up to a certain tolerance ε_NR. This strategy ensured that the fulfilment of the convergence criterion by the global residual R was not due to potentially different scaling of the equations of the system (<ref>)-(<ref>).  Within each iteration of  the Newton-Raphson algorithm, we solved the corresponding linearized system by means of the GMRES <cit.> algorithm with a diagonal preconditioner. We chose a certain tolerance ε_GMRES and a maximum number of iterations as convergence criteria for the GMRES algorithm.



§ SIMULATION STUDY



 §.§ 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 ρ and A within the function m(σ) in (<ref>). In general, tumor growth is usually characterized by a high proliferation index and a low apoptotis index <cit.>, which would translate into choosing |ρ|>|A| within the context of our model. Then, we may represent mild tumors with lower values of these parameters for which the difference |ρ|-|A| is also low, while more aggressive tumors can be modeled with higher values and such that the difference |ρ|-|A| is also larger. For the simulations in this study, we calibrated ρ and A using the proliferation and apoptosis rates previously reported in the literature <cit.> in equations (<ref>)-(<ref>). We further computationally calibrated m_ref in (<ref>)-(<ref>) to ensure |m(σ)-m_refu|<1/3 (see Section <ref>).

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 <cit.>. 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 (<ref>), respectively. Moreover, we analyze the effect of tumor metabolism on its dynamics by choosing a larger or a smaller value of γ_c in (<ref>). 

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 <ref>, 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/m^2 of docetaxel delivered every three weeks <cit.>. We set τ_c=5 days, according to pharmacodynamic studies of this drug <cit.>. 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 <cit.>.  We set τ_a=30 days <cit.>. We computationally calibrated the value of the constants β_c and β_a to produce a meaningful chemotherapeutic response as observed in clinical literature. 

Table <ref> 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 <cit.>.





 §.§ Computational setup


All simulations were run on a 2D square domain with side length of L_d=3000  and 256 isogeometric elements per side. 
We chose a constant time step of Δ t=0.1 days in all simulations.
The convergence of the Newton-Raphson method was set at tolerance ε_NR=10^-3, while for the GMRES algorithm was set at ε_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  and b=200  parallel to the sides of domain. To implement this initial condition, we L^2-projected the hyperbolic tangent function

    ϕ_0(x,y) = 0.5 - 0.5tanh( 10(√((x-L_d/2)^2/a^2 + (y-L_d/2)^2/b^2) - 1) )

over the quadratic B-spline space supporting our spatial discretization. This operation provides the control variables ϕ_0,A of the spline representation of the phase-field initial condition, i.e., ϕ^h_0(x,t)=∑_A=1^n_bϕ_0,A(t)N_A(x) (see Section <ref>).

The initial conditions for the nutrient and the tissue PSA are approximations based on ϕ_0 given by 

    σ_0 = c^0_σ + c^1_σϕ_0

and

    p_0 = c^0_p + c^1_pϕ_0.

The constants c^0_σ, c^1_σ, c^0_p, and c^1_p are computationally calibrated <cit.>, such that σ_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_σ=1 g/L, c^1_σ=-0.8 g/L, c^0_p=0.0625 ng/mL/cc, and c^1_p=0.7975 ng/mL/cc.

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



 §.§ Simulation of untreated prostate cancer growth










Figures <ref> and <ref> 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 <cit.>. 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 <ref>B,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 <ref>C,D). These results also align with the metabolic profile of mild prostatic tumors <cit.>, 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 <ref> and <ref>, 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 <cit.>. This phenomenon has also been observed in previous studies <cit.>. 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 <ref>). 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 <ref>B,E), whereas they become thinner in the opposite scenarios (see Figure <ref>C,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 <ref> displays the plots of tumor volume and serum PSA during each simulation of untreated mild PCa growth, while Figure <ref> 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 <ref>.
In the reference case of untreated aggressive PCa, the tumor volume increase slightly slowed down by the end (see Figure <ref>A).
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 <ref>C,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 <ref> and <ref> 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 <cit.>.
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 <ref>B,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 <cit.>.



 §.§ Simulation of treatment plans


All treatment plans successfully removed mild PCa. Figure <ref> 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 <ref>A). Antiangiogenic therapy took comparatively longer to eliminate the tumor, as this required eight doses instead (see Figure <ref>B). 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 <ref>C). 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 <ref>. Figures <ref> and <ref> 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 <ref>A). 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 <cit.>. 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 <ref>B 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 <ref>B 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 <ref>A).

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 <ref>B). 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 <ref>C 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 <ref>B). 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 <cit.>.

Finally, combined therapy succeeded in eliminating the aggressive tumor right after the tenth cycle. Figure <ref>C 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 <ref>D). We believe that the absence of branching in this simulation may be interpreted as an early sign of successful treatment outcome. Additionally, Figure <ref>D 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.



§ 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 <cit.>. 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 <cit.>, which is an extensively used biomarker of PCa in clinical management of the disease <cit.>.

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 ϕ, representing the phase variable, turns out to lie between 0 (healthy phase) and 1 (tumor phase), and the other components σ 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 <cit.>.
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 <cit.>.
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 <cit.>.
Our model was also able to reproduce the poorer correlation between serum PSA and tumor volume in aggressive cancers exhibiting fast tumor dynamics <cit.>.

Additionally, our results align with the use of cytotoxic chemotherapy as the reference approach to treat advanced PCa <cit.>, which may show further benefit from combination with an antiangiogenic drug in order to treat aggressive tumors <cit.>. 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) <cit.>. 
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. β_c and β_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 <cit.>. 
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 <cit.>.
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 <cit.>.
Microvasculature may be modeled following alternative approaches <cit.>. Hybrid models enable a precise description of the morphology and evolution of microvasculature, but they are more computationally expensive <cit.>. 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 <cit.>. 

Additionally, we plan to explore our mathematical model of PCa growth and cytotoxic and antiangiogenic drug effects within the context of optimal control problems <cit.> 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 <cit.>.
Computational models of tumor growth including the effect of mechanical deformation on cancer development have been found to provide superior predictions of pathological outcome <cit.>. 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 <cit.>. 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 <cit.>. 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 <cit.> and the modeling of the tumor-inhibiting effects of the drugs can be further refined <cit.>.
Radiation effects could be included in the tumor phase-field equation as an additional cytotoxic term <cit.> within our modeling framework. 

Finally, we assumed that proliferation, apoptosis, and the drug effect rates β_c and β_a remained constant during the simulation. However, these parameters may evolve due to treatment action and the phenotypic evolution of the tumor <cit.>. 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. 



abbrv

10

Anderson2008
A. R. Anderson and V. Quaranta.
Integrative mathematical oncology.
Nat. Rev. Cancer, 8(3):227–234, 2008.

Antonarakis2012
E. S. Antonarakis and M. A. Carducci.
Targeting angiogenesis for the treatment of prostate cancer.
Expert Opin. Ther. Targets, 16(4):365–376, 2012.

Baker2004
S. D. Baker, M. Zhao, C. K. K. Lee, J. Verweij, Y. Zabelina, J. R. Brahmer,
  A. C. Wolff, A. Sparreboom, and M. A. Carducci.
Comparative pharmacokinetics of weekly and every-three-weeks
  docetaxel.
Clin. Cancer Res., 10(6):1976–1983, 2004.

Barbu2010
V. Barbu.
Nonlinear differential equations of monotone types in Banach
  spaces.
Springer Science & Business Media, New York, USA, 2010.

Benzekry2013
S. Benzekry and P. Hahnfeldt.
Maximum tolerated dose versus metronomic scheduling in the treatment
  of metastatic cancers.
J. Theor. Biol., 335:235–244, 2013.

Benzekry2012
Benzekry, Sébastien.
Mathematical and numerical analysis of a model for anti-angiogenic
  therapy in metastatic cancers.
ESAIM: M2AN, 46(2):207–237, 2012.

Berges1995
R. R. Berges, J. Vukanovic, J. I. Epstein, M. CarMichel, L. Cisek, D. E.
  Johnson, R. W. Veltri, P. C. Walsh, and J. T. Isaacs.
Implication of cell kinetic changes during the progression of human
  prostatic cancer.
Clin. Cancer Res., 1(5):473–480, 1995.

Bodnar2019
M. Bodnar and M. V. Pérez.
Mathematical and numerical analysis of low-grade gliomas model and
  the effects of chemotherapy.
Commun. Nonlinear Sci. Numer. Simulat., 72:552–564, 2019.

Bogdanska2017
M. Bogdańska, M. Bodnar, J. Belmonte-Beitia, M. Murek, P. Schucht, J. Beck,
  and V. Pérez-García.
A mathematical model of low grade gliomas treated with temozolomide
  and its therapeutical implications.
Math. Biosci., 288:1–13, 2017.

Cavaterra2019
C. Cavaterra, E. Rocca, and H. Wu.
Long-time dynamics and optimal control of a diffuse interface model
  for tumor growth.
Appl. Math. Optim., pages 1–49, 2019.
DOI: https://doi.org/10.1007/s00245-019-09562-5.

Caysa2012
H. Caysa, S. Hoffmann, J. Luetzkendorf, L. P. Mueller, S. Unverzagt, K. Mäder,
  and T. Mueller.
Monitoring of xenograft tumor growth and response to chemotherapy by
  non-invasive in vivo multispectral fluorescence imaging.
PLoS ONE, 7(10):e47927, 2012.

Cha2015
E. K. Cha and J. A. Eastham.
Chemotherapy and novel therapeutics before radical prostatectomy for
  high-risk clinically localized prostate cancer.
Urol. Oncol., 33(5):217–225, 2015.

chung1993time
J. Chung and G. Hulbert.
A time integration algorithm for structural dynamics with improved
  numerical dissipation: the generalized-α method.
J. Appl. Mech., 60(2):371–375, 1993.

Colli2015
P. Colli, G. Gilardi, and D. Hilhorst.
On a Cahn-Hilliard type phase field system related to tumor growth.
Discrete Contin. Dyn. Syst. Ser. A, 35(6):2423–2442, 2015.

Colli2019
P. Colli, G. Gilardi, G. Marinoschi, and E. Rocca.
Sliding mode control for a phase field system related to tumor
  growth.
Appl. Math. Optim., 79(3):647–670, 2019.

Colli2017
P. Colli, G. Gilardi, E. Rocca, and J. Sprekels.
Optimal distributed control of a diffuse interface model of tumor
  growth.
Nonlinearity, 30(6):2518–2546, 2017.

Corwin2013
D. Corwin, C. Holdsworth, R. C. Rockne, A. D. Trister, M. M. Mrugala, J. K.
  Rockhill, R. D. Stewart, M. Phillips, and K. R. Swanson.
Toward patient-specific, biologically optimized radiation therapy
  plans for the treatment of glioblastoma.
PLoS ONE, 8(11):e79115, 2013.

Eisenberger2012
M. A. Eisenberger and E. S. Antonarakis.
The experience with cytotoxic chemotherapy in metastatic
  castration-resistant prostate cancer.
Urol. Clin. North Am., 39(4):573–581, 2012.

Erbersdobler2004
A. Erbersdobler, H. Augustin, T. Schlomm, and R.-P. Henke.
Prostate cancers in the transition zone: Part 1; pathological
  aspects.
BJU Int., 94(9):1221–1225, 2004.

ferlay2015cancer
J. Ferlay, I. Soerjomataram, R. Dikshit, S. Eser, C. Mathers, M. Rebelo, D. M.
  Parkin, D. Forman, and F. Bray.
Cancer incidence and mortality worldwide: sources, methods and major
  patterns in globocan 2012.
Int. J. Cancer, 136(5):E359–E386, 2015.

Ferrara2004
N. Ferrara, K. J. Hillan, H.-P. Gerber, and W. Novotny.
Discovery and development of bevacizumab, an anti-VEGF antibody for
  treating cancer.
Nat. Rev. Drug Discov., 3(5):391–400, 2004.

figg2008angiogenesis
W. D. Figg and J. Folkman.
Angiogenesis: an integrative approach from science to medicine.
Springer Science & Business Media, 2008.

Fraldi2018
M. Fraldi and A. R. Carotenuto.
Cells competition in tumor growth poroelasticity.
J. Mech. Phys. Solids, 112:345–367, 2018.

Frieboes2010
H. B. Frieboes, F. Jin, Y.-L. Chuang, S. M. Wise, J. S. Lowengrub, and
  V. Cristini.
Three-dimensional multispecies nonlinear tumor growth–II: Tumor
  invasion and angiogenesis.
J. Theor. Biol., 264(4):1254–1278, 2010.

Frigeri2015
S. Frigeri, M. Grasselli, and E. Rocca.
On a diffuse interface model of tumour growth.
European J. Appl. Math., 26(2):215––243, 2015.

Gallaher2018
J. A. Gallaher, P. M. Enriquez-Navas, K. A. Luddy, R. A. Gatenby, and A. R.
  Anderson.
Spatial heterogeneity and evolutionary dynamics modulate time to
  recurrence in continuous and adaptive cancer therapies.
Cancer Res., 78(8):2127–2139, 2018.

Garcke2018
H. Garcke, K. F. Lam, and E. Rocca.
Optimal control of treatment time in a diffuse interface model of
  tumor growth.
Appl. Math. Optim., 78(3):495–544, 2018.

Garcke2016
H. Garcke, K. F. Lam, E. Sitka, and V. Styles.
A Cahn-Hilliard-Darcy model for tumour growth with chemotaxis and
  active transport.
Math. Models Methods Appl. Sci., 26(06):1095–1148, 2016.

gomez2018computational
H. Gomez and K. G. van der Zee.
Computational phase-field modeling.
Encyclopedia of Computational Mechanics Second Edition, pages
  1–35, 2018.

Gordon2001
M. S. Gordon, K. Margolin, M. Talpaz, G. W. Sledge, E. Holmgren, R. Benjamin,
  S. Stalter, S. Shak, and D. C. Adelman.
Phase I safety and pharmacokinetic study of recombinant human
  anti-vascular endothelial growth factor in patients with advanced cancer.
J. Clin. Oncol., 19(3):843–850, 2001.

Gorelik2008
B. Gorelik, I. Ziv, R. Shohat, M. Wick, W. D. Hankins, D. Sidransky, and
  Z. Agur.
Efficacy of weekly docetaxel and bevacizumab in mesenchymal
  chondrosarcoma: A new theranostic method combining xenografted biopsies with
  a mathematical model.
Cancer Res., 68(21):9033–9040, 2008.

Hahnfeldt2003
P. Hahnfeldt, J. Folkman, and L. Hlatky.
Minimizing long-term tumor burden: The logic for metronomic
  chemotherapeutic dosing and its antiangiogenic basis.
J. Theor. Biol., 220(4):545–554, 2003.

Hahnfeldt1999
P. Hahnfeldt, D. Panigrahy, J. Folkman, and L. Hlatky.
Tumor development under angiogenic signaling.
Cancer Res., 59(19):4770–4775, 1999.

Hanahan2011
D. Hanahan and R. Weinberg.
Hallmarks of cancer: The next generation.
Cell, 144(5):646–674, 2011.

harma2010comprehensive
V. Härmä, J. Virtanen, R. Mäkelä, A. Happonen, J.-P. Mpindi,
  M. Knuuttila, P. Kohonen, J. Lötjönen, O. Kallioniemi, and M. Nees.
A comprehensive panel of three-dimensional models for studies of
  prostate cancer growth, invasion and drug responses.
PLoS ONE, 5(5):e10431, 2010.

Henares-Molina2017
A. Henares-Molina, S. Benzekry, P. C. Lara, M. García-Rojo, V. M.
  Pérez-García, and A. Martínez-González.
Non-standard radiotherapy fractionations delay the time to malignant
  transformation of low-grade gliomas.
PLoS ONE, 12(6):e0178552, 2017.

Hinow2009
P. Hinow, P. Gerlee, L. J. McCawley, V. Quaranta, M. Ciobanu, S. Wang, J. M.
  Graham, B. P. Ayati, J. Claridge, K. R. Swanson, et al.
A spatial model of tumor-host interaction: application of
  chemotherapy.
Math. Biosci. Eng., 6(3):521–546, 2009.

hughes2005isogeometric
T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs.
Isogeometric analysis: CAD, finite elements, NURBS, exact geometry
  and mesh refinement.
Comput. Methods Appl. Mech. Engrg., 194(39-41):4135–4195,
  2005.

Jackson2000
T. L. Jackson and H. M. Byrne.
A mathematical model to study the effects of drug resistance and
  vasculature on the response of solid tumors to chemotherapy.
Math. Biosci., 164(1):17–38, 2000.

Jain2005
R. K. Jain.
Normalization of tumor vasculature: An emerging concept in
  antiangiogenic therapy.
Science, 307(5706):58–62, 2005.

Jain2014
R. K. Jain, J. D. Martin, and T. Stylianopoulos.
The role of mechanical forces in tumor growth and therapy.
Annu. Rev. Biomed. Eng., 16(1):321–346, 2014.

jansen2000generalized
K. E. Jansen, C. H. Whiting, and G. M. Hulbert.
A generalized-α method for integrating the filtered
  Navier–Stokes equations with a stabilized finite element method.
Comput. Methods Appl. Mech. Engrg., 190(3-4):305–319, 2000.

Kelly2012
W. K. Kelly, S. Halabi, M. Carducci, D. George, J. F. Mahoney, W. M. Stadler,
  M. Morris, P. Kantoff, J. P. Monk, E. Kaplan, N. J. Vogelzang, and E. J.
  Small.
Randomized, double-blind, placebo-controlled phase III trial
  comparing docetaxel and prednisone with or without bevacizumab in men with
  metastatic castration-resistant prostate cancer: CALGB 90401.
J. Clin. Oncol., 30(13):1534–1540, 2012.

Kim2005
J. J. Kim and I. F. Tannock.
Repopulation of cancer cells during therapy: an important cause of
  treatment failure.
Nat. Rev. Cancer, 5(7):516–525, 2005.

Kohandel2007
M. Kohandel, M. Kardar, M. Milosevic, and S. Sivaloganathan.
Dynamics of tumor growth and combination of anti-angiogenic and
  cytotoxic therapies.
Phys. Med. Biol., 52(13):3665–3677, 2007.

Ladyzhenskaia1968
O. A. Ladyženskaja, V. A. Solonnikov, and N. N. Ural'ceva.
Linear and quasi-linear equations of parabolic type, volume 23
  of Translations of Mathematical Monographs.
American Mathematical Society, Providence, USA, 1968.

Lima2017
E. Lima, J. Oden, B. Wohlmuth, A. Shahmoradi, D. Hormuth, T. Yankeelov,
  L. Scarabosio, and T. Horger.
Selection and validation of predictive models of radiation effects on
  tumor growth based on noninvasive imaging data.
Comput. Methods Appl. Mech. Engrg., 327(Supplement C):277–305,
  2017.

Lions1961
J. L. Lions.
Équations Différentielles Opérationnelles et
  Problèmes aux Limites.
Springer-Verlag, Berlin/Göttingen/Heidelberg, Germany, 1961.

lorenzo2019computer
G. Lorenzo, T. J. R. Hughes, P. Dominguez-Frojan, A. Reali, and H. Gomez.
Computer simulations suggest that prostate enlargement due to benign
  prostatic hyperplasia mechanically impedes prostate cancer growth.
Proc. Natl. Acad. Sci. U.S.A., 116(4):1152–1161, 2019.

lorenzo2017hierarchically
G. Lorenzo, M. A. Scott, K. Tew, T. J. R. Hughes, and H. Gomez.
Hierarchically refined and coarsened splines for moving interface
  problems, with particular application to phase-field models of prostate tumor
  growth.
Comput. Methods Appl. Mech. Engrg., 319:515–548, 2017.

lorenzo2016tissue
G. Lorenzo, M. A. Scott, K. Tew, T. J. R. Hughes, Y. J. Zhang, L. Liu,
  G. Vilanova, and H. Gomez.
Tissue-scale, personalized modeling and simulation of prostate cancer
  growth.
Proc. Natl. Acad. Sci. U.S.A., 113(48):E7663–E7671, 2016.

Lu2008
J.-F. Lu, R. Bruno, S. Eppler, W. Novotny, B. Lum, and J. Gaudreault.
Clinical pharmacokinetics of bevacizumab in patients with solid
  tumors.
Cancer Chemother. Pharmacol., 62(5):779–786, 2008.

mehta2001independent
R. Mehta, A. Kyshtoobayeva, T. Kurosaki, E. J. Small, H. Kim, R. Stroup, C. E.
  McLaren, K.-T. Li, and J. P. Fruehauf.
Independent association of angiogenesis index with outcome in
  prostate cancer.
Clin. Cancer Res., 7(1):81–88, 2001.

Miranville2019
A. Miranville, E. Rocca, and G. Schimperna.
On the long time behavior of a tumor growth model.
J. Differential Equations, 267(4):2616–2642, 2019.

Mottet2018
N. Mottet, R. van den Bergh, E. Briers, L. Bourke, P. Cornford, M. D. Santis,
  S. Gillessen, A. Govorov, J. Grummet, A. Henry, T. Lam, M. Mason, H. van der
  Poel, T. van der Kwast, O. Rouviére, T. Wiegel, T. V. den Broeck,
  M. Cumberbatch, N. Fossati, T. Gross, M. Lardas, M. Liew, L. Moris,
  I. Schoots, and P. Willemse.
EAU-ESTRO-ESUR-SIOG Guidelines on Prostate Cancer.
European Association of Urology, 2018.

Mpekris2017
F. Mpekris, J. W. Baish, T. Stylianopoulos, and R. K. Jain.
Role of vascular normalization in benefit from metronomic
  chemotherapy.
Proc. Natl. Acad. Sci. U.S.A., 114(8):1994–1999, 2017.

mukherji2013angiogenesis
D. Mukherji, S. Temraz, D. Wehbe, and A. Shamseddine.
Angiogenesis and anti-angiogenic therapy in prostate cancer.
Crit. Rev. Oncol. Hematol., 87(2):122–131, 2013.

noguchi2000assessment
M. Noguchi, T. A. Stamey, J. E. McNeal, and C. E. Yemoto.
Assessment of morphometric measurements of prostate carcinoma volume.
Cancer, 89(5):1056–1064, 2000.

Perez-Garcia2015
V. M. Pérez-García, M. Bogdanska, A. Martínez-González,
  J. Belmonte-Beitia, P. Schucht, and L. A. Pérez-Romasanta.
Delay effects in the response of low-grade gliomas to radiotherapy: a
  mathematical model and its therapeutical implications.
Math. Med. Biol., 32(3):307–329, 2015.

pettaway2000neoadjuvant
C. A. Pettaway, L. L. Pisters, P. Troncoso, J. Slaton, L. Finn, K. Kamoi, and
  C. J. Logothetis.
Neoadjuvant chemotherapy and hormonal therapy followed by radical
  prostatectomy: feasibility and preliminary results.
J. Clin. Oncol., 18(5):1050–1057, 2000.

picus2011phase
J. Picus, S. Halabi, W. K. Kelly, N. J. Vogelzang, Y. E. Whang, E. B. Kaplan,
  W. M. Stadler, E. J. Small, Cancer, and L. G. B.
A phase 2 study of estramustine, docetaxel, and bevacizumab in men
  with castrate-resistant prostate cancer: results from Cancer and Leukemia
  Group B Study 90006.
Cancer, 117(3):526–533, 2011.

Powathil2007
G. Powathil, M. Kohandel, S. Sivaloganathan, A. Oza, and M. Milosevic.
Mathematical modeling of brain tumors: effects of radiotherapy and
  chemotherapy.
Phys. Med. Biol., 52(11):3291–3306, 2007.

reese2001phase
D. M. Reese, P. Fratesi, M. Corry, W. Novotny, E. Holmgren, and E. J. Small.
A phase II trial of humanized anti-vascular endothelial growth
  factor antibody for the treatment of androgen-independent prostate cancer.
Prostate J., 3(2):65–70, 2001.

Roose2003
T. Roose, P. A. Netti, L. L. Munn, Y. Boucher, and R. K. Jain.
Solid stress generated by spheroid growth estimated using a linear
  poroelasticity model.
Microvasc. Res., 66(3):204–212, 2003.

Saad1986
Y. Saad and M. H. Schultz.
GMRES: A generalized minimal residual algorithm for solving
  nonsymmetric linear systems.
SIAM J. Sci. and Stat. Comput., 7(3):856–869, 1986.

Schmid1993
H.-P. Schmid, J. E. McNeal, and T. A. Stamey.
Observations on the doubling time of prostate cancer. the use of
  serial prostate-specific antigen in patients with untreated disease as a
  measure of increasing cancer volume.
Cancer, 71(6):2031–2040, 1993.

Seruga2011
B. Seruga and I. F. Tannock.
Chemotherapy-based treatment for castration-resistant prostate
  cancer.
J. Clin. Oncol., 29(27):3686–3694, 2011.

Small2012
A. C. Small and W. K. Oh.
Bevacizumab treatment of prostate cancer.
Expert Opin. Biol. Ther., 12(9):1241–1249, 2012.

swanson2001quantitative
K. R. Swanson, L. D. True, D. W. Lin, K. R. Buhler, R. Vessella, and J. D.
  Murray.
A quantitative model for the dynamics of serum prostate-specific
  antigen as a marker for cancerous growth: an explanation for a medical
  anomaly.
Am. J. Pathol., 158(6):2195–2199, 2001.

Tije2005
A. J. ten Tije, J. Verweij, M. A. Carducci, W. Graveland, T. Rogers, T. Pronk,
  M. Verbruggen, F. Dawkins, and S. D. Baker.
Prospective evaluation of the pharmacokinetics and toxicity profile
  of docetaxel in the elderly.
J. Clin. Oncol., 23(6):1070–1077, 2005.

Trock2011
B. J. Trock.
Application of metabolomics to prostate cancer.
Urol. Oncol., 29(5):572–581, 2011.

Vilanova2017
G. Vilanova, I. Colominas, and H. Gomez.
Computational modeling of tumor-induced angiogenesis.
Arch. Computat. Methods Eng., 24(4):1071–1102, 2017.

Vollmer2010
R. T. Vollmer.
Dissecting the dynamics of serum prostate-specific antigen.
Am. J. Clin. Pathol., 133(2):187–193, 2010.

weidner1993tumor
N. Weidner, P. Carroll, J. Flax, W. Blumenfeld, and J. Folkman.
Tumor angiogenesis correlates with metastasis in invasive prostate
  carcinoma.
Am. J. Pathol., 143(2):401–409, 1993.

Weis2015
J. A. Weis, M. I. Miga, L. R. Arlinghaus, X. Li, V. Abramson, A. B.
  Chakravarthy, P. Pendyala, and T. E. Yankeelov.
Predicting the response of breast cancer to neoadjuvant therapy using
  a mechanically coupled reaction–diffusion model.
Cancer Res., 75(22):4697–4707, 2015.

Wise2008
S. Wise, J. Lowengrub, H. Frieboes, and V. Cristini.
Three-dimensional multispecies nonlinear tumor growth-I: Model and
  numerical method.
J. Theor. Biol., 253(3):524–543, 2008.

Xu2016
J. Xu, G. Vilanova, and H. Gomez.
A mathematical model coupling tumor growth and angiogenesis.
PLoS ONE, 11(2):e0149422, 2016.

Yankeelov2013
T. E. Yankeelov, N. Atuegwu, D. Hormuth, J. A. Weis, S. L. Barnes, M. I. Miga,
  E. C. Rericha, and V. Quaranta.
Clinically relevant modeling of tumor growth and treatment response.
Sci. Transl. Med., 5(187):187ps9, 2013.




