Dynamics of Tumor-Immune System with Random Noise

Abstract

With deterministic differential equations, we can understand the dynamics of tumor-immune interactions. Cancer-immune interactions can, however, be greatly disrupted by random factors, such as physiological rhythms, environmental factors, and cell-to-cell communication. The present study introduces a stochastic differential model in infectious diseases and immunology of the dynamics of a tumor-immune system with random noise. Stationary ergodic distribution of positive solutions to the system is investigated in which the solution fluctuates around the equilibrium of the deterministic case and causes the disease to persist stochastically. In some conditions, it may be possible to attain infection-free status, where diseases die out exponentially with a probability of one. Some numerical simulations are conducted with the Euler–Maruyama scheme in order to verify the results. White noise intensity is a key factor in treating infectious diseases.

1. Introduction

A cancerous tumor grows with abnormal hormones, causing it to spread throughout the body and potentially spread to other organs or parts of the body. An initial stage is called a tumor, which does not spread. The highest death rate among all diseases worldwide is cancer because of such risk factors as drinking alcohol, smoking, being overweight, physical inactivity, irregular reproductive patterns, an unhealthy lifestyle, monetary development and urbanization. Tumor cells develop and grow in three stages: first, they escape the immune system; then, they reduce the extracellular matrix; and finally, they help blood vessels survive in their tissue by suppressing the extracellular matrix [1]. Cancer can be treated with a variety of methods, such as chemotherapy [2], radiotherapy [3], surgery or hormone therapy [4] to inhibit and eradicate tumor cell growth. Rihan et al. [5,6] explained the stability, optimal control and bifurcation results of the tumor-immune system with the interaction between tumor cells, natural killer cells and cytotoxic T cells, with time delays and the Michaelis–Menten functional response. The author in [7] discusses multiple delays effect on the interaction between brain tumors and immune components in conjunction with T11 (CD2) target administration as well as the conditions for the local stability of steady states, bifurcation results, and the maximum delay to keep limit cycles stable. Many other researchers have investigated tumor-immune models to understand how the immune system works and examine how immune cells eliminate tumors in different ways [8,9].

Recent studies incorporated different random noises into deterministic tumor-immune models in order to discuss the dynamical behaviors of the corresponding stochastic models (see [10,11]). Basically, the environmental fluctuations stimulate mainly the intrinsic growth rate of host, tumor cells and death rate of effector cells. Moreover, the modeling of such phenomena, stochastic differential equations (SDEs), are more suitable than deterministic models, which offer a more reasonable representation for discussing the long-term kinetics of cell population. Liu et al. [12] studied the dynamical behaviors of tumor-immune responses under chemotherapy treatment; deterministic and stochastic differential equation models were constructed to characterize the dynamical changes in tumor and immune cells. The deterministic model was extended to the stochastic differential equations (SDEs) model and the continuous-time Markov chain (CTMC) model, which accounts for the variability in cellular reproduction, interspecific competition, growth, death, immune response, and chemotherapy. Yang et al. [13] derived the global positive solution and qualitative behaviors of the tumor-immune model with the combination of pulsed immunotherapy, pulsed chemotherapy and white noise effect. Das et al. [14] investigated the deterministic and stochastic modeling of the tumor-immune system under Michaelis–Menten kinetics and also studied the stochastic permanence, global attractivity and weak persistence in mean. The authors in [15] discussed the threshold condition about immune strength for survival, extinction and weak persistence results of a stochastic tumor-immune system.

In this paper, white noise is incorporated into an existing deterministic tumor-immune model to analyze the dynamics of the system. The presence and uniqueness of the global non-negative solution of the stochastic tumor-immune model with a Holling type III functional response is investigated. Using a stochastic Lyapunov function combined with Ito’s formula, we provide a sufficient condition for determining the existing results of stationary distribution, weak persistence, and extinction of tumor cells. The rest of this paper is organized as follows: In Section 2, we formulate the tumor-immune model and study the existence of global positive solution. The stationary distribution and extinction results of this model are derived in Section 3 and Section 4. Some numerical simulations are given in Section 5 to verify the obtained theoretical results. Section 6 contains the conclusion.

2. Stochastic Model for Tumor-Immune Interaction

It worth mentioning here that deterministic models are assumed for tumor-immune interactions; however, there is increasing evidence that better consistency with some phenomena can be provided if the effects of random processes in the system are taken into account. One of the important facts about the impact of the environmental noise is that it can suppress a potential population explosion [16,17,18].

The interaction between cancer and the immune system (IS) has been investigated by many authors using deterministic mathematical models (see [19,20,21,22,23,24,25]). The challenge is to obtain the known biological features without making the mathematics too complex. We include here the following features: The immune system is a multifunctional multipathway. We, therefore, use two immune effectors. The cross-reactivity of the immune system is also taken into account. The model is given by the following:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \dot{T}\left(t\right)=aT-r_{1}T{E}_1-r_{2}T{E}_2,}\hfill \\ {\displaystyle {\dot{E}}_1\left(t\right)=-d_1{E}_1+\frac{T^2{E}_1}{T^2+k_1},}\hfill \\ {\displaystyle {\dot{E}}_2\left(t\right)=-d_2{E}_2+\frac{T^2{E}_2}{T^2+k_2},}\hfill \end{array}\end{array}$$

(1)

