New results for fractional ordinary differential equations in fuzzy metric space

1.   Introduction

In recent ecological research, the intricate dynamics inherent in three-species food chain models have sparked substantial interest, prompting a deeper understanding of complex ecological relationships ^{[2,8,9,12,23,25,33,43,44]}. This fascination has driven investigations into omnivory, a phenomenon where species feed across multiple trophic levels, with particular attention directed towards intraguild predation (IGP) ^{[19,48,50,51]}. IGP, a specialized form of omnivory, involves intricate interactions between a resource (prey), an intermediate predator preying exclusively on the prey and a top predator (referred to as an intraguild predator or IG predator) that consumes both the prey and the intermediate predators ^[14,16]. Given its pronounced impact on species coexistence and the emergence of chaos, the dynamics of IGP have become a central topic in ecological studies. Panja ^[38] analyzed a fractional-order prey-predator model with intraguild predation and observed instability in the system for integer-order derivatives. In contrast, the system remained stable with fractional-order derivatives. In ^[20], the authors explored intraguild predation dynamics, focusing on the impact of interference and predation strength on the coexistence of intraguild predators and prey. Gómez-Hernández et al. ^[13] explored an intraguild predation model and found that interspecific competition consistently reduces infection rates in both species.

Previous research has unveiled various mathematical models to capture the complex dynamics of IGP and its intricate interactions within predator-prey systems. In the realm of ecological modeling, diffusion stands as a fundamental factor closely tied to real environmental processes, influencing a myriad of dynamic behaviors. A wealth of research has delved into the dynamics of diffusive versions of ecological models, shedding light on their intricate mechanisms. Notably, Chen and Wu, in their work ^[4], illuminated the role of attractive predator-taxis in driving Turing bifurcation and the consequent emergence of spatiotemporal patterns within a predator-prey model characterized by nonlinear harvesting. Conversely, they found that the presence of repulsive predator-taxis or the complete absence thereof did not give rise to such phenomena. Moreover, this study provides valuable insights into how diffusion ratios and nonlinear harvesting practices impact the spatial distribution of species and the direction of Turing bifurcation.

In ^[5], the authors undertook a study focusing on the stability of bifurcating solutions within a predator-prey diffusive model with prey-taxis, where they not only confirmed the occurrence of steady state bifurcations but also explored the influence of repulsive prey-taxis in promoting this phenomenon. This research further delved into the determination of bifurcating solution stability. These collective findings offer valuable insights into the realm of ecological dynamics and diffusion-based ecological modeling. Hastings and Powell ^[15] highlighted the potential for chaos within three-species food chain models. On the other hand, Holt and Polis ^[17] introduced a Lotka-Volterra model characterized by linear functional responses, revealing that increased intensity of IGP could induce destabilization within the system ^[26]. Research in the field of intraguild predation (IP) has led to diverse insights into predator-prey dynamics. Hsu et al. ^[18] examined IGP scenarios using Beddington-DeAngelis responses, studying coexistence and competitive exclusion. They highlighted the significance of the interference parameter in shaping stability or oscillation. Kang and Wedekin ^[21] explored type III functional responses for intermediate predators and IG predators, revealing distinctions between generalist and specialist IG predators, and identifying multiple attractors. Freeze et al. ^[10] introduced ratio-dependent functional responses, investigating persistence, permanence and stability through innovative methods. These studies have laid crucial groundwork for understanding predator-prey interactions and set the stage for our exploration of cooperative hunting effects in the context of IGP dynamics. Building on this foundation, Tanabe and Namba ^[47] explored the emergence of chaotic dynamics within IGP scenarios, omitting intraspecific competition among IG prey and predators. Recent advancements by Sen et al. ^[44] and Majumder et al. ^[25] extended these models by incorporating saturating functional responses and intraspecific competition among IG prey and predator, leading to enhanced insights into the rich dynamics such as chaos of predator-prey interactions. However, Majumder et al. ^[25] control chaotic behavior of the system through a memory effect.

Furthermore, cooperative hunting behaviors have emerged as crucial factors shaping ecological stability. A range of adaptive characteristics have evolved among predators, including sharp teeth, claws, fangs and other tools, to facilitate efficient hunting. Predators have developed various strategies such as bait pursuit, collective hunting and other approaches. For instance, carnivores ^[24], lions ^[36], wolves ^[42] and African wild dogs ^[7] frequently collaborate to capture and subdue their prey. Similarly, ants ^[27] and spiders ^[49] exhibit cooperative behaviors in hunting, drawing considerable attention. The concept of hunting cooperation was initially introduced by Cosner et al. ^[6] in a predator-prey model. He proposed a functional response for predators that engage in spatially linear hunting formations and aggregate when encountering significant prey groups. Cooperative foraging has been explored in diverse predator-prey models, revealing scenarios of both destabilization and stabilization ^{[1,3,6,11,31,32,34]}. The incorporation of cooperative hunting behavior among predators in the prey-predator model reveals novel insights into the complex dynamics of ecological systems ^[45,46]. This growing understanding of cooperative behaviors' influence on ecological dynamics provides a strong motivation for our study.

We undertake a comprehensive investigation into the interplay between hunting cooperation among IG predators and IG prey, aiming to elucidate its impact on the intricate dynamics of predator-prey interactions. Our objective is to decipher how this cooperative dimension influences the emergence of chaotic oscillations and its potential to shape the overall stability or instability of ecological communities. By incorporating this cooperative aspect into our model, we aspire to provide novel insights into the intricate interactions between predator and prey, specifically within the context of IGP and the cooperative behaviors that can arise in these ecological systems. The distinctive feature of our study, setting it apart from existing literature ^{[25,41,44,47]} addressing intraguild predation (IGP) scenarios, is the inclusion of cooperative hunting behavior within the model. Through this endeavor, we seek to contribute to the deeper comprehension of the complexities inherent in predator-prey interactions and the factors that govern their dynamics.

2.   Model formulation

We focus on a three-species Lotka-Volterra model that incorporates intraguild predation dynamics. The model involves three major species denoted as $x$, $y$ and $z$, representing a basal prey, an intermediate consumer (IG prey) and an omnivorous top predator (IG predator), respectively. We adopt a linear functional response between these species, as discussed in ^[47]. Additionally, we account for intraspecific competition within the intraguild prey and predator populations, as elaborated in ^[44]. Given the limited availability of resources, competition becomes essential for their sustenance within the community. In our model, we introduce the concept of hunting cooperation, which confers benefits to both intraguild prey and predator. Moreover, the success in capturing a basal prey and intraguild prey increases with the density of intraguild prey and predator, respectively, as highlighted in ^[37]. The model can be expressed as follows:

with non-negative initial conditions. The model parameters $r_{11}, \ d_{11}, \ \alpha_{11}, \ \alpha_{12}, \ c_{11}, \ \alpha_{13}, \ e_{11}, \ \alpha_{22},$ $\alpha_{23}, \ c_{22}, \ d_{22}, \ e_{22}, \ \alpha_{31}, \ \alpha_{33}$ and $d_{33}$ are all positive constants and their ecological descriptions are provided in Table 1. Our model is distinct from previous works like Kang and Wedekin ^[21] and Tanabe and Namba ^[47] due to the incorporation of intraspecific competitions and hunting cooperation behavior among intraguild prey and predator. Furthermore, compared to models studied by Sen et al. ^[44] and Majumder et al. ^[25], our model employs different functional response types between intraguild prey and predator, while also involving cooperative hunting.

3.   Positivity and boundedness of the solution

The positivity of state variables is a fundamental concept in ecological modeling, indicating that the populations of the involved species remain non-negative within the system. Ensuring the positivity of solutions in the proposed model (2.1) is crucial for understanding the persistence of populations.

