Analysis of a fractional-order prey-predator model with prey refuge and predator harvest using the consumption number: Holling type III functional response

Nnaemeka Stanley Aguegboh; Netochukwu Onyiaji; Chibuike Alexander Okeke; Nnaji Ugochukwu Daniel; Okongo Walter; Boubacar Diallo

doi:10.1515/cmb-2025-0023

Artikel Open Access

Analysis of a fractional-order prey-predator model with prey refuge and predator harvest using the consumption number: Holling type III functional response

Nnaemeka Stanley Aguegboh , Netochukwu Onyiaji , Chibuike Alexander Okeke , Nnaji Ugochukwu Daniel , Okongo Walter und Boubacar Diallo

Veröffentlicht/Copyright: 3. Juni 2025

Veröffentlicht von

Veröffentlichen auch Sie bei De Gruyter Brill

Manuskript einreichen Informationen für Autor*innen

Aus der Zeitschrift Computational and Mathematical Biophysics Band 13 Heft 1

Abstract

In this article, we study a fractional-order prey-predator model incorporating prey refuge and predator harvesting, employing a Holling type III functional response. The model’s dynamics are explored using the consumption number C 0 , providing a comprehensive understanding of its ecological implications. We establish the existence and uniqueness of the solution, ensuring the model’s mathematical robustness. Stability analysis is conducted to determine the conditions under which the system remains in equilibrium. Sensitivity analysis is performed to assess the importance of various parameters on the system’s behaviour with respect to the consumption number. Additionally, we investigate the occurrence of transcritical bifurcations, utilizing the prey refuge as a bifurcation parameter. Numerical simulations are carried out to illustrate the theoretical results and to provide deeper insights into the model’s dynamics. Our findings contribute to the understanding of fractional-order systems in ecological modelling and offer valuable perspectives on the management of predator-prey interactions in the presence of prey refuge and predator harvesting.

Keywords: prey-predator; refuge; harvesting; stability; sensitivity analysis; transcritical bifurcation analysis; numerical simulations

MSC 2010: 37A50; 60H10; 60J65; 60J70

1 Introduction

In ecological modelling, understanding the dynamics between prey and predator populations is crucial for managing ecosystems and conserving species. Traditional models often rely on integer-order differential equations to describe these interactions. However, fractional-order models, which incorporate non-integer derivatives, offer a more nuanced approach by capturing complex behaviours that standard models might overlook [25]. In this analysis, we examine a fractional-order prey-predator model that includes two significant factors: prey refuge and predator harvest. Prey refuge refers to safe spaces where prey can avoid predation, while predator harvest addresses the impact of human activities on predator populations [4].

The introduction of fractional derivatives in our model allows us to better represent the real-world dynamics of predator-prey systems. Fractional-order models are advantageous because they can describe memory effects and time delays that are inherent in ecological interactions [25]. By incorporating prey refuge, the model acknowledges that not all prey are equally vulnerable to predation. This adds a layer of realism, as some prey can find sanctuary in environments where predators cannot reach them. Meanwhile, the inclusion of predator harvest reflects the effects of human interventions, such as hunting or fishing, which can significantly alter predator populations and consequently impact the entire ecosystem [4].

Our analysis uses the consumption number, a key parameter that quantifies the rate at which predators consume prey, to evaluate the stability and behaviour of the model. This parameter helps us understand how changes in predator harvest and prey refuge affect the balance between predator and prey populations. By simulating various scenarios with different consumption numbers, we can predict how these ecological changes might influence the stability of the ecosystem. This approach not only provides insights into the dynamics of predator-prey interactions but also helps in developing strategies for sustainable management and conservation efforts.

Fractional calculus is a powerful mathematical tool used to model systems with memory and hereditary properties. Unlike traditional calculus, which deals with derivatives and integrals of integer orders, fractional calculus extends these concepts to non-integer orders. This extension is particularly useful in ecological modelling, where interactions between species are often influenced by past events and experiences [25]. For example, the rate at which predators consume prey can depend not only on the current prey density but also on the history of prey availability and predation pressure. By incorporating fractional calculus into predator-prey models, researchers can more accurately capture these memory effects and predict more complex dynamics [4]. This approach allows for the modelling of systems where the rate of change at any given time is influenced by a weighted history of prior states, providing a more comprehensive understanding of population dynamics.

The concept of prey refuge is essential in understanding how predators and prey interact. Prey refuge refers to safe areas where prey can hide and avoid being eaten by predators. These refuges can be physical places like burrows or thick vegetation, or behavioural changes like being active at night to avoid daytime predators. Having a refuge for prey helps stabilize the predator-prey relationship because it prevents prey populations from being wiped out by too much predation and allows predators and prey to live together more peacefully [12,16]. Mathematical models that include prey refuge show that when prey have places to hide, their populations tend to be more stable, with fewer dramatic swings in numbers [14]. This stability is important for real-world ecosystems, where refuges are common and help keep biodiversity and ecological balance intact [5].

Analysing a fractional-order prey-predator model with prey refuge and predator harvest using the consumption number is crucial for understanding predator-prey dynamics. Unlike traditional models that use whole numbers, fractional-order models use more complex mathematical techniques to better capture how real-world systems work, including delays and memory effects [1,20]. These models can include factors like safe spaces where prey can hide and the impact of human activities on predator populations. The consumption number helps by measuring how much prey is eaten by predators, allowing researchers to see how changes in this rate, along with prey refuge and predator harvest, affect the balance between predators and prey. This approach helps in creating effective conservation strategies and managing ecosystems to maintain their health and stability.

This analysis is a powerful tool for studying the complex relationships between predators and prey. By using advanced math and realistic ecological factors, this model provides valuable insights into how ecosystems function, which is useful for planning conservation efforts and managing habitats [4]. Understanding these details is important for predicting changes in animal populations, protecting biodiversity, and keeping natural environments healthy. Ongoing research and improvements to these models will help scientists better understand the delicate balance in predator-prey relationships and the key factors that influence their stability and resilience [7].

Sharma et al. [26] studied a predator-prey model with a focus on how the ratio of prey to predators and prey harvesting impact the system. They used a mathematical approach called center manifold theory to analyse the stability of the system and identify different types of bifurcations, which are changes in the system’s behaviour. They found several types of bifurcations, including transcritical, Neimark-Sacker, fold, and period-doubling bifurcations, as well as more complex ones like Bogdanov-Takens and Chenciner bifurcations. These require changing two parameters at once for the bifurcation to occur. They used extensive numerical simulations to support their theoretical findings, showing various complex behaviours, including chaotic patterns and specific periods of system behaviour. They also discovered that changing the harvesting parameter can lead to the extinction of the prey population.

Assila et al. [6] studied a fishing environment with two zones: one where fishing is allowed and one where preys are protected. They used a special mathematical model to see how fishing and toxic conditions affect the prey and predator populations. They analysed how stable the populations are locally and globally using different mathematical methods, including Magtinon’s theorem for local stability and Lyapunov’s function for global stability. Their findings were supported by computer simulations to show that their model works.

Ahmed et al. [5] considered a predator-prey model where both predators and prey are harvested. They looked at how the system behaves and changes, focusing on its stability and the conditions under which it becomes chaotic. They used bifurcation theory to show that the system experiences period-doubling and Neimark-Sacker bifurcations, which are specific types of changes in system behaviour. The researchers support their findings with numerical examples to show how changing the harvesting parameters affects the number of stable points and the occurrence of these bifurcations.

