# Modeling and Dynamic Analysis in a Hybrid Stochastic Bioeconomic System with Double Time Delays and Lévy Jumps

Authors: Chao Liu; Longfei Yu; Luping Wang
Journal: Complexity
 (2018)

---

## Abstract

A double delayed hybrid stochastic prey-predator bioeconomic system with Lévy jumps is established and analyzed, where commercial harvesting on prey and environmental stochasticity on population dynamics are considered. Two discrete time delays are utilized to represent the maturation delay of prey and gestation delay of predator, respectively. For a deterministic system, positivity of solutions and uniform persistence of system are discussed. Some sufficient conditions associated with double time delays are derived to discuss asymptotic stability of interior equilibrium. For a stochastic system, existence and uniqueness of a global positive solution are studied. By using the invariant measure theory and singular boundary theory of diffusion process, existence of stochastic Hopf bifurcation and stochastic stability are investigated. By constructing appropriate Lyapunov functions, asymptotic dynamic behavior of the proposed hybrid stochastic system with double time delays and Lévy jumps is discussed. Numerical simulations are provided to show consistency with theoretical analysis.


---

## Body

## 1. Introduction

It is well known that a harvest effort has a strong dynamical impact on the prey-predator system [1–5], which plays a significant role in bioeconomics management among various species in a harvested prey-predator system [6–8]. Furthermore, it is more realistic to investigate the coexistence and interaction mechanism of the harvested prey-predator system by introducing time delays into a model system, such as maturation delay and gestation delay of the population [9–12]. In recent years, combined dynamic effects of the harvest effort and time delay on the population dynamics of prey-predator systems have been widely investigated in [13–18] and references therein, where the asymptotic behavior of the model system around equilibrium is analyzed and stability of bifurcated periodic solutions is studied. Yuan et al. [18] incorporated gestation delay τ>0 as a negative feedback of predator population density in the harvested prey-predator system, which is as follows:
(1)ẋ1t=x1tr1−b1x1t−a1x1tx2tk+x1t−qEx1tm1E+m2x1t,ẋ2t=x2tr2−a2x2t−τk+x1t−τ,where x1t and x2t represent the population density of the prey and predator population, respectively. r1 and r2 stand for the birth rate of the prey and predator population, respectively. b1 represents the intracompetition rate for the prey population, a1 is the maximum value of the per capita reduction rate of the prey population due to predation, and a2 has a similar interpretation to that of a1. k measures the extent to which the surrounding environment provides protection to each population. q denotes the catchability coefficient, E is a constant parameter representing the harvest effort on the prey population, and τ>0 represents the gestation delay of the predator population. By analyzing the dynamical behavior of interior equilibrium and properties of bifurcation phenomena, it reveals that sustainable development of the harvested prey-predator system may be guaranteed by adopting an appropriate harvest effort.It should be noted that the dynamical behavior of the commercially harvested bioeconomic system can be precisely predicted by using stochastic mathematical models [19–23], which can provide an additional degree of realistic reflection in the real world compared to its corresponding deterministic counterpart. Many scholars have incorporated stochastic perturbations into deterministic mathematical models to discuss dynamic effects of environmental noises on population dynamics of the harvested bioeconomic system [19–21], which show that persistence and extinction of population are relevant to time delay and stochastic fluctuations. Combined dynamic effects of time delay and Gaussian white noises on population dynamics of the harvested bioeconomic system as well as optimal harvest control problems are studied in [22, 23]. Recently, it is proved that Lévy jumps can efficiently depict sudden and severe environmental perturbations arising in the real world [24, 25], while these phenomena cannot be described better by Brownian motion.Based on the above analysis, some assumptions are proposed as follows.Assumption 1.
In this paper, we will extend the work in [18] by incorporating commercial harvesting on prey into system (1). Et represents the commercial harvesting effort on prey at time t, w represents the harvesting reward coefficients, c represents the cost per unit harvesting effort for the unit weight of prey, and v is the economic interest of commercial harvesting on prey. Based on the economic theory proposed in [26], an algebraic equation is constructed to study the economic interest of commercial harvesting:
(2)Neteconomicrevenue=totalrevenueTR−totalcostTC.
Based on system (1), TR and TC in (2), it is easy to show that TR=wEtx1t and TC=cEt.Assumption 2.
In this paper, the maturity of the prey population is assumed to be mediated by discrete time delayτ1>0. Furthermore, the reproduction of the predator population after predating the prey population is not instantaneous but will be mediated by some time lag required for gestation of the predator population. τ2>0 represents the gestation delay of the predator population. Hence, we will extend the work in [18] by incorporating two different discrete time delays into system (1), and τ1≠τ2.Assumption 3.
In this paper, the population growth of prey and predator populations affected by environmental stochastic fluctuations is assumed to be a stochastic process. Gaussian white noises and Lévy jumps will be incorporated into system (1) to describe stochastic surrounding environmental factors. σjkj,k=1,2 are nonnegative constants, ξ1t and ξ2t denote multiplicative stochastic excitation and external stochastic excitation related to surrounding environment, respectively. ξjtj=1,2 represents independent Gaussian white noise such that Eξjt=0j=1,2. For any time, t1≠t2, Eξjt1ξjt2=δt2−t1j=1,2, and Eξjt1ξkt2=0j,k=1,2,j≠k, where δ denotes the Dirac delta function. xit− represents the left limit of xit,i=1,2 and γiu>−1,i=1,2. N denotes a Poisson counting measure with characteristic measure λ on a measurable subset Y with λY<+∞ and λ is assumed to be a Lévy measure such that N~dt,du=Ndt,du−λdudt and Y denotes a measurable subset of ℝ+. Throughout this paper, ξiti=1,2 and N are assumed to be independent.In this paper, keeping all these Assumptions1–3 in mind, a double delayed hybrid stochastic prey-predator bioeconomic system with Lévy jumps is established as follows:
(3)dx1t=x1t−τ1r1−b1x1t−τ1−a1x1tx2tk+x1t−qEtx1tm1Et+m2x1tdt+σ11x1tξ1t+σ12ξ2tdt+∫Yγ1ux1t−N~dt,du,dx2t=x2tr2−a2x2t−τ2k+x1t−τ2dt+σ21x2tξ1t+σ22ξ2tdt+∫Yγ2ux2t−N~dt,du,0=Etwx1t−c−v,where the initial conditions for system (3) take the following form:
(4)x1θ≥0,x2θ≥0,θ∈−maxτ1,τ2,0,E0≥0.System (3) can be rewritten in the matrix form as follows:
(5)Ξtdx1tdx2t0=F1x1,x2,EF2x1,x2,EF3x1,x2,E=x1t−τ1r1−b1x1t−τ1−a1x1tx2tk+x1t−qEtx1tm1Et+m2x1tdt+σ11x1tξ1t+σ12ξ2tdt+∫Yγ1ux1t−N~dt,dux2tr2−a2x2t−τ2k+x1t−τ2dt+σ21x2tξ1t+σ22ξ2tdt+∫Yγ2ux2t−N~dt,duEtwx1t−c−v.Remark 1.
Since the algebraic equation in (3) includes no differentiated variables, the third row in matrix
(6)Ξt=100010000
has a corresponding zero row.Recently, some hybrid bioeconomic systems with time delay and stochastic fluctuations are established in [27–30] to investigate the combined dynamic effects of stochastic fluctuation and commercial harvesting on population dynamics. These proposed systems in [27–30] are also constructed by several differential equations with stochastic fluctuations and an algebraic equation. Compared with the previously bioeconomic systems proposed in [19–23, 31] and the references therein, such delayed bioeconomic systems [27–30] can not only discuss coexistence and interaction mechanism of a delayed bioeconomic system under stochastic environmental fluctuations but also investigate population dynamics due to variations of economic interest of commercial harvesting. However, biological characteristics among interacting populations are not considered in [27–29]; time delays such as gestation delay and maturation delay for interacting populations in [27–29] are assumed to be the same discrete value, which contradicts to reality in the real world. Asymptotical stability of interior equilibrium and dynamic effects of Lévy jumps on population dynamics are not discussed in [27–29]. It is proved that Lévy jumps can efficiently depict sudden and severe environmental perturbations arising in the real world [24, 25], while these phenomena can not be described better by Brownian motion. Although the dynamic effects of double time delays have been investigated in [30], the combined dynamic effects of multiple time delays and Lévy jumps on population dynamics have not been investigated in [30]. Dynamic effects of multiple time delays on the hybrid bioeconomic prey-predator system are investigated in [15, 32]. However, the dynamic effect of Lévy jumps and asymptotical stability of solutions of the stochastic system are not studied in [15, 32].To the authors’ best knowledge, population dynamics of the hybrid bioeconomic system with double time delays and Lévy jumps have not been investigated. The rest sections of this paper are organized as follows: in the second section, qualitative analysis of the deterministic system is investigated. Positivity of solutions and uniform persistence of the deterministic system are studied. Asymptotic stability of interior equilibrium of the deterministic system is discussed due to variations of double time delays. In the third section, qualitative analysis of the stochastic system is investigated. Existence and uniqueness of global positive solution are studied. In the absence of double time delays and Lévy jumps, the existence of stochastic Hopf bifurcation and stochastic stability is investigated. In the presence of double time delays and Lévy jumps, by constructing appropriate Lyapunov functions, the asymptotic dynamic behavior of the proposed hybrid stochastic system with double time delays and Lévy jumps is discussed. In the fourth section, numerical simulations are provided to support theoretical findings. Finally, this paper ends with a conclusion.


## 2. Qualitative Analysis of Deterministic System

In the absence of stochastic fluctuations, positivity of solutions and uniform persistence of system (3) with initial conditions (4) will be studied in Lemmas 1 and 2.Lemma 1.
Whenσjk=0j,k=1,2 and γiu=0i=1,2, all solutions of system (3) with initial conditions (4) are positive for all t≥0.Proof 1.
Whenσjk=0j,k=1,2 and γiu=0i=1,2, it is easy to show that Fi:ℝ+3+1→ℝ3 is locally Lipschitz and satisfy the condition, Fi>0, where Fii=1,2,3 can be found in (5).Due to lemma in [33] and Theorem A.4 in [34], all solutions of system (3) with initial conditions (4) exist uniquely and each component of the solution remains within the interval 0,U0 for some U0>0. Standard and simple arguments show that any solution of system (3) with initial conditions (4) always exist and stay positive.Lemma 2.
Whenσjk=0j,k=1,2 and γiu=0i=1,2, if v>0 and τ1 are bounded and satisfy the following inequality:
(7)τ1<min12r1lncb1wr1,12r1lnb1ckmr1+q2+cwkm12r14−b1ckmr1+q8cm1r12,then all solutions of system (3) with initial conditions (4) are uniformly persistent.Proof 2.
Whenσjk=0j,k=1,2 and γiu=0i=1,2, by using Taylor series expansion [35], for x1t, x2t, τ1>0, and τ2>0, it is easy to show that
(8)x1t−τ1=x1t−τ1ddtx1t−τ1ẋ1t+⋯,x2t−τ2=x2t−τ2ddtx2t−τ2ẋ2t+⋯.
Hence, it follows that there existsT1>0 such that
(9)x1t−τ1≤x1t,x2t−τ2≤x2t.According to Lemma1, (9), and the first equation of system (3), it follows from simple computations that ẋ1t<r1x1t−τ1≤r1x1t, which derives that
(10)x1t≤x1t−τ1er1τ1holds for t>T2≔T1+τ1.Fort>T2, it follows from (10) and the first equation of system (3) that
(11)ẋ1t<x1tr1−b1x1t−τ1<x1tr1−b1x1te−2r1τ1,which gives that
(12)limsupt→+∞x1t≤r1e2r1τ1b1≔P1.Ifτ1 is bounded, then P1 is bounded for t>T2.Based on Lemma1 and the second equation of system (3), there exists T3>T2 and it follows from simple computations that ẋ2t>r2x2t1−a2x2t/kr2 holds for t>T3. By using standard comparison arguments, it gives that
(13)liminft→+∞x2t≥kr2a2≔Q2.According to biological interpretations, (12) and the first equation of system (3), it derives that there exists T4>T3 and r1x1t−b1e−2r1τ1x12t−a1x1tx2t/k+r1/b1e2r1τ1>0 holds for t>T4. Furthermore, it follows from Lemma 1 that a1x2t<r1k+r1/b1e2r1τ1 and
(14)limsupt→+∞x2t≤r1a1k+r1b1e2r1τ1≔P2.Ifτ1 is bounded, then P2 is bounded for t>T4.By using the first equation of system (3) and (14), it can be obtained that there exists T5>T4 such that
(15)ẋ1t>r1e−r1τ1x1t−b1x12t−a1P2x1tk−qx1tm1holds for t>T5 and
(16)liminft→+∞x1t≥m1kr12e−2r1τ14b1a1m1P2+qk≔Q1.Furthermore, by using simple computations, it follows from (12) and (16) and the third algebraic equation of system (3) that
(17)liminft→+∞Et≥vwP1−c≔Q3,limsupt→+∞Et≤vwQ1−c≔P3.Ifτ1<1/2r1lncb1/wr1, then it is easy to show that Q3>0 and Q3 is bounded.Ifτ1<1/2r1lnb1ckmr1+q2+cwkm12r14−b1ckmr1+q/8cm1r12, then it is easy to show that P3>0 and P3 are bounded.Based on the above analysis, ifv>0 and τ1 are bounded and satisfy the following inequality
(18)τ1<min12r1lncb1wr1,12r1lnb1ckmr1+q2+cwkm12r14−b1ckmr1+q8cm1r12,then all solutions of system (3) with initial conditions (4) are uniformly persistent.Remark 2.
From the practical and physical perspective of viewpoints, positivity of the solution of the proposed hybrid delayed stochastic bioeconomic system refers to each population survival for a long duration under combined dynamic effects of commercial harvesting, time delay, and stochastic fluctuations. Since natural resources for each population survival are relatively limited within a closed environment, permanence of the hybrid delayed stochastic bioeconomic system interprets that there exist positive finite upper bounds and lower bounds for each population density, which may avoid overpopulation and extinction of each interacting populations. Furthermore, in order to maintain sustainable development of commercially harvested population, commercial harvesting should be constrained within a certain range. It practically interprets there exist positive finite upper bounds and lower bounds for commercial harvesting amount on the predator population.When economic interestv>0, σjk=0j,k=1,2, and γiu=0i=1,2, the interior equilibrium of system (3) can be obtained as follows: M∗x1∗,x2∗,E∗=x1∗,r2k/a2−r2,v/wx1∗−c and x1∗ satisfies the following equation:
(19)x1∗4+A1x1∗3+A2x1∗2+A3x1∗+A4=0,where Ai,i=1,2,3,4 are defined as follows:
(20)A1=m2a2−r2wb1k−r1−cb1wm2b1a2−r2,A2=wm2ka1r2−r1a2−r2+a2−r2b1m1v−cm2b1k−r1wm2b1a2−r2,A3=cm2kr1a2−r2−a1r2k+va2−r2q+m1b1k−r1wm2b1a2−r2,A4=a1r2m1v−qa2−r2q+m1r1wm2b1a2−r2.According to Routh-Hurwitz criterion [35], a sufficient condition for (19) has at least one positive root which is A4<0. Furthermore, M∗ exists provided that x2∗=r2k/a2−r2>0 and E∗=v/wx1∗−c>0. Based on the above analysis, if τ1,τ2,v∈H1, then there exists at least one positive root for (19), and H1 is defined as follows:
(21)H1=τ1,τ2,vτ1≥0,τ2≥0,0<v<minqa2−r2q+m1r1a1r2m1,12r1lncb1wr1,12r1lnb1ckmr1+q2+cwkm12r14−b1ckmr1+q8cm1r12.In the following part, some sufficient conditions associated with double time delays are derived to investigate the local asymptotic stability of the system (3) around M∗. By using the third algebraic equation Et=v/wx1t−c, the system (3) can be rewritten as follows:
(22)ẋ1t=x1t−τ1r1−b1x1t−τ1−a1x1tx2tk+x1t−qvx1tm1v+m2x1twx1t−c,ẋ2t=x2tr2−a2x2t−τ2k+x1t−τ2.For mathematical convenience, some transformationsx1t=x1∗ey1t and x2t=x2∗ey2t are made and system (22) is rewritten as follows:
(23)ẏ1t=−b1x1∗ey1t−τ1ey1tey1t−τ1−1+a1x2∗ey1t−τ1ey1t1k+x1∗−ey2tk+x1t+qvey1t−τ1ey1t1m1v+m2x1∗wx1∗−c−1m1v+m2x1twx1t−c,ẏ2t=a2x2∗x1∗ey1t−τ2−1−k+x1∗ey2t−τ2−1k+x1∗k+x1∗ey1t−τ2.By utilizing the following mathematical relations:(24)ey1t−τ1=ey1t−∫t−τ1tey1sdy1dsds,ey1t−τ2=ey1t−∫t−τ2tey1sdy1dsds,ey2t−τ1=ey2t−∫t−τ1tey2sdy2dsds,ey2t−τ2=ey2t−∫t−τ2tey2sdy2dsds,system (23) can be rewritten as follows:
(25)ẏ1t=−b1x1∗ey1t−τ1ey1t−1ey1t−a1kx2∗ey2t−1k+x1tk+x1∗+qvm2x1∗wx1∗−cey1t−1m1v+m2x1∗wx1∗−cm1v+m2x1twx1∗−c+b1x1∗ey1t−τ1−a1x2∗ey1t−qvey1tm1v+m2x1∗wx1∗−c∫t−τ1tey1sdy1dsds,ẏ2t=a2x1∗x2∗ey1t−1k+x1∗k+x1t−τ2−a2x2∗ey2t−1k+x1t−τ2−a2x1∗x2∗k+x1∗k+x1t−τ2∫t−τ2tey1sdy1dsds+a2x2∗k+x1t−τ2∫t−τ2tey2sdy2dsds.Some sufficient conditions for local asymptotical stability of system (3) around the interior equilibrium M∗ can be concluded as follows.Theorem 1.
Whenσjk=0j,k=1,2 and γiu=0i=1,2, if τ1,τ2,v∈H1∩H2, then system (3) is asymptotically stable around interior equilibrium M∗, where H2 is defined in (A.10).Proof 3.
The proof of Theorem1 can be found in Appendix A of this paper.


