The study on the complex nature of a predator-prey model with fractional-order derivatives incorporating refuge and nonlinear prey harvesting

1.   Introduction

The use of mathematics to model and comprehend biological processes has grown significantly in recent years. The prey-predator framework is an essential tool in mathematical biology that aids in interpreting the dynamics of population interactions in natural environments. In order to provide predictions, these mathematical models are built using data gathered from experiments and observations. The specific dynamics are crucial for comprehending the effects on the environment throughout the modelling phase. Mathematical systems, mostly representing biological events using differential equations, encapsulate their dynamics. Not only do they have extensive uses in mathematics but also fields as diverse as engineering, biology, physics, and meteorology. All of ecology's activities are essentially the scientific investigation and creative analysis of the complex interactions that exist between living things and their surroundings. It includes not only living things but also the communities they build and the inanimate objects in their surroundings, all of which are constantly interacting with one another. The relationship between predators and their prey has always had significance in both the disciplines of ecology and mathematical ecology due to its global relevance ^[1,2,3]. A key component of population dynamics is the relationship between these two entities. Numerous investigators have worked to enhance the fundamental predator-prey system by including various biological factors, such as the prey's fear effect ^[4], prey refuge ^[5], selective and mixed-form harvesting ^[6,7,8], and many more.

There has been a lot of mathematical model building around the effects of prey refuge in recent years ^[9,10,11,12]. Wang et al. ^[13] used both analytical and numerical techniques to investigate an impact on prey refuge as well as predator fear on anti-predator activity. Prey refuges generally manifest in one of two ways: either the number of refuges is directly related to the density of prey population, or a fixed maximal capacity is in place ^{[14,15,16,17,18,19]}. In our model, we focus on the first kind, where $\varepsilon$($0 \leqslant \varepsilon < 1$) is the fixed maximum availability of refuge for prey. Another important aspect of predator-prey dynamics is the harvesting of either group. Fisheries, forests, and wildlife management are common places to find biological resources being harvested. In order to ensure that renewable resources are preserved for future generations, bio-economic models have been created to aid in the scientific management of these assets. The concept of predator-prey interactions is discussed in ^[20], where the predator incorporates a prey refuge and the prey engages in harvesting. In addition to discussing the impact of harvesting, Mukhopadhyay and Bhattacharyya ^[21] examined a scenario in which two predators were vying for only one prey. The predator-prey model with various harvesting scenarios is examined from an economic perspective in ^[22]. Nonlinear harvesting allows our predator-prey system to choose to harvest a single species while protecting the other one. In reality, we limit population growth by harvesting prey or predators. This formulates a predator-prey framework based on the Holling type Ⅱ dynamics, including the effects of nonlinear harvesting.

The above model is predicated on the following assumptions:

(ⅰ) Prey and predator population densities at a given instance of time $t$ are denoted by the variables ${x_1}\, {\text{and}}\, {x_2}$, respectively.

(ⅱ) The inherent growth rate and environmental carrying capacity of prey are $r\, {\text{and}}\, k$, respectively.

(ⅲ) The term $\frac{{b\, {x_1}\, {x_2}}}{{a + {x_1}}}$ denotes Holling type Ⅱ functional response that supports predator species progress, $a$ is the half-saturation constant, $b\left( { > 0} \right)$ is the predator's maximal consumption rate, ${b_1}\left( { > 0} \right)$ is the biomass transformation amount $\left( {0 < {b_1} < b} \right)$, and $c\left( { > 0} \right)$ represents the predator's inherent level of mortality. Here, we can see that if ${b_1} < c$, the right portion of the second equation of (1) never becomes positive. Therefore, we assume that ${b_1} - c > 0$ for model (1).

(ⅳ) The maximum prey harvesting rate is represented by the parameter $h$. With an increase in the harvesting population, the harvesting function $\frac{{h\, {x_1}}}{{h + {x_1}}}$ reaches gradually the saturation limit value $h$. From a biological point of view, we think this harvesting role makes more sense.

