A Rosenzweig–MacArthur Model with Continuous Threshold Harvesting in Predator Involving Fractional Derivatives with Power Law and Mittag–Leffler Kernel

Abstract

The harvesting management is developed to protect the biological resources from over-exploitation such as harvesting and trapping. In this article, we consider a predator–prey interaction that follows the fractional-order Rosenzweig–MacArthur model where the predator is harvested obeying a threshold harvesting policy (THP). The THP is applied to maintain the existence of the population in the prey–predator mechanism. We first consider the Rosenzweig–MacArthur model using the Caputo fractional-order derivative (that is, the operator with the power-law kernel) and perform some dynamical analysis such as the existence and uniqueness, non-negativity, boundedness, local stability, global stability, and the existence of Hopf bifurcation. We then reconsider the same model involving the Atangana–Baleanu fractional derivative with the Mittag–Leffler kernel in the Caputo sense (ABC). The existence and uniqueness of the solution of the model with ABC operator are established. We also explore the dynamics of the model with both fractional derivative operators numerically and confirm the theoretical findings. In particular, it is shown that models with both Caputo operator and ABC operator undergo a Hopf bifurcation that can be controlled by the conversion rate of consumed prey into the predator birth rate or by the order of fractional derivative. However, the bifurcation point of the model with the Caputo operator is different from that of the model with the ABC operator.

1. Introduction

More than 50 years after the model has been proposed, the Rosenzweig–MacArthur predator–prey model [1] has been consistently developed by many scholars to approach the real world phenomena with more realistic mathematical models. The commonsensical modified Rosenzweig–MacArthur models are accomplishable by considering the biological perspectives of ecosystem conditions, for instance the stage structure [2,3], the refuge effect [4,5,6,7,8], the fear effect [9], the Allee effect [10,11], the intraspecific competition [12,13], the cannibalism [14], the infectious disease [15,16,17], and so forth.

On the other hand, the modeling also contemplates the optimal management of bioeconomic resources as in fishery and pest management. Some researchers put an intervention into the predator–prey model such as the harvesting to one or more population [8,18,19,20,21,22]. To protect the population from over-exploitation during the harvesting, some management techniques have been established. One of the famous technique is a continuous threshold harvesting policy (THP) (see [23,24,25,26,27,28,29]). THP is regulated as follows: when the population density above the threshold level, harvesting is permitted; when the population density drops below the threshold level, harvesting is prohibited.

In 2013, Lv et al. [26] proposed the following Rosenzweig–MacArthur model with THP in predator

$$\begin{array}{cc}\hfill \frac{dx}{dt}=& rx\left(1-\frac{x}{K}\right)-\frac{mxy}{a+x},\hfill \\ \hfill \frac{dy}{dt}=& \frac{nxy}{a+x}-dy-H\left(y\right),\hfill \end{array}$$

(1)

where

$$H\left(y\right)=\left\{\begin{array}{ccc}0& ,& \mathrm{if}y<T,\hfill \\ {\displaystyle \frac{h(y-T)}{c+(y-T)}}& ,& \mathrm{if}y\ge T.\hfill \end{array}\right.$$

We describe the biological interpretation of variables and parameters of model (1) in

Table 1. Description of variables and parameters of the model (1).

Variables and Parameters	Description
$x\left(t\right)$	The density of prey
$y\left(t\right)$	The density of predator
r	The intrinsic growth rate of prey
K	The environmental carrying capacity of prey
m	The maximum uptake rate for prey
n	The conversion rate of consumed prey into predator birth
a	The environment protection for prey
d	The natural death rate of predator
h	The harvesting rate
c	The half saturation constant for harvesting
T	The threshold level of harvesting

. Model (1) represents an interaction between two populations with a prey–predator relationship, where THP is only applied for the predator to preserving its populations. Some appealing examples of the ecological model (1) are given by the interaction between Sycanus sp. and Setothosea asigna and between Rhinocoris sp. and Spodoptera litura. Shepard [30] reported that Sycanus sp. and Rhinocoris sp. are the natural predators of the pests such as Setothosea asigna and Spodoptera litura which exist in agricultural lands and plantations. The worrying problem is: How if the density of insects such as Sycanus sp. and Rhinocoris sp. uncontrolled? One solution is applying THP as in model (1).

Lv et al. [26] successfully explored the dynamics of the model (1) including the local stability and the existence of various phenomena (saddle-node, Hopf, cusp, and Bogdanov–Taken bifurcations). Despite their success works, the model with the first-order derivative is limited to its capability of involving all previous conditions to the growth rates of both predator and prey. The growth rates of both populations in the model (1) depend only on the current state. Biologically, the growth rates must be dependent on all of the previous conditions which are known as the memory effects. To account for such memory effects, some researchers proposed to apply the fractional-order derivative instead of the first-order derivative when expressing the growth rate of the population. The fractional-order models are naturally related to systems with memory which exists in most biological models [7,31]. The fractional-order models are also well-liked due to their capability in providing an exact description of different nonlinear phenomena [32]. In recent years, the development of fractional-order models grows rapidly and becomes popular in studying the dynamical behavior of predator–prey interaction [17,33,34,35,36,37,38]. It has been shown that the order of the fractional derivate significantly affects the dynamical behavior of the models, which is in contrast to the first-order derivative models that depend only on the values of parameters.

In this paper, we modify the model of Lv et al. [26] by applying the fractional-order derivative to the left-hand sides of the first-order differential Equation (1). We use two types of fractional-order derivatives, namely the Caputo operator (that is, the operator with the power law kernel) [39] and Atangana–Baleanu operator [40]. The basic difference between these two operators lies on their kernel. Atangana–Baleanu operator has a non-singular and non-local (that is, Mittag–Leffler) kernel while the Caputo operator does not [41,42]. From the biological meanings, a model with Atangana–Baleanu operator gives better results in applying memory effects [43,44,45]. Nevertheless, the Caputo operator has more complex analytical tools in investigating the dynamics of the model such as the local stability [46,47,48,49,50], the global stability [51,52], bifurcation theory [52,53,54], and so on. By considering the deficiencies and advantages, the model with Caputo and Atangana–Baleanu operator are employed in our work. Based on our literature review, the dynamics of the model (1) with Caputo and Atangana–Baleanu operator have never been studied. For this reason, we are interested in investigating the dynamical behavior of model (1) using both Caputo and Atangana–Baleanu fractional-order operators.

If the first-order derivatives $\frac{d}{dt}$ at the left hand sides of model (1) are replaced by the fractional-order derivatives ${\mathcal{D}}_t^{\alpha }$, then we obtain

$$\begin{array}{cc}\hfill {\mathcal{D}}_t^{\alpha }x=& \widehat{r}x\left(1-\frac{x}{K}\right)-\frac{\widehat{m}xy}{a+x},\hfill \\ \hfill {\mathcal{D}}_t^{\alpha }y=& \frac{\widehat{n}xy}{a+x}-\widehat{d}y-H\left(y\right).\hfill \end{array}$$

(2)

Note that the left hand sides of model (2) have the dimension of (time)${}^{-\alpha }$, while the parameters at the right hand sides of model (2) such as $\widehat{r},\widehat{m},\widehat{n},\widehat{d}$, and $\widehat{h}$ have the dimension of (time)${}^{-1}$; this shows the inconsistency of physical dimensions in the model (2) (see [55,56]). To overcome such inconsistency, we rescale the parameters in the model (2) to get the following model

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{D}}_t^{\alpha }x=& {\widehat{r}}^{\alpha }x\left(1-\frac{x}{K}\right)-\frac{{\widehat{m}}^{\alpha }xy}{a+x},\hfill \\ \hfill {\mathcal{D}}_t^{\alpha }y=& \frac{{\widehat{n}}^{\alpha }xy}{a+x}-{\widehat{d}}^{\alpha }y-H\left(y\right),\hfill \end{array}\end{array}$$

(3)

where

$$H\left(y\right)=\left\{\begin{array}{ccc}0& ,& \mathrm{if}y<T,\hfill \\ {\displaystyle \frac{{\widehat{h}}^{\alpha }(y-T)}{c+(y-T)}}& ,& \mathrm{if}y\ge T.\hfill \end{array}\right.$$

By applying new parameters $r={\widehat{r}}^{\alpha },m={\widehat{m}}^{\alpha },n={\widehat{n}}^{\alpha },d={\widehat{d}}^{\alpha }$, and $h={\widehat{h}}^{\alpha }$, we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{D}}_t^{\alpha }x=& rx\left(1-\frac{x}{K}\right)-\frac{mxy}{a+x},\hfill \\ \hfill {\mathcal{D}}_t^{\alpha }y=& \frac{xy}{a+x}-dy-H\left(y\right),\hfill \end{array}\end{array}$$

(4)

where

