# Modeling the Effects of Chemotherapeutic Dose Response on a Stochastic Tumor-Immune Model of Prostate Cancer with Androgen Deprivation Therapy

Authors: Lin Chen; Jin Yang
Journal: Discrete Dynamics in Nature and Society
 (2023)

---

## Abstract

The periodical application of androgen deprivation therapy, immunotherapy, or chemotherapy is an effective method for cancer treatment, but few studies combine them. To explore such comprehensive treatment mechanisms, this paper establishes a pulsed stochastic hybrid dynamics model considering tumor antigenicity and density-dependent mortality. In addition to analyzing the basic properties of solutions such as the tumor-free periodic solution and global attraction of the model, the threshold conditions for the persistence and extinction of prostate cancer cells and effector cells are obtained by using stochastic differential equation theory. Besides, sufficient conditions for the existence of stationary distribution of the system are established. The results reveal that comprehensive therapy or white noise can determine tumor dynamics and suggest that the treatment of prostate cancer should be individualized according to the state of tumor development. Finally, biological significance is discussed and conclusions are given.


---

## Body

## 1. Introduction

Prostate cancer is the second leading cause of cancer death and the most common type of cancer in American men [1, 2]. Prostate specific antigen (PSA) is a protein secreted by prostate cells. It is the main index to detect the presence of malignant tumors [3–5]. Prostate cancer is a hormone-dependent cancer [6]. Prostate cancer cells include androgen dependent (AD) and androgen independent (AI) cells during evolution. Surgery, chemotherapy, and radiotherapy are traditional methods, but these treatments often cannot eliminate cancer cells completely and may have side effects [7]. To this end, androgen deprivation therapy (ADT) and immunotherapy came into being. ADT has become the main treatment for prostate cancer because it affects the proliferation rate and mortality of tumor cells [8, 9]. Besides, ADT includes intermittent androgen deprivation (IAD) and continuous androgen deprivation (CAD) therapy [3]. However, ADT is often accompanied by the emergence of drug-resistant cells, then, an additional treatment is needed to introduce. Since immunotherapy has become the most effective method for the treatment of prostate cancer by stimulating a massive immune response in the target tumor and enhancing the immune system of patients [1, 7, 10–12]. It is reasonable and effective to add immunotherapy to ADT to treat prostate cancer. Moreover, CAD therapy is one of the endocrine therapies for advanced prostate cancer, and there is no clear medical evidence that IAD is better than CAD therapy [3, 13].Mathematical modeling and analysis of ADT and immunotherapy have devoted to the development of mechanisms for tumor progression [1, 8, 9, 14–20]. Ideta et al. established a deterministic model of IAD therapy, illustrating that how the net growth rate of AI cells affects tumor growth and recurrence [16]. Jain and Friedman proposed a mathematical model to simulate prostate cancer response to ADT [17]. In addition, the interaction between tumors and immune cells has been studied by many mathematical models [7, 10–12, 21–24]. Since many studies have shown that single chemotherapy or immunotherapy is far less effective than comprehensive therapy, the therapeutic effect of single chemotherapy is also very limited [1, 10, 25, 26]. Yang et al. pointed out that tumor cells can be eradicated or controlled with comprehensive therapy [11]. Hence, the mathematical model of the combination of ADT and immunotherapy has been studied. Portz and Kuang then developed a model of advanced prostate cancer to test the efficacy of immunotherapy combined with ADT. Their results suggest that comprehensive treatment can stabilize the disease [1]. Moreover, Rutter and Kuang confirmed that IAD therapy can delay the occurrence of drug resistance and improve the quality of life of patients [20].Nevertheless, most of these models are all deterministic [1, 11, 16, 17, 20], but the growth of cancer cells is easily affected by temperature, nutrition, radiation, and other environmental interferences [3, 7, 10, 27–33]. Naturally, the deterministic system of Ideta et al. was developed into the stochastic system of Tanaka and collaborators [16, 27]. In 2018, Hening and Nguyen studied the coexistence and extinction of stochastic Kolmogorov systems and gave sharp conditions under which populations converge exponentially to their unique stationary distributions [32]. Benaim investigated the stochastic persistence and extinction of populations, and the results greatly generalized the results of various stochastic models of population dynamics given by stochastic differential equations or pure jump processes [33]. Zazoua and Wang introduced the stochastic model of prostate cancer and showed the stability of the system that can be determined by white noise [3]. Furthermore, clinical data have shown that immunotherapy can promote chemotherapy [7, 25, 26, 34]. Thereafter, Yang established random tumor-immune models combining immunotherapy and chemotherapy, and observed that all dynamic processes of tumors can be realized by changing environmental noise [7, 10]. In 2021, Chen et al. proposed a mathematical model describing the threshold dynamics of tumors under random interference of tumor antigenicity and white noise [35]. Yang et al. studied stochastic impulsive tumor models combining ADT and immunotherapy to analyze the elimination and persistence of tumor cells [23, 24]. Moreover, Rihan and Alsakaji proposed a stochastic epidemic model with time delay. It is proved that the stochastic delay model is consistent with the physical sensitivity and volatility of the actual observation [29]. Rihan and Rajivganthi introduced a stochastic tumor-immune model with random noise in epidemiology and immunology, and concluded that white noise is a key factor in treating infectious diseases [28]. In 2022, Alsakaji et al. expanded the SEIR epidemic model by combining stochastic perturbation with time delay, established a mathematical model considering vaccination, time delay and random noise, and studied the dynamics of COVID-19 in the UAE [30]. Furthermore, the dynamics of delayed differential model of tumor-immune system with random noise was investigated by Rihan and Alsakaji. Their results showed that in some cases, random noise can completely inhibit the growth of tumors [31]. In addition, the idea of randomness is also widely used in numerical calculation, which is of great significance for the development of stochastic computing procedure [36–42].As our knowledge, immunotherapy combined with chemotherapy or ADT has been widely studied, but there are few scholars combining these three therapies to establish models [1, 7, 10, 20, 34]. Actually, it is necessary to consider the influence of tumor antigenicity on its growth, because it always exists in the whole life cycle of the tumors [11]. In addition, some cancer cells will die due to limited living space and insufficient resources and nutrition [14]. For the purpose, considering both antigenicity and density-dependent mortality of tumors, we propose a pulsed stochastic model combining CAD therapy, immunotherapy, and chemotherapy to explore hybrid dynamics and cancer therapy strategies. We mainly answered the following questions: (1) How do the periods, dosages, and frequencies of in pulse comprehensive therapy affect tumor dynamics? (2) Does the system have a unique ergodic stationary distribution? (3) What is the difference between comprehensive treatment and single treatment and what are the medication strategies for clinical medicine?The rest of this paper is arranged as follows. In Section2, a pulsed stochastic model is proposed and some useful basic knowledge of ISDEs (Impulsive stochastic differential equations) is presented. In Section 3, the expression of a tumor-free solution and the global attraction and boundness of the system are showed. In Section 4, threshold conditions for persistence and extinction of prostate cancer cells and effector cells are obtained. In Section 5, the existence of a unique ergodic stationary distribution of the hybrid system is analyzed. The numerical simulations are given out to verify the theoretical results in Section 6. Finally, this paper is concluded in Section 7.


## 2. Preliminaries and Mathematical Model

### 2.1. Model Formation

In 2017, Rutter and Kuang developed a deterministic model of prostate cancer dendritic cell vaccine combined with hormone therapy to explore its global dynamics [20]. To separate the influence of stochastic noise from immune response, Zazoua and Wang did not include the state variables of effector cells and cytokines related to immunotherapy. Inspired by Tanaka et al. [27], they proposed a stochastic mathematical model to explore the evolution of prostate cancer cells [3, 20]. Let A be the concentration of androgen in the blood. Z1 and Z2 are AD and AI cells, respectively. Their model is as follows:(1)dAdt=−γA−a0−γa0u,dZ1=r1A1−Z1+αZ2K−p1+m11−Aa0Z1dt+σ1Z1dB1t,dZ2=r21−βZ1+Z2KZ2+m11−Aa0Z1dt+σ2Z2dB2t.Since the single treatment can easily make prostate cancer cells develop drug resistance, and AD cells eventually develop into AI cells. It is particularly critical to add new treatments to reduce drug resistance [25, 26]. In 1994, Kuznetsov et al. developed an influential deterministic tumor-immune model [43]. Moreover, some experimental and clinical studies have shown that pulsed immunotherapy combined with chemotherapy is more practical than single therapy for treating cancer [44, 45]. Based on the abovementioned considerations, Yang incorporated the white noise into the model of Kuznetsov et al. [43], and studied a stochastic model of chemotherapy combined with immunotherapy, which showed that the larger white noise can lead to the extinction of tumors and the combination therapy could avoid the shortcomings of a single therapy [34, 43–45]. Then, taking immunotherapy as an auxiliary treatment to prevent drug resistance is reasonable. Let Xt and Yt be the prostate cancer cells and effector cells, respectively. H is the concentration of a chemotherapeutic drug at time t. The system is written as follows [10]:(2)dX=r1−ηXX−aXY−k1HtXdt+δ1XdB1t,dY=bXY1+υX−cXY−dY−k2HtYdt+δ2YdB2t,dH=−μHdt,t≠nT,HnT+=HnT+ι,YnT+=1+enTYnT,t=nT.However, few scholars have explored the dynamic behaviors of complex dynamical systems that combine ADT with immunotherapy and chemotherapy. It is very meaningful to study how dosages, periods (or frequencies) of drug use under pulse comprehensive treatment affect tumor dynamics and provide theoretical guidance for clinical treatment. Inspired by CAD in system (1) and immunotherapy and chemotherapy in system (2), we take tumor antigenicity and white noise into consideration to investigate the evolution of prostate cancer cells under comprehensive treatment [3, 7, 10, 11, 14, 35]. Then, our model is described by(3)dX=rA1−XKX−d11−Aa0+d2X−k1HtX−aXYdt+δ1XdB1t,dY=C−dY−k2HtY−cXYdt+δ2YdB2t,dA=−γA−a0−γa0udt,dH=−μHdt,t≠nT,HnT+=HnT+ι,YnT+=1+enTYnT,t=nT,where Biti=1,2 is an independent Brownian motion. The biological meanings of other parameters are the same as systems (1) and (2), and all of these parameters in system (3) are positive. The biological meanings of other parameters from system (1) to system (3) are shown in Table 1.Table 1
Biological meaning and source of parameters.
ParametersBiological meaningSourceγThe clearance and production rate of androgen[3]uThe efficacy of CAD therapy[3]r1The proliferation rate of AD cells[3]r2The proliferation rate of AI cells[3]αThe competition coefficient of AD cells[3]βThe competition coefficient of AI cells[3]p1The mortality of AD cells[3]a0The normal androgen concentration[3]TThe periods of therapy[3]KThe carrying capacity of cancer cells[3]m1The mutation rate from AD to AI cells[3]σ12The intensities of white noise on the AD cells[3]σ22The intensities of white noise on the AI cells[3]1/ηCarrying capacity of the tumor cells[10]dThe mortality of effector cells[10]bThe maximum accumulation rate of tumors[10]υThe steepness of effector cells[10]cThe inactivation rate of effector cells[10]aThe rate at which the effector cells bind to the tumor cells[10]enTThe net growth rate of effector cells stimulated by immunotherapy[10]μThe degradation rate of chemotherapy[10]ιThe dosages at impulsive point seriesnTn=1,2,3,⋯[10]rThe growth rate of cancer cells[10]d1The death rate of prostate cancer cells[3]d2The density-dependent mortality of prostate cancer cells[46]δ12The intensities of white noise on prostate cancer cells[3]δ22The intensities of white noise on effector cells[10]k1The inhibitory rate of chemotherapy drugs on prostate cancer cells[10]k2The inhibitory rate of chemotherapy drugs on effector cells[10]The evolution of androgen dynamics takes less time to reach equilibrium than cancer cells because androgen dynamics is faster than prostate cancer cells [3]. Thus, we let androgen go into equilibrium A∗=a01−u, and let R1=rA∗=ra01−u,U1=d1u+d2, then system (3) becomes(4)dX=R1X1−XK−U1X−k1HtX−aXYdt+δ1XdB1t,dY=C−dY−k2HtY−cXYdt+δ2YdB2t,dH=−μHdt,t≠nT,HnT+=HnT+ι,YnT+=1+enTYnT,t=nT.Next, we first give the biological meaning and source of the parameters involved in this paper, and then, give the basic knowledge for further research.


### 2.2. Preliminaries

LetΩ,F,Ftt≥0,P be a complete probability space (Bit is defined on this), which conforms to the usual conditions and has a filtration Ftt≥0, where Bit is the independent Brownian motion [7, 10, 47]. Then, we introduce some basic knowledge of the paper.Definition 1.
(see [7, 10, 48]). Zt=Pt,Qt)T is a solution of ISDE (4) with initial values P0≥0 and Q0≥0, t∈R+, if(1)
Zt is absolutely continuous on 0,TandnT,n+1T(2)
For anynT, ZnT−=limt⟶nT−ZtandZnT+=limt⟶nT+ZtandZnT=ZnT−hold true(3)
Ztfollows system (4) for allt∈R+/nTand satisfies the impulsive condition at the impulse pointnTDefinition 2.
For the solutionZt=Pt,Qt)T with initial conditions P0∈R+2 and Q0∈R+2, if for any ϵ∈0,1, there is a solution with a positive constant S such that(5)limsupt⟶∞PZt>S<ϵ.
Then, the solutions of ISDE (4) are stochastically ultimately bounded.Definition 3.
(see [7, 49]). LetZt=Pt,Qt)T be a solution of ISDE (4):(1)
Whenlimt⟶+∞Qt=0, Qt is extinctive.(2)
Whenlimt⟶+∞1/t∫0tQsds=0, Qt is nonpersistent in the mean.(3)
Whenlimt⟶+∞sup1/t∫0tQsds>0, Qt is weakly persistent in the mean.(4)
If for eachε∈0,1, there exist β>0 and δ>0 so that(6)liminft⟶+∞PQt≥β≥1−ε,liminft⟶+∞PQt≤δ≥1−ε.then Yt is stochastically persistent.Definition 4.
(see [7]). Assume Z1t=P1t,Q1t and Z2t=P2t,Q2t are any two solutions of system (4) with Pi0>0, Qi0>0, i=1,2, if(7)limt⟶+∞P1t−P2t=0,and(8)limt⟶+∞Q1t−Q2t=0.
Then, ISDE (4) is called globally attractive.Definition 5.
(see [3]). If(9)gt=1t∫0tgsdsfort>0,then(10)g∗t=liminft⟶∞1t∫0tgsdsfort>0,and(11)g∗t=limsupt⟶∞1t∫0tgsdsfort>0.Lemma 1 (see [7, 10]).
Letgt∈CΩ×R+,R+−0.(1)
If there existη0,t1 and η≥0 such that gt satisfies(12)lngt≤ηt−η0∫0tgsds+∑i=1nβiBit,for any t≥t1, βi is a constant, then(13)limsupt⟶+∞1t∫0tgsds≤ηη0.(2)
If there existη0,t1 and η≥0 such that gt satisfies(14)lngt≥ηt−η0∫0tgsds+∑i=1nβiBit,for any t≥t1, then(15)limsupt⟶+∞1t∫0tgsds≥ηη0.Lemma 2 (see [3]).
AssumeF=Ftt≥0 is a real-valued continuous local martingale vanishing at time zero. If(16)limsupt⟶+∞F,Ftt<∞a.s.,then,(17)limsupt⟶+∞Ftt=0.


## 2.1. Model Formation

In 2017, Rutter and Kuang developed a deterministic model of prostate cancer dendritic cell vaccine combined with hormone therapy to explore its global dynamics [20]. To separate the influence of stochastic noise from immune response, Zazoua and Wang did not include the state variables of effector cells and cytokines related to immunotherapy. Inspired by Tanaka et al. [27], they proposed a stochastic mathematical model to explore the evolution of prostate cancer cells [3, 20]. Let A be the concentration of androgen in the blood. Z1 and Z2 are AD and AI cells, respectively. Their model is as follows:(1)dAdt=−γA−a0−γa0u,dZ1=r1A1−Z1+αZ2K−p1+m11−Aa0Z1dt+σ1Z1dB1t,dZ2=r21−βZ1+Z2KZ2+m11−Aa0Z1dt+σ2Z2dB2t.Since the single treatment can easily make prostate cancer cells develop drug resistance, and AD cells eventually develop into AI cells. It is particularly critical to add new treatments to reduce drug resistance [25, 26]. In 1994, Kuznetsov et al. developed an influential deterministic tumor-immune model [43]. Moreover, some experimental and clinical studies have shown that pulsed immunotherapy combined with chemotherapy is more practical than single therapy for treating cancer [44, 45]. Based on the abovementioned considerations, Yang incorporated the white noise into the model of Kuznetsov et al. [43], and studied a stochastic model of chemotherapy combined with immunotherapy, which showed that the larger white noise can lead to the extinction of tumors and the combination therapy could avoid the shortcomings of a single therapy [34, 43–45]. Then, taking immunotherapy as an auxiliary treatment to prevent drug resistance is reasonable. Let Xt and Yt be the prostate cancer cells and effector cells, respectively. H is the concentration of a chemotherapeutic drug at time t. The system is written as follows [10]:(2)dX=r1−ηXX−aXY−k1HtXdt+δ1XdB1t,dY=bXY1+υX−cXY−dY−k2HtYdt+δ2YdB2t,dH=−μHdt,t≠nT,HnT+=HnT+ι,YnT+=1+enTYnT,t=nT.However, few scholars have explored the dynamic behaviors of complex dynamical systems that combine ADT with immunotherapy and chemotherapy. It is very meaningful to study how dosages, periods (or frequencies) of drug use under pulse comprehensive treatment affect tumor dynamics and provide theoretical guidance for clinical treatment. Inspired by CAD in system (1) and immunotherapy and chemotherapy in system (2), we take tumor antigenicity and white noise into consideration to investigate the evolution of prostate cancer cells under comprehensive treatment [3, 7, 10, 11, 14, 35]. Then, our model is described by(3)dX=rA1−XKX−d11−Aa0+d2X−k1HtX−aXYdt+δ1XdB1t,dY=C−dY−k2HtY−cXYdt+δ2YdB2t,dA=−γA−a0−γa0udt,dH=−μHdt,t≠nT,HnT+=HnT+ι,YnT+=1+enTYnT,t=nT,where Biti=1,2 is an independent Brownian motion. The biological meanings of other parameters are the same as systems (1) and (2), and all of these parameters in system (3) are positive. The biological meanings of other parameters from system (1) to system (3) are shown in Table 1.Table 1
Biological meaning and source of parameters.
ParametersBiological meaningSourceγThe clearance and production rate of androgen[3]uThe efficacy of CAD therapy[3]r1The proliferation rate of AD cells[3]r2The proliferation rate of AI cells[3]αThe competition coefficient of AD cells[3]βThe competition coefficient of AI cells[3]p1The mortality of AD cells[3]a0The normal androgen concentration[3]TThe periods of therapy[3]KThe carrying capacity of cancer cells[3]m1The mutation rate from AD to AI cells[3]σ12The intensities of white noise on the AD cells[3]σ22The intensities of white noise on the AI cells[3]1/ηCarrying capacity of the tumor cells[10]dThe mortality of effector cells[10]bThe maximum accumulation rate of tumors[10]υThe steepness of effector cells[10]cThe inactivation rate of effector cells[10]aThe rate at which the effector cells bind to the tumor cells[10]enTThe net growth rate of effector cells stimulated by immunotherapy[10]μThe degradation rate of chemotherapy[10]ιThe dosages at impulsive point seriesnTn=1,2,3,⋯[10]rThe growth rate of cancer cells[10]d1The death rate of prostate cancer cells[3]d2The density-dependent mortality of prostate cancer cells[46]δ12The intensities of white noise on prostate cancer cells[3]δ22The intensities of white noise on effector cells[10]k1The inhibitory rate of chemotherapy drugs on prostate cancer cells[10]k2The inhibitory rate of chemotherapy drugs on effector cells[10]The evolution of androgen dynamics takes less time to reach equilibrium than cancer cells because androgen dynamics is faster than prostate cancer cells [3]. Thus, we let androgen go into equilibrium A∗=a01−u, and let R1=rA∗=ra01−u,U1=d1u+d2, then system (3) becomes(4)dX=R1X1−XK−U1X−k1HtX−aXYdt+δ1XdB1t,dY=C−dY−k2HtY−cXYdt+δ2YdB2t,dH=−μHdt,t≠nT,HnT+=HnT+ι,YnT+=1+enTYnT,t=nT.Next, we first give the biological meaning and source of the parameters involved in this paper, and then, give the basic knowledge for further research.


## 2.2. Preliminaries

