# Dynamics of a Stochastic Three-Species Food Web Model with Omnivory and Ratio-Dependent Functional Response

Authors: Guirong Liu; Rong Liu
Journal: Complexity
 (2019)

---

## Abstract

This paper is concerned with a stochastic three-species food web model with omnivory which is defined as feeding on more than one trophic level. The model involves a prey, an intermediate predator, and an omnivorous top predator. First, by the stochastic comparison theorem, we show that there is a unique global positive solution to the model. Next, we investigate the asymptotic pathwise behavior of the model. Then, we conclude that the model is persistent in mean and extinct and discuss the stochastic persistence of the model. Further, by constructing a suitable Lyapunov function, we establish sufficient conditions for the existence of an ergodic stationary distribution to the model. Then, we present the application of the main results in some special models. Finally, we introduce some numerical simulations to support the main results obtained. The results in this paper generalize and improve the previous related results.


---

## Body

## 1. Introduction

The dynamic relationship between predators and their preys has long been and will continue to be one of the dominant themes in both ecology and mathematical ecology due to its universal existence and importance [1]. During the past one hundred years, there have been many investigations on predator-prey models. To the best of our knowledge, in the predator-prey interaction, the functional response plays an important role in the population dynamics, and most of the predator-prey models with the functional responses only depend on the prey. However, laboratory experiments show that the ratio-dependent response function is more reasonable in characterizing the relationship between predators and their preys [2]. Arditi and Ginzburg [3] first proposed a ratio-dependent functional response of form αx/x+βy. Kuang and Beretta [4] investigated the following ratio-dependent type predator-prey model:(1)dx1t=x1tr−a1x1t−α12x2tx1t+β12x2tdt,dx2t=x2t−d2+e12α12x1tx1t+β12x2tdt,where x1t and x2t represent population sizes of prey and predator at time t, respectively. r, a1, and d2 stand for the prey intrinsic growth rate, the intraspecific competition rate of the prey, and the predator death rate, respectively; α12, β12, and e12 represent the encounter rate, half capturing saturation constant, and conversion rate, respectively, that predator x2 preys on prey x1.Long-term ecological research studies show that three-species predator-prey models are fundamental building blocks of large scale ecosystems. However, it was only in the 1970s that some scholars began to study the dynamics of three-species predator-prey systems [5]. In particular, Hsu et al. [6] have classified all three-species predator-prey models into five types: two predators competing for one prey, one predator acting on two preys, food chain, food chain with omnivory, and food chain with cycle. Food chain architecture and strengths of species interactions are important determinants of trophic dynamics (see [7]). It is well known that tritrophic food chain model consists of one prey, one intermediate predator, and one top predator. Note that omnivory is a widespread mechanism in interacting populations. In [6], the authors investigated the following three-species predator-prey food chain model with an omnivory top predator:(2)dx1t=x1tr−a11x1t−a12x2t−a13x3tdt,dx2t=x2t−d2+a21x1t−a23x3tdt,dx3t=x3t−d3+a31x1t+a32x2tdt,where x1, x2, and x3 denote the number of prey, intermediate predator, and omnivorous top predator, respectively, r1 is the growth rate of prey, ri is the death rate of species xi i=2,3, a11 is the intraspecific competition rate of prey, a12, a13, and a23 are the capture rates, and a21, a31, and a32 denote the efficiency of food conversion. Model (2) describes that the intermediate predator preys on only the prey and the omnivorous top predator preys on both the prey and the intermediate predator. This is a general part of marine or terrestrial food web ecological systems. Based on model (2), Namba et al. [8] considered the intraspecific competition of the intermediate predator and the intraspecific competition of the top predator. Moreover, the authors demonstrated the stabilizing role of intraspecific competition among intermediate and top predators when the growth rate of prey species is adequate to support both the predator species. Furthermore, Sen et al. [9] investigated the following three-species Lotka–Volterra model with intraguild predation and mixed functional responses:(3)dx1t=x1tr−a11x1t−a12x2t−a13x3tdt,dx2t=x2t−d2+a21x1t−a22x2t−a23x3t1+βx2tdt,dx3t=x3t−d3+a31x1t+a32x2t1+βx2t−a33x3tdt,where x1, x2, and x3 denote the number of prey, intermediate predator, and omnivorous top predator, respectively. Obviously, in [9], the authors considered Holling type-II functional response between the intermediate predator and top predator and other functional responses were assumed to be linear. All meanings of the parameters are exact to or similar as those for (2) except the following. Here, aii is the intraspecific competition rate of species xi i=2,3 and β is the reciprocal of the half-saturation constant.Note that the three-species food web models (2) and (3) with the functional responses only depend on prey density. However, in fact, the predator has to search and compete for food and the ratio-dependent function of the prey and the predator is more suitable to substitute for the model with complicated interaction between the prey and predator. Then, the ratio-dependent type three-species food web model with omnivory is expressed in the form:(4)dx1t=x1tr−a1x1t−α12x2tx1t+β12x2t−α13x3tx1t+β13x3tdt,dx2t=x2t−d2−a2x2t+e12α12x1tx1t+β12x2t−α23x3tx2t+β23x3tdt,dx3t=x3t−d3−a3x3t+e13α13x1tx1t+β13x3t+e23α23x2tx2t+β23x3tdt,where x1t stands for the total number of prey at time t, while x2t and x3t represent the total number of intermediate predators and omnivorous top predators at time t, respectively. Here, r is the intrinsic growth rate of prey; di represents the mortality rate of predator xi i=2,3; ai stands for the intraspecific competition rate of species xi i=1,2,3; α12, β12, and e12 are the encounter rate, half-saturation constant, and conversion rate, respectively, that x2 preys on x1; α13, β13 and e13 stand for the same corresponding denotations that x3 preys on x1; and α23, β23, and e23 represent the same corresponding denotations that x3 preys on x2.As mentioned above, we notice that population models (1)–(4) are described by the deterministic model. This is valid only at the macroscopic scale, that is, the stochastic effects can be neglected or averaged out, in view of the law of large numbers. However, in the real world, populations are actually subject to the environmental fluctuations. Generally speaking, such fluctuations could be modeled by a colored noise. It has been noted that if the colored noise is not strongly correlated, then we can approximate the colored noise by a white noise w˙t, and the approximation works quite well (see [10]). It turns out that the white noise w˙t is formally regarded as the derivative of a Brownian motion wt, i.e., w˙t=dwt/dt (see [11]). As a result, the study of stochastic ecological dynamics model has already become one of the important subjects in biological mathematics.After taking the effect of randomly fluctuating environment into account, many researchers introduced stochastic environmental variation described by the Brownian motion into parameters in the deterministic model to establish the stochastic population model (see [12–15]). Liu and Bai [12] considered the optimal harvesting problem of a stochastic logistic model with time delay. In [13–15], the authors investigated the dynamics of stochastic predator-prey models. Ji et al. [13] discussed a stochastic predator-prey model with modified Leslie–Gowerand Holling-type II schemes. Jovanović and Krstić [14] investigated the extinction of a stochastic predator-prey model with the Allee effect on the prey. Liu and Jiang [15] considered the periodic solution and stationary distribution of stochastic predator-prey models with higher-order perturbation. In [16], considering that fluctuations in the environment would manifest themselves mainly as fluctuations in the intrinsic growth rate of the prey population and in the death rate of the predator population (see [17]), Ji et al. supposed parameters r and d2 in model (1) were perturbed with(5)r⟶r+σ1w˙1t,−d2⟶−d2+σ2w˙2t,where w1t and w2t are mutually independent Brownian motions and σi2 represents the intensity of white noise w˙it i=1,2. Moreover, they investigated the long time behavior of the following stochastic ratio-dependent prey-predator model:(6)dx1t=x1tr−a1x1t−α12x2tx1t+β12x2tdt+σ1x1tdw1t,dx2t=x2t−d2+e12α12x1tx1t+β12x2tdt+σ2x2tdw2t.Based on (6), Wu et al. [18] considered the corresponding nonautonomous stochastic ratio-dependent model. Lv et al. [19] introduced the intraspecific competition of the predator population, denoted by a2, into model (6). Nguyen and Ta [20] considered a corresponding nonautonomous stochastic ratio-dependent prey-predator model, in which the white noise makes the effect on both the growth rates of species and the intraspecific competition coefficient of the species.For the study of stochastic three-species models, consult [21–26] and the references therein. Geng et al. [21] investigated the stability of a stochastic one-predator-two-prey population model with time delay, while Liu et al. [22] studied the stability of a stochastic two-predator one-prey population model with time delay. In [23, 24], the authors discussed the dynamical behaviors of stochastic tri-trophic food-chain models. Li et al. [23] investigated the persistence and nonpersistence of a stochastic food-chain model, while Liu and Bai [24] considered the optimal harvesting problem of a stochastic three species food-chain model. Furthermore, in [25, 26], the stochastic three-species food-chain models with omnivory are discussed. Qiu and Deng [25] investigated the stationary distribution and global asymptotic stability of a stochastic food-web model with omnivory and linear functional response, while R. Liu and G. Liu [26] discussed the persistence in mean and extinction of a stochastic food-web model with intraguild predation and mixed functional responses. In [26], the authors considered Holling type-II functional response between the intermediate predator and the top predator and other functional responses were assumed to be linear.To the best of our knowledge, so far there is no investigation on the dynamics of the stochastic three-species food web model with omnivory and ratio-dependent functional response. The purpose of this paper is to make some contribution in this direction. Recall that parametersr, d2, and d3 in model (4) represent the intrinsic growth rate of the prey population, the death rate of the intermediate predator, and the death rate of the omnivorous top predator, respectively. As done in [16], in this paper, we may replace r, d2, and d3 in model (4), respectively, by(7)r⟶r+σ1w˙1t,−d2⟶−d2+σ2w˙2t,−d3⟶−d3+σ3w˙3t,where w˙it is the white noise and σi2 is the intensity of white noise w˙it i=1,2,3. Then, the stochastic three-species food web model with omnivory and ratio-dependent functional response took the following form:(8)dx1t=x1tr−a1x1t+α12x2tx1t+β12x2t−α13x3tx1t+β13x3tdt+σ1x1tdw1t,dx2t=x2t−d2−a2x2t+e12α12x1tx1t+β12x2t−α23x3tx2t+β23x3tdt+σ2x2tdw2t,dx3t=x3t−d3−a3x3t+e13α13x1tx1t+β13x3t+e23α23x2tx2t+β23x3tdt+σ3x3tdw3t,with x10,x20,x30=x10,x20,x30∈ℝ+3=x1,x2,x3∈ℝ3:xi>0i=1,2,3. All meanings of the parameters are exact to or similar as those for (4) except the following. Here, w=w1t,w2t,w3t:t≥0 represents the three-dimensional standard Brownian motion defined on a complete filtered probability space Ω,ℱ,ℱtt≥0,ℙ satisfying the usual conditions. σi2 represents the intensity of noise wit i=1,2,3. Throughout this paper, unless otherwise specified, we would rather assume that a1>0, a2≥0, a3≥0, α13≥0, α23≥0, β13>0, β23>0, e13>0, and e23>0.


## 2. Existence and Uniqueness of Positive Solution

In this section, we consider the existence of the positive solution for all times. Typically, conditions assuring the nonexplosion of the solution involve local Lipschitz continuity and a linear growth condition. In our case, we miss this last condition, so it is necessary to prove that the solution does not explode at a finite time. To prove the solution is positive and does not explode at a finite time, we use the stochastic comparison theorem. For simplicity, we introduce the following notations:(9)κ1=r−α12β12−α13β13;κ2=e12α12−d2−α23β23;κ3=e13α13+e23α23−d3.Theorem 1.
For any given initial valuex10,x20,x30∈ℝ+3, model (8) has unique global positive solution x1t,x2t,x3t for t≥0, that is, x1t,x2t,x3t∈ℝ+3 with probability one for t∈0,∞.Proof.
Consider the following system:(10)dX1t=r−a1eX1t+α12eX2teX1t+β12eX2t+α13eX3teX1t+β13eX3t−σ122dt+σ1dw1t,dX2t=−d2−a2eX2t+e12α12eX1teX1t+β12eX2t−α23eX3teX2t+β23eX3t−σ222dt+σ2dw2t,dX3t=−d3−a3eX3t+e13α13eX1teX1t+β13eX3t+e23α23eX2teX2t+β23eX3t−σ322dt+σ3dw3t,with initial value X10,X20,X30=lnx10,lnx20,lnx30. Obviously, the coefficients of (10) are locally Lipschitz continuous. Thus, there is a unique maximal local solution X1t,X2t,X3t of (10) for t∈0,τe, where τe denotes the explosion time. Let xit=eXit (i=1,2,3). Using Itô formula, it follows that x1t,x2t,x3t=eX1t,eX2t,eX3t is the unique positive local solution of (8) with initial value x10,x20,x30 for t∈0,τe.
Next, we show thatX1t,X2t,X3t is a global solution of (10), that is, τe=∞. Consider the following two stochastic differential systems:(11)dΦ1t=Φ1tr−a1Φ1tdt+σ1Φ1tdw1t,dΦ2t=Φ2te12α12−a2Φ2tdt+σ2Φ2tdw2t,dΦ3t=Φ3te13α13+e23α23−a3Φ3tdt+σ3Φ3tdw3t,with initial value Φ10,Φ20,Φ30=x10,x20,x30 and(12)dϕ1t=ϕ1tκ1−a1ϕ1tdt+σ1ϕ1tdw1t,dϕ2t=ϕ2tκ2−a2+e12α12β12ϕ1tϕ2tdt+σ2ϕ23tdw2t,dϕ3t=ϕ3tκ3−a3+e13α13β13ϕ1t+e23α23β23ϕ2tϕ3tdt+σ3ϕ3tdw3t,with initial value ϕ10,ϕ20,ϕ30=x10,x20,x30.
Thanks to Lemma 4.2 in [27], systems (11) and (12) can be explicitly solved as follows:(13)Φ1t=expr−σ12/2t+σ1w1t1/x10+a1∫0texpr−σ12/2s+σ1w1sds,Φ2t=expe12α12−σ22/2t+σ2w2t1/x20+a2∫0texpe12α12−σ22/2s+σ2w2sds,Φ3t=expe13α13+e23α23−σ32/2t+σ3w3t1/x30+a3∫0texpe13α13+e23α23−σ32/2s+σ3w3sds,ϕ1t=expκ1−σ12/2t+σ1w1t1/x10+a1∫0texpκ1−σ12/2s+σ1w1sds,ϕ2t=expκ2−σ22/2t+σ2w2t1/x20+∫0ta2+e12α12β12/ϕ1sexpκ2−σ22/2s+σ2w2sds,ϕ3t=expκ3−σ32/2t+σ3w3t1/x30+∫0ta3+e13α13β13/ϕ1s+e23α23β23/ϕ2sexpκ3−σ32/2s+σ3w3sds.
Note that the local solutionx1t,x2t,x3t is positive on 0,τe. Then, from the comparison theorem of stochastic differential equations (see Theorem 3.1 in [28]), it follows that for t∈0,τe(14)0<ϕit≤xit≤Φit,a.s.,i=1,2,3.
Thus, fort∈0,τe(15)lnϕit≤lnxit≤lnΦit,a.s.,i=1,2,3.
Sincelnϕit and lnΦit (i=1,2,3) exist for every t≥0, it follows that τe=∞. Thus, for any initial value X10,X20,X30=lnx10,lnx20,lnx30∈ℝ3, and (10) has a unique global solution X1t,X2t,X3t on 0,∞ a.s. Note that the coefficients of (8) are local Lipschitz continuous. Therefore, for any initial value x10,x20,x30∈ℝ+3, model (8) has a unique global positive solution x1t,x2t,x3t=eX1t,eX2t,eX3t on 0,∞ a.s. The proof is therefore complete.


## 3. Asymptotic Behaviors

Lemma 1 (see [13]).
Consider one-dimensional stochastic differential equation:(16)dxt=xta−bxtdt+σxtdwt,where a, b, and σ are positive constants and wt is standard Brownian motion. For any x0>0, let xt be the solution of equation (16) with initial value x0. If a>σ2/2, then(17)limt⟶∞lnxtt=0,limt⟶∞1t∫0txsds=a−σ2/2b,a.s.Theorem 2.
For anyx10,x20,x30∈ℝ+3, let x1t,x2t,x3t be the solution of model (8) with initial value x10,x20,x30. If a2>0, a3>0, and κi−σi2/2>0 i=1,2,3, then(18)limt⟶∞lnxitt=0,a.s.,i=1,2,3.Proof.
From Theorem 1, it follows that(19)ϕit≤xit≤Φit,a.s.,i=1,2,3.
Note thatϕ1t and Φ1t are the solutions of the following stochastic equations, respectively,(20)dϕ1t=ϕ1tκ1−a1ϕ1tdt+σ1ϕ1tdw1t,dΦ1t=Φ1tr−a1Φ1tdt+σ1Φ1tdw1t.with initial value x10>0. Obviously, from Lemma 1, it follows that if κ1−σ12/2>0:(21)limt⟶∞lnϕ1tt=0,limt⟶∞lnΦ1tt=0a.s.
This, together with (19), yields(22)limt⟶∞lnx1tt=0a.s.
Now, we showlimt⟶∞lnx2t/t=0 a.s. Note that, x2t≤Φ2t a.s. and Φ2t is the solution of equation(23)dΦ2t=Φ2te12α12−a2Φ2tdt+σ2Φ2tdw2t,with initial x20>0. Thus, from Lemma 1 and κ2−σ22/2>0, it follows that(24)limt⟶∞lnΦ2tt=0a.s.
Note thatlimt⟶∞lnϕ1t/t=0, a.s. Thus, for any ε>0, there exists T1>0 such that(25)e−εt≤ϕ1t≤eεt,fort≥T1.
By the strong law of large numbers of local martingales (see Theorem 1.3.4 in [11]), it follows that limt⟶∞σ2w2t/t=0 a.s. Thus, for any ε>0, there exists T2>0 such that(26)−εt≤σ2w2t≤εt,fort≥T2.
From the expression ofϕ2t that for any t≥T=T1∨T2, we have(27)1ϕ2t=1x2Te−κ2−σ22/2t−T−σ2w2t−w2T+a2∫Tte−κ2−σ22/2t−s−σ2w2t−w2sds+∫Tte12α12β12ϕ1se−κ2−σ22/2t−s−σ2w2t−w2sds≤1x2Te−κ2−σ22/2t−T+εt+T+a2∫Tte−κ2−σ22/2t−s+εt+sds+e12α12β12∫Tteεse−κ2−σ22/2t−s+εt+sds.
Hence, fromκ2−σ22/2>0, ε>0 and t≥T, it follows that(28)e−3εt+Tϕ2t≤1x2Te−κ2−σ22/2t−T−2εt+T+a2∫Tte−εt−se−εte−3εTe−κ2−σ22/2t−sds+e12α12β12e−3εT∫Tte−κ2−σ22/2t−s−2εt−sds=1x2Te−κ2−σ22/2t−T−2εt+T+a2κ2+ε−σ22/2e−εte−3εTe−κ2+ε−σ22/2T+e12α12β12κ2+2ε−σ22/2e−3εTe−κ2+2ε−σ22/2T≤1x2T+a2κ2−σ22/2+e12α12β12κ2−σ22/2=˙K1.
That is1/ϕ2t≤K1e3εt+T a.s., for t≥T. Then, −lnϕ2t≤lnK1+3εt+T. Thus, for any ε>0,(29)liminft⟶∞lnϕ2tt≥0a.s.
Consequently,(30)0≤liminft⟶∞lnϕ2tt≤liminft⟶∞lnx2tt≤limsupt⟶∞lnx2tt≤limt⟶∞lnΦ2tt=0a.s.
In addition,(31)lim supt⟶∞lnϕ2tt≤lim  supt⟶∞lnx2tt≤0a.s.
Thus,(32)limt⟶∞lnϕ2tt=0,limt⟶∞lnx2tt=0,a.s.
Similarly, ifκ1−σ12/2>0, κ2−σ22/2>0, and κ3−σ32/2>0, then(33)limt⟶∞lnx3tt=0a.s.
The proof is therefore complete.


## 4. Persistence in Mean and Extinction

In this section, we show that under some conditions, model (8) is persistent in mean and extinct.Theorem 3.
Suppose thata2>0, a3>0, and κi−σi2/2>0 i=1,2,3. Then, for any x10,x20,x30∈ℝ+3, the solution x1t,x2t,x3t of model (8) with initial value x10,x20,x30 obeys(34)lim inft⟶∞1t∫0tx1sds≥κ1−σ12/2a1a.s.,lim  inft⟶∞1t∫0ta2x2s+e12α12β12x2sx1sds≥κ2−σ222a.s.,lim  inft⟶∞1t∫0ta3x3s+e13α13β13x3sx1s+e23α23β23x3sx2sds≥κ3−σ322a.s.Proof.
For preyx1, from Theorem 1, it follows that(35)ϕ1t≤x1ta.s.,and ϕ1t is the solution of the following stochastic equation:(36)dϕ1t=ϕ1tκ1−a1ϕ1tdt+σ1ϕ1tdw1t,with initial value x10>0. Obviously, from Lemma 1, it follows that if κ1−σ12/2>0,(37)limt⟶∞1t∫0tϕ1sds=κ1−σ12/2a1a.s.
This, together with (35), yields(38)lim inft⟶∞1t∫0tx1sds≥limt⟶∞1t∫0tϕ1sds=κ1−σ12/2a1>0a.s.
For intermediate predatorx2, using Itô formula, it follows that(39)lnx2t=∫0t−d2−a2x2s+e12α12x1sx1s+β12x2s−α23x3sx2s+β23x3s−σ222ds+σ2w2t+lnx20≥κ2−σ222t−a2∫0tx2sds−∫0te12α12β12x2sx1s+β12x2sds+σ2w2t+lnx20.
Hence,(40)1t∫0ta2x2s+e12α12β12x2sx1sds≥κ2−σ222+σ2w2tt+lnx20t−lnx2tt.
By the strong law of numbers of local martingales and Theorem 2, we get(41)lim inft⟶∞1t∫0ta2x2s+e12α12β12x2sx1sds≥κ2−σ222a.s.
For the omnivorous top predatorx3, it follows from Itô formula that(42)lnx3t=∫0t−d3−a3x3s+e13α13x1sx1s+β13x3s+e23α23x2sx2s+β23x3s−σ322ds+σ3w3t+lnx30≥κ3−σ322t−∫0ta3x3s+e13α13β13x3sx1s+β13x3s+e23α23β23x3sx2s+β23x3sds+σ3w3t+lnx30.
Hence,(43)1t∫0ta3x3s+e13α13β13x3sx1s+e23α23β23x3sx2sds≥κ3−σ322+σ3w3tt+lnx30t−lnx3tt.
By the strong law of numbers of local martingales and Theorem 2, we get(44)lim inft⟶∞1t∫0ta3x3s+e13α13β13x3sx1s+e23α23β23x3sx2sds≥κ3−σ322a.s.
The proof is therefore complete.Theorem 4.
Suppose thatr−σ12/2<0, e12α12−d2−σ22/2<0 and e13α13+e23α23−d3−σ32/2<0. Then, for any initial value x10,x20,x30∈ℝ+3, model (8) is extinct exponentially with probability one.Proof.
From Itô formula, it follows that(45)lnx1t≤r−σ122t+σ1w1t+lnx10,lnx2t≤e12α12−d2−σ222t+σ2w2t+lnx20,lnx3t≤e13α13+e23α23−d3−σ322t+σ3w3t+lnx30.
Note thatlimt⟶∞σiwit/t+lnxi0/t=0 i=1,2,3 and r−σ12/2<0, e12α12−d2−σ22/2<0, and e13α13+e23α23−d3−σ32/2<0. Then,(46)lim supt⟶∞lnx1tt≤r−σ122<0a.s.,lim supt⟶∞lnx2tt≤e12α12−d2−σ222<0a.s.,lim  supt⟶∞lnx3tt≤e13α13+e23α23−d3−σ322<0a.s.
Therefore, model (8) is extinct exponentially. The proof is complete.


## 5. Stochastic Permanence

In this section, we discuss the stochastic permanence of model (8). The definition of stochastic permanence and stochastically ultimately boundness of model (8) were introduced in the literature [29, 30] as follows.Definition 1 (see [29, 30]).
Model (8) is called stochastically ultimate bounded, if for any ε∈0,1, there exist three positive constants H1=H1ε, H2=H2ε, and H3=H3ε such that the solution x1t,x2t,x3t of model (8) with any initial value x10,x20,x30∈ℝ+3 has the property that(47)lim supt⟶∞ℙxit>Hi<ε,i=1,2,3.Definition 2 (see [29, 30]).
Model (8) is said to be stochastically permanent, if for any ε∈0,1, there exist positive constants δi=δiε, Hi=Hiε, and δi<Hi i=1,2,3, such that the solution x1t,x2t,x3t of model (8) with any initial value x10,x20,x30∈ℝ+3 has the property that(48)lim inft⟶∞ℙxit≤Hi≥1−ε,lim inft⟶∞ℙxit≥δi≥1−ε,i=1,2,3.
It is obvious that if stochastic model (8) is stochastically permanent, its solutions must be stochastically ultimately bounded.


### 5.1. Boundness

