Mathematical insights into neuroendocrine transdifferentiation of human prostate cancer cells

Prostate cancer represents the second most common cancer diagnosed in men and the fifth most common cause of death from cancer worldwide. In this paper, we consider a nonlinear mathematical model exploring the role of neuroendocrine transdifferentiation in human prostate cancer cell dynamics. Sufficient conditions are given for both the biological relevance of the model’s solutions and for the existence of its equilibria. By means of a suitable Liapunov functional the global asymptotic stability of the tumour-free equilibrium is proven, and through the use of sensitivity and bifurcation analyses we identify the parameters responsible for the occurrence of Hopf and saddle-node bifurcations. Numerical simulations are provided highlighting the behaviour discovered, and the results are discussed together with possible improvements to the model.


Introduction
Prostate cancer (PCa) is the second most common cause of cancer among men worldwide [3,18]. 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, 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. [9]. Early mathematical models investigated the link between the concentration of serum prostate specific antigen (PSA) and tumour volume [23, 25,26]. Others explored how different chemical resources would affect tumour growth, e.g. Kuang et al. who proposed the KNE model [15], which applies ideas from ecological stoichiometry to tumour growth, and considers that tumour growth is limited by various physical restraints 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 [8,13,20], 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) [11,12]. 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. [10], who showed that using intermittent ADT can reduce the time to cancer relapse. While this appears counter-intuitive, they still stated 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. [6] 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 [24] by secreting various hormones and proteins [19]. 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 [24]. 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 [21]. Once this switch occurs, neuroendocrine cells are believed to secrete androgen or similar anabolic hormones, which promote tumour growth as androgen receptor signalling can be detected even under ADT.
In 2015, Cerasuolo et al. [5] proposed a discrete delay dynamical system to investigate the theory of 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, as shown in Fig. S1. The mathematical model was inspired by previous work on cell differentiation in hematopoiesis [1,2] 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 [6,17]. In [5] 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. The paper by Cerasuolo et al. did not include any mathematical analysis of the model, which instead has been performed in this paper. Here, we begin by describing the model formulation and finding sufficient conditions for the non-negativity and boundedness of the system solutions (Section 3). In Section 4, after showing the existence of the tumourfree equilibrium, as well as the possible existence of multiple tumour-present equilibria, we study local and global stability of the tumour-free equilibrium. Using sensitivity analysis and bifurcation analysis (Section 5), we identify the bifurcation parameters and the stability switches that can occur. Finally, we illustrate our findings with numerical simulations, and in Section 6, we discuss the ways in which the model can be improved mathematically without losing its biological relevance.

