Stationary distribution and global stability of stochastic predator-prey model with disease in prey population

Abstract

In this paper, a new stochastic four-species predator-prey model with disease in the first prey is proposed and studied. First, we present the stochastic model with some biological assumptions and establish the existence of globally positive solutions. Moreover, a condition for species to be permanent and extinction is provided. The above properties can help to save the dangered population in the ecosystem. Through Lyapunov functions, we discuss the asymptotic stability of a positive equilibrium solution for our model. Furthermore, it is also shown that the system has a stationary distribution and indicating the existence of a stable biotic community. Finally, our results of the proposed model have revealed the effect of random fluctuations on the four species ecosystem when adding the alternative food sources for the predator population. To illustrate our theoretical findings, some numerical simulations are presented.

Keywords:

• Predator-prey system
• stochastically permanent
• extinction
• stochastic stability
• stationary distribution

2020 Mathematics Subject Classifications:

• 60H10
• 60H30
• 92D30
• 93E15

1. Introduction

Mathematical modelling in ecology and epidemiology has become an important tool in analysing the spread and control of infectious diseases. Anderson and May [3] were the first to merge the above two fields and formulate a predator-prey model where some disease infected the prey species. Most infected disease models are created from the classical SIR model of Kermack and McKendrick [18]. Subsequently, many researchers have proposed and studied different predator-prey models in the presence of a disease. Studies in the field of eco-epidemiology have increased enormously in the last two decades, and we only mention a few of them [16, 20, 21, 31, 32, 37]. Qin et al. performed Holling III predator-prey model with an infectious predator. A functional response plays a key role in the predator-prey model and describes the number of prey consumed per predator per unit of time for some given time. Here, we used the Beddington-DeAngelis functional response for the predator and the infected first prey interaction. Mathematically, we may consider both prey-dependent and ratio-dependent models as limiting cases of the general Beddington-DeAngelis type predator-prey model. Much progress [10, 12, 13, 15, 35, 36, 39, 40] have been seen in the study of the predator-prey model with Beddington-DeAngelis functional response. Recently, Arora et al. [1] proposed a predator-prey system highlighting the significance of additional food for predators with a Beddington-DeAngelis functional response.

Another important factor in ecological and eco-epidemiological species interaction is alternative food sources for the predator population. Generally, predators do not depend on a single prey species for food. Predators will change to alternative food when the density of their preferred prey is low. Fryxell and Lundberg [11] has argued that alternative food plays an essential role in promoting the persistence of predator-prey systems. Many authors [1, 9, 45] have considered one predator-two prey systems to show the importance of alternative food for the coexistence of species. Thus, a better understanding of the mechanisms, assumptions and biological hypotheses of predator-prey dynamic models involving alternative food is explicitly needed when the disease is present in the prey population. Here, we consider a model with one predator and two prey populations with the availability of alternative food sources for the predator population when a disease has affected only the first prey population. Srinivasu et al. [38] showed that additional food serves as a tool for the biological conservation of predator-prey systems involving type III functional response.

In reality, a population model can be affected by environmental noise; hence, stochastic population systems play a major role in various branches of applied sciences, including population dynamics and biology. In real life, continuous fluctuations in the environment change the population to a lower or higher extent level. May [30] discovered that all the parameters involved in the population system exhibit random fluctuation as the factors controlling them are not constant. Therefore most authors [5,14,25,41,46] introduce stochastic perturbation into the deterministic models and make the model more realistic for examination; the resulting dynamics solely depend on the intensity of environmental fluctuations. Precisely, enormous efforts have been taken to analyse the dynamical behaviour with environmental fluctuations [8, 21, 28, 34, 39, 40, 42, 44, 45]. Cao et al. [5] transformed the stochastic system with a weak kernel case into an equivalent degenerate system and studied its behaviour. [14] et al. proposed the predator-prey model with Holling II increasing function in the predator. Lu [27] presented Crowley-Martin predator-prey model in a stochastic environment and analysed asymptotical stability in probability by combining the Khasminskii theory of stability with the Lyapunov method. Zhang et al. [43] investigated the influence of anti-predator behaviour due to the fear induced by the predators for the stochastic eco-epidemiological predator-prey model. Recently, Liu et al. [24] developed and analysed a Leslie-Gower stochastic predator-prey model with impulsive toxicant input in polluted environments. Liu [26] carried out several two-species predator-prey models and showed the impact of randomness. A small amount of work has been carried out with the stochastic four-species predator-prey model and disease in the prey.

The present article is organized as follows. Section 2 presents some known definitions, theorems, and lemmas that will be used in the following analysis. In Section 3, we discuss the detailed explanation about the formulation of the stochastic model and analyse the existence and uniqueness of a global positive solution of the system. In Section 4, stochastic permanence and extinction are established under certain parametric restriction. We construct a suitable Lyapunov function to examine the global asymptotic stability for positive equilibrium in Section 5. In Section 6, we establish the conditions for unique stationary distribution. Additionally, we generate some figures to illustrate our main results in Section 7. Finally, we end the paper with a discussion in Section 8.

2. Preliminaries

Consider the d-dimensional stochastic model (SM) of the form (1) $\mathrm{d}X\left(t\right)=f\left(X\right(t),t)\mathrm{d}t+g\left(X\right(t),t)\mathrm{d}B\left(t\right),$(1) with an initial value $X\left(t_0\right)=X_0,t_0\le t\le T<\mathrm{\infty }$, where $f:R^d\times [t_0,T]\to R^d$ and $g:R^d\times [t_0,T]\to R^{d\times m}$ are Borel measurable, $B=\left\{B\right(t){\}}_{t\ge t_0}$ is an $R^m$-valued Wiener process and $X_0$ is an $R^d$-valued random variable.

The differential operator $\mathcal{L}$ associated with the SM (1(1) $\mathrm{d}X\left(t\right)=f\left(X\right(t),t)\mathrm{d}t+g\left(X\right(t),t)\mathrm{d}B\left(t\right),$(1) ) is given by $\mathcal{L}=\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\sum _{i=1}^d{f}_i(x,t)\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}+\frac{1}{2}\sum _{i,j=1}^m(g(x,t)g^T(x,t){)}_{ij}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}_i\mathrm{\partial }{x}_j}.$In addition to the existence and uniqueness assumptions, we assume that f and g satisfy $f(x^{\ast },t)=0$ and $g(x^{\ast },t)=0$ for an equilibrium solution $x^{\ast }$, for $t\ge t_0$. Moreover we assume that $x_0$ is a non-random constant with probability 1.

Definition 2.1

see [22]

The solution $X\left(t\right)$ of (1(1) $\mathrm{d}X\left(t\right)=f\left(X\right(t),t)\mathrm{d}t+g\left(X\right(t),t)\mathrm{d}B\left(t\right),$(1) ) is said to possess stochastic permanent if for any $ϵ\in (0,1),$ there exists a pair of positive constants $\delta =\delta \left(ϵ\right)$ and $\chi =\chi \left(ϵ\right)$ such that for any initial value $X_0\in R_{+}^4,$ the solution $X\left(t\right)$ to (1(1) $\mathrm{d}X\left(t\right)=f\left(X\right(t),t)\mathrm{d}t+g\left(X\right(t),t)\mathrm{d}B\left(t\right),$(1) ) has the properties $\underset{t\to \mathrm{\infty }}{lim inf}\mathbb{P}\left\{|X\right(t)|\ge \delta \}\ge 1-ϵ,\underset{t\to \mathrm{\infty }}{lim inf}P\left\{|X\right(t)|\le \chi \}\ge 1-ϵ.$

Theorem 2.1

see [2]

Assume that f and g satisfy the existence and uniqueness assumptions and that they have continuous coefficients with respect to t.

1. Suppose that there exists a positive definite function $V\in C^{2,1}(U_h\times [t_0,\mathrm{\infty }\left)\right),$ where $U_h=\{x\in R^d:\Vert x-x^{\ast }\Vert <h\},$ for h>0, such that $\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}t\ge t_0,x\in U_h:\mathcal{L}V(x,t)\le 0.$Then, the equilibrium solution $x^{\ast }$ of (1(1) $\mathrm{d}X\left(t\right)=f\left(X\right(t),t)\mathrm{d}t+g\left(X\right(t),t)\mathrm{d}B\left(t\right),$(1) ) is stochastically stable.

2. If, in addition, V is decreasing and there exists a positive definite function $V_1$ such that $V(x,t)\le V_1\left(x\right)\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}t\ge t_0,x\in U_h:\mathcal{L}V(x,t)\le 0,$then the equilibrium solution $x^{\ast }$ is stochastically asymptotically stable.

Let $X\left(t\right)$ be a homogeneous Markov process defined in the l dimensional Euclidean space $E_l$ and described by the following system of stochastic differential equations: $\mathrm{d}X\left(t\right)=b\left(X\right)\mathrm{d}t+\sum _{r=1}^k{f}_r\left(X\right)\mathrm{d}{B}_r\left(t\right).$The diffusion matrix is defined by Hasminskii [19], $A_4\left(x\right)=\left(a_{ij}\right(x\left)\right),a_{ij}\left(x\right)=\sum _{r=1}^{k_1}{f}_r^i\left(x\right)f_r^j\left(x\right).$We assume there exists a bounded domain $U\subset E_l$ with regular boundary Γ, having the following properties:

1. In some neighbourhood of a domain U, the smallest eigenvalue of the diffusion matrix $A_4\left(x\right)$ is bounded away from zero.

2. If $x\in E_l\mathrm{\setminus }U$, the mean time τ at which a path emerging from x reaches the set U is finite, and $\underset{x\in S}{sup}{E}_x\tau <\mathrm{\infty }$ for every compact subset $S\subset E_l$.

Lemma 2.1

If above assumptions hold, then the Markov process $X\left(t\right)$ has a stationary distribution $\phi (\cdot )$. Let $g(\cdot )$ be a function integrable with respect to the measure φ. Then $P_x\left\{\underset{T\to \mathrm{\infty }}{lim}\frac{1}{T}{\int }_0^{T}g\right(X\left(t\right))\mathrm{d}t={\int }_{E_l}g(x\left)\mu \right(\mathrm{d}x\left)\right\}=1,$for all $x\in E_l$.

Lemma 2.2

The following inequality holds $u\le 2(u+1-\log (u\left)\right)-(4-2\log 2),\mathrm{\forall }u>0.$

3. Problem formulation and assumptions

3.1. Mathematical model

In this section, we develop the predator-prey model with two prey and one predator population where the disease infects only the first prey population. Sharma and Samanta [37] pay attention to a similar model and analyse the global stability of the positive equilibrium. We develop the model as follows: (2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S,\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I,\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY,\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ,\end{array}$(2) with initial densities (3) $X_S\left(0\right)=X_{S_0}>0,X_I\left(0\right)=X_{I_0}>0,Y\left(0\right)=Y_0>0,Z\left(0\right)=Z_0>0.$(3) Here, $X_S\left(t\right),X_I\left(t\right),Y\left(t\right)$ and $Z\left(t\right)$ denote the numbers of susceptible first prey, infected first prey, second prey and predator, respectively, and all parameters in the model are positive. We choose the Beverton–Holt functional response to represent the birth rate of the prey population. In the absence of disease and predator population, susceptible first prey population will follow the growth equation $\dot{X_S}=\frac{\mathrm{\Pi }{X}_S}{a+X_S}-b{X}_S,$where Π is the maximum per-capita growth rate of prey and b is per-capita mortality rate; moreover, $\frac{\mathrm{\Pi }{X}_S}{a+X_S}$ can be represented as $\mathrm{\Pi }-\frac{\mathrm{\Pi }a}{a+X_S}$, which approaches Π as $X_S$ goes to infinity.

The parameters used in the system modelling with their meaning are listed in the Table 1.

Table 1. Parameters and meaning – predator-prey model with disease in the first prey.

Parameters	Meaning
α	disease transmission coefficient
β	predation rate of infected first prey
m	half saturation constant for infected first prey
γ	intrinsic growth rate of second prey
δ	predation rate of second prey
c	death rate including natural death rate and disease related death rate of the infected first prey
d, e	natural death rate of second prey and predator

The given model requires the following assumptions:

1. Only the first prey population is infected by the disease and because of the presence of the disease, the first prey population is divided into two sub-classes, namely (i) susceptible first prey population $\left(X_S\right)$ and (ii) infected first prey population $\left(X_I\right)$.

2. The infected first prey population dies before having the capability of reproducing and does not recover from the disease.

3. Disease does not spread to predators from the infected first prey.

4. Predators do not consume the susceptible first prey population. It has more possibilities to consume the infected first prey population as they are easily available to be preyed upon owing to weakening due to the infection, i.e. they lose their power to protect themselves from a predator. The predators also consume the second prey because of their significant availability in the total population. Thus, our assumption is reasonable in such cases. This explains why the infected first prey and second prey increase their availability to predation [21].

5. The interaction between the infected first prey and predator is assumed to be governed by a Beddington-DeAngelis type functional response. This functional response was introduced by Beddington [4] and DeAngelis et al. [7]. The interaction between the second prey and predator is considered a Lotka-Volterra functional response. There is no interaction between the first prey and second prey population.

6. The growth of the second prey population is exponential, so the availability of second prey becomes very large when the predator becomes extinct. Because of this, there is no need to consider the search time for the second prey population.

In reality, environmental noises significantly affect population dynamics, which are essential components in an ecosystem. Taking into account the effect of a randomly fluctuating environment, we included white noise in each equation of the system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S,\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I,\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY,\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ,\end{array}$(2) ). Let us assume that the fluctuations in the environment incorporate themselves mainly as fluctuations in the death rate of the susceptible first prey population, the death rate of the infected first prey population, the growth rate of the second prey population, and the death rate of the predator population, as follows: $\begin{array}{rl}& -b\to -b+{\sigma }_1\dot{B_1}\left(t\right),-c\to -c+{\sigma }_2\dot{B_2}\left(t\right),\gamma \to \gamma +{\sigma }_3\dot{B_3}\left(t\right),\\ & -e\to -e+{\sigma }_4\dot{B_4}\left(t\right),\end{array}$where $B_1\left(t\right),B_2\left(t\right),B_3\left(t\right),B_4\left(t\right)$ represent the independent Brownian motions and ${\sigma }_1,{\sigma }_2,{\sigma }_3,{\sigma }_4$ represent the intensities of the white noises. Therefore the corresponding stochastic system can be described by the Itô equations (4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) In the past several years, no work has been done on the given stochastic predator-prey system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ). Our aim is to analyse the dynamical properties of the stochastic model and show the impact of environmental noise in the population system.

Theorem 3.1

All solutions of the system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S,\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I,\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY,\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ,\end{array}$(2) ) with initial conditions (3(3) $X_S\left(0\right)=X_{S_0}>0,X_I\left(0\right)=X_{I_0}>0,Y\left(0\right)=Y_0>0,Z\left(0\right)=Z_0>0.$(3) ) exist in the interval $[0,\mathrm{\infty })$ and $X_S\left(t\right)>0,X_I\left(t\right)>0,Y\left(t\right)>0,Z\left(t\right)>0,$ for all $t\ge 0$.

Proof.

The proof of this theorem is similar to the proof of Theorem 1 in [37].

3.2. Existence of global positive solution

Theorem 3.2

For $(X_{S_0},X_{I_0},Y_0,Z_0)\in Int\left(R_{+}^4\right)$, the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) possesses a unique positive local solution $\left(X_S\right(t),X_I(t),Y(t),Z(t\left)\right)$ for $t\in [0,{\tau }_e)$ almost surely, where ${\tau }_e$ is the explosion time.

Proof.

Consider the following transformation of variables $p=\log X_S,q=\log X_I,r=\log Y,s=\log Z.$Let us use Itô's formula $\mathcal{L}V=V_t(t,x)+V_x(t,x)f(t,x)+\frac{1}{2}\text{trace}\left(g^T\right(t,x\left)V_{xx}\right(t,x\left)g\right(t,x\left)\right).$Using this, we obtain $\begin{array}{rl}\mathcal{L}p& =(\frac{\mathrm{\Pi }}{a+e^p}-\alpha e^q-b-\frac{1}{2}{\sigma }_1^2)\\ \text{and}\mathrm{d}p& =\mathcal{L}p\mathrm{d}t+g(t,p)\mathrm{d}B\left(t\right).\end{array}$Similarly, from system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ), we obtain (5) $\begin{array}{rl}\mathrm{d}p& =(\frac{\mathrm{\Pi }}{a+e^p}-\alpha e^q-b-\frac{1}{2}{\sigma }_1^2)\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1\left(t\right),\\ \mathrm{d}q& =(\alpha e^p-{\displaystyle \frac{\text{β}{e}^s}{m+\mu e^q+\eta e^s}}-c-\frac{1}{2}{\sigma }_2^2)\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2\left(t\right),\\ \mathrm{d}r& =(\gamma -\delta e^s-d-\frac{1}{2}{\sigma }_3^2)\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3\left(t\right),\\ \mathrm{d}s& =(\delta e^r+\frac{\text{β}{e}^q}{m+\mu e^q+\eta e^s}-e)\mathrm{d}t+{\sigma }_4\mathrm{d}{B}_4\left(t\right),\end{array}$(5) with initial condition $p\left(0\right)=\log X_S\left(0\right),q\left(0\right)=\log X_I\left(0\right),r\left(0\right)=\log Y\left(0\right),S\left(0\right)=\log Z\left(0\right)$. Now, the functions associated with the above system have an initial growth and satisfy a local Lipchitz condition. Therefore there exists a unique local solution $\left(p\right(t),q(t),r(t),s(t\left)\right)$ defined in some interval $[0,{\tau }_e)$. Consequently, $X_S\left(t\right)=e^{p\left(t\right)},X_I\left(t\right)=e^{q\left(t\right)},Y\left(t\right)=e^{r\left(t\right)},Z\left(t\right)=e^{s\left(t\right)}$ is the unique positive local solution of (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) starting from an interior point of the first quadrant.

Theorem 3.3

For any initial condition $(X_{S_0},X_{I_0},Y_0,Z_0)\in Int\left(R_{+}^4\right)$, the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) possesses a unique solution $\left(X_S\right(t),X_I(t),Y(t),Z(t\left)\right)$ for $t\in [0,{\tau }_e)$ and the solution remains in $Int\left(R_{+}^4\right)$ with probability one; therefore, for all $t\ge 0,\left(X_S\right(t),X_I(t),Y(t),Z(t\left)\right)\in Int\left(R_{+}^4\right)$ almost surely.

Proof.

To prove that this solution is global, it is enough to show that $\tau =\mathrm{\infty }$ almost surely. Let $k_0$ be a sufficiently large nonnegative integer such that the closed ball $B\left(k_0\right)\in R_{+}^4$ contains $(X_{S_0},X_{I_0},Y_0,Z_0)$. We choose for any $k\ge k_0$ and define the stop-time by using the ideas in [17] as ${\tau }_k=inf\{t\in [0,{\tau }_e):X_S\notin (\frac{1}{k},k)\mathrm{o}\mathrm{r}{X}_I\notin (\frac{1}{k},k)\mathrm{o}\mathrm{r}Y\notin (\frac{1}{k},k)\mathrm{o}\mathrm{r}Z\notin (\frac{1}{k},k)\},$where $inf\mathrm{\varnothing }=\mathrm{\infty }$ (∅ being the empty set). Clearly ${\tau }_k$ is increasing as $k\to \mathrm{\infty }$.

