# Combined Hormone and Brachy Therapies for the Treatment of Prostate Cancer

Salma Chabbar, Abderrahmane Habbal, Rajae Aboulaich, Nabil Ismaili,
Sanaa El Majjaoui

► To cite this version:

Salma Chabbar, Abderrahmane Habbal, Rajae Aboulaich, Nabil Ismaili, Sanaa El Majjaoui. Combined Hormone and Brachy Therapies for the Treatment of Prostate Cancer. Mathematical Modelling of Natural Phenomena, 2023, 10.1051/mmnp/. [PDF]. hal-03936301

HAL Id: hal-03936301

https://inria.hal.science/hal-03936301

Submitted on 12 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# COMBINED HORMONE AND BRACHY THERAPIES FOR THE TREATMENT OF PROSTATE CANCER *

SALMA CHABBAR<fnref n="1" />,<fnref n="2" />, ABDERRAHMANE HABBAL<fnref n="2" />,<fnref n="3" />, RAJAE ABOULAICH<fnref n="1" />, NABIL ISMAILI<fnref n="4" /> AND SANAA EL MAJJAQUI<fnref n="5" />

**Abstract.** Prostate cancer is a hormone-dependent cancer characterized by two types of cancer cells, androgen-dependent cancer cells and androgen-resistant ones. The objective of this paper is to present a novel mathematical model for the treatment of prostate cancer under combined hormone therapy and brachytherapy. Using a system of partial differential equations, we quantify and study the evolution of the different cell densities involved in prostate cancer and their responses to the two treatments. Numerical simulations of tumor growth under different therapeutic strategies are explored and presented. The numerical simulations are carried out on FreeFem++ using a 2D finite element method.

**Mathematics Subject Classification.** ???, ???

...

## INTRODUCTION

The latest statistics published by the International Agency for Research on Cancer show that in 2018, 18.1 million new cancer cases have been identified and 9.6 million deaths have been recorded worldwide making it the second leading cause of death globally [5]. Prostate cancer ranks third in incidence with 1.28 million cases and represents the second most commonly diagnosed male cancer[5].

Prostate cells need the hormone androgen to survive and function properly. For this to happen, the androgens have to bind to a protein in the prostate cells called Androgen Receptor and activate it [7, 10, 18, 19]. Since androgens act as a growth factor for the cells, one way of treating prostate cancer is through the antihormone therapy that hinder its activity. The Androgen Deprivation Therapy (ADT) aims to either reduce androgen production or to stop the androgens from working through the use of drugs. However, over time, castration-resistant

*Keywords and phrases.* Tumor growth, PDE model, Numerical simulations, Androgen deprivation therapy, Brachytherapy, Combined treatment.

<footnote>* Le Centre National pour la Recherche Scientifique et Technique (CNRST)</footnote>

<footnote n="1">¹ Mohammadia School of Engineering, Mohammed V University in Rabat, LERMA, Avenue Ibn Sina B.P 765, Agdal Rabat, 10090, Morocco</footnote>

<footnote n="2">² University of Côte d'Azur, Inria, LJAD, Parc Valrose 06108 Nice, France</footnote>

<footnote n="3">³ Polytechnic Mohammed VI University, Benguerir, Morocco</footnote>

<footnote n="4">⁴ Department of Medical Oncology, Mohammed VI University of Health Sciences, Cheikh Khalifa International University Hospital, and Mohammed VI International University Hospital, Casablanca, Morocco</footnote>

<footnote n="5">⁵ Department of Radiotherapy, National Institute of Oncology, Avenue Allal El Fassi, Rabat, Morocco</footnote>
cells that are able to sustain growth in a low androgen environment emerge. The castration-resistant cells can either be androgen independent or androgen repressed meaning that they have a negative growth rate when the androgen is abundant in the prostate[16, 17, 19, 24, 33]. In order to delay the development of castration-resistance and reduce its occurrence, the Intermittent Androgen Deprivation Therapy is used. In this case, the antihormone therapy is given in cycles by alternating drug administrations and rest periods. Switching between the different periods is done according to the level of the Prostate Specific Antigen (PSA), a protein produced by prostate cells and present in small quantities in the serum of men with healthy prostates[10, 15, 19].

Several studies have compared the continuous ADT with the intermittent one, and all have concluded that with the continuous ADT, the patient will eventually relapse due to the emergence and proliferation of the castration-resistant prostate cancer cell. In [19], conditions that favor intermittent therapy over continuous treatment in terms of delaying the emergence of castration resistance or eventually achieving remission have been established. In [10], while continuous ADT leads to treatment failure in finite time, the effectiveness of the intermittent treatment depends on the characteristics of the castration-resistant cancer cells and the scheduling of the therapy with the possibility of a scenario where the tumor volume can be bounded. While in [18] the treatment response depends on the competitive interaction between the androgen-dependent and the castration-resistant cancer cells. In the case where the androgen-dependent cells had a competitive advantage over the castration-resistant ones, the intermittent scheduling yielded better results than the continuous treatment. Otherwise, the intermittent ADT resulted in an earlier treatment failure. In addition to the risk of relapse, the androgen deprivation therapy causes some side effects that affect greatly the quality of life of the patient [24, 30]. Combining the ADT with another type of treatment may prove to be more effective and less harmful to the patient.

The brachytherapy is an effective radiation therapy used in the treatment of prostate cancer by placing a sealed radiation source inside the prostate gland[21, 22, 35]. It can be delivered in high dose rates (HDR) or low dose rates (LDR) depending on the radioactive source used and the duration of treatment. In the HDR brachytherapy the source is placed temporarily in the prostate for a few minutes to deliver high dose radiation while for the LDR brachytherapy low radiation dose are delivered from radioactive sources permanently placed in the prostate. The radioactivity of the source decays over time, therefore its presence in the prostate does not cause any long-term concern as its radioactivity disappears eventually[14, 20]. In the present paper, the LDR brachytherapy is investigated for the treatment of prostate cancer using the Iodine-125 as a decaying radioactive source. In practice, brachytherapy is prescribed either as monotherapy, often for localized tumors, or combined with another therapy such as external beam radiation therapy for which the total dose prescribed is divided between internal and external radiation. Brachytherapy can also be prescribed in combination with hormone therapy. However, in the existing literature there is currently no mathematical model that explores this combination of treatments. The objective of this work is to develop a mathematical model to assess the effectiveness of combining androgen deprivation therapy with brachytherapy in the treatment of prostate cancer. The resulting simulations can be used to explore potential therapeutic strategies.

In the next section, we formulate the mathematical model for the treatment of prostate cancer under combined hormone therapy and brachytherapy.

# 1. MATHEMATICAL MODELING OF PROSTATE CANCER UNDER COMBINED THERAPIES

To model the growth and spread of cancer cells in a patient with prostate cancer, we propose a spatial model inspired by the one presented by Friedman A. et al. in [10]. To describe our model, we first note that the volume of the tumor changes with time. We, thus, assume that the tumor is occupying a region $\Sigma(t)$ at time $t$ in a bounded domain $\Omega$ containing it.