## 3. Qualitative Analysis of Stochastic System

Ifτ1,τ2,v∈H1∩H2, then it follows from the third algebraic equation of system (3) that Et=v/wx1t−c. Consequently, system (3) is transformed as follows:
(26)dx1t=x1t−τ1r1−b1x1t−τ1−a1x1tx2tk+x1t−qvx1tm1v+m2x1twx1t−cdt+σ11x1tξ1t+σ12ξ2tdt+∫Yγ1ux1t−N~dt,du,dx2t=x2tr2−a2x2t−τ2k+x1t−τ2dt+σ21x2tξ1t+σ22ξ2tdt+∫Yγ2ux2t−N~dt,du.In the following part, existence and uniqueness of the global positive solution of system (26) are discussed in Theorem 2.Theorem 2.
If∫Yγiu−ln1+γiuλdu≤γ~ii=1,2 and γ~i are positive constants and σjk>0j,k=1,2 is sufficiently small, then system (26) has a unique global positive solution for all t>0 and the solution x1t,x2t∈ℝ+2 for t>0 almost surely.Proof 4.
The proof of Theorem2 can be found in Appendix B of this paper.


### 3.1. Case I: System (3) without Double Time Delays and Lévy Jumps

In this subsection, whenτ1=τ2=0 and γiu=0i=1,2, it follows from simple computations that the linearized form of system (26) around x1∗,x2∗ is as follows:
(27)ẋ1t=α11x1t+α12x2t+σ11x1tξ1t+σ12ξ2t,ẋ2t=α21x1t+α22x2t+σ21x2tξ1t+σ22ξ2t,where α11=r1−2b1x1∗−a1kx2∗/k+x1∗2−qm1v2/m1v+w1x1∗−cm2x1∗2, α12=−a1x1∗/k+x1∗, α21=a2x2∗2/k+x1∗2, and α22=r2−a2x2∗/k+x2∗2.By using Khasminskii transformations,x1t=ωcosθ, x2t=ωsinθ, and system (27) are rewritten as follows:
(28)ω̇t=ωα11cos2θ+α11+α21sinθcosθ+α22sin2θ+ωσ11cos2θ+σ21sin2θξ1t+σ12cosθ+σ22sinθξ2t,θ̇t=α21cos2θ−α12sin2θ+α22−α11cosθsinθ+σ21−σ11sinθcosθξ1t+σ22cosθ−σ12sinθωξ2t.According to Khasminskii limit theorem [36, 37], if τ1=τ2=0, γiu=0i=1,2, and σjk>0j,k=1,2 are sufficiently small, then ωt,θt weakly converges to the two-dimensional Markov diffusion process. Based on the stochastic averaging method, the Itô stochastic differential equation can be obtained as follows:
(29)dω=fωdt+β11dWω+β12dWθ,dθ=fθdt+β21dWω+β22dWθ,where Wω and Wθ represent standard and independent Wiener processes,
(30)fω=ωα11−α222+5ωσ112+σ2128−ωσ11σ214+σ112+σ122+σ212+σ2222ω,fθ=−α122,β112=3ω2σ112+σ2128+σ112+σ122+σ212+σ2222+ω2σ11σ214,β12=0,β21=0,β222=σ21−σ1128+σ112+σ122+σ212+σ2222ω2.Furthermore, some parameter transformations are provided as follows:(31)ϵ1=α11+α212,ϵ2=5σ112+σ212−2σ11σ21,ϵ3=σ112x1∗2+σ212x2∗2+σ122+σ2222,ϵ4=3σ112+σ212+2σ11σ21,ϵ5=α11+α214,ϵ6=σ11−σ212.According toβ12=β21=0 and (31), system (29) can be rewritten as follows:
(32)dω=ω28ϵ1+ϵ2+8ϵ38ωdt+8ϵ3+ϵ4ω28dWω+ϵ5ωdWθ,dθ=ϵ5ωdWω+16ϵ3+ϵ6ω24ωdWθ.It follows from the above analysis andβ12=β21=0 that averaging amplitude ωt refers to the one-dimensional Markov diffusion process, which derives that
(33)dω=ω28ϵ1+ϵ2+8ϵ38ωdt+8ϵ3+ϵ4ω28dWω.Theorem 3.
Whenτ1=τ2=0 and γiu=0i=1,2, stochastic stability of system (26) is discussed due to the variation of stability of averaging amplitude ωt in terms of probability.
(i)
ω=0 is unstable, and system (26) is unstable around x1∗,x2∗ in terms of probability.
(ii)
If28ϵ1+ϵ2−ϵ4>0, then ω=+∞ is attractively natural and system (26) is unstable around x1∗,x2∗ in terms of probability; system (26) does not undergo stochastic Hopf bifurcation around x1∗,x2∗.
(iii)
If28ϵ1+ϵ2−ϵ4=0, then ω=+∞ is strictly natural.
(iv)
If28ϵ1+ϵ2−ϵ4<0, then ω=+∞ is exclusively natural and system (26) is unstable around x1∗,x2∗ in terms of probability.Proof 5.
It follows from the formulation ofϵ3 that ϵ3>0 and β11≠0, which derives that ω=0 is a regular boundary of system (29) (reachable) and ω=0 is unstable. Hence, system (26) is unstable around x1∗,x2∗, which is not relevant to the local stability of the deterministic system around x1∗,x2∗. On the other hand, if ω=+∞, then fω=∞, which derives that ω=+∞ is the second singular boundary of system (29). Based on the singular boundary theory [36, 37], it is easy to show that diffusion exponent ρ1, drifting exponent ρ2, and characteristic value ρ3 of boundary ω=+∞ are computed as follows:
(34)ρ1=2,ρ2=1,ρ3=−limω→+∞2fωωρ1−ρ2β112=−limω→+∞28ϵ3+ω28ϵ1+ϵ28ϵ3+ω2ϵ4=−28ϵ1+ϵ2ϵ4.
If28ϵ1+ϵ2−ϵ4>0, then ω=+∞ is attractively natural and system (26) is unstable around x1∗,x2∗; system (26) does not undergo stochastic Hopf bifurcation around x1∗,x2∗.
If28ϵ1+ϵ2−ϵ4=0, then ω=+∞ is strictly natural, which is just the line of stochastic stability demarcation.
If28ϵ1+ϵ2−ϵ4<0, then ω=+∞ is exclusively natural and system (26) is unstable around x1∗,x2∗ in terms of probability.In the following part, existence of stochastic Hopf bifurcation will be investigated asω=+∞ is exclusively natural. By using Itô stochastic equation of amplitude of ωt, Fokker-Planck equation can be obtained as follows:
(35)∂p∂t=−∂∂ω8ϵ1+ϵ2ω2+8ϵ3p8ω+∂2∂ω28ϵ3+ω2ϵ4p16,with the initial value condition
(36)pω,t∣ω0,t0⟶δω−ω0,t⟶t0,where pω,t∣ω0,t0 denotes the transition probability density of the diffusion process ωt.By virtue of (35), it follows from simple computations that invariant measure of ωt is the steady-state probability density psω, which is as follows:
(37)psω=8ϵ42πϵ3ϵ4Γ2−κπ8ϵ3κΓ1/2−κω2ϵ4ω2+8ϵ3κ−2,where κ=8ϵ1+ϵ2/ϵ4 and Γs=∫0+∞ts−1e−tdt.Theorem 4.
Whenτ1=τ2=0 and γiu=0i=1,2, if 28ϵ1+ϵ2−ϵ4<0 and 8ϵ1+ϵ2−ϵ4<0, then there exists a maximum value ω=ω∗ for psω. Occurrence positions and probabilities of stochastic Hopf bifurcation vary due to values of ϵi,i=1,2,3,4.Proof 6.
By using Namachivaya’s theory [36, 37], it is easy to show that characteristic information of a steady state can be revealed based on the invariant measure. When intension of Gaussian white noises tends to zero, the dynamical behavior of the corresponding deterministic system can be approximated by extreme values of psω. Based on Oseledec ergodic’s theory [36, 37], psω is considered to be the time measurement in the neighbourhood of ωt, which derives that ϵ1<0,, ϵ2>0, ϵ3>0, and ϵ4>0.If there exists a maximum value ofpsω around ζ∗>0, then it is easy to show that ζ∗ is asymptotically stable in terms of probability. On the other hand, if there exists a minimum value of psω (zero), then the opposite case also holds. It is assumed that there exists a maximum value as ω=ω~ for psω; further computations show that
(38)dpsωdωω=ω~=0,d2psωdω2ω=ω~<0,and solutions of the above equation are computed as follows: ω=0 or ω=ω~=−8ϵ3/8ϵ1+ϵ2−ϵ4 as 8ϵ1+ϵ2−ϵ4<0. Furthermore, it can be computed that
(39)d2psωdω2ω=0=27+38ϵ1+ϵ2−4ϵ4ϵ4−1ϵ32+8ϵ1+ϵ2−4ϵ4ϵ4−1>0,d2psωdω2ω=ω~=−8ϵ1+ϵ2−ϵ438ϵ3−8ϵ3ϵ4/8ϵ1+ϵ2−ϵ48ϵ1+ϵ2/ϵ416ϵ328ϵ1+ϵ2−2ϵ43<0.It is easy to show thatpsω=0 when ω=0. Based on the singular boundary theory [36, 37], it shows that system (33) is unstable around ω=0 in terms of probability. Hence, when τ1=τ2=0 and γiu=0i=1,2, system (26) undergoes stochastic Hopf bifurcation around ω=ω~, which derives that
(40)x12t+x22t=−8ϵ38ϵ1+ϵ2−ϵ4,whenω=ω~.Remark 3.
From practical and physical perspective of viewpoints, if the proposed hybrid delayed stochastic bioeconomic system is unstable around the interior equilibrium in terms of probability, then it practically implies that the population density of the prey and predator population will not eventually reach an ecological balance state with a high probability, which is not beneficial to sustainable survival of each population under a commercially harvested ecological environment. If stochastic Hopf bifurcation does not occur and boundary is attractively natural, then it practically interprets that either the prey or the predator population density will sharply increase with a high probability within a short duration, which may result in the corresponding shortage of survival space and resource and beyond environment carrying capacity. The ecological balance state will be destroyed in a high probability.


