Dynamical Behaviors of the Tumor-immune System in a Stochastic Environment

TL;DR

This paper analyzes the stochastic dynamics of the tumor-immune system, focusing on stability, boundary behaviors, and the effects of environmental noise through mathematical and numerical methods.

Contribution

It introduces a stochastic model for tumor-immune interactions and provides rigorous analysis of its stability, boundary behavior, and stationary distributions, supported by simulations.

Findings

Existence and uniqueness of global positive solutions

Conditions for stochastic permanence of the system

Numerical verification of theoretical results

Abstract

This paper investigates dynamic behaviors of the tumor-immune system perturbed by environmental noise. The model describes the response of the cytotoxic T lymphocyte (CTL) to the growth of an immunogenic tumour. The main methods are stochastic Lyapunov analysis, comparison theorem for stochastic differential equations (SDEs) and strong ergodicity theorem. Firstly, we prove the existence and uniqueness of the global positive solution for the tumor-immune system. Then we go a further step to study the boundaries of moments for tumor cells and effector cells and the asymptotic behavior in the boundary equilibrium points. Furthermore, we discuss the existence and uniqueness of stationary distribution and stochastic permanence of the tumor-immune system. Finally, we give several examples and numerical simulations to verify our results.

Equations

$$\left\{\begin{array}[]{lll}\mathrm{d}x(t)&=&\displaystyle\left(\sigma+\frac{\rho{x(t)y(t)}}{\eta+y(t)}-\mu{x(t)y(t)}-\delta{x(t)}\right)\mathrm{d}t,\!\!\!\!\\ \mathrm{d}y(t)&=&\displaystyle\left(\alpha{y(t)}-\beta{y^{2}(t)}-x(t)y(t)\right)\mathrm{d}t\!\!\!\!\end{array}\right.$$

$$-\delta \mathrm{d}t\to -\delta \mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\alpha \mathrm{d}t\to \alpha \mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),$$

