Local and global dynamics of a fractional-order predator-prey system with habitat complexity and the corresponding discretized fractional-order system

TL;DR

This study investigates the stability and complex dynamics of a fractional-order predator-prey model with habitat complexity, revealing conditions for stability, bifurcations, and the effects of discretization on system behavior.

Contribution

It introduces a fractional-order predator-prey model with habitat complexity and analyzes its stability, bifurcations, and discretized dynamics, highlighting fractional-order-dependent phenomena.

Findings

Existence of positivity and boundedness of solutions.

Conditions for local and global stability of equilibrium points.

Presence of Hopf bifurcation and complex dynamics depending on fractional order and step size.

Abstract

This paper is focused on local and global stability of a fractional-order predator-prey model with habitat complexity constructed in the Caputo sense and corresponding discrete fractional-order system. Mathematical results like positivity and boundedness of the solutions in fractional-order model is presented. Conditions for local and global stability of different equilibrium points are proved. It is shown that there may exist fractional-order-dependent instability through Hopf bifurcation for both fractional-order and corresponding discrete systems. Dynamics of the discrete fractional-order model is more complex and depends on both step length and fractional-order. It shows Hopf bifurcation, flip bifurcation and more complex dynamics with respect to the step size. Several examples are presented to substantiate the analytical results.

Equations

$$\frac{dx}{dt}$$

$$\frac{dy}{dt}$$

$${}_0^c{D}_t^{m}x$$

$${}_0^c{D}_t^{m}y$$

$$f(t)=f(a)+\frac{1}{\Gamma (m)}(D_a^{m}f)(\xi ).(t-a{)}^m,$$

$${}_{t_0}^c{D}_t^{m}u(t)\le -\lambda u(t)+\mu ,$$

$$u(t_0)=u_{t_0},$$

$$u(t)\leq\bigg{(}u_{t_{0}}-\frac{\mu}{\lambda}\bigg{)}E_{m}[-\lambda(t-t_{0})^{m}]+\frac{\mu}{\lambda}.$$

$$x(t)$$

$${}_0^c{D}_t^{m}x(t){|}_{t=t_1}=0.$$

$$V(t)=x+\frac{1}{\theta }y,$$

$${}^{c}_{0}D^{m}_{t}V(t)=~{}^{c}_{0}D^{m}_{t}x(t)+~{}^{c}_{0}D^{m}_{t}\frac{1}{\theta}y(t)~{}=~{}rx\bigg{(}1-\frac{x}{K}\bigg{)}-\frac{d}{\theta}y.$$

$${}_0^c{D}_t^{m}V(t)+\eta V(t)=$$

$$=$$

$$\le$$

$${}_0^c{D}_t^{m}V(t)+\eta V(t)\le l,$$

$$V(t)$$

$${}_{t_0}^c{D}_t^{m}x(t)=f(t,x),t>t_0$$

$$H_{1}(X)=rx\bigg{(}1-\frac{x}{K}\bigg{)}-\frac{\alpha(1-c)xy}{1+\alpha(1-c)hx},~{}H_{2}(X)=\frac{\theta\alpha(1-c)xy}{1+\alpha(1-c)hx}-dy.$$

$$\displaystyle\begin{split}\parallel H(X)-H(\bar{X})\parallel=\mid H_{1}(X)-H_{1}(\bar{X})\mid+\mid H_{2}(X)-H_{2}(\bar{X})\mid\\ =\mid rx\bigg{(}1-\frac{x}{K}\bigg{)}-\frac{\alpha(1-c)xy}{1+\alpha(1-c)hx}-r\bar{x}\bigg{(}1-\frac{\bar{x}}{K}\bigg{)}+\frac{\alpha(1-c)\bar{x}\bar{y}}{1+\alpha(1-c)h\bar{x}}\mid\\ +\mid\frac{\theta\alpha(1-c)xy}{1+\alpha(1-c)hx}-dy-\frac{\theta\alpha(1-c)\bar{x}\bar{y}}{1+\alpha(1-c)h\bar{x}}+d\bar{y}\mid\\ =\mid r(x-\bar{x})-\frac{r}{K}(x^{2}-\bar{x}^{2})-\alpha(1-c)\bigg{(}\frac{xy}{1+\alpha(1-c)hx}-\frac{\bar{x}\bar{y}}{1+\alpha(1-c)h\bar{x}}\bigg{)}\mid\\ +\mid\theta\alpha(1-c)\bigg{(}\frac{xy}{1+\alpha(1-c)hx}-\frac{\bar{x}\bar{y}}{1+\alpha(1-c)h\bar{x}}\bigg{)}-d(y-\bar{y})\mid\\ \leq r\mid x-\bar{x}\mid+\frac{r}{K}\mid x^{2}-\bar{x}^{2}\mid\\ +\alpha(1-c)(1+\theta)\mid\bigg{(}\frac{xy}{1+\alpha(1-c)hx}-\frac{\bar{x}\bar{y}}{1+\alpha(1-c)h\bar{x}}\bigg{)}\mid+d\mid y-\bar{y}\mid\\ \leq r\mid x-\bar{x}\mid+\frac{2rM}{K}\mid x-\bar{x}\mid+\alpha(1-c)(1+\theta)\mid xy-\bar{x}\bar{y}\mid\\ +\alpha^{2}(1-c)^{2}(1+\theta)hM^{2}\mid y-\bar{y}\mid+d\mid y-\bar{y}\mid\\ \leq\bigg{(}r+\frac{2rM}{K}+\alpha(1-c)(1+\theta)M\bigg{)}\mid x-\bar{x}\mid\\ +\bigg{(}\alpha(1-c)(1+\theta)M+d+\alpha^{2}(1-c)^{2}(1+\theta)hM^{2}\bigg{)}\mid y-\bar{y}\mid\\ \leq L\parallel(x,y)-(\bar{x},\bar{y})\parallel\\ \leq L\parallel X-\bar{X}\parallel,\end{split}$$

$${}_0^c{D}_t^{m}x(t)=f(x),x(0)=x_0$$

$$\mid arg({\lambda }_i)\mid >\frac{m\pi }{2},i=1,2,.....,n.$$

$$\mathbf{D}(\varphi )=-\left|\begin{array}{ccc}1& a_1& a_2\\ 2& a_1& 0\\ 0& 2& a_1\end{array}\right|=a_1^2-4a_2.$$

$$x^{\ast }$$

$$\mathbf{J}({\mathbf{E}}_1)=\left(\begin{array}{cc}-r& -\frac{\alpha K(1-c)}{1+\alpha K(1-c)h}\\ 0& \frac{\theta \alpha K(1-c)}{1+\alpha K(1-c)h}-d\end{array}\right).$$

$${}^{c}_{t_{0}}D^{m}_{t}\bigg{[}x(t)-x^{*}-x^{*}ln\frac{x(t)}{x^{*}}\bigg{]}\leq\bigg{(}1-\frac{x^{*}}{x(t)}\bigg{)}~{}~{}{{}^{c}_{t_{0}}}D^{m}_{t}x(t),x^{*}\in\Re_{+},\forall m\in(0,1].$$

