Complex dynamics of a predator-prey fishery model: The impact of the Allee effect and bilateral intervention

1.   Introduction

Predator-prey interactions have been an interesting and challenging issue that is frequently discussed in marine ecosystems, especially in fish populations. Predator-prey interactions are the most important component of ecology, determining various factors such as community composition, species behavior and dynamics. Mathematical modeling helps to provide insights into the dynamics of the system, which was investigated in the early pioneering work of Lotka^[1] and Volterra^[2]. In dynamical systems, continuity, equilibrium stability, bifurcation, and control problems are also often studied^[3,4,5,6,7]. Traditionally, predators could be distinguished as specialists or generalists based on whether they ate only one or more types of prey. In general, a predator-prey model can be described by the following ordinary differential equation for both specialists and generalist predators^[8]:

where $P$ and $N$ denote the prey and predator's population densities at moment $t$, respectively. $F(\cdot)$ and $H(\cdot)$ represent the growth of the species in the absence of the other one. $G(\cdot)$ is known as the functional response, which characterizes the average individual of prey consumed by a predator; $\gamma$ represents the conversion rate from prey consumption to predator, $\upsilon(\cdot)$ is a monotonically increasing function. For specialist predators, there is $H(N) = -d < 0$; and for generalist predators, it is required that $\gamma\upsilon(G(0, N))+H(N) > 0$. In literature, different type of functional responses were adopted to model different species, which can be prey dependent^[9,10] or prey and predator dependent^{[11,12,13,14,15]}.

The classical view of population dynamics claims that the higher the population density, the lower the overall growth rate due to the competition for resources. The lower the density, the higher the overall growth rate. However, when the population density is low, Allee^[16] introduced the opposite view that the lower the population density, the lower the overall growth rate, that is, the Allee effect. The Allee effect is a common phenomenon in marine populations^[17,18,19]. When the population density is low, it may affect population development due to pairing restrictions, dispersal, habitat changes, cooperative foraging, cooperative defense, and predator saturation. Therefore, the study of the Allee effect on ecosystems has attracted the attention of many scholars. In general, the Allee effect can be represented by a multiplier of the form $P$-$L$^[20,21,22], where $L$ is the threshold for the Allee effect. When $L < 0$, it is a weak Allee effect, and the Allee effect is always positive no matter how much the prey growth rate decreases. When $L > 0$, it is a strong Allee effect. The strong Allee effect indicates that in order for the population to grow, the population size or density must be higher than $L$; otherwise, the population will die out. Scholars have analyzed the dynamics of the system with the Allee effect of the form $P$-$L$, and discussed the existence of the equilibrium of the system and various bifurcation phenomena such as saddle-node bifurcation and B-T bifurcation^[23,24,25].

Fish is a kind of important ecological resources. In view of the fish resources development problems, scholars studied population behavior by adding harvest items on the basis of continuous systems. However, fishing activities are not continuous, so continuous dynamic systems cannot accurately describe the actual fishing process. In the process of fish harvesting, the periodic harvesting of fish is a kind of human-controlled behavior that can be described by an impulse differential equation, and it has been found that the impulse differential equation is more accurate in describing and portraying the dynamical behavior of the population^[26]. The theory of semi-continuous dynamic systems has been widely used in modeling research on pest management modeling^{[27,28,29,30,31]}. Based on the analysis of the literature, due to the fact that impulsive differential systems (semi-continuous dynamical systems) have the characteristics of both continuous and discrete dynamical systems, there are some studies applying the theory to the development of fish populations in deterministic environments^{[32,33,34,35,36]} and uncertain environment^[37,38,39]. In addition, most of the studies considered fish harvesting or fish protection only (unilateral control); in this study, a bilateral state control^[40,41,42] is considered, and a predator-prey model for conservation and harvesting of two fish species is developed by constructing a semi-continuous dynamical system. When designing the state feedback control strategy, the number of objective fish was used as the state variable for feedback control. On the one hand, when the number of prey fish is small, the Allee effect will lead to their extinction, which will lead to the lack of enough food for prey fish and destroy the balance of the ecosystem. Therefore, it is necessary to release a certain amount of prey fish when the number of prey fish decreases to a certain threshold. When the number of prey fish is high, it is necessary to catch prey fish from the economic point of view. Since the fishing behavior will also lead to the harvesting of some predator fish, in order to maintain the ecological balance and avoid the extinction of predator fish caused by the fishing behavior, it is necessary to release a certain amount of predator fish larvae at the same time. Based on the above two aspects, we propose a bilateral control strategy to maintain the population size of the two species in a suitable range.

This paper considers a predator-prey model in which the predator is a generalist and the growth of the prey is affected by the Allee effect. The structure is as follows: In Section Ⅱ, we describe the fishery model, non-negativity, persistent survivability and discuss the existence and stability of equilibria of the model. In Section Ⅲ, the bifurcation dynamics of the model are discussed using bifurcation theory. In Section Ⅳ, based on this model, we analyze the model of marine fish harvesting and conservation with bilateral controls. The existence and stability of the order-1 and order-2 periodic solutions of the system are analyzed by using the geometry theory of differential equations. In the fifth section, we performed numerical simulations using MATLAB to verify the correctness of the results.

2.   Mathematical model and basic knowledge

2.1.   Fishery model with Allee effect

In this paper, we present a predator-prey model in which predators are generalists, prey growth rates are logical and subject to strong Allee effects, and the functional response is a Holling-Ⅰ type, while the conversion from prey consumption to predator species is saturated,

where $p(T)$, $N(T)$ denote the prey and predator's densities at the moment of $T$; $K$ denotes the prey's environmental holding capacity; $L$ is the threshold of the prey's Allee effect; $r$ and $s$ represent the intrinsic growth rates of prey and predator, respectively; $A$ is the capture rate of prey by the predator; $e$ is the efficiency with which prey is converted to predator; $B$ is the half-saturation constant; $f$ is the intensity of predator density dependence, and $d$ is the predator mortality rate. Since the predators are generalist, it requires that $eA\leq Bd$, $s > d$, and all parameters of model (2.1) are positive.

To facilitate the analysis, let ${x = \frac{P}{K}}$, ${y = N}$, ${t = rKT}$, ${\alpha = \frac{A}{{rK}}}$, ${\beta = \frac{L}{K}}$, ${\gamma = \frac{{eA}}{r}}$, ${\delta = BK}$, ${{s_1} = \frac{s}{{rK}}}$, ${{d_1} = \frac{d}{{rK}}}$. Then system (2.1) is simplified to system (2.2):

and from the biological point of consideration, the model (2.2) is limited in the region

On the other hand, as a renewable resource, fish species are closely related to human life. In order to maintain the balance of the ecosystem during the fishing process, we consider a fish stock control strategy with a combination of fishing and investment. First, to avoid the distinction of prey fish caused by the Allee effect, a quantity ($\eta$) of juvenile prey fish is released when the prey density declines to the level $x = {h_1}$. But at a higher level $x = {h_2}$, a proportion $a$ of prey fish together with a proportion $b$ of predator fish will also be caught for economic purposes, and simultaneously, a quantity $\tau$ of juvenile predator fish is released into the system to maintain the level of predator fish. Based on the control measures, the model can be described as follows:

where $\eta$, $a$, $b$, $\tau$ are all positive, and $a, b\in (0, 1)$.

2.2.   Impulsive semi-continuous system

For a given planar model

Definition 1 (Order-$k$ periodic solution^[32,33,36]) The solution $\tilde{\mathbf{z}}(t) = (\tilde{x}(t), \tilde{y}(t))$ is called periodic if there exists $m (\geqslant 1)$ satisfying $\tilde{\mathbf{z}}_{m} = \tilde{\mathbf{z}}_{0}$. Furthermore, $\tilde{\mathbf{z}}$ is an order-$k$ T-periodic solution with $k\triangleq\min\{j|1\leq j\leq m, \tilde{\mathbf{z}}_{j} = \tilde{\mathbf{z}}_{0}\}$.

Lemma 1 (Analogue of Poincaré Criterion^[32,33,36]). The order-$k$ $T$-periodic solution $\mathbf{z}(t) = (\xi(t), \eta(t))^T$ is orbitally asymptotically stable if $\vert \mu _q \vert < 1$, where

with

$f_1^+ = f_1(\xi(\theta_j^ +), \eta(\theta_j^ +))$, $f_2^+ = f_2(\xi(\theta_j^+), \eta(\theta_j^+))$ and $f_1$, $f_2$, $\frac{\partial {I_1} }{\partial x}$, $\frac{\partial {I_1} }{\partial y}$, $\frac{\partial {I_2} }{\partial x}$, $\frac{\partial {I_2} }{\partial y}$, $\frac{\partial \omega }{\partial x}$, $\frac{\partial \omega }{\partial y}$ are calculated at $(\xi(\theta_j), \eta(\theta_j))$.

3.   Dynamic properties of Model (2.2)

In this section, the bounded-ness of the solution for Model (2.2) is discussed. Moreover, the existence, type, and local stability of the equilibrium as well as the bifurcation properties are verified.

Define

3.1.   Positivity and bounded-ness of the solution

Theorem 1. The solution of Model (2.2) with non-negative initial values will remain non-negative for all time and is bounded on $\mathbb{R}_+^2$.

Proof. By Eq (2.2), it can be obtained that

for all $t\geq0$, as long as $x(0)$ and $y(0)$ are non-negative, then $x(t)$ and $y(t)$ are also non-negative.

Next, we define a function ${u\left(t \right) = \frac{\gamma }{\alpha }x\left(t \right) + y\left(t \right)}$. Then

which implies that

so long as $u_0 = \frac{\gamma}{\alpha}x_0 + y_0$ is bounded, $u(t)$ is bounded in $\Omega$. To sum up, any solution of Model (2.2) starting with a non-negative bounded initial condition is non-negative and bounded in $\Omega$. □

3.2.   Existence and stability of equilibrium

Obviously, four equilibria always exist

Define

and denote $\overline{\gamma}_1\triangleq {d_1}\delta$. Due to the assumptions $eA\leq Bd$ and $s > d$, then $\gamma < \overline{\gamma}_1$, i.e., $y_2(x) < \frac{{{s_1}\delta }}{{\overline{\gamma}_1-\gamma}}$. The positional relationship between ${y_1}(x)$ and ${y_2}(x)$ for different cases is shown in Figure 1.

Let $y_1(x) = y_2(x)$. Then it has

where

Define

Theorem 2. For any of the following cases: $C_1)$ $\alpha < \overline{\alpha}$, $\delta < \overline{\delta}_2$, ${\gamma} < \overline{\gamma}_1$; $C_2)$ $\alpha < \overline{\alpha}$, $\overline{\delta}_2 < \delta < \overline{\delta}_3$, $\overline{\gamma}_2 < \gamma < \overline{\gamma}_1$; $C_3)$ $\alpha > \overline{\alpha}$, $\delta < \overline{\delta}_2$, $0 < {\gamma} < -\overline{\gamma}_2$, if $\rho < 0$ holds, there exists two interior equilibria in Model (2.2); if $\rho = 0$, there exists a unique interior equilibrium in Model (2.2); if $\rho > 0$, there doesn't exists interior equilibrium in Model (2.2).

Proof. Clearly, the existence of an interior equilibrium is equivalent to that of a positive root of Eq (3.3) in the interval (0, 1). Since

and

then for any of the cases $C_1)$–$C_3)$, there is ${a_3} < 0$. Moreover, $g'(0) < 0$, $g'(x_d) = 0$, $g''(x_d) > 0$ and

so that for $x \in (0, {x_d})$, $g'(x) < 0$, for $x \in({x_d}, 0)$, $g'(x) > 0$. If $\rho < 0$, Eq (3.3) has two distinct positive roots $x^*_i\in (0, 1)$, $i = 1, 2$. Denote $y^*_i = y_1(x^*_i)$, $i = 1, 2$. Then two interior equilibria exist in Model (2.2), denoted by ${E^*_{1}}\left({{x^*_{1}}, {y^*_{1}}} \right)$ and ${E^*_{2}}\left({{x^*_{2}}, {y^*_{2}}} \right)$. If $\rho = 0$, Eq (3.3) has a unique positive root $x^* = x_d\in (0, 1)$, then a unique positive equilibrium exists in Model (2.2), denoted by ${E^*}(x^*, y_1(x^*)$; when $\rho > 0$, Eq (3.3) does not have positive root, and thus there does not exist interior equilibrium in Model (2.2). □

3.2.1.   Type and stability of equilibrium

For any equilibrium $\overline{E}(\overline{x}, \overline{y})$, there is

its characteristic equation is

where

1) Boundary equilibria

