# Dynamics and Optimal Harvesting Control for a Stochastic One-Predator-Two-Prey Time Delay System with Jumps

Authors: Tingting Ma; Xinzhu Meng; Zhengbo Chang
Journal: Complexity
 (2019)

---

## Abstract

We consider a stochastic one-predator-two-prey harvesting model with time delays and Lévy jumps in this paper. Using the comparison theorem of stochastic differential equations and asymptotic approaches, sufficient conditions for persistence in mean and extinction of three species are derived. By analyzing the asymptotic invariant distribution, we study the variation of the persistent level of a population. Then we obtain the conditions of global attractivity and stability in distribution. Furthermore, making use of Hessian matrix method and optimal harvesting theory of differential equations, the explicit forms of optimal harvesting effort and maximum expectation of sustainable yield are obtained. Some numerical simulations are given to illustrate the theoretical results.


---

## Body

## 1. Introduction

Many researchers are widely focused on the complex dynamics of biological systems such as delay population systems [1–3], stochastic population systems [4–12], and impulsive population systems [13–15]. Recently, many scholars have investigated two-species models and studied the extinction and persistence [16–18]. However, in the real world, it is a common phenomenon that one predator catches two or more kinds of preys [19]. Consequently, models with three or more species which can explain the dynamical behaviors of the population accurately are investigated [20, 21]. On the other hand, it is necessary and important to consider time delay caused by the competition and predation of species. Delayed differential equations can exhibit much more complex dynamics than differential equations without delay, and stable equilibrium can become unstable with the effects of a time delay. Therefore, many researchers have studied the Lotka-Volterra time delay models with two competitive preys and one predator [22, 23]. Notice that the composite population systems with stochastic effects and time delays present some complex dynamics; thus this causes widespread researchers concern [24–30].In order to rationalize the model, some coefficients should be modified into the existing models. The one-predator-two-prey model with time delays is described by(1)dX1dt=X1ta1-a11X1t-a12X2t-τ12-a13X3t-τ13,dX2dt=X2ta2-a21X1t-τ21-a22X2t-a23X3t-τ23,dX3dt=X3t-a3+a31e-a3τ31X1t-τ31+a32e-a3τ32X2t-τ32-a33X3t,with initial value(2)Xiθ=φiθ,φiθ>0,φiθ∈C-τ0,0,where Xi(t) denotes the size of the i-th species of the prey at time t (i=1, 2), X3(t) denotes the size of the predator at time t, τij≥0 is the time delay, i,j=1,2,3, θ∈[-τ0,0], τ0=max⁡τij. ai>0 (i=1, 2, 3), and ai(i=1, 2) stands for the growth rates of two preys, and a3 stands for the death rate of the predator, respectively. aii>0 is the intraspecific competition rate of species i, i=1,2,3. a31e-a3τ31>0 and a32e-a3τ32>0 are the efficiency of food conversion. a12 and a21 are the interspecific competition rates between 1 and 2, and a13 and a23 are the capture rates.Many systems may suffer environment perturbation, the growth rateai is affected by the white noise [31, 32], B=(B(t),t≥0) is a real-valued Brownian motion defined on a complete probability space (Ω,ϝ,P) with B(0)=0. (3)a^i=ai+σidBit,i=1,2,3.σi>0 (i=1,2,3) denotes the coefficients of the environmental stochastic effects on the preys and the predator populations, respectively.Generally, the dynamical behavior of the species may suffer from the sudden environmental change significantly. Nevertheless, white noise cannot explain huge, occasionally catastrophic disturbances. Therefore, applying the discontinuous stochastic process as Lévy jump to model the abrupt nature phenomenon in ecosystem is necessary [33–35]. In many cases, populations suffer sudden distribution, and this causes changes in the productivity of marine and freshwater species. How to use the discontinuous stochastic process to study these abrupt nature phenomena has been an interesting topic.Recently, some scholars have applied Lévy jump into their models and showed that Lévy jump could describe sudden random environmental perturbations. According to the Lévy decomposition theorem [36], we get that (4)N~t,dv=Nt,dv-λdvtN~(t,dv) is a compensated Poisson process and N is a Poisson counting measure with characteristic measure λ on a measurable subset Y of (0,+∞) with λ(Y)<∞. The distribution of Lévy jump Li(t) can be completely parameterized by (ai,σi,λ) and satisfies the property of infinite divisibility. It is characterized by its characteristic function ΦLi(t); we can get a detailed explanation from the following Lévy-Khintchine formula [37]. There are many other papers about stochastic models with Lévy jump; the readers could refer to [38, 39] and references therein. Considering the inevitable situations in the real world, we assume that the intrinsic growth rates a1 and a2 and the death rate a3 of the model are perturbed by the Lévy jump to signify the sudden climate change, so we introduce the Lévy jump into the underlying stochastic model (1). Taking the economic factors into account, reasonable natural resources management can increase sustainable production and profits. Therefore, harvesting models have been already used to exploited the optimal harvesting policies of renewable resources [40–43]. We only consider to harvest preys x1 and x2 and h1>0 and h2>0 are the harvesting effort rates of x1 and x2, respectively [42, 43]. From system (1), we can obtain the following stochastic harvesting model with Lévy jump: (5)dx1t=x1ta1-h1-a11x1t-a12x2t-τ12-a13x3t-τ13dt+σ1x1tdB1t+x1t∫Yγ1vN~dt,dv,dx2t=x2ta2-h2-a21x1t-τ21-a22x2t-a23x3t-τ23dt+σ2x2tdB2t+x2t∫Yγ2vN~dt,dv,dx3t=x3t-a3+a31e-a3τ31x1t-τ31+a32e-a3τ32x2t-τ32-a33x3tdt+σ3x3tdB3t+x3t∫Yγ3vN~dt,dv.This paper is organized as follows. In Section1, we formulate the model. We show the solution of model (5) is global and positive. We give the main theorems for persistence in mean and extinction under the model (5) and its proof in Section 3. In Section 4, we prove the global attractivity and stability in distribution. Our main aim of this paper is to investigate the optimal strategy of the proposed model; we give the conclusions in Section 5. Finally, we carry out numerical simulations and some figures to support the main conclusions in Section 6.


## 2. Global Positive Solution

In order to explore the dynamical behaviors of ecological population, we first study the positivity of the solutions of system (5).Lemma 1.
There exists an positive constantMmax; we have (6)ln⁡1+γiν≤Mmax,i=1,2,3.Lemma 2 (see [44]).
We assume that, for the eachm>0, there exists an Lm satisfying
D i x , ν - D i y , ν 2 K 0 d ν ≤ L m x - y 2, and Di(x,ν)=γi(ν)xi, i=1,2,3 with x∧y≤m.Lemma 3.
For any given initial valueφ(θ)=(φ1(θ),φ2(θ),φ3(θ))∈C([-τ0,0],R+3), model (5) has a unique global positive solution x(t) on t≥0a.s. Moreover, there is a positive constant K such that (7)limt→∞Exit≤K,i=1,2,3.Proof.
For any given initial value(x1(0),x2(0),x3(0)), there is a unique positive (x1(t),x2(t),x3(t))∈R+3 for t∈[0,τe], where τe is the explosion time. We verify that the positive solution is global, that is τe=∞ as. Let m0≥0 be sufficiently large so x1(t), x2(t), and x3(t) all lie within the interval [1/m0,m0]. For each integer m≥m0, we define the stopping time: (8)τe=inf⁡t∈-ω,τe:xi∉1m0,m0,i=1,2,3.Evidently, τe is strictly increasing when m→∞. Let τ∞=limm→∞τm; thus τ∞≤τm as. If this statement is not true, then there exist pairs of constants T>0, m1≥m0, and 0<ε<1 such that (9)Pτ∞≤T≥ε,m≥m1.Vx1,x2,x3=x1-a-aln⁡x1a+x2-b-bln⁡x2b+x3-c-cln⁡x3c+aa13∫t-τ13tx3rdr+bc23∫t-τ23tx3rdr+aa12∫t-τ12tx2rdr+a21∫t-τ21tx1rdr.Define a C3-function by letting a,b,c be positive constants. UsingItô’s formula, we get (10)dV=LVdt+σ1x1-adB1t+σ2x2-bdB2t+σ3x3-cdB3t+∫Yγ1νx1t--aln⁡1+γ1νx1t+γ2νx2t--bln⁡1+γ2νx2t+γ3νx3t--cln⁡1+γ3νx3tK0dvwhere (11)LV=a1-h1+aa11+bc2a21x1+aa12+a2-h2+ba22x2+aa13+ba23-a3+a31ea3τ31a1+a32ea3τ32a2+a32cx3+-aa1-ah1+12σ12+-ba2-bh2+12σ22+ca3+12σ32-a11x12-a22x22-a33x32+∫Yaγ1ν-aln⁡1+γ1ν+bγ2ν-bln⁡1+γ2ν+cγ3ν-cln⁡1+γ3νN~dt,dv≤K0+3Mmax≤K,where K is a positive constant, (12)Mmax=max⁡∫Yaγ1ν-aln⁡1+γ1ν,∫Ybγ2ν-bln⁡1+γ2ν,∫Ycγ3ν-cln⁡1+γ3ν.Thus, (13)dV≤Kdt+σ1x1-adB1t+σ2x2-bdB2t+σ3x3-cdB3t+∫Yγ1νx1t--aln⁡1+γ1νx1t+γ2νx2t--bln⁡1+γ2νx2t+γ3νx3t--cln⁡1+γ3νx3tN~dt,dv.Taking expectation, yields(14)EVx1τm∧T,x2τm∧T,x3τm∧T≤Vx10,x20,x30+E∫0τm∧TKdt.Set Ωm=τm∧T for m≥m0 and we can obtain P(Ωm)≥ε. For each ω∈Ωm, there are x1(τm,ω), x2(τm,ω), x3(τm,ω) equaling either m or 1/m and yields (15)Vx10,x20,x30+KT≥E1ΩkωVx1τm,ω,x2τm,ω,x3τm,ω≥ε1m-1-ln⁡1m∧m-1-ln⁡m,where 1Ωk(ω) stands for the indicator function of Ωk(ω).
Letm→∞, which implies (16)∞>Vx10,x20,x30+KT=∞is a contradiction. So, we have that τ∞=∞. This completes the proof of Lemma 1.


## 3. Persistence in Mean and Extinction of the Model