In this subsection, we investigate the stochastically ultimate boundness of model (8) in two different ways.Lemma 2 (see [31]).
For any positive constantsp, m, and n, the Bernoulli equation(49)dxtdt=pmxt−pnx1+1/pt,with the initial value x0=x0>0, has the solution(50)xt=mn1−e−mt+m/nx0−1/pe−mtp.Theorem 5.
For anyx10,x20,x30∈ℝ+3, let x1t,x2t,x3t be the solution of model (8) with initial value x10,x20,x30. If a2>0 and a3>0, then for any p≥0,(51)lim supt⟶∞Ex1pt≤r+p/2σ12a1p,lim supt⟶∞Ex2pt≤e12α12+p/2σ22a2p,lim supt⟶∞Ex3pt≤e13α13+e23α23+p/2σ32a3p.That is, the solution of model (8) is uniformly bounded in the pth moment.Proof.
ForΦ1 in system (11), applying Itô formula to Φ1p leads to(52)Φ1pt=x10p+∫0tpΦ1psr+p−12σ12−a1Φ1sds+∫0tpσ1Φ1psdw1s.
Taking the expectation on both sides of the above equation, we have(53)EΦ1pt=x10p+E∫0tpΦ1psr+p−12σ12−a1Φ1sds.
Then, using the Höder inequality, it follows that(54)dEΦ1ptdt=pr+p−12σ12EΦ1pt−pa1EΦ1p+1t≤pr+p2σ12EΦ1pt−pa1EΦ1pt1+1/p=˙pb1EΦ1pt−pa1EΦ1pt1+1/p.
From Lemma 2 and the comparison theorem, it follows that(55)EΦ1pt≤b1a11−e−b1t+b1/a1x10−1e−b1tp.
Note thatb1=r+p/2σ12>0. Thus,(56)limsupt⟶∞EΦ1pt≤r+p/2σ12a1p.
By a similar the discussion as inΦ1t, we also know that(57)limsupt⟶∞EΦ2pt≤e12α12+p−1/2σ22a2p,limsupt⟶∞EΦ3pt≤e13α13+e23α23+p−1/2σ32a3p.
From Theorem 1, it follows that0<xit≤Φit a.s. i=1,2,3. Then, for any p≥0, we have(58)0<Exipt≤EΦipt,i=1,2,3.
Now Theorem 5 follows immediately from the above analysis. The proof is complete.Theorem 6.
For anyx10,x20,x30∈ℝ+3, let x1t,x2t,x3t be the solution of model (8) with initial value x10,x20,x30. Then,(59)limsupt⟶∞Ex1t≤K2dLe12e23+e13,limsupt⟶∞Ex2t≤K2dLe23,limsupt⟶∞Ex3t≤K2dL,where dL=d2∧d3 and K2=e12e23+e13r+dL2/4a1.Proof.
DefineHt=e12e23+e13x1t+e23x2t+x3t. By Itô formula, we have(60)dHt=e12e23+e13rx1−a1x12−α12x1x2x1+β12x2−α13x1x3x1+β13x3dt+σ1x1dw1t+e23−d2x2−a2x22+e12α12x1x2x1+β12x2−α23x2x3x2+β23x3dt+σ2x2dw2t+−d3x3−a3x32+e13α13x1x3x1+β13x3+e23α23x2x3x2+β23x3dt+σ3x3dw3t=−e23d2x2−d3x3+e12e23+e13rx1−e12e23+e13a1x12−e23a2x22−a3x32−e12e23α13x1x3x1+β13x3−e13α12x1x2x1+β12x2dt+e12e23+e13σ1x1dw1t+e23σ2x2dw2t+σ3x3dw3t.
Integrating it from 0 to t and taking expectation yields(61)EHt=H0+E∫0t−e23d2x2−d3x3+e12e23+e13rx1−e12e23+e13a1x12−e23a2x22−a3x32−e12e23α13x1x3x1+β13x3−e13α12x1x2x1+β12x2ds.
Thus, using the Hölder inequality yields(62)dEHtdt=−e23d2Ex2t−d3Ex3t+e12e23+e13rEx1t−e12e23+e13a1Ex12t−e23a2Ex22t−a3Ex32t−e12e23α13Ex1tx3tx1t+β13x3t−e13α12Ex1tx2tx1t+β12x2t≤−dLe12e23+e13Ex1t−dLe23Ex2t−dLEx3t+e12e23+e13r+dLEx1t−e12e23+e13a1Ex1t2=e12e23+e13r+dLEx1t−a1Ex1t2−dLEHt.
It is clear that quadratic functiongx=e12e23+e13r+dLx−a1x2 reaches its maximum value at x=r+dL/2a1>0. Thus, gmax=e12e23+e13r+dL2/4a1=˙K2. Therefore, we have(63)dEHtdt≤K2−dLEHt.
Then, by the comparison theorem, we have(64)0≤lim  supt⟶∞EHt≤K2dL.
Note that the solution of model (8) is positive. Thus,(65)limsupt⟶∞Ex1t≤K2dLe12e23+e13,limsupt⟶∞Ex2t≤K2dLe23,limsupt⟶∞Ex3t≤K2dL.
The proof is therefore complete.Theorem 7.
Model (8) is stochastically ultimate bounded.Proof.
Letx1t,x2t,x3t be solution of (8) with any initial value x10,x20,x30∈ℝ+3. For any ε∈0,1, let H1=K2/dLe12e23+e13ε+1, H2=K2/dLe23ε+1, and H3=K2/dLε+1. Then, by Chebyshev’s inequality(66)ℙxit>Hi≤ExitHi,i=1,2,3.
Hence, from Theorem 6(67)limsupt⟶∞ℙxit>Hi≤limsupt⟶∞ExitHi<ε,i=1,2,3.
The proof is therefore complete.


### 5.2. Stochastic Permanence

In this section, we give some sufficient conditions to guarantee that model (8) is stochastically permanent. Denote γi=˙κi−σi2 i=1,2,3. Define(68)u1t=1ϕ1t,u2t=1ϕ2t,u3t=1ϕ3t.By the Itô formula, we have(69)du1t=a1−γ1u1tdt−σ1u1tdw1t,du2t=a2−γ2u2t+e12α12β12u1tdt−σ2u2tdw2t,du3t=a3−γ3u3t+e13α13β13u1t+e23α23β23u2tdt−σ3u3tdw3t,with initial value u10,u20,u30=1/x10,1/x20,1/x30∈ℝ+3.Lemma 3.
Letx1t,x2t,x3t be the solution of model (8) with any initial value x10,x20,x30∈ℝ+3. If γi>0 i=1,2,3, then(70)limsupt⟶∞E1xit≤Mi,i=1,2,3,where M1=a1/γ1, M2=a2/γ2+a1e12α12β12/γ1γ2, and M3=a3/γ3+a1e13α13β13/γ1γ3++e23α23β23/γ3a2/γ2+a1e12α12β12/γ1γ2.Proof.
First, integrating both sides of the first equation of (69) from 0 to t yields(71)u1t=1x10+∫0ta1−γ1u1sds−∫0tσ1u1sdw1s.
Taking the expectation on both sides of the above equation, we have(72)Eu1t=1x10+E∫0ta1−γ1u1sds.
Thus,(73)dEu1tdt=a1−γ1Eu1t,with initial value Eu10=1/x10. By a simple computation, we can get(74)Eu1t=1x10e−γ1t+a1γ11−e−κ1t.
This, together withγ1>0, yields(75)limt⟶∞Eu1t=a1γ1=M1.
Next, integrating both sides of the second equation of system (69) from 0 to t yields(76)u2t=1x20+∫0ta2−γ2u2s+e12α12β12u1sds−∫0tσ2u2sdw2s.
Taking the expectation on both sides of the above equation, we have(77)Eu2t=1x20+E∫0ta2−γ2u2s+e12α12β12u1sds.
Thus,(78)dEu2tdt=a2−γ2Eu2t+e12α12β12Eu1t,with initial value Eu20=1/x20. By a simple computation, we can get(79)Eu2t=1x20e−γ2t+a2γ21−e−γ2t+e12α12β12∫0te−γ2t−sEu1sds.
From (75), it follows that(80)limt⟶∞∫0te−γ2t−sEu1sds=limt⟶∞∫0teγ2sEu1sdseγ2t=limt⟶∞Eu1tγ2=a1γ1γ2.
This, together with (79) yields(81)limt⟶∞Eu2t≤a2γ2+a1e12α12β12γ1γ2.
At last, integrating both sides of the third equation of system (69) from 0 to t and taking the expectation, we have(82)Eu3t=1x30+E∫0ta3−γ3u3s+e13α13β13u1s+e23α23β23u2sds.
Thus,(83)dEu3tdt=a3−γ3Eu3t+e13α13β13Eu1t+e23α23β23Eu2t,with initial value Eu30=1/x30. By a simple computation, we can get(84)Eu3t=1x30e−γ3t+a3γ31−e−γ3t+e13α13β13∫0te−γ3t−sEu1sds+e23α23β23∫0te−γ3t−sEu2tds.
It follows from (75), (81), and (84) that(85)limt⟶∞Eu3t=a3γ3+e13α13β13limt⟶∞∫0te−γ3t−sEu1sds+e23α23β23limt⟶∞∫0te−γ3t−sEu2tds=a3γ3+a1e13α13β13γ1γ3+e23α23β23γ3a2γ2+a1e12α12β12γ1γ2.
From the comparison theorem of stochastic differential equations, it follows that(86)1xit≤1ϕit=uit,i=1,2,3.
Now, Lemma 3 follows immediately from the above analysis. The proof is complete.Theorem 8.
Ifγi>0 i=1,2,3, then model (8) is stochastically permanent.Proof.
Letx1t,x2t,x3t be solution of (8) with initial value x10,x20,x30∈ℝ+3. For any ε∈0,1, let δi=ε/Mi i=1,2,3, then(87)ℙxit<δi=ℙ1xit>1δi≤E1/xit1/δi=δiE1xit,i=1,2,3.
Thus, from Lemma 3, it follows that(88)lim supt⟶∞ℙxit<δi≤lim supt⟶∞δiE1xit≤ε,i=1,2,3.
This implies(89)liminft⟶∞ℙxit≥δi≥1−ε,i=1,2,3.
Letε∈0,1 be sufficiently small such that δi<Hi. From (67) and Definition 2, model (8) is stochastically permanent. The proof is therefore complete.


## 5.1. Boundness