To establish the positivity of the system (2.1), we integrate its equations and utilize the given initial conditions to obtain the following expressions:

The positivity of these solutions follows from the fact that the right-hand side of each equation is non-negative given the non-negative initial conditions. Consequently, all solutions of the system (2.1) are guaranteed to be non-negative.

While we have established the positivity of solutions, we omit the discussion of the boundedness of solutions, as the boundedness property can be similarly proven and is straightforward in this context.

4.   Equilibria analysis

An equilibrium point in a dynamical system is a state or configuration where the system's variables remain constant over time. In other words, it's a point at which the system's behavior does not change, and the rates of change of the variables are all zero. Equilibrium points are also known as fixed points or steady states. Mathematically, for a dynamical system described by a set of equations, an equilibrium point is found by setting all the derivatives of the variables to zero. Equilibrium points play a crucial role in the analysis of dynamical systems. They serve as important reference states to understand the behavior of the system around those points.

The proposed IGP system has the following feasible equilibria.

$(1)$ Population-free equilibrium $E_0(0, 0, 0)$ exists irrespective of any parametric restriction. From an ecological viewpoint, this equilibrium represents a theoretical scenario where all species involved in the model are entirely absent from the ecosystem. While such a scenario may not directly correspond to a real ecological state, it is a critical reference point.

$(2)$ Basal prey equilibrium $E_1\left(\dfrac{r_{11}-d_{11}}{\alpha_{11}}, 0, 0\right)$ exists under the parametric condition $r_{11} > d_{11}$. From an ecological viewpoint, this equilibrium reflects a scenario where the prey population is the only species present with intermediate and top predator both at zero. Moreover, the existence of $E_1$ suggests that when the prey species experiences a birth rate greater than its death rate, it can maintain its population in the absence of predators.

$(3)$ Basal prey-IG prey equilibrium $E_2(x_2, y_2, 0)$, where $x_2$ and $y_2$ are positive solution(s) of following isoclines:

We can obtain $x$ from first equation and then putting the value of $x$ in the second equation we get a cubic equation for $y$. This is little bit difficult to find the root (i.e., the component of equilibrium point). Therefore, we numerically check the existence $E_2$ by plotting the above isoclines in $x-y$ space as depicted in Figure 1. Figure exhibits unique $E_2$ for the parameter values mentioned in the Table 1. The ecological significance of $E_2$ lies in its representation of a complex ecological interaction where both basal prey and intraguild prey species coexist without the presence of intraguild predators. This equilibrium suggests the existence of resource partitioning mechanisms or the presence of ecological niches that allow both the basal prey and intraguild prey to thrive without direct predation pressure from higher trophic levels.

$(4)$ Basal prey-IG predator equilibrium $E_3\left(\dfrac{\alpha_{13}d_{33}-\alpha_{33}d_{11}+\alpha_{33}r_{11}}{\alpha_{11}\alpha_{33}+\alpha_{13}\alpha_{31}}, 0, {\dfrac{\alpha_{31}r_{11}-\alpha_{11}d_{33}-\alpha_{31}d_{11}}{\alpha_{11}\alpha_{33}+\alpha_{13}\alpha_{31}}}\right).$ $E_3$ exists if the parametric conditions $\alpha_{13}d_{33}+\alpha_{33}r_{11} > \alpha_{33}d_{11}$ and $\alpha_{31}r_{11} > \alpha_{11}d_{33}+\alpha_{31}d_{11}$ hold. Ecologically, this equilibrium illustrates the mechanisms that allow basal prey and intraguild predator to coexist without one driving the intermediate predator to extinction. By understanding the conditions that lead to the coexistence of these two species, ecologists can make more informed decisions regarding the preservation of species and ecosystem health.

$(5)$ Coexistence equilibrium $E_4({x_4, y_4, z_4})$, where the components $x_4, \ y_4$ and $z_4$ can be obtained from following isoclines.

Due to the difficulties of finding parametric values of the components of $E_4$, we find it numerically by plotting above isoclines in $x-y-z$ space as illustrated in Figure 2. It is apparent from the figure that for table parameter, $E_4$ exists uniquely. The coexistence equilibrium $E_4$ represents a critical ecological scenario within a model where all three interacting species coexist within an ecosystem. This equilibrium reflects the richness and stability of biodiversity, highlighting the ability of different species to thrive without the extinction of any species.

5.   Local stability analysis

Mathematically, local stability analysis involves linearizing the system's equations of motion around the equilibrium point. This is typically done by calculating the system's Jacobian matrix at the equilibrium, which describes how small changes in the variables affect the rates of change of the variables themselves. The eigenvalues of this Jacobian matrix are crucial indicators of the local stability properties of the system.

Local stability analysis provides valuable insights into the qualitative behavior of a dynamical system around an equilibrium point. However, it is important to note that local stability analysis only considers behavior near the equilibrium and may not fully capture more complex behaviors, such as oscillations, bifurcations and global attractors, that can arise in nonlinear systems. Here, we check local stability of the proposed system around the system's equilibria. First, we evaluate the Jacobian matrix of the system (2.1) as illustrated below.

where

$(1)$ Eigenvalues of the Jacobian matrix $J(E_0)$ are $-d_{11}+r_{11},$ $-d_{22}$ and $-d_{33}$. From this eigenvalue analysis, we can say that $E_0$ is stable under the parametric condition $d_{11} > r_{11}$. This stability condition is opposite to the existence condition of $E_1$. Thus, $E_0$ and $E_1$ are related via a transcritical bifurcation. Further, we verify the stability of $E_0$ numerically as described in Figure 3. In the figure, we choose $r_{11} = 0.05$ and $d_{11} = 0.1$. It is apparent from this figure that $E_0$ is stable as condition $d_{11} > r_{11}$ holds.

$(2)$ Eigenvalues of the Jacobian matrix $J(E_1)$ are $d_{11}-r_{11}$, $-\dfrac{\alpha_{12}d_{11}e_{11}-\alpha_{12}e_{11}r_{11}+\alpha_{11}d_{22}}{\alpha_{11}}$ and $-\dfrac{\alpha_{11}d_{33}+\alpha_{31}d_{11}-\alpha_{31}r_{11}}{\alpha_{11}}$. Hence, $E_1$ is locally asymptotically stable if the parametric conditions $d_{11} < r_{11}$, $\alpha_{12}d_{11}e_{11}+\alpha_{11}d_{22} > \alpha_{12}e_{11}r_{11}$ and $\alpha_{11}d_{33}+\alpha_{31}d_{11} > \alpha_{31}r_{11}$ hold. $E_1$ loses its stability when either of the conditions fails. We verify its local stability numerically also as described in Figure 4. The parameters used in this figure satisfy stability conditions, as a result, it is stable as depicted in the figure.

$(3)$ Jacobian matrix evaluated at $E_2$ is given by

One eigenvalue of $J(E_2)$ is $e_{22}\alpha_{23}y_2+\alpha_{31}x_2-d_{33}$ and the other two eigenvalues are the roots of the characteristic polynomial of the matrix

Now, eigenvalues of $\widetilde{J}(E_2)$ have negative real part if tr$(\widetilde{J}(E_2)) < 0$ and det$(\widetilde{J}(E_2)) > 0$. Thus, $E_2$ is stable under the following conditions

As it is difficult to follow these conditions for stability of $E_2$, we numerically verify its stability, Figure 5. The parameters used for numerical verification satisfy above stability conditions of $E_2$. As a result $E_2$ is locally asymptotically stable as described in the figure. $(5)$ The Jacobian matrix evaluated at $E_3$ is given by

where

The characteristic equation associated with the matrix $J(E_3)$ is given by

where

By Routh-Hurwitz criteria, $E_3$ is stable if $\Gamma_2 > 0$, $\Gamma_0 > 0$ and $\Gamma_2\Gamma_1 > \Gamma_0$. The expression of $\Gamma_i$, $i = 0, 1, 2$, are large to say about the stability $E_3$. So, we numerically investigate its stability as depicted in Figure 6. It is apparent from the figure that $E_3$ is stable since the parameter used here holds its stability conditions.

$(6)$ The Jacobian matrix evaluated at a representative coexisting equilibrium point

is denoted as

with $J^4_{ij}$ representing the matrix components when $J(x, y, z)$ is computed at

The associated characteristic equation for the matrix $J(E_4)$ is given by

where the coefficients are defined as follows:

The stability of $E_4$ can be determined using the Routh-Hurwitz criteria, which states that $E_4$ is stable if $\Delta_2 > 0$, $\Delta_0 > 0$ and $\Delta_2\Delta_1 > \Delta_0$. The explicit expressions for the components of $E_4$ are not obtainable due to their complexity, leading to the inability to derive parametric restrictions explicitly for the stability of $E_4$. To validate the stability of $E_4$, a numerical example is employed. The chosen parameter values in Figure 7 adhere to the Routh-Hurwitz criteria for the stability of $E_4$. Consequently, $E_4$ is demonstrated to be stable, as depicted in Figure 7.

6.   Bifurcation analysis

Bifurcations in a dynamical system are critical points where the qualitative behavior of the system changes as one or more parameters are varied. In other words, they mark a point at which the stability, periodicity, or overall behavior of the system undergoes a significant transformation. Bifurcations can lead to the emergence of new patterns, such as oscillations, chaos or the creation and destruction of equilibrium points. Bifurcations are fundamental to the study of nonlinear dynamical systems and play a pivotal role in understanding the complex behaviors that can arise in such systems. They are often associated with changes in the number or stability of equilibrium points or the appearance of periodic solutions. Bifurcation analysis involves studying how the system's behavior changes as the parameters are continuously varied. In our proposed IGP system, two types of bifurcations occur: Transcritical and Hopf, which are described in this section.

Theorem 1. The system described by (2.1) undergoes a transcritical bifurcation at the parameter value

concerning the equilibrium point $E_0(0, 0, 0)$.

Proof. Let us consider the Jacobian matrix $J(E_0)|_{r_{11}^{[TB]}}$. One of its eigenvalues is zero, and the other two are negative. Since two eigenvalues are consistently negative and the remaining one is zero, concluding stability-instability is not feasible at this point.

Let

be the eigenvectors of $J{(E_0)}$ and $J(E_0){^T}$ corresponding to the zero eigenvalue, respectively. Next, we define

The conditions outlined by Sotomayor's theorem ^[40] for a transcritical bifurcation are as follows:

$(1)$ $W^TF_{r_{11}}\left(E_0, r_{11}^{[TB]}\right) = 0$,

$(2)$ $W^T\left[DF_{r_{11}}\left(E_0, r_{11}^{[TB]}\right)V\right] = 0$,

$(3)$ $W^T\left[D^2F\left(E_0, r_{11}^{[TB]}\right)(V, V)\right] = -2\alpha_{11}\neq0$.

It's important to note that for a non-degenerate transcritical bifurcation, the second transversality condition should not be zero ^[40]. Hence, the system defined by (2.1) experiences a degenerate transcritical bifurcation ^[22] around $E_0$ when

Consequently, $E_0$ undergoes a transcritical bifurcation as the parameter $r_{11}$ crosses the critical value

and the conditions outlined in the theorem are satisfied. □ Next, we will establish the presence of a Hopf bifurcation concerning the model parameter $r_{11}$. We have designated the critical value of $r_{11}$ as $r_{11}^{[HB]}$.

Theorem 2. The system defined by (2.1) exhibits a Hopf bifurcation if the condition

