Dynamical analysis of a stochastic hybrid predator-prey model with Beddington-DeAngelis functional response and Lévy jumps

1.   Introduction

In ecological mathematics, in order to simulate a variety of different situations (see, e.g., ^[1,2]), many different kinds of predator-prey models ^[3,4,5,6] are studied. In 1975, Beddington-DeAngelis ^[7,8] showed that the predator-dependent functional response can provide a better description for predator-prey models, and introduced the following Beddington-DeAngelis type predator-prey model:

where $x(t)$ denotes the population densities of prey at time $t$, and $y(t)$ represents the population densities of predator at time $t$. $r_1$ denotes the intrinsic growth rate of the prey, $r_2$ denotes the mortality rate of the predator, $a_{11}$ represents the density dependent coefficient of prey, $a_{12}$ represents the capturing rate of predator, $a_{21}$ represents the rate at which nutrients are converted into reproduction of the predator. In addition, $r_i, a_{ij}, m, n$ are all positive constants, $i, j = 1, 2$.

Species are inevitably disturbed by environmental noise (see, e.g., ^{[9,10,11,12,13,14,15]}). Therefore, it is of great significance to study the effect of noise on biological population system. A large number of authors have studied stochastic species models with random factors. Du et al. investigated dynamics of a stochastic Lotka-Volterra model perturbed by white noise ^[10]. Mao et al. got that even small enough noise can have an effect for explosion in species dynamics in ^[12]. From ^[14], we know that Wang et al. studied stationary distribution of a stochastic hybrid phytoplankton-zooplankton model with toxin-producing phytoplankton. Considering the influence of random environment, random disturbances should be introduced into the model to explore the effects of random disturbances. On the one hand, it is assumed that the random interference is white noise. Therefore, we obtain the following stochastic Beddington-DeAngelis type predator-prey model with white noise:

where $\sigma_{i}^{2}(i = 1, 2)$ stands for the intensity of the white noise, $\{B_{1}(t), B_{2}(t)\}_{t\geq 0}$ is a two dimensional Brownian motion. Throughout this paper, the Brownian motion was defined on a complete probability space $(\Omega, \{\mathcal{F}_{t}\}_{t\geq0}, P)$ with a filtration $\{\mathcal{F}_{t}\}_{t \in{R}_+}$ satisfying the usual conditions.

Because growth rates and death rates are especial sensitive to the disturbance of the environment. It is known that the rates often switch randomly because of changes in the environment. For instance, temperature or rain falls ^[16]. A classic example ^[17] is that some species grow at vastly different rates during the dry and rainy seasons. It has been suggested ^{[16,17,18,19]} that these random changes can be simulated using a continuous time finite state Markov chains. Hence, considering this influence, we get the following model:

where $\lambda(t)$ represents a continuous-time Markov chain with a state space $\mathbb{S} = \{1, \; 2, \; 3, \; \cdots, \; n\}$.

On the other hand, white noise is not a good description for some sudden environmental disturbances (such as earthquakes, volcanoes, hurricanes, debris flows, etc.) that are often encountered during the growth of species. Bao et al. ^[20,21] suggested that these phenomena can be described by Lévy jumps, and they studied the stochastic Lotka-Volterra population systems with Lévy jumps. In 2015, Zang et al. ^[22] studied dynamics of a stochastic predator-prey model with Beddington-DeAngelis functional response and Lévy jumps. To the best of our knowledge, no results related to Beddington-DeAngelis predator-prey model with regime-switching and Lévy jumps has been reported. Motivated by this, we will consider the following model:

where $x(t^{-})$ and $y(t^{-})$ are the left limit of $x(t)$ and $y(t)$, respectively. $\tilde{N}(\mathrm{d}t, \mathrm{d}u) \; = N(\mathrm{d}t, \mathrm{d}u)-\lambda(\mathrm{d}u)\mathrm{d}t$, on a measurable subset ${Z}$ of $(0, +\infty)$ with $\lambda({Z}) < +\infty$, $N$ is a Poisson counting measure with characteristic measure $\lambda$, $\gamma_i:{Z}\times \Omega\rightarrow {R}$ is bounded and continuous with respect to $\lambda$, $i = 1, 2$. Let's make the following assumption