In this subsection, we investigate the stochastically ultimate boundness of model (8) in two different ways.Lemma 2 (see [31]).
For any positive constantsp, m, and n, the Bernoulli equation(49)dxtdt=pmxt−pnx1+1/pt,with the initial value x0=x0>0, has the solution(50)xt=mn1−e−mt+m/nx0−1/pe−mtp.Theorem 5.
For anyx10,x20,x30∈ℝ+3, let x1t,x2t,x3t be the solution of model (8) with initial value x10,x20,x30. If a2>0 and a3>0, then for any p≥0,(51)lim supt⟶∞Ex1pt≤r+p/2σ12a1p,lim supt⟶∞Ex2pt≤e12α12+p/2σ22a2p,lim supt⟶∞Ex3pt≤e13α13+e23α23+p/2σ32a3p.That is, the solution of model (8) is uniformly bounded in the pth moment.Proof.
ForΦ1 in system (11), applying Itô formula to Φ1p leads to(52)Φ1pt=x10p+∫0tpΦ1psr+p−12σ12−a1Φ1sds+∫0tpσ1Φ1psdw1s.
Taking the expectation on both sides of the above equation, we have(53)EΦ1pt=x10p+E∫0tpΦ1psr+p−12σ12−a1Φ1sds.
Then, using the Höder inequality, it follows that(54)dEΦ1ptdt=pr+p−12σ12EΦ1pt−pa1EΦ1p+1t≤pr+p2σ12EΦ1pt−pa1EΦ1pt1+1/p=˙pb1EΦ1pt−pa1EΦ1pt1+1/p.
From Lemma 2 and the comparison theorem, it follows that(55)EΦ1pt≤b1a11−e−b1t+b1/a1x10−1e−b1tp.
Note thatb1=r+p/2σ12>0. Thus,(56)limsupt⟶∞EΦ1pt≤r+p/2σ12a1p.
By a similar the discussion as inΦ1t, we also know that(57)limsupt⟶∞EΦ2pt≤e12α12+p−1/2σ22a2p,limsupt⟶∞EΦ3pt≤e13α13+e23α23+p−1/2σ32a3p.
From Theorem 1, it follows that0<xit≤Φit a.s. i=1,2,3. Then, for any p≥0, we have(58)0<Exipt≤EΦipt,i=1,2,3.
Now Theorem 5 follows immediately from the above analysis. The proof is complete.Theorem 6.
For anyx10,x20,x30∈ℝ+3, let x1t,x2t,x3t be the solution of model (8) with initial value x10,x20,x30. Then,(59)limsupt⟶∞Ex1t≤K2dLe12e23+e13,limsupt⟶∞Ex2t≤K2dLe23,limsupt⟶∞Ex3t≤K2dL,where dL=d2∧d3 and K2=e12e23+e13r+dL2/4a1.Proof.
DefineHt=e12e23+e13x1t+e23x2t+x3t. By Itô formula, we have(60)dHt=e12e23+e13rx1−a1x12−α12x1x2x1+β12x2−α13x1x3x1+β13x3dt+σ1x1dw1t+e23−d2x2−a2x22+e12α12x1x2x1+β12x2−α23x2x3x2+β23x3dt+σ2x2dw2t+−d3x3−a3x32+e13α13x1x3x1+β13x3+e23α23x2x3x2+β23x3dt+σ3x3dw3t=−e23d2x2−d3x3+e12e23+e13rx1−e12e23+e13a1x12−e23a2x22−a3x32−e12e23α13x1x3x1+β13x3−e13α12x1x2x1+β12x2dt+e12e23+e13σ1x1dw1t+e23σ2x2dw2t+σ3x3dw3t.
Integrating it from 0 to t and taking expectation yields(61)EHt=H0+E∫0t−e23d2x2−d3x3+e12e23+e13rx1−e12e23+e13a1x12−e23a2x22−a3x32−e12e23α13x1x3x1+β13x3−e13α12x1x2x1+β12x2ds.
Thus, using the Hölder inequality yields(62)dEHtdt=−e23d2Ex2t−d3Ex3t+e12e23+e13rEx1t−e12e23+e13a1Ex12t−e23a2Ex22t−a3Ex32t−e12e23α13Ex1tx3tx1t+β13x3t−e13α12Ex1tx2tx1t+β12x2t≤−dLe12e23+e13Ex1t−dLe23Ex2t−dLEx3t+e12e23+e13r+dLEx1t−e12e23+e13a1Ex1t2=e12e23+e13r+dLEx1t−a1Ex1t2−dLEHt.
It is clear that quadratic functiongx=e12e23+e13r+dLx−a1x2 reaches its maximum value at x=r+dL/2a1>0. Thus, gmax=e12e23+e13r+dL2/4a1=˙K2. Therefore, we have(63)dEHtdt≤K2−dLEHt.
Then, by the comparison theorem, we have(64)0≤lim  supt⟶∞EHt≤K2dL.
Note that the solution of model (8) is positive. Thus,(65)limsupt⟶∞Ex1t≤K2dLe12e23+e13,limsupt⟶∞Ex2t≤K2dLe23,limsupt⟶∞Ex3t≤K2dL.
The proof is therefore complete.Theorem 7.
Model (8) is stochastically ultimate bounded.Proof.
Letx1t,x2t,x3t be solution of (8) with any initial value x10,x20,x30∈ℝ+3. For any ε∈0,1, let H1=K2/dLe12e23+e13ε+1, H2=K2/dLe23ε+1, and H3=K2/dLε+1. Then, by Chebyshev’s inequality(66)ℙxit>Hi≤ExitHi,i=1,2,3.
Hence, from Theorem 6(67)limsupt⟶∞ℙxit>Hi≤limsupt⟶∞ExitHi<ε,i=1,2,3.
The proof is therefore complete.


## 5.2. Stochastic Permanence

In this section, we give some sufficient conditions to guarantee that model (8) is stochastically permanent. Denote γi=˙κi−σi2 i=1,2,3. Define(68)u1t=1ϕ1t,u2t=1ϕ2t,u3t=1ϕ3t.By the Itô formula, we have(69)du1t=a1−γ1u1tdt−σ1u1tdw1t,du2t=a2−γ2u2t+e12α12β12u1tdt−σ2u2tdw2t,du3t=a3−γ3u3t+e13α13β13u1t+e23α23β23u2tdt−σ3u3tdw3t,with initial value u10,u20,u30=1/x10,1/x20,1/x30∈ℝ+3.Lemma 3.
Letx1t,x2t,x3t be the solution of model (8) with any initial value x10,x20,x30∈ℝ+3. If γi>0 i=1,2,3, then(70)limsupt⟶∞E1xit≤Mi,i=1,2,3,where M1=a1/γ1, M2=a2/γ2+a1e12α12β12/γ1γ2, and M3=a3/γ3+a1e13α13β13/γ1γ3++e23α23β23/γ3a2/γ2+a1e12α12β12/γ1γ2.Proof.
First, integrating both sides of the first equation of (69) from 0 to t yields(71)u1t=1x10+∫0ta1−γ1u1sds−∫0tσ1u1sdw1s.
Taking the expectation on both sides of the above equation, we have(72)Eu1t=1x10+E∫0ta1−γ1u1sds.
Thus,(73)dEu1tdt=a1−γ1Eu1t,with initial value Eu10=1/x10. By a simple computation, we can get(74)Eu1t=1x10e−γ1t+a1γ11−e−κ1t.
This, together withγ1>0, yields(75)limt⟶∞Eu1t=a1γ1=M1.
Next, integrating both sides of the second equation of system (69) from 0 to t yields(76)u2t=1x20+∫0ta2−γ2u2s+e12α12β12u1sds−∫0tσ2u2sdw2s.
Taking the expectation on both sides of the above equation, we have(77)Eu2t=1x20+E∫0ta2−γ2u2s+e12α12β12u1sds.
Thus,(78)dEu2tdt=a2−γ2Eu2t+e12α12β12Eu1t,with initial value Eu20=1/x20. By a simple computation, we can get(79)Eu2t=1x20e−γ2t+a2γ21−e−γ2t+e12α12β12∫0te−γ2t−sEu1sds.
From (75), it follows that(80)limt⟶∞∫0te−γ2t−sEu1sds=limt⟶∞∫0teγ2sEu1sdseγ2t=limt⟶∞Eu1tγ2=a1γ1γ2.
This, together with (79) yields(81)limt⟶∞Eu2t≤a2γ2+a1e12α12β12γ1γ2.
At last, integrating both sides of the third equation of system (69) from 0 to t and taking the expectation, we have(82)Eu3t=1x30+E∫0ta3−γ3u3s+e13α13β13u1s+e23α23β23u2sds.
Thus,(83)dEu3tdt=a3−γ3Eu3t+e13α13β13Eu1t+e23α23β23Eu2t,with initial value Eu30=1/x30. By a simple computation, we can get(84)Eu3t=1x30e−γ3t+a3γ31−e−γ3t+e13α13β13∫0te−γ3t−sEu1sds+e23α23β23∫0te−γ3t−sEu2tds.
It follows from (75), (81), and (84) that(85)limt⟶∞Eu3t=a3γ3+e13α13β13limt⟶∞∫0te−γ3t−sEu1sds+e23α23β23limt⟶∞∫0te−γ3t−sEu2tds=a3γ3+a1e13α13β13γ1γ3+e23α23β23γ3a2γ2+a1e12α12β12γ1γ2.
From the comparison theorem of stochastic differential equations, it follows that(86)1xit≤1ϕit=uit,i=1,2,3.
Now, Lemma 3 follows immediately from the above analysis. The proof is complete.Theorem 8.
Ifγi>0 i=1,2,3, then model (8) is stochastically permanent.Proof.
Letx1t,x2t,x3t be solution of (8) with initial value x10,x20,x30∈ℝ+3. For any ε∈0,1, let δi=ε/Mi i=1,2,3, then(87)ℙxit<δi=ℙ1xit>1δi≤E1/xit1/δi=δiE1xit,i=1,2,3.
Thus, from Lemma 3, it follows that(88)lim supt⟶∞ℙxit<δi≤lim supt⟶∞δiE1xit≤ε,i=1,2,3.
This implies(89)liminft⟶∞ℙxit≥δi≥1−ε,i=1,2,3.
Letε∈0,1 be sufficiently small such that δi<Hi. From (67) and Definition 2, model (8) is stochastically permanent. The proof is therefore complete.


## 6. Stationary Distribution and Ergodicity

In this section, we will show that there is an ergodic stationary distribution for the solution of (8). For the completeness of the paper, in this section, we list some theories about stationary distribution (see [32]). Let Xt be a homogeneous Markov process in Ed (denotes d-dimensional Euclidean space), described by the following stochastic differential equation:(90)dXt=bXtdt+gXtdWt,X0=X0.The diffusion matrix of the processXt is defined as JX=gXgTX=aijX.Definition 3 (see [32]).
Letℙt,X,⋅ be the probability measure induced by Xt with initial value X0=X0. That is, ℙt,X0,A=ℙXt∈A∣X0=X0, for any Borel set A∈ℬℝ+d. If there exists a probability measure μ⋅ such that limt⟶∞ℙt,X0,A=μA for all X0∈ℝ+d and A∈ℬℝ+d, then we say that stochastic differential equation (90) has a stationary distribution μ⋅.Lemma 4 (see [33, 34]).
Assume that there exists a bounded domainD⊂Ed with regular boundary Γ and(i)
(A1) There is a positive numberM such that ∑i,j=1daijXξiξj≥Mξ2, X∈D, and ξ∈ℝd;(ii)
(A2) There exists a nonnegativeC2-function V such that there exists a positive constant C, such that(91)LV≤−Cfor  anyX∈EdD.Then, the Markov processXt has a unique ergodic stationary distribution μ⋅. Moreover, if f⋅ is a function integrable with respect to the measure μ, then(92)ℙlimT⟶∞1T∫0TfXtdt=∫Edfxμdx=1.LetXt=x1t,x2t,x3tT, gX=diagσ1x1,σ2x2,σ3x3, Wt=w1t,w2t,w3tT, and(93)bX=x1r−a1x1−α12x2x1+β12x2−α13x3x1+β13x3x2−d2−a2x2+e12α12x1x1+β12x2−α23x3x2+β23x3x3−d3−a3x3+e13α13x1x1+β13x3+e23α23x2x2+β23x3.Then, system (90) reduces to model (8) with diffusion matrix JX=diagσ12x12,σ22x22,σ32x32.Theorem 9.
Ifa2>0, a3>0, κ1−e12α12β12−e13α13β13−σ12>0, κ2−e23α23β23−σ22>0, and κ3−σ32>0, then for any x10,x20,x30∈ℝ+3, model (8) has a stationary distribution and the solutions have ergodic property.Proof.
DefineC2-function V1: ℝ+3⟶ℝ+ by(94)V1X=x1+x2+x3,for X=x1,x2,x3∈ℝ+3. By Itô formula, we have(95)LV1X=x1r−a1x1−α12x2x1+β12x2−α13x3x1+β13x3+x2−d2−a2x2+e12α12x1x1+β12x2−α23x3x2+β23x3+x3−d3−a3x3+e13α13x1x1+β13x3+e23α23x2x2+β23x3≤−a1x12+rx1−a2x22+e12α12−d2x2−a3x32+e13α13+e23α23−d3x3.
DefineC2-function V2: ℝ+3⟶ℝ+ by(96)V2X=x1−1+x2−1+x3−1,for X=x1,x2,x3∈ℝ+3. By Itô formula, we have(97)LV2X=−x1−1r−a1x1−α12x2x1+β12x2−α13x3x1+β13x3+σ12x1−1−x2−1−d2−a2x2+e12α12x1x1+β12x2−α23x3x2+β23x3+σ22x2−1−x3−1−d3−a3x3+e13α13x1x1+β13x3+e23α23x2x2+β23x3+σ32x3−1≤−x1−1κ1−a1x1+σ12x1−1−x2−1κ2−a2x2−e12α12β12x2x1+β12x2+σ22x2−1−x3−1κ3−a3x3−e13α13β13x3x1+β13x3−e23α23β23x3x2+β23x3+σ32x3−1≤−x1−1κ1−a1x1+σ12x1−1−x2−1κ2−a2x2−e12α12β12x2x1+σ22x2−1−x3−1κ3−a3x3−e13α13β13x3x1−e23α23β23x3x2+σ32x3−1=−κ1−e12α12β12−e13α13β13−σ12x1−1+a1−κ2−e23α23β23−σ22x2−1+a2−κ3−σ32x3−1+a3.
LetVX=V1X+V2X. Then,(98)LVX≤−a1x12+rx1−κ1−e12α12β12−e13α13β13−σ12x1−1+a1−a2x22+e12α12−d2x2−κ2−e23α23β23−σ22x2−1+a2−a3x32+e13α13+e23α23−d3x3−κ3−σ32x3−1+a3=fx1+gx2+hx3,where(99)fx1=−a1x12+rx1−κ1−e12α12β12−e13α13β13−σ12x1−1+a1,gx2=−a2x22+e12α12−d2x2−κ2−e23α23β23−σ22x2−1+a2,hx3=−a3x32+e13α13+e23α23−d3x3−κ3−σ32x3−1+a3.
Clearly,fx1, gx2, and hx3 have upper bound on ℝ+. Denote(100)fu=supx1∈ℝ+fx1,gu=supx2∈ℝ+gx2,hu=supx3∈ℝ+hx3.
Fromκ1−e12α12β12−e13α13β13−σ12>0, it follows that(101)LVX≤fx1+gx2+hx3≤fx1+gu+hu⟶−∞,a.s.x1⟶0+orx1⟶+∞.
Similarly, fromκ2−e23α23β23−σ22>0 and κ3−σ32>0, we have(102)LVX≤fx1+gx2+hx3≤fu+gx2+hu⟶−∞,a.s.x2⟶0+orx2⟶+∞,LVX≤fx1+gx2+hx3≤fu+gu+hx3⟶−∞,a.s.x3⟶0+orx3⟶+∞.
Consequently, there existsρ>0 (sufficiently small) such that (103)LVX≤−1,for allx1,x2,x3∈ℝ+3D,where(104)D=x1,x2,x3∈ℝ+3∣ρ<x1<1ρ,ρ<x2<1ρ,ρ<x3<1ρ⊂ℝ+3.
Hence, (A2) in Lemma 4 is satisfied.
Denoteσ2=σ12∧σ22∧σ32. Then, for any X=x1,x2,x3∈D and ξ=ξ1,ξ2,ξ2∈ℝ3, we have(105)∑i,j=13aijXξiξj=σ12x12ξ12+σ22x22ξ22+σ32x32ξ32≥Mξ2,where M=ρ2σ2. Thus, condition (A1) of Lemma 4 holds. According to Lemma 4, model (8) is ergodic and admits a unique stationary distribution. The proof is therefore complete.


## 7. Application of Main Results

In this section, we present the application of the main results in some special models.


### 7.1. Two Species Predator-Prey Model

Letα13=α23=0. Then, the first two equations of (8) form the following closed two-population system:(106)dx1t=x1tr−a1x1t+α12x2tx1t+β12x1tdt+σ1x1tdw1t,dx2t=x2t−d2−a2x2t+e12α12x1tx1t+β12x2tdt+σ2x2tdw2t,with initial value x10,x20∈ℝ+2. This is also the stochastic predator-prey model discussed in [19]. From Theorems 3 and 4, we have the following result.Corollary 1.
Letx1t,x2t be solution of model (106) with initial value x10,x20∈ℝ+2.(i)
Ifa2>0, r−α12/β12−σ12/2>0, and e12α12−d2−σ22/2>0, then(107)liminft⟶∞1t∫0tx1sds≥r−α12/β12−σ12/2a1>0a.s.,liminft⟶∞1t∫0ta2x2s+e12α12β12x2sx1sds≥e12α12−d2−σ222>0a.s.(ii)
Ifr−σ12/2<0 and e12α12−d2−σ22/2<0, then (106) is extinct exponentially with probability one.Remark 1.
It is clear that Corollary 1 is consistent with Theorems 7 and 8 in [19]. Moreover, from Theorems 3 and 4, the persistence in mean and extinction conditions of the three-species model (8) are more complicated. Thus, our work can be seen as the extension of [19].
For model (106), similar to the proof of Theorem 6 (denote H=e12x1+x2), we have(108)limsupt⟶∞Ex1t≤K2′e12d2,limsupt⟶∞Ex2t≤K2′d2,where K2′=e12r+d22/4a1. Furthermore, from Theorems 7–9, for model (106), we have the following result.Corollary 2.
(i)
Model (106) is stochastically ultimate bounded(ii)
Ifr−α12/β12−σ12>0 and e12α12−d2−σ22>0, then model (106) is stochastically permanent(iii)
Ifa2>0, r−α12/β12−e12α12β12−σ12>0, and e12α12−d2−σ22>0, then for any x10,x20∈ℝ+2, model (106) has a stationary distribution and the solutions have ergodic propertyIf we do not consider the intraspecific competition of the predator, i.e.,a2=0 in model (106), then model (6) is available. From Theorems 4, 7, and 8, for model (6), we have the following result.Corollary 3.
(i)
Ifr−σ12/2<0 and e12α12−d2−σ22/2<0, then model (6) is extinct exponentially with probability one(ii)
Model (6) is stochastically ultimate bounded(iii)
Ifr−α12/β12−σ12>0 and e12α12−d2−σ22>0, then model (6) is stochastically permanentRemark 2.
From Theorem 4.11 in [18], it follows that if r−α12/β12−3/2σ12>0 and e12α12−d2−3/2σ32>0, then model (6) is stochastically permanent. Obviously, if conditions of Theorem 4.11 in [18] hold, then conditions in Corollary 3 hold. On the contrary, it is not set up. Therefore, Corollary 3 generalizes and improves Theorem 4.11 in [18].Remark 3.
Ifr−α12/β12−σ12>0 and e12α12−d2−σ32>0, then by Theorem 3.3 in [16], model (6) is persistent in mean, but by Corollary 3, model (6) is stochastically permanent.


### 7.2. Three-Species Food-Chain Model

Letα13=0. Then, we can get the following stochastic three-species food chain model:(109)dx1t=x1tr−a1x1t+α12x2tx1t+β12x2tdt+σ1x1tdw1t,dx2t=x2t−d2−a2x2t+e12α12x1tx1t+β12x2t−α23x3tx2t+β23x3tdt+σ2x2tdw2t,dx3t=x3t−d3−a3x3t+e23α23x2tx2t+β23x3tdt+σ3x3tdw3t,with initial value x10,x20,x30=x10,x20,x30∈ℝ+3. Denote(110)κ1′=r−α12β12;κ2′=e12α12−d2−α23β23;κ3′=e23α23−d3;γi′=κi′−σi2,i=1,2,3.For model (109), we have the following results.Corollary 4.
For anyx10,x20,x30∈ℝ+3, let x1t,x2t,x3t be the solution of model (109) with initial value x10,x20,x30.(i)
Ifa2>0, a3>0, and κI′−σi2/2>0 i=1,2,3, then(111)limt⟶∞lnx1tt=0,limt⟶∞lnx2tt=0,limt⟶∞lnx3tt=0a.s.(ii)
Ifa2>0, a3>0, and κi′−σi2/2>0 i=1,2,3, then(112)liminft⟶∞1t∫0tx1sds≥κ1′−σ12/2a1a.s.,liminft⟶∞1t∫0ta2x2s+e12α12β12x2sx1sds≥κ2′−σ222a.s.,liminft⟶∞1t∫0ta3x3s+e23α23β23x3sx2sds≥κ3′−σ322a.s.(iii)
Ifr−σ12/2<0, e12α12−d2−σ22/2<0, and e23α23−d3−σ32/2<0, then model (109) is extinct exponentially with probability one.Corollary 5.
Model (109) is stochastically ultimate bounded. Furthermore, if γi′>0 i=1,2,3, then model (109) is stochastically permanent.Corollary 6.
Ifa2>0, a3>0, κ1′−e12α12β12−σ12>0, κ2′−e23α23β23−σ22>0, and κ3′−σ32>0, then for any x10,x20,x30∈ℝ+3, model (109) has a stationary distribution and the solutions have ergodic property.If we do not consider the intraspecific competition of the predator, i.e.,a2=a3=0 in model (109), then we obtain the following stochastic three-species food chain model:(113)dx1t=x1tr−a2x1t−α12x2tx1t+β12x2tdt+σ1x1tdw1t,dx2t=x2t−d2+e12α12x1tx1t+β12x2t−α23x3tx2t+β23x3tdt+σ2x2tdw2t,dx3t=x3t−d3+e23α23x2tx2t+β23x3tdt+σ3x3tdw3t,with initial value x10,x20,x30=x10,x20,x30∈ℝ+3. From Theorems 4, 7, and 8, for model (113), we have the following result.Corollary 7.
(i)
Ifr−σ12/2<0, e12α12−d2−σ22/2<0, and e23α23−d3−σ32/2<0, then model (113) is extinct exponentially with probability one(ii)
Model (113) is stochastically ultimate bounded(iii)
Ifγi′>0 i=1,2,3, then model (113) is stochastically permanent


### 7.3. Food-Web Model without Intraspecific Competition of Predators

If we do not consider the intraspecific competition of the predator, i.e.,a2=a3=0 in model (8), then we obtain the following stochastic three-species model:(114)dx1t=x1tr−a1x1t−α12x2tx1t+β12x2t−α13x3tx1t+β13x3tdt+σ1x1tdw1t,dx2t=x2t−d2+e12α12x1tx1t+β12x2t−α23x3tx2t+β23x3tdt+σ2x2tdw2t,dx3t=x3t−d3+e13α13x1tx1t+β13x3t−e23α23x2tx2t+β23x3tdt+σ3x3tdw3t,with initial value x10,x20,x30∈ℝ+3. For model (114), we have the following result.Corollary 8.
(i)
Ifr−σ12/2<0, e12α12−d2−σ22/2<0, and e13α13+e23α23−d3−σ32/2<0, then model (114) is extinct exponentially with probability one(ii)
Model (114) is stochastically ultimate bounded(iii)
Ifγi>0 i=1,2,3, then model (114) is stochastically permanent


## 7.1. Two Species Predator-Prey Model

Letα13=α23=0. Then, the first two equations of (8) form the following closed two-population system:(106)dx1t=x1tr−a1x1t+α12x2tx1t+β12x1tdt+σ1x1tdw1t,dx2t=x2t−d2−a2x2t+e12α12x1tx1t+β12x2tdt+σ2x2tdw2t,with initial value x10,x20∈ℝ+2. This is also the stochastic predator-prey model discussed in [19]. From Theorems 3 and 4, we have the following result.Corollary 1.
Letx1t,x2t be solution of model (106) with initial value x10,x20∈ℝ+2.(i)
Ifa2>0, r−α12/β12−σ12/2>0, and e12α12−d2−σ22/2>0, then(107)liminft⟶∞1t∫0tx1sds≥r−α12/β12−σ12/2a1>0a.s.,liminft⟶∞1t∫0ta2x2s+e12α12β12x2sx1sds≥e12α12−d2−σ222>0a.s.(ii)
Ifr−σ12/2<0 and e12α12−d2−σ22/2<0, then (106) is extinct exponentially with probability one.Remark 1.
It is clear that Corollary 1 is consistent with Theorems 7 and 8 in [19]. Moreover, from Theorems 3 and 4, the persistence in mean and extinction conditions of the three-species model (8) are more complicated. Thus, our work can be seen as the extension of [19].
For model (106), similar to the proof of Theorem 6 (denote H=e12x1+x2), we have(108)limsupt⟶∞Ex1t≤K2′e12d2,limsupt⟶∞Ex2t≤K2′d2,where K2′=e12r+d22/4a1. Furthermore, from Theorems 7–9, for model (106), we have the following result.Corollary 2.
(i)
Model (106) is stochastically ultimate bounded(ii)
Ifr−α12/β12−σ12>0 and e12α12−d2−σ22>0, then model (106) is stochastically permanent(iii)
Ifa2>0, r−α12/β12−e12α12β12−σ12>0, and e12α12−d2−σ22>0, then for any x10,x20∈ℝ+2, model (106) has a stationary distribution and the solutions have ergodic propertyIf we do not consider the intraspecific competition of the predator, i.e.,a2=0 in model (106), then model (6) is available. From Theorems 4, 7, and 8, for model (6), we have the following result.Corollary 3.
(i)
Ifr−σ12/2<0 and e12α12−d2−σ22/2<0, then model (6) is extinct exponentially with probability one(ii)
Model (6) is stochastically ultimate bounded(iii)
Ifr−α12/β12−σ12>0 and e12α12−d2−σ22>0, then model (6) is stochastically permanentRemark 2.
From Theorem 4.11 in [18], it follows that if r−α12/β12−3/2σ12>0 and e12α12−d2−3/2σ32>0, then model (6) is stochastically permanent. Obviously, if conditions of Theorem 4.11 in [18] hold, then conditions in Corollary 3 hold. On the contrary, it is not set up. Therefore, Corollary 3 generalizes and improves Theorem 4.11 in [18].Remark 3.
Ifr−α12/β12−σ12>0 and e12α12−d2−σ32>0, then by Theorem 3.3 in [16], model (6) is persistent in mean, but by Corollary 3, model (6) is stochastically permanent.


## 7.2. Three-Species Food-Chain Model

