Bifurcation analysis of a Leslie-type predator–prey system with prey harvesting and group defense

Abstract

In this paper, we investigate a Leslie-type predator–prey model that incorporates prey harvesting and group defense, leading to a modified functional response. Our analysis focuses on the existence and stability of the system’s equilibria, which are essential for the coexistence of predator and prey populations and the maintenance of ecological balance. We identify the maximum sustainable yield, a critical factor for achieving this balance. Through a thorough examination of positive equilibrium stability, we determine the conditions and initial values that promote the survival of both species. We delve into the system’s dynamics by analyzing saddle-node and Hopf bifurcations, which are crucial for understanding the system transitions between various states. To evaluate the stability of the Hopf bifurcation, we calculate the first Lyapunov exponent and offer a quantitative assessment of the system’s stability. Furthermore, we explore the Bogdanov–Takens (BT) bifurcation, a co-dimension 2 scenario, by employing a universal unfolding technique near the cusp point. This method simplifies the complex dynamics and reveals the conditions that trigger such bifurcations. To substantiate our theoretical findings, we conduct numerical simulations, which serve as a practical validation of the model predictions. These simulations not only confirm the theoretical results but also showcase the potential of the model for predicting real-world ecological scenarios. This in-depth analysis contributes to a nuanced understanding of the dynamics within predator–prey interactions and advances the field of ecological modeling.

1 Introduction

Since the pioneering work of Lotka and Volterra, who introduced a pair of differential equations to describe predator–prey dynamics, predator–prey models have attracted considerable interest from both mathematicians and biologists due to their widespread applicability [1]. Predominantly, predator–prey models have been developed around two main scenarios: (i) the growth function of the predator is directly proportional to its predation function, capturing the functional response of predators to prey availability, and (ii) the growth function of the predator is decoupled from its predation function, suggesting a more intricate relationship between the predator and prey.

Moreover, the growth term of the predator in these models is often portrayed as a function that depends not only on prey density but also on the predator-to-prey ratio, adding another dimension of complexity to the interactions. Among these diverse models, the Leslie-type model has emerged as a particularly influential paradigm. Its sophisticated treatment of predator–prey interactions has solidified its role as a fundamental framework in ecological studies, aiding in the comprehension of the delicate equilibrium within ecosystems.

[2] proposed a Leslie-type model to describe the relationship between predators and their prey:where x and y denote the population densities of the prey and predators, respectively. In the absence of predators, the growth of prey populations is commonly modeled using a logistic growth function, characterized by an intrinsic growth rate r and an environmental carrying capacity K. This model assumes that population growth is limited by the availability of resources as the population size approaches the carrying capacity. The functional response, represented by f(x), describes how the rate of prey consumption by predators varies with the density of the prey population. In contrast, the growth rate of the predator is a Leslie-type growth function expressed as , which depends on both the prey density and the predator density. Here, q denotes the intrinsic per capita growth rate of the predator in the presence of abundant prey, p is a scaling factor that relates the carrying capacity of the predator to that of the prey, and the term reflects the diminishing growth rate as the predator-to-prey ratio approaches one, indicating prey scarcity. This growth function is particularly useful for simulating systems where the success of the predator in capturing the prey is influenced not only by prey abundance but also by the density of the predator population. It captures the competitive pressures and the complex dynamics of resource allocation within an ecosystem, offering a more nuanced understanding of predator–prey interactions. The functional response function f(x) is typically categorized into four classical types: Holling type I [3], Holling type II [4], Holling type III [5, 6], and Holling type IV [7, 8]. These types are defined by specific mathematical expressions that reflect different assumptions about predator behavior and the efficiency of predation. Most predator–prey models incorporate one of these functional responses or their modifications.

Recently, [9, 10] proposed a novel functional response model that replaces the traditional prey density with the square root of prey density in the context of the Holling type II functional response. This innovation is particularly pertinent for prey species that exhibit herding behavior, where only peripheral individuals engage with predators. This revised functional response more accurately reflects how the protective effects of group living affect the rate at which predators encounter and consume the prey. Combining their functional response functions, we propose the following Leslie-type model:where α presents the search efficiency of the predator for the prey and T_h denotes the average handing time for each prey.

Before going into details, by substitutinginto model Eq. 1 and dropping the bars, we obtain