To better simulate the population, a predat-prey model with Bedtington-DeAngelis functional response is studied in this paper. In our model, predator not only has time to find the prey, but also has time to deal with the encounter with other predators. We also consider the influence of environmental factors. The value of this paper is to study the stochastic predator-prey model of Bedtington-DeAngelis function response which is affected by both regime-switching and Lévy jumps. Furthermore, we introduce some numerical simulations to illustrate the main results.

2.   Dynamical behaviors and properties

For the sake of convenience, we define some notations.

Assumption 2.1. For any $j\in\mathbb{S}$, assume that there exists a constant $G_1 > 0$ such that

That is to say, the intensity of Lévy noise is not too large.

Assumption 2.2. For any $j\in\mathbb{S}$, assume that there exist two constants $G_2, G_3 > 0$ such that

Lemma 2.1. Assumptions 2.1 and 2.2 hold, then the model (1.4) with initial value $(x(0), y(0), \lambda(0))$ has a unique positive local solution $(x(t), y(t), \lambda(t))$ on $t\in[0, \tau_e)$ almost surely(a.s.), where $\tau_e$ is the explosion time.

Proof. Consider the following equations:

with initial value $(\varphi(0) = \ln x(0), \psi(0) = \ln y(0), \lambda(0))$, where

According to the Lemma 1 of reference ^[23], we know that the coefficients of model (2.1) are local Lipschitz continuous, hence there is a unique local solution $(\varphi(t), \psi(t), \lambda(t))$ on $t\in[0, \tau_e)$. That is to say, $(x(t) = e^{\varphi(t)}, y(t) = e^{\psi(t)}, \lambda(t))$ is the unique positive local solution of (1.4) with initial value $(x(0), y(0), \lambda(0))$ by $\rm It\hat{o}$'s formula.

Theorem 2.1. Consider the model (1.4), for given initial value$(X(0), \lambda(0))\in{R}_+^{2}\times\mathbb{S}$, where $X(0) = (x(0), y(0))$, there is a unique global solution $(x(t), y(t), \lambda(t))$ on $t\geq0$, and the solution will remain in ${R}_+^{2}$ with probability 1.

Proof. By Lemma 2.1, we just need to prove $\tau_e = +\infty$. Consider the following equations:

By the comparison theorem of stochastic equation (see e.g., ^[24] (Theorem 1 on page 173) or ^[25] (Theorem 3.1)), we get

According to Lemma 3.1 of reference ^[26], (2.2)–(2.4) have the following explicit solutions:

where

Because $M_1(t), M_2(t)$ and $M_3(t)$ are existent on $t\geq0$, then we obtain $\tau_e = +\infty$.

Lemma 2.2. (Lemma 8 in ^[27]) Suppose that $x(t)\in C(\Omega\times[0, +\infty), {R}_+)$, and $\lim_{t\rightarrow +\infty}\frac{F(t)}{t} = 0$ a.s.

(1) If there exist two constants $T > 0$ and $\rho_0 > 0$ such that for all $t\geq T$,

then

(2) If there exist three constants $T > 0$, $\rho > 0$ and $\rho_0 > 0$ such that for all $t\geq T$,

then $\liminf_{t\rightarrow +\infty}t^{-1}\int_0^{t}x(s)\mathrm{d}s\geq \frac{\rho}{\rho_0}$ a.s.

Lemma 2.3. The solution of system (1.4) have the following properties:

Proof. From system (1.4), we get

Set

where $(\bar{x}(t), \bar{y}(t), \lambda(t))$ is a solution of model (2.5) with initial value $(x(0), y(0), \lambda(0))$. According to the comparison theorem of stochastic differential equations, it's easy to figure out $x(t)\leq \bar{x}(t)$ for $t\geq0$.

By using Lemma 2 of reference ^[28], one obtain the following results for $t\geq0$: $\limsup_{t\rightarrow +\infty}\frac{\ln {\bar{x}(t)}}{\ln t}\leq1$, $\limsup_{t\rightarrow +\infty}\frac{\ln {\bar{y}(t)}}{\ln t}\leq1, \; a.s.,$ therefore, we get $\limsup_{t\rightarrow +\infty}\frac{\ln x(t)}{\ln t}\leq1$. That is to say, we have $\limsup_{t\rightarrow +\infty}\frac{\ln x(t)}{t} = \limsup_{t\rightarrow +\infty}\frac{\ln x(t)}{\ln t}\limsup_{t\rightarrow +\infty}\frac{\ln t}{t}\leq\limsup_{t\rightarrow +\infty}\frac{\ln t}{ t} = 0$. Similarly, we have $\limsup_{t\rightarrow +\infty}\frac{\ln y(t)}{t}\leq0$.

Definition 2.1. ^[29] The population $x(t)$ is weakly persistent in the mean if $\langle x\rangle^{*} > 0$, where $\langle x(t)\rangle = \frac{1}{t}\int_0^{t}x(s)\mathrm{d}s$, $(x(t))^{*} = \limsup_{t\rightarrow +\infty}x(t)$.

Theorem 2.2. Assumptions 2.1 and 2.2 hold.

(i) The prey $x(t)$ will go to extinction if $\bar{b}_1 < 0$, i.e., $\lim_{t\rightarrow +\infty}x(t) = 0\; a.s.$.

(ii) The prey $x(t)$ will be weakly persistent in the mean if $\bar{b}_1 > 0$, i.e., $\langle x(t)\rangle^{*} > 0\; a.s.$.

Proof. (i) Applying $\rm It\hat{o}$'s formula to system (1.4), we obtain

where

By the ergodicity of the Markov chain, we have

Hence we have $\bar{b}_2 < 0$. From the first equation of (2.6), we get

so we have

According to Assumption 2.1 and the strong law of large numbers ^[30], we can see that

taking superior limit for both sides of (2.8), one obtain $\limsup_{t\rightarrow +\infty}\frac{\ln x(t)}{t}\leq\bar{b}_1 < 0$. Namely, $\lim_{t\rightarrow +\infty}x(t) = 0$.

(ii) We just need to prove that there is a constant $u_1 > 0$ that makes $\langle x(t)\rangle^{*} = u_1 > 0\; a.s.$ true for the solution $(x(t), y(t), \lambda(t))$ of system (1.4) that has the initial value $(x(0), y(0), \lambda(0))$. If not, for any $\varepsilon_1 > 0$, there exists a solution $(x_1(t), y_1(t), \lambda(t))$ with an initial value $(x(0), y(0), \lambda(0))$ such that $\mathrm{P}\{\langle x_1(t)\rangle^{*} < \varepsilon_1\} > 0$ holds. Suppose $\varepsilon_1$ is small enough, and satisfy the following conditions:

According to the second equation of (2.6), we have

taking superior limit on both sides of (2.10), using virtue of (2.7), (2.9) and (2.10), we have

So, we obtain

By using the first equation of (2.6), we have

therefore, we obtain

By the same method, using virtue of (2.7), (2.9), (2.11) and (2.12), one obtain

Namely, one have $\mathrm{P}\{(t^{-1} \ln x_1(t))^{*} > 0\} > 0$, this is a contradiction to Lemma 2.3. So $\langle x(t)\rangle^{*} > 0$, $x(t)$ will be weakly persistent in the mean a.s.. Complete the proof.

Theorem 2.3. Assumptions 2.1 and 2.2 hold.

(i) The predator $y(t)$ will go to extinction if $a_{11}\bar{b}_2+a_{21}\bar{b}_1 < 0$, i.e., $\lim_{t\rightarrow +\infty}y(t) = 0\; a.s.$.

(ii) The predator $y(t)$ will be weakly persistent in the mean if $\bar{b}_2+\langle\frac{a_{21}\bar{x}(t)}{1+m\bar{x}(t)+n\bar{y}(t)}\rangle^{*} > 0$, i.e., $\langle y(t)\rangle^{*} > 0\; a.s.$, where $(\bar{x}, \bar{y})$ is a solution of system (2.5) with initial value $(x(0), y(0), \lambda(0))$.

Proof. (i) If $\bar{b}_1\leq0$, from the Theorem 2.2, we can know that $\lim_{t\rightarrow \infty}x(t) = 0$. According to (2.10), we have