where $T\equiv T\left(t\right)$ is the tumor cells; $E_1\equiv E_1\left(t\right)$, $E_2\equiv E_2\left(t\right)$ are immune effectors with initial values of $T\left(0\right)>0,E_1\left(0\right)>0,E_2\left(0\right)>0$; a is the growth rate of $T\left(t\right)$; $d_1$ and $d_2$ are the decay rate of $E_1\left(t\right)$ and $E_2\left(t\right)$, respectively; $k_1,k_2$ represent the half saturation parameters; and $r_1$ and $r_2$ are the suppress rate of $T\left(t\right)$ due to presence of $E_1\left(t\right)$ and $E_2\left(t\right)$, respectively. All the above parameters are positive constants. The interaction terms in the second and third equations of system (1) satisfy the cross reactivity property of the immune system.

Lemma 1.

The solution of model (1) with non-negative initial values is non-negative.

Proof. 

The right hand sides of model (1) are continuous functions of dependent variables, and we can easily obtain the following:

$$\begin{array}{cc}\hfill T\left(t\right)& =T\left(0\right)exp\left({\int }_0^t\left[a-r_1{E}_1\left(s\right)-r_2{E}_2\left(s\right)\right]ds\right),\hfill \\ \hfill E_1\left(t\right)& =E_1\left(0\right)exp\left({\int }_0^t\left[-d_1+\frac{T^2\left(s\right)}{T^2\left(s\right)+k_1}\right]ds\right),\hfill \\ \hfill E_2\left(t\right)& =E_2\left(0\right)exp\left({\int }_0^t\left[-d_2+\frac{T^2\left(s\right)}{T^2\left(s\right)+k_2}\right]ds\right).\hfill \end{array}$$

Therefore, it is clear from the above expressions that $T\left(t\right),E_1\left(t\right)$ and $E_2\left(t\right)$ remain positive for all $t>0$ if they initiate from an interior point of the following:

$${\mathbb{R}}_{+}^3=\{T\left(t\right),E_1\left(t\right),E_2\left(t\right):T\left(t\right)>0,E_1\left(t\right)>0,E_2\left(t\right)>0\}.$$

Thus, ${\mathbb{R}}_{+}^3$ is positively invariant for (1). □

The steady states of the system (1) are the following:

$$\begin{array}{c}\hfill (0,0,0);(\sqrt{d_1{k}_1/(1-d_1)},a/r_1,0);(\sqrt{d_2{k}_2/(1-d_2)},0,a/r_2).\end{array}$$

(2)

We assume that in order to avoid the non-biological interior solution in which both immune effectors coexist, we have the following:

$$\begin{array}{c}\hfill (d_1{k}_1/(1-d_1))<<(d_2{k}_2/(1-d_2)).\end{array}$$

(3)

According to the size of the tumor, the first steady state is naive, the second is memory, and the third is endemic. Analyzing the stability of naive state shows that they are unstable. Locally, the memory state is asymptotically stable if the following holds:

$$\begin{array}{c}\hfill d_1<1, \mathrm{and} d_1<d_2.\end{array}$$

(4)

As long as the endemic state remains locally stable, the following holds:

$$\begin{array}{c}\hfill d_2<1, \mathrm{and} d_2<d_1.\end{array}$$

(5)

Hence, we arrive at the following propositions:

Proposition 1.

When (4) is satisfied, then the memory state is asymptotically stable on the local level. Therefore, there is a bifurcation when $d_1=1$.

Growing and interacting cancer cells differ depending on the cell properties, the immune system’s functioning, and the environment where the interaction occurs [26]. An individual’s overall health determines the environment in which the body operates. In order to discuss the impact of body environmental factors on the dynamics of cancer infection, we could extend the deterministic description of the tumor-immune interaction to include stochastic forcing, either additively or multiplicatively. Several researchers have studied mathematical models for studying epidemics in environmental noise, such as [27].

Some authors have introduced randomness into deterministic models to demonstrate the influence of environmental variation on biological systems by taking this into consideration; see [28,29]. For more realistic situation of the development process of the disease, we incorporate the effect of randomization within host by introducing nonlinear perturbation on the natural death rate with white noise into each equation of model (1). The parameters associated with the tumor-immune system are not certain, but the interval to which it belongs can easily be determined. We therefore suggest the following stochastic model:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle dT\left(t\right)=(aT-r_{1}T{E}_1-r_{2}T{E}_2)dt+{\sigma }_{1}T\left(t\right)d{W}_1\left(t\right),}\hfill \\ {\displaystyle d{E}_1\left(t\right)=(-d_1{E}_1+\frac{T^2{E}_1}{T^2+k_1})dt+{\sigma }_2{E}_1\left(t\right)d{W}_2\left(t\right),}\hfill \\ {\displaystyle d{E}_2\left(t\right)=(-d_2{E}_2+\frac{T^2{E}_2}{T^2+k_2})dt+{\sigma }_3{E}_2\left(t\right)d{W}_3\left(t\right),}\hfill \end{array}\end{array}$$

(6)