Let ${\tau }_{\mathrm{\infty }}=\underset{k\to \mathrm{\infty }}{lim}{\tau }_k$; then ${\tau }_{\mathrm{\infty }}\le {\tau }_e$ almost surely. If ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$ almost surely is true, then ${\tau }_e=\mathrm{\infty }$ almost surely. If this statement is not true, that is, ${\tau }_{\mathrm{\infty }}\ne \mathrm{\infty }$, then there exist two constants T>0 and $ϵ\in (0,1)$ such that (6) $\mathbb{P}\{{\tau }_{\mathrm{\infty }}\le T\}>ϵ.$(6) Thus, by denoting ${\mathrm{\Omega }}_k=\{{\tau }_k\le T\},$ there is an integer $k_1\ge k_0$ such that, for all $k\ge k_1$, (7) $\mathbb{P}\{{\tau }_k\le T\}\ge ϵ.$(7) Define a $C^4$ function $V:Int\left(R_{+}^4\right)\to Int\left(R_{+}\right)$ by $\begin{array}{rl}V(X_S,X_I,Y,Z)& =(X_S+1-\log X_S)+(X_I+1-\log X_I)\\ & +(Y+1-\log Y)+(Z+1-\log Z),\end{array}$where the function $V(X_S,X_I,Y,Z)$ is positive definite for all $(X_S,X_I,Y,Z)\in Int\left(R_{+}^4\right)$. Using Itô's formula, we obtain $\begin{array}{rl}\mathcal{L}V& =(1-\frac{1}{X_S})(\frac{\mathrm{\Pi }{X}_S}{a+X_S}-\alpha X_S{X}_I-b{X}_S)+\frac{1}{2}{\sigma }_1^2\\ & +(1-\frac{1}{X_I})(\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I)+\frac{1}{2}{\sigma }_2^2\\ & +(1-\frac{1}{Y})(\gamma Y-\delta YZ-dY)+\frac{1}{2}{\sigma }_3^2\\ & +(1-\frac{1}{Z})(\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ)+\frac{1}{2}{\sigma }_4^2,\\ & =(X_S-1)(\frac{\mathrm{\Pi }}{a+X_S}-\alpha X_I-b)+(X_I-1)(\alpha X_S-{\displaystyle \frac{\text{β}Z}{m+\mu X_I+\eta Z}}-c)\\ & +(Y-1)(\gamma -\delta Z-d)+(Z-1)(\delta Y+{\displaystyle \frac{\text{β}{X}_I}{m+\mu X_I+\eta Z}}-e)\\ & +\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2+\frac{1}{2}{\sigma }_4^2.\end{array}$Taking the differential of $V(S,I,P)$ along the solution trajectories of the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) using above equation, we obtain $\begin{array}{rl}\mathrm{d}V(X_S,X_I,Y,Z)& =\left[\right(X_S-1)(\frac{\mathrm{\Pi }}{a+X_S}-\alpha X_I-b)\\ & +(X_I-1)(\alpha X_S-{\displaystyle \frac{\text{β}Z}{m+\mu X_I+\eta Z}}-c)\\ & +(Y-1)(\gamma -\delta Z-d)+(Z-1)(\delta Y+{\displaystyle \frac{\text{β}{X}_I}{m+\mu X_I+\eta Z}}-e)\\ & +\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2+\frac{1}{2}{\sigma }_4^2]\mathrm{d}t\\ & +{\sigma }_1(X_S-1)\mathrm{d}{B}_1+{\sigma }_2(X_I-1)\mathrm{d}{B}_2\\ & +{\sigma }_3(Y-1)\mathrm{d}{B}_3+{\sigma }_4(Z-1)\mathrm{d}{B}_4,\end{array}$Positivity of $X_S\left(t\right),X_I\left(t\right),Y\left(t\right)$ and $Z\left(t\right)$ implies that $\begin{array}{rl}\mathrm{d}V(X_S,X_I,Y,Z)& \le [\frac{\mathrm{\Pi }{X}_S}{a+X_S}+\alpha X_I+b+\frac{\text{β}Z}{m+\mu X_I+\eta Z}+c+\gamma Y+\delta Z+d+e\\ & +\frac{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2+{\sigma }_4^2}{2}]\mathrm{d}t+{\sigma }_1(X_S-1)\mathrm{d}{B}_1+{\sigma }_2(X_I-1)\mathrm{d}{B}_2\\ & +{\sigma }_3(Y-1)\mathrm{d}{B}_3+{\sigma }_4(Z-1)\mathrm{d}{B}_4,\\ & \le [(\frac{\text{β}}{\eta }+b+c+d+e+\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2+\frac{1}{2}{\sigma }_4^2)+\frac{\mathrm{\Pi }{X}_S}{a}\\ & +\alpha X_I+\gamma Y+\delta Z]\mathrm{d}t+{\sigma }_1(X_S-1)\mathrm{d}{B}_1+{\sigma }_2(X_I-1)\mathrm{d}{B}_2\\ & +{\sigma }_3(Y-1)\mathrm{d}{B}_3+{\sigma }_4(Z-1)\mathrm{d}{B}_4.\end{array}$Define three positive constants $c_1=\frac{\text{β}}{\eta }+b+c+d+e+\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2+\frac{1}{2}{\sigma }_4^2,c_2=2(\frac{\mathrm{\Pi }}{a}+\alpha +\gamma +\delta ).$We utilize Lemma 2.2 (the proof is given in [6]) in the form $x_i\le 2(x_i+1-\log (x_i\left)\right)$, we get that ${\displaystyle \frac{\mathrm{\Pi }}{a}}{X}_S+\alpha X_I+\gamma Y+\delta Z\le c_{2}V(X_S,X_I,Y,Z)$. Therefore $\begin{array}{rl}\mathrm{d}V(X_S,X_I,Y,Z)& \le (c_1+c_{2}V(X_S,X_I,Y,Z\left)\right)\mathrm{d}t+{\sigma }_1(X_S-1)\mathrm{d}{B}_1+{\sigma }_2(X_I-1)\mathrm{d}{B}_2\\ & +{\sigma }_3(Y-1)\mathrm{d}{B}_3+{\sigma }_4(Z-1)\mathrm{d}{B}_4,\\ & \le c_3(1+V(X_S,X_I,Y,Z\left)\right)\mathrm{d}t+{\sigma }_1(X_S-1)\mathrm{d}{B}_1+{\sigma }_2(X_I-1)\mathrm{d}{B}_2\\ & +{\sigma }_3(Y-1)\mathrm{d}{B}_3+{\sigma }_4(Z-1)\mathrm{d}{B}_4,\end{array}$where $c_3=\text{max}(c_1,c_2)$. For $t_1\le T,$ $\begin{array}{rl}{\int }_0^{{\tau }_k\wedge t_1}\mathrm{d}V(X_S,X_I,Y,Z)& \le c_3{\int }_0^{{\tau }_k\wedge t_1}(1+V(X_S,X_I,Y,Z\left)\right)\mathrm{d}t+{\sigma }_1{\int }_0^{{\tau }_k\wedge t_1}(X_S-1)\mathrm{d}{B}_1,\\ & +{\sigma }_2{\int }_0^{{\tau }_k\wedge t_1}(X_I-1)\mathrm{d}{B}_2+{\sigma }_3{\int }_0^{{\tau }_k\wedge t_1}(Y-1)\mathrm{d}{B}_3\\ & +{\sigma }_4{\int }_0^{{\tau }_k\wedge t_1}(Z-1)\mathrm{d}{B}_4,\end{array}$where ${\tau }_k\wedge t_1=min\{{\tau }_k,t_1\}$. Using the properties of Itô's integral and with the above inequality we obtain that $\begin{array}{rl}& V\left(X_S\right({\tau }_k\wedge t_1),X_I({\tau }_k\wedge t_1),Y({\tau }_k\wedge t_1),Z({\tau }_k\wedge t_1\left)\right)\\ & \le V\left(X_S\right(0),X_I(0),Y(0),Z(0\left)\right)+c_3{\int }_0^{{\tau }_k\wedge t_1}(1+V(X_S,X_I,Y,Z\left)\right)\mathrm{d}t.\end{array}$ Taking the expectations of both sides of the above inequality and applying Fubini's theorem, we obtain $\begin{array}{rl}& \mathbb{E}V\left(X_S\right({\tau }_k\wedge t_1),X_I({\tau }_k\wedge t_1),Y({\tau }_k\wedge t_1),Z({\tau }_k\wedge t_1\left)\right)\le V\left(X_S\right(0),X_I(0),Y(0),Z(0\left)\right)\\ & +c_3\mathbb{E}{\int }_0^{{\tau }_k\wedge t_1}(1+V(X_S,X_I,Y,Z\left)\right)\mathrm{d}t,\\ & \le V\left(X_S\right(0),X_I(0),Y(0),Z(0\left)\right)+c_3{t}_1+c_3\mathbb{E}{\int }_0^{{\tau }_k\wedge t_1}V(X_S,X_I,Y,Z)\mathrm{d}t,\\ & \le V\left(X_S\right(0),X_I(0),Y(0),Z(0\left)\right)+c_{3}T\\ & +c_3\mathbb{E}{\int }_0^{t_1}V\left(X_S\right({\tau }_k\wedge t),X_I({\tau }_k\wedge t),Y({\tau }_k\wedge t),Z({\tau }_k\wedge t\left)\right)\mathrm{d}t,\\ & \le V\left(X_S\right(0),X_I(0),Y(0),Z(0\left)\right)+c_{3}T\\ & +c_3{\int }_0^{t_1}\mathbb{E}V\left(X_S\right({\tau }_k\wedge t),X_I({\tau }_k\wedge t),Y({\tau }_k\wedge t),Z({\tau }_k\wedge t\left)\right)\mathrm{d}t.\end{array}$By Grownwall's inequality given in [29] and the above inequality, we obtain (8) $EV\left(X_S\right({\tau }_k\wedge T),X_I({\tau }_k\wedge T),Y({\tau }_k\wedge T),Z({\tau }_k\wedge T\left)\right)\le c_4,$(8) where $c_4=\left(V\right(X_S\left(0\right),X_I\left(0\right),Y\left(0\right),Z\left(0\right))+c_{3}T)e^{c_{3}T}.$From (7(7) $\mathbb{P}\{{\tau }_k\le T\}\ge ϵ.$(7) ), we have $P\left({\mathrm{\Omega }}_k\right)\ge ϵ$. It is noteworthy that there is at least one of $X_S({\tau }_k,\omega ),X_I({\tau }_k,\omega ),Y({\tau }_k,\omega )$ and $Z({\tau }_k,\omega )$ belongs to the set $\{k,\frac{1}{k}\}$, for $\omega \in {\mathrm{\Omega }}_k$ and hence $V\left(X_S\right({\tau }_k),X_I({\tau }_k),Y({\tau }_k),Z({\tau }_k\left)\right)$ is no less than the smallest of $k+1-\log k$ and $\frac{1}{k}+1-\log \left(\frac{1}{k}\right)=\frac{1}{k}+1+\log k$. Consequently, $V\left(X_S\right({\tau }_k),X_I({\tau }_k),Y({\tau }_k),Z({\tau }_k\left)\right)\ge (k+1-\log k)\wedge (\frac{1}{k}+1+\log k).$Hence from (6(6) $\mathbb{P}\{{\tau }_{\mathrm{\infty }}\le T\}>ϵ.$(6) ) and (8(8) $EV\left(X_S\right({\tau }_k\wedge T),X_I({\tau }_k\wedge T),Y({\tau }_k\wedge T),Z({\tau }_k\wedge T\left)\right)\le c_4,$(8) ), it follows that $\begin{array}{rl}{c}_4& \ge \mathbb{E}\left[I_{{\mathrm{\Omega }}_k\left(\omega \right)}V\right(X_S({\tau }_k,\omega ),X_I({\tau }_k,\omega ),Y({\tau }_k,\omega ),Z({\tau }_k,\omega )\left)\right],\\ & \ge ϵ\left[\right(k+1-\log k)\wedge (\frac{1}{k}+1+\log k)],\end{array}$where $I_{{\mathrm{\Omega }}_k\left(\omega \right)}$ is the Indicator function of ${\mathrm{\Omega }}_k$.

Let $k\to \mathrm{\infty }$ lead to a contradiction: $\mathrm{\infty }>c_4=\mathrm{\infty }.$ Therefore ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$ almost surely. Thus $\tau =\mathrm{\infty }$ and $\left(X_S\right(t),X_I(t),Y(t),Z(t\left)\right)\in R_{+}^4$ almost surely. The proof is complete.

Theorem 3.3 tells us that the solution of system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) will remain in $R_{+}^4$. With this property, we continue to discuss how the solution varies in $R_{+}^4$ in more detail.

Theorem 3.4

The solutions of the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) are stochastically ultimately bounded for any initial value $W_0=(X_{S_0},X_{I_0},Y_0,Z_0)\in R_{+}^4$.

Proof.

This proof is similar to the proof of Theorem 3 in [44].

4. Long time behaviour of system

Here, we look at the long-time behaviour of the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ), such as permanence and extinction to determine how the solution varies in $R_{+}^4$ in more detail. First, we impose the following hypotheses, which will be useful later on.

1. $max\{b,c,d,e\}+\frac{1}{2}max\{{\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2,{\sigma }_4^2\}<min\{\frac{\mathrm{\Pi }}{a},\alpha ,\gamma ,\frac{\text{β}}{\mu }\},$

2. $\frac{\mathrm{\Pi }}{a}-b-\frac{{\sigma }_1^2}{2}<0,\alpha -c-\frac{{\sigma }_2^2}{2}<0,\gamma -d-\frac{{\sigma }_3^2}{2}<0,\frac{\text{β}}{\eta }+\delta -e-\frac{{\sigma }_4^2}{2}<0.$

In population dynamics, stochastic permanence is one of the most essential properties for species survival. We discuss this property with the fruitful ideas of [28] as follows:

Theorem 4.1

Suppose that assumption $\left(H1\right)$ holds; then system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) will be stochastically permanent.

Proof.