So we obtain $\limsup_{t\rightarrow +\infty}t^{-1}\ln y(t)\leq\bar{b}_2 < 0$, then $\lim_{t\rightarrow +\infty}y(t) = 0$.

If $\bar{b}_1 > 0$, according to the property of limit and (2.9), there exists a $T$ for sufficiently small $\varepsilon$ such that

on $t > T$. Using Lemma 2.2 and the arbitrariness of $\varepsilon$, one can obtain $\langle x(t)\rangle^{*}\leq\frac{\bar{b}_1}{ a_{11}}$. Substituting the inequality into the second equation of (2.6). By further calculation, we have

Namely, we get $\lim_{t\rightarrow +\infty}y(t) = 0$ a.s..

(ii) We just need to prove that there is a constant $u_2 > 0$ that makes $\langle y(t)\rangle^{*} = u_2 > 0\; a.s.$ true for the solution $(x(t), y(t), \lambda(t))$ of system (1.4) that has the initial value $(x(0), y(0), \lambda(0))$. If not, for any $\varepsilon_2 > 0$, there exists a solution $(x_2(t), y_2(t), \lambda(t))$ with an initial value $(x(0), y(0), \lambda(0))$ such that $\mathrm{P}\{\langle y_2(t)\rangle^{*} < \varepsilon_2\} > 0$ holds. Suppose $\varepsilon_2$ is small enough, and satisfy the following conditions:

From the second equation of (2.6), one can get

where $(\bar{x}, \bar{y})$ is a solution of system (2.5), we have $x_{2}(t)\leq\bar{x}(t)$, $y_{2}(t)\leq\bar{y}(t), \; a.s.$ on $t > 0$.

Due to

thus we have

Next, we define the function $V_1(t): = |\ln \bar{x}(t)-\ln x_2(t)|$, by using $\rm It\hat{o}$'s formula, we obtain

integrating and dividing by $t$ on both sides of the above inequality, one obtain $\frac{V_1(t)-V_1(0)}{t}\leq a_{12}\langle y_2(t)\rangle-a_{11}\langle\bar{x}(t)-x_2(t)\rangle$. Since $\frac{V_1(t)}{t}\geq0$, we have $a_{11}\langle \bar{x}(t)-x_2(t)\rangle\leq a_{12}\langle y_2(t)\rangle+\frac{V_1(0)}{t}$, Note that $V_1(0) = 0$, thus $\langle \bar{x}(t)-x_2(t)\rangle\leq \frac{a_{12}}{a_{11}}\langle y_2(t)\rangle$.

Substituting into (2.14), we get

That is to say, we have

which contradicts Lemma 2.3, so we get $\langle y(t)\rangle^{*} > 0$ a.s., i.e., $y(t)$ is weakly persistent in the mean.

Lemma 2.4. (^[31] Theorem $2.3$) Let Assumptions 2.1 and 2.2 hold. For any $q > 0, \; t\geq0$, there exist constants $K_i(q) > 0, \; i = 1, \; 2$ such that

Assumption 2.3. $a_{11}\geq a_{12}\frac{m}{n}+2a_{21}$, $a_{22}\geq a_{21}\frac{n}{m}+2a_{12}$.

Lemma 2.5. Let $X\big((x_0, y_0, j), t\big) = \big(x((x_0, y_0, j), t), \; y((x_0, y_0, j), t)\big)$ is a solution of model (1.4) with initial data $((x_0, y_0), j)\in D\times \mathbb{S}$, $X\big((\tilde{x}_0, \tilde{y}_0, \tilde{j}), t\big) = \big(x((\tilde{x}_0, \tilde{y}_0, \tilde{j}), t), \; \; y((\tilde{x}_0, \tilde{y}_0, \tilde{j}), t)\big)$ is a solution of model (1.4) with initial data $((\tilde{x}_0, \tilde{y}_0), \tilde{j})\in D\times \mathbb{S}$ respectively, where $D$ is an arbitrary compact subset of ${R}_+^{2}$. If Assumptions 2.1–2.3 hold, then

Proof. For the sake of convenience, we define the following notations:

Define $V_2(t) = \big|\ln x-\ln \tilde{x}\big|+\big|\ln y-\ln \tilde{y}\big|$, then we have

therefore, we get