Over the past three decades, a considerable amount of research has been devoted to understanding the effects of harvesting on the dynamics of predator–prey systems and its implications for the management of renewable resources. [11] explored a predator–prey model that includes Michaelis–Menten-type predator harvesting. In a similar vein, [12] discussed a modified Leslie–Gower model incorporating Michaelis–Menten-type prey harvesting. [13] investigated a model where both predator and prey populations exhibit logistic growth and are subject to nonlinear harvesting. The impact of proportional harvesting on the dynamic behavior of predator–prey models was examined by [14, 15]. [16] provided a detailed bifurcation analysis for a model with Holling-type II functional response and a constant harvesting rate. [17] also studied a model with constant prey harvesting. [18, 19] studied the dynamic behaviors of the predator–prey systems with a Holling- and Leslie-type or Leslie–Gower model with constant-yield prey harvesting, revealing a diverse array of bifurcations. [20] and [21] discussed ratio-dependent predator–prey systems with constant harvesting for both prey and predators, respectively. [22] investigated a system where the prey population forms herds as a defense against predators, and both species are exposed to a constant harvesting rate.

Inspired by these studies, we introduce the impact of a constant harvesting rate for the prey into Eq. 2 and obtain the following results:where h represents the rate of prey harvesting.

It should be noted that system Eq. 2 is a Leslie-type predator–prey model that incorporates prey group defense mechanisms. Furthermore, system Eq. 3 extends this framework by accounting for the constant harvesting of the prey within system Eq. 2. To date, the dynamic behaviors of both systems Eq. 2 and Eq. 3 remain unexplored. Given the significance of biological resources as a sustainable source for human utilization and considering the rapid industrialization and population growth that have led to an inevitable increase in the utilization of various biological assets, it is imperative to investigate the impact of harvesting behavior on the population of both the predator and prey, particularly when prey group defense is involved, for the sake of maintaining ecological balance. In this paper, we study the dynamic behaviors of these two systems. Through rigorous mathematical analysis, we find that the population of the prey will be sustained without the effect of prey harvesting. However, due to the emergence of prey harvesting, system Eq. 3 exhibits complex dynamic behaviors and undergoes saddle-node, Hopf, and Bogdanov–Takens (BT) bifurcations. System Eq. 3 has bi-stable behavior due to saddle-node bifurcation, and whether the prey becomes extinct or not depends on the prey harvesting intensity. Specifically, we determine the maximum sustainable yield as . Once the prey harvesting rate h exceeds h_MSY, both the prey and predator species face extinction. This insight underscores the importance of selecting an appropriate harvesting rate, specifically to ensure the co-existence of the predator and prey, thereby sustaining ecological balance.

The structure of this paper is shown in Figure 1. Section 2 discusses the existence and stability of the equilibria for system Eq. 2. Section 3 delves into the investigation of the equilibria for system Eq. 3, providing both analytical insights and corroborating numerical results. Section 4 focuses on the bifurcation analysis of system Eq. 3, examining both saddle-node and Hopf bifurcations. Additionally, this section explores the Bogdanov–Takens bifurcation within the same system. Concluding remarks are given in Section 5, summarizing the key findings and implications of the study.

FIGURE 1

2 Equilibria of system (3)

2.1 Existence of equilibria

To determine the equilibria of system Eq. 2, we analyze the prey and predator nullclines, which are given by the following equations:

From these nullclines, it is evident that system Eq. 2 has a unique boundary equilibrium at . To identify potential positive equilibria, we solve the following equations:

By substituting , we transform the problem of finding positive solutions to f(x) = 0 into finding positive solutions to g(z) = adz³ + dz² + (1 − ad)z − d = 0. Through algebraic analysis, we can easily determine that g(z) = 0 has a unique positive solution. Consequently, system Eq. 2 has a unique positive equilibrium point, denoted as , where , , and z₀ is the unique positive root of g(z) = 0. Based on this analysis, we can formulate the following theorem regarding the number of equilibria in system Eq. 2.

2.2 Stability of equilibria

Proof. The Jacobian matrix of system Eq. 2 evaluated at is given by

It is obvious that the matrix has a positive eigenvalue b and a negative eigenvalue −1. Therefore, is always a hyperbolic saddle.

The Jacobian matrix of system Eq. 2 evaluated at is given by

According to g(z₀) = 0, we can obtain . Then,

Hence, when b > b₀, , the matrix has two negative part eigenvalues, and is a sink. When b < b₀, then , the matrix has two positive part eigenvalues, and is a source. When b = b₀, then , the matrix has a pair of purely imaginary eigenvalues, and is a center.

3 Equilibria of system (4)

3.1 Existence of equilibria

This section conducts a qualitative analysis of model Eq. 3. With the biological context in mind, we focus on examining the dynamics of system Eq. 3 within the closed first quadrant, denoted as , in the (x, y) coordinate plane.

In order to obtain the equilibria of system Eq. 3, we first study the roots of the following equations:

Regarding the count of boundary equilibria for system Eq. 3, we derived the subsequent theorem.

Proof. When h exceeds , . This implies that the dynamics of system Eq. 3 within are inherently trivial. Consequently, the prey species faces extinction, an event that subsequently triggers the extinction of the predator. It is evident that system Eq. 3 lacks any equilibrium points within the positive quadrant .

We next explore the existence of boundary equilibria under the condition that . Specifically, when y = 0 in Eq. 4, we are confronted with the following equation:

It is straightforward to determine that when , Eq. 5 has two distinct roots, denoted as x₁ and x₂. Conversely, when , the equation yields a unique root, which is .

Note that

Regarding the positive equilibria of system Eq. 3, we have the following theorem.

Proof. When y ≠ 0, then according to Eq. 5, and x is the positive root of the equationwhere and f₁(z) is given in Eq. 6.

Through algebraic analysis, we find that the equation h₁(z) = 0 has a unique positive root, denoted as z₁. When g₁(z₁) < 0, which corresponds to , the equation g₁(z) = 0 admits two distinct positive roots, z₂ and z₃ (with z₂ < z₃). If f₁(z₃) = 0, which occurs when , then the equation f₁(z) = 0 has a unique positive root. Consequently, system Eq. 3 admits a unique positive equilibrium E₃ = (x₃, y₃), where and . If f₁(z₃) < 0, corresponding to , then the equation f₁(z) = 0 has two distinct positive roots, z₄ and z₅ (with z₂ < z₄ < z₃ < z₅). In this case, system Eq. 3 has two distinct positive equilibria, E₄ = (x₄, y₄) and E₅ = (x₅, y₅), where , , , and .

3.2 Stability of equilibrium E₀

FIGURE 2

Proof. The Jacobian matrix J(E₀) of system Eq. 3 evaluated at E₀ is

Given that Det[J(E₀)] = 0, we can conclude that the equilibrium is nonhyperbolic, indicating that it is a degenerate equilibrium. It is evident that the eigenvalues of the Jacobian matrix J(E₀) are λ₁ = 0 and λ₂ = b. To ascertain the stability of E₀, we perform a transformation by shifting the equilibrium to the origin using the transformation . We then expand system Eq. 3 in a power series up to the second order around the origin. Under this transformation, system Eq. 3 becomeswhere P₁(X, Y) and Q₁(X, Y) are smooth functions of at least the third order with respect to (X, Y).

Then, making the transformationand introducing a new time variable τ = bt, for which we retain t to denote τ for notational simplicity, we obtainwhere P₂(x, y) and Q₂(x, y) are smooth functions of at least the third order with respect to (x, y) and

According to Theorem 7.1 from Chapter 2 in [23], the equilibrium E₀ is classified as a saddle-node. This implies that the vicinity of E₀ is bifurcated by two separatrices that approach E₀ from above and below. One region forms a parabolic sector, while the other comprises two hyperbolic sectors. Additionally, the parabolic sector is situated in the right half-plane, and E₀ acts as a repelling saddle-node.

3.3 Stability of equilibria E₁ and E₂

Theorem 3.4. When, system Eq. 3has two different boundary equilibriaE₁andE₂. Moreover,E₁is always a hyperbolic unstable node, andE₂is always a hyperbolic saddle. The corresponding phase portraits are shown inFigure 2B and Figure 3.

FIGURE 3

Proof. The Jacobian matrix J(E₁) of system Eq. 3 evaluated at E₁ is

We can obtain that the eigenvalues are and λ₂ = b > 0. Hence, E₁ is a hyperbolic unstable node.

The Jacobian matrix J(E₂) of system Eq. 3 evaluated at E₂ is

It is easy to observe that the eigenvalues are and λ₂ = b > 0. Hence, E₂ is a hyperbolic saddle.

3.4 Stability of equilibrium E₃

Now, we investigate the stability of equilibrium E₃. The Jacobian matrix J(E₃) of system Eq. 3 evaluated at E₃ is

Straightforward computation shows that, due to ,where

Since f(z₃) = 0, from Theorem 3.2, we can obtain

Substituting g₁(z₃) in the above equation yields . Because z₃ is a positive root of g₁(z) = 0, g₁(z₃) = 0 and l(z₃) = 0. Therefore, Det[J(E₃)] = 0. Therefore, we can establish that equilibrium E₃ is nonhyperbolic, which means it is a degenerate equilibrium. To elaborate further, we present the following theorem:

Proof. Since Tr[J(E₃)] = a₁₀ + b₀₁, we know that the eigenvalues of J(E₃) are λ₁ = 0 and λ₂ = a₁₀ + b₀₁.

Similar to the proof outlined in Theorem 3.3, we perform a transformation that shifts equilibrium E₃ to the origin. We then expand system Eq. 3 in a power series up to the second order around this new origin. This transformation yields a simplified system thatwhere P₁₁(X, Y) and Q₁₁(X, Y) are smooth functions of at least the third order with respect to (X, Y), and

When a₁₀ + b₀₁ ≠ 0, J(E₃) has one zero eigenvalue and a nonzero eigenvalue. To further analyze the stability of this equilibrium, we introduce a transformation thatand introducing a new time variable τ₁ = (a₁₀ + b₀₁)t, for which we retain t to denote τ₁ for notational simplicity, we obtainwhere P₂₁(x, y) and Q₂₁(x, y) are smooth functions of at least the third order with respect to (x, y), and

Straightforward computation shows that

Consequently, as shown by Theorem 7.1 from Chapter 2 in [23], if a₁₀ + b₀₁ ≠ 0, then the equilibrium E₃ qualifies as a saddle-node. Moreover, the characteristics of this saddle-node are determined as follows: (i) if a₁₀ + b₀₁ > 0, the parabolic sector is situated in the upper half-plane, and E₃ acts as a repelling saddle-node and (ii) conversely, if a₁₀ + b₀₁ < 0, the parabolic sector remains in the upper half-plane, but E₃ is an attracting saddle-node.

When a₁₀ + b₀₁ = 0, J(E₃) has one zero eigenvalue with two multiple. Then, making the transformationand introducing a new time variable τ₂ = −b₀₁t = bt. For notational simplicity, we retain t to represent τ₂, then system Eq. 7 becomeswhere P₂₂(x, y) and Q₂₂(x, y) are smooth functions of at least the third order with respect to (x, y) and

To eliminate y² terms in system Eq. 8, we make the following near identity transformation:which transforms system Eq. 8 intowhere P₂₃(X, Y) and Q₂₃(X, Y) are smooth functions of at least the third order with respect to (X, Y). Therefore, when f₂₀ ≠ 0 and f₁₁ ≠ 0 , according to the results presented by [24], E₃ is a cusp of co-dimension 2. If or f₁₁ = 0 , E₃ is a cusp of co-dimension at least 3.

3.5 Stability of equilibria E₄ and E₅

Proof. The Jacobian matrix J(E₄) of system Eq. 3 evaluated at E₄ is

Straightforward computation shows that, due to ,where

Since f(z₄) = 0, from Theorem 3.2, we obtain

Substituting g₁(z₄) by the above equation, then . Because z₂ < z₄ < z₃, z₂ and z₃ are two positive roots of g₁(z) = 0, then g₁(z₄) < 0 and l(z₄) < 0. Therefore, Det[J(E₄)] < 0. Hence, the matrix J(E₄) has a positive eigenvalue and a negative eigenvalue. Hence, E₄ is a hyperbolic saddle.

Similarly, we obtainsince l(z₅) > 0 (z₅ > z₃). Moreover, we can also obtain

Hence, when b > b*, Tr[J(E₅)] < 0, the matrix J(E₅) has two negative part eigenvalues, and E₅ is a sink. When b < b*, then Tr[J(E₅)] > 0, the matrix J(E₅) has two positive part eigenvalues, and E₅ is a source. When b = b*, then Tr[J(E₅)] = 0, the matrix J(E₅) has a pair of pure imaginary eigenvalues, and E₅ is a center.

4 Bifurcation analysis

In this section, we investigate the bifurcations that take place in system Eq. 3.

4.1 Saddle-node bifurcation

According to Theorem 3.1, we provided the conditions for the existence of E₁ and E₂ based on some restrictions. We obtain thatis a saddle-node bifurcation surface. Consequently, we can deduce that on the surface SN₁, system Eq. 3 possesses a unique equilibrium E₀, which is a saddle-node. As the parameters traverse this surface, the quantity of boundary equilibria undergoes a transition from zero to two. Biologically, this bifurcation corresponds to a critical threshold where the maximum sustainable yield (MSY) is at . Beyond this threshold , the prey species faces extinction, leading to the collapse of the system. However, under certain initial conditions where , the prey species can avoid extinction [18, 20]. This phenomenon is shown in the saddle-node bifurcation diagram, as shown in Figure 4B.

FIGURE 4

Similar to the above discussion, we obtainaccording to Theorem 3.2, is another saddle-node bifurcation surface. When the parameters transition from one side of the surface to the other, the number of positive equilibria in system Eq. 3 changes from zero to two. The saddle-node bifurcation is illustrated in Figure 4A.

4.2 Hopf bifurcation

From Theorem 3.6, we know that the positive equilibrium E₅ is a center-type nonhyperbolic equilibrium when b = b* and . Therefore, system (4) may undergo a Hopf bifurcation. To proceed, we first confirm the transversality condition necessary for a Hopf bifurcation,

Therefore, the stability of the positive equilibrium E₅ undergoes a change as the parameters traverse the specific critical surface

To determine the direction of the Hopf bifurcation, we calculate the first Lyapunov number l₁ at the equilibrium E₅ = (x₅, y₅). Initially, we shift the equilibrium to the origin using the transformation . Subsequently, we expand system Eq. 3 into a power series around the origin. Then, system Eq. 3 becomeswhere and are smooth functions of at least the fourth order with respect to and

In accordance with the formula for the first Lyapunov number l₁ presented in [24] [p.353], we proceed with the computationwhereand

Consequently, if l₁ < 0 (which corresponds to M < 0), then E₅ undergoes a supercritical Hopf bifurcation, leading to the emergence of a stable limit cycle in a small neighborhood of E₅ as the parameters across the surface H₁. Conversely, if l₁ > 0 (or M > 0), a subcritical Hopf bifurcation occurs, resulting in an unstable limit cycle appearing in the vicinity of E₅ as the parameters across the surface H₁.

Given the complexity of l₁, which does not readily reveal information about its sign, we turn to a numerical example for further clarification.

Fixing a = 0.4, d = 0.6, and h = 0.04, we obtain b* = 0.3960214438 and M = 0.537259605 and then l₁ > 0, which implies that an unstable limit cycle is created around E₅ = (0.1742665424, 00.2904442373). For b = 0.41 > b*, we obtain M = 0.652746120 and then l₁ > 0, which implies that an unstable limit cycle is created around E₅, where E₅ is asymptotically stable (see Figure 5B). When b = 0.425, a collision occurs between the periodic orbit and saddle E₄, resulting in the formation of a homoclinic orbit. Since M = 0.791862704, then l₁ > 0, which means that periodic orbits remain unstable, as shown in Figure 5A.

FIGURE 5

4.3 Bogdanov–Takens bifurcation

As presented in Theorem 3.5, the positive equilibrium E₃ is identified as a cusp of co-dimension 2. This classification holds under certain conditions where a₁₀ + b₀₁ = 0, , f₂₀ ≠ 0, and f₁₁ ≠ 0. Consequently, system Eq. 3 may undergo a Bogdanov–Takens bifurcation in the immediate vicinity of E₃.

Proof. To obtain the expressions for the saddle-node, Hopf, and homoclinic bifurcation curves, we derive a normal form for the Bogdanov–Takens bifurcation. This derivation uses the techniques outlined by [11, 18, 25].

Let us examine a perturbation of the parameters h and d centered around their BT bifurcation values. We represent this perturbation as h = h_BT + u and d = d_BT + v, where u and v are small deviations from the critical values h_BT and d_BT, respectively. Consider

We begin by translating the equilibrium point E₃ = (x₃, y₃) to the origin using the transformations X = x − x₃ and Y = y − y₃. Subsequently, we expand system Eq. 9 as a power series centered at the origin. Then, we havewhere P₁(X, Y) and Q₁(X, Y) are smooth functions of at least the third order with respect to (X, Y), whose coefficients depend smoothly on u and v, and where

Let x = X and y = m₀₀ + m₁₀X + m₀₁Y + m₁₁XY + m₂₀X² + P₁(X, Y), then system Eq. 10 can be written aswhere Q₂(x, y) are smooth functions of at least the third order with respect to (x, y), whose coefficients vary smoothly on u and v, and

By introducing a new variable τ through the relation dt = (1 − l₀₂x)dτ and subsequently rewriting τ back as t, we can recast system Eq. 11 in a new form:

Let X = x and Y = y(1 − l₀₂x), then system Eq. 12 can be written aswhere Q₃(X, Y) are smooth functions of at least the third order with respect to (X, Y), whose coefficients depend smoothly on u and v, and

Note that d₂₀ is very complex, and we cannot determine the sign of d₂₀ when u and v are small. Hence, we consider two cases: C-1: d₂₀ > 0; C-2: d₂₀ < 0.

C-1: If d₂₀ > 0 for small values of u and v, we perform the following change of variables: x = X, , and . Additionally, we retain variable t to represent τ. With these transformations, system Eq. 13 is transformed into a new formwhere Q₄(x, y) are smooth of at least the third order with respect to (x, y), whose coefficients depend smoothly on u and v, and

Let and Y = y, then system Eq. 14 can be written aswhere Q₅(X, Y) are smooth of at least the third order with respect to (X, Y), whose coefficients depend smoothly on u and v, and

Suppose d₁₁ ≠ 0, then g₁₁ = h₁₁ ≠ 0. Making the change of variables , , and , we then obtain from system Eq. 15 thatwhere Q₆(x, y) are smooth of at least the third order with respect to (x, y), whose coefficients depend smoothly on u and v, and

C-2: If d₂₀ < 0 for small u and v, we perform the following change of variables: x′ = X, , and . Additionally, we retain variable t to represent τ′. With these transformations, system Eq. 12 is transformed into a new form:where are smooth of at least the third order with respect to (x′, y′), whose coefficients depend smoothly on u and v, and

Let and Y′ = y′, then system Eq. 17 can be written aswhere are smooth of at least the third order with respect to (X′, Y′), whose coefficients depend smoothly on u and v, and

Suppose d₁₁ ≠ 0, . Making , , and , we obtain the versal unfolding of system Eq. 18:where are smooth of at least the third order with respect to (x′, y′), whose coefficients depend smoothly on u and v, and

To streamline the analysis, we retain the notation μ₁ and μ₂ for the transformed parameters and , as defined in Equation Eq. 19. This approach reduces the number of cases we must consider. If the Jacobian determinant is non-zero, then the parameter transformations described by Eqs 16 and 19 are homeomorphisms in a small neighborhood of the origin. Consequently, μ₁ and μ₂ can be treated as independent parameters.

Drawing upon the findings obtained by [24] and as supported by [5, 11, 18, 26], we establish that system Eq. 9 experiences a Bogdanov–Takens bifurcation when the parameters (u, v) are in a small neighborhood of the origin. The local representatisons of the bifurcation curves, derived from the analysis, are as follows (“+” for d₂₀ > 0, “-” for d₂₀ < 0):

• (1) The saddle-node bifurcation curve SN = {(u, v): μ₁ = 0, μ₂ ≠ 0};

• (2) The Hopf bifurcation curve ;

• (3) The homoclinic bifurcation curve .

These bifurcation curves H, HL, and SN divide the small neighborhood of the origin in the parameter plane (u, v) into four regions. The bifurcation diagram of system Eq. 9 is given in Figure 6A, from which we obtain the following observations:

FIGURE 6

• (a) When (u, v) = (0, 0), the unique positive equilibrium is a cusp of co-dimension 2.

• (b) When the parameters are situated within region I, the system lacks a positive equilibrium, as shown in Figure 6B. This absence implies that the prey population is inclined toward extinction under nearly all initial conditions.

• (c) When the parameters are positioned along the curve SN, the system exhibits a positive equilibrium that is characterized as a saddle node.

• (d) In region II, system Eq. 9 possesses two positive equilibria. Specifically, E₄ is identified as a saddle, while E₅ acts as a source, as shown in Figure 7A.

• (e) Upon crossing curve H and entering region III, the system undergoes a subcritical Hopf bifurcation, resulting in the emergence of an unstable limit cycle. This phenomenon is depicted in Figure 7B.

• (f) When the parameters are selected along the curve HL, an unstable homoclinic cycle encircles a sink, while the other equilibrium remains a saddle. This configuration is illustrated in Figure 8A.

• (g) In region IV, the system is characterized by the presence of a sink and a saddle, as depicted in Figure 8B.

FIGURE 7

FIGURE 8

5 Conclusion

In this paper, we explore a Leslie-type predator–prey system incorporating prey harvesting. Initially, we examine system Eq. 2 without prey harvesting, establishing that it possesses a unique boundary equilibrium and a unique positive equilibrium. Theorem 2.2 elucidates the stability of both boundary equilibrium and positive equilibrium. Notably, the boundary equilibrium is consistently unstable, indicating that neither the prey nor the predator will ever become extinct without the influence of prey harvesting.

Furthermore, we study the dynamic behaviors of system Eq. 3 to analyze the effect of prey harvesting on the prey and predator. Hence, we first analyze the existence and stability of equilibria for system Eq. 3 in Section 3. As shown in Theorem 3.1, we identify the maximum sustainable yield as . Beyond this threshold , both the prey and predator species face extinction. Our investigation into the stability of the positive equilibria reveals that the coexistence of both species is possible if and the harvesting rate is between and h_MSY, or if b⩾b* and h is less than the minimum of and . Otherwise, either only the prey survives or both species become extinct. This insight underscores the importance of selecting an appropriate harvesting rate, specifically one that is below a critical threshold and less than h_MSY to ensure the coexistence of the predator and prey, thereby sustaining ecological balance.

In Section 4, we chose the constant harvesting rate h as one of the bifurcation parameters to perform the bifurcation analysis of system Eq. 3. We find that system Eq. 3 has complex dynamic behavior and undergoes Hopf bifurcation, Bogdanov–Takens bifurcation of co-dimension 2, and saddle-node bifurcation. Among them, the emergence of saddle-node bifurcation means that system Eq. 3 has bi-stable behavior, and the prey may face extinction or sustain, which depends on the harvesting rate. For the Hopf bifurcation, we calculate the first Lyapunov number to ascertain the direction of the Hopf bifurcation, which indicates whether a stable or unstable limit cycle emerges in the vicinity of the equilibrium E₅. By selecting two parameters from system Eq. 3 as bifurcation parameters, we demonstrate that the system undergoes a Bogdanov–Takens bifurcation of co-dimension 2. This conclusion is drawn from our analysis of the universal unfolding near the cusp.

In conclusion, we constructed and provided a detailed dynamic analysis of a Leslie-type predator–prey system that incorporates prey harvesting. Our findings reveal that prey harvesting introduces complex dynamic behaviors into the system, potentially leading to a series of bifurcations. Furthermore, if the harvesting intensity exceeds a critical threshold, the prey may face extinction. These insights offer valuable contributions to the understanding of predator–prey systems and ecological balance. However, it is worth noting that this paper focuses primarily on qualitative analysis and lacks empirical data to support our theoretical findings. Despite this, significant algorithms have been developed to gather and analyze real-world data. We intend to leverage these algorithms to quantitatively assess the impact of harvesting using the methodologies introduced by [27, 28]. Additionally, this study does not consider the removal of both prey and predator species, and human activity can significantly impact ecosystems, as highlighted by [29]. Drawing inspiration from [30–36], our model can be extended to eco-epidemic systems and reaction–diffusion systems. These areas represent promising directions for future research in this field.

References

• 1.

  NingL-YLuoX-FLiB-LWuY-PSunG-QFengT-CAn effective allee effect may induce the survival of low-density predator. Results Phys (2023) 53:106926.

  • Google Scholar

• 2.

  MayRStability and complexity in model ecosystems. Princeton: Princeton University Press (1973).

  • Google Scholar

• 3.

  SeoGDeAngelisDLA predator-prey model with a holling type i functional response including a predator mutual interference. J Nonlinear Sci (2011) 21:811–33. 10.1007/s00332-011-9101-6

  • CrossRef
  • Google Scholar

• 4.

  HsuSHuangTGlobal stability for a class of predator-prey system. SIAM J Appl Math (1995) 6:763–83.

  • Google Scholar

• 5.

  HuangJRuanSSongJBifurcations in a predator-prey system of leslie type with generalized holling type iii functional response. J Differ Equations (2014) 257:1721–52. 10.1016/j.jde.2014.04.024

  • CrossRef
  • Google Scholar

• 6.

  LiuCChangLHuangYWangZTuring patterns in a predator-prey model on complex networks. Nonlinear Dyn (2020) 99:3313–22. 10.1007/s11071-019-05460-1

  • CrossRef
  • Google Scholar

• 7.

  HuangJXiaoDAnalyses of bifurcations and stability in a predator-prey system with holling type-iv functional response. Acta Math Appl Sin Engl Ser (2004) 20:167–78.

  • Google Scholar

• 8.

  LiYXiaoDBifurcations of a predator-prey system of holling and lesile types. Chaos Soliton Fract (2007) 34:606–20.

  • Google Scholar

• 9.

  BeraSMaitiASamantaGStochastic analysis of a prey-predator model with herd behaviour of prey. Nonlinear Anal.-Model. (2016) 21:345–61. 10.15388/na.2016.3.4

  • CrossRef
  • Google Scholar

• 10.

  BrazaAPredator-prey dynamics with square root functional responses. Nonlinear Anal.-Real (2012) 13:1837–43. 10.1016/j.nonrwa.2011.12.014

  • CrossRef
  • Google Scholar

• 11.

  HuDCaoHStability and bifurcation analysis in a predator–prey system with Michaelis–Menten type predator harvesting. Nonlinear Anal.-Real (2017) 33:58–82. 10.1016/j.nonrwa.2016.05.010

  • CrossRef
  • Google Scholar

• 12.

  GuptaRChandraPBifurcation analysis of modified leslie-gower predator-prey model with michaelis-menten type prey harvesting. J Math Anal Appl (2003) 398:278–95. 10.1016/j.jmaa.2012.08.057

  • CrossRef
  • Google Scholar

• 13.

  DasTMukherjeeRChaudhariKBioeconomic harvesting of a prey-predator fishery. J Biol Dyn (2009) 3:447–62. 10.1080/17513750802560346

  • CrossRef
  • Google Scholar

• 14.

  ChakrabortyKJanaSKarTGlobal dynamics and bifurcation in a stage structured prey-predator fishery model with harvesting. Appl Math Comput (2012) 218:9271–90. 10.1016/j.amc.2012.03.005

  • CrossRef
  • Google Scholar

• 15.

  Rojas-PalmaAGonz’alez-OlivaresEOptimal harvesting in a predator-prey model with allee effect and sigmoid functional response. Appl Math Model (2012) 36:1864–74. 10.1016/j.apm.2011.07.081

  • CrossRef
  • Google Scholar

• 16.

  XiaoDRuanSBogdanov-takens bifurcations in predator-prey systems with constant rate harvesting. Fields Inst Commun (1999) 21:493–506. 10.1090/fic/021/41

  • CrossRef
  • Google Scholar

• 17.

  PengGJiangYLiCBifurcations of a holling-type ii predator-prey system with constant rate harvesting. Int J Bifurcat Chaos (2009) 19:2499–514. 10.1142/s021812740902427x

  • CrossRef
  • Google Scholar

• 18.

  HuangJGongYChenJMultiple bifurcations in a predator-prey system of holling and leslie type with constant-yield prey harvesting. Int J Bifurcat Chaos (2013) 23:1350164. 10.1142/s0218127413501642

  • CrossRef
  • Google Scholar

• 19.

  ZhuCLanKPhase portraits, hopf-bifurcations and limit cycles of leslie-gower predator-prey systems with harvesting rates. Discrete Contin Dyn Syst Ser B (2010) 14:289–306. 10.3934/dcdsb.2010.14.289

  • CrossRef
  • Google Scholar

• 20.

  XiaoDJenningsLSBifurcations of a ratio-dependent predator-prey system with constant rate harvesting. SIAM J Appl Math (2005) 65:737–53. 10.1137/s0036139903428719

  • CrossRef
  • Google Scholar

• 21.

  XiaoDLiWHanMDynamics in a ratio-dependent predator-prey model with predator harvesting. J Math Anal Appl (2006) 324:14–29. 10.1016/j.jmaa.2005.11.048

  • CrossRef
  • Google Scholar

• 22.

  LuoJZhaoYStability and bifurcation analysis in a predator-prey system with constant harvesting and prey group defense. Int J Bifurcat Chaos (2017) 27:1750179. 10.1142/s0218127417501796

  • CrossRef
  • Google Scholar

• 23.

  ZhangZDingTHuangWDongZQualitative theory of differential equation. Beijing: Science Press (1992).

  • Google Scholar

• 24.

  PerkoLDifferential equations and dynamical systems. 3rd ed.New York: Springer-Verlag (2000).

  • Google Scholar

• 25.

  ChenJHuangJRuanSWangJBifurcations of invariant tori in predator-prey models with seasonal prey harvesting. SIAM J Appl Math (2013) 73:1876–905. 10.1137/120895858

  • CrossRef
  • Google Scholar

• 26.

  HuangJGongYRuanSBifurcation analysis in a predator-prey model with constant-yield predator harvesting. Discrete Contin Dynam Syst Ser B (2013) 18:2101–21. 10.3934/dcdsb.2013.18.2101

  • CrossRef
  • Google Scholar

• 27.

  LiHCaoHFengYLiXPeiJOptimization of graph clustering inspired by dynamic belief systems. IEEE Trans Knowledge Data Eng (2023) 1–14. 10.1109/TKDE.2023.3274547

  • CrossRef
  • Google Scholar

• 28.

  LiHFengYXiaCCaoJOverlapping graph clustering in attributed networks via generalized cluster potential game. ACM T Knowl Discov D (2023) 18:1–26. 10.1145/3597436

  • CrossRef
  • Google Scholar

• 29.

  HouL-FSunG-QPercMThe impact of heterogeneous human activity on vegetation patterns in arid environments. Communical Nonlinear SCI (2023) 126:107461. 10.1016/j.cnsns.2023.107461

  • CrossRef
  • Google Scholar

• 30.

  FuSLuoJZhaoYStability and bifurcations analysis in an ecoepidemic system with prey group defense and two infectious routes. Math Comput Simulat (2023) 205:581–99.

  • Google Scholar

• 31.

  GaoSChangLPercMWangZTuring patterns in simplicial complexes. Phys Rev E (2023) 107:014216. 10.1103/physreve.107.014216

  • CrossRef
  • Google Scholar

• 32.

  WangCChangLLiuHSpatial patterns of a predator-prey system of leslie type with time delay. PLoS One (2016) 11:e0150503. 10.1371/journal.pone.0150503

  • CrossRef
  • Google Scholar

• 33.

  LiangJSunG-QEffect of nonlocal delay with strong kernel on vegetation pattern. J Appl Anal Comput (2024) 14:473–505. 10.11948/20230290

  • CrossRef
  • Google Scholar

• 34.

  YangJGongMSunG-QAsymptotical profiles of an age-structured foot-and-mouth disease with nonlocal diffusion on a spatially heterogeneous environment. J Differ Equations (2023) 377:71–112. 10.1016/j.jde.2023.09.001

  • CrossRef
  • Google Scholar

• 35.

  LiuS-MBaiZSunG-QGlobal dynamics of a reaction-diffusion brucellosis model with spatiotemporal heterogeneity and nonlocal delay. Nonlinearity (2023) 36:5699–730. 10.1088/1361-6544/acf6a5

  • CrossRef
  • Google Scholar

• 36.

  HeRLuoXAsamoahJKKZhangYLiYJinZet alA hierarchical intervention scheme based on epidemic severity in a community network. J Math Biol (2023) 87:29. 10.1007/s00285-023-01964-y

  • CrossRef
  • Google Scholar