is met when the parameter $r_{11}$ crosses the critical value $r_{11} = r_{11}^{[HB]}$.

Proof. Considering that the equilibrium point $E_4$ is dependent on the parameter $r_{11}$, we apply the Hopf bifurcation theorem ^[29]. The Hopf bifurcation occurs if the variation matrix $J(E_4)$ has one eigenvalue with a negative real part, where

and the other two eigenvalues are purely imaginary. For purely imaginary eigenvalues, the characteristic Eq (5.2) should satisfy $\Delta_2\Delta_1-\Delta_0 = 0$. We assume that $i\omega$ represents the purely imaginary eigenvalue corresponding to the parameter $r_{11}^{[HB]}$.

By differentiating both sides of the characteristic Eq (5.2) with respect to $r_{11}$, we obtain

Substituting $\lambda = i\omega$ and considering the real part, we arrive at

Therefore, the system described by (2.1) experiences a Hopf bifurcation concerning the parameter $r_{11}$ if the condition

is satisfied at $r_{11} = r_{11}^{[HB]}.$ □

7.   Numerical simulation

In this section, we will delve into our bifurcation analysis results. We will start by examining the bifurcation structure concerning the parameter $r_{11}$, as depicted in Figure 8.

This figure uncovers dynamic behaviors as $r_{11}$ undergoes variations. For lower values of $r_{11}$, the system exclusively features the stable equilibrium $E_0$. The presence of any other equilibria is absent within this parameter range. As $r_{11}$ increases, a critical value is reached, leading to the emergence of the stable equilibrium $E_1$. Simultaneously, the previously stable equilibrium $E_0$ shifts from stability to instability, retaining this behavior for all subsequent higher values of $r_{11}$. With further increments in $r_{11}$, another transcritical point comes into view. At this juncture, the stable equilibrium $E_2$ arises, causing the previously stable $E_1$ to transition into an unstable state. This instability of $E_1$ continues for higher values of $r_{11}$. Finally, the system encounters an additional transcritical point, resulting in the emergence of the stable equilibrium $E_4$, while the previously stable $E_2$ transforms into an unstable state.

Similarly, Figure 9 presents a transcritical bifurcation diagram involving the parameter $e_{11}$. For lower values of $e_{11}$, the system displays a stable equilibrium $E_2$. Beyond a critical value of $e_{11}$, the stable equilibrium $E_4$ comes into play. Consequently, the previously stable $E_2$ shifts into an unstable state. These analyses provide valuable insights into the system's behavior with varying parameters $r_{11}$ and $e_{11}$. They reveal critical points where equilibria undergo shifts in stability and new equilibria emerge, contributing to a deeper comprehension of predator-prey interactions within the context of changing parameter values. Transcritical bifurcation is ecologically significant as it marks a critical point in the dynamics of predator-prey relationships and can lead to shifts in population equilibrium, which, in turn, can impact the stability and biodiversity of ecosystems. Understanding these dynamics is crucial for effective ecological management and conservation efforts.