For the sake of convenience, we introduce some following notations. LetR+3={a=(a1,a2,a3)∈R3∣ai>0,i=1,2,3} and C([-τ0,0],R+3) stand for all continuous functions from [-τ0,0] to R+3. Additionally, we give some notations (17)ηi=∫Yγiv-ln⁡1+γivλdv,bi=ai-hi-12σi2-ηi,i=1,2,b3=ai+12σi2-η3,Mit=∫0t∫Yln⁡1+γivN~ds,dv,A=a11a12a13a21a22a23-a31e-a3τ31-a32e-a3τ32a33,A1=a1-h1a12a13a2-h2a22a23-a3-a32e-a3τ32a33,A2=a11a1-h1a13a21a2-h2a23-a31e-a3τ31-a3a33,A3=a11a12a1-h1a21a22a2-h2-a31e-a3τ31-a32e-a3τ32-a3,A1~=12σ12a12a1312σ22a22a2312σ32-a32e-a3τ32-a33,A2~=a1112σ12a13a2112σ22a23a3112σ32a33,A3~=a11a1212σ12a21a2212σ22-a31e-a3τ31-a32e-a3τ3212σ32,Λ=a11a1-h112σ12a21a2-h212σ22-a31e-a3τ31-a312σ32.Γ1=a22a1-h1-a12a2-h2,Γ2=a11a2-h2-a21a1-h1,Γ3=a31e-a3τ31a1-h1-a11a3,Γ1~=12a22σ12-12a12σ22,Γ2~=12a11σ22-12a21σ12,Γ3~=12a11σ22+12a31σ12.(18)dN1t=a1-h1-12σ12-a11eN1t-a12eN2t-τ12-a13eN3t-τ13-∫Yγ1-ln⁡1+γ1λdvdt+σ1dB1t+∫Yln⁡1+γ1N~dt,dv,dN2t=a2-h2-12σ22-a21eN1t-τ21-a22eN2t-a23eN3t-τ23-∫Yγ2-ln⁡1+γ2λdvdt+σ2dB2t+∫Yln⁡1+γ2N~dt,dv,dN3t=-a3-12σ32+a31e-a3τ31eN1t-τ31+a32e-a3τ32eN2t-τ32-a33eN3t-∫Yγ3-ln⁡1+γ3λdvdt+σ3dB3t+∫Yln⁡1+γ3N~dt,dv,with initial value Ni(t)=ln⁡φi(θ). It is not difficult to know that the coefficients of model (5) satisfy the local Lipschitz condition; therefore, model (5) has a unique local solution x(t) on [0,τe], where τe is the explosion time. Applying Itô’s formula, we can get the following solution: (19)xt=x1t=eN1t,x2t=eN2t,x3t=eN3tT.It is the unique positive local solution to model (5). Now let us prove τe=+∞. Thus, we introduce the following auxiliary model: (20)dY1t=Y1ta1-h1-a11Y1tdt+σ1Y1tdB1t+Y1t∫Yγ1vN~dt,dv,dY2t=Y2ta2-h2-a22Y2tdt+σ2Y2tdB2t+Y2t∫Yγ2vN~dt,dv,dY3t=Y3t-a3+a31e-a3τ31Y1t-τ3+a32e-a3τ32Y2t-τ32-a33Y3tdt+σ3Y3tdB3t+Y3t∫Yγ3vN~dt,dv,with initial value (21)Yiθ=φiθ,i=1,2,3,θ∈-τ0,0.Obviously (22)xiθ≤φiθ,i=1,2,3.Before starting proving, we state several hypotheses. We assume that (23)12σ12+η1>12σ22+η2,which implies that the persistent ability of species 1 is better than that of species 2.Assumption 4.
A > 0, Ai>0,i=1,2,3 express that all the population could coexist under the condition that the model frees from stochastic noises. A33>0, Γi>0,i=1,2, which imply that the two prey populations could coexist when there is no environmental noises and the predators are absent, where Aij is the complement minor of aij in the determinant A,i,j=1,2,3.Assumption 5.
A ~ i > 0, A23>0, A31<0, A32<0.Assumption 6.
a 11 > a 12 + a 13, a22>a21+a23, a33>a31e-a3τ31+a32e-a3τ32.Lemma 7 (see [26, 33]).
Letψ∈C[Ω×[0,+∞),R+].
(i) If there exist three constantsλ, λ0>0, and T≥0, such that (24)lg⁡ψt≤λt-λ0∫0tψsds+Ft,a.s.for all t≥0, F∈[[0,∞)×Ω,R] and limt→∞F(t)/t=0 as.
Then, for allt≥T, where α and δi are constants, then (25)ψt∗≤λλ0a.s.,ifλ≥0,limt→∞ψt=0a.s.,ifλ≤0.(ii) If there exist three constants λ0, λ, and T≥0, such that (26)lg⁡ψt≥λt-λ0∫0tψsds+Ft,a.s.for all t≥0, F∈[[0,∞)×Ω,R] and limt→∞F(t)/t=0 as.
Then, for allt≥T, then 〈ψ(t)〉∗≥λ/λ0, a.s.Lemma 8.
For arbitrarilyτ≥0, (27)limt→+∞∫t-τtYis=0,a.s.,i=1,2.Proof.
(28) l i m t → + ∞ 1 t ∫ 0 t Y i s d s = b i a i i ,  a . s .,   i f   b i ≥ 0 ,   i = 1,2 , l i m t → + ∞ Y i t = 0 ,  a . s .,   i f   b i < 0 ,   i = 1,2 .Therefore, if bi<0, then (29)limt→+∞1t∫t-τtYisds=limt→+∞1t∫0tYisds-∫0t-τYisds=0.If bi≥0, then (30)limt→+∞1t∫t-τtYisds=limt→+∞1t∫0tYisds-∫0t-τYisds=biaii-biaii=0.Lemma 9.
For model (20), we have the following.
(a) If b1,b2<0, then (31)limt→+∞Yit=0,i=1,2,3,a.s.
(b) If b1≥0, b2<0 and b3-a31e-a3τ31b1/a11>0, then (32)limt→+∞Y1t=b1a11,limt→+∞Y2t=0,limt→+∞Y3t=0,a.s.
(c) If b1≥0, b2<0 and b3-a31e-a3τ31b1/a11<0, then (33)limt→+∞Y1t=b1a11,limt→+∞Y2t=0,limt→+∞Y3t=a31e-a3τ31b1-a11b3a11a33,a.s.
(d) If b1<0, b2≥0 and b3-a32e-a3τ32b2/a22>0, then (34)limt→+∞Y1t=0,limt→+∞Y2t=b2a22,limt→+∞Y3t=0,a.s.
(e) If b1≥0, b2≥0 and b3-c32e-a3τ32b2/a22≥0, then (35)limt→+∞Y1t=0,limt→+∞Y2t=b2a22,limt→+∞Y3t=0,a.s.
(f) If b1≥0, b2≥0 and b3-a32e-a3τ32b2/a22-a31e-a3τ31b1/a11≥0, then (36)limt→+∞Y1t=b1a11,limt→+∞Y2t=b2a22,limt→+∞Y3t=0,a.s.
(g) If b1≥0, b2≥0 and b3-a32e-a3τ32b2/a22-a31e-a3τ31b1/a11<0, then (37)limt→+∞Y1t=b1a11,limt→+∞Y2t=b2a22,limt→+∞Y3t=a22a31e-a3τ31b1+a11a32e-a3τ32b2-a11a22b3a11a22a33,a.s.Proof.
First, let us prove(a). Applying Itô’s formula to model (20), we can get that (38)ln⁡Y1t-ln⁡Y10=b1t-a11∫0tY1sds+σ1B1t+M1t,(39)ln⁡Y2t-ln⁡Y20=b2t-a22∫0tY2sds+σ2B2t+M2t,(40)ln⁡Y3t-ln⁡Y30=-b3t+a31e-a3τ31∫0tY1s-τ31ds+a32e-a3τ32∫0tY2s-τ32ds-a33∫0tY3sds+σ3B3t+M3t.Dividing both sides of (38), (39), and (40) by t, we can obtain that (41)t-1ln⁡Y1tY10=b1-a11Y1t+σ1B1tt+M1tt,(42)t-1ln⁡Y2tY20=b2-a22Y2t+σ2B2tt+M2tt,(43)t-1ln⁡Y3tY30=-b3+a31e-a3τ31Y1t+a32e-a3τ32Y2t-a33Y3t+a31t∫t-τ31tY1sds-∫-τ310Y1sds-a32t∫t-τ32tY2sds-∫-τ320Y2sds+σ3B3tt+M3tt.Note that (44)limt→∞Bitt=0,a.s.,i=1,2,3.Firstly, we prove (a). Since b1<0, b2<0, using Lemma 7, yields (45)limt→+∞Y1t=0,a.s.and (46)limt→+∞Y2t=0,a.s.Substituting (45) and (46) into (40) gives (47)ln⁡Y3t-ln⁡Y30≤-b3t+εt-a33∫0tY3sds+σ3B3t+M3t,where ε is small enough satisfying b3+ε<0. Applying (i) in Lemma 7 gives (48)limt→+∞Y3t=0,a.s.Secondly, we prove (b) and (c). Using Lemma 7, since b1≥0 and b2<0, it is easy to obtain that (49)limt→+∞Y1t=b1a11,a.s.and (50)limt→+∞Y2t=0.a.s.As a consequence, we can study and discuss the following model: (51)dY1t=Y1ta1-a11Y1tdt+σ1Y1tdB1t+Y1t∫Yγ1vN~dt,dv,dY3t=Y3t-a3+a31e-a3τ31Y1t-τ31+a32e-a3τ32Y2t-τ32-a33Y3tdt+σ3Y3tdB3t+Y3t∫Yγ3vN~dt,dv.We know that (52)ln⁡Y3t-ln⁡Y30=-b3t+a31e-a3τ31∫0tY1s-τ31ds-a33∫0tY3sds+σ3B3t+M3t.It is not difficult to obtain the system (53)limt→+∞Y3t=0,ifa31e-a3τ31b1a11<0.limt→+∞t-1∫0tY3t=a31e-a3τ31b1-a11b3a11a33,ifa31e-a3τ31b1a11>0.The proofs of (d) and (e) are similar to these of (b) and (c).
Then, we give the proofs of (f) and (g). Using Lemma7, since b1≥0 and b2≥0, it is easy to obtain that (54)limt→+∞Y1t=b1a11,a.s.and (55)limt→+∞Y2t=b2a22,a.s.Using (54) and (55), then combining (44) leads to (56)limt→+∞t-1ln⁡Y10=0,a.s.and(57)limt→+∞t-1ln⁡Y20=0,a.s.Multiplying (41), (42), and (43) by a31e-a3τ31, a32e-a3τ32, and a11a21, respectively, and adding them, we can derive that (58)a31e-a3τ31a22tln⁡Y1tY0t+a11e-a3τ32a32tln⁡Y2tY0t+a11a22tln⁡Y3tY0t=a22a31e-a3τ31b1+a11a32e-a3τ32b2-a11a22b3-a11a22a33Y3t-a11a22a31t∫t-τ31tY1sds-∫-τ310Y1sds-a11a22a32t∫t-τ32tY2sds-∫-τ320Y2sds+a22a31e-a3τ31σ1B1t+M1tt+a11a32e-a3τ32σ2B2t+M2tt+a11a22σ3B3t+M3tt.An application of (56), (57), (58), Lemma 7, and Lemma 8 gives that (59)limt→+∞Y3t=0,a.s.,ifa31e-a3τ31b1a11+a32e-a3τ32b2a22-b3<0,limt→+∞∫0tY3tdst=a22a31e-a3τ31b1+a11a32e-a3τ32b2-a11a22b3a11a22a33,a.s.,ifa31e-a3τ31b1a11+a32e-a3τ32b2a22-b3>0.This completes the proof.Lemma 10.
The solution of model (5) obeys (60)limsupt→+∞ln⁡xitt=0,a.s.,i=1,2,3.Proof.
From Lemma9, we can obtain that either limt→+∞ln⁡Yi(t)=0 or 1/tlimt→+∞Yi(t)ds=c, c is a constant, i=1,2,3. Since Yi(t)≥xi(t), we only need to prove the following results.
Iflimt→+∞ln⁡Yi(t)=0, then(61)limsupt→+∞ln⁡xitt≤limsupt→+∞ln⁡Yitt≤0,a.s.,i=1,2,3.If 1/tlimt→+∞∫0t‍Yi(t)=a constant, then (62)limt→+∞ln⁡Yitt=0,a.s.,i=1,2,3,and hence(63)limsupt→+∞ln⁡xitt≤limsupt→+∞ln⁡Yitt=0,a.s.Theorem 11.
For system (5), define(64)ρ1=a1+h11/2σ12,ρ2=Γ2Γ~2,ρ3=A3A~3,π2=Γ3Γ~3,π3=A2A~2.From Assumptions 4, 5, and 6, we can obtain the following results.(I)
IfΛ>0, then ρ1>ρ2>ρ3; moreover,
(i)
ifρ1<1, then species i, i=1,2,3, go to extinction a.s.; i.e., (65)limt→+∞xit=0,i=1,2,3,a.s.(ii)
Ifρ1>1>ρ2, then the predator goes to extinction a.s., and one prey is persistent in mean, and another goes to extinction; i.e., (66)limt→+∞xit=0,i=2,3,limt→+∞x1t=b1a11,a.s.(iii)
Ifρ2>1>ρ3, then then species 3 goes to extinction and species 1,2 are persistent in mean; i.e., (67)limt→+∞xit=Γi-Γ~iA33,i=1,2,a.s.(iv)
Ifρ3>1, then three species are persistent in mean; i.e., (68)limt→+∞xit=Ai-A~iA,i=1,2,3,a.s.(II)
IfΛ<0, then ρ1>π2>π3; moreover,
(v)
ifρ1<1, then species i, (i=1,2,3) go to extinction, a.s;(vi)
ifρ1>1>π2, then the predator goes to extinction and one prey is persistent in mean, and another goes to extinction, a.s.;(vii)
ifπ2>1>π3, then the species 2 goes to extinction and species 1,3 are persistent in mean, a.s;(viii)
ifπ3>1, then three species are persistent in mean, a.s.Proof.
Now we will give the proof of (I), and the proof of (II) is parallel to (I). We can get(69)a1-h1σ12/2+η1>a2-h2σ22/2+η2,Γ2=a11a2-h2-a21a1-h1>0,and, as a consequence, (70)Γ2~=a1112σ12+η2-a2112σ12+η1>0,ρ1-ρ2>0;we can also compute that (71)Γ2Γ~2-Γ3Γ~3=a11ΛΓ~2Γ~3>0,Γ3Γ~3-A3A~3=a11ΛΓ~3A~3>0.So ρ1>ρ2>ρ3 is established.
Applying Itô’s formula to model (5) yields(72)ln⁡x1t-ln⁡x10=b1t-a11∫0tx1sds-a12∫0tx2s-τ12ds-a13∫0tx3s-τ13ds+σ1B1t+M1t,(73)ln⁡x2t-ln⁡x20=b2t-a21∫0tx1s-τ21ds-a22∫0tx2sds-a23∫0tx3s-τ23ds+σ2B2t+M2t,(74)ln⁡x3t-ln⁡x30=-b3t+a31e-a3τ31∫0tx1s-τ31ds+a32e-a3τ32∫0tx2s-τ32ds-a33∫0tx3sds+σ3B3t+M3t.Firstly, we prove (i). Since a11,a12, and a13 are positive, we can get (75)t-1ln⁡x1tx10=t-1b1-a11∫0tx1sds-a12∫0tx2s-τ12ds-a13∫0tx3s-τ13ds+σ1B1tt+M1tt,(76)t-1ln⁡x2tx20=t-1b2-a21∫0tx1s-τ21ds-a22∫0tx2sds-a23∫0tx3s-τ23ds+σ1B2tt+M2tt,(77)t-1ln⁡x3tx30=t-1-b3+a31e-a3τ31∫0tx1s-τ31ds+a32e-a3τ32∫0tx2s-τ32ds-a33∫0tx3sds+σ1B3tt+M3tt.Note that b1<0, b2<0. By Lemma 9, we get (78)limt→+∞xit=0,a.s.,i=1,2,3.In other words, all the populations go to extinction a.s.
Now we give the proof of (ii). Multiplying both sides of (75) and (76) by (-a21) and a11, respectively, and then adding these two equations, we can get that (79)a11tln⁡x2tln⁡x20-a21tln⁡x1tln⁡x10=Γ2-Γ~2-A33x2t-A32x3t-a12a21t∫t-τ12tx2sds-∫-τ120x2sds-a13a21t∫t-τ13tx3sds-∫-τ130x3sds+a11a21t∫t-τ21tx2sds-∫-τ210x2sds+∫t-τ23tx3sds-∫-τ230x3sds+a11σ2B2t+M2t-a21σ1B1t+M1tt,From Lemma 3.4 of [26], we obtain that either limt→+∞Y3(s)=0 or limt→+∞1/t∫0t‍Y3(s)ds=c, c is a constant; it follows that(80)limt→+∞1t∫t-τtY3sds=c.In view of (22) and Lemma 8, we can obtain that (81)limt→+∞1t∫t-τtxis=0,a.s.,i=1,2,3.Putting (60) and (81) into (79) leads to (82)a11tln⁡x2tln⁡x20≤Γ2-Γ~2+ε-A33x2t+a11σ2B2t+M2t-a21σ1B1t+M1tt,where ε is small enough such that Γ2-Γ~2+ε<0. By Lemma 7, we can get (83)limt→+∞x2t=0,a.s.Taking advantage of the above identity, we will get the following two-species models: (84)dx1t=x1ta1-a11x1t-a13x3t-τ13dt+σ1x1tdB1t+x1t∫Yγ1vN~dt,dv,dx3t=x3t-a3+a31e-a3τ31x1t-τ31-a33x3tdt+σ3x3tdB3t+x3t∫Yγ3vN~dt,dv.It is easy to get the result with the similar proof of that in [40–42](85)limt→+∞1t∫0tx1sds=b1a11,limt→+∞x3t=0.Thirdly, we prove (iii). Denote m, n as the solution of the following equations: (86)a22m+a32e-a3τ31n=a12,a23m+a33n=a13.Consequently, (87)m=A21A11>0,n=-A31A11>0.According to Lemma 7, for arbitrarily given ε>0, there exists a T1>0, for all t>T1, (88)t-1mln⁡x2tx20+nln⁡x3tx30≤ε.Multiplying (75), (76), and (77) by (-1), m, and n, respectively, and adding them, one can observe that for sufficiently large t such that t>T1, (89)t-1ln⁡x1tx10-t-1mln⁡x2tx20+nln⁡x3tx30=A1-A~1A11-AA11x1t+t-1a12∫t-τ12tx2sds-∫-τ120x2sds+t-1a13∫t-τ13tx3sds-∫-τ130x3sds-t-1a21m∫t-τ21tx1sds-∫-τ210x1sds-t-1a23m∫t-τ23tx3sds-∫-τ230x3sds+t-1a31n∫t-τ31tx1sds-∫-τ310x1sds+t-1a32n∫t-τ32tx2sds-∫-τ320x2sds+t-1mσ1B1t-mσ2B2t-nσ3B3t+t-1-M1t+mM2t+nM3t.Substituting (60) and (81) into (89) yields (90)1tln⁡x1tx10≤A1-A~1A11+2ε-AA11x1t+t-1σ1B1t-mσ2B2t-nσ3B3t+t-1mM1t+nM2t-M3t.For sufficiently large t, by A1/A~1>ρ1>ρ2>1, and choosing ε>0 to be sufficiently small, then we have (91)x1t≤A1-A~1A11,a.s.Similarly, denote p,q as the solution of the following equations: (92)a11p-a31e-a3τ31q=a21,a13p+a33q=a23.Then we have (93)p=A12A22>0,q=-A32A22>0.In the same way, we can choose a T2>0, for arbitrarily given ε>0, such that (94)t-1pln⁡x1tx10+qln⁡x3tx30≤ε.It follows that, for any sufficiently small ε, there exists T3 and T4 such that multiplying (75), (76), and (77) by p, (−1) and q, respectively, and adding them, we can obtain that (95)t-1ln⁡x2tx20-t-1pln⁡x1tx10+qln⁡x3tx30=A2-A~2A22-AA22x2t-pt-1a12∫t-τ12tx2sds-∫-τ120x2sds-pt-1a13∫t-τ13tx3sds-∫-τ130x3sds-qt-1a31∫t-τ31tx1sds-∫-τ310x1ds-qt-1a32∫t-τ32tx2sds-∫-τ320x2sds+t-1a21∫t-τ21tx1sds-∫-τ210x1sds+t-1a23∫t-τ23tx3sds-∫-τ230x3sds+t-1pσ1B1t-σ2B2t+qσ3B3t+t-1pM1t-M2t+qM3t.Applying (60) and (81) into (95) yields that (96)t-1ln⁡x2tx20≤A2-A~2A22+2ε-AA22x2t+t-1pσ1B1t-σ2B2t+qσ3B3t+t-1pM1t-M2t+qM3t.Note that A2/A~2>A3/A~3>1. According to the arbitrariness of ε and Lemma 7, we have (97)x2t∗≤A2-A~2A,a.s.Substituting (81), (91), and (97) into (77), (98)1tln⁡x3tx30≤a33A3-A~3A+ε-a33x3t+σ3B3t+M3tt.For sufficiently large t, then according to A3<A~3, Lemma 7, and the arbitrariness of ε, it is not difficult to show (99)limt→+∞x3t=0,a.s.Consequently, model (5) reduces to the following model: (100)dx1t=x1ta1-a11x1t-a13x3t-τ13dt+σ1x1tdB1t+x1t∫Yγ1vN~dt,dv,dx2t=x2ta2-a21x1t-τ21-a22x2tdt+σ2x2tdB2t+x2t∫Yγ2vN~dt,dv,which has already been investigated in [5]. Then similarly to the proof of Theorem 5.1 in [5], the following identities can be derived: (101)limt→+∞x1t=Γ1-Γ~1A33,limt→+∞x2t=Γ2-Γ~2A33,a.s.Fourthly, we prove (iv). Since A3-A~3>0, for arbitrariness of ε, the application of Lemma 7 to (98) yields (102)x3t∗≤A3-A~3A,a.s.(103)a12x2t≤a12x2t∗+ε≤a12A2-A~2A+ε,t>T3,a13x3t≤a13x3t∗+ε≤a13A3-A~3A+ε,t>T4.Substituting (103) into (75), we can obtain that, for sufficiently large t, (104)t-1ln⁡x1tx10≥b1-ε-a11x1t-a12A2-A~2A-a13A3-A~3A-2ε+σ1B1tt+M1tt=a11A1-A~1A-a11x1t-3ε+σ1B1tt+M1tt.According to the arbitrariness of ε and Lemma 7, we have (105)x1t∗≥A1-A~1A,a.s.Analogously, we can prove that 〈x2(t)〉∗≥A2-A~2/A is established.
Substituting (103) and (105) into (77) yields that when t is large enough, (106)t-1ln⁡x3tx30≥-b3+a33x3t-a31e-a3τ31A1-A~1A-a32e-a3τ32A2-A~2A-3ε+σ3B3tt+M3tt=a33A3-A~3A-a33x3t-3ε+σ3B3tt+M3tt.According to the arbitrariness of ε and Lemma 7, we have (107)x3t∗≥A3-A~3A,a.s.Subsequently, we have (108)limt→+∞xit=Ai-A~iA,a.s.,i=1,2,3.This completes the proof.


## 4. Stability in Distribution

For the convenience, we define the following notations:(109)c11=a11,c12=a12,c13=a13,c21=a21,c22=a22,c23=a23,c31=a31e-a3τ31,c32=a32e-a3τ32,c33=a33.Definition 12.
Model (5) is globally attractive if limt→+∞Exi-x~i=0, a.s., i=1,2,3, x(t)=(x1(t),x2(t),x3(t))T and x~(t)=(x~1(t),x~2(t),x~3(t))T are two arbitrary solutions with initial conditions x(θ)∈U and x~θ∈U, respectively.Lemma 13.
For anyp>1, there exists a constant K=K(P) which makes the solution x(t) of model (5) with any given initial value satisfy the property that (110)limsupt→+∞Exipt≤K,i=1,2,3.Proof.
The proof is rather standard and hence is omitted.From Lemma13, there is a T>0 such that, for t≥T, E[xip(t)]≤2K. Note that E[xi(t)] is continuous; thus there is a constant K1>0 such that E[xip(t)]<K1, when -τ0≤t<T. Denote L=max⁡2K,K1; then we have (111)Exipt≤L=Lp,t≥τ0,p>0,i=1,2,3.Lemma 14.
Ifc11>c21+c31, c22>-c12+c32, c33>-c31+c32, then model (5) will be asymptotically stable in distribution; i.e., when t→+∞, there exists a unique probability measure μ(·) such that the transition probability density p(t,ξ,·) of x(t) converges weakly to μ(·) with any given initial value ξ(t)∈C([-τ,0];R+3) [45].Proof.
Denote(112)LA=c12+c13-12-c13-c21c21+c23-c23-c31-c32c31-c32.Define qi as the factor of i-th diagonal element of LA. Then applying Kirchhoff’s Matrix Tree Theorem [46, 47], we can see that qi>0,i=1,2,3.
Define(113)Vit=qiln⁡xit-ln⁡x~it+∑j≠i,j=13qicij∫t-τijtxjr-x~jrdr,i=1,2,3,Vt=∑i=13Vit.Calculating the right differential d+V(t) yields (114)d+Vt=∑i=13qisgnxit-x~itdln⁡xit-ln⁡x~it+∑i=13∑j≠i,j=13qicijxjt-x~jtdt-∑j≠i,j=13qicijxjt-τij-x~jt-τijdt≤∑i=13-qiciixit-x~itdt+∑j≠i,j=13qicijxjt-x~jtdt=-∑i=13∑j=13qiAxjt-x~jtdt.From Theorem 2.3 in [48], (115)∑i=13∑j=13qieijxjt-x~jt=∑i=13∑j=13qieijxit-yit.Then we can obtain (116)d+Vt≤-∑i=13∑j=13qiAxit-x~itdt=-∑i=13qicii-∑j≠i,j=13cijxit-x~itdt.Namely, (117)EVt≤V0-∫0t∑i=13cii-∑j≠i,j=13aijExir-yirdr.Subsequently, (118)EVt+∫0t∑i=13qicii-∑j≠i,j=13cijExir-x~irdr≤V0<∞.Note that V(t)≥0; then Exis-x~is is integrable on [0,+∞).
In other words,Exit-x~it∈L-1[0,+∞),i=1,2,3. Moreover, from model (5), we have (119)Ex1t=x10+∫0tEa1x1s-c11x12r-c12x1rx2r-τ12-c13x1rx3r-τ13dr,Ex2t=x20+∫0tEa2x2r-c21x2rx1r-τ21-c22x22r-c23x2rx3r-τ23dr,Ex3t=x30+∫0tE-a3x3r+c31x3rx1r-τ31+c32x3rx2r-τ32-c33x32rdr.Thus, (120)dExitt≤aiKi,i=1,2.dEx3tt≤-a3K3,where Ki>0. Therefore, E(xi(t),(i=1,2,3) are uniformly continuous. In other words, E(xi(t)),i=1,2,3 are continuously differentiable functions with respect to t. By Lemma 13 and the conclusion of [45], one can observe that (121)limt→+∞Exit-yit=0,i=1,2,3.Suppose that p(t,ξ,dy) is the transition probability density of the process x(t) and P(t,ξ,M) denotes the probability of event x(t;ξ)∈M with the initial value ξ(θ)∈C[-τ0,0];R+3. By Lemma 7 and Chebyshev’s inequality [49], we can see that the family of p(t,ξ,dy) is tight. So we can get a compact subset K∈R+3 such that p(t,ξ,K)≥1-ε∗ for given ε∗.
LetP(C[-τ,0];R+3) be the probability measures on C[-τ0,0];R+3. For arbitrary two measures P1,P2∈P, we define the following Kantorovich metric: (122)dLP1,P2=supl∈L∫R+3lxP1dx-∫R+3lxP2dx,where(123)L=l:C-τ,0;R+3→R:lx-ly≤x-y,l·≤1.For any l∈L and t,u>0, we get (124)Elxt+u;ξ-Elxt;ξ=EElxt+u;ξ∣Fu-Elxt;ξ=∫R+3Elxt;ξpu,ξ,dϕ-Elxt;ξ≤∫R+3Elxt;ϕ-Elxt;ξpu,ξ,dϕ.From (121), there exists a T>0 such that, for t≥T, we have (125)supl∈LElxt;ϕ-Elxt;ξ≤ϵ∗3.Obviously (126)Elxt+u;ξ-Elxt;ξ≤ϵ∗.By the arbitrariness of l, we have (127)supg∈LElxt+u;ξ-Elxt;ξ≤ϵ∗.Thus, (128)dLpt+u,ξ,·,pt,ξ,·≤ϵ∗,∀t≥T,u>0.Therefore, pt,0,·:t≥0 is Cauchy in P with metric dL.
There exists a uniqueμ(·)∈PC-τ0,0;R+3 such that limt→∞dL(p(t,0,·),μ(·))=0. From (132), we can obtain (129)limt→∞dLpt,ξ,·,pt,0,·=0.Consequently,(130)limt→∞dLpt,ξ,·,μ·≤limt→∞dLpt,ξ,·,pt,0,·+limt→∞dLpt,0,·,μ·=0.This completes the proof of Lemma 14.


## 5. Optimal Harvesting

We give the following extra notions to get the optimal harvesting policy:(131)D=a11a12a21a22,Λ1=λ1,λ2H=DD-1H+I-1Q,where Q=(Q1,Q2)T, Qi=ai-1/2σi2-ηi,i=1,2, and I is the unit matrix.Theorem 15.
Suppose thata11>a21,a22>a12, A>0, and D-1+D-1T is positive definite. (i)
Ifλi≥0 and when hi=λi,i=1,2, we have Θ1>0,Θ2>0; then the optimal harvesting effort is H∗=Λ1=[D(D-1)T+I]-1Q and the maximum of ESY is(132)Y∗=Λ1TD-1Q-Λ1,whereΘi=Ai-A~i,i=1,2.(ii)
Whenhi=λi,i=1,2, there is Θi≤0 or λi<0,i=1,2; then the optimal harvesting policy does not exist.Proof.
LetG={H=(h1,h2)T∈R2∣Θi+1>0,hi>0,i=1,2}. When (iv) of Theorem 11 holds, we find that, for every H∈G, and if the optimal harvesting effort H∗ exists, then it must belong to G.
Firstly we prove (i). It is easy to see thatΛ1∈G, so G is not empty. By (iv) of Theorem 11, for any H∈G, we have (133)limt→+∞t-1∫0tHTxsds=∑i=12hilimt→+∞t-1∫0txisds=HTD-1Q-H.Applying to Lemma 14, model (18) has a unique invariant measure μ(·). According to Corollary 3.4.3 in [46], we can get that μ(·) is strong mixing. At the same time, it is ergodic by Theorem 3.2.6 in [46]. Hence, it then follows from (3.3.2) in [46] that (134)limt→+∞t-1∫0tHTxsds=∫R+2HTxμdx.Let ϱ(x) represent the stationary probability density of model (18); we obtain (135)YH=limt→+∞∑i=12Ehixit=limt→+∞EHTxt=∫R+2HTxϱxdx.Since the invariant measure of model (18) is unique, then, according to the one-to-one correspondence between ϱ(x) and its corresponding invariant measure μ(·), we obtain (136)∫R+2HTxϱxdx=∫R+2HTxμdx.That is to say, (137)YH=HTD-1Q-H.Assume that Λ1=(λ1,λ2)T is the unique stagnation point of the following equation: (138)dYHdH=dHTdHD-1Q-H+ddHQ-HTD-1TH=D-1Q-D-1+D-1TH=0.It holds that Λ1=[D(D-1)T+I]-1Q. We can take use of the following Hessian matrix [42, 43]: (139)ddHTdYHdH=ddHdYHdHTT=ddHQTD-1T-HTD-1+P-1TT=-D-1-D-1Tis negatively defined, so Λ1 is the unique extreme point of Y(H). In other words, if Λ1∈G, i.e., λi≥0 and Θi>0,i=1,2, then the optimal harvesting effort is H∗=Λ1 and Y∗ is the maximum value of ESY. We are now in the position to prove (ii). Suppose that the optimal harvesting effort H-∗=(h-1∗,h-2∗)T exists. So H-∗∈G; i.e., Θihi=h-i∗>0,h-i∗≥0,i=1,2. That is to say, if H-∗ is the optimal harvesting effort, then H-∗ must be the unique solution of (138). However, Λ1=(λ1,λ2)T is also the solution of (138). Hence, λi=h-i∗≥0, and Θihi=λi=Θihi=h-i∗>0. It contradicts with the condition.
This completes the proof of Theorem15.


## 6. Numerical Simulations

In this section, we carry out extensive numerical simulations using MATLAB by choosing the following parameters to check the model. And the phase diagrams are given too. Some figures are as follows:a1=4,a2=2.5,a3=0.7,a11=0.8,a12=0.3,a13=0.45,a21=0.1,a22=0.5,a23=0.5,a31=0.4,a32=0.3,a33=0.5,γ1=1,γ2=1,γ3=1,τ12=0.6,τ13=0.7,τ21=0.5,τ23=1,τ31=0.4,τ32=0.6. From Assumptions 4, 5, and 6, we can obtain the following results.(A) Persistent and Extinction. If Λ>0, then ρ1>ρ2>ρ3.(i)
Ifρ1<1, then species i, i=1,2,3, go to extinction.(ii)
Ifρ1>1>ρ2, then the predator goes to extinction and one prey is persistent in mean, and another goes to extinction.(iii)
Ifρ2>1>ρ3, then the species 3 goes to extinction and species 1,2 are persistent in mean.(iv)
Ifρ3>1, then three species are persistent in mean.In Figure1, Λ=46.3305>0, ρ1=0.3061<1; it shows that all the populations are extinct.Figure 1
The extinction of three species is given in Theorem11. (b) stands for the phase portrait of (a). The density of white noises is taken: σ1=7, σ2=10, σ3=6.
(a) (b)In Figure2, Λ=668.1985>0, ρ1=2.7500>1. We can find that only the prey is persistent in mean and another prey and predator are extinct.Figure 2
One prey persistent in mean and others are extinct as given in Theorem11. (b) stands for the phase portrait of (a). The density of white noises is taken: σ1=0.5, σ2=22, σ3=30.
(a) (b)In Figure3, Λ=114.4603>0, ρ1=2.7500>ρ2=1.4194>1. We can find that only the two preys are persistent in mean and the predator is extinct.Figure 3
The persistence in mean of two preys is given in Theorem11. (b) stands for the phase portrait of (a). The density of white noises is taken: σ1=0.5,σ2=0.9,σ3=12.
(a) (b)In Figure4, Λ=0.2990>0, ρ1=7.0400>ρ2=4.0146>ρ3=2.5178>1, and the noise is small.Figure 4
The persistence in mean of two preys and of three species are given in Theorem11. (b) stands for the phase portrait of (a). The density of white noises is taken: σ1=0.8,σ2=0.65,σ3=0.42.
(a) (b)(B) Optimal Harvesting. Regarding the optimal harvesting effort, when a11=0.8>a21=0.3,a22=0.5>a12=0.1, it is not difficult to estimate that D-1+(D-1)T is positive definite. Note that Λ=[D(D-1)T+I]-1Q; we can observe Λ=(λ1,λ2)T=(0.9778,0.9243)T. Then we can find Γ2>0,Γ3>0. Therefore, by Theorem 15, we can observe h1=λ1=0.9778,h2=λ2=0.9243,Y∗=ΛTD-1(Q-Λ)=3.3388. Thus the optimal harvesting policy exists (see Figure 3).In Figure5, in order to check the conclusion, we choose another two types of harvesting data. It is obvious that the optimal harvesting policy leads to the maximum of expectation of sustainable yield.Figure 5
The optimal harvesting effort and the maximum of ESY [45]. Red line is with h1=0.9778,h2=0.9243, blue line is with h1=0.5,h2=0.4, and green line is with h1=0.2,h2=0.1.


## 7. Conclusion

This paper is about a stochastic three-species population model with time delays and Lévy jumps [49]. We also consider the optimal harvesting of preys [50]. To begin with, we establish the modified model. In Theorem 11, we obtain the sufficient criteria for extinction and persistence in mean of each species. Lévy jump is important to study. When ρ1<1, all species go to extinction. When ρ1>1>ρ2, then the predator x3 goes to extinction and x2 is persistent in mean, x1 goes to extinction. When ρ2>1>ρ3, then x3 goes to extinction and x1, x2 are persistent in mean. When ρ3>1, then three species are persistent in mean. The main purpose of this paper is to study the optimal harvesting. After discussing the stability of distribution, we study the optimal harvesting and obtain the maximum yield of two preys.From the numerical simulations, we list the following biological meanings.( 1 ) We find the noise can cause the variation of species. When the noise is large, it in reality can suppress the increase of population, then it dies out.( 2 ) Time delay and Lévy jump have important effects on the persistence in mean and the harvesting yield.In traditional papers, scholars consider two species or three species without optimal harvesting. We consider a three-species model with Lévy jump and optimal harvesting. Time delay is ineluctable in the ecological environment and is necessary to consider delays. Recently, stochastic models with the telephone noise have been studied by many authors [51, 52]. In the future research, we hope to add more realistic conditions and study more interesting topics, for example, pulse process, Markov Chain, telephone noise, and partial differential system [53–55].



---

*Source: 5342031-2019-03-11.xml*