Model formulation
The model considers two cell populations, L(t), the LNCaP androgen-dependent cells and the transdifferentiated non-malignant neuroendocrine androgen-independent cells, N (t). The growth of L(t) cells is affected by changes in the androgen concentration in the environment, A(t), as it influences their ability to proliferate. In the experiments the androgen was introduced into the Petri dishes through a charcoal stripped serum, and as this was the only external source of androgen, the serum could be considered equivalent to the initial androgen concentration. The N (t) cells are assumed to be androgen independent and post-mitotic; they are unable to proliferate and can only be produced by transdifferentiation of L(t) cells. The L(t) cells go through asymmetric cell division, which implies that upon undergoing proliferation a proportion of the daughter cells will be differentiated. The L(t) cells are divided into three compartments: mature/resting cells, L(t), proliferating cells, P (t), and transdifferentiating cells, T (t). Note that the quantities P (t) and T (t) denote different phases of the L(t) cells' life cycle, and do not therefore represent new cell types. Two discrete delays are considered: τ 1 , the average cell cycle duration for L(t) cells, and τ 2 , representing the duration of transdifferentiation from L(t) to N (t) cells.
For brevity in the following, we will denote x(t) as x and x(t − τ i ) as x τi for x ∈ {A, L, N, P, T } and i = 1, 2. All parameters are assumed to be positive.
Androgen is depleted with constant rate ϕ, and is secreted by N cells with secretion rate, κ: The L cell death occurs at a constant per capita mortality rate, δ. The differentiation to N happens with the differentiation efficiency k t , and depends on the concentration of androgen in the environment and, as suggested by experimental evidence [5], is regulated by the Ricker function α(A) = rAe −aA , where r is the gradient of the differentiation increase, and a the inverse of the maximum differentiation rate. Resting cells become proliferating with an introduction rate β(L, A), a continuous function that is zero in the absence of androgen in the medium, increases as the androgen concentration increases, and decreases as the L-cell concentration increases. β(L, A) is expressed as the product of two terms, the first term represents the inhibition of the mitotic reentry rate (the rate at which L cells become proliferating) due to cell density, and is described by a Hill function; while the second term represents the dependence of the introduction rate on the amount of androgen in the medium and is represented by a Michaelis-Menten function. The β function is based on the β functions used by Adimy et al. in [1] and [2]: Here, β 0 is the maximum rate of cell transfer from the resting phase to the proliferating phase, b is the half-saturation constant for androgen concentration, and θ and n have similar roles to the Hill coefficients and represent the response to LNCaP population changes. Finally, we note that β(L, A) is a decreasing function in L such that 0 β(L, A) < β 0 for all values of L and A.
The L cells are generated through asymmetrical cell division. This means that rather than mitosis always resulting in two L cells, a proportion of daughter cells are generated as N cells. This proportion is regulated by the function α(A τ1 ), where the delay related to cell cycle duration is taken into account, and by k p , i.e. the proliferation rate of L. The equation for L cells is given by The N -cells die at a constant per capita mortality rate µ, they are generated through the asymmetrical cell division of P -cells and by the differentiation of L-cells. The equation for N cells is given by The transdifferentiating and proliferating cell populations satisfy the equations where τ 2 represents the time T -cells take to become fully differentiated into N cells, and γ is the death rate of the proliferating cells. We observe that equations (1), (2) and (3) are decoupled from (4) and (5) as neither T nor P appear in the equations of A, L and N . Moreover, equations (4) and (5) can be solved explicitly as follows: For biological relevance, we consider the initial conditions as positive continuous func- and N (0) 0. We can therefore restrict the study to the model with only three equations: https://www.journals.vu.lt/nonlinear-analysis By observing that k p α(A τ1 ) measures the loss of second generation cells during mitosis by differentiation, it is clear that we need to assume 0 1 − k p α(A τ1 ) 1.

Remark 1.
We can see that k p α(A τ1 ) 0 holds as long as A 0 for all t, and that a sufficient condition for k p α(A τ1 ) 1 is that where α max = α(1/a) = re −1 /a.

System properties
In this section, we investigate some basic properties of the solutions to system (7), in particular non-negativity and boundedness. The initial conditions for system (7) are the positive continuous functions φ(ω) . To understand the parameter values for which model (7) is biologically relevant, we will study non-negativity and boundedness of the solutions. We prove the following results.
Lemma 1. If k p ae/r, then any solution of the system given by equations (6) and (7) remains non-negative whenever it exists.
Proof. It is possible to prove that A is non-negative by contradiction.
We will now determine the sign of N (t 0 ). From (3), at t 0 , we have and by using the variation of constant formula we get that is Hence, as all initial conditions are non-negative and as the integrals are strictly positive and continuous functions, we have that N (t 0 ) 0 on the closed interval [0, t 0 ], which means that dA(t 0 )/dt 0. Therefore, on the interval (t 0 , t 0 + ), we have that A(t) is negative and increasing, which is in contradiction to the initial hypothesis A(t 0 ) = 0, and therefore A is non-negative.
By contradiction and using the non-negativity of A it is possible to prove that L is non-negative.
Suppose there exists t 0 > 0 and 0 This is a contradiction of L(t) < 0 for t ∈ (t 0 , t 0 + ), and therefore L is non-negative. The non-negativity of N can be proved using the same reasoning as for L.

