PCa dynamics with neuroendocrine differentiation and distributed delay

1.   Introduction

Prostate cancer (PCa) is the second most common cause of cancer among men worldwide ^[1,2]. Much work has been done to understand the development of this disease. Over the last few decades many biological models, such as TRAMP and LADY ^[3] have been created using various strains of genetically engineered mice to simulate PCa growth observed in humans. Due to the limitations of these mouse models, in vitro experiments have been designed using cells taken from specific strains of human prostate cancer, such as the LNCaP cell line that was established in 1980 by Horoszewicz et al. ^[4]. In the last few years, many mathematical models have been proposed to describe the rich dynamics of prostate cancer with or without treatment, and one the major challenges has been the development of dynamical systems that would represent in vivo conditions, while being mathematically tractable ^[5]. Early mathematical models investigated the link between the concentration of serum prostate specific antigen (PSA) and tumour volume ^[6,7,8]. Others explored how different chemical resources would affect tumour growth, e.g. Kuang et al. proposed the KNE model ^[9], in which they apply ideas from ecological stoichiometry to tumour growth, and consider that tumour growth is limited by various physical constraints such as space, nutrients, and vasculature. However, as the link between androgen, the male sex hormones testosterone and dihydrotestosterone, and prostate cancer was explored further ^[10,11,12], mathematical models began to focus on the response of PCa tumour growth to androgen concentration. In 2004 Jackson modelled the tumour as two populations, one androgen dependent (AD) and the other androgen independent (AI) ^[13,14]. Using these models Jackson predicted that androgen deprivation therapy (ADT) would successfully control tumour growth in a well-defined region of the parameter space for all time, but that cancer could recur through AI mechanisms after undergoing ADT. These concepts were further explored by Ideta et al. ^[15], who showed that using intermittent ADT can reduce the time to cancer relapse. They stated, counter intuitively, that the reduced time before cancer reappears is an acceptable trade off when considering the known side effects of ADT and other possible improvements in the quality of the patient's life. Later, in 2010, Eikenberry et al. ^[16] proposed two mathematical models that consider the intracellular kinetics of the androgen and its receptors. The authors found that decreasing androgen levels could increase the PCa cell mutation rate, resulting in a more heterogeneous population.

Recent work has considered the role of neuroendocrine cells in the re-emergence of PCa tumours. Neuroendocrine cells are specialised secretion cells, with a cell structure similar to neurons. They are found throughout the human body, including in glands such as the prostate, and usually contribute to the homeostasis of the surrounding tissues ^[17] by secreting various hormones and proteins ^[18]. Whilst there have been multiple observations of neuroendocrine cells being present in prostate tumours, the current theories for their role in cancer development are still considered controversial. One of these theories concerns the role of neuroendocrine transdifferentiation, which is believed to be caused by the reduction in androgen levels ^[17]. This theory proposes that after the tumour has been under castrate conditions for 16 to 18 months, such as those caused by ADT, a proportion of the PCa cells undergo transdifferentiation and become neuroendocrine cells. Transdifferentiation is the irreversible switch of one type of cell to another ^[19]. Once this switch occurs, neuroendocrine cells are believed to secrete androgen or similar anabolic hormones which promote tumour growth.

In 2015 Cerasuolo et al. ^[20] proposed a discrete delay dynamical system to investigate neuroendocrine transdifferentiation in PCa, based on in vitro experiments of androgen-deprived conditions on LNCaP cells growing in Petri dishes. In these experiments the LNCaP cells were first grown in an androgen-rich environment, before being transferred to the androgen-deprived condition of the Petri dishes. The mathematical model was inspired by previous work on cell differentiation in hematopoiesis ^[21,22] and on PCa, where the cancer cell population is divided in androgen-dependent and androgen-independent cells that are able to proliferate under in vivo conditions ^[16,23]. In ^[20] the model represents two cancer cell populations, one with androgen-dependent cells and the other with neuroendocrine androgen-independent cells. The model also considers the androgen concentration.

The model was parameterised against experimental data and was used to forecast tumour growth in the long term. In the first instance, the authors showed agreement between the in vitro experiments and the simulated growth curve. They then simulated the long-term behaviour over 400 days, and observed that at the beginning androgen-dependent cells would become almost extinguished while the neuroendocrine cells remained nearly constant for the first 150 days. This was followed by an increase in both cell populations, with the system reaching an equilibrium after another 150 days. This predictive analysis led to the assertion that androgen-dependent cells react to hormone deprivation by favouring the establishment of a neuroendocrine cell population, which leads to the development of androgen-resistant PCa in most patients. The authors concluded that the main novelty of the model stemmed from considering cell transdifferentiation as a consequence of androgen concentration, and that the most interesting result was the active role that differentiated cells can play in sustaining the tumour in castrate conditions. A mathematical and numerical analysis of the system proposed in ^[20] was performed by Turner et al. ^[24]. Here, the authors found sufficient conditions that would ensure the existence of tumour-present equilibria, but it was not possible to locate or explore the equilibria stability analytically. In ^[24], several amendments were suggested that would aid the analysis of the model, as well as move the model from in vitro to in vivo conditions.

