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

: The current research presents a predator-prey model that incorporates both a Gilpin-Ayala growth function and a Holling type III functional response. Two Lyapunov functions are established to conﬁrm the global asymptotic stability of the positive equilibrium P ∗ and the predator extinction equilibrium P k . Considering ecological protection and commercial incentives, we also incorporated a weighted harvesting strategy and pulse control into the model. We investigated intricate dynamical problems instigated by the weighting harvesting and pulse e ﬀ ects, and a ﬃ rmed the existence and local asymptotic stability of both predator-extinction periodic solution and positive order-1 periodic solution. In the end, a suite of numerical simulations were carried out using MATLAB, aiming to corroborate the theoretical ﬁndings and deliver conclusions rooted in a biological context.


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-II 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 predatorprey 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 α 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 α ∈ (0, 1), they transformed the model into a nonlinear system, thereby increasing its analytical complexity [22].Nevertheless, this exponential parameter α 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 šlez-Olivares and Rojas-Palma [9] furthered this area of study by proposing a model incorporating the Holling III 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 III 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 III 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.

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 III 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 α is a positive parameter to modify the classical logistic model and represents a nonlinear measure of interspecific interference.The term τ 1 A 2 A 2 +β denotes the functional response of the predator, commonly referred to as the Holling type III response function.The variables τ 1 , τ 2 , β and ω 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 µ and 1 − µ, respectively.
The model driven by this weighted harvesting strategy can be expressed as follows: where c 1 and c 2 represent capture rates, while Λ 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 M and N, respectively (see Figure 1).Let O 0 be the intersection point of N and the A-axis.For any point V 1 ∈ N, we denote the trajectory that originates from V 1 by Γ 1 (V 1 ).If the trajectory Γ 1 (V 1 ) intersects with the phase set N and pulse set M at points V 2 and V 3 , and point V 3 is mapped to a point V 4 ∈ N through pulse effect, then The schematic diagram of the successor function.
Lemma 1. (Lasalle invariance principle) Consider the following system of differential equations where 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 Ḡ and possesses a continuous first-order partial derivative, V(A, B) is referred to as a Lyapunov function on G if Assuming M ⊂ S is the largest invariant subset within S , and given a solution x(t, t 0 , (A 0 , B 0 )) ⊂ G of the system, it follows that x(t, , t 0 , (A 0 , B 0 )) → M as t → ∞.

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 dA dt | A=0 = 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, where f (A) ≥ 0 for A ∈ (0, K], and f (A) < 0 for A > K. Thus, we obtain A(t) ≤ max{A(0), K}.Subsequently, we denote Through direct calculation, we can deduce that the maximum point of f (A) is A m = K α ω+r r(α+1) .Denote the maximum value by M m = f (A m ), we have Define and we can easily know that Ω 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 τ 2 K 2 K 2 +β > ω, 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 > ω, 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.
The preceding discussion establishes that a positive equilibrium 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 τ 2 K 2 K 2 +β > ω, if one of the following conditions holds: (I) ξ 1 (A) = 0 has no or only one positive root.(II) ξ 1 (A) = 0 has two positive roots A 1 and A 2 (where A 1 < A 2 ).If A * ∈ (0, A 1 ), then F(A) > B * for all A ∈ (0, A * ) and F(A) < B * for all A ∈ (A * , K) holds; if A * ∈ (A 2 , K), then F(A) > B * for all A ∈ (0, A * ) and F(A) < B * for all A ∈ (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 then we have Clearly, 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 then it follows, by differentiate with respect to A, that In the special case of α = 1, we obtain the maximum point For the case α ∈ (0, 1), when 1+α and ξ 1 (A) = 0, we have When 0 < α < 1, we observe that ξ 2 (0) 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) = ξ 2 (A) and y 1 (A) = A K α either have no intersection points, or they intersect at a single point A c such that √ β < A c < K, or they intersect at two points, A 1 and A 2 , where Using an analogous analysis for α > 1, the distribution of the roots for ξ 1 (A) = 0 exhibits a similar behavior.Therefore, for any α > 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 dV dt ≤ 0 for all (A, B) ∈ Ω 0 , and dV dt = 0 if, and only if, (A, B) = (A * , B * ) (see Figure 2 (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 ∈ (0, A * ) and F(A) < B * for A ∈ (A * , K) when A * ∈ (0, A 1 ); (b) F(A) > B * for A ∈ (0, A * ) and F(A) < B * for A ∈ (A * , K) when A * ∈ (A 2 , K).

. The geometric profile of function F(A).
Thus, we can also claim that dV dt ≤ 0 for all (A, B) ∈ Ω 0 , and 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.

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) ∈ 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).
Assume H is a point on the phase set with dB dA . 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: (3.9) 2) has a positive order-1 periodic solution, if one of the conditions (i)-(iii) holds: Proof.(i) Referring to Theorem 1, the boundary periodic solution of (2.2) is unstable under this condition.For any P 1 ∈ N ∩ U(L, δ), where δ is sufficiently small, the trajectory originating from P 1 intersects the impulse set M at point P − 1 and then is instantly pulsed to the point P + 1 in the phase set 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 ) ∈ N satisfies function (3.6), the trajectory of system (2.2) starting from P 2 (A 2 , B 2 ) first intersects the impulse set M at P − 2 (A − 2 , B − 2 ) and then pulses to the phase set N at P + 2 .There are three possible situations: (1)P + 2 = P 2 , (2)P + 2 = P + 2 and (3)P + 2 = 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 = P + 2 , then we obtain that S 1 (P 2 ) = d(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 ∈ 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 Since the trajectory of system (2.2) starting from P 2 (A 2 , B 2 ) is tangent to the phase set N, the trajectory of system (2.2) starting from P + 2 crosses N once more at P + 2 before reaching 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 S 2 (P In addition, due to the continuity of solution in the system (2.1) and for ε = d(P 2 , P + 2 )/2, we know that there must exist some δ < ε such that D ∈ U(P 2 , δ), d(D − , P − 2 ) < ε, and d(D Thus, there must exist a point E ∈ P 2 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 M. Let the intersections of 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 µK < I < I m , the system doesn't yield an extinction periodic solution, and all orbits starting from the phase set N reach the pulse set M. Similar to the proof in (i), choosing P 1 ∈ N with B P 1 = (1 − τ 2 h)B Z , we can confirm the existence of an order-1 periodic solution of system (2.2).
(iii) Choose P 1 ∈ N with d(P 1 , L) = Λ, 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.

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 T .Let's consider 2 ) and B c := A c to represent the target densities of prey and predator, respectively.As h relies on A c , and both h and Λ depend on µ and A c , we assume that . We can also infer that T depends on µ 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., max where A 2 < A c < A − 2 and 0 ≤ µ ≤ 1.By solving the optimization problem (3.11), we can deduce the optimal harvest level A * c and weight µ * , as well as the optimal predator release amount Λ * , the optimal harvesting intensity h * and the optimal harvesting period T * .Note that these results depend on n i , where i = 1, 2, 3, 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: K 2 +β ≈ 0.193 > ω.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 ω = 0.3, then τ 2 K 2 K 2 +β ≈ 0.193 < ω, 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.

Simulation for system (2.2)
For the complex pulse system (2.2) we selected the following harvesting parameters: With Λ = 0, I = 1.5, µ = 0.5 satisfying I ≤ µK = 5 and c 2 = 0.3 satisfying c 2 > c2 ≈ 0.061, a predator-extinction periodic solution exists for system (2.2).As c 2 > c2 , this period solution is orbitally I = 1.5, µ = 0.4 > µ * , 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, µ = 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 , µ = 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 Λ 0, with fixed parameters I = 1.5 < I m , µ = 0.5 and choosing Λ = 0.8 < Λ m , the system (2.2) exhibits a positive order-1 periodic solution (see Figure 7(a), Theorem 4 (iii)).If Λ=1.32 > Λ 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, Λ = 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   When the release amount of predator cubs is reduced to Λ = 1.3 (which is less than or equal to Λ 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.

Conclusions
This study presents and examines a Gilpin-Ayala predator-prey system with state feedback weighted harvest strategy and a Holling III 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 α, while its stability is predominantly affected by τ 2 , ω and β (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 (Λ = 0), a predator-extinction periodic solution occurs when I < µK and B(0) = 0. Further, this solution is asymptotically stable if c 2 > c2 , 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 (Λ = 0) and keeping c 2 < c2 , 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 ≤ I m and Λ < Λ 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.