LetΩ,F,Ftt≥0,P be a complete probability space (Bit is defined on this), which conforms to the usual conditions and has a filtration Ftt≥0, where Bit is the independent Brownian motion [7, 10, 47]. Then, we introduce some basic knowledge of the paper.Definition 1.
(see [7, 10, 48]). Zt=Pt,Qt)T is a solution of ISDE (4) with initial values P0≥0 and Q0≥0, t∈R+, if(1)
Zt is absolutely continuous on 0,TandnT,n+1T(2)
For anynT, ZnT−=limt⟶nT−ZtandZnT+=limt⟶nT+ZtandZnT=ZnT−hold true(3)
Ztfollows system (4) for allt∈R+/nTand satisfies the impulsive condition at the impulse pointnTDefinition 2.
For the solutionZt=Pt,Qt)T with initial conditions P0∈R+2 and Q0∈R+2, if for any ϵ∈0,1, there is a solution with a positive constant S such that(5)limsupt⟶∞PZt>S<ϵ.
Then, the solutions of ISDE (4) are stochastically ultimately bounded.Definition 3.
(see [7, 49]). LetZt=Pt,Qt)T be a solution of ISDE (4):(1)
Whenlimt⟶+∞Qt=0, Qt is extinctive.(2)
Whenlimt⟶+∞1/t∫0tQsds=0, Qt is nonpersistent in the mean.(3)
Whenlimt⟶+∞sup1/t∫0tQsds>0, Qt is weakly persistent in the mean.(4)
If for eachε∈0,1, there exist β>0 and δ>0 so that(6)liminft⟶+∞PQt≥β≥1−ε,liminft⟶+∞PQt≤δ≥1−ε.then Yt is stochastically persistent.Definition 4.
(see [7]). Assume Z1t=P1t,Q1t and Z2t=P2t,Q2t are any two solutions of system (4) with Pi0>0, Qi0>0, i=1,2, if(7)limt⟶+∞P1t−P2t=0,and(8)limt⟶+∞Q1t−Q2t=0.
Then, ISDE (4) is called globally attractive.Definition 5.
(see [3]). If(9)gt=1t∫0tgsdsfort>0,then(10)g∗t=liminft⟶∞1t∫0tgsdsfort>0,and(11)g∗t=limsupt⟶∞1t∫0tgsdsfort>0.Lemma 1 (see [7, 10]).
Letgt∈CΩ×R+,R+−0.(1)
If there existη0,t1 and η≥0 such that gt satisfies(12)lngt≤ηt−η0∫0tgsds+∑i=1nβiBit,for any t≥t1, βi is a constant, then(13)limsupt⟶+∞1t∫0tgsds≤ηη0.(2)
If there existη0,t1 and η≥0 such that gt satisfies(14)lngt≥ηt−η0∫0tgsds+∑i=1nβiBit,for any t≥t1, then(15)limsupt⟶+∞1t∫0tgsds≥ηη0.Lemma 2 (see [3]).
AssumeF=Ftt≥0 is a real-valued continuous local martingale vanishing at time zero. If(16)limsupt⟶+∞F,Ftt<∞a.s.,then,(17)limsupt⟶+∞Ftt=0.


## 3. Global Positive Solution

### 3.1. Tumor-free Solution

It is obvious that the dynamics of chemotherapeutic drugs can be described by(18)dHt=−μH,t≠nT,HnT+=HnT+ι,t=nT.Let us compute the expression for theT periodic solution HTt as the following:(19)HTt=ιe−μt−nT1−e−μT,where t∈nT,n+1T, HTnT+=ι/1−e−μT.Lemma 3 (see [10, 50, 51]).
A unique positive periodic solution of system (5) can be expressed by HTt, and it holds that limt⟶∞Ht=HTt, then for any ϵ>0 yields(20)HTt−ϵ<Ht<HTt+ϵandlimt⟶∞1t∫0tHTsds=ιμT.For the convenience of writing, we useHt instead of HTt. Suppose that the prostate cancer cells can be eradicated, then let Xt=0, and system (4) becomes(21)dY=C−dY−k2HtYdt+δ2YdB2t,dH=−μH,t≠nT,HnT+=HnT+ι,YnT+=1+enTYnT,t=nT.Since the explicit expression ofHt is given, then(22)dY=C−dY−k2HtYdt+δ2YdB2t,t≠nT,YnT+=1+enTYnT,t=nT.Theorem 1.
There exists a unique global positive solutionYt of system (8) for any initial value Y0+=Y0, where(23)Yt=∏0<nT<t1+enTY0expC−d−k2Ht−δ222t+δ2B2t.Proof.
For anyt∈nT,n+1T, define a Lyapunov function Vt=lnYt. Applying Itô’s formula gives(24)dlnYt=C−d−k2Ht−0.5δ22dt+δ2dB2t.
We integrate the above equation fromnT to t, then(25)lnYt−lnYnT=C−d−k2Ht−0.5δ22t−nT+δ2B2t−B2nT.
Therefore,(26)Yt=YnTexpC−d−k2Ht−0.5δ22t−nT+δ2B2t−B2nT.
Whent=nT+, after a single immunotherapy, we have(27)Yt=∏0<nT<t1+enTYnTexpC−d−k2Ht−δ222t−nT+δ2B2t−B2nT.
An application of mathematical induction yields(28)Yt=∏0<nT<t1+enTY0expC−d−k2Ht−δ222t+δ2B2t.
This completes the proof.


### 3.2. Global Positive Solution of System (4)

Now, we will study the global dynamics of ISDE (4). Based on reference [7], we just need to explore the following system (29), which is equivalent to the original system (4). Since the dynamics of CAD drugs and chemotherapy drugs have been studied, then(29)dX=R1X1−XK−k1HtX−aXY−U1Xdt+δ1XdB1t,t≠nT,dY=C−dY−k2HtY−cXYdt+δ2YdB2t,YnT+=1+enTYnT,t=nT.It is expedient to define a SDE (Stochastic Differential Equation) without pulsed effects:(30)dX1=X1R11−X1K−k1Ht+U1−a∏0<nT<t1+enTY1dt+δ1X1dB1t,dY1=Y1C−d−k2Ht−cX1dt+δ2Y1dB2t,with X10,Y10=X0,Y0∈R+2.Theorem 2.
For any initial valueX0,Y0, system (9) has a unique positive solution Xt,Yt and the solution Xt,Yt will remain in R+2.Proof.
According to Yang et al. [10, 52], system (10) has a unique global positive solution X1t,Y1t. Setting(31)Xt,Yt=X1t,∏0<nT<t1+enTY1t,thereafter, since X1t,Y1t is absolute continuous, which results in the absolute continuity of Xt,Yt for any t∈nT,n+1T⊂0,+∞,n∈ℵ.
Whent≠nT, the derivatives of (31) with respect to (30) lead to the following:(32)dXt=dX1t=X1R11−XK−a∏0<nT<t1+enTY1t−k1Ht+U1dt+δ1X1dB1t,dYt=∏0<nT<t1+enTdY1t=Y−cX1+C−d−k2Htdt+δ2Y1dB2t=−cXY+C−dY−k2HtYdt+δ2YdB2t.
Oncet=nT,(33)XnT+=limt⟶nT+X1t=X1nT=XnT,YnT+=limt⟶nT+∏0<iT<t1+aiTY1t=1+enT∏0<iT<nT1+aiTY1nT=1+enTYnT.
Consequently, system (29) has a global unique positive solution. This completes the proof of Theorem 2.Theorem 3.
IfC−d−k2Ht<0, then the positive solutions Zt=Xt,Yt of system (4) are stochastically ultimately bounded.Proof.
We also define equations without respect to the pulsed immunotherapy, then(34)dX1=R11−X1K−k1Ht−a∏0<nT<t1+enTY1−U1X1dt+δ1XdB1t,dY1=C−dY1−k2HtY1−cX1Y1dt+δ2Y1dB2t.
To prove the boundedness of the solutions of system (4), we just need to illustrate that the solutions of system (34) are stochastically ultimately bounded.
ForX1, we let Vt,X1=etX1q where q>0. Applying Itô’s formula yields(35)dVX1=et1+qR11−X1R1+Ka∏0<nT<t1+enTY1R1K−k1Ht+U1+q−12δ12X1qdt+qetδ1X1qdB1t≤et1+qR1+q−12δ12X1q−qR1KX1q+1dt+qetδ1XqdB1t≤etX1q−11+qR1+q−12δ12X1−qR1KX12dt+qetδ1X1qdB1t≤M1etdt+qetδ1X1qdB1t,where M1 is a positive constant. By some calculations,(36)EetX1qt≤X1q0+M1et−1.
Then,(37)EX1qt≤X1q0et+M11−e−t.
Consequently,(38)limsupt⟶+∞EX1qt≤M1.
ForY1, note that(39)dY1≤C−d−k2HtY1dt+δ2Y1dB2t.
Letφt be the solution of(40)dφt=C−d−k2Htφdt+δ2φdB2t,forallt≥0.
Then, forφt, using the same method as previous yields(41)Eφt=φ0+12EC−d−k2HtY2.
Due toC−d−k2Ht<0, (42)EY1≤Eφt≤φ0≤N1,where N1 is a positive constant, and then(43)EY1≤Eφt≤φ0≤N1.
That is to say,(44)limsupt⟶+∞EY1qt≤N1p≐M2,where M2 is a positive constant. Note that(45)X1t)2+Y1t)2)q/2≤2q/2X1t)q+Y1t)q.
It follows from systems (38) and (44) that(46)limsupt⟶+∞E∣Z∣q≤2q/2M1+M2<∞.
From Chebychev’s inequality, this proof is completed.Theorem 4.
The solution of system (4) is globally attractive.Proof.
Without loss of generality, we letZ1t=P1t,Q1t and Z2t=P2t,Q2t be any two solutions of system (9), which have initial values Pi0>0, Qi0>0, i=1,2. First of all, define the Lyapunov function as follows:(47)Vt=lnP1t−lnP2t+lnQ1t−lnY2t,where t>0 and t≠nT. Then, we compute the upper right derivative d+Vt of Vt and apply the Itô’s formula to system (29),(48)d+Vt=signP1t−P2tdlnP1t−lnP2t+signQ1t−Q2tdlnQ1t−lnQ2t=signP1t−P2t−R1KP1t−P2t−aQ1t−Q2tdt+signQ1t−Q2t−cP1t−P2tdt=−R1K+cP1t−P2t−aQ1t−Q2tdt≤−πP1t−P2t+Q1t−Q2tdt≐−πVtdt,where π=minR1/K+c,a. Second, when t=nT,(49)VnT+=lnP1nT+−lnP2nT++lnQ1nT+−lnQ2nT+=lnP1nT−lnP2nT+ln1+enTQ1nT−ln1+enTQ2nT=lnX1nT−lnP2nT+lnQ1nT−lnQ2nT=VnT.
From 0 tot, we integrate system (48) and compute the expectation of both sides, and we have(50)Vt≤V0−π∫0tVsds.
Hence,(51)Vt+π∫0tVsds≤V0<∞.
Vt>0 is always valid which gives rise to limt⟶+∞Vt=0.
That is to say,(52)limt⟶∞P1t−P2t=0andlimt⟶∞Q1t−Q2t=0.
This completes the proof.Theorem 5.
For0≤s<t, if ∏s≤nT<t1+enT≤MM>0 and 1+C−d−k2ιe−μT/1−e−μT≤0, then, the solution xt=Xt,Yt of system (9) satisfies(53)limt⟶∞EXt≤MN,where N=cK/4R1R1+1−k1ιe−μT/1−e−μT2.Proof.
DefineV1t=e1Xt+e2Yt, here e1=c>0 and e2=−a<0. For any t∈n−1T,nT, by use of Itô’s formula for system (9), we have(54)dV1t=e1dX+e2dY=LVtdt+e1δ1XdB1t+e2δ2YdB2t,with(55)LVt=e1XR11−XK−aY−k1Ht−U1+e2Y−cX+C−d−k2Ht.
Moreover, we define Lyapunov functionV2t=etV1t and employ Itô’s formula, so one can get that(56)dV2t=etV1tdt+etdV1t=etV1tdt+etLVtdt+e1δ1XtdB1t+e2δ2YtdB2t.
Integrating the abovementioned equation fromn−1T to t and computing the expectations, (57)EetV1t=en−1TV1n−1T+E∫n−1TtesV1s+LVsds.
Note that(58)LV+V1=e1R1X1−XK−U1X−aXY−k1HtX+e2C−dY−cXY−k2HtY+e1X+e2Y≤e1R1X−R1X2K−k1HtX+X+−ae1−ce2XY+e2Y1+C−d−k2Ht≤cK4R1R+1−k1ιe−μT1−e−μT2−aY1+C−d−k2ιe−μT1−e−μT≤N.
From systems (57) and (58), we deduce that(59)etEV1t≤en−1TV1n−1T−N+etN.
Taking both sides of the derivative of system (59) leads to the following:(60)dEV1t≤N−EV1tdt.
Whent=nT,(61)EV1nT+=e1EXnT++e2EYnT+=e1EXnT+e21+enTEYnT≤1+enTe1EXnT+e2EYnT=1+enTEV1nT.
By systems (60) and (61),(62)dEV1t≤N−EV1tdt,t≠nT,EV1nT+≤1+enTEV1nT,t=nT.
To proveEV1t is bounded, consider the following impulsive system:(63)dwt=N−wtdt,t≠nT,wnT+=1+enTwnT,t=nT.
The unique solution of system (63) is deduced by(64)wt=w0nt,0+N∫0tnt,sds,where nt,s=∏s≤nT<t1+enTe−t−s. The methods of references [7, 10] give rise to limt⟶+∞wt=MN. Based on the comparison theorem of impulsive differential equations [7, 10, 53, 54],(65)limt⟶+∞EV1t≤limt⟶+∞wt=MN.
This completes the proof.Remark 1.
Under certain conditions, Theorem5 indicates that there is an upper bound on the expectation of the solution of ISDE (4). Biologically, when the pulse immunotherapy or pulsed perturbations are bounded, the number of prostate cancer cells is controllable.