Note that $EV_2(t)\geq 0$, so we have

Hence $E|x-\tilde{x}|\in L^{1}[0, +\infty), \; E|y-\tilde{y}|\in L^{1}[0, +\infty)$. By Barb$\check{\rm a}$ lat's conclusions (see, e.g., ^[32]), now we just proof that $E(x(t))$ and $E(y(t))$ are uniformly continuous with respect to $t$. As a matter of fact, thanks to system (1.4), one get

So, $E(x(t))$ and $E(y(t))$ are continuously differentiable. In addition, according to above formulas and Lemma 2.4, we have

where $K_1 > 0, \; K_2 > 0$ are two constants, $r_1^{u} = \max_{i\in\mathbb{S}}r_1(i)$. Namely, $E(x(t))$ and $E(y(t))$ are uniformly continuous.

Theorem 2.4. Suppose that Assumptions 2.1–2.3 hold, if $\bar{b}_1 > 0$, $\bar{b}_2+\langle\frac{a_{21}\bar{x}(t)}{1+m\bar{x}(t)+n\bar{y}(t)}\rangle^{*} > 0$, then system (1.4) has a unique stationary measure $\eta(\cdot\times\cdot)$ which is ergodic.

To prove this theorem we need to introduce more notations. Let ${B}({R}_+^{2}\times\mathbb{S})$ represent all the probability measures defined on ${R}_+^{2}\times\mathbb{S}$. For any two measures $p_1, \; p_2\in{B}$, define the metric $d_H$ as follows

where

Let us now present the following lemmas.

Lemma 2.6. For every $q > 0$ and any compact subset $D$ of ${R}_+^{2}$, $\sup_{(X(0), j)\in D\times\mathbb{S}}E\big[\sup _{0\leq s\leq t}\big|X^{X(0), j}(s)\big|^{q}\big] < +\infty, \; \forall t\geq0$, where $X^{X(0), j} = (x^{X(0), j}, y^{X(0), j})$.

Proof. From (1.4), we have

According to the H$\ddot{\rm o}$lder inequality and the moment inequality of stochastic integrals, exists $k = 1, 2, \cdot\cdot\cdot$, such that

By Lemma $2.4$, exist a positive constant $K(q)$ such that $E\big|x(t)\big|^{q}\leq K(q)$, $t\in[0, +\infty)$. Thus we can obtain that

In view of the Burkholder-Davis-Gundy inequality (see Theorem 1.7.3 of reference ^[33]), we have

where $\sigma_1^{u} = \max_{i\in\mathbb{S}}\sigma_1(i)$. Make use of Assumption 2.1 and Kunita's first inequality (see Theorem 4.4.23 in reference ^[34]), one obtain

According to (2.15)–(2.18), we have

Similarly, we obtain

Therefore, we have

Lemma 2.7. Assumption 2.1–2.3 hold, then for any compact subset $D$ of ${R}_+^{2}$, we have

uniformly in $X(0), \tilde{X}(0)\in D$ and $i, j\in \mathbb{S}$.

Proof. The proof is similar to Lemma $5.6$ of reference ^[13]. By using Lemmas 2.5 and 2.6, the proof is easy to prove, hence it is omitted.

Lemma 2.8. Assumption 2.1–2.3 hold. Then for any $(X(0), i)\in{R}_+^{2}\times\mathbb{S}$, $\{p(t, X(0), i, \; \cdot\times\cdot)|t\geq 0\}$ is Cauchy in the space ${B}({R}_+^{2}\times\mathbb{S})$ with metric $d_H$.

Proof. Fix any $(X(0), i)\in{R}_+^{2}\times\mathbb{S}$, we just need to proof that for any $\varepsilon_3 > 0$, there is a $T > 0$ such that

This is equivalent to

For any $h\in H$ and $t, s > 0$, we have

where $\bar{D}_R = \{X\in {R}_+^{2}\big||X|\leq R\}$, $\bar{D}_R^{C} = {R}_+^{2}-\bar{D}_R$. According to the well-known Chebyshev inequality, the family of transition probabilities $\{p(t, X(0), i, \mathrm{d}z_0\times\{l\}|t\geq 0)\}$ is tight. That is to say, for any $\varepsilon > 0$ there is a compact subset $D = D(\varepsilon, X(0), i)$ of ${R}_+^{2}$ such that $p(t, X(0), i, D\times \mathbb{S})\geq 1-\varepsilon, \; \forall t\geq 0$, there is a positive number $R$ sufficiently large for

Notes that by Lemma 2.7, there is a $T > 0$ such that

Substituting (2.22) and (2.23) into (2.21), we get

Because $h$ is arbitrary, the inequality (2.20) must hold.

Proof of Theorem 2.4. First, we prove that there is a probability measure $\eta(\cdot\times\cdot)\in{B}$ such that for any $(X(0), j)\in{R}_+^{2}\times\mathbb{S}$, the transition probability $p(t, X(0), j, \cdot\times\cdot)$ of $X((X(0), j), t)$ converges weakly to $\eta(\cdot\times\cdot)$. Based on proposition 2.5 in reference ^[35], the weak convergence of probability measures is a metric concept, namely, $p(t, X(0), j, \cdot\times\cdot)$ converges weakly to $\eta(\cdot\times\cdot) \; \Leftrightarrow$ there is a metric $d$ such that $\lim_{t\rightarrow +\infty}d(p(t, X(0), j, \cdot\times\cdot), \eta(\cdot\times\cdot)) = 0$. Hence, we just need to show that for any $(X(0), i)\in {R}_+^{2}\times\mathbb{S}$,

By using Lemma 2.8, $\{p(t, 0, 1, \cdot\times\cdot)|t\geq 0\}$ is Cauchy in the space $B({R}_+^{2}\times\mathbb{S})$ with metric $d_H$. Therefore, there is a unique $\eta(\cdot\times\cdot)\in B$ such that

By Lemma 2.7, we get

Namely, the distribution of $(X(t), \lambda(t))$ converges weakly to $\eta$. Due to the Kolmogorov-Chapman equation, we know that $\eta$ is invariant. Applying the Corollary 3.43 in reference ^[36], we obtain that $\eta$ is strong mixing. Therefore, $\eta$ is ergodic by Theorem 3.2.6 and (3.3.2) in reference ^[36].

Corollary 2.1. Consider model (1.4), Assumptions 2.1–2.3 hold, then

(i) if $\bar{b}_1 > 0$, $\bar{b}_2+\langle\frac{a_{21}\bar{x}(t)}{1+m\bar{x}(t)+n\bar{y}(t)}\rangle^{*} > 0$, then model (1.4) has a unique stationary distribution $\eta(\cdot\times\cdot)$.

(ii) if $\bar{b}_1 > 0$, $a_{11}\bar{b}_2+a_{21}\bar{b}_1 < 0$, then $x(t)$ has a unique ergodic stationary distribution, $y(t)$ goes to extinction.

(iii) if $\bar{b}_1 < 0$, $a_{11}\bar{b}_2+a_{21}\bar{b}_1 < 0$, then both $x(t)$ and $y(t)$ will go to extinction.

3.   Numerical simulations

In this section, we will validate our theoretical results with the help of numerical simulation examples taking parameters. Our results show that the existence of stationary distribution has close relations with random disturbances. Let's consider the model (1.4) with $\mathbb{S} = \{1, 2\}$ to make it easier to understand. Therefore, according to the law of the Markov chain, system (1.4) can be regarded as a hybrid system which switches between the following two subsystems:

(i) Firstly, let us consider the effects about the distribution of Markov chain. We know the following facts: if the above two subsystems have a stationary distributions, then the hybrid system (1.4) still has a stationary distribution due to regime switching; if one of the two subsystems has a stationary distribution and the other does not, then the hybrid system (1.4) may have a stationary distribution, or not. To see the above more clearly, let's use several simulations to illustrate the impacts. Here, we only present the second case by letting the distribution of the Markov chain change(i.e., let $\pi$ change).

Example 1. Choose $a_{11} = 0.75, \; a_{12} = 0.25, \; a_{21} = 0.15, \; a_{22} = 0.95, \; {Z} = (0, +\infty), \; m = 0.15, \; n = 0.2, \; \lambda({Z}) = 1$. We have $a_{11}-a_{12}\frac{m}{n}-2a_{21} = 0.26\geq 0$, $a_{22}-a_{21}\frac{n}{m}-2a_{12} = 0.25\geq 0$.