First, we show that for any initial value $W\left(0\right)=\left(X_S\right(0),X_I(0),Y(0),Z(0\left)\right)\in R_{+}^4$, the solution $W\left(t\right)=\left(X_S\right(t),X_I(t),Y(t),Z(t\left)\right)$ satisfies $\underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}\left(\frac{1}{|W\left(t\right){|}^{\kappa }}\right)\le M,$where κ is an arbitrary positive constant satisfying (9) $max\{b,c,d,e\}+{\displaystyle \frac{(\kappa +1)}{2}}max\{{\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2,{\sigma }_4^2\}<min\{{\displaystyle \frac{\mathrm{\Pi }}{a}},\alpha ,\gamma ,{\displaystyle \frac{\text{β}}{\mu }}\}.$(9) By (9(9) $max\{b,c,d,e\}+{\displaystyle \frac{(\kappa +1)}{2}}max\{{\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2,{\sigma }_4^2\}<min\{{\displaystyle \frac{\mathrm{\Pi }}{a}},\alpha ,\gamma ,{\displaystyle \frac{\text{β}}{\mu }}\}.$(9) ), there exists an arbitrary positive constant $\rho >0$ satisfying (10) $\kappa min\{{\displaystyle \frac{\mathrm{\Pi }}{a}},\alpha ,\gamma ,{\displaystyle \frac{\text{β}}{\mu }}\}-\rho -\kappa max\{b,c,d,e\}-{\displaystyle \frac{\kappa (\kappa +1)}{2}}max\{{\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2,{\sigma }_4^2\}>0.$(10) Define $V(X_S,X_I,Y,Z)=X_S+X_I+Y+Z$ for $(X_S,X_I,Y,Z)\in R_{+}^4$ and $U(X_S,X_I,Y,Z)={\displaystyle \frac{1}{V(X_S,X_I,Y,Z)}}$; by Itô's formula, we have $\begin{array}{rl}\mathrm{d}V(X_S,X_I,Y,Z)& =\{X_S(\frac{\mathrm{\Pi }}{a+X_S}-\alpha X_I-b)+X_I(\alpha X_S-\frac{\text{β}Z}{m+\mu X_I+\eta Z}-c)\\ & +Y(\gamma -\delta Z-d)+Z(\delta Y+\frac{\text{β}{X}_I}{m+\mu X_I+\eta Z}-e)\}\mathrm{d}t\\ & +\left({\sigma }_1{X}_S\mathrm{d}{B}_1\right(t)+{\sigma }_2{X}_I\mathrm{d}{B}_2(t)+{\sigma }_{3}Y\mathrm{d}{B}_3(t)+{\sigma }_4\mathrm{d}{B}_4(t\left)\right),\end{array}$and $\begin{array}{rl}\mathrm{d}U\left(W\right)& =\{-U^2(W)[X_S(\frac{\mathrm{\Pi }}{a+X_S}-\alpha X_I-b)+X_I(\alpha X_S-\frac{\text{β}Z}{m+\mu X_I+\eta Z}-c)\\ & +Y(\gamma -\delta Z-d)+Z(\delta Y+\frac{\text{β}{X}_I}{m+\mu X_I+\eta Z}-e)]\}\mathrm{d}t+U^3\left(W\right)\left[\right({\sigma }_1{X}_S{)}^2\\ & +({\sigma }_2{X}_I{)}^2+({\sigma }_{3}Y{)}^2+\left({\sigma }_{4}Z{)}^2\right]\mathrm{d}t-U^2\left(W\right)\left({\sigma }_1{X}_S\mathrm{d}{B}_1\right(t)+{\sigma }_2{X}_I\mathrm{d}{B}_2(t)\\ & +{\sigma }_{3}Y\mathrm{d}{B}_3\left(t\right)+{\sigma }_{4}Y\mathrm{d}{B}_4\left(t\right)),\\ & =\mathcal{L}U\left(W\right)\mathrm{d}t-U^2\left(W\right)\left({\sigma }_1{X}_S\mathrm{d}{B}_1\right(t)+{\sigma }_2{X}_I\mathrm{d}{B}_2(t)+{\sigma }_{3}Y\mathrm{d}{B}_3(t)+{\sigma }_{4}Z\mathrm{d}{B}_4(t\left)\right).\end{array}$Under Assumption $\left(H1\right)$, we can certainly choose a positive constant κ satisfying (9(9) $max\{b,c,d,e\}+{\displaystyle \frac{(\kappa +1)}{2}}max\{{\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2,{\sigma }_4^2\}<min\{{\displaystyle \frac{\mathrm{\Pi }}{a}},\alpha ,\gamma ,{\displaystyle \frac{\text{β}}{\mu }}\}.$(9) ). Again, applying Itô's formula, we obtain $\begin{array}{rl}\mathcal{L}(1+U(w){)}^{\kappa }& =\kappa (1+U(W\left){)}^{\kappa -1}\mathcal{L}U\right(W)+\frac{1}{2}\kappa (\kappa -1\left)\right(1+U\left(W\right){)}^{\kappa -2}{U}^4\left(W\right)\\ & \times \left[\right({\sigma }_1{X}_S{)}^2+({\sigma }_2{X}_I{)}^2+({\sigma }_{3}Y{)}^2+\left({\sigma }_{4}Z{)}^2\right].\end{array}$Then, we can choose $\rho >0$ sufficiently small such that it satisfies (10(10) $\kappa min\{{\displaystyle \frac{\mathrm{\Pi }}{a}},\alpha ,\gamma ,{\displaystyle \frac{\text{β}}{\mu }}\}-\rho -\kappa max\{b,c,d,e\}-{\displaystyle \frac{\kappa (\kappa +1)}{2}}max\{{\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2,{\sigma }_4^2\}>0.$(10) ); consequently $\begin{array}{rl}\mathcal{L}{e}^{\rho t}(1+U(W){)}^{\kappa }& =\rho e^{\rho t}(1+U(W){)}^{\kappa }+e^{\rho t}\mathcal{L}(1+U\left(W\right){)}^{\kappa },\\ & =e^{\rho t}(1+U(W\left){)}^{\kappa -2}[\rho \right(1+U\left(W\right){)}^2+G],\end{array}$where $\begin{array}{rl}G& =-\kappa U^2\left(W\right)\{X_S(\frac{\mathrm{\Pi }}{a+X_S}-\alpha X_I-b)+X_I(\alpha X_S-\frac{\text{β}Z}{m+\mu X_I+\eta Z}-c)\\ & +Y(\gamma -\delta Z-d)+Z(\delta Y+\frac{\text{β}{X}_I}{m+\mu X_I+\eta Z}-e)\}\\ & -\kappa U^3\left(W\right)\{X_S(\frac{\mathrm{\Pi }}{a+X_S}-\alpha X_I-b)+X_I(\alpha X_S-\frac{\text{β}Z}{m+\mu X_I+\eta Z}-c)\\ & +Y(\gamma -\delta Z-d)+Z(\delta Y+\frac{\text{β}{X}_I}{m+\mu X_I+\eta Z}-e)\}\\ & +\kappa U^3\left(W\right)[({\sigma }_1{X}_S{)}^2+({\sigma }_2{X}_I{)}^2+({\sigma }_{3}Y{)}^2\\ & +({\sigma }_{4}Z{)}^2]+\frac{\kappa (\kappa +1)}{2}{U}^4(W\left)[\right({\sigma }_1{X}_S{)}^2+({\sigma }_2{X}_I{)}^2+({\sigma }_{3}Y{)}^2+({\sigma }_{4}Z{)}^2].\end{array}$Let us now discuss the upper boundedness of the function $(1+U(W\left){)}^{\kappa -2}\right[\rho (1+U(W){)}^2+G]$. It is easy to see that $\kappa U^3\left(W\right)[({\sigma }_1{X}_S{)}^2+({\sigma }_2{X}_I{)}^2+({\sigma }_{3}Y{)}^2+({\sigma }_{4}Z{)}^2]\le \kappa U\left(W\right)max\{{\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2,{\sigma }_4^2\},$and $\begin{array}{rl}& \frac{\kappa (\kappa +1)}{2}{U}^4\left(W\right)[({\sigma }_1{X}_S{)}^2+({\sigma }_2{X}_I{)}^2+({\sigma }_{3}Y{)}^2+({\sigma }_{4}Z{)}^2]\\ & \le \frac{\kappa (\kappa +1)}{2}{U}^2\left(W\right)max\{{\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2,{\sigma }_4^2\}.\end{array}$ Hence $\begin{array}{rl}\mathcal{L}{e}^{\rho t}(1+U(W){)}^{\kappa }& =e^{\rho t}(1+U(W\left){)}^{\kappa -2}[\rho \right(1+U\left(W\right){)}^2+G],\\ & \le e^{\rho t}(1+U(W){)}^{\kappa -2}\{\rho +[2\rho -\kappa min\{\frac{\mathrm{\Pi }}{a},\alpha ,\gamma ,\frac{\text{β}}{\mu }\}\\ & +\kappa max\{b,c,d,e\}\phantom{\frac{\mathrm{\Pi }}{a}}+\kappa max\{{\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2,{\sigma }_4^2\}]U\left(W\right)\\ & +[\rho -\kappa min\{\frac{\mathrm{\Pi }}{a},\alpha ,\gamma ,\frac{\text{β}}{\mu }\}\\ & +\kappa max\{b,c,d,e\}+\frac{\kappa (\kappa +1)}{2}max\{{\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2,{\sigma }_4^2\}]U^2\left(W\right)\}.\end{array}$From (9(9) $max\{b,c,d,e\}+{\displaystyle \frac{(\kappa +1)}{2}}max\{{\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2,{\sigma }_4^2\}<min\{{\displaystyle \frac{\mathrm{\Pi }}{a}},\alpha ,\gamma ,{\displaystyle \frac{\text{β}}{\mu }}\}.$(9) ) and (10(10) $\kappa min\{{\displaystyle \frac{\mathrm{\Pi }}{a}},\alpha ,\gamma ,{\displaystyle \frac{\text{β}}{\mu }}\}-\rho -\kappa max\{b,c,d,e\}-{\displaystyle \frac{\kappa (\kappa +1)}{2}}max\{{\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2,{\sigma }_4^2\}>0.$(10) ), we know that there exists a positive constant Q such that $\mathcal{L}{e}^{\rho t}(1+U(W){)}^{\kappa }\le Q{e}^{\rho t}$ which implies that $\mathbb{E}[e^{\rho t}(1+U(W){)}^{\kappa }]\le (1+U\left(0\right){)}^{\kappa }+\frac{Q(e^{\rho t}-1)}{\rho }.$Therefore $\underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}[U^{\kappa }\left(W\right(t\left)\right)]\le \underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}[(1+U(W\left(t\right)){)}^{\kappa }]\le \frac{Q}{\rho }.$Note that $(X_S+X_I+Y+Z{)}^{\kappa }\le 4^{\kappa }({X_S}^4+{X_I}^4+Y^4+Z^4{)}^{\frac{\kappa }{4}}=4^{\kappa }|W{|}^{\kappa },$where $W=(X_S,X_I,Y,Z)\in R_{+}^4$. Accordingly $\underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}\left[\frac{1}{|W\left(t\right){|}^{\kappa }}\right]\le 4^{\kappa }\underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}[U^{\kappa }\left(W\right(t\left)\right)]\le 4^{\kappa }\frac{Q}{\rho }:=M.$Thus, for any $ϵ>0$, let $\nu =(\frac{ϵ}{M}{)}^{\frac{1}{\kappa }}.$ By Chebyshev's inequality, we obtain $P\{|W\left(t\right)|<\nu \}=P\{|W\left(t\right){|}^{-\kappa }>{\nu }^{-\kappa }\}\le \mathbb{E}[|W\left(t\right){|}^{-\kappa }]/{\nu }^{-\kappa }={\nu }^{\kappa }\mathbb{E}[|W\left(t\right){|}^{-\kappa }],$that is, $\underset{t\to \mathrm{\infty }}{lim inf}P\left\{|W\right(t)|\ge \nu \}\ge 1-ϵ,$Similarly, we obtain that for any $ϵ>0$, there exists $\chi >0$ such that $\underset{t\to \mathrm{\infty }}{lim inf}P\left\{|W\right(t)|\le \chi \}\ge 1-ϵ$. Hence the theorem is proved.

In ecology, the worse occurs when a species disappears completely. Here we prove that if the noise is sufficiently large, the solution will become extinct with probability one.

Theorem 4.2

Assume that $\left(H2\right)$ holds, then for any given initial value $W\left(0\right)=\left(X_S\right(0),X_I(0),Y(0),Z(0\left)\right)\in R_{+}^4$, the solution $W\left(t\right)=\left(X_S\right(t),X_I(t),Y(t),Z(t\left)\right)$ of system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) will be extinct with probability one.

Proof.

Define the Lyapunov function $V_1=\ln X_S$. Then, by Itô's formula, we have $d(\ln X_S)=(\frac{\mathrm{\Pi }}{a+X_S}-\alpha X_I-b-\frac{1}{2}{\sigma }_1^2)\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1\left(t\right).$Integrating it from 0 to t yields $\begin{array}{rl}\ln X_S\left(t\right)& \le \ln X_S\left(0\right)+(\frac{\mathrm{\Pi }}{a}-b-\frac{{\sigma }_1^2}{2})t-\alpha {\int }_0^t{X}_I\left(u\right)\mathrm{d}u\\ & +{\sigma }_1{\int }_0^t\mathrm{d}{B}_1\left(u\right),\end{array}$which implies that $\ln X_S\left(t\right)\le \ln X_S\left(0\right)+(\frac{\mathrm{\Pi }}{a}-b-\frac{{\sigma }_1^2}{2})t+{\sigma }_1{B}_1\left(t\right).$Dividing by t on both sides and letting $t\to \mathrm{\infty }$, we obtain $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln X_S\left(t\right)}{t}\le (\frac{\mathrm{\Pi }}{a}-b-\frac{{\sigma }_1^2}{2})<0\text{almost surely}.$Furthermore, define the Lyapunov function $V_2=\ln X_I$; by Itô's formula, we have $d(\ln X_I)=(\alpha X_S-\frac{\text{β}Z}{m+\mu X_I+\eta Z}-c-\frac{1}{2}{\sigma }_2^2)\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2\left(t\right).$Integrating it from 0 to t, we obtain $\begin{array}{rl}\ln X_I\left(t\right)& =\ln X_I\left(0\right)-(c+\frac{{\sigma }_2^2}{2})t+\alpha {\int }_0^t{X}_S\left(u\right)\mathrm{d}u-\text{β}{\int }_0^t\frac{Z\left(u\right)}{m+\eta X_I\left(u\right)+\mu Z\left(u\right)}\mathrm{d}u\\ & +{\sigma }_2{\int }_0^t\mathrm{d}{B}_2\left(u\right).\end{array}$Consequently $\ln X_I\left(t\right)\le \ln X_I\left(0\right)+(\alpha -c-\frac{{\sigma }_2^2}{2})t+{\sigma }_2{B}_2\left(t\right).$Dividing by t and letting $t\to \mathrm{\infty }$, we obtain $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln X_I\left(t\right)}{t}\le (\alpha -c-\frac{{\sigma }_2^2}{2})<0,\text{almost surely}.$Define the Lyapunov function $V_3=\ln Y$; by Itô's formula, we arrive at $d(\ln Y)=(\gamma -\delta Z-d-\frac{1}{2}{\sigma }_3^2)\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3\left(t\right).$Integrating it, we get $\ln Y\left(t\right)=\ln Y\left(0\right)+(\gamma -d-\frac{{\sigma }_3^2}{2})t-\delta {\int }_0^{t}Z\left(u\right)\mathrm{d}u+{\sigma }_3{\int }_0^t\mathrm{d}{B}_3\left(u\right).$Accordingly, $\ln Y\left(t\right)\le \ln Y\left(0\right)+(\gamma -d-\frac{{\sigma }_3^2}{2})t+{\sigma }_3{B}_3\left(t\right).$Then, dividing it by t and letting $t\to \mathrm{\infty }$, we obtain $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln Y\left(t\right)}{t}\le (\gamma -d-\frac{{\sigma }_3^2}{2})<0,\text{almost surely}.$Similarly, we define the Lyapunov function $V_4=\ln Z$ and by Itô's formula, we get $d(\ln Z)=(\delta Y+\frac{\text{β}{X}_I}{m+\mu X_I+\eta Z}-e-\frac{1}{2}{\sigma }_4^2)\mathrm{d}t+{\sigma }_4\mathrm{d}{B}_4\left(t\right).$Integrating it, we have $\begin{array}{rl}\ln Z\left(t\right)& =\ln Z\left(0\right)-(e+\frac{{\sigma }_4^2}{2})t+\text{β}{\int }_0^t\frac{X_I\left(u\right)}{m+\eta X_I\left(u\right)+\mu Z\left(u\right)}\mathrm{d}u\\ & +\delta {\int }_0^{t}Y\left(u\right)\mathrm{d}u+{\sigma }_4{\int }_0^t\mathrm{d}{B}_4\left(u\right).\end{array}$Then $\ln Z\left(t\right)\le \ln Z\left(0\right)+(\frac{\text{β}}{\eta }+\delta -e-\frac{{\sigma }_4^2}{2})t+{\sigma }_4{B}_4\left(t\right).$Dividing it by t and letting $t\to \mathrm{\infty }$, we have $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln Z\left(t\right)}{t}\le (\frac{\text{β}}{\eta }+\delta -e-\frac{{\sigma }_4^2}{2})<0,\text{almost surely}.$The desired assertion is thus proved.

5. Stochastic asymptotic stability

Using the ideas of Pitchaimani and Rajaji [33], we shall analyse the stochastic stability, which will lead to certain conditions. Here, the discussion is entirely about the mathematical theory on stochastic asymptotic stability of the positive equilibrium solution $(X_S^{\ast },X_I^{\ast },Y^{\ast },Z^{\ast })$ of the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ).

The positive equilibrium point $E^{\ast }=({X_S}^{\ast },{X_I}^{\ast },Y^{\ast },Z^{\ast })$ where $Z^{\ast }=\frac{\gamma -d}{\delta }$ and others are obtained as follows:

Let ${X_I}^{\ast }$ be a positive root of the following quadratic equation $g\left({X_I}^{\ast }\right)={{X_I}^{\ast }}^2+C{X_I}^{\ast }+D=0,$where $\begin{array}{rl}C& =\frac{1}{2\alpha \delta \mu (c+a\alpha )}(\alpha \text{β}\gamma -d\alpha \text{β}+cm\alpha \delta +am{\alpha }^2\delta -cd\alpha \eta -ad{\alpha }^2\eta \\ & +c\alpha \gamma \eta +a{\alpha }^2\gamma \eta +bc\delta \mu +ab\alpha \delta \mu -\mathrm{\Pi }\alpha \delta \mu )\end{array}$and $\begin{array}{rl}D& =\frac{1}{2\alpha \delta \mu (c+a\alpha )}(b\text{β}\gamma -bd\text{β}+bmc\delta +abm\alpha \delta -\mathrm{\Pi }m\alpha \delta \\ & -bcd\eta -abd\alpha \eta +\mathrm{\Pi }d\alpha \eta +bc\gamma \eta +ab\alpha \gamma \eta -\mathrm{\Pi }\alpha \gamma \eta ).\end{array}$The roots of the above quadratic equation are ${X_I}^{\ast }=\frac{-C\pm \sqrt{C^2-4D}}{2}.$When any one of the following cases is satisfied, then the equilibrium points ${X_I}^{\ast }$ can have one or two positive values.

1. C<0 and D<0,

2. C<0, D>0 and $C^2-4D>0$,

3. $C>0,D<0$,

$\text{where}{X_S}^{\ast }=\frac{\mathrm{\Pi }-a(a{X_I}^{\ast }+b)}{a{X_I}^{\ast }+b},Y^{\ast }=\frac{e(m+\mu {X_I}^{\ast }+\eta Z^{\ast })-\text{β}{X_I}^{\ast }}{\delta (m+\mu {X_I}^{\ast }+\eta Z^{\ast })}.$The following relations must hold for the positiveness of ${X_S}^{\ast }$ and $Y^{\ast }$ $\mathrm{\Pi }>a^2{X_I}^{\ast }+ab\text{and}e(m+e\mu {X_I}^{\ast }+\eta Z^{\ast })>\text{β}{X_I}^{\ast }\text{or}e\mu >\text{β}.$We have used the criteria in Theorem 2.1 and the ideas of [23] to prove the stochastic system stability.

Theorem 5.1

The equilibrium solution $E^{\ast }=({X_S}^{\ast },{X_I}^{\ast },Y^{\ast },Z^{\ast })$ of the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) is stochastically asymptotically stable, if (11) $\begin{array}{rl}& {X_I}^{\ast }>max\{\frac{Z^{\ast }(\eta +2\mu )}{\mu },\frac{Z^{\ast }\eta }{2\eta +\mu }\}\text{and}\\ & {\displaystyle \frac{1}{2}}{\displaystyle \frac{{\sigma }_1^2{X_S}^{\ast }}{{X_S}^2}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\sigma }_2^2{X_I}^{\ast }}{{X_I}^2}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\sigma }_3^2{Y}^{\ast }}{Y^2}}+\frac{1}{2}{\displaystyle \frac{{\sigma }_4^2{Z}^{\ast }}{Z^2}}\le \varphi (X_S,X_I,Y,Z),\end{array}$(11) $\begin{array}{rl}\text{where}\varphi (X_S,X_I,Y,Z)& =\frac{\mathrm{\Pi }(X_S-{X_S}^{\ast }{)}^2}{(\alpha +X_S)(\alpha +X_S^{\ast })}+\frac{\text{β}(\mu {X_I}^{\ast }-Z^{\ast }(\eta +2\mu \left)\right)}{2(m+\mu {X_I}^{\ast }+\eta Z^{\ast })}(X_I-{X_I}^{\ast }{)}^2\\ & +\frac{\text{β}\left({X_I}^{\ast }\right(2\eta +\mu )-\eta Z^{\ast })}{2(m+\mu {X_I}^{\ast }+\eta Z^{\ast })}(Z-Z^{\ast }{)}^2.\end{array}$