Remark 2.
We observe that the positivity of A and L imply the positivity of both P and T by considering recurrence arguments applied to the integral forms (6).
Given the non-negativity of the solutions (Lemma 1), in order to prove their boundedness, we will show that they are bounded above.
and Lemma 1 is satisfied, then the solutions of system (7) are bounded.
Proof. Let We can observe that Z(t) is a differentiable function in [0, +∞), and its derivative iṡ which, with substitution ofL from (7), becomeṡ https://www.journals.vu.lt/nonlinear-analysis Therefore from 0 α(A) α max we have thaṫ If (9) holds and as L is positive and β(L, A) is a decreasing function in L such that β(L, A) < β 0 , it follows thatŻ(t) 0 for t 0 and therefore Z(t) is decreasing. As Z(t) is a non-negative function, as the integral is of a strictly positive and continuous function on the closed interval of [t, t + τ 1 ], then it follows that Z(t) is bounded, which implies that L is bounded on [0, ∞).
Let us now consider equation (3) for N . Using the variation of constant formula, we obtain That is where L * = sup [0,+∞) L(t). We conclude that N is bounded for all t ∈ [0, +∞).
We will now prove that A is bounded. As we did for N , by considering the variation of constant formula applied to equation (1) we get that

Remark 3.
It is important to note that if 2e −γτ1 − 1 0, then condition (9) is automatically satisfied. Also, by solving (9) for τ 1 we are able to find a lower bound for the first delay, which ensures that the model is biologically relevant, giving 4 Equilibria and their stability analysis

Existence of equilibria
The general equilibrium of system (7)Ē = (Ā,L,N ) can be found by solving the following equations: From the second equation of (10) we obtain that eitherL = 0, which leads us to the tumour-free equilibrium E 0 = (0, 0, 0), or that the following must be satisfied: This can be solved for β(L,Ā), giving By observing that by definition β(L, A) 0 for all L, A, then in order to be able to solve equation (11), we must assume that . We denote the right-hand side of (11) by β * (Ā), then (10) yields the other components ofĒ as functions solely dependent onĀ. Solving the third equation of (10) forN gives us which, via the first equation of (10), yields tō .
From (11) we have that at an internal equilibrium β(L,Ā) = β * (Ā), and by substituting expression (12) obtained forL we can generate an equation inĀ, which when solved will give us the equilibria of (7). This leads to the transcendental equation We observe that which implies that, like the Ricker function, β * (Ā) has a maximum at 1/a. If we now analyse the left-hand side of (13), we can see that .
From the expression of the derivative forL and since β * has a maximum at 1/a, the following result holds.
thenL is monotonically increasing in A.
From Proposition 1 and by observing that f (0) = 0 we can see that the function β(L,Ā) will initially increase and then tend to zero as A → ∞, which ensures the existence of at least one maximum A * A 1 for β(L,Ā). Given these observations, a sufficient condition for the existence of internal equilibria (tumour-present) would be g(L(1/a))f (1/a) > β * (1/a). However, this would be far from an optimal sufficient condition as will be shown with a numerical example below. The above results are summarised in the following proposition. Proposition 2. System (7) admits the following equilibria: (i) For all parameter values, system (7) admits the tumour-free equilibrium E 0 = (0, 0, 0). If inequality (9) is satisfied, then E 0 is the only equilibrium.
then (7) will admit at least two tumour-present equilibria provided that and condition (8) holds.
In Fig. S2, we observe, by plotting the left-hand side, which we relabelβ(Ā), against the right-hand side β * (Ā) whilst varyingĀ, that there can be up to 4 intersections between the two functions. As β * (0) = 0 ( Fig. S2 inset), each of the intersections represent distinct tumour-present equilibria.

Remark 4.
It is possible to observe that by changing the parameter value κ while keeping the other parameters fixed, the number of tumour-present equilibria changes through 0, 1, 2, 3, and 4 (not shown).
In Section 5, we explore in detail how the number of equilibria and their stability properties change depending on parameter values.

