# Dynamical Analysis and Optimal Harvesting Strategy for a Stochastic Delayed Predator-Prey Competitive System with Lévy Jumps

Authors: Guodong Liu; Kaiyuan Liu; Xinzhu Meng
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

This paper develops a theoretical framework to investigate optimal harvesting control for stochastic delay differential systems. We first propose a novel stochastic two-predator and one-prey competitive system subject to time delays and Lévy jumps. Then we obtain sufficient conditions for persistence in mean and extinction of three species by using the stochastic qualitative analysis method. Finally, the optimal harvesting effort and the maximum of expectation of sustainable yield (ESY) are derived from Hessian matrix method and optimal harvesting theory of delay differential equations. Moreover, some numerical simulations are given to illustrate the theoretical results.


---

## Body

## 1. Introduction

Optimal control problem in the field of biological mathematics has been widely concerned by researchers. Resource exploitation always aims to maximum sustainable yield (MSY) or the profit associated with the maximum economic yield (MEY) [1]. How to obtain MSY or MEY involves the optimal harvesting control problem. Therefore, it is interesting and meaningful to investigate optimal harvesting control strategies for biological population models, especially stochastic population models. A number of researchers have investigated single-species or two-species population models [2–4]. However, the above two classes of models can not describe some natural phenomena completely and it is believed that models with three or more species can explain the dynamical behaviors of the population accurately [5–8]. Predator-prey models are arguably known as the most fundamental building blocks of the biosystems and ecosystems as the biomasses are grown out of their resource, which have been widely investigated [9–17]. As we all know, after a predator catches and feeds on a prey, the number of the predators will not increase at once, which needs the processes of digestion and absorption. Therefore, the time delays were considered in many differential systems [18–21], especially stochastic delay models [22–29].The basic model we consider is based on the delay Lotka-Volterra model with two competitive predators and one prey. We propose the following model by assuming that the random perturbations of intrinsic growth rate are subjected to Gaussian white noise [30]:(1)dxt=xtr-1-a11xt-a12y1t-τ12-a13y2t-τ13dt+σ-1xtdB1t,dy1t=y1tr-2-a21e-d21τ21xt-τ21-a22y1t-a23e-d23τ23y2t-τ23dt+σ-2y1tdB2t,dy2t=y2tr-3-a31e-d31τ31xt-τ31-a32e-d32τ32y1t-τ32-a33y2tdt+σ-3y2tdB3t,with initial value(2)xθ=φ0θ,y1θ=φ1θ,y2θ=φ2θ,θ∈-τ0,0,τ0=maxi,j=1,2,3⁡τij,where x(t) denotes the number of the prey species at time t and yi(t)(i=1, 2) denote the predator species at time t. r¯1>0 stands for the growth rate of prey. r¯2<0 and r¯3<0 stand for the death rates of the two predators, respectively. aii>0 is the interspecific competition rate of species i, i=1,2,3, and a23e-d23τ23>0 and a32e-d32τ32>0 are the competition rates of two predators. a12 and a13 are positive representing capture rates, and a21e-d21τ21 and a31e-d31τ31 are negative representing the growth rates from prey relatively. τij≥0 signifies the time delay, and φi(θ)>0(i=0,1,2) are continuous functions on [-τ0,0]. Bi(t) (i=1,2,3) are mutually independent Brownian motion with Bi(0)=0 defined on a complete probability space (Ω,F,P) with a filtration Ftt≥0 that satisfies that it is right continuous and increasing with F0 that contains all P-null sets.Various harvesting models have been used to investigate the optimal harvesting policies of renewable resources (e.g., Beddington and May [30], Neubert [31]). Recently, many authors explored the optimal harvesting models [32–37]. According to the catch-per-unit-effort (CPUE) hypothesis, we consider that predators y1 and y2 are subject to exploitation with harvesting effort rates h1>0 and h2>0, respectively.However, in some cases, models just perturbed by the Gaussian white noise can not effectively and efficiently describe the circumstance when the species suffer sudden catastrophic disturbance in nature. The sudden environmental change can affect the dynamical behavior of the species significantly. Therefore, it is necessary to use the discontinuous stochastic process (e.g., Lévy jump) to model the abrupt nature phenomenon in ecosystem [38–43].To introduce the Lévy jump into the underlying stochastic model (1), we first give some facts about the Lévy jump [33]. Generally, a Lévy process Li(t) can be decomposed into the sum of a linear drift, a Brownian motion Bi(t), and a superposition of centered Poisson processes with different jump sizes λ(dv) which is the rate of arrival of the Poisson process with jump of size v. According to the Lévy decomposition theorem [44], we know that (3)Lit=a~it+σ~iBit+∫YvN~t,dv,i=1,2,…,n,where a~i∈R, σ~i>0, N~(t,dv)=N(t,dv)-λ(dv)t is a compensated Poisson process, and N is a Poisson measure with characteristic measure λ on a measurable subset Y of (0,+∞) with λ(Y)<∞. The distribution of a Lévy process Li(t) has the property of infinite divisibility and is characterized by its characteristic function ϕLi(t), which is given by the following Lévy-Khintchine formula [45]: (4)ϕLit≔EeιtLi=expa~iιt-12σ~i2t2+∫Yeιtv-1-ιtvIv<1λdv, where IG is the indicator function of set G and ι is the imaginary unit. The distribution of Lévy jump Li(t) can be completely parameterized by (a~i,σ~i,λ).Motivated by the above discussion, we can assume that the intrinsic growth rater¯1 and the death rates r¯2 and r¯3 of model (1) perturbed by the Lévy jump to signify the sudden climate change, r¯i→r¯i+γ¯idLi(t) [42, 46], and then we can obtain the following stochastic model incorporating Lévy jump:(5)dxt=xtr1-a11xt-a12y1t-τ12-a13y2t-τ13dt+σ1xtdB1t+xt∫Yγ1vN~dt,dv,dx2t=y1tr2-h1-a21e-d21τ21xt-τ21-a22y1t-a23e-d23τ23y2t-τ23dt+σ2y1tdB2t+y1t∫Yγ2vN~dt,dv,dx3t=y2tr3-h2-a31e-d31τ31xt-τ31-a32e-d32τ32y1t-τ32-a33y2tdt+σ3y2tdB3t+y2t∫Yγ3vN~dt,dv,here ri=r¯i+γ¯ia~i, σi=σ¯i+γ¯iσ~i, γi(v)=γ¯i(v), i=1,2,3, and with initial value (2).In this paper, we devote our main attention to obtain the optimal harvesting control strategy of system (5). To this end, we firstly investigate the dynamical behavior of the three species including persistence in mean and extinction and asymptotically stable distribution. Then we explore how the time delay, sudden environmental shock expressed by Lévy jump, and other factors affect the optimal harvesting policy and the maximum expectation of sustainable yield.This paper is organized as follows. We discuss the persistence in mean and extinction of the three species in Section2. Based on the conclusion of Section 2, we consider the optimal harvesting policy in Section 3. Finally, we conclude our results by numerical simulations and discussions in the last section.


## 2. Persistence in Mean and Extinction

For convenience, we give some notations which will be used for analyzing our main results. LetR+3=a=a1,a2,a3∈R3∣ai>0,i=1,2,3 and let C([-τ0,0];R+3) stand for all continuous functions from [-τ0,0] to R+3, where R3 is a 3-dimensional Euclidean space and R+3 is the positive cone in R3. Correspondingly, R+=(0,+∞). Denote d1=r1, d2=r2-h1, and d3=r3-h2. To make it more convenient, let c11=a11, c12=a12, c13=a13, c21=a21e-d21τ21, c22=a22, c23=a23e-d23τ23, c31=a31e-d31τ31, c32=a32e-d32τ32, and c33=a33. Additionally, denote α1=(c11,c21,c31)T, α2=(c12,c22,c32)T, α3=(c13,c23,c33)T and β=(d1,d2,d3)T, and γ=1/2(σ12,σ22,σ32)T. Accordingly, we define(6)Φ=det⁡α1,α2,α3,Υ=det⁡α1,β,γ,Φ1=det⁡β,α2,α3,Φ~1=detγ,α2,α3,Φ2=det⁡α1,β,α3,Φ~2=det⁡α1,γ,α3,Φ3=det⁡α1,α2,β,Φ~3=detα1,α2,γ,Δ1=c22d1-c12d2,Δ2=c11d2-c21d1,Δ3=c11d3-c31d1,Δ~1=c22σ122-c12σ222,Δ~2=c11σ222-c21σ122,Δ~3=c11σ322-c31σ222,ft=t-1∫0tfsds,f∗=lim⁡supt→+∞⁡t-1∫0tfsds,f∗=lim⁡inft→+∞⁡t-1∫0tfsds.Φij is the complement minor of cij in the deterministic Φ, i,j=1,2,3. Now we give a fundamental assumption of Lévy jump.Assumption 1.
Assume thatγi(v)>-1, and there exists a constant K such that (7)∫Yln⁡1+γiv2λdv≤K<∞.Assumption1 implies that the intensity of Lévy jump can not be sufficiently large.Lemma 2.
For any given initial valueφ(θ)=(φ0(θ),φ1(θ),φ2(θ))T∈C([-τ0,0];R+3), model (5) has a unique global positive solution X(t)=(x(t),y1(t),y2(t))T on t≥0a.s.Proof.
Firstly, we consider the following stochastic model:(8)dln⁡xt=d1-12σ12-c11xt-c12y1t-τ12-c13y2t-τ13-∫Yγ1-ln⁡1+γ1λdvdt+σ1dB1t+∫Yln⁡1+γ1N~dt,dv,dln⁡y1t=d2-12σ22-c21xt-τ21-c22y1t-c23y2t-τ23-∫Yγ2-ln1+γ2λdvdt+σ2dB2t+∫Yln⁡1+γ2N~dt,dv,dln⁡y2t=d3-12σ32-c31xt-τ31-c32y1t-τ32-c33y2t-∫Yγ3-ln⁡1+γ3λdvdt+σ1dB1t+∫Yln⁡1+γ1N~dt,dvwith initial value (2). It is easy to see that the coefficients of model (8) satisfy the local Lipschitz condition; therefore model (8) has a unique local solution X(t) on [0,τe], where τe is the explosion time. According to Itô’s formula, we can find that(9)Xt=xt,y1t,y2tT is the unique positive local solution to model (5). Now let us prove τe=+∞. Thus, we introduce the following auxiliary model:(10)dut=utd1-c11utdt+σ1utdB1t+ut∫Yγ1vN~dt,dv,dv1t=v1td2-c21ut-τ21-c22v1tdt+σ2v1tdB2t+v1t∫Yγ2vN~dt,dv,dv2t=v2td3-c31ut-τ31-c33v2tdt+σ3v2tdB3t+v2t∫Yγ3vN~dt,dvwith initial value(11)uθ=φ0θ,v1θ=φ1θ,v2θ=φ2θ,θ∈-τ0,0.Taking advantage of the comparison theorem for stochastic equation [47] yields that, for t∈[0,τe],(12)xθ<uθ,yiθ<viθa.s.,i=1,2.According to Theorem 4.2 in Jiang and Shi [48], the explicit solution of model (10) is(13)ut=expb1t+σ1B1t+M1tu-10+c11∫0texp⁡b1s+σ1B1s+M1sds,v1t=exp⁡b2t+c21∫0tut-τ21ds+σ2B2t+M2tv1-10+c22∫0texp⁡b2s-∫0sc21uς-τ21dς+σ2B2s+M2s,v2t=exp⁡b3t+c31∫0tut-τ31ds+σ3B3t+M3tv2-10+c33∫0texp⁡b3s-∫0sc31uς-τ31dς+σ3B3s+M2s,where ηi=∫Y[γi(v)-ln(1+γi(v))]λ(dv), bi=di-(1/2)σi2-ηi, and Mit=∫0t∫Yln⁡(1+γi(v))N~(ds,dv), i=1,2,3 [33]. It is not difficult to see that y1(t), y2(t), and y3(t) are existent on t≥0; thereby τe=+∞. This completes the proof of Lemma 2.Lemma 3 (see [22]).
Suppose thatϕ(t)∈C[Ω×[0,+∞);R+].
(i) If there exist three constantsλ0, λ, and T≥0, such that(14)ln⁡ϕt≤λt-λ0∫0tϕsds+α∑i=13Bit+∑i=13δi∫0t∫Yln⁡1+γivN~ds,dv,a.s.,for all t≥T, where α and δi are constants; then(15)ϕ∗≤λλ0a.s.,ifλ≥0,limt→∞ϕt=0a.s.,ifλ≤0.(ii) If there exist three constants λ0, λ, and T≥0, such that(16)ln⁡ϕt≥λt-λ0∫0tϕsds+α∑i=13Bit+∑i=13δi∫0t∫Yln⁡1+γivN~ds,dv,a.s.,for all t≥T; then ϕ∗≥λ/λ0 a.s.Lemma 4.
For model (10), consider the following:
(a) Ifb1<0, then(17)limt→+∞ut=0,limt→+∞vit=0,a.s.,i=1,2.(b) If b1≥0, b2-c21b1/c11≥0, and b3-c31b1/c11≥0, then(18)limt→+∞⁡ut=b1c11,limt→+∞⁡v1t=b2-c21b1/c11c22,limt→+∞⁡v2t=b3-c31b1/c11c33,a.s.;(c) If b1≥0, b2-c21b1/c11<0, and b3-c31b1/c11<0, then(19)limt→+∞⁡ut=b1c11,limt→+∞⁡v1t=0,limt→+∞⁡v2t=0,a.s..(d) If b1≥0, b2-c21b1/c11≥0, and b3-c31b1/c11<0, then(20)limt→+∞⁡ut=b1c11,limt→+∞⁡v1t=b2-c21b1/c11c22,limt→+∞⁡v2t=0,a.s..(e) If b1≥0, b2-c21b1/c11<0, and b3-c31b1/c11≥0, then(21)limt→+∞ut=b1c11,limt→+∞v1t=0,limt→+∞v2t=b3-c31b1/c11c33,a.s.Proof.
Firstly, we prove(a). Applying Itô’s formula to model (10), we can get that(22)ln⁡ut-ln⁡u0=b1t-c11∫0tusds+σ1B1t+M1t,(23)ln⁡v1t-ln⁡v10=b2t-c21∫0tus-τ21ds-c22∫0tv1sds+σ2B2t+M2t,(24)ln⁡v2t-ln⁡v20=b3t-c31∫0tus-τ31ds-c33∫0tv2sds+σ3B3t+M3t.It follows from (22) that(25)t-1ln⁡utu0≤b1+σ1B1tt+M1tt. By Assumption 1, the quadratic variation of M1(t) is (26)M1t,M1tt=∫0t∫Yln⁡1+γ1v2λdvds≤t∫Yln⁡1+γ1v2λdv<ct, where M1t,M1t(t) is Meyer’s angle bracket process. Utilizing the strong law of large numbers for local martingales gives that limt→∞B1t/t=0 and limt→∞M1(t)/t=0. Since b1<0,(27)limt→+∞ut=0,a.s.Substituting (27) into (23) gives that(28)ln⁡v1t-ln⁡v10≤b2t+εt-c22∫0tv1sds+σ2B2t+M2t,where ε is small enough satisfying b2+ε<0. Applying (i) in Lemma 3 gives(29)limt→+∞⁡v1t=0,a.s. Similarly, we can get that(30)limt→+∞⁡v2t=0,a.s.
Secondly, we prove (b). Sinceb1≥0, applying Lemma 3 to (22) yields that(31)limt→+∞⁡ut=b1c11,a.s.We know that(32)t-1ln⁡utu0=b1-c11ut+σ1B1tt+M1tt,(33)t-1ln⁡v1tv10-t-1c21∫t-τ21tusds-∫-τ210usds=b2-c21ut-c22v1t+σ2B2tt+M2tt,(34)t-1ln⁡v2tv20-t-1c31∫t-τ31tusds-∫-τ310usds=b3-c31ut-c33v2t+σ3B3tt+M3tt.Using (31) in (32) gives that(35)limt→+∞t-1ln⁡u0=0,a.s.Multiplying (32) and (33) by -c21 and c11, respectively, and, adding them, we can derive that(36)t-1c11ln⁡v1tv10-t-1c11c21∫t-τ21tusds-∫-τ210usds-t-1ln⁡utu0=b2c11-c11c22v1t+t-1c11σ2B2t-b1c21-t-1c21σ1B1t-t-1c21M1t+t-1c11M2t.An application of (31) gives that(37)limt→+∞⁡t-1∫t-τ21tusds=limt→+∞⁡t-1∫0tusds-t-1∫0t-τ21usds=0,a.s.Consequently, utilizing (35), (36), (37), and Lemma 3 yields that (38)limt→+∞v1t=b2-c21b1/c11c22,a.s. Similarly, we can get when b3-c31b1/c11≥0, (39)limt→+∞v2t=b3-c31b1/c11c33,a.s.
Thirdly, we prove (c). Ifb2-c21b1/c11<0, according to (35), (36), (37), and Lemma 3, one can obtain that (40)limt→+∞v1t=0,a.s. Similarly, if b3-c31b1/c11<0, we can obtain (41)limt→+∞⁡v2t=0,a.s.
The proofs of (d) and (e) are similar to that of (b) and (c), respectively, and hence are omitted.Lemma 5.
The solution of model (5) obeys(42)lim⁡supt→+∞⁡ln⁡xtt≤0,lim⁡supt→+∞ln⁡yitt≤0,a.s.,i=1,2.Proof.
Sincex(t)≤u(t), yi(t)≤vi(t), i=1,2, we only need to prove(43)lim⁡supt→+∞ln⁡utt≤0,lim⁡supt→+∞ln⁡vitt≤0,a.s.,i=1,2.According to the proof of Lemma 4, with the same method of [22], we can get (43). Therefore, (42) is obtained.Theorem 6.
For system (5), if Φ>0, Φi>0, i=1,2,3, and Υ>0, we set ρ1=2d1/σ12, ρ2=Δ2/Δ~2, and ρ3=Φ3/Φ~3, and then ρ1>ρ2>ρ3; moreover,(i)
ifρ1<1, then prey species and predator species go to extinction a.s.; i.e.,(44)limt→+∞xt=0,limt→+∞yit=0,i=1,2,a.s.;(ii)
ifρ1>1>ρ2, then two predators go to extinction a.s., and prey is persistent in mean; i.e.,(45)limt→+∞xt=b1c11,a.s.;(iii)
ifρ2>1>ρ3, then predator species y2 goes to extinction, while prey species x and predator species y1 are persistent in mean; i.e.,(46)limt→+∞⁡xt=Δ1-Δ~1Φ33,limt→+∞⁡y1t=Δ2-Δ~2Φ33a.s.;(iv)
ifρ3>1, then three species are persistent in mean; i.e.,(47)limt→+∞xt=Φ1-Φ~1Φ,limt→+∞y1t=Φ2-Φ~2Φ,limt→+∞y2t=Φ3-Φ~3Φ,a.s.Proof.
It is easy to proveρ1>ρ2>ρ3. Applying Itô’s formula to model (5) yields that(48)ln⁡xt-ln⁡x0=b1t-c11∫0txsds-c12∫0ty1s-τ12ds-c13∫0ty2s-τ13ds+σ1B1t+M1t,(49)ln⁡y1t-ln⁡y10=b2t-c22∫0ty1sds-c21∫0txs-τ21ds-c23∫0ty2s-τ23ds+σ2B2t+M2t,(50)ln⁡y2t-ln⁡y20=b3t-c33∫0ty2sds-c31∫0txs-τ31ds-c32∫0ty1s-τ32ds+σ3B3t+M3t.
Firstly, we prove (i). Sincec11,c12, and c13 are positive, we can get(51)t-1ln⁡xtx0≤b1+σ1B1tt+M1tt.Note that ρ1=2d1/σ12<1; then b1<0. By Lemma 3, we get(52)limt→+∞xt=0,a.s.Substituting this identity into (49), we can observe that, for sufficiently large t,(53)ln⁡y1t-ln⁡y10≤b2t+εt-c22∫0ty1sds+σ2B2t+M2t,where ε is small enough such that b2+ε<0. By Lemma 3, we can get(54)limt→+∞⁡y1t=0,a.s. Similarly, we have(55)limt→+∞⁡y2t=0,a.s.
Secondly, we prove (ii). Since1>ρ2>Δ3/Δ~3, we know(56)Δ3Δ~3<1,i.e.,c11d3-c31d1c11σ32/2-c31σ12/2<1. Simplifying the above inequality gives that c11b3<c31b1, which means b3-c31b1/c11<0. Furthermore, from ρ2<1, we have b2-c21b1/c11<0. According to (c) in Lemma 3 and (12), one can observe that(57)limt→+∞⁡yit=0,a.s.,i=1,2. Substituting the above identity into (48) and using Lemma 3, we can get(58)limt→+∞⁡xt=b1c11.
Thirdly, we prove (iii). Dividing (48), (49), and (50) by t, we can derive the following equations:(59)t-1ln⁡xtx0=b1-c11xt-c12y1t-c13y2t+t-1c12∫t-τ12ty1sds-∫-τ120y1sds+t-1c13∫t-τ13ty2sds-∫-τ130y2sds+t-1σ1B1t+t-1M1t,(60)t-1ln⁡y1ty10=b2-c22y1t-c21xt-c23y2t+t-1c21∫t-τ21txsds-∫-τ210xsds+t-1c23∫t-τ23ty2sds-∫-τ230y2sds+t-1σ2B2t+t-1M2t,(61)t-1ln⁡y2ty20=b3-c33y2t-c31xt-c32y1t+t-1c31∫t-τ31txsds-∫-τ310xsds+t-1c32∫t-τ32ty1sds-∫-τ320y1sds+t-1σ3B3t+t-1M3t.Denote m, n as the solution of the following equations: (62)c11m+c21n=c31,c12m+c22n=c32.Consequently,(63)m=-Φ13Φ33>0,n=Φ23Φ33>0. By (12), (43), and Lemma 5, we have(64)lim⁡supt→+∞⁡ln⁡xt≤0,lim⁡supt→+∞⁡ln⁡yit≤0,i=1,2.In addition, (65)lim⁡supt→+∞⁡t-1∫t-τ21txsds=lim⁡supt→+∞⁡t-1∫t-τ31txsds=0,lim⁡supt→+∞⁡t-1∫t-τ12ty1sds=lim⁡supt→+∞⁡t-1∫t-τ32ty1sds=0,lim⁡supt→+∞⁡t-1∫t-τ13ty2sds=lim⁡supt→+∞⁡t-1∫t-τ23ty2sds=0.According to Lemma 5, for arbitrarily given ε>0, there exists a T1>0 such that when t>T1(66)t-1mln⁡xtx0+nln⁡y1ty10≤ε.Multiplying (59), (60), and (61) by -m, -n, and 1, respectively, and adding them, one can observe that, for sufficiently large t such that t>T1,(67)t-1ln⁡y2ty20-t-1mln⁡xtx0+ln⁡y1ty10=Φ3-Φ~3Φ33-ΦΦ33y2t-mt-1c12∫t-τ12ty1sds-∫-τ120y1sds-mt-1c13∫t-τ13ty2sds-∫-τ130y2sds-nt-1c21∫t-τ21txsds-∫-τ210xsds-nt-1c23∫t-τ23ty2sds-∫-τ230y2sds+t-1c31∫t-τ31txsds-∫-τ310xsds+t-1c32∫t-τ32ty1sds-∫-τ320y1sds-t-1mσ1B1t+nσ2B2t-σ3B3t-t-1mM1t+nM2t-M3t.Using (43) in (81) yields that(68)t-1ln⁡y2ty20≤Φ3-Φ~3Φ33+2ε-ΦΦ33y2t-t-1mσ1B1t+nσ2B2t-σ3B3t-t-1mM1t+nM2t-M3t.Since ρ3=Φ3/Φ~3<1, we can choose ε>0 to be sufficiently small such that Φ3-Φ~3/Φ33+2ε<0. Making use of the arbitrariness ε and Lemma 3 gives that (69)limt→+∞y2t=0,a.s. Consequently, model (5) reduces to the following model:(70)dxt=xtd1-c11xt-c12y1t-τ12dt+σ1xtdB1t+xt∫Yγ1vN~dt,dv,dy1t=y1td2-c21xt-τ21-c22y1tdt+σ2y1tdB2t+y1t∫Yγ2vN~dt,dv.For system (70), similarly to the proof of Theorem 1 in [49], the following identities can be derived:(71)limt→+∞⁡xt=Δ1-Δ~1Φ33,limt→+∞⁡y1t=Δ2-Δ~2Φ33,a.s.
Fourthly, we prove (iv). By (68), since ρ3>1, we know from the arbitrariness of ε and Lemma 3 that(72)y2∗≤Φ3-Φ~3Φ,a.s.Denote p,q as the solution of the following equations: (73)c11p+c31q=c21,c13p+c33q=c23.Then we have(74)p=-Φ12Φ22>0,q=Φ32Φ22>0. According to Lemma 5, for arbitrarily given ε>0, there exists a T2>0, such that(75)t-1pln⁡xtx0+qln⁡y2ty20≤ε. Multiplying (59), (60), and (61) by -p, 1, and –q, respectively, and, adding them, we can obtain that, for t>T2,(76)t-1ln⁡y1ty10-t-1pln⁡xtx0+qln⁡y2ty20=Φ2-Φ~2Φ22-ΦΦ22y1t-pt-1c12∫t-τ12ty1sds-∫-τ120y1sds-pt-1c13∫t-τ13ty2sds-∫-τ130y2sds-qt-1c31∫t-τ31txsds-∫-τ310xsds-qt-1c32∫t-τ32ty1sds-∫-τ320y1sds+t-1c23∫t-τ23ty2sds-∫-τ230y2sds+t-1c21∫t-τ21txsds-∫-τ210xsds-t-1pσ1B1t-σ2B2t+qσ3B3t-t-1pM1t-M2t+qM3t.Using (43) in (76) yields that (77)t-1ln⁡y1ty10≤Φ2-Φ~2Φ22+2ε-ΦΦ22x2t-t-1pσ1B1t-σ2B2t+qσ3B3t-t-1pM1t-M2t+qM3t. Note that Φ2/Φ~2>Φ3/Φ~3>1. According to the arbitrariness of ε and Lemma 3, we have(78)y1∗≤Φ2-Φ~2Φ,a.s.It follows that, for any sufficiently small ε, there exist T3 and T4 such that(79)c12y1t≤c12y1∗+ε≤c12Φ2-Φ~2Φ+ε,t>T3,c13y2t≤c13y2∗+ε≤c13Φ3-Φ~3Φ+ε,t>T4.Substituting (79) into (59) results in that, for sufficiently large t, (80)t-1ln⁡xtx0≥b1-c11xt-c12Φ2-Φ~2Φ-c13Φ3-Φ~3Φ-3ε+t-1σ1B1t+t-1M1t=c11Φ1-Φ~1Φ-c11xt-3ε+t-1σ1B1t+t-1M1t. According to the arbitrariness of ε and Lemma 3, we have(81)x∗≥Φ1-Φ~1Φ,a.s.By c31<0, c32>0, one can observe that, for every sufficiently small ε, there exist T5 and T6, such that(82)c31xt≤c31x∗+ε≤c31Φ1-Φ~1Φ+ε,t>T5,c32y1t≤c32y1∗+ε≤c32Φ3-Φ~3Φ+ε,t>T6.Using (82) in (61) yields for t large enough (83)t-1ln⁡y2ty20≥b3-c33y2t-c31Φ1-Φ~1Φ-c32Φ2-Φ~2Φ-3ε+t-1σ3B3t+t-1M3t=c33Φ3-Φ~3Φ-c33y2t-3ε+t-1σ3B3t+t-1M3t. According to the arbitrariness of ε and Lemma 3, we have(84)y2∗≥Φ3-Φ~3Φ,a.s.Combining (72) with (84), one can observe that(85)limt→+∞y2t=Φ3-Φ~3Φ,a.s.In the similar way, using (72), (81), and (60) and then combining them with (78) yield that(86)limt→+∞y1t=Φ2-Φ~2Φ,a.s.Subsequently, we can observe that (87)limt→+∞xt=Φ1-Φ~1Φ,a.s.This completes the proof of Theorem 6.


## 3. Optimal Harvesting

Lemma 7.
For anyp>1, there exists a constant K=K(p) which makes the solution X(t) of model (5) satisfy the property that (88)lim⁡supt→+∞⁡Expt≤K,lim⁡supt→+∞⁡Eyipt≤K,i=1,2.Proof.
The proof is rather standard and hence is omitted.From Lemma7, there is a T>0 such that, for t≥T, E[xp(t)]≤2K and E[yip(t)]≤2K. Note that E[x(t)] and E[yi(t)](i=1,2) are continuous; thus there is a constant K1>0 such that E[xp(t)]<K1 and E[yip(t)]<K1 when -τ0≤t<T. Denote L=max2K,K1; then we have (89)Expt≤L=LpandEyipt≤L=Lp,t≥τ0,p>0,i=1,2.Lemma 8.
Ifa11>a12+a13,a22>-a21e-d21τ21+a23e-d23τ23,anda33>-a31e-d31τ31+a32e-d32τ32, then model (5) will be asymptotically stable in distribution; i.e., when t→+∞, there is a unique probability measure μ(·) such that the transition probability density p(t,ξ,·) of X(t) converges weakly to μ(·) with any given initial value ξ(θ)∈C([-τ,0];R+3).Proof.
Since the proof of Lemma8 is rather standard and hence is omitted. The similar proof can be found in [22].We give the following extra notions to get the optimal harvesting policy:(90)P=a22a23e-d23τ23a32e-d32τ32a33,Λ=λ1,λ2T=PP-1T+I-1Q,where Q=(Q1,Q2)T, Qi=ri+1-1/2σi+12-ηi+1, i=1,2, and I is the unit matrix.Theorem 9.
Suppose thata22>a23e-d23τ23,a33>a32e-d32τ32, Φ>0, and P-1+P-1T is positive definite. (i)
Ifλi≥0 and when hi=λi, i=1,2, we have Γ2>0, Γ3>0, and then the optimal harvesting effort is H∗=Λ=[P(P-1)T+I]-1Q and the maximum of ESY is (91)Y∗=ΛTP-1Q-Λ,whereΓi=Φi-Φ~i, i=2,3.(ii)
Whenhi=λi, i=1,2, there is a Γi≤0 or λi<0, and then the optimal harvesting strategy does not exist.Proof.
DefineG=H=h1,h2T∈R2∣Γi+1>0,hi>0,i=1,2. Therefore, the harvesting effort H∈G. If the optimal harvesting effort H∗ exists, it must belong to G.
Firstly, we prove (i). Obviously,G is not empty, since Λ∈G. By (iv) of Theorem 6, for any H∈G, we have(92)limt→+∞⁡t-1∫0tHTXsds=∑i=12hilimt→+∞⁡t-1∫0tyisds=HTP-1Q-H.According to Lemma 8, we obtain that model (5) has a unique invariant measure μ(·) which is strong mixing and ergodic by Corollary 3.4.3 and Theorem 3.2.6 in [50], respectively. Hence, it can be derived from (3.3.2) in [50] that (93)limt→+∞t-1∫0tHTXsds=∫R+2HTXμdX. Let ϱ(X) represent the stationary probability density of model (5). Since the invariant measure of model (5) is unique, and then, by the one-to-one correspondence between ϱ(x) and its corresponding invariant measure μ(·), we obtain(94)YH=limt→+∞⁡∑i=12Ehiyit=limt→+∞EHTXt=∫R+2HTXϱXdX=∫R+2HTXμdX.Combining (92) with (94), we can get (95)YH=HTP-1Q-H. Let Λ=(λ1,λ2)T be the unique solution of the following equation:(96)dYHdH=dHTdHP-1Q-H+ddHQ-HTP-1TH=P-1Q-P-1+P-1TH=0.Hence Λ=[P(P-1)T+I]-1Q. Obviously, the following Hessian matrix (97)ddHTdYHdH=ddHdYHdHTT=ddHQTP-1T-HTP-1+P-1TT=-P-1-P-1Tis negative definite, so Λ is the global maximum point of Y(H). In other words, if Λ∈G, i.e., λi≥0 and Γi+1>0, i=1,2, then the optimal harvesting effort is H∗=Λ and Y∗ is the maximum value of ESY.
Secondly, we prove (ii). Obviously, if there is ani(i=1,2) such that hi<0, the optimal harvesting strategy does not exist. Then we suppose that the optimal harvesting effort H¯∗=(h¯1∗,h¯2∗)T exists. So H¯∗∈G, i.e., Γi+1hi=h¯i∗>0, h¯i∗≥0, i=1,2. That is to say, H¯∗ is the unique solution of (96). On the other hand, Λ=(λ1,λ2)T is also the solution of (96). Hence, λi=h¯i∗≥0, and Γi+1hi=λi=Γi+1hi=h¯i∗>0, i=1,2. It is in contradiction with the condition, which implies that the optimal harvesting strategy does not exist.
This completes the proof of Theorem9.