In Bourega’s study [9], they explored how predators and prey interact using a model inspired by the Kolmogorov model. Specifically, they looked at a variant that includes a feature called the allee effect and a specific way predators hunt known as the Holling type II functional response. Instead of using traditional differential equations, they applied fractional-order differential equations to model the system. They analysed the existence and uniqueness of solutions that are always non-negative and identified conditions under which stable points exist. They used a method called linearization to study how stable these points are locally. Additionally, this study addressed an optimization problem related to controlling this predator-prey system. They utilized Pontryagin’s maximum principle, a mathematical tool, to determine the best way to manage the prey’s refuge area, which protects some prey from being hunted, ensuring the system operates optimally under specified conditions.

Debnath et al. [12] analysed a predator-prey model where the prey and predators are harvested. The prey is harvested linearly, while the predators are harvested nonlinearly. They used uniform persistence theory to examine how both species can survive over a long time. They analysed the stability of various points in the model using the centre manifold theorem. They also looked at the global behaviour of the system, focusing on how both species coexist. The study includes different types of bifurcations (changes in system behaviour) of one and two parameters to understand how the system changes with varying conditions. Their findings show that the system behaves more realistically and complexly with harvesting than without. They used numerical simulations to support their theoretical analysis.

Xie and Zhang [30] investigated how fear of predators affects a system where prey can hide, using a Holling type III prey-predator model. They examined the system’s stability and found that fear can stabilize the system by stopping cycles and reducing predator numbers without causing their extinction. They also showed that prey shelters help keep predators around. They used Hopf bifurcation analysis to understand the changes in the system’s behaviour and confirmed their findings with computer simulations.

In their study, Barman et al. [7] developed a predator-prey model using fractional-order differential equations, considering two significant factors: predator-induced fear and prey refuge. They used fractional calculus to explore how these factors influence the dynamics of prey populations whose life cycles involve memory. They analytically validated the biological relevance and essential properties of their model, such as positivity and boundedness of solutions. Stability analysis of equilibrium points was systematically conducted, revealing that the system exhibits Hopf bifurcation around the interior equilibrium point concerning parameters related to predator fear and prey refuge. Through numerical simulations, they demonstrated that the system’s stability improves under fractional-order derivatives compared to integer-order derivatives. They observed that higher levels of predator-induced fear and prey refuge independently contribute to system stability in the integer-order model. Conversely, decreasing the model order shifts the system dynamics from instability to stability, especially under lower levels of fear and refuge. Overall, their findings underscore the role of memory effects in promoting stable coexistence within predator-prey interactions.

Selvam et al. [18] looked at a predator-prey system using a special type of math called fractional calculus. They focused on how prey hiding spots (refuges) affect this system over time, using a model that takes specific time steps. They analysed the system’s stability and found conditions under which the populations of both predators and prey remain balanced. Their theoretical findings were supported by numerical examples and visual aids such as graphs and diagrams to help explain the dynamics of the system.

Yang [31] developed a model that looks at how prey and predators interact, taking into account prey hiding places and how predators hunt. The study used the method of integrated semigroup and Hopf bifurcation theory to show that when a certain condition is met, the populations start to cycle in a regular pattern. This happens when a key parameter changes. The results showed that prey hiding places help keep both prey and predator populations balanced. This finding was supported by computer simulations.

Haque et al. [14] explored how prey safety zones affect predator-prey relationships using two mathematical models. They introduced a new idea where the number of prey in refuges is not just affected by predator encounters but also follows a Holling type II response. The second model included competition among predators. The researchers thoroughly examined both models, investigating various equilibrium states where predators and prey can coexist, and assessed their stability locally and globally. They used saddle-node bifurcation analysis across different parameter ranges to understand how the systems behave. Their findings showed that the dynamics of these ecosystems are heavily influenced by two key factors: the constant of prey refuge and intra-predator competition.

The motivation for this research is to better understand how predators and prey interact in nature, especially using new mathematical methods to make more accurate predictions [22,27]. The goal is to use advanced mathematical techniques to analyse how factors such as safe spaces for prey and human actions that affect predator populations impact these interactions. By applying fractional calculus to a prey-predator model with prey refuge and predator harvest [3,28], this research aims to fill a gap left by traditional models that may not fully capture complex real-world dynamics. The primary aim is to create a detailed model of how these factors influence the balance between predators and prey. Objectives include improving predictions of population changes, enhancing conservation strategies, and providing better insights for managing natural ecosystems.

This research on a fractional-order prey-predator model with prey refuge and predator harvest using the consumption number is crucial for understanding and managing the complex interactions between predators and prey in natural ecosystems [29]. Using advanced math called fractional calculus, the study aims to capture the intricate behaviours and delayed responses that traditional models might miss. The model also considers safe spaces where prey can hide and the effects of human activities on predator populations. Additionally, the prey refuge and predator harvest using the consumption number realistically show how these factors influence the balance between predator and prey populations. This research aims to develop a more accurate and detailed model of predator-prey dynamics. The goal is to provide better insights for effective conservation strategies and ecosystem management.

2 Model formulation and assumptions

We extend the classical prey-predator model with prey refuge used in [17]. The prey-predator model with prey refuge is given as

(2.1) d x d t = α x 1 − x k − β ( 1 − m ) x y 1 + a ( 1 − m ) x , d y d t = − γ y + c β ( 1 − m ) x y 1 + a ( 1 − m ) x ,

where x is the population density of the prey, y is the population density of the predator, α is the growth rate of the prey, k is the carrying capacity of prey, β is the removal rate of preys by predators, a is the removal rate of x when y is halved, c is the biomass conversion rate of x into y and γ is the reduction of y due to other factors. The constant proportion m ϵ [0, 1) of the prey can take refuge to avoid predation.

2.1 Assumptions of the modified model

The fractional prey-predator mathematical model is formulated based on the following assumptions:

Given a time t , let x ( t ) and y ( t ) represent the densities of prey and predator, respectively.
Assume that, in the absence of a predator, the prey population grows logistically to its carrying capacity K at an inherent growth rate τ .
Let α represent the predator’s conversion efficiency and ω the food-independent mortality rate.
The harvesting rate is denoted by γ y , and the percentage of the predator population that is harvested in a given amount of time is represented by the constant γ .
To evade predators, a constant proportion c ϵ [0, 1) of the prey can seek refuge. This exposes ( 1 − c ) x of the prey to potential predators.
β represents the conversion of prey biomass into predator biomass
η represents x when y is half.

We will examine how prey refuge and predator harvesting affects the predator-prey model with a Holling type III response function, which is the Lotka-Volterra type. The modified fractional model is given as:

(2.2) D t α c x ( t ) = τ x 1 − x K − α ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 , D t α c y ( t ) = α β ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 − ω y − γ y ,

where D α denotes the Caputo fractional derivative and 0 < α ≤ 1 . When α = 1 , the model reduces to the classical case. We give the initial value of model (2.2) as follows x ( t 0 ) = x 0 ≥ 0 and y ( t 0 ) = y 0 ≥ 0 .

2.2 Preliminaries [19,23]

Definition 1

The fractional derivative of order α in the Caputo sense of a function f : ℜ + → ℜ is given by

(2.3) D t α f ( t ) = 1 Γ ( α − n ) ∫ α t f ( n ) ( τ ) d τ ( t − τ ) α + 1 − n ( n − 1 < α ≤ n ) .

Definition 2

The Laplace transform of the Caputo derivative is given by

(2.4) ∫ 0 ∞ e − p t { D t α f ( t ) } d t = p α F ( p ) − ∑ k = 0 n − 1 p α − k − 1 f ( k ) ( 0 ) , ( n − 1 < α ≤ n ) .

Definition 3

The fractional integral of order α of a function f : ℜ + → ℜ is given by

(2.5) J α ( f ( x ) ) = 1 Γ ( α ) ∫ 0 x ( x − t ) α − 1 f ( t ) d t , α > 0 , x > 0 ,