$$H\left(y\right)=\left\{\begin{array}{ccc}0& ,& \mathrm{if}y<T,\hfill \\ {\displaystyle \frac{h(y-T)}{c+(y-T)}}& ,& \mathrm{if}y\ge T.\hfill \end{array}\right.$$

This paper is organized as follows. In Section 2, we investigate the dynamics of model (4) with the Caputo operator. We identify the existence, uniqueness, non-negativity, and boundedness of solutions. Furthermore, we explore the dynamics of the model by examining the existence of the equilibrium points, their local and global stability, and the existence of Hopf bifurcation. In Section 3, the existence and uniqueness of solutions of the model with Atangana–Baleanu operator are verified. In Section 4, we explore the dynamics of the model using both operators numerically. We demonstrate numerically the stability of the equilibrium point, and the occurrence of forward and Hopf bifurcations. We end our works with a brief conclusion in Section 5.

2. The Caputo Fractional-Order Rosenzweig–MacArthur Model with THP in Predator

2.1. Model Formulation

The operator of Caputo fractional-order derivative is defined as follows

Definition 1.

[48] Let $\alpha \in (0,1]$, $f\in C^n([0,+\infty ),\mathbb{R})$, and $\Gamma (·)$ is the Gamma function. The Caputo fractional derivative of order-α is defined by

$${}^C{\mathcal{D}}_t^{\alpha }f\left(t\right)=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t{(t-s)}^{-\alpha }{f}^{\prime }\left(s\right)ds,t\ge 0.$$

(5)

The kernel of Caputo operator is known as the power law kernel. By applying Definition 1 to model (4), we get the Caputo fractional order Rosenzweig–MacArthur model with THP in predator

$$\begin{array}{cc}\hfill {}^C{\mathcal{D}}_t^{\alpha }x=& rx\left(1-\frac{x}{K}\right)-\frac{mxy}{a+x}\equiv F_1,\hfill \\ \hfill {}^C{\mathcal{D}}_t^{\alpha }y=& \frac{nxy}{a+x}-dy-H\left(y\right)\equiv F_2.\hfill \end{array}$$

(6)

2.2. Existence and Uniqueness

In this part, we study the existence and uniqueness of model (6).

Lemma 1.

[57] Consider a Caputo fractional-order system

$${}^C{\mathcal{D}}_t^{\alpha }x\left(t\right)=f(t,x\left(t\right)),t>0,x\left(0\right)\ge 0,\alpha \in (0,1],$$

(7)

where $f:(0,\infty )\times \Theta \to {\mathbb{R}}^n,\Theta \subseteq {\mathbb{R}}^n$. A unique solution of Equation (7) on $(0,\infty )\times \Theta$ exists if $f(t,x(t\left)\right)$ satisfies the locally Lipschitz condition with respect to x.

Since the right hand-side of model (6) is a piecewise function which is switched when the number of predators passes through the threshold level, we divide the existence and uniqueness of the solution into two cases, namely $y\ge T$ and $y<T$. We start from $y\ge T$. Consider the region $\Theta \times [0,T_{+}]$ where $\Theta :=\left\{(x,y)\in {\mathbb{R}}^2:max\left(\left|x\right|,\left|y\right|\right)\le \gamma ,y\ge T\right\}$, $T_{+}<+\infty$, and a mapping $F(\Lambda )=(F_1(\Lambda ),F_2(\Lambda ))$. For any $\Lambda =\left(x_{,}y\right)\in \Theta$ and $\overline{\Lambda }=\left(\overline{x},\overline{y}\right)\in \Theta$, we obtain

$$\begin{array}{cc}\hfill ∥F(\Lambda )-F\left(\overline{\Lambda }\right)∥=& \left|F_1(\Lambda )-F_1\left(\overline{\Lambda }\right)\right|+\left|F_2(\Lambda )-F_2\left(\overline{\Lambda }\right)\right|\hfill \\ \hfill =& \left|rx\left(1-\frac{x}{K}\right)-\frac{mxy}{a+x}-\left(r\overline{x}\left(1-\frac{\overline{x}}{K}\right)-\frac{m\overline{x}\overline{y}}{a+\overline{x}}\right)\right|+\hfill \\ \hfill & \left|\frac{nxy}{a+x}-dy-H\left(y\right)-\left(\frac{n\overline{x}\overline{y}}{a+\overline{x}}-d\overline{y}-H\left(\overline{y}\right)\right)\right|\hfill \\ \hfill =& \left|r(x-\overline{x})-\frac{r}{K}\left(x^2-{\overline{x}}^2\right)-m\left(\frac{xy}{a+x}-\frac{\overline{x}\overline{y}}{a+\overline{x}}\right)\right|+\hfill \\ \hfill & \left|n\left(\frac{xy}{a+x}-\frac{\overline{x}\overline{y}}{a+\overline{x}}\right)-d\left(y-\overline{y}\right)-\left(H\left(y\right)-H\left(\overline{y}\right)\right)\right|\hfill \\ \hfill =& r\left|x-\overline{x}\right|+\frac{r}{K}\left|x+\overline{x}\right|\left|x-\overline{x}\right|+(m+n)\left|\frac{xy}{a+x}-\frac{\overline{x}\overline{y}}{a+\overline{x}}\right|+\hfill \\ \hfill & d\left|y-\overline{y}\right|+\left|\frac{h(y-T)}{c+(y-T)}-\frac{h(\overline{y}-T)}{c+(\overline{y}-T)}\right|\hfill \\ \hfill =& r\left|x-\overline{x}\right|+\frac{r}{K}\left|x+\overline{x}\right|\left|x-\overline{x}\right|+(m+n)\left|\frac{ay(x-\overline{x})+(a\overline{x}+\overline{x}x)(y-\overline{y})}{(a+x)(a+\overline{x})}\right|+\hfill \\ \hfill & d\left|y-\overline{y}\right|+\left|\frac{ch(y-\overline{y})}{(c+y-T)(c+\overline{y}-T)}\right|\hfill \\ \hfill \le & r\left|x-\overline{x}\right|+\frac{2\gamma r}{K}\left|x-\overline{x}\right|+\frac{m+n}{a^2}\left|ay(x-\overline{x})+(a\overline{x}+\overline{x}x)(y-\overline{y})\right|+\hfill \\ \hfill & d\left|y-\overline{y}\right|+\frac{h}{c}\left|y-\overline{y}\right|\hfill \\ \hfill \le & r\left|x-\overline{x}\right|+\frac{2\gamma r}{K}\left|x-\overline{x}\right|+\frac{(m+n)\gamma }{a}\left|x-\overline{x}\right|+\hfill \\ & \frac{(m+n)(a\gamma +{\gamma }^2)}{a^2}\left|y-\overline{y}\right|+d\left|y-\overline{y}\right|+\frac{h}{c}\left|y-\overline{y}\right|\hfill \\ \hfill =& \left(r+\frac{2\gamma r}{K}+\frac{(m+n)\gamma }{a}\right)\left|x-\overline{x}\right|+\left(\frac{(m+n)(a\gamma +{\gamma }^2)}{a^2}+d+\frac{h}{c}\right)\left|y-\overline{y}\right|\hfill \\ \hfill \le & M_1∥\Lambda -\overline{\Lambda }∥,\hfill \end{array}$$

where $M_1=max\left\{r+{\displaystyle \frac{2\gamma r}{K}}+{\displaystyle \frac{(m+n)\gamma }{a}},{\displaystyle \frac{(m+n)(a\gamma +{\gamma }^2)}{a^2}}+d+{\displaystyle \frac{h}{c}}\right\}$. Hence, $F(\Lambda )$ satisfies the Lipschitz condition for $y\ge T$. By using similar approaches, when $y<T$, it is easy to check that $∥F(\Lambda )-F\left(\overline{\Lambda }\right)∥\le M_2∥\Lambda -\overline{\Lambda }∥$, where $M_2=max\left\{r+{\displaystyle \frac{2\gamma r}{K}}+{\displaystyle \frac{(m+n)\gamma }{a}},{\displaystyle \frac{(m+n)(a\gamma +{\gamma }^2)}{a^2}}+d\right\}$ and hence the Lipschitz condition is also satisfied. According to Lemma 1, model (6) with non-negative initial condition has a unique solution $\Lambda \left(t\right)=\left(x\right(t),y(t\left)\right)\in \Theta$. Thus, we establish the following theorem.

Theorem 1.

For each non-negative initial condition $(x_0,y_0)\in \Theta$, there exists a unique solution $\left(x\right(t),y(t\left)\right)\in \Theta$ of model (6), which is defined for all $t\ge 0$.

2.3. Non-Negativity and Boundedness

The solution of model (6) is required to be nonnegative and bounded to establish a biologically well-behaved model. To determine the non-negativity and boundedness of the solution of model (6), the following lemmas are needed.

Lemma 2.

[58] Let $0<\alpha \le 1$. Suppose that $f\left(t\right)\in C[a,b]$ and ${}^C{\mathcal{D}}_t^{\alpha }f\left(t\right)\in C[a,b]$. If ${}^C{\mathcal{D}}_t^{\alpha }f\left(t\right)\ge 0,\forall t\in (a,b)$, then $f\left(t\right)$ is a non-decreasing function for each $t\in [a,b]$. If ${}^C{\mathcal{D}}_t^{\alpha }f\left(t\right)\le 0,\forall t\in (a,b)$, then $f\left(t\right)$ is a non-increasing function for each $t\in [a,b]$.

Lemma 3.

(Standard comparison theorem for Caputo fractional-order derivative [31]). Let $x\left(t\right)\in C\left([0,+\infty )\right)$. If $x\left(t\right)$ satisfies ${}^C{\mathcal{D}}_t^{\alpha }x\left(t\right)+\lambda x\left(t\right)\le \mu ,x\left(0\right)=x_0$, where $\alpha \in (0,1],(\lambda ,\mu )\in {\mathbb{R}}^2$ and $\lambda \ne 0$, then $x\left(t\right)\le \left(x_0-{\displaystyle \frac{\mu }{\lambda }}\right)E_{\alpha }[-\lambda t^{\alpha }]+{\displaystyle \frac{\mu }{\lambda }}$.

In the following theorem, we prove the non-negativity and boundedness of solutions using the above lemmas.

Theorem 2.

All solutions of model (6) with non-negative initial conditions are non-negative and uniformly bounded.

Proof. 

We start by proving that, if the initial condition is non-negative, then $x\left(t\right)\ge 0$ for all $t\to \infty$. Suppose that it is incorrect; then, we can find $t_1>0$ such that

$$\left\{\begin{array}{ccc}x\left(t\right)\hfill & >& 0,0\le t<t_1,\hfill \\ x\left(t_1\right)\hfill & =& 0,\hfill \\ x\left(t_1^{+}\right)\hfill & <& 0\hfill \end{array}\right.$$

(8)

By employing (8) and the first equation of model (6), we obtain

$${\left.{}^C{\mathcal{D}}_t^{\alpha }x\left(t_1\right)\right|}_{x\left(t_1\right)=0}=0.$$

(9)

Based on Lemma 2, we have $x\left(t_1^{+}\right)=0$, which contradicts the fact that $x\left(t_1^{+}\right)<0$. Thus, $x\left(t\right)\ge 0$ for all $t\ge 0$. Using a similar procedure, we conclude $y\left(t\right)\ge 0$ for all $t>0$.

Now, we adopt the similar manner as in [34] to show the boundedness of solutions. By setting up a function $\mathcal{V}\left(t\right)=x+{\displaystyle \frac{my}{n}}$, we get

$$\begin{array}{cc}\hfill {}^C{\mathcal{D}}_t^{\alpha }\mathcal{V}\left(t\right)+d\mathcal{V}\left(t\right)=& {}^C{\mathcal{D}}_t^{\alpha }x+\frac{m}{n}{}^C{\mathcal{D}}_t^{\alpha }y+dx+\frac{dmy}{n}\hfill \\ \hfill =& rx\left(1-\frac{x}{K}\right)-\frac{mxy}{a+x}+\frac{m}{n}\left(\frac{nxy}{a+x}-dy-H\left(y\right)\right)+dx+\frac{dmy}{n}\hfill \\ \hfill =& (d+r)x-\frac{r{x}^2}{K}-\frac{mH\left(y\right)}{n}\hfill \\ \hfill =& -\frac{r}{K}{\left(x-\frac{(d+r)K}{2r}\right)}^2+\frac{{(d+r)}^{2}K}{4r}-\frac{mH\left(y\right)}{n}\hfill \\ \hfill \le & \frac{{(d+r)}^{2}K}{4r}.\hfill \end{array}$$

According to the standard comparison theorem for Caputo fractional-order derivative in Lemma 3, we achieve the following inequality

$$\mathcal{V}\left(t\right)\le \left(\mathcal{V}\left(0\right)-\frac{{(d+r)}^{2}K}{4r}\right)E_{\alpha }\left[-d{\left(t\right)}^{\alpha }\right]+\frac{{(d+r)}^{2}K}{4r},$$

where $E_{\alpha }$ is the one-parameter Mittag–Leffler function. Since $E_{\alpha }\left[-d{\left(t\right)}^{\alpha }\right]\to 0$ as $t\to 0$, we acquire $\mathcal{V}\left(t\right)\to {\displaystyle \frac{{(d+r)}^{2}K}{4r}}$ for $t\to \infty$. Hence, with non-negative initial condition, all solutions are restricted to the region ${\Theta }_M$ where

$${\Theta }_M:=\left\{(x,y)\in {\mathbb{R}}_{+}^2:x+\frac{my}{n}\le {\displaystyle \frac{{(d+r)}^{2}K}{4r}}+\epsilon ,\epsilon >0\right\}.$$

Consequently, all solutions of model (6) are uniformly bounded. □

2.4. The Existence of Equilibrium Point

The first commonplace technique in studying the dynamical behavior of a fractional-order system is investigating the existence of the equilibrium point. Due to the biological nature, we give the following definition.

Definition 2.

Consider a Caputo fractional-order system

$${}^C{\mathcal{D}}_t^{\alpha }\overrightarrow{x}=\overrightarrow{f}(\overrightarrow{x});\overrightarrow{x}\left(0\right)={\overrightarrow{x}}_0;\alpha \in (0,1].$$

(10)

A point ${\overrightarrow{x}}^{*}$ is called an equilibrium point of Equation (10) if it satisfies $\overrightarrow{f}\left({\overrightarrow{x}}^{*}\right)=0$. Particularly, it is called the biological equilibrium point of Equation (10) if it satisfies ${\overrightarrow{x}}^{*}\ge 0$.

Based on Definition 2, the equilibrium point of model (6) is obtained by solving

$$\begin{array}{cc}\hfill \left[r-\frac{rx}{K}-\frac{my}{a+x}\right]x=& 0,\hfill \\ \hfill \frac{nxy}{a+x}-dy-H\left(y\right)=& 0.\hfill \end{array}$$

(11)

Thus, we get four feasible biological equilibrium points as follows.

(i)

The equilibrium points when $y<T$ are

(i.a)

the origin point $E_0=(0,0)$ which always exists;

(i.b)

the predator extinction point $E_1=(K,0)$ which always exists; and

(i.c)

the first co-existence point $\widehat{E}=\left(\widehat{x},{\displaystyle \frac{\left(K-\widehat{x}\right)(a+\widehat{x})r}{mK}}\right)$, with $\widehat{x}={\displaystyle \frac{ad}{n-d}}$, which exists if $n>{\displaystyle \frac{ad}{K}}+d$ and $\left(K-\widehat{x}\right)(a+\widehat{x})<{\displaystyle \frac{TmK}{r}}$.

(ii)

The equilibrium point when $y\ge T$ is the second co-existence point $E^{*}=\left(x^{*},y^{*}\right)$ where $y^{*}={\displaystyle \frac{\left(K-x^{*}\right)(a+x^{*})r}{mK}}$ and $x^{*}$ is the positive roots of polynomial ${\beta }_1{x}^4+{\beta }_2{x}^3+{\beta }_3{x}^2+{\beta }_{4}x+{\beta }_5=0$ where

$$\begin{array}{cc}\hfill {\beta }_1=& (n-d)r^2,\hfill \\ \hfill {\beta }_2=& \left[(an+2dK)-2(ad+nK)\right]r^2,\hfill \\ \hfill {\beta }_3=& (nrK+4adr+cdm+mnT+hm)rK\hfill \\ \hfill & -((drK+2anr+cmn+dmT)K+a^{2}dr)r,\hfill \\ \hfill {\beta }_4=& \left(\right(anr+cmn+dmT)K+(2adr+hm+cdm\left)a\right)rK\hfill \\ \hfill & -\left(\right(2adr+hm+cdm+mnT)K+admT)rK,\hfill \\ \hfill {\beta }_5=& \left[(adr+hm)mT-(adr+hm+cdm)ar\right]K^2.\hfill \end{array}$$

$E^{*}$ exists if $0<x^{*}<K$ and $\left(K-x^{*}\right)(a+x^{*})\ge {\displaystyle \frac{TmK}{r}}$.

2.5. Local Asymptotic Stability

In this part, we discuss the local stability of $E_0$, $E_1$, $\widehat{E}$, and $E^{*}$. For this aim, we need the following theorem.

Theorem 3.

(Matignon condition [48,59]) The equilibrium point ${\overrightarrow{x}}^{*}$ of system (10) is locally asymptotically stable if all eigenvalues ${\lambda }_j$ of the Jacobian matrix $J=\partial \overrightarrow{f}/\partial \overrightarrow{x}$ evaluated at ${\overrightarrow{x}}^{*}$ satisfy $|arg({\lambda }_j\left)\right|>\alpha \pi /2$. If there exists at least one eigenvalue satisfying $|arg({\lambda }_k\left)\right|>\alpha \pi /2$ and $|arg({\lambda }_l\left)\right|<\alpha \pi /2$, $k\ne l$, then ${\overrightarrow{x}}^{*}$ is a saddle-point.

Now, we present Theorems 4–7 to show the local stability properties of $E_0$, $E_1$, $\widehat{E}$, and $E^{*}$.

Theorem 4.

The origin point $E_0=(0,0)$ is always a saddle point.

Proof. 

When $E_0=(0,0)$, the model (6) has Jacobian matrix $J\left(E_0\right)=\left[\begin{array}{cc}r& 0\\ 0& -d\end{array}\right]$, where its eigenvalues are ${\lambda }_1=r>0$ and ${\lambda }_2=-d<0$. It is clear that $\left|arg\left({\lambda }_1\right)\right|=0<\alpha \pi /2$ and $\left|arg\left({\lambda }_2\right)\right|=\pi >\alpha \pi /2$. Therefore, based on Theorem 3, $E_0$ is always a saddle point. □

Theorem 5.

The predator extinction point $E_1=(K,0)$ is locally asymptotically stable if $n<{\displaystyle \frac{ad}{K}}+d$. Otherwise, if $n>{\displaystyle \frac{ad}{K}}+d$, then $E_1=(K,0)$ is a saddle point.

Proof. 

The Jacobian matrix of model (6) evaluated at $E_1$ is $J\left(E_1\right)=\left[\begin{array}{cc}-r& -\frac{mK}{a+K}\\ 0& \frac{nK}{a+K}-d\end{array}\right]$. The eigenvalues of $J\left(E_1\right)$ are ${\lambda }_1=-r<0$ and ${\lambda }_2={\displaystyle \frac{nK}{a+K}}-d$. Clearly, $\left|arg\left({\lambda }_1\right)\right|=\pi >\alpha \pi /2$ and $\left|arg\left({\lambda }_2\right)\right|=\pi >\alpha \pi /2$ if $n<{\displaystyle \frac{ad}{K}}+d$ and $\left|arg\left({\lambda }_2\right)\right|=0<\alpha \pi /2$ if $n>{\displaystyle \frac{ad}{K}}+d$. Hence, we have the theorem. □

Remark 1.

It is noted that the existence condition for the first co-existence point $\widehat{E}$ contradicts the stability condition of $E_1$. Consequently, if $E_1$ is locally asymptotically stable, then $\widehat{E}$ does not exist. This condition also indicates the existence of forward bifurcation, which is confirmed numerically in the next section.

Theorem 6.

Let $\Delta ={\displaystyle \frac{4\left(K-\widehat{x}\right)anr\widehat{x}}{{(a+\widehat{x})}^{2}K}}-{\left({\displaystyle \frac{K-\widehat{x}}{a+\widehat{x}}}-1\right)}^2{\displaystyle \frac{r^2{\widehat{x}}^2}{K^2}}$ and $\widehat{\alpha }={\displaystyle \frac{2}{\pi }}{tan}^{-1}\left({\displaystyle \frac{\sqrt{\Delta }(a+\widehat{x})K}{\left(K-a-2\widehat{x}\right)r\widehat{x}}}\right)$. Suppose that one of the following statements is obeyed.

(i)

$\widehat{x}>{\displaystyle \frac{K-a}{2}}$; or

(ii)

$\widehat{x}<{\displaystyle \frac{K-a}{2}}$, $\Delta >0$, and $\alpha <\widehat{\alpha }$.

Then, the first co-existence point $\widehat{E}=\left(\widehat{x},{\displaystyle \frac{\left(K-\widehat{x}\right)(a+\widehat{x})r}{mK}}\right)$ is locally asymptotically stable.

Proof. 

We first observe that the Jacobian matrix of model (6) evaluated at $\widehat{E}$ is

$$J\left(\widehat{E}\right)=\left[\begin{array}{cc}\left({\displaystyle \frac{K-\widehat{x}}{a+\widehat{x}}}-1\right){\displaystyle \frac{r\widehat{x}}{K}}& -{\displaystyle \frac{m\widehat{x}}{a+\widehat{x}}}\\ {\displaystyle \frac{\left(K-\widehat{x}\right)anr}{(a+\widehat{x})mK}}& 0\end{array}\right].$$

(12)

The eigenvalues of the Jacobian matrix (12) are the solutions of the characteristic equation

$${\lambda }^2-\left(\frac{K-\widehat{x}}{a+\widehat{x}}-1\right)\frac{r\widehat{x}}{K}\lambda +\frac{\left(K-\widehat{x}\right)anr\widehat{x}}{{(a+\widehat{x})}^{2}K}=0,$$

namely

$$\begin{array}{cc}\hfill {\lambda }_1=& \left(\frac{K-\widehat{x}}{a+\widehat{x}}-1\right)\frac{r\widehat{x}}{2K}+\frac{i\sqrt{\Delta }}{2},\hfill \\ \hfill {\lambda }_2=& \left(\frac{K-\widehat{x}}{a+\widehat{x}}-1\right)\frac{r\widehat{x}}{2K}-\frac{i\sqrt{\Delta }}{2}.\hfill \end{array}$$

(13)

When $\widehat{x}>{\displaystyle \frac{K-a}{2}}$, the real parts of ${\lambda }_{1,2}$ are negative. Thus, the eigenvalues (13) always satisfy $\left|arg\left({\lambda }_{1,2}\right)\right|>\alpha \pi /2$ for any $\Delta$. If $\widehat{x}<{\displaystyle \frac{K-a}{2}}$ and $\Delta \le 0$, then ${\lambda }_1{\lambda }_2={\displaystyle \frac{\left(K-\widehat{x}\right)anr\widehat{x}}{{(a+\widehat{x})}^{2}K}}>0$ and ${\lambda }_1+{\lambda }_2=\left({\displaystyle \frac{K-\widehat{x}}{a+\widehat{x}}}-1\right){\displaystyle \frac{r\widehat{x}}{K}}>0$, meaning that $\left|arg\left({\lambda }_{1,2}\right)\right|=0<\alpha \pi /2$. If $\widehat{x}<{\displaystyle \frac{K-a}{2}}$ and $\Delta >0$, then $\left|arg\left({\lambda }_{1,2}\right)\right|>\alpha \pi /2$ for $\alpha <\widehat{\alpha }$. Hence, we prove the theorem. □

Theorem 7.

Let

$$\begin{array}{cc}\hfill \xi =& {\displaystyle \frac{({(y^{*}-T)}^2-cT)h}{y^{*}{(c+y^{*}-T)}^2}}+\left({\displaystyle \frac{m{y}^{*}}{{(a+x^{*})}^2}}-{\displaystyle \frac{r}{K}}\right)x^{*},\hfill \\ \hfill \theta =& \left({\displaystyle \frac{m{y}^{*}}{{(a+x^{*})}^2}}-{\displaystyle \frac{r}{K}}\right){\displaystyle \frac{({(y^{*}-T)}^2-cT)h}{y^{*}{(c+y^{*}-T)}^2}}{x}^{*}+\left(1-{\displaystyle \frac{x^{*}}{a+x}}\right){\displaystyle \frac{mn{x}^{*}{y}^{*}}{{(a+x^{*})}^2}}.\hfill \end{array}$$

If one of the following statements is satisfied, then the second co-existence point $E^{*}=(x^{*},y^{*})$ is locally asymptotically stable:

(i)

$\theta >0$ and $\xi <0$; or

(ii)

${\xi }^2<4\theta$, $\xi >0$, and $\alpha <{\alpha }^{*}$.

Proof. 

The Jacobian matrix of model (6) evaluated at $E^{*}$ is given by

$$J\left(E^{*}\right)=\left[\begin{array}{cc}\left({\displaystyle \frac{m{y}^{*}}{{(a+x^{*})}^2}}-{\displaystyle \frac{r}{K}}\right)x^{*}& -{\displaystyle \frac{m{x}^{*}}{a+x^{*}}}\\ \left(1-{\displaystyle \frac{x^{*}}{a+x^{*}}}\right){\displaystyle \frac{n{y}^{*}}{a+x^{*}}}& {\displaystyle \frac{({(y^{*}-T)}^2-cT)h}{y^{*}{(c+y^{*}-T)}^2}}\end{array}\right].$$

(14)

The eigenvalues of (14) are obtained by solving the characteristic equation ${\lambda }^2-\xi \lambda +\theta =0$. Thus, we have ${\lambda }_{1,2}={\displaystyle \frac{\xi }{2}}\pm {\displaystyle \frac{\sqrt{{\xi }^2-4\theta }}{2}}$. If $\theta >0$ and $\xi <0$, then $\left|arg\left({\lambda }_{1,2}\right)\right|>\alpha \pi /2$. If ${\xi }^2<4\theta$ and $\xi >0$, then $\left|arg\left({\lambda }_{1,2}\right)\right|>\alpha \pi /2$ for $\alpha <{\alpha }^{*}$. Using Theorem 3, the local stability of $E^{*}$ is completely proven. □

2.6. Global Asymptotic stability

To study the global stability of equilibrium points, we need the following lemmas.

Lemma 4.

[51] Let $x\left(t\right)\in C\left({\mathbb{R}}_{+}\right)$, $x^{*}\in {\mathbb{R}}_{+}$, and its Caputo fractional derivative of order-α exists for any $\alpha \in (0,1]$. Then, for any $t>0$, we have ${}^C{\mathcal{D}}_t^{\alpha }\left[x\left(t\right)-x^{*}-x^{*}ln{\displaystyle \frac{x\left(t\right)}{x^{*}}}\right]\le \left(1-{\displaystyle \frac{x^{*}}{x\left(t\right)}}\right){{}^C\mathcal{D}}_t^{\alpha }x\left(t\right)$.

Lemma 5.

(Generalized Lasalle Invariance Principle [52]). Suppose Ω is a bounded closed set and every solution of system

$${}^C{\mathcal{D}}_t^{\alpha }x\left(t\right)=f\left(x\left(t\right)\right),$$

(15)

which starts from a point in Ω remains in Ω for all time. Let $V\left(x\right):\Omega \to \mathbb{R}$ be a continuously differentiable function such that

$${}^C{\mathcal{D}}_t^{\alpha }{V|}_{\left(15\right)}\le 0.$$

Let $E:=\left\{x|{}^C{\mathcal{D}}_t^{\alpha }{V|}_{\left(15\right)}=0\right\}$ and M be the largest invariant set of E. Then, every solution $x\left(t\right)$ originating in Ω tends to M as $t\to \infty$.

Since model (6) is basically a piecewise fractional-order differential equations that depends on T, the analysis of the global stability is split into two regions defined by ${\Omega }_1:=\left\{(x,y):x\ge 0,y<T\right\}$ and ${\Omega }_2:=\left\{(x,y):x\ge 0,y\ge T\right\}$. Therefore, the global stabilities of $E_1,\widehat{E}$, and $E^{*}$ are investigated as follows.

Theorem 8.

If $n<{\displaystyle \frac{ad}{K}}$, then the predator extinction point $E_1=(K,0)$ is globally asymptotically stable in the region ${\Omega }_1$.

Proof. 

By specifying a positive Lyapunov function ${\mathcal{L}}_1(x,y)=\left(x-K-Kln{\displaystyle \frac{x}{K}}\right)+{\displaystyle \frac{my}{n}}$, and conforming to Lemma 4, we obtain

$$\begin{array}{cc}\hfill {}^C{\mathcal{D}}_t^{\alpha }{\mathcal{L}}_1(x,y)\le & \left(\frac{x-K}{x}\right){}^C{\mathcal{D}}_t^{\alpha }x+\frac{m}{n}{}^C{\mathcal{D}}_t^{\alpha }y\hfill \\ \hfill =& (x-K)\left(r-\frac{rx}{K}-\frac{my}{a+x}\right)+\frac{m}{n}\left(\frac{nxy}{a+x}-dy\right)\hfill \\ \hfill =& 2rx-rK-\frac{r{x}^2}{K}+\frac{mKy}{a+x}-\frac{dmy}{n}\hfill \\ \hfill =& -\frac{r}{K}{\left(x-K\right)}^2+\frac{mKy}{a+x}-\frac{dmy}{n}\hfill \\ \hfill \le & -\frac{r}{K}{\left(x-K\right)}^2+\frac{mKy}{a}-\frac{dmy}{n}\hfill \\ \hfill \le & -\left(\frac{d}{n}-\frac{K}{a}\right)my.\hfill \end{array}$$

Owing to the fact that $n<{\displaystyle \frac{ad}{K}}$, we have ${}^C{\mathcal{D}}_t^{\alpha }{\mathcal{L}}_1(x,y)\le 0$. In consequence of Lemma 5, $E_1$ is globally asymptotically stable in the region ${\Omega }_1$. □

Remark 2.

Notice $E_1$ is locally asymptotically stable if $n<{\displaystyle \frac{ad}{K}}+d$ and is globally asymptotically stable if $n<{\displaystyle \frac{ad}{K}}$. Hence, if the global stability condition is fulfilled, then the local stability is also achieved but not vice versa. This fact reinforces that the global stability condition has a larger attracting region than that of the local stability condition.

Theorem 9.

If $\left(K-\widehat{x}\right)(a+\widehat{x})<a^2$, then the first co-existence point $\widehat{E}=\left(\widehat{x},{\displaystyle \frac{\left(K-\widehat{x}\right)(a+\widehat{x})r}{mK}}\right)$ is globally asymptotically stable in the region ${\Omega }_1$.

Proof. 

Let $\widehat{E}=\left(\widehat{x},\widehat{y}\right)$ where $\widehat{y}={\displaystyle \frac{\left(K-\widehat{x}\right)(a+\widehat{x})r}{mK}}$ and $\phi$ is a positive constant. By considering a Lyapunov function ${\mathcal{L}}_2(x,y)=\left[x-\widehat{x}-\widehat{x}ln{\displaystyle \frac{x}{\widehat{x}}}\right]+\phi \left[y-\widehat{y}-\widehat{y}ln{\displaystyle \frac{y}{\widehat{y}}}\right]$, and using Lemma 4, we get

$$\begin{array}{cc}\hfill {}^C{\mathcal{D}}_t^{\alpha }{\mathcal{L}}_2(x,y)\le & \left(\frac{x-\widehat{x}}{x}\right){}^C{\mathcal{D}}_t^{\alpha }x+\phi \left(\frac{y-\widehat{y}}{y}\right){}^C{\mathcal{D}}_t^{\alpha }y\hfill \\ \hfill =& (x-\widehat{x})\left(r-\frac{rx}{K}-\frac{my}{a+x}\right)+\phi (y-\widehat{y})\left(\frac{nx}{a+x}-d\right)\hfill \\ \hfill =& (x-\widehat{x})\left(\frac{r\widehat{x}}{K}+\frac{m\widehat{y}}{a+\widehat{x}}-\frac{rx}{K}-\frac{my}{a+x}\right)+\phi (y-\widehat{y})\left(\frac{nx}{a+x}-\frac{n\widehat{x}}{a+\widehat{x}}\right)\hfill \\ \hfill =& -(x-\widehat{x})\left(\left(x-\widehat{x}\right)\frac{r}{K}+\left(\frac{y}{a+x}-\frac{\widehat{y}}{a+\widehat{x}}\right)m\right)+n\phi (y-\widehat{y})\left(\frac{x}{a+x}-\frac{\widehat{x}}{a+\widehat{x}}\right)\hfill \\ \hfill =& -{\left(x-\widehat{x}\right)}^2\frac{r}{K}-(x-\widehat{x})\left(\frac{(y-\widehat{y})(a+\widehat{x})+(\widehat{x}-x)\widehat{y}}{(a+x)(a+\widehat{x})}\right)m+\left(\frac{(x-\widehat{x})(y-\widehat{y})}{(a+x)(a+\widehat{x})}\right)an\phi \hfill \\ \hfill =& -{\left(x-\widehat{x}\right)}^2\frac{r}{K}-\left(\frac{(x-\widehat{x})(y-\widehat{y})(a+\widehat{x})}{(a+x)(a+\widehat{x})}\right)m+\left(\frac{{(x-\widehat{x})}^2\widehat{y}}{(a+x)(a+\widehat{x})}\right)m+\hfill \\ \hfill & \left(\frac{(x-\widehat{x})(y-\widehat{y})}{(a+x)(a+\widehat{x})}\right)an\phi \hfill \\ \hfill \le & -{\left(x-\widehat{x}\right)}^2\left(\frac{r}{K}-\frac{m\widehat{y}}{a^2}\right)+\left(\frac{(x-\widehat{x})(y-\widehat{y})}{(a+x)(a+\widehat{x})}\right)(an\phi -(a+\widehat{x})m).\hfill \end{array}$$

By choosing $\phi ={\displaystyle \frac{(a+\widehat{x})m}{an}}$ and substituting the value of $\widehat{y}$, we get

$${}^C{\mathcal{D}}_t^{\alpha }{\mathcal{L}}_2(x,y)\le -{\left(x-\widehat{x}\right)}^2\left(a^2-\left(K-\widehat{x}\right)(a+\widehat{x})\right)\frac{r}{a^{2}K}.$$

Therefore, if $\left(K-\widehat{x}\right)(a+\widehat{x})<a^2$, then ${}^C{\mathcal{D}}_t^{\alpha }{\mathcal{L}}_2(x,y)\le 0$. Using Lemma 5, we conclude that $\widehat{E}$ is globally asymptotically stable in the region ${\Omega }_1$. □

Based on Theorem 2, $x\left(t\right)$ and $y\left(t\right)$ are bounded. Let $\Psi$ be the upper bound of $y\left(t\right)$ such that $0<T\le y\left(t\right)\le \Psi$. The global stability of $E^{*}$ is stated in the following theorem.

Theorem 10.

If $y^{*}<min\left\{{\displaystyle \frac{a^{2}r}{mK}},{\displaystyle \frac{cT}{\Psi -T}}+T\right\}$, then the second co-existence point $E^{*}=(x^{*},y^{*})$ is globally asymptotically stable in the region ${\Omega }_2$.

Proof. 

We consider a positive Lyapunov function ${\mathcal{L}}_3\left(E^{*}\right)=\left[x-x^{*}-x^{*}ln{\displaystyle \frac{x}{x^{*}}}\right]+{\phi }^{*}\left[y-y^{*}-y^{*}ln{\displaystyle \frac{y}{y^{*}}}\right]$. According to Lemma 4, the fractional derivative of ${\mathcal{L}}_3\left(E^{*}\right)$ satisfies

$$\begin{array}{cc}\hfill {}^C{\mathcal{D}}_t^{\alpha }{\mathcal{L}}_2(x,y)\le & \left(\frac{x-x^{*}}{x}\right){}^C{\mathcal{D}}_t^{\alpha }x+{\phi }^{*}\left(\frac{y-y^{*}}{y}\right){}^C{\mathcal{D}}_t^{\alpha }y\hfill \\ \hfill =& (x-x^{*})\left(r-\frac{rx}{K}-\frac{my}{a+x}\right)+{\phi }^{*}(y-y^{*})\left(\frac{nx}{a+x}-d-\frac{h(y-T)}{(c+y-T)y}\right)\hfill \\ \hfill =& (x-x^{*})\left(\frac{r{x}^{*}}{K}+\frac{m{y}^{*}}{a+x^{*}}-\frac{rx}{K}-\frac{my}{a+x}\right)\hfill \\ \hfill & +{\phi }^{*}(y-y^{*})\left(\frac{nx}{a+x}-\frac{n{x}^{*}}{a+x^{*}}+\frac{h(y^{*}-T)}{(c+y^{*}-T)y^{*}}-\frac{h(y-T)}{(c+y-T)y}\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill =& -{(x-x^{*})}^2\frac{r}{K}-(x-x^{*})\left(\frac{(y-y^{*})(a+x^{*})m-(x-x^{*})m{y}^{*}}{(a+x^{*})(a+x)}\right)\hfill \\ \hfill & +{\phi }^{*}(y-y^{*})\left(\frac{(x-x^{*})an}{(a+x^{*})(a+x)}-\frac{(y-y^{*})chT-(y-y^{*})(y^{*}-T)(y-T)h}{(c+y^{*}-T)(c+y-T)y^{*}y}\right)\hfill \\ \hfill =& -{(x-x^{*})}^2\frac{r}{K}-\frac{(x-x^{*})(y-y^{*})(a+x^{*})m}{(a+x^{*})(a+x)}+\frac{{(x-x^{*})}^{2}m{y}^{*}}{(a+x^{*})(a+x)}\hfill \\ \hfill & +\frac{(x-x^{*})(y-y^{*})an{\phi }^{*}}{(a+x^{*})(a+x)}-\frac{(cT-(y^{*}-T)(y-T)){(y-y^{*})}^{2}h{\phi }^{*}}{(c+y^{*}-T)(c+y-T)y^{*}y}\hfill \\ \hfill \le & -{(x-x^{*})}^2\frac{r}{K}-\frac{(x-x^{*})(y-y^{*})(a+x^{*})m}{(a+x^{*})(a+x)}+{(x-x^{*})}^2\frac{m{y}^{*}}{a^2}\hfill \\ \hfill & +\frac{(x-x^{*})(y-y^{*})an{\phi }^{*}}{(a+x^{*})(a+x)}-\frac{(cT-(y^{*}-T)(y-T)){(y-y^{*})}^{2}h{\phi }^{*}}{(c+y^{*}-T)(c+y-T)y^{*}y}\hfill \end{array}$$

By taking ${\phi }^{*}={\displaystyle \frac{(a+x^{*})m}{an}}$ and remembering that $y\left(t\right)<\Psi$, we obtain

$$\begin{array}{cc}\hfill {}^C{\mathcal{D}}_t^{\alpha }{\mathcal{L}}_2\left(E^{*}\right)\le & -{(x-x^{*})}^2\left(\frac{r}{K}-\frac{m{y}^{*}}{a^2}\right)-\frac{(cT-(y^{*}-T)(\Psi -T)){(y-y^{*})}^2(a+x^{*})mh}{(c+y^{*}-T)(c+y-T)an{y}^{*}y}.\hfill \end{array}$$

It is easily confirmed that, if $y^{*}<min\left\{{\displaystyle \frac{a^{2}r}{mK}},{\displaystyle \frac{cT}{\Psi -T}}+T\right\}$, then ${}^C{\mathcal{D}}_t^{\alpha }{\mathcal{L}}_2(x,y)\le 0$. Based on Lemma 5, $\widehat{E}$ is globally asymptotically stable in the region ${\Omega }_2$. □

2.7. The Existence of Hopf Bifurcation

The Hopf bifurcation is a local phenomenon when a stable equilibrium point loses its stability and all nearby solutions converge to a periodic solution namely limit-cycle if a parameter is varied [54,60,61]. It is shown that many fractional-order models involving the Caputo operator undergo a Hopf bifurcation which is driven by the order of the derivative (see [2,17,34,53,62]). The difference between the Hopf bifurcation in the integer-order model and that in the fractional-order model lies on the property of the limit-cycle. In the integer-order model, the limit-cycle is a periodic orbit which does not exist in the fractional-order model [63]. In the fractional-order model, the limit-cycle is not a periodic solution, but all nearby solutions converge to a limit-cycle [56,62].

Adapted from Theorem 3 in [62], for two dimensional Caputo fractional-order system, a Hopf bifurcation occurs when the eigenvalues ${\lambda }_{1,2}$ of the Jacobian matrix evaluated at the equilibrium point satisfy the following conditions:

(i)

${\lambda }_{1,2}=\psi \pm \omega i$ where $\psi >0$;

(ii)

$m\left({\alpha }^{*}\right)={\alpha }^{*}\pi /2-{min}_{1\le i\le 2}\left|arg\left({\lambda }_i\right)\right|=0$; and

(iii)

${\left.{\displaystyle \frac{dm\left(\alpha \right)}{d\alpha }}\right|}_{\alpha ={\alpha }^{*}}\ne 0$.

Now, consider the stability condition in Theorems 6 and 7. For $y<T$, the Jacobian matrix of model (6) evaluated at $\widehat{E}$ has a pair of complex eigenvalues if $\Delta >0$. We can easily confirm that, if $\widehat{x}<{\displaystyle \frac{K-a}{2}}$, then the real part of the eigenvalues are positive. We also have that $m\left(\widehat{\alpha }\right)=0$ and ${\left.{\displaystyle \frac{dm\left(\alpha \right)}{d\alpha }}\right|}_{\alpha =\widehat{\alpha }}\ne 0$. Therefore, $\widehat{E}$ undergoes a Hopf bifurcation when $\alpha$ crosses $\widehat{\alpha }$. A similar circumstance also occurs when $y\ge T$. When ${\xi }^2<4\theta$, the Jacobian matrix of model (6) has a pair of complex eigenvalues. The real part of the eigenvalues are also positive when $\xi >0$. We can also check that $m\left({\alpha }^{*}\right)=0$ and ${\left.{\displaystyle \frac{dm\left(\alpha \right)}{d\alpha }}\right|}_{\alpha ={\alpha }^{*}}\ne 0$. This means the Hopf bifurcation also occurs when $y\ge T$. Therefore, we have the following theorem.

Theorem 11.

(i) 

Let $\Delta >0$ and $\widehat{x}<{\displaystyle \frac{K-a}{2}}$. The first co-existence point $\widehat{E}$ undergoes a Hopf bifurcation when α passes through $\widehat{\alpha }$ in the region ${\Omega }_1$.

(ii) 

Let ${\xi }^2<4\theta$ and $\xi >0$. The second co-existence point $E^{*}$ undergoes a Hopf bifurcation when α passes through ${\alpha }^{*}$ in the region ${\Omega }_2$.

3. The Atangana–Baleanu Fractional-Order Rosenzweig–MacArthur Model with THP in Predator

3.1. Model Formulation

In this section, we consider a fractional-order Rosenzweig–MacArthur model with THP in predator involving the Atangana–Baleanu operator. Specifically, we consider the Atangana–Baleanu operator in Caputo sense of order-$\alpha$ which is defined as follows.

Definition 3.

[40] Suppose $0<\alpha \le 1$. The Atangana–Baleanu fractional integral and derivative in Caputo sense of order-α (ABC) is defined by

$$\begin{array}{cc}\hfill {}^{ABC}{\mathcal{I}}_t^{\alpha }f\left(t\right)=& \frac{1-\alpha }{B\left(\alpha \right)}f\left(t\right)+\frac{\alpha }{\Gamma \left(\alpha \right)B\left(\alpha \right)}{\int }_0^t{(t-s)}^{\alpha -1}f\left(s\right)ds,\hfill \\ \hfill {}^{ABC}{\mathcal{D}}_t^{\alpha }f\left(t\right)=& \frac{B\left(\alpha \right)}{1-\alpha }{\int }_0^t{E}_{\alpha }\left[-\frac{\alpha }{1-\alpha }{(t-s)}^{\alpha }\right]f^{\prime }\left(s\right)ds,\hfill \end{array}$$

where $t\ge 0,$ $f\in C^n([0,+\infty ),\mathbb{R})$, $E_{\alpha }\left(t\right)={\sum }_{k=0}^{\infty }{\displaystyle \frac{t^k}{\Gamma (\alpha k+1)}}$ is the Mittag–Leffler function and $B\left(\alpha \right)$ is a normalization function with $B\left(0\right)=B\left(1\right)=1$. In this paper, we define $B\left(\alpha \right)=1-\alpha +{\displaystyle \frac{\alpha }{\Gamma \left(\alpha \right)}}$.

By applying Definition 3 to model (4), we get the following fractional-order model with ABC operator

$$\begin{array}{cc}\hfill {}^{ABC}{\mathcal{D}}_t^{\alpha }x=& rx\left(1-\frac{x}{K}\right)-\frac{mxy}{a+x}\equiv G_1,\hfill \\ \hfill {}^{ABC}{\mathcal{D}}_t^{\alpha }y=& \frac{nxy}{a+x}-dy-H\left(y\right)\equiv G_2.\hfill \end{array}$$

(16)

Due to the lack of analytical theory, model (16) is investigated numerically. However, we first show the existence and uniqueness of the solution of the model (16).

3.2. Existence and Uniqueness

We start this work by representing the Lipschitz condition of the kernels of model (16). Since the harvesting is performed by obeying threshold harvesting policy, we give the proof into two cases i.e., $y<T$ and $y\ge T$.

We start for $y\ge T$. Let $x,\overline{x},y$, and $\overline{y}$ be functions satisfying $∥x∥\le a_1$, $∥\overline{x}∥\le a_2$, $∥y∥\le b_1$, and $∥\overline{y}∥\le b_2$. Suppose that

$$\begin{array}{cc}\hfill g_1=& r+(a_1+a_2)\frac{r}{K}+\frac{my}{a},\hfill \\ \hfill g_2=& n+d+{\displaystyle \frac{h}{c}}.\hfill \end{array}$$

Therefore, we get

$$\begin{array}{cc}\hfill ∥G_1\left(x\right)-G_1\left(\overline{x}\right)∥=& ∥\left(rx\left(1-\frac{x}{K}\right)-\frac{mxy}{a+x}\right)-\left(r\overline{x}\left(1-\frac{\overline{x}}{K}\right)-\frac{m\overline{x}y}{a+\overline{x}}\right)∥\hfill \\ \hfill =& ∥rx-\frac{r{x}^2}{K}-\frac{mxy}{a+x}-r\overline{x}+\frac{r{\overline{x}}^2}{K}+\frac{m\overline{x}y}{a+\overline{x}}∥\hfill \\ \hfill =& ∥r(x-\overline{x})-\frac{r}{K}(x^2-{\overline{x}}^2)-my\left(\frac{x}{a+x}-\frac{\overline{x}}{a+\overline{x}}\right)∥\hfill \\ \hfill =& ∥r(x-\overline{x})-\frac{r}{K}(x+\overline{x})(x-\overline{x})-\frac{amy(x-\overline{x})}{(a+x)(a+\overline{x})}∥\hfill \\ \hfill \le & r∥x-\overline{x}∥+(a_1+a_2)\frac{r}{K}∥x-\overline{x}∥+\frac{my}{a}∥x-\overline{x}∥\hfill \\ \hfill =& \left(r+(a_1+a_2)\frac{r}{K}+\frac{my}{a}\right)∥x-\overline{x}∥\hfill \\ \hfill =& g_1∥x-\overline{x}∥,\hfill \end{array}$$

(17)

and

$$\begin{array}{cc}\hfill ∥G_2\left(y\right)-G_2\left(\overline{y}\right)∥=& ∥\left(\frac{nxy}{a+x}-dy-{\displaystyle \frac{h(y-T)}{c+(y-T)}}\right)-\left(\frac{nx\overline{y}}{a+x}-d\overline{y}-{\displaystyle \frac{h(\overline{y}-T)}{c+(\overline{y}-T)}}\right)∥\hfill \\ \hfill =& ∥\frac{nxy}{a+x}-dy-{\displaystyle \frac{h(y-T)}{c+(y-T)}}-\frac{nx\overline{y}}{a+x}+d\overline{y}+{\displaystyle \frac{h(\overline{y}-T)}{c+(\overline{y}-T)}}∥\hfill \\ \hfill =& ∥\frac{nx(y-\overline{y})}{a+x}-d(y-\overline{y})-{\displaystyle \frac{ch(y-\overline{y})}{(c+y-T)(c+\overline{y}-T)}}∥\hfill \\ \hfill \le & n∥y-\overline{y}∥+d∥y-\overline{y}∥+{\displaystyle \frac{h}{c}}∥y-\overline{y}∥\hfill \\ \hfill =& \left(n+d+{\displaystyle \frac{h}{c}}\right)∥y-\overline{y}∥\hfill \\ \hfill =& g_2∥y-\overline{y}∥.\hfill \end{array}$$

(18)

When $y<T$, by utilizing the similar manner, we achieve $∥G_2\left(y\right)-G_2\left(\overline{y}\right)∥\le g_3∥y-\overline{y}∥$ where $g_3=n+d$. Accordingly, the kernel of (16) satisfies the Lipschitz condition. Furthermore, if $0\le g_1<1$ and $0\le g_3<g_2<1$, then $G_1$ and $G_2$ are contracted.

Now, by employing the fixed-point theorem, the solution of model (16) is investigated. By utilizing the Atangana–Baleanu fractional integral operator, model (16) is transformed into the following Volterra-type integral equation.

$$\begin{array}{cc}\hfill x\left(t\right)-x\left(0\right)=& \frac{1-\alpha }{B\left(\alpha \right)}{G}_1(t,x)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}{\int }_0^t{(t-s)}^{\alpha -1}{G}_1(s,x)ds,\hfill \\ \hfill y\left(t\right)-y\left(0\right)=& \frac{1-\alpha }{B\left(\alpha \right)}{G}_2(t,y)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}{\int }_0^t{(t-s)}^{\alpha -1}{G}_2(s,y)ds.\hfill \end{array}$$

(19)

In a recursive form, Equation (19) is written as

$$\begin{array}{cc}\hfill x_n\left(t\right)=& \frac{1-\alpha }{B\left(\alpha \right)}{G}_1(t,x_{n-1})+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}{\int }_0^t{(t-s)}^{\alpha -1}{G}_1(s,x_{n-1})ds,\hfill \\ \hfill y_n\left(t\right)=& \frac{1-\alpha }{B\left(\alpha \right)}{G}_2(t,y_{n-1})+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}{\int }_0^t{(t-s)}^{\alpha -1}{G}_2(s,y_{n-1})ds.\hfill \end{array}$$

(20)

The associated initial conditions along with Equation (20) are $x_0\left(t\right)=x\left(0\right)$ and $y_0\left(t\right)=y\left(0\right)$. By considering Equation (20), we have the difference expression of successive terms as follows.

$$\begin{array}{cc}\hfill {\Phi }_{1,n}\left(t\right)=& x_n\left(t\right)-x_{n-1}\left(t\right)\hfill \\ \hfill =& \frac{1-\alpha }{B\left(\alpha \right)}\left(G_1(t,x_{n-1})-G_1(t,x_{n-2})\right)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}{\int }_0^t{(t-s)}^{\alpha -1}\left(G_1(t,x_{n-1})-G_1(t,x_{n-2})\right)ds,\hfill \\ \hfill {\Phi }_{2,n}\left(t\right)=& y_n\left(t\right)-y_{n-1}\left(t\right)\hfill \\ \hfill =& \frac{1-\alpha }{B\left(\alpha \right)}\left(G_2(t,y_{n-1})-G_2(t,y_{n-2})\right)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}{\int }_0^t{(t-s)}^{\alpha -1}\left(G_2(t,y_{n-1})-G_2(t,y_{n-2})\right)ds.\hfill \end{array}$$

(21)

According to Equation (21), we get $x_n\left(t\right)={\sum }_{i=1}^n{\Phi }_{1,i}\left(t\right)$ and $y_n\left(t\right)={\sum }_{i=1}^n{\Phi }_{2,i}\left(t\right)$. By applying Equations (17), (18) and (21), we have that