In regime 1, choose $r_1(1) = 0.08, \; r_2(1) = 0.20, \; \sigma_1(1) = 0.53, \; \gamma_1(1) = 0.15, \; \sigma_2(1) = 0.59, \; \gamma_2(1) = 0.10\;$, thus $b_1(1) = -0.07 < 0, \; b_2(1) = -0.38 < 0, \; a_{11}b_2(1)+a_{21}b_1(1) = -0.30 < 0$. Therefore, according to Corollary 2.1, in subsystem (3.1), prey $x(t)$ and predator $y(t)$ go to extinction, see Figure 1.

In regime 2, choose $r_1(2) = 0.45, \; r_2(2) = 0.01, \; \sigma_1(2) = 0.08, \; \gamma_1(2) = 0.09, \; \sigma_2(2) = 0.11, \; \gamma_2(2) = 0.11$, therefore $b_1(2) = 0.44 > 0, \; b_2(2) = -0.02 < 0, \; a_{11}b_2(1)+a_{21}b_1(1) = 0.05 > 0$, but $\bar{b}_2+\langle\frac{a_{21}\bar{x}(t)}{1+m\bar{x}(t)+n\bar{y}(t)}\rangle^{*} > 0$ is difficult to verify. In subsystem (3.2), Figure 2 illustrates the situation that the model (1.4) has a unique stationary distribution $\eta(\cdot\times\cdot)$.

And then we're going to choose different $\pi$.

Case (a): We choose $\pi = (0.05\; 0.95)$. Therefore $\bar{b}_1 = 0.05\times(-0.07)+0.95\times0.44 = 0.41 > 0$; $\bar{b}_2 = 0.05\times(-0.38)+0.95\times(-0.02) = -0.04 < 0$, $a_{11}\bar{b}_2+a_{21}\bar{b}_1 = 0.75\times(-0.04)+0.15\times0.41 = 0.03 > 0$, but $\bar{b}_2+\langle\frac{a_{21}\bar{x}(t)}{1+m\bar{x}(t)+n\bar{y}(t)}\rangle^{*} \; > 0$ is difficult to verify. Figure 3 illustrates that the model (1.4) has a unique stationary distribution $\eta(\cdot\times\cdot)$.

Case (b): We choose $\pi = (0.3\; 0.7)$. Therefore $\bar{b}_1 = 0.3\times(-0.07)+0.7\times0.44 = 0.29 > 0$; $\bar{b}_2 = 0.3\times(-0.38)+0.7\times(-0.02) = -0.13 < 0$, and $a_{11}\bar{b}_2+a_{21}\bar{b}_1 = 0.75\times(-0.13)+0.15\times0.29 = -0.05 < 0$. Thus according to Corollary 2.1, the prey $x(t)$ has a unique ergodic stationary distribution, $y(t)$ goes to extinction, see Figure 4.

Case (c): We choose $\pi = (0.9\; 0.1)$. Therefore $\bar{b}_1 = 0.9\times(-0.07)+0.1\times0.44 = -0.02 < 0$; $\bar{b}_2 = 0.9\times(-0.38)+0.1\times(-0.02) = -0.34 < 0$, and $a_{11}\bar{b}_2+a_{21}\bar{b}_1 = 0.75\times(-0.34)+0.15\times(-0.02) = -0.26 < 0$. Therefore, according to Corollary 2.1, prey $x(t)$ and predator $y(t)$ go to extinction, see Figure 5.

From the above examples, we can know the following facts: one subsystem has a stationary distribution while the other doesn't have, different things can happen when we choose different $\pi$. When $\pi = (0.9\; 0.1)$, both prey $x(t)$ and predator $y(t)$ go to extinction; when $\pi = (0.3\; 0.7)$, the prey $x(t)$ has a unique ergodic stationary distribution, $y(t)$ goes to extinction; when $\pi = (0.05\; 0.95)$, the model (1.4) has a unique ergodic stationary distribution.

(ii) Next, we consider the effect of white noise on the population. For convenience, we will consider only subsystem (3.1). One can easily calculate the following facts