Definition 4

The fractional integral of the Caputo fractional derivative of order α of a function f : ℜ + → ℜ is given by

(2.6) J α { D α f ( t ) } = f ( t ) − ∑ k = 0 n − 1 f ( k ) ( 0 ) t k k ! , t > 0 .

Definition 5

A two-parameter function of the Mittag-Leffler type is defined by the series expansion:

(2.7) E α , β ( z ) = ∑ k = 0 ∞ z k Γ ( α k + β ) , ( α , β > 0 ) .

2.3 Existence, uniqueness, and boundedness of solution

Theorem 2.1

All solutions of equation (2.2) are non-negative for any initial values ( x ( 0 ) , y ( 0 ) ) ϵ R + 2 .

Proof

Since D t α c x ( t ) = τ x 1 − x K − α ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 , D t α c x ( 0 ) = 0 if x ( 0 ) = 0 . D t α c x = 0 implies that the population density of the prey does not change. Then, we show that if x ( 0 ) > 0 , then x ( t ) ≥ 0 for every t > 0 . Suppose that the assertion is wrong, so there exists t * > 0 such that

(2.8) x ( t ) > 0 , 0 ≤ t < t * , x ( t ) = 0 , t = t * , x ( t ) < 0 , t ≥ t * .

From (2.2) and (2.8), we obtain that D t α c x = 0 , t = t * . Therefore, the population density of the prey does not change when t = t * . Following the previous assertion, x ( t ) = 0 , t = t * , so that x ( t ) = 0 , t > t * . This contradicts the assertion that x ( t ) < 0 for t > t * . Hence, x ( t ) ≤ 0 for all t > 0 . is true. In a similar manner, it can be proved that y ( t ) ≤ 0 for every t > 0 .

Next, we will apply a fixed-point theorem, such as Banach’s fixed-point theorem (also called the contraction mapping theorem), to demonstrate the existence and uniqueness of the solution of system (2.2) [24]. Rewriting equation (2.2), we obtain

(2.9) D t α c x ( t ) = τ x 1 − x K − α ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 , D t α c y ( t ) = α β ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 − ω y − γ y .

Introducing a function f on the right-hand side of the system,

(2.10) f 1 ( x , y ) = τ x 1 − x K − α ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 , f 2 ( x , y ) = α β ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 − ω y − γ y .

Thus, system (2.2) can be expressed as

(2.11) D t α c x ( t ) = f 1 ( x ( t ) , y ( t ) ) , D t α c y ( t ) = f 2 ( x ( t ) , y ( t ) ) .

Theorem 2.2

If a function F maps a complete metric space X into itself and is a contraction (i.e., there exists a constant 0 < k < 1 such that d ( T ( x ) , T ( y ) ) ≤ k d ( x , y ) for all x , y ∈ X ), then T has a unique fixed point.

The fractional integral of order α of a function g is defined as

(2.12) I α g ( t ) = 1 Γ ( α ) ∫ 0 t ( t − s ) α − 1 g ( s ) d s .

Given D t α c x ( t ) = f 1 ( x ( t ) , y ( t ) ) , this can be expressed as

(2.13) x ( t ) = x ( 0 ) + I α f 1 ( x ( t ) , y ( t ) ) ,

and in the same manner,

(2.14) y ( t ) = y ( 0 ) + I α f 2 ( x ( t ) , y ( t ) ) .

Let a ( t ) = ( x ( t ) , y ( t ) ) . Define the operator T as

(2.15) T ( z ) ( t ) = a ( 0 ) + I α f ( a ( t ) ) ,

where f ( a ) = ( f 1 ( x , y ) , f 2 ( x , y ) ) .

We must show that T is a contraction on an appropriate function space, such as C ( [ 0 , T ] , R 2 ) , in order to use Banach’s theorem with the sup norm.

Consider the space X = C ( [ 0 , T ] , R 2 ) with norm ‖ a ‖ = max t ∈ [ 0 , T ] ‖ a ( t ) ‖ .

If a 1 ( t ) = ( x 1 ( t ) , y 1 ( t ) ) and a 2 ( t ) = ( x 2 ( t ) , y 2 ( t ) ) , we need to show that

(2.16) ‖ T ( a 1 ) − T ( a 2 ) ‖ ≤ k ‖ a 1 − a 2 ‖ ,

for some 0 < k < 1 .

Proof

(2.17) f 1 ( a ) = τ x 1 − x K − α ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 ,

(2.18) f 2 ( a ) = α β ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 − ω y − γ y .

For all a = ( x , y ) , a ¯ = ( x ¯ , y ¯ ) ϵ Ω ,

(2.19) ∣ ∣ f ( a ) − f ( a ¯ ) ∣ ∣ ≤ ∣ f 1 ( a ) − f 1 ( a ¯ ) ∣ + ∣ f 2 ( a ) − f 2 ( a ¯ ) ∣ = τ x 1 − x K − α ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 − τ x ¯ 1 − x ¯ K − α ( 1 − c ) x ¯ 2 y ¯ 1 + η ( 1 − c ) x ¯ 2 + α β ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 − ω y − γ y − α β ( 1 − c ) x ¯ 2 y ¯ 1 + η ( 1 − c ) x ¯ 2 − ω y ¯ − γ y ¯ ≤ ∣ τ x − x ¯ ∣ + x 2 − x ¯ 2 K + α ( 1 − c ) x 2 y ( 1 + η ( 1 − c ) x ¯ 2 ) − α ( 1 − c ) x ¯ 2 y ¯ ( 1 + η ( 1 − c ) x 2 ) ( 1 + η ( 1 − c ) x 2 ) ( 1 + η ( 1 − c ) x ¯ 2 ) + α β ( 1 − c ) x 2 y ( 1 + η ( 1 − c ) x ¯ 2 ) − α β ( 1 − c ) x ¯ 2 y ¯ ( 1 + η ( 1 − c ) x 2 ) ( 1 + η ( 1 − c ) x 2 ) ( 1 + η ( 1 − c ) x ¯ 2 ) + ∣ ( ω + γ ) ( y − y ¯ ) ∣ ≤ τ + τ ( x + x ¯ ) K + ( x α ( 1 − c ) + β ) η ( 1 − c ) y ( 1 + η ( 1 − c ) x 2 ) ( 1 + η ( 1 − c ) x ¯ 2 ) ∣ x − x ¯ ∣ + ( ω + γ ) ( ∣ y − y ¯ ∣ ) .

It can be demonstrated in the discussion that the system solution (2.2) is bounded in Ω . There is a positive constant M = max { x , y } , ∀ t ≥ 0 . Therefore, we have

(2.20) ∣ ∣ f ( a ) − f ( a ¯ ) ∣ ∣ ≤ τ + 2 M K + ( M α ( 1 − c ) + β ) η ( 1 − c ) M ( 1 + η ( 1 − c ) ) ( 1 + η ( 1 − c ) ) ∣ x − x ¯ ∣ + ( ω + γ ) ∣ ( y − y ¯ ) ∣ = L 1 ∣ x − x ¯ ∣ + L 2 ∣ y − y ¯ ∣ ,

with

L 1 = τ + 2 M K + ( M α ( 1 − c ) + β ) η ( 1 − c ) M ( 1 + η ( 1 − c ) ) ( 1 + η ( 1 − c ) ) , L 2 = ( ω + γ ) .

By choosing a positive constant L = max { L 1 , L 2 } , we obtain

(2.21) ∣ ∣ f ( a ) − f ( a ¯ ) ∣ ∣ ≤ L ∣ ∣ a − a ¯ ∣ ∣ .

Given the Lipschitz continuity of f 1 and f 2 in their arguments (which can be shown under suitable regularity conditions), we can bound the difference as

(2.22) ‖ I α ( f 1 ( a ) − f 1 ( a ¯ ) ) ‖ ≤ L 1 ‖ a − a ¯ ‖ ,

(2.23) ‖ I α ( f 2 ( a ) − f 2 ( a ¯ ) ) ‖ ≤ L 2 ‖ a − a ¯ ‖ .

Since I α is a bounded linear operator and under the assumption that f 1 and f 2 are Lipschitz continuous with constants L 1 and L 2 , we have

(2.24) ‖ T ( a ) − T ( a ¯ ) ‖ ≤ L ‖ a − a ¯ ‖ .

Banach’s fixed-point theorem ensures a unique fixed point in the case when L < 1 , indicating the existence and uniqueness of the solution if T is a contraction.

To show boundedness, we examine the behaviour of the solutions x ( t ) and y ( t ) over time.

Consider system (2.2)

(2.25) D t α c x ( t ) = τ x 1 − x K − α ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 , D t α c y ( t ) = α β ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 − ω y − γ y .

For boundedness, we construct the following Lyapunov function:

(2.26) V ( x , y ) = x + y .

Taking the Caputo derivative of V ,

(2.27) D t α c V ( x , y ) = D t α c x ( t ) + D t α c y ( t ) .

Substituting the equations,

(2.28) D t α c V ( x , y ) = τ x 1 − x K − α ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 + α β ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 − ω y − γ y .

Simplifying (2.28) gives

(2.29) D t α c V ( x , y ) = τ x 1 − x K − ( ω + γ ) y .

For the solution to be bounded, D t α c V ( x , y ) ≤ 0 should hold in the region where x and y are large. Note that:

(2.30) c D t α V ( x , y ) ≤ τ K − ( ω + γ ) y .

As y → ∞ , − ( ω + γ ) y becomes larger, ensuring V ( x , y ) (and thus x and y ) remains bounded. This implies that there exists some constant M such that x ( t ) , y ( t ) ≤ M for all t ≥ 0 .□

3 Model analysis

3.1 Equilibria

System (2.2) has the following equilibrium points:

the prey-predator free equilibrium E 1 ( x , y ) = ( 0,0 ) ,
the predator free equilibrium E 2 ( x , y ) = ( K , 0 ) , and
the coexistence equilibrium E 3 ( x , y ) = ( x * , y * ) ,

where x * = ω + γ α ( 1 − c ) − ( ω + γ ) + η ( 1 − c ) and y * = β τ ω + γ ( 1 − c ) ( α β − η ( ω + γ ) ) 1 − ω + γ ( 1 − c ) ( α β − η ( ω + γ ) ) K ω + γ

Remark

For each of the equilibrium points E i (for i = 1 , 2, 3) to exist, the inequalities E 1 < E 2 < E 3 , 0 ≤ x ≤ K , and α ( 1 − c ) + η ( 1 − c ) > ( ω + γ ) must be satisfied.

3.2 Consumption number, C 0

The consumption number (also known as predation number) often denoted as C 0 can be defined as the rate at which predators convert consumed prey into new predators relative to the prey density. Based on the availability and consumption of the prey, it indicates the potential increase of the predator population.

We use the consumption number, C 0 , as a threshold parameter that provides a criterion for the system’s stability around the equilibrium points. In epidemiological models, this quantity can also be linked to the basic reproduction number, R 0 , [13]. The next-generation matrix method is used to calculate it in a similar way [11]:

(3.1) C 0 = α β ( 1 − c ) K 2 ( 1 + η ( 1 − c ) K 2 ) ( ω + γ ) .

In ecological terms, we can express C 0 as the combination of parameters that quantifies the prey consumption required for survival. So, C 0 = 1 implies that the predator uses prey biomass at a rate almost corresponding to their own biomass loss. For C 0 < 1 , less prey is consumed per unit of predator biomass loss. For C 0 > 1 , more prey is consumed per unit of predator biomass loss.

3.2.1 Stability analysis

We linearize the system around the equilibrium points and look at the Jacobian matrix’s eigenvalues in order to analyse stability.

The Jacobian matrix J of system (2.2) at any point ( x , y ) is given as

(3.2) J = ∂ f 1 ∂ x ∂ f 1 ∂ y ∂ f 2 ∂ x ∂ f 2 ∂ y ,

where

f 1 ( x , y ) = τ x 1 − x K − α ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 ,

f 2 ( x , y ) = α β ( 1 − c ) x 2 y 1 + η ( 1 − c ) x 2 − ω y − γ y .

The partial derivatives are

∂ f 1 ∂ x = τ 1 − 2 x K − 2 α ( 1 − c ) x y ( 1 + η ( 1 − c ) x 2 ) 2 , ∂ f 1 ∂ y = − α ( 1 − c ) x 2 1 + η ( 1 − c ) x 2 , ∂ f 2 ∂ x = 2 α β ( 1 − c ) x y ( 1 + η ( 1 − c ) x 2 ) 2 , ∂ f 2 ∂ y = α β ( 1 − c ) x 2 1 + η ( 1 − c ) x 2 − ( ω + γ ) .

At ( x , y ) = (0, 0),
J = τ 0 0 − ( ω + γ ) .
The eigenvalues are λ 1 = τ and λ = − ( ω + γ )
1. λ 1 = τ : If τ > 0 , then the equilibrium point (0, 0) is unstable
2. λ 2 = − ( ω + γ ) : Since ( ω + γ ) > 0 , this eigenvalue is negative, implying stability in the y -direction.
Therefore, (0, 0) is a saddle point and unstable.
At ( x , y ) = ( K , 0 ) ,
J = − τ − α ( 1 − c ) K 2 1 + η ( 1 − c ) K 2 0 α β ( 1 − c ) K 2 1 + η ( 1 − c ) K 2 − ( ω + γ ) .
The eigenvalues are
λ 1 = − τ ,

λ 2 = α β ( 1 − c ) K 2 1 + η ( 1 − c ) K 2 − ( ω + γ ) = ( ω + γ ) ( C 0 − 1 ) .
Hence, the equilibrium point E 2 ( x , y ) = ( K , 0 ) is stable if C 0 < 1 and unstable if C 0 > 1 .
At ( x , y ) = ( x * , y * ) ,
J = τ 1 − 2 x * K − 2 α ( 1 − c ) x * y * ( 1 + η ( 1 − c ) x 2 ) 2 − α ( 1 − c ) ( x * ) 2 1 + η ( 1 − c ) ( x * ) 2 2 α β ( 1 − c ) x * y * ( 1 + η ( 1 − c ) ( x * ) 2 ) 2 α β ( 1 − c ) ( x * ) 2 1 + η ( 1 − c ) ( x * ) 2 − ( ω + γ ) .
The characteristic equation is given as
τ − 2 x * τ K − 2 α ( 1 − c ) x * y * ( 1 + η ( 1 − c ) x 2 ) 2 − λ α β ( 1 − c ) ( x * ) 2 1 + η ( 1 − c ) ( x * ) 2 − ( ω + γ ) − λ + α ( 1 − c ) ( x * ) 2 1 + η ( 1 − c ) ( x * ) 2 2 α β ( 1 − c ) x * y * ( 1 + η ( 1 − c ) ( x * ) 2 ) 2 = 0 .
We rewrite the equation in the standard form:
( a 11 − λ ) ( a 22 − λ ) + a 12 a 21 = 0 ,
where
a 11 = τ − 2 x * τ K − 2 α ( 1 − c ) x * y * ( 1 + η ( 1 − c ) ( x * ) 2 ) 2 , a 22 = α β ( 1 − c ) ( x * ) 2 1 + η ( 1 − c ) ( x * ) 2 − ( ω + γ ) , a 12 = α ( 1 − c ) ( x * ) 2 1 + η ( 1 − c ) ( x * ) 2 , a 21 = 2 α β ( 1 − c ) x * y * ( 1 + η ( 1 − c ) ( x * ) 2 ) 2 .
The characteristic polynomial is
Δ ( λ ) = λ 2 + ( a 11 + a 22 ) λ + ( a 11 a 22 − a 12 a 21 ) .
Thus, the coefficients of the characteristic polynomial are
λ 2 + b 1 λ + b 0 = 0 ,
where
b 1 = a 11 + a 22 ,