## 3.1. Tumor-free Solution

It is obvious that the dynamics of chemotherapeutic drugs can be described by(18)dHt=−μH,t≠nT,HnT+=HnT+ι,t=nT.Let us compute the expression for theT periodic solution HTt as the following:(19)HTt=ιe−μt−nT1−e−μT,where t∈nT,n+1T, HTnT+=ι/1−e−μT.Lemma 3 (see [10, 50, 51]).
A unique positive periodic solution of system (5) can be expressed by HTt, and it holds that limt⟶∞Ht=HTt, then for any ϵ>0 yields(20)HTt−ϵ<Ht<HTt+ϵandlimt⟶∞1t∫0tHTsds=ιμT.For the convenience of writing, we useHt instead of HTt. Suppose that the prostate cancer cells can be eradicated, then let Xt=0, and system (4) becomes(21)dY=C−dY−k2HtYdt+δ2YdB2t,dH=−μH,t≠nT,HnT+=HnT+ι,YnT+=1+enTYnT,t=nT.Since the explicit expression ofHt is given, then(22)dY=C−dY−k2HtYdt+δ2YdB2t,t≠nT,YnT+=1+enTYnT,t=nT.Theorem 1.
There exists a unique global positive solutionYt of system (8) for any initial value Y0+=Y0, where(23)Yt=∏0<nT<t1+enTY0expC−d−k2Ht−δ222t+δ2B2t.Proof.
For anyt∈nT,n+1T, define a Lyapunov function Vt=lnYt. Applying Itô’s formula gives(24)dlnYt=C−d−k2Ht−0.5δ22dt+δ2dB2t.
We integrate the above equation fromnT to t, then(25)lnYt−lnYnT=C−d−k2Ht−0.5δ22t−nT+δ2B2t−B2nT.
Therefore,(26)Yt=YnTexpC−d−k2Ht−0.5δ22t−nT+δ2B2t−B2nT.
Whent=nT+, after a single immunotherapy, we have(27)Yt=∏0<nT<t1+enTYnTexpC−d−k2Ht−δ222t−nT+δ2B2t−B2nT.
An application of mathematical induction yields(28)Yt=∏0<nT<t1+enTY0expC−d−k2Ht−δ222t+δ2B2t.
This completes the proof.


## 3.2. Global Positive Solution of System (4)

Now, we will study the global dynamics of ISDE (4). Based on reference [7], we just need to explore the following system (29), which is equivalent to the original system (4). Since the dynamics of CAD drugs and chemotherapy drugs have been studied, then(29)dX=R1X1−XK−k1HtX−aXY−U1Xdt+δ1XdB1t,t≠nT,dY=C−dY−k2HtY−cXYdt+δ2YdB2t,YnT+=1+enTYnT,t=nT.It is expedient to define a SDE (Stochastic Differential Equation) without pulsed effects:(30)dX1=X1R11−X1K−k1Ht+U1−a∏0<nT<t1+enTY1dt+δ1X1dB1t,dY1=Y1C−d−k2Ht−cX1dt+δ2Y1dB2t,with X10,Y10=X0,Y0∈R+2.Theorem 2.
For any initial valueX0,Y0, system (9) has a unique positive solution Xt,Yt and the solution Xt,Yt will remain in R+2.Proof.
According to Yang et al. [10, 52], system (10) has a unique global positive solution X1t,Y1t. Setting(31)Xt,Yt=X1t,∏0<nT<t1+enTY1t,thereafter, since X1t,Y1t is absolute continuous, which results in the absolute continuity of Xt,Yt for any t∈nT,n+1T⊂0,+∞,n∈ℵ.
Whent≠nT, the derivatives of (31) with respect to (30) lead to the following:(32)dXt=dX1t=X1R11−XK−a∏0<nT<t1+enTY1t−k1Ht+U1dt+δ1X1dB1t,dYt=∏0<nT<t1+enTdY1t=Y−cX1+C−d−k2Htdt+δ2Y1dB2t=−cXY+C−dY−k2HtYdt+δ2YdB2t.
Oncet=nT,(33)XnT+=limt⟶nT+X1t=X1nT=XnT,YnT+=limt⟶nT+∏0<iT<t1+aiTY1t=1+enT∏0<iT<nT1+aiTY1nT=1+enTYnT.
Consequently, system (29) has a global unique positive solution. This completes the proof of Theorem 2.Theorem 3.
IfC−d−k2Ht<0, then the positive solutions Zt=Xt,Yt of system (4) are stochastically ultimately bounded.Proof.
We also define equations without respect to the pulsed immunotherapy, then(34)dX1=R11−X1K−k1Ht−a∏0<nT<t1+enTY1−U1X1dt+δ1XdB1t,dY1=C−dY1−k2HtY1−cX1Y1dt+δ2Y1dB2t.
To prove the boundedness of the solutions of system (4), we just need to illustrate that the solutions of system (34) are stochastically ultimately bounded.
ForX1, we let Vt,X1=etX1q where q>0. Applying Itô’s formula yields(35)dVX1=et1+qR11−X1R1+Ka∏0<nT<t1+enTY1R1K−k1Ht+U1+q−12δ12X1qdt+qetδ1X1qdB1t≤et1+qR1+q−12δ12X1q−qR1KX1q+1dt+qetδ1XqdB1t≤etX1q−11+qR1+q−12δ12X1−qR1KX12dt+qetδ1X1qdB1t≤M1etdt+qetδ1X1qdB1t,where M1 is a positive constant. By some calculations,(36)EetX1qt≤X1q0+M1et−1.
Then,(37)EX1qt≤X1q0et+M11−e−t.
Consequently,(38)limsupt⟶+∞EX1qt≤M1.
ForY1, note that(39)dY1≤C−d−k2HtY1dt+δ2Y1dB2t.
Letφt be the solution of(40)dφt=C−d−k2Htφdt+δ2φdB2t,forallt≥0.
Then, forφt, using the same method as previous yields(41)Eφt=φ0+12EC−d−k2HtY2.
Due toC−d−k2Ht<0, (42)EY1≤Eφt≤φ0≤N1,where N1 is a positive constant, and then(43)EY1≤Eφt≤φ0≤N1.
That is to say,(44)limsupt⟶+∞EY1qt≤N1p≐M2,where M2 is a positive constant. Note that(45)X1t)2+Y1t)2)q/2≤2q/2X1t)q+Y1t)q.
It follows from systems (38) and (44) that(46)limsupt⟶+∞E∣Z∣q≤2q/2M1+M2<∞.
From Chebychev’s inequality, this proof is completed.Theorem 4.
The solution of system (4) is globally attractive.Proof.
Without loss of generality, we letZ1t=P1t,Q1t and Z2t=P2t,Q2t be any two solutions of system (9), which have initial values Pi0>0, Qi0>0, i=1,2. First of all, define the Lyapunov function as follows:(47)Vt=lnP1t−lnP2t+lnQ1t−lnY2t,where t>0 and t≠nT. Then, we compute the upper right derivative d+Vt of Vt and apply the Itô’s formula to system (29),(48)d+Vt=signP1t−P2tdlnP1t−lnP2t+signQ1t−Q2tdlnQ1t−lnQ2t=signP1t−P2t−R1KP1t−P2t−aQ1t−Q2tdt+signQ1t−Q2t−cP1t−P2tdt=−R1K+cP1t−P2t−aQ1t−Q2tdt≤−πP1t−P2t+Q1t−Q2tdt≐−πVtdt,where π=minR1/K+c,a. Second, when t=nT,(49)VnT+=lnP1nT+−lnP2nT++lnQ1nT+−lnQ2nT+=lnP1nT−lnP2nT+ln1+enTQ1nT−ln1+enTQ2nT=lnX1nT−lnP2nT+lnQ1nT−lnQ2nT=VnT.
From 0 tot, we integrate system (48) and compute the expectation of both sides, and we have(50)Vt≤V0−π∫0tVsds.
Hence,(51)Vt+π∫0tVsds≤V0<∞.
Vt>0 is always valid which gives rise to limt⟶+∞Vt=0.
That is to say,(52)limt⟶∞P1t−P2t=0andlimt⟶∞Q1t−Q2t=0.
This completes the proof.Theorem 5.
For0≤s<t, if ∏s≤nT<t1+enT≤MM>0 and 1+C−d−k2ιe−μT/1−e−μT≤0, then, the solution xt=Xt,Yt of system (9) satisfies(53)limt⟶∞EXt≤MN,where N=cK/4R1R1+1−k1ιe−μT/1−e−μT2.Proof.
DefineV1t=e1Xt+e2Yt, here e1=c>0 and e2=−a<0. For any t∈n−1T,nT, by use of Itô’s formula for system (9), we have(54)dV1t=e1dX+e2dY=LVtdt+e1δ1XdB1t+e2δ2YdB2t,with(55)LVt=e1XR11−XK−aY−k1Ht−U1+e2Y−cX+C−d−k2Ht.
Moreover, we define Lyapunov functionV2t=etV1t and employ Itô’s formula, so one can get that(56)dV2t=etV1tdt+etdV1t=etV1tdt+etLVtdt+e1δ1XtdB1t+e2δ2YtdB2t.
Integrating the abovementioned equation fromn−1T to t and computing the expectations, (57)EetV1t=en−1TV1n−1T+E∫n−1TtesV1s+LVsds.
Note that(58)LV+V1=e1R1X1−XK−U1X−aXY−k1HtX+e2C−dY−cXY−k2HtY+e1X+e2Y≤e1R1X−R1X2K−k1HtX+X+−ae1−ce2XY+e2Y1+C−d−k2Ht≤cK4R1R+1−k1ιe−μT1−e−μT2−aY1+C−d−k2ιe−μT1−e−μT≤N.
From systems (57) and (58), we deduce that(59)etEV1t≤en−1TV1n−1T−N+etN.
Taking both sides of the derivative of system (59) leads to the following:(60)dEV1t≤N−EV1tdt.
Whent=nT,(61)EV1nT+=e1EXnT++e2EYnT+=e1EXnT+e21+enTEYnT≤1+enTe1EXnT+e2EYnT=1+enTEV1nT.
By systems (60) and (61),(62)dEV1t≤N−EV1tdt,t≠nT,EV1nT+≤1+enTEV1nT,t=nT.
To proveEV1t is bounded, consider the following impulsive system:(63)dwt=N−wtdt,t≠nT,wnT+=1+enTwnT,t=nT.
The unique solution of system (63) is deduced by(64)wt=w0nt,0+N∫0tnt,sds,where nt,s=∏s≤nT<t1+enTe−t−s. The methods of references [7, 10] give rise to limt⟶+∞wt=MN. Based on the comparison theorem of impulsive differential equations [7, 10, 53, 54],(65)limt⟶+∞EV1t≤limt⟶+∞wt=MN.
This completes the proof.Remark 1.
Under certain conditions, Theorem5 indicates that there is an upper bound on the expectation of the solution of ISDE (4). Biologically, when the pulse immunotherapy or pulsed perturbations are bounded, the number of prostate cancer cells is controllable.