The first amendment concerned the proliferation of PCa cells. A condition of in vivo models is that inter-cellular competition for space and resources have a large effect on the growth rate of the population. This led us to consider a logistic growth instead of the cell-doubling mechanism used in ^[20]. The second major change concerned the two delays in the model. In ^[24] the authors showed that both delays are not highly significant parameters for the model and are not involved in the local stability of the trivial steady state. In the revised new model, we therefore no longer consider a delay in the transdifferentiation from LNCaP cells to NE cells, nor in the mitosis rate. However, like Hutchinson ^[25], we assume that the per capita growth rate is negatively affected by the population density after a time $\tau$ representing the development and maturation of the species. Following Cassidy and Humphries ^[26], who consider a heterogeneous cell-cycle duration more realistic than a discrete delay, we consider $\tau$ to be a distributed delay, as discrete delays assume uniform and constant behaviour across all cells, while a distributed delay will allow for variation in a small interval. In this way we take into account the fact that the dependence of the growth limitation on population density is distributed over a past time interval and not concentrated in a single instant. With this distributed delay, our logistic term will be similar to that in Volterra's proposed predator-prey system introduced in 1931, which was translated by Scudo and Ziegler's ^[27].

The paper is organised as follows. Section 2 is devoted to the model formulation and the proof of its basic properties. The model is a distributed-delay system that considers the dynamics of two cell populations, the androgen dependent LNCaP cells, denoted by $L(t)$, and the neuroendocrine androgen independent cells, denoted by $N(t)$, and the dynamics of the androgen concentration in the environment, denoted by $A(t)$. Sufficient conditions that ensure non-negativity and boundedness of the solutions are found. In Section 3 conditions for the existence of a tumour-free and a tumour-present equilibrium are given. The local stability analysis of the tumour-free equilibrium is performed. Furthermore, sufficient conditions for the existence of bifurcations for the tumour-present equilibrium are given in the case of a uniformly distributed delay. In Section 4 we illustrate the systems behaviour with numerical simulations in the case of uniformly and exponentially distributed delays, and finally a discussion of our results and further research ideas are in given in Section 5.

2.   The model and its basic properties

For brevity in the following, we will denote $x(t)$ as $x$ for $x \in \{A \, ,L \, ,N\}$. The equation for the androgen dynamics follows the modelling by Baez and Kuang ^[28], where $A_{min}$ represents the minimum androgen concentration required to sustain the tumour, and $A_{max} > A_{min}$, the maximum androgen concentration of the system. Here we consider that androgen can be produced from other endocrine glands such as the adrenal glands and kidneys, with production rate $\gamma$ and diffusion rate proportional to the difference between $A_{max}$ and $A$. Androgen depletion depends on its maximum rate, $\mu_A$, and assumes the existence of a minimum threshold of androgen required to sustain the cancer. Androgen is secreted by $N$ cells with secretion rate, $\kappa$, which leads to the equation

In the model we assume that $L$ cells are generated through asymmetric cell division, undergo apoptosis at rate $\delta_L$, and transdifferentiate into $N$ cells with a maximum rate $k_p$ during growth, and a maximum rate $k_t$ when fully mature. The process of transdifferentiation depends on the androgen concentration through the Ricker function $\alpha(A)$, which is

The proliferation of $L$ cells is governed by $A$ through the Droop equation ^[29]

where $\beta_P$ is the maximum proliferation rate of $L$ and $A_{min}$ is the minimum androgen concentration required for tumour growth, as in (2.1). Also, we assume that the growth of $L$ cells is limited by the population density, and that such limitation is distributed over a past time interval $\tau$ ^[26]. Therefore, the equation for $L$ is

where $\eta_k$ is the carrying capacity of the environment, $\omega(\theta)$, the delay kernel or delay distribution function and finally $L_{\theta} = L(t+\theta)$ is the delay term for LNCaP cells.

The growth of $N$-cells depends on the production through asymmetrical cell division of $L$ cells and on the transdifferentiation of $L$ cells; and based on experimental evidence, the death rate contains a linear term for the apoptosis, with maximum rate $\mu_N$, and a quadratic term, $-\delta_N N^2$, representing the intraspecific competition for space and resources ^[20,28]:

The final model is

For biological relevance we consider that all parameters are positive, and that the initial conditions are non-negative continuous functions defined as

where $\phi = (\phi_1, \phi_2, \phi_3)^T \in C$ is such that $\phi_i (\theta) \geq 0, \, (-\tau\leq \theta \leq 0, \, i = 1, 2, 3)$. $C$ denotes the Banach space $C([-\tau, 0], \mathbb{R}^3_+)$ of continuous functions mapping the interval $[-\tau, 0]$ into $\mathbb{R}^3_+$ with the supremum norm

where $|\cdot|$ is any norm in $\mathbb{R}^3_+$ ^[30].

Also, by observing that the proportion of $L$ cells produced by asymmetric division must satisfy

and that $\alpha_{max} = \frac{r e^{-1}}{a}$, then we must impose the condition

To prove the boundedness of the solutions to (2.2) we use Theorem 2.1 from Faria and Liz ^[31], which is

Theorem 2.1. Consider

where the following hold.

(H1) We have

(i) $b: \mathrm{R} \to \mathrm{R}$ is a continuous function, and there are constants $\beta_0$, $\beta^0$ such that $0 < \beta_0\leq b(t)\leq \beta^0$ for all $t \in \mathrm{R}$;

(ii) $M:C \to \mathrm{R}$ is a positive linear operator with $l: = ||M|| > 0$.

Define $x_* = l^{-1}$ and consider the IVP for (2.4) with $x_0 = \varphi$, with solution $x(\varphi)(t)$. Then the following results hold.

(i) $x(\varphi)(t)$ is defined for $t\geq 0$, bounded below from zero and bounded on $[0, \infty)$.

(ii) If $x(\varphi)(t)$ is non-oscillatory around $x_*$, then $x(\varphi)(t)\to x_*$ as $t \to \infty$.

(iii) Moreover, if $x(\varphi)(t)$ is oscillatory around $x_*$ and

then $x(\varphi)(t) \to x_*$ as $t \to \infty$.

Lemma 2.2. Let $\phi_i(s) \geq 0$ and $\phi_1(s) \geq A_{min}$, with $s \in [-\tau, 0]$ and $i = 2, 3$. If either of the following conditions hold,

where $L_B = \sup_{[0, +\infty)}L(t)$ and $\Delta = \int^0_{-\tau}\omega(\theta)\mathrm{d}\theta$, then the solutions of (2.2) are non-negative and bounded. In particular, $A(t) \geq A_{min}\, \, \forall t\geq 0$.

Proof. By considering the second equation of (2.2)

with $L(0) = \phi_2(0) > 0$, then

If $L(0) = 0$ then $L(t) = 0 \; \forall\, t\geq 0$.

Let $t_0 > 0$ be such that $A(t_0) = A_{min}$, we wish to prove that $A(t) \geq A_{min} \, \forall t > 0$. Let us assume, on the contrary, that $A(t) < A_{min}$ for all $t \in (t_0, t_0+\epsilon)$ with $0 < \epsilon < \tau$ and for $t < t_0, \, A(t) > A_{min}$. Then the first equation of (2.2) at $t_0$ becomes

We will determine the sign of $N(t_0)$. From the third equation of (2.2) at $t_0$

where $f_N(N(t_0)) : = (\delta_NN(t_0)+\mu_N)$ and $B(A(t_0), L(t_0)) = k_t\alpha(A(t_0))L(t_0)$. Then by the variation of constant formula, we obtain

and as all initial conditions are non-negative, $A(t_0) > 0$ and $L(t)\geq 0 \, \forall\, t\geq 0$, then $N(t_0) \geq 0$, and therefore $\frac{dA(t_0)}{dt}\geq 0$. This is a contradiction of the initial hypothesis, and therefore $A(t) \geq A_{min}\, \, \forall\, t\geq 0$.

Before proving the positivity of $N$, we need to prove the boundedness of $L$. Following Faria and Liz ^[31], and as $L(\theta) = \phi_1(\theta)\geq 0$ with $\theta \in [-\tau, 0]$ and $L(t) \geq 0, \, A(t)\geq A_{min} \, \forall t \geq 0$, from the second equation of (2.2) we obtain

that is,

where $M: C\to \mathrm{R}$ is a non-zero positive integral linear operator, i.e. $M(\phi_2)\geq 0 \, \forall\, \phi_2\in C, \, \text{s.t.}\, \phi_2\geq 0.$

We can observe that all hypotheses of Theorem 2.1 are satisfied for the equation

This result, together with Gronwall's Inequality, ensures that $L$ is bounded on $[0, +\infty)$.

We can now prove the positivity of $N$. Let us consider the third equation of (2.2)

Let there exist a $t_0$ such that $N(t_0) = 0$, $N(t) > 0$ for $t < t_0$, and $N(t) < 0$ for $t\in(t_0, t_0+\epsilon]$, then

so if either hypothesis (2.5) or (2.6) holds then $\frac{d N(t_0)}{dt}$ is greater than zero, which is a contradiction of $N(t) < 0$ for $t\in (t_0, t_0+\epsilon]$. Therefore, $N(t) \geq 0$ for all $t \geq 0$.

Having shown non negativity of $N$ we now prove its boundedness. Let us consider the third equation of (2.2)

As $L$ is bounded, $\alpha(A) \leq \alpha_{max} \, \forall A$ and from (2.3), our equation becomes

where $L_B: = \sup_{t \in [0, +\infty)}L(t)$. This can be solved to give

and therefore $N(t)$ is bounded above for all $t \geq 0$.

Finally, from the first equation of (2.2) and using the variation of constant formula and Gronwall's Inequality we get

where $N_B: = \sup_{[0, +\infty)}N(t)$. Therefore $A(t)$ is bounded for all $t \geq 0$, which completes the proof of non-negativity and boundedness of the solutions to (2.2).

3.   Equilibria and their stability properties

The equilibria of (2.2) solve the system

and

where $\Delta = \int^0_{-\tau}\omega(\theta)\mathrm{d}\theta$, and $X^*$ is an equilibrium point with $X = (A, L, N)$. From (3.1) we get that

and by considering that (3.2) has a common factor $L^*$ it becomes apparent that there is a tumour-free equilibrium $E_1$ at

It is also possible to prove the existence of a tumour-present equilibrium (all non-zero coordinates). From (3.2) we get

This can be solved for $L^*$ to give

Since, for biological reasons, $L^* > 0$, we obtain

Summing (3.2) and (3.3) we get that

and substituting $F(A)\left(1-\frac{L^*\Delta}{\eta_k}\right) = \frac{\delta_L+k_t\alpha(A^*)}{1-k_p\alpha(A^*)}$ into (3.8) gives

Therefore, for the equilibrium to exist (3.7) and (3.9) must be satisfied. Observing that $L^*$ and $N^*$ can be expressed as functions of $A^*$, by (3.4) and (3.6) we can rewrite (3.9) as a function of $A^*$. If we consider a tumour-present equilibrium, then the right-hand side becomes

and the left-hand side becomes

Note that $G(A^*)$ is a quadratic function whose minimum can be obtained by solving

which yields

Therefore, if we define

then

which means that $\min_{A > A_1} G(A) < 0$. As $A_2 < A_1$ and $G(A_1) = 0$ we have that $G(A^*) > 0$ for all $A^* > A_1$, which is also the condition that ensures the existence of an equilibrium with positive $N^*$. If we now focus our attention on $H(A^*)$, defined in (3.10), from its expression we can see that

and that

Therefore, if we can show that $H(A_1) \geq 0$ then there exist $\bar A > A_1$ such that $H(\bar A) = G(\bar A)$, which implies the existence of a tumour-present equilibrium. Since,

for $H(A_1)$ to be non-negative, we require that either the factor inside both brackets are greater than zero, both are less than zero or that one or both are equal to zero. As the second bracket is positive due to (2.3) then we only need to consider the first bracket being greater or equal to zero. Therefore, if we define

we require that

The following lemma can be proved easily,

Lemma 3.1. If $A^* > A_1 \geq \dfrac{1}{a}$ and $R_0 \leq 1$ then condition (3.7) holds.

Proof. If $A^* > A_1$ and $A_1 \geq \frac{1}{a}$ then we obtain that $\alpha(A_1) > \alpha(A^*)$ and since the function $h(A) = \frac{A}{A-A_{min}}$ is strictly decreasing for all $A \in \mathbb R^+$ we obtain that

Using Lemma 3.1 we can state the following:

Theorem 3.2. System $\mathrm{(2.2)}$ always admits the tumour-free equilibrium

and if the conditions

and

hold, then there exists a tumour-present equilibrium $E^* = (A^*, L^*, N^*)$, whose coordinates are given by

and where $A^* \, (>A_1)$ solves the equation

3.1.   Local stability analysis

Given the generic equilibrium $\bar{E} = (\bar{A}, \bar{L}, \bar{N})$ and by using the following change of variables $x_1 = A-\bar{A}, \, x_2 = L-\bar{L}$ and $x_3 = N-\bar{N}$ we can rewrite system (2.2) as follows

The stability analysis of $\bar{E}$ can be performed by analysing the characteristic equation associated with the linearised system of (3.13), which is calculated by using a truncated Taylor expansion about $(0, 0, 0)$, that is:

where

and

Remark 3.1. We can represent system $\mathrm{(3.14)}$ as a sum of two terms as follows:

where

and

The characteristic equation is then given by

which can be written as:

where

At $E_1 = (A_1, \, 0, \, 0)$ equation (3.15) yields

which admits the two negative roots $\lambda_1 = -(\gamma+\mu_A)$ and $\lambda_2 = -\mu_N,$ and

Therefore $E_1$ is locally asymptotically stable if

that is,

This result is summarised in the next lemma.

Lemma 3.3. The tumour-free equilibrium $E_1$ is locally asymptotically stable if

and is unstable otherwise.

Remark 3.2. We can observe that the condition for the instability of $E_1$ ensures the existence of a tumour-present equilibrium (see Theorem 3.2).

If we consider $E^* = (A^*, L^*, N^*)$, the tumour-present equilibrium, then from (3.5) we obtain $b_2 = 0$, and from (3.3) we get that

and therefore

Using this we can rewrite $b_4 = N^*(\delta_N N^*+ \mu_N) > 0$. Then the associated characteristic equation is

For $\tau = 0$, this becomes

From the conditions (3.11) and (3.12) that ensure the existence of a tumour-present equilibrium, and observing that $A^* > 1/a$, we obtain $\alpha'(A^*) < 0$ and $L^* < \frac{\eta_k}{\Delta}$ which imply $b_1 > 0$. Hence, by Descarte's rule we have that at least one solution to (3.17) has a positive real part, and therefore the tumour-present equilibrium is unstable when $\tau = 0$, and we thus have the following theorem.

Theorem 3.4. If $\tau = 0$ then whenever $E^*$ exists, it is unstable.

To understand how the stability properties of $E^*$ change with increasing $\tau$ we need to specify the distribution function $\omega(\theta)$. We will consider two delay distributions, the uniform and the exponential distributions, ^[32,33]. The uniform distribution uses the kernel

and therefore gives equal weight to all of the history incorporated by the delay. From a biological point of view, this type of distribution assumes that the development and maturation time are bounded above by $\tau$ days, and that an individual cell will mature in $t$ days, with $t \in[0, \tau]$. The uniform distribution is typically used if there is no a priori information on the behaviour of the delay, as is the case considered in this paper. On the other hand, the exponential distribution is the one mostly used in the literature and is represented by

In this case the history of the $L$ cells is assumed to have a much greater weight closer to the present time. Firstly, using the uniform distribution (3.18), we have that when $\tau \neq 0$

Then for $\lambda \neq 0$ and $\tau \neq 0$ the characteristic equation (3.16) becomes

We can observe that the characteristic equation is of the type:

where $P_n(\lambda, \tau) = \sum_{k = 0}^{n} p_k(\tau) \lambda^k$ and $Q_m(\lambda, \tau) = \sum_{k = 0}^{m} q_k(\tau) \lambda^k$, with $n, m\in \mathbb{N}, \, n > m$ and $p_k(\tau)$ and $q_k(\tau)$ are continuous and differentiable functions of $\tau$. In this case $n = 4$ and $m = 2$. For such characteristic equations, where the coefficients depend upon the delay, we can follow the geometric stability switch criteria introduced by Beretta and Kuang ^[34] to define functions whose zeroes are the critical values at which Hopf bifurcations occur.

Let us consider $\lambda = i\omega$. We separate $P_4(i \omega, \tau)$ and $Q_2(i \omega, \tau)$ into real part

and imaginary part

where

We define the function $F(\omega, \tau)$ as:

which is such that its zeroes are the Hopf frequencies associated with the imaginary roots of the characteristic equation (3.20). Factorising $\omega^2$ out we obtain the equation:

Note that $\omega = 0$ cannot be considered as a solution to $F(\omega, \tau)$ for $\lambda \neq 0$.

Lemma 3.5. If $A^* > \frac{1}{a}+A_{min}$ and $\tau < \tau^* \left( = \frac{2 (b_1 b_4 c_1 \kappa z_2 + z_1 z_3)}{b_1^2 b_4^2 \kappa^2}\right)$ in (3.21), then $F(\omega, \tau)$ has at least one positive real root $\omega^*$.

Proof. Taking $A^* > \frac{1}{a}+A_{min}$ implies that $b_3 < 0$ and therefore $z_1 > 0$, which ensures that $\tau^* > 0$. Since $F(\omega, \tau)$ is a continuous function of $\omega$ over the interval $[0, +\infty)$, ${\lim_{\omega \to \infty}} F(\omega, \tau) = +\infty,$ and $F(0, \tau) = -2 z_1 z_3 - 2 b_1 b_4 c_1 z_2 \kappa + b_1^2 b_4^2 \kappa^2 \tau < 0$ when $\tau < \tau^*$, then there exists at least one positive solution to the equation $F(\omega, \tau) = 0$.

