Open AccessArticle

Dynamics of a Prey–Predator Model with Group Defense for Prey, Cooperative Hunting for Predator, and Lévy Jump

Hengfei Chen

¹,

Ming Liu

^1,*

and

Xiaofeng Xu

^2,3,4

Department of Mathematics, Northeast Forestry University, Harbin 150040, China

Engineering Research Center of Agricultural Microbiology Technology, Ministry of Education, Harbin 150080, China

Heilongjiang Provincial Key Laboratory of Ecological Restoration and Resource Utilization for Cold Region, Harbin 150080, China

⁴

School of Mathematical Science, Heilongjiang University, Harbin 150080, China

Author to whom correspondence should be addressed.

Axioms 2023, 12(9), 878; https://doi.org/10.3390/axioms12090878

Submission received: 29 July 2023 / Revised: 24 August 2023 / Accepted: 12 September 2023 / Published: 14 September 2023

(This article belongs to the Special Issue Recent Advances in Applied Mathematics and Artificial Intelligence)

Download

Browse Figures

Versions Notes

Abstract

A stochastic predator–prey system with group cooperative behavior, white noise, and Lévy noise is considered. In group cooperation, we introduce the Holling IV interaction term to reflect group defense of prey, and cooperative hunting to reflect group attack of predator. Firstly, it is proved that the system has a globally unique positive solution. Secondly, we obtain the conditions of persistence and extinction of the system in the sense of time average. Under the condition that the environment does not change dramatically, the intensity of cooperative hunting and group defense needs to meet certain conditions to make both predators and preys persist. In addition, considering the system without Lévy jump, it is proved that the system has a stationary distribution. Finally, the validity of the theoretical results is verified by numerical simulation.

Keywords:

predator–prey; Holling IV; Lévy jump; cooperative hunting

MSC:

65C99; 92D25

1. Introduction

The predator–prey model is one of the most interesting and important topics in mathematical ecology, and has attracted many mathematicians and ecologists to study it. In the mid-20th century, Leslie and Gower [1,2] proposed a Leslie–Gower-type predator–prey model, which is characterized by the decrease in the number of predators being inversely proportional to its per capita preference for food availability. Considering that the number of predators was limited by the most important food, AzizAlaoui and Okiye [3] proposed the modified Leslie–Gower-type model. As far as we know, there are many articles on the Leslie–Gower model, most of which are deterministic equations. There are not many studies on Leslie–Gower with white noise and Lévy noise, and almost no research on the Leslie–Gower model with cooperative hunting and group defense functions.

In the predator–prey model, the predation rate intuitively reflects the relationship between the two populations. Scholars have proposed a series of functional responses to describe the predation rate [4,5,6,7,8]. Ali et al. [9] numerically studied the effects of nonlinear reaction–diffusion equations on the dynamics of prey–predator interactions. In nature, cooperation between species of the same species is common, and many scholars have introduced cooperation items into the modeling of functional responses for predation rates. Among them, there is cooperative hunting, such as lions hunting faster animals [10,11] and wolves hunting collectively against animals larger than themselves [12], and many scholars have studied the significance of hunting cooperation [13,14]. Chow and Jang [15] studied a predator–prey system with

\frac{a y}{1 + a y}

as the cooperative hunting term, where a is the coefficient related to cooperation strength. The results show that large-scale cooperative hunting may promote the persistence of predators. If there is no cooperation, predators will become extinct. Of course, species at the bottom of the food chain usually carry out group defense through mass reproduction. Andrews [16] proposed the so-called Monod–Haldane function and used it to model the inhibitory effect at high concentrations. Then, Sokol and Howell [17] proposed a simplified Monod–Haldane function of the form

\frac{K_{1} S}{K_{2} + S^{2}}

. After that, Shen [18] used a simple Holling IV to describe the group defense of prey for research. Bai [19] studied the global dynamics of a predator–prey system with cooperative hunting, and found that cooperative hunting is beneficial to the coexistence of prey and predator. Du [20] considered the cooperative hunting and group defense of the system without random perturbations, and the effect of predator cooperative hunting and predation aggregation on the stability of coexistence and system dynamics. Recently, many scholars have studied the models with group defense [21,22,23,24].

Inspired by the above articles, we add the Holling IV type function and cooperative hunting term into the Leslie–Gower model to show the cooperation of prey and predator.