In short, subsystem (3.1) has a stationary distribution if for some $\sigma_1^{2}(1)$ and $\sigma_2^{2}(1)$, then the stationary distribution could disappear with the increase of $\sigma_1^{2}(1)$ or $\sigma_2^{2}(1)$. To see this more clearly, let's look at the following example.

Example 2. Consider the subsystem (3.1), the values of $\sigma_1(1)$ and $\sigma_2(1)$ are given in the table below, and the values of other parameters are the same with Example 1.

Therefore, in subsystem (3.1), according to Corollary 2.1, the prey $x(t)$ has a unique ergodic stationary distribution with $\sigma_1(1) = 0.04$, $x(t)$ is extinct at all other values. $y(t)$ goes to extinction due to $a_{11}b_2+a_{21}b_1 < 0$, see Figure 6. We know that when $\sigma_1(1)$ and $\sigma_2(1)$ increase, the stationary distribution of prey $x(t)$ could disappear, and the greater the intensity of random perturbation, the faster the species dies out.

(iii) Finally, consider the Lévy jumps. To keep things simple, we let $\gamma_i(u) = \varepsilon_i(i = 1, 2)$. We have

we know that $\varepsilon_i-\ln(1+\varepsilon_i)\geq0$, $\varepsilon_i\geq-1, i = 1, 2$. Therefore, the effect of $\varepsilon_i$ on the population is similar to that of $\sigma_i^{2}$, so it is omitted here.

After the above numerical simulation, there are similar examples in reality. In ^[37], the authors show that global warming has a profound bottom-up impact upon marine top-predators. Marine pollution via heavy metals, organochlorides, oil products and plastics is a recurrent threat to seabirds on a worldwide scale. All threats mentioned above cause substantial disturbance to seabird populations. Most seabirds feed on fish, and fish stocks are overexploited by industrial fishing, which can lead to large numbers of seabirds going hungry, all threats that cause significant disruption to bird populations. In ^[37], lesser sandeels Ammodytes marinus, which used to be the food-base of a vast seabird community around the British Isles, have been depleted due to the combined effects of overfishing and climate change. Seabirds, in particular kittiwakes Rissa tridactyla, now feed increasingly on snake pipefish. We know that not all seabirds are geographically malleable ^[38]. Seabirds may face extinction as some endemic diseases are trapped in restricted areas due to the effects of climate change. This is most likely the case for the Gal$\acute{\rm a}$pagos penguin Spheniscus mendiculus ^[39] and marblefinch bracketail fern ^[40]. Above examples show that under the influence of environmental disturbance caused by climate warming and Marine pollution, on the one hand, the population will die out under some certain conditions; on the other hand, under some certain conditions, the predator changes its foraging ecology so that it can survive, which indicates that the predator and the new prey will not die out.

4.   Conclusions

For the predator-prey model with Beddington-DeAngelis functional response, it has important theoretical and practical significance in real life, and has received extensive attention. However, a stochastic hybrid predator-prey model with Beddington-DeAngelis functional response and Lévy jumps has not been developed. In our work, our method easy to understand, but the sufficient conditions that make the conclusion hold may be a little strict, and the weakening of the conditions needs further study. Our main result is Theorem 2.4, which establishes sufficient conditions for the existence and uniqueness of an ergodic stationary distribution. Corollary 2.1 indicates that the existence of stationary distribution and extinction of model (1.4) depends on the signs of $\bar{b}_1$, $a_{11}\bar{b}_2+a_{21}\bar{b}_1$ and $\bar{b}_2+\langle\frac{a_{21}\bar{x}(t)}{1+m\bar{x}(t)+n\bar{y}(t)}\rangle^{*}$. Our results show that the existence of stationary distribution has close relations with the disturbance of environment.

There are still some interesting questions that deserve further consideration. For example, maybe one can give the threshold between the weak persistence in the mean and extinction for predator in future, further weakening the conditions in Theorem 2.4. Another interesting question is to consider what happens if other parameters are perturbed by noises. It is also interesting to consider other population systems, such as, food chain models, competition model and so on.

Acknowledgments

We would like to thank the the Fundamental Research Funds for the Central Universities (Grant number 3122021073).

Conflict of interest

All authors declare no conflicts of interest in this paper.