Lemma 3.5 states that if $\tau \in (0, \tau^*) ( = I)$ then the tumour-present equilibrium can experience a stability switch in which system (2.2) undergoes a Hopf bifurcation. However, that does not give us information about the critical value $\tau_c$ at which the bifurcation occurs. For this, let us substitute $\lambda = i \omega$ in $D(\lambda, \tau) = 0$. In this way we obtain:

By considering real and imaginary parts we obtain:

which gives

Dividing $\cos(\tau\omega)$ by $\sin(\tau\omega)$ we obtain that

where

with

from which we obtain the map

where $\omega(\tau)$ is a positive real solution to (3.21). We now define the functions

which are continuous and differentiable in $\tau$ as proved in ^[34], and whose solutions represent the values of the time delay at which the Hopf bifurcation occurs.

Let us now assume that $\omega(\tau)$ is a positive real root of (3.21) defined for $\tau \in I$ and that $\tau_c$ is a root of (3.22) for some $n \in \mathbb{N}$. Then Theorem 2.2. of ^[34] ensures that there exists a pair of roots that crosses the imaginary axis from left to right at the pure imaginary values $\lambda_\pm(\tau_c) = \pm i \omega(\tau_c)$, with

and from right to left otherwise, where $\lambda(\tau)$ is a solution to the characteristic equation (3.20) for $\tau$ in a neighbourhood of $\tau_c$.

Therefore, from Lemma 3.5 and applying Theorem 2.2. of ^[34] we can obtain the following theorem about the stability of the tumour-present equilibrium of system (2.2) and the occurrence of a Hopf bifurcation.

Theorem 3.6. Let $\omega(\theta)$ be defined as in (3.18). If (i) $A^* > \frac{1}{a}+A_{min}$, (ii) (3.21) has at least one positive real root $\omega(\tau)$ defined for $\tau \in I$, and (iii) at some $\tau_c \in I$, $S_k(\tau_c) = 0$ for some $k \in \mathbb{N}$ and $F_\omega'(\omega(\tau_c), \tau_c) \neq 0$, then the tumour-present equilibrium undergoes a Hopf bifurcation at a critical value of the time delay $\tau = \tau_c$ if $\rm{sign}\{ F_\omega'(\omega(\tau_c), \tau_c) \} \cdot \rm{sign} \left\{ \left.\frac{d S_n(\tau)}{d\tau}\right|_{\tau = \tau_c} \right\} > 0$.

The results stated above, together with Theorem 3.6, will be used to calculate the values $\tau_c$ where the stability switches occur once the parameter values are chosen.

On the ground of analytical tractability, the exponential distribution will only be used to run numerical simulations and to emphasise how strongly the choice of the delay kernel can influence the dynamics of the system.

4.   Numerical simulations

In order to perform numerical simulations of (2.2) we transform the distributed-delay into an equivalent discrete-delay system using the linear chain trick ^[35]. This transformation retains the same equilibria, properties and behaviour of the original system, and it is specific to the distribution function used. Details of the transformation in the case of uniform and exponential distributions can be found in the Appendix.

4.1.   Parameter values and initial conditions

Using (A.1), (A.3) and the dde23 solver in Matlab R2020a we were able to produce numerical simulations that highlight the difference between the dynamics produced by the two distributions. For the following figures we define the history $\phi(\theta) = (\phi_1, \phi_2, \phi_3, \phi_4)^T$, such that $\phi_i(\theta) \geq 0$ for $-\tau \leq \theta < 0$ and $i = 1, 2, 3, 4.$ To mimic in vivo conditions we choose the history to reflect normal physiological conditions for the prostate, and the initial conditions to represent the start of the ADT. Therefore, we have that $\phi(\theta) = [10, \, 3.1, \, 0, \, U_0]^T$ with initial conditions $[1, \, 3.1, \, 0, \, U_0]^T$ ^[20], where for the uniform distribution $U_0 = 3.1$, and for the exponential distribution $U_0 = \frac{3.1\zeta}{\psi}(1-e^{-\psi\tau})$. It should be noted that all simulations are run for a time span of 200 days, and for clarity we plot only the dynamics of $A, \, L$ and $N$.

Most of the parameter values are taken from ^[20] (see Table 1). However, the value for $\kappa$ has been decreased to $0.009$ (from $0.017$ in ^[20]) to account for the fact that we are now considering external sources of androgen secretion and a higher value for this parameter would have resulted in an over-production of androgen; and $\mu_N$ ($\mu$ in ^[20]) has been increased from $0.006$ to $0.08$, as we now consider an in vivo situation in which cells can become apoptotic more easily than in in vitro conditions, such as the one described in ^[20] where all resources for cell growth except for androgen were provided in abundance.

