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

1.   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:

where $r_1(\cdot) > 0$ is the growth rate of prey, $r_2(\cdot) < 0$ is the death rate of predator, $a_{ii} > 0$ is the intra-specific 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 $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\geq0}, P)$, where $\{\mathcal{F}_t\}_{t\geq0}$ is a filtration. It is right continuous and $\mathcal{F}_0$ contains all P-null sets; $\sigma_i^2(\cdot)$ is the intensity of stochastic noise, $i = 1, 2.$ The regime of switching $\alpha(t)$ is a Markovian chain in a finite state space $\mathbb{S} = \{1, 2, \ldots, N\}$. The generator of $\alpha(t)$ is defined as $\Gamma = (\gamma_{ij})_{N\times N}$ with

where $\vartriangle t > 0, \gamma_{ij}$ is the transition rate from the $i$th stage to the $j$th stage and $\gamma_{ij}\geq0$ if $i\neq j$ while $\gamma_{ii} = -\sum_{i\neq j}\gamma_{ij}, i, j\in \mathbb{S}$.

It is often assumed that every sample of $\alpha(t)$ is a right continuous step function and irreducible with a finite simple jumps in any finite subinterval of $R_+ = [0, \infty)$. It obeys a unique stationary distribution $\pi = (\pi_1, \pi_2, \ldots, \pi_N)$ satisfying $\pi\Gamma = 0$ and $\sum_{k = 1}^N\pi_k = 1, \pi_k > 0, k\in \mathbb{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

with initial value

where $x_i(t^-)$ stands for the left limit of $x_i(t)$; $\int_{-\tau_{ij}}^0x_i(t+\theta)d\mu_{ij}(\theta)$ is Lebesgue-Stieltjes integral; $\tau_{ij} > 0$ is time delay; $\mu_{ij}(\theta)$ is bounded, nondecreasing variation function defined on $[-\tau, 0]$ with $\tau = \max_{i, j = 1, 2}\{\tau_{ij}\}, \varsigma\in\mathbb{S}, i, j = 1, 2;$ $R_+^2 = \{x = (x_1, x_2)\in R^2, x_i > 0, i = 1, 2\}$ with $|x(t)| = \sqrt{\sum_{i = 1}^2x_i^2};$ $N$ is a Poisson counting measure, $\widetilde{N}(dt, du) = N(dt, du)-\lambda(du)$ is the component of $N,$ where $\lambda$ is the characteristic measure on a measurable subset $\mathbb{Z}\subset R_+ = [0, \infty)$ such that $\lambda(\mathbb{Z}) < \infty$. 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.

2.   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\; \mbox{and}\; g: R^2\times R_+\times \mathbb{S}\to R^2, h: R^2\times R_+\times \mathbb{S}\times\mathbb{Z}\to R^2$ are measurable functions. Let $V\in C^{2, 1}(R^2\times R_+\times \mathbb{S}, R^2)$. Define the operator $L$ as follows:

where $V_t(x, t, k) = \frac{\partial V(x, t, k)}{\partial t}, V_x(x, t, k) = (\frac{\partial V(x, t, k)}{\partial x_1}, \frac{\partial V(x, t, k)}{\partial x_2}), V_{xx}(x, t, k) = (\frac{\partial^2 V(x, t, k)}{\partial x_i\partial x_j})_{2\times 2}, i, j = 1, 2.$ 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.

Definition 2.1 (^[27]). Let $x(t) = (x_1(t), x_2(t))^T\in R_+^2$ be a solution of system (1.2). Then,

$(a)$ the population $x_i(t)$ is said to be extinct if $\lim\limits_{t\rightarrow \infty}x_i(t) = 0, i = 1, 2$;

$(b)$ the population $x_i(t)$ is said to be persistent in the mean if $\lim\limits_{t\to\infty} \frac{1}{t}\int_0^t x_i(s)ds = K$ a.s., where $K$ is a positive constant, $i = 1, 2.$

Lemma 2.1 (^[4]). Suppose that $Z(t)\in C[\Omega\times[0, +\infty), R_+]$ and $\lim\limits_{t\rightarrow \infty}F(t)/t = 0, a.s.$

${\rm (a)}$ If there exist two positive constants $T$ and $\lambda_0$ such that, for all $t > T$,

then

${\rm (b)}$ If there exist three positive constants $T, \lambda_0$ and $\lambda$ such that, for all $t > T$,

then

Further, for the need of our discussion, we give some technical assumptions.

Assumption 1. For any $\alpha\in \mathbb{S}$ and $i = 1, 2,$ we assume $\gamma_i(\alpha, 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.

Lemma 2.2. Let Assumptions 1 and 2 hold. For any given initial value $x(\theta)\in C([-\tau, 0];R_+^2), \alpha(\theta) = \varsigma,$ system (1.2) has a unique solution$x (t)$ on $t\in[-\tau, \infty),$ and the solution remains in $R_+^2$ with probability one. Moreover, for any $p > 0$, there is a $K(p) > 0$such that $\limsup_{t\rightarrow \infty}\mathbb{ E}|x(t)|^p\leq K(p).$

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\sum_{i = 1}^2 x_i^p(t).$ The process is similar to reference ^[14] and is omitted.

For simplicity, we introduce some notations to end this section.

3.   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.

${(i)}$ If $\sum_{k = 1}^N\pi_k\xi(k) < 0,$ then $\lim\limits_{t\to \infty}W(t) = 0.$

${(ii)}$ If $\sum_{k = 1}^N\pi_k\xi(k) > 0,$ then$\lim\limits_{t\to \infty} \frac{1}{t} \int_0^tW(s)ds = \frac{\sum_{k = 1}^N\pi_k\xi(k)}{a+\int_{-\tau}^0d\mu(\theta)},$

where $\xi(k) = r(k)-\frac{\sigma^2(k)}{2} -\int_Z[\gamma(k, u)-\ln(1+\gamma(k, u))]\lambda(du).$

Proof. 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

Since

it follows that

Then

Therefore (3.3) leads to (3.4) as follows.

On the other hand, by the ergodicity of Markovian chain, we have

Now we give the proof of (ⅰ) and (ⅱ).

(ⅰ) By comparison method, using Lemma 2.1, we can easily derive from (3.4) that

If $\sum_{k = 1}^N\pi_k\xi(k) < 0$, it is clear that $\lim\limits_{t\to \infty}W(t) = 0.$

(ⅱ) If $\sum_{k = 1}^N\pi_k\xi(k) > 0,$ by Lemma 2.1 again, we have

Therefore,

The proof is completed.

Next, we consider the following comparison system of (1.2):

By Lemma 3.1, we derive from (3.5) that

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\leq y_1$ and $y_2\leq \hat{y}_2.$ Consequently, we have the following results.

Lemma 3.2. For system (3.5), the following statements of Table 1 hold.

Now we give the main result on the extinction and persistence in the mean of the species of system (1.2).

Theorem 3.1. For system (1.2), the following statements hold (see the Table 2).

Proof. For system (1.2), by using the Itô's formula to $\ln x_i(t), i = 1, 2,$ we have

and

Similar with the proof of Lemma 3.1, we have

and

Integrating both sides of (3.6) from 0 to $t$, and combining (3.8) and (3.9), then

By the same argumentation, we have

By the ergodicity of Markovian chain, then

Now we prove the conclusion of Theorem 3.1.

(ⅰ) If $\bar{\eta}_1 < 0,$ then by Lemma 3.2, we have $\lim\limits_{t\to \infty}x_1(t) = 0.$ Since $\bar{\eta}_2 < 0,$ it follows from (3.11) that $\lim\limits_{t\to \infty}x_2(t) = 0.$ Hence $\lim\limits_{t\to \infty}x_i(t) = 0, i = 1, 2.$

(ⅱ) If $\bar{\eta}_1 > 0$, then we derive from Lemma 3.2 that $\limsup\limits_{t\to \infty} \frac{1}{t}\int_0^tx_1(s)ds\leq\frac{\bar{\eta}_1}{A_{11}}.$ Substituting it into (3.11) and using Lemma 2.1 gives

By the condition $\Delta_2 < 0$, then $\lim\limits_{t\to \infty}x_2(t) = 0.$

From (3.11) again, by the comparison theorem ^[4,14,27], we have $\liminf\limits_{t\to \infty} \frac{1}{t}\int_0^tx_1(s)ds\geq\frac{\bar{\eta}_1}{A_{11}}.$ Therefore,

(ⅲ) If $\Delta_2 > 0$, it is clear that $\bar{\eta}_1 > 0$. We compute $(3.10)\times A_{21}+(3.11)\times A_{11},$ then

where $M_i(t) = \int_0^t\sigma_i(\alpha(s))dB_i(s) + \int_0^t\int_\mathbb{Z}\ln(1+\gamma_i(u, \alpha(s)))\tilde{N}(ds, du), i = 1, 2.$ By the strong law of large numbers ^[9], we have

On the other hand, by (3.10) and (ⅱ), we can derive that $\lim\limits_{t\rightarrow \infty} \frac{\ln x_1(t)}{t} = 0.$ By comparison method (using Lemma 2.1 again), then

Substituting (3.12) into (3.10) and using Lemma 2.1 again, we can obtain that

Consequently, we have $\lim_{t\rightarrow \infty} \frac{1}{t}\int_0^tx_i(s)ds = \frac{\Delta_i}{\Delta}, i = 1, 2$. The proof is completed.

Remark 3.1. Compared with (1.1), model (1.2) is more popular and contains (1.1) as its special case. Another difference between (1.1) and (1.2) is that the Lévy jump is considered in (1.2) while it is not considered in (1.1).

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.

4.   Asymptotically stable in distribution of (1.2)

For simplicity, we denote the solution of (1.2) with initial data $x(\theta) = \phi(\theta), \alpha(\theta) = \varsigma$ by $x(t, \phi, \varsigma)$. Let $x(t, \phi, \varsigma)$ be the $R_+^2\times \mathbb{S}-$valued stochastic Markovian process. Let $\mathbb{B}\subseteq R_+^2$ be a Borel measurable set and $\mathbb{D}\subseteq \mathbb{S}$, and we denote the transmission probability of the event $\{x(t, \phi, \varsigma)\in \mathbb{B}\times \mathbb{D}\}$ by $P(t, \phi, \varsigma, \mathbb{B}\times \mathbb{D}),$ that is, $P(t, \phi, \varsigma, \mathbb{B}\times \mathbb{D}) = \sum_{u\in \mathbb{D}}\int_\mathbb{B}P(t, \phi, \varsigma, dx\times\{u\}).$ Denote by $P(R_+^2, \mathbb{S})$ all the probability measures on $R_+^2\times \mathbb{S}$, and for any two $P_1, P_2\in P(R_+^2, \mathbb{S})$, we define the metric $d_L$ as follows:

where $L = \{f:C([-\tau, 0], R_+^2)\times \mathbb{S}\rightarrow R:|f(x_1, k)-f(x_2, \tilde{k})\leq |x_1-x_2|+|k-\tilde{k}|, |f(\cdot, \cdot)|\leq 1\}.$

Definition 4.1 (^[9]). The process $x(t, \phi, \varsigma)$ is said to be asymptotically stable in distribution if there exists a probability measure $u(\cdot\times \cdot)$ on $R_+^2\times \mathbb{S}$ such that the transmission probability $P(t, \phi, \varsigma, dx\times\{k\})$ of $x(t, \phi, \varsigma)$ converges weakly to $u(dx\times\{k\})$ as $t\rightarrow \infty$ for every $(\phi, \varsigma)\in C([-\tau, 0], R_+^2)\times \mathbb{S}$. System (1.2) is said to be asymptotically stable in distribution if the solution $x(t, \phi, \varsigma)$ of (1.2) is asymptotically stable in distribution.

Lemma 4.1 (^[30]). Let $f(t)$ be a nonnegative function defined on $[0, \infty)$ such that $f(t)$ is integrable on $[0, \infty)$ and is uniformly continuous on $[0, \infty)$, then $\lim_{t\rightarrow \infty}f(t) = 0.$

For the need of discussion, we give the following technical assumption.

Assumption 3. $a_{ii}-\int_{-\tau_{ii}}^0d\mu_{ii}(\theta)-\int_{-\tau_{ji}}^0d\mu_{ji}(\theta) > a_{ji}, i, j = 1, 2, i\neq j.$

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, \phi, \varsigma)$ and $\widetilde{x}(t) = \widetilde{x}(t, \varphi, \varsigma)$ be two solutions of (1.2) withinitial value $x(\theta) = \phi, \widetilde{x}(\theta) = \varphi$ and $\alpha(\theta) = \varsigma$, respectively.If Assumption 3 holds, then we have

Proof. Define Lyapunov function $V_i(t) = |\ln x_i(t)-\ln \widetilde{x}_i(t)|, t\geq0, i = 1, 2.$ 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:

Let $V(t) = V_1(t)+V_2(t)+V_3(t)$, then by Assumption 3, after computation, we have

Therefore,

Hence,

for $i, j = 1, 2, i\neq j.$

On the other hand, by (1.2) we have,

where

Taking expectation from both sides of (4.2), then

Consequently, $\mathbb{E}(x_i(t))$ is continuously differentiable with respect to t. Moreover,

That is to say, $\mathbb{E}(x_i(t))$ is uniformly continuous. An application of Lemma 4.1 gives $\lim_{t\rightarrow \infty}\mathbb{E}|x_i(t)-\widetilde{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.

Lemma 4.2. For any compact subset $\mathbb{B}\subseteq R_+^2$ and $(\phi, \varsigma)\in C([-\tau, 0], \mathbb{B})\times \mathbb{S}$, the family of transmission probability of the solution$P(t, \phi, \varsigma, dx\times\{u\})$ is tight.

Proof. For (4.2), by use of the Hölder inequality and the moment inequality of stochastic integrals, there exist $k = 1, 2, ...$ such that

By Lemma 2.2 (i.e., $\mathbb{E}|x_i(t)|^p\leq K(p)$), then

where $M_1(p)$ is a constant. On the other hand, using the Burkholder-Davis-Gundy inequality ^[19], there exists a constant $C_p$ such that

As to the semimartingale part, by applying the Kunita's first inequality ^[19], there exists a constant $D(p)$ such that

Therefore, for any $s\in[0, t]$, we have

That is, the probability measure set $P(t, \phi, \varsigma, dx\times\{u\})$ is tight for $(\phi, \varsigma)\in C([-\tau, 0], R_+^2)\times \mathbb{S}$. This completes the proof.

Lemma 4.3. Let Assumptions 1–3 hold, then for any $(\phi, \varsigma)\in C([-\tau, 0], R_+^2)\times \mathbb{S}$, the transmission probability$P(t, \phi, \varsigma, \cdot\times\cdot:t\geq0)$ of the solution of (1.2) is Cauchy in the space $P(R_+^2\times \mathbb{S})$with the metric $d_L$ defined as before.

Proof. Let $\mathbb{B} = \{x\in C([-\tau, 0], R_+^2):\sigma\leq|x|\leq\varrho\},$ where $\varrho$ is a sufficiently large positive number, and $\sigma$ is a sufficiently small positive number. By the tightness of the transmission probability of solutions of (1.2), for any $\varepsilon > 0,$ we have $P (t, x, \alpha, \mathbb{B}^C\times \mathbb{S})\leq\varepsilon, t\geq0,$ where $\mathbb{B}^C = R_+^2/\mathbb{B}.$

For any $f\in L$, $t, s > 0$, and initial value $(\phi, \varsigma)\in C([-\tau, 0], R_+^2)\times\mathbb{ S}$, by the characteristic of condition expectation, we have

where

By the proof of Theorem 4.6 in ^[23] (4.22 of page 106), for any $\varepsilon > 0$ and sufficiently large $t$, we have

Then

Hence, $d_L(P(t+s, \phi, \varsigma, \cdot\times\cdot), P(t, \phi, \varsigma, \cdot\times\cdot))\leq3\varepsilon$ for $s > 0$ and sufficiently large $t$, that is, $P(t, \phi, \varsigma, \cdot\times\cdot:t\geq0)$ is Cauchy in the space $P(C([-\tau, 0], R_+^2)\times \mathbb{S})$. The proof is completed.

Theorem 4.2. Under the conditions of Lemma 4.3, the solution process $x(t)$ of (1.2) is asymptotically stable in distribution.

Proof.By Lemma 4.3, the transmission probability $P(t, \phi, \varsigma, \cdot\times\cdot:t\geq0)$ is Cauchy in the space $P(R_+^2\times \mathbb{S})$ with the metric $d_L$. Hence, for any fixed $\varphi\in C([-\tau, 0], R_+^2)$, $P(t, \widetilde{x}, {\varsigma}, \cdot\times\cdot:t\geq0)$ is Cauchy in $P(R_+^2\times \mathbb{S})$. Then there exists a probability measure $u(\cdot\times\cdot)$ such that

On the other hand, by Theorem 4.1, we have

Therefore, by the triangle inequality, we can derive from (4.10) and (4.11) that

Since the weak convergence of probability measure is a metric concept, (4.12) shows that for any initial value $(\phi, {\varsigma})\in C([-\tau, 0], R_+^2)\times \mathbb{S}$, the probability measure $P(t, \phi, {\varsigma}, \cdot\times\cdot:t\geq0)$ of the solution of (1.2) converges weakly to the probability measure $u(\cdot\times\cdot).$ By Definition 4.1, then the solution process $x(t)$ of (1.2) is asymptotically stable in distribution. This completes the proof.

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.

5.   Optimal harvesting effort

In this section, we consider the optimal harvesting effort (OHE) of (1.2). By ^[8,27,28,29], the OHE problem is to find a constant $h^* = (h_1^*, h_2^*)^T$ such that both $x_1$ and $x_2$ survive, and $Y(h^*) = \lim_{t\rightarrow \infty}\mathbb{E}(h_i^*x_i(t))$ is maximum.

Theorem 5.1. Let $\lambda = (\lambda_1, \lambda_2)^T = (A(A^{-1})^T+I)^{-1}\xi,$ where $I$ is the unit matrix, $\xi = (\bar{\xi}_1, \bar{\xi}_2)^T.$

(1) If $A+(A^{-1})^T$ is positive semi-definite, and $\Delta_2|_{h_i = \lambda_i} > 0, \lambda_i\geq0, i = 1, 2,$ then the OHE is $h^* = \lambda$, and the maximum of the sustainable yield is

(2) If the conditions of (1) fail, then there exists no OHE.

Proof. We define $\Xi = \{h = (h_1, h_2)^T\in R_+^2, \Delta_2 > 0\}$, then for any $h\in \Xi,$ by Theorem 3.1, we have $\lim\limits_{t\rightarrow \infty} \frac{1}{t}\int_0^tx_i(t)dt = \frac{\Delta_i}{\Delta}, i = 1, 2.$ Conversely, if the $h^*$ of OHE exists, obviously $h^*\in\Xi.$

(1) By the given conditions, it is clear that $\lambda\in\Xi,$ that is, the set $\Xi$ is not empty. For any $h\in\Xi,$ 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 $\rho,$ 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 $\rho$ and the invariant measure $u$ (see, e.g.^[31], page 105), we have

By (5.1)–(5.3), we observe that

Clearly, by computing the derivative of the variable $h$, we have $\frac{d(Y(h))}{dh} = A^{-1}\xi-(A^{-1}+(A^{-1})^T)h.$ Denote the solution of $\frac{d(Y(h))}{dh} = 0$ by $\lambda$, then $\lambda = (AA^{-1}+I)\xi.$ Further, after computation, 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 $\lambda$ is a global maximum point of $Y(h)$. That is, if the conditions $\Delta_2|_{h_i = \lambda_i} > 0, \lambda_i\geq0 (i = 1, 2)$ hold, then the OHE is $h^* = \lambda$, and $Y(h^*) = \lambda^TA^{-1}(\xi-\lambda).$

(2) Firstly we prove that there is no OHE if $\Delta_2|_{h_i = \lambda_i, i = 1, 2} > 0,$ $\lambda_1\geq0$ and $\lambda_2\geq0$, but $A^{-1}+(A^{-1})^T$ is not positive semi-definite. By above conditions, clearly the set $\Xi$ is not empty. Let $(\iota_{ij})_{2\times2} = A^{-1}+(A^{-1})^T$, then $\iota_{11} = \frac{2A_{22}}{\Delta} > 0.$ It implies that $A^{-1}+(A^{-1})^T$ is not negative semi-definite. Noting that $A^{-1}+(A^{-1})^T$ is not positive semi-definite, hence $A^{-1}+(A^{-1})^T$ is indefinite. The extreme value theory ^[32] shows that $Y(h)$ has no extreme points. Therefore, there is no OHE.

Secondly, under the conditions $\Delta_2|_{h_i = \lambda_i, i = 1, 2} < 0,$ or $\lambda_1\leq0$, or $\lambda_2\leq0$, we prove there is also no OHE. Otherwise, suppose the OHE is $\tilde{h}^* = (\tilde{h}_1^*, \tilde{h}_2^*)^T$, then $\tilde{h}^*\in\Xi.$ Thus $\Delta_2|_{h_i^* = \tilde{h}_i^*} > 0, \tilde{h}_i^*\geq0, i = 1, 2.$ That means $\tilde{h}_i^*$ is also the solution of $\frac{d(Y(h))}{dh} = 0,$ which contradicts with the uniqueness of the solution. This completes the proof.

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.

6.   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.

For simplicity, we assume that the continuous-time discrete state Markovian chain $\alpha(t)$ takes value in the space $\mathbb{S} = \{1, 2\}$, then system (1.2) reduces to the following subsystem (6.1) and (6.2):

and

For (6.1) and (6.2), unless otherwise stated, we always take $a_{11} = 0.8, a_{12} = 0.6, a_{21} = 0.4, a_{22} = 0.7$, $h_i = 0, \sigma_i^2 = 0.2, \int_{-\tau_{ii}}^0d\mu_{ii}(\theta) = 0.3,$ $\int_{-\tau_{ij}}^0d\mu_{ij}(\theta) = 0.2, \gamma_i(1) = \gamma_i(2) = 0.4,$ $i = 1, 2, i\neq j, \mathbb{Z} = [0, \infty),$ and $\lambda(\mathbb{Z}) = 1.$ Then $A_{11} = 1.1, A_{12} = 0.8, -A_{21} = -0.6, A_{22} = 1, \Delta = 1.58.$

For (6.1), let $r_1(1) = 0.06, r_2(1) = -0.02.$ It is easy to verify that Assumptions 1–3 hold, and ${\eta}_{1}(1) = -0.1035 < 0.$ Theorem 3.1 implies both $x_1$ and $x_2$ are extinct (see Figure 1 (a)).

For (6.2), let $r_1(2) = 0.7, r_2(2) = -0.05.$ Similarly Assumptions 1–3 hold, and $\eta_{1}(2) = 0.5365 > 0, \eta_{2}(2) = -0.2135,$ $\Delta_1 = 0.7073, \Delta_2 = 0.087 > 0.$

By Theorem 3.1 again, we can obtain that both $x_1$ and $x_2$ are persistent in the mean, and

Figure 1 (b) reveals it clearly.

Next we will study the effect of regime switching, time delays and Lévy jumps on the extinction and persistence in the mean of all species, respectively. We discuss in three cases.

● Case (ⅰ) The effect of switching $\pi$.

(1) Let the stationary distribution $\pi = (0.1, 0.9).$ By computation, then

and

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 $\alpha(t)$ with stationary distribution $\pi = (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 $\lambda = (0.3258, 0.0532)^T,$ the green line represents the yield $Y_1 = 0.0484$ with $\lambda_1 = (0.1558, 0.0632)^T,$ the blue line represents the yield $Y_2 = 0.0601$ with $\lambda_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.

(2) Let $\pi = (0.9, 0.1).$ Then $\bar{\eta}_1 = 0.9*(-0.1035)+0.1*0.5365 = -0.0395 < 0.$ Theorem 3.1 implies that both $x_1(t)$ and $x_2(t)$ are to be extinct. That is, the switching $\pi$ leads to the extinction of $x_1(t)$ and $x_2(t)$ (see Figure 4 (a)).

(3) Let $\pi = (0.4, 0.6).$ Then by computing, we have

By Theorem 3.1 again, then $x_2$ is extinct and $x_1$ is persistent in the mean, and

It shows the switching $\pi$ leads to the extinction of $x_2(t)$ and the different persistence in the mean of $x_1(t)$ (see Figure 4 (b)).

● Case (ⅱ) The effect of time-delays.

For the persistent case, if we take $\int_{-\tau_{ii}}^0d\mu_{ii}(\theta) = 0.87, \int_{-\tau_{ij}}^0d\mu_{ij}(\theta) = 0.3,$ other parameters are same as before, then $\bar{\eta}_1 = 0.4725 > 0, \Delta_2 = -0.0208 < 0,$ which leads to the extinction of $x_2(t)$ and

That is, time delays destroy the persistent case, see Figure 4 (c).

● Case (ⅲ) The effect of Lévy jumps.

For the persistent case, if we take $\gamma_i(1) = \gamma_i(2) = 1, i = 1, 2, i\neq j,$ then $\overline{{\eta}}_{1} = -0.4189 < 0.$ By Theorem 3.1, both $x_1$ and $x_2$ go to be extinct. It implies that too large Lévy jumps leads to the extinction of $x_1(t)$ and $x_2(t)$, see Figure 4 (d).

7.   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 $\alpha(t) = 1, h_i = 0, \gamma_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 $\mu_{ii}(\theta)$ are constant on $[-\tau, 0], a_{ij} = 0, i\neq 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 $i$th species at time $t$ is often affected by the abundance of the $j$th species on the interval $[t-\tau_i, t]$, rather than only on the time $t-\tau_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.

Acknowledgements

The author is very grateful to the anonymous referees and the editor for their helpful and suggestive comments. This research was supported by the National Natural Science Foundation of China (Nos. 11861027 and 11871475).

Conflict of interest

The authors declare that they have no competing interests in this paper.