\{\begin{matrix} d x (t) = & x (t) [a - b x (t) - \frac{q y^{2} (t)}{(m + y (t)) (c + x^{2} (t))}] d t, \\ d y (t) = & y (t) [δ - \frac{h y (t)}{l + x (t)}] d t, \end{matrix}

(1)

with initial values

x (0) = x_{0}

and

y (0) = y_{0}

x (t)

and

y (t)

denote the number of prey and predators at time t. It should be noted that

a, δ, b, q, m, c, h

, and l are all positive numbers. The parameters a and

δ

are the growth rate of the prey and the predator, respectively. b represents the intensity of competition between individuals of the prey. l describes the degree of protection provided by the environment to the predator. h denotes the maximum per capita reduction rate of the prey.

\frac{x (t)}{c + x^{2} (t)}

is the simplified Holling IV functional response function, where c describes the degree of group defense.

\frac{q y (t)}{m + y (t)}

is the cooperative hunting term, where q and m are cooperative coefficients, reflecting the intensity of cooperative hunting between predators.

What we know is that there are endless noise disturbances in nature. In the past, random biological models [25,26,27,28] have been a hot topic for biologists. In the study of stochastic interference, there are both white noise [29,30] and Lévy jump [31,32] considering sudden environmental disturbance. Random perturbations described by Brownian motion describe continuous effects, and Lévy jumps are considered to describe sudden and violent environmental changes well. Therefore, considering both white noise and Lévy noise is more realistic, so we obtain the following stochastic model:

\{\begin{matrix} d x (t) = & x (t) [a - b x (t) - \frac{q y^{2} (t)}{(m + y (t)) (c + x^{2} (t))}] d t \\ + σ_{1} x (t) d B_{1} (t) + \int_{Γ} x (t^{-}) γ_{1} (u) \tilde{N} (d u, d t), \\ d y (t) = & y (t) [δ - \frac{h y (t)}{l + x (t)}] d t + σ_{2} y (t) d B_{2} (t) \\ + \int_{Γ} y (t^{-}) γ_{2} (u) \tilde{N} (d u, d t), \end{matrix}

(2)

where

B_{1} (t)

and

B_{2} (t)

are mutually independent standard Brownian motions, and

x (t^{-}),

y (t^{-})

denote the left limits of

x (t)

and

y (t)

N (d u, d t)

is a Poisson counting measure defined on

λ (d u)

. The characteristic measure

λ

on the measure subset

Γ

[0, + \infty]

is such that

λ (Γ) < \infty

N (d u, d t)

is defined on

R_{+} \times (R - {0})

\tilde{N} (d t, d u) = N (d t, d u) - λ (d u) d t

is the corresponding martingale measure. In addition,

γ_{i} (u) (γ_{i} (u) > - 1, \forall u \in Γ)

represents the intensity of Lévy noise changing the prey and predator populations.

Throughout this paper, we give an Assumption 1 that there is a positive constant

P > 0

such that

Assumption 1.

\begin{matrix} \int_{Γ} {[ln (1 + γ_{i} (u))]}^{2} λ (d u) < P, \int_{Γ} γ_{i} {(u)}^{2} λ (d u) < P, i = 1, 2 . \end{matrix}

(3)

We will frequently use the following inequality:

k - 1 - ln k \geq 0, k > 0 .

The rest of the paper is organized as follows. In Section 2, we prove that system 2 has a unique global positive solution. In Section 3, we give the conditions for the extinction and persistence of the prey and predator in the sense of time average. In Section 4, we prove that a systemwithout Lévy noise has ergodic stationary distribution. In Section 5, we give appropriate parameters for numerical simulation to verify the correctness of our theorem. In the end, we summarize this article.

2. Existence and Uniqueness of a Global Positive Solution

The existence and uniqueness of global positive solutions is the basis for studying the dynamic properties of stochastic differential systems. In this section, we will first prove the existence and uniqueness of the local positive solution; then, we prove the global existence and uniqueness of the positive solution of the system.

Lemma 1

([33]). Denote by

Σ (t)

a local martingale vanishing at

t = 0

. Define

\bar{Σ} (t) = \int_{0}^{t} \frac{d 〈 Σ 〉 (s)}{{(1 + s)}^{2}}, t \geq 0,

where

〈 Σ 〉 (t) = 〈 Σ, Σ 〉 (t)

stands for the Meyer’s angle bracket process. If

\underset{t \to + \infty}{lim sup} \bar{Σ} < + \infty

, then

lim_{t \to + \infty} t^{- 1} Σ (t) = 0

, a.s.

Lemma 2.

For

t \in [0, τ_{e})

(

τ_{e}

is the explosion time), model (2) has a unique solution for any initial value

(x_{0}, y_{0}) \in R_{+}^{2}

Proof.

Consider the following equation:

\{\begin{matrix} d M (t) = d ln x (t) = & [a - b e^{M (t)} - \frac{q e^{2 N (t)}}{(m + e^{N (t)}) (c + e^{2 M (t)})} - \frac{1}{2} σ_{1}^{2} \\ + \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u)] d t + σ_{1} d B_{1} (t) + \int_{Γ} ln (1 + γ_{1} (u)) \tilde{N} (d t, d u), \\ d N (t) = d ln y (t) = & [δ - \frac{h e^{N (t)}}{l + e^{M (t)}} - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)] d t \\ + σ_{2} d B_{2} (t) + \int_{Γ} ln (1 + γ_{2} (u)) \tilde{N} (d t, d u), \end{matrix}

(4)

with initial values

M (0) = ln x_{0}

N (0) = ln y_{0}

. Obviously, (4) satisfies the local Lipschitz condition. Therefore, for

t \in [0, τ_{e})

, there is a local solution

(ln x (t), ln y (t))

to (4) with the initial value

(ln x_{0}, ln y_{0})

. It can be known from

I t \hat{o} ’ s

formula that

(x (t), y (t))

is the unique local positive solution of (2) with initial value

(x_{0}, y_{0})

. □

Theorem 1.

For

\forall t \in [0, + \infty)

, model (2) has a unique global positive solution for any initial value

(x_{0}, y_{0}) \in R_{+}^{2}

Proof.

From Lemma 2, we can prove that

τ_{e} = \infty

a . s .

to prove that the solution is global. If

k_{0} > 0

is sufficiently large such that

x_{0}, y_{0} \in [\frac{1}{k_{0}}, k_{0}]

, the following stopping time sequence is defined for each integer

k > k_{0}

τ_{k} = inf \{t \in [0, τ_{e}) : min {x (t), y (t)} \leq \frac{1}{k} o r max {x (t), y (t)} \geq k\} .

Denote

τ_{\infty} = lim_{k \to \infty} τ_{k}

. Then, we can obtain

τ_{\infty} \leq τ_{e}

, a.s. We can obtain

τ_{e} = \infty

a.s. by proving

τ_{\infty} = \infty

a.s. If

τ_{\infty} \neq \infty

, there are constants

T \geq 0

and

ε \in (0, 1)

such that

P (τ_{e} \leq T) > ε

. Then,

\exists k_{1} \geq k_{0}

such that the following holds:

P (τ_{e} \leq T) \geq ε for all k > k_{1} .

Define a

C^{2}

-function

V : R_{+}^{2} \to R_{+}

V (x, y) = (x - 1 - ln x) + (y - 1 - ln y) .

Since

s - 1 - ln s \geq 0

, for all

s > 0

V (x, y)

, is nonnegative. Applying

I t {\hat{o}}^{'} s

formula to

V (x, y)

, we can obtain

\begin{matrix} d V (x, y) = & L V (x, y) d t + σ_{1} (x - 1) d B_{1} (t) + σ_{2} (y - 1) d B_{2} (t) \\ + \int_{Γ} [x (t^{-}) γ_{1} (u) - ln (1 + γ_{1} (u))] \tilde{N} (d t, d u) \\ + \int_{Γ} [y (t^{-}) γ_{2} (u) - ln (1 + γ_{2} (u))] \tilde{N} (d t, d u), \end{matrix}

(5)

where

\begin{matrix} L V (x, y) = & (x - 1) (a - b x (t) - \frac{q y^{2} (t)}{(m + y (t)) (c + x^{2} (t))}) + (y - 1) (δ - \frac{h y (t)}{l + x (t)}) \\ + \frac{1}{2} σ_{1}^{2} + \frac{1}{2} σ_{2}^{2} + \int_{Γ} [γ_{1} (u) - ln (1 + γ_{1} (u))] λ (d u) \\ + \int_{Γ} [γ_{2} (u) - ln (1 + γ_{2} (u))] λ (d u) . \end{matrix}

Applying Assumption 1, we obtain

\begin{matrix} L V (x, y) = & a x - b x^{2} - \frac{q y^{2}}{(m + y) (c + x^{2})} - (a - b x - \frac{q y^{2}}{(m + y) (c + x^{2})}) + \\ δ y - \frac{h y^{2}}{l + x} - (δ - \frac{h y}{l + x}) + \int_{Γ} [γ_{1} (u) - ln (1 + γ_{1} (u))] λ (d u) \\ + \int_{Γ} [γ_{2} (u) - ln (1 + γ_{2} (u))] λ (d u) \\ \leq & - b x^{2} + (a + b) x + \frac{h y}{l} + (\frac{q}{c} + δ) y + \frac{1}{2} σ_{1}^{2} + \frac{1}{2} σ_{2}^{2} - a - δ \\ + \int_{Γ} [γ_{1} (u) - ln (1 + γ_{1} (u))] λ (d u) + \int_{Γ} [γ_{2} (u) - ln (1 + γ_{2} (u))] λ (d u) \\ = & M_{1} + M_{2} y, \end{matrix}

where

\begin{matrix} M_{1} = & max_{x \geq 0} {- b x^{2} + (a + b) x + \frac{h y}{l} + (\frac{q}{c} + δ) y + \frac{1}{2} σ_{1}^{2} \\ + \frac{1}{2} σ_{2}^{2} - a - δ + \int_{Γ} [γ_{1} (u) - ln (1 + γ_{1} (u))] λ (d u) \\ + \int_{Γ} [γ_{2} (u) - ln (1 + γ_{2} (u))] λ (d u)}, \end{matrix}

\begin{matrix} M_{2} = \frac{a}{c} + δ + \frac{h}{l}, \end{matrix}

are positive constants. And because

y \leq 2 (y - 1 - ln y) + 2 ln 2 \leq 2 V (x, y) + 2 ln 2

, the following equation is true:

\begin{matrix} L V (x, y) & \leq M_{1} + M_{2} (2 V (x, y) + 2 ln 2) \\ = M_{1} + 2 ln 2 + 2 M_{2} V (x, y) \\ \leq M (1 + V (x, y)), \end{matrix}

(6)

where

M = m a x {M_{1} + 2 M_{2} ln 2, 2 M_{2}}

. Combining (5) and (6), we can obtain that

\begin{matrix} d V (x, y) \leq & M (1 + V (x, y)) d t + σ_{1} (x - 1) d B_{1} (t) + σ_{2} (y - 1) d B_{2} (t) \\ + \int_{Γ} [x (t^{-}) γ_{1} (u) - ln (1 + γ_{1} (u))] \tilde{N} (d t, d u) \\ + \int_{Γ} [y (t^{-}) γ_{2} (u) - ln (1 + γ_{2} (u))] \tilde{N} (d t, d u) . \end{matrix}

(7)

Integrating the two sides of (7) from 0 to

τ_{k} \land T

, and then taking the expectation, we can obtain

\begin{matrix} E V (x (τ_{k} \land T), y (τ_{k} \land T)) & \leq V (x_{0}, y_{0}) + M E \int_{0}^{τ_{k} \land T} (1 + V (x, y)) d t \\ \leq V (x_{0}, y_{0}) + M T + M E \int_{0}^{τ_{k} \land T} V (x, y) d t . \end{matrix}

Applying the Gronwall’s inequality to the above equation, we obtain

\begin{matrix} E V (x (τ_{k} \land T), y (τ_{k} \land T)) \leq (V (x_{0}, y_{0}) + M T) e^{M} T . \end{matrix}

Let

Ω_{k} = τ_{k} \leq T

; we have

P (Ω_{k}) \geq ε

. Therefore, for

\forall ω \in Ω_{k}

, there is at least one value equal to k or

\frac{1}{k}

x (τ_{k}, Ω)

y (τ_{k}, Ω)

. Note that

V (x (τ_{k}), y (τ_{k})) \geq (k - 1 - ln m) \land (\frac{1}{k} - 1 - ln \frac{1}{k})

. Consequently,

\begin{matrix} (V (x_{0}, y_{0}) + M T) e^{M T} & \geq E (I_{Ω_{k} (ω)}, V (x (τ_{k}), y (τ_{k}))) \\ \geq ε (k - 1 - ln m) \land (\frac{1}{k} - 1 - ln \frac{1}{k}), \end{matrix}

where

I_{Ω_{k} (ω)}

is the characteristic function of

ω_{k}

. When

k \to \infty

, we infer that

\begin{matrix} ε (k - 1 - ln k) \land (\frac{1}{k} - 1 - ln \frac{1}{k}) \to + \infty . \end{matrix}

We derive the contradiction, so we can obtain

τ_{\infty} = \infty

, and the theorem is proved. □

3. Existence and Demise of Biological Populations

In this section, we derive the conditions for the extinction and persistence of prey and predators.

Theorem 2.

Suppose that

(x (t), y (t))

is a positive solution of (2) with initial value

(x_{0}, y_{0})

. If Assumption 1 is satisfied and

a > \frac{1}{2} σ_{1}^{2} - \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u),

δ < \frac{1}{2} σ_{2}^{2} - \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u),

then the prey population x persists in the sense of time average, and the predator population y is extinct, i.e.,

\underset{t \to + \infty}{lim inf} t^{- 1} \int_{0}^{t} x (s) d s > 0, lim_{t \to + \infty} y (t) = 0 .

Proof.

Using

I t {\hat{o}}^{'} s

formula in (2) gives

\begin{matrix} d ln x (t) = & [a - b x (t) - \frac{q y^{2} (t)}{(m + y (t)) (c + x^{2} (t))} - \frac{1}{2} σ_{1}^{2} \\ + \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u)] d t \\ + σ_{1} d B_{1} (t) + \int_{Γ} ln (1 + γ_{1} (u)) \tilde{N} (d t, d u), \end{matrix}

(8)

\begin{matrix} d ln y (t) = & [δ - \frac{h y (t)}{l + x (t)} - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)] d t \\ + σ_{2} d B_{2} (t) + \int_{Γ} ln (1 + γ_{2} (u)) \tilde{N} (d s, d u) . \end{matrix}

(9)

It follows from (9) that

\begin{matrix} d ln y (t) \leq & [δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)] d t + σ_{2} d B_{2} (t) \\ + \int_{Γ} ln (1 + γ_{2} (u)) \tilde{N} (d s, d u) . \end{matrix}

(10)

We take the integral of 0 to t on both sides of (10), and the following holds:

\begin{matrix} ln y (t) \leq & ln y_{0} + \int_{0}^{t} [δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)] d t + σ_{2} B (t) \\ + \int_{0}^{t} \int_{Γ} ln (1 + γ_{2} (u)) \tilde{N} (d s, d u) . \end{matrix}

Then, we obtain

\begin{matrix} \frac{ln y (t)}{t} \leq & [δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)] + \frac{σ_{2} B (t)}{t} \\ + \frac{Σ_{1} (t)}{t} + \frac{ln y_{0}}{t} . \end{matrix}

where

Σ_{1} (t) = \int_{0}^{t} \int_{Γ} ln (1 + γ_{2} (u)) \tilde{N} (d s, d u)

, according to Lemma 1.

〈 Σ_{1}, Σ_{1} 〉 (t) = t \int_{Γ} {[ln (1 + γ_{i} (u))]}^{2} λ (d u) < K t i = 1, 2 .

In view of Lemma 1, we obtain

lim_{t \to + \infty} \frac{Σ_{1} (t)}{t} = 0 .

In other words, we have shown that

\begin{matrix} lim_{t \to + \infty} t^{- 1} ln y (t) \leq δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u) < 0 . \end{matrix}

According to the conclusion of the third section in [34], we can obtain

lim_{t \to + \infty} y (t) = 0 .

Therefore, for

\forall ε > 0

, there exist

t_{0}

and

Ω_{ε}

such that

P (Ω_{ε}) \geq 1 - ε

\frac{q y^{2} (t)}{(m + y) (c + x^{2} (t))} \geq ε

when

t > t_{0}

and

ω \in Ω_{ε}

; thus,

\begin{matrix} x (t) (a - b x (t) - ε) d t + σ_{1} x (t) d B_{1} (t) + \int_{Γ} x (t^{-}) γ_{1} (u) \tilde{N} (d u, d t) \leq d (x) \\ \leq x (t) (a - b x (t)) d t + σ_{1} x (t) d B_{1} (t) + \int_{Γ} x (t^{-}) γ_{1} (u) \tilde{N} (d u, d t) . \end{matrix}

By using the comparison theorem, we can obtain

\begin{matrix} \underset{t \to \infty}{lim inf} t^{- 1} \int_{0}^{t} x (s) d s \geq \frac{a - ε - \frac{1}{2} σ_{1}^{2} + \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u)}{b}, \\ \underset{t \to \infty}{lim sup} t^{- 1} \int_{0}^{t} x (s) d s \leq \frac{a - \frac{1}{2} σ_{1}^{2} + \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u)}{b} . \end{matrix}

By the arbitrariness of

ε

lim_{t \to \infty} t^{- 1} \int_{0}^{t} x (s) d s = \frac{a - \frac{1}{2} σ_{1}^{2} + \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u)}{b} > 0 .

□

Theorem 3.

Suppose that

(x (t), y (t))

is a positive solution of (2) with initial value

(x_{0}, y_{0})

. If Assumption 1 is satisfied and

a < \frac{1}{2} σ_{1}^{2} - \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u),

δ > \frac{1}{2} σ_{2}^{2} - \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u),

then x is extinct, and y persists in the sense of time average, i.e.,

lim_{t \to + \infty} x (t) = 0, \underset{t \to + \infty}{lim inf} t^{- 1} \int_{0}^{t} y (s) d s > 0 .

Proof.

Similar to the proof of Theorem 2, we can obtain

\begin{matrix} lim_{t \to + \infty} t^{- 1} ln x (t) \leq a - \frac{1}{2} σ_{1}^{2} + \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u) < 0 . \end{matrix}

Then, we have

lim_{t \to + \infty} x (t) = 0 .

Therefore, for

\forall ε > 0

, there exist

t_{0}

and

Ω_{ε}

such that

P (Ω_{ε}) \geq 1 - ε

x (t) \leq ε

when

t > t_{0}

and

ω \in Ω_{ε}

; thus,

\begin{matrix} y (t) [δ - \frac{h y (t)}{l}] d t + σ_{2} y (t) d B_{2} (t) + \int_{Γ} y (t^{-}) γ_{2} (u) \tilde{N} (d u, d t) \leq d y (t) \\ \leq y (t) [δ - \frac{h y (t)}{l + ε}] d t + σ_{2} y (t) d B_{2} (t) + \int_{Γ} y (t^{-}) γ_{2} (u) \tilde{N} (d u, d t) . \end{matrix}

(11)

Using the comparison theorem, we obtain

\begin{matrix} \underset{t \to \infty}{lim inf} t^{- 1} \int_{0}^{t} y (s) d s \geq \frac{l (δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u))}{h}, \\ \underset{t \to \infty}{lim sup} t^{- 1} \int_{0}^{t} y (s) d s \leq \frac{l (δ - ε - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u))}{h} . \end{matrix}

From the arbitrariness of

ε

, the following equation can be given:

lim_{t \to \infty} t^{- 1} \int_{0}^{t} y (s) d s = \frac{l}{h} {δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)} > 0 .

□

Theorem 4.

Suppose that

(x (t), y (t))

is a positive solution of (2) with initial value

(x_{0}, y_{0})

. If Assumption 1 is satisfied and

a < \frac{1}{2} σ_{1}^{2} - \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u),

δ < \frac{1}{2} σ_{2}^{2} - \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u),

then x and y are all extinct, i.e.,

lim_{t \to + \infty} x (t) = 0, lim_{t \to + \infty} y (t) = 0 .

Proof.

By the proof of Theorems 2–3, this is obviously true. □

Theorem 5.

Suppose that

(x (t), y (t))

is a positive solution of (2) with initial value

(x_{0}, y_{0})

. If Assumption 1 is satisfied and

δ > \frac{1}{2} σ_{2}^{2} - \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u), q^{2} - 4 (c - q l) \leq 0,

h [a - \frac{1}{2} σ_{1}^{2} + \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u)] > [δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)],

then x and y both persist in the sense of time average, i.e.,

\underset{t \to + \infty}{lim inf} t^{- 1} \int_{0}^{t} x (s) d s > 0, \underset{t \to + \infty}{lim inf} t^{- 1} \int_{0}^{t} y (s) d s > 0 .

Proof.

Integrating both sides of (9), we obtain

\begin{matrix} ln y (t) = & ln y_{0} + [δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)] t + σ_{2} B (t) \\ - \int_{0}^{t} \frac{h y (s)}{l + x (s)} d s + \int_{0}^{t} \int_{Γ} ln (1 + γ_{2} (u)) \tilde{N} (d s, d u), \end{matrix}

(12)

and

\begin{matrix} ln y (t) \geq & ln y_{0} + [δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)] t + σ_{2} B (t) \\ - \int_{0}^{t} \frac{h y (s)}{l} d s + \int_{0}^{t} \int_{Γ} ln (1 + γ_{2} (u)) \tilde{N} (d s, d u) . \end{matrix}

(13)

Define

\tilde{y} (t)

as the solution of the following equation:

\{\begin{matrix} d \tilde{y} (t) = \tilde{y} (t) [δ - \frac{h \tilde{y} (t)}{l}] d t + σ_{2} \tilde{y} (t) d B_{2} (t) + \int_{Γ} \tilde{y} (t^{-}) γ_{2} (u) \tilde{N} (d u, d t), \\ \tilde{y} (0) = y_{0} . \end{matrix}

Then, we can obtain

\begin{matrix} lim_{t \to \infty} t^{- 1} \int_{0}^{t} \tilde{y} (s) d s = \frac{l (δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u))}{h} > 0 . \end{matrix}

Using comparison theorem in (11), we obtain

y (t) \geq \tilde{y} (t)

. Then, there is

\begin{matrix} lim_{t \to \infty} t^{- 1} ln y (t) \geq lim_{t \to \infty} t^{- 1} ln \tilde{y} (t) = 0 . \end{matrix}

(14)

Divide t on both sides of (12) and take the inferior limit, yielding

\begin{matrix} \underset{t \to \infty}{lim inf} t^{- 1} ln y (t) = & δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u) \\ - \underset{t \to \infty}{lim inf} t^{- 1} \int_{0}^{t} \frac{h y (s)}{l + x (s)} d s . \end{matrix}

Bringing (14) into the above formula, we observe that

\begin{matrix} \underset{t \to \infty}{lim inf} t^{- 1} \int_{0}^{t} \frac{y (s)}{l + x (s)} d s \leq \frac{δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)}{h} . \end{matrix}

(15)

Similarly, for (13) we can obtain

\begin{matrix} \underset{t \to \infty}{lim inf} t^{- 1} \int_{0}^{t} y (s) d s \geq \frac{l (δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u))}{h} . \end{matrix}

(16)

Because

q^{2} - 4 (c - q l) \leq 0

, we can deduce that

\begin{matrix} \frac{q y^{2}}{(m + y (t)) (c + x^{2} (t))} \leq \frac{q y (t)}{c + x^{2} (t)} \leq \frac{y (t)}{l + x (t)} . \end{matrix}

(17)

Integrating both sides of (8) and using (17), then we can obtain

\begin{matrix} ln x (t) \geq & ln x_{0} + [a - \frac{1}{2} σ_{1}^{2} + \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u)] t + σ_{1} B_{1} (t) \\ - \int_{0}^{t} \frac{y (s)}{l + x (s)} d s - b \int_{0}^{t} x (s) d s + \int_{0}^{t} \int_{Γ} ln (1 + γ_{1} (u)) \tilde{N} (d s, d u) . \end{matrix}

(18)

From (15), we can obtain

\begin{matrix} t^{- 1} \int_{0}^{t} \frac{y (s)}{l + x (s)} d s \leq \frac{δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)}{h} + ε, \end{matrix}

(19)

where

ε

is an any small positive real number. Similar to the proof of inequality (14), we can obtain

\begin{matrix} lim_{t \to \infty} t^{- 1} ln x (t) \leq 0 . \end{matrix}

Then divide t on both sides of inequality (18); then, from (19), we have

\begin{matrix} \frac{ln x (t)}{t} \geq & \frac{ln x_{0}}{t} + [a - \frac{1}{2} σ_{1}^{2} + \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u)] + \frac{σ_{1} B_{1} (t)}{t} \\ - \frac{δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)}{h} - ε - \frac{b}{t} \int_{0}^{t} x (s) d s \\ + \frac{1}{t} \int_{0}^{t} \int_{Γ} ln (1 + γ_{1} (u)) \tilde{N} (d s, d u) . \end{matrix}

(20)

Take the inferior limit on both sides of (19) at the same time, then using the arbitrariness of

ε

results in

\begin{matrix} lim_{t \to \infty} t^{- 1} \int_{0}^{t} x (s) d s \geq & \frac{1}{b} [a - \frac{1}{2} σ_{1}^{2} + \int_{Γ} [ln (1 + γ_{1} (u)) - γ_{1} (u)] λ (d u)] \\ - \frac{1}{b h} [δ - \frac{1}{2} σ_{2}^{2} + \int_{Γ} [ln (1 + γ_{2} (u)) - γ_{2} (u)] λ (d u)] > 0 . \end{matrix}

□

4. Stationary Distribution without Lévy Noise

Furthermore, when

γ_{i} (u) = 0 (i = 1, 2)

, this means excluding drastic environmental changes. System (2) becomes the following system:

\{\begin{matrix} d x (t) = & x (t) [a - b x (t) - \frac{q y^{2} (t)}{(m + y (t)) (c + x^{2} (t))}] d t + σ_{1} x (t) d B_{1} (t), \\ d y (t) = & y (t) [δ - \frac{h y (t)}{l + x (t)}] d t + σ_{2} y (t) d B_{2} (t) . \end{matrix}

(21)

We consider the stationary distribution and ergodicity of this system. Let

X (t)

be a homogeneous Markov process in

E_{k}

(

E_{k}

denotes the k-dimensional Euclidean space), which can be described by the following stochastic process:

d X (t) = h (X) d t + \sum_{η = 1}^{m} σ_{η} d B_{η} (t) .

Its diffusion matrix is as follows:

Ψ (x) = (a_{i j} (x)), a_{i j} (x) = \sum_{η = 1}^{m} σ_{η}^{i} σ_{η}^{j} .

Assume that there exists a bounded domain

U \subset E_{k}

with regular boundary; then, according to the conclusion in the second section of [35], we can prove that there exists a neighborhood U and a nonnegative function

g (x, y)

such that

L g

is negative for any

E_{h} / U

, which can be used as a sufficient condition to prove that the system has a stationary distribution.

Theorem 6.

For any initial value

(x_{0}, y_{0})

, if the following conditions hold, then System (21) has an ergodic stationary distribution.

δ - \frac{(φ + 1) σ_{2}^{2}}{2} > 0, a + δ - \frac{1}{2} σ_{1}^{2} - \frac{1}{2} σ_{2}^{2} \neq 0 .

Proof.

We define a

C^{2} -

function

g (x, y)

as follows:

g (x, y) = x - F ln x + y - F ln y + y^{- φ} - G (x^{*}, y^{*}), 0 < φ < 1,

where

G (x, y) = x - F ln x + y - F ln y + y^{- φ}

, and

G (x^{*}, y^{*})

is the minimum value of

G (x, y)

F

is a undetermined constant; we will determine its value in the proof process. According to

I t \hat{o} ’ s

formula,

\begin{matrix} L g (x, y) = & x [a - b x - \frac{q y^{2}}{(m + y) (c + x^{2})}] - F [a - b x - \frac{q y^{2}}{(m + y) (c + x^{2})}] \\ + y (δ - \frac{h y}{m + x}) - F (δ - \frac{h y}{m + x}) - φ y^{- φ - 1} \cdot y (δ - \frac{h y}{m + x}) \\ + \frac{F}{2} σ_{1}^{2} + \frac{F}{2} σ_{2}^{2} + \frac{φ (φ + 1)}{2} σ_{2}^{2} y^{- φ} \\ = & - b x^{2} + (a + F b) x + \frac{q F y^{2}}{(m + y) (c + x^{2})} - \frac{q x y^{2}}{(m + y) (c + x^{2})} \\ + δ y - F a - F δ + \frac{- h y^{2} + F h y + φ h y^{- φ + 1}}{m + x} + \frac{F}{2} σ_{1}^{2} + \frac{F}{2} σ_{2}^{2} \\ - & φ (δ - \frac{(φ + 1) σ_{2}^{2}}{2}) y^{- φ} . \end{matrix}

Consider the following bounded regions:

U = \{(x, y) \in R_{+}^{2} | ε < x < \frac{1}{ε}, ε < y < \frac{1}{ε}\} .

We have

R_{+}^{2} / U = U_{1} \cup U_{2} \cup U_{3} \cup U_{4}

, where

\begin{matrix} U_{1} = \{(x, y) \in R_{+}^{2}, x > \frac{1}{ε}\}, & U_{2} = \{(x, y) \in R_{+}^{2}, ε > x > 0\}, \\ U_{3} = \{(x, y) \in R_{+}^{2}, x > ε, y > \frac{1}{ε}\}, & U_{4} = \{(x, y) \in R_{+}^{2}, ε > y > 0\} . \end{matrix}

And

ε

is a small positive number satisfying the following conditions:

\begin{matrix} - \frac{b}{2 ε^{2}} + ▵_{1} + \frac{▵_{2}}{m} < - 1, \\ - \frac{| ▵_{2} |}{m} + (a + F) ε < - 1, \\ ▵_{1} + \frac{▵_{2}}{m} - \frac{h}{2 m ε^{2}} < - 1, \\ ▵_{1} + \frac{▵_{2}}{m} - φ (δ - \frac{(φ + 1) σ_{2}^{2}}{2}) \frac{1}{ε^{φ}} < - 1 . \end{matrix}

where

▵_{1} = sup_{x \in (0, + \infty)} \{- \frac{b}{2} x^{2} + (a + F b) x\} - F a - F δ + \frac{F}{2} σ_{1}^{2} + \frac{F}{2} σ_{2}^{2} < \infty,

▵_{2} = sup_{y \in (0, + \infty)} \{\frac{- h y^{2}}{2} + F h y + m δ y + \frac{q m F y}{c} + φ h y^{- φ + 1}\} < \infty .

It is easy to prove that

▵_{2}

is greater than zero. Next, we will discuss it in four parts.

(1): When $(x, y) \in U_{1}$ , we have

$\begin{matrix} L g (x, y) = & - b x^{2} + (a + F b) x + \frac{q F y^{2}}{(m + y) (c + x^{2})} - \frac{q x y^{2}}{(m + y) (c + x^{2})} \\ + δ y - F a - F δ + \frac{- h y^{2} + F h y + φ h y^{- φ + 1}}{m + x} + \frac{F}{2} σ_{1}^{2} + \frac{F}{2} σ_{2}^{2} \\ - φ (δ - \frac{(φ + 1) σ_{2}^{2}}{2}) y^{- φ} \\ = & - \frac{b}{2} x^{2} - \frac{b}{2} x^{2} + (a + F b) x - F a - F δ + \frac{F}{2} σ_{1}^{2} + \frac{F}{2} σ_{2}^{2} \\ - \frac{q x y^{2}}{(m + y) (c + x^{2})} + \frac{q F y^{2}}{(m + y) (c + x^{2})} + δ y \\ + \frac{\frac{- h y^{2}}{2} + F h y + φ h y^{- φ + 1}}{m + x} - \frac{h y^{2}}{2 (m + x)} - φ (δ - \frac{(φ + 1) σ_{2}^{2}}{2}) y^{- φ} \\ \leq & - \frac{b}{2} x^{2} + ▵_{1} + \frac{▵_{2}}{m} \\ < & - \frac{b}{2 ε^{2}} + ▵_{1} + \frac{▵_{2}}{m} \\ < & - 1 . \end{matrix}$
(2): When $(x, y) \in U_{2}$ , we have

$\begin{matrix} L g (x, y) \leq & \frac{▵_{2}}{m} + (a + F) x - F (a + δ - \frac{1}{2} σ_{1}^{2} + \frac{1}{2} σ_{2}^{2}) . \end{matrix}$

Choose

$F = \frac{2 ▵_{2}}{m (a + δ - \frac{1}{2} σ_{1}^{2} - \frac{1}{2} σ_{2}^{2})} .$

Therefore, we have

$\begin{matrix} L g (x, y) < - \frac{▵_{2}}{m} + (a + F) ε < - 1 . \end{matrix}$
(3): When $(x, y) \in U_{3}$ , we have

$\begin{matrix} L g (x, y) \leq & ▵_{1} + \frac{▵_{2}}{m} - \frac{h}{2 m} y^{2} \\ < & ▵_{1} + \frac{▵_{2}}{m} - \frac{h}{2 m ε^{2}} \\ < & - 1 . \end{matrix}$
(4): When $(x, y) \in U_{4}$ , we have

$\begin{matrix} L g (x, y) \leq & ▵_{1} + \frac{▵_{2}}{m} - φ (δ - \frac{(φ + 1) σ_{2}^{2}}{2}) y^{- φ} \\ < & ▵_{1} + \frac{▵_{2}}{m} - φ (δ - \frac{(φ + 1) σ_{2}^{2}}{2}) \frac{1}{ε^{φ}} \\ < & - 1 . \end{matrix}$

In summary, for

\forall (x, y) \in R_{+}^{2} / U

, we have

L g < - 1

. Moreover, we can also find a constant

M = min_{(x, y) \in U \subset R_{2}^{+}} {σ_{1}^{2} x^{2}, σ_{2}^{2} y^{2}}

that satisfies

σ_{1}^{2} x^{2} ζ_{1}^{2} + σ_{2}^{2} y^{2} ζ_{2}^{2} \geq M | ζ^{2} |, for all (x, y) \in U, ζ = (ζ_{1}, ζ_{2}) \in R^{2} .

At this point, we prove the sufficient condition for the existence of stationary distribution, and then System (21) has a unique stationary distribution. □

5. Numerical Simulations

In this section, we choose the parameters that satisfy the conditions of the theorem to carry out numerical simulation to verify the correctness of the theorem. We take the determined initial values

(x_{0}, y_{0}) = (2.9, 1.4)

and let

σ_{1} = σ_{2} = 0.6

γ_{1} = γ_{2} = 0.06

Example 1.

In order to verify Theorem 2, we choose

a = 4, b = 2, q = 5, m = 2, c = 5, δ = 0.1,

h = 9, l = 1

; then, we have

4 - \frac{{0.6}^{2}}{2} + \int_{Γ} [ln (1 + 0.06) - 0.06] \approx 3.81 > 0

and

0.1 - \frac{{0.6}^{2}}{2} + \int_{Γ} [ln (1 + 0.06) - 0.06] \approx - 0.08 < 0

. Therefore, by Theorem 2, the predator dies out and the prey persists in the sense of time average. This is consistent with Figure 1.

In order to verify Theorem 3, we choose

a = 0.1, b = 3, q = 6, m = 5, c = 10,

δ = 5, h = 5, l = 1

; then, we have

0.1 - \frac{{0.6}^{2}}{2} + \int_{Γ} [ln (1 + 0.06) - 0.06] \approx - 0.08 < 0

and

5 - \frac{{0.6}^{2}}{2} + \int_{Γ} [ln (1 + 0.06) - 0.06] \approx 4.81 > 0

. Therefore, by Theorem 3, the predator persists in the sense of time average and the prey dies out. This is consistent with Figure 2.

In order to verify Theorem 4, we choose

a = 0.1, b = 2, q = 5, m = 2, c = 5, δ = 0.1,

h = 9, l = 1

; then, we have

0.1 - \frac{{0.6}^{2}}{2} + \int_{Γ} [ln (1 + 0.06) - 0.06] \approx - 0.08 < 0

and

0.1 - \frac{{0.6}^{2}}{2} + \int_{Γ} [ln (1 + 0.06) - 0.06] \approx - 0.08 < 0

. Therefore, by Theorem 4, both the predator and the prey die out. This is consistent with Figure 3.

In order to verify Theorem 5, we choose

a = 30, b = 3, q = 6, m = 5, c = 20, δ = 5,

h = 5, l = 1

; then, we have

5 - \frac{{0.6}^{2}}{2} + \int_{Γ} [ln (1 + 0.06) - 0.06] \approx 4.81 > 0

6^{2} - 4 (20 - 6 \times 1) = - 20 < 0

and

5 (30 - \frac{{0.6}^{2}}{2} + \int_{Γ} [ln (1 + 0.06) - 0.06]) - (5 - \frac{{0.6}^{2}}{2} + \int_{Γ} [ln (1 + 0.06) - 0.06]) \approx 144.24 > 0

. Therefore, by Theorem 5, both the predator and the prey persist in the sense of time average. This is consistent with Figure 4. In order to further explore the influence of cooperative hunting and group defense on existence, we only change q and c in the above parameters. Only changing q = 6 to q = 15 makes the cooperative hunting intensity increase and does not meet the theorem conditions; Figure 5 shows that the prey perishes. Only changing c = 30 to c = 50 reduces the population defense strength and does not satisfy the theorem condition; Figure 6 shows that the prey perishes. This tells us that when the intensity of cooperative hunting is too large, or the intensity of group defense is too low, it will lead to the demise of the prey, which is also in line with our common sense.

To verify the condition that System (2) without Lévy jump has an ergodic stationary distribution, we assume that

(x_{0}, y_{0}) = (0.31, 1.725)

, and we choose the following parameter values:

σ_{1} = σ_{2} = 0.1, a = 2, b = 0.7, q = 0.6, m = 0.12, c = 0.1, δ = 2.5, h = 0.75, l = 0.21

. Then, we have

δ - \frac{(φ + 1) σ_{2}^{2}}{2} > 2.5 - 0.01 > 0

and

a + δ - \frac{1}{2} σ_{1}^{2} - \frac{1}{2} σ_{2}^{2} = 0.31 + 1.725 - 0.01 \neq 0

. As shown in Figure 7, the sample paths are concentrated in the elliptical region, indicating that the system is stochastically stable.

6. Conclusions

This paper discusses the dynamic properties of a predator–prey model with cooperative hunting, group defense, and Lévy noise. First, we prove the uniqueness of the global positive solution of System (2) by constructing a suitable equation. Then, we obtain the theorem that the intensity of white noise and Lévy noise more than a certain extent may lead to the extinction of the population; when the intensity of white noise and Lévy noise is small, if the intensity of cooperative hunting and group defense satisfies an inequality, the two populations can continue to coexist. Furthermore, if the intensity of cooperative hunting is too large or the intensity of group defense is too small, the prey will become extinct. Finally, we conclude that the system has a stationary distribution when there is no Lévy noise.

Author Contributions

H.C.: Writing-Original Draft, Methodology, Formal analysis M.L.: Methodology, Writing-Review and Editing X.X.: Data Analysis, Visualization. All authors have read and agreed to the published version of the manuscript.

Funding

This work is supported by College Students Innovations Special Project funded by Northeast Forestry University (202310225237), the Natural Science Foundation of Heilongjiang Province (No. LH2022A002) and the National Natural Science Foundation of China (No. 12071115).

Data Availability Statement

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

Conflicts of Interest

The authors declare that they have no competing interests.

References

Leslie, P. Some further notes on the use of matrices in population mathematic. Biometrica 1948, 35, 213–245. [Google Scholar] [CrossRef]
Leslie, P.; Gower, J.C. The properties of a stochastic model for the predator-prey type of interaction between two species. Biometrika 1960, 47, 219–234. [Google Scholar] [CrossRef]
Aziz-Alaoui, M.; Okiye, M. Boundedness and global stability for a predator-prey model with modified Leslie-Gower and Holling-type II schemes. Appl. Math. Lett. 2003, 16, 1069–1075. [Google Scholar] [CrossRef]
Ble, G.; Castellanos, V. Stable limit cycles in an intraguild predation model with general functional responses. Math. Methods Appl. Sci. 2021, 45, 2219–2233. [Google Scholar] [CrossRef]
Islam, Y.; Shah, F. Functional Response of Harmonia axyridis to the Larvae of Spodoptera litura: The Combined Effect of Temperatures and Prey Instars. Front. Plant Sci. 2022, 13, 849574. [Google Scholar] [CrossRef]
Ble, G.; Guzman-Arellano, C. Coexistence in a four-species food web model with general functional responses. Chaos Solitons Fractals 2021, 153, 111555. [Google Scholar] [CrossRef]
Fu, C.; Xu, T.; Wen, S.; Zhang, S.; Liu, F.; Du, C.H.; Zhang, L. Predation Behaviors and Functional Responses of Picromerus lewisi to the Larvae of Ostrinia furnacalis. Chin. J. Biol. Control 2022, 37, 956–962. [Google Scholar]
Nisal, A.; Diwekar, U.; Hanumante, N.; Shastri, Y.; Cabezas, H. Integrated model for food-energy-water (FEW) nexus to study global sustainability: The main generalized global sustainability model (GGSM). PLoS ONE 2022, 17, e0267403. [Google Scholar] [CrossRef]
Ali, I.; Rasool, G.; Alrashed, S. Numerical simulations of reaction–diffusion equations modeling prey–predator interaction with delay. Int. J. Biomath. 2018, 11, 1850054. [Google Scholar] [CrossRef]
Scheel, D.; Packer, C. Group hunting behaviour of lions: A search for cooperation. Anim. Behav. 1991, 41, 697–709. [Google Scholar] [CrossRef]
Heinsohn, R.; Packer, C. Complex cooperative strategies in group-territorial African lions. Science 1995, 269, 1260–1262. [Google Scholar] [CrossRef]
Schmidt, P.; Mech, L. Wolf pack size and food acquisition. Am. Nat. 1997, 269, 513–517. [Google Scholar] [CrossRef]
Bowman, R. Apparent cooperative hunting in Florida Scrub-Jays. Wilson Bull. 2003, 115, 197–199. [Google Scholar] [CrossRef]
Hannah, K.C. An apparent case of cooperative hunting in immature Northern Shrikes. Wilson Bull. 2005, 117, 407–409. [Google Scholar] [CrossRef]
Chow, Y.; Jang, S. Cooperative hunting in a discrete predator-prey system. J. Biol. Dyn. 2019, 13, 247–264. [Google Scholar] [CrossRef]
Andrews, J. A mathematical model for the continuous culture of microorganisms utilizing inhibitory substrates. Biotechnol. Bioeng. 1968, 10, 707–723. [Google Scholar] [CrossRef]
Sokol, W.; Howell, J. Kinetics of phenol oxidation by washed cells. Biotechnol. Bioeng. 1980, 23, 2039–2049. [Google Scholar] [CrossRef]
Shen, C. Permanence and global attractivity of the food-chain system with Holling IV type functional response. Appl. Math. Comput. 2007, 194, 179–185. [Google Scholar] [CrossRef]
Bai, D.; Tang, J. Global Dynamics of a Predator–Prey System with Cooperative Hunting. Appl. Sci. 2023, 13, 8178. [Google Scholar] [CrossRef]
Du, Y.; Niu, B.; Wei, J. A predator-prey model with cooperative hunting in the predator and group defense in the prey. Am. Inst. Math. Sci. 2022, 27, 5845–5881. [Google Scholar] [CrossRef]
Yao, W.; Li, X. Complicate bifurcation behaviors of a discrete predator–prey model with group defense and nonlinear harvesting in prey. Appl. Anal. 2022, 102, 2567–2582. [Google Scholar] [CrossRef]
Fu, S.; Zhang, H. Effect of hunting cooperation on the dynamic behavior for a diffusive Holling type II predator-prey model. Commun. Nonlinear Sci. Numer. Simul. 2021, 99, 105807. [Google Scholar] [CrossRef]
Pal, S.; Pal, N.; Samanta, S. Chattopadhyay, J. Fear effect in prey and hunting cooperation among predators in a Leslie-Gower model. Ecol. Complex. 2019, 16, 5146–5179. [Google Scholar]
Qi, H.; Meng, X.; Hayat, T.; Hobiny, A. Stationary distribution of a stochastic predator-prey model with hunting cooperation. Appl. Math. Lett. 2022, 124, 107662. [Google Scholar] [CrossRef]
Liu, M.; Du, C.; Deng, M. Persistence and extinction of a modified Leslie–Gower Holling-type II stochastic predator–prey model with impulsive toxicant input in polluted environments. Nonlinear Anal. Hybrid Syst. 2018, 27, 177–190. [Google Scholar] [CrossRef]
Lin, Y.; Jiang, D. Long-time behavior of a stochastic predator-prey model with modified Leslie-Gower and Holling-type II schemes. J. Math. 2016, 9, 121–138. [Google Scholar] [CrossRef]
Yu, M.; Lo, W. Dynamics of microorganism cultivation with delay and stochastic perturbation. Nonlinear Dyn. 2020, 101, 501–519. [Google Scholar]
Higham, D. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Rev. 2001, 433, 525–546. [Google Scholar] [CrossRef]
Slimani, S.; de Fitte, P. Dynamics of a prey-predator system with modified Leslie-Gower and holling type II schemes incorporating a prey refuge. Discret. Contin. Dyn. Syst.-Ser. B 2019, 24, 5003–5039. [Google Scholar]
Han, B.; Jiang, D. Stationary distribution, extinction and density function of a stochastic prey-predator system with general anti-predator behavior and fear effect. Chaos Solitons Fractals 2022, 162, 112458. [Google Scholar] [CrossRef]
Qiao, M.; Yuan, S. Analysis of a stochastic predator-prey model with prey subject to disease and Lévy noise. Stochastics Dyn. 2019, 19, 1950038. [Google Scholar] [CrossRef]
He, X.; Liu, M.; Xu, X. Analysis of stochastic disease including predator-prey model with fear factor and Lévy jump. Math. Biosci. Eng. 2023, 20, 1750–1773. [Google Scholar] [CrossRef] [PubMed]
Lipster, R. A strong law of large numbers for local martingales. Stochastics 1980, 3, 217–228. [Google Scholar]
Khasminskii, R.; Klebaner, F. Long term behavior of solutions of the Lotka-Volterra system under small random perturbations. Ann. Appl. Probab. 2001, 11, 952–963. [Google Scholar] [CrossRef]
Khasminskii, R. Stochastic Stability of Differential Equations. In Stochastic Modeling and Applied Probability; Springer: Berlin/Heidelberg, Germany, 2012; pp. 99–136. [Google Scholar]

Figure 1. We select the following parameter values:

a = 4, b = 2, q = 5, m = 2, c = 5, δ = 0.1,

h = 9, l = 1

. The predator y dies out and the prey x persists in the sense of time average.

Figure 1. We select the following parameter values:

a = 4, b = 2, q = 5, m = 2, c = 5, δ = 0.1,

h = 9, l = 1

. The predator y dies out and the prey x persists in the sense of time average.

Figure 2.

a = 0.1, b = 3, q = 6, m = 5, c = 10, δ = 5, h = 5, l = 1

. The predator y persists in the sense of time average and the prey x dies out.

Figure 2.

a = 0.1, b = 3, q = 6, m = 5, c = 10, δ = 5, h = 5, l = 1

. The predator y persists in the sense of time average and the prey x dies out.

Figure 3.

a = 0.1, b = 2, q = 5, m = 2, c = 5, δ = 0.1, h = 9, l = 1

. Then, both the predator y and the prey x die out.

Figure 3.

a = 0.1, b = 2, q = 5, m = 2, c = 5, δ = 0.1, h = 9, l = 1

. Then, both the predator y and the prey x die out.

Figure 4. We select the following parameter values:

a = 30, b = 3, q = 6, m = 5, c = 50, δ = 5,

h = 5, l = 1

. Then, both the predator y and the prey x persist in the sense of time average.

Figure 4. We select the following parameter values:

a = 30, b = 3, q = 6, m = 5, c = 50, δ = 5,

h = 5, l = 1

. Then, both the predator y and the prey x persist in the sense of time average.

Figure 5.

a = 30, b = 3, q = 6, m = 5, c = 20, δ = 5, h = 5, l = 1

. When letting c decrease so that the group defense becomes larger, the predator perishes.

Figure 5.

a = 30, b = 3, q = 6, m = 5, c = 20, δ = 5, h = 5, l = 1

. When letting c decrease so that the group defense becomes larger, the predator perishes.

Figure 6.

a = 30, b = 3, q = 15, m = 5, c = 50, δ = 5, h = 5, l = 1

. Increasing q makes the cooperative hunting intensity increase, and the prey perishes.

Figure 6.

a = 30, b = 3, q = 15, m = 5, c = 50, δ = 5, h = 5, l = 1

. Increasing q makes the cooperative hunting intensity increase, and the prey perishes.

Figure 7.

σ_{1} = σ_{2} = 0.1, a = 2, b = 0.7, q = 0.6, m = 0.12, c = 0.1, δ = 2.5, h = 0.75, l = 0.21

. The distribution of the sample path in the phase space.

Figure 7.

σ_{1} = σ_{2} = 0.1, a = 2, b = 0.7, q = 0.6, m = 0.12, c = 0.1, δ = 2.5, h = 0.75, l = 0.21

. The distribution of the sample path in the phase space.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Chen, H.; Liu, M.; Xu, X. Dynamics of a Prey–Predator Model with Group Defense for Prey, Cooperative Hunting for Predator, and Lévy Jump. Axioms 2023, 12, 878. https://doi.org/10.3390/axioms12090878

AMA Style

Chen H, Liu M, Xu X. Dynamics of a Prey–Predator Model with Group Defense for Prey, Cooperative Hunting for Predator, and Lévy Jump. Axioms. 2023; 12(9):878. https://doi.org/10.3390/axioms12090878

Chicago/Turabian Style

Chen, Hengfei, Ming Liu, and Xiaofeng Xu. 2023. "Dynamics of a Prey–Predator Model with Group Defense for Prey, Cooperative Hunting for Predator, and Lévy Jump" Axioms 12, no. 9: 878. https://doi.org/10.3390/axioms12090878

APA Style

Chen, H., Liu, M., & Xu, X. (2023). Dynamics of a Prey–Predator Model with Group Defense for Prey, Cooperative Hunting for Predator, and Lévy Jump. Axioms, 12(9), 878. https://doi.org/10.3390/axioms12090878

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Dynamics of a Prey–Predator Model with Group Defense for Prey, Cooperative Hunting for Predator, and Lévy Jump

Abstract

1. Introduction

2. Existence and Uniqueness of a Global Positive Solution

3. Existence and Demise of Biological Populations

4. Stationary Distribution without Lévy Noise

5. Numerical Simulations

6. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI