Abstract
Introduction:
This study investigates a three-dimensional fractional predator-prey model with adaptive refuge and predator cannibalism. The model combines logistic prey growth, Holling type-II predation, saturating predator self-regulation, and a dynamically evolving refuge variable interpreted as a global hunting-inhibition index.
Methods:
Memory effects are introduced through the Caputo derivative of order 0 < α ≤ 1. Positivity, invariance of the biologically feasible region, existence, uniqueness, and boundedness are established under transparent dissipativity assumptions. Boundary and coexistence equilibria are derived, and local stability is analyzed through the Jacobian matrix and Matignon's criterion. A fractional Hopf-type stability boundary is formulated in terms of eigenvalue arguments. Numerical simulations are performed using the Adams-Bashforth-Moulton predictor-corrector method.
Results:
The positive coexistence equilibrium is reduced to explicit feasibility conditions. For the Hopf-type parameter set, the critical fractional order is α* ≈ 0.8410, separating locally stable and unstable fractional-order regimes. The simulations show that fractional memory, refuge activation, and refuge decay substantially modify transient prey-predator dynamics.
Discussion:
The results show that adaptive refuge and memory can reshape boom-bust transients and shift the stability threshold of coexistence states. The model provides a compact framework for studying the combined effects of behavioral inhibition, cannibalistic self-regulation, and fractional memory in ecological systems.
1 Introduction
Predator–prey models remain one of the central mathematical frameworks for describing trophic interactions. The classical Lotka–Volterra equations introduced self-sustained oscillations through the coupling of prey growth and predator consumption [1, 2]. Subsequent ecological models introduced carrying capacity, saturating predation, harvesting, fear, refuge, disease, stage structure, and predator self-regulation in order to obtain mechanisms closer to observed population dynamics [3–5].
Two mechanisms are particularly relevant for the present work. The first is refuge or avoidance. Refuge may represent spatial shelter, behavioral evasion, reduction of predator activity, or fear-induced reduction of exposure. In predator–prey models, refuge can change extinction thresholds, delay predator growth, and alter the stability of coexistence states [4–6]. The second mechanism is predator cannibalism. Cannibalism is observed in many species and may function as a nonlinear self-limiting process when predator density becomes large. Mathematical studies show that cannibalism can stabilize coexistence, create additional equilibria, or induce bifurcations depending on its strength and saturation structure [7–9].
Fractional-order derivatives provide an additional modeling layer when the present rate of change depends not only on the current state but also on past ecological history. In biological populations, delayed behavioral responses, accumulated environmental stress, learning, territorial memory, and non-instantaneous recovery processes motivate nonlocal-in-time formulations. The Caputo derivative is especially convenient because it permits standard initial conditions and has become common in fractional ecological models [9–12].
The main mathematical difference between integer-order and fractional-order stability is captured by Matignon's criterion. For a commensurate fractional system of order 0 < α ≤ 1, an equilibrium is locally asymptotically stable when every eigenvalue λi of the Jacobian satisfies
Thus the fractional order α is not a passive numerical parameter. It changes the stability wedge and may create a stability transition even when all ecological parameters are fixed [13, 14]. These perspectives are complemented by Volterra-type Lyapunov methods for fractional systems [15], fractional Adams algorithms [16, 17], recent predator-prey bifurcation studies [18], global stability criteria for fractional differential equations [19], and discrete Rosenzweig-MacArthur predator-prey dynamics [20].
The model proposed here extends two-dimensional predator–prey systems with refuge and cannibalism by adding a dynamic refuge state. Instead of assuming a fixed refuge fraction, the refuge level evolves in response to intrinsic refuge growth, predator-induced activation, and natural decay. This produces a three-dimensional system in which population densities and interaction inhibition coevolve.
The novelty of the present study is threefold:
a dynamically evolving refuge variable is coupled simultaneously to predation and cannibalism through a global hunting-inhibition interpretation;
the positive equilibrium is reduced to a scalar feasibility condition, clarifying when coexistence is biologically meaningful;
the fractional stability boundary is computed and illustrated numerically, showing how the order α can act as a bifurcation parameter.
The paper is organized as follows. Section 2 introduces the model and gives a parameter table. Section 3 proves positivity, invariance, well-posedness, and boundedness. Section 4 derives boundary and coexistence equilibria. Section 5 presents local stability analysis. Section 6 develops the fractional Hopf-type condition. Section 7 gives the numerical method and simulations. Section 8 discusses the ecological interpretation, and Section 9 concludes the paper.
2 Model formulation and biological interpretation
Let N(t) denote prey density, P(t) predator density, and m(t) ∈ [0, 1] an adaptive refuge or hunting-inhibition level. The quantity m(t) is not interpreted as a literal physical shelter used identically by prey and predators. Rather, it is a dimensionless global inhibition index measuring the reduction of effective risky encounters. This index may combine spatial hiding, reduction of predator search efficiency, avoidance behavior, and reduction of predator conflict under environmental structure. Under this interpretation, the factor 1 − m(t) reduces both interspecific predation and intraspecific cannibalistic encounters.
The proposed Caputo fractional-order model is
where 0 < α ≤ 1 and
The Caputo derivative is
For α = 1, system (Equation 1) reduces to the corresponding classical system. The variables, parameters, biological meanings, and assumptions used in the model are summarized in Table 1.
Table 1
| Symbol | Biological meaning | Assumption |
|---|---|---|
| N(t) | Prey density | N ≥ 0 |
| P(t) | Predator density | P ≥ 0 |
| m(t) | Global refuge/hunting-inhibition level | 0 ≤ m ≤ 1 |
| r | Intrinsic prey growth rate | r > 0 |
| K | Prey carrying capacity | K > 0 |
| b1 | Maximum predation rate | b1 > 0 |
| k1 | Half-saturation constant for predation | k1 > 0 |
| c1 | Conversion rate from consumed prey to predator growth | c1 > 0 |
| c2 | Additional predator recruitment or background gain | c2 ≥ 0 |
| e | Predator mortality rate | e > 0 |
| b2 | Maximum cannibalism/self-limitation rate | b2 > 0 |
| k2 | Half-saturation constant for cannibalism | k2 > 0 |
| s | Intrinsic refuge/adaptation growth rate | s > 0 |
| η | Predator-induced refuge activation rate | η > 0 |
| h3 | Half-saturation constant for refuge activation | h3 > 0 |
| δ | Refuge decay rate | δ > 0 |
| α | Caputo fractional order/memory intensity | 0 < α ≤ 1 |
Parameters of system (Equation 1).
3 Basic analytical properties
Let
The following results show that system (Equation 1) is biologically meaningful in Ω.
Theorem 3.1 (Positive invariance). If (N0, P0, m0) ∈ Ω, then every solution of system (Equation 1) remains in Ω for all t ≥ 0 for which the solution exists.
Proof. The right-hand side of the first equation satisfies
Therefore the vector field is tangent to the boundary N = 0. Similarly,
so the boundary P = 0 is invariant.
For the refuge variable, the two relevant boundary signs are
and
Thus the vector field points inward or is tangent on the boundary faces m = 0 and m = 1. The fractional Nagumo-type invariance principle for Caputo systems implies that the closed region Ω is positively invariant. □
Theorem 3.2 (Existence and uniqueness). For every initial condition in Ω, system (Equation 1) admits a unique local solution.
Proof. The right-hand side of system (Equation 1) is continuously differentiable on every compact subset of Ω because all denominators are strictly positive:
Hence the vector field is locally Lipschitz on Ω. The standard existence and uniqueness theorem for Caputo fractional differential equations gives a unique local solution. □
Theorem 3.3 (A sufficient boundedness condition). Assume c2 < e. Then N(t) and P(t) are bounded on [0, ∞), and m(t) ∈ [0, 1] for all t ≥ 0.
Proof. Define
Using system (Equation 1), the predation transfer terms cancel and
Because the last term is nonpositive and c2 < e, there exist constants μ > 0 and C > 0 such that
The fractional comparison lemma gives
where Eα is the Mittag–Leffler function. Hence V(t) is bounded, and therefore both N(t) and P(t) are bounded. The invariance theorem already gives 0 ≤ m(t) ≤ 1. □
Remark 3.4. The condition c2 < e is a transparent sufficient dissipativity condition. More general boundedness conditions may be derived by imposing a uniform lower bound on the effective cannibalism term, but the above assumption is adequate for the parameter set used in the main simulations.
4 Equilibrium points
Equilibria of a Caputo autonomous system are obtained by setting the right-hand sides of Equation 1 equal to zero. Therefore they do not depend on α.
4.1 Boundary equilibria
The trivial equilibrium is
If P = 0, the prey equation gives N = 0 or N = K, and the refuge equation becomes
Thus the prey-only equilibria are
and, when s > δ,
A prey-free predator equilibrium has the form
It must satisfy
and
Consequently, a necessary feasibility condition is
where mP ∈ (0, 1) is determined by the refuge equation above.
4.2 Positive coexistence equilibrium
Let
denote a positive coexistence equilibrium. The notation E* will be used only for this positive equilibrium.
Set
From the prey equation with N > 0, one obtains
which implies
Substituting qP = A(N) into the predator equation gives
Hence
The corresponding refuge and predator coordinates are
It remains to satisfy the refuge equation. Therefore the positive equilibrium is obtained from the scalar equation
Theorem 4.1 (Feasibility of coexistence). System (Equation 1) has a positive coexistence equilibrium if there exists Nc ∈ (0, K) such that
where q(N) and G(N) are defined by Equations 2, 3. In that case,
This formulation makes explicit what is meant by the positive equilibrium and gives directly testable positivity conditions.
5 Local stability analysis
Let denote the vector field of Equation 1. The Jacobian matrix J(N, P, m) = DF(N, P, m) is
where, with q = 1 − m and D = k2 + qP,
Let λ1, λ2, λ3 be the eigenvalues of J(E). For 0 < α ≤ 1, Matignon's criterion states that E is locally asymptotically stable if and only if
5.1 Boundary equilibria
At E0 = (0, 0, 0),
so
Since r > 0, the trivial equilibrium is unstable.
At E1 = (K, 0, 0), the eigenvalues are
Thus E1 is locally asymptotically stable when
The first inequality means that predators cannot invade the prey-only state.
At E2 = (K, 0, 1 − δ/s), which exists for s > δ, the eigenvalues are
Therefore E2 is locally asymptotically stable when
Again, the second condition is a predator-invasion condition, now reduced by the refuge factor δ/s.
5.2 Positive equilibrium
At the coexistence equilibrium Ec, the characteristic polynomial is
where
For the integer-order system, the Routh–Hurwitz conditions
are sufficient and necessary for local asymptotic stability. For the fractional system, however, the decisive condition is Equation 4. The coefficients ai are useful, but the eigenvalue arguments must be checked directly.
6 Fractional Hopf-type stability boundary
Fractional systems may lose or gain local stability as the order α changes. Define
By Matignon's criterion, the equilibrium E is stable when
If a pair of complex eigenvalues λ1,2 = u ± iv with v ≠ 0 determines the minimum angle and u > 0, then
The use of atan2 is important because it gives the correct quadrant of the eigenvalue argument.
Theorem 6.1 (Fractional Hopf-type transition). Assume that a coexistence equilibrium Ec(μ) depends smoothly on a parameter μ, and that the Jacobian has a simple conjugate pair λ1,2(μ) = u(μ) ± iv(μ) with v(μ) ≠ 0. If
and
then the equilibrium crosses the fractional stability boundary at μ = μ0. Under the usual nondegeneracy hypotheses, this crossing corresponds to a Hopf-type transition in the fractional-order system.
The biological interpretation is direct. Smaller values of α enlarge the fractional stability wedge and may suppress oscillations. Larger values of α approach the integer-order stability requirement and may allow oscillatory growth.
7 Numerical method and simulations
The nonlinear fractional system is solved by the fractional Adams–Bashforth–Moulton predictor–corrector method. For
the equivalent Volterra integral equation is
Let tn = nh. The predictor is
and the corrector is
where
and, for 1 ≤ j ≤ n − 1,
7.1 Baseline parameter set
Unless otherwise stated, the numerical simulations use
The initial condition is
with h = 0.05 and T = 20.
For this parameter set, the positive equilibrium is
and the eigenvalues are
Thus
Consequently, this coexistence equilibrium is locally stable only for orders below this threshold. For α = 0.8, the simulated trajectories should therefore be interpreted as transient dynamics rather than as evidence of local asymptotic stability of Ec. For the baseline parameter set, Figures 1–4 show the main numerical responses. Figure 1 displays the prey density for different fractional orders, Figure 2 shows the corresponding predator dynamics, Figure 3 reports the refuge activation variable, and Figure 4 presents the three-dimensional phase-space trajectory for α = 0.8.
Figure 1
Figure 2
Figure 3
Figure 4
7.2 Effects of refuge activation and refuge decay
The refuge activation parameter η controls how strongly predator density increases the hunting-inhibition level. Figure 5 shows that changing η modifies prey recovery after the initial predator expansion.
Figure 5
The decay parameter δ controls how quickly refuge disappears in the absence of sustained activation. Figure 6 shows that increasing δ lowers the refuge level and therefore changes the long-term interaction intensity.
Figure 6
7.3 Numerical illustration of the fractional Hopf-type threshold
To illustrate the fractional stability boundary, consider the parameter set
The coexistence equilibrium is
and the eigenvalues are
The corresponding critical fractional order is
Hence the coexistence equilibrium is locally stable for α < α* and loses fractional stability for α > α*. This transition is illustrated by the spectral-margin curve in Figure 7.
Figure 7
8 Discussion
The revised analysis clarifies the ecological and mathematical role of the refuge variable. A single variable m(t) cannot literally represent all details of prey hiding and predator cannibalism avoidance. It should instead be interpreted as a global inhibition index that reduces the encounter intensity of risky interactions. This modeling choice keeps the dimension low while allowing feedback between predator pressure and interaction suppression.
The equilibrium analysis also shows why it is important to distinguish between transient numerical patterns and asymptotic stability. For the baseline parameter set, the positive equilibrium exists but its local fractional stability threshold is low. Therefore simulations at α = 0.8 describe a biologically meaningful transient response, but they do not by themselves prove convergence to the positive equilibrium. This correction strengthens the interpretation of the numerical results.
The fractional order modifies the stability wedge. In the Hopf-type example, the critical value α* ≈ 0.8410 shows that memory may stabilize the coexistence state for smaller orders and destabilize it as the system approaches the integer-order limit. This gives a precise mathematical meaning to the statement that memory can suppress oscillatory tendencies.
9 Conclusions
A three-dimensional fractional predator–prey model with adaptive refuge and predator cannibalism has been analyzed. The model is well posed in the biologically feasible region, and a sufficient boundedness condition was provided. Boundary equilibria were classified, and the positive coexistence equilibrium was expressed through a scalar feasibility equation.
The local stability analysis was rewritten using the correct Jacobian matrix and Matignon's criterion. The notation E* was fixed to denote the positive coexistence equilibrium only. The Hopf section was strengthened by introducing the fractional stability boundary
Numerical simulations now compare different fractional orders, include a three-dimensional phase-space trajectory, and examine the effects of refuge activation and refuge decay. A separate Hopf-type parameter set demonstrates the crossing of the fractional stability threshold.
Future work may consider separate refuge variables for prey evasion and cannibalism avoidance, spatial diffusion, noncommensurate fractional orders, parameter estimation from ecological data, and optimal control strategies based on refuge modulation.
Statements
Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.
Author contributions
ÁS: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. LM: Formal analysis, Investigation, Methodology, Validation, Writing – original draft, Writing – review & editing. FG: Conceptualization, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing – original draft, Writing – review & editing.
Funding
The author(s) declared that financial support was not received for this work and/or its publication.
Acknowledgments
The authors thank the reviewers for their careful reading and constructive recommendations.
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that Generative AI was used in the creation of this manuscript.
The author(s) used Generative AI tools only for language editing, proofreading, formatting assistance, and preparation of responses to production queries. All mathematical derivations, numerical results, interpretations, and conclusions were checked and approved by the authors, who take full responsibility for the content of the published article.
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.
References
1.
LotkaAJ. Elements of Physical Biology. Baltimore, MD: Williams & Wilkins (1925).
2.
VolterraV. Fluctuations in the abundance of a species considered mathematically. Nature. (1926) 118:558–60. doi: 10.1038/118558a0
3.
HollingCS. Some characteristics of simple types of predation and parasitism. Can Entomol. (1959) 91:385–98. doi: 10.4039/Ent91385-7
4.
MondalSSamantaG. Dynamics of predator-prey system with prey refuge and harvesting. Phys A. (2019) 534:122301. doi: 10.1016/j.physa.2019.122301
5.
YaseenRMHelalMMDehingiaKMohsenAA. Effect of the fear factor and prey refuge in an asymmetric predator-prey model. Braz J Phys. (2024) 54:214. doi: 10.1007/s13538-024-01594-9
6.
HameedAAl-HusseinHF. The impact of fear and refuge on the dynamics of predator-prey model: Stability and simulation. Partial Differ Equ Appl Math. (2025) 13:101029. doi: 10.1016/j.padiff.2024.101029
7.
DengHChenFZhuZLiZ. Dynamic behaviors of Lotka-Volterra predator-prey model incorporating predator cannibalism. Adv Differ Equ. (2019) 2019:359. doi: 10.1186/s13662-019-2289-8
8.
ZhangFChenYLiJ. Dynamical analysis of a stage-structured predator-prey model with cannibalism. Math Biosci. (2019) 307:33–41. doi: 10.1016/j.mbs.2018.11.004
9.
RayungsariMSuryantoAKusumawinahyuWMDartiI. Dynamics analysis of a predator-prey fractional-order model incorporating predator cannibalism and refuge. Front Appl Math Stat. (2023) 9:1122330. doi: 10.3389/fams.2023.1122330
10.
PodlubnyI. Fractional Differential Equations. San Diego, CA: Academic Press (1999).
11.
KilbasAASrivastavaHMTrujilloJJ. Theory and Applications of Fractional Differential Equations. Amsterdam: Elsevier (2006).
12.
PanigoroHSSuryantoADartiIKilicmanA. A fractional-order predator-prey model with ratio-dependent functional response and linear harvesting. Mathematics. (2019) 7:1100. doi: 10.3390/math7111100
13.
MatignonD. Stability results for fractional differential equations with applications to control processing. in IMACS-SMC proceedings, Vol. 2, Lille. (1996), p. 963–8.
14.
LiCWuG. Stability of fractional differential equations and applications to fractional-order systems. Nonlinear Dyn. (2010) 62:895–904.
15.
Vargas-De-LeónC. Volterra-type Lyapunov functions for fractional-order epidemic systems. Commun Nonlinear Sci Numer Simul. (2015) 24:75–85. doi: 10.1016/j.cnsns.2014.12.013
16.
DiethelmKFordNJFreedAD. Detailed error analysis for a fractional Adams method. Numer Algorithms. (2004) 36:31–52. doi: 10.1023/B:NUMA.0000027736.85078.be
17.
DiethelmKFordNJFreedAD. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn. (2002) 29:3–22. doi: 10.1023/A:1016592219341
18.
HassanSSGoldarSMohsenAAPatiR. Dynamics and bifurcation analysis of a discrete predator-prey model with dual Allee effects. Int J Bifurcat Chaos. (2026) 36:2650038. doi: 10.1142/S0218127426500380
19.
HuoHFZhaoXQZhuL. The global stability of fractional differential equations. Appl Math Lett. (2015) 45:14–9.
20.
GoldarSSardarPBiswasSHassanSSPatiRMohsenAAet al. Nonlinear dynamics and bifurcation analysis in a modified discrete-time Rosenzweig-MacArthur predator-prey model. Comput Math Model. (2025) 36:414–48. doi: 10.1007/s10598-025-09626-y
Summary
Keywords
Adams–Bashforth–Moulton method, adaptive refuge, Caputo derivative, fractional Hopf bifurcation, fractional predator–prey model, matignon stability, predator cannibalism
Citation
Salas S. ÁH, Martínez H. LJ and Gallego L. FA (2026) Adaptive refuge and memory effects in a three-dimensional fractional predator–prey system with predator cannibalism. Front. Appl. Math. Stat. 12:1860657. doi: 10.3389/fams.2026.1860657
Received
20 April 2026
Revised
10 June 2026
Accepted
15 June 2026
Published
10 July 2026
Volume
12 - 2026
Edited by
Khalid Hattaf, Centre Régional des Métiers de l'Education et de la Formation (CRMEF), Morocco
Reviewed by
Ahmed Mohsen, University of Baghdad, Iraq
Maya Rayungsari, University of Brawijaya, Indonesia
Updates
Copyright
© 2026 Salas S., Martínez H. and Gallego L..
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Álvaro H. Salas S., ahsalass@unal.edu.co
Disclaimer
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.