Complex dynamics in a slow–fast modified Leslie–Gower model: the effects of predator migration under central place foraging framework

1 Introduction

Ecological systems are inherently complex due to nonlinear interactions between species and continuous adaptation to changing environmental conditions. For years, mathematical modeling has been considered one of the important tools for understanding complex processes such as population persistence, extinction, and oscillation, which are highly influenced by interspecific interactions and environmental heterogeneity (Lotka, 1925; Volterra, 1926). The fundamental predator–prey models rely on simple assumptions, homogenous habitats, and instantaneous biological responses, which may not efficiently reflect the ecosystem dynamics in the real world.

The natural ecosystem is often divided into a number of interconnected habitats or zones, in which species move between them as a result of the availability of resources, predation, or the impact of the environment. As a result, patchy predator–prey models have attracted considerable attention as spatial movement can significantly impact population stability and coexistence. According to several studies, migration alone can either stabilize or destabilize the population dynamics, depending upon whether it is density dependent, asymmetric, or occurs over a range of timescales (Cantrell and Cosner, 2004; Cosner, 2009). In particular, two-patch systems provide a minimal, yet powerful enough, framework for investigating how spatial structure influences ecological outcomes.

There is a crucial component of predator–prey interactions that governs the rate at which predators consume their prey, known as functional response. Even though all Holling-type functional responses are quite effective mathematically, they fail to include biological mechanisms such as predator interference among themselves during predation. The Beddington–DeAngelis functional response addresses this limitation by accounting for mutual interference among predators, and it has been shown to generate richer dynamical behaviors compared to classical formulations (Beddington, 1975; DeAngelis et al., 1975). Furthermore, the growth of predators is often restricted by environmental constraints such as alternative food availability. This is why the modified Leslie–Gower type formulations are used to describe predator dynamics under prey scarcity (Leslie, 1948; Aziz-Alaoui and Okiye, 2003).

In recent years, the ecological models incorporating multiple time scales have captured significant attention. It is often observed in biological systems that certain biological mechanisms, like migration, behavioral switching, etc., occur much faster than demographic processes such as reproduction and mortality. This separation of timescales leads to slow–fast dynamical systems where standard analytical approaches may fail to capture the long-term behavior of the system. Through the aggregation method, fast variables are replaced with their quasi-equilibrium distributions, which results in a lower dimensional system in which slow dynamics can be accurately described (Auger and Roussarie, 1987; Mchich et al., 2005; Auger et al., 2012; Manoharan and Rani, 2026). This method has been successfully applied to a variety of predator–prey models with fast migration that gives both analytical tractability and biological interpretation (Chiorino et al., 1999; Auger et al., 2000; Mchich et al., 2005; El Abdllaoui et al., 2007; Rani et al., 2019; Manoharan et al., 2025a; Manoharan et al., 2025b).

In earlier studies, such as that of Chiorino et al. (1999), the studied models focused on prey migration between heterogeneous zones while the predator remains stationary in one zone. They had used the Lotka–Volterra functional response to show the interaction between species. Notably, they have taken asymmetric prey migration between two spatially heterogeneous zones, and later, other researchers have extended these models by including factors like density-dependent migration (Mchich et al., 2005) and reinforced repulsion mechanism (Auger et al., 2000).

In recent studies, many authors have studied the movement of species between spatially heterogeneous zones. X. Huang et al. (2023) have studied artificially controlled migration, which helps to control two populations in a stable co-existence state. In another study, K. Sato et al. (2015) have studied the pair approximation (PA) theory that while migration can help predators to persist over a wide range of parameters, prey migration can cause trouble for the prey.

Although there are advances occurring in modeling with migration, many studies depend on simplified predator–prey interactions even in a heterogeneous environment. In this context, functional response plays a crucial role in advancing predator–prey dynamics. Tahani Al-Karkhi et al (Al-Karkhi and Abukoglu, 2024) have studied a spatiotemporal model with Beddington–DeAngelis functional response showing an interaction between zooplankton and phytoplankton. S. Khajanchi et al (Khajanchi, 2014) investigated a stage-structured model with Beddington–DeAngelis-type functional response, where the model admits periodic solutions.

Another important aspect of predator–prey interaction is the dependence of predator growth on the availability of prey. Some models assume that the carrying capacity of the predator is directly proportional to the prey population. Under conditions like severe scarcity of food, predators may switch to alternative food; consequently, due to the limited supply of their preferred food, their growth becomes restricted. Thus, researchers have proposed a modified Leslie–Gower model to capture the ecological dynamics more efficiently (Gakkhar and Singh, 2012; Singh and Gakkhar, 2014; Kaur et al., 2021; Al-Momen and Naji, 2022).

In recent years, several studies have examined the dynamics of complex predator–prey systems incorporating patch structure, modified functional responses, and two-timescale analysis. There has been evidence that studies of predator–prey models with density-dependent migration, fear effects, and environmental heterogeneity lead to the emergence of a Hopf, Hopf–Hopf, and Bogdanov–Takens bifurcations. These findings emphasize the importance of spatial behavioral mechanisms in shaping ecosystem dynamics (Manoharan et al., 2025a; Kumari and Menon, 2025; Pal et al., 2025). Frequently, populations are observed to undergo oscillations, multi-stability, and sudden regime shifts, which are caused by these bifurcations.

Even though recent advances in spatial ecology have extensively examined the movement of species between heterogeneous zones, the majority of these mathematical models tend to focus on prey migration (e.g (Chiorino et al., 1999)) or assume that both species can freely migrate across patches (Manoharan et al., 2025a). However, there is a critical spatial and behavioral mechanism overlooked by these classical frameworks: the central place foraging (CPF) (Burke and Montevecchi, 2009; Bouchet et al., 2026). During critical life stages, such as the breeding season, many marine predators (e.g., seabirds like puffins) cannot freely migrate. They are instead geographically anchored to a terrestrial breeding colony and must travel rapidly, periodically, to geographically distinct feeding grounds. On the other hand, although the CPF is deeply used in empirical ecology, it is rarely incorporated into continuous dynamical systems to evaluate how these fast-scale, daily foraging decisions affect slow-scale macro-population shifts. This behavior of hunting and resting of predators is explained in Adimy and Banerjee, 2025.

To bridge this significant interdisciplinary gap, the present study proposes a three-dimensional slow–fast predator–prey system specifically designed using the CPF framework. This framework fundamentally differs from existing slow–fast spatial models. While traditional models typically assume either that both species migrate freely between patches (Manoharan et al., 2025a) or that prey migrate while predators remain stationary (Chiorino et al., 1999), our CPF formulation imposes a strict asymmetric spatial constraint: the prey are stationary, but the predators are geographically anchored to a breeding colony and must commute. Furthermore, this study introduces a unique prey-density-dependent return migration rate (m^(x)=αx). In this way, the biological reality that a central place forager’s provisioning time is influenced by local prey availability is effectively translated into a mathematical framework. By mathematically translating these micro-scale foraging constraints into a continuous slow–fast dynamical system, this study reveals fundamentally new ecological understandings. Specifically, it demonstrates that the interplay between environmental prey growth (r₁) and the predator’s outward foraging effort (m) under these strict commuting constraints acts as the primary driver of macro-scale population collapse and bistability trap phenomena that previous continuous models have primarily attributed to free-ranging spatial avoidance. As a result, this framework demonstrates that the application of time-scale separation in this context is not merely a mathematical assumption, but rather a fundamental biological reality of central place foragers, which yields novel, predictive dynamical insight. This model provides a realistic framework for understanding how these fast-scale spatial commuting constraints couple with a Beddington–DeAngelis functional response and modified Leslie–Gower growth to trigger complex bifurcations and sudden regime shifts.

The structure of the remaining part is as follows. Section 2 presents the complete slow–fast model, and its biological assumptions are formulated. In addition, this section contains key biological feasibility conditions of the aggregated model, such as positivity, boundedness, and permanence of the solutions. Section 3 analyzes the existence and stability of all biologically feasible equilibrium points of the aggregated system. In Section 4, a global sensitivity analysis is performed to find the most influential parameters for the aggregated model. Section 5 is devoted to the analysis of the codimension-one and codimension-two bifurcations. The concluding remarks are discussed in Section 6.

2 Model formulation

This section presents the biological assumptions along with the mathematical formulation for the predator–prey model. This section begins with the classical predator–prey model and subsequently incorporates spatial structure, predator migration, non-linear functional responses, and time-scale separation.

Basic predator–prey framework: In general, a predator–prey dynamics can be described by the following system of autonomous differential equations:

(1) dxdt=xG(x)−cp(x)y,dydt=ep(x)y−dy,

where x(t) and y(t) denote the population densities of prey and predator, respectively. The function G(x) denotes the growth of prey in the absence of predation, and d denotes the natural mortality rate of the predator. The term cp(x)y, denoted as the functional response, represents the rate of prey consumption by the predator. Similarly, ep(x)y represents the numeric response, which accounts for the growth of the predator by consuming prey.

Extension to spatial structure with predator migration: The system (Equation 1) can be further developed by the assumption that one of the species (say, predator) is migrating between two spatially heterogeneous zones while the prey remains only in one zone, and thus predation is occurring in only one zone (Chiorino et al., 1999).

(2) dxdt=xG1(x)−cp(x)y1,dy1dt=(My2−M^y1)+ep(x)y1−dy1,dy2dt=(M^y1−My2)+G2(y2).

Here y₁(t) and y₂(t) represent predator densities in zones 1 and 2, respectively, and y(t) = y₁(t) + y₂(t) as total predator density. In addition, x(t) denotes the predator species in zone 1. Subsequently, M represents the migration of predators from zone 1 to 2, and M^represents the migration from zone 2 to 1. Since the predators are not migrating, the predation takes place only in zone 1 and is represented in the first equation of Equation 2.

Slow–fast system: Unlike other studies [e.g (Chiorino et al., 1999)], which specifically evaluate the effect of prey migration between two heterogeneous zones, whereas predators do not migrate, the present study incorporates another ecological case. The other scenario in this study is that the predators are migrating between zones, such as between the foraging zone (for food) and the resting zone (for shelter), whereas the prey do not migrate and remain confined to their resource-rich zone. This kind of behavior is observed near bird–fish ecology, where predators migrate to their foraging zone for food, and after successfully capturing fish, they return to their breeding colony (Corkhill, 1973; Elliott et al., 1977; du P Bothma and Le Riche, 1984; Scheel and Packer, 1991; Seyrling et al., 2022). A schematic representation of this phenomena is represented in Figure 1. The reasons for returning to the resting zone can be assumed as follows:

• Load-size constraints and provisioning: Once a sufficient load of prey is successfully captured (often multiple capelin held crosswise in the beak), the puffin reaches its maximum capacity to carry food and must return to provision its chick (Corkhill, 1973).
• Avoidance of kleptoparasitism: The pelagic foraging zone exposes the seabird predators, such as puffins, to high risks of kleptoparasitism from larger avian predators, such as skuas and gulls. Returning to respective burrows provides a safe refuge for delivering food (Corkhill, 1973).
• Obligate terrestrial breeding: Although seabirds like puffins are adapted to live in the open ocean, the biological processes of mating, egg incubation, and chick rearing require a stable nesting location on land (Corkhill, 1973).
• Burrow maintenance and defense: The breeding colony requires constant presence not only for feeding but also to maintain the structural integrity of the burrow and defend it from intraspecific competition.

The interaction between prey and predator is modeled using a Beddington–DeAngelis functional response (Khajanchi, 2014), given by:

(3) cp(x,y)=cx1+ax+by,

where a, b, and c are non-negative parameters. Here c describes the capture rate of the predator and a describes the handling time of the predator per prey. The intensity of mutual interference among the individuals of the predators is reflected by the term by(b ≥ 0). Furthermore, the assumptions are as follows:

• The prey population is assumed to grow logistically with carrying capacity K and is given by the growth function G(x):
• The net demographic change of the predator population follows a modified Leslie–Gower formulation (Singh and Gakkhar, 2014; Al-Momen and Naji, 2022; Manoharan et al., 2025a),. In this spatially divided system, it is important to distinguish between local reproduction and overall demographic recruitment. While local reproduction (e.g., hatching eggs) occurs strictly in the breeding colony (resting zone), the predator populations in both zones are dynamically supplemented by regional immigration from outside the local system—for example, seabirds frequently exhibit “prospecting” behaviors, where non-breeding birds travel long distances to join active foraging flocks (in the foraging zone) or an alternative breeding colony (resting zone) in search of food and shelter, respectively (Corkhill, 1973; Péron and Grémillet, 2013). Therefore, the recruitment parameters r₂ and r₃ represent composite demographic influxes (regional immigration to zone 1; local birth, along with immigration in zone 2). To capture this distributed demographic dynamic, the formulation incorporates alternative food resources K₁ (representing survival fallbacks when the preferred prey is scarce) denoted by:
• Considering migrations between zones occurring on a faster time scale compared to species interaction (Chiorino et al., 1999; Auger et al., 2000; Dao Duc et al., 2008), the model is given as: where ϵ is a dimensionless parameter with 0 < ϵ << 1.

Considering all of the above assumptions and the Equations 3, 4 and 5, the following system is obtained:

(6) dxdt=r1x(1−xK)−cxy11+ax+by1,dy1dt=(my2−m^(x)y1)ϵ+y1(r2−a1y1x+K1),dy2dt=(m^(x)y1−my2)ϵ+y2(r3−a2y2x+K1),

with initial conditions:

(x(0),y1(0),y2(0))∈ℝ+3.

Here y₁(t) and y₂(t) denote the predator (e.g., seabird) densities in zone 1 (foraging grounds) and zone 2 (breeding colony), respectively. Other biological interpretations of parameters are given in Table 1.

Table 1. Biological representation and nominal values of parameters used in the model (Equation 13).

Furthermore, the prey denoted by x(t) is assumed to be the primary, high-energy schooling fish (e.g., capelin) preferred by the seabird population, denoted by y(t) = y₁(t)+y₂(t). This preference is driven by the high lipid content and dense schooling nature of the prey, which optimizes the predator’s foraging efficiency (Burke and Montevecchi, 2009). However, the predator also relies on alternative food resources provided by the pelagic environment (such as zooplankton or less optimal forage fish) denoted by K₁ (Corkhill, 1973). The prey exists entirely within the pelagic foraging grounds and does not migrate, representing a localized resource patch, while the seabirds act as central place foragers commuting between the foraging grounds and the breeding colony.

The outward migration rate of predators from the breeding colony to the foraging grounds, denoted by m, is assumed to be constant. This outward flight is driven intrinsically by adult metabolic needs and the continuous provisioning demands of the chicks, operating independently of the actual prey density at sea. Therefore, m is considered independent of x. However, the return migration rate from the foraging grounds back to the colony, denoted by m^, is highly dependent on prey availability. As central place foragers must secure a sufficient load of prey (e.g., a full beak load of fish) before returning to provision offspring, a higher prey density (x) significantly reduces search and handling time. Consequently, the more abundant the prey, the quicker the seabird can complete its foraging trip and return to the colony (Corkhill, 1973). Hence, the return migration rate is functionally dependent on prey density and can be defined as m^=m^(x)=αx, where > 0 is a constant. However, one can take other functional forms such as m^(x)=α1x+α2orm^(x)=α1x2+α2, where α₁ and α2∈ℝ, etc.

Furthermore, under extreme environmental conditions where the primary pelagic food stock completely collapses (i.e., x(t) = 0), seabirds are known to abandon their breeding attempts to ensure adult survival. In this model, this biological reality is captured by assuming that the predators abandon the terrestrial colony and remain entirely at sea searching for alternative sustenance. Thus, m^(x)=0 at x=0. Consequently, as the prey grows logistically with carrying capacity K, the return migration m^(x) is assumed to be bounded above by αK. Considering these assumptions, the return migration rate of the predator from the foraging grounds to the colony is given as follows:

m^(x)={αx,0<x≤K,0,x=0.

Let τ be the fast time scale associated with the rapid daily commuting flights, given by τ=tϵ, where t represents the slow demographic time scale of population growth. Since ϵ is a small dimensionless parameter (0<ϵ≪1), taking the limit as ϵ→0 isolates the fast behavioral dynamics. The system (Equation 6) can thus be rewritten in terms of the fast time scale as:

(7) dxdτ=ϵ(r1x(1−xK)−cxy11+ax+by1),dy1dτ=(my2−m^(x)y1)+ϵy1(r2−a1y1x+K1),dy2dτ=(m^(x)y1−my2)+ϵy2(r3−a2y2x+K1),

with initial conditions:

(x(0),y1(0),y2(0))∈ℝ+3.

2.1 Aggregated model

Consider that migration between two zones is conserved, and the total predator density is given by y = y₁ + y₂. Considering the non-perturbed part of Equation 7 with ϵ = 0, thus the fast subsystem is given by:

(8) dxdτ=0,dy1dτ=my2−m^(x)y1dy2dτ=m^(x)y1−my2.

From the last two equations of Equation 8, dy1dτ=0 and dy2dτ=0 give the fast equilibrium in the following Equation 9:

(9) y1=mm+m^(x)y, y2=m^(x)m+m^(x)y,

and the critical manifold is defined in the following Equation 10:

(10) M0={(y1,y2)|my2=m^(x)y1},

where x(τ) is fixed and acts as a constant in fast timescale and thus in the critical manifold. It can be observed that the critical manifold M₀ is a straight line in the y₁y₂− plane. To verify the hyperbolicity of M₀, consider the Jacobian of this fast-subsystem as follows:

(11) Jfast=[−m^(x)mm^(x)−m]

The eigenvalues of Equation 11 are given by:

(12) λ1=0 & λ2=−(m+m^(x))<0.

From Equation 12, eigenvector corresponding to eigenvalue λ1=0 is Vλ1=[mm^(x)1], which is a tangent to critical manifold M₀. This makes λ₁ a trivial eigenvalue, and thus it does not affect the hyperbolicity of the manifold M₀ (Fenichel, 1979).

Again, the eigenvector corresponding to eigenvalue λ2=−(m+m^(x)) is Vλ2=[−11], which is normal to critical manifold M0. This makes λ2 a non-trivial eigenvalue, and thus the hyperbolicity of M0 is dependent on λ2. Since λ2=−(m+m^(x))<0, M0 is always normally attracting and never loses its hyperbolicity. Thus, there exists no contact point on the critical manifold M₀, and hence dynamics like a slow–fast Hopf point or slow–fast saddle point cannot be discussed for this particular model (Yao et al., 2024).

However, with the help of center manifold theorem, the system can be approximated to an aggregated system, where the errors between the solutions of the initial system (Equation 7) and the new aggregate system are of order ϵ (see Supplementary Appendix A), i.e., for every combination of the slow variable values x and y, the fast-subsystem (Equation 8) tends toward an equilibrium state. Thus, the dynamics of the slow–fast system (Equation 7) can be observed through the method of aggregation (Chiorino et al., 1999; Auger et al., 2000). Using these fast equilibrium points y1∗ and y2∗ in place of y₁ and y₂ in Equation 7, the aggregated model is given by:

(13) dxdt=r1x(1−xK)−cxy1+B(x)x+by=xf(x,y),dydt=y(R(x)−A(x)yx+K1)=yg(x,y).

x(0)>0,y(0)>0,

where:

t=ϵτ, B(x)=α+am+aαxm, R(x)=r2f1(x)+r3f2(x), A(x)=a1f12(x)+a2f22(x),

f1(x)=mm+m^(x), f2(x)=m^(x)m+m^(x).

It is important to note that the biological implications of the fast-scale time separation do not vanish during this quasi-equilibrium reduction. The micro-scale behavioral parameters governing the daily commute, specifically the outward foraging effort (m) and the return rate parameter (α), survive the aggregation and become permanently embedded instead within the nonlinear structure of the reduced two-dimensional system. Consequently, the reduced model fundamentally retains the essential consequences of the central place foraging constraint, allowing us to evaluate how daily micro-scale decisions dictate macro-scale demographic survival.

2.2 Positivity

Theorem 1. All of the solutions x(t) and y(t) of system (Equation 13) are positive ∀t > 0, for x(0) > 0 and y(0) > 0.

Proof.

(14) x(t)=x(0)exp ∫0t(r1−r1xK−cy1+B(x)x+by)dt,

(15) y(t)=y(0)exp ∫0t(R(x)−A(x)yx+K1)dt.

From the above Equations 14 and 15, it can be verified that all of the solutions of system (Equation 13) with initial conditions x(0) > 0 and y(0) > 0 are positive. Therefore, any solution from the first quadrant of xy−plane is positive.

2.3 Boundedness

Theorem 2. All solutions of system (Equation 13) that start in ℝ+2 are ultimately uniformly bounded and eventually enter the region:

(16) S(x,y)={(x,y)∈ℝ+2:0<x≤K,0<y≤L}

Proof. From the first equation of Equation 13,

dxdt≤r1x(1−xK)

Standard comparison implies that:

lim supt→∞ x(t)≤K.

Next, considering the last equation of Equation 13:

(17) dydt=y(R(x)−yA(x)x+K1)dydt≤y(L1−yL2),

where:

L1=supx∈[0,K]R(x)=Rmax , L2=infx∈[0,K]A(x)x+K1=Amin K+K1.

In addition to determining L₁, R′(x) is given by:

R′(x)=(r3−r2)mα(m+αx)2,

where R′(x) denotes the differentiation of R(x) with respect to x. Depending on the values of r₂ and r₃, the maximum of R(x) is summarized in Table 2.

Table 2. Maximum of R(x) on (0, K).

Hence, L1=max {r2,r2m+r3αKm+αK} is a constant. Subsequently, in order to determine L2, A′(x) is given by:

A′(x)=2mα(αxa2−ma1)(m+αx)3,

where A′(x) denotes the differentiation with respect to x, which vanishes at:

xcritical=ma1αa2.

Since A′(x) changes sign from negative to positive at xcritical, the function A(x) attains a minimum at this point. If 0 < x_critical < K, then:

Amin =a1a2a1+a2.

Otherwise, A(x) is monotone on [0, K], and thus:

Amin=min{A(0),A(K)}=min {a1,a1m2+a2α2K2(m+αK)2}

Hence, L2=Amin K+K1 is a positive constant. Next, solving inequality (Equation 17) yields:

y(t)≤L1L2(CeL1t1+CeL1t),

for some constant C > 0. Therefore,

lim supt→∞y(t)≤L, whereL=L1L2.

Combining the bounds for x(t) and y(t), it follows that all solutions of system (Equation 13) with positive initial conditions eventually enter and remain in the region:

S(x,y)={(x,y)∈ℝ+2:0<x≤K,0<y≤L}.

2.4 Permanence

The system Equation 13 is said to be permanent if there exists positive constants m1,m2,M1, and M2, such that every positive solution (x(t),y(t)) of (13) with initial conditions (x0,y0)∈ℝ+2 satisfies:

m1≤lim t→∞inf x(t)≤lim t→∞sup x(t)≤M1,

m2≤lim t→∞inf y(t)≤lim t→∞sup y(t)≤M2.

Theorem 3. The system (Equation 13) is said to be permanent if:

r1>cb

Proof. From the first equation of Equation 13 and Equation 16, one can obtain that:

dxdt≤r1x(1−xK),

lim t→∞sup x(t)≤max {x(0),K}≡M1

Similarly, from the second equation of Equation 13,

dydt≤y(Rmax −yAmin K+K1),

lim t→∞sup y(t)≤max {y(0),L}≡M2,

Next, from the first equation of Equation 13 to establish the lower bounds, note from the first equation of Equation 13 that:

dxdt≥x(r1−r1Kx−cb)=x((r1−cb)−r1Kx).

If r1>cb, then comparing it with a logistic equation implies:

lim inft→∞x(t)≥K(br1−c)br1≡m1>0.

Similarly, from the second equation:

dydt≥y(Rmin −yAmax 0+K1)

lim t→∞inf y(t)≥K1RminAmax ≡m2,

where:

Rmin =minx∈(0,K)R(x), Amax =maxx∈(0,K)A(x).

Since R(x) is monotonic on (0, K) and A(x) attains its extremum at the endpoints, the constants R_min and A_max are well defined and positive. Hence,

lim inft→∞y(t)≥K1Rmin Amax ≡m2>0.

Therefore, it can be observed that whenever the condition r1>cb holds, all solutions of system (Equation 13) with positive initial conditions remain bounded away from zero and infinity. Ecologically, (cb) represents the maximum per-capita predation rate that the prey can experience due to mutual interference among predators. Thus, whenever the intrinsic growth rate of prey (r₁) exceeds this maximum predation rate (cb), the prey will persist in the ecosystem alongside the predator.

3 Equilibrium points and stability

This section thoroughly discusses the existence of equilibrium points and their nature (stable, unstable, etc.) with the help of the Jacobian matrix.

The system (Equation 13) admits the following four biologically feasible equilibrium points:

• The trivial equilibrium point P₀(0, 0) always exists.
• P1(x¯,0) = (K, 0) is the predator-free equilibrium point.
• P2(0,y¯)=(0,r2K1a1) is the prey-free equilibrium point.
• The interior equilibrium point P3(x,y) exists, where x* and y* can be determined from the following equations:

(18) y=r1(K−x)(1+B(x)x)cK−r1b(K−x*),

(19) A1x5+A2x4+A3x3+A4x2+A5x*+A6=0.

where Ai's are given by:

A1=r1a2α2(α+am)m,

A2=(1−aK)r1a2α3+(aa2−br3)r1mα2m,

A3=r1a1aαm(1−K)−Kr1a2α2+(r2+r3)(c−br1)mαK−r1r3bα2KK1−r1r2bm2−(r2+r3)r1mαbK1,

A4=(1−aK)r1a1αm+(aa1−br2)r1m2+(cr3K1−r1a2−r1r3bmK1)Kα2+m(r2+r3)(cαK−r1bαK−r1bK1),

A5=(1−aK)r1a1m2−(1−cr3K1)mαK+(c−br1)Kr2m2+(r2+r3K1)r1mαbK−r1r2bm2KK1,

A6=cr2m2KK1−a1r1m2K−br1r2m2KK1.

The sufficient condition for the existence of at least one real positive value of x=x* is A6<0. Consequently, because x<K,it follows that y>0. Hence, the positive interior equilibrium P3(x,y) from Equation 18 exists if:

(20) cr2K1a1+br2K1<r1

Furthermore, observe that:

(21) cr2K1a1+br2K1<cb

Hence, r1>cb (from Equations 20 and 21) provides a sufficient condition for the existence of at least one positive interior equilibrium P3(x,y). Ecologically, whenever the intrinsic growth rate of the prey (r1) exceeds the maximum per-capita predation rate (cb) that the prey experience due to mutual interference among predators, a biologically feasible coexistence equilibrium exists.

With the help of Descartes’ Rule of Signs (Perko, 2013), the possible number of positive roots of Equation 19 can be predicted, and this will further provide the number of existing interior equilibrium points as represented in Tables 3, 4. In Figures 2, 3, the nullclines of system (Equation 13) show the existence of interior equilibria.

Table 3. Possibilities of solutions using Descarte's Rule of sign changes when A₆<0.

Table 4. Possibilities of solutions using Descarte's Rule of sign changes when A₆>=0.

The Jacobian matrix of system (Equation 13) at any point (x, y) in the following Equation 22:

(22) J|(x,y)=[j11j12j21j22]

where:

j11=f(x,y)+x(−r1K+cmy(α+am+2aαx)(m+αx+amx+aαx2+bmy)2),j12=−cmx(1+ax)(m+αx)(m+αx+amx+αax2+bmy)2,

j21=y(αm(r3−r2)(m+αx)2)−y2(x+K1)2(2αm(a2αx−a1m)(x+K1)(m+αx)3−a1m2+a2α2x2(m+αx)2),

j22=g(x,y)−yx+K1(a1m2+a2α2x2(m+αx)2).

The Jacobian matrix of system (Equation 13) at trivial equilibrium (0,0) is given by:

(23) J|(0,0)=[r100r2]

The eigenvalues of the Jacobian matrix (Equation 23) are r₁ > 0 and r₂ > 0, which makes the equilibrium point P₀(0, 0) an unstable node. Thus, in the neighborhood of this equilibrium, there are unstable manifolds along both the x− and y− directions. This indicates that for any small perturbation in the densities of prey x(t) or predator y(t), the trajectories of system (Equation 13) will move away from the equilibrium P₀. In other words, this implies that the extinction of both species simultaneously is unstable, and the system may tend to evolve toward persistence of at least one species or both species, depending on the parametric conditions.

The Jacobian matrix of system (Equation 13) at predator free equilibrium point (x¯,0)=(K,0) is given by:

(24) J|(x¯,0)=[−r1−cmKm+(α+am+aαK)K0r2m+r3αKm+αK].

The eigenvalues of the above Jacobian matrix from Equation 24 are −r1<0 and r2m+r3αKm+αK>0, which makes the equilibrium P1(x¯,0) a saddle point having a stable manifold in the x− direction and an unstable manifold in the y− direction.

Therefore, for any small perturbation in prey density x(t), the trajectories of system (Equation 13) will return toward equilibrium P1(x¯,0). However, for a small positive increment in predator density y(t), the trajectories of system (Equation 13) will move into the first quadrant of the xy− plane away from the predator-free equilibrium state. Biologically, this indicates that, in system (Equation 13), the extinction of the predator species is unstable in the presence of abundant prey.

The Jacobian matrix of system (Equation 13) at prey free equilibrium at (0,y¯)=(0,r2K1a1) is:

(25) J|(0,y¯)=[r1−cr2K1a1+br2K10mr22+(r2+r3)αr2K1ma1−r2].

From Equation 25, The prey free equilibrium point P2(0,y¯) is asymptotically stable if:

(26) r1<cr2K1a1+br2K1,

(27) or K1>r1a1r2(c−br1), provided r1<cb.

Therefore, near a prey-free equilibrium, the prey population will go extinct, and the predator population will survive with some alternate food. These dynamics occur if the growth rate of prey is sufficiently low according to Equation 26. Equivalently, if the availability of alternative food (K₁) exceeds the critical threshold and the intrinsic prey growth rate (r₁) remains sufficiently low according to Equation 27, the predator population will persist near P₂ and suppress any prey growth, leading to prey extinction in system (Equation 13).

The Jacobian matrix of system (Equation 13) at existing interior equilibrium point P3(x,y) is given by:

(28) J|(x,y)=[j11j12j21j22]

where:

j11=−r1xK+cmxy(α+am+2aαx)(m+αx+amx+aαx2+bmy)2,

j12=−cmx(1+ax)(m+αx)(m+αx+amx+αax2+bmy*)2,

j21=y(αm(r3−r2)(m+αx)2)−y2(x+K1)2(2αm(a2αx−a1m)(x+K1)(m+αx)3−a1m2+a2α2x2(m+αx*)2),

j22=−yx+K1(a1m2+a2α2x2(m+αx)2).

The interior equilibrium point is locally asymptotically stable if tr(J (x^∗, y^∗)) < 0 and det(J(x^∗, y^∗)) >:

0, and these are given by:

tr(J)=−r1xK+cmxy(α+am+2aαx)(m+αx+amx+aαx2+bmy)2−yx+K1(a1m2+a2α2x2(m+αx)2),

Det(J)=(−r1xK+cmxy(α+am+2aαx)(m+αx+amx+aαx2+bmy)2)×(−yx+K1(a1m2+a2α2x2(m+αx)2))+(cmxy(1+ax)(m+αx+amx+aαx2+bmy)2)(αm(r3−r2)m+αx−y(x+K1)2(2αm(a2αx−a1m)(x+K1)(m+αx)2−a1m2+a2α2x2m+αx*)).

The analytical expression for det(J) > 0 and tr(J) < 0 is too complex to provide here. Thus, a sufficient condition is discussed, which gives det(J) > 0 and tr(J) < 0 and consequently ensures the stability of P₃(x^∗, y^∗).

From Equation 28, if j₁₁ < 0 and j₂₁ > 0, then tr(J) < 0 and det(J) > 0, respectively. Hence, the sufficient conditions for P₃(x^∗, y^∗) to be locally asymptotically stable (LAS) are as follows:

• r1>cKL(α+am+2aαK)m
• r3>r2+2αa2KL(m+αK)K1m3

The first condition gives a sufficient intrinsic growth rate for the prey population. Subsequently, the second condition specifies how much the predator’s growth in the breeding colony should exceed its growth in the foraging zone. With these sufficient growth rates, both populations get stabilized in the ecosystem (Equation 13).

Example 1. With the help of the parameter set provided in dataset (Equation 29), the stability of the interior equilibrium P₃ is analyzed, and the interior equilibrium point is calculated as P₃(x^∗, y^∗) = (6.9439,2.8352). The numerical simulation in Figures 4, 5 illustrates the same.

The parametric values in the dataset (Equation 29) are illustrative and satisfy all the biologically valid and necessary conditions for the existence of interior equilibrium P₃(x^∗, y^∗). Evaluating the Jacobian matrix (Equation 28) at P₃ yields:

(29) r1=1.5, r2=0.1, r3=0.2, K=25, a=0.04, b=0.1, c=0.8, m=0.5, α=0.03, K1=10, a1=1.5, a2=0.3.

J|(P3)=[−0.1898−1.22450.0940−0.2350],

and tr(J) and det(J) are given as follows:

tr(J)=−0.4248<0,

det(J)=0.1597>0.

Since tr(J) < 0 and det(J) > 0, the interior equilibrium P₃ is a locally asymptotically stable equilibrium.

4 Sensitivity analysis

To analyze the influence of parameters on the dynamics of system (Equation 13), a global sensitivity analysis is carried out by using two well-established statistical techniques, called partial rank correlation coefficient (PRCCs) and Latin hypercube sampling (LHS). This approach has been widely used as a reliable and efficient sampling-based method for global sensitivity analysis in ecological modeling (Biswas et al., 2020; Mondal et al., 2021; Chatterjee and Meng, 2024; Chatterjee et al., 2024). Compared with Monte Carlo methods, LHS offers a highly efficient sampling framework that provides an unbiased coverage of multidimensional parameter space with relatively few computational iterations. On the other hand, PRCC is explicitly chosen because it is highly robust for evaluating nonlinear but monotonic relationships, which frequently characterize continuous predator–prey dynamical systems. As a result, PRCC isolates the independent linear effect of each parameter on the population density through rank-transformed correlation while excluding the nonlinear confounding effects of other variables (Biswas et al., 2020; Mondal et al., 2021; Chatterjee and Meng, 2024; Chatterjee et al., 2024). The values of PRCCs lie in the interval [−1, 1]. The sign represents the type of correlation (positive or negative), while the value indicates the strength of the correlation. The parameters in Table 1 are considered input parameters, while the prey population (x(t)) and predator population (y(t)) are considered output variables. This study considers that a uniform distribution is assigned to each parameter, and the sampling range for each parameter is taken as ±25% of their respective nominal value listed in Equation 29 (Mondal et al., 2021). This range ±25% provides a substantial biological variance while ensuring that all parameters remain strictly positive. To ensure the statistical robustness and convergence of the LHS matrix, a large sample size of N = 1,000 runs is generated. The PRCC was calculated using the population density output evaluated after a sufficiently long integration time (t = 1,000) to ensure that transient dynamics had settled and the system had reached its final state. The significance of the PRCC values was determined using Spearman partial correlation. Uncertainty and statistical significance were reported using a strict p-value threshold, where ∗∗ denotes p < 0.01.

The PRCC values corresponding to the steady-state prey population density (x) are illustrated in Figure 6a. The result indicates that prey population dynamics are strongly and positively correlated with parameters like r₁, K, a, b, m, and a₂. This indicates that for an increment in these parameters, the prey density increases significantly. On the other hand, parameters like r₂, r₃, c, α, and K₁ share a strong negative correlation with the prey dynamics. This indicates that a small increment in these parameter values leads to a decline in prey population density. Figure 6b illustrates PRCC values corresponding to the steady-state of predator population density (y). The parameters r₁, K, a, b, α, and a₁ are strongly and positively correlated with predator population density, whereas the parameters c and m have a strong negative correlation with predator population density.

From this steady-state sensitivity analysis of system (Equation 13), the study establishes a robust criterion for selecting the most influential parameters for the subsequent local bifurcation analysis. First, evaluating the symmetric magnitude of the PRCC indices reveals that the intrinsic growth rate (r₁) of prey and the rate of predation (c) are the most dominant drivers of the entire ecosystem. They are the only parameters that exert a highly significant, high-magnitude influence (|PRCC| > 0.50) across both species simultaneously. Other parameters exert asymmetrical influences—for example, a₁ strongly influences the prey population but weakly influences the predator population. Furthermore, to investigate qualitative shifts caused by the non-linear functional response proposed in this model, the study isolates the aggregation parameters m and α. While their PRCC magnitudes are moderate, the analysis confirms that they remain highly statistically significant (p < 0.01) drivers of both populations. Therefore, since local bifurcation analysis investigates qualitative shifts in global stability, this study chooses r₁, c, m, and α as they represent a good analytical balance of the system’s most dominant global drivers r₁ and c, along with its core behavioral mechanism m and α.

5 Bifurcation analysis

In this section, various types of local (transcritical, saddle-node, and Hopf) and global (Bogdanov–Takens, Bautin, and cusp) bifurcations are analyzed for system (Equation 13).

To highlight the mathematical and ecological advantages of the proposed CPF framework, it is crucial to compare these bifurcation structures with previously published migration-based predator–prey models. In earlier spatial models utilizing aggregation methods, patch selection and density-dependent migration generally act as stabilizing forces or introduce standard codimension-one limit cycles via Hopf bifurcations (Chiorino et al., 1999). While more recent models incorporating fast, two-way density-dependent migrations for both species have successfully reported codimension-two Bautin bifurcations leading to bistability (El Abdllaoui et al., 2007), these frameworks rely on the assumption that populations can freely adjust their spatial distributions to avoid predation or seek resources (El Abdllaoui et al., 2007). In contrast, the rigid CPF constraints in system (Equation 13) impose a strict spatial anchor on the predator’s breeding colony. A distinct mathematical advantage of the CPF framework is demonstrating that even without bidirectional spatial escape, the foraging commuting parameters (m and α) alone are sufficient to drive the system into highly degenerate codimension-two structures, including Bogdanov–Takens, generalized Hopf, and, specifically, cusp bifurcations. Ecologically, this captures complex phenomena such as bistability traps, proving that recovery from a regime shift in a central place foraging ecosystem requires much more than simply restoring pre-collapse conditions.

5.1 Transcritical bifurcation

A transcritical bifurcation arises when two equilibrium points collide and exchange their stability as a system parameter is varied (Perko, 2013). In the system (13), this phenomenon appears around the prey-free equilibrium P2(0,y¯) when the parameter r1 is varied.

This occurs when one of the eigenvalues of the Jacobian matrix becomes zero. Since the eigenvalues of the Jacobian matrix (25) evaluated at the prey-free equilibrium P2(0,y¯) are:

λ1=r1−cr2K1a1+br2K1, λ2=−r2,

the eigenvalue λ₁ vanishes when:

r1=r1tc=cr2K1a1+br2K1

Theorem 4. The system (Equation 13) undergoes a transcritical bifurcation around the prey free equilibrium P2(0,y¯) at,

(30) r1tc=cr2K1a1+br2K1

provided that the non-degeneracy conditions of Sotomayor’s theorem are satisfied, i.e.,

Δ1=WTFr1(P2,r1tc)=0,

Δ2=WT[DFr1(P2,r1tc)V]≠0,

Δ3=WT[D2F(P2,r1tc)(V,V)]≠0.

where F = [xf,yg]^T, and V and W denote the right and left eigenvectors corresponding to the zero eigenvalue.

Proof. The eigenvectors of J(0,y¯) and (J(0,y¯))T corresponding to zero eigenvalues are given by:

V=(a2mK1(αr3+r2m+r2α),1)T and W = (1,0)^T respectively.

For r1tc=cr2K1a1+br2K1, from Equation 30 with the help of Sotomayor’s theorem (Perko, 2013), the values of Δ₁, Δ₂, and Δ₃ are given as follows:

Δ1=WTFr1(P2,r1tc)=0,whereF=[xf,yg]T.

Δ2=WT[DFr1(P2,r1tc)V]=a2mK1(αr3+r2m+r2α)≠0,

Δ3=WT[D2F(P2,r1tc)(V,V)]=2ca2m(r2a2(Kaa1m+Kαa1−ma1−mbr2K1)−Kma12(αr3+r2m+r2α))KK1(a1+br2K1)2(αr3+r2m+r2α)2≠0.

Due to algebraic complexity, the proof of Δ₃ ≠;= 0 is demonstrated through the following example:

The above analysis reveals that the threshold r1=r1tc represents the minimum intrinsic growth rate required for the persistence of the prey population. As observed in Figure 7, if environmental degradation (e.g., climate change, resource depletion, etc.) causes the decline of the intrinsic growth rate of prey below r1tc=0.5, the seabird fish ecosystem will face an ecological shift. As the prey population lacks a sufficient growth rate, it cannot sustain the predation pressure and eventually goes extinct from the ecosystem. Furthermore, this leads the predator population (e.g., seabirds) to rely entirely on alternative food resources (K₁) for survival.

Here BP indicates branching point or transcritical bifurcation point. It can be clearly observed that before the transcritical line r1tc=0.5, the system has only a stable prey-free equilibrium P2(0,y¯), denoted by the green line. This changes to a blue line after r1=0.5, indicating its unstable nature due to a transcritical bifurcation. Hence, in this region, the stability is exchanged with the interior equilibrium P3(x,y).

5.1.1 Numerical example

Example 2. For the dataset from Equation 29, the Jacobian matrix J(P3,a1tc) has zero eigenvalues at r1tc=0.5. The Jacobian matrix is given by,

(31) J|(P2,r1tc)=[000.0187−0.1]

From Equation 31, The eigenvector corresponding to zero eigenvalues (i.e., a₁₁ = 0) are V = (0.03165, 1)^T and W = (1, 0)^T. Now using Sotomayor’s theorem (Perko, 2013), Δ₁ = 0, Δ₂ ≈ 0.2542 ≠ 0, and Δ₃ ≈−0.3540 ≠ 0.

5.2 Hopf bifurcation

The characteristic equation of the Jacobian matrix (Equation 28) is given by:

(32) λ2−tr(J(x,y))λ+det (J(x,y))=0

For some critical value of c, obtained from tr(J(x^∗, y^∗)) = 0, given by:

c=cH=(m+αx+amx+aαx2+bmy)2mxy(α+am+2aαx)×(r1xK+yx+K1(a1m2+a2α2x2(m+αx)2)),

the above equation (Equation 32) will have purely imaginary root tr(J(x,y))=0 and det(J(x,y))>0. Then, the interior equilibrium P3(x,y) of system (Equation 13) changes its behavior from spiral sink to spiral source or vice-versa through Hopf bifurcation.

Theorem 5. The necessary and sufficient condition for system (Equation 13) to exhibit Hopf bifurcation around interior equilibrium point P₃(x^∗, y^∗) is that there exist c = c^Hsuch that:

• tr(J(x,y))|cH=0, and det (J(x,y))|cH>0,
• (dRe(λ)dc)cH≠0.

Proof. With the definition of critical value c=cH at interior equilibrium point P3(x,y), it can be stated that tr(J(x,y))|cH=0 and det(J(x,y))|cH>0. Furthermore, when c∈(cH−ϵ,cH+ϵ),ϵ>0, in the ϵ− neighborhood of cH, the roots of Equation 32 are given as:

λ1,2=M±iN

where M=Re(λ)=tr(J(x,y))2, N=4det (J(x,y))−(tr(J(x,y)))22.

Next, for the transversality condition:

Re(λ)=12tr(J(x,y))

Differentiating both sides with respect to c at c = c^H:

(dRe(λ)dc)|c=cH=12(dtr(J(x,y))dc)|c=cH

⇒(dRe(λ)dc)|c=cH=12(∂tr(J)∂xdxdc+∂tr(J)∂ydydc+∂tr(J)∂c)|c=cH

⇒(dRe(λ)dc)|c=cH≠0

Due to the presence of highly non-linear terms, the proof of condition (dRe(λ)dc)|c=cH≠0 is observed numerically.

5.2.1 Direction of limit cycle

The direction and stability properties of the periodic trajectories originating from P₃(x^∗, y^∗) through Hopf bifurcation can be calculated with the first Lyapunov coefficient l₁. Using the transformation z₁ = x − x^∗ and z₂ = y − y^∗, system (Equation 13) becomes (Manoharan et al., 2025a):

dz1dt=r1(z1+x)(1−z1+xK)−c(z1+x)(z2+y)1+(B(z1+x))(z1+x)+b(z2+y*)

dz2dt=(r2m+r2α(z1+x)m+α(z1+x))(z2+y)−(z2+y)2z1+x+K1(a1m2+a2α2(z1+x)2(m+α(z1+x*))2)

Expanding Taylor’s series of the above system at (z₁,z₂) = (0, 0) up to terms of order 3 produces the following system Equations 33 and 34:

(33) dz1dt=p10z1+p01z2+p20z12+p11z1z2+p02z22+p30z13+p21z12z2+p12z1z22+p03z23+O(|z|4)

(34) dz2dt=q10z1+q01z2+q20z12+q11z1z2+q02z22+q30z13+q21z12z2+q12z1z22+q03z23+O(|z|4)

where the coefficients p_ij and q_ij are given in Supplementary Appendix A.

By neglecting the higher-order terms of order 4 above from system (Equation 34), the first Lyapunov coefficient is given by:

l1=−3π2p01Δ3/2{(p10q10(p112+p11q02+p02q11)+p10p01(p112+p20q11+p11q02)+q102(p11p02+2p02q02) −2p10q10(q022−p20p02)−2p10p01(p202−q20q02)−p012(2p20q20+q11q20)+(p01q10−2p102)(q11q02−p11p20))−(p102+p01q10)(3(q10q03−p01p30)+2p10(p21+q12)+(q10p12−p01q21))}

where:

Δ=(det(J(x,y)))cH

Hence, if l₁ < 0, then a stable limit cycle emerges from the equilibrium point, indicating a supercritical Hopf bifurcation. Similarly, if l₁ > 0, then the equilibrium point leads to an unstable limit cycle which is again overlapped with a stable limit cycle, indicating a sub-critical Hopf bifurcation.

Example 3. The Jacobian matrix J(P₃, c^H) of P₃(4.937133, 2.020082) have purely imaginary eigenvalues at c^H = 1.045345, where the other parameters are from Equation 29. The trace of the Jacobian matrix is given by:

tr(J(P3))cH=−r1xK+cHmxy(α+am+2aαx)(m+αx+amx+aαx2+bmy)2−yx+K1(a1m2+a2α2x2(m+αx)2),

tr(J(P3))cH≈−2.393841×10−7≈0.

det(J(x,y))=(−r1xK+cmxy(α+am+2aαx)(m+αx+amx+aαx2+bmy)2)×(−yx+K1(a1m2+a2α2x2(m+αx)2))+(cmx(1+ax)(m+αx)(m+αx+amx+aαx2+bmy)2)(y(αm(r3−r2)(m+αx)2)−y2(x+K1)2(2αm(a2αx−a1m)(x+K1)(m+αx)3−a1m2+a2α2x2(m+αx)2))

and the transversality condition is given by:

(d(Re(λ))dc)cH=12dtr(J(x,y))dc|cH=12(∂tr(J)∂xdxdc+∂tr(J)∂ydydc+∂tr(J)alc)|c=cH.

The calculations above and in Figure 8 are performed to illustrate the continuation of codimension-one parametric bifurcation. While Figure 9 illustrates the time series and phase portrait of system (Equation 13) at the Hopf bifurcation point P₃(4.937133, 2.020082), Figure 8 shows the emergence of stable limit cycles (green curves) from the Hopf point P₃(4.937133, 2.020082).

5.3 Stability switching through Hopf bifurcation

This subsection investigates the stability switching at the interior equilibrium induced by two successive Hopf bifurcations under the variance of a single-system parameter. The system (Equation 13) experiences a stability switch at the interior equilibrium through two distinct Hopf bifurcations as the capture rate (c) of the predator varies. In this phenomenon, the interior equilibrium of the system first loses its stability through one Hopf bifurcation and again regains stability via another Hopf bifurcation. The following numerical example in (5.3.1) demonstrates this phenomenon.

5.3.1 Numerical example

Example 4. To illustrate the existence of two Hopf bifurcations numerically, the parameter set from Equation 29 is considered with threshold parameter c^H:

The first Hopf point occurs at c=cH1=1.045345, where P₃ = (4.937133, 2.020082), and the Jacobian matrix admits a pair of purely imaginary eigenvalues with vanishing real part. The transversality condition ddcRe (λ)|c=cH1≠0 is satisfied, which confirms the Hopf bifurcation.

Similarly, a second Hopf point is detected at c=cH2=1.627138 for which the interior equilibrium becomes P₃ = (2.062776, 1.123972). Between these two Hopf points, the interior equilibrium remains unstable and oscillates periodically (Volterra, 1926), while stability returns after cH2. This phenomenon is illustrated in Figures 10, 11.

Figure 12 depicts the variation of the real part of the eigenvalues with respect to c, showing two crossings of the imaginary axis. The one-parameter bifurcation diagram in Figure 10 confirms the analytical predictions.

The one-parameter bifurcation diagram with respect to the capture rate of predator (c) (Figure 11) reveals a complex, non-monotonic behavioral response within the seabird–fish ecosystem, culminating in a transcritical extinction threshold. The dynamics can be categorized into four distinct ecological phases:

• Sustainable foraging (c<cH1): At a lower capture rate, the predator–prey interaction is stable. Seabirds catch the fish at a sustainable rate, allowing the prey population to revive itself through logistic growth.
• Predator-driven periodic expansion and collapse (cH1<c<cH2): The coexistence equilibrium loses stability as the capture rate of the predator crosses the first Hopf threshold (cH1≈1.045), and the system oscillates in a stable periodic oscillation. Biologically, the seabirds become too efficient and rapidly over-exploit the prey. Consequently, food demands of the breeding colony increase, and the predator population eventually suffers a demographic collapse due to starvation. This further helps relieve predation pressure and allows the prey to revive, creating perpetual periodic expansion and collapse of population dynamics.
• Depleted stability (cH2<c<cBP): At an extremely high capture rate, the system undergoes another Hopf bifurcation (cH2≈1.627). The amplitude of oscillations dampens, and the system returns to a stable, but severely depleted, equilibrium. Here the primary prey is suppressed to a low density. The system avoids immediate collapse only because seabirds rely heavily on alternative food resources (K₁) to supplement their diet while continuously hunting the primary prey.
• Transcritical extinction (c = c_BP): If the predator capture rate continues to increase, the system reaches a critical branch point (c_BP ≈ 1.944). At this threshold, the coexistence equilibrium collides with the prey-free equilibrium and exchanges stability. Biologically, predation pressure becomes absolute; the primary fish population (x) can no longer be sustained and is at risk of extinction. Beyond this point, seabirds must survive solely on alternative food resources (K₁).

5.4 Role of migration parameters (m and α)

5.4.1 Effect of migration parameter m

The effect of the migration parameter m, which governs the movement of the predator from the breeding colony to the foraging zone, is studied with numerical simulation as shown in Figure 13. Table 5 summarizes as well the detected bifurcation points and limit point of cycles (LPC).

Table 5. Codimension-one bifurcation points of system (Equation 13) with respect to the migration parameter m.

From Figure 13a, system (Equation 13) exhibits a rich dynamical behavior with respect to parameter m, where it can be observed that there exist two Hopf bifurcations (H₁ and H₂) and two saddle-node bifurcations (LP₁ and LP₂). From Figure 13c, in R₁, there exists only one stable interior equilibrium with high population densities. As the value of m crosses LP₂ and enters R₂, the system exhibits three interior equilibria (stable, saddle, and unstable). Furthermore, as the parametric value of m increases, system (Equation 13) enters the subcritical region R₃, where there exist three similar equilibria along with stable and unstable limit cycles. In region R₄, the unstable limit cycles shrink, and only a stable limit cycle exists along with three unstable equilibria. At the threshold LP₁, the unstable equilibrium and saddle equilibrium merge and disappear as the parametric value of m gradually increases. Thus, in R₅, there exists only one unstable interior equilibrium along with a stable limit cycle. This stable limit cycle shrinks near H₂, and the low-population-density unstable equilibrium gains stability through super-critical Hopf bifurcation (see Figures 13a, b).

From the analysis, it is observed that the subcritical Hopf bifurcation gives a warning for the ecosystem, predicting large-amplitude oscillations which possess the risk of extinction of at least one species (e.g., the stable limit cycle in 13(b) becomes extremely close to the x and y− axes)—that is, whenever unstable oscillations are observed due to a change in the outward migration rate of predators, the system becomes fragile and can shift regime from R₃ to R₄ or R₅ for a minor increase in values of m (e.g., due to an increase in foraging efforts). Hence, the coexistence equilibrium loses stability, and the population oscillates, which can cause risk in a real-life scenario.

5.4.2 Effect of migration parameter α

A similar dynamical behavior is observed when the migration parameter α is varied, which controls predator movement from the foraging ground to the breeding colony. The corresponding bifurcation points are listed in Table 6 and the dynamics can be observed in Figures 14a, b.

Table 6. Codimension-one bifurcation points of system (Equation 13) with respect to the migration parameter α.

As α increases, the prey population temporarily declines until the threshold H₂. However, as soon as the predators start to leave the foraging ground, the predation pressure drops, and the prey gets a chance to revive. Both populations oscillate in a stable periodic oscillation. A further increase in the value of α gives a saddle-node bifurcation at LP₁, where saddle and unstable equilibria emerge. The unstable interior equilibrium undergoes a sub-critical Hopf bifurcation at H₁, followed by a saddle-node bifurcation at LP₁, where one saddle and one unstable interior equilibria merge. Whenever the parametric value of α crosses H₁, the unstable equilibrium becomes stable with high population density. Interestingly, whenever the system exhibits an unstable limit cycle and the α value increases, the interior equilibrium becomes stable.

Interestingly, from Tables 5, 6, the Hopf bifurcation points H₁ and H₂ of both parameters m and α are occurring at the same interior equilibria (11.522553, 15.871092) and (6.280760, 3.036407) respectively. Furthermore, the saddle–node bifurcation points L₁ and L₂ are also coinciding for both m and α at (9.197642, 14.222238) and (2.419143, 7.5558), respectively.

5.5 Saddle-node bifurcation

If at the interior equilibrium det(J(x^∗, y^∗)) = 0 and tr(J(x^∗, y^∗)) ≠ 0, then one of the eigenvalues of J(x^∗, y^∗) is zero for m = m^SN, which further creates a possibility for a saddle-node bifurcation at P₃(x^∗, y^∗).

Theorem 6. The system (Equation 13) at the interior equilibrium P₃(x^∗, y^∗) exhibits a saddle-node bifurcation for a critical value m = m^SN, where m^SN satisfies:

(35) det(J(P3))|m=mSN=0, tr(J(P3))|m=mSN≠0.

Furthermore, the non-degeneracy of the bifurcation is ensured through the transversality conditions of the Sotomayor theorem (Perko, 2013).

Proof. The existence of a saddle-node bifurcation in system (Equation 13) is established by using Sotomayor’s theorem. According to the theorem, the Jacobian J(x^∗, y^∗) possesses a simple zero eigenvalue, while all other eigenvalues have a non-zero real part at m = m^SN. The eigenvectors of J(x^∗, y^∗) and its transpose (J(x^∗, y^∗))^T associated with the zero eigenvalue are given by:

V=[Kcmx(1+ax)(m+αx)−r1x(m+αx+amx+aαx2+bmy)2+Kcmxy(α+am+2aαx*)1] and W=[−1j11j21]

Now compute:

Δ1=WTFm(P3;mSN), where Fm(P3;mSN)=∂F∂m|(P3;mSN) and F=[xf,yg]T,=−cαx2y(1+ax)(mSN+αx+amSNx+aαx2+bmSNy)2−j11j21×y(αx(r2−r3)(mSN+αx)2−2αxy(a1mSN−a2αx)(x+K1)(mSN+αx)3)≠0,

Δ2*=WT[D2F(P3,mSN)(V,V)]=(1)·(∂2(xf(x,y))∂x2v12+2∂2(xf(x,y))∂x∂yv1v2+∂2(xf(x,y))∂y2v22)|(P3,mSN)+{(−j11j21)·(∂2(yg(x,y))∂x2v12+2∂2(yg(x,y))∂x∂yv1v2+∂2(yg(x,y))∂y2v22)}|(P3,mSN)≠0.

where:

∂2∂y2(xf(x,y))=2m2bcx(m+(α+am+aαx)x*+bmy)2−2m3b2cxym+(α+am+aαx)x+bmy

∂2∂x2(yg(x,y))=−2α2my(r3−r2)(m+αx)3−2y2α2m(a2m+3a1m−2a2αx)(x+K1)(m+αx)4+4y2αm(a2αx−a1m)(x+K1)2(m+αx)3−2y2(a1m2+a2α2x2)(m+αx)2(x+K1)3

∂2∂x∂y(yg(x,y))=αm(r3−r2)(m+αx)2−2y(2α2a2x(m+αx)2(x+K1)−(a1m2+a2α2x2)(m+3αx+2αK1)(m+αx)3(x+K1)2)

∂2∂y2(yg(x,y))=−2(a1m2+a2α2x2(x+K1)(m+αx)2)

Due to the presence of highly non-linear algebraic terms, the non-zero conditions of Δ1* and Δ2* have been demonstrated using the following example.

5.5.1 Numerical example

Example 5. Using the parameter values listed in dataset (Equation 29), system (Equation 13) admits an interior equilibrium P₃(2.419179, 7.555952) at the critical migration parameter m^SN = 0.030507. The Jacobian matrix evaluated at (P₃, m^SN) is given by:

(36) J|(P3,mSN)=[0.7465−0.36030.3530−0.1704]

From Equation 36, the determinant is det(J) ≈ 1.3768 × 10⁻⁷ ≈ 0 and the trace is tr (J) ≈ 0.5761 ≠ 0, which confirms the condition of eigenvalues (from Equation 35). The transversality conditions are Δ1≈−21.020615≠0 and Delta2≈0.068965≠0. Since both Δ1* and Δ2* are non-zero, all of the conditions of Sotomayor’s theorem are satisfied. Hence, the existence of saddle-node bifurcation is confirmed at m = m^SN. Figure 13a illustrates the saddle-node bifurcation, as the existence of LP₁ and LP₂ points can be observed. Similarly, Figure 15 shows these two saddle-node bifurcations as the interior equilibrium P_3a splits into P_3a1 and P_3a2. Furthermore, the two equilibria P_3a2 and P_3b collide and merge into P_3c.

Ecologically, the saddle-node bifurcation represents the sudden creation or annihilation of feasible coexistence states between the seabirds and their prey. Since the outward migration m is driven by intrinsic colony demands (e.g., chick provisioning, adult metabolic needs, etc.), the parameter m serves as the foraging effort of the breeding colony. The saddle-node bifurcation demonstrates that this foraging effort must remain within a strict, narrow bandwidth. If m crosses the critical threshold m1SN≈0.030507, an interior equilibrium is created, and if m crosses the threshold m2SN≈0.040982264, then two interior equilibria collide and vanish.

5.6 Generalized Hopf (Bautin) bifurcation

In this subsection, the codimension-two bifurcation known as Bautin or generalized Hopf (GH) bifurcation is investigated. This bifurcation appears when the Jacobian (Equation 28) admits purely imaginary eigenvalues at the non-saddle interior equilibrium P₃(x^∗, y^∗) and the first Lyapunov coefficient vanishes.

Additionally, considering the bifurcation parameters as r₁ and m, the critical point for the generalized Hopf (GH) bifurcation will lie on the Hopf bifurcation curve in the r₁m − parametric plane. The following theorem will show that system (Equation 13) at P₃(x^∗, y^∗) exhibits a GH bifurcation for the threshold parameter ((r1GH,mGH).

Theorem 7. The system (Equation 13) undergoes codimension-two GH bifurcation at parametric threshold values (r1GH,mGH), when the following conditions hold:

• det(J(P₃; r1=r1GH,m = m^GH)) > 0
• tr(J(P₃; r1=r1GH,m = m^GH)) = 0
• l₁(P₃; r1=r1GH,m = m^GH) = 0
• l₂(P₃; r1=r1GH,m = m^GH) ≠ 0

where I₁ and l₂ denote the first and second Lyapunov coefficients.

Proof. The proof of this has been extended in Supplementary Appendix C.

5.6.1 Numerical example

Example 6. The system (Equation 13) undergoes GH bifurcation at interior equilibrium P₃(11.125944, 15.318847) for threshold parameter r1GH=2.3475494 and mGH=0.075013686, i.e.,

• det(J(P₃; r1=r1GH,m = m^GH)) ≈ 9.486185 × 10⁻² > 0.
• tr(J(P₃; r1=r1GH,m = m^GH)) ≈ 9.097402 × 10⁻⁹ ≈ 0.
• l₁(P₃; r1=r1GH,m = m^GH) ≈ 0.
• l₂(P₃; r1=r1GH,m = m^GH) ≈−2.019507 × 10⁻⁴ ≠ 0

Hence, GH bifurcation takes place at P₃(x^∗, y^∗) for threshold parametric value (r1GH,mGH). To numerically observe the GH bifurcation, the study utilized the numerical continuation software MATCONT. The computations were performed using the following exact continuation settings and tolerances: VarTolerance = 10⁻⁶, FunTolerance = 10⁻⁶, TestTolerance = 10⁻⁵, InitStepsize = 0.01, MaxStepsize = 3, and MaxNumPoints = 300.

The points GH and BT in Figure 16a are generalized Hopf and Bogdanov–Takens bifurcation points, respectively. Figure 16b illustrates how the Hopf bifurcation curve (cyan curve) arises from two Hopf points H₁ and H₂, and the GH point occurs on it.

Ecological implications (ecosystem resilience and regime shifts): Biologically, the generalized Hopf (Bautin) bifurcation acts as a critical threshold for stability in ecosystem resilience. This bifurcation separates the parameter regimes of a predictable, stable population coexistence from highly fragile, oscillating states. As observed in Figure 17b, when the intrinsic growth rate of the prey (r₁) and the outward migration of the predator reach the critical threshold (r1GH,mGH), the first Lyapunov coefficient vanishes, and the coexistence equilibrium becomes a weak focus. Thus, the black trajectories appear not to reach the equilibrium P₃.

Furthermore, this bifurcation shows potential regime shifts in the parameter space. If the prey’s intrinsic growth rate (r₁) decreases below the critical threshold r1GH or if the outward migration of predators to the foraging zone (m) exceeds the threshold m^GH (e.g., due to depletion of resources in the foraging zone for prey, increase in food demand at breeding colony, etc.), the ecosystem loses its resilience. This forces a regime shift, where the system enters into high-amplitude oscillations (see Figures 17a, d).

Conversely, if ecological conditions improve, such that the intrinsic growth rate of prey (r₁) exceeds r1GH (overcoming predation pressure) or the outward migration of predators drops below m^GH, the ecosystem regains its resilience. The high-amplitude oscillations dampen entirely, and the predator–prey dynamics stabilize around the coexistence equilibrium P₃ (see Figures 17c, e).

5.7 Bogdanov–Takens bifurcation

This subsection investigates the appearance of the Bogdanov–Takens (BT) bifurcation of system (Equation 13) in the neighborhood of the interior equilibrium point P₃(x^∗, y^∗). It is a two-parametric bifurcation of an equilibrium point in an autonomous system. In a two-dimensional parametric plane, such a bifurcation arises at the intersection of a Hopf bifurcation curve and a saddle-node bifurcation curve.

Equivalently, a BT bifurcation may appear when the Jacobian matrix (Equation 28) of system (Equation 13), evaluated at the interior equilibrium P₃(x^∗, y^∗), possesses a zero eigenvalue with an algebraic multiplicity of two. The existence of this BT bifurcation can further be ensured with the transversality and nondegeneracy conditions.

Since it is observed that for parameter m both Hopf and saddle-node bifurcations are taking place, parameter m is chosen as one bifurcation parameter for BT bifurcation, and another parameter is taken as r₁ since BT bifurcation requires two parameters to take place. Assuming that the model undergoes BT bifurcation at P₃(x^∗, y^∗), the parametric thresholds are considered (r1,m)=(r1BT,mBT). From general bifurcation theory (Perko 2013), it is known that for interior equilibrium P₃(x^∗, y^∗) the Jacobian matrix (Equation 28) is similar to [0100] at (r1BT,mBT)) (Perko, 2013; Pal et al., 2019; Sen et al., 2020; Sen et al., 2022; Kumar et al., 2025; Vishwakarma and Mondal, 2026; Hu et al., 2025). The following theorem provides a condition for system (Equation 13) to undergo BT bifurcation.

Theorem 8. The system undergoes a codimension-two BT bifurcation at parametric threshold values (r1BT,mBT) when the following conditions hold:

• tr(J(P₃; r1=r1BT, m = m^BT)) = 0
• det(J(P₃; r1=r1BT, m = m^BT)) = 0
• and the transversality condition is given by:

∂2F1∂x2−∂F1∂x∂2F1∂x∂y∂F1∂y+∂2F2∂x∂y≠0, 12∂F1∂x∂2F1∂x2−(∂F1∂x)2∂2F1∂x∂y∂F1∂y+12∂F1∂y∂2F2∂y2−∂F1∂x∂2F2∂x∂y≠0,

where F₁ = xf(x, y) and F₂ = yg(x, y) from (Equation 13) (Sen et al., 2020).

5.7.1 Numerical example

Example 7. All of the values of parameters for the numerical example are the same as in (Equation 29), except threshold parameters r1BT and mBT. System (Equation 13) undergoes BT bifurcation at interior equilibrium P₃(10.580850, 14.717666) at these threshold parameters r1BT=0.769236 and mBT=0.019108, i.e.,

• tr(J(P₃; r1=r1BT,m = m^BT)) = 1.405364 × 10⁻⁶ ≈ 0.
• det(J(P₃; r1=r1BT,m = m^BT)) = −1.649798 × 10⁻⁷ ≈ 0.

The transversality condition is:

∂2F1∂x2−∂F1∂x∂2F1∂x∂y∂F1∂y+∂2F2∂x∂y≈0.199378≠0,

12∂F1∂x∂2F1∂x2−(∂F1∂x)2∂2F1∂x∂y∂F1∂y+12∂F1∂y∂2F2∂y2−∂F1∂x∂2F2∂x∂y≈0.018540≠0,

Therefore, BT bifurcation takes place at P₃(x^∗, y^∗) for threshold parametric value (r1BT,mBT).

An illustration of this phenomenon of BT bifurcation is shown in Figure 16. Figure 16a, b show how the Hopf bifurcation curve (cyan curve) intersects with the saddle-node bifurcation curve (magenta curve), forming a Bogdanov–Takens (BT) point at P₃(10.580850, 14.717666). The numerical continuation software MATCONT is used to capture the BT bifurcation in the r₁m– parametric plane as shown in Figure 16a. The computations were performed using the following exact continuation settings and tolerances: VarTolerance = 10⁻⁶, FunTolerance = 10⁻⁶, TestTolerance = 10⁻⁵, InitStepsize = 0.01, MaxStepsize = 3, and MaxNumPoints = 300.

Ecological implications (ecosystem resilience and regime shifts): Biologically, the Bogdanov–Takens (BT) bifurcation serves as a critical threshold of sensitivity for the seabird–fish ecosystem (Equation 13). It represents an extreme tipping point, where the intrinsic growth rate of the prey (r₁) and the predators’ outward foraging migration (m) reach a highly sensitive, unbalanced state (see Figure 18b). The resilience of the ecosystem is completely compromised at this critical threshold point.

Since this point organizes multiple bifurcation curves (see Figure 16), any slight environmental fluctuation results in a regime shift—for example, if due to environmental degradation, the prey growth rate drops (r1<r1BT) or the foraging pressure increases (m > m^BT), and the coexistence equilibrium P₃_b completely vanishes. The system (Equation 13) is pushed into a stable limit cycle with extreme, large-amplitude periodic oscillations (see Figures 18a, d), leaving the populations highly vulnerable to stochastic extinction.

Conversely, if conditions favor prey survival either through increased intrinsic growth (r1>r1BT) or a reduction in the predators’ outward migration (m<mBT), the ecosystem undergoes a stabilizing regime shift. With a decrease in the intensity of foraging pressure, the prey population revives. As a result, the degenerate equilibrium splits and gives birth to a new, highly stable coexistence state (P₃_b₂) alongside a saddle point (P₃_b₁) (see Figures 18c, e). Through this transition, the ecosystem’s resilience is effectively restored, allowing both species to stably persist without dangerous cyclical collapses.

5.8 Cusp bifurcation

This subsection will discuss the codimension-two cusp bifurcation. This bifurcation occurs at a critical point on a saddle-node bifurcation curve where the quadratic terms of the system’s normal form (Kuznetsov, 1998) on the center manifold vanish, leading to a higher-order degeneracy. Geometrically, this point marks the tip of a wedge-shaped region in the two-parametric plane, inside which the system has three coexisting equilibria, whereas the system possesses only one equilibrium outside this region. The Cusp point acts as an organizing center for two branches of the saddle-node bifurcation curve, which meet tangentially at this point (Kuznetsov, 1998).

Let us consider system (Equation 13) is given as x˙=F(x), where x=[xy], and F(x)=[xf(x,y)yg(x,y)], then the Taylor expansion of the system is given by:

F(x)≈Jx+12!B(x,x)+13!C(x,x,x)+…

Here J is the Jacobian (Equation 28) of system (Equation 13). The term B(u, v) is a symmetric bilinear form whose components are given by the second-order partial derivatives of F(x). Similarly, C(u, v, w) is a symmetric trilinear form involving the third-order partial derivatives of F(x). Using u = v = w = x:

B(x,x)=[∂2(xf)∂x2x2+2∂2(xf)∂x∂yxy+∂2(xf)∂y2y2∂2(yg)∂x2x2+2∂2(yg)∂x∂yxy+∂2(yg)∂y2y2]and

C(x,x,x)=[∂3(xf)∂x3x3+3∂3(xf)∂x2∂yx2y+3∂3(xf)∂x∂y2xy2+∂3(xf)∂y3y3∂3(yg)∂x3x3+3∂3(yg)∂x2∂yx2y+3∂3(yg)∂x∂y2xy2+∂3(yg)∂y3y3]

Theorem 9. The system (Equation 13) undergoes a codimension two cusp bifurcation at the interior equilibrium point P3(x,y) for a critical set of parameter values (r1CP,mCP), when the following conditions hold:

• Det(J(P3;r1=r1CP,m=mCP))=0.
• tr(J(P3;r1=r1CP,m=mCP))≠0.
• The quadratic normal form coefficient on the center manifold is zero at P3(x,y) for (r1CP,mCP), which is given by where, V=[−j12j111] and W=[−1j11j21] are the eigenvectors of Jacobian J (Equation 28) and J^T with respect to zero eigenvalue, and p and q are normalized eigenvectors of V &W such that 〈p, q〉 = 1.
• The cubic normal form coefficient is non-zero, given by where, h solves J(P3;r1=r1CP,m=mCP)h=−B(p,p)|(P3;r1=r1CP,m=mCP).

Proof. For detailed proof, one can follow (Kuznetsov, 1998; Lessard and Pugliese, 2024). The illustration of the above phenomena is performed through the following example.

Example 8. The numerical simulations are performed in Figure 19 to illustrate the phenomenon of cusp bifurcation. For Figure 19(a), all MATCONT settings, tolerance, and step size are the same as used in GH and BT bifurcations.

Again, it is observed that, for some parametric values of (r1CP,mCP) inside the cusp region of Figure 19a, there can be the existence of two stable interiors. This phenomenon is illustrated in Figure 20. Here, the red dots denote the stable equilibrium, whereas the cyan dots denotes saddle equilibrium. The green and black curves denote the hetero-clinic orbits of a saddle equilibrium. The figure clearly demonstrates that if the initial condition of the predator–prey system lies on the right side of the black curves, then the system will attain the stable interior equilibrium P₃_a with high population density. If the system starts before the black curves and outside the magenta limit cycle, then the system will attain the axial saddle equilibrium P₂, and if it starts inside the magenta limit cycle, the system will attain stable interior equilibrium P₃_c with low population density.

It can be observed that the trajectories of system (Equation 13) move more quickly toward the equilibrium P₃_a compared to P₃_c in their respective basin of attractors. In other words, the high population density equilibrium (P₃_a) gets stabilized more quickly as compared to the low population density equilibrium (P₃_c) (See Figure 21).

Ecological implications (ecosystem resilience and regime shifts): Biologically, the cusp bifurcation reveals the existence of a dangerous bistability trap within the ecosystem. It can be observed from Figures 20 and 21 that when parameters fall inside the cusp region (inside the bounded region of Figure 19a), the ecosystem loses its predictable, single-state resilience. It instead exhibits alternative stable states: a high-population-density coexistence (P₃_a) and a severely depleted, low-population-density survival state (P₃_c), separated by the saddle equilibrium (P₃_b), which acts as a tipping point.

Crucially, this bifurcation explains how history and initial conditions guide the future ecological state. If a severe environmental shock (e.g., sudden resource depletion or overpredation) forces the population densities across the saddle boundary (the black separatrix curves in Figure 20), a catastrophic regime shift occurs, trapping the ecosystem in the depleted state (P₃_c).

The resilience of the ecosystem is fundamentally compromised due to this bistability. This means that merely restoring the environmental parameters (such as prey’s intrinsic growth r₁ or outward foraging migration m of the predator) to their previously favorable levels is not significant enough to recover the populations. In order to reverse this regime shift and restore the ecosystem to the high-population-density state (P₃_a), conservation efforts must be introduced forcefully to push the population densities back across the separatrix threshold.

6 Conclusion

This study investigated the dynamics of a three-dimensional slow–fast predator–prey system modeled around the constraints of the central place foraging (CPF) framework. Considering the fact that central place foragers, like marine seabirds, commute between the foraging zone (e.g., sea) and their breeding zone (e.g., island with breeding colony) on a much faster timescale than demographic changes, the study incorporates the idea of slow–fast dynamics. Furthermore, using the method of aggregation, the fast commuting dynamics are reduced to their quasi-equilibrium distributions. This yielded a two-dimensional aggregated system on the slow demographic timescale, incorporating a Beddington–DeAngelis functional response and a modified Leslie–Gower growth formulation.

To explicitly distinguish the contributions of this CPF framework from standard modified Leslie–Gower formulations, it is crucial to recognize what drives the observed complex bifurcations. In generic, non-spatial predator–prey models, population collapses and regime shifts (such as those indicated by saddle-node or cusp bifurcations) are typically triggered by macro-environmental degradation, such as a drop in prey carrying capacity or increased direct predator mortality. However, the introduction of the present slow–fast CPF model reveals a fundamentally different mechanism: our numerical simulations suggest that even in a biologically robust environment, catastrophic regime shifts can be induced solely by micro-scale behavioral constraints—specifically, the commuting effort (m). If this spatial framework were not introduced, the extreme vulnerability of central place foragers to their own time-budget constraints would remain mathematically invisible. This model uniquely demonstrates that, within this theoretical framework, stable environmental conditions may not guarantee population survival if the physical and temporal costs of commuting cross critical bifurcation thresholds.

Under specific parametric conditions, the analytical results confirm the fundamental well-posedness of the aggregated model by establishing properties such as positivity, boundedness, and permanence. The local stability analysis of the trivial equilibrium reveals that whenever the initial population density is in close proximity to P₀, both species cannot go extinct simultaneously. Similarly, the analytical results of the predator-free equilibrium show that whenever the initial population density is in the neighborhood of P₁, predator extinction is unstable due to prey abundance. Furthermore, the analysis of the prey-free equilibrium (P₂) established the baseline survival constraints, revealing that alternative food resources (K₁) allow predators to sustain even in the absence of the primary prey population. In spite of the system’s rich non-linear dynamics and ecological regime shifts, specific behavioral and demographic parameters play a crucial role: the prey’s intrinsic growth rate (r₁), the predator’s capture rate (c), the outward migration rate (m), and the prey-dependent migration coefficient (α).

A one-parametric numerical bifurcation analysis with respect to the prey’s intrinsic growth rate (r₁) reveals a critical transcritical bifurcation, establishing a theoretical threshold for the sustainable prey population. If, due to any environmental degradation, the growth rate of prey declines below this threshold, the model indicates that the predator population completely suppresses and drives the prey to extinction. However, once r₁ crosses the threshold (r1tc), the prey successfully overcomes the predation pressure and establishes a stable coexistence.

On the other hand, numerical continuation reveals that the predator’s capture rate (c) drives a complex stability switch of the coexistence equilibrium through two successive Hopf bifurcations. As the foraging efficiency of the predator increases, a stable coexistence equilibrium with high population density loses stability and oscillates periodically, but with a further increase in foraging efficiency, it dampens this oscillation due to top-down predation pressure, and the system returns to a stable coexistence with lower population density.

Additionally, the one-parameter bifurcation analysis of m and α illustrates the profound theoretical impact of the predator’s behavioral commuting on the structural stability of the ecosystem. The outward migration rate (m) and the prey-dependent return rate (α) reveal a strict bandwidth for sustainable foraging. Fluctuations in these parameters trigger a complex sequence of saddle-node and Hopf bifurcations, forcing the system through regions of multiple stable states (e.g., stable coexistence along with stable periodic oscillation). Ecologically, this suggests that the predator’s spatial flexibility highly influences the ecosystem’s resilience; an imbalance between outward foraging effort and density-dependent return forces the system into highly fragile states. Within the model, an early warning signal is the emergence of stable limit cycles before saddle-node collapse, which indicates that the ecosystem is dangerously close to a catastrophic regime shift.

This study further analyzed two-parameter numerical continuation to identify critical codimension two bifurcations: the Bogdanov–Takens (BT), generalized Hopf (GH), and cusp bifurcation points. A double-zero eigenvalue characterizes the BT point, which serves as a fundamental organizing center from which saddle-node and Hopf bifurcation curves emerge. On the other hand, the GH point acts as a sensitive boundary, separating supercritical regimes of stable, predictable population cycling from subcritical domains of extreme fragility. In these subcritical regions, minor parameter drifts can cause high-amplitude oscillations that increase the extinction risk of at least one species. Furthermore, the detection of a cusp bifurcation allows the presence of alternative stable states. Thus, the cusp point theoretically implies that the long-term behavior of the ecosystem can be entirely dictated by initial conditions if the respective parameters are well within the cusp region.

The mathematical thresholds derived in this study offer theoretical insights—rather than empirically validated predictions—for real-world ecosystem management and conservation, particularly for central place foragers like marine seabirds. Drawing upon the ecological implications of our bifurcation analysis, we highlight two potential conservation strategies:

• Preventing deterministic extinction: The transcritical threshold (r1tc) determines a theoretical absolute minimum prey growth rate required before total ecosystem collapse. This conceptually emphasizes the need for strict, data-driven fishing quotas in the foraging zones of seabirds. The model suggests that the management of seabird colonies could be supported by maintaining a sustainable baseline prey population by ensuring that human harvesting does not impact the intrinsic growth of the species.
• Avoiding bistability traps and regime shifts: As revealed by the higher-order cusp and generalized Hopf bifurcations, simply restoring favorable environmental conditions after a population collapse is often insufficient due to a bistability trap. This mathematical hysteresis conceptually mirrors observed ecological transitions in marine environments, such as sudden seabird colony abandonments where predator populations fail to recover even after local fish stocks are artificially replenished. Therefore, our theoretical results suggest that conservation strategies must be proactive rather than reactive. Management efforts could benefit from monitoring the predators’ outward foraging effort (m) and ensuring sufficient local prey density to optimize return rates (α). Keeping these commuting parameters within the safe operational bandwidths identified in our numerical analysis could theoretically serve as an early warning system to prevent sudden, irreversible ecological regime shifts.

Overall, this aggregated slow–fast model analytically demonstrates that macro-scale demographic survival is closely related to micro-scale behavioral and foraging parameters. Importantly, these theoretical outcomes differ substantially from standard non-spatial frameworks by mathematically demonstrating that catastrophic ecosystem collapses do not strictly require massive macro-environmental degradation; they can theoretically be triggered solely by the behavioral time-budget constraints of the foragers. Our numerical simulations show that even minor environmental degradation may push this ecosystem toward critical threshold points due to the extreme sensitivity of central place foragers to changes in commuting cost, prey revival, and foraging efficiency. Using this mathematical framework, future ecological management and conservation efforts will be able to gain a theoretical understanding of the vulnerability of the spatially distributed predator–prey system.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2026.1896471/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Jiaxin Zhang, Huazhong University of Science and Technology, China

Reviewed by: Srilena Kundu, Helmholtz Institute for Functional Marine Biodiversity (HIFMB), Germany

Anand Singh, Birla Institute of Technology and Science, India