$$\begin{array}{cc}\hfill ∥{\Phi }_{1,n}\left(t\right)∥\le & \frac{1-\alpha }{B\left(\alpha \right)}{g}_1∥{\Phi }_{1,n-1}∥+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}{g}_1{\int }_0^t∥{\Phi }_{1,n-1}\left(s\right)∥{(t-s)}^{\alpha -1}ds,\hfill \\ \hfill ∥{\Phi }_{2,n}\left(t\right)∥\le & \frac{1-\alpha }{B\left(\alpha \right)}{g}_2∥{\Phi }_{2,n-1}∥+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}{g}_2{\int }_0^t∥{\Phi }_{2,n-1}\left(s\right)∥{(t-s)}^{\alpha -1}ds.\hfill \end{array}$$

(22)

Therefore, by using (22), the existence and uniqueness of model (16) is presented as follows.

Theorem 12.

Model (16) has a unique solution if we can find $t_0$ such that

$$\frac{(1-\alpha )g_i}{B\left(\alpha \right)}+\frac{t_0^{\alpha }{g}_i}{B\left(\alpha \right)\Gamma \left(\alpha \right)}<1,i=1,2,3.$$

(23)

Proof. 

Let $x\left(t\right)$ and $y\left(t\right)$ be bounded functions, and hence the Lipschitz condition is satisfied. Thus, according to Equation (22), we obtain the following inequalities.

$$\begin{array}{cc}\hfill ∥{\Phi }_{1,n}\left(t\right)∥\le & ∥x_0∥{\left(\frac{(1-\alpha )g_1}{B\left(\alpha \right)}+\frac{t^{\alpha }{g}_1}{B\left(\alpha \right)\Gamma \left(\alpha \right)}\right)}^n,\hfill \\ \hfill ∥{\Phi }_{2,n}\left(t\right)∥\le & ∥y_0∥{\left(\frac{(1-\alpha )g_2}{B\left(\alpha \right)}+\frac{t^{\alpha }{g}_2}{B\left(\alpha \right)\Gamma \left(\alpha \right)}\right)}^n.\hfill \end{array}$$

(24)

Therefore, the continuity and existence of solution are proved since $∥{\Phi }_{1,n}\left(t\right)∥\to 0$ and $∥{\Phi }_{2,n}\left(t\right)∥\to 0$ as $n\to \infty$ and $t=t_0$. Now, suppose that

$$\begin{array}{cc}\hfill x\left(t\right)-x\left(0\right)=& x_n\left(t\right)-{\mathrm{Y}}_{1,n}\left(t\right),\hfill \\ \hfill y\left(t\right)-y\left(0\right)=& y_n\left(t\right)-{\mathrm{Y}}_{2,n}\left(t\right).\hfill \end{array}$$

(25)

We confirm that

$$∥{\mathrm{Y}}_{1,n}\left(t\right)∥\le {\left(\frac{(1-\alpha )}{B\left(\alpha \right)}+\frac{t^{\alpha }}{B\left(\alpha \right)\Gamma \left(\alpha \right)}\right)}^{n+1}{g}_1^{n+1}$$

(26)

It is clear that $∥{\mathrm{Y}}_{1,n}\left(t\right)∥\to 0$ when $n\to \infty$. By using a similar manner, we acquire $∥{\mathrm{Y}}_{2,n}\left(t\right)∥\to 0,n\to \infty$. Finally, the uniqueness of solution for the model is proven. Suppose that there exists different set of solutions denote by $\tilde{x}\left(t\right)$ and $\tilde{y}\left(t\right)$; then, we have

$$x\left(t\right)-\tilde{x}\left(t\right)=\frac{1-\alpha }{B\left(\alpha \right)}(G_1(t,x)-G_1(t,\tilde{x}))+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}{\int }_0^t(G_1(s,x)-G_1(s,\tilde{x})){(t-s)}^{\alpha -1}ds.$$

(27)

Taking the norm for both sides and using a simplification as in (22) and (24), we obtain

$$∥x\left(t\right)-\tilde{x}\left(t\right)∥\left(1-\frac{(1-\alpha )g_1}{B\left(\alpha \right)}-\frac{t^{\alpha }{g}_1}{B\left(\alpha \right)\Gamma \left(\alpha \right)}\right)\le 0.$$

