Dynamics and optimal harvesting of a stochastic predator-prey system with regime switching, S-type distributed time delays and L´evy jumps

: This work is concerned with a stochastic predator-prey system with S-type distributed time delays, regime switching and L´evy jumps. By use of the stochastic di ﬀ erential comparison theory and some inequality techniques, we study the extinction and persistence in the mean for each species, asymptotic stability in distribution and the optimal harvesting e ﬀ ort of the model. Then we present some simulation examples to illustrate the theoretical results and explore the e ﬀ ects of regime switching, distributed time delays and L´evy jumps on the dynamical behaviors, respectively.


Introduction
The natural world is full of environmental perturbations almost everywhere, which play an important role in the ecological system [1][2][3][4][5]. As May [6] said that the growth rates in ecological systems should be stochastic due to the influences of random noises, so the solution would fluctuate around some average values, not being a steady positive point. Therefore, it is more rational by considering a stochastic environmental perturbation (white noise) to better investigate the real population systems. In mathematical modelling, the white noise has been extensively introduced to reveal the stochastic population dynamics, see e.g. [3][4][5]7,8] and references cited therein.
In real world, population systems may inevitably suffer some abrupt changes. It is well known that the growth rates of some species are affected by seasonal factors, which are much different in summer from those in winter. These phenomena can be described by a continuous-time Markovian process with a finite state space, i.e. colorful noise in mathematical modelling. The colorful noise may take several values and switch among different regimes of environment, which is memoryless, and the waiting time for the next switching follows an exponential distribution. The effect of colorful noise on population dynamics has attracted many researchers [9][10][11][12][13][14]. For example, by using a Markovian switching process to model the colorful noise in environment, Liu, He and Yu [14] proposed the following stochastic system with harvesting and regime-switching: x 2 (t + s)dµ 1 (s) dt + σ 1 (α(t))x 1 (t)dB 1 (t), dx 2 (t) = x 2 (t) r 2 (α(t)) − h 2 + a 21 0 −τ 2 x 1 (t + s)dµ 2 (s) − a 22 x 2 (t) dt + σ 2 (α(t))x 2 (t)dB 2 (t), (1.1) where r 1 (·) > 0 is the growth rate of prey, r 2 (·) < 0 is the death rate of predator, a ii > 0 is the intraspecific competition rate, a 12 > 0 is the capture rate and a 21 > 0 is the conversion rate of food; h i is the harvesting constant; B 1 (t) and B 2 (t) are standard independent Brownian motions defined on a complete probability space (Ω, F , {F t } t≥0 , P), where {F t } t≥0 is a filtration. It is right continuous and F 0 contains all P-null sets; σ 2 i (·) is the intensity of stochastic noise, i = 1, 2. The regime of switching α(t) is a Markovian chain in a finite state space S = {1, 2, . . . , N}. The generator of α(t) is defined as where t > 0, γ i j is the transition rate from the ith stage to the jth stage and γ i j ≥ 0 if i j while γ ii = − i j γ i j , i, j ∈ S. It is often assumed that every sample of α(t) is a right continuous step function and irreducible with a finite simple jumps in any finite subinterval of R + = [0, ∞). It obeys a unique stationary distribution π = (π 1 , π 2 , . . . , π N ) satisfying πΓ = 0 and N k=1 π k = 1, π k > 0, k ∈ S. The switching mechanism of the hybrid system is referred to [9].
In the study of ecological dynamics, the current growth of populations is usually influenced by its past history, that is, time delay is often inevitable in the natural ecosystems [15]. S-type distributed time delay is such a distributed delay that the integral is Lebesgue-Stieltjes. Just as the authors [16] said, "systems with discrete time delays and those with continuously distributed time delays do not contain each other. However, systems with S-type distributed time delays contain both." So it is interesting to consider the impact of S-type distributed delay on the ecological dynamics. Models with S-type distributed time delays have been studied by many authors [17,18]. On the other hand, earthquake, harvesting and epidemics often happen in natural world. These sudden environmental perturbations are so strong and can change the population size in a very short time, which can not be described by white noise [19,20]. Many experiments show that, due to the influence of environmental disturbance, the distribution of many species exhibits a scale-free characteristic, and hence the biologists introduce a non-Gaussian Lévy jump to characterize it. Models with Lévy jump are studied by many researchers and many nice results have been obtained [21][22][23].
Considering the effects of S-type time delays and Lévy jumps on system (1.1), we establish the following stochastic model The Markov chain, Brownian motion and Lévy jumps are mutually independent.
Our main goal of this paper is as follows. First, since the study of dynamical behaviors of predator-prey system is an important topic [7,24], we establish some sufficient conditions assuring the extinction, persistence in the mean for all species of system (1.2).
Second, in the study of long-term behaviors of species, the existence of a unique probability measure plays an important role in stochastic models with Lévy jumps (see [25][26][27]). Hence, it is very interesting to analyze the asymptotic stability in distribution of (1.2).
Third, in mathematical biology, it is valuable to keep the species persistent to maintain the biological balance. Consequently, the optimal harvesting strategy of renewable resources becomes more and more important [8,[27][28][29]. By ergodic method, we will study the optimal harvesting strategy of system (1.2).
The rest of this paper is organized as follows. Section 2 begins with some definitions, important lemmas and notations. Section 3 is devoted to the extinction and persistence in the mean for species. Section 4 and Section 5 focus on the asymptotic stability in distribution and the optimal harvesting strategy of (1.2), respectively. Section 6 gives some numerical simulations to verify the main results. Finally, Section 7 presents a brief discussion to conclude this paper.

Preliminaries
To begin with this section, we introduce some notations of the Itô's integral for a stochastic differential equation with Markovian switching and Lévy jumps [19,23]. Let where f and g : . Define the operator L as follows: The generalized Itô's formula with jumps (see, for example [19,23]) is defined as Now we give the definitions of extinction and persistence in the mean of each species, and an important comparison theorem.
If there exist three positive constants T, λ 0 and λ such that, for all t > T , Further, for the need of our discussion, we give some technical assumptions.
Assumption 1. For any α ∈ S and i = 1, 2, we assume γ i (α, u) > −1 and where K > 0 denotes a positive and finite constant unless otherwise stated, which may be different in different places.
Assumption 2. For any p > 0, we assume Assumptions 1 and 2 imply that the intensity of Lévy noise cannot be too strong, otherwise, the solution of system (1.2) may explode in some finite time [23]. About the existence of nontrivial positive solutions of (1.2), we have the following lemma.
Proof. The proof of the existence of solutions is straightforward, one may refer to [14,23]. As to the proof of the boundedness of expectation, one can get it by using the generalized Itô's formular with jumps to function e t 2 i=1 x p i (t). The process is similar to reference [14] and is omitted. For simplicity, we introduce some notations to end this section.

Extinction and persistence in the mean of species
In this section, we study the long-term behaviors of (1.2). Consider the following system, Lemma 3.1. For system (3.1), the following statements hold.
Using the Itô's formula to ln W(t) and computing the stochastic differential along the solution of (3.1), we have Integrating both sides of (3.2) from 0 to t, and divided by t from both sides, then On the other hand, by the ergodicity of Markovian chain, we have Now we give the proof of (i) and (ii).
(i) By comparison method, using Lemma 2.1, we can easily derive from (3.4) that lim sup . .
The proof is completed. Next, we consider the following comparison system of (1.2): Consider the following comparison system of the second equation of (3.5), Similar with the previous reasoning, using Lemma 3.1 again, we can derive that By the comparison theory [1,9], it is clear that x 1 ≤ y 1 and y 2 ≤ŷ 2 . Consequently, we have the following results.
Lemma 3.2. For system (3.5), the following statements of Table 1 hold.

Conditions
Species y 1 Species y 2 η 1 < 0 Extinction Extinction Now we give the main result on the extinction and persistence in the mean of the species of system (1.2).

Cases
Conditions Species x 1 Species x 2 iη 1 < 0 Extinction Extinction and Similar with the proof of Lemma 3.1, we have Integrating both sides of (3.6) from 0 to t, and combining (3.8) and (3.9), then (3.10) By the same argumentation, we have (ii) Ifη 1 > 0, then we derive from Lemma 3.2 that lim sup . Substituting it into (3.11) and using Lemma 2.1 gives lim sup By the condition ∆ 2 < 0, then lim t→∞ x 2 (t) = 0.
From (3.11) again, by the comparison theorem [4,14,27], we have lim inf x 2 (s)ds where M i (t) = On the other hand, by (3.10) and (ii), we can derive that lim x 2 (s)ds = Substituting (3.12) into (3.10) and using Lemma 2.1 again, we can obtain that Consequently, we have lim The proof is completed. Remark 3.2. Theorem 3.1 implies that the regime switching, time delays and Lévy jumps will bring large influence to the dynamics of (1.2). By changing their values, the species of (1.2) may change from persistence in the mean to extinction, and vice versa, which is analyzed and verified well by numerical simulations in Section 6.

Asymptotically stable in distribution of (1.2)
For simplicity, we denote the solution of (1.2) with initial data x(θ) = φ(θ), α(θ) = ς by x(t, φ, ς). Let x(t, φ, ς) be the R 2 + ×S−valued stochastic Markovian process. Let B ⊆ R 2 + be a Borel measurable set and D ⊆ S, and we denote the transmission probability of the event {x(t, φ, ς) ∈ B × D} by P(t, φ, ς, B × D), that is, P(t, φ, ς, B × D) = u∈D B P(t, φ, ς, dx × {u}). Denote by P(R 2 + , S) all the probability measures on R 2 + × S, and for any two P 1 , P 2 ∈ P(R 2 + , S), we define the metric d L as follows: where The process x(t, φ, ς) is said to be asymptotically stable in distribution if there exists a probability measure u(·×·) on R 2 + ×S such that the transmission probability P(t, φ, ς, dx×{k}) of 2) is said to be asymptotically stable in distribution if the solution x(t, φ, ς) of (1.2) is asymptotically stable in distribution. For the need of discussion, we give the following technical assumption.
Assumption 3 means that under the effect of time delays, the intraspecific competition rate is still greater than the interaction rate between different species. Theorem 4.1. Let x(t) = x(t, φ, ς) and x(t) = x(t, ϕ, ς) be two solutions of (1.2) with initial value x(θ) = φ, x(θ) = ϕ and α(θ) = ς, respectively. If Assumption 3 holds, then we have We calculate the right differential of V 1 (t) along the solutions of (1.2), then Similarly, we have Define V 3 (t) as follows:

then by Assumption 3, after computation, we have
Hence, On the other hand, by (1.2) we have, Taking expectation from both sides of (4.2), then Consequently, E(x i (t)) is continuously differentiable with respect to t. Moreover, That is to say, E(x i (t)) is uniformly continuous. An application of Lemma 4.1 gives lim t→∞ E|x i (t) − x i (t)| = 0. This completes the proof.
Remark 4.1. Theorem 4.1 will be applied later to prove the stability in distribution of (1.2). Assumption 3 is necessary in our proof. If without S-type time delays and Lévy jumps, Reference [14] implies that Assumption 3 is unnecessary and may be dropped.
Remark 4.2. The asymptotic stability in distribution of species reveals the existence and uniqueness of an invariant probability measure, which is the basis of discussing the optimal harvesting effort in Section 5.
(1) By the given conditions, it is clear that λ ∈ Ξ, that is, the set Ξ is not empty. For any h ∈ Ξ, by Theorem 3.1, after computation, we have Theorem 4.1 implies that the distribution of solutions of (1.2) is asymptotically stable, and hence suppose the stationary probability density is ρ, then The asymptotical stability in distribution means that there exists a unique invariant measure u. In view of the one-to-one correspondence between the stationary probability density ρ and the invariant measure u (see, e.g. [31], page 105), we have which is negative semi-definite for all h by the given conditions, then we derive from Theorem 4.1.5 of [32] that λ is a global maximum point of Y(h). That is, if the conditions ∆ 2 | h i =λ i > 0, λ i ≥ 0(i = 1, 2) hold, then the OHE is h * = λ, and Y(h * ) = λ T A −1 (ξ − λ).
Remark 5.1. The existence and uniqueness of an invariant probability measure plays a key role to derive (5.4). With the invariant measure, by using the extremum theory, we find the maximum of Y(h), which is very popular and can be applied to resolve some similar problems.

Examples and simulations
In this section, by numerical analysis, we give some examples and apply the Milstein method [33] to illustrate our theoretical results, and explore the effects of regime switching, time delays and Lévy jumps on the system dynamics.
• Case (i) The effect of switching π.
(1) Let the stationary distribution π = (0.1, 0.9). By computation, then That is, both x 1 (t) and x 2 (t) are persistent in the mean. We call it the "persistent case". Figure  2 (a) is the long-term behaviours and Figure 2 (b) is the probability densities of x 1 (t) and x 2 (t), respectively. Figure 3 gives the time series and histogram of Markov chain α(t) with stationary distribution π = (0.1, 0.9).  For the persistent case, Theorem 5.1 shows that there exists OHE, and The maximum of the sustainable yield is By numerical simulation, we get Figure 2 (c). In Figure 2 (c), the red line represents the yield Y = 0.0669 with λ = (0.3258, 0.0532) T , the green line represents the yield Y 1 = 0.0484 with λ 1 = (0.1558, 0.0632) T , the blue line represents the yield Y 2 = 0.0601 with λ 2 = (0.4258, 0.0332) T , respectively. Obviously, the maximum of sustainable yield is Y(h * ) = 0.0669. Figure 2 (c) verify it well. For the interpretation of the references to color, readers are referred to the web version of the article.
• Case (iii)The effect of Lévy jumps.

Conclusions and discussions
In this paper, we study the dynamics of a stochastic predator-prey system with S-type distributed time delays, regime switching and Lévy jumps. Theorem 3.1 gives the sufficient conditions assuring the extinction and persistence in the mean of each species. Theorem 4.2 shows that (1.2) is asymptotically stable in distribution. Theorem 5.1 gives the optimal harvesting effort. Finally, some examples are given and the effects of regime switching, distributed time delays and Lévy jumps are discussed by numerical analysis. Figure 4 (a) and (b) imply that different regime switching may lead to very different long-term behaviours of species. Figure 2 (c) shows that the OHE is based on the suitable regime switching. Figure 4 (c) and (d) show that too large delays or Lévy jumps will destroy the persistent case of (1.2), respectively. All these show that regime switching, distributed time delays and Lévy jumps play key role in the dynamics of (1.2).
Further, (1.2) is very popular and contains many researched model as its special cases. Firstly, (1.1) is contained in (1.2), and hence, one can obtain the sufficient conditions of the extinction and persistence in the mean for species of (1.1), that is, Theorems 1 and 2 of Reference [14] are contained in our results. Secondly, if α(t) = 1, h i = 0, γ i = 0(i = 1, 2), then we get the model studied by Wang et.al. [18]. Our result coincides with Theorem 2.2 of [18]. Thirdly, if µ ii (θ) are constant on [−τ, 0], a i j = 0, i j, and then we get the discrete time delays model proposed in Reference [5,8]. Similarly we can obtain Theorem 2.2 in [8] and Lemma 3 in [5]. As Liu et. al. stated in [14], the growth of the ith species at time t is often affected by the abundance of the jth species on the interval [t − τ i , t], rather than only on the time t − τ i . Hence the S-type distributed delays can fit with some real biological systems better. In above sense, we improve and generalize the obtained conclusions of [5,8,14,18]. However, the switching does not appear in all parameters and the control inputs are all constants, if they are dependent on time t, then the persistence in the mean and the optimal harvesting strategy can not be established and is still unknown. On the other hand, for predator-prey system, if the predator is provided with additional food [34], or the fear of prey induced by predator appears [35], what will happen is very interesting. All these will be our research work in the future.