## 4. Numerical Simulations and Discussions

It is imperative to understand the influence of environmental perturbations on the coexistence and extinction of species. In this paper, we consider a stochastic competitive delay model of two predators and one prey, taking both Gaussian white noise and Lévy jump into account. Compared with [24], we consider a model where the prey species predated by two different predators and also we introduce harvesting efforts into the two predator species, which is novel and more realistic. This relationship can be found in Yangtze river. There are two kinds of fishes, sturgeon and siniperca chuatsi, which feed on some shrimps. The two fishes are harvested on spring and our model reflects this phenomenon. Additionally, our main purpose is not only to investigate the dynamical behavior, but also to obtain the optimal harvesting strategy of model (5). How to deal with time delay and jump are two key points. The work [51] used the graph-theoretic approach to deal with the delay and jump in network. In this paper, we utilize variable substitution to eliminate the delay. Quardratic variation and the strong law of large numbers for local martingales are applied to deal with the jump.Theorem6 describes sufficient conditions for persistence in mean and extinction of three species, which are derived from the comparison theorem of stochastic differential equations and limit superior theory. When it comes to the harvesting of predators, it is essential to consider the optimal harvesting policy and the maximum expectation of sustainable yield (ESY). Theorem 9 gives the optimal harvesting effort and the maximum of ESY by using Hessian Matrix method and optimal harvesting theory of differential equations.The authors of [33] have investigated the influence of competition of two species on the optimal harvesting strategy. In contrast to [33], the dynamical behaviors of a three species model are analyzed more sufficiently. Not only do we consider the competition of the species, we also take the predation into account. Our results show that the capture rates can affect the persistence in mean and extinction of the species. Correspondingly, the capture rates also influence the optimal harvesting strategy and the maximum expectation of sustainable yield. To sum up, our main results are summarized as follows:(I)
Extinction and persistenceDefineρ1=2d1/σ12, ρ2=Δ2/Δ~2, and ρ3=Φ3/Φ~3.(i) Ifρ1<1, then three species go to extinction a.s.; i.e.,(98)limt→+∞⁡xt=0,limt→+∞⁡yit=0,i=1,2,a.s..(ii) Ifρ1>1>ρ2, then two predators go to extinction a.s., and prey is persistent in mean; i.e.,(99)limt→+∞⁡xt=b1c11,a.s..(iii) Ifρ2>1>ρ3, then predator species y2 goes to extinction, while prey species x and predator species y1 are persistent in mean; i.e.,(100)limt→+∞xt=Δ1-Δ~1Φ33,limt→+∞y1t=Δ2-Δ~2Φ33a.s..(iv) Ifρ3>1, then three species are persistent in mean; i.e.,(101)limt→+∞⁡xt=Φ1-Φ~1Φ,limt→+∞⁡y1t=Φ2-Φ~2Φ,limt→+∞⁡y2t=Φ3-Φ~3Φ,a.s.(II)
Optimal harvesting strategyDefineΛ=(λ1,λ2)T=[P(P-1)T+I]-1Q.(i) Ifλi≥0 and when hi=λi, i=1,2, we have Γ2>0,Γ3>0, then the optimal harvesting effort is H∗=Λ=[P(P-1)T+I]-1Q, and the maximum of ESY is (102)Y∗=ΛTP-1Q-Λ,whereΓi=Φi-Φ~i, i=2,3.(ii) Ifλi<0, then the optimal harvesting strategy does not exist.Next, we give some numerical simulations to illustrate the biological significance of the results. We chooser1=1.2, r2=-0.05, r3=-0.005, h1=0.1, h2=0.005, a11=1.6, a12=1.2, a13=0.3, a21e-d21τ21=-0.85, a22=1.9, a23e-d23τ23=0.6, a31e-d31τ31=-0.4, a32e-d32τ32=1, a33=2.1, τ12=3, τ13=7, τ21=1, τ23=5, τ31=4, τ32=10, and γ1=γ2=γ3=1. Additionally, we denote x(θ)=0.5+0.1sin⁡θ, y1(θ)=0.1+0.05sin⁡θ, and y2(θ)=0.05+0.03sin⁡θ. The densities of white noises will be given in Figures 1(a)–1(d).Figure 1
The persistence in mean and extinction of three species are given in Theorem6. The densities of white noises are taken: (a) σ1=0.77, σ2=0.32, and σ3=0.32, (b) σ1=0.77, σ2=0.32, and σ3=5, (c) σ1=0.6, σ2=10, and σ3=8, (d) σ1=15, σ2=8, and σ3=10.
(a) (b) (c) (d)Figure1(a) shows that all three species are persistent in mean when ρ3>1. Figure 1(b) shows that the prey and one predator are persistent in mean while another predator is extinct when ρ2>1>ρ3. We can find that only the prey is persistent in mean and the two predators are extinct when ρ1>1>ρ2; see Figure 1(c). It is obvious that all three species are extinct when ρ1<1; see Figure 1(d).Regarding the optimal harvesting effort, whena22=1.9>a23e-d23τ23=0.6, a33=2.1>a32e-d32τ32=1, it is not difficult to estimate that P-1+(P-1)T is positive definite. Note that Λ=[P(P-1)T+I]-1Q, and we can observe Λ=(λ1,λ2)T=(0.4452,0.3307)T. Then we can find Γ2>0, Γ3>0. The conditions in Theorem 9 hold; therefore, we have h1=λ1=0.4452, h2=λ2=0.3307, and Y∗=ΛTP-1(Q-Λ)=0.31. Thus the optimal harvesting policy exists (see Figure 2).Figure 2
The optimal harvesting effort and the maximum of ESY are given in Theorem9. Red line is with h1=0.4452, h2=0.3307, blue line is with h1=0.35, h2=0.4 and green line is with h1=0.1, h2=0.6.In Figure2, we illustrate not only the optimal harvesting policy but also another two harvesting policies. It is conspicuous that the optimal harvesting policy leads to the maximum of expectation of sustainable yield.Based on the theoretical analysis and numerical simulations, we present the main biological and ecological meanings of our results.( 1 ) The perturbations we considered are not only Gaussian white noises but also Lévy jump. The results reveal that the Lévy jump may significantly affect the optimal harvesting effort and the maximum of ESY.( 2 ) Time delay is imperative in the ecological environment. It has significant relationship with the persistence in mean and optimal harvesting policy of model (5).( 3 ) We have investigated not only environmental disturbance on the species but harvesting effort affected by human and social factors. The results provide theoretical references for some modern fields, such as fishery management. It is beneficial for people to make a rational exploitation and derive maximum profit.As a matter of fact, with the environmental pollution continually becoming worse, it is significant to consider the species in the polluted environment [4, 33]. There is no exaggeration that plenty of interesting topics deserve further investigations, for example, nonautonomous system, Markov process, impulsive effect, and partial differential system [40, 52–56]. We leave these for future investigations.



---

*Source: 2187274-2019-02-12.xml*