### 3.2. Case II: System (3) with Double Time Delays and Lévy Jumps

It should be noted thatx1∗,x2∗ is not the equilibrium of the stochastic system (26); we will show that solution of system (26) is going around x1∗,x2∗ under certain conditions.Theorem 5.
Forτ1,τ2,v∈H1∩H2, when ∫Yγiu−ln1+γiuλdu≤γ~ii=1,2 and γ~i are positive constants and σjk>0j,k=1,2 is sufficiently small, then
(41)limsupt→∞1tE∫0tx1s−x1∗2ds≤B1,limsupt→∞1tE∫0tx2s−x2∗2ds≤B2,where B1 and B2 are defined as follows:
(42)B1=a1m1kP1x1∗P1+x2∗1+P1+k+x1∗qkP+x1∗+a1m1P2x1∗km1b1k+x1∗+x1∗Q1x1∗+P12Q1+2qvQ12P12+x1∗2+m2x1∗wQ1−cx1∗1+x1∗σ112+2+Q12σ122+2Q121+x1∗γ~12b1m2x1∗Q12wQ1−c,B2=P2k+P1a2P2+x2∗k+P1+P1P2+x1∗x2∗+a2P2kP21+x2∗+x2∗P1+P2a2k+Q1k1+Q2+Q1+x2∗σ212σ222+2x2∗+σ222+4γ~22a2k1+Q2+Q1.Proof 7.
The proof of Theorem5 can be found in Appendix C of this paper.


## 3.1. Case I: System (3) without Double Time Delays and Lévy Jumps