(28)

For $t=t_0$, we have (23), thus $∥x\left(t\right)-\tilde{x}\left(t\right)∥=0$ and hence $x\left(t\right)=\tilde{x}\left(t\right)$. Applying the same algebraic procedures, we can show that $y\left(t\right)=\tilde{y}\left(t\right)$. Therefore, the solution is unique. □

4. Numerical Simulations

In this section, the numerical simulations of Caputo model (6) and ABC model (16) are performed to illustrate the previous theoretical results. In the literature, there exist many numerical methods to solve a system of fractional differential equations, such as the Monte Carlo method [64], the Grünwald–Letnikov method [65,66] and the predictor–corrector method [67,68,69]. We apply the predictor–corrector scheme proposed by Diethelm et al. [67] to solve the Caputo fractional-order model and the predictor–corrector scheme proposed by Baleanu et al. [69] to solve the Atangana–Baleanu in Caputo sense model (ABC). Due to the limitation of field data, we use hypothetical parameter values that correspond to the theoretical results. The initial parameter values are given as follows.

$$r=0.5,K=1,m=0.3,a=0.2,d=0.1,h=0.1,T=0.9,c=0.1,\mathrm{and} \alpha =0.9.$$

(29)

In Figure 1, we plot a bifurcation diagram by varying the conversion rate of consumed prey into predator birth n in interval $[0.08,0.2]$. We notice that the parameter values (29) and the interval $0.08\le n\le 0.2$ lead to the non-existence of equilibrium point in ${\Omega }_2$. Therefore, the first numerical simulations are focused on the dynamics in ${\Omega }_1$.

For $0.08\le n<n_1^{*}=0.12$, Theorem 5 states that the predator extinction point $E_1=(1,0)$ is the only equilibrium point which is asymptotically stable. To see this behavior, we take $n=0.1$ and plot the phase portrait and the time series as in Figure 2. It is seen that all solutions with initial values in both ${\Omega }_1$ and ${\Omega }_2$ are convergent to $E_1$. When the initial value is in ${\Omega }_2$, then the solution is oscillating when it crosses the threshold harvesting level and eventually converges to $E_1$.

When n passes through $n_1^{*}$, the equilibrium point $E_1=(1,0)$ undergoes a forward bifurcation. In this case, there appear two equilibrium points, namely the unstable $E_1$ and an asymptotically stable $\widehat{E}$. Figure 1 shows that $\widehat{E}$ is asymptotically stable if $0.12<n\lesssim n_2^{*}=0.1557$. In Figure 3, we show the phase portrait and time series for the case of $n=0.14$. We see that $E_0=(0,0)$ and $E_1=(1,0)$ are saddle points, while $\widehat{E}=(0.5,0.5833)$ is asymptotically stable. This circumstance corresponds to Theorems 4–6 and 9.

Furthermore, if we increase the value of n such that $n>n_2^{*}$, then the equilibrium point $\widehat{E}$ loses its stability, and all solutions converge to a limit-cycle. This situation confirms the occurrence of a Hopf bifurcation driven by n. For example, if we select $n=0.2$ then all equilibrium points are unstable and all solutions eventually converge to limit cycle (see Figure 4).

Now, we perform simulation using the same parameter values as in Figure 4, but with a lower threshold value, namely $T=0.5$. In this case, there is no equilibrium point $\widehat{E}$ in ${\Omega }_1$, and $E^{*}=(0.5954,0.5364)$ occurs in ${\Omega }_2$. According to Theorem 7, $E^{*}$ is asymptotically stable. Such dynamics can be clearly seen in Figure 5. This simulation shows that by applying the THP when the interior equilibrium point is stable, we can choose a suitable constant of threshold so that the existence of both prey and predator are maintained.

Notice that, in Figure 2, Figure 3, Figure 4 and Figure 5, we see that both model with Caputo operator (6) and Atangana–Baleanu operator (16) have similar dynamical behavior. The noticeable difference between them is the orbit of solutions and the diameter of the limit-cycle. In Figure 4, the diameter of the limit-cycle of the model with ABC operator looks shorter than that of the Caputo operator, which may give different dynamics when a Hopf bifurcation occurs. To get more detail view, we plot a bifurcation diagram by varying the order of the fractional derivative $\left(\alpha \right)$ (see Figure 6). In this simulation, we use parameter values as in Figure 4 and vary the order-$\alpha$ in the interval $[0.6,0.9]$. It is clearly seen that, besides the diameter of the limit-cycle, the bifurcation points of Caputo model and ABC model are also different. The Caputo model has an earlier bifurcation point than that of the ABC model. To show this situation, we show some numerical simulations using different values of $\alpha$ (see Figure 7). For $\alpha =0.7$, the equilibrium point $\widehat{E}$ of both model are asymptotically stable. For $\alpha =0.772$, the equilibrium point $\widehat{E}$ of the Caputo model loses its stability, while that of the ABC model remains asymtotically stable. For $\alpha =0.83$, the equilibrium point $\widehat{E}$ of both models are unstable.

From the biological point of view, all previous numerical simulations show that the dynamical properties of both Caputo model and ABC model are similar except when the eigenvalues of the Jacobian matrix evaluated at the interior equilibrium point $\widehat{E}$ are a pair of complex conjugate with positive real part. There is a biological condition such that the prey and predator densities are eventually periodic for the Caputo model, while for ABC model, the densities of both predator and prey are eventually constant.

5. Conclusions

The dynamics of a Rosenzweig–MacArthur model with continuous threshold harvesting in predator involving the Caputo fractional-order derivative and ABC fractional-order derivative are studied. We prove the existence and uniqueness of solutions of both Caputo and ABC models. Particularly, we completely investigate the dynamics of the Caputo model including the non-negativity, boundedness, local stability, global stability, and the existence of Hopf bifurcation. From the biological meanings, the extinction of both populations never occurs since the origin point $\left(E_0\right)$ is a saddle point. Some of the situations that might occur are: (1) the predator goes extinct while the prey still survives, which is indicated by the stability of $E_1$; (2) both predator and prey co-exist and converge to a constant population density, which happens if the interior point $\widehat{E}$ or $E^{*}$ are asymptotically stable; and (3) both predator and prey co-exist where both population densities change periodically, namely when a Hopf bifurcation occurs. We show numerically that our model may undergo a forward bifurcation or a Hopf bifurcation. The Hopf bifurcation in models with both Caputo operator and ABC operator can be driven by the conversion rate of consumed prey into the predator birth rate or by the order of fractional derivative. Our numerical simulations show that the Hopf bifurcation point of both models are different.