1. Introduction
Understanding the mechanisms that enable complex, self-organized population dynamics in natural systems remains a central focus of theoretical ecology [1–3]. Omnivory, in which species feed across multiple trophic levels, is a widespread feature of both aquatic and terrestrial ecosystems [4–6], generating cross-trophic interactions that can profoundly influence population fluctuations, community persistence, and ecosystem stability [7]. A key form of omnivory is intraguild predation (IGP), in which a top predator consumes both a basal resource and a competing intermediate consumer. By combining competition and predation within a single trophic module, IGP systems can exhibit complex dynamical behaviors [8–10]. Previous studies have shown that IGP can stabilize populations through mechanisms such as adaptive foraging and morph switching [11–14], yet it can also destabilize them, producing large-amplitude oscillations or chaotic fluctuations [15–19].
In IGP systems, chaotic population dynamics are often attributed to period-doubling bifurcations [20, 21]. Yet, other mechanisms capable of generating chaos remain largely unexplored. One such mechanism is the Shilnikov homoclinic bifurcation, in which a trajectory departing from a saddle-focus equilibrium returns along its stable manifold, forming a homoclinic loop that gives rise to a Smale horseshoe structure and a countable set of periodic orbits, leading to chaotic dynamics [22, 23]. While Shilnikov-type chaos has been extensively studied in simple food-chain models [24–27], its potential role in IGP systems has not been examined.
Here, we develop a minimal model of intraguild predation [28, 29], and demonstrate that Shilnikov-type homoclinic orbits and chaotic dynamics can arise in IGP systems, giving rise to complex, self-organized population fluctuations. By combining numerical simulations, local bifurcation analysis, and numerical evaluation of Lyapunov exponents, we identify the conditions under which distinct coexistence modes emerge independently of stochastic forcing or environmental variability. The model is applicable to natural IGP communities. In particular, our results show that the simulated self-organized oscillations quantitatively reproduce empirical field data from a host-parasite system [30] across communities with differing productivity levels. Overall, this study provides new insights into how intraguild predation can drive the emergence of complex, self-organized dynamics in natural ecosystems.
2. Results
2.1. A minimal model of intraguild predation
To investigate the effects of omnivory on complex population dynamics in simple communities, we developed a minimal model of intraguild predation (equation (1 )) involving three interacting populations: a basal resource R, an intermediate predator C1, and a top predator C2 (figure 1). This three-species configuration, without additional trophic levels or interaction types, represents the simplest ecological framework in which both predation and competition coexist, with a minimal number of variables and a simple structure, while retaining the basic intraguild predation observed in complex food webs.
Figure 1. Schematic of a minimal intraguild predation model. An intermediate predator C1 feeds on the basal resource R, while the top predator C2 preys upon both R and C1, generating a dual interaction of competition and predation. Arrows indicate the direction of biomass flow. |
Following the classical MacArthur consumer-resource framework [31–33], the resource R exhibits logistic growth, $\dot{R}=rR\left(1-\frac{R}{{K}_{0}}\right)$, reflecting its intrinsic reproduction limited by the environmental carrying capacity. Meanwhile, the consumer populations C1 and C2 increase solely through predation and decline in the absence of prey. Specifically, predation by Ci reduces the prey S at a per-capita rate fi(S), and increases the biomass of Ci proportionally, scaled by conversion efficiency wi,S.
The system involves three feeding interactions, ${f}_{1}\left(R\right)$, ${f}_{2}\left(R\right)$, and f2(C1), corresponding to the arrows in figure 1. We assume Holling type-II functional responses [34] for all feeding processes, ${f}_{i}(S)=\frac{{a}_{i,S}S}{1+{b}_{i,S}S}$, which describes saturation effects in consumer feeding rates commonly observed across aquatic and terrestrial systems [35]. The population dynamics of the system can thus be described by the following equations:
$\begin{eqnarray}\left\{\begin{array}{l}\frac{{\rm{d}}R}{{\rm{d}}t}=rR\left(1-\frac{R}{{K}_{0}}\right)-\frac{{a}_{1,R}\,R}{1+{b}_{1,R}\,R}\ {C}_{1}-\frac{{a}_{2,R}\,R}{1+{b}_{2,R}\,R}\ {C}_{2},\\ \frac{{\rm{d}}{C}_{1}}{{\rm{d}}t}={w}_{1,R}\,\frac{{a}_{1,R}\,R}{1+{b}_{1,R}\,R}\ {C}_{1}-\frac{{a}_{2,{C}_{1}}\,{C}_{1}}{1+{b}_{2,{C}_{1}}\,{C}_{1}}\ {C}_{2}-{m}_{1}\,{C}_{1},\\ \frac{{\rm{d}}{C}_{2}}{{\rm{d}}t}={w}_{2,R}\,\frac{{a}_{2,R}\,R}{1+{b}_{2,R}\,R}\ {C}_{2}+{w}_{2,{C}_{1}}\,\frac{{a}_{2,{C}_{1}}\,{C}_{1}}{1+{b}_{2,{C}_{1}}\,{C}_{1}}\ {C}_{2}-{m}_{2}\,{C}_{2}.\end{array}\right.\end{eqnarray}$
Here, r and K0 denote the intrinsic growth rate and carrying capacity of the basal biotic resource. ai,S and 1/bi,S represent the attack rate and saturation coefficient characterizing consumer i's functional response to prey species S, where $S\in \left\{R,{C}_{1}\right\}$ and $i\in \left\{1,2\right\}$. The term wi,S denotes the biomass conversion efficiency, and mi represents the mortality rates of consumer i. All state variables remain non-negative under biologically feasible conditions.2.2. Emergence of a Shilnikov homoclinic orbit
The IGP system described by equation (1 ) admits five biologically relevant equilibria. The null state, ${E}_{0}=\left(0,0,0\right)$, corresponds to the extinction of all populations. The basal-resource state, ${E}_{R}=\left({K}_{0},0,0\right)$, represents the condition where only the resource persists at its carrying capacity. Two single-consumer equilibria also exist, ${E}_{1}=\left({\rho }_{1},{\xi }_{1},0\right)$ and ${E}_{2}=\left({\rho }_{2},0,{\xi }_{2}\right)$, where ρ1 = $\frac{{m}_{1}}{{w}_{1,R}\,{a}_{1,R}-{b}_{1,R}\,{m}_{1}}$, ξ1 = ${\rho }_{1}\frac{r\,{w}_{1,R}}{{m}_{1}}$ $\times \left(1-\frac{{\rho }_{1}}{{K}_{0}}\right)$, ρ2 = $\frac{{m}_{2}}{{w}_{2,R}\,{a}_{2,R}-{b}_{2,R}\,{m}_{2}}$, ξ2 = ${\rho }_{2}\frac{r\,{w}_{2,R}}{{m}_{2}}$ $\left(1-\frac{{\rho }_{2}}{{K}_{0}}\right)$, with C1 and C2 coexisting with the resource, respectively. Finally, the coexistence equilibrium, ${E}^{* }=\left({R}^{* },{C}_{1}^{* },{C}_{2}^{* }\right)$, represents the state in which all three populations persist simultaneously.
Next, we assess the local stability of the coexistence equilibrium E* by analyzing the eigenvalues of its Jacobian matrix $J\left({E}^{* }\right)$. The eigenvalues $\left({{\Lambda }}_{1},{{\Lambda }}_{2},{{\Lambda }}_{3}\right)$ satisfy the characteristic polynomial: Λ3 + σ1Λ2 + σ2Λ + σ3 = 0, where σ1, σ2, and σ3 are real coefficients determined by the system parameters. According to the Routh–Hurwitz criterion, E* is locally stable if σ1, σ2 > 0 and σ1σ2 > σ3. Violation of any of these conditions destabilizes the equilibrium, potentially leading to a saddle-node or Hopf bifurcation. A saddle-focus equilibrium arises when $J\left({E}^{* }\right)$ has one real eigenvalue ${{\Lambda }}_{1}\in {\mathbb{R}}$ and a complex conjugate pair Λ2,3 = μ ± iδ, with $\mu ,\delta \in {\mathbb{R}}$ and δ ≠ 0, where μ and δ denote the real and imaginary parts, respectively. A saddle-focus satisfying the Shilnikov conditions μ Λ1 < 0 and ∣μ∣ < ∣Λ1∣ is a prerequisite for a Shilnikov homoclinic orbit, a classical route to chaos.
Consistent with these criteria, numerical simulations confirm ${E}^{* }=\left(7.35,7.36,3.29\right)$ as a saddle-focus. The corresponding eigenvalues are Λ1 = –3.03 and Λ2,3 = 0.108 ± 0.098i, satisfying the Shilnikov conditions for a homoclinic orbit. Trajectories diverge from E* and return along its stable manifold, forming a closed homoclinic loop, consistent with a Shilnikov homoclinic orbit [36] (figure 2(b)). Increasing resource productivity r induces a transition from limit-cycle oscillations to a Shilnikov homoclinic orbit (figures 2(c)–(d)). These shifts indicate that resource enrichment destabilizes the coexistence equilibrium, generating self-sustained oscillation via the Shilnikov mechanism.
Figure 2. Emergence of Shilnikov homoclinic attractors in a minimal intraguild predation model as resource productivity r increases across different ecosystems. (a)–(b) Time series of population abundances, showing the transition from regular oscillations to Shilnikov homoclinic oscillations. (c) Oscillating coexistence trajectories from two distinct initial conditions, converging to the same attracting limit-cycle orbit under the same parameters as panel (a). (d) The corresponding 3D phase-space trajectories associated with panel (b), demonstrating the transition from a stable limit cycle to a homoclinic orbit. See SM section III for simulation details of figures 2–4. |
2.3. Shilnikov-type chaos through homoclinic bifurcations
A Shilnikov-type homoclinic bifurcation can induce chaos near the saddle-focus equilibrium E*. To characterize the chaotic dynamics of the IGP system, we performed sensitivity analysis and computed the Lyapunov spectrum (figure 3). The prey (R) and predators (C1, C2) exhibit aperiodic, irregular oscillations, providing direct evidence of chaos (figure 3(a)). A small perturbation in the initial predator abundance (ΔC1 = 10−4) causes trajectories to diverge rapidly (figure 3(b)), demonstrating sensitivity to initial conditions, a hallmark of chaotic dynamics. In the three-dimensional phase space $\left(R,{C}_{1},{C}_{2}\right)$, a single trajectory traces a horseshoe-shaped strange attractor (figure 3(c)), consistent with Shilnikov-type chaos. Lyapunov exponents reveal a positive largest exponent $\left({\lambda }_{1}=0.04\gt 0\right)$, and a negative sum of all exponents (Σiλi = –0.89 < 0), confirming exponential divergence of nearby trajectories and overall dissipative behavior (figure 3(d)). These spectral properties indicate chaos confined to a bounded strange attractor.
Figure 3. Chaotic dynamics in a minimal intraguild predation model. (a) Time series of predator and prey abundances showing irregular, aperiodic oscillations. (b) Divergence of trajectories under small perturbations (ΔC1 = 10−4), demonstrating sensitivity to initial conditions. (c) Chaotic attractor in the $\left(R,{C}_{1},{C}_{2}\right)$ phase space, showing horseshoe-type strange attractor. (d) Lyapunov spectrum $({\lambda }_{1},{\lambda }_{2},{\lambda }_{3})=\left(0.04,0,\,\unicode{x02013}\,0.93\right)$, computed using the Benettin algorithm from the time series in (a); a positive λ1 and negative sum of exponents confirm chaos. (e) Lyapunov spectrum (λ1, λ2 and λ3 ) as r varies, illustrating the onset and persistence of chaos. (f) Bifurcation diagram of local maxima and minima of C1 and C2 across the same r range as in (e), showing periodic (isolated points) and chaotic (dense vertical bands) regimes. |
To further examine the dependence of chaotic dynamics on resource productivity r, we computed the Lyapunov spectrum (Figure 3(e)), and the corresponding bifurcation diagrams of the local maxima and minima of C1 and C2 (figure 3(f)). In regions where the largest Lyapunov exponent becomes positive (λ1 > 0), the bifurcation diagrams show dense vertical bands, indicating the onset and persistence of chaos. By contrast, isolated points in the bifurcation diagrams correspond to periodic oscillations, where λ1 = 0. These observations demonstrate that the system transitions between periodic and chaotic regimes, across the range of r, with chaotic windows consistently identified by both Lyapunov exponents and bifurcation patterns (figures 3(e)–(f)).
The saddle-focus E* serves as the organizing center of the chaotic dynamics. Homoclinic bifurcations destabilize nearby orbits, producing recurrent, irregular fluctuations across trophic levels while all three species coexist. These results demonstrate that even a minimal intraguild predation model can generate chaos via a Shilnikov-type mechanism, revealing a novel pathway for complex population dynamics induced by omnivory.
2.4. Comparison of model predictions with field observations
We evaluate the model against field data from a host-parasite community [30], a representative IGP system. This community includes the specialist parasitoid Encarsia perniciosi (intermediate predator) and the facultative parasitoid Aphytis melinus (top predator), both exploiting the common host Aonidiella aurantii (California red scale) [30]. Model dynamics are examined across three productivity levels (grapefruit: low, citrus: medium, and lemon: high) by varying the resource growth rate (r), reflecting observed differences in host egg densities among cultivars [30, 46].
Across all levels, communities exhibit sustained self-organized oscillations (figures 4(a)–(c)). At low productivity, the intermediate predator (C1) dominates, whereas high productivity strengthens direct resource-omnivore interactions, elevating the mean abundance of the top predator (C2). Meanwhile, the time-averaged relative abundance of the basal resource (R) increases from 0.87 to 0.92, highlighting productivity-dependent shifts in species dominance.
Figure 4. Comparison of model results with field data from a host-parasite community [30]. (a)–(c) Simulated time series of species abundances under increasing resource productivity (modeled via parameter r). Markers indicate time-averaged abundances, and error bars represent oscillation amplitudes. (d) Comparison of model results and experimental data [30] based on time-averaged relative abundances across productivity gradients (grapefruit: low (L); citrus: medium (M); lemon: high (H)). Model-experiment correspondence is quantified using Bray–Curtis similarity, with values of 0.93 (low), 0.97 (medium) and 0.97 (high). Species are color-coded by row: R (Aonidiella aurantii) in orange, C1 (Encarsia perniciosi) in blue, and C2 (Aphytis melinus) in green, with columns indicating productivity levels marked as: low (L, triangle), medium (M, star), and high (H, sphere). See SM section IV. for parameter settings and empirically supported data [37–45] for figure 4 of A. melinus and E. perniciosi. |
Furthermore, we quantify model performance using Bray–Curtis similarities between simulated and observed time-averaged relative abundances, yielding values above 0.9 across all three productivity levels (figure 4(d)), exceeding the widely accepted ecological threshold of 0.8 [47]. This strong correspondence demonstrates that the model reliably captures community-level patterns and indicates that omnivory interactions within the intraguild predation framework are sufficient to generate ecologically realistic dynamics.
3. Discussion
Omnivory through intraguild predation is pervasive in natural ecosystems and has been recognized as a key driver of complex population dynamics [9, 29, 48]. Although mechanisms such as adaptive foraging and morphological switching that stabilize IGP systems have been extensively studied [49, 50], the role of Shilnikov homoclinic bifurcations in shaping omnivory dynamics remains largely unexplored. Here, we show that a minimal intraguild predation model can facilitate Shilnikov homoclinic bifurcations and give rise to intrinsic chaotic dynamics near a saddle-focus equilibrium. By combining numerical simulations and Lyapunov spectrum analysis, we reveal self-organized coexistence modes ranging from regular oscillations to Shilnikov homoclinic orbits and Shilnikov-type chaotic fluctuations. Furthermore, our model quantitatively explains the species coexistence patterns of a host-parasite community, which constitutes a natural IGP community observed in field studies [30].
Previous studies of IGP systems have emphasized the route to chaos via period-doubling bifurcations [20, 21, 51]. In contrast, our findings reveal Shilnikov-type bifurcations leading to chaos in IGP systems. This mechanism involves a global bifurcation in which trajectories depart from a saddle-focus coexistence equilibrium and return along its stable manifold, forming a homoclinic loop that generates a Smale horseshoe structure and a countable set of unstable periodic orbits [22, 23, 52]. It produces irregular bursting, long excursions near the equilibrium, and sensitive dependence on initial conditions, yielding transient and spectral signatures distinct from those of period-doubling cascades. Although Shilnikov-type chaos has been reported in simple food-chain systems lacking omnivory [24–27], omnivory and intraguild predation are widespread features of natural ecosystems [4]. Therefore, our results highlight the ecological relevance of Shilnikov-type dynamics in IGP systems and provide a mechanistic explanation for intrinsic population fluctuations observed in field data.
Finally, from an ecological perspective, intraguild predation plays a key role in promoting species coexistence [53]. A well-known constraint on species diversity in natural ecosystems is the competitive exclusion principle (CEP), which states that two consumer species competing for a single type of resource cannot coexist at steady state [54–58]. Although our previous studies have identified mechanisms that can break CEP, such as pack hunting [59] and intraspecific predation [60, 61], intraguild predation promotes species coexistence through a food-web mechanism that alleviates the constraints of CEP by creating contextual differences [62]. Disruption of omnivory, whether due to behavioral or environmental constraints in an IGP system, can impair species coexistence that would otherwise be facilitated (figure S1). Our results provide an explanation for sustained species coexistence under resource enrichment and offer new insights into the paradox of enrichment [63].