Stability analysis
Here, we derive results concerning the local stability of the tumour-free equilibrium that will additionally be used in Section 4.2.1 in order to prove its global stability. Local stability of the equilibria is determined by the characteristic equation, i.e.
where I ∈ M 3×3 (R) is the identity matrix, and X, Y and Z are the Jacobian matrices with respect to the non-delayed, τ 1 -delayed and τ 2 -delayed variables, respectively (the explicit expression of the characteristic equation can be found in the Supplemental material).
Theorem 1. The tumour-free state E 0 is always locally asymptotically stable.
Remark 5. We observe that the stability of the tumour-free state E 0 is not dependent on either of the delays of system (7).
Because of the complexity of the explicit form of the characteristic equation and the impossibility to get explicit expressions of the tumour-present states, the stability of the latter will be analysed numerically.

Global stability of E 0
We will now prove that under the sufficient condition that ensures the boundedness of solutions and the uniqueness of the tumour-free equilibrium, it is possible to prove that E 0 is globally asymptotically stable.
Theorem 2. If condition (9) holds, then all solutions (A, L, N ) of (7) converge to the tumour-free state E 0 as t → ∞.
From the results of local stability in Theorem 1 and the global attractivity in 2 we obtain the following corollary. Corollary 1. If condition (9) holds, then the tumour-free equilibrium E 0 = (0, 0, 0) of system (7) is globally asymptotically stable in R 3 + .

Numerical results
In this section, we provide numerical results that allow us to explore system (7) further. First, we use bifurcation analysis to analyse different system behaviours with respect to the parameters that cannot be estimated experimentally, then we run a sensitivity analysis to establish which parameters contribute most to the uncertainty in the model outputs, and finally, we present some numerical simulations under various parameter conditions.

Bifurcation analysis
In [5] the authors show that only three of the parameters in system (7) cannot be estimated by designing suitable experiments, these are k t , k p and κ. Therefore, such parameters could prove to be mathematically interesting, and we explored the behaviour of (7) as k t , k p and κ vary (k t is not shown as it exhibits the same bifurcation behaviour as k p ). In order to investigate the parameters κ and k p for bifurcations, we used the software DDE Biftool v3.1.1 [7] with the parameter values as in Table 1 and with the ranges [0.00001, 0.02] and [0, 0.7] for κ and k p , respectively. Results from the bifurcation analysis are shown in Figs. 2(a) and S3, with Fig. 2(a) showing a shorter range as the behaviour of κ does not change for κ > 0.01. The k p range ensures that condition (8) was satisfied.    Table 1. Red indicates an unstable equilibrium, an unstable equilibrium with a stable limit cycle, and green is a stable equilibrium. The pink circle indicates the presence of a saddle-node bifurcation, and the black square indicates a Hopf bifurcation. (b) shows the unstable tumour-present equilibrium near the tumor free state E 0 .
As we can see from Fig. 2(a), varying κ in the range [0.00001, 0.005] greatly affects the number of tumour-present equilibria as well as their stability properties. For very small values, κ < 0.00024, only the tumour-free state E 0 exists, but increasing κ leads to a saddle node bifurcation, which results in two tumour-present equilibria, one stable and one unstable. Figure 2(b) shows that the unstable tumour-present equilibrium heads very close to the tumour-free state E 0 . When κ ≈ 0.00055, the stable tumour-present equilibrium undergoes a Hopf bifurcation becoming unstable, and a limit cycle appears. When κ ≈ 0.00077, a new saddle-node bifurcation occurs, again producing two additional tumourpresent equilibria, one stable and one unstable. As κ increases, the new stable tumourpresent equilibrium undergoes a Hopf bifurcation, leading to the appearance of a limit cycle (κ ≈ 0.00096), but continuing to increase κ leads to a second Hopf bifurcation, returning to a stable tumour-present equilibrium when κ ≈ 0.0047. The second unstable tumour-present equilibrium does not change behaviour, while κ increases. Figure S3(a) shows that changing k p does not change the stability or behaviour of the largest stable tumour-present equilibrium or the smallest unstable tumour-present equilibrium (which is shown in Fig. S3(b)). As k p increases through 0.347, a saddle node bifurcation occurs. The stable tumour-present equilibrium undergoes a Hopf bifurcation as k p increases, leading to the appearance of an unstable tumour-present equilibrium and a limit cycle. Increasing k p further does not change the stability or behaviour of any of the tumour-present equilibria.
Both Figs. 2(a) and S3(a) confirm that the tumour-free equilibrium E 0 is always stable.
A bifurcation analysis was also performed for r, a and θ (not shown). The bifurcation diagrams for a and r, Figs. S3(c) and S3(d) (see Supplemental material), respectively, were obtained with parameter ranges such that condition (8) was satisfied. While with r the system exhibits the same bifurcations and stability changes as with k p , increasing a decreases the number of tumour-present equilibria with a saddle node bifurcation at a ≈ 1.691.

Sensitivity analysis
A global sensitivity analysis (SA) was performed to examine the response of system (7) to parameter variation following the approach by Marino et al. [16]. Since the relationship, including monotonicity, between the parameters and the outputs is not known a priori, two methods were used: (i) Latin Hypercube Sampling/Partial Rank Correlation Coefficient (LHS/PRCC) sensitivity analysis; and (ii) an extended Fourier amplitude sensitivity test (eFAST) [16]. The two methods provide different rankings to the same parameters, so a dual method approach helps to highlight additional parameters of interest. In this way, we were able to reduce some of the uncertainty caused by performing the analysis on a large number of parameters simultaneously.
We considered all parameters taking values in R. Therefore, n, which is the exponent in the β(L, A) function and is assumed to be an integer, was excluded and was taken to be n = 4. The parameter values used for the SA were either derived as means of parameter ranges obtained from experimental data (see [5]) or were taken from literature (Table S1). Parameter ranges were constructed to reflect a 10% and a 20% variation in the data, and the parameters were assumed to be random variables with uniform distributions. We took the history to be the constant values φ(ω) = (10, 0.1, 0) when ω ∈ [− max{τ 1 , τ 2 }, 0), and the initial conditions to be φ(0) = (3.74, 1.95, 17.58), close to the upper stable tumour-present equilibrium shown in Fig. 2. This was done to reflect the experimental conditions [5] and to start close to a stable equilibrium to avoid possible unstable manifolds that could disrupt the results. It should be noted that condition (8) was met for all simulations.
The sensitivity analysis was run against L and N outputs in order to highlight the effect of parameter values on the outcome of the tumour. The PRCC were generated using a Monte Carlo method called the Latin Hyperspace Sampling (LHS) technique [22]. The analysis was performed using the Matlab R2019b solver dde23 over 200 runs, simulating 500 days of PCa growth to allow solutions to converge to a final behaviour. The calculated PRCC values can be seen in Table 2. The parameters with large PRCC values (> 0.5 or < −0.5) statistically have the most influence [27]. However, it should be noted that any parameter with a PRCC value with modulus larger than that of the dummy value may be significant [16].
The eFAST method was used to generate the first order sensitivity (S i ) and total order sensitivity (S Ti ) of each parameter. It works by calculating the variance of the model output based on the input parameters variation. This variance is partitioned to determine what fraction can be explained by each input parameter. Parameters with total order sensitivity less than or equal to the dummy parameter's total order sensitivity are considered to have no significant influence on the outcome of the system [16]. The values in Table 2  From Table 2 and the bar graphs in Figs. 3 we observe that in the 10% range, there are five incidents of non significance from PRCC with k p , γ, τ 1 and τ 2 being not significant for L, and b for N . Moving to the 20% range, seen in Figs. S5(a) and S5(b), we find that N has no not-significant parameters, and k p and τ 2 become significant for L. We should note that increasing the parameter range only caused δ for L to become highly significant with no parameters becoming less significant, and that for N , all of k p , γ and β 0 for the 10% range and γ, b and τ 1 for the 20% range have small enough PRCC values for their significance to be investigated further.
The eFAST rankings from Table 2 and the bar graphs in Figs. 4, S7(a) and S7(b) show that κ and θ are consistently sensitive across both parameter ranges. Increasing the parameter range to 20% causes all other parameters to become insignificant for L, while k p , γ, µ, b and τ 2 change their significance for N with µ becoming the only insignificant parameter for N .
Comparing the PRCC and eFAST results we see that eFAST finds that more parameters are insensitive on both parameter ranges, and that κ, b, θ and τ 1 are the parameters for which both methods record the same results across both variables and parameter ranges. Finally, r is an interesting parameter as for both methods, the significance for L and N does not change when increasing the parameter range, even though the methods do not give the same significance.
From Table 2 it is clear that κ and θ are the most significant parameters with both methods, with ϕ, µ and a consistently highly significant from the PRCC. It is notable that these parameters regulate the reintroduction of mature cells in the proliferating phase (θ),  how much testosterone is made available through N -cells to L-cells for their growth (κ), the androgen degradation rate (ϕ), the death rate of N -cells (µ) and the maximum rate of the differentiation from L to N cells (a).