$$\left\{\!\!\!\begin{array}[]{lll}\mathrm{d}x(t)&=&\displaystyle\left(\sigma+\frac{\rho{x(t)y(t)}}{\eta+y(t)}-\mu{x(t)y(t)}-\delta{x(t)}\right)dt+\sigma_{1}{x(t)}\mathrm{d}B_{1}(t),\!\!\!\!\\ \mathrm{d}y(t)&=&\displaystyle\left(\alpha{y(t)}-\beta{y^{2}(t)}-x(t)y(t)\right)dt+\sigma_{2}{y(t)}\mathrm{d}B_{2}(t),\!\!\!\!\end{array}\right.$$

$$\begin{split}{\cal{L}}V(x,y,t)&=V_{t}(x,y,t)+V_{x}(x,y,t)\Big{(}\sigma+\frac{\rho{xy}}{\eta+y}-\mu{xy}-\delta x\Big{)}\\ &~{}~{}+V_{y}(x,y,t)\Big{(}\alpha{y}-\beta{y^{2}}-xy\Big{)}+\frac{1}{2}\Big{(}V_{xx}(x,y,t)\sigma^{2}_{1}x^{2}+V_{yy}(x,y,t)\sigma^{2}_{2}y^{2}\Big{)}.\end{split}$$

$${\tau }_m=\inf \left\{t\in [0,{\tau }_e):\min \{x(t),y(t)\}\le \frac{1}{m}\text{mbox}or\max \{x(t),y(t)\}\ge m\right\},$$

$$\mathbb{P}({\tau }_{\infty }\le T)>\epsilon .$$

$$\mathbb{P}({\tau }_m\le T)\ge \frac{\epsilon }{2}.$$

$$V(x,y)=(x+1-\log x)+(y+1-\log y),\forall (x,y)\in R_{+}^2.$$

$$\mathbb{E}[V(x({\tau }_m\wedge T),y({\tau }_m\wedge T))]=V(x_0,y_0)+\mathbb{E}{\int }_0^{{\tau }_m\wedge T}\mathcal{L}V(x(t),y(t))\mathrm{d}t,$$

$$\begin{array}{rl}\mathcal{L}V(x,y)=& (\sigma +\delta -\alpha +\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2)+\frac{\rho xy}{\eta +y}+x+(\mu +\alpha +\beta )y\\ & -\mu xy-\delta x-\frac{\sigma }{x}-\frac{\rho y}{\eta +y}-\beta y^2-xy\\ \le & (\sigma +\delta +\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2)+(\rho +1)x+(\mu +\alpha +\beta )y.\\ \le & v_1+2(\rho +1)(x+1-\log x)+2(\mu +\alpha +\beta )(y+1-\log y)\\ \le & v_1+v_{2}V(x,y),\end{array}$$

$$\mathbb{E}[V(x({\tau }_m\wedge T),y({\tau }_m\wedge T))]\le V(x_0,y_0)+v_{1}T+v_2\mathbb{E}{\int }_0^{{\tau }_m\wedge T}V(x(t),y(t))\mathrm{d}t.$$

$$\mathbb{E}[V(x({\tau }_m\wedge T),y({\tau }_m\wedge T))]\le (V(x_0,y_0)+v_{1}T)e^{v_{2}T}.$$

$$(m+1-\log m)\wedge (\frac{1}{m}+1+\log m)\le V(x({\tau }_m\wedge T),y({\tau }_m\wedge T)).$$

$$\begin{array}{rl}\frac{\epsilon }{2}(m+1-\log m)\wedge (\frac{1}{m}+1+\log m)& \le \mathbb{E}[{\mathbf{I}}_{{\Omega }_m}(\omega )V(x({\tau }_m\wedge T),y({\tau }_m\wedge T))]\\ & \le (V(x_0,y_0)+v_{1}T)e^{v_{2}T},\end{array}$$

$$\underset{m\to \infty }{\lim }(m+1-\log m)\wedge (\frac{1}{m}+1+\log m)=\infty ,$$

$$\infty \le (V(x_0,y_0)+v_{1}T)e^{v_{2}T}<\infty ,$$

$$\left\{\begin{array}[]{lll}\mathrm{d}\psi(t)=\psi(t)\left[\alpha-\beta\psi(t)\right]\mathrm{d}t+\sigma_{2}\psi(t)\mathrm{d}B_{2}(t),\\ \psi(0)=y_{0}>0,\end{array}\right.$$

$$\mathbb{E}{\psi }^k(t)\le {\left[\frac{1}{x_0}{e}^{-(\alpha +\frac{k-1}{2}{\sigma }_2^2)t}+\frac{2\beta }{2\alpha +(k-1){\sigma }_2^2}\left(1-e^{-(\alpha +\frac{k-1}{2}{\sigma }_2^2)t}\right)\right]}^{-k}.$$

$$\limsup\limits_{t\rightarrow\infty}\mathbb{E}\psi^{k}(t)\leq\varrho_{k}:=\Big{(}\frac{2\alpha+(k-1)\sigma^{2}_{2}}{2\beta}\Big{)}^{k},\,\forall k>1.$$

$$\underset{t\to +\infty }{\mathrm{limsup}}\mathbb{E}[y^k(t)]\le {\varrho }_k;$$

$$\underset{t\to +\infty }{\mathrm{limsup}}\mathbb{E}[y^k(t)]\le ({\varrho }_2{)}^{\frac{k}{2}}.$$

$$\underset{t\to +\infty }{\mathrm{limsup}}\mathbb{E}[y^k(t)]\le \underset{t\to +\infty }{\mathrm{limsup}}\mathbb{E}[{\psi }^k(t)]\le {\varrho }_k.$$

$$\underset{t\to +\infty }{\mathrm{limsup}}\mathbb{E}[y^k(t)]\le \underset{t\to +\infty }{\mathrm{limsup}}{\left[\mathbb{E}[y^2(t)]\right]}^{\frac{k}{2}}\le ({\varrho }_2{)}^{\frac{k}{2}},$$

$$\underset{t\to \infty }{\mathrm{limsup}}\mathbb{E}[(1+x(t)+cy(t){)}^{\theta }]\le L(c,\theta ),$$

$$f_1^{\prime }(y)=\frac{-(\mu +c)y^2-2(\mu +c)\eta y+(\rho -\eta (\mu +c))\eta }{(\eta +y{)}^2}.$$

$$y_1=-\eta -\sqrt{\frac{\eta \rho }{\mu +c}}<0,y_2=-\eta +\sqrt{\frac{\eta \rho }{\mu +c}}.$$

$$V_1(x,y)=(1+x+cy{)}^{\theta },\forall (x,y)\in {\mathbb{R}}_{+}^2.$$

$$\displaystyle\begin{split}&\displaystyle\mathcal{L}V_{1}(x,y)\\ =&\displaystyle\theta(1+x+cy)^{\theta-2}\Bigg{\{}\bigg{[}\frac{\rho y}{\eta+y}-(\mu+c)y-\delta+\frac{\theta-1}{2}\sigma^{2}_{1}\bigg{]}x^{2}\\ &~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\Big{[}-c\big{(}\beta+\mu+c\big{)}y^{2}+\big{(}c\alpha+\frac{c\rho y}{\eta+y}-c\delta-\mu-c\big{)}y+\frac{\rho y}{\eta+y}-\delta+\sigma\Big{]}x\\ &~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-c^{2}\beta y^{3}+\Big{(}c^{2}\alpha-c\beta+\frac{\theta-1}{2}c^{2}\sigma^{2}_{2}\Big{)}y^{2}+c(\alpha+\sigma)y+\sigma\Bigg{\}}.\end{split}$$

$$L_1(\theta ):=\delta +\frac{(1-\theta )}{2}{\sigma }_1^2-\frac{\kappa }{\theta }>0.$$

$$M_{v_1}(t):=e^{\kappa t}{V}_1(x(t),y(t))-V_1(x_0,y_0)-{\int }_0^t\mathcal{L}[e^{\kappa t}{V}_1(x(s),y(s))]\mathrm{d}s$$

Full text

Keywords. Tumor-immune system; Stochastic permanence; Comparison theorem; Invariant measure; Ergodicity.

1 Introduction

At present, cancer is considered to be one of the most complicated diseases to be treated clinically and one of the most dreadful killers in the world today. Keeping in mind its devastating nature, a great deal of human and economic resources are devoted to the research on cancer biology and subsequent development of proper therapeutic measures. Surgery, radiation therapy, and chemotherapy are the three traditional therapy procedures that are practised for treatment of cancer. However, all these procedures are characterized by a relatively low efficacy and high toxicity for the patient. Therefore, compared with traditional treatment methods, emerging immunotherapy has great development prospects. Immunotherapy, also known as biological therapy, usually refers to the use of cytokines, a protein hormone that mediate both natural and specific immunity to induce antitumor responses of immune system.

Mathematical models of tumour-immune system and their dynamical behaviors [1, 3], help us to understand better how host immune cells and cancerous cells evolve and interact. In order to get closer to reality more and more tumour-immune models have been studied, for instance, [5, 7, 9, 18, 21, 27, 32, 33, 34] and reference therein. It’s worth noticing that a classical mathematical simplified tumour-immune model

[TABLE]

is proposed to simulate the interaction of the CTL with immunogenic tumor cells and took into account the inactivation of effector cells as well as the penetration of effector cells into tumor cells by Kuznetsov and Taylor [19], where $x$ represents non-dimensional local concentration of effector cells (EC),  $y$ represents the non-dimensional local concentration of  tumor cells (TC). Their model can be applied to describe two different mechanisms of the tumor: tumor dormancy and sneaking through. Yafia [34] studied the stability of the equilibriums and proved the existence of a family of periodic solutions bifurcating from the nontrivial steady state of the Kuznetsov-Taylor model with a delay. More complete bibliography about the evolution of cells and the relevant role of cellular phenomena in directing the body toward recovery or toward illness can be found in [6, 26, 28].

In the tumor tissue, the growth rate and cytotoxic parameters are influenced by many environmental factors, e.g. the degree of vascularization of tissues, the supply of oxygen, the supply of nutrients, the immunological state of the host, chemical agents, temperature, radiations, gene expression, protein synthesis and antigen shedding from the cell surface, etc. Due to the complexity, it is unavoidable that in the course of time the parameters of the system undergo random variations which give them a stochastic character [11, 13, 14, 22]. Inclusion of randomness in mathematical models of biological and biochemical processes is thus necessary for better understanding of mechanisms which govern the biological systems. Considering the impact of the stochastic volatility of environment, we assume that environmental fluctuations mainly affect the culling rate of effector cells  $\delta$ and the intrinsic growth rate of tumor cells  $\alpha$.

[TABLE]

where $B_{1}(t)$ and $B_{2}(t)$ are the $1$-dimensional Brown motion and independent, and $\sigma_{1}$ and $\sigma_{2}$ denote the intensity of white noises. Thus the stochastic tumor-immune model is described by the following SDE

[TABLE]

with an initial value $x(0)=x_{0}>0$, $y(0)=y_{0}>0$. Based on the actual background of the model, we assume that $\sigma_{1},\sigma_{2}$ and all other parameters are non-negative real numbers. Obviously, the model  $(\ref{eq1.5})$ degenerates into  $(\ref{eq1.4})$ if  $\sigma_{1}=0$, $\sigma_{2}=0$.

In the last years, stochastic growth models for cancer cells have been developed, one can see [2, 12, 31] and reference therein. Lyapunov exponent method and Fokker-Planck method are used to investigate the stability of the stochastic models by numerical simulations. Mukhopadhyay and Bhattacharyya [25] analyzed the stochastic stability for a stochastic virus-tumor-immune model. Riccardo, Dumitru and Oana [29] studied the stochastic stability of the stochastic  Kuznetsov-Taylor model by constructing the Lyapunov function. Li and Cheng [20] established the tumor-immune model describing the interaction and competition between the tumor cells and immune system based on the Michaelis-Menten enzyme kinetics, and gave the threshold conditions for extinction, weak persistence and stochastic persistence of tumor cells by the rigorous theoretical proofs, to name a few.

In this paper our main aim is to investigate the stochastic Kuznetsov-Taylor tumor-immune model (1.2), which describes the response of the CTL to the growth of immunogenic tumor cells. Combing the stochastic Lyapunov analysis with the comparison principle for SDEs and making use of the strong ergodicity theorem, we discuss the asymptotic behaviors including the stochastic ultimately boundedness in moment, the limit distribution as well as the ergodicity. Especially, it is pointed out that when tumor cells subject to strong stochastic perturbations, the density of tumor cells is exponentially decreasing while the density of effector cells tends to the stationary distribution. On the other hand, due to weak noises, existence and uniqueness of the stationary distribution with the support set in $\mathbb{R}^{2}_{+}$ is yielded, which implies that tumor cells and immune cells are stochastically permanence. These obtained judgement criteria on extinction and permanence will provide us some inspirations on how to make more effective and precise therapeutic schedule to eliminate tumor cells and improve the treatment of cancer.

The rest of the paper is arranged as follow. Section 2 gives some notation and proves the existence of the unique global positive solution. Section 3 obtains the ultimate moment boundedness of the global positive solution. Section 4 yields the ergodicity of tumor cells and effector cells in the stochastic tumor-immune model which implies the stochastic permanence of cells. Section 5 presents a couple of examples and numerical simulations to illustrate our results. Section 6 concludes this paper.

2 Global positive solution

Throughout this paper, let $(\Omega,{\cal F},\{{\cal F}_{t}\}_{t\geq 0},\mathbb{P})$ be a complete filtered probability space with $\{{\cal F}_{t}\}_{t\geq 0}$ satisfying the usual conditions (that is, it is right continuous and ${\cal F}_{0}$ contains all $\mathbb{P}$-null sets). Let $B(t)=(B_{1}(t),B_{2}(t))^{T}$ be an $2$-dimensional Brownian motion defined on the probability space. Let $|\cdot|$ denote both the Euclidean norm in $\mathbb{R}^{2}$. Also let $\mathbb{R}_{+}=\{x\in\mathbb{R}|x>0\}$ and $\mathbb{R}_{+}^{2}=\{(x,y)\in\mathbb{R}^{2}|x>0,y>0\}$. Also let $C$ denote a generic positive constant whose value may change in different appearances. Moreover, let $C^{2,1}(\mathbb{R}^{2}\times\mathbb{R}_{+};{{\mathbb{R}}}_{+})$ denote the family of all nonnegative functions $V(x,t)$ on $\mathbb{R}^{2}\times\mathbb{R}_{+}$ which are continuously twice differentiable in $x$ and once differentiable in $t$. For each $V\in C^{2,1}(\mathbb{R}^{2}\times\mathbb{R}_{+};\mathbb{R}_{+})$, define an operator ${\cal{L}}$ such that ${\cal{L}}V$: $\mathbb{R}^{2}\times\mathbb{R}_{+}\rightarrow\mathbb{R}$ satisfying

[TABLE]

Since $x(t)$ represents the density of  EC,  $y(t)$ represents the density of  TC, both $x(t)$ and $y(t)$ in (1.2) should be positive. The theorem below gives an affirmative answer.

Theorem 2.1

For any initial value $(x_{0},y_{0})\in{\mathbb{R}^{2}_{+}},$ the equation $(\ref{eq1.5})$ has a unique global positive solution $(x(t),y(t))$ for all $t\geq 0$ with probability one.

Proof. Note that the coefficients of $(\ref{eq1.5})$ are locally Lipschitz continuous, so for any given initial value $(x_{0},y_{0})\in{\mathbb{R}^{2}_{+}}$, there is a unique positive local solution $(x(t),y(t))~{}(t\in[0,\tau_{e}))$, where $\tau_{e}$ is the explosion time . To show this solution is global, we need to prove $\tau_{e}=\infty~{}a.s.$ Choose an $m_{0}\geq 1$ such that $x_{0}\in({1}/{m_{0}},m_{0})$, $y_{0}\in({1}/{m_{0}},m_{0})$. For any positive $m\geq m_{0}$, define the stopping time as follows

[TABLE]

We set $\inf\emptyset=\infty$, clearly,  $\tau_{m}\leq\tau_{e}$, and $\tau_{m}$ is increasing as $m\rightarrow\infty$. Let $\tau_{\infty}=\lim\limits_{m\rightarrow\infty}\tau_{m}$,  $\tau_{\infty}\leq\tau_{e}$ a.s. If we can prove that $\tau_{\infty}=\infty~{}a.s.$ then  $\tau_{e}=\infty~{}a.s.$

Here we use the proof method of contradiction. Suppose that  $\tau_{\infty}=\infty~{}a.s.$ doesn’t hold, then there exist constants $T>0$ and $\varepsilon\in(0,1)$ such that

[TABLE]

This implies that exists an $m_{1}\geq m_{0}$ such that for all $m\geq m_{1}$

[TABLE]

Define

[TABLE]

Using the Itô formula, we have

[TABLE]

where

[TABLE]

with  $v_{1}=\sigma+\delta+\frac{1}{2}\sigma_{1}^{2}+\frac{1}{2}\sigma_{2}^{2}$,  $v_{2}=2(\rho+1+\mu+\alpha+\beta)$. This, together with (2.4), implies

[TABLE]

The Gronwall inequality yields that

[TABLE]

Let $\Omega_{m}=\{\omega:\tau_{m}\leq T\}$, then $\forall\omega\in\Omega_{m}$, at least one of $x(\tau_{m}(\omega)\wedge T)$ and $y(\tau_{m}(\omega)\wedge T)$ is equal to $1/m$ or  ${m}$. Hence, we have

[TABLE]

Due to (2.3) and (2.5), we arrive at

[TABLE]

where $\mathbf{I}_{\Omega_{m}}(\omega)$ is the indicate function of  $\Omega_{m}$. On the other hand, one observes

[TABLE]

Taking $m\rightarrow\infty$ in (2.6), we obtain

[TABLE]

which results in a contradiction. The proof is therefore complete. $\quad\Box$

3 Moment boundedness

Based on the existence result of positive solutions, this section focuses on the moment estimation of the processes $x(t)$ and $y(t)$. In order to discuss the uniform boundary of $\mathbb{E}[y^{k}(t)]$, we introduce an auxiliary process $\psi(t)$ described by

[TABLE]

where $B_{2}(t)$ is the Brownian motion defined in $(\ref{eq1.5})$. By utilizing a comparison theorem, one observes that  $0<y(t)\leq\psi(t)$ for all $t\geq 0$ a.s. The following result is taken from [4], we cite it as a lemma.

Lemma 3.1

${}^{\hbox{\cite[cite]{[\@@bibref{}{H2}{}{}]}}}$* Let $\psi(t)$ be the solution of (3.1). The it holds that for any $k>1,$*

[TABLE]

Therefore, we have

[TABLE]

We now investigate the asymptotic properties of the moments of $y(t)$.

Theorem 3.1

For any $k>1$, we have

[TABLE]

For any $0<k\leq 1$, we have

[TABLE]

Proof. Since $0<y(t)\leq\psi(t)$, If $k>1$, by lemma 3.1, we have

[TABLE]

Additionally, If $0<k\leq 1$, by Hölder inequality, we have

[TABLE]

as required. $\quad\Box$

Next, we continue to consider the asymptotic property of the moments of $x(t)$. By virtue of the interaction between $x(t)$ and $y(t)$ and the positivity of $y(t)$ we provide the following sufficient result for the moment boundedness of $x(t)$.

Theorem 3.2

For any $\theta\in(0,1+2\delta/\sigma_{1}^{2})$ and $c>[\rho/\eta-\mu]^{+}$,

[TABLE]

where $L(c,\theta)$ is a positive constant dependent on $\theta$ and $c$, which is defined by (3.7) below.

Proof. Define the function $f_{1}(y)=\displaystyle\frac{\rho y}{\eta+y}-(\mu+c)y$ for any $y\geq 0$, then

[TABLE]

Solving $f_{1}^{\prime}(y)=0$, we obtain two roots

[TABLE]

Therefore, if $y_{2}\leq 0$, namely $\rho\leq(\mu+c)\eta$, we then have   $f_{1}^{\prime}(y)<0$, $\forall y>0$, and  $f_{1}(y)<f_{1}(0)=0$, $\forall y>0$. For a fixed  $\theta\in(0,1+2\delta/\sigma_{1}^{2})$, define

[TABLE]

By virtue of (2.1) computing $\mathcal{L}V_{1}(x,y)$ leads to

[TABLE]

Since $0<\theta<1+2\delta/\sigma_{1}^{2}$, that is  $\displaystyle\delta+\frac{(1-\theta)}{2}\sigma_{1}^{2}>0$. Now, choose a positive constant $\kappa:=\kappa(\theta)$ sufficiently small such that

[TABLE]

By the Itô formula,

[TABLE]

is a local martingale. And

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

Let $n_{0}>0$ be sufficiently large for $x_{0},~{}y_{0}$ lying within the interval $({1}/{n_{0}},n_{0})$. For any constant $n\geq n_{0}$, define the stopping time

[TABLE]

Note $\xi_{n}$ is monotonically increasing and hence has a (finite or infinite) limit. Denote the limit by $\xi_{\infty}$. For any $n$ sufficiently large, we have $\xi_{n}\geq\tau_{n}$, where $\tau_{n}$ is defined by $(\ref{eq3.1})$. By Theorem 2.1, we have  $\tau_{\infty}=\infty$, then  $\xi_{\infty}=\infty$. The local martingale property implies that $\mathbb{E}[M_{v_{1}}(t\wedge\xi_{n})]=0$. That is, for any $t\geq 0$

[TABLE]

From the definition of $\xi_{n}$, we have  $e^{\kappa(t\wedge\xi_{n})}(1+x(t\wedge\xi_{n})+cy(t\wedge\xi_{n}))^{\theta}$ is monotonically increasing. Let $n\rightarrow\infty$, we obtain

[TABLE]

By the monotone convergence theorem,

[TABLE]

By the dominated convergence theorem,

[TABLE]

Therefore, letting $n\rightarrow\infty$ in (3.5) yields

[TABLE]

This implies

[TABLE]

Then

[TABLE]

Letting $t\rightarrow\infty$, we have

[TABLE]

The proof is complete. $\quad\Box$

The positivity of $y(t)$ implies the follow result directly.

Corollary 3.1

For any $\theta\in(0,1+2\delta/\sigma_{1}^{2})$ and $c>[\rho/\eta-\mu]^{+}$,

[TABLE]

where $L(c,\theta)$ is defined in Theorem 3.2.

Due to the inequality direction in stochastic analysis it is difficult to find the lower bound of the moment of $x(t)$. Alternatively, we try to look for the upper bound of the moment of $1/x(t)$. Thus we get the following result.

Lemma 3.2

If  $\theta\in(0,2)$, then there exists an $L>0$ such that

[TABLE]

Proof. Let

[TABLE]

Choosing a positive constant $\kappa$ and applying the Itô formula lead to

[TABLE]

is a local martingale, where

[TABLE]

Using the Young inequality yields

[TABLE]

where

[TABLE]

Let $n_{0}>0$ be sufficiently large for the initial value  ${x}_{0}$ lying within the interval $({1}/{n_{0}},n_{0})$. For any $n\geq n_{0}$, define the stopping time

[TABLE]

Note $\tilde{\xi}_{n}$ is monotonically increasing and hence has a (finite or infinite) limit. Denote the limit by $\tilde{\xi}_{\infty}$. For any $n$ sufficiently large, we have $\tilde{\xi}_{n}\geq\tau_{n}$, where $\tau_{n}$ is defined by $(\ref{eq3.1})$. By Theorem 2.1, we have $\tau_{\infty}=\infty$, so  $\tilde{\xi}_{\infty}=\infty$. The local martingale property implies that $\mathbb{E}[M_{v_{2}}(t\wedge\tilde{\xi}_{n})]=0$. That is, for any $t\geq 0$

[TABLE]

From the definition of $\tilde{\xi}_{n}$, we have  $e^{\kappa(t\wedge\tilde{\xi})}(1+\frac{1}{x(t\wedge\tilde{\xi}_{n})})^{\theta}$ is monotonically increasing. Letting $n\rightarrow\infty$ yields

[TABLE]

By the monotone convergence theorem one notices that as $n\rightarrow\infty$

[TABLE]

Noting that  $\mathbb{E}[y^{5}(t)]$ and $\mathbb{E}\left[y^{2}(t)\right]$ are bounded uniformly with respect to $t\in(0,\infty)$, by the Fubini theorem and (3.9), we obtain

[TABLE]

Using the dominated convergence theorem implies that as $n\rightarrow\infty$

[TABLE]

Therefore, letting $n\rightarrow\infty$ yields

[TABLE]

This together with Theorem 3.1 implies

[TABLE]

where \displaystyle L_{7}:=L_{6}+\frac{\mu}{5}\sup_{t\geq 0}\mathbb{E}\big{(}y^{5}(t)\big{)}+\frac{\mu}{2}\sup_{t\geq 0}\mathbb{E}\big{(}y^{2}(t)\big{)}<\infty, hence

[TABLE]

We therefore obtain

[TABLE]

The proof is complete. $\quad\Box$

4 Existence and uniqueness of invariant measure

This section is devoted to analyze the invariant measure. Define function

[TABLE]

Similar to the analysis of the function $f_{1}(y)$ in Theorem 3.2, we obtain

• $(\mathrm{i})$

If  $\rho\leq\mu\eta$,  $f_{2}(y)<0$,  $\forall y>0$.

• $(\mathrm{ii})$

If $\rho>\mu\eta$,  $f_{2}(y)\leq(\sqrt{\rho}-\sqrt{\mu\eta})^{2},$  $\forall y>0$.

This implies that for any $y>0$, $f_{2}(y)\leq[(\sqrt{\rho}-\sqrt{\mu\eta})\vee 0]^{2}.$ We now introduce a new auxiliary process $\varphi(t)$ with respect to $x(t)$ described by

[TABLE]

where

[TABLE]

If $\delta-h^{2}>0$, by solving the Fokker-Planck equation (see details in [8]), the process $\varphi(t)$ has a unique stationary distribution $\nu(\cdot)$ which is the inverse Gamma distribution with parameter

[TABLE]

with a notation abuse slightly, we write $\phi\sim IG(a_{1},~{}b_{1})$, with probability density

[TABLE]

For any $p>0$, by the strong law of large numbers we deduce that

[TABLE]

Especially, if $p=1$,  $M_{1}=\displaystyle\frac{\sigma}{\delta-h^{2}}$. Moreover,  the stationary distribution of $\displaystyle\frac{1}{\varphi(t)}$ is the Gamma distribution with parameter $a_{1}$ and $b_{1}$, see details in [10]. Therefore, by the Itô formula and the strong law of large numbers, noting that the mean of Gamma distribution is $a_{1}/b_{1},$ we arrive at

[TABLE]

By virtue of the comparison theorem it follows that  $0<x(t)\leq\varphi(t)$ for all $t\geq 0$ a.s. This implies,

[TABLE]

Furthermore, we derive the following result from (4.3).

Lemma 4.1

If $\delta-h^{2}>0$,

[TABLE]

Moreover,

[TABLE]

Now, we consider the auxiliary process  $\psi(t)$ defined by (3.1). If $\displaystyle 2\alpha<\sigma_{2}^{2}$, it can be easily verified that $\lim\limits_{t\rightarrow\infty}\psi(t)=0~{}~{}a.s$. If $\displaystyle 2\alpha>\sigma_{2}^{2}$, by solving the Fokker-Planck equation (see details in [10]), the process $\psi(t)$ has a unique stationary distribution $\lambda(\cdot)$, and obeys the Gamma distribution with parameter

[TABLE]

with a notation abuse slightly, we write $\psi\sim G(a_{2},~{}b_{2})$, with density

[TABLE]

For any $p>0$, by the strong law of large numbers we derive that

[TABLE]

In particular, if $p=1$, we have $\displaystyle\bar{M}_{1}=\frac{1}{\beta}(\alpha-\frac{\sigma_{2}^{2}}{2})$. Therefore, using the Itô formula implies

[TABLE]

By virtue of the comparison theorem it follows that  $0<y(t)\leq\psi(t)$ for all $t\geq 0$ a.s. One observes that

[TABLE]

Furthermore, we yield the following result from (4.7).

Lemma 4.2

If $\displaystyle 2\alpha>\sigma_{2}^{2}$,

[TABLE]

Moreover,

[TABLE]

To obtain more properties of the solution, we go a further step to consider the equation on the boundary

[TABLE]

where  $t_{0}\geq 0$ will be chosen latter. By solving the Fokker-Planck equation (see details in [8]), the process $z(t)$ has a unique stationary distribution $\mu(\cdot)$, and obeys the inverse Gamma distribution with parameter

[TABLE]

With a notation abuse slightly, we write $z\sim IG(a_{3},~{}b_{3}),$ with probability density

[TABLE]

In the following, we will reveal the long-time behavior of the tumor cells and the effector cells, if the intensity of the noise $\sigma_{2}$ is large sufficiently.

Theorem 4.1

If $\lambda_{1}:=\displaystyle\frac{\sigma_{2}^{2}}{2}-\alpha>0$, then we have

[TABLE]

and the distribution of  $x(t)$ converges weakly to a unique invariant probability measure $\pi_{1}(\cdot)$.

Proof. If $\displaystyle\alpha<\frac{\sigma_{2}^{2}}{2}$, by virtue of the Itô formula, it follows from (3.1) that

[TABLE]

which implies that

[TABLE]

For any $\varepsilon>0$, let $t_{0}>0$ be sufficiently large such that  $\mathbb{P}(\Omega_{\varepsilon})>1-\varepsilon$, where

[TABLE]

and

[TABLE]

Case $(1)$. If $\rho\leq\mu\eta$,  $\displaystyle f_{2}(y)=\frac{\rho y}{\eta+y}-\mu y\leq 0$, by the comparison theorem, we have $\mathbb{P}\{z(t)\geq x(t),~{}~{}\forall t\geq t_{0}\}=1$. By the Itô formula, we deduce that for almost all $\omega\in\Omega_{\varepsilon}$, $\forall t\geq t_{0}$,

[TABLE]

Case $(2)$. If $\rho>\mu\eta$,  $f_{2}(0)>0$. Due to (4.12) and the continuity of $f_{2}(y)$ at $y=0$, one may choose $t_{1}\geq t_{0}$ such that for all $t\geq t_{1}$, $y(t)$ is sufficiently small such that $\displaystyle f_{2}(y(t))=\frac{\rho y(t)}{\eta+y(t)}-\mu y(t)>0$. By the comparison theorem, we have $\mathbb{P}\{x(t)\geq z(t),~{}~{}\forall t\geq t_{1}\}=1$. By the Itô formula we deduce that, for almost all $\omega\in\Omega_{\varepsilon}$, $\forall t\geq t_{1}$,

[TABLE]

Therefore

[TABLE]

Let $\pi_{1}^{*}(\cdot)$ be the invariant measure of $\ln z(t)$. In order to show that the distribution of $x(t)$ converges weakly to a probability measure $\pi_{1}(\cdot)$, we only need to prove that the distribution of $\ln x(t)$ converges weakly to $\pi_{1}^{*}(\cdot)$. Let $\mathcal{\mathbb{missing}}{P}(\mathbb{R})$ represents the family of all probability measures on $\mathbb{R}$. For any $\mathbb{P}_{1},\mathbb{P}_{2}\in{\mathbb{P}(\mathbb{R})}$, define the distance as in [23]

[TABLE]

where

[TABLE]

By the Portmanteau theorem, we need to prove that for any $f\in\mathbb{L},$

[TABLE]

Since  the diffusion is nondegenerate, it is well known that as $t\rightarrow\infty$ the distribution of $\ln z(t)$ converges weakly to the unique stationary distribution $\pi_{1}^{*}(\cdot)$, namely,

[TABLE]

We now compute

[TABLE]

This, together with  $(\ref{eq3.25})$ and $(\ref{eq3.26})$, yields

[TABLE]

The proof is therefore complete. $\quad\Box$

In order to investigate the probability law for the small noises we prepare two lemmas.

Lemma 4.3

If $\displaystyle\frac{\sigma_{2}^{2}}{2}<\alpha$, then the property

[TABLE]

holds, where \lambda_{2}:=\displaystyle\frac{1}{\sigma}\Big{[}\frac{\mu}{\beta}\Big{(}\alpha-\frac{\sigma_{2}^{2}}{2}\Big{)}+\delta+\frac{\sigma_{1}^{2}}{2}\Big{]}.

Proof. For any $(x_{0},y_{0})\in\mathbb{R}_{+}^{2}$, using the fact $0<y(t)\leq\psi(t)$ a.s. and the Itô formula, we have

[TABLE]

Letting $t\rightarrow\infty$, by the strong law of large numbers, $(\ref{eq3.38})$ and $(\ref{eq3.21})$ we deduce that

[TABLE]

The proof is complete. $\quad\Box$

Lemma 4.4

If $\delta>h^{2}$ and $\displaystyle\alpha-\frac{\sigma_{2}^{2}}{2}-\frac{\sigma}{\delta-h^{2}}>0$, then the inequality

[TABLE]

holds, where $\lambda_{3}:=\displaystyle\frac{1}{\beta}\left(\alpha-\frac{\sigma_{2}^{2}}{2}-\frac{\sigma}{\delta-h^{2}}\right)$.

Proof. For any $(x_{0},y_{0})\in\mathbb{R}_{+}^{2}$, since $0<x(t)\leq\varphi(t)$ a.s., an application of the Itô formula yields

[TABLE]

Taking  $t\rightarrow\infty$, by $(\ref{eq3.19})$ and $(\ref{eq3.39})$ we have

[TABLE]

The proof is complete. $\quad\Box$

Now, let us prove the existence of the invariant measure of the equation  $(\ref{eq1.5})$.

Theorem 4.2

If $\delta>h^{2}$ and  $\displaystyle\alpha-\frac{\sigma_{2}^{2}}{2}-\frac{\sigma}{\delta-h^{2}}>0$, then the process $(x(t),y(t))$ has an invariant probability measure on $\mathbb{R}_{+}^{2}$.

Proof. Let  $\hbar$ and $H$ be two positive constants such that $\hbar<\min\{\lambda_{1},\lambda_{2}\},~{}H>\max\{\lambda_{1},\lambda_{2}\}$, where $\lambda_{1}$ and $\lambda_{2}$ defined in Theorem 4.1 and Lemma 4.3, respectively. By the Hölder inequality, we have

[TABLE]

therefore

[TABLE]

Moreover, we have

[TABLE]

Hence, combing (4.18) with (4.19) yields

[TABLE]

It follows from Lemma 4.3 that

[TABLE]

Similarly, by Lemma 4.1 and 4.2, one observes that

[TABLE]

and

[TABLE]

Let $A=\{(x,y):\hbar\leq x\leq H,~{}\hbar\leq y\leq H\}$. Choose more precise $\hbar$ and $H$ such that

[TABLE]

From $(\ref{eq3.30})$-$(\ref{eq3.33})$, we derive

[TABLE]

Taking expectation on both sides, we have

[TABLE]

Using the Fatou lemma and the Fubini theorem yields

[TABLE]

where $\mathbb{P}(t,(x,y),\cdot)$ is the transition probability of $(x(t),y(t))$. Obviously, the Markov process $(x(t),y(t))$ on the state space $\mathbb{R}_{+}^{2}=\{x>0,y>0\}$ has the Feller property. Thus, $(\ref{eq3.34})$ and [24] imply that $(x(t),y(t))$ has an invariant probability measure $\pi_{2}(\cdot).$ $\quad\Box$

We now present the uniqueness of the invariant measure of  $(x(t),y(t))$.

Theorem 4.3

Under the conditions of Theorem 4.2 the solution $(x(t),y(t))$ of (1.2) has a unique invariant measure.

Proof.   For convenience, let $2\zeta:=\displaystyle(\delta-h^{2})\left(\alpha-\frac{\sigma_{2}^{2}}{2}\right)-\sigma.$ Obviously, the given conditions imply $\zeta>0.$ Furthermore, choose a constant $c>0$ small sufficiently such that

[TABLE]

Define $U:\mathbb{R}_{+}^{2}\rightarrow\mathbb{R}_{+}$ by

[TABLE]

Computing ${\cal L}U(x,y)$ yields

[TABLE]

Noting that the definition of $h$ in (4.2) and $\delta>h^{2},~{}\displaystyle\frac{\mu y}{x}\leq\frac{\sigma}{2x^{2}}+\frac{\mu^{2}}{2\sigma}$, one observes that

[TABLE]

This yields that there exist positive constants $h,H,\tilde{h},\tilde{H}$ such that

[TABLE]

By [17, pp.106-122], $(x(t),y(t))$ is positive recurrent with respect to $D$. Then the desired assertion follows. $\quad\Box$

Moreover, by [16] and [17], we have the following ergodicity result.

Theorem 4.4

Under the conditions of Theorem 4.2, the model (1.2) has a unique invariant probability measure $\pi_{2}$ with support $\mathbb{R}_{+}^{2}$. Moreover,

• $(\mathrm{i})$

For any $\pi_{2}$-integrable $f(x,y):\mathbb{R}_{+}^{2}\rightarrow\mathbb{R}$, we have

[TABLE]

• $(\mathrm{ii})$

Let $\|\cdot\|_{\mbox{var}}$ denote the total variation norm, for $(x,y)\in\mathbb{R}_{+}^{2}$, we have

[TABLE]

• $(\mathrm{iii})$

For any $\varepsilon>0$, there is  $\zeta\in(0,1)$ such that

[TABLE]

For a biological system the property $(\mathrm{iii})$ of Theorem 4.2 is also called stochastic strong permanence.

5 Examples and numerical simulations

In this section, we mainly illustrate the effects of noise intensity on effector cells and tumor cells. We select the data in [19] and [30] , see Table  4.1 below.

**Table 4.1: The Significance and value of the parameters

parameter Real value/unit Biological significance

     $a$ 0.18 /day the intrinsic growth rate of the TC

     $b$ $2.0\times 10^{-9}$ /day Reciprocal of environmental capacity of TC

     $s$ $1.3\times 10^{4}$ /piece$\times$day the normal rate of inflow into the tumor site for EC

     $d$ 0.0412 /day the coefficient of destruction and migration of EC

     $g$ $2.019\times 10^{7}$ piece the positive constant in response functional

     $q$ 0.1245 /day $q=fK,~{}\displaystyle K=\frac{k_{1}}{k_{2}+k_{3}+k_{-1}}$

     $r_{1}$ $2.422\times 10^{-10}$ /day$\times$piece $r_{1}=Kk_{3}$

     $r_{2}$ $1.101\times 10^{-7}$ day$\times$piece $r_{2}=Kk_{2}$

**

where $f$ is the positive constant of response function, $k_{1}$ and $k_{-1}$ describe the rates of binding of EC to TC and detachment of EC from TC without damaging cells, $k_{2}$ is the rate at which EC-TC interactions irreversibly program TC for lysis, and $k_{3}$ is the rate at which EC-TC interactions inactivate EC. The non-dimensional treatment of the equation is done by selecting the order of magnitude scales $E_{0}$ and $T_{0}$ for the $E$ and $T$ cell populations, respectively, where  $E_{0}=T_{0}=10^{6}$ cells [19]. Using the non-dimensionalization method in [19] yields coefficients in the model (1.2) as follows

[TABLE]

In addition, applying the Milstein method in Higham [15], we obtain the discrete equation as follows:

[TABLE]

where $\xi_{k}$, $\eta_{k}$ ($k=1,2,\ldots$) are two independent Gaussian random variables, and both obey the normal distribution with mean  0 and variance  1.

Example 5.1

Choose the noise intensities $\sigma_{1}=0.2$, $\sigma_{2}=2$ in the stochastic tumor-immune model $(\ref{eq1.5})$. Then we have

[TABLE]

Theorem 4.1 tell us that the density of tumor cells  $y(t)$ is exponentially decreasing, see the right side of Figure 1. Meanwhile, Theorem 4.1 also shows that the distribution of effector cells $x(t)$ weakly converges to the unique invariant probability measure $\pi_{1}(\cdot)$, the inverse gamma distribution with $a_{3}=19.715$ and $b_{3}=5.905$. To further illustrate the result of Theorem 4.1, we use the $\mathrm{K}$-$\mathrm{S}$ test with a significance level of $0.05$ to check if the stationary distribution of $x^{-1}(t)$ is the gamma distribution. At this level of significance, by $\mathrm{Matlab}$ we do confirm that the stationary distribution of $x^{-1}(t)$ is the Gamma distribution. And because $x^{-1}(t)\sim G(19.715,~{}5.905)$ is equivalent to $x(t)\sim IG(19.715,~{}5.905)$, we know that the stationary distribution of $x(t)$ is the inverse gamma distribution.

Furthermore, to more intuitively illustrate the result of Theorem 4.1, we plot the empirical density function of $x(t)$ and the density function of the inverse gamma distribution $IG(19.715,~{}5.905)$ in Figure 3. One observes obviously from the Figure 3 that as $t\rightarrow\infty$, the distribution of $x(t)$ weakly converges to $\pi_{1}(\cdot)$. Thus, this example illustrates the significance of the result of Theorem 4.1. On the other hand, we compare the simulations of the stochastic model with the deterministic for the same parameter values. Figure 2 depicts that the path of $y(t)$ in the deterministic model tends to a positive equilibrium with frequency vibration, namely, the tumor cells are not extinct. However, in Figure 2, one observes that when the noise intensity is large such that $2\alpha<\sigma_{2}^{2}$, the tumor cells are extinct. It is revealed that stochastic factors cannot be ignored, and their existence plays a key role in the permanence and extinction of the tumor cells.**

Example 5.2

In the stochastic tumor-immune model $(\ref{eq1.5})$, let $\sigma_{2}=0.25$, which implies that the stochastic environment has a weak effect on the intrinsic growth rate of tumor cells. At the same time, the binding rate of $\mathrm{EC}$ to $\mathrm{TC}$ will be decreased when the immune response of the effector cells to the tumor cells is weak or the tumor cells are less irritating to the effector cells. Therefore, in this example we reduce the binding rate $k_{1}$ of $\mathrm{EC}$ and $\mathrm{TC}$ in the literature [19], let $\rho=0.613$. Compute

[TABLE]

These imply that the conditions of Theorems 4.2 and 4.4 hold. By virtue of Theorem 4.2 and 4.4 the solution $(x(t),y(t))$ of the stochastic tumor-immune model $(\ref{eq1.5})$ has a unique invariant probability measure $\pi^{*}$, and the system is stochastically permanent. Figure  4 depicts the trajectories of the effector cells $x(t)$ and the tumor cells $y(t)$ in $(\ref{eq1.5})$. Figure 5 is the phase diagram with respect to the model $(\ref{eq1.5})$. Figure 6 depicts the empirical density of the invariant measure $\pi^{*}$ of the stochastic model $(\ref{eq1.5})$. Therefore, this example verifies the theoretical results of Theorems  4.2 and  4.4 well. **

6 Concluding remarks

This paper mainly studies the dynamical behaviors of the tumor-immune model proposed by Kuznetsov and Taylor [19] perturbed by the environment noise. Firstly, we prove the existence and uniqueness of the global positive solution for the tumor-immune system by the method of stochastic Lyapunov analysis. Next, by constructing appropriate comparison equations, we obtain the asymptotic moment boundedness of the effector cells and the tumor cells. Regarded the boundary equation $(\ref{eq3.24})$ as a bridge, it is pointed out that when tumor cells are subject to the strong noise, the density of the tumor cells decays to zero at an exponential rate while the density of effector cells tends to a stationary distribution. Furthermore, when the noise intensity of tumor cells is small relatively, by analyzing the upper and lower limits of the density of tumor cells and effector cells at time-average, we prove the existence and uniqueness of the stationary distribution of the stochastic tumor-immune model $(\ref{eq1.5})$. Moreover, the ergodicity and the stochastic permanence is obtained . Finally, all of our main results are illustrated and verified by numerical simulations. Overall, the fact is revealed that the intensity of stochastic noise for tumor cells plays a key role in the permanence and extinction of tumor cells. As for strong noise, compared with the deterministic, the dynamical behaviors of the stochastic model are different, and even richer.