Loading [MathJax]/jax/output/SVG/jax.js
Research article

PCa dynamics with neuroendocrine differentiation and distributed delay


  • Received: 19 July 2021 Accepted: 24 September 2021 Published: 08 October 2021
  • Prostate cancer is the fifth most common cause of death from cancer, and the second most common diagnosed cancer in men. In the last few years many mathematical models have been proposed to describe the dynamics of prostate cancer under treatment. So far one of the major challenges has been the development of mathematical models that would represent in vivo conditions and therefore be suitable for clinical applications, while being mathematically treatable. In this paper, we take a step in this direction, by proposing a nonlinear distributed-delay dynamical system that explores neuroendocrine transdifferentiation in human prostate cancer in vivo. Sufficient conditions for the existence and the stability of a tumour-present equilibrium are given, and the occurrence of a Hopf bifurcation is proven for a uniform delay distribution. Numerical simulations are provided to explore differences in behaviour for uniform and exponential delay distributions. The results suggest that the choice of the delay distribution is key in defining the dynamics of the system and in determining the conditions for the onset of oscillations following a switch in the stability of the tumour-present equilibrium.

    Citation: Leo Turner, Andrew Burbanks, Marianna Cerasuolo. PCa dynamics with neuroendocrine differentiation and distributed delay[J]. Mathematical Biosciences and Engineering, 2021, 18(6): 8577-8602. doi: 10.3934/mbe.2021425

    Related Papers:

    [1] Alacia M. Voth, John G. Alford, Edward W. Swim . Mathematical modeling of continuous and intermittent androgen suppression for the treatment of advanced prostate cancer. Mathematical Biosciences and Engineering, 2017, 14(3): 777-804. doi: 10.3934/mbe.2017043
    [2] Paolo Fergola, Marianna Cerasuolo, Edoardo Beretta . An allelopathic competition model with quorum sensing and delayed toxicant production. Mathematical Biosciences and Engineering, 2006, 3(1): 37-50. doi: 10.3934/mbe.2006.3.37
    [3] Tin Phan, Changhan He, Alejandro Martinez, Yang Kuang . Dynamics and implications of models for intermittent androgen suppression therapy. Mathematical Biosciences and Engineering, 2019, 16(1): 187-204. doi: 10.3934/mbe.2019010
    [4] Avner Friedman, Harsh Vardhan Jain . A partial differential equation model of metastasized prostatic cancer. Mathematical Biosciences and Engineering, 2013, 10(3): 591-608. doi: 10.3934/mbe.2013.10.591
    [5] Shunyi Li . Hopf bifurcation, stability switches and chaos in a prey-predator system with three stage structure and two time delays. Mathematical Biosciences and Engineering, 2019, 16(6): 6934-6961. doi: 10.3934/mbe.2019348
    [6] Tyler Cassidy, Morgan Craig, Antony R. Humphries . Equivalences between age structured models and state dependent distributed delay differential equations. Mathematical Biosciences and Engineering, 2019, 16(5): 5419-5450. doi: 10.3934/mbe.2019270
    [7] Nancy Azer, P. van den Driessche . Competition and Dispersal Delays in Patchy Environments. Mathematical Biosciences and Engineering, 2006, 3(2): 283-296. doi: 10.3934/mbe.2006.3.283
    [8] Xinran Zhou, Long Zhang, Tao Zheng, Hong-li Li, Zhidong Teng . Global stability for a class of HIV virus-to-cell dynamical model with Beddington-DeAngelis functional response and distributed time delay. Mathematical Biosciences and Engineering, 2020, 17(5): 4527-4543. doi: 10.3934/mbe.2020250
    [9] Edoardo Beretta, Dimitri Breda . Discrete or distributed delay? Effects on stability of population growth. Mathematical Biosciences and Engineering, 2016, 13(1): 19-41. doi: 10.3934/mbe.2016.13.19
    [10] Mostafa Adimy, Abdennasser Chekroun, Claudia Pio Ferreira . Global dynamics of a differential-difference system: a case of Kermack-McKendrick SIR model with age-structured protection phase. Mathematical Biosciences and Engineering, 2020, 17(2): 1329-1354. doi: 10.3934/mbe.2020067
  • Prostate cancer is the fifth most common cause of death from cancer, and the second most common diagnosed cancer in men. In the last few years many mathematical models have been proposed to describe the dynamics of prostate cancer under treatment. So far one of the major challenges has been the development of mathematical models that would represent in vivo conditions and therefore be suitable for clinical applications, while being mathematically treatable. In this paper, we take a step in this direction, by proposing a nonlinear distributed-delay dynamical system that explores neuroendocrine transdifferentiation in human prostate cancer in vivo. Sufficient conditions for the existence and the stability of a tumour-present equilibrium are given, and the occurrence of a Hopf bifurcation is proven for a uniform delay distribution. Numerical simulations are provided to explore differences in behaviour for uniform and exponential delay distributions. The results suggest that the choice of the delay distribution is key in defining the dynamics of the system and in determining the conditions for the onset of oscillations following a switch in the stability of the tumour-present equilibrium.



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

    For brevity in the following, we will denote x(t) as x for x{A,L,N}. The equation for the androgen dynamics follows the modelling by Baez and Kuang [28], where Amin represents the minimum androgen concentration required to sustain the tumour, and Amax>Amin, 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 γ and diffusion rate proportional to the difference between Amax and A. Androgen depletion depends on its maximum rate, μ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, κ, which leads to the equation

    dAdt=γ(AmaxA)μA(AAmin)+κN. (2.1)

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

    α(A)=rAeaA.

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

    F(A)=βP(1AminA),

    where βP is the maximum proliferation rate of L and Amin 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 τ [26]. Therefore, the equation for L is

    dLdt=(1kpα(A))F(A)L(10τLθηkω(θ)dθ)δLLktα(A)L,

    where ηk is the carrying capacity of the environment, ω(θ), the delay kernel or delay distribution function and finally Lθ=L(t+θ) 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 μN, and a quadratic term, δNN2, representing the intraspecific competition for space and resources [20,28]:

    dNdt=kpα(A)F(A)L(10τLθηkω(θ)dθ)+ktα(A)LδNN2μNN.

    The final model is

    dAdt=γ(AmaxA)μA(AAmin)+κN,dLdt=(1kpα(A))F(A)L(10τLθηkω(θ)dθ)δLLktα(A)L,dNdt=kpα(A)F(A)L(10τLθηkω(θ)dθ)+ktα(A)LδNN2μNN. (2.2)

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

    A(θ)=ϕ1(θ);L(θ)=ϕ2(θ);N(θ)=ϕ3(θ),

    where ϕ=(ϕ1,ϕ2,ϕ3)TC is such that ϕi(θ)0,(τθ0,i=1,2,3). C denotes the Banach space C([τ,0],R3+) of continuous functions mapping the interval [τ,0] into R3+ with the supremum norm

    ||ϕ||=supθ[τ,0]|ϕ(θ)|

    where || is any norm in R3+ [30].

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

    0<1kpα(A)1A,

    and that αmax=re1a, then we must impose the condition

    kpαmax1. (2.3)

    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

    ˙x(t)=b(t)x(t)[1M(x(t))], (2.4)

    where the following hold.

    (H1) We have

    (i) b:RR is a continuous function, and there are constants β0, β0 such that 0<β0b(t)β0 for all tR;

    (ii) M:CR is a positive linear operator with l:=||M||>0.

    Define x=l1 and consider the IVP for (2.4) with x0=φ, with solution x(φ)(t). Then the following results hold.

    (i) x(φ)(t) is defined for t0, bounded below from zero and bounded on [0,).

    (ii) If x(φ)(t) is non-oscillatory around x, then x(φ)(t)x as t.

    (iii) Moreover, if x(φ)(t) is oscillatory around x and

    ttτb(τ)dτ32for large t,

    then x(φ)(t)x as t.

    Lemma 2.2. Let ϕi(s)0 and ϕ1(s)Amin, with s[τ,0] and i=2,3. If either of the following conditions hold,

    LBΔηk<1,or (2.5)
    LBΔηk1andkt>(LBΔηk1)kpβP, (2.6)

    where LB=sup[0,+)L(t) and Δ=0τω(θ)dθ, then the solutions of (2.2) are non-negative and bounded. In particular, A(t)Amint0.

    Proof. By considering the second equation of (2.2)

    dLdt=L((1kpα(A))F(A)(10τLθηkω(θ)dθ)δLktα(A)),

    with L(0)=ϕ2(0)>0, then

    L(t)=L(0)exp[(1kpα(A))F(A)(10τLθηkω(θ)dθ)δLktα(A)]>0t0.

    If L(0)=0 then L(t)=0t0.

    Let t0>0 be such that A(t0)=Amin, we wish to prove that A(t)Amint>0. Let us assume, on the contrary, that A(t)<Amin for all t(t0,t0+ϵ) with 0<ϵ<τ and for t<t0,A(t)>Amin. Then the first equation of (2.2) at t0 becomes

    dA(t0)dt=γ(AmaxAmin)+κN(t0).

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

    dN(t0)dt=ktα(A(t0))L(t0)δNN(t0)2μNN(t0),dN(t0)dt+N(t0)fN(N(t0))=:B(A(t0),L(t0)),

    where fN(N(t0)):=(δNN(t0)+μN) and B(A(t0),L(t0))=ktα(A(t0))L(t0). Then by the variation of constant formula, we obtain

    N(t0)=N(0)exp[t00fN(N(s))ds]+t00B(A(s),L(s))ds

    and as all initial conditions are non-negative, A(t0)>0 and L(t)0t0, then N(t0)0, and therefore dA(t0)dt0. This is a contradiction of the initial hypothesis, and therefore A(t)Amint0.

    Before proving the positivity of N, we need to prove the boundedness of L. Following Faria and Liz [31], and as L(θ)=ϕ1(θ)0 with θ[τ,0] and L(t)0,A(t)Amint0, from the second equation of (2.2) we obtain

    ˙L<F(A)L(11ηk0τLθω(θ)dθ),

    that is,

    ˙L<βPL(1M(Lθ)),

    where M:CR is a non-zero positive integral linear operator, i.e. M(ϕ2)0ϕ2C,s.t.ϕ20.

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

    ˙L=βPL(1M(Lθ)).

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

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

    dNdt=kpα(A)F(A)L(10τLθηkω(θ)dθ)+ktα(A)LδNN2μNN.

    Let there exist a t0 such that N(t0)=0, N(t)>0 for t<t0, and N(t)<0 for t(t0,t0+ϵ], then

    dN(t0)dt=L(t0)α(A(t0))(kpF(A(t0))+ktkpF(A(t0))ηk0τL(t0+θ)ω(θ)dθ),

    so if either hypothesis (2.5) or (2.6) holds then dN(t0)dt is greater than zero, which is a contradiction of N(t)<0 for t(t0,t0+ϵ]. Therefore, N(t)0 for all t0.

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

    dNdt=kpα(A)F(A)L(10τLθηkω(θ)dθ)+ktα(A)LδNN2μNN.

    As L is bounded, α(A)αmaxA and from (2.3), our equation becomes

    dNdt<βPLB+ktαmaxLBμNN,

    where LB:=supt[0,+)L(t). This can be solved to give

    N(t)<N(0)eμNt+LB(βP+ktαmax)μN(1eμNt),

    and therefore N(t) is bounded above for all t0.

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

    ˙A<γAmax+μAAminA(γ+μA)+κNB,A(t)<e(γ+μA)tA(0)+(κNB+γAmax+μAAmin)e(γ+μA)tt0e(γ+μA)sds,A(t)<e(γ+μA)tA(0)+(κNB+γAmax+μAAminγ+μA)(1e(γ+μA)t),

    where NB:=sup[0,+)N(t). Therefore A(t) is bounded for all t0, which completes the proof of non-negativity and boundedness of the solutions to (2.2).

    The equilibria of (2.2) solve the system

    (γ+μA)A(γAmax+μAAmin)=κN, (3.1)
    0=L((1kpα(A))F(A)(1LηkΔ)δLktα(A)), (3.2)

    and

    N(δNN+μN)=kpα(A)F(A)L(1LηkΔ)+ktα(A)L, (3.3)

    where Δ=0τω(θ)dθ, and X is an equilibrium point with X=(A,L,N). From (3.1) we get that

    N=(γ+μA)A(γAmax+μAAmin)κ, (3.4)

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

    (A1,L1,N1)=(γAmax+μAAminγ+μA,0,0).

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

    (1kpα(A))F(A)(1LΔηk)δLktα(A)=0. (3.5)

    This can be solved for L to give

    L=ηkΔ(1A(δL+ktα(A))βP(AAmin)(1kpα(A))). (3.6)

    Since, for biological reasons, L>0, we obtain

    A(δL+ktα(A))βP(AAmin)(1kpα(A))<1. (3.7)

    Summing (3.2) and (3.3) we get that

    N(δNN+μN)=LF(A)(1LΔηk)δLL, (3.8)

    and substituting F(A)(1LΔηk)=δL+ktα(A)1kpα(A) into (3.8) gives

    N(δNN+μN)=L(δL+ktα(A)1kpα(A)δL). (3.9)

    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

    H(A)=ηkΔ(1A(δL+ktα(A)βp(AAmin)(1kpα(A)))(δL+ktα(A)1kpα(A)δL), (3.10)

    and the left-hand side becomes

    G(A):=δNκ2(((γ+μA)A(γAmax+μAAmin))2+κμNδN((γ+μA)A(γAmax+μAAmin))).

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

    dGdA=δNκ2(2(γ+μA)((γ+μA)A(γAmax+μAAmin))+κμN(γ+μA)δN)=0,

    which yields

    A=(γAmax+μAAmin)γ+μAκμN(γ+μA)2δN(γ+μA)2.

    Therefore, if we define

    A2=2δN(γAmax+μAAmin)κμN2δN(γ+μA),

    then

    G(A2)=μ2N4δ2N,

    which means that minA>A1G(A)<0. As A2<A1 and G(A1)=0 we have that G(A)>0 for all A>A1, 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

    limAAminH(A)+andlimAA+minH(A),

    and that

    limA+H(A)0.

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

    H(A1)=ηkΔδN(1A1(δL+ktα(A1))βP(A1Amin)(1kpα(A1)))(δL+ktα(A1)1kpα(A1)δL),

    for H(A1) 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

    R0=A1(δL+ktα(A1))βP(A1Amin)(1kpα(A1)),

    we require that

    R01.

    The following lemma can be proved easily,

    Lemma 3.1. If A>A11a and R01 then condition (3.7) holds.

    Proof. If A>A1 and A11a then we obtain that α(A1)>α(A) and since the function h(A)=AAAmin is strictly decreasing for all AR+ we obtain that

    A(δL+ktα(A))βP(AAmin)(1kpα(A))<A1(δL+ktα(A1))βP(A1Amin)(1kpα(A1))1.

    Using Lemma 3.1 we can state the following:

    Theorem 3.2. System (2.2) always admits the tumour-free equilibrium

    E1=(A1,L1,N1)=(γAmax+μAAminγ+μA,0,0),

    and if the conditions

    R01, (3.11)

    and

    A11a, (3.12)

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

    L=ηkΔ(1A(δL+ktα(A))βP(AAmin)(1kpα(A))),
    N=(γ+μA)A(γAmax+μAAmin)κ,

    and where A(>A1) solves the equation

    δNκ2(((γ+μA)A(γAmax+μAAmin))2+κμNδN((γ+μA)A(γAmax+μAAmin)))=ηkΔ(1A(δL+ktα(A)βp(AAmin)(1kpα(A)))(δL+ktα(A)1kpα(A)δL).

    Given the generic equilibrium ˉE=(ˉA,ˉL,ˉN) and by using the following change of variables x1=AˉA,x2=LˉL and x3=NˉN we can rewrite system (2.2) as follows

    dx1dt=(γ+μA)x1+κx3+γ(AmaxˉA)μA(ˉAAmin)+κˉN,dx2dt=(1kpα(x1+ˉA))F(x1+ˉA)(x2+ˉL)(11ηk0τ(x2θ+ˉL)ω(θ)dθ)(δL+ktα(x1+ˉA))(x2+ˉL),dx3dt=kpα(x1+ˉA)F(x1+ˉA)(x2+ˉL)(11ηk0τ(x2θ+ˉL)ω(θ)dθ)+ktα(x1+ˉA)(x2+ˉL)δN(x3+ˉN)2μN(x3+ˉN). (3.13)

    The stability analysis of ˉ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:

    dx1dt=(γ+μA)x1+κx3,dx2dt=b1x1+b2x2c10τx2θω(θ)dθ,dx3dt=b3x1+b4x2(2δNˉN+μN)x3c20τx2θω(θ)dθ, (3.14)

    where

    b1=kpα(ˉA)F(ˉA)ˉL(1ˉLΔηk)+F(ˉA)(1kpα(ˉA))ˉL(1ˉLΔηk)ktα(ˉA)ˉL,b2=(1kpα(ˉA))F(ˉA)(1ˉLΔηk)(δL+ktα(ˉA)),b3=kpα(ˉA)F(ˉA)ˉL(1ˉLΔηk)+ktα(ˉA)ˉL+kpα(ˉA)F(ˉA)ˉL(1ˉLΔηk),b4=kpα(ˉA)F(ˉA)(1ˉLΔηk)+ktα(ˉA),

    and

    c1=1ηk(1kpα(ˉA))F(ˉA)ˉL,c2=1ηkkpα(ˉA)F(ˉA)ˉL.

    Remark 3.1. We can represent system (3.14) as a sum of two terms as follows:

    ˙x_=Bx_+0τC(θ)x_θdθ

    where

    B=((γ+μA)0κb1b20b3b4(2δNˉN+μN))

    and

    C=(0000c100c20).

    The characteristic equation is then given by

    det[λIB0τCω(θ)eλθdθ]=0,

    which can be written as:

    ||λ+γ+μA0κb1λb2+c1Φ(λ)0b3b4+c2Φ(λ)λ+2δNˉN+μN||=0, (3.15)

    where

    Φ(λ)=0τω(θ)eλθdθ.

    At E1=(A1,0,0) equation (3.15) yields

    (λ+γ+μA)(λ(1kpα(A1)F(A1)+δL+ktα(A1))(λ+μN)=0,

    which admits the two negative roots λ1=(γ+μA) and λ2=μN, and

    λ3=(1kpα(A1))F(A1)δLktα(A1).

    Therefore E1 is locally asymptotically stable if

    δL+ktα(A1)>(1kpα(A1))F(A1),

    that is,

    A1(δL+ktα(A1))βP(A1Amin)(1kpα(A1))>1.

    This result is summarised in the next lemma.

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

    A1(δL+ktα(A1))βP(A1Amin)(1kpα(A1))>1

    and is unstable otherwise.

    Remark 3.2. We can observe that the condition for the instability of E1 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 b2=0, and from (3.3) we get that

    α(A)(kpF(A)(1LΔηk)+kt)=N(δNN+μN),

    and therefore

    (kpF(A)(1LΔηk)+kt)>0.

    Using this we can rewrite b4=N(δNN+μN)>0. Then the associated characteristic equation is

    0=λ3+λ2(γ+2δNN+μA+μN)+λ(2γδNN+2μAδNN+γμN+μAμNκb3)κb1b4+Φ(λ)[c1λ2+λc1(γ+μA+2δNN+μN)+c1(2γδNNκb3+2μAδNN+γμN+μAμN)+κb1c2]. (3.16)

    For τ=0, this becomes

    0=λ3+λ2(γ+2δNN+μA+μN)+λ(2γδNN+2μAδNN+γμN+μaμNκb3)κb1b4. (3.17)

    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 α(A)<0 and L<ηkΔ which imply b1>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 τ=0, and we thus have the following theorem.

    Theorem 3.4. If τ=0 then whenever E exists, it is unstable.

    To understand how the stability properties of E change with increasing τ we need to specify the distribution function ω(θ). We will consider two delay distributions, the uniform and the exponential distributions, [32,33]. The uniform distribution uses the kernel

    ω(θ)={1τifτθ0,0otherwise, (3.18)

    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 τ days, and that an individual cell will mature in t days, with t[0,τ]. 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

    ω(θ)=ζeψθ. (3.19)

    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 τ0

    Φ(λ)={1λ=0,1eτλτλλ0.

    Then for λ0 and τ0 the characteristic equation (3.16) becomes

    τλ4+λ3(γ+2δNN+μA+μN)τ+λ2(c1+2τδN(γ+μA)N+μNτ(γ+μA)κτb3)+λ(c1(γ+2δNN+μA+μN)κτb1b4)+c1(2δN(γ+μA)N+μN(γ+μA)κb3)+κb1c2eτλ(c1λ2+c1λ(γ+2δNN+μA+μN)+c1(2δN(γ+μA)N+μN(γ+μA)κb3)+κb1c2)=0. (3.20)

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

    D(λ,τ)=Pn(λ,τ)+eτλQm(λ,τ)

    where Pn(λ,τ)=nk=0pk(τ)λk and Qm(λ,τ)=mk=0qk(τ)λk, with n,mN,n>m and pk(τ) and qk(τ) are continuous and differentiable functions of τ. 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 λ=iω. We separate P4(iω,τ) and Q2(iω,τ) into real part

    PR(iω,τ)=τω4ω2(c1+τz1)+z3,QR(iω,τ)=z3+c1ω2,

    and imaginary part

    PI(iω,τ)=ω(c1z2κτb1b4)τz2ω3,QI(iω,τ)=c1z2ω,

    where

    z1=(2δNN+μN)(γ+μA)κb3,z2=γ+2δNN+μA+μN,z3=c1z1+κb1c2.

    We define the function F(ω,τ) as:

    F(ω,τ)=|P4(iω,τ)|2|Q2(iω,τ)|2,

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

    F(ω,τ)=τω6+ω4(τz222(c1+τz1))+ω4(2c1z1+τz21+2z32c1z2(z2κτb1b4))2z1z32b1b4c1z2κ+b21b24κ2τ=0. (3.21)

    Note that ω=0 cannot be considered as a solution to F(ω,τ) for λ0.

    Lemma 3.5. If A>1a+Amin and τ<τ(=2(b1b4c1κz2+z1z3)b21b24κ2) in (3.21), then F(ω,τ) has at least one positive real root ω.

    Proof. Taking A>1a+Amin implies that b3<0 and therefore z1>0, which ensures that τ>0. Since F(ω,τ) is a continuous function of ω over the interval [0,+), limωF(ω,τ)=+, and F(0,τ)=2z1z32b1b4c1z2κ+b21b24κ2τ<0 when τ<τ, then there exists at least one positive solution to the equation F(ω,τ)=0.

    Lemma 3.5 states that if τ(0,τ)(=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 τc at which the bifurcation occurs. For this, let us substitute λ=iω in D(λ,τ)=0. In this way we obtain:

    PR(iω,τ)+iPI(iω,τ)+eiωτ(QR(iω,τ)+iQI(iω,τ))=0.

    By considering real and imaginary parts we obtain:

    τω4ω2(c1+τz1)+z3=cos(τω)(z3c1ω2)+c1ωz2sin(τω),ω(c1z2κτb1b4)τz2ω3=c1z2ωcos(τω)sin(τω)(z3c1ω2),

    which gives

    sin(τω)=τω(z2(c1z1+z3)ω2+b1b4κ(z3c1ω2))z232c1z3ω2+c21ω2(z22+ω2),cos(τω)=c1z2ω2(c1z2b1b4κτz2τω2)+(z3c1ω2)(z3c1ω2z1τω2+τω4)z232c1z3ω2+c21ω2(z22+ω2).

    Dividing cos(τω) by sin(τω) we obtain that

    τω(τ)=θ(τ)+kπ,kN,τI

    where

    θ(τ)=ATanifsin(τω)>0,cos(τω))>0;θ(τ)=π/2ifsin(τω)=1,cos(τω))=0;θ(τ)=π+ATanifcos(τω))<0;θ(τ)=3π/2ifsin(τω)=1,cos(τω))=0;θ(τ)=2π+ATanifsin(τω)<0,cos(τω))>0.

    with

    ATan=arctanc1z2ω2(c1z2b1b4κτz2τω2)+(z3c1ω2)(z3c1ω2z1τω2+τω4)τω(z2(c1z1+z3)ω2+b1b4κ(z3c1ω2)),

    from which we obtain the map

    τk(τ)=θ(τ)+kπω(τ),kN,τI,

    where ω(τ) is a positive real solution to (3.21). We now define the functions

    Sk(τ)=ττk(τ),kN, (3.22)

    which are continuous and differentiable in τ 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 ω(τ) is a positive real root of (3.21) defined for τI and that τc is a root of (3.22) for some nN. 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 λ±(τc)=±iω(τc), with

    sign{dRe(λ)dτ|λ=iω(τc)}=sign{Fω(ω(τc),τc)}sign{dSn(τ)dτ|τ=τc}>0,

    and from right to left otherwise, where λ(τ) is a solution to the characteristic equation (3.20) for τ in a neighbourhood of τ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 ω(θ) be defined as in (3.18). If (i) A>1a+Amin, (ii) (3.21) has at least one positive real root ω(τ) defined for τI, and (iii) at some τcI, Sk(τc)=0 for some kN and Fω(ω(τc),τc)0, then the tumour-present equilibrium undergoes a Hopf bifurcation at a critical value of the time delay τ=τc if sign{Fω(ω(τc),τc)}sign{dSn(τ)dτ|τ=τc}>0.

    The results stated above, together with Theorem 3.6, will be used to calculate the values τ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.

    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.

    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 ϕ(θ)=(ϕ1,ϕ2,ϕ3,ϕ4)T, such that ϕi(θ)0 for τθ<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 ϕ(θ)=[10,3.1,0,U0]T with initial conditions [1,3.1,0,U0]T [20], where for the uniform distribution U0=3.1, and for the exponential distribution U0=3.1ζψ(1eψτ). 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 κ 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 μN (μ 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.

    Table 1.  Parameter values taken from [20] or amended from the first model.
    Parameter Definition Value Units
    γ External androgen production rate 0.013 day1
    Amax Maximum androgen concentration of the system 6 %
    μA Androgen depletion rate 0.08 day1
    Amin Minimum androgen concentration for tumour growth 0.1 %
    κ Secretion rate from N cells 0.009 day1
    δL L cell apoptosis rate 0.013 day1
    kt Maximum transdifferentiation rate for mature L 0.52 day1
    kp Maximum transdifferentiation rate during growth 0.41 day1
    r Gradient of the differentiation increase 3.67 -
    a Inverse of the maximum differentiation rate 1.5 -
    βP Maximum proliferation rate of L 1.4 day1
    ηk Carrying capacity of the environment 3 ×106 cells/l
    δN Intraspecifc competition death rate for N 0.013 day1
    μN Apoptosis death rate for N 0.08 day1
    τ History of the distributed delay 1.42 day
    ψ Scale parameter for the exponential distribution 1 -
    ζ Shape parameter for the exponential distribution 1.318 -
    t Numerical simulation time span 200 day

     | Show Table
    DownLoad: CSV

    Remark 4.1. As shown in Figure 1, by using the parameter values in Table 1 we obtain that the function S0(τ) has two roots. The first root, (τ1c0) is such that sign{dRe(λ)dτ|λ=iω(τ1c)}<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 (τ2c7.6) is such that sign{dRe(λ)dτ|λ=iω(τ2c)}>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 Sk(τ),k>0, suggesting that the only stability switch is the one observed with S0(τ).

    Figure 1.  S0(τ) (solid line) and S1(τ) (dashed line) obtained with parameter values defined in Table 1.

    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 E2, which is to be expected as in this case A11a and condition (3.7) is satisfied, which imply the existence of a (stable) tumour-present equilibrium (Theorem 3.2).

    Figure 2.  Plot of system (A.1) and (A.3) with parameter values as in Table 1 in (a), and with a=1 in (b).

    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 E1. As proved in Theorem 3.6, the stability of the tumour-present equilibrium is affected by τ. In Figure 3 different τ values are considered, τ=0.01,10,20, to demonstrate the effect on the system dynamics.

    Figure 3.  Plot of system (A.1) with parameter values as in Table 1, and different τ values illustrating the stability switch of the tumour-present equilibrium in the case of uniform delay distribution.

    Decreasing τ (Figures 3a) does not affect the solutions of the system (A.1) with uniform distribution. Increasing τ to 10 days introduces oscillatory behaviour to system (A.1), as seen in Figure 3b. Increasing further to τ=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×104. In the case of the exponential distribution, the dynamics of the solutions does not change by altering τ 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.

    Table 2.  Parameter ranges used for searching stability switch in the case of the exponential distribution.
    Parameter Range Parameter Range
    γ [0,0.1] δN [0,0.1]
    Amax [3,20] μN [0,0.1]
    μA [0,0.1] r [3,4.4]
    Amin [0.001,1] a [1,2]
    κ [0.0001,0.01] βP [0,2]
    δL [0,0.1] ηk [2,12]
    kt [0,1] τ [0.01,20]
    kp [0,0.61] ψ [0.005,1]
    ζ [0.052,100.5] t 500

     | Show Table
    DownLoad: CSV

    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.

    Figure 4.  Plots showing different oscillatory behaviours in system (3.19) with parameter values for (a) (or (b) respectively) given by γ=0.04(0.05), δN=0.07(0.005), Amax=4.3(19), μN=0.03(0.06), μA=0.01(0.002), r=4.3(3.7), Amin=1(0.7), a=1.6(1.5), κ=0.003(0.002), βP=1.5(1.5), δL=0.06(0.007), ηk=11.7(4.7), kt=0.2(0.8), τ=15(4.8), kp=0.3(0.4), ψ=0.055(0.12) and ζ=0.097(0.27).

    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 γ, which is smaller in Table 1, and κ, 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 ψ=0.055, and keeping all other parameters but the delay equal to the baseline values (Table 1), a Hopf bifurcation occurs when τ11.8 (Not shown).

    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 ηk and the maximum proliferation rate βp, and the transdifferentiation rates kt and kp. The stability properties of the tumour-free equilibrium E1 were explored fully through a local stability analysis, and the delay τ 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, Amax and Amin, and the external androgen production rate, γ, allowed us to explore further the effect of androgen concentration on the existence of the equilibria. As the stability of E1 is reliant on A1, 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 τ is a bifurcation parameter, and the numerical simulations confirmed this finding and showed that the system undergoes a supercritical Hopf bifurcation as the delay τ increases. Also, we observed that as τ 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 τ is a bifurcation parameter for the stability of the tumour-present equilibrium only when changed together with the scale parameter ψ. Exploring the parameter ranges further, we were able to find examples of oscillatory behaviour with the exponential delay distribution, with indication that also γ and κ 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, α(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 α(A) while retaining biological plausibility in the process.

    The authors declare there is no conflicts of interest.

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

    U(t)=1τttτL(s)ds,

    where s=t+θ. Finding the derivative with respect to t gives

    dUdt=1τ(L(t)L(tτ))

    Hence, system (2.2) becomes

    dAdt=γ(AmaxA)μA(AAmin)+κN,dLdt=(1kpα(A))F(A)L(1Uηk)δLLktα(A)L,dNdt=kpα(A)F(A)L(1Uηk)+ktα(A)LδNN2μNN,dUdt=1τ(L(t)L(tτ)), (A.1)

    with the following initial conditions for the original state variables

    {A(t),L(t),N(t)}={ϕ1(t),ϕ2(t),ϕ3(t)}fort[τ,0], (A.2)

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

    U(0)=1τ0τϕ2(s)ds.

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

    U(t)=ttτζeψ(st)L(s)ds,

    where s=t+θ. Therefore,

    dUdt=ζLζeψτL(tτ)ψU.

    Then system (2.2) becomes

    dAdt=γ(AmaxA)μA(AAmin)+κN,dLdt=(1kpα(A))F(A)L(1Uηk)δLLktα(A)L,dNdt=kpα(A)F(A)L(1Uηk)+ktα(A)LδNN2μNN,dUdt=ζLζeψτL(tτ)ψU. (A.3)

    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

    U(0)=0τζeψsϕ2(s)ds.

    Since distribution kernels can be assumed, without loss of generality, to be positive-definite and normalised to unity [37], we take Δ=1. We then find the following relationship between ζ,ψ and τ.

    0τζeψθdθ=1,

    thus,

    ζψ(1eτψ)=1,

    hence

    ζ=ψ(1eτψ).

    This relationship is used in 4.1 to define the parameter ζ in the simulations.



    [1] 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 (2012), 1133–1145. doi: 10.1002/ijc.27711. doi: 10.1002/ijc.27711
    [2] W. H. Organisation, Cancer Today, available from: https://gco.iarc.fr/today, Last Accessed: 2020-06-05.
    [3] P. J. Hensley, N. Kyprianou, Modeling prostate cancer in mice: limitations and opportunities, J. Androl., 33 (2012), 133–144. doi: 10.2164/jandrol.111.013987. doi: 10.2164/jandrol.111.013987
    [4] J. Horoszewicz, S. Leong, T. Ming-Chu, Z. Wajsman, M. Friedman, L. Papsidero, et al., The LNCaP cell line - a new model for studies on human prostatic carcinoma, Prog. Clin. Biol. Res., 37 (1980), 115–132.
    [5] T. Phan, S. Crook, A. Bryce, C. Maley, E. Kostelich, Y. Kuang, Mathematical modeling of prostate cancer and clinical application, Appl. Sci., 10 (2020), 2721. doi: 10.3390/app10082721. doi: 10.3390/app10082721
    [6] K. Swanson, L. True, D. Lin, K. Buhler, R. Vessella, J. 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 (2001), 2513–2522. doi: 10.1016/S0002-9440(10)64691-3. doi: 10.1016/S0002-9440(10)64691-3
    [7] R. 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 (2002), 1692–1698. doi: 10.1002/cncr.10443. doi: 10.1002/cncr.10443
    [8] R. Vollmer, P. Humphrey, Tumor volume in prostate cancer and serum prostate-specific antigen: Analysis from a kinetic viewpoint, Am. J. Pathol., 119 (2003), 80–89. doi: 10.1309/UNAQ-JTFP-B1RQ-BQD4. doi: 10.1309/UNAQ-JTFP-B1RQ-BQD4
    [9] Y. Kuang, J. Nagy, J. Elser, Biological stoichiometry of tumor dynamics: mathematical models and analysis, Disc. Cont. Dyn. Sys. B, 4 (2004), 221–240. doi: 10.3934/dcdsb.2004.4.221. doi: 10.3934/dcdsb.2004.4.221
    [10] C. Heinlein, C. Chang, Androgen receptor in prostate cancer, Endocr. Rev., 25 (2004), 276–308. doi: 10.1210/er.2002-0032. doi: 10.1210/er.2002-0032
    [11] P. Koivisto, M. Kolmer, T. Visakorpi, O. Kallioniemi, Androgen receptor gene and hormonal therapy failure of prostate cancer, Am. J. Pathol., 152 (1998), 1–9.
    [12] R. Rittmaster, A. Manning, A. Wright, L. Thomas, S. Whitefield, R. Norman, et al., Evidence for atrophy and apoptosis in the ventral prostate of rats given the 5 alpha-reductase inhibitor finasteride, Endocrinology, 136 (1995), 741–748. doi: 10.1210/endo.136.2.7835306. doi: 10.1210/endo.136.2.7835306
    [13] T. Jackson, A mathematical investigation of the multiple pathways to recurrent prostate cancer: Comparison with experimental data, Neoplasia, 6 (2004), 697–704. doi: 10.1593/neo.04259. doi: 10.1593/neo.04259
    [14] T. Jackson, A mathematical model of prostate tumor growth and androgen-independent relapse, Disc. Cont. Dyn. Sys. B, 4 (2004), 187–201. doi: 10.3934/dcdsb.2004.4.187. doi: 10.3934/dcdsb.2004.4.187
    [15] A. Ideta, G. Tanaka, T. Takeuchi, K. Aihara, A mathematical model of intermittent androgen suppression for prostate cancer, J. Nonlinear Sci., 18 (2008), 593–614. doi: 10.1007/s00332-008-9031-0. doi: 10.1007/s00332-008-9031-0
    [16] S. Eikenberry, J. Nagy, Y. Kuang, The evolutionary impact of androgen levels on prostate cancer in a multi-scale mathematical model, Biol. Direct, 5. doi: 10.1186/1745-6150-5-24.
    [17] S. Terry, H. Beltran, The many faces of neuroeondocrine differentiation in prostate cancer progression, Front. Oncol., 4 (2014), 1–9. doi: 10.3389/fonc.2014.00060. doi: 10.3389/fonc.2014.00060
    [18] V. Perrot, Neuroendocrine differentiation in the progression of prostate cancer: an update to recent developments, Open J. Urol., 2 (2012), 173–182. doi: 10.4236/oju.2012.223032. doi: 10.4236/oju.2012.223032
    [19] C. Shen, Z. Burke, D. Tosh, Transdifferentiation, metaplasia and tissue regeneration, Organogenesis, 1 (2004), 36–44. doi: 10.4161/org.1.2.1409. doi: 10.4161/org.1.2.1409
    [20] 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 (2015), 2975–2986. doi: 10.1158/0008-5472.CAN-14-3830. doi: 10.1158/0008-5472.CAN-14-3830
    [21] 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 (2010), 2913–2929. doi: 10.1016/j.nonrwa.2009.11.001. doi: 10.1016/j.nonrwa.2009.11.001
    [22] M. Adimy, F. Crauste, S. Ruan, Modelling hematopoiesis mediated by growth factors with applications to periodic hematological diseases, Bull. Math. Biol., 68 (2006), 2321–2351. doi: 10.1007/s11538-006-9121-9. doi: 10.1007/s11538-006-9121-9
    [23] J. Morken, A. Packer, R. Everett, J. Nagy, Y. Kuang, Mechanisms of resistance to intermittent androgen deprivation in patients with prostate cancer identified by a novel computational method, Cancer Res., 74 (2014), 3673–3683. doi: 10.1158/0008-5472.CAN-13-3162. doi: 10.1158/0008-5472.CAN-13-3162
    [24] L. Turner, A. Burbanks, M. Cerasuolo, Mathematical insights into neuroendocrine transdifferentiation of human prostate cancer cells, Nonlinear Anal. Model. Control, 26 (2021), 884–913. doi: 10.15388/namc.2021.26.24441. doi: 10.15388/namc.2021.26.24441
    [25] G. Hutchinson, Circular causal systems in ecology., Ann. N. Y. Acad. Sci., 50 (1948), 221–246. doi: 10.1111/j.1749-6632.1948.tb39854.x. doi: 10.1111/j.1749-6632.1948.tb39854.x
    [26] T. Cassidy, A. R. Humphries, A mathematical model of viral oncology as an immuno-oncology instigator, Math. Med. Biol., 37 (2020), 117–151. doi: 10.1093/imammb/dqz008. doi: 10.1093/imammb/dqz008
    [27] F. Scudo, J. Ziegler, The golden age of theoretical ecology, 1923-1940: A collection of the works of V.Volterra, V.A. Kostitzin, A.J. Lotka, and A.N. Kolmogoroff, vol. 22, Springer, 1978.
    [28] J. Baez, Y. Kuang, Mathematical models of androgen resistance in prostate cancer patients under intermittent androgen suppression therapy, Appl. Sci., 6 (2016), 352. doi: 10.3390/app6110352. doi: 10.3390/app6110352
    [29] M. Droop, Vitamin B12 and marine ecology, IV: The kinetics of uptake, growth and inhibition in Monochrysis lutheri, J. Mar. Biol. Assoc, UK, 48 (1968), 689–733. doi: 10.1017/S0025315400019238. doi: 10.1017/S0025315400019238
    [30] B. Buonomo, M. Cerasuolo, The effect of time delay in plant-pathogen interactions with host demography, Math. Biosci. Eng., 12 (2015), 473–490. doi: 10.3934/mbe.2015.12.473. doi: 10.3934/mbe.2015.12.473
    [31] T. Faria, E. Liz, Boundedness and asymptotic stability for delayed equations of logistic type., Proc. Math. Roy. Soc. Edinb., 133 (2003), 1057–1073. doi: 10.1017/S030821050000281X. doi: 10.1017/S030821050000281X
    [32] A. Ahmadian, M. Bin Suleiman, F. Ismail, Numerical simulation of tumor development stages using artificial neural network., Trends Appl. Sci. Res., 7 (2012), 132–141. doi: 10.3923/tasr.2012.132.141. doi: 10.3923/tasr.2012.132.141
    [33] M. Piotrowska, M. Bodnar, Influence of distributed delays on the dynamics of a generalized immune system cancerous cells interactions model, Commun. Nonlinear Sci. Numer. Simulat., 54 (2018), 38. doi: 10.1016/j.cnsns.2017.06.003. doi: 10.1016/j.cnsns.2017.06.003
    [34] E. Beretta, Y. Kuang, Geometric stability switch criteria in delay differential systems with delay dependent parameters, SIAM J. Math. Anal., 33 (2002), 1144–1165. doi: 10.1137/S0036141000376086. doi: 10.1137/S0036141000376086
    [35] H. Smith, An Introduction to Delay Differential Equations with Applications to the Life Sciences, Springer, New York, 2010. doi: 10.1007/978-1-4419-7646-8.
    [36] P.-L. Clermont, X. Ci, H. Pandha, Y. Wang, F. Crea, Treatment-emergent neuroendocrine prostate cancer: molecularly driven clinical guidelines, Int. J. Endocr. Oncol., 6 (2019), IJE20. doi: 10.2217/ije-2019-0008. doi: 10.2217/ije-2019-0008
    [37] B. Rahman, K. Blyuss, Y. Kyrychko, Dynamics of neural systems with discrete and distributed time delays, SIAM J. Appl. Dyn., 14 (2015), 2069–2095. doi: 10.1137/15M1006398. doi: 10.1137/15M1006398
  • This article has been cited by:

    1. Andrew Burbanks, Marianna Cerasuolo, Roberto Ronca, Leo Turner, A hybrid spatiotemporal model of PCa dynamics and insights into optimal therapeutic strategies, 2023, 355, 00255564, 108940, 10.1016/j.mbs.2022.108940
    2. William Meade, Allison Weber, Tin Phan, Emily Hampston, Laura Figueroa Resa, John Nagy, Yang Kuang, High Accuracy Indicators of Androgen Suppression Therapy Failure for Prostate Cancer—A Modeling Study, 2022, 14, 2072-6694, 4033, 10.3390/cancers14164033
    3. Zheng Fu, Benzhong Jia, Advances in the role of heat shock protein 90 in prostate cancer, 2022, 54, 0303-4569, 10.1111/and.14376
    4. Tin Phan, Allison Weber, Alan H. Bryce, Yang Kuang, The prognostic value of androgen to PSA ratio in predictive modeling of prostate cancer, 2023, 176, 03069877, 111084, 10.1016/j.mehy.2023.111084
  • Reader Comments
  • © 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2924) PDF downloads(73) Cited by(4)

Figures and Tables

Figures(4)  /  Tables(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog