Qualitative analysis of a prey–predator model with prey refuge and intraspecific competition among predators

In this study, we consider a prey–predator model with prey refuge and intraspecific competition between predators using the Crowley–Martin functional response and investigate the dynamic characteristics of spatial and nonspatial prey–predator systems via both analytical and numerical methods. The local stability of nontrivial interior equilibrium, the existence of a Hopf bifurcation, and the stability of bifurcating periodic solutions are obtained in the absence of diffusion. For the spatial system, the Turing and non-Turing patterns are evaluated for some set of parametric belief functions, and we obtain some interesting results in terms of prey and predator inhabitants. We present the results of numerical simulations that demonstrate that both prey and predator populations do not converge to a stationary equilibrium state at any foreseeable future time when the parametric values are processed in the Turing domain.


Introduction
Understanding the dynamics of interacting prey-predator models is of paramount importance for gaining insights into the long-term behavior of species in ecological systems.Over the years, researchers have extensively investigated various nonlinear coupled ordinary differential equations to model the complex interactions among prey and their potential predators.Among the classical models, the Lotka-Volterra model, introduced independently by Lotka and Volterra in the early twentieth century, has served as a fundamental basis for studying predator-prey interactions.Subsequent advancements by Holling [18] in the form of the Rosenzweig-MacArthur model [33] incorporated densitydependent prey growth and functional response.
While these classical models have provided valuable insights into the general characteristics of prey-predator populations, recent research has shown that incorporating additional factors can better reflect the intricacies of natural ecosystems.Consequently, researchers have introduced various functional responses to represent the feeding rates of predators [5,49,50].Notably, the Lotka-Volterra and Holling types have been widely used in this context.
Additionally, the most crucial concept in prey-predator models is the 'functional response' , which describes the rate at which a predator attacks a given number of prey.The Lotka-Volterra and Holling types are the most commonly used functional responses to represent the average feeding rate of a predator.
Lotka-Volterra type : f 1 (u) = δu, Holling Type-II : f 2 (u) = δu/(1 + αu), Holling Type-III : f 3 (u) = δu 2 / 1 + αu 2 , Holling Type-IV : f 4 (u) = δu/ 1 + αu 2 , where u denotes the density of the prey population, and δ and α are positive constants that respectively represent the effects of capture rate and handling time.Notably, these four functional response features are prey-dependent and generally not influenced by predator inhabitants.However, a crucial aspect often overlooked in these models is the consideration of mutual interference among potential predators.To address this limitation, Beddington [2] and DeAngelis [11] proposed a functional response f (u, v) = δu/(1 + αu + βv) that accounts for predator (v) interference β, while Crowley and Martin [10] introduced a similar response, taking into account both searching and handling interference: which is referred to as the Crowley-Martin functional response.These modifications have proven essential in understanding the dynamics of real-world predator-prey systems, and they offer a more accurate representation of their behavior [38].
In light of these considerations, the present study focuses on formulating a preypredator model that accounts for prey refuge and intraspecific competition among predators using the Crowley-Martin functional response.The inclusion of prey refuge and competition among predators in the model is crucial as it addresses important ecological aspects that have significant implications for species coexistence and ecosystem stability.
The rest of this paper is organized as follows: Sect. 2 presents the mathematical formulation of the prey-predator model with prey refuge and intraspecific competition among predators using the Crowley-Martin functional response.In Sect.3, we analyze the existence and local and global dynamics of equilibria.Subsequently, in Sect.4, we discuss the presence of Hopf bifurcation and its periodic solution behaviors.Section 5 focuses on the occurrence of Turing instability in the interior equilibrium point under the influence of diffusion.Theoretical results are verified through numerical simulations in Sect.6, where we also uncover additional meaningful phenomena.Finally, Sect.7 offers concluding remarks and highlights the significance of our findings within the broader context of ecological research.
By comprehensively investigating the dynamic behavior of this extended prey-predator model, we aim to contribute to a deeper understanding of ecological interactions and enhance the applicability of predator-prey models in real-world scenarios.

Mathematical model and analysis
Bazykin [1] introduced a prey-predator model based on the Crowley-Martin (C-M) functional response, which considers the simultaneous involvement of potential predators in various activities, including searching and handling prey as well as interacting with other predators.The C-M functional response is described by the total instantaneous per capita feeding rate equation Here, X denotes the prey population, Y denotes the predator population, and c (units: 1/predator) represents the magnitude of interference among predators.The term Y -1 is used to account for the mechanistic nature of predator dependence, indicating that a predator does not interfere with itself in these models.When Y = 1, the equation reduces to the traditional Holling Type II response.
Prey refuge, a significant concept in prey-predator models, plays a crucial role in supporting a constant proportion of prey through predation [16,17].It profoundly influences the coexistence of predators and prey.In theoretical ecology, understanding the impact of prey refuge on the dynamics of prey-predator interactions has been a topic of interest.Many authors [9,12,21,34] have found that prey refuge stabilizes prey-predator dynamics, preserving prey biomass and preventing extinction due to predation.Therefore, the strength of prey refuges affects the foraging efficiency of predators.
To incorporate the effect of prey refuge, we integrate it into the C-M functional response term.As a result, we arrive at the nonautonomous C-M type prey-predator model, accounting for prey refuge and intraspecific competition among predators for prey: (2.1) In this model, X(t) represents the biomass of the prey population, and Y (t) denotes the biomass of the predator population at time t.The refuge protects a population of prey denoted by mX, where m ∈ [0, 1) is a constant, making (1m)X the only prey available to predators.The parameters r 1 , r 2 , ρ, δ 1 , δ 2 , δ 3 , χ 1 , χ 2 , χ 3 , and χ 4 are all positive constants with ecological meaning as follows in Table 1.
For simplicity, we proceed to nondimensionalize (2.1) using the following scaling: This leads us to the following form of the model: In this formulation, we define dimensionless constants as follows: , and σ = δ 3 ρ δ 2 .Our physical world exists within a spatial universe, where environmental interactions play a vital role in shaping ecological communities.Various aspects, such as the atmosphere, geological materials, and biological frameworks, differ significantly across the planet.Consequently, population dynamics depend on both space and time, incorporating spatial movement.To analyze the spatial dynamics of prey-predator models, we consider the following system of partial differential equations with homogeneous Neumann boundary conditions: In this spatial extension, u(x, t) and v(x, t) represent the densities of prey and predators, respectively, in the fixed open bounded domain ⊂ R N at time t.The diffusion coefficients D 1 and D 2 account for the spatial movement of prey and predators, respectively.The Laplacian operator ∇ 2 = ∂ 2 ∂u 2 + ∂ 2 ∂v 2 describes spatial interactions in a two-dimensional space.
The boundary ∂ of the domain is smooth, and ν represents the outward unit normal vector of the boundary.The homogeneous Neumann boundary conditions indicate that the prey-predator system is self-contained with no population flux across the boundary.

Equilibrium points and stability analysis 3.1 Equilibria
System (2.2) has trivial and axial equilibrium points where E 0 , E 1 , E 2 always exist, where the following occur: (i) Trivial equilibrium point: E 0 = (0, 0), both prey and predator are extinct.
(ii) Axial equilibrium points: E 1 = (r, 0), only prey survives; + can be obtained by solving the following system of equations: Solving these nullclines, we obtain a cubic equation in u: where (H): provided the cubic equation (3.1) has at most one positive solution u * (say).Using the value of u * , we obtain the value of v * as follows: Below we assume that all E * (u * , v * ) ensure the above conditions (H) and ζ > σ .

Local stability analysis
The dynamic behavior of the equilibrium points can be studied by computing the eigenvalues of the Jacobian matrix J of system (2.2), namely, The existence and local stability of the equilibrium solutions can be stated as follows.
(ii) The Jacobian matrix of system (2.2) evaluated at the axial equilibrium point E 1 = (r, 0) is given by tr Therefore, E 1 is locally asymptotically stable when ζ < σ .
(iii) The Jacobian matrix of system (2.2) evaluated at the axial equilibrium point Then the eigenvalues of the above Jacobi matrix , both eigenvalues are negative.Hence, the axial equilibrium point E 2 is locally asymptotically stable when and because of our assumption that r > 0. This shows that the locally asymptotically stable situation never occurs.Therefore, E 2 is always unstable.Also note that ζ < σ ⇒ E 2 < 0 (i.e., when E 2 / ∈ {{0} ∪ R} 2 + , the predators also die out).

Interior equilibrium qualitative behaviors
The Jacobian matrix evaluated at the coexistence equilibrium point E * (u * , v * ) is where Then the trace and determinant of the Jacobian matrix (3.3) is Therefore, the characteristic equation of the linearized system of (2.2) at The qualitative behaviors of the interior equilibrium point E * (u * , v * ) may be stated as follows.

Case (i):
If T < 0 and D > 0, then the characteristic roots of (3.4) are either both negative reals or complex conjugates with a negative real part.Therefore, E * is either a stable node (T 2 > 4D) or a stable spiral (T 2 < 4D).

Case (ii):
If T > 0 and D > 0, then the characteristic roots of (3.4) are either both positive reals or complex conjugates with a positive real part.Hence, E * is either an unstable node (T 2 > 4D) or an unstable spiral (T 2 < 4D).Case (iii): If D < 0 and T = 0, then the characteristic roots of (3.4) are both are real with opposite signs.Therefore E * is a saddle point.Case (iv): If D < 0 and T = 0, then the characteristic roots of (3.4) both are real numbers with same magnitude and opposite signs.Hence, E * is a saddle node.Case (v): If D > 0 and T = 0, then the characteristic roots of (3.4) are purely complex conjugate.Therefore, E * is a center.

Global stability analysis
Here, we provide the result to attain global stability in the nontrivial interior equilibrium 2) and spatial (2.3) systems.
Proof The proof is usually developed by applying the Lyapunov function.We consider the subsequent positive definite Lyapunov function in R about the equilibrium E * : This kind of Lyapunov function was first utilized in [13], and later it was broadly exploited by several researchers.We can simply verify that ∂V 1 ∂u > 0 for u > u * and ∂V 1 ∂u < 0 for 0 < u < u * , and 2), we obtain where Similarly, where We define the Lyapunov function Generating dV dt via (3.5) and (3.6) gives The coefficient of (uu * ) 2 is and (1m)u * (βα) < 1 hold consequently from (3.7), then we obtain dV dt < 0. Hence, by Lyapunov's asymptotic stability theorem, the interior equilibrium E * of system (2.2) is globally asymptotically stable.Now, we select the Lyapunov function for the diffusion system (2.3) Hence, differentiating E(t) with respect to t along the solutions of system (2.3), we obtain Considering the zero-flux boundary conditions ∂ ν u = ∂ ν v = 0, x ∈ ∂ and applying Green's first identity in the plane, we obtain Therefore, the equilibrium E * of the spatial system (2.3) is globally asymptotically stable.

Stability behavior of Hopf bifurcation
Here, we discuss the behavior of Hopf bifurcation.For this specific purpose, we propose the perturbation u = u 1 + u 1 in our local system (2.2).Then, broadening in Taylor series, we obtain where a 10 = r -2u - η+(1-m)uσ ).Therefore, a 10 + b 01 = 0 and a 10 b 01a 01 b 10 > 0. Various other coefficients are to be determined as given below. 3, 3  , 4 , where , By [32], the first Lyapunov number to conclude the dynamics (stable or unstable) of the limit cycle arising through Hopf bifurcation has Theorem 4.4 When < 0, the direction of Hopf bifurcation is supercritical and the bifurcated periodic solutions are stable; when > 0, the direction of Hopf bifurcation is subcritical and the bifurcated periodic solutions are unstable.

Diffusion-driven instability
Through this section, we concentrate on the prey-predator system with self-diffusion and examine the occurrence of Turing instability in the equilibrium point under diffusion effects (diffusion-driven instability).Theorem (4.3) states that whenever ξ > ξ 0 , the nontrivial interior equilibrium E * is locally asymptotically stable for the nondiffusion system (2.2).
We consider the influences of diffusion on the stable nontrivial interior equilibrium E * of (2.3) under the supposition ξ > ξ 0 .Subsequently, for the diffusion system (2.3), we should consider the one-dimensional space = (0, π) with a smooth boundary ∂ : (5.1) This is actually the notable operator u → -u xx with Neumann boundary conditions.The analogous eigenvalues and normalized eigenfunctions are where k = 1, 2, 3, . . . .
Linearizing the above diffusion system (5.1) at E * , we obtain where and J is the Jacobian matrix pointed out in Sect.3.3.L indicates a linear operator whose domain is (5.2) Thus, all eigenvalues of L are obtained by the eigenvalues of where and Simply by examining the distribution of characteristic roots of J k , we obtain the imminent conclusion.
Theorem 5.1 Assume that ξ > ξ 0 and D(ξ ) > 0. Then (i) the equilibrium E * = (u * , v * ) of the nondiffused system (2.2) is locally asymptotically stable; (ii) the equilibrium E * = (u * , v * ) of the diffused system (5.1) is locally asymptotically stable if and only if the following conditions hold: Furthermore, for the diffused system (5.1), the solution E * is unstable (that is, Turing instability occurs) if Proof Because ξ > ξ 0 (T < 0) signifies that T k < 0 for all k ≥ 0.Moreover, by the meaning of T k , we have the relationship that for every k ≥ 0, Consequently, the real part of the characteristic values signs of (5.3) are discovered by the sign of D(k 2 ), separately.The symmetric axis of the graph We know that when D k < 0 (D(k 2 ) < 0), the characteristic roots of J k (5.3) are two real roots with opposite signs.Notice that in D(k 2 ), d 1 d 2 > 0 and k 2 > 0. Therefore, whenever ξ d 1ξ 0 d 2 < 0, D(k 2 ) attains the bare minimum at k 2 = k 2 min .Therefore, when (C3) holds, D(k 2 min ) is negative.However, this means that any one of the characteristic roots of J k have a positive real part, that is, that E * is the unstable solution of (5.1).As a result, we consider that whenever (C3) holds, Turing instability occurs.(C1) implies that D k > 0 for all k ≥ 0 (because D 0 = D > 0), and (C2) implies that D(k 2 min ) is positive, and therefore that all the characteristic roots of J k have negative real parts.In this manner, any of conditions (C1) and (C2) guarantee that the characteristic roots of J k have negative real parts.Consequently, if any one of the conditions of (C1) and (C2) holds, then E * is the stable equilibrium solution of (5.1).
By fixing the same set of above parameters other than ζ and σ with condition ζ = 1.4 < σ = 1.7, our resultant phase plane is shown in Fig.
Moreover, when ξ goes through ξ 0 from the right-hand side of ξ 0 , the equilibrium point E * loses its steadiness and a Hopf bifurcation occurs, as shown in Fig. 5 and Fig. 6.For the above set of parameters, we obtain = -9.01885< 0. Therefore, from Theorem 4.4, Figure 4 Stability behavior with respect to time t and phase portraits of system (2.2) with ξ > ξ 0 and initial data (u 0 , v 0 ) = (0.5, 0.3) Figure 5 Periodic behavior with respect to time t and phase portraits of system (2.2) with ξ = ξ 0 and initial data (u 0 , v 0 ) = (0.32, 0.36) Figure 6 Unstable behavior with respect to time t and phase portraits of system (2.2) with ξ < ξ 0 and initial data (u 0 , v 0 ) = (0.32, 0.36) the direction of the Hopf bifurcation at ξ = ξ 0 is supercritical, and the bifurcating periodic solutions are asymptotically stable, as shown in phase diagrams of Fig. 5 and Fig. 6.
The formation of a limit cycle around the interior equilibrium point with initial data (0.4, 0.3) and (0.5, 0.5) inside and outside of the limit cycle, respectively, is shown in Fig. 7.The diagram shows the limit cycle is stable.
To verify the occurrence of Turing of the diffusive prey-predator system (5.rium solution E * of system (5.1) is unstable.That is, Turing instability occurs.The Turing instability and the corresponding contour diagram are shown in Fig. 9.

Conclusion
Our study focused on investigating the dynamic behavior of a prey-predator model that incorporates prey refuge and interference among predators.The interaction between prey and predators was governed by the Crowley-Martin response function.Our key contribution lies in the thorough analysis of the refuge function during intraspecific competition among predators for prey.Throughout our investigation, we successfully established the existence criteria for biologically meaningful axial and interior equilibrium points, and we assessed their stability.Particularly, we found that when the growth rate of predators, based on the conversion coefficients from individual prey to individual predators, exceeds the death rate of predators (ζ = r 2 ρ/δ 2 > σ = δ 3 ρ/δ 2 ), only predators will survive (v * > 0).Moreover, we identified the influential role of the parameter ξ , independent of E * , in causing a Hopf bifurcation around E * .Theorem 4.3 reveals that when D > 0 and the proportion between the maximal per capita predator consumption rate and the intensity of competition among individuals of the prey species (ξ = δ 2 /ρ) exceeds the critical value ξ 0 , the prey-predator inhabitants are stable for any initial interior population.This emphasizes the significance of prey refuge and interference among predators in maintaining stable populations.However, when the ratio between the maximal per capita predator consumption rate and the intensity of competition among individuals of the prey species falls below the critical value ξ 0 , the population size becomes unstable and cannot be precisely determined.When this ratio precisely equals the critical value ξ 0 , the prey-predator inhabitant dynamics exhibit periodic changes due to the occurrence of Hopf bifurcation, as presented in Theorem 4.4.
Furthermore, we delved into the dynamics of the prey-predator population, considering both spatial and temporal motion.To capture spatial movements, we analyzed a diffusion system (5.1) and explored the diffusion-driven instability of the spatial system in detail.We observed that the stability of the interior equilibria E * could vary from stable to unstable, even when ξ > ξ 0 was satisfied due to the occurrence of diffusion-driven effects.The corresponding results are presented in Theorem 5.1.To validate our analytical findings, we provided numerical examples in Sect.6, ensuring the robustness of our results.
As we conclude our study, we recognize that an exciting avenue for future research would involve exploring the impact of stochastic noise on the model, particularly in terms of habitat-dependent parameters.Incorporating stochastic elements could offer valuable insights into the system's resilience and adaptability under varying environmental conditions.
In summary, our work contributes to a deeper understanding of the complex dynamics within prey-predator ecosystems, shedding light on how refuge and interactions among predators play pivotal roles in governing population stability and behavior.Our findings open the door to further investigations and offer potential applications in ecological management and conservation efforts.

Figure 7 Figure 8
Figure 7 Existence of a stable limit cycle around the interior equilibrium E * = (0.315572, 0.354658)

Table 1
Ecological meaning