b 0 = a 11 a 22 − a 12 a 21 .
To apply the Routh-Hurwitz criterion, we need to check the following conditions for stability:

1. b 1 > 0 and 2. b 0 > 0 .

Let us compute these coefficients:
b 1 = τ − 2 x * τ K − 2 α ( 1 − c ) x * y * ( 1 + η ( 1 − c ) ( x * ) 2 ) 2 + α β ( 1 − c ) ( x * ) 2 1 + η ( 1 − c ) ( x * ) 2 − ( ω + γ ) , b 0 = τ − 2 x * τ K − 2 α ( 1 − c ) x * y * ( 1 + η ( 1 − c ) ( x * ) 2 ) 2 α β ( 1 − c ) ( x * ) 2 1 + η ( 1 − c ) ( x * ) 2 − ( ω + γ ) − α ( 1 − c ) ( x * ) 2 1 + η ( 1 − c ) ( x * ) 2 2 α β ( 1 − c ) x * y * ( 1 + η ( 1 − c ) ( x * ) 2 ) 2 .
Stability conditions
1. Condition 1: b 1 > 0
  τ − 2 x * τ K − 2 α ( 1 − c ) x * y * ( 1 + η ( 1 − c ) ( x * ) 2 ) 2 + α β ( 1 − c ) ( x * ) 2 1 + η ( 1 − c ) ( x * ) 2 − ( ω + γ ) > 0 .
2. Condition 2: b 0 > 0
  τ − 2 x * τ K − 2 α ( 1 − c ) x * y * ( 1 + η ( 1 − c ) ( x * ) 2 ) 2 α β ( 1 − c ) ( x * ) 2 1 + η ( 1 − c ) ( x * ) 2 − ( ω + γ ) > α ( 1 − c ) ( x * ) 2 1 + η ( 1 − c ) ( x * ) 2 2 α β ( 1 − c ) x * y * ( 1 + η ( 1 − c ) ( x * ) 2 ) 2 .
  By analysing the signs of b 1 and b 0 using the Routh-Hurwitz criterion, we can determine the stability of the equilibrium point. If both conditions are met, the system is stable. If either condition is not met, the system is unstable. For a Matlab code that computes this stability easily, see Appendix.

3.2.2 Sensitivity analysis

This section shows how varying the parameter values affects the consumption number, R 0 . Finding the important parameter that might be an essential threshold for preserving the prey (by reducing the predation number) in an ecosystem is imperative.

The following are the mathematical representations of the consumption number, C 0 , sensitivity index towards the parameters α , β , c , K , η , ω , and γ :

∂ C 0 ∂ α = β ( 1 − c ) K 2 ( 1 + η ( 1 − c ) K 2 ) ( ω + γ ) , ∂ C 0 ∂ β = α ( 1 − c ) K 2 ( 1 + η ( 1 − c ) K 2 ) ( ω + γ ) , ∂ C 0 ∂ c = − α β K 2 ( 1 + η ( 1 − c ) K 2 ) ( ω + γ ) 2 , ∂ C 0 ∂ K = 2 α β ( 1 − c ) K ( 1 + η ( 1 − c ) K 2 ) ( ω + γ ) , ∂ C 0 ∂ η = − α β ( 1 − c ) 2 K 4 ( 1 + η ( 1 − c ) K 2 ) 2 ( ω + γ ) , ∂ C 0 ∂ ω = α β ( 1 − c ) K 2 ( ω + γ ) 2 ⋅ ( 1 + η ( 1 − c ) K 2 ) ∂ C 0 ∂ γ = − α β ( 1 − c ) K 2 ( 1 + η ( 1 − c ) K 2 ) ( ω + γ ) 2 .

It can be inferred that some derivatives seem positive and that the consumption number, C 0 , increases as any of the aforementioned positive value factors, α , β , and K increase. The elasticity for each parameter is shown in the following:

E C 0 , α = ∂ C 0 ∂ α × α C 0 = 1 , E C 0 , β = ∂ C 0 ∂ β × β C 0 = 1 , E C 0 , c = ∂ C 0 ∂ c = − c 1 − c , E C 0 , K = ∂ C 0 ∂ K = 2 [ ( 1 + η ( 1 − c ) K 2 ) ( ω + γ ) − η ( 1 − c ) K ( ω + γ ) ] K ( 1 − K ) , E C 0 , η = ∂ C 0 ∂ η = − 2 η ( 1 − c ) K 2 ( ω + γ ) 1 + η ( 1 − c ) K 2 , E C 0 , ω = ∂ C 0 ∂ ω = − ω ω + γ , E C 0 , γ = ∂ C 0 ∂ γ = − γ ω + γ .

It is evident that E C 0 , α , E C 0 , β , and E C 0 , K are all positive as a result. This suggests that the consumption number will increase in value when the values of these parameters α , β , and K are increased. Table 1 shows the sensitivity indices of C 0 with respect to each of the parameter while Figure 1 shows the mesh plot that agrees with the sensitivity indices. Each of the plots (a)–(d) corresponds to each of the parameter with positive sensitivity index.

Table 1

C 0 sensitivity indices

Parameters	Indices
α	1.0000
β	1.0000
c	− 0.0002325
K	0.00030997
η	− 0.9998
ω	− 0.77157
γ	− 0.2284

$Figure 1 Simulation of C 0 {C}_{0} with respect to the parameters. Source: Created by the authors.$

Figure 1

Simulation of C 0 with respect to the parameters. Source: Created by the authors.

3.2.3 Global sensitivity analysis

This sub-section contains the sensitivity analysis of system (2.2) using the Latin hypercube sampling (LHS) method, which is described in [8,21] with 1,000 simulations per run. The partial rank correlation coefficient (PRCC) is a sampling-dependent, efficient, and reliable method that performs global sensitivity analysis, allowing us to examine the monotonicity between model parameters and the model output, or in this case, when the consumption number, C 0 , is varied. Generally, the following formula can be used to obtain the standard correlation coefficient, or ρ , for any two variables, x and y :

ρ = ∑ i ( x i − x ¯ ) ( y i − y ¯ ) ∑ i ( x i − x ¯ ) 2 ∑ i ( y i − y ¯ ) 2 ,

where { ( x i , y i ) : x i ϵ x , y i ϵ y } are a set of paired sample data [8]. The tornado plot of the PRCCs for the parameters in our suggested model is displayed in Figure 2. When parameters with negative PRCCs are increased, the consumption decreases, and when parameters with positive PRCCs are increased, the basic consumption number, C 0 , increases.

Figure 2

Tonado plot displaying the PRCCs of system (2.2). Source: Created by the authors.

3.2.4 Transcritical bifurcation analysis

By varying the parameter around each equilibrium point, we examine the dynamic behaviour of system (2.2). Sotomayor’s theorem [15,32] applications for local bifurcation are being modified. This theorem states that the given system’s Jacobian matrix has transcritical bifurcation at the predator-free equilibrium point E 2 ( x , y ) = ( K , 0 ) . We choose the bifurcation point c = c * and C 0 = 1 . It is easy to verify that the system’s Jacobian matrix at ( E 2 , c * ) can be computed as J = D f ( E 2 , c * )

(3.3) J = − τ − α ( 1 − c * ) K 2 1 + η ( 1 − c * ) K 2 0 α β ( 1 − c * ) K 2 1 + η ( 1 − c * ) K 2 − ( ω + γ ) ,

where c * = 1 − ω + γ K 2 ( α β − η ( ω + γ ) ) .

From the Jacobian matrix (3.3), one of the eigenvalues is zero in the direction of y (since C 0 = 1 ) and the other is negative in the direction of x ( λ x = − τ ) .

Furthermore, the eigenvector W = ( w 1 , w 2 , w 3 ) T , which corresponds to λ y meet the condition J w = λ w then J w = 0 ,

− τ − α ( 1 − c * ) K 2 1 + η ( 1 − c * ) K 2 0 α β ( 1 − c * ) K 2 1 + η ( 1 − c * ) K 2 − ( ω + γ ) w 1 w 2 = 0 0 ,

which yields

(3.4) − τ w 1 − α ( 1 − c * ) K 2 1 + η ( 1 − c * ) K 2 w 2 = 0 ,

(3.5) α β ( 1 − c * ) K 2 1 + η ( 1 − c * ) K 2 − ( ω + γ ) w 2 = 0 .

Solving 3.4 gives

(3.6) w 1 = − b 1 w 2 τ ,

where

b 1 = α ( 1 − c * ) K 2 1 + η ( 1 − c * ) K 2 .

Therefore,

W = b 11 w 2 w 2 = w 2 b 11 1 ,

where b 11 = − b 1 τ

In the same vein, the eigenvector Z = ( z 1 , z 2 , z 3 ) T , which corresponds to λ x of J T , can be written as

0 0 − b 1 τ z 1 z 2 = 0 0 .

Hence, we obtain

(3.7) z 1 = − τ z 2 b 11 .

Therefore,

Z = b 12 z 2 z 2 = z 2 b 12 1 ,

where b 12 = − τ b 11 .

System (2.2) could be written in the form

d x d t = f ( x ) ,

where X = ( x , y ) T and f ( f 1 , f 2 ) T , then we find d f d c = f c .

We obtain that

f c = − α x 2 y + 2 α η ( 1 − c ) x 4 y ( 1 + η ( 1 − c ) x 2 ) 2 − α β x 2 y ( 1 + η ( 1 − c ) x 2 ) 2 .

Then,

f c ( E 2 , c * ) = 0 0 .

Z T ⋅ f c ( E 2 , c * ) = 0 ,

D f c ( E 2 , c * ) = 0 α ( 2 K ( 1 + η ( 1 − c ) K 2 ) − 4 η ( 1 − c ) K 3 ) ( ( 1 + η ( 1 − c ) x 2 ) ) 3 0 − α K 2 + 2 α η ( 1 − c ) K 4 ( 1 + η ( 1 − c ) K 2 ) 2 ,

Z T ⋅ [ D f c ( E 2 , c * ) ⋅ W ] = z 2 w 2 α K 2 + 2 α η ( 1 − c ) K 4 ( 1 + η ( 1 − c ) K 2 ) 2 ≠ 0 .

The transcritical bifurcation occurs in predator-free equilibrium when the parameter c passes through the bifurcation value c * , according to the Sotomayor theorem.

3.3 Numerical simulation

The Adams-Bashforth-Moulton strategy is the most widely used numerical technique for resolving fractional-order initial value problems [2]. Examine the following fractional differential equation:

(3.8) D t α c F i ( t ) = g j ( t , F j ( t ) ) , F j r ( 0 ) = F j 0 r ,

r = 0 , 1 , 2 , … , α , j ∈ N ,

where F j 0 r is the arbitrary real number, α > 0 , and the fractional differential operator D t α is linked to the well-known Volterra integral equation in the Caputo sense:

(3.9) F j ( t ) = ∑ n = 0 α − 1 F j 0 r t n n ! + 1 Γ ( α ) ∫ 0 t ( t − u ) α − 1 f j ( u , F j ( u ) ) d u , j ∈ N .

In this article, we study the numerical solution of system (2.2) model using the Adam’s-Bashforth-Moulton predictor-corrector scheme. The algorithm is shown in the following:

Let h = T m ˆ , t n = n h , n = 0 , 1 , 2 , … , m ˆ .

Corrector formulae:

x n + 1 = x 0 + h α Γ ( α + 2 ) τ x n + 1 p 1 − x n + 1 p K − α ( 1 − c ) x n + 1 2 p y n + 1 p 1 + η ( 1 − c ) x n + 1 2 p + h α Γ ( α + 2 ) ∑ j = 0 n α j , n + 1 τ x j 1 − x j K − α ( 1 − c ) x j 2 y j 1 + η ( 1 − c ) x j 2 ,

y n + 1 = y 0 + h α Γ ( α + 2 ) α β ( 1 − c ) x n + 1 p y n + 1 p 1 + η ( 1 − c ) x n + 1 2 p − ω y n + 1 p − γ y n + 1 p + h α Γ ( α + 2 ) ∑ j = 0 n α j , n + 1 α β ( 1 − c ) x j y j 1 + η ( 1 − c ) x j 2 − ω y j − γ y j .

Predictor formulae:

x n + 1 p = x 0 + 1 Γ ( α ) ∑ j = 0 n ζ j , n + 1 τ x j ( 1 − x j K ) − α ( 1 − c ) x j 2 y j 1 + η ( 1 − c ) x j 2 , y n + 1 p = y 0 + 1 Γ ( α ) ∑ j = 0 n ζ j , n + 1 α β ( 1 − c ) x j y j 1 + η ( 1 − c ) x j 2 − ω y j − γ y j ,

where

α j , n + 1 = n α + 1 − ( n − α ) ( n + 1 ) α , if j = 0 ( n − j + 2 ) α + 1 + ( n − j ) α + 1 − 2 ( n − j + 1 ) α + 1 , 0 ≤ j ≤ n , 1 , if j = 1 ,

and

ζ = h α α [ ( n + 1 − j ) α − ( n − j ) α ] , 0 ≤ j ≤ n .

Using the parameter values given in Table 2, we carried out the numerical simulations of the model system (2.2) in the Caputo sense to validate our analytical findings with x = 1,000 and y = 2,000 .

Table 2

Parameter values

Parameters	Value	Source
τ	0.18	[17]
K	898	[10]
α	0.6	Estimated
η	0.02	[17]
c	0–0.9	Estimated
β	0.02	Estimated
ω	1.52	Estimated
γ	0.3-0.45	[17]

3.3.1 Discussion of results

The dynamics of the prey-predator model are depicted in Figure 3, where the predator hunting is represented by γ = 0.9 , the prey refuge by c = 0.3 , and the consumption number by C 0 = 0.914 < 1 . It is clear that there is a substantial decline in the prey species at a low-prey refuge and a significant increase in the predator population, both of which lead to a significant decline in the prey species eventually due to increased predator hunting and the necessity of the prey species being eaten by the predators.

$Figure 3 Dynamics of prey-predator model for C 0 < 1 {C}_{0}\lt 1 with c = 0.3 c=0.3 and γ = 0.9 \gamma =0.9 . Source: Created by the authors.$

Figure 3

Dynamics of prey-predator model for C 0 < 1 with c = 0.3 and γ = 0.9 . Source: Created by the authors.

The behaviour of the prey-predator model is depicted in Figure 4, where the prey refuge increases from c = 0.34 to c = 0.9 with γ = 0.9 and C 0 < 1 . A gradual fall in the prey population is likely to result in a corresponding considerable decrease in the predator population. It also demonstrates that the prey population will not go extinct if C 0 < 1 with those values for c and γ .

$Figure 4 Dynamics of prey-predator model for C 0 < 1 {C}_{0}\lt 1 with c = 0.9 c=0.9 and γ = 0.9 \gamma =0.9 . Source: Created by the authors.$

Figure 4

Dynamics of prey-predator model for C 0 < 1 with c = 0.9 and γ = 0.9 . Source: Created by the authors.

The prey-predator model is plotted in Figure 5 with a consumption number greater than one ( C 0 = 1.12 ). This was accomplished by drastically reducing the prey refuge c = 0.3 , reducing the refuge hunt γ = 0.5 , and increasing the predator conversion efficiency, α = 1.5 . This graph demonstrates how the population of the prey will rapidly go extinct when the consumption number is more than 1, as they would be consumed by the predator, and how the predator population will gradually decline as a result of the preys’ disappearance from the ecosystem.

$Figure 5 Dynamics of prey-predator model for C 0 > 1 {C}_{0}\gt 1 with c = 0.4 c=0.4 and γ = 0.2 \gamma =0.2 . Source: Created by the authors.$

Figure 5

Dynamics of prey-predator model for C 0 > 1 with c = 0.4 and γ = 0.2 . Source: Created by the authors.

The prey-predator model’s dynamics when C 0 < 1 , with c = 0.1 and γ = 0.2 , are depicted in Figure 6. The graphs indicate that there will be a reduction in the populations of prey and predator if the consumption number is less than one, with prey refuge being 0.1 and predator harvest being 0.2. However, none of these species will die extinct; instead, they will coexist.

$Figure 6 Dynamics of prey population for C 0 < 1 {C}_{0}\lt 1 with c = 0.1 c=0.1 and γ = 0.2 \gamma =0.2 , different values of α \alpha . Source: Created by the authors.$

Figure 6

Dynamics of prey population for C 0 < 1 with c = 0.1 and γ = 0.2 , different values of α . Source: Created by the authors.

The variation of the predator and prey species, respectively, with varying values of α , where C 0 < 1 , is depicted in Figures 7 and 8. The main purpose of these graphs is to demonstrate how the dynamics of the population are impacted by variations in the fractional derivative’s order. It was evident that a lower derivative order value at any given time resulted in higher biomass levels in both the prey and predator populations.

$Figure 7 Dynamics of prey population for C 0 < 1 {C}_{0}\lt 1 at different values of α \alpha . Source: Created by the authors.$

Figure 7

Dynamics of prey population for C 0 < 1 at different values of α . Source: Created by the authors.

$Figure 8 Dynamics of predator population for C 0 < 1 {C}_{0}\lt 1 at different values of α \alpha . Source: Created by the authors.$

Figure 8

Dynamics of predator population for C 0 < 1 at different values of α . Source: Created by the authors.

3.3.2 Conclusion

This article presents a thorough analysis of a fractional-order prey-predator model with the incorporation of prey refuge and predator harvesting, using a Holling type III functional response. Our study confirms the existence and uniqueness of solutions, demonstrating the mathematical soundness of the proposed model. The stability analysis reveals the conditions under which the system achieves equilibrium, while the sensitivity analysis highlights the critical parameters affecting the system’s dynamics. By employing the consumption number as a bifurcation parameter, we successfully identify the occurrence of transcritical bifurcations, providing significant insights into the ecological behaviour of the system. The numerical simulations corroborate our theoretical findings, offering a detailed visual representation of the model’s dynamics and also shows that the order of derivative plays a vital role in modelling. Our work emphasizes the importance of fractional-order models in ecological studies, showing their potential to capture more complex and realistic behaviours compared to classical integer-order models. The insights gained from this analysis can aid in the development of more effective strategies for the management and conservation of ecological systems, particularly in scenarios involving prey refuge and predator harvesting. The study provides valuable insights for wildlife conservation, sustainable resource management, and ecological forecasting, advocating for protected habitats and controlled harvesting strategies. Future research should focus on validating the model with real-world ecological data, incorporating climate effects, and extending it to multi-species interactions to improve conservation efforts and ecosystem sustainability.

Acknowledgements

Not applicable.

Funding information: The authors state no funding involved.
Author contributions: NSA: conceptualization; formal analysis; validation; visualization; writing – original draft; and writing. NO: conceptualization; validation; review and editing. CAO: conceptualization; validation; review and editing. NUD: conceptualization; validation; review and editing. WO: conceptualization; supervision; review and editing. BD: conceptualization; validation; review and editing.
Conflict of interest: The authors state no conflict of interest.
Ethical approval: The conducted research is not related to either human or animals use.
Data availability statement: Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.

Appendix

References

[1] Aguegboh N. S., Okongo, W., Boubacar, D., Dasumani, M., Nnamaga, K. C., … Nwachukwu, A. U. (2025). A novel approach to modeling malaria with treatment and vaccination as control strategies in Africa using the Atangana-Baleanu derivative. Modeling Earth Systems and Environment, 11(2), 110, DOI: https://doi.org/10.1007/s40808-024-02273-z. Suche in Google Scholar

[2] Aguegboh, N. S., Phineas, K. R., Felix, M., & Diallo, B. (2023). Modeling and control of hepatitis B virus transmission dynamics using fractional order differential equations. Communications in Mathematical Biology and Neuroscience, 2023, Art. ID 107, DOI: http://doi.org/10.28919/cmbn/8174. Suche in Google Scholar

[3] Ahmed, R., Akhtar, S., Farooq, U., & Ali, S. (2023). Stability, bifurcation, and chaos control of predator-prey system with additive Allee effect. Communications in Mathematical Biology and Neuroscience, 2023, Article ID 9, DOI: https://doi.org/10.28919/cmbn/7824.Suche in Google Scholar

[4] Ahmed, R. Mushtaq, J., Saher, S., & Saeed, H. M. A. (2022). Dynamic analysis of a predator-prey model with Holling type-II functional response and prey refuge by using a NSFD scheme. Communications in Mathematical Biology and Neuroscience, 2022, 111, DOI: http://doi.org/10.289/cmbn/7735. Suche in Google Scholar

[5] Ahmed, R., Yazdani, M. S., & Saher, S. (2023). Complex dynamics of a predator-prey model with harvesting effects on both predator and prey. International Journal of Nonlinear Analysis and Applications, 14(10), 95–106, DOI: https://doi.org/10.22075/ijnaa.2023.25915.3161. Suche in Google Scholar

[6] Assila, C., Benazza, H., Lemnaouar, M. R., & Louartassi, Y. (2023). Harvesting of a fractional-order prey-predator model with Holling type ii functional response and reserved area for prey in the presence of toxicity. International Journal of Mathematical Modelling and Numerical Optimisation, 13(3), 275–296, DOI: https://doi.org/10.1504/IJMMNO.2023.132296. Suche in Google Scholar

[7] Barman, D., Roy, J., Alrabaiah, H., Panja, P., Mondal, S. P., & Alam, S. (2021). Impact of predator incited fear and prey refuge in a fractional order prey predator model. Chaos, Solitons & Fractals, 142, 110420, DOI: https://doi.org/10.1016/j.chaos.2020.110420. Suche in Google Scholar

[8] Blower, S. M., & Dowlatabadi, H. (1994). Sensitivity and uncertainty analysis of complex models of disease transmission: an HIV model, as an example. International Statistical Review/Revue Internationale de Statistique, 62(2), 229–243, DOI: https://doi.org/10.2307/1403510. Suche in Google Scholar

[9] Bourega, A. (2023). Stability analysis and optimization problem fractional of a predator-prey system with Holling II functional response. Asian-European Journal of Mathematics, 16(3), 2350048, DOI: https://doi.org/10.1142/S1793557123500481. Suche in Google Scholar

[10] Chakraborty, S., Pal, S., & Bairagi, N. (2012). Predator-prey interaction with harvesting: Mathematical study with biological ramifications. Applied Mathematical Modelling, 36(9), 4044–4059, DOI: https://doi.org/10.1016/j.apm.2011.11.029. Suche in Google Scholar

[11] Collins, O. C., & Duffy, K. J. (2016). Consumption threshold used to investigate stability and ecological dominance in consumer-resource dynamics. Ecological Modelling, 319, 155–162, DOI: https://doi.org/10.1016/j.ecolmodel.2015.03.021. Suche in Google Scholar

[12] Debnath, S., Majumdar, P., Sarkar, S., & Ghosh, U. (2022). Global dynamics of a prey-predator model with Holling type III functional response in the presence of harvesting. Journal of Biological Systems, 30(1), 225–260, DOI: https://doi.org/10.1142/S0218339022500073. Suche in Google Scholar

[13] den Driessche, P. V., & Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences, 180(1–2), 29–48, DOI: https://doi.org/10.1016/S0025-5564(02)00108-6. Suche in Google Scholar PubMed

[14] Haque, M., Rahman, M. D., Venturino, Z., & Li, B. L. (2014). Effect of a functional response-dependent prey refuge in a predator-prey model. Ecological Complexity, 20, 248–256, DOI: https://doi.org/10.1016/j.ecocom.2014.04.001. Suche in Google Scholar

[15] Khalaf K. Q., Majeed, A. A., & Naji, R. K. (2015). The local bifurcation and the hopf bifurcation for eco-epidemiological system with one infectious disease. Gen, 31(1), 18–41, http://www.geman.in. Suche in Google Scholar

[16] Ma, Z., Li, W., Zhao, Y., Wang, W., Zhang, H., & Li, Z. (2009). Effects of prey refuges on a predator-prey model with a class of functional responses: the role of refuges. Mathematical Biosciences, 218(2), 73–79, DOI: https://doi.org/10.1016/j.mbs.2008.12.008. Suche in Google Scholar PubMed

[17] Madueme, P.-G.U., Eze, V. O., & Aguegboh, N. S. (2021). Dynamics of prey predator model with prey refuge using a threshold parameter. Journal of Mathematics and Computer Science, 11(5), 5937–5946, DOI: https://doi.org/10.28919/jmcs/6184. Suche in Google Scholar

[18] Maria Selvam, A. G., Janagaraj, R., & Jacintha, M. (2020). Dynamical analysis of a discrete fractional order prey-predator system incorporating a prey refuge with Holling type II response. Journal of Physics: Conference Series, 1597, 012008, DOI: https://doi.org/10.1088/1742-6596/1597/1/012008. Suche in Google Scholar

[19] Nnaemeka, S., Ofe, U., Onyiaji, N. E., & Lovelyn, U. O. (2021). Fractional model on the dynamics of chicken pox with vaccination. International Journal of Mathematics Trends and Technology-IJMTT, 67, Article ID IJMTT-V67I8P501, DOI: https://doi.org/10.14445/22315373/IJMTT-V67I8P501. Suche in Google Scholar

[20] Nnaemeka, S. A., & Amanso, O. (2021). Analysis of a model on the transmission dynamics (with prevention and control) of Hepatitis B. Journal of Fractional Calculus and Applications, 12(1), 76–89, DOI: https://doi.org/10.21608/jfca.2021.308741. Suche in Google Scholar

[21] Okongo, W., Okelo, J. A., Gathungu, D. K., Moore, S. E., & Nnaemeka, S. A. (2024). Transmission dynamics of monkeypox virus with age-structured human population: A mathematical modeling approach. Journal of Applied Mathematics, 2024(1), 9173910, DOI: https://doi.org/10.1155/2024/9173910. Suche in Google Scholar

[22] Parwaliya, A., Singh, A., & Kumar, A. (2024). Hopf bifurcation in a delayed prey-predator model with prey refuge involving fear effect. International Journal of Biomathematics, 17(5), 2350042, DOI: https://doi.org/10.1142/S1793524523500420. Suche in Google Scholar

[23] Podlubny, I. (1998). Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications (Vol. 198). Elsevier. Suche in Google Scholar

[24] Rayungsari, M., Suryanto, A., Kusumawinahyu, W. M., & Darti, I. (2023). Dynamics analysis of a predator-prey fractional-order model incorporating predator cannibalism and refuge. Frontiers in Applied Mathematics and Statistics, 9, 1122330, DOI: https://doi.org/10.3389/fams.2023.1122330. Suche in Google Scholar

[25] Rihan, F. A., Lakshmanan, S., Hashish, A. H., Rakkiyappan, R., & Ahmed E. (2015). Fractional-order delayed predator-prey systems with Holling type-ii functional response. Nonlinear Dynamics, 80, 777–789, DOI: https://doi.org/10.1007/s11071-015-1905-8. Suche in Google Scholar

[26] Sharma, V. S., Singh, A., & Malik, P. (2024). Bifurcation patterns in a discrete predator-prey model incorporating ratio-dependent functional response and prey harvesting. Qualitative Theory of Dynamical Systems, 23(2), 74, DOI: https://doi.org/10.1007/s12346-023-00929-2. Suche in Google Scholar

[27] Singh, A., & Sharma, V. S. (2023). Bifurcations and chaos control in a discrete-time prey-predator model with Holling type-ii functional response and prey refuge. Journal of Computational and Applied Mathematics, 418, 114666, DOI: https://doi.org/10.1016/j.cam.2022.114666. Suche in Google Scholar

[28] Singh, A., & Sharma, V. S. (2023). Codimension-2 bifurcation in a discrete predator-prey system with constant yield predator harvesting. International Journal of Biomathematics, 16(5), 2250109, DOI: https://doi.org/10.1142/S1793524522501091. Suche in Google Scholar

[29] Suleman, A., Ahmed, R., Alshammari, F. S., & Shah, N. A. (2023). Dynamic complexity of a slow-fast predator-prey model with herd behavior. AIMS Mathematics, 8(10), 24446–24472, DOI: https://doi.org/10.3934/math.20231247. Suche in Google Scholar

[30] Xie, B., & Zhang, N. (2022). Influence of fear effect on a Holling type III prey-predator system with the prey refuge. AIMS Mathematics, 7(2), 1811–1830, DOI: https://doi.org/10.3934/math.2022104. Suche in Google Scholar

[31] Yang, P. (2019). Hopf bifurcation of an age-structured prey-predator model with Holling type II functional response incorporating a prey refuge. Nonlinear Analysis: Real World Applications, 49, 368–385, DOI: https://doi.org/10.1016/j.nonrwa.2019.03.014. Suche in Google Scholar

[32] Yu, X., Zhu, Z., Lai, L., & Chen, F. (2020). Stability and bifurcation analysis in a single-species stage structure system with Michaelis-Menten-type harvesting. Advances in Difference Equations, 2020, 1–18, DOI: https://doi.org/10.1186/s13662-020-02652-7. Suche in Google Scholar

Received: 2024-12-20

Revised: 2025-03-05

Accepted: 2025-03-21

Published Online: 2025-06-03

This work is licensed under the Creative Commons Attribution 4.0 International License.

Artikel in diesem Heft

https://doi.org/10.1515/cmb-2025-0023

Schlagwörter für diesen Artikel

prey-predator; refuge; harvesting; stability; sensitivity analysis; transcritical bifurcation analysis; numerical simulations

Creative Commons

BY 4.0