Dynamics of a Gilpin-Ayala predator-prey system with state feedback weighted harvest strategy

1.   Introduction

Biological populations in nature exhibit intricate interactions, including competition, cooperation and predation. The predator-prey relationship is one of the most prevalent ecological relationships in the natural world. Both predator and prey play pivotal roles as vital biological resources, impacting human life and economic growth. To preserve the renewal potential of these biological resources, rational development and scientific management are imperative^[12,16,24]. Employing limited renewable resources for sustainable development and utilization has emerged as a shared interest among economic management scientists, mathematicians and ecologists ^[4,5].

Italian mathematician V. Volterra introduced the Lotka-Volterra model in 1926 to depict variations in fish population. This model still serves as a foundational reference for predator-prey interaction models today ^{[18,19,28,34]}. However, the classic Lotka-Volterra model overlooks the limited capacity for predation in its portrayal of predator-prey relationships, thereby oversimplifying the functional response of the predator to the prey. Holling addressed this oversight in 1965 by introducing three distinct functional response functions ^[11]. Subsequent studies, such as that by M. Zhang and L. Chen ^[32], delved into a state feedback impulsive controlled predator-prey system with the Holling-Ⅱ functional response. They examined its equilibrium and presented the state feedback model, which assumes the harvest of predators and supplementation of preys as impulsive disturbances. This model aims to find the dynamical balance between prey and predator populations. Further research by M. Huang et al. ^[13,14,15] investigated a class of time dependent impulsive switching systems and discussed a series of population control strategies. Nevertheless, state feedback impulsive control proves to be more appropriate for studying predator-prey systems ^[29,33]. In contrast to time dependent harvesting strategies, state feedback harvesting strategies consider the current status of species^[25,30], and reduce potential adverse effects on species' sustainability. The weighted sum of harvesting methods, which take into account the states of both predator and prey populations ^[7,24], promote the sustained development of natural populations.

Biomathematical researchers have previously used logistic growth models to characterize predator-prey relationships. However, such models tend to depict population growth in a linear manner, which is not always accurate. To rectify this shortcoming, Gilpin-Ayala introduced the exponential parameter $\alpha$ to the classical logistic population growth term. In 1973, Ayala and Gilpin improved the logistic growth model, proposing the following model ^[1]:

By introducing a parameter $\alpha \in (0, 1)$, they transformed the model into a nonlinear system, thereby increasing its analytical complexity ^[22]. Nevertheless, this exponential parameter $\alpha$ overcomes the limitations of the classical logistic population growth theory, specifically the assumption that each addition to the population reduces the population growth rate by a constant. This refined model offers a more nuanced and accurate representation of various population dynamics. Building on this idea, Jiang et al. ^[17,26,27] further developed this concept by allowing the Gilpin-Ayala parameter to vary with state under Markov switching and pulse interference.

Gonz$\check{s}$lez-Olivares and Rojas-Palma ^[9] furthered this area of study by proposing a model incorporating the Holling Ⅲ functional response and the Allee effect:

Further research on this model was presented in ^[30] where the model construction incorporated nonlinear impulse disturbances. Both the predation mortality and predator release rates depend on their densities, adding realism to the model.

The Gilpin-Ayala growth model provides a more precise representation of population density changes. By introducing the exponential parameter, it overcomes the inherent limitations of the classic logistic population growth theory. This theory traditionally posits that each addition to the population uniformly reduces the growth rate. In reality, the Gilpin-Ayala model offers a more insightful description of the dynamic behaviors of population organisms and stands as one of the most critical numerical models in biological research ^[2]. Moreover, compared to the initial two Holling functional response functions ^[32], the Holling Ⅲ type ^[8,21] more authentically mirrors consumption rate variations at different resource densities. Based on the above literature analysis and research from ^[20], this paper formulates a predator-prey model that incorporates the Gilpin-Ayala growth, the Holling Ⅲ type functional response and state feedback impulsive control.

This paper is organized as follows: In Section 2 we develop a Gilpin-Ayala predator-prey system with a discontinuous weighted harvest strategy and provide the lemmas and definitions to be used in subsequent sections. Section 3 presents proof of positivity and boundedness for the uncontrolled system, as well as the global stability of equilibrium points. We also discuss the existence and stability for the predator extinction periodic solution and positive order-1 periodic solution for the controlled system. In Section 4 we provide some numerical simulations to validate these theoretical results, and, finally, a brief conclusion is given in Section 5.

2.   Model formulation and preliminaries

In this paper, we consider a sufficiently large prey population and thus neglect the Allee effect, such as in ^[3,9,20]. Consequently, we introduce a predator-prey model featuring both Gilpin-Ayala growth and Holling Ⅲ functional response as follows:

where $A(t)$ and $B(t)$ denote the densities of prey and predator at time $t$, respectively. The variable $r$ stands for the intrinsic growth rate of the prey population, while $K$ represents the maximum environmental carrying capacity of the prey population in the absence of predators and harvesting. The symbol $\alpha$ is a positive parameter to modify the classical logistic model and represents a nonlinear measure of interspecific interference. The term $\frac{\tau_{1}A^{2}}{A^{2}+\beta}$ denotes the functional response of the predator, commonly referred to as the Holling type Ⅲ response function. The variables $\tau_{1}$, $\tau_{2}$, $\beta$ and $\omega$ correspond to the predator capture rate, the conversion coefficient specifying the ratio of newborn predators to captured prey, the half-saturation constant, and the death rate of the predator, respectively.

Considering commercial interests and the sustainability of biological resources, reasonable harvesting is permitted. The authors in ^[30] focused only on monitoring the predator population state to trigger pulse harvesting, examining a nonlinear ecological control model with a complex discrete map. Similar harvesting mechanisms have been very common in previous studies ^{[1,17,22,26,27,30]}. In ^[24], a different state dependent harvesting method was proposed that permitted the harvesting only when the combined densities of both prey and predator, weighted by specific factors, reached a predefined threshold. This approach not only ensures the sustainable development of the two biological resources but also safeguards commercial interests. Drawing on their approach, we also consider a weighted state pulse harvest for the above system (2.1). Denote the predefined harvesting threshold for the weighted sum of the densities of prey and predator as $I$. The intensity of this harvesting is measured by $h$ and the prey and predator populations are assigned weights of $\mu$ and $1-\mu$, respectively. The model driven by this weighted harvesting strategy can be expressed as follows:

where $c_{1}$ and $c_{2}$ represent capture rates, while $\Lambda$ describes the constant number of predator released at time $t$.

The practice of harvesting two types of populations simultaneously and releasing a certain number of predators in model (2.2) is also in line with many practical application scenarios. For example, in fishing operations, when the combined count of phytoplankton or shrimp (representing the prey) and fish (acting as the predator) hits a certain limit, there's an initiative to harvest the fish. However, to optimize salvage costs and uphold ecological equilibrium, phytoplankton is simultaneously harvested. This not only boosts the economic returns, but also enhances fish yield per unit area. Subsequently, to guarantee the yield and maintain a balanced ecosystem, a specific quantity of fish fry is introduced into the environment. This approach mirrors the prevalent fishery practice known as "rotational catch and stocking."

The foundational concepts and theories related to state feedback pulse dynamic systems are important for our subsequent discussion. Descriptions of such dynamic systems, along with definitions of pulse set, phase set, successor function, order $k$-periodic solution and its stability, are commonly found in various literature. Given this prevalence, we will exclude certain rudimentary explanations here (interested readers can refer to ^[6]). Instead, we'll focus on detailing the definition for a distinct type of successor function and the criteria to determine the order-k periodic solution.

Definition 1. (Successor function^[25,31]) Assume that the pulse set and the phase set of a state feedback pulse dynamic system are two lines denoted by $\bar{M}$ and $\bar{N}$, respectively (see Figure 1). Let $O_{0}$ be the intersection point of $\bar{N}$ and the A-axis. For any point $V_{1}\in\bar{N}$, we denote the trajectory that originates from $V_{1}$ by $\Gamma_{1}(V_{1})$. If the trajectory $\Gamma_{1}(V_{1})$ intersects with the phase set $\bar{N}$ and pulse set $\bar{M}$ at points $V_{2}$ and $V_{3}$, and point $V_{3}$ is mapped to a point $V_{4}\in \bar{N}$ through pulse effect, then $S_{1}(V_{2}) = d(V_{4}, O_{0})-d(V_{2}, O_{0})$ is called the type-1 successor function of point $V_{1}$, while $S_{2}(V_{1}) = d(V_{4}, O_{0})-d(V_{1}, O_{0})$ is called the type-2 successor function of point $V_{1}$, where the function $d(\cdot, \cdot)$ represents the directed distance between two points.

Lemma 1. (Lasalle invariance principle) Consider the following system of differential equations

where $(A, B)\in G\subset\mathbb{R}^{2}$ with $G$ is a bounded set. $P(A, B)$ and $Q(A, B)$ are continuous and satisfy the Lipschitz condition. Given a scalar function $V(A, B)$ that is continuous over $\bar{G}$ and possesses a continuous first-order partial derivative, $V(A, B)$ is referred to as a Lyapunov function on $G$ if

Define

Assuming $M\subset S$ is the largest invariant subset within $S$, and given a solution $x(t, t_{0}, (A_{0}, B_{0}))\subset G$ of the system, it follows that $x(t, ,t_{0}, (A_{0}, B_{0}))\rightarrow M$ as $t\rightarrow \infty$.

Lemma 2. (Analogue of the Poincare criterion^[10,23,25]) The $T$-periodic solution $\tilde{\rho}_{L}(t) = (\tilde{A}(t), \tilde{B}(t))$ of the system

is orbitally asymptotically stable and exhibits asymptotic phase behavior if the Floquet multiplier $\mu_{2}$ satisfies $|\mu_{2}| < 1$. Here,

and

where $P^{+} = P(\tilde{A}(\tau_{k}^{+}), \tilde{N}(\tau_{k}^{+}))$ and $Q^{+} = Q(\tilde{A}(\tau_{k}^{+}), \tilde{N}(\tau_{k}^{+}))$. The values of $P, \; Q, \; \frac{\partial a}{\partial A}, \frac{\partial a}{\partial B}, \frac{\partial b}{\partial A}, \frac{\partial b}{\partial B}, \frac{\partial\phi}{\partial A}, \frac{\partial\phi}{\partial B}$ are computed at the point $(\tilde{A}(\tau_{k}), \tilde{B}(\tau_{k}))$.

3.   Main works

3.1.   Positive and boundedness of solutions for system (2.1)

Let's suppose that $A(0) = A_{0} > 0$ and $B(0) = B_{0} > 0$ are positive initial conditions. Initially, we assess the system (2.1) with these conditions for the positivity and boundedness of the solutions. Considering that the right-hand sides of (2.1) are smooth functions with respect to both $A$ and $B$ and $\frac{dA}{dt}|_{A = 0} = \frac{dB}{dt}|_{B = 0} = 0$, the state space for the system (2.1) is the positive quadrant $\{(A, B):A > 0, B > 0\}$, which is an invariant set. Consequently, the solutions of the system (2.1) are positive. Next, we verify the boundedness of the solutions for the system (2.1).

Undoubtedly,

Let

where $\bar{f}(A)\geq0$ for $A\in(0, K]$, and $\bar{f}(A) < 0$ for $A > K$. Thus, we obtain $A(t)\leq\max\{A(0), K\}$.

Subsequently, we denote $N: = \frac{\tau_{2}}{\tau_{1}}A+B$, then we have

Through direct calculation, we can deduce that the maximum point of $\tilde{f}(A)$ is $A_{m} = K\sqrt[^{\alpha}]{\frac{\omega+r}{r(\alpha+1)}}$. Denote the maximum value by $M_{m} = \tilde{f}(A_{m})$, we have

and

Define

and we can easily know that $\Omega_{0}$ is a positive invariant set for the system (2.1).

3.2.   Equilibria and stability analysis of system (2.1)

Model (2.1) always has two equilibria: a saddle point $P_{o}(0, 0)$ and a boundary equilibrium $P_{k}(K, 0)$. Moreover, a positive equilibrium $P^{*}(A^{*}, B^{*})$ exists when $\frac{\tau_{2}K^{2}}{K^{2}+\beta} > \omega$, where

To examine the local stability of the equilibrium points $P_{k}(K, 0)$ and $P^{*}(A^{*}, B^{*})$, we first determine the Jacobian matrix for the system (2.1) at any point $(A, B)$. It is given by

The stability of the equilibrium points will be determined by the eigenvalues of the Jacobian matrix calculated at the respective equilibrium points.

By a simple calculation, the trivial equilibrium point $P_{o}(0, 0)$ is unstable as long as it exists. Substituting $(A, B) = (K, 0)$ into (3.2), we get the characteristic values $\lambda_{k1} = -r\alpha$, $\lambda_{k2} = \frac{\tau_{2}K^{2}}{K^{2}+\beta}-\omega$. Clearly, if $\frac{\tau_{2}K^{2}}{K^{2}+\beta} < \omega$, then $P^{*}$ does not exist and $P_{k}(K, 0)$ is locally stable; if $\frac{\tau_{2}K^{2}}{K^{2}+\beta} > \omega$, i.e., $P^{*}$ exists, then $P_{k}$ is unstable and both $P_{k}$ and $P_{o}$ are saddle points.

From the above discussion, both $P_{o}(0, 0)$ and $P_{k}(K, 0)$ are unstable when $P^{*}$ exists. To investigate the local asymptotic stability of the positive equilibrium point $P^{*}(A^{*}, B^{*})$, we calculate

from which we obtain

The trace of $P^{*}$ is given by

This implies that $P^{*}$ is locally asymptotically stable if the following condition holds

In order to illustrate the global dynamic behavior of $P_{k}(K, 0)$ and $P^{*}(A^{*}, B^{*})$, we construct two Lyapunov functions in the following two theorems.

Theorem 1. The boundary equilibrium point $P_{k}(K, 0)$ of system (2.1) is globally asymptotically stable when $\frac{\tau_{2}K^{2}}{K^{2}+\beta} < \omega$.

Proof. Consider a function $V(A, B)$ defined by

By taking the derivative of system (2.1), it can be obtained that

Since $\frac{\tau_{2}K^{2}}{K^{2}+\beta} < \omega$ and $A\leq K$, we have $(\tau_{2}-\omega)K^{2}-\beta\omega < 0$ and $(\tau_{2}-\omega)A^{2}-\beta\omega < 0$. This leads to the conclusion that $\frac{dV}{dt}\leq0$. Moreover, $\frac{dV}{dt} = 0$ if, and only if, $(A, B) = (K, 0)$. By the Lasalle invariance principle (Lemma 1), the boundary equilibrium point $P_{k}(K, 0)$ of system (2.1) is globally asymptotically stable. The proof is completed. □

The preceding discussion establishes that a positive equilibrium $P^{*}$ exists if $\frac{\tau_{2}K^{2}}{K^{2}+\beta} > \omega$. To assess the global stability of $P^{*}$, we let

Theorem 2. The positive equilibrium point $P^{*}(A^{*}, B^{*})$ of system (2.1) is globally asymptotically stable when $\frac{\tau_{2}K^{2}}{K^{2}+\beta} > \omega$, if one of the following conditions holds:

(Ⅰ) $\xi_{1}(A) = 0$ has no or only one positive root.

(Ⅱ) $\xi_{1}(A) = 0$ has two positive roots $A_{1}$ and $A_{2}$ (where $A_{1} < A_{2}$). If $A^{*} \in (0, A_{1}')$, then $F(A) > B^{*}$ for all $A \in (0, A^{*})$ and $F(A) < B^{*}$ for all $A \in (A^{*}, K)$ holds; if $A^{*} \in (A_{2}', K)$, then $F(A) > B^{*}$ for all $A \in(0, A^{*})$ and $F(A) < B^{*}$ for all $A \in (A^{*}, K)$ holds, where $A_{1}', A_{2}'$ are two points that must exist in this case such that $F(A'_{1}) = F(A_{2})$ and $F(A'_{2}) = F(A_{1})$ with the ordering $0<A'_{1}<A_{1}<A_{2}<A'_{2}<K$.

Proof. Let us define a function

then we have

Clearly, $\frac{dV}{dt} = 0$ if, and only if, $(A, B) = (A^{*}, B^{*})$.

Next, we consider the function

and compute its derivative with respect to $A$, we have

Let

then it follows, by differentiate with respect to $A$, that

In the special case of $\alpha = 1$, we obtain the maximum point $\xi_{1\max}(A) = \frac{K^{2}}{27}-\beta$ of $\xi_{1}(A) = 0$. For the roots of $\xi_{1}(A) = 0$ we note that: $(i)$ If $\frac{K^{2}}{27} < \beta$, then there is no root. $(ii)$ If $\frac{K^{2}}{27} = \beta$, then there is a unique root $A = A_{c}$ satisfying $\sqrt{\beta} < A_{c} < K$. $(iii)$ If $\frac{K^{2}}{27} > \beta$, then there are two roots $A_{1}, A_{2}$ satisfying $\sqrt{\beta} < A_{1} < A_{2} < K$.

For the case $\alpha\in(0, 1)$, when $\beta(1-\alpha)-A^{2}(1+\alpha) = 0$, i.e., $A = \sqrt{\frac{\beta(1-\alpha)}{1+\alpha}}$, we have $\xi_{1}(A) = A^{2}-\beta < 0$. For $A\neq\sqrt{\frac{\beta(1-\alpha)}{1+\alpha}}$ and $\xi_{1}(A) = 0$, we have

When $0 < \alpha < 1$, we observe that $\xi_{2}(0) = \frac{1}{1-\alpha} > 1, \lim_{A\rightarrow \{\sqrt{\frac{1-\alpha}{1+\alpha}\beta}\}^- }\xi_{2}(A) = +\infty,$ $\lim_{A\rightarrow \{\sqrt{\frac{1-\alpha}{1+\alpha}\beta}\}^+ }\xi_{2}(A) = -\infty, 0 < \xi_{2}(\sqrt{\beta}) = 0, \xi_{2}(K) < 1$ and $y_1(A) = \xi_{2}(A)$ is strictly monotonically increasing on intervals $[0, \sqrt{\frac{1-\alpha}{1+\alpha}\beta})$ and $(\sqrt{\frac{1-\alpha}{1+\alpha}\beta}, K]$. While $y_2(A) = \left(\frac{A}{K}\right)^{\alpha}$ is strictly monotonically increasing on interval $[0, K]$ with $y_2(0) = 0$ and $y_2(K) = 1$. The graphs of the functions $y_2(A) = \xi_{2}(A)$ and $y_1(A) = \left(\frac{A}{K}\right)^{\alpha}$ either have no intersection points, or they intersect at a single point $A_{c}$ such that $\sqrt{\beta} < A_{c} < K$, or they intersect at two points, $A_{1}$ and $A_{2}$, where $\sqrt{\beta} < A_{1} < A_{2} < K$.

Using an analogous analysis for $\alpha > 1$, the distribution of the roots for $\xi_{1}(A) = 0$ exhibits a similar behavior. Therefore, for any $\alpha > 0$, there are three cases for the root of the function $F'(A) = 0$. Let's summarize these aspects in two cases:

$(I)$ For the cases (i) and (ii), we obtain that $F(A)-B^{*} > 0$ when $A < A^{*}$ and $F(A)-B^{*} < 0$ when $A > A^{*}$. This implies that $\frac{dV}{dt}\leq 0$ for all $(A, B)\in \Omega_0$, and $\frac{dV}{dt} = 0$ if, and only if, $(A, B) = (A^{*}, B^{*})$ (see Figure 2(a) and (b)).

$(II)$ For the case (iii), we identify two positive points, $A_{1}$ and $A_{2}$, such that ${F'}(A_{i}) = 0$ for $i = 1, 2$. Additionally, two distinct points, $A'_{1}$ and $A'_{2}$, must exist satisfying $F(A'_{1}) = F(A_{2})$ and $F(A'_{2}) = F(A_{1})$, with the ordering $0<A'_{1}<A_{1}<A_{2}<A'_{2}<K$ (see Figure 2(c)). Based on this, we deduce the following:

(a) $F(A)>B^{*}$ for $A \in (0, A^{*})$ and $F(A)<B^{*}$ for $A \in (A^{*}, K)$ when $A^{*} \in (0, A'_{1})$;

(b) $F(A)>B^{*}$ for $A \in (0, A^{*})$ and $F(A)<B^{*}$ for $A \in (A^{*}, K)$ when $A^{*} \in (A'_{2}, K)$.

Thus, we can also claim that $\frac{dV}{dt}\leq 0$ for all $(A, B)\in \Omega_0$, and $\frac{dV}{dt} = 0$ if, and only if, $(A, B) = (A^{*}, B^{*})$ According to Lasalle invariance principle (Lemma 1), we affirm that the positive equilibrium point $P^{*}(A^{*}, B^{*})$ of system (2.1) is globally asymptotically stable if one of the conditions listed in Theorem 2 holds. The proof is completed. □

3.3.   Complex dynamic behavior of impulsive system (2.2)

To discuss the dynamic behavior of system (2.2), we first define the impulse set and phase set as follows:

A trajectory of system (2.2) that is tangent to the phase set at $(A, B)\in \bar{N}$ must satisfy the following conditions:

Next, we will concentrate on the dynamic behavior of both predator-extinction periodic solution and positive periodic solution of the system (2.2).

For $\Lambda = 0$, we obtain that $B\equiv 0$ if $B(0) = 0$. At this point system (2.2) can be transformed into the following form:

Assuming that $\bar{N}$ intersects the $A$-axis at the point $L(A_{0}, 0)$, we denote $A = \tilde{A}(t)$ as a solution of the equation $A'(t) = rA(1-(A/K)^{\alpha})$ with the initial condition $A(0) = A_{0}$. We denote a solution of the system (3.7) as $\tilde{\rho}_{L}(t) = (\tilde{A}(t, A_{0}), 0)$, where

Let $\tilde{T} = \int^{I/\mu}_{A_{0}}K^{\alpha}/rA(K^{\alpha}-A^{\alpha})dA = \tilde{T}(A_{0})$, which satisfies $\tilde{A}(\tilde{T}) = I/\mu$ and $\tilde{A}(\tilde{T}^{+}) = A_{0}$. This illustrates that $A(t) = \tilde{A}(t, A_{0})$ is a periodic solution of (3.7) with period $\tilde{T}$.

Theorem Under the conditions $\Lambda = 0$ and $I\leq\mu K$, if $c_{2} > \tilde{c}_{2}$, then the system (2.2) admits an orbitally asymptotically stable predator-extinction periodic solution $\tilde{\rho}_{L}(t) = (\tilde{A}(t, A_{0}), 0)$. Here, $\tilde{c}_{2}$ is defined as

Proof. Assume $I > \mu K$, we notice that the trajectory of system (3.7) will not reach $A = I/\mu$ for any initial value $(A(t_{0}), 0)$. Therefore, there is no periodic solution in the system (3.7) and the predator-extinction equilibrium $(K, 0)$ is asymptotically stable. Conversely, if $I/\mu \leq K$ and $\Lambda = 0$, according to Definition 1, we find that $S_{1}(L) = d(L, L)-d(L, L) = 0$, and $\tilde{\rho}_{L}(t) = (\tilde{A}(t, A_{0}), 0)$ is a predator-extinction periodic solution of system (2.2).

According to Lemma 2, we have

By direct calculations, we get

Whereupon, we can obtain

Moreover,

where $\tilde{A} = \tilde{A}(t, A_{0})$, and then

Thus, when $c_{2} > \tilde{c}_{2}$, we have $\mu_{2} < 1$ and the periodic solution $\tilde{\rho}_{L}(t) = (\tilde{A}(t, A_{0}), 0)$ is orbitally asymptotically stable. On the other hand, if $0 < c_{2} < \tilde{c}_{2}$, we find that $\mu_{2} > 1$, implying the instability of the periodic solution $\tilde{\rho}_{L}(t) = (\tilde{A}(t, A_{0}), 0)$. The proof is completed. □

Assume $H$ is a point on the phase set with ${\frac{dB}{dA}}|_{H} = \frac{-\mu(1-c_{2}h)}{(1-c_{1}h)(1-\mu)}$. The trajectory of the system (2.2) starting from $H$ intersects the isoclinic line $f_{1}(A, B) = 0$ and the impulse set at points $H'$ and $H''$, respectively. Define some thresholds for the system (2.2) as follows:

Theorem 4. System (2.2) has a positive order-1 periodic solution, if one of the conditions $(i)$–$(iii)$ holds: $(i)\; \Lambda = 0, I\leq \min\{I^{*}, \mu K\}$, $0 < c_{2} < \tilde{c}_{2}$; $(ii)\; \Lambda = 0, \mu < \mu^{*}, \mu K < I < I_{m}$; $(iii)\; 0 < \Lambda < \Lambda_{m}, I\leq I_{m}$, where $\Lambda_{m} = \min \{\frac{c_2h}{1-\mu}, \frac{c_1h(1-c_2h)\Lambda}{(1-\mu)(1-c_1h)}\}$.

Proof. $(i)$ Referring to Theorem 1, the boundary periodic solution of (2.2) is unstable under this condition. For any $P_{1}\in {\bar{N}}\cap U(L, \delta)$, where $\delta$ is sufficiently small, the trajectory originating from $P_{1}$ intersects the impulse set $\bar{M}$ at point $P_{1}^{-}$ and then is instantly pulsed to the point $P_{1}^{+}$ in the phase set $\bar{N}$. Here, $P_{1}^{+}$ is above $P_{1}$ as shown in Figure 3. Following Definition 1, we have

Assuming that the point $P_{2}(A_{2}, B_{2})\in \bar{N}$ satisfies function (3.6), the trajectory of system (2.2) starting from $P_{2}(A_{2}, B_{2})$ first intersects the impulse set $\bar{M}$ at $P_{2}^{-}(A_{2}^{-}, B_{2}^{-})$ and then pulses to the phase set $\bar{N}$ at $P_{2}^{+}$. There are three possible situations: $(1)P_{2}^{+} = P_{2}$, $(2)P_{2}^{+} = \underline{P}_{2}^{+}$ and $(3)P_{2}^{+} = \overline{P}_{2}^{+}$ (shown in Figure 3).

Firstly, for the case of $P_{2}^{+} = P_{2}$, by Definition 1 we have $S_{1}(P_{2}) = d(P_{2}^{+}, L)-d(P_{2}, L) = 0$. This implies the solution starting from $P_{2}$ forms an order-1 periodic solution of system (2.2). Suppose $P_{2}^{+} = \underline{P}_{2}^{+}$, then we obtain that $S_{1}(P_{2}) = d(\underline{P}_{2}^{+}, L)-d(P_{2}, L) < 0$. As we have already established that $S_{1}(P_{1}) = d(P_{1}^{+}, L)-d(P_{1}, L) > 0$, a point $E\in \overline{P_{1}P_{2}}$ must exist such that $S_{1}(E) = 0$. Consequently, there exists a positive order-1 periodic solution starting from point $E$.

In the scenario $P_{2}^{+} = \overline{P}_{2}^{+}$, we have $S_{1}(P_{2}) = d(\overline{P}_{2}^{+}, L)-d(P_{2}, L) > 0$. Since the trajectory of system (2.2) starting from $P_{2}(A_{2}, B_{2})$ is tangent to the phase set $\bar{N}$, the trajectory of system (2.2) starting from $P_{2}^{+}$ crosses $\bar{N}$ once more at $\overline{P}_{2}^{+}$ before reaching $\bar{M}$ at $P_{2}^{+-}$ and then pulsing back to the phase set at $P_{2}^{++}$. According to the properties of system (2.1), $P_{2}^{+-}$ is below $P_{1}^{-}$ and $P_{2}^{++}$ is below $P_{1}^{+}$, which leads to

In addition, due to the continuity of solution in the system (2.1) and for $\varepsilon = d(P_{2}, P_{2}^{+})/2$, we know that there must exist some $\delta < \varepsilon$ such that $D\in U(P_{2}, \delta)$, $d(D^{-}, P_{2}^{-}) < \varepsilon$, and $d(D^{+}, \overline{P}_{2}^{+})\leq \max\{(1-\tau_{1}), (1-\tau_{2})\}d(D^{-}, P_{2}^{-}) < 0$.

Consequently,

Thus, there must exist a point $E\in \overline{P_{2}\overline{P}_{2}^{+}}$ such that $S_{2}(E) = 0$, and the trajectory from $E$ forms a positive order-1 periodic solution.

$(ii)$ For this case, we can easily have that the positive equilibrium $(A^{*}, B^{*})$ is above the pulse line $\bar{M}$. Let the intersections of $\bar{M}$ with $f_{1}(A, B) = 0$ and $A$-axis be the points $Z$ and $Z_1$, respectively, where point $Z_1$ is on the right of point $(K, 0)$. Therefore, we can obtain

Because $\mu K < I < I_{m}$, the system doesn't yield an extinction periodic solution, and all orbits starting from the phase set $\bar{N}$ reach the pulse set $\bar{M}$. Similar to the proof in $(i)$, choosing $P_{1}\in \bar{N}$ with $B_{P_{1}} = (1-\tau_{2}h)B_{Z}$, we can confirm the existence of an order-1 periodic solution of system (2.2).

$(iii)$ Choose $P_{1}\in \bar{N}$ with $d(P_{1}, L) = \Lambda$, and similar to the discussion in case (i), we can obtain the existence of a positive order-1 periodic solution for (2.2). The proof is completed. □

In Theorem 4, we've already proven the existence of positive order-1 periodic solution of system (2.2). In the upcoming part, we aim to study the stability of the periodic solution. For simplicity, we introduce some notations in the following.

Let's assume a positive order-1 periodic solution of (2.2) as $\bar{\rho}_{E}(t) = (\bar{A}(t), \bar{B}(t))$ with period $\bar{T}_{1}$, and the orbit $\Gamma (E) = \widehat{EE^{-}E}$, where $E = (\bar{A}_{0}, \bar{B}_{0})\in \bar{N}$, $E^{-} = (\bar{A}_{1}, \bar{B}_{1})\in \bar{M}$, $\bar{A}_{1}\triangleq \bar{A}(\bar{T})$, $\bar{B}_{1}\triangleq \bar{B}(\bar{T})$, $\bar{A}_{0}\triangleq (1-c_{1}h)\bar{A_{1}}$, $\bar{B}_{0}\triangleq (1-c_{2}h)\bar{B}_{1}+\Lambda$. Furthermore, we mark $\bar{P}_{0}\triangleq P(\bar{A}_{0}, \bar{B}_{0})$, $\bar{P}_{1}\triangleq P(\bar{A}_{1}, \bar{B}_{1})$, $\bar{Q}_{0}\triangleq Q(\bar{A}_{0}, \bar{B}_{0})$, $\bar{Q}_{1}\triangleq Q(\bar{A}_{1}, \bar{B}_{1})$.

Denote

Theorem 5. The positive order-1 periodic solution $\bar{\rho}_{E}(t) = (\bar{A}(t), \bar{B}(t))$ of system (2.2) is orbitally asymptotically stable if

Proof. According to Lemma 2, we have

and

When condition (3.10) is satisfied, we can obtain

which implies that the positive order-1 periodic solution of system (2.2) is orbitally asymptotically stable. □

3.4.   Optimal harvest control strategy

The existence of a positive periodic solution for the system (2.2), as shown in Theorem 4, justifies the possibility of implementing a periodic harvesting strategy with a period of $\bar{T}$. Let's consider $A_{c}(A_{2} < A_{c} < A_{2}^{-})$ and $B_{c}: = \varrho A_{c}$ to represent the target densities of prey and predator, respectively. As $h$ relies on $A_{c}$, and both $h$ and $\Lambda$ depend on $\mu$ and $A_{c}$, we assume that

We can also infer that $\bar{T}$ depends on $\mu$ and $A_{c}$. Let $n_{1}$ and $n_{2}$ denote the sale prices of unit prey and predator, respectively, while $n_{3}$ signifies the cost per unit capture strength and $n_{4}$ represents the cost of feeding predators per unit.

To assure the sustainable renewal of resources while maximizing commercial profits, it is essential to optimize the harvesting process. The profit from harvesting can be formulated as

The primary objective here is to optimize the profit cycle, i.e.,

where $A_{2} < A_{c} < A_{2}^{-}$ and $0\leq\mu\leq1$.

By solving the optimization problem (3.11), we can deduce the optimal harvest level $A_{c}^{*}$ and weight $\mu^{*}$, as well as the optimal predator release amount $\Lambda^{*}$, the optimal harvesting intensity $h^{*}$ and the optimal harvesting period $\bar{T}^{*}$. Note that these results depend on $n_{i}$, where $i = 1, 2, 3, 4$.

4.   Numerical simulation

This section presents some numerical simulations to validate the conclusions drawn in section 3. We adopt the following model parameters for numerical simulations:

4.1.   Simulation for system (2.1)

From (3.1) we can directly compute the positive equilibrium $P^{*}(A^{*}, B^{*})\approx(1.87, 1.61)$ with $\tau_{2} = 0.2, \omega = 0.1, \beta = 3.5$, $\frac{\tau_{2}K^{2}}{K^{2}+\beta}\approx 0.193 > \omega$. Under these conditions, the positive equilibrium point $P^{*}$ of system (2.1) is globally asymptotically stable, as illustrated in Figure 4(a)–(c). This observation aligns with the conclusion in Theorem 2. However, if $\omega = 0.3$, then $\frac{\tau_{2}K^{2}}{K^{2}+\beta}\approx 0.193 < \omega$, the positive equilibrium point disappears, and the boundary equilibrium point $P_{k} = (K, 0)$ becomes globally asymptotically stable, as shown in Figure 4(d)–(f). This finding is consistent with Theorem 1.

4.2.   Simulation for system (2.2)

For the complex pulse system (2.2) we selected the following harvesting parameters:

With $\Lambda = 0, I = 1.5, \mu = 0.5$ satisfying $I\leq\mu K = 5$ and $c_{2} = 0.3$ satisfying $c_{2} > \tilde{c}_{2}\approx 0.061$, a predator-extinction periodic solution exists for system (2.2). As $c_{2} > \tilde{c}_{2}$, this period solution is orbitally asymptotically stable, as depicted in Figure 5(a)–(c). This aligns with the conclusion of Theorem 3. If, however, $0 < c_{2} < \tilde{c}_{2}$ (for instance $c_{2} = 0.04$), the predator-extinction periodic solution becomes unstable, and a positive periodic solution emerges, as shown in Figure 5(d)–(f). This corresponds to the positive periodic solution for the case $(i)$ in Theorem 4.

From (3.9), we derive $\mu^{*}\approx 0.17, I_{m}\approx 1.95$. Choosing $\mu = 0.1, I = 1.5$ satisfies $\mu < \mu^{*}$ and $\mu K < I < I_{m}$. Under these conditions, the predator-extinction periodic solution does not exist. Figure 6(a) shows the existence of a positive periodic solution, which is orbitally asymptotically stable. When $I = 1.5, \mu = 0.4 > \mu^{*}$, the predator-extinction periodic solution emerges while the positive periodic solution disappears (see Figure 6(b)). It's worth mentioning that a positive period solution also exists when $I = 1.5, \mu = 0.2$ (see Figure 6(c)), suggesting that the conditions in Theorem 4$(ii)$ are necessary but not sufficient. If $I = 1.98 > I_{m}, \mu = 0.1$, the system's periodic behavior ceases, and the trajectory approaches the positive equilibrium point $P^{*}(1.87, 1.61)$ closely (see Figure 6(d)).

For $\Lambda\neq0$, with fixed parameters $I = 1.5 < I_{m}, \; \mu = 0.5$ and choosing $\Lambda = 0.8 < \Lambda_{m}$, the system (2.2) exhibits a positive order-1 periodic solution (see Figure 7(a), Theorem 4 $(iii)$). If $\Lambda$ = $1.32 > \Lambda_{m}$ while keeping other parameters unchanged, the predator population surpasses the threshold and the periodic solution becomes unstable, causing the prey population to become extinct after a finite number of pulses (Figure 7(c) and (e)). In order to maintain symbiosis between two populations, we can decrease prey harvesting and predator release while increasing the capture of adult predators. Setting $c_{1} = 0.3, \Lambda = 1, c_{2} = 0.6$ results in the disappearance of the positive periodic solution for system (2.2) and the trajectory converges to the positive equilibrium point $P^{*}(1.87, 1.61)$ (see Figure 7 (b)).

When the release amount of predator cubs is reduced to $\Lambda = 1.3$ (which is less than or equal to $\Lambda_{m}$) and the effective capture rate of the prey is set to $c_{1} = 0.4$, system $(2.2)$ exhibits a bistable phenomenon, as seen in Figure 8(a) and (b). This system features both a locally asymptotically stable equilibrium point and a positive periodic solution. This implies that in the system $(2.2)$, a pulse effect will not manifest if the weighted sum of the two populations is insufficient to activate it. However, a positive equilibrium does exist.

5.   Conclusions

This study presents and examines a Gilpin-Ayala predator-prey system with state feedback weighted harvest strategy and a Holling Ⅲ type functional response. Through the construction of two Lyapunov functions, we have primarily studied the global stability of both the coexistence equilibrium point $P^{*}$ and the boundary equilibrium point $P_{k}$. The position of the positive equilibrium point is influenced by the magnitude of the exponential factor $\alpha$, while its stability is predominantly affected by $\tau_{2}$, $\omega$ and $\beta$ (refer to Theorem 1, Figure 4). Aiming to optimize harvest economic value and support sustainable ecological development, we introduced a weighted harvest strategy and initiating harvesting once the weighted sum of the two populations reached a predetermined threshold, causing complex dynamic effects on the predator-prey model.

For the complex pulse harvesting system, our focus lay on analyzing the existence and stability of predator-extinction periodic solutions and positive order-1 periodic solutions. In the absence of predator cub release ($\Lambda = 0$), a predator-extinction periodic solution occurs when $I < \mu K$ and $B(0) = 0$. Further, this solution is asymptotically stable if $c_{2} > \tilde{c}_{2}$, as demonstrated in Theorem 3 (see Figure 5); that is, when no predator seedlings are placed and harvesting parameter $c_{2}$ meets the threshold condition, the predator population will become extinct. From the perspective of ecological balance and commercial value, we do not want any species to become extinct. To this end, options include introducing predator pups into the environment or reducing the predator capture rate $c_{2}$. Initially, without releasing predator cubs ($\Lambda = 0$) and keeping $c_{2} < \tilde{c}_{2}$, the harvesting system (2.2) yields an asymptotically stable positive order-1 periodic solution (refer to Theorem 4, Figure 5). Moreover, when predator cubs are released, a positive order-1 periodic solution exists provided $I\leq I_{m}$ and $\Lambda < \Lambda_{m}$ (see Theorem 4, Figures 6 and 7). It indicates that, due to human harvesting action, the number of individuals in the area will decrease, which is not conducive to population growth and harvesting. Proper supplementation of predator seedlings can help improve population development and increase economic benefits.

A clear demonstration of bistability in system (2.2) is evident from Figures 8. This crucial observation underlines the system's sensitivity to initial conditions. In the future, we hope to prove it theoretically through mathematical methods. Due to the complexity of the switching system, only the local stability of the positive order-1 periodic solution was proven in Theorem 5, but its global stability cannot be determined using our current methods. In the future, we will try appropriate methods to study it, such as constructing Lyapunov function methods. In addition, due to a more comprehensive exploration of the predator-prey relationship in nature, we will incorporate the Allee effect, time delay effect and optimal control into the model in future work.

Use of AI tools declaration

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

Acknowledgements

This work was supported by the Program for Science & Technology Innovation Talents in Universities of Henan Province (No. 21HASTIT026) and the Scientific Research Foundation of Graduate School of Xinyang Normal University (No. 2022KYJJ055).

Conflict of interest

The authors declare no conflict of interest regarding the publication of this paper.