Numerical simulations
Using the Matlab R2019b solver dde23 with step size tolerance 10 −6 , we simulated the various system behaviours disclosed by the bifurcation analysis as well as the global stability property of the tumour-free state when condition (9) holds. In fact, as Fig. 2(a) illustrates, by varying κ in the interval [0, 0.0047] the system undergoes five bifurcations with changes in both the number of equilibria and their stability properties. In this subsection, we will explore numerically the dynamics of system (7) for κ = 0.00017, κ = 0.0009, κ = 0.003 and κ = 0.005. Unless otherwise stated, we use the parameter values from Table 1 with final time at 1000 days. We consider the following constant history values φ(ω) = (10, 0.1, 0), when ω ∈ [− max{τ 1 , τ 2 }, 0), to replicate the high androgen concentration environment used to prepare the cancer cells before the experiment [5], and  If κ = 0.00017, then condition (9) is satisfied and the tumour-free equilibrium, E 0 , is globally asymptotically stable. From Fig. 5 we can see that the L component approaches the trivial equilibrium faster than both A and N , taking less than 300 days for all initial values, while A and N can take up to 900 days for the largest initial value. A plausible explanation is that N -cells mortality rate, µ, is much lower than the L-cells one, γ, and A is produced by N -cells. If condition (9) is not satisfied, then the system converges to a stable tumour-present equilibrium for all φ(0) considered (Fig. S8).
For κ = 0.0009, it is possible to distinguish three different behaviours of the system as shown in Fig. 6, with one initial value heading to a stable tumour-present equilibrium, one going to the tumour-free state E 0 and the others tending towards a limit cycle, consistent with Fig. 2(a).
When κ = 0.003, we can see that the solutions approach one of two limit cycles, and the initial condition close to the tumour-free state E 0 is still in its basin of attraction. It  is interesting to note that the limit cycle corresponding to higher androgen levels (HAL) has a shorter period, approximately 8 days, compared to the limit cycles corresponding to lower androgen levels (LAL), over 400 days, which also has a longer period compared to the corresponding limit cycle in the case when κ = 0.0009. We observe that in the case of the HAL limit cycle the oscillation only happens along the L-axis and with a period comparable to the τ 2 delay value from Table 1, while for the LAL limit cycle, all components oscillate as shown in Fig. 7. It is worthy of note that when the initial values are such that the trajectories are attracted by the second limit cycle, there are long periods of time (∼ 200 days) in which the tumour cells are almost extinct. In a clinical setting, this would correspond to considering the patient disease-free as N -cells are non-tumorigenic and cannot be detected through the blood-based PCa tumour markers [5].
Finally, for κ = 0.005, all of the initial conditions, including the one close to the tumour-free state E 0 , give rise to trajectories that tend to the stable tumour-present   equilibrium (Fig. 8). Interestingly, for initial conditions such that L < 0.9, we observe a transition phase with large oscillations in L disappearing into the tumour-present equilibrium, suggesting that the limit cycle could have become unstable ( Fig. 2(a)). Increasing κ removes such behaviour as shown in Figs. S4, S5 obtained with the baseline parameters.