$$V(x,y)=\bigg{(}x-K-Kln\frac{x}{K}\bigg{)}+\frac{y}{\theta}.$$

$$\displaystyle\begin{split}{}^{c}_{0}D^{m}_{t}V(x,y)\leq&\frac{(x-K)}{x}{{}^{c}_{0}}D^{m}_{t}x(t)+\frac{1}{\theta}{{}^{c}_{0}}D^{m}_{t}y(t)\\ =&(x-K)\bigg{[}r(1-\frac{x}{K})-\frac{\alpha(1-c)y}{1+\alpha(1-c)hx}\bigg{]}+\frac{\alpha(1-c)xy}{1+\alpha(1-c)hx}-\frac{dy}{\theta}\\ =&(x-K)\bigg{[}-\frac{r}{K}(x-K)-\frac{\alpha(1-c)y}{1+\alpha(1-c)hx}\bigg{]}+\frac{\alpha(1-c)xy}{1+\alpha(1-c)hx}-\frac{dy}{\theta}\\ =&-\frac{r}{K}(x-K)^{2}+\frac{\alpha K(1-c)y}{1+\alpha(1-c)hx}-\frac{dy}{\theta}\\ \leq&-\frac{r}{K}(x-K)^{2}+\bigg{[}\alpha K(1-c)-\frac{d}{\theta}\bigg{]}y.\\ \end{split}$$

$$\mathbf{J}({\mathbf{E}}^{\ast })=\left(\begin{array}{cc}r(1-\frac{2x^{\ast }}{K})-\frac{\alpha (1-c)y^{\ast }}{(1+\alpha (1-c)h{x}^{\ast }{)}^2}& -\frac{\alpha (1-c)x^{\ast }}{(1+\alpha (1-c)h{x}^{\ast })}\\ \frac{\theta \alpha (1-c)y^{\ast }}{(1+\alpha (1-c)h{x}^{\ast }{)}^2}& 0\end{array}\right).$$

$${\xi }^2-trace(J^{\ast })\xi +det(J^{\ast })=0,$$

$${\xi }_{1,2}=\frac{1}{2}[trace(J^{\ast })\pm \sqrt{trace(J^{\ast }{)}^2-4det(J^{\ast })}].$$

$${\xi }_i=\frac{1}{2}[trace(J^{\ast })\pm i\sqrt{4det(J^{\ast })-trace(J^{\ast }{)}^2}],i=1,2,$$

Full text

††journal:

1 Introduction

Fractional calculus is the area of mathematics where derivatives and integrals can be extended to an arbitrary order. There are different approaches to study the dynamical behaviors of population models, e.g. ordinary differential equations (ODE), partial differential equations (PDE), difference equations (DE), fractional-order differential equations (FDE) etc. The first three techniques are being extensively used for a long time. However, the fractional-order differential equations have gained considerable importance only in the recent past due to their ability of providing an exact or approximate description of different nonlinear phenomena. The main advantage of fractional-order system is that they allow greater degrees of freedom than an integer order system [1]. FDE are preferably used since they are naturally related to systems with memory which exists in most biological phenomena [2]. Moreover, FDE has close relations to fractals which has wide applications in mathematical biology. Recently, some authors have investigated the importance of fractional-order differential equations in several biological systems, e.g. ecological system with delay [3, 4], control based epidemiological system [5], ecological system with diffusion [6] etc. It also has applications in other fields of science and engineering [7, 8, 9, 10, 11]. Some recent studies discuss about the approximate solution of nonlinear fractional-order differential population models [1, 12] and some others study the qualitative behavior of nonlinear interactions of biological systems [13, 14, 15, 16, 17]. However, existence and proof of Hopf bifurcation, that causes oscillations in population densities due to fractional-order, is rare in the fractional-order population model. We address this issue along with others in a fractional-order predator-prey model considered in Caputo sense. Lot of discrete models on biological systems have been proposed and analyzed. However, discretization of a fractional-order population model is rare. Elsadany and Matouk [18] recently studied a fractional-order Lotka-Volterra predator prey model with its discretization. They showed complex dynamics even in a simpler prey-predator model. In the second phase of this paper, we construct the discrete version of the continuous fractional-order system and reveal its dynamics.

All most every habitat, whether it is aquatic or terrestrial, contains some kind of complexity. For example, sea grass, aquatic weeds, salt marshes, littoral zone vegetation, mangroves, coral reefs etc. make aquatic habitat complex. Both field and laboratory experiments confirm that habitat complexity increases persistency of interacting species [19, 20, 21, 22, 23, 24]. A general hypothesis is that habitat complexity reduces predation rates by decreasing predator-prey interaction and thereby increases population persistency. A Rosenzweig-MacArthur predator-prey model [25] that incorporates the effect of habitat complexity can be represented by the following coupled nonlinear system:

[TABLE]

This model says that the prey population $x$ grows logistically with intrinsic growth rate $r$ to its carrying capacity $K$. Predator $y$ captures the prey at a maximum rate $\alpha$ in absence of any habitat complexity $(c=0)$. In presence of complexity, predation rate decreases to $\alpha(1-c)$, where the dimensionless parameter $c$ is called the degree or strength of complexity. The value of $c$ ranges from [math] to $1$. In particular, $c=0.4$ implies that predation rate decreases by $40\%$ due to habitat complexity. If $c=0$, i.e. if there is no habitat complexity then the system (1) reduces to well known Rosenzweig-MacArthur model [25]. However, if $c=1$ then $y\rightarrow 0$ as $t\rightarrow\infty$ and the prey population grows logistically to its maximum value $K$. The parameter $\theta$ $(0<\theta<1)$ is the conversion efficiency, measuring the number of newly born predators for each captured prey and $d$ is the death rate of predator. All parameters are assumed to be positive. For construction and more explanation of the model, readers are referred to [26].

Considering the fractional derivatives in the sense of Caputo, we have the following fractional-order model corresponding to the integer order model (1):

[TABLE]

where ${}^{c}_{0}D^{m}_{t}$ is the Caputo fractional derivative with fractional-order $m$ $(0<m\leq 1)$. The main advantage of Caputo’s approach is that the initial conditions for the fractional differential equations with Caputo derivatives takes the similar form as for integer-order differential equations [27, 28], and thus takes the advantage of defining integer order initial conditions for fractional-order differential equations. We analyze system (2) with the initial conditions $x(0)>0,~{}y(0)>0.$ In this paper, we prove different mathematical results, like existence, non-negativity and boundedness of the solutions of fractional-order system (2). We establish conditions for local and global stabilities of different equilibrium points. It is shown that the interior equilibrium may switch its stability through Hopf bifurcation for some critical value of the fractional-order when the degree of complexity is low. A discrete system generally produces more complex dynamics than its continuous counterpart [18]. Here we construct a discrete version of the fractional-order prey-predator model (2). We prove local stability of different fixed points of the discrete system along with the existence conditions of Hopf and flip bifurcations. Numerical examples are presented for both systems in support of the analytical results. He and Lai (2011) has discretized a continuous type predator-prey model by Euler method. Using center manifold theorem, it is shown that the system undergoes flip and Neimark-Sacker bifurcations. Period doubling bifurcation leading to chaos was also shown through numerical simulations. However, they have not studied the dynamics of fractional-order discrete system. Abdelaziz et al. (2018) transformed an integer order SI-type epidemic model to a fraction-order discrete epidemic model and analyzed it to show flip and Neimark-Sacker bifurcations. But did not analyze the qualitative behavior of the fractional-order system. Here we study both the fractional-order and discretized fractional-order predator-prey systems. We compare the qualitative behavior of integer order system with the fractional-order and discretized fractional-order systems.

The rest of the paper is organized as follows. The next section contains well-posedness, existence and uniqueness of the solutions of the fractional-order system. Qualitative behavior of different equilibrium points are also presented here. Section 3 deals with stability and hopf bifurcation of fractional-order discrete system. Different examples are presented to illustrate the observed dynamics in Section 4. The paper ends with a brief discussion in Section 5.

2 Well-posedness

2.1 Nonnegativity and boundedness

Considering the biological significance of the model, we are only interested in solutions that are nonnegative and bounded in the region $\Re^{2}{}_{+}=\{z\in\Re^{2}|z\geq 0\}$ and $z(t)=(x(t),y(t))^{T}$. To prove the nonnegativity and uniform boundedness of our system, we shall use the following results.

Lemma 2.1 [29] Suppose that $f(t)\in C[a,b]$ and $D^{m}_{a}f(t)\in C(a,b]$ with $0<m\leq 1$. The Generalized Mean Value Theorem states that

[TABLE]

*where $a\leq\xi\leq t$, $\forall t\in(a,b]$.

From this lemma, one can easily prove the following result.

Corollary 2.1 [13, 29] * Suppose $f(t)\in C[a,b]$ and ${}^{c}_{t_{0}}D^{m}_{t}f(t)\in C[a,b]$, $0<m\leq 1$. If ${}^{c}_{t_{0}}D^{m}_{t}f(t)\geq 0,\forall t\in(a,b)$ then $f(t)$ is a non decreasing function for each $t\in[a,b]$ and if ${}^{c}_{t_{0}}D^{m}_{t}f(t)\leq 0,\forall t\in(a,b)$ then $f(t)$ is a non-increasing function for each $t\in[a,b]$.

Lemma 2.2 [13] * Let $u(t)$ be a continuous function on $[t_{0},\infty)$ and satisfying*

[TABLE]

[TABLE]

where $0<m\leq 1$, $(\lambda,\mu)\in\Re^{2}$, $\lambda\neq 0$ and $t_{0}\geq 0$ is the initial time. Then its solution has the form

[TABLE]

Theorem 2.1 * All solutions of system (2) which start in $\Re^{2}_{+}$ are nonnegative and uniformly bounded.

Proof First we show that the solutions $x(t)\in\Re^{2}_{+}$ are nonnegative if it start with positive initial values. If not, then there exists a $t_{1}>0$ such that

[TABLE]

Using (2.1) in the first equation of (2), we have

[TABLE]

According to Corollary $2.1$, we have $x(t_{1}^{+})=0$, which contradicts the fact $x(t_{1}^{+})<0$. Therefore, we have $x(t)\geq 0,~{}\forall~{}t\geq 0$. Using similar arguments, we can prove $y(t)\geq 0,\forall t\geq 0$.

Next we show that all solutions of system (2) which initiate in $\Re^{2}_{+}$ are uniformly bounded. Define a function

[TABLE]

Taking fractional time derivative, we have

[TABLE]

Now, for each $\eta>0$, we have

[TABLE]

If we take $\eta<d$ then right hand side of (2.1) is bounded in $\Re^{2}_{+}$ and there exist a constant $l>0$ (say) such that

[TABLE]

where $l=\frac{K}{4r}(r+\eta)^{2}$.

Applying Lemma $2.2$, we then have

[TABLE]

For $t\rightarrow\infty$, we thus have $V(t)\rightarrow\frac{l}{\eta}$. Therefore, $0<V(t)\leq\frac{l}{\eta}$. Hence all solutions of the system (2) that starts from $\Re^{2}_{+}$ are confined in the region $B=\{(x,y)\in\Re^{2}_{+}|0<V(t)\leq\frac{l}{\eta}+\epsilon$, for any $\epsilon>0,0<\eta<d,l=\frac{K}{4r}(r+\eta)^{2}\}$. Hence the theorem.

2.2 Existence and uniqueness

Here we study the existence and uniqueness of the solution of our system (2). We have the following Lemma due to Li et al [30].

Lemma 2.3 Consider the system

[TABLE]

with initial condition $x_{t_{0}}$, where $0<m\leq 1$, $f:[t_{0},\infty)\times\Omega\rightarrow\Re^{n}$, $\Omega\in\Re^{n}$. If $f(t,x)$ satisfies the locally Lipschitz condition with respect to $x$ then there exists a unique solution of the above system on $[t_{0},\infty)\times\Omega$.

We study the existence and uniqueness of the solution of system (2) in the region $\Omega\times[0,T]$, where $\Omega=\{(x,y)\in\Re^{2}|~{}max\{|x|,|y|\}\leq M\}$, $T<\infty$ and $M$ is large. Denote $X=(x,y)$, $\bar{X}=(\bar{x},\bar{y})$. Consider a mapping $H:\Omega\rightarrow\Re^{2}$ such that $H(X)=(H_{1}(X),H_{2}(X))$, where

[TABLE]

For any $X,\bar{X}\in\Omega$, it follows from (9) that

[TABLE]

where $L=max\{r+\frac{2rM}{K}+\alpha(1-c)(1+\theta)M,\alpha(1-c)(1+\theta)M+d+\alpha^{2}(1-c)^{2}(1+\theta)hM^{2}\}$. Thus $H(X)$ satisfies Lipschitz condition with respect to $X$ and following Lemma $2.3$, there exists a unique solution $X(t)$ of system (2) with initial condition $X(0)=(x(0),y(0))$.

2.3 Stability of equilibrium points

We have the following stability result on fractional-order differential equations.

Theorem 2.2 [31] Consider the following fractional-order system

[TABLE]

with $0<m\leq 1,x\in\Re^{n}$ and $f:\Re^{n}\rightarrow\Re^{n}$. The equilibrium points of the above system are calculated by solving the equation $f(x)=0$. These equilibrium points are locally asymptotically stable if all eigenvalues $\lambda_{i}$ of the jacobian matrix $J=\frac{\partial f}{\partial x}$ evaluated at the equilibrium points satisfy

[TABLE]

For any quadratic polynomial $\phi(x)=x^{2}+a_{1}x+a_{2}$, the discriminant $D(\phi)$ of the polynomial $\phi$ is given by

[TABLE]

The generalized Routh-Hurwitz stability conditions for fractional-order systems are then given by the following proposition [2, 32, 33].

Proposition 2.1

• (i)

If $D(\phi)\geq 0$, $a_{1}>0$ and $a_{2}>0$, then the equilibrium $E^{*}$ is locally asymptotically stable for $0<m\leq 1$.

• (ii)

If $D(\phi)<0$, $a_{1}<0$ and $\mid\tan^{-1}(\frac{\sqrt{4a_{2}-a_{1}^{2}}}{a_{1}})\mid>\frac{m\pi}{2}$, $0<m<1$, then the equilibrium $E^{*}$ is locally asymptotically stable.

The system (2) has three equilibrium points: $(i)$ $E_{0}=(0,0)$ as the trivial equilibrium, $(ii)$ $E_{1}=(K,0)$ as the predator-free equilibrium and $(iii)$ $E^{*}=(x^{*},y^{*})$ as the interior equilibrium, where

[TABLE]

Note that the equilibria $E_{0}$ and $E_{1}$ always exist. The interior equilibrium $E^{*}$ exists if $0<c<c_{1}$ and $\theta>\theta_{1}$, where $c_{1}=1-\frac{d}{\alpha K(\theta-hd)}$, $\theta_{1}=hd+\frac{d}{\alpha K}$.

Theorem 2.3 *(a) The trivial equilibrium point $E_{0}$ is a saddle point. (b) The predator-free equilibrium point $E_{1}$ is locally asymptotically stable if $c>c_{1}$, $\theta>\theta_{1}$ and a saddle if $c<c_{1}.$ *

Proof The proof of part (a) is straightforward and omitted. The Jacobian matrix corresponding to $E_{1}$ is given by

[TABLE]

The corresponding eigenvalues are $\xi_{1}=-r~{}(<0)$, $\xi_{2}=\frac{\theta\alpha K(1-c)}{1+\alpha K(1-c)h}-d$. If $c<c_{1}$, then $\xi_{2}>0$ and $\mid arg(\xi_{2})\mid=0$. In this case, $E_{1}=(K,0)$ is a saddle point.

If $c>c_{1}$ and $\theta>\theta_{1}$ then $\xi_{2}<0$. Consequently, $\mid arg(\xi_{i})\mid=\pi>\frac{m\pi}{2},\forall m\in(0,1]$, $i=1,2,$ and the equilibrium $E_{1}=(K,0)$ is locally asymptotically stable. In other words, when the degree of complexity is high and the conversion efficiency of predator exceeds some lower threshold value, then the predator-free equilibrium becomes locally asymptotically stable.

To prove the global stability of $E_{1}$, we use the following Lemma.

Lemma 2.4 [14] Let $x(t)\in\Re_{+}$ be a continuous and derivable function. Then for any time instant $t>t_{0}$

[TABLE]

Theorem 2.4 *The predator-free equilibrium $E_{1}$ is globally asymptotically stable for any $m\in(0,1]$ if $c>c_{1}$, $\theta>\theta_{1}$, where $c_{1}=1-\frac{d}{\alpha K(\theta-hd)}$, $\theta_{1}=hd+\frac{d}{\alpha K}$.

Proof Consider the Lyapunov function

[TABLE]

Here $V(x,y)>0$ for all values of $x(t),y(t)>0$ and $V=0$ only at $E_{1}=(K,0)$. Calculating the $mth$ order fractional derivative of $V(x,y)$ along the solution of (2) and using Lemma $2.4$ when $t_{0}=0$, we have

[TABLE]

One can note that ${}^{c}_{0}D^{m}_{t}V(x,y)\leq 0,\forall(x,y)\in R^{2}_{+}$ if $\alpha K(1-c)-\frac{d}{\theta}<0$, i.e., if $\frac{d}{\theta}>\alpha K(1-c)>\frac{\alpha K(1-c)}{1+\alpha(1-c)hK}$. This implies ${}^{c}_{0}D^{m}_{t}V(x,y)\leq 0,\forall(x,y)\in R^{2}_{+}$ if $c>c_{1}$, $\theta>\theta_{1}$ and ${}^{c}_{0}D^{m}_{t}V(x,y)=0$ at $E_{1}$. Therefore, the only invariant set on which ${}^{c}_{0}D^{m}_{t}V(x,y)=0$ is the singleton $\{E_{1}\}$. Then by Lemma $4.6$ in [15], it follows that the predator-free equilibrium $E_{1}$ is globally asymptotically stable if $c>c_{1}$ and $\theta>\theta_{1}$. This completes the proof.

Remark 2.1 It is to be noted that stability of the predator-free equilibrium does not depend on the fractional-order $m$.

Theorem 2.5 The following statements are true for the stability of the interior equilibrium point $E^{*}$ of system (2).

• (a)

If  $trace(J^{*})<0$, i.e. if   $c_{2}<c<c_{1}$ with $\theta>\theta_{2}$, $\alpha>\frac{1}{Kh}$ then the interior equilibrium $E^{*}$ is locally asymptotically stable for $0<m\leq 1$, where $c_{2}=1-\frac{\theta+hd}{\alpha Kh(\theta-hd)}$, $c_{1}=1-\frac{d}{\alpha K(\theta-hd)}$ and $\theta_{2}=\frac{hd(\alpha Kh+1)}{\alpha Kh-1}.$

• (b)

If  $0<trace(J^{*})<2\sqrt{det(J^{*})}$, i.e. if   $0<c<c_{2}$ with $\theta>\theta_{2}$, $\alpha>\frac{1}{Kh}$ then for any $m\in(0,m^{*})$, the interior equilibrium $E^{*}$ is locally asymptotically stable and unstable for any $m\in(m^{*},1]$. A Hopf bifurcation occurs at $m=m^{*}$, where $m^{*}=\frac{2}{\pi}\mid\cos^{-1}(\frac{trace(J^{*})}{2\sqrt{det(J^{*})}})\mid$.

• (c)

If  $trace(J^{*})\geq 2\sqrt{det(J^{*})}$, then the interior equilibrium $E^{*}$ is unstable for any  $m\in(0,1]$.

Proof For the interior equilibrium $E^{*}$, the Jacobian matrix is given by

[TABLE]

The corresponding characteristic equation is given by

[TABLE]

where $trace(J^{*})=\frac{rd\{\alpha hk(1-c)-1-2\alpha h(1-c)x^{*}\}}{\alpha\theta K(1-c)}$ and $det(J^{*})=\frac{rd\{\alpha K(1-c)(\theta-hd)-d\}}{K\{\alpha(\theta-hd)+\alpha hd\}(1-c)}$.

Therefore, the roots of this equation are given by

[TABLE]

(a) Note that $trace(J^{*})$ will be negative if $\alpha hk(1-c)-1-2\alpha h(1-c)x^{*}<0$, i.e., if $c>1-\frac{\theta+hd}{\alpha Kh(\theta-hd)}=c_{2}$ with $\theta>\frac{hd(1+\alpha Kh)}{\alpha Kh-1}$, $\alpha>\frac{1}{Kh}$. Since $det(J^{*})>0$, both roots of (12) are negative real or complex conjugate with negative real parts. Hence $\mid arg(\xi_{1,2})\mid>\frac{m\pi}{2},\forall m\in(0,1]$. So the positive interior equilibrium $E^{*}$ is locally asymptotically stable for $0<m\leq 1$ if $c_{2}<c<c_{1}$ with $\theta>max[hd+\frac{d}{\alpha K},\frac{hd(1+\alpha Kh)}{\alpha Kh-1}]$, $\alpha>\frac{1}{Kh}$. Noting that $max[hd+\frac{d}{\alpha K},\frac{hd(1+\alpha Kh)}{\alpha Kh-1}]=\frac{hd(1+\alpha Kh)}{\alpha Kh-1}=\theta_{2}$, one gets the required stability result. This completes the proof of (a).

(b) The condition $0<trace(J^{*})<2\sqrt{det(J^{*})}$ will hold if $0<c<c_{2}$ with $\theta>\theta_{2}$, $\alpha>\frac{1}{Kh}$, where $c_{2}=1-\frac{\theta+hd}{\alpha Kh(\theta-hd)}$. Since $0<trace(J^{*})<2\sqrt{det(J^{*})}$, the equation (12) has two complex conjugate roots with positive real part given by

[TABLE]

with $\mid arg(\xi_{1,2})\mid=\mid\cos^{-1}(\frac{trace(J^{*})}{2\sqrt{det(J^{*})}})\mid$. Assume that there exists a $m^{*}\in(0,1]$ such that $\mid\cos^{-1}(\frac{trace(J^{*})}{2\sqrt{det(J^{*})}})\mid=\frac{m^{*}\pi}{2}$. Then, following Theorem $3.1$, we have $\mid arg(\xi_{i})\mid>\frac{m\pi}{2}$ for all $m\in(0,m^{*})$ and $\mid arg(\xi_{i})\mid<\frac{m\pi}{2}$ for all $m\in(m^{*},1]$. Therefore, the positive interior equilibrium $E^{*}$ is locally asymptotically stable for  $0<m<m^{*}$ and unstable for $m^{*}<m\leq 1$ when $0<trace(J^{*})<2\sqrt{det(J^{*})}$, i.e., when $0<c<c_{2}$. A Hopf bifurcation will occur at $m=m^{*}$ under the following conditions [38, 39]:

[TABLE]

Note that $Real(\xi_{i})=\frac{1}{2}trace(J^{*})>0$ and $\mid arg(\xi_{i})\mid=\frac{m^{*}\pi}{2},~{}i=1,2,$ by assumption. Also, $\frac{d}{dm}[Real{\xi_{i}}]\mid_{m=m^{*}}=\frac{\pi}{2}\neq 0.$ Therefore, a Hopf bifurcation exists as $m$ crosses the critical value $m^{*}$. The equilibrium $E^{*}$ is thus stable for all $m\in(0,m^{*})$ and unstable for all $m\in(m^{*},1].$ This completes the proof.

(c) As $trace(J^{*})\geq 2\sqrt{det(J^{*})}$, the equation (12) has two real roots given by $\xi_{1,2}=\frac{1}{2}[trace(J^{*})\pm\sqrt{trace(J^{*})^{2}-4det(J^{*})}]$. Now for the positive root $\xi_{1}$, we note that $\mid arg(\xi_{1})\mid=0$. Since the eigenvalue $\xi_{1}$ does not satisfy $\mid arg(\xi_{1})\mid>\frac{m\pi}{2},\forall m\in(0,1]$, therefore $E^{*}$ is unstable for any $m\in(0,1]$. This completes the proof of (c).

Theorem 2.6 The interior equilibrium $E^{*}$ is globally asymptotically stable for any $m\in(0,1]$ if $c_{2}<c<c_{1}$ with $\theta>\theta_{2},\alpha>\frac{1}{Kh}$, where $c_{2}=1-\frac{\theta+hd}{\alpha Kh(\theta-hd)},~{}c_{1}=1-\frac{d}{\alpha K(\theta-hd)}~{}and~{}\theta_{2}=\frac{hd(1+\alpha Kh)}{\alpha Kh-1}$.

Proof Let us consider the Lyapunov function

[TABLE]

Here $V(x,y)>0$ for all values of $x,y>0$ and $V=0$ only at $E^{*}=(x^{*},y^{*})$. Considering the $mth$ order fractional derivative of $V(x,y)$ along the solution of (2) and using Lemma $2.4$, we have

[TABLE]

One can note that ${}^{c}_{0}D^{m}_{t}V(x,y)\leq 0,\forall(x,y)\in R^{2}_{+}$ if $\alpha h(1-c)(K-x^{*})-1<0$, i.e., if $\alpha hK(1-c)<\frac{\theta}{\theta-hd}<\frac{\theta+hd}{\theta-hd}$. This implies ${}^{c}_{0}D^{m}_{t}V(x,y)\leq 0,\forall(x,y)\in R^{2}_{+}$ if $c_{2}<c<c_{1}$, $\theta>\theta_{2}$, $\alpha>\frac{1}{Kh}$ and ${}^{c}_{0}D^{m}_{t}V(x,y)=0$ implies that $(x,y)=(x^{*},y^{*})$. Therefore, the only invariant set on which ${}^{c}_{0}D^{m}_{t}V(x,y)=0$ is the singleton $\{E^{*}\}$. Then, following Lemma $4.6$ in [15], the interior equilibrium $E^{*}$ is globally asymptotically stable if the conditions in the theorem are satisfied. This completes the proof.

3 Discretized fractional-order model and its analysis

We first construct the discrete fractional-order model corresponding to the system (2). Following Elsadany and Matouk [18], discretization of the model system (2) with piecewise constant arguments can be done in the following manner:

[TABLE]

with initial condition $x(0)=x_{0}>0$ and $y(0)=y_{0}>0$.

Let $t\in[0,s)$, so that $t/s\in[0,1)$. In this case, we have

[TABLE]

and the solution of this fractional differential equation can be written as

[TABLE]

In the second step, we assume $t\in[s,2s)$ so that $t/s\in[1,2)$ and obtain

[TABLE]

The solution of this equation reads

[TABLE]

where $J^{m}_{s}=\frac{1}{\Gamma(m)}\int^{t}_{s}(t-\tau)^{(m-1)}d\tau,m>0$.

Repeating the discretization process $n$ times, we have

[TABLE]

where $t\in[ns,(n+1)s)$.

Making $t\rightarrow(n+1)s$, we obtain the corresponding fractional discrete model of the continuous fractional model (2) as

[TABLE]

It is noticeable that Euler discrete model is a special case of this generalized discrete model when $m\rightarrow 1$.

3.1 Existence and stability of fixed of points

In the following, we investigate the dynamics of the discretized fractional-order model (3). At the fixed point, we have $x_{n+1}=x_{n}=x$ and $y_{n+1}=y_{n}=y$. One can easily compute that (3) has the same fixed points as in the fractional-order system (2) given by $E_{0}=(0,0)$, $E_{1}=(K,0)$ and $E^{*}=(x^{*},y^{*})$, where

[TABLE]

The fixed point $E^{*}$ exists if $0<c<c_{1}$ and $\theta>\theta_{1}$, where $c_{1}=1-\frac{d}{\alpha K(\theta-hd)}$, $\theta_{1}=hd+\frac{d}{\alpha K}$.

The Jacobian matrix of system (3) at any arbitrary fixed point point $(x,y)$ reads

[TABLE]

where

[TABLE]

Let $\xi_{1}$ and $\xi_{2}$ be the eigenvalues of the Jacobian matrix (15). Then we have the following definition and lemma.

Definition 3.1 [41, 42] *A fixed point $(x,y)$ of system (3) is called stable if $\mid\xi_{1}\mid<1$, $\mid\xi_{2}\mid<1$ and a source if $\mid\xi_{1}\mid>1$, $\mid\xi_{2}\mid>1$. It is called a saddle if $\mid\xi_{1}\mid<1$, $\mid\xi_{2}\mid>1$ or $\mid\xi_{1}\mid>1$, $\mid\xi_{2}\mid<1$ and a nonhyperbolic fixed point if either $\mid\xi_{1}\mid=1$ or $\mid\xi_{2}\mid=1$. It is called a spiral source if $\xi_{1,2}=\alpha\pm i\beta,\beta\neq 0,\alpha,\beta\in R$ and $\mid\xi_{1,2}\mid>1$.

Lemma 3.1 [41] Let $\xi_{1}$ and $\xi_{2}$ be the eigenvalues of Jacobian matrix (15). Then $\mid\xi_{1}\mid<1$ and $\mid\xi_{2}\mid<1$ if the following condition holds:

[TABLE]

Theorem 3.1 * (a) The fixed point $E_{0}$ is always unstable for $0<m\leq 1$. It will be a saddle point if $0<s<s_{1}$ and a source if $s>s_{1}$. If $s=s_{1}$, then $E_{0}$ is nonhyperbolic, where $s_{1}=\sqrt[m]{\frac{2m\Gamma(m)}{d}}$. *

(b) The fixed point $E_{1}$ is stable for $0<m\leq 1$ if $c>c_{1}$, $s<min\{s_{2},s_{3}\}$, where $s_{2}=\sqrt[m]{\frac{2m\Gamma(m)}{r}}$, $s_{3}=\sqrt[m]{\frac{2m\Gamma(m)\{1+\alpha Kh(1-c)\}}{d-K\alpha(1-c)(\theta-hd)}}$. It is a saddle point if $c>c_{1}$, $s_{3}<s<s_{2}$ ; or $c>c_{1}$, $s_{2}<s<s_{3}$; and a source if $c>c_{1}$, $s>max\{s_{2},s_{3}\}$.

*(c) The fixed point $E^{*}$ is locally asymptotically stable for $0<m\leq 1$ if $c_{2}<c<c_{1}$ with $\theta>\theta_{2},\alpha>\frac{1}{Kh}$ and $s<min\{s_{4},s_{5}\}$, where $s_{4}=\sqrt[m]{\frac{m\Gamma(m)G}{H}}$, $s_{5}=\sqrt[m]{\frac{2m\Gamma(m)}{G}}$,

$G=\frac{rx^{*}}{K\theta}[\theta+hd-\alpha hK(1-c)(\theta-hd)]$ and $H=\frac{rx^{*}(\theta-hd)}{K\theta}[\alpha K(1-c)(\theta-hd)-d]$*.

Proof At the fixed point $E_{0}$, the eigenvalues are $\xi_{1}=1+r\frac{s^{m}}{m\Gamma(m)}$ and $\xi_{2}=1-d\frac{s^{m}}{m\Gamma(m)}$. Since $|\xi_{1}|>1$, $E_{0}$ is always unstable for $0<m\leq 1$. In fact, it is a saddle point if $0<s<\sqrt[m]{\frac{2m\Gamma(m)}{d}}$ for which $|\xi_{2}|<1$ and a source if $s>\sqrt[m]{\frac{2m\Gamma(m)}{d}}$ for which $|\xi_{2}|>1$. Again, it becomes nonhyperbolic if $s=\sqrt[m]{\frac{2m\Gamma(m)}{d}}$ for any $m\in(0,1]$.

The eigenvalues evaluated at the fixed point $E_{1}$ are evaluated as

[TABLE]

Note that for $0<m\leq 1$, $|\xi_{1,2}|<1$ hold if

[TABLE]

Therefore, $E_{1}$ is locally asymptotically stable for $0<m\leq 1$ if $c>c_{1}$ and $s<min\{s_{2},s_{3}\}$. However, $|\xi_{1}|>1$ if $s>s_{2}$ and $|\xi_{2}|>1$ if $s>s_{3}$ with $c>c_{1}$. Thus, $E_{1}$ will be a source if $c>c_{1}$ and $s>max\{s_{2},s_{3}\}$. The fixed point $E_{1}$ will be a saddle point if either of the conditions $(i)~{}s_{3}<s<s_{2}$, $c>c_{1}$ or $(ii)~{}s_{2}<s<s_{3}$, $c>c_{1}$ holds.

At the interior fixed point $E^{*}$, the Jacobian matrix is evaluated as

[TABLE]

where $a_{11}=1-\frac{s^{m}}{m\Gamma(m)}G,$ $a_{12}=-\frac{s^{m}}{m\Gamma(m)}\frac{\alpha(1-c)(\theta-hd)x^{*}}{\theta},$ $a_{21}=\frac{s^{m}}{m\Gamma(m)}\frac{r(\theta-hd)(K-x^{*})}{K},$ $a_{22}=1$ and $a_{12}a_{21}=-{(\frac{s^{m}}{m\Gamma(m)})}^{2}H$ with $G=\frac{rx^{*}}{K\theta}[\theta+hd-\alpha hK(1-c)(\theta-hd)]$ and

$H=\frac{rx^{*}(\theta-hd)}{K\theta}[\alpha K(1-c)(\theta-hd)-d]$.

Note that $H>0$ if $c<c_{1}$ and $G>0$ if $c>c_{2},\theta>\theta_{2},\alpha>\frac{1}{Kh}$. After some algebraic manipulations, we have

[TABLE]

Thus, $1-trace(J)+det(J)={(\frac{s^{m}}{m\Gamma(m)})}^{2}H>0$ if $c<c_{1}$. Also, $1-det(J)=(\frac{s^{m}}{m\Gamma(m)})(G-(\frac{s^{m}}{m\Gamma(m)})H)$ is positive if $s<s_{4}$, where $s_{4}=\sqrt[m]{\frac{m\Gamma(m)G}{H}}$ and $c>c_{2}$ with $\theta>\theta_{2},\alpha>\frac{1}{Kh}$.

One can compute that

[TABLE]

This expression will be positive if $0<s<s_{5}$, where $s_{5}=\sqrt[m]{\frac{2m\Gamma(m)}{G}}$. Therefore, the fixed point $E^{*}$ is stable if $c_{2}<c<c_{1}$ and $s<min\{s_{4},s_{5}\}$ for any $m\in(0,1]$ and unstable otherwise. Hence the theorem.

Remark 3.1 Here we also observe that the predator-free fixed point $E_{1}$ looses its stability through transcritical bifurcation (a real eigenvalue that passes through $+1$) when $1-trace(J)+det(J)=0$ at $c=c_{1}$ for any $m$ [18]. Again our model system (3) undergoes a flip bifurcation (a real eigenvalue becomes equal to $-1$) when $1+trace(J)+det(J)=0$ at the predator-free fixed point $E_{1}$ for $c=c_{1}$ and $s=s_{5}=\sqrt[m]{\frac{2m\Gamma(m)}{G}}$.

Remark 3.2 Note that the eigenvalues of $J(x^{*},y^{*})$ are

[TABLE]

Therefore, $\xi_{1,2}$ are complex conjugate if $G^{2}-4H<0$, i.e., if $-2\sqrt{H}<G<2\sqrt{H}$. Now,

[TABLE]

and this modulus is equal to unity if $det(J)=1$, i.e., if $s=\sqrt[m]{\frac{m\Gamma(m)G}{H}}=s_{4}$. Since $G=\frac{s^{m}}{m\Gamma(m)}H>0,$ the previous inequality becomes $0<G<2\sqrt{H}$. Therefore, we can conclude that $J(x^{*},y^{*})$ has complex conjugate roots with unit modulus if parameters belong to the set

[TABLE]

Therefore, if the parameter $s$ varies in the neighborhood of $s_{4}$ and $(m,s,r,K,\alpha,\theta,h,c,d)\in U$, the system (3) may undergo a Hopf bifurcation around the equilibrium $E^{*}$.

3.2 Hopf Bifurcation and its stability

Here we prove the existence of Hopf bifurcation around $E^{*}=(x^{*},y^{*})$ and its stability. Let $S_{1}=\frac{s^{m}}{m\Gamma(m)}$ and $S^{*}$ be a perturbation in the bifurcation parameter $S_{1}$, where $|S^{*}|<<1$. Then a perturbation form of model (3) can be represented as [35]

[TABLE]

Let $X_{n}=x_{n}-x^{*},Y_{n}=y_{n}-y^{*}$ so that the fixed point $E^{*}=(x^{*},y^{*})$ of the map (3.2) is transformed into the origin. The transformed system reads

[TABLE]

where $c_{11}=1-(S_{1}+S^{*})G,$ $c_{12}=-(S_{1}+S^{*})\frac{\alpha(1-c)(\theta-hd)x^{*}}{\theta},$ $c_{21}=(S_{1}+S^{*})\frac{r(\theta-hd)(K-x^{*})}{K},$ $c_{22}=1,$ $c_{13}=-\frac{\alpha(1-c)(S_{1}+S^{*})}{2(1+\alpha(1-c)hx^{*})^{2}},$ $c_{23}=\frac{\theta\alpha(1-c)(S_{1}+S^{*})}{2(1+\alpha(1-c)hx^{*})^{2}},$ $c_{12}c_{21}=-(S_{1}+S^{*})^{2}H$ and $S_{1}=\frac{s_{4}^{m}}{m\Gamma(m)}$ with $G=\frac{rx^{*}}{K\theta}[\theta+hd-\alpha hK(1-c)(\theta-hd)]$, $H=\frac{rx^{*}(\theta-hd)}{K\theta}[\alpha K(1-c)(\theta-hd)-d]$.

The characteristic equation associated with the linearization of the model (3.2) at $(X_{n},Y_{n})=(0,0)$ is given by

[TABLE]

where

[TABLE]

Since the parameters $(m,s,r,K,\alpha,\theta,h,c,d)\in U$ and $S^{*}$ varies in a small neighborhood of $S^{*}=0$, and the roots of (18) are pair of complex conjugate numbers $\lambda_{1}$ and $\lambda_{2}$ denoted by

[TABLE]

Therefore, $|\lambda_{1,2}|=\sqrt{q(S^{*})}$. Since $q(S^{*})=1$ at $S^{*}=0$, when $s=s_{4}=\sqrt[m]{\frac{m\Gamma(m)G}{H}}$, then $|\lambda_{1,2}|=1$ at $S^{*}=0$ for $s=s_{4}$.

Consequently for $s=s_{4}$,

[TABLE]

Also, at $S^{*}=0$, $\lambda_{1,2}^{n}\neq 1$ for $n=1,2,3,4$ (nonresonance conditions), which is equivalent to

[TABLE]

Since $p^{2}(0)-4q(0)<0$ and $q(0)=1$, we have $p^{2}(0)<4$; then $p(0)\neq\pm 2$. It is only require that $p(0)\neq 0,1$, which leads to

[TABLE]

for $s=s_{4}.$ Next, we study the normal form of the model (3.2) at $S^{*}=0$. Let $\delta=Re(\lambda_{1,2})$ and $\beta=Im(\lambda_{1,2})$. We construct an invertible matrix

[TABLE]

and consider the translation

[TABLE]

Thus, the map (3.2) becomes

[TABLE]

where

[TABLE]

In order to undergo Hopf bifurcation, we require that the following discriminatory quantity $\gamma$ be nonzero

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

From the above analysis and Theorem $3.2$ in [34], following theorem can be stated.

Theorem 3.2 *If conditions (21) and (24) hold, then the system (3) undergoes Hopf bifurcation at the positive fixed point $E^{*}=(x^{*},y^{*})$ when the parameter $s$ varies in the small neighborhood of $s_{4}$. Furthermore if $\gamma<0$ (respectively $\gamma>0$), then an attracting (respectively repelling) invariant closed curve bifurcates from the fixed point $E^{*}=(x^{*},y^{*})$ for $s>s_{4}$ (respectively $s<s_{4}$), where $s_{4}=\sqrt[m]{\frac{m\Gamma(m)G}{H}}$.

A comparison table on dynamical behaviors of system (1) with corresponding fractional-order and discretized fractional-order versions has been given in Table $3.1$.

Table 3.1. Comparison of dynamical behaviors of three systems.

[TABLE]

4 Numerical Simulations

In this section, we perform extensive numerical computations of fractional-order differential equations (FDE) system (2) for different fractional values of $m$ $(0<m\leq 1)$ as well as the fractional-order discrete system (3). We use Adams-type predictor corrector method for the numerical solution of FDE system (2). It is an effective method to give numerical solutions of both linear and nonlinear FDE [36, 37]. We first replace our system (2) by the following equivalent fractional integral equations:

[TABLE]

and then apply the PECE (Predict, Evaluate, Correct, Evaluate) method.

Three examples are presented to illustrate the analytical results of FDE system obtained in the previous section. To explore the effect of habitat complexity and fractional-order, we varied $c$ and $m$ in their respective ranges $0<c<1$ and $0<m<1$. We also plotted the solutions for $m=1$, whenever necessary, to compare the solutions of fractional-order system with that of integer order system.

Example 1: We considered the parameter values as $r=2.65$, $K=898$, $\alpha=0.045$, $h=0.0437$, $d=1.06$ and initial point $x(0)=10,y(0)=5$ from [26]. Step size in all simulations is considered as $0.05$. Note that the condition $\alpha>\frac{1}{Kh}$ is always satisfied by this parameter set. Following Theorem $2.3(b)$, we compute $c_{1}=0.8445$, $\theta_{1}=0.0726$ and select $c=0.86(>c_{1})$, $\theta=0.215(>\theta_{1})$ to show that the predator-free equilibrium $E_{1}$ of the system (2) is asymptotically stable for all $m\in(0,1]$ (Fig. 1). It is noticeable that the solutions reach to the equilibrium more slowly as the value of $m$ gets smaller. The phase planes presented in Fig. 2 show that the solution trajectory with different initial conditions (denoted by stars) reach to the equilibrium point (red circle) in each case, following Theorem $2.4$, depicting the global stability of the predator-free equilibrium $E_{1}$ for different values of $m$.

Example 2: For the same parameter values as in Example 1, we compute $\theta_{2}=0.1673$, $c_{2}=0.1227$, $c_{1}=0.8445$ and $trace(J(E^{*}))=-0.3398<0$. Thus, following Theorem $2.5(a)$, if we choose $\theta=0.215~{}(>\theta_{2})$ and $c=0.45~{}(c_{2}<c<c_{1})$ then solutions for all $m$ eventually converge to the equilibrium point $E^{*}$ where both the prey and predator populations coexist in the form of a stable equilibrium (Fig. 3).

If we choose $c<c_{2}=0.1227,$ say $c=0.05,$ and $\theta=0.215(<\theta_{2})$ as before then we obtain $trace(J^{*})=0.0437$, $2\sqrt{det(J^{*})}=2.7152$. Therefore, from Theorem $2.5(b)$, there exists a critical value $m=m^{*}=0.9898$ below which $E^{*}$ is stable and above which it is unstable. The stable behavior of the system (2) for $m=0.95(<0.9898)$ is presented in Fig. 4a and the unstable behavior of the system for $m=0.995(>0.9898)$ in Fig. 4b. A Hopf bifurcation occurs at $m=m^{*}$. One can obtain a critical value $m^{*}$ for each $c\in(0,c_{2})$, following Theorem 2.5(b), and can draw a stability region of $E^{*}$ in $c-m$ plane. The bifurcation curve separates the stable and unstable region (see Fig. 4c).

Example 3: Global stable behavior of system (2) around the interior equilibrium $E^{*}$ is presented in Fig. 5. This figure shows that solutions with different initial conditions converge to the coexisting equilibrium $E^{*}$ for all values of $m$ when the conditions of Theorem $2.6$ are fulfilled.

Example 4: To illustrate the corresponding discrete system (3) of the fractional-order system (2), we consider the same parameter set as in Example 1. Stability of fixed points depends on the step size $s$ for different fractional-order $m$ (see Theorem 3.1). Assigning $c=0.86$ for $E_{1}$ and $c=0.45$ for $E^{*}$, we present the ranges of step size for the stability of the corresponding fixed point for different values of fractional-order $m$ in Table 5.1.

Table 5.1. Restriction on the step size, following Theorem 3.1, for the stability of fixed points $E_{1}$ and $E^{*}$ for different fractional-order $m$.

[TABLE]

For example, when $m=0.95$ then the step size $s$ should be less than $0.7279$ and $0.1940$ for $E_{1}$ and $E^{*}$, respectively, to be stable whenever they exist and unstable if it exceeds. In Fig. 6a, we have plotted the stable behavior of the fixed point $E_{1}$ for $s=0.68$ and the unstable oscillatory behavior is presented in Fig. 6b for $s=0.8$. Similar stable and unstable behaviors of the fixed point $E^{*}$ are presented in Figs. 6c-6d for $s=0.12$ and $s=0.22$, respectively. Following Theorem 3.2, we obtain a pair of complex conjugate eigenvalues as $\lambda_{1,2}=0.9635\pm 0.2678i$, where \mid\lambda_{1,2}\mid=1,~{}\frac{d|\lambda_{1,2}|}{dS^{*}}\bigg{|}_{S^{*}=0}=0.1699>0 and $\gamma=-0.000000019961<0$ at $(m,s,r,K,\alpha,\theta,h,c,d)=(0.95,0.194,2.65,898,0.045,0.215,0.0437,0.45,1.06)\in U$. This implies that the system (3) undergoes a Hopf bifurcation at the fixed point $E^{*}=(253.9056,97.8867)$ for $s=0.194~{}(=s_{4})$. The bifurcation diagrams (Figs. 7a, 7c) represents it succiently. Also, the bifurcating closed curve is stable as $\gamma$ is negative. The system shows period doubling bifurcations leading to chaos as the step-size is further increased. Clear maginifed pictures of period doubling are presented in Figs. 7b & 7d. Phase diagrams of the fractional-order discrete system (3) for some particular values of $s$ are presented in Fig. 8.

One can also compute $1+trace(J)+det(J)=0$ and $\xi_{1}=1-r\frac{s^{m}}{m\Gamma(m)}=-1,$ \xi_{2}=1+\frac{s^{m}}{m\Gamma(m)}\bigg{(}\frac{\theta\alpha(1-c)K}{1+\alpha(1-c)hK}-d\bigg{)}=1 at the predator-free fixed point $E_{1}$ for $c=c_{1}$ and $s=s_{5}=\sqrt[m]{\frac{2m\Gamma(m)}{G}}=0.7279$. Thus, following Remark 3.1, the predator-free fixed point $E_{1}$ undergoes a flip bifurcation (Fig. 9) at the critical step size $s=s_{5}$.

5 Discussion

This paper generalizes the results of continuous system predator-prey model that considers the effect of habitat complexity. This generalization has been accomplished in two phases. In the first phase, we constructed a fractional-order predator-prey model considering the fractional derivatives in Caputo sense. In the second phase, the fractional-order predator-prey model was discretized. Rigorous mathematical and computational results in relation to the stability of both the systems was presented. Proving existence of Hopf bifurcation with respect to the fractional-order of the derivatives in both discrete and fractional systems is rare in contemporary studies. We have presented it both theoretically and numerically, showing the novelty of this study. For the fractional-order system, we proved different mathematical results like positivity and boundedness, local and global stability of different equilibrium points. It is shown that the trivial equilibrium $E_{0}$ is always unstable saddle and the predator-free equilibrium is globally asymptotically stable for any value of fractional-order $m\in(0,1]$ if the degree of habitat complexity exceeds some upper threshold value $c_{1}$. The solution, however, takes more time to reach the predator-free equilibrium as the value of fraction order is reduced. At the intermediate level of habitat complexity $(c_{2}<c<c_{1})$, the system becomes both locally and globally asymptotically stable around $E^{*}$ for any value of $m$. These dynamics are consistent with the integer order system $m=1$. Stability of the interior equilibrium, however, depends on the fractional-order $m$ if the strength of habitat complexity is very low and the system shows order-dependent instability. If $0<c<c_{2}$ then there exists a critical value $m^{*}$ of the fractional-order $m\in(0,1)$ such that the coexistence equilibrium is stable if $m<m^{*}$ and unstable if $m$ crosses $m^{*}$. In case of integer order system ($m=1$), the coexistence equilibrium is, however, unstable for all $c\in(0,c_{2})$. Simulation results also agree perfectly with the analytical results. Discretization of the fractional-order system was done with piecewise constant arguments and the dynamics of this discrete model was explored. It is observed that the dynamics of the discrete system depends on both the step-size and fractional-order. Existence of Hopf and flip bifurcations have been shown both theoretically and numerically. It is also observed that the discrete fractional-order system shows more complex dynamics as the step size becomes larger. Our simulation results revealed that the discrete system shows period doubling route to chaos for larger step-size.