Proof.

Consider the Lyapunov function $\begin{array}{rl}V(X_S,X_I,Y,Z)& =(X_S-{X_S}^{\ast }-{X_S}^{\ast }\ln \frac{X_S}{{X_S}^{\ast }})+(X_I-{X_I}^{\ast }-{X_I}^{\ast }\ln \frac{X_I}{{X_I}^{\ast }})\\ & +(Y-Y^{\ast }-Y^{\ast }\ln \frac{Y}{Y^{\ast }})+(Z-Z^{\ast }-Z^{\ast }\ln \frac{Z}{Z^{\ast }}).\end{array}$The infinitesimal generator $\mathcal{L}$ acting on the Lyapunov function can be written as follows $\begin{array}{rl}\mathcal{L}V& =(X_S-{X_S}^{\ast })(\frac{\mathrm{\Pi }}{a+X_S}-\alpha X_I-b)+(X_I-{X_I}^{\ast })(\alpha X_S-\frac{\text{β}Z}{m+\mu X_I+\eta Z}-c)\\ & +(Y-Y^{\ast })(\gamma -\delta Z-d)+(Z-Z^{\ast })(\delta Y+\frac{\text{β}{X}_I}{m+\mu X_I+\eta Z}-e)\\ & +\frac{1}{2}\frac{{\sigma }_1^2{X_S}^{\ast }}{S^2}+\frac{1}{2}\frac{{\sigma }_2^2{X_I}^{\ast }}{{X_I}^2}+\frac{1}{2}\frac{{\sigma }_3^2{Y}^{\ast }}{Y^2}+\frac{1}{2}\frac{{\sigma }_4^2{Z}^{\ast }}{Z^2},\\ & =(X_S-{X_S}^{\ast })(\frac{\mathrm{\Pi }}{a+X_S}-\frac{\mathrm{\Pi }}{a+{X_S}^{\ast }}-\alpha (X_I-{X_I}^{\ast }\left)\right)+(X_I-{X_I}^{\ast })\left(\alpha \right(X_S-{X_S}^{\ast })\\ & +\frac{\text{β}{Z}^{\ast }}{m+\mu {X_I}^{\ast }+\eta Z^{\ast }}-\frac{\text{β}Z}{m+\mu X_I+\eta Z})+(Y-Y^{\ast })(-\delta (Z-Z^{\ast }\left)\right)\\ & +(Z-Z^{\ast })\left(\delta \right(Y-Y^{\ast })+\frac{\text{β}{X}_I}{m+\mu X_I+\eta Z}-\frac{\text{β}{X_I}^{\ast }}{m+\mu {X_I}^{\ast }+\eta Z^{\ast }})\\ & +\frac{1}{2}\frac{{\sigma }_1^2{X_S}^{\ast }}{{X_S}^2}+\frac{1}{2}\frac{{\sigma }_2^2{X_I}^{\ast }}{{X_I}^2}+\frac{1}{2}\frac{{\sigma }_3^2{Y}^{\ast }}{Y^2}+\frac{1}{2}\frac{{\sigma }_4^2{Z}^{\ast }}{Z^2},\\ & =-\frac{\mathrm{\Pi }(X_S-{X_S}^{\ast }{)}^2}{(\alpha +X_S)(\alpha +X_S^{\ast })}+\left(\frac{\text{β}\mu Z^{\ast }}{(m+\mu X_I+\eta Z)(m+\mu {X_I}^{\ast }+\eta Z^{\ast })}\right)(X_I-{X_I}^{\ast }{)}^2\\ & -\left(\frac{\text{β}\eta {X_I}^{\ast }}{(m+\mu X_I+\eta Z)(m+\mu {X_I}^{\ast }+\eta Z^{\ast })}\right)(Z-Z^{\ast }{)}^2\\ & -\left(\frac{\text{β}(\mu {X_I}^{\ast }-\eta Z^{\ast })}{(m+\mu X_I+\eta Z)(m+\mu {X_I}^{\ast }+\eta Z^{\ast })}\right)(X_I-{X_I}^{\ast })(Z-Z^{\ast })\\ & +\frac{1}{2}\frac{{\sigma }_1^2{X_S}^{\ast }}{{X_S}^2}+\frac{1}{2}\frac{{\sigma }_2^2{X_I}^{\ast }}{{X_I}^2}+\frac{1}{2}\frac{{\sigma }_3^2{Y}^{\ast }}{Y^2}+\frac{1}{2}\frac{{\sigma }_4^2{Z}^{\ast }}{Z^2},\\ & \le -\frac{\mathrm{\Pi }(X_S-{X_S}^{\ast }{)}^2}{(\alpha +X_S)(\alpha +X_S^{\ast })}+\left(\frac{\text{β}\mu Z^{\ast }}{m+\mu {X_I}^{\ast }+\eta Z^{\ast }}\right)(X_I-{X_I}^{\ast }{)}^2-\left(\frac{\text{β}\eta {X_I}^{\ast }}{m+\mu {X_I}^{\ast }+\eta Z^{\ast }}\right)\\ & \times (Z-Z^{\ast }{)}^2-\left(\frac{\text{β}(\mu {X_I}^{\ast }-\eta Z^{\ast })}{m+\mu {X_I}^{\ast }+\eta Z^{\ast }}\right)\left[\frac{(X_I-{X_I}^{\ast }{)}^2+(Z-Z^{\ast }{)}^2}{2}\right]\\ & +\frac{1}{2}\frac{{\sigma }_1^2{X_S}^{\ast }}{{X_S}^2}+\frac{1}{2}\frac{{\sigma }_2^2{X_I}^{\ast }}{{X_I}^2}+\frac{1}{2}\frac{{\sigma }_3^2{Y}^{\ast }}{Y^2}+\frac{1}{2}\frac{{\sigma }_4^2{Z}^{\ast }}{Z^2},\\ & =-\frac{\mathrm{\Pi }(X_S-{X_S}^{\ast }{)}^2}{(\alpha +X_S)(\alpha +X_S^{\ast })}-\left(\frac{\text{β}(\mu {X_I}^{\ast }-Z^{\ast }(\eta +2\mu )}{2(m+\mu {X_I}^{\ast }+\eta Z^{\ast })}\right)(X_I-{X_I}^{\ast }{)}^2\\ & -\left(\frac{\text{β}\left({X_I}^{\ast }\right(2\eta +\mu )-\eta Z^{\ast })}{2(m+\mu {X_I}^{\ast }+\eta Z^{\ast })}\right)(Z-Z^{\ast }{)}^2+\frac{1}{2}\frac{{\sigma }_1^2{X_S}^{\ast }}{{X_S}^2}+\frac{1}{2}\frac{{\sigma }_2^2{X_I}^{\ast }}{{X_I}^2}\\ & +\frac{1}{2}\frac{{\sigma }_3^2{Y}^{\ast }}{Y^2}+\frac{1}{2}\frac{{\sigma }_4^2{Z}^{\ast }}{Z^2},\\ & =-\varphi (X_S,X_I,Y,Z)+\frac{1}{2}\frac{{\sigma }_1^2{X_S}^{\ast }}{{X_S}^2}+\frac{1}{2}\frac{{\sigma }_2^2{X_I}^{\ast }}{{X_I}^2}+\frac{1}{2}\frac{{\sigma }_3^2{Y}^{\ast }}{Y^2}+\frac{1}{2}\frac{{\sigma }_4^2{Z}^{\ast }}{Z^2}.\end{array}$One can observe that $\varphi (X_S,X_I,Y,Z)\ge 0$, if ${X_I}^{\ast }>max\{\frac{Z^{\ast }(\eta +2\mu )}{\mu },\frac{Z^{\ast }\eta }{2\eta +\mu }\}$. By assumption, $\mathcal{L}V(X_S,X_I,Y,Z)$ becomes negative definite; therefore the equilibrium solution $E^{\ast }$ of the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) is stochastically asymptotically stable.