In this subsection, whenτ1=τ2=0 and γiu=0i=1,2, it follows from simple computations that the linearized form of system (26) around x1∗,x2∗ is as follows:
(27)ẋ1t=α11x1t+α12x2t+σ11x1tξ1t+σ12ξ2t,ẋ2t=α21x1t+α22x2t+σ21x2tξ1t+σ22ξ2t,where α11=r1−2b1x1∗−a1kx2∗/k+x1∗2−qm1v2/m1v+w1x1∗−cm2x1∗2, α12=−a1x1∗/k+x1∗, α21=a2x2∗2/k+x1∗2, and α22=r2−a2x2∗/k+x2∗2.By using Khasminskii transformations,x1t=ωcosθ, x2t=ωsinθ, and system (27) are rewritten as follows:
(28)ω̇t=ωα11cos2θ+α11+α21sinθcosθ+α22sin2θ+ωσ11cos2θ+σ21sin2θξ1t+σ12cosθ+σ22sinθξ2t,θ̇t=α21cos2θ−α12sin2θ+α22−α11cosθsinθ+σ21−σ11sinθcosθξ1t+σ22cosθ−σ12sinθωξ2t.According to Khasminskii limit theorem [36, 37], if τ1=τ2=0, γiu=0i=1,2, and σjk>0j,k=1,2 are sufficiently small, then ωt,θt weakly converges to the two-dimensional Markov diffusion process. Based on the stochastic averaging method, the Itô stochastic differential equation can be obtained as follows:
(29)dω=fωdt+β11dWω+β12dWθ,dθ=fθdt+β21dWω+β22dWθ,where Wω and Wθ represent standard and independent Wiener processes,
(30)fω=ωα11−α222+5ωσ112+σ2128−ωσ11σ214+σ112+σ122+σ212+σ2222ω,fθ=−α122,β112=3ω2σ112+σ2128+σ112+σ122+σ212+σ2222+ω2σ11σ214,β12=0,β21=0,β222=σ21−σ1128+σ112+σ122+σ212+σ2222ω2.Furthermore, some parameter transformations are provided as follows:(31)ϵ1=α11+α212,ϵ2=5σ112+σ212−2σ11σ21,ϵ3=σ112x1∗2+σ212x2∗2+σ122+σ2222,ϵ4=3σ112+σ212+2σ11σ21,ϵ5=α11+α214,ϵ6=σ11−σ212.According toβ12=β21=0 and (31), system (29) can be rewritten as follows:
(32)dω=ω28ϵ1+ϵ2+8ϵ38ωdt+8ϵ3+ϵ4ω28dWω+ϵ5ωdWθ,dθ=ϵ5ωdWω+16ϵ3+ϵ6ω24ωdWθ.It follows from the above analysis andβ12=β21=0 that averaging amplitude ωt refers to the one-dimensional Markov diffusion process, which derives that
(33)dω=ω28ϵ1+ϵ2+8ϵ38ωdt+8ϵ3+ϵ4ω28dWω.Theorem 3.
Whenτ1=τ2=0 and γiu=0i=1,2, stochastic stability of system (26) is discussed due to the variation of stability of averaging amplitude ωt in terms of probability.
(i)
ω=0 is unstable, and system (26) is unstable around x1∗,x2∗ in terms of probability.
(ii)
If28ϵ1+ϵ2−ϵ4>0, then ω=+∞ is attractively natural and system (26) is unstable around x1∗,x2∗ in terms of probability; system (26) does not undergo stochastic Hopf bifurcation around x1∗,x2∗.
(iii)
If28ϵ1+ϵ2−ϵ4=0, then ω=+∞ is strictly natural.
(iv)
If28ϵ1+ϵ2−ϵ4<0, then ω=+∞ is exclusively natural and system (26) is unstable around x1∗,x2∗ in terms of probability.Proof 5.
It follows from the formulation ofϵ3 that ϵ3>0 and β11≠0, which derives that ω=0 is a regular boundary of system (29) (reachable) and ω=0 is unstable. Hence, system (26) is unstable around x1∗,x2∗, which is not relevant to the local stability of the deterministic system around x1∗,x2∗. On the other hand, if ω=+∞, then fω=∞, which derives that ω=+∞ is the second singular boundary of system (29). Based on the singular boundary theory [36, 37], it is easy to show that diffusion exponent ρ1, drifting exponent ρ2, and characteristic value ρ3 of boundary ω=+∞ are computed as follows:
(34)ρ1=2,ρ2=1,ρ3=−limω→+∞2fωωρ1−ρ2β112=−limω→+∞28ϵ3+ω28ϵ1+ϵ28ϵ3+ω2ϵ4=−28ϵ1+ϵ2ϵ4.
If28ϵ1+ϵ2−ϵ4>0, then ω=+∞ is attractively natural and system (26) is unstable around x1∗,x2∗; system (26) does not undergo stochastic Hopf bifurcation around x1∗,x2∗.
If28ϵ1+ϵ2−ϵ4=0, then ω=+∞ is strictly natural, which is just the line of stochastic stability demarcation.
If28ϵ1+ϵ2−ϵ4<0, then ω=+∞ is exclusively natural and system (26) is unstable around x1∗,x2∗ in terms of probability.In the following part, existence of stochastic Hopf bifurcation will be investigated asω=+∞ is exclusively natural. By using Itô stochastic equation of amplitude of ωt, Fokker-Planck equation can be obtained as follows:
(35)∂p∂t=−∂∂ω8ϵ1+ϵ2ω2+8ϵ3p8ω+∂2∂ω28ϵ3+ω2ϵ4p16,with the initial value condition
(36)pω,t∣ω0,t0⟶δω−ω0,t⟶t0,where pω,t∣ω0,t0 denotes the transition probability density of the diffusion process ωt.By virtue of (35), it follows from simple computations that invariant measure of ωt is the steady-state probability density psω, which is as follows:
(37)psω=8ϵ42πϵ3ϵ4Γ2−κπ8ϵ3κΓ1/2−κω2ϵ4ω2+8ϵ3κ−2,where κ=8ϵ1+ϵ2/ϵ4 and Γs=∫0+∞ts−1e−tdt.Theorem 4.
Whenτ1=τ2=0 and γiu=0i=1,2, if 28ϵ1+ϵ2−ϵ4<0 and 8ϵ1+ϵ2−ϵ4<0, then there exists a maximum value ω=ω∗ for psω. Occurrence positions and probabilities of stochastic Hopf bifurcation vary due to values of ϵi,i=1,2,3,4.Proof 6.
By using Namachivaya’s theory [36, 37], it is easy to show that characteristic information of a steady state can be revealed based on the invariant measure. When intension of Gaussian white noises tends to zero, the dynamical behavior of the corresponding deterministic system can be approximated by extreme values of psω. Based on Oseledec ergodic’s theory [36, 37], psω is considered to be the time measurement in the neighbourhood of ωt, which derives that ϵ1<0,, ϵ2>0, ϵ3>0, and ϵ4>0.If there exists a maximum value ofpsω around ζ∗>0, then it is easy to show that ζ∗ is asymptotically stable in terms of probability. On the other hand, if there exists a minimum value of psω (zero), then the opposite case also holds. It is assumed that there exists a maximum value as ω=ω~ for psω; further computations show that
(38)dpsωdωω=ω~=0,d2psωdω2ω=ω~<0,and solutions of the above equation are computed as follows: ω=0 or ω=ω~=−8ϵ3/8ϵ1+ϵ2−ϵ4 as 8ϵ1+ϵ2−ϵ4<0. Furthermore, it can be computed that
(39)d2psωdω2ω=0=27+38ϵ1+ϵ2−4ϵ4ϵ4−1ϵ32+8ϵ1+ϵ2−4ϵ4ϵ4−1>0,d2psωdω2ω=ω~=−8ϵ1+ϵ2−ϵ438ϵ3−8ϵ3ϵ4/8ϵ1+ϵ2−ϵ48ϵ1+ϵ2/ϵ416ϵ328ϵ1+ϵ2−2ϵ43<0.It is easy to show thatpsω=0 when ω=0. Based on the singular boundary theory [36, 37], it shows that system (33) is unstable around ω=0 in terms of probability. Hence, when τ1=τ2=0 and γiu=0i=1,2, system (26) undergoes stochastic Hopf bifurcation around ω=ω~, which derives that
(40)x12t+x22t=−8ϵ38ϵ1+ϵ2−ϵ4,whenω=ω~.Remark 3.
From practical and physical perspective of viewpoints, if the proposed hybrid delayed stochastic bioeconomic system is unstable around the interior equilibrium in terms of probability, then it practically implies that the population density of the prey and predator population will not eventually reach an ecological balance state with a high probability, which is not beneficial to sustainable survival of each population under a commercially harvested ecological environment. If stochastic Hopf bifurcation does not occur and boundary is attractively natural, then it practically interprets that either the prey or the predator population density will sharply increase with a high probability within a short duration, which may result in the corresponding shortage of survival space and resource and beyond environment carrying capacity. The ecological balance state will be destroyed in a high probability.