Letα13=0. Then, we can get the following stochastic three-species food chain model:(109)dx1t=x1tr−a1x1t+α12x2tx1t+β12x2tdt+σ1x1tdw1t,dx2t=x2t−d2−a2x2t+e12α12x1tx1t+β12x2t−α23x3tx2t+β23x3tdt+σ2x2tdw2t,dx3t=x3t−d3−a3x3t+e23α23x2tx2t+β23x3tdt+σ3x3tdw3t,with initial value x10,x20,x30=x10,x20,x30∈ℝ+3. Denote(110)κ1′=r−α12β12;κ2′=e12α12−d2−α23β23;κ3′=e23α23−d3;γi′=κi′−σi2,i=1,2,3.For model (109), we have the following results.Corollary 4.
For anyx10,x20,x30∈ℝ+3, let x1t,x2t,x3t be the solution of model (109) with initial value x10,x20,x30.(i)
Ifa2>0, a3>0, and κI′−σi2/2>0 i=1,2,3, then(111)limt⟶∞lnx1tt=0,limt⟶∞lnx2tt=0,limt⟶∞lnx3tt=0a.s.(ii)
Ifa2>0, a3>0, and κi′−σi2/2>0 i=1,2,3, then(112)liminft⟶∞1t∫0tx1sds≥κ1′−σ12/2a1a.s.,liminft⟶∞1t∫0ta2x2s+e12α12β12x2sx1sds≥κ2′−σ222a.s.,liminft⟶∞1t∫0ta3x3s+e23α23β23x3sx2sds≥κ3′−σ322a.s.(iii)
Ifr−σ12/2<0, e12α12−d2−σ22/2<0, and e23α23−d3−σ32/2<0, then model (109) is extinct exponentially with probability one.Corollary 5.
Model (109) is stochastically ultimate bounded. Furthermore, if γi′>0 i=1,2,3, then model (109) is stochastically permanent.Corollary 6.
Ifa2>0, a3>0, κ1′−e12α12β12−σ12>0, κ2′−e23α23β23−σ22>0, and κ3′−σ32>0, then for any x10,x20,x30∈ℝ+3, model (109) has a stationary distribution and the solutions have ergodic property.If we do not consider the intraspecific competition of the predator, i.e.,a2=a3=0 in model (109), then we obtain the following stochastic three-species food chain model:(113)dx1t=x1tr−a2x1t−α12x2tx1t+β12x2tdt+σ1x1tdw1t,dx2t=x2t−d2+e12α12x1tx1t+β12x2t−α23x3tx2t+β23x3tdt+σ2x2tdw2t,dx3t=x3t−d3+e23α23x2tx2t+β23x3tdt+σ3x3tdw3t,with initial value x10,x20,x30=x10,x20,x30∈ℝ+3. From Theorems 4, 7, and 8, for model (113), we have the following result.Corollary 7.
(i)
Ifr−σ12/2<0, e12α12−d2−σ22/2<0, and e23α23−d3−σ32/2<0, then model (113) is extinct exponentially with probability one(ii)
Model (113) is stochastically ultimate bounded(iii)
Ifγi′>0 i=1,2,3, then model (113) is stochastically permanent


## 7.3. Food-Web Model without Intraspecific Competition of Predators

If we do not consider the intraspecific competition of the predator, i.e.,a2=a3=0 in model (8), then we obtain the following stochastic three-species model:(114)dx1t=x1tr−a1x1t−α12x2tx1t+β12x2t−α13x3tx1t+β13x3tdt+σ1x1tdw1t,dx2t=x2t−d2+e12α12x1tx1t+β12x2t−α23x3tx2t+β23x3tdt+σ2x2tdw2t,dx3t=x3t−d3+e13α13x1tx1t+β13x3t−e23α23x2tx2t+β23x3tdt+σ3x3tdw3t,with initial value x10,x20,x30∈ℝ+3. For model (114), we have the following result.Corollary 8.
(i)
Ifr−σ12/2<0, e12α12−d2−σ22/2<0, and e13α13+e23α23−d3−σ32/2<0, then model (114) is extinct exponentially with probability one(ii)
Model (114) is stochastically ultimate bounded(iii)
Ifγi>0 i=1,2,3, then model (114) is stochastically permanent


## 8. Numerical Simulations

In this section, we use the Milstein method (see [35]) to substantiate our main results. The numerical simulations of population dynamics are carried out for the academic tests with the arbitrary values of the vital rates and other parameters which do not correspond to some specific biological populations and exhibit only the theoretical properties of numerical solutions of the considered model. To illustrate the theoretical results, we take the parameter values as following with different noise intensities:(115)r=0.62,d2=0.0005,d3=0.0006,a1=0.0005,a2=0.0002,a3=0.0003,α12=0.2,α13=0.15,α23=0.08,e12=0.8,e13=0.5,e23=0.6,β12=β13=β23=1,x10=800,x20=300,x30=200.In Figure1, we choose σi=0 i=1,2,3 and get the solutions of the corresponding deterministic model.Figure 1
The trajectories of model (8) with σ1=σ2=σ3=0.(i)
Assume thatσ12=1.4, σ22=0.4, and σ32=0.3. By a simple computation, r−σ12/2=−0.08<0, e12α12−d2−σ22/2=−0.0405<0, and e13α13+e23α23−d3−σ32/2=−0.0276<0. Thus, the condition of Theorem 4 holds. From Theorem 4, model (8) will become extinct with probability one. As can be seen from Figure 2, all the population becomes extinct.
(ii)
Assume thatσ12=0.44, σ22=0.08, and σ32=0.08. Then, κ1=r−α12/β12−α13/β13=0.27, κ2=e12α12−d2−α23/β23=0.0795, and κ3=e13α13+e23α23−d3=0.1224. Thus, κ1−σ12/2=0.05>0, κ2−σ22/2=0.0395>0, and κ3−σ32/2=0.0824>0. That is, the conditions of Theorem 3 hold. In view of Theorem 3, model (8) is persistent in mean. As can be seen from Figure 3, all the populations are permanent in mean. This is consistent to Theorem 3.
(iii)
Assume thatσ12=0.02, σ22=0.01, σ32=0.01. From (ii), it follows that κ1=0.27, κ2=0.0795 and κ3=0.1224. Thus, κ1−σ12=0.25>0, κ2−σ22=0.0695>0 and κ3−σ32=0.1124>0. Hence, the conditions of Theorem 8 hold. In view of Theorem 8, model (8) is stochastically permanent. From Figure 4 that all the populations are stochastically permanent. This is consistent to Theorem 8.
(iv)
Assume thatσ12=0.02, σ22=0.01, and σ32=0.01. From (ii), it follows that κ1=0.27, κ2=0.0795, and κ3=0.1224. Furthermore, κ1−e12α12β12−e13α13β13−σ12=0.015>0, κ2−e23α23β23−σ22=0.0215>0, and κ3−σ32=0.1244>0. Thus, the conditions of Theorem 9 hold. Therefore, model (8) has a stationary distribution according to Theorem 9 (see Figures 5 and 6).Figure 2
The trajectories of stochastic model (8) with σ12=1.4, σ22=0.4, and σ32=0.3.Figure 3
The trajectories of stochastic model (8) with σ1=0.44, σ2=0.08, and σ3=0.08.Figure 4
The trajectories of stochastic model (8) with σ1=0.02, σ2=0.01, and σ3=0.01.Figure 5
The density for each population at timet=10000 with σ1=0.02, σ2=0.01, and σ3=0.01.Figure 6
The density for each population at timet=50000 with σ1=0.02, σ2=0.01, and σ3=0.01.


## 9. Conclusions and Discussions

This paper is concerned with a stochastic three-species predator-prey food web model with omnivory and ratio-dependent functional response. First, by the comparison theorem of stochastic differential equations, we prove the existence and uniqueness of global positive solution of the model. Next, we investigate an important asymptotic property of the solution, which is crucial to the study of the dynamic behavior of the model. Then, under some conditions, we conclude that the model is persistent in mean and extinct. Moreover, we discuss the stochastic persistence of the model. Furthermore, by constructing a suitable Lyapunov function, we establish sufficient conditions for the existence of an ergodic stationary distribution to the model. Then, we present the application of the main results in some special models. Finally, some numerical simulations are introduced to support the main results.In Section4, we prove that there are two typical phenomena arising in accordance with the relative values of the parameters of the model. In Theorem 3, we give the conditions on the parameters that informally can be stated by saying that the noise intensities σi2 (i = 1, 2, 3) are small compared to the other parameters, such that the species in model (8) are persistent in mean. From Theorem 4, it follows that in the case that the noise intensities σi2 (i = 1, 2, 3) are large with respect to the other parameters, then the solution of model (8) tends to extinction almost surely.Later, in Section5, we discuss on the stochastic permanence of the solution. This concept, which can be paraphrased by saying that the species in model (8) will survive forever, is one of the most important and interesting topics in the analysis of the model. From Theorem 8, if the noise intensities σi2 (i = 1, 2, 3) are small compared to the other parameters, such that κi−σi2>0 (i=1,2,3), then model (8) is stochastically permanent.Moreover, in Section6, by constructing a suitable Lyapunov function, we show that there is an ergodic stationary distribution for the solution of model (8). In Theorem 9, we give the conditions on the parameters that can be stated by saying that the intensity σi2 of white noise w˙it is sufficiently small, such that the solution model (8) has an ergodic stationary distribution.The results in this paper generalize and improve the previous related results. From Remark 1, we know that our work can be seen as the extension of [19]. From Remark 2, we know that Theorem 8 generalizes and improves Theorem 4.11 in [18].



---

*Source: 4876165-2019-11-13.xml*