# Dynamic Analysis of a Tumor-Immune System under Allee Effect

Authors: Chunmei Zeng; Shaojuan Ma
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, we develop a definite tumor-immune model considering Allee effect. The deterministic model is studied qualitatively by mathematical analysis method, including the positivity, boundness, and local stability of the solution. In addition, we explore the effect of random factors on the transition of the tumor-immune system from a stable coexistence equilibrium point to a stable tumor-free equilibrium point. Based on the method of stochastic averaging, we obtain the expressions of the steady-state probability density and the mean first-passage time. And we find that the Allee effect has the greatest impact on the number of cells in the system when the Allee threshold value is within a certain range; the intensity of random factors could affect the likelihood of the system crossing from the coexistence equilibrium to the tumor-free equilibrium.


---

## Body

## 1. Introduction

With the continuous improvement of living standards, individuals pay more attention to their own health. As one of the most harmful diseases to human health and life, malignant tumor has aroused wide concern in the world. At the same time, it has also become a research hotspot in academia. Thompson [1] discussed strategies to target YAP/TAZ activity in cancer and suggested that YAP and TAZ are important drivers of solid tumor growth, metastasis, and resistance to therapy. Du et al. [2] represented a first-time report of NSCLC-intrinsic PD-1 expression and a potential mechanism by which PD-1 blockade may promote cancer growth. Studies on ascorbic acid, which chemosensitize colorectal cancer cells and synergistically inhibit tumor growth, were presented in [3]. Chaoul et al. [4] found that the lack of MHC class II molecules was beneficial to the infiltration of T-cells into tumors, thus enhancing the control of tumor growth. Villalobos et al. [5] demonstrated the anticancer effect of SOCS1 in prostate cancer and identified its potential inhibitory mechanism. Volz et al. [6] studied a new family of potent TLR9 agonists, EnanDIM, and it can inhibit tumor growth.Although there are some achievements on finding factors which can effectively inhibit or control tumor growth. What we should not ignore is that the human body itself has a mechanism that can resist tumor cells—the immune system. It can not only protect the body from various internal and external factors against the normal immune system but also defend the body from tumor interference so that the normal operation of the body tissues. This function of the immune system has urged mathematicians and biologists exert to the research on its dynamic properties. Adriana et al. [7] discussed the regulation of angiogenesis by innate immune cells in the tumor microenvironment, specific features, and roles of major players. Key steps for visualizing the complex interaction between malignant tumors and the immune system were proposed and discussed in [8]. Jin et al. [9] summarized the current understanding of autophagy and its regulation in cancer. Medzhitov and Janeway [10] introduced the role of innate immune mechanism in the development of organisms and illustrated that it is the first line of defense of human immunity. The dynamics between tumor and immune cells in microenvironment can be referred to [11–14]. Khajanchi [15] described how tumor cells evolve and survive after encounter with the immune system. Terry et al. [16] presented new insights on the role of EMT in tumor-immune escape.At present, the research on tumor-immune system has obtained above results, and Katrin et al. [17] showed that tumors have the Allee effect. According to the study of literature [18], the Allee effect is a correlation between any component of individual fitness and the number or density of the same individual. Since the Allee effect has been shown to change the optimal control decisions and control costs and estimate the potential risks in ecology, Katrin et al. [17] believed that it was also crucial for tumor growth control. However, so far, it has been neglected in tumor growth and persistence. Therefore, under the excitation of the above research, we attempt to add the Allee effect on the basis of Pillis [19] model, study its dynamic behavior, and explore the influence of Allee effect on tumor cells.The remainder of the paper is organized as follows. In Section2, we develop a definite tumor-immune model with the Allee effect and analyze the property of the deterministic model, including the positivity, boundness, and local stability. Section 3 is devoted to discuss the tumor-immune model with the Allee effect, which is disturbed by white noise. We carry out numerical simulations in Section 4. In Section 5, we take a discussion and draw a conclusion.


## 2. The Deterministic Tumor-Immune Model with Allee Effect

### 2.1. Establishment of Deterministic Model

In 2003, Pillis and Radunskaya [20, 21] proposed a mathematical model to describe the relationship between the tumor and immune system:(1)dEdt=s+ρETα+T−β1ET−d1E,dNdt=r1N1−b1N−β2NT,dTdt=r2T1−b2T−β3ET−β4NT.Yet, for the tumor itself, it has its own optimal density in the process of growth and reproduction, whether it is too dense or too thin will have an inhibitory effect on its growth. Based on this, we study the relationship between tumor cells, normal cells, and immune cells under the background of the Allee effect. The mathematical model is as follows:(2)dEdt=s+ρETα+T−β1ET−d1E,dNdt=r1N1−b1N−β2NT,dTdt=r2T1−b2TT−m−β3ET−β4NT,where E represents the number of immune cells, T is the amount of tumor cells, N can be thought of the quantity of normal cells, s and d1 represent the resource constant and mortality of effector immune cells, respectively, and m is Allee threshold and satisfies the condition of m⋅b2<1. Since the presence of tumor cells would activate effector immune cells, this term is expressed as ρET/α+T by Michaelis–Menten expression. Competition between effector immune cells and tumor cells would lead to the death of two types of cells, which are expressed as −β1ET and −β3ET. The maximum growth rate and environmental carrying capacity of normal cells are represented by r1 and b1−1, respectively. Similarly, the tumor cells are expressed by r2 and b2−1. β2andβ4 stand for the cell death rate caused by the competition between normal cells and tumor cells, respectively.In this study, we assume that the tumor is homogeneous. And according to the physical meaning of each variable, the initial conditions of model (2) are(3)E0=E0>0,N0=N0>0,T0=T0>0.


### 2.2. Positiveness and Boundedness of Solutions

Theorem 1.
For anyt greater than 0, the solution of system (2) is positive and bounded under the initial conditions (3), when d1>ρ.Proof.
(1)
Positiveness: integrating both sides of formula (2), we have the following relations:(4)Et=E0exp∫0tρTξα+Tξ−β1Tξ−d1dξ+s∫0texp∫0tρTξα+Tξ−β1Tξ−d1dξdθ,Nt=N0exp∫0tr11−b1Nξ−β2Tξdξ,Tt=T0exp∫0tr21−b2TξTξ−m−β3Eξ−β4Nξdξ.The initial condition isE0>0,N0>0,T0>0. When t>0, it can be obtained from the above equation:(5)Et>0,Nt>0,Tt>0.(2)
Boundedness: considering a functionV=Et+Nt+Tt when d1>ρ, according to equation (2), we can obtain that(6)V˙=s+ρETα+T−β1ET−d1E+r1N1−b1N−β2NT+r2T1−b2TT−m−β3ET−β4NT≤s+r14b1+R+η−minβ1,β2,β3,β4V.Then, the above equation can be written as follows:(7)dVdt+τV≤φ,whereτ=minβ1,β2,β3,β4 and φ=s+r1/4b1+R+η, where R is the maximum value obtained by r2T1−b2TT−m within the range of Tt>0 and η is a positive number selected which is greater than the difference between E+N+T and E+NT. We then have(8)0<VE,N,T<φτ1−exp−τt+VE0,N0,T0exp−τt.Whent⟶∞, 0<V<φ/τ. This shows that, under initial conditions (3), the solution of system (2) is bounded.


### 2.3. Classification of Equilibrium Points

According to equation (2), we get the following three kinds of equilibrium points:(1)
Tumor-free equilibrium point:P0=s/d1,1/b1,0.(2)
Coexisting equilibrium point:P1=sa+sT∗/βT∗+d1a+T∗−ρT∗,r1−β2T∗/r1b1,T∗, where T∗ satisfies the equation(9)fT∗=T4+BT3+CT2+DT+E=0,and the coefficients of equation (9) are(10)B=α−m−1b2+d1−ρβ1+β2β4r1r2b1b2,C=−mα+αd1+mρ−md1β1+m−αb2+ρ−r2d1r2b2β1+β1β4r1−αβ2+β2β4ρ−d1r1r2b1b2β1,D=αd1r1r2mb1+r1b1β3sα+r1αd1β4r1r2b1b2β1.(3)
Abnormal equilibrium point:p2=s/d1,0,0,P3=sα+sT∗/1T∗+d1α+T∗−ρT∗,0,T∗, where T∗ is determined by the equation(11)fT∗=A1T4+B1T3+C1T2+D1T+E1=0,and the coefficients of above equation are(12)A1=r2b2β1,B1=r2b2αβ1+r2b2d1−ρr2b2−β1r2−mr2b2β1,C1=r2b2αd1−r2αβ1−r2d1+r2ρ−mr2b2d1+mρr2b2+mβ1r2−mr2b2αβ1,D1=mr2αβ1+mr2d1−mr2ρ−αr2d1−αmd1r2b2+β3s,E1=αsβ3+αmr2d1.


### 2.4. Stability of Equilibrium Points

In this section, we analyze the local stability of the equilibrium points. The calculation shows that the number of normal cells at the abnormal equilibrium point is zero, which is not consistent with the practical significance of biology, so we only discuss the local asymptotic stability of tumor-free equilibrium point and coexistence equilibrium point.The Jacobian matrixJ of system (2) at the equilibrium point is as follows:(13)ρTα+T−β1T−d10αρEα+T2−β1E0−2r1b1N+r1−β2T−β2N−β3T−β4T−3r2b2T2+2mb2r2+2r2T−mr2−β3E−β4N.


#### 2.4.1. Stability of Tumor-Free Equilibrium PointP0

The eigenvalues of Jacobian matrix atP0 are(14)λ1=−d1<0,λ2=−r1<0,λ3=−mr2−sβ3d1−β4b1.And by calculation, the equilibrium solutionP0 is locally asymptotically stable if and only if −sβ3/r2d1−β4/r2b1<m<1/b2, while P0 is the saddle point with index two when m<−sβ3/r2d1−β4/r2b1.


#### 2.4.2. Stability of Coexistence Equilibrium PointP1

When the equilibrium point isP1, the characteristic equation of matrix (13) at this point is(15)J−λE=λ3+r1−C2−A2+β4N∗−β2T∗λ2+D2λ+E2=0,where T∗ satisfies equation (9), and(16)A2=ρT∗α+T∗−β1T∗−d1,B2=αρα+T∗2−β1sα+sT∗βT∗+d1a+T∗−ρT∗,C2=−3r2b2T∗+2mr2b2+2r2T∗−mr2−sβ3a+T∗βT∗+d1a+T∗−ρT∗,D2=A2−r1C2−β4N∗−r1A2+A2β2+B2β3+C2β2T∗−2β2β4N∗T∗,E2=r1A2C2−B2β2β3T∗−A2C2β2T∗−A2β4r1N∗+B2r1β3T∗+2A2β2β4N∗T∗,N∗=r1−β2T∗r1b1.According to Routh–Hurwitz criterion, we know that(17)a0=1,a1=r1+mr2+d1+3−2mr2b2+β1−2r2−ρα+T∗T∗+r1−β2T∗r1b1β4+β3−αρT∗α+T∗2+β1T∗sα+sT∗βT∗+d1a+T∗−ρT∗,a2=A2−r1C2−r1−β2T∗r1b1β4−r1A2+A2β2+β3B2+C2β2T∗−2β2β4T∗r1−β2T∗r1b1,a3=r1−β2T∗A2C2+r1−β2B2β3T∗+2β2T∗−r1r1−β2T∗r1b1A2β4.The determinant of Hurwitz can be calculated according to the above expressions:(18)Δ1=a1,Δ2=a1a0a3a2,Δ3=a1a10a3a2a100a3.Thus, whenΔ1>0,Δ2>0,andΔ3>0 are satisfied at the same time, coexistence balance P1 is locally asymptotically stable.By equation (9), we know that coexistence equilibrium expression T∗ is the root of quartic equation, namely, there are at most four coexistence equilibrium points. Considering the meaningful number of equilibrium and its stability is related to parameter selection, and in order to make the results more intuitive, we carry out numerical simulation according to effective range of variables. As a result, the number and stability of the coexistence equilibrium point are shown in Figure 1. In this numerical simulation, the horizontal axis is ρ value from 0 to 2, while the vertical axis is m value from −2 to 1.Figure 1
The number of coexistence equilibrium points and their stability.In order to reveal the role of the Allee effect in the system, we only change the value ofm. By referring to [20, 21], we make ρ=0.01; then, the time history diagrams of three cells in different m are as follows.Observing Figure2, we know when m∈−∞,−1.2∪−0.5,1, the number of the three types of cells do not change significantly. When m belongs to −1.2,−0.5, the cell change trend is the most significant. This indicates that when the Allee threshold of tumor cells is within the range of −1.2 to −0.5, the Allee effect has the greatest influence on the number of cells in the system, and this also proves that the Allee effect has an impact on the tumor-immune system.Figure 2
Time history diagram of differentm (ρ=0.01). (a) m = 0.5. (b) m = 0.1. (c) m = 0.0. (d) m = −0.5. (e) m = −0.7. (f) m = −0.8. (g) m = −1.0. (h) m = −1.2. (i) m = −1.5.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 2.1. Establishment of Deterministic Model

In 2003, Pillis and Radunskaya [20, 21] proposed a mathematical model to describe the relationship between the tumor and immune system:(1)dEdt=s+ρETα+T−β1ET−d1E,dNdt=r1N1−b1N−β2NT,dTdt=r2T1−b2T−β3ET−β4NT.Yet, for the tumor itself, it has its own optimal density in the process of growth and reproduction, whether it is too dense or too thin will have an inhibitory effect on its growth. Based on this, we study the relationship between tumor cells, normal cells, and immune cells under the background of the Allee effect. The mathematical model is as follows:(2)dEdt=s+ρETα+T−β1ET−d1E,dNdt=r1N1−b1N−β2NT,dTdt=r2T1−b2TT−m−β3ET−β4NT,where E represents the number of immune cells, T is the amount of tumor cells, N can be thought of the quantity of normal cells, s and d1 represent the resource constant and mortality of effector immune cells, respectively, and m is Allee threshold and satisfies the condition of m⋅b2<1. Since the presence of tumor cells would activate effector immune cells, this term is expressed as ρET/α+T by Michaelis–Menten expression. Competition between effector immune cells and tumor cells would lead to the death of two types of cells, which are expressed as −β1ET and −β3ET. The maximum growth rate and environmental carrying capacity of normal cells are represented by r1 and b1−1, respectively. Similarly, the tumor cells are expressed by r2 and b2−1. β2andβ4 stand for the cell death rate caused by the competition between normal cells and tumor cells, respectively.In this study, we assume that the tumor is homogeneous. And according to the physical meaning of each variable, the initial conditions of model (2) are(3)E0=E0>0,N0=N0>0,T0=T0>0.


## 2.2. Positiveness and Boundedness of Solutions

Theorem 1.
For anyt greater than 0, the solution of system (2) is positive and bounded under the initial conditions (3), when d1>ρ.Proof.
(1)
Positiveness: integrating both sides of formula (2), we have the following relations:(4)Et=E0exp∫0tρTξα+Tξ−β1Tξ−d1dξ+s∫0texp∫0tρTξα+Tξ−β1Tξ−d1dξdθ,Nt=N0exp∫0tr11−b1Nξ−β2Tξdξ,Tt=T0exp∫0tr21−b2TξTξ−m−β3Eξ−β4Nξdξ.The initial condition isE0>0,N0>0,T0>0. When t>0, it can be obtained from the above equation:(5)Et>0,Nt>0,Tt>0.(2)
Boundedness: considering a functionV=Et+Nt+Tt when d1>ρ, according to equation (2), we can obtain that(6)V˙=s+ρETα+T−β1ET−d1E+r1N1−b1N−β2NT+r2T1−b2TT−m−β3ET−β4NT≤s+r14b1+R+η−minβ1,β2,β3,β4V.Then, the above equation can be written as follows:(7)dVdt+τV≤φ,whereτ=minβ1,β2,β3,β4 and φ=s+r1/4b1+R+η, where R is the maximum value obtained by r2T1−b2TT−m within the range of Tt>0 and η is a positive number selected which is greater than the difference between E+N+T and E+NT. We then have(8)0<VE,N,T<φτ1−exp−τt+VE0,N0,T0exp−τt.Whent⟶∞, 0<V<φ/τ. This shows that, under initial conditions (3), the solution of system (2) is bounded.


## 2.3. Classification of Equilibrium Points

According to equation (2), we get the following three kinds of equilibrium points:(1)
Tumor-free equilibrium point:P0=s/d1,1/b1,0.(2)
Coexisting equilibrium point:P1=sa+sT∗/βT∗+d1a+T∗−ρT∗,r1−β2T∗/r1b1,T∗, where T∗ satisfies the equation(9)fT∗=T4+BT3+CT2+DT+E=0,and the coefficients of equation (9) are(10)B=α−m−1b2+d1−ρβ1+β2β4r1r2b1b2,C=−mα+αd1+mρ−md1β1+m−αb2+ρ−r2d1r2b2β1+β1β4r1−αβ2+β2β4ρ−d1r1r2b1b2β1,D=αd1r1r2mb1+r1b1β3sα+r1αd1β4r1r2b1b2β1.(3)
Abnormal equilibrium point:p2=s/d1,0,0,P3=sα+sT∗/1T∗+d1α+T∗−ρT∗,0,T∗, where T∗ is determined by the equation(11)fT∗=A1T4+B1T3+C1T2+D1T+E1=0,and the coefficients of above equation are(12)A1=r2b2β1,B1=r2b2αβ1+r2b2d1−ρr2b2−β1r2−mr2b2β1,C1=r2b2αd1−r2αβ1−r2d1+r2ρ−mr2b2d1+mρr2b2+mβ1r2−mr2b2αβ1,D1=mr2αβ1+mr2d1−mr2ρ−αr2d1−αmd1r2b2+β3s,E1=αsβ3+αmr2d1.


## 2.4. Stability of Equilibrium Points

In this section, we analyze the local stability of the equilibrium points. The calculation shows that the number of normal cells at the abnormal equilibrium point is zero, which is not consistent with the practical significance of biology, so we only discuss the local asymptotic stability of tumor-free equilibrium point and coexistence equilibrium point.The Jacobian matrixJ of system (2) at the equilibrium point is as follows:(13)ρTα+T−β1T−d10αρEα+T2−β1E0−2r1b1N+r1−β2T−β2N−β3T−β4T−3r2b2T2+2mb2r2+2r2T−mr2−β3E−β4N.


### 2.4.1. Stability of Tumor-Free Equilibrium PointP0

The eigenvalues of Jacobian matrix atP0 are(14)λ1=−d1<0,λ2=−r1<0,λ3=−mr2−sβ3d1−β4b1.And by calculation, the equilibrium solutionP0 is locally asymptotically stable if and only if −sβ3/r2d1−β4/r2b1<m<1/b2, while P0 is the saddle point with index two when m<−sβ3/r2d1−β4/r2b1.


### 2.4.2. Stability of Coexistence Equilibrium PointP1

When the equilibrium point isP1, the characteristic equation of matrix (13) at this point is(15)J−λE=λ3+r1−C2−A2+β4N∗−β2T∗λ2+D2λ+E2=0,where T∗ satisfies equation (9), and(16)A2=ρT∗α+T∗−β1T∗−d1,B2=αρα+T∗2−β1sα+sT∗βT∗+d1a+T∗−ρT∗,C2=−3r2b2T∗+2mr2b2+2r2T∗−mr2−sβ3a+T∗βT∗+d1a+T∗−ρT∗,D2=A2−r1C2−β4N∗−r1A2+A2β2+B2β3+C2β2T∗−2β2β4N∗T∗,E2=r1A2C2−B2β2β3T∗−A2C2β2T∗−A2β4r1N∗+B2r1β3T∗+2A2β2β4N∗T∗,N∗=r1−β2T∗r1b1.According to Routh–Hurwitz criterion, we know that(17)a0=1,a1=r1+mr2+d1+3−2mr2b2+β1−2r2−ρα+T∗T∗+r1−β2T∗r1b1β4+β3−αρT∗α+T∗2+β1T∗sα+sT∗βT∗+d1a+T∗−ρT∗,a2=A2−r1C2−r1−β2T∗r1b1β4−r1A2+A2β2+β3B2+C2β2T∗−2β2β4T∗r1−β2T∗r1b1,a3=r1−β2T∗A2C2+r1−β2B2β3T∗+2β2T∗−r1r1−β2T∗r1b1A2β4.The determinant of Hurwitz can be calculated according to the above expressions:(18)Δ1=a1,Δ2=a1a0a3a2,Δ3=a1a10a3a2a100a3.Thus, whenΔ1>0,Δ2>0,andΔ3>0 are satisfied at the same time, coexistence balance P1 is locally asymptotically stable.By equation (9), we know that coexistence equilibrium expression T∗ is the root of quartic equation, namely, there are at most four coexistence equilibrium points. Considering the meaningful number of equilibrium and its stability is related to parameter selection, and in order to make the results more intuitive, we carry out numerical simulation according to effective range of variables. As a result, the number and stability of the coexistence equilibrium point are shown in Figure 1. In this numerical simulation, the horizontal axis is ρ value from 0 to 2, while the vertical axis is m value from −2 to 1.Figure 1
The number of coexistence equilibrium points and their stability.In order to reveal the role of the Allee effect in the system, we only change the value ofm. By referring to [20, 21], we make ρ=0.01; then, the time history diagrams of three cells in different m are as follows.Observing Figure2, we know when m∈−∞,−1.2∪−0.5,1, the number of the three types of cells do not change significantly. When m belongs to −1.2,−0.5, the cell change trend is the most significant. This indicates that when the Allee threshold of tumor cells is within the range of −1.2 to −0.5, the Allee effect has the greatest influence on the number of cells in the system, and this also proves that the Allee effect has an impact on the tumor-immune system.Figure 2
Time history diagram of differentm (ρ=0.01). (a) m = 0.5. (b) m = 0.1. (c) m = 0.0. (d) m = −0.5. (e) m = −0.7. (f) m = −0.8. (g) m = −1.0. (h) m = −1.2. (i) m = −1.5.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 2.4.1. Stability of Tumor-Free Equilibrium PointP0

The eigenvalues of Jacobian matrix atP0 are(14)λ1=−d1<0,λ2=−r1<0,λ3=−mr2−sβ3d1−β4b1.And by calculation, the equilibrium solutionP0 is locally asymptotically stable if and only if −sβ3/r2d1−β4/r2b1<m<1/b2, while P0 is the saddle point with index two when m<−sβ3/r2d1−β4/r2b1.


## 2.4.2. Stability of Coexistence Equilibrium PointP1

When the equilibrium point isP1, the characteristic equation of matrix (13) at this point is(15)J−λE=λ3+r1−C2−A2+β4N∗−β2T∗λ2+D2λ+E2=0,where T∗ satisfies equation (9), and(16)A2=ρT∗α+T∗−β1T∗−d1,B2=αρα+T∗2−β1sα+sT∗βT∗+d1a+T∗−ρT∗,C2=−3r2b2T∗+2mr2b2+2r2T∗−mr2−sβ3a+T∗βT∗+d1a+T∗−ρT∗,D2=A2−r1C2−β4N∗−r1A2+A2β2+B2β3+C2β2T∗−2β2β4N∗T∗,E2=r1A2C2−B2β2β3T∗−A2C2β2T∗−A2β4r1N∗+B2r1β3T∗+2A2β2β4N∗T∗,N∗=r1−β2T∗r1b1.According to Routh–Hurwitz criterion, we know that(17)a0=1,a1=r1+mr2+d1+3−2mr2b2+β1−2r2−ρα+T∗T∗+r1−β2T∗r1b1β4+β3−αρT∗α+T∗2+β1T∗sα+sT∗βT∗+d1a+T∗−ρT∗,a2=A2−r1C2−r1−β2T∗r1b1β4−r1A2+A2β2+β3B2+C2β2T∗−2β2β4T∗r1−β2T∗r1b1,a3=r1−β2T∗A2C2+r1−β2B2β3T∗+2β2T∗−r1r1−β2T∗r1b1A2β4.The determinant of Hurwitz can be calculated according to the above expressions:(18)Δ1=a1,Δ2=a1a0a3a2,Δ3=a1a10a3a2a100a3.Thus, whenΔ1>0,Δ2>0,andΔ3>0 are satisfied at the same time, coexistence balance P1 is locally asymptotically stable.By equation (9), we know that coexistence equilibrium expression T∗ is the root of quartic equation, namely, there are at most four coexistence equilibrium points. Considering the meaningful number of equilibrium and its stability is related to parameter selection, and in order to make the results more intuitive, we carry out numerical simulation according to effective range of variables. As a result, the number and stability of the coexistence equilibrium point are shown in Figure 1. In this numerical simulation, the horizontal axis is ρ value from 0 to 2, while the vertical axis is m value from −2 to 1.Figure 1
The number of coexistence equilibrium points and their stability.In order to reveal the role of the Allee effect in the system, we only change the value ofm. By referring to [20, 21], we make ρ=0.01; then, the time history diagrams of three cells in different m are as follows.Observing Figure2, we know when m∈−∞,−1.2∪−0.5,1, the number of the three types of cells do not change significantly. When m belongs to −1.2,−0.5, the cell change trend is the most significant. This indicates that when the Allee threshold of tumor cells is within the range of −1.2 to −0.5, the Allee effect has the greatest influence on the number of cells in the system, and this also proves that the Allee effect has an impact on the tumor-immune system.Figure 2
Time history diagram of differentm (ρ=0.01). (a) m = 0.5. (b) m = 0.1. (c) m = 0.0. (d) m = −0.5. (e) m = −0.7. (f) m = −0.8. (g) m = −1.0. (h) m = −1.2. (i) m = −1.5.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 3. Stochastic Tumor-Immune Model with Allee Effect

### 3.1. Establishment of Stochastic Model

Considering that the growth of tumor cells can be affected by various uncertainties, such as vibration and temperature, and these can be simulated with noise, therefore, based on the deterministic model (2), we propose the following stochastic model:(19)E.1=s+ρE1T1α+T1−β1E1T1−d1E,N.1=r1N11−b1N1−β2N1T1+N1W1t,T.1=r2T11−b2T1T1−m−β3E1T1+β4N1T1+T1W2t.Among them,E1t, N1t, and T1t are three random processes, which represent the number of immune cells, normal cells, and tumor cells, respectively. Wit is the independent Gaussian white noise with zero mean intensity Dii=1,2, which acts on the natural growth rate of normal cells ρ and the nature growth rate of tumor r2, and we assume that both damping and excitation are small. The specific meanings of other parameters are the same as those in model (2).


### 3.2. The Stochastic Averaging Method

In order to facilitate the application of the stochastic average method, we first simplified model (19), setting(20)f=−r2m,g=r1b,G3=−γ1b1N1,G1E1,N1,T1=sE1+ρT1a+T1−β1T1,G2=rT1−b2T1+b2m−βE.Then, model (19) can be written as follows:(21)E˙1=E1−d1+G1E1,N1,T1,N˙1=N1r1−β2T1+G3E1,N1,T1+N1W1t,T˙1=T1−β4N1+f+G2E1,N1,T1+T1W2t.The above equation (21) is converted into the following Ito^ stochastic differential equation by Ito^ formula:(22)dE1=E1−d1+G1E1,N1,T1dt,dN1=N1r1−β2T1+G3E1,N1,T1+12D1dt+D1N1dB1t,dT1=T1f−β4N1+G2E1,N1,T1+12D2dt+D2T1dB2t,where B1t and B2t are independent Wiener process and 1/2D1E1dt and 1/2D2T1dt are the Wong–Zakai amendments. Defining a random process,(23)RN1,T1=β2T1−β4N1+f−r1+flnβ4N1f−r1lnβ2T1r1.Among them,Rf/β4,r1/β2=0. β2, β4, and r1 are positive, m<0. Then, the Ito^ stochastic differential equation of random process Rt is(24)dR=f−β4N1G3+β2T1−r1G2−12β4N1D1+12β2T1D2dt+D1f−β4N1dB1+D2β2T1−γ1dB2.Using the random average principle of Stratorovich-Khasminskii, the averageIto^ stochastic differential equation of Rt is(25)dR=mRdt+σRdBt.Whence,(26)mR=f−β4N1G3+β2T1−r1G2t+12β2D2T1t−12β4D1N1t,(27)σ2R=D1f−β4N12t+D2β2T1−r12t,and according to the formula,(28)⋅t=1TR∮⋅dt.We have the following expressions:(29)N1t=fβ4,T1t=r1β2,(30)r1f−β4N12t=fr1−β2T12t=gRTR,where(31)gR=r1∮f−β4N1T1dT1.Then, with the help of equations (26)–(30), the drift coefficient and diffusion coefficient can be simplified as follows:(32)mR=12r1D2−12fD1−b1gRTR,σ2R=fD1−r1D2r1fgRTR.


#### 3.2.1. Stationary Probability Density

Base on equation (23), we obtain that(33)lnβ4N1f≈1+β4N1−ff−12β4N1−ff2,(34)lnβ4N1r1≈1+β2T1−r1r1−12β2T1−r1r12.Substitute (33) and (34) into (23); then, it can be written as follows:(35)r1β4N1−f2+fβ2T1−r12≈2r1fR.Due to equations,(36)β4N2−f2t≈fR,r1−β2T12t≈r1R,thus, we obtain that(37)mR≈12r1D2−12fD1−b1r1fR,σ2R≈fD1−r1D2R.WhenR⟶4.3782, the boundary is nonsingular; when R⟶∞, it is the second kind of singular boundary. The diffusion index a=1, the drift index β=1,m+∞<0. Therefore, the nontrivial stationary probability density exists. And its steady-state probability density function pr satisfies the following FPK equation:(38)ddrmrpr−12d2dr2σ2rpr=0.Whence,(39)mr=mRR=r,σ2r=σ2RR=r.Solving the FPK equation, we then have(40)pr=CTrgrexp∫0rr1fr1D2−fD1Tx−2r1b1fgxfD1−r1D2gxdx,where C is a normalization constant, and(41)C=C1r1ffD1−r1D2,dgrdr=r1∮β4N1T1β4N1−fdT1.Hence, equation (40) can be simplified as follows:(42)pr=CTrexpυr,υ=−2r1b1ffD1−r1D2.SupposeR=r, the conditional probability density of random process N1 is pN1/r; then, the joint probability density of R and N1 can be expressed as follows:(43)pr1,N1=prN1r1−β2T1Tr.As a result, we obtain the stationary probability density ofN1 and T1:(44)pN1,T1=CN1T1expυr.Thus, we get the steady-state probability densities ofN1 and T1, respectively. They are(45)pN1=∫0∞pN1,T1dT1=∫0∞CN1T1expυrdT1,pT1=∫0∞pN1,T1dN1=∫0∞CN1T1expυrdN1.


#### 3.2.2. Mean First-Passage Time

When people are suffering from tumor, we hope that the number of tumor cells go to zero. However, according to the calculation in the above section, we know that the steady-state probability density of tumor cells is present, which indicates that tumor cells will not be extinct under natural state, but the effects of random factors may change this. Therefore, it is practical to study the average first-passage timeβ from the coexistence equilibrium point to the tumor-free equilibrium point, and the Pontryagin equation is satisfied that(46)1+mx0dβdx0+12σ2x0d2βdx02=0,where(47)mx0=mRR=x0,σ2x0=σ2RR=x0.From the above studies, we know thatxl is nonsingular and xs is natural. According to the constant variation method, β can be written as follows:(48)βx0=2∫x0xsψ−1μ∫xlμ1σ2υψυ+C1dμ,where C1 is the integral constant, and(49)ψυ=exp∫2mυσ2υdυ.


## 3.1. Establishment of Stochastic Model

Considering that the growth of tumor cells can be affected by various uncertainties, such as vibration and temperature, and these can be simulated with noise, therefore, based on the deterministic model (2), we propose the following stochastic model:(19)E.1=s+ρE1T1α+T1−β1E1T1−d1E,N.1=r1N11−b1N1−β2N1T1+N1W1t,T.1=r2T11−b2T1T1−m−β3E1T1+β4N1T1+T1W2t.Among them,E1t, N1t, and T1t are three random processes, which represent the number of immune cells, normal cells, and tumor cells, respectively. Wit is the independent Gaussian white noise with zero mean intensity Dii=1,2, which acts on the natural growth rate of normal cells ρ and the nature growth rate of tumor r2, and we assume that both damping and excitation are small. The specific meanings of other parameters are the same as those in model (2).


## 3.2. The Stochastic Averaging Method

In order to facilitate the application of the stochastic average method, we first simplified model (19), setting(20)f=−r2m,g=r1b,G3=−γ1b1N1,G1E1,N1,T1=sE1+ρT1a+T1−β1T1,G2=rT1−b2T1+b2m−βE.Then, model (19) can be written as follows:(21)E˙1=E1−d1+G1E1,N1,T1,N˙1=N1r1−β2T1+G3E1,N1,T1+N1W1t,T˙1=T1−β4N1+f+G2E1,N1,T1+T1W2t.The above equation (21) is converted into the following Ito^ stochastic differential equation by Ito^ formula:(22)dE1=E1−d1+G1E1,N1,T1dt,dN1=N1r1−β2T1+G3E1,N1,T1+12D1dt+D1N1dB1t,dT1=T1f−β4N1+G2E1,N1,T1+12D2dt+D2T1dB2t,where B1t and B2t are independent Wiener process and 1/2D1E1dt and 1/2D2T1dt are the Wong–Zakai amendments. Defining a random process,(23)RN1,T1=β2T1−β4N1+f−r1+flnβ4N1f−r1lnβ2T1r1.Among them,Rf/β4,r1/β2=0. β2, β4, and r1 are positive, m<0. Then, the Ito^ stochastic differential equation of random process Rt is(24)dR=f−β4N1G3+β2T1−r1G2−12β4N1D1+12β2T1D2dt+D1f−β4N1dB1+D2β2T1−γ1dB2.Using the random average principle of Stratorovich-Khasminskii, the averageIto^ stochastic differential equation of Rt is(25)dR=mRdt+σRdBt.Whence,(26)mR=f−β4N1G3+β2T1−r1G2t+12β2D2T1t−12β4D1N1t,(27)σ2R=D1f−β4N12t+D2β2T1−r12t,and according to the formula,(28)⋅t=1TR∮⋅dt.We have the following expressions:(29)N1t=fβ4,T1t=r1β2,(30)r1f−β4N12t=fr1−β2T12t=gRTR,where(31)gR=r1∮f−β4N1T1dT1.Then, with the help of equations (26)–(30), the drift coefficient and diffusion coefficient can be simplified as follows:(32)mR=12r1D2−12fD1−b1gRTR,σ2R=fD1−r1D2r1fgRTR.


### 3.2.1. Stationary Probability Density

Base on equation (23), we obtain that(33)lnβ4N1f≈1+β4N1−ff−12β4N1−ff2,(34)lnβ4N1r1≈1+β2T1−r1r1−12β2T1−r1r12.Substitute (33) and (34) into (23); then, it can be written as follows:(35)r1β4N1−f2+fβ2T1−r12≈2r1fR.Due to equations,(36)β4N2−f2t≈fR,r1−β2T12t≈r1R,thus, we obtain that(37)mR≈12r1D2−12fD1−b1r1fR,σ2R≈fD1−r1D2R.WhenR⟶4.3782, the boundary is nonsingular; when R⟶∞, it is the second kind of singular boundary. The diffusion index a=1, the drift index β=1,m+∞<0. Therefore, the nontrivial stationary probability density exists. And its steady-state probability density function pr satisfies the following FPK equation:(38)ddrmrpr−12d2dr2σ2rpr=0.Whence,(39)mr=mRR=r,σ2r=σ2RR=r.Solving the FPK equation, we then have(40)pr=CTrgrexp∫0rr1fr1D2−fD1Tx−2r1b1fgxfD1−r1D2gxdx,where C is a normalization constant, and(41)C=C1r1ffD1−r1D2,dgrdr=r1∮β4N1T1β4N1−fdT1.Hence, equation (40) can be simplified as follows:(42)pr=CTrexpυr,υ=−2r1b1ffD1−r1D2.SupposeR=r, the conditional probability density of random process N1 is pN1/r; then, the joint probability density of R and N1 can be expressed as follows:(43)pr1,N1=prN1r1−β2T1Tr.As a result, we obtain the stationary probability density ofN1 and T1:(44)pN1,T1=CN1T1expυr.Thus, we get the steady-state probability densities ofN1 and T1, respectively. They are(45)pN1=∫0∞pN1,T1dT1=∫0∞CN1T1expυrdT1,pT1=∫0∞pN1,T1dN1=∫0∞CN1T1expυrdN1.


### 3.2.2. Mean First-Passage Time

When people are suffering from tumor, we hope that the number of tumor cells go to zero. However, according to the calculation in the above section, we know that the steady-state probability density of tumor cells is present, which indicates that tumor cells will not be extinct under natural state, but the effects of random factors may change this. Therefore, it is practical to study the average first-passage timeβ from the coexistence equilibrium point to the tumor-free equilibrium point, and the Pontryagin equation is satisfied that(46)1+mx0dβdx0+12σ2x0d2βdx02=0,where(47)mx0=mRR=x0,σ2x0=σ2RR=x0.From the above studies, we know thatxl is nonsingular and xs is natural. According to the constant variation method, β can be written as follows:(48)βx0=2∫x0xsψ−1μ∫xlμ1σ2υψυ+C1dμ,where C1 is the integral constant, and(49)ψυ=exp∫2mυσ2υdυ.


## 3.2.1. Stationary Probability Density

Base on equation (23), we obtain that(33)lnβ4N1f≈1+β4N1−ff−12β4N1−ff2,(34)lnβ4N1r1≈1+β2T1−r1r1−12β2T1−r1r12.Substitute (33) and (34) into (23); then, it can be written as follows:(35)r1β4N1−f2+fβ2T1−r12≈2r1fR.Due to equations,(36)β4N2−f2t≈fR,r1−β2T12t≈r1R,thus, we obtain that(37)mR≈12r1D2−12fD1−b1r1fR,σ2R≈fD1−r1D2R.WhenR⟶4.3782, the boundary is nonsingular; when R⟶∞, it is the second kind of singular boundary. The diffusion index a=1, the drift index β=1,m+∞<0. Therefore, the nontrivial stationary probability density exists. And its steady-state probability density function pr satisfies the following FPK equation:(38)ddrmrpr−12d2dr2σ2rpr=0.Whence,(39)mr=mRR=r,σ2r=σ2RR=r.Solving the FPK equation, we then have(40)pr=CTrgrexp∫0rr1fr1D2−fD1Tx−2r1b1fgxfD1−r1D2gxdx,where C is a normalization constant, and(41)C=C1r1ffD1−r1D2,dgrdr=r1∮β4N1T1β4N1−fdT1.Hence, equation (40) can be simplified as follows:(42)pr=CTrexpυr,υ=−2r1b1ffD1−r1D2.SupposeR=r, the conditional probability density of random process N1 is pN1/r; then, the joint probability density of R and N1 can be expressed as follows:(43)pr1,N1=prN1r1−β2T1Tr.As a result, we obtain the stationary probability density ofN1 and T1:(44)pN1,T1=CN1T1expυr.Thus, we get the steady-state probability densities ofN1 and T1, respectively. They are(45)pN1=∫0∞pN1,T1dT1=∫0∞CN1T1expυrdT1,pT1=∫0∞pN1,T1dN1=∫0∞CN1T1expυrdN1.


## 3.2.2. Mean First-Passage Time

When people are suffering from tumor, we hope that the number of tumor cells go to zero. However, according to the calculation in the above section, we know that the steady-state probability density of tumor cells is present, which indicates that tumor cells will not be extinct under natural state, but the effects of random factors may change this. Therefore, it is practical to study the average first-passage timeβ from the coexistence equilibrium point to the tumor-free equilibrium point, and the Pontryagin equation is satisfied that(46)1+mx0dβdx0+12σ2x0d2βdx02=0,where(47)mx0=mRR=x0,σ2x0=σ2RR=x0.From the above studies, we know thatxl is nonsingular and xs is natural. According to the constant variation method, β can be written as follows:(48)βx0=2∫x0xsψ−1μ∫xlμ1σ2υψυ+C1dμ,where C1 is the integral constant, and(49)ψυ=exp∫2mυσ2υdυ.


## 4. Numerical Simulation

Refer to the experimental data of article [20–29] (Table 1) and the analysis of stability of equilibria in the previous section, and we obtain that system (19) of the tumor-free equilibrium is stable, and the corresponding critical state is denoted by xs. Calculate based on Table 1 and reference to the Figure 1, we know the result of coexisting equilibrium belongs to the red area when the m value is −0.8, −0.9, and −1.0. Namely, system (19) has only one stable coexistence equilibrium, and the corresponding critical state is denoted by xl. As a result, xl is nonsingular, xs is natural, and the expression of the first-passage time from xl to xs satisfies equation (46). In order to investigate the effect of the Allee threshold and noise intensity on the first-passage time, we numerically simulate the above results and obtain Figures 3 and 4.Table 1
The definitions and corresponding value of variables in models (2) and (19).
ParameterDefinitionValueSourcesConstant source of effector cells0.33[21, 23]ρMaximum effector cells recruitment rate by tumor cells0.01[20, 21]αSteepness coefficient of the effector cells recruitment0.3[20, 22]β1Effector cells inactivation rate by tumor cells1.0[23, 24]d1Death rate of effector cells0.2[20, 29]r1Normal cells’ growth rate1.0[20, 25]r2Tumor growth rate1.5[26, 27]b2Reciprocal of tumor carrying capacity1.0[20, 28]b1Reciprocal of normal cells’ carrying capacity1.0[20, 28]β2Fractional normal cells killed by tumor cells1.0[23, 24]β3Fractional tumor cells killed by effector cells0.5[23, 24]β4Inactivation rate of tumor cells owing competition with normal cells1.0[20, 23]mThreshold value of Allee effectm<1Figure 3
Mean first-passage time varies withD1. (a) m = −1.0. (b) m = −0.9. (c) m = −0.8.
(a)(b)(c)Figure 4
Mean first-passage time varies withD2. (a) m = −1.0. (b) m = −0.9. (c) m = −0.8.
(a)(b)(c)As can be seen from Figures3 and 4, the time of first passage is limited only when the random intensity is within a certain range. This indicates that the existence of random factors does not cause the first-passage to occur, and only when the random intensity is within a certain range will the system cross from the coexistence equilibrium to the tumor-free equilibrium.


## 5. Conclusions

In this paper, we analyze the property of the deterministic tumor-immune model under the Allee effect, including the positivity and boundness. The results show that, in the condition ofd1>ρ, the solution of system (2) is positive and bounded. At the same time, we know that the Allee effect affects the growth and reproduction of cells in the tumor-immune system according to Figure 2. And the influence is most significant when the Allee threshold of tumor cells is in the range of −1.2 to −0.5.And in the meantime, we are devoted to explore the effect of random factors on the transition of the tumor-immune system from the stable coexistence equilibrium to the stable tumor-free equilibrium. Based on the method of stochastic averaging, we obtain the expressions of the steady-state probability density and the mean first-passage time. Through the numerical simulation about the expression of first-passage time, we know that the existence of random factors will not make the first-passage must happen, in other words, the random strength only within a certain range; the system will pass from the stable coexistence equilibrium point to the tumor-free equilibrium point. Our conclusions provide a possible method to change the number of tumor cells, which is expected to be helpful for the treatment of tumors.



---

*Source: 4892938-2020-09-02.xml*