## 3.2. Case II: System (3) with Double Time Delays and Lévy Jumps

It should be noted thatx1∗,x2∗ is not the equilibrium of the stochastic system (26); we will show that solution of system (26) is going around x1∗,x2∗ under certain conditions.Theorem 5.
Forτ1,τ2,v∈H1∩H2, when ∫Yγiu−ln1+γiuλdu≤γ~ii=1,2 and γ~i are positive constants and σjk>0j,k=1,2 is sufficiently small, then
(41)limsupt→∞1tE∫0tx1s−x1∗2ds≤B1,limsupt→∞1tE∫0tx2s−x2∗2ds≤B2,where B1 and B2 are defined as follows:
(42)B1=a1m1kP1x1∗P1+x2∗1+P1+k+x1∗qkP+x1∗+a1m1P2x1∗km1b1k+x1∗+x1∗Q1x1∗+P12Q1+2qvQ12P12+x1∗2+m2x1∗wQ1−cx1∗1+x1∗σ112+2+Q12σ122+2Q121+x1∗γ~12b1m2x1∗Q12wQ1−c,B2=P2k+P1a2P2+x2∗k+P1+P1P2+x1∗x2∗+a2P2kP21+x2∗+x2∗P1+P2a2k+Q1k1+Q2+Q1+x2∗σ212σ222+2x2∗+σ222+4γ~22a2k1+Q2+Q1.Proof 7.
The proof of Theorem5 can be found in Appendix C of this paper.


## 4. Numerical Simulation

Numerical simulations are carried out to show combined dynamic effects of double time delays and Lévy jumps on population dynamics. Parameters are partially taken from numerical simulations in [18], r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, and c=1 with appropriate units.


### 4.1. Numerical Simulation for Deterministic System

By utilizing the above parameter values, it follows from (21) that there exists an interior equilibrium when 0<v<2.38. In order to facilitate the following analysis, v=0.5 is arbitrarily selected within 0,2.38 which is enough to merit the corresponding analysis in this paper, and an interior equilibrium M∗4.5384,1.7647,0.0231 can be obtained as v=0.5. Based on Theorem 1 and (A.10), if 0<τ1<τ1∗=1.4451 and 0<τ2<τ2∗=6.4223, then it is easy to show that interior equilibrium M∗4.5384,1.7647,0.0231 of system (22) is asymptotically stable. Based on the above analysis, some numerical simulations are supported. System (22) with τ1=0.75<τ1∗ and τ2=2.92 is locally asymptotically stable around M∗4.5384,1.7647,0.0231, whose dynamical responses are plotted in Figure 1(a). System (22) with τ1=τ1∗=1.4451 and τ2=2.92 is unstable around M∗4.5384,1.7647,0.0231, whose dynamical responses are plotted in Figure 1(b). Furthermore, system (22) with τ1=0.81 and τ2=3.65<τ2∗ is locally asymptotically stable around M∗4.5384,1.7647,0.0231, whose dynamical responses are plotted in Figure 2(a). System (22) with τ1=0.81 and τ2=τ2∗=6.4223 is unstable around M∗4.5384,1.7647,0.0231, whose dynamical responses are plotted in Figure 2(b).Figure 1
Parameter values are given as follows:r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, c=1, and v=0.5. When σjk=0j,k=1,2 and γiu=0i=1,2, dynamical responses of system (21) are plotted, where (a) τ1=0.75, τ2=2.92, and system (21) is locally asymptotically stable around 4.5384,1.7647,0.0231 and (b) τ1=1.4451, τ2=2.92, and system (21) is unstable stable around 4.5384,1.7647,0.0231.
(a)
(b)Figure 2
Parameter values are given as follows:r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, c=1, and v=0.5. When σjk=0j,k=1,2 and γiu=0i=1,2, dynamical responses of system (21) are plotted, where (a) τ1=0.81, τ2=3.65, and system (21) is locally asymptotically stable around 4.5384,1.7647,0.0231 and (b) τ1=0.81, τ2=6.4223, and system (21) is unstable stable around 4.5384,1.7647,0.0231.
(a)
(b)


### 4.2. Numerical Simulation for Stochastic System

In the absence of double time delays and Lévy jumps, by takingσ11=0.1 and σ21=0.5 as fixed values, it follows from simple computations that ϵ1=−0.8982, ϵ2=1.2, and ϵ4=0.88, which derive 28ϵ1+ϵ2−ϵ4=−12.8512<0 and 8ϵ1+ϵ2−ϵ4=−6.8656<0. Based on Theorems 3 and 4 of this paper, when v=0.5, probabilities and occurrence positions of stochastic Hopf bifurcation of system (27) with σ11=0.1 and σ21=0.5 and different values of σ12=σ22=0.1,0.2,0.3,0.5 are computed in Table 1. Corresponding to the computation results in Table 1, the probability density curve psω and occurrence positions of stochastic Hopf bifurcation of system (27) due to variations of ω are plotted in Figure 3. According to (40), the corresponding phase portraits of probability densities psω with respect to x1t and x2t are plotted in Figure 4.Table 1
Probabilities and occurrence positions of stochastic Hopf bifurcation of system (26) with σ11=0.1 and σ21=0.5 and different values of σ12 and σ22.
σ12=σ22=0.1
σ12=σ22=0.2
σ12=σ22=0.3
σ12=σ22=0.5
ϵ1
−0.8982
−0.8982
−0.8982
−0.8982
ϵ2
1.2
1.2
1.2
1.2
ϵ3
0.5023
0.5323
0.5823
0.7423
ϵ4
0.88
0.88
0.88
0.88
ω~
0.7650
0.7876
0.8237
0.93
psω~
0.9705
0.9427
0.9013
0.7983Figure 3
Parameter values for numerical simulations are given as follows:r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, c=1, and v=0.5. When γiu=0i=1,2 and τ1=τ2=0, the probability densities psω and occurrence positions of stochastic Hopf bifurcation of system (27) with σ11=0.1 and σ21=0.5 and different values of σ12 and σ22, where (a) σ12=σ22=0.1, (b) σ12=σ22=0.2, (c) σ12=σ22=0.3, and (d) σ12=σ22=0.5.Figure 4
Parameter values for numerical simulations are given as follows:r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, c=1, and v=0.5. The phase portraits of probability densities psω corresponding to Figure 3 in the x1−x2−psω space, where (a) σ12=σ22=0.1, (b) σ12=σ22=0.2, (c) σ12=σ22=0.3, and (d) σ12=σ22=0.5.
(a)
(b)
(c)
(d)In the presence of double time delays and Lévy jumps, if0<τ1<τ1∗ and 0<τ2<τ2∗, then it follows from Theorem 5 that limsupt→∞1/tE∫0tx1s−x1∗2ds≤B1 and limsupt→∞1/tE∫0tx2s−x2∗2ds≤B2, where B1 and B2 have been defined in Theorem 5. It follows from the mathematical form of B1 and B2 that σij and γii,j=1,2 affect asymptotic behaviors of system (3) around interior equilibrium M∗ as well as fluctuation intensity of population densities of x1t and x2t. By assuming τ1=0.8, τ2=3.5, v=0.5, σ11=0.1, and σ21=0.5, dynamical responses of system (3) with σ12=0.1 and σ22=0.1 and different values of γ1 and γ2 are obtained based on Euler-Maruyama method and plotted in Figure 5, where γ1=γ2=0.1 and γ1=γ2=0.5. Similarly, by assuming τ1=0.8, τ2=3.5, v=0.5, σ11=0.1, and σ21=0.5, dynamical responses of system (3) with σ12=0.5 and σ22=0.5 and different values of γ1 and γ2 are obtained based on Euler-Maruyama method and plotted in Figure 6, where γ1=γ2=0.1 and γ1=γ2=0.5.Figure 5
Parameter values for numerical simulations are given as follows:r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, c=1, and v=0.5. Dynamical responses of system (3) with σ11=0.1, σ21=0.5, σ12=0.1, and σ22=0.1 and different values of γ1 and γ2, where (a) γ1=γ2=0.1 and (b) γ1=γ2=0.5.
(a)
(b)Figure 6
Parameter values for numerical simulations are given as follows:r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, c=1, and v=0.5. Dynamical responses of system (3) with σ11=0.1, σ21=0.5, σ12=0.5, and σ22=0.5 and different values of γ1 and γ2, where (a) γ1=γ2=0.1 and (b) γ1=γ2=0.5.
(a)
(b)Based on the numerical simulations in Figures1 and 2, it follows that local asymptotic stability of system (22) around interior equilibrium switches due to variations of time delays. System (22) shows unstable dynamical responses when time delays cross critical values, which can be computed based on Theorem 1. It follows from the mathematical formulation of ϵ3 and numerical computational results in Table 1 that the value of ϵ3 increases as intensities of external excitations σ12 and σ22 values increase. Furthermore, it follows from Figure 3 the occurrence positions where system (27) undergoes stochastic Hopf bifurcation switch higher as the intensities of external excitations σ12 and σ22 values increase. On the other hand, the corresponding stochastic Hopf bifurcation will decrease in a higher probability, which is observed in Figure 4. Based on the numerical simulations in Figures 5 and 6, it can be concluded that the magnitude of environmental external stochastic excitation σ12 and σ22 values and Lévy jump γ1 and γ2 values plays significant roles to determine the magnitude of oscillation of population dynamics. It follows from the detailed numerical comparison results between Figures 5 and 6 that oscillation magnitude of dynamical responses increases as Lévy jump γ1 and γ2 values increase for the fixed environmental external stochastic excitation σ12 and σ22 values. On the other hand, oscillation magnitudes of dynamical responses increase as external stochastic excitation σ12 and σ22 values increase for the fixed environmental Lévy jump γ1 and γ2 values.


## 4.1. Numerical Simulation for Deterministic System

By utilizing the above parameter values, it follows from (21) that there exists an interior equilibrium when 0<v<2.38. In order to facilitate the following analysis, v=0.5 is arbitrarily selected within 0,2.38 which is enough to merit the corresponding analysis in this paper, and an interior equilibrium M∗4.5384,1.7647,0.0231 can be obtained as v=0.5. Based on Theorem 1 and (A.10), if 0<τ1<τ1∗=1.4451 and 0<τ2<τ2∗=6.4223, then it is easy to show that interior equilibrium M∗4.5384,1.7647,0.0231 of system (22) is asymptotically stable. Based on the above analysis, some numerical simulations are supported. System (22) with τ1=0.75<τ1∗ and τ2=2.92 is locally asymptotically stable around M∗4.5384,1.7647,0.0231, whose dynamical responses are plotted in Figure 1(a). System (22) with τ1=τ1∗=1.4451 and τ2=2.92 is unstable around M∗4.5384,1.7647,0.0231, whose dynamical responses are plotted in Figure 1(b). Furthermore, system (22) with τ1=0.81 and τ2=3.65<τ2∗ is locally asymptotically stable around M∗4.5384,1.7647,0.0231, whose dynamical responses are plotted in Figure 2(a). System (22) with τ1=0.81 and τ2=τ2∗=6.4223 is unstable around M∗4.5384,1.7647,0.0231, whose dynamical responses are plotted in Figure 2(b).Figure 1
Parameter values are given as follows:r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, c=1, and v=0.5. When σjk=0j,k=1,2 and γiu=0i=1,2, dynamical responses of system (21) are plotted, where (a) τ1=0.75, τ2=2.92, and system (21) is locally asymptotically stable around 4.5384,1.7647,0.0231 and (b) τ1=1.4451, τ2=2.92, and system (21) is unstable stable around 4.5384,1.7647,0.0231.
(a)
(b)Figure 2
Parameter values are given as follows:r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, c=1, and v=0.5. When σjk=0j,k=1,2 and γiu=0i=1,2, dynamical responses of system (21) are plotted, where (a) τ1=0.81, τ2=3.65, and system (21) is locally asymptotically stable around 4.5384,1.7647,0.0231 and (b) τ1=0.81, τ2=6.4223, and system (21) is unstable stable around 4.5384,1.7647,0.0231.
(a)
(b)


## 4.2. Numerical Simulation for Stochastic System

In the absence of double time delays and Lévy jumps, by takingσ11=0.1 and σ21=0.5 as fixed values, it follows from simple computations that ϵ1=−0.8982, ϵ2=1.2, and ϵ4=0.88, which derive 28ϵ1+ϵ2−ϵ4=−12.8512<0 and 8ϵ1+ϵ2−ϵ4=−6.8656<0. Based on Theorems 3 and 4 of this paper, when v=0.5, probabilities and occurrence positions of stochastic Hopf bifurcation of system (27) with σ11=0.1 and σ21=0.5 and different values of σ12=σ22=0.1,0.2,0.3,0.5 are computed in Table 1. Corresponding to the computation results in Table 1, the probability density curve psω and occurrence positions of stochastic Hopf bifurcation of system (27) due to variations of ω are plotted in Figure 3. According to (40), the corresponding phase portraits of probability densities psω with respect to x1t and x2t are plotted in Figure 4.Table 1
Probabilities and occurrence positions of stochastic Hopf bifurcation of system (26) with σ11=0.1 and σ21=0.5 and different values of σ12 and σ22.
σ12=σ22=0.1
σ12=σ22=0.2
σ12=σ22=0.3
σ12=σ22=0.5
ϵ1
−0.8982
−0.8982
−0.8982
−0.8982
ϵ2
1.2
1.2
1.2
1.2
ϵ3
0.5023
0.5323
0.5823
0.7423
ϵ4
0.88
0.88
0.88
0.88
ω~
0.7650
0.7876
0.8237
0.93
psω~
0.9705
0.9427
0.9013
0.7983Figure 3
Parameter values for numerical simulations are given as follows:r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, c=1, and v=0.5. When γiu=0i=1,2 and τ1=τ2=0, the probability densities psω and occurrence positions of stochastic Hopf bifurcation of system (27) with σ11=0.1 and σ21=0.5 and different values of σ12 and σ22, where (a) σ12=σ22=0.1, (b) σ12=σ22=0.2, (c) σ12=σ22=0.3, and (d) σ12=σ22=0.5.Figure 4
Parameter values for numerical simulations are given as follows:r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, c=1, and v=0.5. The phase portraits of probability densities psω corresponding to Figure 3 in the x1−x2−psω space, where (a) σ12=σ22=0.1, (b) σ12=σ22=0.2, (c) σ12=σ22=0.3, and (d) σ12=σ22=0.5.
(a)
(b)
(c)
(d)In the presence of double time delays and Lévy jumps, if0<τ1<τ1∗ and 0<τ2<τ2∗, then it follows from Theorem 5 that limsupt→∞1/tE∫0tx1s−x1∗2ds≤B1 and limsupt→∞1/tE∫0tx2s−x2∗2ds≤B2, where B1 and B2 have been defined in Theorem 5. It follows from the mathematical form of B1 and B2 that σij and γii,j=1,2 affect asymptotic behaviors of system (3) around interior equilibrium M∗ as well as fluctuation intensity of population densities of x1t and x2t. By assuming τ1=0.8, τ2=3.5, v=0.5, σ11=0.1, and σ21=0.5, dynamical responses of system (3) with σ12=0.1 and σ22=0.1 and different values of γ1 and γ2 are obtained based on Euler-Maruyama method and plotted in Figure 5, where γ1=γ2=0.1 and γ1=γ2=0.5. Similarly, by assuming τ1=0.8, τ2=3.5, v=0.5, σ11=0.1, and σ21=0.5, dynamical responses of system (3) with σ12=0.5 and σ22=0.5 and different values of γ1 and γ2 are obtained based on Euler-Maruyama method and plotted in Figure 6, where γ1=γ2=0.1 and γ1=γ2=0.5.Figure 5
Parameter values for numerical simulations are given as follows:r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, c=1, and v=0.5. Dynamical responses of system (3) with σ11=0.1, σ21=0.5, σ12=0.1, and σ22=0.1 and different values of γ1 and γ2, where (a) γ1=γ2=0.1 and (b) γ1=γ2=0.5.
(a)
(b)Figure 6
Parameter values for numerical simulations are given as follows:r1=2, b1=0.4, a1=0.5, k=6, q=0.1, m1=1, m2=2, r2=0.3, a2=2, w=5, c=1, and v=0.5. Dynamical responses of system (3) with σ11=0.1, σ21=0.5, σ12=0.5, and σ22=0.5 and different values of γ1 and γ2, where (a) γ1=γ2=0.1 and (b) γ1=γ2=0.5.
(a)
(b)Based on the numerical simulations in Figures1 and 2, it follows that local asymptotic stability of system (22) around interior equilibrium switches due to variations of time delays. System (22) shows unstable dynamical responses when time delays cross critical values, which can be computed based on Theorem 1. It follows from the mathematical formulation of ϵ3 and numerical computational results in Table 1 that the value of ϵ3 increases as intensities of external excitations σ12 and σ22 values increase. Furthermore, it follows from Figure 3 the occurrence positions where system (27) undergoes stochastic Hopf bifurcation switch higher as the intensities of external excitations σ12 and σ22 values increase. On the other hand, the corresponding stochastic Hopf bifurcation will decrease in a higher probability, which is observed in Figure 4. Based on the numerical simulations in Figures 5 and 6, it can be concluded that the magnitude of environmental external stochastic excitation σ12 and σ22 values and Lévy jump γ1 and γ2 values plays significant roles to determine the magnitude of oscillation of population dynamics. It follows from the detailed numerical comparison results between Figures 5 and 6 that oscillation magnitude of dynamical responses increases as Lévy jump γ1 and γ2 values increase for the fixed environmental external stochastic excitation σ12 and σ22 values. On the other hand, oscillation magnitudes of dynamical responses increase as external stochastic excitation σ12 and σ22 values increase for the fixed environmental Lévy jump γ1 and γ2 values.


## 5. Conclusion

In this paper, a double delayed hybrid stochastic bioeconomic system with commercial harvesting and Lévy jumps is established, which extends work done in [18] by incorporating double time delays and stochastic fluctuations. The dynamical model proposed in [18] is composed of ordinary differential equations, which are utilized to study the interaction mechanism of the prey-predator system with single time delay. Compared with the system established in [18], an algebraic equation is introduced into system (3), which concentrates on the dynamic effect of the economic interest of commercial harvesting on population dynamics and provides a straightforward way to investigate complex dynamics due to the variation of economic interests. Furthermore, two discrete time delays, which represent the maturation delay of the prey and gestation delay of the predator, are incorporated into system (3). Population growth of the prey and predator population affected by environmental stochastic fluctuations is assumed to be a stochastic process. Gaussian white noises and Lévy jumps are incorporated into system (3) to describe stochastic surrounding environmental factors. For the deterministic system, positivity of solutions and uniform persistence of system are discussed in Lemmas 1 and 2, respectively. Some sufficient conditions associated with double time delays are derived to discuss the asymptotic stability of interior equilibrium in Theorem 1. For the stochastic system, existence and uniqueness of global positive solution are studied in Theorem 2. By using the invariant measure theory and singular boundary theory of the diffusion process, the existence of stochastic Hopf bifurcation and stochastic stability is investigated in the absence of double time delays and Lévy jumps, which can be found in Theorems 3 and 4, respectively. By constructing appropriate Lyapunov functions, the asymptotic dynamic behavior of the proposed hybrid stochastic system with double time delays and Lévy jumps is discussed in Theorem 5.From the practical perspective of viewpoints, positivity of the proposed hybrid stochastic double delayed bioeconomic system is relevant to the prey and predator population survival for a long duration under a commercially harvested environment. Generally speaking, the natural resources for the prey and predator population survival are under sever intraspecies competition. Therefore, permanence of the proposed hybrid stochastic double delayed bioeconomic system biologically means that there are certain positively finite upper constraints and positively finite lower constraints for the prey and predator population, which may, to a great extent, avoid population extinction and overpopulation. Furthermore, permanence of the proposed hybrid stochastic double delayed bioeconomic system biologically also means that there exist positively finite upper constraints and positively finite lower constraints for commercial harvesting amount on the prey population. The practical interpretations are introduced as follows: with the purpose of maintaining sustainable development of commercially harvested prey population resources, commercial harvesting should be regulated within certain harvesting range by formulating some constructive policies.From the practical perspective of viewpoints, if the proposed hybrid stochastic double delayed bioeconomic system is unstable in terms of probability, then it practically interprets that the prey and predator population density will not eventually arrive at an ecological balance level with a high probability. The biological interpretations show that such phenomenon is not advantageous for a sustainable development of interacting populations under a commercial harvest effect. If the boundary is attractively natural and stochastic Hopf bifurcation does not occur, then it practically means that the interacting population densities may dramatically increase with a high probability during short time and it may reduce the corresponding survival resources and space, which are beyond the environment carrying capacity and destroy the ecological balance state with a high probability.From the resource management of viewpoints, the numerical simulations reveal that population density may remain an ideal level by controlling double time delays within certain constraints and the critical values of time delays may increase as economic interests of commercial harvesting increase. The sustainable development of commercially harvested population can be indirectly achieved by formulating relevant policy to regulate an economic interest within some appropriate ranges.Recently, some hybrid bioeconomic systems with time delay and stochastic fluctuations are established in [27–30] to investigate the combined dynamic effects of stochastic fluctuation and commercial harvesting on population dynamics. These proposed systems in [27–30] are also constructed by several differential equations with stochastic fluctuations and an algebraic equation. Compared with the previously bioeconomic systems proposed in [19–23, 31] and the references therein, such delayed bioeconomic systems [27–30] can not only discuss coexistence and interaction mechanism of delayed bioeconomic system under stochastic environmental fluctuations but also investigate population dynamics due to variations of the economic interest of commercial harvesting. However, biological characteristics among interacting populations are not considered in [27–29]; time delays such as gestation delay and maturation delay for interacting populations in [27–29] are assumed to be the same discrete value, which contradicts to the reality in the real world. Asymptotical stability of interior equilibrium and dynamic effects of Lévy jumps on population dynamics are not discussed in [27–29]. It is proved that Lévy noise can efficiently depict sudden and severe environmental perturbations arising in the real world [24, 25], while these phenomena can not be described better by Brownian motion. Although the dynamic effects of double time delays have been investigated in [30], the combined dynamic effects of multiple time delays and Lévy jumps on population dynamics have not been investigated in [30]. Dynamic effects of multiple time delays on the hybrid bioeconomic prey-predator system are investigated in [15, 32]. However, the dynamic effect of Lévy jumps and asymptotical stability of solutions of the stochastic system are not studied in [15, 32]. To the authors’ best knowledge, population dynamics of the hybrid bioeconomic system with double time delays and Lévy jumps have not been investigated. Compared with the related work, we can investigate combined dynamic effects of double time delays and Lévy jumps on population dynamics by analyzing stability analysis and stochastic dynamical behavior of system (3) in this paper; these analytical findings make this paper have some positive and new features.



---

*Source: 8294153-2018-08-13.xml*