At $O(0, 0)$, there is

since $\lambda_1 = {{-\beta}} < 0$, $\lambda_2 = s_1-d_1 > 0$, then ${O}$ is an unstable higher-order singularity.

At ${E^1}\left(0, \frac{1}{f}\left({\frac{{{s_1} - {d_1}}}{{{d_1}}}} \right)\right)$, there is

since $\lambda_1 = 0$, ${\lambda _2} = - \frac{\alpha }{f}\left({\frac{{{s_1} - {d_1}}}{{{d_1}}}} \right) < 0$, so that ${E^1}$ is locally stable.

At ${E^2}(\beta, 0)$, there is

since $\lambda_1 = {{\beta(1-\beta)}} > 0$, ${\lambda _2} = \frac{{\gamma \beta }}{{1 + \delta \beta }} + {s_1} - {d_1} > 0$, then ${E^2}$ is unstable.

At ${E^3}(1, 0)$, there is

since $\lambda_1 = {{\beta-1}} < 0$, $\lambda _2 = \frac{\gamma }{{1 + \delta }} + {s_1} - {d_1} > 0$, then ${E^3}$ is unstable.

2) Interior equilibrium

At $E^*(x^*, y^*)$, there are $g_1(x^*, y^*) = 0$ and $g_2(x^*, y^*) = 0$, then

thus,

where $x^*y^*\frac{{\partial {g_1}}}{{\partial y}}\frac{{\partial {g_2}}}{{\partial y}}|_{x^*, y^*} > 0$, then the sign of $\text{DET}(J(E^*))$ is identical to that of $\frac{{\text{d}{y_1}}}{{\text{d}x}}|_{x = x^*} - \frac{{\text{d}{y_2}}}{{\text{d}x}}|_{x = x^*}$. Next, it will discuss the sign of $\frac{{\text{d}{y_1}}}{{\text{d}x}}|_{x = x^*} - \frac{{\text{d}{y_2}}}{{\text{d}x}}|_{x = x^*}$ for different cases in Theorem 2.

ⅰ) When $\rho < 0$ holds, two interior equilibria ${E^*_{1}}\left({{x^*_{1}}, {y^*_{1}}} \right)$ and ${E^*_{2}}\left({{x^*_{2}}, {y^*_{2}}} \right)$ with $0 < {x^*_{1}} < {x^*_{2}} < 1$ exists in Model (2.2), as illustrated in Figure 2(a). It can be easily checked that

Besides, at ${E^*_{{1}}}$, there is $\frac{{\text{d}{y_1}}}{{\text{d}x}}|_{x = x^*_1} > \frac{{\text{d}{y_2}}}{{\text{d}x}}|_{x = x^*_1}$, thus

i.e., ${E^*_{1}}$ is unstable. Similarly, at ${E^*_{2}}$, there is $\frac{{\text{d}{y_1}}}{{\text{d}x}}|_{x = x^*_2} < \frac{{\text{d}{y_2}}}{{\text{d}x}}|_{x = x^*_2}$, thus

i.e., ${E^*_{2}}$ is a locally asymptotically stable node.

ⅱ) When $\rho = 0$ holds, system (2.2) has a unique interior equilibrium ${E^*}\left({{x^*}, {y^*}} \right)$, as shown in Figure 2(b). Let $X = x - {x^*}, Y = y - {y^*}$, then ${E^*}$ is converted to the origin $O(0, 0)$, and the model is written as

where,

and ${P_3}(X, Y)$ is a function of $(X, Y)$ with degree of three or higher. The Jacobian matrix at ${E^*}$ is

thus

a) If $\text{TR}(J({E^*})) = {a_{11}} + {a_{22}} = 0$, then ${\lambda _1} = {\lambda _2} = 0$. Make the transformation

then Model (3.5) is converted into the following standard form:

where

and ${\bar P_3}(X, Y)$ is a function of three or more degrees about $(X, Y)$.

Let $\tau = {\bar a_{12}}t$. For convenience, $t$ is still used to represent $\tau$, then it has

where

and ${\tilde P_3}(X, Y)$ is a function of $(X, Y)$ with three or higher degree. Model (3.6) can be transformed into the following form^[43]:

and if ${\tilde B_1} \ne 0$, then ${E^*}$ is a cusp point.

Meanwhile, if ${\tilde B_2} + 2{\tilde A_1} \ne 0$, ${E^*}$ is a cusp of codimension two. If ${\tilde B_2} + 2{\tilde A_1} = 0$, ${E^*}$ is a cusp with at least codimension three.

b) If $\text{TR}(J({E^*})) = {a_{11}} + {a_{22}}\ne 0$, then ${\lambda _1} = 0$, ${\lambda _2} \ne 0$. Make the transformation

system (3.5) is converted into the following standard form

where

Next, introduce a new variable $\tau = {a'_{22}}t$ (for convenience, is represented by $t$), then

where, ${\hat A_i} = \frac{{{{A'}_i}}}{{{{a'}_{22}}}}, {\hat B_i} = \frac{{{{B'}_i}}}{{{{a'}_{22}}}}, i = 1, 2, 3$.

If ${\hat A_1} \ne 0$ then ${E^*}$ is unstable. Meanwhile, ${E^*}$ is a saddle node of attraction.

To sum up, the following result hold.

Theorem 3. For Model (2.2), 1) ${O}(0, 0)$ is unstable, ${E^1}(0, \frac{1}{f}\left({\frac{{{s_1} - {d_1}}}{{{d_1}}}} \right))$ is locally stable, ${E^2}(\beta, 0)$, ${E^3}(1, 0)$ is unstable; 2) For the interior equilibrium, when $\rho > 0$, ${E^*_{{1}}}\left({{x^*_{{1}}}, {y^*_{{1}}}} \right)$ is a saddle (unstable), and ${E^*_{{2}}}\left({{x^*_{{2}}}, {y^*_{{2}}}} \right)$ is a stable node; when $\rho = 0$, if $\text{TR}(J({E^*})) = 0$ and ${\tilde B_1} \ne 0$, ${E^*}\left({{x^*}, {y^*}} \right)$ is a cusp of codimension two in case of ${\tilde B_2} + 2{\tilde A_1} \ne 0$, and a cusp with at least codimension three in case of ${\tilde B_2} + 2{\tilde A_1} = 0$; If $\text{TR}(J({E^*})) \ne 0$ and ${\hat A_1} \ne 0$, ${E^*}\left({{x^*}, {y^*}} \right)$ is an attractive saddle node.

4.   Parameter's influence

4.1.   Saddle-node bifurcation

Let $\alpha = {\alpha _0}$ satisfy

Denote

and define

Theorem 4. Let the parameters of Model (2.2) satisfy $\text{TR}(J({E^*})) \ne 0$ and ${\hat A_1} \ne 0$. If ${\Phi_2} \ne 0$, system (2.2) undergoes a saddle node bifurcation near ${E^*}\left({{x^*}, {y^*}} \right)$ when $\alpha = {\alpha _0}$.

Proof. At ${E^*}\left({{x^*}, {y^*}} \right)$, there is

Let $V = {\left({{V_1}, {V_2}} \right)^T}$ ($W = {\left({{W_1}, {W_2}} \right)^T}$) be the eigenvector corresponding to the zero eigenvalue of $J\left({{E^*}} \right)$ (${\left({J({E^*})} \right)^T}$). Then $V = {\left({{\omega_1}, - {\xi_1}} \right)^T}$, $W = {\left(- {\xi_2}, {{\xi_1}} \right)^T}$. Let $g = {\left({{g_1}, {g_2}} \right)^T}$. Then

and

Since

and

then according to Sotomayor's theorem^[44], Model (2.2) undergoes a saddle-node bifurcation near ${E^*}\left({{x^*}, {y^*}} \right)$ when $\alpha = {\alpha_0}$. □

Remark 1. For Model (2.2), it can be concluded that for a lager $\alpha$, the interior equilibrium does not exist. When $\alpha$ decreases to $\alpha = \alpha_0$, a unique equilibrium exists in the system. As $\alpha$ decreases below $\alpha_0$, the system undergoes a saddle-node bifurcation at the interior equilibrium ${E^*}$, giving rise to two interior equilibria ${E^*_{{1}}}$ and ${E^*_{{2}}}$.

4.2.   Bogdanov-Takens bifurcation

According to Theorem 3, when $\rho = 0$, $\text{TR}(J({E^*})) = 0$, ${\tilde B_1} \ne 0$ and ${\tilde B_2} + 2{\tilde A_1} \ne 0$, Model (2.2) to undergo a Bogdanov-Takens bifurcation of codimension two near ${E^*}$ when $\left({\alpha, \beta } \right) = \left({{\alpha _0}, {\beta _0}} \right)$. Next, it will show the universal unfolding of the Bodmanov-Takens bifurcation of codimension two under parameter perturbation when $\alpha$ and $\beta$ are taken as bifurcation parameters.

Define

Let $\left({{\epsilon_1}, {\epsilon_2}} \right)$ be a parameter vector in a small neighborhood of $(0, 0)$. Then

Theorem 5. Let the parameters of Model (2.2) satisfy $\rho = 0$, $\text{TR}(J({E^*})) = 0$, ${\tilde B_1} \ne 0$, ${\tilde B_2} + 2{\tilde A_1} \ne 0$ and $\left| {\frac{{\partial \left({{O_1}, {O_2}} \right)}}{{\partial \left({{\epsilon_1}, {\epsilon_2}} \right)}}} \right| \ne 0$. When $\left({\alpha, \beta } \right)$ varies in the neighborhood of $\left({{\alpha _0}, {\beta _0}} \right)$, Model (2.2) changes in the small neighborhood of $E^*({x^*}, {y^*})$, and a codimensional $2$ Bogdanov-Takens branching occurs.

Proof. Consider the perturbation system

For $\left({\alpha, \beta } \right) = \left({{\alpha _0}, {\beta _0}} \right)$, there is $\text{DET}(J({E^*})) = 0, \text{TR}(J({E^*})) = 0$. With the transformation ${x_1} = x - {x^*}$, ${y_1} = y - {y^*}$, we can obtain

where ${N_3}({x_1}, {y_1}, {\epsilon_1}, {\epsilon_2})\in{C^\infty }$.

Make the following transformation:

then Model (4.1) is converted to

where ${N'_3}({x_2}, {y_2}, {\epsilon_1}, {\epsilon_2})\in{C^\infty }$ with coefficients depending smoothly on ${\epsilon_1}, {\epsilon_2}$.

Next, introduce the variable $\tau$, denoted as $\text{d}t = \left({1-{H_6}{x_2}} \right)\text{d}\tau$ (still denote $\tau$ as $t$), then

Let ${x_3} = {x_2}, {y_3} = \left({1 - {H_6}{x_2}} \right){y_2}$, then the above system of equations is transformed into

where ${\tilde N_3}({x_3}, {y_3}, {\epsilon_1}, {\epsilon_2})\in{C^\infty }$ with coefficients depending smoothly on ${\epsilon_1}, {\epsilon_2}$.

(ⅰ) When ${J_4} < 0$, the following transformations are applied to the variables:

still denote $\tau$ as $t$, there is

where ${M_3}({x_4}, {y_4}, {\epsilon_1}, {\epsilon_2})\in{C^\infty }$ with at least third order.

Let ${x_5} = {x_4} - \frac{{K_2}}{2}$, ${y_5} = {y_4}$, and obtain

where ${M'_3}({x_5}, {y_5}, {\epsilon_1}, {\epsilon_2})\in{C^\infty }$ with at least third order.

If set ${J_5} \ne 0$, then ${L_3} \ne 0$. Define new variables: ${x_6} = - L_3^2{x_5}$, ${y_6} = L_3^3{y_5}$, $\tau = - \frac{t}{{L{}_3}}$, and denote ${x_6}$ by $x$, ${y_6}$ by $y$, and $\tau$ by $t$, which yields that

where ${\tilde N_3}({x_2}, {y_2}, {\epsilon_1}, {\epsilon_2})\in{C^\infty }$ with at least third order, and ${O_1}$, ${O_2}$ can be represented by $\epsilon_1$ and $\epsilon_2$.

(ⅱ) When ${J_4} > 0$, the following transformations are applied

still denote $\tau$ by $t$, it has

where ${N_3}({x_4}, {y_4}, {\epsilon_1}, {\epsilon_2})\in{C^\infty }$ with at least third order, and

Let ${x_5} = {x_4} + {{{{K'_2}}}}/{2}, {y_5} = {y_4}$. Then it has

where ${N'_3}({x_5}, {y_5}, {\epsilon_1}, {\epsilon_2})\in{C^\infty }$ with at least third order, and

If set ${J_5} \ne 0$, then ${L'_3} \ne 0$. Define new variables: ${x_6} = {L'_3}^2{x_5}$, ${y_6} = {L'_3}^3{y_5}$, $\tau = {t}/{{{{L'_3}}}}$, and still denote ${x_6}$ by $x$, ${y_6}$ by $y$, $\tau$ by $t$, which yields that

where ${O'_1} = {L'_1}{L'_3}^4, {O'_2} = {L'_2}{L'_3}$, and ${\epsilon_1}$, ${\epsilon_2}$ can be represented by ${O'_1}$, ${O'_2}$.

For convenience of discussion, ${O'_1}$, ${O'_2}$ is still denoted by ${O_1}, {O_2}$. When $\left| {\frac{{\partial \left({{O_1}, {O_2}} \right)}}{{\partial \left({{\epsilon_1}, {\epsilon_2}} \right)}}} \right| \ne 0$, Models (4.2) and (4.3) are the cardinal folds of the Bogdanov-Takens bifurcation^[44], when $\left({\alpha, \beta } \right)$ varies in the vicinity of $\left({{\alpha_0}, {\beta_0}} \right)$, Model (2.2) undergoes a codimension 2 Bogdanov-Takens bifurcation in a small neighborhood of ${E^*}\left({{x^*}, {y^*}} \right)$. □

5.   Complex dynamics of control-Model (2.3)

It only focuses on the case of $\rho < 0$ in Theorem 3. For Model (2.3), there are

5.1.   Order-1 periodic solution

Denote ${l_1} = \frac{1}{\alpha }\left({1 - x} \right)\left({x - \beta } \right)$ as the prey isocline, ${l_2} = \frac{{\left({{s_1} - {d_1}} \right)\left({1 + \delta x} \right) - \gamma {x}}}{{f\left[{{d_1} + \left({{d_1}\delta - \gamma } \right){x}} \right]}}$ as the predator's isocline, ${l_3}$, ${l_4}$, ${l_5}$, ${l_6}$ as the saddle point separatrix of $E^*_1$ in different directions. The intersection point between $l_1$ and $\mathcal{N}_2$ is denoted by $A_3$. The intersection point between $l_4$ and $\mathcal{N}_2$ is denoted by $M$. The trajectory from $A_3$ intersects $\mathcal{M}_2$ at the point ${B_3}$. Define $\overline{\tau}_2\triangleq y_{M}-(1-b)y_{B_3}$. Denote

Definition 2 (Successor function). For a point $S\in \mathcal{N}_1$, if the trajectory from $S$ intersects $\mathcal{M}_1$, then denote the intersection point by $S^{-}\in \mathcal{M}_1$. Under the impulse effect, the point $S^{-}$ is mapped to $S^{+}\in \mathcal{N}_1$. In such a case, we can define $f^{\text{I}}_{\text{SOR}_1}$: $\mathcal{N}_1\rightarrow R$, $S\rightarrow f^{\text{I}}_{\text{SOR}_1}(S)\triangleq y_{S^{+}}-y_{S}$. If the trajectory from $S$ intersects $\mathcal{M}_2$, then denote the intersection point by $S^{-}$. Under the impulse effect, the point $S^{-}$ is mapped to $S^{+}\in \mathcal{N}_2$. If the trajectory from $S^{+}$ intersects $\mathcal{M}_1$, denote the intersection point by $S^{+-}\in \mathcal{M}_1$. Under the impulse effect, the point $S^{+-}$ is mapped to $S^{++}\in \mathcal{N}_1$. In such cases, we can define $f^{\text{II}}_{\text{SOR}_1}$: $\mathcal{N}_1\rightarrow R$, $S\rightarrow f^{\text{II}}_{\text{SOR}_1}(S)\triangleq y_{S^{++}}-y_{S}$. Similarly, For a point $S'\in \mathcal{N}_2$, define $f^{\text{I}}_{\text{SOR}_2}$: $\mathcal{N}_2\rightarrow R$, $S'\rightarrow f^{\text{I}}_{\text{SOR}_2}(S')\triangleq y_{S'^{+}}-y_{S'}$.

Theorem 6. For system (2.3) with model parameters satisfying any one of $C_1)$–$C_3)$ and $\rho > 0$, if D-1) ${h_1} + \eta < {x^*_{1}} < \left({1 - a} \right){h_2} < {h_2} < {x^*_{2}}$, $y_2(h_1)\geq y_1(h_1+\eta)$ and $\tau < \overline{\tau}_2$, an order-1 periodic solution exists in each $\Omega_i$, $i = 1, 2$; if D-2) ${h_1} < \left({1 - a} \right){h_2} < {h_1} + \eta < {x^*_{1}} < {h_2} < {x^*_{2}}$ and $y_2(h_1)\geq y_1(h_1+\eta)$, an order-1 periodic solution exists in $\Omega_1$; if D-3) ${h_1} < {x^*_{1}} < {h_1} + \eta < \left({1 - a} \right){h_2} < {h_2} < {x^*_{2}}$ and $\tau < \overline{\tau}_2$, an order-1 periodic solution exists in $\Omega_2$.