The suggested model (1) was built using a couple of first-order derivative nonlinear differential equations, in which the instantaneous population density determines how the population density changes over time. However, in actuality, a phenomenon known as the memory effect occurs, wherein the history of all past situations also influences the present state ^[23]. A fractional differential system may depict phenomena or systems with memory and genetic properties ^[24]. Drawing on the work of L'Hôpital and Leibniz, who provided a rigorous evaluation of derivatives of order 1/2, Liouville ^[25] developed the notion of fractional-order derivative. Riemann reinterpreted Liouville's formulation as the Riemann-Liouville fractional derivative operator ^[26], which he achieved by directly generalising the Cauchy formula. Euler's investigation of non-integer integration particularly inspired him to develop the Gamma function to be a generalisation of the factorial notion for non-integer numbers, which is used to explain the non-integer order derivative idea proposed by Liouville and Riemann ^[27]. When solving differential equations, Michele Caputo changed the Riemann-Liouville operator in 1967 to eliminate the need for beginning conditions. The modified operator's definition has become known as the Caputo fractional-order derivative operator. This derivative takes into account the memory effects and long-range dependencies present in many real-world systems, allowing for a more accurate representation of their behaviour. The Caputo fractional derivative can handle both smooth and non-smooth functions, making it applicable to a wide range of systems and phenomena. In addition, it possesses desirable mathematical properties, including linearity, compatibility with initial conditions, and the ability to preserve the order of differentiability. Many studies have been conducted on predator-prey frameworks involving non-integer derivatives of the Caputo type ^{[28,29,30,31,32,33,34,35]}. So, using the Caputo non-integer order derivative operator, we apply several modifications and analyses within the predator-prey framework (1) that include selective nonlinear harvesting and refuge for prey.

The following is the outline of the article. Section 2 outlines the process of developing the model and its fundamental attributes. One of the fundamental characteristics is the ability to confirm that the model's solutions are non-negative, exist, unique, and bounded. Section 3 presents the outcomes of the dynamical analysis. The findings include the presence and consistency of equilibrium points. The Hopf bifurcation is the one that is studied, and both local and global stability are examined. Verification of analytical findings is accomplished in Section 4 using numerical simulations and interpretations. In Section 5, we summarise our findings.

2.   Formulation and fundamental characteristics of the model

The model is created using the Caputo non-integer order derivative operator to the left side of the framework (1), which incorporates selective nonlinear harvesting and refuge for the prey.

with $\xi \in \mathbb{R}$, $0 < \xi \leqslant 1,$ and $^CD_0^\xi$ indicates the ${\xi ^{{\text{th}}}}$ order of Caputo non-integer derivative operator given by $D_t^\xi f\left( t \right) = \frac{1}{{\Gamma \left( {1 - \xi } \right)}}{\int\limits_0^t {\left( {t - \zeta } \right)} ^{ - \xi }}f\left( \zeta \right)\, d\zeta .$ All of system (2)'s parameters are defined in terms of system (1)'s assumptions.

Theorem 1. ^[36] Let $\frac{{{d^\xi }u}}{{d{t^\xi }}} = Au;\, \, u\left( 0 \right) = {u_0}$ be the linear and non-integer order differential equation, where $0 < \xi \leqslant 1, \, \, x \in {\mathbb{R}^n}$ and $A$ is any arbitrary matrix of size $n.$ In addition, there are:

(ⅰ) For $u = 0$ to be asymptotically resilient, it is necessary and sufficient that every latent value ${\lambda _i}$ of $A$ fulfil $\left| {\arg \left( {{\lambda _i}} \right)} \right| > \frac{{\xi \pi }}{2},$

(ⅱ) for $u = 0$ to be a stable solution, each latent value ${\lambda _i}$ of $A$ must satisfy $\left| {\arg \left( {{\lambda _i}} \right)} \right| \geqslant \frac{{\xi \pi }}{2}$, and all latent values that fulfil $\left| {\arg \left( {{\lambda _i}} \right)} \right| = \frac{{\xi \pi }}{2}$ must have the equal geometric and algebraic multiplicity.

