Complex Dynamics of a Predator-Prey System With Gompertz Growth and Herd Behavior

. The complex dynamics of a predator-prey system in discrete time are studied. In this system, we consider the prey’s Gompertz growth and the square-root functional response. The existence of ﬁxed points and stability are examined. Using the center manifold and bifurcation theory, we found that the system undergoes transcritical bifurcation, period-doubling bifurcation, and Neimark-Sacker bifurcation. In addition, numerical examples are presented to illustrate the consistency of the analytical ﬁndings. The bifurcation diagrams show that the positive ﬁxed point is stable if the death rate of the predator is greater than a threshold value. Biologically, it means that to prevent the predator population from growing uncontrollably and stability of the positive ﬁxed point, the predator’s death rate should be greater than the threshold value.


Introduction
In mathematical ecology, a predator-prey interaction is essential due to its significance and universal existence.One of the foremost imperative subjects in mathematical ecology is the dynamic mutual action between predator and prey populations, which aids in conserving species in a habitat.The first predator-prey interaction was derived independently by Lotka [1], and Volterra [2], and it is known as the Lotka-Volterra predator-prey system.Since then, numerous researchers have made much progress considering different biological facts.
The functional response is an essential component of predator-prey interactions within population dynamics.It illustrates the link between the predator's consumption rate and the density of prey.
It denotes how much prey each predator consumes.Holling (1965) proposed three distinct types of functional responses [3].Subsequently, a number of functional responses were proposed by scholars such as Crowley-Martin [4] and Beddington-DeAngelis [5,6].Later, several research studies investigated models centered on predator-prey interactions, encompassing diverse categories of functional responses [7][8][9][10][11][12].Some prey populations exhibit herd behavior, where interactions between predators and prey occur primarily at the boundaries of the prey population.The nature of this interaction is not sufficiently explicable through Holling-type responses.It is noteworthy that a particular population of prey exhibits collective behavior known as herd behavior.Therefore, the rate at which a predator captures its prey differs from conventional models.As an example, the predation rate of zooplankton by fish in marine ecosystems exceeds the predation rate of phytoplankton.The observed phenomenon entails the manifestation of herding behavior among phytoplankton.The utilization of the square root of the prey population was employed by Ajraldi et al. [13] as a method to investigate the herd behavior of the prey population.This strategy facilitated the predator's ability to interact with the prey in the boundary region of the group.Numerous scholars have investigated the mechanisms of predator-prey relationships by utilizing square-root functional responses [14][15][16][17][18][19][20][21][22][23].
We examine a predator-prey interaction when the prey exhibits group defense, which is represented by the set of ordinary differential equations given as:    dx dt = xg(x, K) − y p(x), dy dt = −dy + q(x)y , (1.1) where x, y , K > 0, and d > 0 are the prey population density, predator population density, carrying capacity, and predator's death rate, respectively.In the absence of a predator, the function g(x, K) indicates the particular growth rate of the prey.We employ Gompertz growth [24] of prey g(x, K) = r ln( K x ), with a natural growth rate r > 0. The predator response function is represented by p(x), and we assume it is of square root type, i.e, p(x) = m √ x, where m > 0 represents the predator's search efficiency.Moreover, q(x) describes the conversion rate of prey.For some positive constant c, q(x) = cp(x) in Gauss' model.Considering all assumptions, system (1.1) takes the following form: Due to the square root term, the Jacobian of system (1.2) possesses a singularity.We utilize the transformations x(t) = u 2 (t) and y (t) = v (t) to understand it better.After applying this transformation to system (1.2), we obtain In cases where a population has no overlapping generations, discrete systems that are governed by difference equations are comparatively more appropriate than continuous systems.Several species, such as monocarpic plants and semelparous animals, exhibit distinct generations that remain different from one another, and reproductive events occur at anticipated intervals during specific mating seasons.Discrete-time mappings or difference equations are employed to depict their interrelationships.
Discrete-time dynamical systems have complicated and diversified dynamical properties [25][26][27][28][29][30][31].The system to be analyzed in this work is then obtained by using the forward Euler technique on the system (1.3) as follows: where h > 0 denotes the step size.The main contributions of this study are listed as follows: • A discrete-time predator-prey system with Gompertz growth and herd behavior is proposed.
• The study investigates the existence and topological categorization of fixed points.
The following outlines the format of the paper: Section 2 investigates the existence and stability of fixed points in the system (1.4).In Section 3, we employ the center manifold theorem and bifurcation theory to analyze local bifurcation analysis at fixed points of the system (1.4).Section 4 presents some examples to verify our theoretical results.Some closing observations are added in Section 5.

Existence and stability of fixed points
The long-term behavior of dynamical system depends on fixed points.Stable fixed points cause the system to converge, while unstable ones cause oscillation or chaos.The existence and stability conditions for the fixed points of system (1.4) are investigated in this section.By simple algebraic computations, the system (1.4) is found to have two fixed points: ).The first fixed point P 1 is a boundary point.Biologically, it indicates that when there are no predators around, the number of prey approaches the square root of the carrying capacity.The second fixed point, P 2 is the unique positive fixed point in the system (1.4) if K > d 2 c 2 m 2 .Now, we will look at the local stability of the fixed points in the system (1.4).The eigenvalues of the Jacobian matrix computed at the fixed points define the local stability of the fixed points.At the point ( ū, v ), the Jacobian matrix of system (1.4) is We apply the following results to examine the stability of the fixed points: Lemma 2.1.[32] Let P (q) = q2 + Sq + T be the characteristic polynomial of J( ū, v ), and q 1 , q 2 be the two roots of P (q) = 0.
Lemma 2.2.[32] Let P (q) = q 2 + Sq + T and P (1) > 0. If q 1 , q 2 are roots of P (q) = 0, then The Jacobian matrix computed at the fixed point P 1 ( √ K, 0) is The following result describes the topological classification of P 1 ( √ K, 0): (ii) a source if any of the following is true: (iii) a SP if any of the following is true: (iv) NHP if any of the following is true: It is evident that if d = cm √ K, then one of the eigenvalues of J(P 1 ) is 1.Consequently, a transcritical bifurcation may occur if the parameters change in a close neighborhood of Γ 1 .
a PD bifurcation may occur if the parameters change in a close neighborhood of Γ 2 or Γ 3 , where and The Jacobian matrix at .
The characteristic polynomial of J(P 2 ) is where and Through simple calculations, we obtain Using the lemma (2.1), we acquire the local dynamics of the fixed point (2) source if any of the following is true: (3) SP if any of the following is true: and ) NHP if any of the following is true: The condition for J(P 2 ) eigenvalues to be unit-modulus complex is established when Thus, the NS bifurcation occurs at point P 2 in the system (1.4) when the parameters are altered in the vicinity of Ω 1 .
, one of the eigenvalues of the matrix J(P 2 ) is −1, while the other eigenvalue λ meets the condition that |λ| = 1.If the parameters fluctuate slightly around Ω 2 or Ω 3 , a PD bifurcation may occur, where and

Local Bifurcation Analysis
The present section delves into various types of fixed-point bifurcations that could potentially manifest in the system (1.4).Bifurcations in predator-prey systems manifest as a result of alterations in the parameters.A slight modification in the parameter leads to a bifurcation.Bifurcations in predator-prey systems play a crucial role in forecasting the dynamics of wild populations and formulating viable approaches for sustainable management.Improper management of bifurcations can lead to disturbances in population dynamics and the consequent destruction of ecosystems.For a detailed bifurcation analysis, we recommend that readers refer to [31,[33][34][35][36][37][38][39].
This section discusses transcritical bifurcation at the boundary fixed point P 1 ( √ K, 0) using center manifold theory.Assuming that (h, r, K, m, c, d) ∈ Γ 1 , and γ be small perturbation in d, the subsequent perturbation of the system (1.4) is taken into consideration: In order to translate fixed point P 1 ( √ K, 0) to (0, 0), we define the translation map as follows: As a result of this translation map, the system (3.1)transforms to .
The following theorem provides the parametric conditions for the presence and direction of transcritical bifurcation for system (1.4) at its boundary fixed point P 1 ( √ K, 0).

PD Bifurcation at
The PD bifurcation at the fixed point P 1 ( √ K, 0) in Γ 2 is discussed in this section.For the domain Γ 3 , a similar investigation can be done.Suppose that (h, r, K, m, c, d) ∈ Γ 2 , and γ be minimal change in h, the subsequent perturbation of the system (1.4) is taken into consideration: In order to translate the fixed point P 1 ( √ K, 0) to (0, 0), we define the translation map as follows: As a result of this translation map, the system (3.5)transforms to Under the following transformation where Subsequently, the center manifold W C for the system (3.8) is computed and characterized as follows: where System (3.8) when confined to the center manifold is expressed as follows: F : For map (3.9) to experience PD bifurcation, we need the value of the following two expressions to be non-zero: Fenenen .
From simple computations, we obtain The following theorem provides the parametric conditions for the presence and direction of PD bifurcation for system (1.4) at its boundary fixed point P 1 ( √ K, 0).

PD Bifurcation at
).The PD bifurcation at the fixed point P 2 ( d cm , r d cm 2 ln( kc 2 m 2 d 2 )) for the domain Ω 3 is discussed in this section.For the domain Ω 2 , a similar investigation can be done.Suppose that (h, r, K, m, c, d) ∈ Ω 3 , and γ be minimal change in h, the subsequent perturbation of the system (1.4) is taken into consideration: (3.10) In order to translate fixed point P 2 ( d cm , r d cm 2 ln( kc 2 m 2 d 2 )) to (0, 0), we define the translation map as follows: ).As a result of this translation map, the system (3.10)transforms to where , the eigenvalues of J(P 2 ) are λ 1 = −1 and Under the following transformation where .
System (3.13) restricted to the center manifold is For map (3.14) to experience PD bifurcation, we need the value of the following two expressions to be non-zero: Fenenen .
From simple computations, we obtain The following result is drawn from the calculations mentioned above: Theorem 3.3.Suppose that (h, r, K, m, c, d) ∈ Ω 3 .The system (1.4) undergoes PD bifurcation at the fixed point . Moreover, if l 2 > 0 (respectively l 2 < 0), then the period-2 orbits that bifurcate from ) are stable (respectively, unstable).

NS Bifurcation at
).This section discusses the NS bifurcation at the fixed point P 2 ( d cm , r d cm 2 ln( kc 2 m 2 d 2 )) for the domain Ω 1 .Suppose that (h, r, K, m, c, d) ∈ Ω 1 , and γ be minimal change in h 3 , the subsequent perturbation of the system (1.4) is taken into consideration: We define ) to (0, 0).As a result of this translation map, the system (3.15)transforms to where The characteristic equation corresponding to the linearized part of the system (3.16) at the origin is where The equation (3.17) roots are complex that have the property |λ 1,2 | = 1, which are given by By computations, we obtain and Moreover, it is required that A 2 , therefore α(0) = ±2.We only need to require that α(0) = 0, 1, which leads to A 2 1 = 2A 2 , A 2 .The canonical form at γ = 0 of the linear part in (3.16) is achieved through the utilization of the subsequent transformation: Upon application of the mapping (3.18), the system (3.16)transforms as follows: where The below-mentioned number L explains how the invariant curve appears in a system going through the NS bifurcation.
, where From the calculations shown above, we can obtain the following result: Theorem 3.4.Consider (h, r, K, m, c, d) ∈ Ω 1 and A 2 1 = 2A 2 , A 2 .When the parameter h differs in a small neighborhood of h 3 = A 1 A 2 , the system (1.4) undergoes NS bifurcation at the fixed point P 2 if L = 0. Additionally, if L > 0, a repelling invariant closed curve bifurcates from P 2 , whereas if L < 0, an attracting invariant closed curve bifurcates from P 2 .

Numerical examples
This section will strengthen our theoretical investigations of the system by presenting some numerical examples that illustrate its numerous qualitative properties.
4.1.PD bifurcation of the system (1.4) at P 2 by using h as bifurcation parameter.
Consider the following parametric values: c = 0.9, m = 0.5, d = 0.15, r = 1.5, k = 0.12, and initial conditions u 0 = 0.5, v 0 = 0.05.For these values, the bifurcation value is h = 1.39251, and the positive fixed point is obtained as P 2 (0.333333, 0.076961).The eigenvalues of J(P 2 ) are obtained as λ 1 = −1, λ 2 = 0.991606, confirming that the system (1.4) undergoes period doubling bifurcation at P 2 (0.333333, 0.076961) as h passes through h 0 = 1.39251.Figures (1a,  2d-2k) depict phase portraits of the system (1.4) for some values of h.From the figures, P 2 is a LAS for h < 2.19813 but loses stability at h = 2.19813, where the system (1.4) goes through NS bifurcation.A smooth curve that is invariant emerges for values of h ≥ 2.19813, with its radius increasing proportionally to the increase in h.The sudden disappearance of the invariant curve and subsequent emergence of a  Sensitive dependency on initial conditions implies that small changes in the initial populations of predator and prey may result in drastically different long-term population dynamics.This phenomenon highlights the importance of accurate and precise measurements of the initial conditions in ecological investigations and the requirement for strong modeling tools that can account for data uncertainty and variability.It also highlights the difficulties in forecasting the long-term dynamics of ecological systems since even minor errors in the initial conditions may lead to considerable errors in predictions.

Conclusion
The nonlinear dynamics of a novel discrete-time predator-prey system with square root functional response and Gompertz growth of prey, generated using the forward Euler discretization method, were examined in this study.The system has a boundary fixed point ( √ K, 0) which always exists, and an interior fixed point ( d cm , r d cm 2 ln( kc 2 m 2 d 2 )) which exists only if the carrying capacity of prey is greater than a threshold value.Using bifurcation theory and the center manifold theorem, it is shown that the system's fixed points have transcritical bifurcation, PD bifurcation, and NS bifurcation.The less integral step size h can stabilize the positive fixed point.However, the large integral step size can destabilize the positive fixed point, resulting in more complex dynamical behaviors, as depicted in the figures.Moreover, we see that if the death rate of the predator is greater than a threshold value, then the positive fixed point is stable; otherwise, it is unstable.This threshold value shows the lowest death rate necessary to prevent the predator population from exploding and to maintain a stable equilibrium between predator and prey populations.Moreover, it is demonstrated that the trajectories of the system (1.4) are sensitive to the initial values.Slight variations in initial conditions can significantly affect the long-term dynamics of ecological systems.The sensitive dependence on initial conditions emphasizes the importance of accurate and precise measurements and the need for robust modeling techniques in ecological studies.

Conflicts of Interest:
The authors declare that there are no conflicts of interest regarding the publication of this paper.
figures for h ∈ [1.1, 1.7].The MLE is plotted in figure (1c).The fixed point is a LAS for these parametric values if 0 < h < 1.39251.Figures(1d,1e,1f) depict phase portraits of the system (1.4) for distinct values of h.From the figures, the P 2 (0.333333, 0.076961) is a LAS for 0 < h < 1.39251, but fails to retain stability at h = 1.39251,where the system (1.4) experiences PD bifurcation.

Figure ( 4 )
Figure(4) depicts 2 perturbed trajectories in blue and red colors to highlight the sensitivity of the system (1.4) to initial conditions.The two trajectories are initially overlapping and identical, but the divergence between them grows fast after a few repetitions.Figure (4) plots the u− and v −coordinates of the two trajectories for the system (1.4) against the number of iterations, revealing a reactive reliance on the initial conditions.The initial perturbation of 2 trajectories is 0.0001.The 2 trajectories with initial points (u 0 , v 0 ) = (0.4,1.4) (trajectory in blue color) and (u 0 , v 0 ) = (0.4001, 1.4001) (trajectory in red color) are calculated and plotted in figure(4), respectively.The trajectories of the system (1.4) are depicted in figure (4) to be sensitively dependent on the initial conditions, i.e., complex dynamic behavior happens with an initial perturbation.

Figure 4 .
Figure 4. Sensitive dependence on initial conditions of the system (1.4) 15, and the positive fixed point is obtained as P 2 (0.333333, 1.50408).The eigenvalues of J(P 2 ) are λ 1 = 0.591212 − 0.806518i , λ 2 = 0.591212 + 0.806518i with |λ 1,2 | = 1, indicating that the system (1.4) is experiencing NS bifurcation at P 2 (0.333333, 1.50408) as the bifurcation parameter d passes through d 0 = 0.15.Figures (3a,3b) depict bifurcation diagrams for d ∈ [0.14, 0.17].The MLE is plotted in figure3c.The bifurcation diagrams show that the positive fixed point is stable if the death rate of the predator is greater than a threshold value.Biologically, it means that to prevent the predator population from growing uncontrollably and stability of the positive fixed point, the predator's death rate should be greater than the threshold value.