6. Existence of stationary distribution

Here, we shall give the existence of the stationary distribution of prey and predator populations. For this purpose, we must find the stationary distribution for solutions of the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ), which implies stability in the stochastic sense. We use lemma 2.1, (that is given in [19]) to prove the main theorem on the stationary distribution.

Theorem 6.1

Assume $\begin{array}{rl}& \left(\mathrm{i}\right){X_I}^{\ast }>max\{\frac{Z^{\ast }(\eta +2\mu )}{\mu },\frac{Z^{\ast }\eta }{2\eta +\mu }\}\\ & \left(\mathrm{i}\mathrm{i}\right)\upsilon <min\left\{2m\mathrm{\Pi }A\right({X_S}^{\ast }{)}^2,\text{β}(\mu {X_I}^{\ast }-Z^{\ast }(2\mu +\eta \left)\right)\\ & \times {({X_I}^{\ast }+\frac{M_2\mu }{2\text{β}(\mu {X_I}^{\ast }-Z^{\ast }(2\mu +\eta \left)\right)})}^2,\\ & m\delta A(Y^{\ast }{)}^2,(m\delta A+\text{β}\left({X_I}^{\ast }\right(\mu +2\eta )-\eta Z^{\ast }))\\ & \times {(Z^{\ast }+\frac{M_2}{2(m\delta A+\text{β}({X_I}^{\ast }(\mu +2\eta )-\eta Z^{\ast }\left)\right)})}^2\}\end{array}$where $({X_S}^{\ast },{X_I}^{\ast },Y^{\ast },Z^{\ast })$ is the equilibrium of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S,\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I,\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY,\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ,\end{array}$(2) ), $A=m+\mu {X_I}^{\ast }+\eta Z^{\ast },\upsilon =M_{2}A+M_2\eta Z^{\ast }+\frac{M_2^2{\mu }^2}{4\left(\text{β}\right(\mu {X_I}^{\ast }-Z^{\ast }(2\mu +\eta )\left)\right)}+\frac{M_2^2{\eta }^2}{4(m\delta A+\text{β}({X_I}^{\ast }(\mu +2\eta )-\eta Z^{\ast }\left)\right)}\text{and}{M}_2=A({X_S}^{\ast }{\sigma }_1^2+{X_I}^{\ast }{\sigma }_2^2+2Y^{\ast }{\sigma }_3^2+Z^{\ast }{\sigma }_4^2)$ simultaneously. Then there exists a stationary distribution $\phi (\cdot )$ for the stochastic system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ).

Proof.

The proof of this theorem is based on the following positive definite function $V:E_4\to R_{+}$, where $E_4=\text{Int}\left(R_{+}^4\right)$, $\begin{array}{rl}V(X_S,X_I,Y,Z)& =(m+\mu {X_I}^{\ast }+\eta Z^{\ast })(X_S-{X_S}^{\ast }-{X_S}^{\ast }\ln \frac{X_S}{{X_S}^{\ast }})+(m+\mu {X_I}^{\ast }+\eta Z^{\ast })\\ & \times (X_I-{X_I}^{\ast }-{X_I}^{\ast }\ln \frac{X_I}{{X_I}^{\ast }})\\ & +2(m+\mu {X_I}^{\ast }+\eta Z^{\ast })(Y-Y^{\ast }-Y^{\ast }\ln \frac{Y}{Y^{\ast }})\\ & +(m+\mu {X_I}^{\ast }+\eta Z^{\ast })(Z-Z^{\ast }-Z^{\ast }\ln \frac{Z}{Z^{\ast }}).\end{array}$Applying Itô's formula, we obtain $\begin{array}{rl}\mathrm{d}{V}_1& =A(1-\frac{{X_S}^{\ast }}{X_S})\mathrm{d}{X}_S+\frac{A}{2}\frac{{X_S}^{\ast }}{{X_S}^2}(d{X}_S{)}^2,\\ & =\left[A\right(X_S-{X_S}^{\ast })(\frac{\mathrm{\Pi }}{a+X_S}-\alpha X_I-b)+\frac{A}{2}{X_S}^{\ast }{\sigma }_1^2]\mathrm{d}t+{\sigma }_{1}A(X_S-{X_S}^{\ast })\mathrm{d}{B}_1,\\ & =\left[A\right(X_S-{X_S}^{\ast })(-\mathrm{\Pi }(X_S-{X_S}^{\ast })-\alpha (X_I-{X_I}^{\ast }\left)\right)+\frac{A}{2}{X_S}^{\ast }{\sigma }_1^2]\mathrm{d}t\\ & +{\sigma }_{1}A(X_S-{X_S}^{\ast })\mathrm{d}{B}_1,\\ & =[-\mathrm{\Pi }A(X_S-{X_S}^{\ast }{)}^2-\alpha A(X_I-{X_I}^{\ast })(X_S-{X_S}^{\ast })+\frac{A}{2}{X_S}^{\ast }{\sigma }_1^2]\mathrm{d}t\\ & +{\sigma }_{1}A(X_S-{X_S}^{\ast })\mathrm{d}{B}_1,\end{array}$ $\begin{array}{rl}\mathrm{d}{V}_2& =A(1-\frac{{X_I}^{\ast }}{X_I})\mathrm{d}{X}_I+\frac{A}{2}\frac{{X_I}^{\ast }}{{X_I}^2}(\mathrm{d}{X}_I{)}^2,\\ & =\left[A\right(X_I-{X_I}^{\ast })(\alpha X_S-\frac{\text{β}Z}{m+\mu X_I+\eta Z}-c)+\frac{A}{2}{X_I}^{\ast }{\sigma }_2^2]\mathrm{d}t\\ & +{\sigma }_{2}A(X_I-{X_I}^{\ast })\mathrm{d}{B}_2,\\ & =\left[A\right(X_I-{X_I}^{\ast })\left(\alpha \right(X_S-{X_S}^{\ast })+\frac{\text{β}{Z}^{\ast }}{m+\mu {X_I}^{\ast }+\eta Z^{\ast }}-\frac{\text{β}Z}{m+\mu X_I+\eta Z})\\ & +\frac{A}{2}{X_I}^{\ast }{\sigma }_2^2]\mathrm{d}t+{\sigma }_{2}A(X_I-{X_I}^{\ast })\mathrm{d}{B}_2,\\ & =\left[\alpha A\right(X_S-{X_S}^{\ast }\left)\right(X_I-{X_I}^{\ast })-\frac{\text{β}(m+\mu {X_I}^{\ast })}{m+\mu X_I+\eta Z}(X_I-{X_I}^{\ast }\left)\right(Z-Z^{\ast })\\ & +\frac{\text{β}\mu Z^{\ast }}{m+\mu X_I+\eta Z}(X_I-{X_I}^{\ast }{)}^2+\frac{A}{2}{X_I}^{\ast }{\sigma }_2^2]\mathrm{d}t+{\sigma }_{2}A(X_I-{X_I}^{\ast })\mathrm{d}{B}_2,\\ \mathrm{d}{V}_3& =2A(1-\frac{Y^{\ast }}{Y})\mathrm{d}Y+A\frac{Y^{\ast }}{Y^2}(\mathrm{d}Y{)}^2,\\ & =\left[2A\right(Y-Y^{\ast }\left)\right(\gamma -d-\delta Z)+A{Y}^{\ast }{{\sigma }_3}^2]\mathrm{d}t+2{\sigma }_{3}A(Y-Y^{\ast })\mathrm{d}{B}_3,\\ & =[-2A\delta (Y-Y^{\ast }\left)\right(Z-Z^{\ast })+A{Y}^{\ast }{{\sigma }_3}^2]\mathrm{d}t+2{\sigma }_{3}A(Y-Y^{\ast })\mathrm{d}{B}_3,\end{array}$and $\begin{array}{rl}\mathrm{d}{V}_4& =A(1-\frac{Z^{\ast }}{Z})\mathrm{d}Z+\frac{A}{2}\frac{Z^{\ast }}{Z^2}(\mathrm{d}Z{)}^2,\\ & =\left[A\right(Z-Z^{\ast })(\delta Y+\frac{\text{β}{X}_I}{m+\mu X_I+\eta Z}-e)+\frac{A}{2}{Z}^{\ast }{\sigma }_4^2]\mathrm{d}t+{\sigma }_{4}A(Z-Z^{\ast })\mathrm{d}{B}_4,\\ & =\left[A\right(Z-Z^{\ast })\left(\delta \right(Y-Y^{\ast })+\frac{\text{β}{X}_I}{m+\mu X_I+\eta Z}-\frac{\text{β}{X_I}^{\ast }}{m+\mu {X_I}^{\ast }+\eta Z^{\ast }})+\frac{A}{2}{Z}^{\ast }{\sigma }_4^2]\mathrm{d}t\\ & +{\sigma }_{4}A(Z-Z^{\ast })\mathrm{d}{B}_4,\\ & =\left[A\delta \right(Y-Y^{\ast }\left)\right(Z-Z^{\ast })+\frac{\text{β}(m+\eta Z^{\ast })}{m+\mu X_I+\eta Z}(X_I-{X_I}^{\ast }\left)\right(Z-Z^{\ast })\\ & -\frac{\text{β}\eta {X_I}^{\ast }}{m+\mu X_I+\eta Z}(Z-Z^{\ast }{)}^2+\frac{A}{2}{Z}^{\ast }{\sigma }_4^2]\mathrm{d}t+{\sigma }_{4}A(Z-Z^{\ast })\mathrm{d}{B}_4,\end{array}$Therefore $\begin{array}{rl}\mathrm{d}V& =\mathrm{d}{V}_1+\mathrm{d}{V}_2+\mathrm{d}{V}_3+\mathrm{d}{V}_4,\\ & =\mathcal{L}V\mathrm{d}t+{\sigma }_{1}A(X_S-{X_S}^{\ast })\mathrm{d}{B}_1+{\sigma }_{2}A(X_I-{X_I}^{\ast })\mathrm{d}{B}_2+2{\sigma }_{3}A(Y-Y^{\ast })\mathrm{d}{B}_3\\ & +{\sigma }_{4}A(Z-Z^{\ast })\mathrm{d}{B}_4,\end{array}$where $\begin{array}{rl}\mathcal{L}V& =-\mathrm{\Pi }A(X_S-{X_S}^{\ast }{)}^2-\frac{\text{β}(m+\mu {X_I}^{\ast })}{m+\mu X_I+\eta Z}(X_I-{X_I}^{\ast }\left)\right(Z-Z^{\ast })+\frac{\text{β}\mu Z^{\ast }}{m+\mu X_I+\eta Z}\\ & \times (X_I-{X_I}^{\ast }{)}^2-2A\delta (Y-Y^{\ast }\left)\right(Z-Z^{\ast })+A\delta (Y-Y^{\ast }\left)\right(Z-Z^{\ast })\\ & +\frac{\text{β}(m+\eta Z^{\ast })}{m+\mu X_I+\eta Z}(X_I-{X_I}^{\ast })(Z-Z^{\ast })-\frac{\text{β}\eta {X_I}^{\ast }}{m+\mu X_I+\eta Z}(Z-Z^{\ast }{)}^2+\frac{M_2}{2},\\ & =-\mathrm{\Pi }A(X_S-{X_S}^{\ast }{)}^2+\frac{\text{β}\mu Z^{\ast }}{m+\mu X_I+\eta Z}(X_I-{X_I}^{\ast }{)}^2-A\delta (Y-Y^{\ast })(Z-Z^{\ast })\\ & -\frac{\text{β}\eta {X_I}^{\ast }}{m+\mu X_I+\eta Z}(Z-Z^{\ast }{)}^2-\frac{\text{β}(\mu {X_I}^{\ast }-\eta Z^{\ast })}{m+\mu X_I+\eta Z}(X_I-{X_I}^{\ast }\left)\right(Z-Z^{\ast })+\frac{M_2}{2},\\ & \le -\frac{\mathrm{\Pi }A(m+\mu X_I+\eta Z)(X_S-{X_S}^{\ast }{)}^2}{m+\mu X_I+\eta Z}+\frac{\text{β}\mu Z^{\ast }}{m+\mu X_I+\eta Z}(X_I-{X_I}^{\ast }{)}^2\\ & -A\delta \left[\frac{(Y-Y^{\ast }{)}^2+(Z-Z^{\ast }{)}^2}{2}\right]-\frac{\text{β}\eta {X_I}^{\ast }}{m+\mu X_I+\eta Z}(Z-Z^{\ast }{)}^2\\ & -\frac{\text{β}(\mu {X_I}^{\ast }-\eta Z^{\ast })}{m+\mu X_I+\eta Z}\left[\frac{(X_I-{X_I}^{\ast }{)}^2+(Z-Z^{\ast }{)}^2}{2}\right]+\frac{M_2}{2},\\ & \le \frac{1}{2(m+\mu X_I+\eta Z)}[-2\mathrm{\Pi }Am(X_S-{X_S}^{\ast }{)}^2-\text{β}(\mu {X_I}^{\ast }-Z^{\ast }(2\mu +\eta \left)\right)(X_I-{X_I}^{\ast }{)}^2\\ & -Am\delta (Y-Y^{\ast }{)}^2-(Am\delta +\text{β}\left({X_I}^{\ast }\right(\mu +2\eta )-\eta Z^{\ast })\left)\right(Z-Z^{\ast }{)}^2]+M_2.\end{array}$Note that $\begin{array}{rl}2(m+\mu X_I+\eta Z)\mathcal{L}V& =-2\mathrm{\Pi }Am(X_S-{X_S}^{\ast }{)}^2-\text{β}(\mu {X_I}^{\ast }-Z^{\ast }(2\mu +\eta )\left)\right(X_I-{X_I}^{\ast }{)}^2\\ & -Am\delta (Y-Y^{\ast }{)}^2-(Am\delta +\text{β}\left({X_I}^{\ast }\right(\mu +2\eta )-\eta Z^{\ast })\left)\right(Z-Z^{\ast }{)}^2\\ & +M_2(m+\mu X_I+\eta Z),\\ & \le -2\mathrm{\Pi }Am(X_S-{X_S}^{\ast }{)}^2-Am\delta (Y-Y^{\ast }{)}^2\\ & -\text{β}(\mu {X_I}^{\ast }-Z^{\ast }(2\mu +\eta \left)\right)\\ & \times {(X_I-({X_I}^{\ast }+\frac{M_2\mu }{2\text{β}(\mu {X_I}^{\ast }-Z^{\ast }(2\mu +\eta \left)\right)}))}^2\\ & -(Am\delta +\text{β}({X_I}^{\ast }(\mu +2\eta )-\eta Z^{\ast }\left)\right)\\ & \times {(Z-(Z^{\ast }+\frac{M_2}{2(m\delta A+\text{β}({X_I}^{\ast }(\mu +2\eta )-\eta Z^{\ast }\left)\right)}))}^2+\upsilon .\end{array}$Now if υ satisfies the following condition, $\begin{array}{rl}\upsilon & <min\left\{2m\mathrm{\Pi }A\right({X_S}^{\ast }{)}^2,\text{β}(\mu {X_I}^{\ast }-Z^{\ast }(2\mu +\eta \left)\right){({X_I}^{\ast }+\frac{M_2\mu }{2\text{β}(\mu {X_I}^{\ast }-Z^{\ast }(2\mu +\eta \left)\right)})}^2,\\ & m\delta A(Y^{\ast }{)}^2,(m\delta A+\text{β}\left({X_I}^{\ast }\right(\mu +2\eta )-\eta Z^{\ast }))\\ & \times {(Z^{\ast }+\frac{M_2}{2(m\delta A+\text{β}({X_I}^{\ast }(\mu +2\eta )-\eta Z^{\ast }\left)\right)})}^2\}.\end{array}$Here, the ellipsoid $\begin{array}{rl}& 2\mathrm{\Pi }Am(X_S-{X_S}^{\ast }{)}^2+Am\delta (Y-Y^{\ast }{)}^2+\text{β}(\mu {X_I}^{\ast }-Z^{\ast }(2\mu +\eta \left)\right)\\ & \times {(X_I-({X_I}^{\ast }+\frac{M_2\mu }{2\text{β}(\mu {X_I}^{\ast }-Z^{\ast }(2\mu +\eta \left)\right)}))}^2+(Am\delta +\text{β}({X_I}^{\ast }(\mu +2\eta )-\eta Z^{\ast }\left)\right)\\ & \times {(Z-(Z^{\ast }+\frac{M_2}{2(m\delta A+\text{β}({X_I}^{\ast }(\mu +2\eta )-\eta Z^{\ast }\left)\right)}))}^2=\upsilon ,\end{array}$lies entirely in Int$\left(R_4^{+}\right)$. Let U to be a neighbourhood of the ellipsoid with Ū $\subseteq E_4=\text{Int}{R}_4^{+}$, where Ū is the closure of U. Then, for every $x\in U\mathrm{\setminus }{E}_4$, we have $\mathcal{L}V<0$, which implies condition $\left(P2\right)$ in Lemma 2.1 is satisfied.

In addition, there exist $M=min\{{{\sigma }_1}^2{X_S}^2,{{\sigma }_2}^2{X_I}^2,{{\sigma }_3}^2{Y}^2,{{\sigma }_4}^2{Z}^2,\mathrm{\forall }(X_S,X_I,Y,Z)\in$Ū$\}>0$ such that $\sum _{i,j=1}^4{a}_{ij}{\xi }_i{\xi }_j={{\sigma }_1}^2{X_S}^2{{\xi }_1}^2+{{\sigma }_2}^2{X_I}^2{{\xi }_2}^2+{{\sigma }_3}^2{Y}^2{{\xi }_3}^2+{{\sigma }_4}^2{Z}^2{{\xi }_4}^2\ge M|\xi {|}^2,$for all $(X_S,X_I,Y,Z)\in$Ū, $\xi \in R^4$, which implies that condition $\left(P1\right)$ in above Lemma 2.1 is satisfied. Therefore we can conclude from the above discussion that the stochastic system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) has a stationary distribution $\phi (\cdot )$.

7. Numerical simulation

In this section, let us consider some examples to visualize our theoretical findings by simulating the stochastic system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ). To find the strong solution of system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) with a given initial value, we use the Milstein method. Consider the following discretization equations of model (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) $\begin{array}{rl}{X}_{S_{k+1}}& =X_{S_k}+X_{S_k}(\frac{\mathrm{\Pi }}{a+X_{S_k}}-\alpha X_{I_k}-b)\mathrm{\Delta }t+{\sigma }_1{X}_{S_k}\sqrt{\mathrm{\Delta }t}{\xi }_k+\frac{{\sigma }_1^2}{2}{X}_{S_k}({\xi }_k^2-1)\mathrm{\Delta }t,\\ X_{I_{k+1}}& =X_{I_k}+X_{I_k}(\alpha X_{S_k}-\frac{\text{β}{Z}_k}{m+\mu X_{I_k}+\eta Z_k}-c)\mathrm{\Delta }t\\ & +{\sigma }_2{X}_{I_k}\sqrt{\mathrm{\Delta }t}{\eta }_k+\frac{{\sigma }_2^2}{2}{X}_{I_k}({\eta }_k^2-1)\mathrm{\Delta }t,\\ Y_{k+1}& =Y_k+Y_k(\gamma -\delta Z_k-d)\mathrm{\Delta }t+{\sigma }_3{Y}_k\sqrt{\mathrm{\Delta }t}{\zeta }_k+\frac{{\sigma }_3^2}{2}{Y}_k({\zeta }_k^2-1)\mathrm{\Delta }t,\\ Z_{k+1}& =Z_k+Z_k(\delta Y_k+\frac{\text{β}{X}_{I_k}}{m+\mu X_{I_k}+\eta Z_k}-e)\mathrm{\Delta }t+{\sigma }_4{Z}_k\sqrt{\mathrm{\Delta }t}{\varrho }_k+\frac{{\sigma }_4^2}{2}{Z}_k({\varrho }_k^2-1)\mathrm{\Delta }t,\end{array}$where time increment $\mathrm{\Delta }t>0$ and ${\xi }_k,{\eta }_k,{\zeta }_k,{\varrho }_k,k=1,2,3,\dots ,n$, are independent random Gaussian random variables $N(0,1)$ that can be generated numerically by pseudo random number generators.

We can construct suitable simulations for the stochastic system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) and the corresponding deterministic system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S,\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I,\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY,\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ,\end{array}$(2) ) with the help of MATLAB software and possible parameters. Choose the initial value as $X_S\left(0\right)=0.6,X_I\left(0\right)=0.5,Y\left(0\right)=0.5,Z\left(0\right)=0.8$ for all figures that are going to be discussed in this platform. (12) $\begin{array}{rl}& \mathrm{\Pi }=3,a=0.8,\alpha =2,b=1,\text{β}=2.5,m=1,\mu =1,\eta =2,\\ & c=1.2,\gamma =2,\delta =1.5,d=0.1,e=\mathrm{0.5.}\end{array}$(12) The only difference between Figure 1(a,b) is the intensity of white noise. From the figures, we see that the solution curve of the stochastic model (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) always fluctuates around the solution curves of the deterministic system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S,\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I,\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY,\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ,\end{array}$(2) ). We observe that with the decreasing of ${\sigma }_1,{\sigma }_2,{\sigma }_3,{\sigma }_4$ values, the dynamics of the stochastic system become increasingly similar to the deterministic system.

Figure 1. (a) This graph shows that the solution curve of deterministic (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S,\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I,\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY,\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ,\end{array}$(2) ) and stochastic system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) with large noise ${\sigma }_1={\sigma }_2=0.06,{\sigma }_3={\sigma }_4=0.03$ and (b) small noise ${\sigma }_1={\sigma }_2=0.03,{\sigma }_3={\sigma }_4=0.01$.

For the parametric $\mathrm{\Pi }=3,a=3,\alpha =1,b=1,\text{β}=1.5,m=0.5,\mu =1,\eta =1,c=0.5,\gamma =1,\delta =1.5,d=1.3,e=1,$ and intensity values ${\sigma }_1=0.7,{\sigma }_2=0.6,{\sigma }_3={\sigma }_4=0.5$ condition $\left(H2\right)$ in Theorem 4.2 is verified by the observer; hence, Figure 2(a) shows that all the four species of the stochastic system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ), as well as the corresponding deterministic system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S,\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I,\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY,\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ,\end{array}$(2) ), will die out. Furthermore, in Figure 2(b), the parameters chosen are the same as those used in (12(12) $\begin{array}{rl}& \mathrm{\Pi }=3,a=0.8,\alpha =2,b=1,\text{β}=2.5,m=1,\mu =1,\eta =2,\\ & c=1.2,\gamma =2,\delta =1.5,d=0.1,e=\mathrm{0.5.}\end{array}$(12) ); moreover, the intensities ${\sigma }_1={\sigma }_2=0.03,{\sigma }_3={\sigma }_4=0.01$ satisfy condition $\left(H1\right)$ in Theorem 4.1, and simulations also confirm that all three species of both systems are permanent.

Figure 2. (a) This graph shows that the species of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S,\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I,\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY,\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ,\end{array}$(2) ) and (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) goes to extinction and (b) each species in both system goes to permanent.

Figure 3 illustrates Theorem 5.1. We choose the parameters values as in (12(12) $\begin{array}{rl}& \mathrm{\Pi }=3,a=0.8,\alpha =2,b=1,\text{β}=2.5,m=1,\mu =1,\eta =2,\\ & c=1.2,\gamma =2,\delta =1.5,d=0.1,e=\mathrm{0.5.}\end{array}$(12) ) with different initial values $(0.3,0.2,0.4,0.5),(0.2,0.5,0.4,0.7),(0.5,0.2,0.6,0.8)$ and $(0.2,0.3,0.3,0.4)$. Our model admits a stochastically asymptotically stable solution around the positive equilibrium for the above parameters. From this, we can conclude that the positive equilibrium agrees with Theorem 5.1 under the assumed conditions. The stochastic stability shown for each of the four species in the phase-portrait plane is presented in Figure 4.

Figure 3. This graph shows the stochastic stability of the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) around the positive equilibrium with different initial values and ${\sigma }_1=0.04,{\sigma }_2=0.03,{\sigma }_3=0.02,{\sigma }_4=0.01$.

Figure 4. The phase trajectories clearly shows the stochastic stability of individual species in the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) with different initial values and with the noise ${\sigma }_1=0.04,{\sigma }_2=0.03,{\sigma }_3=0.02,{\sigma }_4=0.01$ converges to the region where the positive equilibrium occur.

The parameter values were chosen (12(12) $\begin{array}{rl}& \mathrm{\Pi }=3,a=0.8,\alpha =2,b=1,\text{β}=2.5,m=1,\mu =1,\eta =2,\\ & c=1.2,\gamma =2,\delta =1.5,d=0.1,e=\mathrm{0.5.}\end{array}$(12) ), and the choice of ${\sigma }_1,{\sigma }_2,{\sigma }_3$ and ${\sigma }_4$ was based on the conditions needed for the existence of a stationary distribution (see Theorem 6.1). We kept all parameter values fixed and repeated the simulation 10, 000 times but never observed any extinction scenario appear at any future time. Figure 5 shows that all species are distributed normally around the mean values $(0.30,0.85,0.07,1.27)$. After that, we increased the strengths of environmental forcing to ${\sigma }_1=0.06,{\sigma }_2=0.06,{\sigma }_3=0.03,{\sigma }_4=0.03$ and, again, we focussed on the species distribution fluctuate around the deterministic steady-state value; however, the amplitude of fluctuation is more comparable to that of the earlier case given in Figure 6. In this case, the susceptible first prey population is distributed within $(0.17,0.42)$; the infected first prey population occupies within $(0.7,1.1)$, the second prey lies within $(0,0.2)$, and the predator population remains within the range $(0.79,1.78)$. In the earlier case, the distribution of each species was within $(0.22,0.37),(0.75,0.99),(0,0.09)$ and $(1,1.5)$. Figures 5 and 6 show that in stochastic stability, all the solution curves will converge to one particular region. In deterministic stability, all the solution curves will converge to a single point. In Figures 5 and 6, the left-hand side shows that the mean solution of the stochastic system coincided with the solution curves of the deterministic system. Figure 7(a,b) present a clear representation of the stationary distribution for all species in a single graph, similar to Figures 5 and 6.

Figure 5. The left panel shows the solution trajectories of both deterministic and stochastic systems from one simulation run; the right panel shows the stationary distribution of each species in the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) separately from 10,000 simulation runs with intensity of noise ${\sigma }_1={\sigma }_2=0.03,{\sigma }_3={\sigma }_4=0.01$.

Figure 6. The left panel represents the solution trajectories of both system from one simulation run; the right panel represents the stationary distribution of each species in the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) separately from 10,000 simulation run with intensity of noise ${\sigma }_1={\sigma }_2=0.06,{\sigma }_3={\sigma }_4=0.03$.

Figure 7. (a) This graph presents the distribution of all species of system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) in one picture with intensity of small noise ${\sigma }_1={\sigma }_2=0.03,{\sigma }_3={\sigma }_4=0.01$ and (b) large noise ${\sigma }_1={\sigma }_2=0.06,{\sigma }_3={\sigma }_4=0.03$.

8. Discussion

From an applications viewpoint, the population model has become more interesting and studied extensively by several researchers. Our work is the first attempt to consider the stochastic four species predator-prey model with disease in the first prey (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ). The main objective of the present paper is to study the model regarding how the intensities of environmental noises affect each population in the system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) where an alternative food source is considered for predator population and disease has affected only the first prey population. Our analytical findings have pointed out that the intensity of noises can play a significant role in the survival of all species. These result can be verified with numerical support.

First, we explain the rigour of our model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S,\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I,\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY,\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ,\end{array}$(2) ) under several assumptions. We introduce the stochastic term to the system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S,\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I,\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY,\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ,\end{array}$(2) ) and show that the stochastic system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) has a local and global positive solution for any positive initial value. The property of permanence and extinction is more desirable since it means the long-time and short time survival in a population dynamics. Under assumption $\left(H1\right)$, our stochastic system (4(4) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{X}_S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\Pi }{X}_S}{a+X_S}}-\alpha X_S{X}_I-b{X}_S+{\sigma }_1{X}_S{B}_1\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{X}_I}{\mathrm{d}t}}=\alpha X_S{X}_I-{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-c{X}_I+{\sigma }_2{X}_I{B}_2\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}t}}=\gamma Y-\delta YZ-dY+{\sigma }_{3}Y{B}_3\left(t\right),\\ {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=\delta YZ+{\displaystyle \frac{\text{β}{X}_{I}Z}{m+\mu X_I+\eta Z}}-eZ+{\sigma }_{4}Z{B}_4\left(t\right).\end{array}$(4) ) attains the permanence. The stochastic system goes to the extinction state if assumption $\left(H2\right)$ holds. In real life, this may happen when a serious disease or severe weather occurs. Furthermore, we discuss the stochastic asymptotic stability of positive equilibrium with the help of the invariance principle and Lyapunov's second method that are given in Theorem 5.1. Thereafter, we have shown the existence of a stationary distribution for all species under certain parameter restriction discussed in Theorem 6.1. This parametric restriction follows the idea that a large intensity of noise can destabilize the system, and in this particular situation we cannot find a stationary distribution.

In a stochastic system, we see that larger intensities of the white noise lead to increasingly larger fluctuations of the solutions. The numerical result shows that strong noise may make the species extinct, which can occur in reality. In permanence, if we reduce the intensity of noise, solutions become stable. From this aspect, one can see the major difference when comparing the deterministic with stochastic models: that is, in the deterministic system the parameters are assumed to be constant, which is not true biologically in real life. Thus, to make our model more fruitful there arises a need to include the environmental variations while modelling the real world scenario mathematically. The use of additional food has been widely recognized by experimental scientists as one of the important tools for biological balance. The present study gives the importance of additional food for predator and the impact of randomness into the system. With the proper choice of parameters satisfying conditions $\left(H1\right)$ and $\left(H2\right)$, remedial measures could be taken to maintain the balance of the ecosystem.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Our paper contains numerical experimental results, and values for these experiments are included in the paper. The data is freely available.

Additional information

Funding

The work of first author is supported by the DST-INSPIRE Fellowship [grant number DST/INSPIRE Fellowship/2017/IF170244], Department of Science and Technology, New Delhi. The second author is thankful to the DST-FIST [grant number SR/FST/MSI-115/2016(Level-I)], Department of Science and Technology, New Delhi for providing financial support. The last author was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) [grant number 2021R1F1A1048937].