The evolution of the prostate tumor is described through the evolution of the densities of the different cells composing it. The growth and spread of prostate cells are presented in the following section.
## 1.1. Dynamics of the prostate cells

In this section, the considered dynamics governing cell behavior are discussed and modeled. First regarding the movement of the cells and second regarding their growth.

### 1.1.1. Movements of the prostate cells

The cells are surrounded and supported by the extracellular matrix (ECM) which, as an adhesive substrate, plays an important role in cell adhesion, movement and motility [2, 9, 25].

In our work, we consider two different types of cell movement: a random movement of the cells and a directional cell movement induced by the adhesion gradient of the ECM. In this case the cells move in the direction of a positive gradient of the chemicals bound to the ECM. This directional motility also known as haptotaxis is involved in several biological processes such as wound healing and tumor cell invasion [29].

Metastasis or tumor cell invasion is one of the most dangerous aspects of cancer, as the tumor is no longer localized but rather spreads to other locations. This phenomenon occurs due to the degradation of the ECM by cancer cells and was widely studied and taken into consideration while modeling tumor evolution [2, 11, 25, 34, 37].

By degrading their support, cancer cells are able to migrate and spread. To do so, they produce the matrix degradative enzymes (MDEs) that erode the ECM locally by altering cell adhesion. This also enables the cells to make a space where they can move by diffusion (random motility).

In the present work, we used the simplified method provided by Alina Toma et al. in [34] to model tumor growth which does not explicitly model the MDEs and where the degraded ECM can be simulated directly using the cancer cell distribution. The ECM $f$ is, thus, directly degraded by cancer cells with a coefficient $\alpha_f$ and is regenerating with a rate $\beta_f$.

### 1.1.2. Growth of the prostate cells

Since androgen acts as a growth factor for the cells, the evolution of prostate cells depends on the level of androgen available in their environment. In the present work, we used the same growth rate profiles for androgen-dependent and castration-resistant cells as presented in [10, 15, 17, 33].

The growth rate $K_{A_d}$ of the androgen dependent cells $A_d$ is an increasing and saturating function of the androgen level $A$ defined as follows:

$$K_{A_d}(A) = \underbrace{\lambda_{A_d} \left( \lambda_1 + (1 - \lambda_1) \frac{A}{A + k_1} \right)}_{\text{Proliferation rate}} - \underbrace{\mu_{A_d} \left( \mu_1 + (1 - \mu_1) \frac{A}{A + k_1} \right)}_{\text{Death rate}} \quad (1)$$

The growth rate $K_{A_r}$ of castration-resistant cells $A_r$ is a decreasing function of androgen defined as follows:

$$K_{A_r}(A) = \underbrace{\lambda_{A_r}}_{\text{Proliferation rate}} - \underbrace{\mu_{A_r} \left( \mu_2 + (1 - \mu_2) \frac{A}{A + k_2} \right)}_{\text{Death rate}} \quad (2)$$

In an environment full of androgen the androgen dependent cells $A_d$ proliferate at a maximum rate $\lambda_{A_d}$ and die at a rate $\mu_{A_d}$ while castration-resistant cells are considered androgen-repressed and have a negative growth rate in the presence of androgen. As in [17], their proliferation rate is considered androgen independent and has a constant value of $\lambda_{A_r}$ while their death rate is an increasing and saturating function of androgen. $\lambda_1$ quantifies the impact of androgen concentration on $A_d$ proliferation. It represents the $A_d$ proliferation rate in the absence of androgens compared to its value under abundant androgens. A value of 0 as in [15] means that $A_d$ cells do not multiply without androgens while a value of 0.8 as in [17] means that $A_d$ cells lose 20% of their maximum proliferation rate when there is no more androgen in the medium. Likewise, $\mu_1$ and $\mu_2$ quantify the impact of androgens concentration on $A_d$ and $A_r$ death respectively. They represent their death rate in the absence of androgens compared to its value under abundant androgens such that $\mu_1 > 1$ to increase $A_d$ death in a low
androgen environment and $\mu_2 < 1$ to lower $A_r$ death in the unfavorable environment.

Moreover, in order to survive in an androgen-depleted environment, $A_d$ cells mutate into $A_r$ cells with a rate $K_m$. The more the environment is unfavorable for the $A_d$ cells, the more likely they will mutate into $A_r$ to survive [8]. The mutation rate $K_m$ is therefore a decreasing function of androgen. The mutation is considered irreversible for simplicity [10, 15, 33] and the mutation rate is defined as:

$$K_m(A) = K_{max} \left(1 - \frac{A}{A_0}\right) \qquad (3)$$

where $A_0$ is the normal level of androgen concentration in healthy condition and $K_{max}$ is the maximum mutation rate.

Normal cell proliferation is regulated through contact inhibition. This process allows a healthy cell to stop its own proliferation when it comes into contact with another cell. Cancer cells do not have this characteristic, resulting in uncontrolled and excessive cancer cell proliferation and solid tumor development [26]. Compared to cancer cells, healthy cells $H$ are assumed to multiply at a constant growth rate $K_H$ while considering the maximum cell carrying capacity $K_0$ [10].

## 1.2. Mathematical modeling of both therapies

In this section, we are interested in modeling the effects of the hormone therapy and brachytherapy and the response of prostate cells to both therapies.

### 1.2.1. Hormone therapy

Androgen deprivation therapy (ADT) aims to deplete the prostate of androgens or to prevent them from getting into prostate cells in order to kill the cancer cells that rely on it to survive and multiply. We model the case when drugs that directly interfere with the production of androgens in the prostate are used.

Throughout the treatment period, the ADT can be delivered continuously or intermittently, alternating periods of treatment with periods of rest. The intermittent androgen suppression (IAS) method depends on the PSA levels to control the treatment schedule.

The prostate-specific antigen (PSA) is a protein produced by prostate cells and is normally found at very low levels in the blood of men with healthy prostate and is used as a biomarker for monitoring prostate cancer [1]. The IAS is switched on when the PSA level is above a certain critical level $P_0$ and is switched off when the PSA level goes down below the critical level $P_0$. This results in depletion of androgen production and this decrease is modeled by reduction in the normal rate of androgen production by a factor of $1 + \lambda_u u(P)$. Here, $u(\cdot)$ is a smooth function and $P$ represents the PSA level such that $u(P) = 1$ when the IAS is switched on and $u(P) = 0$ when the IAS is switched off. For simplicity, and as in [10, 15, 28], the PSA concentration $P$ is considered proportional to the total number of cancer cells as follows:

$$P = P(t) = \int_{\Omega} \beta_{PSA}(A_d + A_r) dx \qquad (4)$$

The diffusion of the androgen concentration $A$ in the prostate under hormone therapy is therefore defined as follows:

$$\frac{\partial A}{\partial t} - D_A \Delta A = \underbrace{\frac{\lambda_A}{1 + \lambda_u u(P, t)}}_{\text{drug-inhibited production}} - \underbrace{\mu_A A}_{\text{natural decay}} - \underbrace{\lambda_0 (A_d + A_r) A}_{\text{consumption by } A_d \text{ and } A_r} \qquad (5)$$

Eq.(5) can easily be adapted to model the second type of drug used in hormone therapy that prevents androgens from working instead of interfering with their production. This can be modeled by targeting the existing androgens, thereby increasing their decay rate instead of reducing new androgen production.
### 1.2.2. Brachytherapy

In radiation therapy, X-rays or photons travel through the patient's tissue to reach the tumor. In brachytherapy the decaying radioactive sources are placed inside the organ and release energetic beams whose energy depends on the nature of the source used. As the photons travel through the medium, their energy decreases with an attenuation coefficient $k$[20]. $k$ quantifies the energy loss caused by the distance between the irradiated area and the radioactive source. Its value depends on the energy of the photons and the medium crossed[3, 27]. Moreover, the radioactivity of the source declines over time with a decay constant $\lambda$.

The linear quadratic LQ model is generally used to evaluate the biological effect of a radiation dose $D$ on the irradiated tissue [6, 23, 36]. The survival probability $S$ of cells exposed to a radiation dose $D$ is given by:

$$S = \exp(-(\alpha D + \beta GD^2)) \qquad (6)$$

where $\alpha$ and $\beta$ describe the radiosensitivity of the tissue. $\alpha$ is the yield rate for lethal lesions, and $\beta$ for sublethal-reparable lesions. The parameter $G$ is a dose-rate factor which depends on the radiation course.

We consider $M$ radiative source locations and denote by $x_m$ the locations of the sources, $1 \le m \le M$.

During the radiation course $[0, \tau]$, the irradiation by the source $m$ is given at a dose rate $R_m(t) = dD_m(t)/dt$ such that:

$$R_m(t) = R_{m0} \exp(-\lambda t)$$

where $R_{m0}$ is the source $m$ initial dose rate and $\lambda$ its radioactive decay constant.

The radiation dose rate $R_S(x, t)$ absorbed at a given $x \in \Omega$ at time $t$ from the source $m$ is given by:

$$R_S(x, t) = R_m(t) \exp(-k \|x - x_m\|)$$

where $k$ is the attenuation factor.

The total absorbed dose rate $R_T$ is given by:

$$R_T(x, t) = \sum_m R_m(t) \exp(-k \|x - x_m\|). \qquad (7)$$

Radiotherapy does not distinguish between healthy and cancerous cells. The target cells in brachytherapy are therefore $A_d$, $A_r$ and $H$. However, normal healthy cells are less sensitive to radiation than cancer cells. Their DNA repair capacity from sublethal lesions is generally greater. Therefore, we mitigate the brachytherapy killing rate of healthy cells by a factor $\varepsilon$ compared to cancer cells.

The brachytherapy destroys the cells with a space and time dependent killing rate, $\mu_B(x, t)$ such that the death rate from a total dose rate delivered at time $t$ is:

$$\mu_B(x, t) = 1 - S = 1 - \exp(-(\alpha R_T(x, t) + \beta GR_T(x, t)^2)) \qquad (8)$$

## 2. A PDE MODEL FOR PROSTATE CANCER UNDER COMBINED TREATMENT

We summarize the dynamics of cancer cells, healthy cells and androgens in the context of prostate cancer under hormone therapy and brachytherapy along with the different variables used in the following section.

### 2.1. Prostate cancer mathematical model

The variables used to describe the evolution of prostate cancer are listed as follows:

(1) $A_d$ = Density of androgen-dependent cancer cells,

(2) $A_r$ = Density of androgen-repressed cancer cells,

(3) $H$ = Density of healthy cells,

(4) $f$ = Density of the ECM,

(5) $A$ = Concentration of androgens,
(6) $P$ = Concentration of prostate-specific antigen (PSA).

We have the following system of partial differential equations modeling the dynamics of cells, ECM and the androgen concentration:

$$ 
\left\{
\begin{array}{ll}
\partial_t A_d + \chi \nabla \cdot (A_d \nabla f) - D \Delta A_d = K_{A_d}(A) A_d - K_m(A) A_d - \mu_B(x,t) A_d & \text{on } \Omega \times [0, T] \\
\partial_t A_r + \chi \nabla \cdot (A_r \nabla f) - D \Delta A_r = K_{A_r}(A) A_r + K_m(A) A_d - \mu_B(x,t) A_r & \text{on } \Omega \times [0, T] \\
\partial_t H + \chi \nabla \cdot (H \nabla f) - D \Delta H = K_H H \left(1 - \frac{A_d + A_r + H}{K_0}\right) - \varepsilon \mu_B(x,t) H & \text{on } \Omega \times [0, T] \\
\partial_t f = -\alpha_f (A_d + A_r) f + \beta_f f & \text{on } \Omega \times [0, T] \\
\partial_t A - D_A \Delta A = \frac{\lambda_A}{1 + \lambda_u u(P,t)} - \mu_A A - \lambda_0 (A_d + A_r) A & \text{on } \Omega \times [0, T]
\end{array}
\right.
\quad (9) \\
\text{with boundary conditions :}
\quad
\begin{array}{ll}
A_d = 0 & \text{on } \partial \Omega \times [0, T] \\
A_r = 0 & \text{on } \partial \Omega \times [0, T] \\
H = 0 & \text{on } \partial \Omega \times [0, T] \\
f = f_0 & \text{on } \partial \Omega \times [0, T] \\
A = A_0. & \text{on } \partial \Omega \times [0, T]
\end{array}
$$

Boundary equations for $A_d$, $A_r$, $f$ and $A$ are derived from the observation that the outer medium is assumed to be healthy, thus containing neither type of cancer cell, and with a constant androgen concentration and ECM density. The boundary condition on the healthy prostate cells $H$ is derived from the prostate being enclosed in $\Omega$. The known parameters in the above PDE model are described in Table.1.

## 2.2. Numerical evaluation of the PDE system parameters

In order to carry out the simulations, the numerical values listed in Table.2 were used. In addition, for the parameters related to the androgen evolution, the value of the production and natural decay rates of androgen are selected so as to respect the saturating level of androgens $A_0$ which represents the natural level of androgen in healthy conditions and is given by $A_0 = \lambda_A/\mu_A[10]$. With $\mu_A$ being of about $10^{-2}day^{-1}$ [33] we used $\lambda_A = 0.2nMday^{-1}$ and $\mu_A = 10^{-2}day^{-1}$. To have a reference value for $\lambda_0$ we took into consideration the carrying capacity $K_0$ ($\sim 10^7 \#cells/cm^3$) and the natural decay rate of androgen ($\sim 10^{-2}day^{-1}$). The lambda used in the simulations is therefore $\lambda_0 = 10^{-9} \#cells^{-1}cm^3day^{-1}$.

To compute the rate of PSA production by cancer cells, we relied on the results in [32] who quantified the PSA production rate per tumor volume: $\beta_{PSA} = \beta_c * Volume_{cells}/Volume_{prostate}$ where $Volume_{cells} = 1.5 10^{-8}cm^3[2]$.

Normal prostate cells have a rate of division almost equal to their rate of apoptosis, so their growth rate is extremely low[4]. In this work we used $K_H = 0.003day^{-1}$.

For the LDR brachytherapy, the source used in this work is the Iodine 125 (I-125). It has a half-life of ($\sim 60$ days), which is the time the source needs to lose about half of its radioactive activity. The decay constant of I-125 is thus $\lambda = 0.011day^{-1}$ [21]. In terms of photon absorption, the human tissues are comparable to the "water" medium and the I-125 source produces photons with a low energy of ($\sim 28keV$), the attenuation coefficient $k$ used is thus $k = 0.409cm^{-1}$[7]. The attenuation coefficient is a decreasing function of photon energy, the photons produced in low-dose brachytherapy therefore have a high attenuation coefficient which allows for a better preservation of the surrounding healthy tissues.
<table><caption>TABLE 1. Parameters description</caption><thead><tr><th>Parameter</th><th>Description</th><th>Comment</th></tr></thead><tbody><tr><td>χ</td><td>Haptotaxis coefficient</td><td>Constant related to cell-ECM adhesion</td></tr><tr><td>D</td><td>Cell diffusion coefficient</td><td>Responsible for the non motile part of the cell migration by inducing a progressive wavefront</td></tr><tr><td>K<sub>A<sub>d</sub></sub></td><td>A<sub>d</sub> growth rate</td><td>Increasing and saturating function of A. Negative when A is small and tends to a positive constant when A is large.</td></tr><tr><td>K<sub>A<sub>r</sub></sub></td><td>A<sub>r</sub> growth rate</td><td>Decreasing function of A. Negative when A is large.</td></tr><tr><td>K<sub>m</sub></td><td>A<sub>d</sub> mutation rate</td><td>Irreversible change in genetic sequence.</td></tr><tr><td>μ<sub>B</sub></td><td>Brachytherapy killing rate</td><td>Biological effect of brachytherapy on prostate cells</td></tr><tr><td>K<sub>0</sub></td><td>Maximum cell carrying capacity</td><td>Capacity in number of cells per cm<sup>3</sup></td></tr><tr><td>K<sub>H</sub></td><td>H growth rate</td><td>Positive constant</td></tr><tr><td>α<sub>f</sub></td><td>ECM degradation coefficient</td><td>Quantifies the effect of cancer cells on ECM degradation</td></tr><tr><td>β<sub>f</sub></td><td>ECM remodeling rate</td><td>Natural regeneration of ECM</td></tr><tr><td>D<sub>A</sub></td><td>Androgen diffusion coefficient</td><td>Responsible for the diffusion of androgen in the cells environment</td></tr><tr><td>u</td><td>Hormone therapy switch</td><td>The function that activates/deactivates the hormone treatment in case of intermittent schedule</td></tr><tr><td>λ<sub>u</sub></td><td>Hormone therapy intensity</td><td>Coefficient to control the intensity of the androgen deprivation</td></tr><tr><td>λ<sub>A</sub></td><td>Androgens production rate</td><td>Normal androgen concentration for adult males</td></tr><tr><td>μ<sub>A</sub></td><td>Androgens natural decay rate</td><td>Androgen dynamics time constant</td></tr><tr><td>λ<sub>0</sub></td><td>Androgens consumption rate</td><td>Consumption rate by cancer cells</td></tr><tr><td>P</td><td>PSA concentration</td><td>Prostate Specific Antigen concentration in blood</td></tr></tbody></table>

### 3. NUMERICAL RESOLUTION OF THE PDE SYSTEM

The PDE model (9) is solved using the general purpose finite element package FreeFem++ [12, 13]. To do so, we need to write the weak formulation for each variable. We consider the discretization in time:

$$t_0 = 0; t_1 = t_0 + dt = dt; \dots; t_n = ndt; \dots; t_N = Ndt = T.$$

We use an implicit scheme for the discretization in time of all the equations of except for the healthy cells PDE for which a semi-implicit scheme is used. The system (9) is therefore discretized in time as follows:
<table><caption>TABLE 2. Numerical values of the parameters and their sources</caption><thead><tr><th>Parameter</th><th>Value</th><th>Unit</th><th>Source</th></tr></thead><tbody><tr><td>χ</td><td>2600</td><td>cm<sup>2</sup>M<sup>-1</sup>s<sup>-1</sup></td><td>[2]</td></tr><tr><td>D</td><td>10<sup>-9</sup></td><td>cm<sup>2</sup>s<sup>-1</sup></td><td>[2]</td></tr><tr><td>λ<sub>A<sub>d</sub></sub></td><td>0.4621</td><td>day<sup>-1</sup></td><td>[17]</td></tr><tr><td>λ<sub>A<sub>r</sub></sub></td><td>0.4621</td><td>day<sup>-1</sup></td><td>[17]</td></tr><tr><td>μ<sub>A<sub>d</sub></sub></td><td>0.3812</td><td>day<sup>-1</sup></td><td>[17]</td></tr><tr><td>μ<sub>A<sub>r</sub></sub></td><td>0.4765</td><td>day<sup>-1</sup></td><td>[17]</td></tr><tr><td>λ<sub>1</sub></td><td>0 - 1</td><td>-</td><td>[33]</td></tr><tr><td>μ<sub>1</sub></td><td>1.35</td><td>-</td><td>[17]</td></tr><tr><td>μ<sub>2</sub></td><td>0.25 - 1.0</td><td>-</td><td>[17]</td></tr><tr><td>K<sub>max</sub></td><td>5 × 10<sup>-5</sup></td><td>day<sup>-1</sup></td><td>[33]</td></tr><tr><td>A<sub>0</sub></td><td>20</td><td>nM</td><td>[33]</td></tr><tr><td>K<sub>0</sub></td><td>6.7 × 10<sup>7</sup></td><td>#cells/cm<sup>3</sup></td><td>[2]</td></tr><tr><td>α<sub>f</sub></td><td>1.49 × 10<sup>-14</sup></td><td>#cells<sup>-1</sup>cm<sup>3</sup>s<sup>-1</sup></td><td>[34]</td></tr><tr><td>β<sub>f</sub></td><td>10<sup>-7</sup></td><td>s<sup>-1</sup></td><td>[34]</td></tr><tr><td>D<sub>A</sub></td><td>10<sup>-11</sup> - 10<sup>-13</sup></td><td>m<sup>2</sup>s<sup>-1</sup></td><td>[31]</td></tr><tr><td>P<sub>0</sub></td><td>10</td><td>ng/ml</td><td>[33]</td></tr><tr><td>β<sub>c</sub></td><td>1.7210 - 6.9722</td><td>ng/mm<sup>3</sup>/day</td><td>[32]</td></tr><tr><td>λ</td><td>0.011</td><td>day<sup>-1</sup></td><td>[36]</td></tr><tr><td>α</td><td>0.15</td><td>Gy<sup>-1</sup></td><td>[36]</td></tr><tr><td>β</td><td>0.048</td><td>Gy<sup>-2</sup></td><td>[36]</td></tr></tbody></table>

$$
\left\{
\begin{aligned}
& \frac{A_d^{n+1} - A_d^n}{dt} + \chi \nabla \cdot (A_d^{n+1} \nabla f^n) - D \Delta A_d^{n+1} = K_{A_d}(A^n) A_d^{n+1} - K_m(A^n) A_d^{n+1} \\
& - \mu_B A_d^{n+1}, \\
& \frac{A_r^{n+1} - A_r^n}{dt} + \chi \nabla \cdot (A_r^{n+1} \nabla f^n) - D \Delta A_r^{n+1} = K_{A_r}(A^n) A_r^{n+1} + K_m(A^n) A_d^{n+1} \\
& - \mu_B A_r^{n+1}, \\
& \frac{H^{n+1} - H^n}{dt} + \chi \nabla \cdot (H^{n+1} \nabla f^n) - D \Delta H^{n+1} = K_H H^{n+1} \left(1 - \frac{A_d^{n+1} + A_r^{n+1} + H^n}{K_0}\right) \\
& - \varepsilon \mu_B(x, t) H^{n+1}, \\
& \frac{f^{n+1} - f^n}{dt} = -\alpha_f (A_d^{n+1} + A_r^{n+1}) f^{n+1} + \beta_f f^{n+1}, \\
& \frac{A^{n+1} - A^n}{dt} - D_A \Delta A^{n+1} = \frac{\lambda_A}{1 + \lambda_u u} - \mu_A A^{n+1} - \lambda_0 (A_d^{n+1} + A_r^{n+1}) A^{n+1}, \\
& + B.C. + I.C.
\end{aligned}
\right.
\quad (10)
$$

In the brachytherapy killing rate expression Eq.(8) we recall that the parameter *G* represents the dose protraction factor. For the external-beam radiation therapy where the total radiation dose is split into *n* individual dose over the treatment period, *G* = 1/*n* [36]. In the present work, although permanent implants are used, the discretization in time leads to the splitting of the total dose into (*T*/*dt*) doses where *dt* is the time step chosen in the simulations. For this reason we used *G* = 1/(*T*/*dt*) in the simulations. The weak formulations of the problem (10) are written as follows:
$$ \left\{ \begin{array}{l} \text{Find } A_d^{n+1} \in H_0^1(\Omega) \text{ such that } \forall \varphi \in H_0^1(\Omega) \\ \int_{\Omega} \frac{A_d^{n+1} - A_d^n}{dt} \varphi - \chi \int_{\Omega} A_d^{n+1} \nabla f^n \cdot \nabla \varphi + D \int_{\Omega} \nabla A_d^{n+1} \nabla \varphi - \int_{\Omega} K_{A_d}(A^n) A_d^{n+1} \varphi \quad (11) \\ + \int_{\Omega} K_m(A^n) A_d^{n+1} \varphi + \int_{\Omega} \mu_B A_d^{n+1} \varphi = 0, \end{array} \right. $$

$$ \left\{ \begin{array}{l} \text{Find } A_r^{n+1} \in H_0^1(\Omega) \text{ such that } \forall \varphi \in H_0^1(\Omega) \\ \int_{\Omega} \frac{A_r^{n+1} - A_r^n}{dt} \varphi - \chi \int_{\Omega} A_r^{n+1} \nabla f^n \cdot \nabla \varphi + D \int_{\Omega} \nabla A_r^{n+1} \nabla \varphi - \int_{\Omega} K_{A_r}(A^n) A_r^{n+1} \varphi \quad (12) \\ - \int_{\Omega} K_m(A^n) A_d^{n+1} \varphi + \int_{\Omega} \mu_B A_r^{n+1} \varphi = 0, \end{array} \right. $$

$$ \left\{ \begin{array}{l} \text{Find } H^{n+1} \in H_0^1(\Omega) \text{ such that } \forall \varphi \in H_0^1(\Omega) \\ \int_{\Omega} \frac{H^{n+1} - H^n}{dt} \varphi - \chi \int_{\Omega} H^{n+1} \nabla f^n \cdot \nabla \varphi + D \int_{\Omega} \nabla H^{n+1} \nabla \varphi \quad (13) \\ - \int_{\Omega} K_H H^{n+1} \left( 1 - \frac{A_d^{n+1} + A_r^{n+1} + H^n}{K_0} \right) + \varepsilon \int_{\Omega} \mu_B H^{n+1} \varphi = 0, \end{array} \right. $$

$$ \left\{ \begin{array}{l} \text{Find } A^{n+1} \in H^1(\Omega) \text{ such that } \forall \varphi \in H_0^1(\Omega) \\ \int_{\Omega} \frac{A^{n+1} - A^n}{dt} \varphi + D_A \int_{\Omega} \nabla A^{n+1} \nabla \varphi - \int_{\Omega} \frac{\lambda_A}{1 + \lambda_u u} \varphi + \mu_A \int_{\Omega} A^{n+1} \varphi \quad (14) \\ + \lambda_0 \int_{\Omega} (A_d^{n+1} + A_r^{n+1}) A^{n+1} \varphi = 0, \end{array} \right. $$

The computation of $f^{n+1}$ is done as follows:

$$ f^{n+1} = \frac{1}{1 + dt(\alpha_f(A_d^{n+1} + A_r^{n+1}) - \beta_f)} f^n \quad (15) $$

## 3.1. Simulation results of the prostate PDE system

Androgen deprivation results in side effects related to the male characteristics of the patients that greatly influence their quality of life. Therefore, the evolution of the androgen level during the treatment should be monitored. This is also one of the reasons why intermittent treatment may be considered. Discontinuous androgen deprivation not only delays the appearance of castration-resistant cancer cells and reduces their number, but might also help to preserve the body's androgen reserves.

With brachytherapy, the concern is the preservation of the organs near the prostate that will also be exposed to the radiation depending on their proximity to the radioactive sources. Particularly, the urethra, which crosses the prostate at its center. Therefore, the dose absorbed by the urethra needs to be monitored [35].

The outputs of our model are thus the evolution of the different cell densities and of the androgen concentration as well as the distributions of the isodose curves received by the prostate and the urethra.

For the simulations, we considered a transverse section of the prostate. The prostate appears as an ellipse with a semi-major axis of 4cm and semi-minor axis of 3cm and is crossed at its center by the urethra which is shaped as a circle with a diameter of 10mm.

As initial conditions of our model, the tumor mass is considered localized on the left lobe of the prostate and the cancer cells represent 14% of the maximum cell carrying capacity $K_0$. The radioactive sources are placed inside the tumor mass. The objective is to explore different therapeutic strategies using the androgen deprivation therapy and the brachytherapy. The initial distributions of the prostate cells is shown in Figure 1.

We first consider the case of each therapy administered as a monotherapy before investigating the combination of the androgen deprivation therapy with the brachytherapy. The analysis of the results is supported, in each scenario, by clinical information.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Initial cells distribution in black and white printing (black indicates areas with no cell density) for: (A) Healthy prostate cells (B) Androgen dependent cancer cells (C) Androgen repressed cancer cells</figcaption></figure>

### 3.1.1. Androgen deprivation therapy simulations

At first, we deactivate the brachytherapy and deliver only the androgen deprivation therapy. To compare the effect of using continuous androgen deprivation and an intermittent treatment schedule, we observe the results of six months of treatment over a one-year period. For the continuous androgen suppression therapy (CAS), the treatment is applied over the first six months while for the intermittent androgen suppression therapy (IAS) the treatment is applied throughout the year depending on the PSA level for a maximum total duration of drug administration not exceeding six months. The numerical results are shown in Figure 2, Figure 3, Figure 4 and Figure 5.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Prostate cancer cell distribution in percentage of the maximum cell carrying capacity of the prostate at $t = 364$ day after CAS monotherapy.</figcaption></figure>

The prostate tumor is composed of androgen dependent and androgen resistant cancer cells. In Figure 4(A) where the CAS is administered over the first six month, the treatment is effective at first and the tumor volume comprising essentially of androgen-dependent cancer cells decreases quickly.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Prostate cancer cell distribution in percentage of the maximum cell carrying capacity of the prostate at $t = 364$ day after IAS monotherapy.</figcaption></figure>

However, the percentage of the androgen repressed cells starts to rise steadily and the treatment is no longer effective. From Figure 5 we notice that this happens when the environment is consistently depleted of androgen and the androgen concentration remains almost constant at a low level in the case of the CAS therapy. In Figure 5, at the end of the CAS therapy ($t = 182$ days) the androgen level starts to increase towards its saturation level causing the tumor volume composed mainly of the androgen-repressed cancer cells to start decreasing. However, the tumor had reached a percentage of almost 40% of the prostate volume which is a critical volume. In the case of the IAS therapy, the evolution of the androgen-repressed cancer cells, which do not respond to the androgen deprivation therapy, is contained and the tumor volume is bounded as shown in Figure 4 (B). The androgen level does not drop below 60% of its saturation level upon restarting the intermittent treatment, unlike in the continuous treatment where the androgen level is kept at 45% throughout the six months of treatment. This could represent a significant benefit for patients who struggle with the side effects of androgen deprivation.

From a clinical point of view, continuous hormone therapy necessarily leads to treatment failure because of the loss of sensitivity of cancer cells to androgen. Moreover, hormone therapy does not cure prostate cancer, but at best stops the proliferation of cancer cells and reduces tumor volume. This supports the results we obtained using a monotherapy based on androgen deprivation therapy. Furthermore, intermittent hormone therapy on its own is used for advanced late-stage cancer patients who need lifelong treatment. In this case, IAS helps in controlling the tumor volume, but most importantly, it offers the patient the opportunity to rest during the treatment breaks.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Cancer cells and tumor evolution under androgen deprivation therapy in percentage of the maximum cell carrying capacity of the prostate: (A) Continuous deprivation, (B) Intermittent deprivation.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Androgen deprivation therapy outputs: (A) Androgen concentration evolution in percentage of the androgen saturation level $A_0$, (B) ADT treatment schedule.</figcaption></figure>

### 3.1.2. Brachytherapy simulations

In order to observe the impact of brachytherapy on the different cell densities, we deactivate the hormone therapy to keep only the brachytherapy. Brachytherapy aims to kill cells through radiation emitted from a radioactive source placed inside the tumor mass. $I-125$ with its half-life of ~60 days retains only 50% of its radioactivity after 60 days, 12.5% after 180 days and < 2% after 360 days. The intensity of the treatment therefore gradually decreases until it disappears. Radiation does not distinguish between healthy and cancerous tissues. It is therefore necessary to monitor the dose absorbed by the organs at risk surrounding the prostate. The minimum dose to treat the tumor initially considered by brachytherapy alone for our model is 120Gy. We are interested in combining brachytherapy with another therapy to decrease the radiation dose while still curing the tumor. This way we can better protect the patients from the side effects of radiation. For a total dose of 120Gy, the isodose curves at the beginning of the treatment when the radioactivity is at its maximum are
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Isodose curves at $t = 0$ for a total dose of 120Gy.</figcaption></figure>

shown in Figure 6.

The difference between the total dose of 120Gy and the maximum absorbed dose of 106.11Gy in Figure 6 is due to the distance between the radioactive sources. The distribution of the absorbed doses in the prostate and urethra are shown in Figure 7 and the biological effect of these doses on the cancer cells is shown in Figure 8. The results of a treatment with a total dose of 120Gy and 80Gy are presented. The dose of 120Gy represents the minimum dose necessary to fully treat the initial considered tumor. The 80Gy dose represents a dose at which we can clearly observe a cancer relapse and therefore the dose to which we combine hormone therapy afterwards.

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Organ absorbed doses at $t = 0$ for a total dose of 120Gy.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Cancer cells evolution under brachytherapy in percentage of the maximum cell carrying capacity of the prostate.</figcaption></figure>

When comparing Figure 8 and Figure 4, brachytherapy is more effective than hormonal therapy since it allows to treat the tumor unlike hormonal therapy which, in the same time frame, allows at best to bound the tumor volume. And even in that case the bounded tumor volume is more important than that of the relapsed tumor with brachytherapy. However, brachytherapy does not spare healthy tissue as shown in Figure 9.

<figure><img src="image_10.png" /><figcaption>FIGURE 9. Evolution of healthy prostate cells in percentage of the maximum cell carrying capacity of the prostate in the case of different treatments.</figcaption></figure>

During brachytherapy treatment, the loss in healthy tissue is more important than in hormone therapy. We notice the positive growth of the healthy cells when the radioactivity of the sources has almost disappeared and the tumor has been treated. By reducing the administered dose from 120Gy to 80Gy, we can better preserve healthy cells. Nevertheless, the damage to healthy tissue remains significant compared to hormone therapy. From a clinical point of view, brachytherapy is recommended when the tumor is localized. In this case, it treats fully the tumor. It represents the equivalent of surgery since it kills the cells, unlike hormone therapy which stops their proliferation and aims to control their growth. In our case, given the considered initial conditions, the results obtained are in good accordance with the reality.
### 3.1.3. Combined treatment simulations

In order to find the best way to combine hormone therapy with brachytherapy we investigate two aspects. First, the timetable of the combination, either by initiating both treatments simultaneously or by delaying the start of one therapy compared to the other. Second, the impact of the intensity of each treatment on the effectiveness of the combination, by using continuous or intermittent for hormone therapy and by varying the dosage for brachytherapy.

#### Treatment schedule

The dose considered for brachytherapy is 80Gy for which we try to prevent the relapse or at least delay it by using a combined treatment. As hormone therapy has given better results in its intermittent version, we consider, at first, the combination of intermittent androgen deprivation therapy with a brachytherapy of 80Gy. In the first scenario, both therapies start simultaneously at ($t = 0$). As shown in Figure 10(A), the tumor relapse is delayed compared to the case where only 80Gy brachytherapy is used along with a much slower tumor regrowth. However, it is not enough to prevent the relapse of the disease. This is explained by the fact that the hormone therapy switched off at an early stage (after 13 days of treatment according to Figure 10(B)) and did not reactivate afterwards. In fact, at the beginning of the treatment, the curves for brachytherapy alone and in combination almost overlap except for a slightly more rapid decrease in tumor volume in the case of combined treatment, the brachytherapy dominates the hormone therapy and the resulting tumor volume does not allow the IAS to resume. To remedy this, a second scenario where the start of brachytherapy is delayed compared to that of hormone therapy was studied. In the case where only IAS was administered, the hormone therapy needed 26 days for the PSA level to fall below the PSA threshold to switch off the treatment (see Figure 5(B)). Therefore, for the second scenario, initially only the IAS is delivered at ($t=0$). Brachytherapy is then activated after 26 days at 80Gy. In this case the treatments are not simultaneous but occur one after the other: androgen deprivation therapy first then brachytherapy. Figure 10(A) shows that this combination leads to a slower decrease of the tumor volume at the beginning but better result in fine.

<figure><img src="image_11.png" /><figcaption>FIGURE 10. Different combining schedules: (A) Tumor evolution in percentage of the maximum cell carrying capacity under different combining protocols of androgen deprivation therapy and brachytherapy, (B) IAS schedule in the simultaneous combined treatment.</figcaption></figure>
## Treatment intensity

To study the impact of the intensity of each therapy on the combination of treatments, we first replace the intermittent hormone therapy with its continuous version and combine it with brachytherapy to assess the impact of hormone therapy on the treatment combination. As shown in Figure 11, at the beginning of treatment, the two curves for brachytherapy with intermittent deprivation therapy and brachytherapy with continuous androgen deprivation overlap. At the end of the observation period, we note a slight reduction of the tumor volume due to the use of continuous deprivation instead of the intermittent one. However, Figure 12 shows a significant loss in androgen due to CAS compared to the slight gain in tumor volume reduction. Using hormone therapy for nearly a month before stopping it and then switching to brachytherapy seems to be the best alternative.

<figure><img src="image_12.png" /><figcaption>FIGURE 11. Combined brachytherapy with both forms of hormone therapy: (A) Tumor evolution under brachytherapy at 80Gy combined with both forms of hormonal therapy in percentage of the maximum cell carrying capacity, (B) Androgen deprivation schedule in the case of a combined treatment of ADT therapy with brachytherapy at 80Gy.</figcaption></figure>

Secondly, we evaluate the impact of the dose used in brachytherapy on the effectiveness of the treatment. We first activate only the androgen deprivation therapy for 26 days before switching to brachytherapy at 80Gy for the first scenario and 70Gy for the second. As shown in Figure 13, the graphs of the two scenarios are identical for the first 26 days when only hormone therapy is activated. Thereafter, a slightly slower decrease in tumor volume is observed for the scenario with a dose of 70Gy. Furthermore, tumor relapse occurs earlier with a lower dose of brachytherapy. The radiation dose therefore plays an important role in the remission of the tumor.

Clinically, hormone therapy is combined with brachytherapy for advanced and metastatic cancers. Hormone therapy is then administered either at the same time than radiation therapy as initial treatment or before radiation therapy to reduce tumor volume and make the radiation treatment more effective. In our case, even if the tumor initially considered does not necessarily require a combined treatment, the latter allowed us to reduce the radiation dose as well as the androgen deprivation time and thus preserve the androgen level in the prostate.
<figure><img src="image_13.png" /><figcaption>FIGURE 12. Androgen concentration evolution in percentage of the androgen saturation level $A_0$ for a combined treatment of brachytherapy at 80Gy with both hormonal therapy.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 13. Tumor evolution in percentage of the maximum cell carrying capacity of the prostate under combined androgen deprivation therapy and brachytherapy.</figcaption></figure>

## 4. CONCLUSION

Prostate cancer is a heterogeneous cancer with a good prognosis. The treatment is multimodal, most often combining radiotherapy and hormone therapy. Clinically, brachytherapy finds its interest in the localized stages and its role in the advanced stages remains complementary. Whereas, hormone therapy is of interest in the advanced stages. Hence the interest in developing new therapeutic strategies with the help of mathematical modeling to explore the role of hormone therapy in the early stages but also to improve its combination with brachytherapy depending on the progress of the disease. The most relevant strategies are: neoadjuvant hormone therapy to brachytherapy and the combination of hormone therapy and brachytherapy.
In this paper, we developed a novel mathematical model for prostate cancer under hormone therapy and brachytherapy. The dynamics of the prostate cells helped us to write the partial differential equations describing the evolution of the different cell densities composing the prostate. Androgens are diffused in the prostate and act as a growth factor for androgen-dependent cancer cell. The first treatment implemented is the hormone therapy that affects the production of androgens in order to reduce the proliferation of the cells that depend on it. When we used continuous androgen suppression, the treatment became ineffective when the androgen-repressed cancer cells started proliferating and led to treatment failure. However, its intermittent form has given better results. While it did not treat the tumor, intermittent androgen suppression managed to control its volume. The second treatment implemented is internal radiotherapy for which we studied the effect of different doses on the tumor. The effect of radiation on the healthy tissue of the prostate and on the organ at risk, the urethra, was also assessed. Our results show that, unlike hormone therapy, brachytherapy was able to treat the tumor.

Combining hormone therapy with brachytherapy allowed us to reduce the dose used from 120Gy to 80Gy. When the treatments are given at the same time, the final tumor volume is significantly reduced compared to using each therapy separately. However, starting with hormone therapy gave better results. When using intermittent hormone therapy combined with brachytherapy, we found that once the PSA level drops below the critical level, it stays at reasonable levels and therefore the hormone therapy does not reactivate. When we use continuous hormone therapy instead, the prostate is unnecessarily deprived of androgen for an almost non-existent reduction in tumor volume compared to intermittent deprivation. The use of hormone therapy over a short period of time is therefore sufficient to give good results. The results also showed that the dose used in the combined treatments affects the tumor relapse.

The model we presented in this work is a first attempt to model and study the combination of hormone therapy with brachytherapy for prostate cancer. As discussed with clinicians, it would be interesting to investigate different initial tumor conditions. The initial conditions considered in this paper, represent the case where brachytherapy alone can treat the tumor and hormone therapy alone is not routinely used. We chose reasonably simple initial conditions to experiment with the model and get some initial results. The placement of the radioactive sources is therefore intuitive and simple since the tumor is located on one lobe of the prostate. The model could be extended to take into consideration the optimization of the seeds locations in order to minimize the dose absorbed by the organs at risk while minimizing the tumor volume. In the case presented here, we found that we do not need to optimize the treatment schedule of hormone therapy when used in combination. However, by changing the initial conditions and considering a more advanced stage of the tumor, hormone therapy could play a more important role and therefore also require an optimization of its parameters.

## REFERENCES

[1] Adhyam, M. and Gupta, A. K. (2012). A review on the clinical utility of psa in cancer prostate. *Indian journal of surgical oncology*, 3(2):120-129.

[2] Anderson, A. R. (2005). A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion. *Mathematical medicine and biology: a journal of the IMA*, 22(2):163-186.

[3] Attix, F. H. (2008). *Introduction to radiological physics and radiation dosimetry*. John Wiley & Sons.

[4] Berges, R. R., Vukanovic, J., Epstein, J. I., CarMichel, M., Cisek, L., Johnson, D. E., Veltri, R. W., Walsh, P. C., and Isaacs, J. T. (1995). Implication of cell kinetic changes during the progression of human prostatic cancer. *Clinical cancer research*, 1(5):473-480.

[5] Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., and Jemal, A. (2018). Global cancer statistics 2018: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries. *CA: a cancer journal for clinicians*, 68(6):394-424.

[6] Dale, R. G. (1985). The application of the linear-quadratic dose-effect equation to fractionated and pro-tracted radiotherapy. *The British journal of radiology*, 58(690):515-528.

[7] Deloule, S. (2014). *Développement d'une méthode de caractérisation spectrale des faisceaux de photons d'énergies inférieures à 150 keV utilisés en dosimétrie*. PhD thesis, Université Paris Sud-Paris XI.
[8] Eikenberry, S. E., Nagy, J. D., and Kuang, Y. (2010). The evolutionary impact of androgen levels on prostate cancer in a multi-scale mathematical model. *Biology direct*, 5(1):1–28.

[9] Friedman, A. (2012). Cancer as multifaceted disease. *Mathematical Modelling of Natural Phenomena*, 7(1):3–28.

[10] Friedman, A. and Jain, H. V. (2013). A partial differential equation model of metastasized prostatic cancer. *Mathematical biosciences and engineering: MBE*, 10(3):591–608.

[11] Fritz, M., Lima, E. A., Nikolić, V., Oden, J. T., and Wohlmuth, B. (2019). Local and nonlocal phase-field models of tumor growth and invasion due to ecm degradation. *Mathematical Models and Methods in Applied Sciences*, 29(13):2433–2468.

[12] Hecht, F. (2012). New development in freefem++. *Journal of numerical mathematics*, 20(3-4):251-266.

[13] Hecht, F., Pironneau, O., Le Hyaric, A., and Ohtsuka, K. (2005). Freefem++ manual.

[14] Holm, Å. (2013). *Mathematical optimization of HDR brachytherapy*. PhD thesis, Linköping University Electronic Press.

[15] Ideta, A. M., Tanaka, G., Takeuchi, T., and Aihara, K. (2008). A mathematical model of intermittent androgen suppression for prostate cancer. *Journal of nonlinear science*, 18(6):593–614.

[16] Jackson, T. L. (2004a). A mathematical investigation of the multiple pathways to recurrent prostate cancer: comparison with experimental data. *Neoplasia*, 6(6):697–704.

[17] Jackson, T. L. (2004b). A mathematical model of prostate tumor growth and androgen-independent relapse. *Discrete & Continuous Dynamical Systems-B*, 4(1):187.

[18] Jain, H. V., Clinton, S. K., Bhinder, A., and Friedman, A. (2011). Mathematical modeling of prostate cancer progression in response to androgen ablation therapy. *Proceedings of the National Academy of Sciences*, 108(49):19701–19706.

[19] Jain, H. V. and Friedman, A. (2013). Modeling prostate cancer response to continuous versus intermittent androgen ablation therapy. *Discrete & Continuous Dynamical Systems-B*, 18(4):945.

[20] Jones, B., Dale, R., Bleasdale, C., Tan, L., and Davies, M. (1994). A mathematical model of intraluminal and intracavitary brachytherapy. *The British journal of radiology*, 67(800):805–812.

[21] Kanikowski, M., Skowronek, J., Kubaszewska, M., CHICHEŁ, A., and Milecki, P. (2008). Permanent implants in treatment of prostate cancer. *Reports of Practical Oncology & Radiotherapy*, 13(3):150–167.

[22] King, M. T., Keyes, M., Frank, S. J., Crook, J. M., Butler, W. M., and Rossi, P. J. e. a. (2021). Low dose rate brachytherapy for primary treatment of localized prostate cancer: A systemic review and executive summary of an evidence-based consensus statement. *Brachytherapy*, 20(6):1114–1129.

[23] Mi, H., Petitjean, C., Dubray, B., Vera, P., and Ruan, S. (2014). Prediction of lung tumor evolution during radiotherapy in individual patients with pet. *IEEE transactions on medical imaging*, 33(4):995–1003.

[24] Miyamoto, H., Messing, E. M., and Chang, C. (2004). Androgen deprivation therapy for prostate cancer: current status and future prospects. *The Prostate*, 61(4):332–353.

[25] Nguyen Edalgo, Y. T. and Ford Versypt, A. N. (2018). Mathematical modeling of metastatic cancer migration through a remodeling extracellular matrix. *Processes*, 6(5):58.

[26] Pavel, M., Renna, M., Park, S. J., Menzies, F. M., Ricketts, T., Füllgrabe, J., Ashkenazi, A., Frake, R. A., Lombarte, A. C., Bento, C. F., et al. (2018). Contact inhibition controls cell survival and proliferation via yap/taz-autophagy axis. *Nature communications*, 9(1):1–18.

[27] Podgorsak, E. B. (2005). Basic radiation physics. *Radiation oncology physics: a handbook for teachers and students*. Vienna: IAEA, 7.

[28] Portz, T., Kuang, Y., and Nagy, J. D. (2012). A clinical data validated mathematical model of prostate cancer growth under intermittent androgen suppression therapy. *Aip Advances*, 2(1):011002.

[29] Rikitake, Y. and Takai, Y. (2011). Chapter three - directional cell migration: Regulation by small g proteins, nectin-like molecule-5, and afadin. *International Review of Cell and Molecular Biology*, 287:97–143.

[30] Sharifi, N., Gulley, J. L., and Dahut, W. L. (2005). Androgen deprivation therapy for prostate cancer. *Jama*, 294(2):238–244.
[31] Sieminska, L., Ferguson, M., Zerda, T. W., and Couch, E. (1997). Diffusion of steroids in porous sol-gel glass: application in slow drug delivery. *Journal of Sol-Gel Science and Technology*, 8(1):1105-1109.

[32] Swanson, K. R., True, L. D., Lin, D. W., Buhler, K. R., Vessella, R., and Murray, J. D. (2001). A quantitative model for the dynamics of serum prostate-specific antigen as a marker for cancerous growth: an explanation for a medical anomaly. *The American journal of pathology*, 158(6):2195-2199.

[33] Tanaka, G., Hirata, Y., Goldenberg, S. L., Bruchovsky, N., and Aihara, K. (2010). Mathematical modelling of prostate cancer growth and its application to hormone therapy. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 368(1930):5029-5044.

[34] Toma, A., Mang, A., Schuetz, T. A., Becker, S., and Buzug, T. M. (2012). A novel method for simulating the extracellular matrix in models of tumour growth. *Computational and mathematical methods in medicine*, 2012.

[35] Verry, C. (2011). Traitement du cancer localisé de la prostate par curiethérapie à l'iode 125: à propos de 200 patients traités de juillet 2001 à janvier 2011 au chu de grenoble. *Thèse de doctorat*.

[36] Wang, J. Z., Guerrero, M., and Li, X. A. (2003). How low is the $α/β$ ratio for prostate cancer? *International Journal of Radiation Oncology*\* *Biology*\* *Physics*, 55(1):194-203.

[37] Wong, H. and Tang, W. (2011). Effects of ecm degradation rate, adhesion, and drag on cell migration in 3d. In *5th Kuala Lumpur International Conference on Biomedical Engineering 2011*, pages 428-431. Springer.