Remark 4.1. As shown in Figure 1, by using the parameter values in Table 1 we obtain that the function $S_0(\tau)$ has two roots. The first root, ($\tau_c^1 \approx 0$) is such that sign$\left\{ \left.\frac{dRe (\lambda)}{d\tau}\right|_{\lambda = i \omega(\tau_c^1)} \right\} < 0$, therefore there exists a pair of pure imaginary roots that crosses the imaginary axis from the right to the left half-plane, which implies that the tumour-present equilibrium becomes stable. The second root ($\tau_c^2 \approx 7.6$) is such that sign$\left\{ \left.\frac{dRe (\lambda)}{d\tau}\right|_{\lambda = i \omega(\tau_c^2)} \right\} > 0$, which means that a pair of pure imaginary roots crosses the imaginary axis from left to right, and therefore a Hopf bifurcation occurs (Theorem 3.6). The same behaviour is observed for all $S_k(\tau), \, k > 0$, suggesting that the only stability switch is the one observed with $S_0(\tau)$.

Figure 2a shows the system dynamics with uniform and exponential distributed delay obtained using the parameters in Table 1. With this choice of parameter values, both transformed systems (A.1) and (A.3) exhibit the same dynamics with the trajectories tending towards the tumour-present equilibrium $E_2$, which is to be expected as in this case $A_1 \geq \frac{1}{a}$ and condition (3.7) is satisfied, which imply the existence of a (stable) tumour-present equilibrium (Theorem 3.2).

By choosing $a = 1$ neither conditions (3.11) nor (3.12) for the existence of the tumour-present equilibrium is satisfied. This case is represented in Figure (2b), where for both distributions the system solutions tend towards the tumour-free equilibrium $E_1$. As proved in Theorem 3.6, the stability of the tumour-present equilibrium is affected by $\tau$. In Figure 3 different $\tau$ values are considered, $\tau = 0.01, 10, 20$, to demonstrate the effect on the system dynamics.

Decreasing $\tau$ (Figures 3a) does not affect the solutions of the system (A.1) with uniform distribution. Increasing $\tau$ to $10$ days introduces oscillatory behaviour to system (A.1), as seen in Figure 3b. Increasing further to $\tau = 20$, we observe an increase in the amplitude and wavelength of the oscillations, Figure 3c. We can also note that the $L$ population appears to get closer to zero cyclically, with a value of approximately $1.1\times 10^{-4}$. In the case of the exponential distribution, the dynamics of the solutions does not change by altering $\tau$ without also modifying the other parameters.

Since, in the case of exponential distributed delay, it was not possible to prove analytical results on the stability properties of the tumour-present equilibrium, we explored the dynamics of the system numerically, by considering the parameter ranges in Table 2. Simulations were run selecting parameter values at random for each simulation.

We should observe that the proposed parameter ranges far exceed the ranges found in ^[20], and that the parameter combinations were chosen so that they would satisfy condition (2.3), which ensures the biological relevance of the solutions. Using this method we were able to find oscillatory behaviours also in the case of exponential delay distribution, as observed in Figure 4.

In Figure 4a we can see that both $L$ and $N$ populations oscillate with similar frequency (roughly 1 cycle every 23 days), but with $L$ having a much larger amplitude than $N$. The $A$ concentration does also oscillate, but with a very small amplitude (less than $0.02\%$). Interestingly, it is possible to find a parameter combination such that an oscillatory behaviour, with a large frequency (roughly one cycle every 8 days) but small amplitude, can only be observed in $L$ (Figure 4b), while $A$ and $N$ appear to have reached a steady state.

As we have varied all the parameters to produce Figure 4 isolating the key parameters is very difficult, but we can observe that those parameters that are similar in Figures 4a and 4b, but different from the baseline parameters are $\gamma$, which is smaller in Table 1, and $\kappa$, which is much larger in 1. These parameters regulate the production of androgen, either from other endocrine glands or from $N$ cells, further emphasising that the androgen dynamics plays a critical role in the survival and growth of the PCa tumour. A further investigation showed that by choosing the scale parameter $\psi = 0.055$, and keeping all other parameters but the delay equal to the baseline values (Table 1), a Hopf bifurcation occurs when $\tau \simeq 11.8$ (Not shown).

5.   Conclusions

In this paper we proposed a new nonlinear distributed-delay dynamical system to model neuroendocrine transdifferentiation and the dynamics of human prostate cancer in in vivo conditions. The main aim of this work was to analyse the effect of the delay in the development and maturation of the human prostate cancer cell population, and of the delay distribution, on prostate cancer growth. The model was developed as an extension of the work by Cerasuolo et al. ^[20], whose strengths and weaknesses were discussed in detail by Turner et al. ^[24].

Analytically, we showed that key parameters for the biological relevance of the system are those regulating the final size of the $L$-cell population, the carrying capacity $\eta_k$ and the maximum proliferation rate $\beta_p$, and the transdifferentiation rates $k_t$ and $k_p$. The stability properties of the tumour-free equilibrium $E_1$ were explored fully through a local stability analysis, and the delay $\tau$ was found to be a bifurcation parameter when the kernel of the distributed delay was modelled as a uniform distribution, with the critical value of the bifurcation parameter identified. Interestingly, it was observed that the presence of the delay initially stabilises the system, but as the delay increases the tumour-present equilibrium loses its newly acquired stability properties. Finally, numerical examples were produced to highlight the different possible behaviours that can be obtained with different delay distributions. We found that also in the case of exponential distribution changes in the delay together with changes in the scale parameter cause a Hopf bifurcation.

The modifications introduced in this model to the original system studied in ^[24] have significantly changed the behaviour of the system. The most evident changes are the reduction of the number of possible equilibria from five to two, the loss of a trivial equilibrium (all zero coordinates), and the existence of a tumour-free equilibrium whose stability properties depend on the parameter values, but it is independent of the delay. The changes to the system equations in (2.1), to incorporate the maximum androgen concentration of the system, the minimum concentration for tumour growth, $A_{max}$ and $A_{min}$, and the external androgen production rate, $\gamma$, allowed us to explore further the effect of androgen concentration on the existence of the equilibria. As the stability of $E_1$ is reliant on $A_1$, this suggests that in androgen-depleted conditions the extinction of the PCa tumour is linked to the underlying androgen dynamics, and in particular to the production of androgen from the other endocrine glands.

The stability of the tumour-present equilibrium was found to be dependent on the choice of the distribution function. In the case of a uniform distribution (system (3.18)), the local stability analysis showed that the delay $\tau$ is a bifurcation parameter, and the numerical simulations confirmed this finding and showed that the system undergoes a supercritical Hopf bifurcation as the delay $\tau$ increases. Also, we observed that as $\tau$ increases the frequency of the oscillations decreases and the amplitude increases, which indicates that the longer the LNCaP cells take to fully mature, the more volatile the PCa tumour's growth becomes.

The analytical study of the system with exponential distribution (3.19) was very involved, didn't lead us to any conclusive result and was therefore omitted. However, numerical simulations showed that $\tau$ is a bifurcation parameter for the stability of the tumour-present equilibrium only when changed together with the scale parameter $\psi$. Exploring the parameter ranges further, we were able to find examples of oscillatory behaviour with the exponential delay distribution, with indication that also $\gamma$ and $\kappa$ can potentially be bifurcation parameters (not shown).

While the model proposed in this work successfully extends the one developed in ^[20] to the case of in vivo conditions, the new system does have its limitations. Whilst the introduction of a logistic growth term for LNCaP cells and the contribution of inter-cellular resource competition to the death rate for NE cells has brought a concept of limited resources (space, nutrients, etc.) to the system, it does not consider the spatial distribution of the cells and the effect this has on resource availability. Also, currently there is no experimental evidence of the oscillatory behaviours displayed in the system by neuroendocrine and LNCaP tumour cells, as there are no studies that produced time series of the dynamic evolution of the transdifferentiation in patients (or clonal cells) over an extended period of time. This is probably due to the fact that neuroendocrine tumours are difficult to diagnose as they usually occur in patients with multiple metastases, and in this condition clinicians are discouraged from performing biopsies ^[36]. Another limitation is that while replacing the proliferation rate function with $F(A)$ and introducing the logistic growth function has helped in the simplification of the model, $\alpha(A)$ has remained unchanged, which means that also in this case as in ^[24] its transcendental nature caused complications in the analysis of the system. Further studies could focus on gaining additional insights into the androgen-dependent mechanisms regulating LNCaP cell transdifferentiation, and look for ways to simplify $\alpha(A)$ while retaining biological plausibility in the process.

Conflict of interests

The authors declare there is no conflicts of interest.

A.   Appendix

A.1.   Uniform Distribution

In the case of the uniform distribution (3.18), we augment the system by introducing the new state variable

where $s = t+\theta$. Finding the derivative with respect to $t$ gives

Hence, system (2.2) becomes

with the following initial conditions for the original state variables

and following ^[33] we define at $t = 0$ the initial condition for $U$, which is

A.2.   Exponential Distribution

For the exponential distribution (3.19) we introduce the new state variable

where $s = t+\theta$. Therefore,

Then system (2.2) becomes

with the initial conditions (A.2) for the original state variables and following ^[33] we define at $t = 0$ the initial condition for $U$, which is

Since distribution kernels can be assumed, without loss of generality, to be positive-definite and normalised to unity ^[37], we take $\Delta = 1$. We then find the following relationship between $\zeta, \, \psi$ and $\tau$.

thus,

hence

This relationship is used in 4.1 to define the parameter $\zeta$ in the simulations.