Proof. For case D-1) ${h_1} + \eta < {x^*_{1}} < \left({1 - a} \right){h_2} < {h_2} < {x^*_{2}}$ and $y_2(h_1)\geq y_1(h_1+\eta)$, denote the intersection point between $l_2$ and $\mathcal{N}_1$ by $A_0$. Select a point $A_1\in \mathcal{N}_1$ above $A_0$; the trajectory from $A_1$ intersects $\mathcal{M}_1$ at $B_1$. Under the impulse effect, the point $B_1$ is mapped to $A^{+}_1\in \mathcal{N}_1$. Since $g_1(x_{A_1}, y_{A_1}) < 0$, $g_2(x_{A_1}, y_{A_1}) < 0$, then we have $f^{\text{I}}_{\text{SOR}_1}(A_1) = y_{A^{+}_1}-y_{A_1} = y_{B_1}-y_{A_1} < 0$. Besides, denote the intersection between $l_1$ and $\mathcal{N}_1$ by $A_2$. Since $g_1(x_{A_2}, y_{A_2}) = 0$, $g_2(x_{A_2}, y_{A_2}) > 0$ and $g_1(x_S, y_S) < 0$ for $S\in U(A_1, \epsilon)$ with $x_S < h_1+\eta$ and $y_S < y_{A_1}$, then we have $f^{\text{I}}_{\text{SOR}_1}(A_2) = y_{A^{+}_2}-y_{A_2} = y_{B_2}-y_{A_2} > 0$ due to $y_2(h_1)\geq y_1(h_1+\eta)$. The continuity of $f^{\text{I}}_{\text{SOR}_1}$ implies that a point $S\in \overline{A_1A_2}$ exists so that $f^{\text{I}}_{\text{SOR}_1}(S) = 0$, i.e., an order-1 periodic solution exists in $\Omega_1$ (Figure 3(a)). Similarly, it can be proved that an order-1 periodic solution exists in $\Omega_1$ for case D-2).

Under the impulse effect, ${B_3}$ is mapped to ${A_3}^{+}\in{\mathcal{N}_2}$, where ${y_{{A_3}^ + }} = \left({1 - b} \right){y_{{B_3}}} + \tau$. Define ${\overline{\tau}_1} \triangleq {y_{{A_3}}} - \left({1 - b} \right){y_{{B_3}}}$.

ⅰ) If $\tau = \overline{\tau}_1$, then $f^{\text{I}}_{\text{SOR}_2}({A_3}) = {y_{{A_3}^ + }} - {y_{{A_3}}} = 0$, i.e., an order-1 periodic solution exists in $\Omega_1$ (Figure 3(a)).

ⅱ) If $\tau < \overline{\tau}_1$, then $f^{\text{I}}_{\text{SOR}_2}({A_3}) = {y_{{A_3}^ + }} - {y_{{A_3}}} < 0$. On the other hand, for ${A_4}((1-a)h_2, \tau)$, there is $f^{\text{I}}_{\text{SOR}_2}({A_4}) = {y_{{A_4}^ + }} - {y_{{A_4}}} > 0$. The continuity of $f^{\text{I}}_{\text{SOR}_2}$ implies that a point $S'\in \overline{A_3A_4}$ exists so that $f^{\text{I}}_{\text{SOR}_2}(S') = 0$, i.e., an order-1 periodic solution exists in $\Omega_2$ (Figure 3(b)).

ⅲ) If $\overline{\tau}_1 < \tau < \overline{\tau}_2$, then $f^{\text{I}}_{\text{SOR}_2}({A_3}) = {y_{{A_3}^ + }} - {y_{{A_3}}} > 0$. According to the trend of the trajectory and the fact that any two trajectories cannot be intersected, select $D\in \mathcal{N}_2$ above and sufficiently close to the ${A_3}$, then $f^{\text{I}}_{\text{SOR}_2}(D) = {y_{{D^ + }}} - {y_D} > 0$. On the other hand, for $C = {A^{+}_3}$, there is $f^{\text{I}}_{\text{SOR}_2}(C) = {y_{{C^ + }}} - {y_C} < {y_{{A^{+}_{3}}}} - {y_C} = 0$. The continuity of $f^{\text{I}}_{\text{SOR}_2}$ implies that a point $S'\in \overline{CD}$ exists so that $f^{\text{I}}_{\text{SOR}_2}(S') = 0$, i.e., an order-1 periodic solution exists in $\Omega_2$ (Figure 3(c)).

Similarly, system (2.3) has an order-1 periodic solution in $\Omega_2$ for case $D$-$3$).

To sum up, an order-1 periodic solution exists in Model (2.3) for case D-1). □

Let $\mathbf{z}_i(t) = (\phi_i(t), \varphi_i(t))$ $(k-1)T_i\leq t\leq kT_i$, $k\in \mathbb{N}$ be the order-1 periodic solution in $\Omega_i$, $i = 1, 2$. For $\mathbf{z}_1(t) = (\phi_1(t), \varphi_1(t))$ $(k-1)T_1\leq t\leq kT_1$, denote

For $\mathbf{z}_2(t) = (\phi_2(t), \varphi_2(t))$ $(k-1)T_2\leq t\leq kT_2$, denote

Theorem 7. For the model parameters with any one of $C_1)$–$C_3)$, $\rho > 0$ and D-1), $\mathbf{z}_i(t) = (\phi_i(t), \varphi_i(t))$ $(k-1)T_i\leq t\leq kT_i$ is orbitally asymptotically stable if $\mu_i < 1$, where

with $\Lambda\left(t \right) = \exp \left\{ {\int_{{0^ + }}^{{T_2}} {\left[{\phi_2(t)\left({\beta + 1 - 2\phi_2(t)} \right) -\frac{{{s_1}f\varphi_2(t)}}{{{{\left({1 + f\varphi_2(t)} \right)}^2}}}} \right]} dt} \right\}$.

Proof. For system (2.3), there are

and ${\omega_1}\left({x, y} \right) = x - {h_1}, {I_{11}}\left({x, y} \right) = \eta, {I_{21}}\left({x, y} \right) = 0$, So it can be concluded that

Denote ${f_{{1^ + }}}\triangleq {f_1}\left({\phi_1\left({T_1^ + } \right), \varphi_1\left({T_1^ + } \right)} \right), {f_{{2^ + }}}\triangleq {f_2}\left({\phi_1\left({T_1^ + } \right), \varphi_1 \left({T_1^ + } \right)} \right)$. Then by Lemma 1, there are

and

If $\mu_1 < 1$, the order-1 periodic solution $\mathbf{z}_1(t) = (\phi_1(t), \varphi_1(t))$ $(k-1)T_1\leq t\leq kT_1$ is orbitally asymptotically stable.

Similarly, for the order-1 periodic solution $\mathbf{z}_2(t) = (\phi_2(t), \varphi_2(t))$ $(k-1)T_2\leq t\leq kT_2$, there is

where $\Lambda\left(t \right) \triangleq \exp \left\{ {\int_{{0^ + }}^{{T_2}} {\left[{\phi_2(t)\left({\beta + 1 - 2\phi_2(t)} \right) - \frac{{{s_1}f\varphi_2(t)}}{{{{\left({1 + f\varphi_2(t)} \right)}^2}}}} \right]} dt} \right\}$. By Lemma 1, if $\mu_2 < 1$, the order-1 periodic solution $\mathbf{z}_2(t) = (\phi_2(t), \varphi_2(t))$ $(k-1)T_2\leq t\leq kT_2$ is orbitally asymptotically stable. The proof is completed. □

5.2.   Order-2 periodic solution

Here only the case for ${h_1} < \left({1 - a} \right){h_2} < {x^*_{1}} < {h_1} + \eta < {h_2} < {x^*_{2}}$ is considered. Let ${G_1}\left({{h_1} + \eta, {y_{{G_1}}}} \right)$ be the intersection point between ${l_1}$ and ${\mathcal{N}_1}$, ${G_2}\left({{h_1}, {y_{{G_2}}}} \right)$ be the intersection point between ${l_2}$ and ${\mathcal{M}_1}$. The trajectory that starts with ${G_1}$ intersects ${\mathcal{M}_2}$ at ${B_1}\left({{h_2}, {y_{{B_1}}}} \right)$. Define $\overline{\tau}_3\triangleq y_{G_1}-(1-b)y_{B_1}$.

Theorem 8. For the model parameters with any one of $C_1)$-$C_3)$, $\rho > 0$ and D-4) ${h_1} < \left({1 - a} \right){h_2} < {x^*_{1}} < {h_1} + \eta < {h_2} < {x^*_{2}}$, if $by_{G_2}\leq \tau\leq \min\{y_{G_2}, \overline{\tau}_3\}$ holds, an order-2 periodic solution exists in system (2.3).

Proof. For ${G_1}\in \mathcal{N}_1$, the trajectory starting from ${G_1}$ intersects with ${\mathcal{M}_2}$ at ${B_1}$. Then ${B_1}$ is mapped to ${B''_1}\in{\mathcal{N}_2}$, and next intersects with $\mathcal{M}_1$ at $\hat{B}_1$, and then $\hat{B}_1$ is mapped to $G^{+}_1$ under impulse effect. Since $\tau < \overline{\tau}_3$, then ${y_{B''_1}} = \left({1 - b} \right){y_{{B_1}}} + \tau < {y_{{G_1}}}$. Since $g_1(B''_1) < 0$ and $g_2(B''_1) < 0$, it has ${y_{{{\hat B}_1}}} < {y_{{{B''}_1}}}$, then ${y_{G_1^ + }} = y_{\hat{B}_1} < y_{B''_1} = \left({1 - b} \right){y_{{B_1}}} + \tau < {y_{{G_1}}}$, i.e., $f^{\text{II}}_{\text{SOR}1}({G_1}) = {y_{G_1^ + }} - {y_{{G_1}}} < 0$.

On the other hand, since $\tau\leq y_{G_2}$, then $G_2\in \mathcal{N}_1$. The trajectory starting from $A_1\in \mathcal{N}_1$ with $y_{A_1} = y_{G_2}$ intersects with $\mathcal{M}_2$ at $A^{-}_1$, and then $A^{-}_1$ is mapped to $A^{+}_1$. Then it intersects with $\mathcal{M}_1$ at $A^{+-}_1$, and next $A^{+-}_1$ is mapped to $A^{++}_1$. Since $y_{A^{+}_1} = (1-b)y_{A^{-}_1}+\tau > (1-b)y_{G_2}+\tau > y_{G_2}$, so $f^{\text{II}}_{\text{SOR}1}({A_1}) = y_{A^{++}_1} - y_{A_1} = y_{A^{+-}_1}-y_{G_2} > 0$.

The continuity of $f^{\text{II}}_{\text{SOR}_1}$ implies that a point $S\in \overline{A_1G_1}\subset \mathcal{N}_1$ exists so that $f^{\text{II}}_{\text{SOR}_1}(S) = 0$, i.e., an order-2 periodic solution exists (Figure 4(a)).

Similarly, a point $S'\in \overline{A_2B2}\subset \mathcal{N}_2$ exists so that $f^{\text{II}}_{\text{SOR}_2}(S') = 0$, i.e., an order-2 periodic solution exists (Figure 4(b)). □

Let $\mathbf{z}_{3}(t) = (\phi_3(t), \varphi_3(t))$ $(k-1)(T_1+T_2)\leq t\leq k(T_1+T_2)$, $k\in \mathbb{N}$ be the order-2 periodic solution. Denote

Define

where

Theorem 9. For the model parameters with any one of $C_1)$–$C_3)$, $\rho > 0$ and D-4) ${h_1} < \left({1 - a} \right){h_2} < {x^*_{1}} < {h_1} + \eta < {h_2} < {x^*_{2}}$, and $by_{G_2}\leq \tau\leq \min\{y_{G_2}, \overline{\tau}_3\}$, if

then $\mathbf{z}_{3}(t) = (\phi_3(t), \varphi_3(t))$ $(k-1)(T_1+T_2)\leq t\leq k(T_1+T_2)$ is orbitally asymptotically stable.

Proof. For convenience, denote the intersection point between $\mathbf{z}_{3}(t) = (\phi_3(t), \varphi_3(t))$ and $\mathcal{N}_1$ ($\mathcal{M}_2$, $\mathcal{N}_2$, $\mathcal{M}_1$) by $L_1(h_1+\eta, \delta_3)$ ($L_2(h_2, \delta_4)$, $L_3((1-a)h_2, (1-b)\delta_4+\tau)$, $L_4(h_1, \delta_3)$). Then, according to analogue of Poincaré Criterion, there are

and

Then

If $\mu < 1$, then $\mathbf{z}_{3}(t) = (\phi_3(t), \varphi_3(t))$ $(k-1)(T_1+T_2)\leq t\leq k(T_1+T_2)$ is orbitally asymptotically stable. □

6.   Numerical simulations and verification

6.1.   Numerical simulations for system (2.2)

For system (2.2) with model parameters $\alpha = 0.09$, $\beta = 0.15$, $\gamma = 0.6$, $\delta = 0.9$, ${s_1} = 1.12$, ${d_1} = 1.1$, $f = 0.21$, there is $\rho < 0$, then two interior equilibria exist in the system, and the phase diagram is presented in Figure 5(a); For the model parameters $\alpha = 0.219$, $\beta = 0.051$, $\gamma = 0.6$, $\delta = 0.9$, ${s_1} = 1.12$, ${d_1} = 1.1$, $f = 0.21$, there is $\rho = 0$, then a unique equilibrium exists in system (2.2), which is a sharp point. The phase diagram of the system (2.2) for such a case is presented in Figure 5(b). Moreover, system (2.2) undergoes a Bogdanov-Takens bifurcation of codimension two in a very small neighborhood of the unique interior equilibrium.

While for system (2.2) with model parameters $\alpha = 0.09$, $\beta = 0.29$, $\gamma = 0.6$, $\delta = 0.9$, ${s_1} = 1.12$, ${d_1} = 1.1$, $f = 0.21$, there is also $\rho = 0$, in such case a unique equilibrium exists in system (2.2), which is a saddle node, and the phase diagram is presented in Figure 5(c). System (2.2) undergoes a saddle-node bifurcation of codimension one in a very small neighborhood of the interior equilibrium.

For system (2.2) with model parameters $\alpha = 0.09$, $\beta = 0.35, \gamma = 0.6$, $\delta = 0.9$, ${s_1} = 1.12$, ${d_1} = 1.1$, $f = 0.21$, there is $\rho > 0$. In such a case, system (2.2) does not have interior equilibrium, and the phase diagram of the system (2.2) is presented in Figure 5(d).

For system (2.2) with parameters $\beta = 0.051$, $\gamma = 0.6$, $\delta = 0.9$, ${s_1} = 1.12$, ${d_1} = 1.1$, $f = 0.21$, the bifurcation diagrams of the residual dimension 1 for the system (2.2) are shown in Figure 6(a), (b) when $\alpha$ is selected as the bifurcation parameter. The result shows that for larger $\alpha$, the interior equilibrium does not exist. When $\alpha$ decreases to $\alpha = \alpha_0$, a unique equilibrium exists in the system. As $\alpha$ decreases below $\alpha_0$, the system undergoes a saddle-node bifurcation at the interior equilibrium ${E^*}$, giving rise to two interior equilibria ${E^*_{{1}}}$ and ${E^*_{{2}}}$.

6.2.   Simulations for system (2.3)

For system (2.3) with model parameters $\alpha = 0.1$, $\beta = 0.15$, $\gamma = 0.6$, $\delta = 0.9$, ${s_1} = 1.12$, ${d_1} = 1.1$, $f = 0.21$, the control parameters: ${h_1} = 0.07$, $\eta = 0.1$, ${h_2} = 0.6$, $a = 0.35$, $b = 0.38$, $\tau = 0.51$, an order-1 periodic solution can be formed in both $\Omega_1$ and $\Omega_2$ (Theorem 6), as presented in Figure 7(a); For model parameters $\alpha = 0.1$, $\beta = 0.15$, $\gamma = 0.6$, $\delta = 0.9$, ${s_1} = 1.12$, ${d_1} = 1.1$, $f = 0.21$ and control parameters ${h_1} = 0.07$, $\eta = 0.128$, ${h_2} = 0.44$, $a = 0.7$, $b = 0.08$, $\tau = 0.51$, an order-1 periodic solution can be formed in $\Omega_1$ (Theorem 6), as presented in Figure 7(b); while for control parameters ${h_1} = 0.14$, $\eta = 0.16$, ${h_2} = 0.6$, $a = 0.25$, $b = 0.38$, $\tau = 0.51$, an order-1 periodic solution can be formed (Theorem 6), as presented in Figure 7(c).

For system (2.3) with given model parameters $\alpha = 0.09$, $\beta = 0.215$, $\gamma = 0.6$, $\delta = 0.9$, ${s_1} = 1.12$, ${d_1} = 1.1$, $f = 0.21$ and control parameters ${h_1} = 0.07$, $\eta = 0.1$, ${h_2} = 0.4$, $a = 0.65$, $b = 2$, $\tau = 3$, an order-2 periodic solution can be formed in system (2.3) (Theorem 3.2), as presented in Figure 8(a); while for the control parameters ${h_1} = 0.07$, $\eta = 0.4$, ${h_2} = 0.68$, $a = 0.65$, $b = 0.85$, $\tau = 0.5$, a different order-2 periodic solution can be formed in system (2.3) (Theorem 8), as presented in Figure 8(b).

7.   Conclusions

Considering that the Allee effect is an important mechanism in ecosystems and a realistic description of the interaction between species, we presented a model of prey-predator system with prey's Allee effect and generalist predator in the context of fish resources (Models (2.1) or (2.2)). We investigated the dynamic properties of Model (2.2) such as the type and stability of the boundary equilibria as well as the existence and stability of the interior equilibrium in detail (Theorems 1–3, Figure 5).

To show the influence of the parameters on the dynamics of Model (2.2), we analyzed the bifurcations in the predation system by selecting the capture rate of prey by predator and Allee threshold as key parameters. We showed that Model (2.2) will undergo a saddle-node bifurcation as changing of the capture rate $\alpha$ (Figure 6), and undergo a Bogdanov-Takens bifurcation of codimension at least 2 and 3 as changing of $(\alpha, \beta)$.

To achieve sustainable and efficient exploitation of fish stocks, we adopted a bilateral intervention strategy, i.e., to avoid the distinction of prey fish caused by the Allee effect, releasing juvenile prey fish is adopted at a lower-level $x = h_1$; while for economic purposes, capturing both prey and predator fish is adopted at a higher-level $x = h_2$. We obtained the conditions for the existence and stability of the order-1 periodic solution (Theorems 6, 7, Figure 7) and order-2 periodic solution (Theorems 8, 9, Figure 8) of the control system (2.3). The results showed that the extinction can be prevented by control even when the prey density is low, while in the case of the prey density increasing to a certain extent, fishing activities can be taken in a periodic way (periodic solution) to obtain the fish resources. Therefore, as long as the fish stocks are properly managed, the number of fish stocks can be controlled within an appropriate range, and the sustainable development and utilization of biological resources can be realized.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The work was supported by the National Natural Science Foundation of China (No. 11401068).

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.