Theorem 2. ^[37] Let $\frac{{{d^\xi }u}}{{d{t^\xi }}} = f\left( u \right);\, \, u\left( 0 \right) = {u_0},$ be the fractional order differential equation with $0 < \xi \leqslant 1, \, \, x \in {\mathbb{R}^n}$. The aforementioned system's equilibrium points may be found by solving $f\left( u \right) = 0.$ If the equilibrium points meets the criteria $\left| {\arg \left( {{\lambda _i}} \right)} \right| > \frac{{\xi \pi }}{2},$ then the system is locally asymptotically stable for each latent value $\lambda {\, _i}$ of the community matrix $J = \frac{{\partial f}}{{\partial u}}$.

Lemma 1. ^[38] Assume that $\xi$ is a member of $\left( {0, \, \, 1} \right].$ Let us take a function $f\left( t \right)$ that is a member of $C\left[ {0, \, \, a} \right]$ such that $D_t^\xi f\left( t \right) \in C\left[ {0, \, \, a} \right]$. If $D_t^\xi f\left( t \right) \geqslant 0$, for every $t \in \left( {0, \, \, a} \right),$ then it is non-decreasing function for all $t \in \left[ {0, \, \, a} \right]$. If $D_t^\xi f\left( t \right) \leqslant 0$, then it is non-increasing function for all $t \in \left[ {0, \, \, a} \right]$.

2.1.   Existence, uniqueness, non-negativity, and boundedness

Theorem 3. The framework (2) has a unique solution to every non-negative starting condition.

Proof. Establishing this theorem requires finding a unique solution to the dynamical framework within the domain $\left[ {0, \, \infty } \right) \times \, \chi , \,$ where $\chi = \left\{ {\left( {{x_1}, \, \, {x_2}} \right) \in {\mathbb{R}^2}:\max \left\{ {\left| {{x_1}} \right|, \, \left| {{x_2}} \right|} \right\} \leqslant T} \right\}.$

Let us take $V = \left( {{x_1}, \, \, {x_2}} \right), \, \, \bar V = \left( {{{\bar x}_1}, \, \, {{\bar x}_2}} \right)$ and now consider a mapping $H\left( V \right) = \left( {{H_1}\left( V \right), \, {H_2}\left( V \right)} \right),$ where

Then, for any $V, \, \, \bar V\, \in \chi$ we have

where $L = \max \left\{ {\left( {r + \frac{{2rT}}{k} + \frac{{{h^2}}}{{{{\left( {h + T} \right)}^2}}} + \frac{{a\left( {1 - \varepsilon } \right)T}}{{{{\left( {a + \left( {1 - \varepsilon } \right)T} \right)}^2}}}\left( {b + {b_1}} \right)} \right), \, \, \left( {c + \frac{{\left( {1 - \varepsilon } \right)T}}{{a + \left( {1 - \varepsilon } \right)T}}\left( {b + {b_1}} \right)} \right)\, } \right\}$.

Consequently, the Lipschitz criteria is satisfied by the mapping $H\left( V \right)$. Therefore, it may be concluded that, given the starting condition ${V_0} = \left( {{x_{{1_0}}}, \, \, {x_{{2_0}}}} \right)$, where ${x_{{1_0}}} \geqslant 0$ and ${x_{{2_0}}} \geqslant 0$, the existence and uniqueness theorem of ^[39] states that there is a single solution to the system of differential equations (2).

Theorem 4. There exist non-negative and uniformly bounded solutions to a pair of differential equations (2) that begin in the area $\mathbb{R}_ + ^2$.

Proof. We begin by determining if the solution, which originates in the domain $\mathbb{R}_ + ^2$, is non-negativity. Hence, it is necessary to establish that, for any $t \geqslant 0$, ${x_1}\left( t \right) \geqslant 0$, and ${x_2}\left( t \right) \geqslant 0$. Assume that not each $t \geqslant 0$ achieve the criteria of ${x_1}\left( t \right) \geqslant 0$. Then, there exist ${t_1} > 0$ such that ${x_1}\left( t \right) > 0$ for $t \in \left[ {0, \, \, {t_1}} \right), \, \, {x_1}\left( {{t_1}} \right) = 0$ and ${x_1}\left( {{t^ + }_1} \right) < 0.$ Subsequently, we may deduce from the initial equation in (2) that $D_t^\xi {x_1}\left( t \right){|_{t = {t_1}}} = 0.$ We find that ${x_1}\left( t \right) \geqslant 0$ for every $t \geqslant 0$ because Lemma 1 suggests that ${x_1}\left( {{t^ + }_1} \right) = 0$, which goes against the assumption ${x_1}\left( {{t^ + }_1} \right) < 0$. It is simple to establish that the remaining equation of (2) has a non-negative solution that begins within the domain of $\mathbb{R}_ + ^2$ using a similar procedure. Thus, for any $t \geqslant 0$, it provides ${x_2}\left( t \right) \geqslant 0$ in the same approach. Following that, we have to establish the uniform boundedness of the (2) solutions, which begin in the domain $\mathbb{R}_ + ^2$.

Considering $S\left( t \right) = {x_1}\left( t \right) + \frac{b}{{{b_1}}}{x_2}\left( t \right)$, we can use Eq (2) to get

It is deduced from Lemma 3 of ^[40] that

where ${E_\xi }$ represents the Mittag-Leffler function. By applying Lemma 5 and Corollary 6 from ^[41], we obtain $S\left( t \right) \leqslant \frac{{k{{\left( {r + c} \right)}^2}}}{{4r}}, \, \, t \to \infty$. Consequently, any solution to (2) that begins in $\mathbb{R}_ + ^2$ is confined to the region $\Omega = \left\{ {\left( {{x_1}\, , \, {x_2}} \right) \in \mathbb{R}_ + ^2:S \leqslant \frac{{k{{\left( {r + c} \right)}^2}}}{{4r}} + \delta , \, \, \delta > 0} \right\}.$

3.   Dynamical analysis

3.1.   Existence of equilibrium points

The steady state condition is solved at the equilibrium points of the framework (2), which are found by making $D_0^\xi {x_1}\left( t \right) = 0$ and $D_0^\xi {x_2}\left( t \right) = 0$. We can obtain the following equilibrium points:

(ⅰ) The extinction equilibrium ${E_0} = \left( {0, 0} \right),$

(ⅱ) the auxiliary equilibrium point ${E_1} = \left( {{x_1}, 0} \right),$ where $rx_1^2 + r\left( {h - k} \right){x_1} + kh\left( {1 - r} \right) = 0, \,$ ${x_1} = \frac{{r\left( {k - h} \right) + \sqrt {{r^2}{{\left( {h + k} \right)}^2} - 4rkh} }}{{2r}}$ and ${x_1}$ is positive for $k > h$, $\frac{{r{{\left( {h + k} \right)}^2}}}{{4kh}} > 1$,

(ⅲ) the interior equilibrium point ${E_2} = \left( {x_1^*, \, \, x_2^*} \right),$ where $x_1^* = \frac{{ac}}{{\left( {{b_1} - c} \right)\left( {1 - \varepsilon } \right)}}$ which is positive for $\varepsilon < 1$ and $x_2^* = \frac{{a{b_1}\left[ {kh\left( {r - 1} \right){{\left( {{b_1} - c} \right)}^2}{{\left( {1 - \varepsilon } \right)}^2} + rac\left( {k - h} \right)\left( {{b_1} - c} \right)\left( {1 - \varepsilon } \right) - r{a^2}{c^2}} \right]}}{{kb{{\left( {{b_1} - c} \right)}^2}{{\left( {1 - \varepsilon } \right)}^2}\left[ {h\left( {{b_1} - c} \right)\left( {1 - \varepsilon } \right) + ac} \right]}}$ which is positive for $r > 1$ and $\frac{{\left( {k - h} \right)\left( {{b_1} - c} \right)\left( {1 - \varepsilon } \right)}}{{ac}} > 1$.

3.2.   Local stability

The local stability of the equilibrium points of the framework (2) is shown by the findings of the community matrix's latent values and the application of the stability criteria theorem in Petras ^[37]. Assuming that ${E^*}$ is an equilibrium point of the framework (2). According to the Matignon stability criteria theorem, ${E^*}$ is local asymptotically stable provided each of the latent values ${\lambda _i}$ of a community matrix,

that satisfies $\left| {\arg \left( {{\lambda _i}} \right)} \right| > \frac{{\xi \pi }}{2}.$

Theorem 5. The origin point ${E_0}\left( {0, 0} \right)$ is locally asymptotically stable if $r < 1$ and saddle point if $r > 1.$

Proof. The Jacobian matrix for ${E_0} = \left( {0, 0} \right)$ is

The latent values of $J\left( {{E_0}} \right)$ are obtained as ${\lambda _1} = r - 1$ and ${\lambda _2} = - c < 0.$ When $r < 1$, then the eigenvalue ${\lambda _1}$ is not positive. Consequently, we get $\left| {\arg \left( {{\lambda _1}} \right)} \right| = \left| {\arg \left( {{\lambda _2}} \right)} \right| = \pi$, indicating that the criterion $\left| {\arg \left( {{\lambda _i}} \right)} \right| > \frac{{\xi \pi }}{2}$, $i = 1, \, \, 2$ has been achieved. Thus, ${E_0}$ is asymptotically stable. If $r > 1$, then one latent value is positive and another is negative, with $\left| {\arg \left( {{\lambda _1}} \right)} \right| = 0$ and $\left| {\arg \left( {{\lambda _2}} \right)} \right| = \pi$. Thus, we have the situation where one latent value meets the requirement $\left| {\arg \left( {{\lambda _i}} \right)} \right| > \frac{{\xi \pi }}{2}$, another meets the requirement $\left| {\arg \left( {{\lambda _i}} \right)} \right| < \frac{{\xi \pi }}{2}$. ${E_0}$ is the saddle point as a result.

Remark 1. According to the aforementioned study, the extinction equilibrium ${E_0}$ remains unstable $h = 0$ (the system has no harvesting).

Theorem 6. The equilibrium point ${E_1} = \left( {x_1^*, 0} \right)$ is locally asymptotically stable if $\sqrt {{{hk} / r}} - h < x_1^* < \frac{{ac}}{{\left( {1 - \varepsilon } \right)\left( {{b_1} - c} \right)}}$.

Proof. The community matrix for ${E_1}$ is

The latent values of $J\left( {{E_1}} \right)$ are ${\lambda _1} = - \left[ {\frac{{r{x_1}}}{k} - \frac{{h{x_1}}}{{{{\left( {h + {x_1}} \right)}^2}}}} \right]$, which is negative if ${x_1} > \left( {\sqrt {{{hk} / r}} - h} \right)$ and ${\lambda _2} = \frac{{{b_1}\left( {1 - \varepsilon } \right){x_1}}}{{a + \left( {1 - \varepsilon } \right){x_1}}} - c$, which is negative if ${x_1} < \frac{{ac}}{{\left( {1 - \varepsilon } \right)\left( {{b_1} - c} \right)}}$. Consequently, we get $\left| {\arg \left( {{\lambda _1}} \right)} \right| = \left| {\arg \left( {{\lambda _2}} \right)} \right| = \pi$, indicating that the criterion $\left| {\arg \left( {{\lambda _i}} \right)} \right| > \frac{{\xi \pi }}{2}$, $i = 1, \, \, 2$ has been achieved. Thus, ${E_1}$ is asymptotically stable if $\sqrt {{{hk} / r}} - h < x_1^* < \frac{{ac}}{{\left( {1 - \varepsilon } \right)\left( {{b_1} - c} \right)}}$.

Theorem 7. The positive equilibrium point ${E^*} = \left( {x_1^*, x_2^*} \right)$ is locally asymptotically stable if any of the conditions listed below are fulfilled:

ⅰ) ${\text{trace}}\left( {J\left( {{E^*}} \right)} \right) \leqslant 0$,

ⅱ) ${\text{trace}}\left( {J\left( {{E^*}} \right)} \right) > 0, \, \, {\text{trac}}{{\text{e}}^{\text{2}}}\left( {J\left( {{E^*}} \right)} \right) - 4\det \left( {J\left( {{E^*}} \right)} \right) < 0\,$ and

${\text{sqrt}}\left( {\left| {{\text{trac}}{{\text{e}}^{\text{2}}}\left( {J\left( {{E^*}} \right)} \right) - 4\det \left( {J\left( {{E^*}} \right)} \right)} \right|} \right) > {\text{trace}}\left( {J\left( {{E^*}} \right)} \right)\tan \left( {\frac{{\xi \pi }}{2}} \right) .$

Proof. The community matrix for ${E^*} = \left( {x_1^*, x_2^*} \right)$ is given by

Then, ${\text{trace}}\left( {J\left( {{E^*}} \right)} \right) = - \frac{{r{x_1}}}{k} + \frac{{b{{\left( {1 - \varepsilon } \right)}^2}{x_1}{x_2}}}{{{{\left( {a + \left( {1 - \varepsilon } \right){x_1}} \right)}^2}}} + \frac{{h{x_1}}}{{{{\left( {h + {x_1}} \right)}^2}}}$ and $\det \left( {J\left( {{E^*}} \right)} \right) = \frac{{ab{b_1}{{\left( {1 - \varepsilon } \right)}^2}{x_1}{x_2}}}{{{{\left( {a + \left( {1 - \varepsilon } \right){x_1}} \right)}^3}}}$

(1) If ${\text{trace}}\left( {J\left( {{E^*}} \right)} \right) \leqslant 0$, then we have to consider three cases:

Case ⅰ) For ${\text{trace}}\left( {J\left( {{E^*}} \right)} \right) = 0$, we will get $\left| {\arg \left( {{\lambda _i}} \right)} \right| = \frac{\pi }{2}, \, \, \, i = 1, \, \, 2$ and ${E^*}$ is locally asymptotically stable.

Case ⅱ) For ${\text{trace}}\left( {J\left( {{E^*}} \right)} \right) < 0$ and ${\text{trac}}{{\text{e}}^{\text{2}}}\left( {J\left( {{E^*}} \right)} \right) - 4\det \left( {J\left( {{E^*}} \right)} \right) \geqslant 0$, then both latent roots are negative and $\left| {\arg \left( {{\lambda _i}} \right)} \right| = \pi , \, i = 1, 2$. Thus, ${E^*}$ is locally asymptotically stable.

Case ⅲ) For ${\text{trace}}\left( {J\left( {{E^*}} \right)} \right) < 0$ and ${\text{trac}}{{\text{e}}^{\text{2}}}\left( {J\left( {{E^*}} \right)} \right) - 4\det \left( {J\left( {{E^*}} \right)} \right) < 0$, then both latent roots are complex conjugates with non-positive real part and $\left| {\arg \left( {{\lambda _i}} \right)} \right| > \frac{{\xi \pi }}{2}, \, i = 1, 2$. Therefore, ${E^*}$ is locally asymptotically stable for any $\xi$.

(2) If ${\text{trace}}\left( {J\left( {{E^*}} \right)} \right) > 0$, ${\text{trac}}{{\text{e}}^2}\left( {J\left( {{E^*}} \right)} \right) - 4\det \left( {J\left( {{E^*}} \right)} \right) < 0$ and ${\text{sqrt}}\left( {\left| {{\text{trac}}{{\text{e}}^{\text{2}}}\left( {J\left( {{E^*}} \right)} \right) - 4\det \left( {J\left( {{E^*}} \right)} \right)} \right|} \right) > {\text{trace}}\left( {J\left( {{E^*}} \right)} \right)\tan \left( {\frac{{\xi \pi }}{2}} \right)$, then we obtain ${Im} g\left( {{\lambda _1}} \right) = - {Im} g\left( {{\lambda _2}} \right) = 4\det \left( {J\left( {{E^*}} \right)} \right) - {\text{trac}}{{\text{e}}^{\text{2}}}\left( {J\left( {{E^*}} \right)} \right)$ and ${Re} \left( {{\lambda _i}} \right) = {\text{trace}}\left( {J\left( {{E^*}} \right)} \right) > 0, \, i = 1, 2$. From these we can conclude that $\left| {{Im} g\left( {{\lambda _i}} \right)} \right| > \tan \left( {\frac{{\xi \pi }}{2}} \right)\left( {{Re} \left( {{\lambda _i}} \right)} \right), \, i = 1, 2$. Consequently, ${E^*}$ is locally asymptotically stable if $\left| {\arg \left( {{\lambda _i}} \right)} \right| > \frac{{\xi \pi }}{2}, \, i = 1, 2$.

3.3.   Global stability

We next proceed to examine ${E^*}$'s global stability. To achieve this objective, Lemma 3.1 of ^[42] and the generalised LaSalle's invariance principle by Huo et al. ^[43] are employed.

Now, consider a function $V\left( {{x_1}, \, {x_2}} \right) = \left( {{x_1} - x_1^* - x_1^*\log \frac{{{x_1}}}{{x_1^*}}} \right) + {a_1}\left( {{x_2} - x_2^* - x_2^*\log \frac{{{x_2}}}{{x_2^*}}} \right)$, where ${a_1} > 0$ is a constant. Lemma 3.1 of above reference and the $\xi$-order non-integer derivative of $V\left( {{x_1}, \, {x_2}} \right)$ along the solution of (2) give us

the following can be derived by assuming $a{a_1}{b_1} = \left( {ab + b\left( {1 - \varepsilon } \right)x_1^*} \right)$,

Therefore, the positive equilibrium point ${E^*}$ is globally asymptotically stable.

3.4.   Existence of Hopf bifurcation

The transition from stable spiral to unstable behaviour at the equilibrium point causes Hopf bifurcation to occur when the bifurcation parameter reaches the critical value. It is obvious from Lemma 1 that the order of the derivative $\xi$ determines the stability of a dynamical system. Therefore, $\xi$ may be considered a bifurcation parameter, even if the non-integer order system does not possess the identical Hopf bifurcation condition as the integer-order framework.

Let us examine a system with fractional order

where $\xi \in \left[ {0, 1} \right), \, x \in {\Re ^2}$ and let ${E^*}$ is an equilibrium point. Based on the analysis of a non-integer order Hopf bifurcation in ^[44,45], it can be concluded that, under certain circumstances, system (10) experiences a Hopf bifurcation given that the value $m$ passes the threshold ${m^*}$ around the equilibria ${E^*}$.

ⅰ) The characteristic equation has two complex conjugate eigenvalues

ⅱ) ${\nu _{1, 2}}\left( {\xi , {m^*}} \right) = 0, \, {\text{where}}\, \, {\nu _{1, 2}} = \frac{{\xi \pi }}{2} - \left| {\arg \left( {{\lambda _i}\left( m \right)} \right)} \right|$,

ⅲ) ${\left. {\frac{{\partial {\nu _{1, 2}}}}{{\partial m}}} \right|_{m = {m^*}}} \ne 0$.

4.   Numerical simulations

Theoretical results have been validated through numerical simulations employing the predictor-corrector technique on non-integer differential equations. In order to resolve fractional order nonlinear differential equations, numerical methods are employed. This research explores the impact of altering several factors on the non-integer order predator-prey dynamical framework (2).

Figure 1 illustrates the phase portraits of the framework (2) about ${E_2}$ for parameter values as $r = 2.5, k = 0.899, \beta = 0.945, \varepsilon = 0.058, \alpha = 0.559, h = 0.291,$ $d = 0.059, \, \, {\beta _1} = 0.895$ and three different values of $\xi = 0.98, 0.92\, , {\text{and}}\, 0.85$ in order to study the influence of the derivative order ($\xi$) on the predator-prey dynamics. Based on these numbers, we know that system (2) goes from unstable to stable and that at $\xi = 0.92$, it shows evidence of Hopf bifurcation. Hence, we deduce that the dynamics of the system under consideration are significantly affected by the non-integer order derivative.

The influence of prey refuge rate $\varepsilon$ has been visually explored for $\xi = 1, 0.92$, with the values of $\varepsilon = 0.08, 0.6$, and 0.8, as well as the set of parameter values $r = 2.5, k = 0.899,$ $\beta = 0.945, \varepsilon = 0.058, h = 0.291, d = 0.059, \, {\beta _1} = 0.895$ shown in Figures 2 and 3. The equilibrium point ${E_2}$, as seen in Figure 2, is not stable at $\xi = 1$ and $\varepsilon = 0.08$ but, it should be emphasised, becomes stable at $\xi = 0.92$ and $\varepsilon = 0.08$. This leads us to the conclusion that the existence of prey refuge affects the prey-predator system, and fractional order $\xi$ contributes to species persistence by stabilising the system.

Considering the identical selection of parameter values as illustrated in Figure 1, the influence of prey harvesting rate $h$ has been visually explored for $\xi = 0.98$ and 1. Figures 4 and 5 illustrate that, when $\xi = 0.98$, the system (2) is stable at $h = 0.5$ then not stable for $h = 0.1$. In connection with the parameter $h$, bifurcation diagrams are provided in Figure 6.

5.   Conclusions

In the present study, we explored a fractional-order Holling type Ⅱ predator-prey framework that includes selective nonlinear harvesting and refuge in prey. As soon as the density of the harvested population reaches a certain limit, the harvesting grows smoothly to that value. The basic objective of this research is to understand how fractional order, refuge, and nonlinear harvesting affect the model's dynamic behaviour. The existence, uniqueness, non-negativity, and boundedness of the system's solution were defined first (2). Next, the local resilience of each feasible equilibrium point in the predator-prey framework was investigated using Matignon's criterion. Through the subsequent development of a suitable Lyapunov function, a global stability analysis of system (2) has been explored. Next, using $\xi$ and $h$ as parameters, the presence of Hopf bifurcation was investigated. Furthermore, the theoretical conclusions were numerically confirmed by employing the Adams-Bashforth PECE method (time-domain).

The system is stabilised by the prey refuge rate, according to numerical simulations. The framework exhibits stability in both integer-order as well as fractional orders when the refuge rates of prey populations are reduced. When the rate of prey refuge increases at a certain level of predation and harvesting, the system becomes unstable. Furthermore, it has been shown that nonlinear harvesting plays a vital part in maintaining system stability. The system shows stable behaviour at larger values of prey harvesting, whereas fractional order $\xi = 0.98$ shows stability at lower prey population harvesting. Therefore, given a constant value of the prey refuge amount, the system shows stable behaviour when there is a significant quantity of prey harvesting. Analysing the influence of spatial organisation on the overall behaviour of the system in the presence of a nonlinear harvesting term might provide valuable insights. Although this model can include spatial structure inpattern-forming Turing dynamics, it may be challenging to conduct comparable analyses and simulations for nonlinear harvesting techniques in this particular scenario. Therefore, more investigation into this matter is necessary. Also, the incorporation of prey-refuging strategy switching based on the prey-to-predator density ratio in the proposed model offers a more realistic ecological representation. These results underscore the importance of threshold-based refuge seeking for preserving species and maintaining ecosystem stability.

Use of AI tools declaration

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

Acknowledgments

The authors extend their appreciation to Prince Sattam bin Abdulaziz University for funding this research work through the project number (PSAU/2024/01/78916).

Conflict of interest

The authors state that none of the work presented in this study may have been influenced by any known conflicting financial interests or personal relationships.