**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–8], the fear effect [9], the Allee effect [10,11], the intraspecific competition [12,13], the cannibalism [14], the infectious disease [15–17], and so forth.

On the other hand, the modeling also contemplates the optimal managemen<sup>t</sup> 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–22]. To protect the population from over-exploitation during the harvesting, some managemen<sup>t</sup> techniques have

been established. One of the famous technique is a continuous threshold harvesting policy (THP) (see [23–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{aligned} \frac{dx}{dt} &= rx\left(1 - \frac{x}{K}\right) - \frac{mxy}{a+x},\\ \frac{dy}{dt} &= \frac{nxy}{a+x} - dy - H(y), \end{aligned} \tag{1}$$

where

$$H(y) = \begin{cases} \begin{array}{c} 0 \\ \frac{h(y-T)}{c+(y-T)} \end{array} & \text{if } y < T. \end{cases}$$

We describe the biological interpretation of variables and parameters of model (1) in Table 1. 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).


**Table 1.** Description of variables and parameters of the 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–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–45]. Nevertheless, the Caputo operator has more complex analytical tools in investigating the dynamics of the model such as the local stability [46–50], the global stability [51,52], bifurcation theory [52–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 *ddt* at the left hand sides of model (1) are replaced by the fractional-order derivatives <sup>D</sup>*<sup>α</sup>t*, then we obtain

$$\begin{aligned} \mathcal{D}\_t^a \mathbf{x} &= \mathbb{f} \mathbf{x} \left( 1 - \frac{\mathbf{x}}{K} \right) - \frac{\hat{m} \mathbf{x} y}{a + \mathbf{x}}, \\ \mathcal{D}\_t^a y &= \frac{\hat{n} \mathbf{x} y}{a + \mathbf{x}} - \hat{d} y - H(y). \end{aligned} \tag{2}$$

Note that the left hand sides of model (2) have the dimension of (time)−*α*, while the parameters at the right hand sides of model (2) such as *r*ˆ, *m*<sup>ˆ</sup> , *n*ˆ, ˆ *d*, and ˆ *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 ge<sup>t</sup> the following model

$$\begin{aligned} \mathcal{D}\_t^\alpha \mathbf{x} &= \mathfrak{f}^\alpha \mathbf{x} \left( 1 - \frac{\mathbf{x}}{K} \right) - \frac{\mathfrak{M}^\alpha \mathbf{x} \mathfrak{y}}{a + \mathbf{x}}, \\ \mathcal{D}\_t^\alpha \mathbf{y} &= \frac{\hat{\mathfrak{n}}^a \mathfrak{x} \mathfrak{y}}{a + \mathbf{x}} - \hat{d}^\alpha \mathfrak{y} - H(\mathfrak{y}), \end{aligned} \tag{3}$$

where

$$H(y) = \begin{cases} \begin{array}{c} 0 \\ \hat{h}^a(y-T) \end{array} & \text{if } y < T. \\\ \begin{array}{c} \hat{h}^a(y-T) \\ \hline c+(y-T) \end{array} & \text{if } y \ge T. \end{cases}$$

By applying new parameters *r* = *<sup>r</sup>*<sup>ˆ</sup>*α*, *m* = *m*<sup>ˆ</sup> *α*, *n* = *n*ˆ *α*, *d* = ˆ *d<sup>α</sup>*, and *h* = ˆ *h<sup>α</sup>*, we obtain

$$\begin{aligned} \mathcal{D}\_t^a \mathbf{x} &= r \mathbf{x} \left( 1 - \frac{\mathbf{x}}{K} \right) - \frac{m \mathbf{x} \mathbf{y}}{a + \mathbf{x}}, \\ \mathcal{D}\_t^a \mathbf{y} &= \frac{\mathbf{x} \mathbf{y}}{a + \mathbf{x}} - dy - H(\mathbf{y}), \end{aligned} \tag{4}$$

where

$$H(y) = \begin{cases} \begin{array}{c} 0 \\ \frac{h(y-T)}{c+(y-T)} \end{array} & \text{if } y < T. \end{cases}$$

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 α* ∈ (0, 1]*, f* ∈ *Cn*([0, <sup>+</sup>∞), <sup>R</sup>)*, and* <sup>Γ</sup>(·) *is the Gamma function. The Caputo fractional derivative of order-α is defined by*

$$\, \, ^C \mathcal{D}\_t^a f(t) = \frac{1}{\Gamma(1-a)} \int\_0^t (t-s)^{-a} f'(s)ds, t \ge 0. \tag{5}$$

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

$$\begin{aligned} \, ^C \mathcal{D}\_t^a \mathbf{x} &= r \, \times \left( 1 - \frac{\mathbf{x}}{K} \right) - \frac{m \mathbf{x} \, y}{a + \mathbf{x}} \equiv \mathbf{F}\_1, \\ ^C \mathcal{D}\_t^a \mathbf{y} &= \frac{n \mathbf{x} \, y}{a + \mathbf{x}} - dy - H(y) \equiv \mathbf{F}\_2. \end{aligned} \tag{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*

$$\prescript{C}{}{\mathcal{D}}\_t^{\mathfrak{a}} \mathbf{x}(t) = f(t, \mathbf{x}(t)),\\ t > 0, \mathbf{x}(0) \ge 0, \mathfrak{a} \in (0, 1], \tag{7}$$

*where f* : (0, ∞) × Θ → R*<sup>n</sup>*, Θ ⊆ R*n. A unique solution of Equation (7) on* (0, ∞) × Θ *exists if f*(*<sup>t</sup>*, *x*(*t*)) *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* ≥ *T* and *y* < *T*. We start from *y* ≥ *T*. Consider the region Θ × [0, *T*+] where Θ := &(*<sup>x</sup>*, *y*) ∈ R<sup>2</sup> : max (|*x*|, |*y*|) ≤ *γ*, *y* ≥ *<sup>T</sup>*', *T*+ < +<sup>∞</sup>, and a mapping *<sup>F</sup>*(Λ)=(*<sup>F</sup>*1(Λ), *<sup>F</sup>*2(Λ)). For any Λ = (*<sup>x</sup>*,*y*) ∈ Θ and Λ¯ = (*x*¯, *y*¯) ∈ Θ, we obtain

99*F*(Λ) − *F*(Λ¯ )99 = *<sup>F</sup>*1(Λ) − *<sup>F</sup>*1(Λ¯ ) + *<sup>F</sup>*2(Λ) − *<sup>F</sup>*2(Λ¯ ) = *rx* 1 − *xK* − *mxy a* + *x* − *rx*¯ 1 − *x*¯*K* − *mx*¯*y*¯ *a* + *x*¯ + *nxy a* + *x* − *dy* − *<sup>H</sup>*(*y*) − *nx*¯*y*¯ *a* + *x*¯ − *dy*¯ − *<sup>H</sup>*(*y*¯) = *r*(*x* − *x*¯) − *rK x*2 − *x*¯2 − *m xy a* + *x* − *x*¯*y*¯ *a* + *x*¯ + *n xy a* + *x* − *x*¯*y*¯ *a* + *x*¯ − *d* (*y* − *y*¯) − (*H*(*y*) − *<sup>H</sup>*(*y*¯)) =*r* |*x* − *x*¯| + *rK* |*x* + *x*¯| |*x* − *x*¯| + (*m* + *n*) *xy a* + *x* − *x*¯*y*¯ *a* + *x*¯ + *d* |*y* − *y*¯| + *h*(*y* − *T*) *c* + (*y* − *T*) − *h*(*y*¯ − *T*) *c* + (*y*¯ − *T*) =*r* |*x* − *x*¯| + *rK* |*x* + *x*¯| |*x* − *x*¯| + (*m* + *n*) *ay*(*x* − *x*¯)+(*ax*¯ + *xx*¯ )(*y* − *y*¯) (*a* + *x*)(*a* + *x*¯) + *d* |*y* − *y*¯| + *ch*(*y* − *y*¯) (*c* + *y* − *<sup>T</sup>*)(*c* + *y*¯ − *T*) ≤*r* |*x* − *x*¯| + 2*γr K* |*x* − *x*¯| + *m* + *n a*2 |*ay*(*x* − *x*¯)+(*ax*¯ + *xx*¯ )(*y* − *y*¯)| + *d* |*y* − *y*¯| + *hc* |*y* − *y*¯| ≤*r* |*x* − *x*¯| + 2*γr K* |*x* − *x*¯| + (*m* + *n*)*γ a* |*x* − *x*¯| + (*m* + *n*)(*aγ* + *γ*<sup>2</sup>) *a*2 |*y* − *y*¯| + *d* |*y* − *y*¯| + *hc* |*y* − *y*¯| = *r* + 2*γr K* + (*m* + *n*)*γ a* |*x* − *x*¯| + (*m* + *n*)(*aγ* + *γ*<sup>2</sup>) *a*2 + *d* + *hc* |*y* − *y*¯| ≤*M*1 99Λ − Λ¯ 99 ,

where *M*1 = max *r* + 2*γr K* + (*m* + *n*)*γ a* , (*m* + *n*)(*aγ* + *γ*<sup>2</sup>) *a*2 + *d* + *hc* . Hence, *F*(Λ) satisfies the Lipschitz condition for *y* ≥ *T*. By using similar approaches, when *y* < *T*, it is easy to check that 99*F*(Λ) − *F*(Λ¯ )99 ≤ *M*2 99Λ − Λ¯ 99, where *M*2 = max *r* + 2*γr K* + (*m* + *n*)*γ a* , (*m* + *n*)(*aγ* + *γ*<sup>2</sup>) *a*2 + *d* and hence the Lipschitz condition is also satisfied. According to Lemma 1, model (6) with non-negative initial condition has a unique solution <sup>Λ</sup>(*t*)=(*x*(*t*), *y*(*t*)) ∈ Θ. Thus, we establish the following theorem.

**Theorem 1.** *For each non-negative initial condition* (*<sup>x</sup>*0, *y*0) ∈ Θ*, there exists a unique solution* (*x*(*t*), *y*(*t*)) ∈ Θ *of model* (6)*, which is defined for all t* ≥ 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 < *α* ≤ 1*. Suppose that f*(*t*) ∈ *<sup>C</sup>*[*<sup>a</sup>*, *b*] *and <sup>C</sup>*D*<sup>α</sup>t f*(*t*) ∈ *<sup>C</sup>*[*<sup>a</sup>*, *b*]*. If <sup>C</sup>*D*<sup>α</sup>t f*(*t*) ≥ 0, ∀ *t* ∈ (*a*, *b*)*, then f*(*t*) *is a non-decreasing function for each t* ∈ [*a*, *b*]*. If <sup>C</sup>*D*<sup>α</sup>t f*(*t*) ≤ 0, ∀*t* ∈ (*a*, *b*)*, then f*(*t*) *is a non-increasing function for each t* ∈ [*a*, *b*]*.*

**Lemma 3.** *(Standard comparison theorem for Caputo fractional-order derivative [31]). Let x*(*t*) ∈ *C* ([0, <sup>+</sup>∞))*. If x*(*t*) *satisfies <sup>C</sup>*D*<sup>α</sup>t x*(*t*) + *<sup>λ</sup>x*(*t*) ≤ *μ*, *x*(0) = *x*0*, where α* ∈ (0, <sup>1</sup>],(*<sup>λ</sup>*, *μ*) ∈ R<sup>2</sup> *and λ* = 0*, then x*(*t*) ≤ *x*0 − *μλ <sup>E</sup>α*[−*λt<sup>α</sup>*] + *μλ.*

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*(*t*) ≥ 0 for all *t* → ∞. Suppose that it is incorrect; then, we can find *t*1 > 0 such that

$$\begin{cases} \begin{array}{rcl} \mathbf{x}(t) &> & \mathbf{0}, \mathbf{0} \le t < t\_{1\prime} \\ \mathbf{x}(t\_{1}) &=& \mathbf{0}, \\ \mathbf{x}(t\_{1}^{+}) &< & \mathbf{0} \end{array} \end{cases} \tag{8}$$

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

$$\left. \, \prescript{C}{}{\mathcal{D}}\_t^{\alpha} \mathbf{x}(t\_1) \right|\_{\mathbf{x}(t\_1) = 0} = 0. \tag{9}$$

Based on Lemma 2, we have *<sup>x</sup>*(*t*+1 ) = 0, which contradicts the fact that *<sup>x</sup>*(*t*+1 ) < 0. Thus, *x*(*t*) ≥ 0 for all *t* ≥ 0. Using a similar procedure, we conclude *y*(*t*) ≥ 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 V(*t*) = *x* + *myn*, we ge<sup>t</sup>

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

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

$$\mathcal{V}(t) \le \left(\mathcal{V}(0) - \frac{(d+r)^2 K}{4r}\right) E\_a \left[-d(t)^a\right] + \frac{(d+r)^2 K}{4r}.$$

where *Eα* is the one-parameter Mittag–Leffler function. Since *Eα* [−*d*(*t*)*α*] → 0 as *t* → 0, we acquire V(*t*) → (*d* + *r*)<sup>2</sup>*<sup>K</sup>* 4*r* for *t* → ∞. Hence, with non-negative initial condition, all solutions are restricted to the region Θ*M* where

$$\Theta\_M := \left\{ (\mathbf{x}, \mathbf{y}) \in \mathbb{R}\_+^2 : \mathbf{x} + \frac{my}{n} \le \frac{(d+r)^2 K}{4r} + \varepsilon, \varepsilon > 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^\mathfrak{a}\vec{\mathfrak{x}} = \vec{f}(\vec{\mathfrak{x}}); \vec{\mathfrak{x}}(0) = \vec{\mathfrak{x}}\_0; \mathfrak{a} \in (0, 1]. \tag{10}$$

*A point x*<sup>∗</sup> *is called an equilibrium point of Equation* (10) *if it satisfies f*(*x*<sup>∗</sup>) = 0*. Particularly, it is called the biological equilibrium point of Equation* (10) *if it satisfies x*<sup>∗</sup> ≥ 0*.*

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

$$\begin{cases} r - \frac{r\chi}{K} - \frac{my}{a+x} \Big] \ge -0, \\ \frac{nxy}{a+x} - dy - H(y) = 0. \end{cases} \tag{11}$$

Thus, we ge<sup>t</sup> four feasible biological equilibrium points as follows.

(i) The equilibrium points when *y* < *T* are


$$\begin{aligned} \text{(i.c)} \quad &\text{the first co-existence point } \mathcal{E} = \left( \pounds, \frac{\left( \mathcal{K} - \hat{\mathfrak{x}} \right) \left( a + \hat{\mathfrak{x}} \right) r}{m \mathcal{K}} \right), \text{with } \mathfrak{x} = \frac{ad}{n - d}, \text{ which exists if } \\ &n > \frac{ad}{K} + d \text{ and } \left( \mathcal{K} - \hat{\mathfrak{x}} \right) \left( a + \hat{\mathfrak{x}} \right) < \frac{TmK}{r}. \end{aligned}$$

(ii) The equilibrium point when *y* ≥ *T* is the second co-existence point *E*∗ = (*x*<sup>∗</sup>, *y*<sup>∗</sup>) where *y*∗ = (*K* − *x*<sup>∗</sup>)(*a* + *x*<sup>∗</sup>)*r mK* and *x*<sup>∗</sup> is the positive roots of polynomial *β*1*x*<sup>4</sup> + *β*2*x*<sup>3</sup> + *β*3*x*<sup>2</sup> + *β*4*<sup>x</sup>* + *β*5 = 0 where

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

*E*∗ exists if 0 < *x*<sup>∗</sup> < *K* and (*K* − *x*<sup>∗</sup>)(*a* + *x*<sup>∗</sup>) ≥ *TmKr* .

#### *2.5. Local Asymptotic Stability*

In this part, we discuss the local stability of *E*0, *E*1, *E* ˆ , and *E*<sup>∗</sup>. For this aim, we need the following theorem.

**Theorem 3.** *(Matignon condition [48,59]) The equilibrium point x*<sup>∗</sup> *of system* (10) *is locally asymptotically stable if all eigenvalues λj of the Jacobian matrix J* = *∂ f* /*∂x evaluated at x*<sup>∗</sup> *satisfy* | arg(*<sup>λ</sup>j*)| > *απ*/2*. If there exists at least one eigenvalue satisfying* | arg(*<sup>λ</sup>k*)| > *απ*/2 *and* | arg(*<sup>λ</sup>l*)| < *απ*/2*, k* = *l, then x*<sup>∗</sup> *is a saddle-point.*

> Now, we present Theorems 4–7 to show the local stability properties of *E*0, *E*1, *E* ˆ , and *E*<sup>∗</sup>.

**Theorem 4.** *The origin point E*0 = (0, 0) *is always a saddle point.*

**Proof.** When *E*0 = (0, <sup>0</sup>), the model (6) has Jacobian matrix *J*(*<sup>E</sup>*0) = # *r* 0 0 −*d* \$, where its eigenvalues are *λ*1 = *r* > 0 and *λ*2 = −*d* < 0. It is clear that |arg(*<sup>λ</sup>*1)| = 0 < *απ*/2 and |arg(*<sup>λ</sup>*2)| = *π* > *απ*/2. Therefore, based on Theorem 3, *E*0 is always a saddle point.

**Theorem 5.** *The predator extinction point E*1 = (*<sup>K</sup>*, 0) *is locally asymptotically stable if n* < *adK* + *d. Otherwise, if n* > *adK* + *d, then E*1 = (*<sup>K</sup>*, 0) *is a saddle point.*

**Proof.** The Jacobian matrix of model (6) evaluated at *E*1 is *J*(*<sup>E</sup>*1) = # −*r* − *mK a*+*K* 0 *nK a*+*K* − *d* \$. The eigenvalues of *J*(*<sup>E</sup>*1) are *λ*1 = −*r* < 0 and *λ*2 = *nK a* + *K* − *d*. Clearly, |arg(*<sup>λ</sup>*1)| = *π* > *απ*/2 and |arg(*<sup>λ</sup>*2)| = *π* > *απ*/2 if *n* < *adK* + *d* and |arg(*<sup>λ</sup>*2)| = 0 < *απ*/2 if *n* > *adK* + *d*. Hence, we have the theorem.

**Remark 1.** *It is noted that the existence condition for the first co-existence point E*ˆ *contradicts the stability condition of E*1*. Consequently, if E*1 *is locally asymptotically stable, then E*ˆ *does not exist. This condition also indicates the existence of forward bifurcation, which is confirmed numerically in the next section.*

**Theorem 6.** *Let* Δ = 4 (*K* − *x*ˆ) *anrx*<sup>ˆ</sup> (*a* + *x*<sup>ˆ</sup>)<sup>2</sup>*<sup>K</sup>* − *K* − *x*ˆ *a* + *x*ˆ − 12 *r*2*x*ˆ2 *K*<sup>2</sup> *and α*ˆ = 2*π* tan−<sup>1</sup> ! √Δ(*a* + *x*<sup>ˆ</sup>)*<sup>K</sup>* (*K* − *a* − <sup>2</sup>*x*<sup>ˆ</sup>)*rx*<sup>ˆ</sup> "*. Suppose that one of the following statements is obeyed.*

$$\begin{aligned} \text{(i)} \qquad &\quad \text{f} > \frac{K-a}{2}; \text{or} \\ \text{(ii)} \qquad &\quad \text{f} < \frac{K-a}{2}, \text{\(\Lambda > 0\)} \end{aligned}$$

2 *,*  *and α* < *α*ˆ*.*

*Then, the first co-existence point E*ˆ = *x*ˆ, (*K* − *x*<sup>ˆ</sup>)(*a* + *x*<sup>ˆ</sup>)*r mK is locally asymptotically stable.*

**Proof.** We first observe that the Jacobian matrix of model (6) evaluated at *E* ˆ is

$$f(\hat{E}) = \begin{bmatrix} \left(\frac{K-\hat{\mathfrak{x}}}{a+\hat{\mathfrak{x}}}-1\right)\frac{r\hat{\mathfrak{x}}}{K} & -\frac{m\hat{\mathfrak{x}}}{a+\hat{\mathfrak{x}}}\\ \frac{(K-\hat{\mathfrak{x}})\,amr}{(a+\hat{\mathfrak{x}})\,mK} & 0 \end{bmatrix}.\tag{12}$$

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

$$
\lambda^2 - \left(\frac{K-\mathfrak{k}}{a+\mathfrak{k}} - 1\right) \frac{r\mathfrak{k}}{K} \lambda + \frac{(K-\mathfrak{k})\operatorname{amr}\mathfrak{k}}{(a+\mathfrak{k})^2 K} = 0,
$$

namely

$$\begin{split} \lambda\_1 &= \left( \frac{K - \pounds}{a + \pounds} - 1 \right) \frac{r\p}{2K} + \frac{i\sqrt{\Delta}}{2}, \\ \lambda\_2 &= \left( \frac{K - \pounds}{a + \pounds} - 1 \right) \frac{r\p}{2K} - \frac{i\sqrt{\Delta}}{2}. \end{split} \tag{13}$$

When *x*ˆ > *K* − *a* 2 , the real parts of *λ*1,2 are negative. Thus, the eigenvalues (13) always satisfy |arg(*<sup>λ</sup>*1,2)| > *απ*/2 for any Δ. If *x*ˆ < *K* − *a* 2 and Δ ≤ 0, then *λ*1*λ*2 = (*K* − *x*ˆ) *anrx*<sup>ˆ</sup> (*a* + *x*<sup>ˆ</sup>)<sup>2</sup>*<sup>K</sup>* > 0 and *λ*1 + *λ*2 = *K* − *x*ˆ *a* + *x*ˆ − 1 *rx*<sup>ˆ</sup>*K* > 0, meaning that |arg(*<sup>λ</sup>*1,2)| = 0 < *απ*/2. If *x*ˆ < *K* − *a* 2 and Δ > 0, then |arg(*<sup>λ</sup>*1,2)| > *απ*/2 for *α* < *α*ˆ. Hence, we prove the theorem.

#### **Theorem 7.** *Let*

$$\begin{aligned} \xi &= \frac{((y^\*-T)^2 - cT)h}{y^\*(c+y^\*-T)^2} + \left(\frac{my^\*}{(a+x^\*)^2} - \frac{r}{K}\right)x^\*,\\ \theta &= \left(\frac{my^\*}{(a+x^\*)^2} - \frac{r}{K}\right)\frac{((y^\*-T)^2 - cT)h}{y^\*(c+y^\*-T)^2}x^\* + \left(1 - \frac{x^\*}{a+x}\right)\frac{mmx^\*y^\*}{(a+x^\*)^2}.\end{aligned}$$

*If one of the following statements is satisfied, then the second co-existence point E*∗ = (*x*<sup>∗</sup>, *y*<sup>∗</sup>) *is locally asymptotically stable:*


**Proof.** The Jacobian matrix of model (6) evaluated at *E*∗ is given by

$$J(E^{\*}) = \begin{bmatrix} \left(\frac{my^{\*}}{(a+\mathbf{x}^{\*})^{2}} - \frac{r}{K}\right)\mathbf{x}^{\*} & -\frac{mx^{\*}}{a+\mathbf{x}^{\*}}\\ \left(1 - \frac{\mathbf{x}^{\*}}{a+\mathbf{x}^{\*}}\right)\frac{ny^{\*}}{a+\mathbf{x}^{\*}} & \frac{((y^{\*}-T)^{2}-cT)h}{y^{\*}(c+y^{\*}-T)^{2}} \end{bmatrix}.\tag{14}$$

The eigenvalues of (14) are obtained by solving the characteristic equation *λ*<sup>2</sup> − *ξλ* + *θ* = 0. Thus, we have *λ*1,2 = *ξ*2 ± :*ξ*<sup>2</sup> − 4*θ* 2 . If *θ* > 0 and *ξ* < 0, then |arg(*<sup>λ</sup>*1,2)| > *απ*/2. If *ξ*2 < 4*θ* and *ξ* > 0, then |arg(*<sup>λ</sup>*1,2)| > *απ*/2 for *α* < *<sup>α</sup>*<sup>∗</sup>. 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*(*t*) ∈ *C* (R+)*, x*<sup>∗</sup> ∈ R+*, and its Caputo fractional derivative of order-α exists for any α* ∈ (0, 1]*. Then, for any t* > 0*, we have <sup>C</sup>*D*<sup>α</sup>t* .*x*(*t*) − *x*<sup>∗</sup> − *x*<sup>∗</sup> ln *x*(*t*) *x*<sup>∗</sup> / ≤ 1 − *x*<sup>∗</sup> *x*(*t*) *<sup>C</sup>*D*<sup>α</sup>t <sup>x</sup>*(*t*)*.*

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

$$\prescript{\mathsf{C}}{}{\mathcal{D}}\_{t}^{\mathsf{u}}\mathbf{x}(t) = f(\mathbf{x}(t)),\tag{15}$$

*which starts from a point in* Ω *remains in* Ω *for all time. Let <sup>V</sup>*(*x*) : Ω → R *be a continuously differentiable function such that*

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

*Let E* := &*x*|*<sup>C</sup>*D*<sup>α</sup>t <sup>V</sup>*|(15) = 0' *and M be the largest invariant set of E. Then, every solution x*(*t*) *originating in* Ω *tends to M as t* → ∞*.*

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 Ω1 := {(*<sup>x</sup>*, *y*) : *x* ≥ 0, *y* < *T*} and Ω2 := {(*<sup>x</sup>*, *y*) : *x* ≥ 0, *y* ≥ *<sup>T</sup>*}. Therefore, the global stabilities of *E*1, *E*ˆ, and *E*∗ are investigated as follows.

**Theorem 8.** *If n* < *adK , then the predator extinction point E*1 = (*<sup>K</sup>*, 0) *is globally asymptotically stable in the region* Ω1*.*

**Proof.** By specifying a positive Lyapunov function L1(*<sup>x</sup>*, *y*) = *x* − *K* − *K* ln *xK* + *myn* , and conforming to Lemma 4, we obtain

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

Owing to the fact that *n* < *adK* , we have *<sup>C</sup>*D*<sup>α</sup>t* L1(*<sup>x</sup>*, *y*) ≤ 0. In consequence of Lemma 5, *E*1 is globally asymptotically stable in the region Ω1.

**Remark 2.** *Notice E*1 *is locally asymptotically stable if n* < *adK* + *d and is globally asymptotically stable if n* < *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* (*K* − *x*<sup>ˆ</sup>)(*a* + *x*ˆ) < *a*2 *, then the first co-existence point E*ˆ = *x*ˆ, (*K* − *x*<sup>ˆ</sup>)(*a* + *x*<sup>ˆ</sup>)*r mK is globally asymptotically stable in the region* Ω1*.*

**Proof.** Let *E* ˆ = (*x*ˆ, *y*ˆ) where *y*ˆ = (*K* − *x*<sup>ˆ</sup>)(*a* + *x*<sup>ˆ</sup>)*r mK* and *ϕ* is a positive constant. By considering a Lyapunov function L2(*<sup>x</sup>*, *y*) = *x* − *x*ˆ − *x*ˆ ln *xx*ˆ <sup>+</sup> *ϕ* .*y* − *y*ˆ − *y*ˆ ln *yy*ˆ /, and using Lemma 4, we ge<sup>t</sup>

*<sup>C</sup>*D*<sup>α</sup>t* L2(*<sup>x</sup>*, *y*) ≤ *x* − *x*ˆ *x <sup>C</sup>*D*<sup>α</sup>t x* + *ϕ y* − *y*ˆ *y <sup>C</sup>*D*<sup>α</sup>t y* =(*x* − *x*ˆ) *r* − *rxK* − *my a* + *x* + *ϕ*(*y* − *y*ˆ) *nx a* + *x* − *d* =(*x* − *x*ˆ) *rx*<sup>ˆ</sup>*K* + *my*<sup>ˆ</sup> *a* + *x*ˆ − *rxK* − *my a* + *x* + *ϕ*(*y* − *y*ˆ) *nx a* + *x* − *nx*<sup>ˆ</sup> *a* + *x*ˆ = − (*x* − *x*ˆ) (*x* − *x*ˆ) *rK* + *y a* + *x* − *y*ˆ *a* + *x*ˆ *m* + *nϕ*(*y* − *y*ˆ) *x a* + *x* − *x*ˆ *a* + *x*ˆ = − (*x* − *x*<sup>ˆ</sup>)<sup>2</sup> *rK* − (*x* − *x*ˆ) (*y* − *y*<sup>ˆ</sup>)(*a* + *x*<sup>ˆ</sup>)+(*x*<sup>ˆ</sup> − *x*)*y*<sup>ˆ</sup> (*a* + *x*)(*a* + *x*ˆ) *m* + (*x* − *x*<sup>ˆ</sup>)(*y* − *y*ˆ) (*a* + *x*)(*a* + *x*ˆ) *anϕ* = − (*x* − *x*<sup>ˆ</sup>)<sup>2</sup> *rK* − (*x* − *x*<sup>ˆ</sup>)(*y* − *y*<sup>ˆ</sup>)(*a* + *x*ˆ) (*a* + *x*)(*a* + *x*ˆ) *m* + (*x* − *x*<sup>ˆ</sup>)<sup>2</sup>*y*<sup>ˆ</sup> (*a* + *x*)(*a* + *x*ˆ) *m*+ (*x* − *x*<sup>ˆ</sup>)(*y* − *y*ˆ) (*a* + *x*)(*a* + *x*ˆ) *anϕ* ≤ − (*x* − *x*<sup>ˆ</sup>)<sup>2</sup> *rK* − *my*<sup>ˆ</sup> *a*2 + (*x* − *x*<sup>ˆ</sup>)(*y* − *y*ˆ) (*a* + *x*)(*a* + *x*ˆ) (*anϕ* − (*a* + *<sup>x</sup>*<sup>ˆ</sup>)*m*).

By choosing *ϕ* = (*a* + *x*<sup>ˆ</sup>)*m an* and substituting the value of *y*ˆ, we ge<sup>t</sup>

$${}^{C}\mathcal{D}\_{1}^{a}\mathcal{L}\_{2}(\mathfrak{x},\mathfrak{y}) \leq -\left(\mathfrak{x}-\mathfrak{x}\right)^{2}\left(a^{2}-\left(K-\mathfrak{x}\right)\left(a+\mathfrak{x}\right)\right)\frac{r}{a^{2}K}.$$

Therefore, if (*K* − *x*<sup>ˆ</sup>)(*a* + *x*ˆ) < *a*2, then *<sup>C</sup>*D*<sup>α</sup>t* L2(*<sup>x</sup>*, *y*) ≤ 0. Using Lemma 5, we conclude that *E*ˆ is globally asymptotically stable in the region Ω1.

Based on Theorem 2, *x*(*t*) and *y*(*t*) are bounded. Let Ψ be the upper bound of *y*(*t*) such that 0 < *T* ≤ *y*(*t*) ≤ Ψ. The global stability of *E*∗ is stated in the following theorem.

**Theorem 10.** *If y*∗ < min *a*2*r mK*, *cT* Ψ − *T* + *T, then the second co-existence point E*∗ = (*x*<sup>∗</sup>, *y*<sup>∗</sup>) *is globally asymptotically stable in the region* Ω2*.*

**Proof.** We consider a positive Lyapunov function L3(*E*<sup>∗</sup>) = *x* − *x*<sup>∗</sup> − *x*<sup>∗</sup> ln *xx*∗ <sup>+</sup> *ϕ*∗ .*y* − *y*∗ − *y*∗ ln *yy*∗ /. According to Lemma 4, the fractional derivative of L3(*E*<sup>∗</sup>) satisfies

$$\begin{split} \, ^C \mathcal{D}\_t^\mathbf{z} \mathcal{L}\_2(\mathbf{x}, y) &\leq \left( \frac{\mathbf{x} - \mathbf{x}^\*}{\mathbf{x}} \right) \, ^C \mathcal{D}\_t^\mathbf{z} \mathbf{x} + \boldsymbol{\varrho}^\* \left( \frac{\mathbf{y} - \mathbf{y}^\*}{\mathbf{y}} \right) \, ^C \mathcal{D}\_t^\mathbf{x} \boldsymbol{y} \\ &= (\mathbf{x} - \mathbf{x}^\*) \left( r - \frac{r\mathbf{x}}{K} - \frac{my}{a + \mathbf{x}} \right) + \boldsymbol{\varrho}^\* (\mathbf{y} - \mathbf{y}^\*) \left( \frac{n\mathbf{x}}{a + \mathbf{x}} - d - \frac{h(\mathbf{y} - T)}{(c + \mathbf{y} - T)\mathbf{y}} \right) \, ^C \mathcal{D}\_t^\mathbf{x} \boldsymbol{y} \\ &= (\mathbf{x} - \mathbf{x}^\*) \left( \frac{rx^\*}{K} + \frac{my^\*}{a + \mathbf{x}^\*} - \frac{rx}{K} - \frac{my}{a + \mathbf{x}} \right) \\ &\quad + \boldsymbol{\varrho}^\* (\mathbf{y} - \mathbf{y}^\*) \left( \frac{nx}{a + \mathbf{x}} - \frac{nx^\*}{a + \mathbf{x}^\*} + \frac{h(\mathbf{y}^\* - T)}{(c + \mathbf{y}^\* - T)\mathbf{y}^\*} - \frac{h(\mathbf{y} - T)}{(c + \mathbf{y} - T)\mathbf{y}} \right) \end{split}$$

$$\begin{split} &= -\left(x-x^{\*}\right)^{2}\frac{r}{K} - \left(x-x^{\*}\right)\left(\frac{(y-y^{\*})(a+x^{\*})m - (x-x^{\*})my^{\*}}{(a+x^{\*})(a+x)}\right) \\ &\quad + \begin{split} &\quad + \mathfrak{g}^{\*}(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) \\ &= -\left(x-x^{\*}\right)^{2}\frac{r}{K} - \frac{(x-x^{\*})(y-y^{\*})(a+x^{\*})m}{(a+x^{\*})(a+x)} + \frac{(x-x^{\*})^{2}my^{\*}}{(a+x^{\*})(a+x)} \\ &\quad + \frac{(x-x^{\*})(y-y^{\*})anq^{\*}}{(a+x^{\*})(a+x)} - \frac{(cT-(y^{\*}-T)(y-T))(y-y^{\*})^{2}hq^{\*}}{(c+y^{\*}-T)(c+y-T)y^{\*}y} \\ &\leq -\left(x-x^{\*}\right)^{2}\frac{r}{K} - \frac{(x-x^{\*})(y-y^{\*})(a+x^{\*})m}{(a+x^{\*})(a+x)} + (x-x^{\*})^{2}\frac{my^{\*}}{a^{2}} \\ &\quad + \frac{(x-x^{\*})(y-y^{\*})anq^{\*}}{(a+x^{\*})(a+x)} - \frac{(cT-(y^{\*}-T)(y-T))(y-y^{\*})^{2}hq^{\*}}{(c+y^{\*}-T)(c+y-T)y^{\*}y} \end{split}$$

By taking *ϕ*∗ = (*a* + *x*<sup>∗</sup>)*m an* and remembering that *y*(*t*) < Ψ, we obtain

$$\prescript{C}{}{\mathcal{D}}\_{t}^{a}\mathcal{L}\_{2}(\mathcal{E}^{\*}) \leq - \left(\mathbf{x} - \mathbf{x}^{\*}\right)^{2} \left(\frac{r}{K} - \frac{my^{\*}}{a^{2}}\right) - \frac{(cT - (y^{\*} - T)(\Psi - T))(y - y^{\*})^{2}(a + \mathbf{x}^{\*})mh}{(c + y^{\*} - T)(c + y - T)any^{\*}y}.$$

It is easily confirmed that, if *y*∗ < min *a*2*r mK*, *cT* Ψ − *T* + *T*, then *<sup>C</sup>*D*<sup>α</sup>t* L2(*<sup>x</sup>*, *y*) ≤ 0. Based on Lemma 5, *E* ˆ is globally asymptotically stable in the region Ω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 *λ*1,2 of the Jacobian matrix evaluated at the equilibrium point satisfy the following conditions:


$$\left(\overleftrightarrow{\mathrm{iii}}\right)\_{\alpha} \qquad \frac{d\dot{m}(\vec{a})}{d\alpha}\Big|\_{\alpha=\alpha^\*} \neq 0.$$

Now, consider the stability condition in Theorems 6 and 7. For *y* < *T*, the Jacobian matrix of model (6) evaluated at *E* ˆ has a pair of complex eigenvalues if Δ > 0. We can easily confirm that, if *x*ˆ < *K* − *a* 2 , then the real part of the eigenvalues are positive. We also have that *m*(*α*<sup>ˆ</sup>) = 0 and *dm*(*α*) *dα α*=*α*ˆ = 0. Therefore, *E*ˆ undergoes a Hopf bifurcation when *α* crosses *α*ˆ. A similar circumstance also occurs when *y* ≥ *T*. When *ξ*2 < 4*θ*, the Jacobian matrix of model (6) has a pair of complex eigenvalues. The real part of the eigenvalues are also positive when *ξ* > 0. We can also check that *m*(*α*<sup>∗</sup>) = 0 and *dm*(*α*) *dα α*=*α*∗ = 0. This means the Hopf bifurcation also occurs when *y* ≥ *T*. Therefore, we have the following theorem.

**Theorem 11.** *(i) Let* Δ > 0 *and x*ˆ < *K* − *a* 2 *. The first co-existence point E*ˆ *undergoes a Hopf bifurcation when α passes through α*ˆ *in the region* Ω1*.*

*(ii) Let ξ*2 < 4*θ and ξ* > 0*. The second co-existence point E*∗ *undergoes a Hopf bifurcation when α passes through α*<sup>∗</sup> *in the region* Ω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-*α* which is defined as follows.

**Definition 3.** *[40] Suppose* 0 < *α* ≤ 1*. The Atangana–Baleanu fractional integral and derivative in Caputo sense of order-α (ABC) is defined by*

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

*where t* ≥ 0, *f* ∈ *Cn*([0, <sup>+</sup>∞), <sup>R</sup>)*, <sup>E</sup>α*(*t*) = ∑∞*<sup>k</sup>*=<sup>0</sup> *tk* <sup>Γ</sup>(*αk* + 1) *is the Mittag–Leffler function and <sup>B</sup>*(*α*) *is a normalization function with B*(0) = *B*(1) = 1*. In this paper, we define <sup>B</sup>*(*α*) = 1 − *α* + *α* <sup>Γ</sup>(*α*)*.*

By applying Definition 3 to model (4), we ge<sup>t</sup> the following fractional-order model with ABC operator

$$\begin{aligned} \, ^{ABC}\mathcal{D}\_t^a \mathbf{x} &= r \mathbf{x} \left( 1 - \frac{\mathbf{x}}{K} \right) - \frac{m \mathbf{x} \mathbf{y}}{a + \mathbf{x}} \equiv \mathcal{G}\_1, \\ ^{ABC}\mathcal{D}\_t^a \mathbf{y} &= \frac{n \mathbf{x} \mathbf{y}}{a + \mathbf{x}} - dy - H(\mathbf{y}) \equiv \mathcal{G}\_2. \end{aligned} \tag{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* ≥ *T*.

We start for *y* ≥ *T*. Let *x*, *x*¯, *y*, and *y*¯ be functions satisfying *x* ≤ *a*1, *x*¯ ≤ *a*2, *y* ≤ *b*1, and *y*¯ ≤ *b*2. Suppose that

$$\begin{aligned} \mathbf{g}\_1 &= r + (a\_1 + a\_2)\frac{r}{K} + \frac{my}{a}, \\ \mathbf{g}\_2 &= n + d + \frac{h}{c}. \end{aligned}$$

Therefore, we ge<sup>t</sup>

*<sup>G</sup>*1(*x*) − *<sup>G</sup>*1(*x*¯) = 9 9 9 9 *rx* 1 − *x K* − *mxy a* + *x* − *rx*¯ 1 − *x*¯ *K* − *mxy*¯ *a* + *x*¯ 9 9 9 9 = 9 9 9 9 *rx* − *rx*<sup>2</sup> *K* − *mxy a* + *x* − *rx*¯ + *rx*¯2 *K* + *mxy*¯ *a* + *x*¯ 9 9 9 9 = 9 9 9 9*<sup>r</sup>*(*<sup>x</sup>* − *x*¯) − *r K*(*x*<sup>2</sup> − *x*¯ 2) − *my x a* + *x* − *x*¯ *a* + *x*¯ 9 9 9 9 = 9 9 9 9*<sup>r</sup>*(*<sup>x</sup>* − *x*¯) − *r K*(*<sup>x</sup>* + *x*¯)(*x* − *x*¯) − *amy*(*x* − *x*¯) (*a* + *x*)(*a* + *x*¯) 9 9 9 9 ≤*r x* − *x*¯ + (*<sup>a</sup>*1 + *<sup>a</sup>*2) *r K x* − *x*¯ + *my a x* − *x*¯ = *r* + (*<sup>a</sup>*1 + *<sup>a</sup>*2) *r K* + *my a x* − *x*¯ =*g*1 *x* − *x*¯ , (17)

and

$$\begin{split} \left\|(G\_{2}(y)-G\_{2}(\bar{y}))\right\| &= \left\| \left(\frac{nxy}{a+x}-dy-\frac{h(y-T)}{c+(y-T)}\right)-\left(\frac{nx\bar{y}}{a+x}-d\bar{y}-\frac{h(\bar{y}-T)}{c+(\bar{y}-T)}\right) \right\| \\ &= \left\|\frac{nxy}{a+x}-dy-\frac{h(y-T)}{c+(y-T)}-\frac{nx\bar{y}}{a+x}+d\bar{y}+\frac{h(\bar{y}-T)}{c+(\bar{y}-T)}\right\| \\ &= \left\|\frac{nx(y-\bar{y})}{a+x}-d(y-\bar{y})-\frac{ch(y-\bar{y})}{(c+y-T)(c+\bar{y}-T)}\right\| \\ &\leq n\left\|y-\bar{y}\right\|+d\left\|y-\bar{y}\right\|+\frac{h}{c}\left\|y-\bar{y}\right\| \\ &= \left(n+d+\frac{h}{c}\right)\left\|y-\bar{y}\right\| \end{split} \tag{18}$$

When *y* < *T*, by utilizing the similar manner, we achieve *<sup>G</sup>*2(*y*) − *<sup>G</sup>*2(*y*¯) ≤ *g*3 *y* − *y*¯ where *g*3 = *n* + *d*. Accordingly, the kernel of (16) satisfies the Lipschitz condition. Furthermore, if 0 ≤ *g*1 < 1 and 0 ≤ *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{split} \mathbf{x}(t) - \mathbf{x}(0) &= \frac{1-a}{B(a)} \mathbf{G}\_1(t, \mathbf{x}) + \frac{a}{B(a)\Gamma(a)} \int\_0^t (t-s)^{a-1} \mathbf{G}\_1(s, \mathbf{x}) ds, \\ \mathbf{y}(t) - \mathbf{y}(0) &= \frac{1-a}{B(a)} \mathbf{G}\_2(t, \mathbf{y}) + \frac{a}{B(a)\Gamma(a)} \int\_0^t (t-s)^{a-1} \mathbf{G}\_2(s, \mathbf{y}) ds. \end{split} \tag{19}$$

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

$$\begin{split} \mathbf{x}\_{n}(t) &= \frac{1-\mathfrak{a}}{\mathcal{B}(a)} \mathbf{G}\_{1}(t, \mathbf{x}\_{n-1}) + \frac{\mathfrak{a}}{\mathcal{B}(a)\Gamma(a)} \int\_{0}^{t} (t-s)^{a-1} \mathbf{G}\_{1}(s, \mathbf{x}\_{n-1}) ds, \\ \mathbf{y}\_{n}(t) &= \frac{1-\mathfrak{a}}{\mathcal{B}(a)} \mathbf{G}\_{2}(t, \mathbf{y}\_{n-1}) + \frac{\mathfrak{a}}{\mathcal{B}(a)\Gamma(a)} \int\_{0}^{t} (t-s)^{a-1} \mathbf{G}\_{2}(s, \mathbf{y}\_{n-1}) ds. \end{split} \tag{20}$$

The associated initial conditions along with Equation (20) are *<sup>x</sup>*0(*t*) = *x*(0) and *y*0(*t*) = *y*(0). By considering Equation (20), we have the difference expression of successive terms as follows.

$$\begin{split} \boldsymbol{\Phi}\_{1,n}(t) &= \mathbf{x}\_{n}(t) - \mathbf{x}\_{n-1}(t) \\ &= \frac{1-\mathfrak{a}}{\mathcal{B}(a)} \left( \mathbf{G}\_{1}(t, \mathbf{x}\_{n-1}) - \mathbf{G}\_{1}(t, \mathbf{x}\_{n-2}) \right) + \frac{\mathfrak{a}}{\mathcal{B}(a)\Gamma(a)} \int\_{0}^{t} (t-s)^{\mathfrak{a}-1} \left( \mathbf{G}\_{1}(t, \mathbf{x}\_{n-1}) - \mathbf{G}\_{1}(t, \mathbf{x}\_{n-2}) \right) ds, \\ \boldsymbol{\Phi}\_{2,n}(t) &= \boldsymbol{y}\_{n}(t) - \mathbf{y}\_{n-1}(t) \\ &= \frac{1-\mathfrak{a}}{\mathcal{B}(a)} \left( \mathbf{G}\_{2}(t, \mathbf{y}\_{n-1}) - \mathbf{G}\_{2}(t, \mathbf{y}\_{n-2}) \right) + \frac{\mathfrak{a}}{\mathcal{B}(a)\Gamma(a)} \int\_{0}^{t} (t-s)^{\mathfrak{a}-1} \left( \mathbf{G}\_{2}(t, \mathbf{y}\_{n-1}) - \mathbf{G}\_{2}(t, \mathbf{y}\_{n-2}) \right) ds. \end{split} \tag{21}$$

According to Equation (21), we ge<sup>t</sup> *xn*(*t*) = ∑*ni*=<sup>1</sup> *<sup>Φ</sup>*1,*i*(*t*) and *yn*(*t*) = ∑*ni*=<sup>1</sup> *<sup>Φ</sup>*2,*<sup>i</sup>*(*t*). By applying Equations (17), (18) and (21), we have that

$$\begin{split} \|\|\Phi\_{1,n}(t)\|\| \leq & \frac{1-a}{B(a)} \mathfrak{g}\_{1} \|\|\Phi\_{1,n-1}\|\| + \frac{a}{B(a)\Gamma(a)} \mathfrak{g}\_{1} \int\_{0}^{t} \|\Phi\_{1,n-1}(s)\|\| (t-s)^{a-1} ds, \\ \|\|\Phi\_{2,n}(t)\|\| \leq & \frac{1-a}{B(a)} \mathfrak{g}\_{2} \|\|\Phi\_{2,n-1}\|\| + \frac{a}{B(a)\Gamma(a)} \mathfrak{g}\_{2} \int\_{0}^{t} \|\Phi\_{2,n-1}(s)\|\| (t-s)^{a-1} ds. \end{split} \tag{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-a)\mathcal{g}\_i}{B(a)} + \frac{t\_0^a \mathcal{g}\_i}{B(a)\Gamma(a)} < 1, i = 1, 2, 3. \tag{23}$$

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

$$\begin{split} \|\|\Phi\_{1,n}(t)\|\| \leq & \|\ge\_{0}\|\left(\frac{(1-a)\mathcal{G}\_{1}}{B(a)} + \frac{t^{a}\mathcal{G}\_{1}}{B(a)\Gamma(a)}\right)^{n}\| \\ \|\|\Phi\_{2,n}(t)\|\| \leq & \|\|y\_{0}\|\| \left(\frac{(1-a)\mathcal{G}\_{2}}{B(a)} + \frac{t^{a}\mathcal{G}\_{2}}{B(a)\Gamma(a)}\right)^{n} . \end{split} \tag{24}$$

Therefore, the continuity and existence of solution are proved since *<sup>Φ</sup>*1,*n*(*t*) → 0 and *<sup>Φ</sup>*2,*n*(*t*) → 0 as *n* → ∞ and *t* = *t*0. Now, suppose that

$$\begin{aligned} \mathbf{x}(t) - \mathbf{x}(0) &= \mathbf{x}\_n(t) - \mathbf{Y}\_{1,n}(t), \\ \mathbf{y}(t) - \mathbf{y}(0) &= \mathbf{y}\_n(t) - \mathbf{Y}\_{2,n}(t). \end{aligned} \tag{25}$$

We confirm that

$$\|\|Y\_{1,u}(t)\|\| \le \left(\frac{(1-a)}{B(a)} + \frac{t^a}{B(a)\Gamma(a)}\right)^{n+1} g\_1^{n+1} \tag{26}$$

It is clear that <sup>Υ</sup>1,*n*(*t*) → 0 when *n* → ∞. By using a similar manner, we acquire <sup>Υ</sup>2,*n*(*t*) → 0, *n* → ∞. Finally, the uniqueness of solution for the model is proven. Suppose that there exists different set of solutions denote by *x*˜(*t*) and *y*˜(*t*); then, we have

$$\mathbf{x}(t) - \tilde{\mathbf{x}}(t) = \frac{1 - a}{B(a)} (\mathbf{G}\_1(t, \mathbf{x}) - \mathbf{G}\_1(t, \tilde{\mathbf{x}})) + \frac{a}{B(a)\Gamma(a)} \int\_0^t (\mathbf{G}\_1(s, \mathbf{x}) - \mathbf{G}\_1(s, \tilde{\mathbf{x}}))(t - s)^{a - 1} ds. \tag{27}$$

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

$$\|\|\mathbf{x}(t) - \tilde{\mathbf{x}}(t)\|\left(1 - \frac{(1 - a)g\_1}{B(a)} - \frac{t^a g\_1}{B(a)\Gamma(a)}\right) \le 0. \tag{28}$$

For *t* = *t*0, we have (23), thus *x*(*t*) − *x*˜(*t*) = 0 and hence *x*(*t*) = *<sup>x</sup>*˜(*t*). Applying the same algebraic procedures, we can show that *y*(*t*) = *y*˜(*t*). 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–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.

$$
\sigma = 0.5, K = 1, m = 0.3, a = 0.2, d = 0.1, h = 0.1, T = 0.9, c = 0.1, \text{ and } a = 0.9. \tag{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 ≤ *n* ≤ 0.2 lead to the non-existence of equilibrium point in Ω2. Therefore, the first numerical simulations are focused on the dynamics in Ω1.

For 0.08 ≤ *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 Ω1 and Ω2 are convergen<sup>t</sup> to *E*1. When the initial value is in Ω2, then the solution is oscillating when it crosses the threshold harvesting level and eventually converges to *E*1.

When *n* passes through *<sup>n</sup>*<sup>∗</sup>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 *E* ˆ . Figure 1 shows that *E* ˆ is asymptotically stable if 0.12 < *n* - *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 *E* ˆ = (0.5, 0.5833) is asymptotically stable. This circumstance corresponds to Theorems 4–6 and 9.

**Figure 1.** Bifurcation diagram of Caputo model (6) and ABC model (16) driven by the conversion rate of consumed prey into predator birth (*n*) with constant parameter values (29). There exists two bifurcations namely a forward bifurcation which occurs when *n* passes through *n*∗1 ≈ 0.12, and a Hopf bifurcation which occurs when *n* passes through *n*∗2≈ 0.1557.

**Figure 2.** Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29) and *n* = 0.1: Figure (**a**) phase portrait; and Figure (**b**) time series.

Furthermore, if we increase the value of *n* such that *n* > *<sup>n</sup>*<sup>∗</sup>2, then the equilibrium point *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 *E* ˆ in Ω1, and *E*∗ = (0.5954, 0.5364) occurs in Ω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.

**Figure 3.** Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29) and *n* = 0.14: Figure (**a**) phase portrait; and Figure (**b**) time series.

Notice that, in Figures 2–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 ge<sup>t</sup> more detail view, we plot a bifurcation diagram by varying the order of the fractional derivative (*α*) (see Figure 6). In this simulation, we use parameter values as in Figure 4 and vary the order-*α* 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 *α* (see Figure 7). For *α* = 0.7, the equilibrium point *E* ˆ of both model are asymptotically stable. For *α* = 0.772, the equilibrium point *E* ˆ of the Caputo model loses its stability, while that of the ABC model remains asymtotically stable. For *α* = 0.83, the equilibrium point *E* ˆ of both models are unstable.

**Figure 4.** Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29) and *n* = 0.2: Figure (**a**) phase portrait; and Figure (**b**) time series.

(**a**) Phase portrait (**b**) *n* = 0.14 **Figure 5.** Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29),*n*=0.2and*T*=0.5.

**Figure 6.** Bifurcation diagram of Caputo model (6) and ABC model (16) driven by the order of the fractional-derivative (*α*) with constant parameter values (29) and *n* = 0.2. There exists a Hopf bifurcation where the bifurcation points of the Caputo model and ABC model are different.

**Figure 7.** *Cont.*

**Figure 7.** Numerical simulations of Caputo model (6) and ABC model (16) with parameter values (29), *n* = 0.2, and *α* = {0.7, 0.772, 0.83}: Figure (**a**) phase portrait; and Figure (**b**) time series.

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 *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.