Continuing our exploration, we now delve into the system's behavior surrounding the coexistence equilibrium $E_4$. This investigation involves the examination of phase portraits and time series solutions for various values of $r_{11}$, as showcased in Figure 10. As the parameter $r_{11}$ undergoes increments, the system's dynamics around the coexistence equilibrium $E_4$ progress through distinct phases, leading to diverse patterns of behavior. For lower values of $r_{11}$, the system exhibits a stable behavior around the coexistence equilibrium, manifesting as confined trajectories in the phase portrait, Figure 10a. As $r_{11}$ increases, the system's behavior shifts to one-period oscillations. Trajectories trace out regular, repetitive patterns in the phase portrait, Figure 10b. Further increments in $r_{11}$ lead to the emergence of two-period oscillations, Figure 10c. Trajectories now complete two cycles before repeating themselves in the phase portrait. The system continues to evolve with higher values of $r_{11}$, resulting in four-period oscillations, Figure 10d. Trajectories now complete four cycles within the same phase portrait. As $r_{11}$ further increases, the system transitions into chaotic behavior, Figure 10e. Trajectories become highly sensitive to initial conditions, leading to unpredictable and complex patterns in the phase portrait. Moreover, we observed that the chaotic attractors size increases on increaing $r_{11}$ (see Figure 11a). Finally, for larger values of $r_{11}$ both the predator populations extinct making the axial equilibrium point $E_1$ stable (see Figure 11b). In essence, Figure 10 vividly portrays the evolution of the system's dynamics around the coexistence equilibrium $E_4$ in response to changing values of $r_{11}$. This comprehensive visualization allows us to discern the progression from stable behavior to chaotic complexity and highlights the intricate nature of predator-prey interactions in the context of varying parameter settings. Chaotic dynamics in predator-prey systems can lead to irregular and unpredictable population cycles. This unpredictability can help maintain biodiversity by preventing any one species from dominating the ecosystem for extended periods, fostering coexistence and preventing extinctions.

Moving forward, our analysis takes us to the exploration of species' maximum and minimum densities while systematically varying a key model parameter, as detailed in Figure 12. The insights gleaned from this visualization provide a comprehensive understanding of how the system's behavior evolves in response to specific parameter alterations. Figure 12a provides a revealing depiction of the system's response to changes in the parameter $r_{11}$. Initially, the system sustains a stable coexistence state. This phase occurs after a transcritical bifurcation point where a switch occurs between equilibrium states. Then the system enters a phase of one-period oscillations, a behavior characterized by regular and predictable cycles. As $r_{11}$ continues to increase, the system progresses into a phase of many-period oscillations. This pattern is marked by an increasing number of cycles within each oscillatory pattern, indicating a higher level of complexity in the population dynamics. Ultimately, the system experiences a transition into a chaotic state. In this phase, trajectories in the system become highly sensitive to initial conditions, leading to unpredictable and irregular behavior. This chaotic behavior underscores the intricate interactions that can arise within an ecosystem. Similar patterns of behavior are observed for the parameter $\alpha_{23}$ (Figure 12f), with the primary distinction being the absence of an initial transcritical bifurcation. Furthermore, we observe that parameter $c_{22}$ (Figure 12b) influences the system's dynamics by inducing two cycles, leading to an increase in the number of cycles and finally chaos. The presence of chaos and transitions between stable states and periodic states can reflect the system's ability to adapt to changing environmental conditions. Ecosystems that exhibit such behavior may be more resilient, can support biodiversity, are capable of adjusting to disturbances and evolving in response to environmental shifts. On the other hand, parameters $\alpha_{11}$ (Figure 12c), $\alpha_{22}$ (Figure 12d), $d_{11}$ (Figure 12e) and $\alpha_{12}$ (Figure 12g) exert control over the chaotic behavior by modulating the transition from many-period oscillations to one-period oscillations, bringing about a more ordered pattern and ultimately stable coexistence. The ability to control chaotic behavior through specific parameters is ecologically significant because it offers opportunities for managing, restoring and preserving ecosystems. Lastly, Figure 12h demonstrates that the parameter $c_{11}$ is capable of driving the system towards more orderly behavior from an irregular behavior, as it produces a limit cycle through multiple cycles.

Finally, we present the two-parameter bifurcation structure, varying the IG prey cooperation during hunting basal prey parameter $(c_{11})$ and Consumption rate of IG predator to IG prey parameter $(\alpha_{23}$, as shown in Figure 13. This figure reveals the division of the entire parameter space into five subregions, each exhibiting distinct dynamics. The two-parameter bifurcation structure provide a comprehensive view of the model's behavior under varying two parameters, highlighting the occurrence of stable regions, limit cycles, two period, four period and chaotic dynamics. We note that for low values of consumption rate of IG predator, the cooperation of IG predator helps the system to keep in a stable state position. However, increasing in the consumption rate of IG predators push on the system into chaotic oscillations through periodic-doubling oscillations.

8.   Conclusions

In conclusion, we delve into a three-species Lotka-Volterra model embodying intraguild predation dynamics. The model incorporates complex interactions between a basal prey, an intraguild prey and an omnivorous top predator. These interactions are defined by linear functional responses, alongside considerations for intraspecific competition and cooperative hunting, both of which play essential roles in shaping population dynamics within the ecological community. The proposed model diverges from prior works through its integration of intraspecific competition and cooperative hunting. We analyze equilibria and parametric conditions under which they are stable. Numerical stability of all steady states are justified. Our exploration encompasses a comprehensive bifurcation analysis that reveals intricate behaviors driven by parameter variations.

The examination of bifurcation structures, particularly concerning parameters such as $r_{11}$ and $e_{11}$, showcases critical shifts in equilibria stability and the emergence of new equilibrium states. We uncover transitions from stable to unstable states and vice versa, offering insights into the interplay between predator-prey interactions and varying parameter values. Moreover, our investigation into the dynamics around the coexistence equilibrium $E_4$ illuminates a diverse range of behaviors as parameter $r_{11}$ evolves. These behaviors include stable coexistence, one-period oscillations, multi-period oscillations and ultimately chaotic dynamics. This detailed analysis provides a nuanced understanding of the system's response to changing parameter values.

Furthering our comprehension, the exploration of species' maximum and minimum densities in response to parameter variations demonstrates distinct patterns. These patterns range from stable coexistence to chaotic behavior, each induced by different parameter influences. Parameters such as $\alpha_{12}$, $c_{11}$ and $d_{11}$ exert significant control over the system's chaotic dynamics, promoting stable or orderly behavior. Strong hunting cooperation behavior of intraguild predator produces chaos in the system whereas the same behavior of intraguild prey can control the chaotic nature of the system. In contrast to the findings of Sen et al. ^[44], our research demonstrates that a chaotic regime can emerge from a state of stable coexistence when considering the consumption rate of intraguild predators towards intraguild prey it may be due to the assumption of hunting cooperation. However, Vandermeer ^[50] observed chaotic oscillations in an intraguild predation (IGP) system in the absence of self-regulation for the same parameter. A recent study by Namba et al. ^[35] highlighted the stabilizing influence of strong self-regulation on a similar tri-trophic food web featuring type-I functional responses. Sen et al. ^[44] reported that robust self-regulation within both IG prey and predator populations can drive the dynamics of the system from a chaotic state to stable steady-state coexistence due to the saturated consumption of IG prey by IG predators. In our research, we extend this concept by revealing that strong self-regulation has the potential to guide the system's dynamics from a chaotic regime to stable coexistence. This stabilization is facilitated by the phenomenon of hunting cooperation. Also for low values of consumption rate of IG predator, the cooperation of IG predator helps the system to keep in a stable state position but increasing in the consumption rate of IG predators push on the system into chaotic oscillations through periodic-doubling oscillations.

In summary, this study offers a comprehensive understanding of complex predator-prey interactions in the context of intraguild predation dynamics. The integrated effects of intraspecific competition and cooperative hunting, coupled with systematic parameter variations, provide valuable insights into the intricate behaviors exhibited by ecological systems. Our findings converge as well as diverge from the existing literature. By shedding light on the intricate interplay between species and parameters, this research contributes to the broader field of ecological modeling and enhances our understanding of the complex dynamics that characterize natural ecosystems.

Use of AI tools declaration

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

Acknowledgments

This work has been supported by the Deanship of Scientific Research at King Khalid University through the Research Group program with Grant Number RGP2/202/44. The authors are very grateful to the editor and reviewers for their careful reading of the manuscript.

Conflict of interest

The authors declare that they have no conflicts of interest.