## 4. Extinction and Persistence

To investigate the threshold conditions for extinction and persistence of effector cells and prostate cancer cells. LetVt=lnX1t and set Xt=X1t,Yt=∏0<nT<t1+enTY1t, then using Itô’s formula to the first equation of system (30) gives(66)dlnX1t=R1−R1KX1−a∏0<nT<t1+enTY1−k1Ht−U1−δ122dt+δ1dB1t=R1−R1KX−aY−k1Ht−U1−δ122dt+δ1dB1t.Lemma 4.
DefineMit=∫0tδidBiti=1,2, then integrating the abovementioned equation from 0 to t yields(67)lnXt−lnX0=R1−U1−12δ12t−R1K∫0tXsds−k1∫0tHsds−a∫0tYsds+M1t.
Similarly,(68)lnYt−lnY0=∑0<nT<tln1+enT−d−C+12δ22t−c∫0tXsds−k2∫0tHsds+M2t.Theorem 6.
(1)
If(69)R1−U1−12δ12−k1ιμT<0,then the prostate cancer cells Xt are extinct.(2)
If(70)R1−U1−12δ12−k1ιμT=0,then the prostate cancer cells Xt become nonpersistent in the mean.(3)
If(71)R1−U1−12δ12−k1ιμT>0,(72)limt⟶+∞sup∑0<nT<tln1+enTt<d−C+12δ22+k2ιμT,then the prostate cancer cells Xt are weakly persistent in the mean.Proof.
(1)
From system (25),(73)1tlnXtX0=R1−U1−12δ12−R1K∫0tXsdst−k1∫0tHsdst−a∫0tYsdst+M1tt.Due to<Mit,Mit>=∫0tδi2ds, and in the light of Lemma 2, we have(74)limt⟶+∞Mitt=0.Taking the superior limit of equality (27) yields(75)limt⟶+∞suplnXtt≤R1−U1−12δ12−k1ιμT−R1KX∗t−aY∗t<0.It indicates thatlimt⟶+∞Xt=0, and the prostate cancer cells become extinct.(2)
For any fixedε>0, there is a constant t1 so that for all t≥t1,(76)k1∫0tHsdst>k1ιμT−ε2,M1tt≤ε2.Based on system (73),(77)1tlnXtX0=R1−U−12δ12−R1K∫0tXsdst−k1∫0tHsdst−a∫0tYsdst+M1tt≤R1−U−12δ12−k1ιμT−RK∫0tXsdst+M1tt+ε.Owing to Lemma1, if ε is small enough, then(78)limsupt⟶+∞1t∫0tXsds≤KR1−U1−1/2δ12−k1ι/μTR=0.SinceX∗t≥0, it indicates that limt⟶+∞sup1/t∫0tXsds≥0.Then,limt⟶∞sup1/t∫0tXsds=0; thus, the prostate cancer cells become nonpersistent in the mean.(3)
Whenlimt⟶∞supXt/t<0, the superior limit of system (27) gives rise to(79)R1KX∗t+aY∗t≥R1−U1−12δ12−k1ιμT−limsupt⟶∞lnXtt>0.
As a result,X∗t≥0. Else, for any t∗∈X∗t,t∗=0, then Y∗t,t∗>0. If X∗t,t∗=0, we use the same method mentioned previously,(80)limsupt⟶+∞lnYt,t∗t≤limsupt⟶+∞∑0<nT<tln1+enTt−d−C+12δ22−k2ιμT<0,which means that Y∗t,t∗=0; it is a contradiction that Y∗t,t∗>0, i.e., limt⟶∞Xt/t>0. Then, the prostate cancer cells are weakly persistent in the mean. This completes the proof.Theorem 7.
(1)
If(81)R1−U1−12δ12−k1ιμT<0(82)limsupt⟶∞∑0<nT<tln1+enTt<d−C+12δ22+k2ιμT,then the effector cellsYt are extinct.(2)
If(83)R1−U1−12δ12−k1ιμT=0,then the effector cellsYt become nonpersistent in the mean.(3)
If(84)limsupt⟶+∞∑0<nT<tln1+enTt−d−C+12δ22+k2ιμT−cKR1R1−U1−12δ12−k1ιμT>0,then the effector cellsYt are weakly persistent in the mean.Proof.
(1)
According to system (26), we can deduce(85)1tlnYtY0=∑0<nT<tln1+enTt−d−C+12δ22−c1t∫0tXsds−k21t∫0tHsds+M2tt.Taking the superior limit of system (85) gives(86)limsupt⟶+∞lnYtt≤limsupt⟶∞∑0<nT<tln1+enTt−d−C+12δ22+k2ιμT−cX∗t.Then,(87)limsupt⟶∞lnYtt≤0.Hence,limt⟶+∞Yt=0; then, the effector cells are extinct.(2)
For any fixedε>0, there is a t2 so that for all t≥t2, one can obtain(88)k1∫0tHsdst>k1ιμT−ε2,M1tt≤ε2,in the light of system (73) yields(89)a1t∫0tYsds=−1tlnXtX0+R1−U1−12δ12−R1K∫0tXsdst−k1∫0tHsdst+M1tt≤R1−U1−12δ12−k1ιμT+ε.Asε is small enough, we take the superior limit that leads to the following:(90)limsupt⟶+∞1t∫0tYsds≤0.Thus,limsupt⟶∞1/t∫0tYsds=0; then, the effector cells become nonpersistent in the mean.(3)
From system (30), we know(91)1tlnYtY0=∑0<nT<tln1+enTt−d−C+12δ22−c1t∫0tXsds−k21t∫0tHsds+M2tt.
Adding systems (73) and (85) and taking the superior limit, there is a t3>0 so that(92)lnYtt≥limt⟶+∞sup∑0<nT<tln1+enTt−d−C+12δ22+k2ιμT+R1−U1−12δ12−k1ιμT−R1K+cX∗t−aY∗t.
It follows from system (78), and then(93)aY∗t≥limt⟶+∞sup∑0<nT<tln1+enTt−d−C+12δ22+k2ιμT−cKR1R1−U1−12δ12−k1ιμT>0.
Hence,Y∗t=limsupt⟶+∞1/t∫0tYsds≥0; then, the effector cells are weakly persistent in the mean. This completes the proof.Assumption 1.
There are two positive constantsm1 and M1 such that m1≤∏0<nT<t1+enT≤M1.Theorem 8.
Based on Assumption1, if ς=mint≥0R1−U1−0.5δ12−k1ι/μT−aM1y0>0, the prostate cancer cells become stochastically permanent.Proof.
First, we want to show that there are two constantsβ>0 and ϱ>0 so that liminft⟶+∞PXt≥β≥1−ε and liminft⟶+∞PXt≤ϱ≥1−ε for any ε∈0,1. For one thing, defining a Lyapunov function V1x=1/X1X1>0 and using Itô’s formula to the first equation of system (10) lead to the following:(94)dV1X1=−V1X1R1−U1−R1KX1−k1Ht−a∏0<nT<t1+enTY1dt+V1X1δ12dt−V1X1δ1dB1t.
Select a constantω>0 which satisfies ς>0.5ωδ12 and establish V2X1=1+V1X1)ω. An application of Itô’s formula yields(95)dV2X1=ω1+V1X1ω−1dV1X1+0.5ωω−11+V1X1ω−2dV1X12=ω1+V1X1ω−2−V1X1−V1X1R1−U1−R1KX1−k1Ht−a∏0<nT<t1+enTY1+V1X1+V1X12δ12+0.5ω−1V1X12δ12dt−ω1+V1X1ω−1V1X1δ1dB1t=ω1+V1X1ω−2−V1X12R1−U1−0.5δ12−k1Ht−0.5ωδ12−a∏0<nT<t1+enTY1+V1X1−R1+U1+R1K+a∏0<nT<t1+enTY1+k1Ht+δ12+R1Kdt−ω1+V1X1ω−1V1X1δ1dB1t≤ω1+V1X1ω−2−V1X12ς−0.5ωδ12+V1X1U1+R1K+aM1Y0+δ12+k1ιμT+R1Kdt−ω1+V1X1ω−1V1X1δ1dB1t.
Thereafter, choose a sufficiently smallϱ such that(96)ς−0.5ωδ12>ϱω>0.
DefiningV3X1=expϱtV2X1 and employing Itô’s formula give(97)dV3X1=ϱexpϱtV2X1dt+expϱtdV2X1≤ωexpϱt1+V1X1)ω−2ϱ1+V1X1)2ω−V1X1)2ς−0.5ωδ12+V1X1U1+R1K+aM1Y0+δ12+k1ιμT+R1Kdt−ωexpϱt1+V1X1)ω−1V1X1δ1dB1t≐expϱtPX1dt−ωexpϱt1+V1X1)ω−1V1X1δ1dB1t,where(98)PX1=ω1+V1X1ω−2−ς−0.5ωδ12−ϱωV1X1)2+U1+R1K+aM1Y0+δ12+k1ιμT+2ϱωV1X1+R1K+ϱω.
LetE1=ς−0.5ωδ12−ϱ/ω, E2=U1+R1/K+aM1Y0+δ12+k1ι/μT+2ϱ/ω, and E3=R1/K+ϱ/ω; from (96), we know that E1>0, E2>0, and E3>0. Thus, we can rewrite PX1 as follows:(99)PX1=ω1+1X1ω−2−E1X12+E2X1+E3≐P1X1.
We show thatPX1 is upper bounded when X1>0. If 1/X1≥E2+E22+4E1E3/2E1≐λ1, then PX1≤0. If 0<1/X1≤λ1, then P1X1≤4E1E3+E22/4E1. Besides, if ω≥2, then ω1+1/X1ω−2≤ω1+λ1ω−2; if ω<2, then ω1+1/X1ω−2≤ω. Therefore, for X1>0, we always have PX1≤P0=λ24E1E3+E22/4E1, where λ2=maxω,ω1+λ1ω−2. Clearly, PX1 is always upper bounded. Furthermore,(100)dV3X1≤expϱtPX1dt−ωexpϱt1+V1X1)ω−1V1X1δ1dB1t≤P0expϱtdt−ωexpϱt1+V1X1)ω−1V1X1δ1dB1t.
From 0 tot, integrating the previous inequality and taking the expectation,(101)EV3X1t≤V3X10+P0ϱexpϱt.
SinceV3X1t=expϱt1+V1X1t)ω, then we have(102)EV3X1t=Eexpϱt1+V1X1t)ω≤V3X10+P0ϱexpϱt=1+V1X10)ω+P0ϱexpϱt.
Taking the superior limit, then(103)limsupt⟶+∞E1X1tω=limsupt⟶+∞EV1X1t)ω≤limsupt⟶+∞E1+V1X1t)ω≤P0ϱ.
From system (31), one can obtain Xt=X1t,(104)limsupt⟶+∞E1Xtω=limsupt⟶+∞E1X1tω≤P0ϱ≐PM.
Furthermore, using the same method as in references [7, 10], then(105)limsupt⟶+∞E1Xtω≤P0ϱ≐PM.
For arbitraryε>0, we denote β=ε1/ω/PM1/ω. Making use of Chebyshev’s inequality gives(106)limsupt⟶+∞PXt<β=limsupt⟶+∞P1Xωt>1βω≤limsupt⟶+∞E1/Xωtβ−ω=limsupt⟶+∞βωE1Xωt=ε.
Thus,(107)liminft⟶+∞PXt≥β≥1−ε.
For another, defining a Lyapunov functionV3X1t=X1qtX1>0 and in terms of Itô’s formula with the first equation of system (30), it indicates that(108)dV3X1t=qV3X1tR1−U1−R1KX1−k1Ht−a∏0<nT<t1+enTY1t+0.5q−1δ12dt+qδ1V3X1tdB1t≤qV3X1tR1−U1−R1KX1−k1ιμT+0.5q−1δ12dt+qδ1V3X1tdB1t.
Integrating it from 0 tot and then taking the expectation lead to(109)EV3X1t−EV3X10≤q∫0tEV3X1sR1−U1−R1KX1−k1ιμT+0.5q−1δ12ds.
The derivative of the upper formula gives(110)dEV3X1tdt≤qEV3X1tR1−U1−k1ιμT+0.5q−1δ12−qR1KEX1q+1t.
On account of Hölder’s inequality, it can be obtained that(111)dEV3X1tdt≤qEV3X1tR1−U1−k1ιμT+0.5q−1δ12−qR1KEX1qtq+1q.
Letmt=EV3X1t, one can obtain that(112)dmtdt≤qmtR1−U1−k1ιμT+0.5q−1δ12−R1Km1qt≤qmtR1−U1−k1ιμT+0.5qδ12−R1Km1qt.
From the standard comparison theorem, we conclude that(113)limsupt⟶+∞EX1qt=limsupt⟶+∞EV3X1t=limsupt⟶+∞mt≤R1−U1−k1ι/μT+0.5qδ12R1/Kq.
SinceXt=X1t,(114)limsupt⟶+∞EXqt=limsupt⟶+∞EX1qt≤KR1−U1−k1ι/μT+0.5qδ12R1q.
Similarly, it follows from Chebyshev’s inequality that(115)liminft⟶+∞PXt≤ϱ≥1−ε.
As a consequence, it follows from system (4) of Definition 3. This completes the proof of Theorem 8.


## 5. Stationary Distribution and Ergodicity of System (4)

Here, we wonder the existence of a unique ergodic stationary distribution of ISDE (4). If g is a bounded function on R+, we define gu=supt∈R+gt. For more details, readers can refer to references [3, 46, 55].Define a SDE that equivalents to the original system (4) without pulsed effects:(116)dX=XR11−X1K−k1Ht+U1−a∏0<nT<t1+enTYdt+δ1XdB1t,dY=YC−d−k2Ht−cXdt+δ2YdB2t.WithXt,Yt∈R+2.Theorem 9.
Assume(117)R1−U1−k1Ht−12δ12>0,C−d−k2Ht−12δ22>0.Then, system (4) has a unique ergodic stationary distribution.Proof.
Let(118)V1X,Y=lnX+lnY.
Then,(119)LV1X,Y=R1−U1−k1Ht−12δ12−R1XK−a∏0<nT<t1+enTY+C−d−k2Ht−12δ22−cX.
Define(120)V2X,Y=1Xθ+1Yθ,0<θ<1.
By some calculations, we obtain(121)LV2X,Y=−θX−θR1−U1−k1Ht−a∏0<nT<t1+enTY−θ+12δ12+R1θKX1−θ−θY−θC−d−k2Ht−cX−θ+12δ22≤−θX−θR1−U1−k1Ht−θ+12δ12+R1θKX1−θ−θY−θC−d−k2Ht−θ+12δ22.
Let(122)VX,Y=V1X,Y+V2X,Y=lnX+lnY+1Xθ+1Yθ.
We obtain(123)LVX,Y=LV1+LV2≤−R1K+cX+R1−U1−k1Ht−12δ12−a∏0<nT<t1+enTY+C−d−k2Ht−12δ22≤−θX−θR1−U1−k1Ht−θ+12δ12+R1θKX1−θ−θY−θC−d−k2Ht−θ+12δ22.
It is easy to obtain that(124)LVX,Y≤ϕ1X+ϕ2Y,where(125)ϕ1X=−R1K+cX+R1−U1−k1Ht−12δ12−θX−θR1−U1−k1Ht−θ+12δ12+R1θKX1−θ,ϕ2Y=−a∏0<nT<t1+enTY+C−d−k2Ht−12δ22−θY−θC−d−k2Ht−θ+12δ22.Case 1.
IfX⟶0+, then(126)LV=ϕ1X+ϕ2Y≤ϕ1X+ϕ2u⟶−∞.
IfY⟶0+, then(127)LV=ϕ1X+ϕ2Y≤ϕ1u+ϕ2Y⟶−∞.Case 2.
IfX⟶+∞, then(128)LV=ϕ1X+ϕ2Y≤ϕ1X+ϕ2u⟶−∞.
IfY⟶+∞, then(129)LV=ϕ1X+ϕ2Y≤ϕ1u+ϕ2Y⟶−∞.
Using the analysis method of Theorem 5.1 in reference [46], we conclude that when X⟶0+ or Y⟶0+ or X⟶+∞ or Y⟶+∞, we can get LV⟶−∞. Hence, we take sufficiently small ν>0 and let U:=ν,1/ν×ν,1/ν, then(130)LVX,Y≤−1forallX,Y∈IntR+2U.
Moreover,(131)∑i,j=12bijX,Yξiξj=δ12X2ξ12+δ22Y2ξ22≥minX,Y∈Uδ12X2,δ22Y2ξ2forallX,Y∈U,ξ1,ξ2∈R2.
It follows from Theorem9, which completes the proof.


## 6. Numerical Simulations

Next, numerical simulations are conducted to verify our results by the Milstein higher order method [7, 35, 56]. The approximate solution system (4) is obtained and the discretization equations of system (4) are as follows:(132)dXk+1=Xk+XkrAk1−XkK−d11−Aka0+d2−k1Hk−aYkΔt+δ1XkΔtξk+δ122Xkξk2−1Δt,dYk+1=Yk+YkC−d−k2Hk−cXkΔt+δ2YkΔtηk+δ222Ykηk2−1Δt,dAk+1=Ak+−γAk−a0−γa0uΔt,dHk+1=Hk−μHkΔt,and when t=nT system (4) implements pulsed therapy, i.e., when modk,T=0, then(133)Yk+1=1+akYk,Hk+1=Hk+bk,where ξk and ηkk=1,2,3,⋯ have a distribution N0,1, which represents independent Gaussian random variables. We set time increment Δt=0.01. We replace the periodic solution DTt with the maximum D∗, where D∗=ι/1−e−μT.


### 6.1. Extinction and Persistence of Tumor Cells and Effector Cells

The standard parameter values of the pulseless random system are selected from the classical references [3, 7, 10, 20], which are scientific and reliable to a certain extent.In Figure1(a), we fix parameters as in Figure 1(a). By a simple calculation,(134)R1−U1−12δ12−k1ιμT=−0.926<0,and by Theorem 6, the prostate cancer cells are extinct (Figure 1(a)). Set r=1.663, and keep other parameters as in Figure 1(a); it is easy to see(135)R1−U1−12δ12−k1ιμT=0.Figure 1
(a–c) The extinction, nonpersistence in the mean, and weakly persistence in the mean of prostate cancer cells. (d–f) The extinction, nonpersistence in the mean, and weakly persistence in the mean of effector cells. (a)r=1.2, C=0.47, δ1=2.5, and enT=0.05; (b) C=0.47, δ1=2.5, and enT=0.05; (c) δ1=2.5 and enT=0.05; (d) r=1.2, δ1=2.5, and enT=0.1; (e) r=1.663, δ1=2.5, and enT=0.3; (f) C=0.47 and enT=0.3. The initial values were fixed as X0,Y0=0.1,0.5, red for X0,Y0=10,0.5 and blue for X0,Y0=300,5, and other parameters were fixed as δ2=0.5, u=0.5, ι=0.1, T=100, A0=5, H0=0.2, a0=4, a=1, d1=0.2, d2=0.1, d=0.3, k1=0.5, k2=0.5, γ=0.08, μ=0.5, K=1000, and c=0.00311.
(a)(b)(c)(d)(e)(f)In terms of Theorem6, the prostate cancer cells are nonpersistent in the mean (Figure 1(b)). Let r=1.85andC=0.4, then(136)R1−U1−12δ12−k1ιμT=0.374>0,and(137)limsupt⟶∞∑0<nT<tln1+enTt≈0.01<d−C+12δ22+k2ιμT=0.026.From Theorem6, the prostate cancer cells are weakly persistent in the mean (Figure 1(c)). It is easy to see that as the growth rate increases, the dynamic behaviors of prostate cancer cells changes from extinction to persistence, which means that the growth rate of tumors is positively correlated with their persistence under certain conditions (Figures 1(a)–1(c)).In Figure1(d), set C=0.4, by computing,(138)R1−U1−12δ12−k1ιμT=−0.926<0,and(139)limsupt⟶∞∑0<nT<tln1+enTt≈0.01<d−C+12δ22+k2ιμT=0.026.According to Theorem7, the effector cells are extinct. In Figure 1(e), we set C=0.7andak=0.3 and keep all other parameters as shown in Figure 1(b). One can get that(140)R1−U1−12δ12−k1ιμT=0.From Theorem7, the effector cells become nonpersistent in the mean. In Figure 1(f), r=1.663,δ1=3,andak=0.3; keeping all other parameters as shown in Figure 1(a), we obtain(141)limt⟶+∞sup∑0<nT<tln1+enTt−d−C+12δ22+k2ιμT−cKR1R1−U1−12δ12−k1ιμT≈1.3297>0.Theorem7 suggests that the effector cells are weakly persistent in the mean.In Figure2(a), setting r=3andδ1=2 and keeping all other parameters as shown in Figure 1(c), then(142)ς=mint≥0R1−U1−12δ12−k1ιμT−aM1y0=2.694>0.Figure 2
Stochastic permanence of prostate cancer cells. (a)r=3andδ1=2; (b) r=3andδ1=0.5. Other parameters were fixed as δ2=0.5, u=0.5, ι=0.1, C=0.4, enT=0.05, and T=100.
(a)(b)Theorem8 means that the prostate cancer cells become stochastically permanent (Figure 2(a)). In Figure 2(b), let δ1=0.5 and it is easy to see that the amplitude becomes smaller (Figure 2(b)), and then(143)ς=mint≥0R1−U1−12δ12−k1ιμT−aM1y0=4.569.It is concluded that with the increase ofς, the stronger the random persistence, the smaller the corresponding amplitude (Figure 2).


### 6.2. Effects of Random Perturbation on the Dynamics of Tumor Cells

To investigate how these disturbances affect the dynamical behaviors of prostate cancer cells, we do the following things.In Figure3(a), we set δ1=1 and fix all other parameters as Figure 1(b) and then(144)ι=mint≥0R1−U1−12δ12−k1ιμT−aM1y0=1.52>0.Figure 3
The effects of white noise on the evolution of prostate cancer cells. (a)r=1.663andδ1=1; (b) r=1.663andδ1=2.2; (c) r=1.663andδ1=2.5; (d) r=1.663andδ1=3. Other parameters were fixed as δ2=0.5, u=0.5, ι=0.1, C=0.47, enT=0.05, and T=100.
(a)(b)(c)(d)By Theorem8, the prostate cancer cells are stochastically persistent (Figure 3(a)). If we set δ1=2.2, by some computations,(145)ι=mint≥0R1−U1−12δ12−k1ιμT−aM1y0=−0.4<0,(146)R1−U1−12δ12−k1ιμT=0.705>0,and(147)limsupt⟶∞∑0<nT<tln1+enTt≈0.01<d−C+12δ22+k2ιμT=0.026.In the light of Theorem6, the dynamics of prostate cancer cells changes from stochastically persistent to weakly persistent in the mean (Figure 3(b)). If we set δ1=2.5, we have(148)R1−U1−12δ12−k1ιμT=0.Theorem6 explains that the dynamic behaviors of prostate cancer cells change from weakly persistent in the mean to nonpersistent in the mean (Figure 3(c)). If we set δ1=3, then(149)R1−U1−12δ12−k1ιμT=−1.375<0.By Theorem6, the dynamics of prostate cancer cells change from nonpersistent in the mean to extinct (Figure 3(d)). It can be seen that white noise can determine all the dynamic behaviors of prostate cancer cells. Increasing the intensities of white noise can accelerate the death of cancer cells, which indicates that white noise has a great influence on the evolution of tumors.


### 6.3. Monotherapy and Comprehensive Therapy

In this subsection, the initial values were fixed asX0,Y0=10,0.5, and other parameters are fixed as shown in (Figure 1(c)).If only CAD therapy is used, we just need to adjust the value ofu, and the feasible method of this therapy is to increase its intensities (Figure 4). If only chemotherapy is applied (Figures 5), the feasible methods of chemotherapy include increasing the dosages of chemotherapy (Figures 5(a) and 5(b)) and decreasing the pulsed periods (Figures 5(c) and 5(d)). If only immunotherapy is initiated, the feasible options of immunotherapy are to increase the dosages of immunotherapy (Figures 6(a) and 6(b)) or reduce the duration of immunotherapy (Figures 6(c) and 6(d)).Figure 4
The effects of CAD alone on the evolution of prostate cancer cells. (a)u=0.1; (b) u=0.2; (c) u=0.3; (d) u=0.4; (e) u=0.5; (f) u=0.6; (g) u=0.7; (h) u=0.8; (i) u=0.9. Other parameters were fixed as r=1.85, δ1=2.5, δ2=0.5, T=100, ι=0, C=0.4, and enT=0.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 5
The effects of chemotherapy alone on the evolution of prostate cancer cells. (a)ι=0.2,T=50; (b) ι=2,T=50; (c) ι=0.2,T=20; (d) ι=0.2,T=80 (red for (c), (d)). Other parameters were fixed as r=1.85, δ1=2.5, δ2=0.5, u=0, C=0.4, and enT=0.
(a)(b)(c)(d)Figure 6
The effects of immunotherapy alone on the evolution of prostate cancer cells. (a)enT=0.5,T=50; (b) enT=1.7,T=50; (c) enT=0.5,T=20; (d) enT=0.5,T=80 (red for (c), (d)). Other parameters were fixed as r=1.85, δ1=2.5, δ2=0.5, ι=0, C=0.4, and u=0.
(a)(b)(c)(d)From (Figures4–6), it has been observed that although the extinction of prostate cancer cells can be achieved by any of the previous single treatments under large stochastic fluctuations, these three kinds of treatments have their drawbacks such as resistance and toxic reaction [44, 45]. Hence, we will show how combination therapy of immunotherapy and chemotherapy together with CAD affects the evolution of prostate cancer cells. It can be seen from (Figures 4–6) that appropriately increasing the intensities of CAD therapy (or chemotherapy), increasing the dosages of chemotherapy (or immunotherapy), or shortening the treatment periods of immunotherapy (or increasing the frequencies of treatment) is more conducive to cancer treatment. Furthermore, compared with (Figures 4–7), we found that the tiny changes in hybrid therapy may cause big changes, such as to shorten the death time of prostate cancer cells (Figure 7).Figure 7
The effects of comprehensive therapy on the evolution of prostate cancer cells. (a)T=50,ι=2,enT=1.7,andu=0.6; (b) T=20,ι=0.2,enT=0.5,andu=0.2; (c) T=40,ι=2,enT=1.7,andu=0.6. Other parameters were fixed as r=1.85, δ1=2.5, δ2=0.5, and C=0.4.
(a)(b)(c)Next, we will explore the development of prostate cancer cells with different initial values under the same parameters, we setak=1.3, u=0.5, bk=0.1, and T=30, and other parameters were fixed as shown in Figure 1(c).In Figure8, it is revealed that the treatment of prostate cancer is associated with the initial tumor state of each patient. Therefore, the medical profession should formulate corresponding treatment strategies according to the development of patients’ condition.Figure 8
The dynamics of prostate cancer cells with different initial values under optimal strategy. (a)X0,Y0=0.1,0.5; (b) X0,Y0=10,0.5; (c) X0,Y0=300,5; other parameters were fixed as r=1.85, δ1=2.5, δ2=0.5, T=30, ι=0.1, C=0.4, u=0.5, and enT=1.3.
(a)(b)(c)


### 6.4. Existence of Stationary Distribution of System (4)

With the parameters as shown in Figure9, one can see that all conditions of Theorem 9 are satisfied. If we set δ1=0.1,δ2=0.1,u=0.1,ak=0.1,bk=0.1, by a calculation, then(150)R1−U1−k1Ht−12δ12≈6.485>0,C−d−k2Ht−12δ22≈0.045>0.Figure 9
Stationary distribution of deterministic model and stochastic model. (a)X0,Y0=0.1,0.5; (b) X0,Y0=10,0.5; (c) X0,Y0=300,5; other parameters were fixed as r=1.85, C=0.4, and T=100.
(a)(b)(c)By Theorem9, system (4) has a unique stationary distribution (Figure 9). The existence of stationary distribution shows that the dynamic behaviors of prostate cancer cells fluctuate within a certain range when the white noise is small.


## 6.1. Extinction and Persistence of Tumor Cells and Effector Cells

The standard parameter values of the pulseless random system are selected from the classical references [3, 7, 10, 20], which are scientific and reliable to a certain extent.In Figure1(a), we fix parameters as in Figure 1(a). By a simple calculation,(134)R1−U1−12δ12−k1ιμT=−0.926<0,and by Theorem 6, the prostate cancer cells are extinct (Figure 1(a)). Set r=1.663, and keep other parameters as in Figure 1(a); it is easy to see(135)R1−U1−12δ12−k1ιμT=0.Figure 1
(a–c) The extinction, nonpersistence in the mean, and weakly persistence in the mean of prostate cancer cells. (d–f) The extinction, nonpersistence in the mean, and weakly persistence in the mean of effector cells. (a)r=1.2, C=0.47, δ1=2.5, and enT=0.05; (b) C=0.47, δ1=2.5, and enT=0.05; (c) δ1=2.5 and enT=0.05; (d) r=1.2, δ1=2.5, and enT=0.1; (e) r=1.663, δ1=2.5, and enT=0.3; (f) C=0.47 and enT=0.3. The initial values were fixed as X0,Y0=0.1,0.5, red for X0,Y0=10,0.5 and blue for X0,Y0=300,5, and other parameters were fixed as δ2=0.5, u=0.5, ι=0.1, T=100, A0=5, H0=0.2, a0=4, a=1, d1=0.2, d2=0.1, d=0.3, k1=0.5, k2=0.5, γ=0.08, μ=0.5, K=1000, and c=0.00311.
(a)(b)(c)(d)(e)(f)In terms of Theorem6, the prostate cancer cells are nonpersistent in the mean (Figure 1(b)). Let r=1.85andC=0.4, then(136)R1−U1−12δ12−k1ιμT=0.374>0,and(137)limsupt⟶∞∑0<nT<tln1+enTt≈0.01<d−C+12δ22+k2ιμT=0.026.From Theorem6, the prostate cancer cells are weakly persistent in the mean (Figure 1(c)). It is easy to see that as the growth rate increases, the dynamic behaviors of prostate cancer cells changes from extinction to persistence, which means that the growth rate of tumors is positively correlated with their persistence under certain conditions (Figures 1(a)–1(c)).In Figure1(d), set C=0.4, by computing,(138)R1−U1−12δ12−k1ιμT=−0.926<0,and(139)limsupt⟶∞∑0<nT<tln1+enTt≈0.01<d−C+12δ22+k2ιμT=0.026.According to Theorem7, the effector cells are extinct. In Figure 1(e), we set C=0.7andak=0.3 and keep all other parameters as shown in Figure 1(b). One can get that(140)R1−U1−12δ12−k1ιμT=0.From Theorem7, the effector cells become nonpersistent in the mean. In Figure 1(f), r=1.663,δ1=3,andak=0.3; keeping all other parameters as shown in Figure 1(a), we obtain(141)limt⟶+∞sup∑0<nT<tln1+enTt−d−C+12δ22+k2ιμT−cKR1R1−U1−12δ12−k1ιμT≈1.3297>0.Theorem7 suggests that the effector cells are weakly persistent in the mean.In Figure2(a), setting r=3andδ1=2 and keeping all other parameters as shown in Figure 1(c), then(142)ς=mint≥0R1−U1−12δ12−k1ιμT−aM1y0=2.694>0.Figure 2
Stochastic permanence of prostate cancer cells. (a)r=3andδ1=2; (b) r=3andδ1=0.5. Other parameters were fixed as δ2=0.5, u=0.5, ι=0.1, C=0.4, enT=0.05, and T=100.
(a)(b)Theorem8 means that the prostate cancer cells become stochastically permanent (Figure 2(a)). In Figure 2(b), let δ1=0.5 and it is easy to see that the amplitude becomes smaller (Figure 2(b)), and then(143)ς=mint≥0R1−U1−12δ12−k1ιμT−aM1y0=4.569.It is concluded that with the increase ofς, the stronger the random persistence, the smaller the corresponding amplitude (Figure 2).


## 6.2. Effects of Random Perturbation on the Dynamics of Tumor Cells

To investigate how these disturbances affect the dynamical behaviors of prostate cancer cells, we do the following things.In Figure3(a), we set δ1=1 and fix all other parameters as Figure 1(b) and then(144)ι=mint≥0R1−U1−12δ12−k1ιμT−aM1y0=1.52>0.Figure 3
The effects of white noise on the evolution of prostate cancer cells. (a)r=1.663andδ1=1; (b) r=1.663andδ1=2.2; (c) r=1.663andδ1=2.5; (d) r=1.663andδ1=3. Other parameters were fixed as δ2=0.5, u=0.5, ι=0.1, C=0.47, enT=0.05, and T=100.
(a)(b)(c)(d)By Theorem8, the prostate cancer cells are stochastically persistent (Figure 3(a)). If we set δ1=2.2, by some computations,(145)ι=mint≥0R1−U1−12δ12−k1ιμT−aM1y0=−0.4<0,(146)R1−U1−12δ12−k1ιμT=0.705>0,and(147)limsupt⟶∞∑0<nT<tln1+enTt≈0.01<d−C+12δ22+k2ιμT=0.026.In the light of Theorem6, the dynamics of prostate cancer cells changes from stochastically persistent to weakly persistent in the mean (Figure 3(b)). If we set δ1=2.5, we have(148)R1−U1−12δ12−k1ιμT=0.Theorem6 explains that the dynamic behaviors of prostate cancer cells change from weakly persistent in the mean to nonpersistent in the mean (Figure 3(c)). If we set δ1=3, then(149)R1−U1−12δ12−k1ιμT=−1.375<0.By Theorem6, the dynamics of prostate cancer cells change from nonpersistent in the mean to extinct (Figure 3(d)). It can be seen that white noise can determine all the dynamic behaviors of prostate cancer cells. Increasing the intensities of white noise can accelerate the death of cancer cells, which indicates that white noise has a great influence on the evolution of tumors.


## 6.3. Monotherapy and Comprehensive Therapy

In this subsection, the initial values were fixed asX0,Y0=10,0.5, and other parameters are fixed as shown in (Figure 1(c)).If only CAD therapy is used, we just need to adjust the value ofu, and the feasible method of this therapy is to increase its intensities (Figure 4). If only chemotherapy is applied (Figures 5), the feasible methods of chemotherapy include increasing the dosages of chemotherapy (Figures 5(a) and 5(b)) and decreasing the pulsed periods (Figures 5(c) and 5(d)). If only immunotherapy is initiated, the feasible options of immunotherapy are to increase the dosages of immunotherapy (Figures 6(a) and 6(b)) or reduce the duration of immunotherapy (Figures 6(c) and 6(d)).Figure 4
The effects of CAD alone on the evolution of prostate cancer cells. (a)u=0.1; (b) u=0.2; (c) u=0.3; (d) u=0.4; (e) u=0.5; (f) u=0.6; (g) u=0.7; (h) u=0.8; (i) u=0.9. Other parameters were fixed as r=1.85, δ1=2.5, δ2=0.5, T=100, ι=0, C=0.4, and enT=0.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 5
The effects of chemotherapy alone on the evolution of prostate cancer cells. (a)ι=0.2,T=50; (b) ι=2,T=50; (c) ι=0.2,T=20; (d) ι=0.2,T=80 (red for (c), (d)). Other parameters were fixed as r=1.85, δ1=2.5, δ2=0.5, u=0, C=0.4, and enT=0.
(a)(b)(c)(d)Figure 6
The effects of immunotherapy alone on the evolution of prostate cancer cells. (a)enT=0.5,T=50; (b) enT=1.7,T=50; (c) enT=0.5,T=20; (d) enT=0.5,T=80 (red for (c), (d)). Other parameters were fixed as r=1.85, δ1=2.5, δ2=0.5, ι=0, C=0.4, and u=0.
(a)(b)(c)(d)From (Figures4–6), it has been observed that although the extinction of prostate cancer cells can be achieved by any of the previous single treatments under large stochastic fluctuations, these three kinds of treatments have their drawbacks such as resistance and toxic reaction [44, 45]. Hence, we will show how combination therapy of immunotherapy and chemotherapy together with CAD affects the evolution of prostate cancer cells. It can be seen from (Figures 4–6) that appropriately increasing the intensities of CAD therapy (or chemotherapy), increasing the dosages of chemotherapy (or immunotherapy), or shortening the treatment periods of immunotherapy (or increasing the frequencies of treatment) is more conducive to cancer treatment. Furthermore, compared with (Figures 4–7), we found that the tiny changes in hybrid therapy may cause big changes, such as to shorten the death time of prostate cancer cells (Figure 7).Figure 7
The effects of comprehensive therapy on the evolution of prostate cancer cells. (a)T=50,ι=2,enT=1.7,andu=0.6; (b) T=20,ι=0.2,enT=0.5,andu=0.2; (c) T=40,ι=2,enT=1.7,andu=0.6. Other parameters were fixed as r=1.85, δ1=2.5, δ2=0.5, and C=0.4.
(a)(b)(c)Next, we will explore the development of prostate cancer cells with different initial values under the same parameters, we setak=1.3, u=0.5, bk=0.1, and T=30, and other parameters were fixed as shown in Figure 1(c).In Figure8, it is revealed that the treatment of prostate cancer is associated with the initial tumor state of each patient. Therefore, the medical profession should formulate corresponding treatment strategies according to the development of patients’ condition.Figure 8
The dynamics of prostate cancer cells with different initial values under optimal strategy. (a)X0,Y0=0.1,0.5; (b) X0,Y0=10,0.5; (c) X0,Y0=300,5; other parameters were fixed as r=1.85, δ1=2.5, δ2=0.5, T=30, ι=0.1, C=0.4, u=0.5, and enT=1.3.
(a)(b)(c)


## 6.4. Existence of Stationary Distribution of System (4)

With the parameters as shown in Figure9, one can see that all conditions of Theorem 9 are satisfied. If we set δ1=0.1,δ2=0.1,u=0.1,ak=0.1,bk=0.1, by a calculation, then(150)R1−U1−k1Ht−12δ12≈6.485>0,C−d−k2Ht−12δ22≈0.045>0.Figure 9
Stationary distribution of deterministic model and stochastic model. (a)X0,Y0=0.1,0.5; (b) X0,Y0=10,0.5; (c) X0,Y0=300,5; other parameters were fixed as r=1.85, C=0.4, and T=100.
(a)(b)(c)By Theorem9, system (4) has a unique stationary distribution (Figure 9). The existence of stationary distribution shows that the dynamic behaviors of prostate cancer cells fluctuate within a certain range when the white noise is small.


## 7. Conclusions

Many studies have pointed out that the development of prostate cancer cells is inevitably affected by environmental disturbances such as nutrients and temperature [3, 7, 10, 27]. To our knowledge, there is no conclusive evidence in the medical field that IAD therapy is better than CAD therapy, and experiments have shown that immunotherapy can promote chemotherapy [3, 7, 10, 13]. However, there are few studies that combine CAD, immunotherapy and chemotherapy to study the dynamics of prostate cancer. Based on these considerations, we proposed a pulsed stochastic model, which combined the three treatments mentioned previously by incorporating tumor antigenicity and density-dependent mortality.We first explore the pharmacokinetics of chemotherapy and obtain the expression of the tumor-free solution, indicating that the system has a unique global positive solution. Then, the global attraction of solution and the boundness of expectation are proved, which indicates that prostate cancer cells cannot grow indefinitely and can be controlled under limited pulse immunotherapy. Furthermore, the threshold conditions of extinction and persistence for prostate cancer cells and effector cells are provided by using the theorems of ISDEs and Itô’s formula. Moreover, the sufficient conditions of stochastically permanence of prostate cancer cells and the existence of ergodic stationary distribution of the system are established. Finally, numerical simulations are carried out to confirm the theoretical results and provide guidance for treatment. The details of the theoretical results are shown in Table2.Table 2
The main theoretical results and corresponding conditions.
ConditionsTheoretical resultsR1−U1−1/2δ12−k1ι/μT<0Xt: extinctR1−U1−1/2δ12−k1ι/μT=0Xt: nonpersistent in the meanR1−U1−1/2δ12−k1ι/μT>0 and limt⟶+∞sup∑0<nT<tln1+enT/t<d−C+1/2δ22+k2ι/μTXt: weakly persistent in the meanR1−U1−1/2δ12−k1ι/μT<0 and limt⟶∞sup∑0<nT<tln1+enT/t<d−C+1/2δ22+k2ι/μTYt: extinctR1−U1−1/2δ12−k1ι/μT=0Yt: nonpersistent in the meanlimt⟶+∞sup∑0<nT<tln1+enT/t−d−C+1/2δ22+k2ι/μT−cK/R1R1−U1−1/2δ12−k1ι/μT>0Yt: weakly persistent in the meanς=mint≥0R1−U1−0.5δ12−k1ι/μT−aM1y0>0Xt: stochastically permanentR1−U1−k1Ht−1/2δ12>0 and C−d−k2Ht−1/2δ22>0Xt: stationary distributionBiologically, we obtain the following conclusions: (1) White noise can change tumor dynamics and has a negative effect on the evolution of prostate cancer cells. (2) Compared with single treatment, comprehensive therapy can significantly reduce the time of tumor regression and prevent tumor recurrence. The results show that increasing the intensities of CAD therapy (or chemotherapy), increasing the dosages of chemotherapy (or immunotherapy), or shortening the treatment periods of immunotherapy (or increasing the frequencies of treatment) are feasible treatments for prostate cancer (Figures4–7). (3) This paper investigates the development of prostate cancer cells with different initial values under the same parameter conditions, declaring that the treatment of prostate cancer can be adjusted according to the initial tumor state (Figure 8). (4) The existence of stationary distribution shows that small noise means stochastic stability, while large noise is destructive to stability of the system and leads to cure (Figure 9).Some interesting questions deserve further investigation, and more effective but complex models can be studied. Firstly, without using the traditional Lyapunov function method, how to analyze the properties of impulsive-free boundary stochastic differential systems [32, 33]? Secondly, Markov chain is often used to simulate random factors in ecosystems. Based on the original model, what is the dynamic behaviours after considering such random factors? Besides, it is interesting to investigate IAD with pulses. We leave these questions for future work.



---

*Source: 6887913-2023-03-09.xml*