where ${\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2$ are intensities of the environmental white noises. $W_1\left(t\right),W_2\left(t\right),W_3\left(t\right)$ are mutually independent standard Brownian motions with $W_i\left(0\right)=0$ ($i=1,2,3$). We define the basic concepts of probability theory and SDEs. Let $(\Omega ,\mathfrak{F},{\left\{{\mathfrak{F}}_t\right\}}_{t\ge 0},P)$ be complete probability space with filtration ${\left\{{\mathfrak{F}}_t\right\}}_{t\ge 0}$ satisfying the usual conditions. See more details about Ito’s formula (see [16,30,31]).

Let $y\left(t\right)$ be a regular time-homogeneous Markov process in ${\mathbb{R}}^n$ defined by SDE:

$$\begin{array}{c}\hfill dy\left(t\right)=f\left(y\right(t\left)\right)dt+g\left(y\right(t\left)\right)dW\left(t\right).\end{array}$$

(7)

The diffusion matrix of the process $y\left(t\right)$ is described as $\mathcal{A}\left(y\right)=\left(b_{ij}\left(y\right)\right)$, $b_{ij}\left(y\right)=g^i\left(y\right)g^j\left(y\right).$

Lemma 2

(Ref. [32]). The Markov process $y\left(t\right)$ has a unique ergodic stationary distribution $\pi (·)$ if there exists a bounded open domain $\mathcal{D}\subset {\mathbb{R}}^n$ with regular boundary $\Gamma ,$ having the following properties:

i.

The diffusion matrix $\mathcal{A}\left(y\right)$ is strictly positive definite for all $y\in \mathcal{D}.$

ii.

There exists a non-negative $C^2-$ function $\mathcal{V}$ such that $L\mathcal{V}$ is negative for any ${\mathbb{R}}^n\setminus \mathcal{D}.$

Theorem 1.

Model (6) has a unique positive solution $(T\left(t\right),E_1\left(t\right),E_2\left(t\right))$ on $t\ge 0$ with$(T\left(0\right),E_1\left(0\right),E_2\left(0\right))\in R_{+}^3$, and the solution remains in ${\mathbb{R}}_{+}^3$ with probability 1.

Proof. 

Basically, the coefficients of system (6) satisfy the local Lipschitz condition. Then, (6) has a unique local solution $(T\left(t\right),E_1\left(t\right),E_2\left(t\right))$ on $[0,{\tau }_e]$, where ${\tau }_e$ is an exposure time. Then, we prove that ${\tau }_e=+\infty .$ Let us follow the similar proof of Theorem 3.1 in [16]. The major step is to describe a non-negative $C^2$ function $V:{\mathbb{R}}_{+}^3\to {\mathbb{R}}_{+}$ such that ${lim}_{h\to \infty ,(T,E_1,E_2)\in {\mathbb{R}}_{+}^3\setminus D_h}infV(T,E_1,E_2)=+\infty ,andLV(T,E_1,E_2)\le \mathcal{K},$ where $D_h=(\frac{1}{h},h)\times (\frac{1}{h},h)\times (\frac{1}{h},h),$ and $\mathcal{K}$ is a non-negative constant. Define a function $V:{\mathbb{R}}_{+}^3\to {\mathbb{R}}_{+}$ as follows

$$V(T,E_1,E_2)=(T-1-lnT)+(E_1-1-ln{E}_1)+(E_2-1-ln{E}_2).$$

The non-negativity of this function can be seen from $\varsigma -1-ln\varsigma \ge 0$, $\forall \varsigma >0.$ By applying the Ito’s formula, we can obtain the following:

$$\begin{array}{cc}\hfill dV(T,E_1,E_2)=& LV(T,E_1,E_2)dt+{\sigma }_1(T-1)d{W}_1\left(t\right)\hfill \\ & +{\sigma }_2(E_1-1)d{W}_2\left(t\right)+{\sigma }_3(E_2-1)d{W}_3\left(t\right),\hfill \\ \hfill LV(T,E_1,E_2)=& (1-\frac{1}{T})(aT-r_{1}T{E}_1-r_{2}T{E}_2)+(1-\frac{1}{E_1})(-d_1{E}_1+\frac{T^2{E}_1}{T^2+k_1})\hfill \\ & +(1-\frac{1}{E_2})(-d_2{E}_2+\frac{T^2{E}_2}{T^2+k_2})+\frac{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2}{2},\hfill \\ \hfill \le & aT+r_1{E}_1+r_2{E}_2+T^2(E_1+E_2)+d_1+d_2+\frac{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2}{2},\hfill \end{array}$$

Applying the superior of the co-efficient of above inequality and using the positiveness of $T,E_1$ and $E_2$, there exists a positive constant $\mathcal{K}$ such that $LV\le \mathcal{K}.$ The rest of the proof follows that of [16,31] and hence, it is omitted. □

3. Existence of Ergodic Stationary Distribution

Here, we discuss the stationary distribution and extinction results for model (6), which helps to identify whether the disease is persistent or can be eradicated.

Theorem 2.

If $d-\frac{{\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2}{2}>0,$ then model (6) has a unique ergodic stationary distribution for any $(T\left(0\right),E_1\left(0\right),E_2\left(0\right))\in {\mathbb{R}}_{+}^3$.

Proof. 

The diffusion matrix of (6) is calculated as follows:

$$\mathcal{A}=\left(\begin{array}{ccc}{\sigma }_1^2{T}^2& 0& 0\\ 0& {\sigma }_2^2{E}_1^2& 0\\ 0& 0& {\sigma }_3^2{E}_2^2\end{array}\right),$$

which is positive definite for any compact subset of ${\mathbb{R}}_{+}^3$. Condition (i) in Lemma 2 is verified.

Define the $C^2-$ function $V:{\mathbb{R}}_{+}^3\to \mathbb{R}$ as follows:

$$\begin{array}{cc}\hfill V(T,E_1,E_2)& =-lnT-ln{E}_1-ln{E}_2+T^{-\theta }+E_1^{-\theta }+\frac{1}{\vartheta +1}{(T+E_1+E_2)}^{\vartheta +1},\hfill \\ \hfill & =V_1+V_2+V_3.\hfill \end{array}$$

Further, $V(T,E_1,E_2)$ is not only continuous, but also goes to $+\infty$, as $(T,E_1,E_2)$ tends to the boundary of ${\mathbb{R}}_{+}^3$ and $\parallel (T,E_1,E_2)\parallel \to \infty .$ V must have a minimum point $(T\left(0\right),E_1\left(0\right),E_2\left(0\right))$ in the interior of ${\mathbb{R}}_{+}^3.$ Define a $C^2-$function $\tilde{V}:{\mathbb{R}}_{+}^3\to {\mathbb{R}}_{+}$ as

$$\tilde{V}(T,E_1,E_2)=V(T,E_1,E_2)-V(T\left(0\right),E_1\left(0\right),E_2\left(0\right)).$$

Let $V_1=-lnT-ln{E}_1-ln{E}_2,V_2=T^{-\theta }+E_1^{-\theta },V_3=\frac{1}{\vartheta +1}{(T+E_1+E_2)}^{\vartheta +1}$, $\theta >0$, $d=max\{d_1,d_2\}$ and constant $\vartheta >1$ satisfies $d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)>0,$ and the following:

$$\begin{array}{cc}\hfill M_1& =\underset{(T,E_1,E_2)\in R_{+}^3}{sup}\{((a+d)T+(E_1+E_2)T^2){(T+E_1+E_2)}^{\vartheta }-\frac{1}{2}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)]\hfill \\ \hfill & \times (T^{\vartheta +1}+E_1^{\vartheta +1}+E_2^{\vartheta +1})\}<\infty ,\hfill \\ \hfill M_2& =\underset{(T,E_1,E_2)\in R_{+}^3}{sup}\{\theta T^{-\theta }{r}_1{E}_1+M_1-\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)](T^{\vartheta +1}+E_1^{\vartheta +1}+E_2^{\vartheta +1})\hfill \\ & +\theta E_1^{-\theta }(d_1+\frac{(\theta +1)}{2}{\sigma }_2^2)+r_1{E}_1+r_2{E}_2+H\}.\hfill \end{array}$$

Applying Ito’s formula L on the functions $V_1,V_2$ and $V_3$ give the following:

$$\begin{array}{cc}\hfill L{V}_1& \le r_1{E}_1+r_2{E}_2+d_1+d_2+\frac{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2}{2},\hfill \end{array}$$

(8)

$$\begin{array}{cc}L{V}_2\hfill & \le -\theta T^{-\theta }(a-r_1{E}_1-r_2{E}_2)+\frac{\theta (\theta +1)}{2}{\sigma }_1^2{T}^{-\theta }-\theta E_1^{-\theta }(-d_1+\frac{T^2}{T^2+k_1})+\frac{\theta (\theta +1)}{2}{\sigma }_2^2{E}_1^{-\theta },\hfill \\ \hfill & \le -\theta T^{-\theta }(a-\frac{(\theta +1)}{2}{\sigma }_1^2)-\theta E_1^{-\theta }{T}^2+\theta T^{-\theta }(r_1{E}_1+r_2{E}_2)+\theta E_1^{-\theta }(d_1+\frac{(\theta +1)}{2}{\sigma }_2^2),\hfill \end{array}$$

(9)

$$\begin{array}{cc}L{V}_3\hfill & \le {(T+E_1+E_2)}^{\vartheta }\left\{(a+d)T+(E_1+E_2)T^2-d(T+E_1+E_2)\right\}+\frac{\vartheta }{2}{(T+E_1+E_2)}^{\vartheta -1}\hfill \\ & \times [{\sigma }_1^2{T}^2+{\sigma }_2^2{E}_1^2+{\sigma }_3^2{E}_2^2],\hfill \\ \hfill & \le {(T+E_1+E_2)}^{\vartheta }\left\{(a+d)T+(E_1+E_2)T^2-d(T+E_1+E_2)\right\}+\frac{\vartheta }{2}{(T+E_1+E_2)}^{\vartheta +1}\hfill \\ & \times ({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2),\hfill \\ \hfill & \le \left\{(a+d)T+(E_1+E_2)T^2\right\}{(T+E_1+E_2)}^{\vartheta }-{(T+E_1+E_2)}^{\vartheta +1}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)],\hfill \\ \hfill & \le M_1-\frac{1}{2}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)](T^{\vartheta +1}+E_1^{\vartheta +1}+E_2^{\vartheta +1}).\hfill \end{array}$$

(10)

From Equations (8)–(10), we have the following:

$$\begin{array}{cc}L\tilde{V}\hfill & \le -\theta T^{-\theta }(a-\frac{(\theta +1)}{2}{\sigma }_1^2)-\theta E_1^{-\theta }{T}^2+\theta T^{-\theta }(r_1{E}_1+r_2{E}_2)+\theta E_1^{-\theta }(d_1+\frac{(\theta +1)}{2}{\sigma }_2^2)+r_1{E}_1\hfill \\ & +r_2{E}_2+d_1+d_2+\frac{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2}{2}+M_1-\frac{1}{2}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)](T^{\vartheta +1}+E_1^{\vartheta +1}+E_2^{\vartheta +1}),\hfill \\ \hfill & \le -\theta T^{-\theta }(a-\frac{(\theta +1)}{2}{\sigma }_1^2)-\theta E_1^{-\theta }{T}^2+\theta T^{-\theta }{r}_2{E}_2-\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)](T^{\vartheta +1}\hfill \\ \hfill & +E_1^{\vartheta +1}+E_2^{\vartheta +1})+\theta T^{-\theta }{r}_1{E}_1+M_1-\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)](T^{\vartheta +1}+E_1^{\vartheta +1}+E_2^{\vartheta +1})\hfill \\ \hfill & +\theta E_1^{-\theta }(d_1+\frac{(\theta +1)}{2}{\sigma }_2^2)+r_1{E}_1+r_2{E}_2+H.\hfill \end{array}$$

Define a bounded closed set as follows:

$$D=\left\{(T,E_1,E_2)\in {\mathbb{R}}_{+}^3:ϵ\le T\le \frac{1}{ϵ},{ϵ}^2\le E_1\le \frac{1}{{ϵ}^2},{ϵ}^2\le E_2\le \frac{1}{{ϵ}^2}\right\}.$$

In the set ${\mathbb{R}}_{+}^3\setminus D,$ let us choose $ϵ>0$ satisfies the following:

$$\begin{array}{cc}\hfill & -\frac{\theta }{{ϵ}^{\theta }}(a-\frac{(\theta +1)}{2}{\sigma }_1^2)+M_2\le -1,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & -\frac{\theta {ϵ}^2}{{ϵ}^{2\theta }}+M_2\le -1,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & 0<ϵ<\frac{1}{r_2},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & -\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)]\frac{1}{{ϵ}^{\vartheta +1}}+M_2\le -1,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & -\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)]\frac{1}{{ϵ}^{2(\vartheta +1)}}+M_2\le -1,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & -\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)]\frac{1}{{ϵ}^{2(\vartheta +1)}}+M_2\le -1,\hfill \end{array}$$

(16)

where $M_2={sup}_{(T,E_1,E_2)\in {\mathbb{R}}_{+}^3}\left\{\theta T^{-\theta }{r}_1{E}_1+M_1-\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)](T^{\vartheta +1}+E_1^{\vartheta +1}+E_2^{\vartheta +1})+\theta E_1^{-\theta }(d_1+\frac{(\theta +1)}{2}{\sigma }_2^2)+r_1{E}_1+r_2{E}_2+H\right\}$ and $H=d_1+d_2+\frac{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2}{2}.$

Let us prove condition (ii) of Lemma 2 to show that $L\tilde{V}\le -1$ for $(T,E_1,E_2)\in {\mathbb{R}}_{+}^3\setminus D$ and ${\mathbb{R}}_{+}^3\setminus D={\bigcup }_{i=1}^6{D}_i,$ where the following holds:

$$\begin{array}{cc}\hfill D_1& =\{(T,E_1,E_2)\in {\mathbb{R}}_{+}^3;0<T<ϵ\},\hfill \\ \hfill D_2& =\{(T,E_1,E_2)\in {\mathbb{R}}_{+}^3;0<E_1<{ϵ}^2,T\ge ϵ\},\hfill \\ \hfill D_3& =\{(T,E_1,E_2)\in {\mathbb{R}}_{+}^3;0<E_2<{ϵ}^2,E_1\ge {ϵ}^2,T\ge ϵ\},\hfill \\ \hfill D_4& =\{(T,E_1,E_2)\in {\mathbb{R}}_{+}^3;T>\frac{1}{ϵ}\},\hfill \\ \hfill D_5& =\{(T,E_1,E_2)\in {\mathbb{R}}_{+}^3;E_1>\frac{1}{{ϵ}^2}\},\hfill \\ \hfill D_6& =\{(T,E_1,E_2)\in {\mathbb{R}}_{+}^3;E_2>\frac{1}{{ϵ}^2}\}.\hfill \end{array}$$

Case (i) For $(T,E_1,E_2)\in D_1$, we obtain the following:

$$\begin{array}{cc}L\tilde{V}\hfill & \le -\theta T^{-\theta }(a-\frac{(\theta +1)}{2}{\sigma }_1^2)+\theta T^{-\theta }{r}_1{E}_1+M_1-\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)](T^{\vartheta +1}+E_1^{\vartheta +1}+E_2^{\vartheta +1})\hfill \\ & +\theta E_1^{-\theta }(d_1+\frac{(\theta +1)}{2}{\sigma }_2^2)+r_1{E}_1+r_2{E}_2+H,\hfill \\ \hfill & \le -\theta T^{-\theta }(a-\frac{(\theta +1)}{2}{\sigma }_1^2)+M_2,\hfill \\ \hfill & \le -\frac{\theta }{{ϵ}^{\theta }}(a-\frac{(\theta +1)}{2}{\sigma }_1^2)+M_2\le -1,\hfill \end{array}$$

which is obtained from (11).

Case (ii). For $(T,E_1,E_2)\in D_2,$ we have the following:

$$\begin{array}{cc}\hfill L\tilde{V}& \le -\theta E_1^{-\theta }{T}^2+M_2,\hfill \\ \hfill & \le -\theta \frac{{ϵ}^2}{{ϵ}^{2\theta }}+M_2\le -1,\hfill \end{array}$$

which is obtained from (12).

Case (iii). For $(T,E_1,E_2)\in D_3,$ we have the following:

$$\begin{array}{cc}\hfill L\tilde{V}& \le -\theta E_1^{-\theta }{T}^2+r_2{E}_2+M_2,\hfill \\ \hfill & \le -\theta \frac{{ϵ}^2}{{ϵ}^{2\theta }}+r_2{ϵ}^2+M_2<-1,\hfill \end{array}$$

which is obtained from (13) and $-\theta \frac{{ϵ}^2}{{ϵ}^{2\theta }}+M_2\le -2.$

Case (iv). For $(T,E_1,E_2)\in D_4,$ we have the following:

$$\begin{array}{cc}\hfill L\tilde{V}& \le -\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)]T^{\vartheta +1}+M_2,\hfill \\ \hfill & \le -\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)]\frac{1}{{ϵ}^{(\vartheta +1)}}+M_2\le -1,\hfill \end{array}$$

which is obtained from (14).

Case (v). For $(T,E_1,E_2)\in D_5,$ we obtain the following:

$$\begin{array}{cc}\hfill L\tilde{V}& \le -\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)]E_1^{\vartheta +1}+M_2,\hfill \\ \hfill & \le -\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)]\frac{1}{{ϵ}^{2(\vartheta +1)}}+M_2\le -1,\hfill \end{array}$$

which is obtained from (15).

Case (vi). For $(T,E_1,E_2)\in D_6,$ it yields the following:

$$\begin{array}{cc}\hfill L\tilde{V}& \le -\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)]E_2^{\vartheta +1}+M_2,\hfill \\ \hfill & \le -\frac{1}{4}[d-\frac{\vartheta }{2}({\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2)]\frac{1}{{ϵ}^{2(\vartheta +1)}}+M_2\le -1,\hfill \end{array}$$

which is obtained from (16).

Therefore, following the above discussion, there exists a $ϵ>0,$ such that $L\tilde{V}(T,E_1,E_2)\le -1,$ for all $(T,E_1,E_2)\in {\mathbb{R}}_{+}^3\setminus D.$ Based on Lemma 2, the model (6) has unique ergodic stationary distribution. □

4. Extinction

Theorem 3.

Let $(T\left(t\right),E_1\left(t\right),E_2\left(t\right))$ be the solution of (6) with $(T\left(0\right),E_1\left(0\right),E_2\left(0\right))\in {\mathbb{R}}_{+}^3$. If $a<\frac{{\sigma }_1^2}{2}$, then the tumor cell $T\left(t\right)$ populations will die out, i.e., ${lim}_{t\to \infty }T\left(t\right)=0.$

Proof. 

Applying Ito’s formula to the first equation of (6), one can obtain the following:

$$\begin{array}{c}\hfill d\left(lnT\left(t\right)\right)=(a-r_1{E}_1-r_2{E}_2-\frac{{\sigma }_1^2}{2})dt+{\sigma }_{1}d{W}_1\left(t\right).\end{array}$$

Taking integration from 0 to t on both sides and dividing by t, we have the following:

$$\begin{array}{cc}\hfill \frac{lnT\left(t\right)-lnT\left(0\right)}{t}& =a-\frac{r_1}{t}{\int }_0^t{E}_1\left(s\right)ds-\frac{r_2}{t}{\int }_0^t{E}_2\left(s\right)ds-\frac{{\sigma }_1^2}{2}+\frac{{\sigma }_1{W}_1\left(t\right)}{t},\hfill \\ \hfill & \le a-\frac{{\sigma }_1^2}{2}+\frac{{\sigma }_1{W}_1\left(t\right)}{t}.\hfill \end{array}$$

By using the strong law of large numbers for local martingales, ${lim}_{t\to \infty }\frac{W_1\left(t\right)}{t}=0,a.s$

$$\begin{array}{cc}\hfill \underset{t\to \infty }{lim}\sup \frac{lnT\left(t\right)}{t}& \le a-\frac{{\sigma }_1^2}{2}<0,\hfill \\ \hfill \mathrm{Moreover},& \underset{t\to \infty }{lim}T\left(t\right)=0.\hfill \end{array}$$

Defining $ln(E_1\left(t\right)+E_2\left(t\right))$ and applying Ito’s formula, we obtain the following:

$$\begin{array}{cc}d\left(ln(E_1\left(t\right)+E_2\left(t\right))\right)\hfill & =\frac{1}{E_1\left(t\right)+E_2\left(t\right)}\left\{-d_1{E}_1\left(t\right)-d_2{E}_2\left(t\right)+\frac{T^2\left(t\right)E_1\left(t\right)}{T^2\left(t\right)+k_1}++\frac{T^2\left(t\right)E_2\left(t\right)}{T^2\left(t\right)+k_2}\right\}dt\hfill \\ \hfill & -\frac{{\sigma }_2^2{E}_1^2\left(t\right)+{\sigma }_3^2{E}_2^2\left(t\right)}{2{(E_1\left(t\right)+E_2\left(t\right))}^2}dt+\frac{{\sigma }_2{E}_1\left(t\right)}{E_1\left(t\right)+E_2\left(t\right)}d{W}_2\left(t\right)+\frac{{\sigma }_3{E}_2\left(t\right)}{E_1\left(t\right)+E_2\left(t\right)}d{W}_3\left(t\right).\hfill \end{array}$$

Based on ${lim}_{t\to \infty }T\left(t\right)=0$, there exists $t_1>0$ such that $T\left(t\right)<ϵ$ when $t>t_1$ and $k=min\{k_1,k_2\}$ and $d_0=min\{d_1,d_2\}$.

$$\begin{array}{cc}\hfill d\left(ln(E_1\left(t\right)+E_2\left(t\right))\right)& \le \left\{\frac{1}{k}-d_0\right\}dt+\frac{{\sigma }_2{E}_1\left(t\right)}{E_1\left(t\right)+E_2\left(t\right)}d{W}_2\left(t\right)+\frac{{\sigma }_3{E}_2\left(t\right)}{E_1\left(t\right)+E_2\left(t\right)}d{W}_3\left(t\right).\hfill \end{array}$$

Let $P_1\left(t\right)={\int }_0^t\frac{{\sigma }_2{E}_1\left(t\right)}{E_1\left(t\right)+E_2\left(t\right)}d{W}_2\left(s\right)$ and $P_2\left(t\right)={\int }_0^t\frac{{\sigma }_3{E}_2\left(t\right)}{E_1\left(t\right)+E_2\left(t\right)}d{W}_3\left(s\right)$ be local martingales with quadratic variations as follows:

$$\begin{array}{c}\hfill <P_1\left(t\right),P_1\left(t\right){>}_t={\sigma }_2^2{\int }_0^t{\left(\frac{E_1\left(t\right)}{E_1\left(t\right)+E_2\left(t\right)}\right)}^{2}ds\le {\sigma }_2^{2}t,\\ \hfill <P_2\left(t\right),P_2\left(t\right){>}_t={\sigma }_3^2{\int }_0^t{\left(\frac{E_2\left(t\right)}{E_1\left(t\right)+E_2\left(t\right)}\right)}^{2}ds\le {\sigma }_3^{2}t.\end{array}$$

Using the strong law of large numbers for the local martingales, ${lim}_{t\to \infty }\frac{P_1\left(t\right)}{t}=0$, ${lim}_{t\to \infty }\frac{P_2\left(t\right)}{t}=0,a.s.$ Taking integration from 0 to t on both sides and dividing by t, we have the following:

$$\begin{array}{cc}\hfill \frac{ln(E_1\left(t\right)+E_2\left(t\right))-ln(E_1\left(0\right)+E_2\left(0\right))}{t}& \le \left\{\frac{1}{k}-d_0\right\}+\frac{1}{t}{\int }_0^t\frac{{\sigma }_2{E}_1\left(s\right)}{E_1\left(s\right)+E_2\left(s\right)}d{W}_2\left(s\right)\hfill \\ & +\frac{1}{t}{\int }_0^t\frac{{\sigma }_3{E}_2\left(s\right)}{E_1\left(s\right)+E_2\left(s\right)}d{W}_3\left(s\right),\hfill \\ \hfill & \le \left\{\frac{1}{k}-d_0\right\}+\frac{P_1\left(t\right)}{t}+\frac{P_2\left(t\right)}{t},\hfill \\ \hfill \underset{t\to \infty }{lim}\sup \frac{ln(E_1\left(t\right)+E_2\left(t\right))}{t}& \le \left\{\frac{1}{k}-d_0\right\}.\hfill \end{array}$$

□

We arrive at the following remarks:

Remark 1.

If $a>\frac{{\sigma }_1^2}{2}$ and $\left\{\frac{1}{k}-d_0\right\}<0$, we can obtain results, such as ${lim}_{t\to \infty }SupT\left(t\right)>0$, ${lim}_{t\to \infty }{E}_1\left(t\right)=0$ and ${lim}_{t\to \infty }{E}_2\left(t\right)=0.$ Clearly, the tumor cells $T\left(t\right)$ are weakly persistent in the mean a.s.

Remark 2.

Theorem 2 shows that under small white noises, the tumor cell $T\left(t\right)$ and effector cells $E_1\left(t\right)$ and $E_2\left(t\right)$ distribution approaches to an invariant measure as $t\to \infty$. That is, the tumor cell $T\left(t\right)$ tends to a dormant steady state, stochastic in nature.

Remark 3.

Theorem 3 shows that when the stochastic perturbation for tumor cells $T\left(t\right)$ is strong enough, the tumor goes to extinction, while the effector cells $E_1\left(t\right)$ and $E_2\left(t\right)$ distribution converges to a steady state $\frac{1}{k}-d_0$. We can easily see that ${\sigma }_1$ is a critical parameter to eradicate the tumor cells $T\left(t\right)$, and the effector cells $E_1\left(t\right)$ and $E_2\left(t\right)$ approach a steady state stochastic in nature.

5. Numerical Simulations

In this section, we use Euler–Maruyama method for solving SDEs discussed in detail in Refs. [16,33], to obtain the discretization transformation of (6) as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill T_{j+1}=& T_j+[a{T}_j-r_1{T}_j{E}_{1,j}-r_2{T}_j{E}_{2,j}]\Delta t+{\sigma }_1{T}_j\sqrt{\Delta t}{\zeta }_{1,j},\hfill \\ \hfill E_{1,j+1}=& E_{1,j}+[-d_1{E}_{1,j}+\frac{T_j^2{E}_{1,j}}{T_j^2+k_1}]\Delta t+{\sigma }_2{E}_{1,j}\sqrt{\Delta t}{\zeta }_{2,j},\hfill \\ \hfill E_{2,j+1}=& E_{2,j}+[-d_2{E}_{2,j}+\frac{T_j^2{E}_{2,j}}{T_j^2+k_2}]\Delta t+{\sigma }_3{E}_{2,j}\sqrt{\Delta t}{\zeta }_{3,j},\hfill \end{array}\end{array}$$

(17)

where $\Delta t>0$ is the time increment, ${\zeta }_{i,j},(i=1,2,3)$ are independent Gaussian random variables, which follow the distribution $N(0,1)$.

We assign the following parameter values: $a=1,$$r_1=1$, $r_2=1$, $k_1=0.3$, $k_2=0.7$, $d_1=0.3$, $d_2=0.7$.

The deterministic model (1), with long-range temporal memory has endemic steady states, which are locally asymptotically stable, according to conditions (4) and (5): $d_1=0.3<1$, $d_2=0.7$ and $d_2=0.3<1,d_1=0.7$, respectively. Figure 1 and Figure 2 show the numerical simulations of the model with stable memory and endemic steady states. With two cases of parameter values $d_1=d_2=0.92<1$ and $d_1=d_2=1.02>1$, the stability and instability conditions of model (1) are presented in Figure 3.

Now, we incorporate white noise in the model to show the dynamics of the stochastic model (6). Initially, we consider the white noise values ${\sigma }_1=0.5,{\sigma }_2=0.6,{\sigma }_3=0.8$, $k_1=0.3$, $k_2=0.7,d_1=0.3,d_2=0.7$ (right banner of Figure 4) and $k_1=0.7,k_2=0.3,d_1=0.7$, $d_2=0.3$ (left banner of Figure 4); the threshold conditions of unique stationary distribution $d-\frac{{\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2}{2}=0.075>0$ are satisfied. The left and right banners of Figure 4 show that the tumor cells $T\left(t\right)$, effector cells $E_1\left(t\right)$ and $E_2\left(t\right)$ fluctuate randomly. We then slightly increase the white noise values: ${\sigma }_1=0.9,{\sigma }_2=1.1,{\sigma }_3=1.3,k_1=0.3,k_2=0.7$, $d_1=0.3,d_2=0.7$ (right banner of Figure 5) and $k_1=0.7$, $k_2=0.3$, $d_1=0.7$, $d_2=0.3$ (left banner of Figure 5). When the condition of weak persistence $a-\frac{{\sigma }_1^2}{2}=0.595>0$ is satisfied, we can see in the left and right banners of Figure 5 a weak persistence in the mean of $T\left(t\right)$. Tumor cell load $T\left(t\right)$ gradually decreases and fluctuates in the neighborhood of zero, defining weak persistence in the mean. The intensities of white noise can reduce the certain degree of cancer cells and suppress tumor growth, but not completely eliminate the cancer cells. Thus, the mutation and diffusion of tumor cells can be controlled by varying the strength of noise.

By increasing the white noise values ${\sigma }_1=1.5,{\sigma }_2=1.6,{\sigma }_3=1.8,k_1=0.3,k_2=0.7,d_1=0.3,d_2=0.7$ (right banner of Figure 6) and $k_1=0.7,k_2=0.3,d_1=0.7,d_2=0.3$ (left banner of Figure 6), the condition of extinction results in Theorem 3$a-\frac{{\sigma }_1^2}{2}=-0.125<0$ is satisfied. Figure 6 shows that the tumor cells $T\left(t\right)$ can die out if the white noise value increases. This serves as evidence that the intensity of environmental noise plays an important role to accelerate tumor extinction.

Remark 4.

As compared to Figure 4, Figure 5 and Figure 6, white noise perturbations on growth and destruction rates promote the elimination of tumor cells $T\left(t\right)$. However, increasing the white noise values alters the intra-specific competition rate, which affects the persistence–extinction properties of tumor cells.

6. Conclusions

As enzymatic reactions of proteins are sensitive to environmental changes in the bio-chemical system, in the present paper, we incorporated the noise effect into a tumor-immune model with Holling type III response functions to discuss the fluctuations in cell dynamics. We investigated the existence of global non-negative solutions to such a model. We derived sufficient conditions for unique stationary distribution/stochastic persistence, extinction, and weak persistence of tumor cells, using Ito’s formula and the Lyapunov function. As a result of the stochastic tumor immune model, the following theoretical conditions are obtained:

(i)

If $d-\frac{{\sigma }_1^2\vee {\sigma }_2^2\vee {\sigma }_3^2}{2}>0,$ then the tumor cells $T\left(t\right)$ are stochastic persistent and permanent in the tissue a.s.

(ii)

If $a>\frac{{\sigma }_1^2}{2}$, then $T\left(t\right)$ shows weak persistence/non-persistence in the mean a.s.

(iii)

If $a<\frac{{\sigma }_1^2}{2}$, then the tumor cells $T\left(t\right)$ are eradicated/eliminated a.s.

The numerical simulations are verified based on these sufficient conditions. Furthermore, theoretical and numerical analyses show that the intensity of the noise effect is a prominent factor in controlling and suppressing tumor cell growth in the presence of immune effectors. Dynamical behavior in deterministic and stochastic tumor-immune models is useful to an oncologist and treatment protocols. Treatment protocols can be adapted to cancer patients with the help of environmental fluctuations for better outcomes during public interventions.

Stochastic epidemic models with Markovian switching and time delays are the focus of our future research.