Discussion and conclusion
In this paper, we explored the dynamics of the discrete-delay dynamical system proposed by Cerasuolo et al. [5] to investigate the neuroendocrine transdifferentiation of PCa cells grown in vitro in androgen-deprived conditions. The model, based on experimental evidence, had already been compared to biological data, but the authors had not performed an exhaustive analytical and/or numerical study of this dynamical system. Here, analytical properties of system (7) solutions, such as boundedness and positivity, were proved showing that the model has the necessary features to be biologically significant. The local and global stability analysis of the tumour-free equilibrium E 0 , together with the sensitivity and bifurcation analyses provide further insights into the system dynamics. Finally, through the use of the latter numerical methods, we were able to determine which parameters are the most influential on the outcome of tumour growth. Via numerical examples we highlighted the different behaviours of the model.
Interestingly, the most significant parameters from both methods of sensitivity analysis were κ and θ, the production of testosterone from N -cells and the reintroduction of mature cells in the proliferating phase, respectively. Other highly significant parameters were the androgen degradation rate, ϕ, the N -cells death rate, µ, and the maximum transdifferentiation rate, a, which suggest that in androgen-depleted conditions the PCa tumour cells survival is heavily dependent on the rate at which L-cells transform into Ncells, as well as the efficiency of the N -cells at surviving and producing testosterone. However, the bifurcation analysis showed that the most interesting dynamics can be observed by changing κ. In particular, from Fig. 2(a) we can see a very rich dynamics emerging as κ takes values in the range (0, 0.0047). Even a small increase in κ allows to move from a situation of global stability of the trivial equilibrium to one in which the L-cells will almost always survive (κ > 0.00024) depending on the initial conditions. The sensitivity analysis demonstrates that both τ 1 and τ 2 are never highly significant for system (7).
The analytical and numerical results show that when subjected to ADT, tumour L-cells that fail to divide or, equivalently, to move from the mature to the proliferating phase quickly enough, eventually become extinct (Fig. 5). Bi-stability is observed for κ ∈ (0.00024, 0.00055) (Fig. S6), for larger values, oscillatory behaviours appear, until finally when κ 0.0047, solutions will either go towards the tumour-free equilibrium E 0 or towards the stable tumour-present equilibrium. However, what is of note in all of these cases is that as soon as there is enough (even if minimal) release of testosterone from N -cells, the tumour is able to sustain itself, regardless of the lack of androgen imposed by the therapy. It is interesting to note that most of the cyclic behaviours can be observed only for a small window of κ values and that even if counter-intuitive such a survival strategy could be the most successful one for the tumour cells. Certainly, from a therapeutic point of view it is of extreme importance to consider the possibility of being unable to detect the tumour when in fact it persists under a different guise.
From the sensitivity analysis we can conclude that the delays are neither highly significant for both cell populations L and N at the same time, nor that they are consistently significant using either method. This seems to suggest that modifications to system (7) could involve removing one or both delays as well as giving additional emphasis to the sensitive parameters κ, θ, ϕ, µ, a and associated functions. Hence, further studies could focus on the feedback loop regulating the release of testosterone and on gaining additional insights into the androgen-dependent mechanisms regulating LNCaP cell transdifferentiation.
Acknowledgment. The authors would like to thank Professor Edoardo Beretta for his comments on the modelling aspects of this problem, and Dr. Alessia Ligresti for her kind permission to print Fig. S1. They would also like to thank the referees for their very useful and insightful comments that helped improve the manuscript further.

Supplemental material
Linear stability: Characteristic equation Equation (16) can be written as

Scheme of LNCaP neuroendocrine trans-differentiation
Mature NED cells Cell seeding (GM) characterized by multipolar showing the cell morphology changes during the transdifferentiation process [5] (right panel is reprinted from [5]). Figure S2. A plot ofβ(Ā) (blue line) and of β * (Ā) (red dotted line) using parameters from [5]. The parameter values satisfy the necessary condition (15) for the existence of tumour-present equilibria, but not the sufficient one (14).  Table 1. Red indicates an unstable equilibrium, blue an unstable equilibrium with a stable limit cycle, and green is a stable equilibrium. The pink circle indicates the presence of a saddle-node bifurcation, and the black square indicates a Hopf bifurcation; (b) shows the unstable tumour-present equilibrium near the tumor free state E 0 . https://www.journals.vu.lt/nonlinear-analysis