HTML

ePub

PDF

1 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 de- signed 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 tumour- free 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.

2 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: τ₁, the average cell cycle duration for L(t) cells, and τ₂, 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 as xτ_i for x and . All parameters are assumed to be positive.

Androgen is depleted with constant rate , and is secreted by 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 , 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 , 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. 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, 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 is a decreasing function in L such that 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 , 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 τ₂ 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 functions and for , where with and satisfying

and

We can therefore restrict the study to the model with only three equations:

By observing that measures the loss of second generation cells during mitosis by differentiation, it is clear that we need to assume .

Remark 1. We can see that holds as long as for all t, and that a sufficient condition for is that

where

3 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 with . Hance, C denotes the Banach space of continuous functions mapping the interval into with the supremum norm and is any norm in [4].

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

Let be such that for with and and for . Then (1) at becomes

We will now determine the sign of . From (3), at t0, 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 on the closed interval , which means that . Therefore, on the interval , we have that is negative and increasing, which is in contradiction to the initial hypothesis , 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 and such that for , and for . Therefore (2) at becomes.

This is a contradiction of for, 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 bounded- ness, we will show that they are bounded above.

Lemma 2.Provided that

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 , and its derivative is

which, with substitution of from (7), becomes

Therefore from we have that

If (9) holds and as L is positive and is a decreasing function in L such that it follows that for 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 , then it follows that Z(t) is bounded, which implies that L is bounded on .

Let us now consider equation (3) for N . Using the variation of constant formula, we obtain

That. Is

which gives

where . We conclude that N is bounded for all .

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

That is

which is

where . Therefore, we see that l, and that A is bounded for all . This concludes the proof.

Remark 3. It is important to note that if , then condition (9) is automatically satisfied. Also, by solving (9) for τ₁ 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

4.1 Existence of equilibria

The general equilibrium of system (7) can be found by solving the following equations:

From the second equation of (10) we obtain that either , which leads us to the tumour-free equilibrium , or that the following must be satisfied:

This can be solved for , giving

By observing that by definition for all L, A, then in order to be able to solve equation (11), we must assume that , which is satisfied if . We denote the right-hand side of (11) by , then yields the other components of as functions solely dependent on . Solving the third equation of (10) for gives us

which, via the first equation of (10), yields to

From (11) we have that at an internal equilibrium , and by substituting expression (12) obtained for 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

and

which implies that, like the Ricker function, has a maximum at . If we now analyse the left-hand side of (13), we can see that

Also, if we write

then

From the expression of the derivative for

and since has a maximum at , the following result holds.

Proposition 1. If , then there exists a value such that the function

then is monotonically increasing in A.

From Proposition 1 and by observing that we can see that the function will initially increase and then tend to zero as , which ensures the existence of at least one maximum for . Given these observations, a sufficient condition for the existence of internal equilibria (tumour-present) would be . 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 . If inequality (9) is satisfied, then E₀ is the only equilibrium.

(ii) If

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

4.2 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 is the identity matrix, and X, Y and Z are the Jacobian matrices with respect to the non-delayed, τ₁-delayed and τ₂-delayed variables, respectively (the explicit expression of the characteristic equation can be found in the Supplemental material).

Theorem 1.The tumour-free state E₀ is always locally asymptotically stable.

Proof. The characteristic equation (16) computed at E₀ reduces to

which gives

Hence, since , the tumour-free state E₀ is always stable.

Remark 5. We observe that the stability of the tumour-free state E₀ 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.

4.2.1 Global stability of E₀

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₀ 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₀ as .

Proof. Consider the function

Along the solutions of (7), we obtain

We then consider the functional

Finally, we define a Liapunov functional

Since

then

giving

Since Lemma 1 ensures the positivity of maximals, when condition (9) holds, then we obtain that . Therefore, the maximum invariant set in V(t)7=0 is E₀, and by applying the Liapunov-LaSalle-type theorem for delayed systems [14] we get as (t), thus E₀ is globally attactive.

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 of system (7) is globally asymptotically stable in

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

5.1 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 (kt 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 v₃.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 kp, respectively. Results from the bifurcation analy- sis 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.

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₀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 E0. When κ 0.00055, the stable tumour-present equilibrium undergoes a Hopf bifurcation becoming unstable, and a limit cycle appears. When κ ≈ Importar tabla ≈ Importar tabla 0.00077, a new saddle-node bifurcation occurs, again producing two additional tumour- present equilibria, one stable and one unstable. As κ increases, the new stable tumour- present 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 kp 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 kp increases through 0.347, a saddle node bifurcation occurs. The stable tumour-present equilibrium undergoes a Hopf bifurcation as kp increases, leading to the appearance of an unstable tumour-present equilibrium and a limit cycle. Increasing kp 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₀ 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 kp, increasing a decreases the number of tumour-present equilibria with a saddle node bifurcation at a 1.691.

5.2 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 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 when , and the initial conditions to be , 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 (Si) and total order sensitivity (STi ) 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 (eFAST columns) represent the ranks of STi in decreasing order (1 is the most sensitive).

The PRCC of the parameters were also computed dynamically and plotted against time as shown in Figs. S4(a), S4(b), S6(a) and 5(d) (see Supplemental material), which demonstrate the variation in PRCC over time with respect to the 10% and 20% parameter ranges. Across both parameter ranges the greatest change occured in the first 100 days with all parameters reaching a constant PRCC at around 150 days. The time plots of L and N across the 200 runs, Figs. S4(c), S4(d), S6(c) and S6(d) show that all the solutions head to a stable tumour-present equilibrium.

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, γ, τ₁ and τ₂ being not significant for L, and b for N . Moving to the 20% range, seen in Figs. S 5 (a) and S 5 (b), we find that N has no not-significant parameters, and k_p and τ₂ 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 kp, γ and β₀ for the 10% range and γ, b and τ₁ 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, S 7 (a) and S 7 (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 kp, γ, µ, b and τ₂ 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 τ₁ 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).

5.3 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 bifurcationswith changes in both the number of equilibria and their stability properties. In this sub-section, we will explore numerically the dynamics of system (7) for κ= 0.00017,κ=0.0009,κ= 0.003and κ= 0.005. Unless otherwise stated, we use the parameter valuesfrom Table 1 with final time at 1000 days. We consider the following constant history values φ(ω) = (10,0.1,0), when , to replicate the high androgenconcentration environment used to prepare the cancer cells before the experiment [5], and

the constant initial values . In each figure, we present solutions for A. L and N with different initial values and the corresponding phase portrait.

If κ = 0.00017, then condition (9) is satisfied and the tumour-free equilibrium, E₀, 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- cell. If condition (9) is not satisfied, then the system converges to a stable tumour-present equilibrium for all 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₀ 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₀ is still in its basin of attraction.

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 τ₂ 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₀, 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.

6 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₀, 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 anal- ysis 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 N - cells, 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 τ₁and τ₂ 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₀ 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 counterintuitive 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 sig- nificant 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

Acknowledgments

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.

References

1. M. Adimy, F. Crauste, C. Marquet, Asymptotic behaviour and stability switch for a mature- immature model of cell differentiation, Nonlinear Anal., Real World Appl., 11(4):2913–2929, 2010, https://doi.org/10.1016/j.nonrwa.2009.11.001.

2. M. Adimy, F. Crauste, S. Ruan, Modelling hematopoiesis mediated by growth factors with applications to periodic hematological diseases, Bull. Math. Biol., 68(8):2321–2351, 2006, https://doi.org/10.1007/s11538-006-9121-9

3. F. Bray, J. Ren, E. Masuyer, J. Ferlay, Global estimates of cancer prevalence for 27 sites in the adult population in 2008, Int. J. Cancer, 132(5):1133–1145, 2012, https://doi.org/10.1002/ijc.27711.

4. B. Buonomo, M. Cerasuolo, The effect of time delay in plant-pathogen interactions with host demography, Math. Biosci. Eng., 12(3):473–490, 2015, https://doi.org/10.3934/mbe.2015.12.473

5. M. Cerasuolo, D. Paris, F.A. Iannotti, D. Melck, R. Verde, E. Mazzarella, A. Motta, A. Ligresti, Neuroendocrine transdifferentiation in human prostate cancer cells: An integrated approach, Cancer Res., 75(15):2975–2986, 2015, https://doi.org/10.1158/0008-5472. CAN-14-3830.

6. S.E. Eikenberry, J.D. Nagy, Y. Kuang, The evolutionary impact of androgen levels on prostate cancer in a multi-scale mathematical model, Biol. Direct, .:24, 2010, https://doi.org/10.1186/17456150-5-24.

7. K. Engelborghs, T. Luzyanina, G. Samaey, G. Roose, K. Verheyden, http://twr.cs.kuleuven.be/research/software/delay/ddebiftool.shtml.

8. C.A. Heinlein, C. Chang, Androgen receptor in prostate cancer, Endocr. Rev, 25:276–308, 2004, https://doi.org/10.1210/er.2002-0032.

J.S. Horoszewicz, S.S. Leong, T. Ming-Chu, Z.L. Wajsman, M. Friedman, L. Papsidero, U. Kim, L.S. Chai, S.K. KakatiS., Arya, A.A. Sandberg, The LNCaP cell line – A new model for studies on human prostatic carcinoma, Prog. Clin. Biol. Res., 37:115–132, 1980.

10. A. Ideta, G. Tanaka, T. Takeuchi, K. Aihara, A mathematical model of intermittent androgen suppression for prostate cancer, J. Nonlinear Sci., 18:593–614, 2008, https://doi.org/10.1007/s00332-008-9031-0.

11. T.L. Jackson, A mathematical investigation of the multiple pathways to recurrent prostate cancer: Comparison with experimental data, Neoplasia, .:697–704, 2004, https://doi.org/10.1593/neo.04259.

12. T.L. Jackson, A mathematical model of prostate tumor growth and androgen-independent relapse, Discrete Contin. Dyn. Syst., Ser. B, .:187–201, 2004, https://doi.org/10.3934/dcdsb.2004.4.187.

13. P. Koivisto, M. Kolmer, T. Visakorpi, O.P. Kallioniemi, Androgen receptor gene and hormonal therapy failure of prostate cancer, Am. J. Pathol., 152:1–9, 1998.

14. Y. Kuang, J.D. Nagy, S.E. Eikenberry, Introduction to Mathematical Oncology, Chapman & Hall/CRC, New York, 2016.

15. Y. Kuang, J.D. Nagy, J.J. Elser, Biological stoichiometry of tumor dynamics: Mathematical models and analysis, Discrete Contin. Dyn. Syst., Ser. B, .:221–240, 2004, https://doi.org/10.3934/dcdsb.2004.4.221.

16. S. Marino, I.B. Hogue, C.J. Ray, D.E. Kirschner, A methodology for performing global uncertainty and sensitivity analysis in systems biology, J. Theor. Biol., 254(1):178–198, 2008, https://doi.org/10.1016/j.jtbi.2008.04.011.

17. J.D. Morken, A. Packer, R.A. Everett, J.D. Nagy, Y. Kuang, Mechanisms of resistance to intermittent androgen deprivation in patients with prostate cancer identified by a novel computational method, Cancer Res., 74:3673–3683, 2014, https://doi.org/10. 1158/0008-5472.CAN-13-3162.

18. World Health Organisation, https://gco.iarc.fr/today.

19. V. Perrot, Neuroendocrine differentiation in the progression of prostate cancer; an update to recent developments, Open J. Urol., .:173–182, 2012, https://doi.org/10.4236/oju.2012.223032

R.S. Rittmaster, A.P. Manning, A.S. Wright, L.N. Thomas, S. Whitefield, R.W. Norman, C.B. Lazier, G. Rowden, Evidence for atrophy and apoptosis in the ventral prostate of rats given the 5 alpha-reductase inhibitor finasteride, Endocrinology, 136:741–748, 1995, https://doi.org/10.1210/endo.136.2.7835306.

C.N. Shen, Z.D. Burke, D. Tosh, Transdifferentiation, metaplasia and tissue regeneration, Organogenesis, 1(2):36–44, 2004, https://doi.org/10.4161/org.1.2.1409.

22. M. Stein, Large sample properties of simulations using latin hypercube sampling, Technometrics, 29(2):143–151, 1987, https://doi.org/10.1080/00401706. 1987.10488205.

23. K.R. Swanson, L.D. True, D.W. Lin, K.R. Buhler, R. Vessella, J.D. Murray, A quantitative model for the dynamics of serum prostate-specific antigen as a marker for cancerous growth: An explanation for a medic anomaly, Am. J. Pathol., 163:2513–2522, 2001, https://doi. org/10.1016/S0002-9440 (10)64691-3.

24. S. Terry, H. Beltran, The many faces of neuroendocrine differentiation in prostate cancer progression, Front. Oncol., .:1–9, 2014, https://doi.org/10.3389/fonc.2014. 00060.

25. R.T. Vollmer, S. Egaqa, S. Kuwao, S. Baba, The dynamics of prostate antigen during watchful waiting of prostate carcinoma: A study of 94 japanese men, Cancer, 94:1692–1698, 2002, https://doi.org/10.1002/cncr.10443.

26. R.T. Vollmer, P.A. Humphrey, Tumor volume in prostate cancer and serum prostate-specific antigen: Analysis from a kinetic viewpoint, Am. J. Pathol., 119:80–89, 2003, https://doi.org/10.1309/UNAQ-JTFP-B1RQ-BQD4.

27. J. Wu, R. Dhingra, M. Gambhir, J.V. Remais, Sensitivity analysis of infectious disease models: Methods, advances and their application, J. R. Soc. Interface, 10(86):2012–2018, 2013, https://doi.org/10.1098/rsif.2012.1018.

IR