ANALYSIS OF A DISCRETIZED FRACTIONAL-ORDER PREY–PREDATOR MODEL UNDER WIND EFFECT

Abstract

While constructing mathematical models, scientists usually consider biotic factors, but it is crystal-clear that abiotic factors, such as wind, are also important as biotic factors. From this point of view, this paper is devoted to the investigation of some bifurcation properties of a fractional-order prey–predator model under the effect of wind. Using fractional calculus is very popular in modelling, since it is more effective than classical calculus in predicting the system’s future state and also discretization is one of the most powerful tools to study the behaviour of the models. In this paper, first of all, the model is discretized by using a piecewise discretization approach. Then, the local stability of fixed points is considered. We show using the centre manifold theorem and bifurcation theory that the system experiences a flip bifurcation and a Neimark–Sacker bifurcation at a positive fixed point. Finally, numerical simulations are given to demonstrate our results.

1. Introduction

Differential equations with fractional orders are widely used in science, and this extensive use has recently sparked much research interest. The so-called fractional-order integral or derivative operators are the main producers of fractional-order differential equations. Due to the nonlocal nature of their integral and derivative operators, fractional-order operators are more effective than other classical deterministic operators at predicting the future state of the system, because they depend not only on the current state, but also on all of its past states. The most popular and widely used fractional-order operators in applications of fractional-order differential equations are the Riemann–Liouville, Caputo and Grünwald–Letnikov kinds [28]. Now, let us give the definition of fractional-order integration and differentiation [44, 46].

Definition 1.1. The fractional integral (or the Riemann Liouville integral) of order $\beta \in \mathbb {R}^+$ of the function $f(t), \hspace {0.5em} t>0$ , is defined by

$$ \begin{align*} I_{a}^{\beta}f(t)=\int_a^t \frac{(t-s)^{\beta-1}}{\Gamma(\beta)} f(s)\,ds. \end{align*} $$

The fractional derivative of order $\alpha \in (n-1,n)$ is defined by the following two (nonequivalent) ways.

1. (1) Riemann–Liouville fractional derivative: take the fractional integral of order $(n-\alpha )$ , and then take the nth derivative as follows:

   $$ \begin{align*} D^\alpha f(t)= D^{n} I_{a}^{(n-\alpha)} f(t), \quad D =\frac{d}{dt}, \quad n=1,2,3,\ldots. \end{align*} $$

2. (2) Caputo-fractional derivative: take the nth derivative, and then take a fractional integral of order $(n-\alpha )$

   $$ \begin{align*} D^\alpha f(t)=I_{a}^{(n-\alpha)} D^{n} f(t). \end{align*} $$

Even though it is more restrictive than the Riemann–Liouville, we consider the fractional derivative provided by Caputo in this study [11]. This is because this definition is more suitable for problems consisting of fractional differential equations with initial conditions. The following result contains the main properties from fractional calculus [45].

Remark 1.2. Let $L^1=L^1[a,b]$ be the class of Lebesgue integrable functions on $[a,b], a<b<\infty $ , $\beta , \gamma \in \mathbb {R^+}$ and $\alpha \in (0,1)$ . Then:

1. (1) if $I_{a}^{\beta }:L^1\rightarrow L^1$ and $f(x) \in L^1$ , then $I_{a}^{\gamma }I_{a}^{\beta }f(x)=I_{a}^{\gamma +\beta } f(x)$ ;

2. (2) $\lim _{\beta \to n }I_{a}^{\beta }f(x){\kern-1pt}={\kern-2pt} I_{a}^{n}f(x)$ uniformly on $[a,b], n{\kern-1pt}={\kern-1pt}1,2,3,\ldots ,$ where $I_{a}^{1}f(t){\kern-1pt}={\kern-2.5pt}\int _0^t f(s) ds$ ;

3. (3) $\lim _{\beta \to 0 }I_{a}^{\beta }f(t)=f(t)$ weakly;

4. (4) if $f(t)$ is absolutely continuous on $[a,b]$ , then $\lim _{\alpha \to 1} D^\alpha f(t)={df(t)}/{dt}$ .

In applied mathematics, the process of converting continuous models defined by differential equations into their discrete equivalent is known as discretization. Analysing the dynamics of the discretization is crucial for understanding how well the discrete model approximates the original continuous system. Discretization is often applied to systems to perform numerical simulations or to design algorithms that run on digital platforms. As the step size approaches zero, the discrete system converges to the original continuous system. In this limit, the discrete-time model approaches the behaviour of the differential equations, ensuring that any discrepancies between the two systems vanish. In the literature, there are a large number of nonlinear fractional differential equations without an analytical solution, so they need to be solved with the help of some numerical and discretization techniques such as the Adomian decomposition method [1, 22], differential transform method [40, 51], Euler method [34], extrapolation method [18], Grünwald–Letnikov method [30] and variational iteration method [55, 56].

In the last decade, we see a new type of discretization technique for fractional equations, called piecewise discretization, which is constructed with the help of the piecewise constant arguments [2, 19, 20, 23, 32, 33, 39, 42, 57, 58, 60, 61]. The systematic study of problems involving piecewise constant arguments began in the early 1980s with the work of Shah and Wiener [50], who introduced the term “differential equations with piecewise constant argument” (DEPCA). A comprehensive source on this class of equations is provided in [52]. Busenberg and Cooke [10] were pioneers in applying such deviating arguments to mathematical models, specifically in the study of vertically transmitted diseases, by reducing their analysis to discrete equations. The main source for DEPCA theory are the papers [17, 52]. In 1991, Györi was the first who used the piecewise constant argument for approximation [24]. After this work, there has been a significant interest in the usage of the piecewise constant argument for a numerical approach [16, 25–27]. In the late 2010s, Akhmet [4] introduced the equation

$$ \begin{align*} x'(t)=f(t, x(t), x\gamma(t))), \end{align*} $$

where $\gamma (t)$ is a piecewise constant argument of generalized type, that is, given $(t_{k})_{k \in \mathbb {Z}}$ and $(\zeta _{k})_{k \in \mathbb {Z}}$ such that $t_{k}<t_{k+1}$ for all $ k \in \mathbb {Z}$ with $\lim _{ k\rightarrow \pm \infty }t_{k}=\pm \infty $ and $t_{k} \leq \zeta _{k}\leq t_{k+1}$ , if $t \in I_{k}=[t_{k}, t_{k+1})$ , then $\gamma (t)=\zeta _{k}$ . These type of equations are called “differential equations with piecewise constant argument of generalized” (DEPCAG) type. They have continuous solutions, even when $\gamma (t)$ is not, producing a discrete equation. Several aspects of this equation have garnered much interest, see [3, 5, 13–15] and the references therein.

In this paper, we propose to use a piecewise discretization method. For this aim, let us first explain the usage of the method for any fractional-order differential equation. Agarwal et al. [2] considered the equation

(1.1) $$ \begin{align} D^{\alpha} x ( t)=f \bigg(x \bigg(r \bigg[\frac{t}{r} \bigg] \bigg)\bigg),\quad x (t )=x_0,\quad t\leq 0, \end{align} $$

as a discretized version of the equation

$$ \begin{align*} \begin{cases} D^{\alpha} x ( t) = f (x (t) ), & t>0, \\ x (0)=x_{0}, & t \leq 0, \end{cases} \end{align*} $$

where r is the discretization parameter. Let $t \in [0,r)$ , then ${t}/{r} \in [0,1) $ . So from (1.1),

$$ \begin{align*} D^{\alpha} x ( t) = f (x_{0}), \: t \in [0,r). \end{align*} $$

Thus,

$$ \begin{align*} x_{1} ( t) = x_{0}+ \frac{t^{\alpha}}{\Gamma (1+\alpha)}f (x_{0}). \end{align*} $$

Let $t \in [r,2r)$ , so ${t}/{r} \in [1,2) $ . Thus,

$$ \begin{align*} D^{\alpha} x ( t)=f (x_{1}(r)), \quad t \in [r,2r), \end{align*} $$

whose solution is

$$ \begin{align*} x_{2} ( t)=x_{1}(r)+\frac{(t-r)^{\alpha}}{\Gamma(1+\alpha)}f(x_{1}(r)). \end{align*} $$

Let $t \in [2r,3r)$ , then ${t}/{r} \in [2,3) $ . Thus, we obtain

$$ \begin{align*} D^{\alpha} x ( t)=f (x_{2}(r)), \quad t \in [2r,3r). \end{align*} $$

So,

$$ \begin{align*} x_{3} ( t)=x_{2}(2r)+\frac{(t-2r)^{\alpha}}{\Gamma(1+\alpha)}f(x_{2}(2r)). \end{align*} $$

Repeating the process, when $t \in [nr,(n+1)r)$ ,

$$ \begin{align*} D^{\alpha} x ( t)=f (x_{n}(nr)), \quad t \in [nr,(n+1)r), \end{align*} $$

which yields

$$ \begin{align*} x_{n+1} ( t)=x_{n}(nr)+\frac{(t-nr)^{\alpha}}{\Gamma(1+\alpha)}f(x_{n}(nr)). \end{align*} $$

Letting $t \rightarrow (n+1)r$ in the above, the corresponding difference equation is obtained as

(1.2) $$ \begin{align} y_{n+1}=y_n+\frac{r^\alpha}{\Gamma (\alpha+1)} f(y_n) \end{align} $$

with $x(nr)=y_n$ .

Prey and predator live in the same ecosystem, and both of them have their own, but mostly the same, biotic factors, such as speed, stealth and camouflage (to hide while pursuing the prey or to hide from the predator). They also possess good senses of smell, sight and hearing (to locate the prey or to detect the predator), and poison (to kill the prey or to spray when approached or bitten), among other traits. In mathematical biology, these factors are usually considered as the main subjects that affect the model. However, recently, there has been a significant increase in the number of studies exploring the impact of abiotic effects. These include the influence of predator-induced fear [6, 8, 36, 37], climate change [48, 49, 54], seasonal variations [12, 41, 47] and wind [7, 43]. These elements can of course be increased.

Now, let us turn the wheels to our problem. Barman et al. [7] investigated the effect of wind in a prey–predator model. They constructed two types of functional responses to describe the whole dynamics of the considered system under the fact that wind may either decrease or increase the predation rate of predators. They showed positivity and boundedness of the solutions of the systems and examined the stability of the equilibrium points. Also for particular cases, it is demonstrated that wind can change the stability of the equilibrium point through Hopf bifurcation. Furthermore, considering the fact that wind flow cannot be constant for a time period, they developed the discussed system and investigated this system with the help of numerical simulations. In their work, they considered

(1.3) $$ \begin{align} \begin{cases} \displaystyle\frac{dx}{dt}=rx \bigg( 1-\frac{x}{k}\bigg) - \displaystyle\frac{c_{1}xy}{1+\omega+bx+{b \omega x}/({1+\omega})}, \\[6pt] \displaystyle\frac{dy}{dt}= \displaystyle\frac{c_{2} x y}{1+ \omega+ b x+{b \omega x}/({1+\omega})}-d y , \end{cases} \end{align} $$

with the initial conditions

(1.4) $$ \begin{align} x(0)=x_0 \quad \textrm{and} \quad y(0)=y_0. \end{align} $$

The variables and all positive parameters in this model are defined in Table 1.

Table 1 Description of the variables and parameters in (1.3).

Our aim is to investigate the fractional-order discrete version of this problem to see the effect of both the fractional order and the discretization parameter in the occurrence of flip and Neimark–Sacker bifurcation.

In recent years, there has been a considerable amount of work on flip and Neimark–Sacker bifurcations: in [31], a discrete-time predator–prey system is considered, and the chaos control is investigated. The effect of prey refuge proportional to the predator in a discrete-time prey–predator model is studied in [38]. The authors consider the dynamical behaviour of a discrete-time prey–predator system with Leslie type in [9], and in [59], a discrete Lotka–Volterra-type predator–prey system with refuge effect is analysed. The reader is also referred to the references therein.

Before stating our main problem, let us first apply nondimensionalization to (1.3) to use one of the most important benefits: reduction of the numbers of the parameters. For this aim, let us define the variables

(1.5) $$ \begin{align} x=k P,\quad y=\frac{r Q}{c_1},\quad t=\frac{\tau}{r}. \end{align} $$

Then, the derivatives become

(1.6) $$ \begin{align} \frac{dx}{dt}=k r\frac{dP}{d\tau}, \quad \frac{dy}{dt}=\frac{r^2}{c_1}\frac{dQ}{d\tau}. \end{align} $$

Substituting (1.5) and (1.6) in the initial value problem (1.3)–(1.4) gives the nondimensionalized form

$$ \begin{align*} \begin{cases} \displaystyle\frac{dP}{d\tau}=P (1-P ) - \displaystyle\frac{P Q}{a_1+a_2 P},\\[6pt] \displaystyle\frac{dQ}{d\tau}= \displaystyle\frac{b_{1} P Q}{a_1+a_2 P} -b_{2} Q, \end{cases} \end{align*} $$

with

$$ \begin{align*} P(0)=P_0,\quad Q(0)=Q_0, \end{align*} $$

where

$$ \begin{align*} a_{1}=1+\omega,\quad a_{2}=bk+\frac{b \omega k}{1+\omega},\quad b_{1}=\frac{c_2 k}{r},\quad b_{2}=\frac{d}{r}, \quad P_{0}=\frac{x_0}{k},\quad Q_{0}=\frac{y_0 c_1}{r}. \end{align*} $$

Now, let us consider this initial value problem with fractional derivatives and apply the piecewise discretization process. For $\alpha \in (0,1)$ and $h>0$ as the discretization parameter, it can be written as

(1.7) $$ \begin{align} \begin{cases} D^\alpha P( \tau )= P \bigg( \bigg[ \displaystyle\frac{\tau}{h}\bigg] h \bigg) \bigg(1-P\bigg( \bigg[ \displaystyle\frac{\tau}{h}\bigg] h \bigg)\bigg)-\displaystyle\frac{P\bigg( \bigg[ \displaystyle\frac{\tau}{h}\bigg] h \bigg) Q \bigg( \bigg[ \displaystyle\frac{\tau}{h}\bigg] h \bigg)}{a_1+a_2 P\bigg( \bigg[ \displaystyle\frac{\tau}{h}\bigg] h \bigg)} , \\ D^\alpha Q( \tau )= \frac{b_1 P \bigg( \displaystyle\bigg[ \displaystyle\frac{\tau}{h}\bigg] h \bigg) Q\bigg( \bigg[ \displaystyle\frac{\tau}{h}\bigg] h \bigg)}{a_1+a_2 P \bigg( \bigg[ \displaystyle\frac{\tau}{h}\bigg] h \bigg)}-b_2 Q \bigg( \bigg[ \displaystyle\frac{\tau}{h}\bigg] h \bigg) , \end{cases} \end{align} $$

with

$$ \begin{align*} P(0)=P_0,\quad Q(0)=Q_0. \end{align*} $$

The equations in (1.7) are both in the form of (1.1), so the corresponding difference equation, from (1.2), is

(1.8) $$ \begin{align} \begin{cases} M_{n+1}=M_{n}+ \displaystyle\frac{\rho} {\Gamma (\alpha +1)}\bigg(M_{n} (1-M_{n} )-\displaystyle\frac{M_{n} N_{n} }{a_1+a_2 M_{n}}\bigg), \\ N_{n+1}=N_{n}+ \displaystyle\frac{\rho} {\Gamma (\alpha +1)} \bigg(\displaystyle\frac{b_1 M_{n} N_{n} }{a_1+a_2 M_{n}}-b_2 N_{n}\bigg), \end{cases} \end{align} $$

where $\rho =h^{\alpha }$ ,

$$ \begin{align*} P ( nh )=M_n, \quad Q ( nh )=N_n \end{align*} $$

with the initial conditions

$$ \begin{align*} M_0=P_0 \quad \text{and} \quad N_0=Q_0. \end{align*} $$

The presentation in this paper is organized as follows. Section 1 is devoted to the introduction of the stated problem. Preliminaries needed for the paper are given in Section 2. Stability and bifurcation analysis are investigated in Sections 3 and 4, respectively. Numerical examples are stated in Section 5 and, finally, Section 6 presents the conclusion.

2. Preliminaries

In this section, we recall some definitions and lemmas that will be useful in the later sections.

Consider the difference system

$$ \begin{align*} \begin{cases} M_{n+1}=f(M_{n}, N{n}), \\ N_{n+1}=g(M_{n}, N{n}). \end{cases} \end{align*} $$

The Jacobian of this system at any fixed point $(\bar {M},\bar {N})$ is

$$\begin{align*}J = \begin{pmatrix} \dfrac{\partial f}{\partial M_n} & \dfrac{\partial f}{\partial N_n} \\[7pt] \dfrac{\partial g}{\partial M_n} & \dfrac{\partial g}{\partial N_n} \end{pmatrix}\bigg|_{(\bar{M},\bar{N})} , \end{align*}$$

and the corresponding characteristic equation is

(2.1) $$ \begin{align} \lambda^2 - \text{Tr}(J) \cdot \lambda + \det(J) = 0. \end{align} $$

Definition 2.1 [31].

Let $\lambda _1$ and $\lambda _2$ be the characteristic roots of (2.1). A fixed point $(\bar {M},\bar {N})$ is called:

1. (i) sink if $\lvert \lambda _1 \rvert <1$ and $\lvert \lambda _2 \rvert <1$ , and it is locally asymptotically stable;

2. (ii) source if $\lvert \lambda _1 \rvert>1$ and $\lvert \lambda _2 \rvert>1$ , and it is locally unstable;

3. (iii) saddle if $\lvert \lambda _1 \rvert>1$ and $\lvert \lambda _2 \rvert <1$ or $\lvert \lambda _1 \rvert <1$ and $\lvert \lambda _2 \rvert>1$ ;

4. (iv) nonhyperbolic if either $\lvert \lambda _1 \rvert =1$ or $\lvert \lambda _2 \rvert =1$ .

Lemma 2.2 [29].

Let $F(\lambda ) =\lambda ^2 + B\lambda + C,$ where B and C are constants. Suppose $F(1)> 0$ , and $\lambda _{1}$ and $\lambda _{2}$ are two roots of $F(\lambda ) = 0.$ Then:

1. (i) $\lvert \lambda _{1} \rvert < 1$ and $\lvert \lambda _{2} \rvert < 1$ if and only if $F(-1)>0$ and $C<1$ ;

2. (ii) $\lambda _{1}=-1$ and $\lvert \lambda _{2} \rvert \neq 1$ if and only if $F(-1)=0$ , and $B\neq 0,2$ ;

3. (iii) $\lvert \lambda _{1} \rvert> 1$ and $\lvert \lambda _{2} \rvert < 1$ if and only if $F(-1)<0$ ;

4. (iv) $\lvert \lambda _{1} \rvert> 1$ and $\lvert \lambda _{2} \rvert> 1$ if and only if $F(-1)>0$ and $C>1$ ;

5. (v) $\lambda _{1} $ and $ \lambda _{2} $ are the conjugate complex roots and $\lvert \lambda _{1} \rvert = \lvert \lambda _{2} \rvert =1$ if and only if ${B^2-4C<0}$ and $C = 1.$

Definition 2.3 [35].

The bifurcation associated with the appearance of $\mu _{1}=-1$ is called a flip (or period-doubling) bifurcation, where $\mu _{1}$ is an eigenvalue of the nonhyperbolic fixed point of the the system

(2.2) $$ \begin{align} x_{n+1}=f( x_n, \alpha ), \quad x \in \mathbb{R}^{n}, \quad \alpha \in \mathbb{R} \end{align} $$

with smooth f.

Theorem 2.4 [35].

Suppose that a one-dimensional system

$$ \begin{align*} x_{n+1}=f( x_n, \alpha ), \quad x \in \mathbb{R}, \quad \alpha \in \mathbb{R}, \end{align*} $$

with smooth f, has at $\alpha = 0$ the fixed point $ x_0= 0$ and let $\mu =f_x (0,0)=-1.$ Assume that the following nondegeneracy conditions are satisfied:

1. (i) $\tfrac {1}{2}( f_{xx}(0,0))^{2} + \tfrac {1}{3} f_{xxx}(0,0)\neq 0$ ;

2. (ii) $f_{x \alpha }(0,0)\neq 0$ .

Then, there is a flip bifurcation in the system.

Definition 2.5 [35].

The bifurcation corresponding to the presence of $\mu _{1,2}=e^{\pm i \theta _{0}}$ , $0<\theta _{0}<\pi $ , is called a Neimark–Sacker bifurcation, where $\mu _{1,2}$ are the eigenvalues of the nonhyperbolic fixed point of the system (2.2).

Theorem 2.6 [35].

Suppose that the system

$$ \begin{align*} \begin{pmatrix} X_{n+1}\\ Y_{n+1} \end{pmatrix} = \begin{pmatrix} f(X_n, Y_n, \delta ) \\ g(X_n, Y_n, \delta ) \end{pmatrix}, \end{align*} $$

with a parameter $\delta $ , has a pair of complex conjugate eigenvalues $\lambda _{1,2}=r(\delta )e^{\pm i \theta (\delta )}$ at the fixed point $(0,0)$ , where $\delta _c$ is the critical parameter value, $r(\delta _c)=1$ and $\theta (\delta _c)=\theta _0$ . Assume that the following conditions are satisfied:

1. (i) ${d(r(\delta ))}/{d \delta } \rvert _{\delta =\delta _c} \neq 0$ ;

2. (ii) $e^{i k \theta _0} \neq 1$ for $k=1,2,3,4$ and $a(\delta _c)= Re (e^{i\theta _0}c_{1}(\delta _c)) \neq 0$ .

Then, there is a neighbourhood of origin in which a unique closed curve bifurcates from the origin as $\delta $ passes through the critical value $\delta _c$ . It should be mentioned here that the sign of $a(\delta _c)$ determines the direction of the appearance of the invariant curve in a generic system exhibiting the Neimark–Sacker bifurcation: if $a(\delta _c)<0$ , then there is a supercritical Neimark–Sacker bifurcation which is stable; if $a(\delta _c)>0$ , then there is a subcritical Neimark–Sacker bifurcation which is unstable.

3. Stability analysis

The difference system (1.8) has three equilibrium points given by

(3.1) $$ \begin{align} E_0=(0,0), \quad E_1=(1,0) \quad \text{and} \quad E_{*}=( M_{*}, N_{*} ) \end{align} $$

with

(3.2) $$ \begin{align} M_{*}=\frac{a_1 b_2}{b_1-a_2 b_2}, \quad N_{*}=\frac{a_1 b_1 (b_1- b_2(a_1+a_2))}{(b_1-a_2 b_2)^{2}}. \end{align} $$

When simulating biological systems, positivity ensures that the model can make meaningful predictions about real-world behaviour. Negative values in simulations would lead to outcomes that are not interpretable or applicable to the biological phenomenon being studied. To ensure the positivity of the point $E_{*}=( M_{*}, N_{*} )$ , it is assumed that

(3.3) $$ \begin{align} b_1> b_2(a_1+a_2). \end{align} $$

Now, let us investigate the stability of each point.

3.1. Stability of $E_0$

Theorem 3.1. $E_0=(0,0)$ is:

1. (i) saddle if $0<\rho < ({2\Gamma (\alpha +1)})/{b_{2}}$ ;

2. (ii) source if $\rho> ({2\Gamma (\alpha +1)})/{b_{2}}$ ,

where $\rho =h^\alpha $ , h is the discretization parameter and $\alpha $ is the order of the equation.

Proof. The Jacobian matrix at $E_{0}$ is

(3.4) $$ \begin{align} J(E_{0})= \begin{array}{cc} \left(\begin{array}{cc} \displaystyle\frac{\rho}{\Gamma (\alpha +1)}+1 & 0 \\ 0 & 1-\displaystyle\frac{\rho b_2}{\Gamma (\alpha +1)} \end{array}\right) \end{array}\!\!. \end{align} $$

The characteristic polynomial corresponding to (3.4) is calculated as

$$ \begin{align*} F(\lambda)=\lambda ^2+ \bigg(\frac{(b_2-1) \rho}{\Gamma (\alpha +1)}-2 \bigg)\lambda + \frac{(\Gamma (\alpha +1)+\rho ) (\Gamma (\alpha +1)-b_2\rho)}{\Gamma (\alpha +1)^2}, \end{align*} $$

whose eigenvalues are

$$ \begin{align*} \lambda_1= 1-\frac{b_2 \rho}{\Gamma (\alpha +1)},\quad \lambda_2=1+\frac{\rho}{\Gamma (\alpha +1)}>1. \end{align*} $$

Hence, we have

$$ \begin{align*} 1-\frac{b_2 \rho}{\Gamma (\alpha +1)} \begin{cases} \in (-1,1)&if \ 0<\rho <\dfrac{2\Gamma (\alpha +1)}{b_{2}}, \\[6pt] <-1&if \ \rho>\dfrac{2\Gamma (\alpha +1)}{b_{2}}. \end{cases} \end{align*} $$

From Definition 2.1 and Lemma 2.2, the proof is complete.

3.2. Stability of $E_1$

Theorem 3.2. Assume that $b_1<(a_1+a_2) b_2.$ Then, $E_{1}=(1,0)$ is:

1. (i) sink if $0<\rho <\min \{\delta _{1}, \delta _{2}\}$ ;

2. (ii) saddle if $\min \{\delta _{1}, \delta _{2}\}<\rho <\max \{\delta _{1}, \delta _{2}\}$ ;

3. (iii) source if $\rho>\max \{\delta _{1}, \delta _{2}\},$

where

$$ \begin{align*} \delta_1=2\Gamma (\alpha +1),\quad \delta_2=\frac{2 (a_1+a_2) \Gamma (\alpha +1)}{(a_1+a_2) b_2-b_1}. \end{align*} $$

Proof. The Jacobian matrix at $E_{1}$ is obtained as

$$ \begin{align*} J(E_{1})=\left( \begin{array}{@{}cc@{}} 1-\dfrac{\rho}{\Gamma (\alpha +1)} & - \dfrac{\rho}{\Gamma (\alpha +1) (a_1+a_2)} \\[10pt] 0 & 1 + \dfrac{(b_1-(a_1+a_2)b_2) \rho}{(a_1+a_2)\Gamma (\alpha +1)} \\ \end{array} \right)\!, \end{align*} $$

which has the characteristic polynomial

$$ \begin{align*} F(\lambda)&=\lambda ^2+ \bigg(\frac{((a_1+a_2) b_2+a_1+a_2-b_1)\rho }{(a_1+a_2) \Gamma (\alpha +1)}-2\bigg)\lambda\\ &\quad+\frac{(\Gamma (\alpha +1)-\rho ) ((a_1+a_2) (\Gamma (\alpha +1)-b_2 \rho )+b_1 \rho )}{(a_1+a_2) \Gamma (\alpha +1)^2}. \end{align*} $$

Hence, the following eigenvalues are found:

$$ \begin{align*} \lambda_{1}=1-\frac{\rho}{\Gamma (\alpha +1)},\quad \lambda_{2}= 1-\frac{((a_1+a_2) b_2-b_1) \rho}{(a_1+a_2) \Gamma (\alpha +1)}. \end{align*} $$

Since $b_1<(a_1+a_2) b_2$ ,

$$ \begin{align*} 1-\frac{\rho}{\Gamma (\alpha +1)}\begin{cases} \in (-1,1)& \text{if } 0<\rho <2\Gamma (\alpha +1), \\ <-1& \text{if } \rho>2\Gamma (\alpha +1), \end{cases} \end{align*} $$

and

$$ \begin{align*} 1-\dfrac{((a_1+a_2) b_2-b_1) \rho}{(a_1+a_2) \Gamma (\alpha +1)} \begin{cases} \in (-1,1)& \text{if } 0<\rho <\dfrac {2(a_1+a_2) \Gamma (\alpha +1)}{((a_1+a_2) b_2-b_1) \rho}, \\[10pt] <-1&\text{if } \rho>\dfrac {2(a_1+a_2) \Gamma (\alpha +1)}{((a_1+a_2) b_2-b_1) \rho}. \end{cases} \end{align*} $$

So, Definition 2.1 and Lemma 2.2 complete the proof.

3.3. Stability of $E_*$

Theorem 3.3. Assume that $ b_1(a_1 - a_2) + a_2 b_2(a_1 + a_2)> 0$ and (3.3) is true. Then, the positive equilibrium point $E_{*}$ of system (1.8) is:

1. (i) sink fixed point if one of the following conditions holds:

   $$ \begin{align*} 0<\rho<\rho_1\quad \textit{when}\ \xi_2> 4\xi_1;\\ 0<\rho<\frac{1}{\xi_1}\quad \textit{when} \ \xi_2\leq 4\xi_1, \end{align*} $$

2. (ii) source fixed point if one of the following conditions holds:

   $$ \begin{align*} \rho>\rho_2\quad \textit{when} \ \xi_2> 4\xi_1;\\ \rho>\frac{1}{\xi_1}\quad \textit{when} \ \xi_2\leq 4\xi_1, \end{align*} $$

3. (iii) saddle fixed point if the following condition holds:

   $$ \begin{align*} \rho_1<\rho<\rho_2\quad \textit{when} \ \xi_2> 4\xi_1, \end{align*} $$

where $\rho =h^\alpha $ and

(3.5) $$ \begin{align} \begin{aligned} \xi_1&=\frac{({b_1}-{a_2} {b_2}) ({b_1}-{b_2} ({a_1}+{a_2}))}{\Gamma (\alpha +1) ({b_1} ({a_1}-{a_2})+{a_2} {b_2} ({a_1}+{a_2}))},\\ \xi_2&=\frac{b_2 (b_1 (a_1-a_2)+a_2 b_2 (a_1+a_2))}{ \Gamma (\alpha +1)b_1 (b_1-a_2 b_2)},\\ \rho_1&=\frac{\xi_2-\sqrt{\xi_2^{2}-4\xi_1\xi_2}}{\xi_1\xi_2},\\ \rho_2&=\frac{\xi_2+\sqrt{\xi_2^{2}-4\xi_1\xi_2}}{\xi_1\xi_2}. \end{aligned} \end{align} $$

Proof. The Jacobian matrix at $E_{*}$ , given in (3.1) and (3.2), is

(3.6) $$ \begin{align} J_*=\left( \begin{array}{@{}cccc@{}} 1-\dfrac{\rho b_2 ((a_1-a_2) b_1+a_2 (a_1+a_2) b_2)}{\Gamma (\alpha +1) b_1 (b_1-a_2 b_2)} & -\dfrac{\rho b_2}{\Gamma (\alpha +1) b_1} \\[8pt] \dfrac{\rho (b_1-(a_1+a_2) b_2)}{\Gamma (\alpha +1)} & 1 \\ \end{array} \right), \end{align} $$

whose characteristic equation is

(3.7) $$ \begin{align} \lambda^2+\Phi \lambda+\Psi=0, \end{align} $$

where

$$ \begin{align*} \Phi&=\bigg(\frac{b_2 ((a_1-a_2) b_1+a_2 (a_1+a_2) b_2) \rho}{b_1 (b_1-a_2 b_2) \Gamma (\alpha +1)}-2\bigg), \\ \Psi&=1+\frac{b_2 \rho (\rho (b_1-a_2 b_2) (b_1-b_2 (a_1+a_2))-\Gamma (\alpha +1) (b_1 (a_1-a_2)+a_2 b_2 (a_1+a_2)))}{b_1 \Gamma (\alpha +1)^2 (b_1-a_2 b_2)}. \end{align*} $$

To apply Lemma 2.2(i), consider the left-hand side of (3.7) as a function of $\lambda $ , that is,

(3.8) $$ \begin{align} F(\lambda)=\lambda^2+\Phi \lambda+\Psi, \end{align} $$

or

(3.9) $$ \begin{align} F(\lambda)=\lambda^2+( \rho \xi_2-2)\lambda+(\rho^2 \xi_1 \xi_2- \rho \xi_2 +1) \end{align} $$

with $h^\alpha =\rho $ , and $\xi _1\: and \: \xi _2 $ are given in (3.5). Equation (3.8) is called the characteristic polynomial corresponding to the positive equilibrium point $E_{*}$ .

Taking $\lambda =1$ and $\lambda =-1$ above, and considering $\Phi $ and $ \Psi $ , $F(1)$ and $F(-1)$ are calculated as

$$ \begin{align*} F(1)&=\frac{b_2 (b_1-(a_1+a_2) b_2) \rho^{2}}{b_1 \Gamma (\alpha +1)^2}, \\ F(-1)&=\frac{b_2(b_1-(a_1+a_2)b_2)}{b_1 \Gamma^2 (\alpha +1)}\rho^{2} - \frac{b_2((a_1-a_2)b_1+a_2(a_1+a_2)b_2)}{b_1(b_1-a_2b_2)\Gamma (\alpha +1)}\rho+4. \end{align*} $$

From (3.3), it is clear that $F(1)>0$ . To ensure the positivity of $F(-1)$ , it is written as

$$ \begin{align*} F(-1)&=\xi_1 \xi_2\rho^2-2\xi_2\rho+4\\ &=(\rho-\rho_1)(\rho-\rho_2), \end{align*} $$

where $\xi _1, \: \xi _2, \: \rho _1, \: \rho _2$ are given in (3.5). So, the roots of $F(-1)$ are obtained as

$$ \begin{align*} \rho=\rho_1=\frac{\xi_2-\sqrt{\xi_2^{2}-4\xi_1\xi_2}}{\xi_1\xi_2} \quad \text{and} \quad \rho=\rho_2=\frac{\xi_2+\sqrt{\xi_2^{2}-4\xi_1\xi_2}}{\xi_1\xi_2}. \end{align*} $$

First, we assume that $\xi _2>4\xi _1$ , then $\rho _1$ and $\rho _2$ are real and distinct. Under conditions (3.3) and $a_1>a_2$ , it is obtained that $\rho _1$ and $\rho _2$ are positive with $0<\rho <\rho _1$ or $\rho>\rho _2$ , and this implies that $F(-1)>0$ . Moreover, when $0<\rho <\rho _1,$ we get $\rho <\rho _1<{1}/{\xi _1}<\rho _2$ and $C=\xi _1\xi _2\rho ^2-\xi _2\rho +1<1.$ When $\rho>\rho _2,$ it is found $\rho>{1}/{\xi _1}$ and then $C>1$ .

Second, if $\:\xi _2=4\xi _1$ and $\:\xi _2<4\xi _1$ , then $F(-1)$ has real roots $ \rho _{1}= \rho _{2}={1}/{\xi _1}$ and does not have real roots, respectively, and for each case, $F(-1)$ becomes positive. Also, if $0<\rho <{1}/{\xi _1},$ then $C<1$ or if $\rho>{1}/{\xi _1},$ then $C>1.$

So, from Definition 2.1 and Lemma 2.2, the proof is complete.

4. Bifurcation analysis

This section is devoted to the investigation of the flip and Neimark–Sacker bifurcation analysis of the equilibrium point $E_*$ given in (3.1)–(3.2). Here, we choose $\rho $ as a bifurcation parameter.

4.1. Flip bifurcation

In this part, we study flip bifurcation for the system (1.8). We prove that under certain conditions, the flip bifurcation arises from the positive fixed point $E_*$ .

Theorem 4.1. If the following conditions are satisfied by the characteristic polynomial (3.8), then one may conclude that there is a flip bifurcation at $E_*$ . When $\xi _2> 4\xi _1$ , then $F(-1)$ has the roots $ \rho =\rho _1$ and $\rho =\rho _2$ such that

$$ \begin{align*} \rho \neq \frac{2}{\xi_2} \quad \textrm{and} \quad \rho \neq \frac{4}{\xi_2}, \end{align*} $$

where $\rho _1$ and $\rho _2$ are given in (3.5).

Proof. From Lemma 2.2(ii), it is known that if

$$ \begin{align*} F(-1)=0 \quad \text{and} \quad \rho \xi_2-2 \neq \{0,2\}, \end{align*} $$

then there may be flip bifurcation at the equilibrium point $E_*$ . From Theorem 3.3, it is known that when $\xi _2> 4\xi _1$ , then $F(-1)$ has two real roots, which means that if we take $\rho =\rho _1$ or $\rho _2$ , then we get $F(-1)=0$ . Moreover, considering $ \rho \neq {2}/{\xi _2}$ and $ \rho \neq {4}/{\xi _2}$ leads us to $\rho \xi _2-2 \neq \{0,2\},$ and thus the proof is complete.

However, let us take

(4.1) $$ \begin{align} \rho=\rho_1=\frac{\xi_2-\sqrt{\xi_2^{2}-4\xi_1\xi_2}}{\xi_1\xi_2} , \end{align} $$

then the characteristic roots of (3.9) are found as $\lambda _1=-1$ and $\lambda _2=3-\rho _1 \xi _2$ , with $\lvert \lambda _{2} \rvert \neq 1$ .

We will now focus on the analytic construction of the flip bifurcation, or, to put it another way, use Theorem 2.4 to establish the following theorem (Theorem 4.2) at the point $E_*$ . However, to apply Theorem 2.4, we must first apply to our system (1.8) the centre manifold theorem [9, 35, 53].

Theorem 4.2. If the following conditions are satisfied, then there is a flip bifurcation for the system (1.8):

1. (i) $\xi _2>4 \xi _1$ ;

2. (ii) $2 m_1+2( h_1 m_2 + m_5 ) \neq 0$ ;

3. (iii) $m_3 \neq 0$ .

Furthermore, if $2 m_1+2( h_1 m_2 + m_5 )$ is positive, then the flip bifurcation is supercritical, otherwise it becomes subcritical. Here, $\xi _1$ and $\xi _2$ are given in (3.5), and $h_1, m_1, m_2, m_3, m_5$ are given in the following proof in (4.6) and (4.8).

Proof. For the application of the centre manifold theorem, let us define the right-hand side of the equations in system (1.8) as

$$ \begin{align*} \begin{cases} F(M_{n},N_{n})=M_{n}+ \dfrac{\rho} {\Gamma (\alpha +1)}\bigg(M_{n} (1-M_{n} )-\dfrac{M_{n} N_{n} }{a_1+a_2 M_{n}}\bigg),\\ G(M_{n},N_{n})=N_{n}+ \dfrac{\rho} {\Gamma (\alpha +1)}\bigg(\dfrac{b_1 M_{n} N_{n} }{a_1+a_2 M_{n}}-b_2 N_{n}\bigg). \end{cases} \end{align*} $$

Expanding the Taylor series of these functions around the fixed point

$$ \begin{align*} E_{*}=( M_{*}, N_{*} )= \bigg(\frac{a_1 b_2}{b_1-a_2 b_2}, \frac{a_1 b_1 (b_1- b_2(a_1+a_2))}{(b_1-a_2 b_2)^{2}}\bigg) \end{align*} $$

gives us the following equivalent system of (1.8):

(4.2) $$ \begin{align} \begin{pmatrix} X_n\\ Y_n \end{pmatrix} =\begin{pmatrix} B_{11} & B_{12} \\ B_{21} & 1 \end{pmatrix} \begin{pmatrix} X_n\\ Y_n \end{pmatrix} + \begin{pmatrix} f_1(X_n,Y_n,\tilde{\rho},\bar{\rho})\\ g_1(X_n,Y_n,\tilde{\rho},\bar{\rho}) \end{pmatrix} , \end{align} $$

where

(4.3) $$ \begin{align} \begin{aligned} B_{11}&=\frac{b_{1} \Gamma (\alpha +1) (b_{1}-a_{2} b_{2})-b_{2} \rho (a_{1} (a_{2} b_{2}+b_{1})+a_{2} (a_{2} b_{2}-b_{1}))}{b_{1} \Gamma (\alpha +1) (b_{1}-a_{2} b_{2})},\\ B_{12}&=-\frac{b_{2} \rho }{b_{1} \Gamma (\alpha +1)},\\ B_{21}&=\frac{\rho (b_{1}-b_{2} (a_{1}+a_{2}))}{\Gamma (\alpha +1)}, \end{aligned} \end{align} $$

and

(4.4) $$ \begin{align} \begin{aligned} f_1(X_n,Y_n,\tilde{\rho},\bar{\rho})&= \frac{(\bar{\rho}+\tilde{\rho }) (b_{1}-a_{2} b_{2})^2}{a_{1} b_{1}^2 \Gamma (\alpha +1)} X_{n} Y_{n}-\frac{2 a_{2} (\bar{\rho}+\tilde{\rho }) (a_{2} b_{2}-b_{1})^3}{a_{1}^2 b_{1}^3 \Gamma (\alpha +1)} X_{n}^{2}Y_{n}\\ & \quad -\frac{2 (\bar{\rho}+\tilde{\rho }) (a_{1} (-a_{2}^2 b_{2}^2+a_{2} b_{1} b_{2}+b_{1}^2)-a_{2} (b_{1}-a_{2} b_{2})^2)}{a_{1} b_{1}^2 \Gamma (\alpha +1)} X_{n}^{2}\\ & \quad -\frac{6 a_{2}^2 (\bar{\rho}+\tilde{\rho }) (b_{1}-a_{2} b_{2})^2 (b_{2} (a_{1}+a_{2})-b_{1})}{a_{1}^2 b_{1}^3 \Gamma (\alpha +1)} X_{n}^{3}+H.O.T.,\\ g_1(X_n,Y_n,\tilde{\rho},\bar{\rho})&=\frac{(\bar{\rho}+\tilde{\rho }) (b_{1}-a_{2} b_{2})^2}{a_{1} b_{1} \Gamma (\alpha +1)}X_{n}Y_{n} +\frac{2 a_{2} (\bar{\rho}+\tilde{\rho }) (a_{2} b_{2}-b_{1})^3}{a_{1}^2 b_{1}^2 \Gamma (\alpha +1)} X_{n}^{2}Y_{n}\\ & \quad +\frac{2 a_{2} (\bar{\rho}+\tilde{\rho }) (b_{1}-a_{2} b_{2}) (b_{2} (a_{1}+a_{2})-b_{1})}{a_{1} b_{1} \Gamma (\alpha +1)}X_{n}^{2}\\ & \quad -\frac{6 a_{2}^2 (\bar{\rho}+\tilde{\rho }) (b_{1}-a_{2} b_{2})^2 (b_{2} (a_{1}+a_{2})-b_{1})}{a_{1}^2 b_{1}^2 \Gamma (\alpha +1)} X_{n}^{3}+H.O.T., \end{aligned} \end{align} $$

with $\rho =\rho _{1}$ as in (4.1), where $\xi _1$ , $\xi _2$ are given in (3.5), and H.O.T. represents the higher order terms.

In this expansion, we take $\bar {\rho }=\rho _1$ and $M_{n}-M_{*}=X_{n}, \: N_{n}-N_{*}=Y_{n}, \: \rho -\bar {\rho }=\tilde {\rho }$ , which transform the fixed point $E_*$ to the origin and the bifurcation parameter’s critical value to zero, that is, $\rho _c=0$ . Next, we construct the matrix

$$ \begin{align*} V= \begin{pmatrix} -2 && 2-\bar{\rho} \xi_2 \\ B_{21} && B_{21} \end{pmatrix} , \end{align*} $$

whose columns are the eigenvectors corresponding to the eigenvalues of $J_*$ , given in (3.6), when $\rho =\bar {\rho }=\rho _1$ , which are $ \lambda _1=-1 \quad \text {and} \quad \lambda _2 = 3-\rho _1 \xi _2.$ Here, $\rho _{1}$ is given in (4.1) with $\xi _1$ and $\xi _2$ as in (3.5), and $B_{21}$ is given in (4.3). Using the transformation

$$ \begin{align*} \begin{pmatrix} X_n\\ Y_n \end{pmatrix}=V \begin{pmatrix} U_n\\ V_n \end{pmatrix} \end{align*} $$

in (4.2), the following system is obtained:

(4.5) $$ \begin{align} \begin{pmatrix} U_{n+1}\\ V_{n+1} \end{pmatrix}= \begin{pmatrix} -1 && 0\\ 0 && 3-\rho_1 \xi_2\end{pmatrix} \begin{pmatrix} U_n\\ V_n \end{pmatrix}+ \begin{pmatrix} f_2(U_n,V_n)\\ g_2(U_n,V_n)\end{pmatrix} \quad \text{with} \end{align} $$

$$ \begin{align*} f_2(U_n,V_n)&=m_1 U_n^{3}+m_2 U_n V_n^{2} + m_3 U_n^{2} V_n+m_4 V_n^{3}+m_5 U_n^{2}\\&\quad+m_6 U_n V_n+m_7 V_n^{2}+H.O.T.,\\ g_2(U_n,V_n)&=n_1 U_n^{3}+n_2 U_n V_n^{2} + n_3 U_n^{2} V_n+n_4 V_n^{3}+n_5 U_n^{2}+n_6 U_n V_n+n_7 V_n^{2}+H.O.T., \end{align*} $$

where

$$ \begin{align*} m_1&=(a_1^2 b_1^3 \Gamma (\alpha +1))^{-1}8 a_2 (\bar{\rho}+\tilde{\rho }) (b_1-a_2 b_2)^2 (b_1 T_{12}-T_{11})(a_2 (b_2 (6 a_1+B_{21})-6 b_1)\nonumber\\ &\quad +6 a_2^2 b_2-b_1 B_{21}), \nonumber\\ m_2&=(a_1^2 b_1^3 \Gamma (\alpha +1))^{-1}(a_1^2 b_1^3 \Gamma (\alpha +1))^{-1}2 a_2 (\bar{\rho }+\tilde{\rho }) (r_2 \bar{\rho }-2) (b_1-a_2 b_2)^2 (b_1 T_{12}-T_{11})\nonumber\\ &\quad \times (r_2 \bar{\rho } (a_2 (18 a_1 b_2-18 b_1+b_2 B_{21})+18 a_2^2 b_2-b_1 B_{21}+2 a_2 (-18 a_1 b_2+18 b_1\nonumber\\ &\quad +b_2 B_{21})-36 a_2^2 b_2-2 b_1 B_{21})a_1+B_{21}))-18 a_2^2 b_2+b_1 B_{21}), \nonumber\\ m_3&=(a_1^2 b_1^3 \Gamma (\alpha +1))^{-1}8 a_2 (\bar{\rho }+\tilde{\rho }) (b_1-a_2 b_2)^2 (b_1 T_{12}-T_{11}) (\bar{\rho }+\tilde{\rho }) (r_2 \bar{\rho } (a_2(b_2 (9 a_1\nonumber\\ &\quad +B_{21})-9 b_1)+9 a_2^2 b_2-b_1 B_21)+a_2(18 b_1-b_2 (18 a_1+B_{21})) -18 a_{2}^2 b_2+b_1 B_{21}), \nonumber\\ m_4&=(a_1^2 b_1^3 \Gamma (\alpha +1))^{-1}(2 a_2 (\bar{\rho }+\tilde{\rho }) (r_2 \bar{\rho }-2)^2 (b_1-a_2 b_2)^2 (b_1 T_{12}-T_{11})\nonumber\\ &\quad \times (3 a_2 r_2 \bar{\rho } (b_2 (a_1+a_2)-b_1)+a_2 (-6 a_1 b_2+6 b_1+b_2 B_{21})-6 a_2^2b_2-b_1 B_{21})),\nonumber\\ m_5&=-(a_{1} b_{1}^2 \Gamma (\alpha +1))(2 (\bar{\rho }+\tilde{\rho }) (4a_{1}(-a_{2}^2 b_{2}^2 T_{11}+b_{1}^2 (T_{11}-a_{2} b_{2} T_{12})\nonumber\\ &\quad +a_{2} b_{1} b_{2} (a_{2} b_{2} T_{12}+T_{11}))+(4 a_{2}+B_{21}) (b_{1}-a_{2} b_{2})^2 (b_{1} T_{12}-T_{11}))), \nonumber\\ m_6&=-(a_1 b_1^2 \Gamma (\alpha +1))^{-1}((\bar{\rho }+\tilde{\rho }) (r_2 \bar{\rho } (8 a_1 (-a_2^2 b_2^2 T_{11}+b_1^2 (T_{11}-a_2 b_2 T_{12})\nonumber\\ &\quad +a_2 b_1 b_2 (a_2 b_2 T_{12}+T_{11}))+(8 a_2+B_{21})(b_1-a_2 b_2)^2 (b_1 T_{12}-T_{11}))\nonumber\\ &\quad -16 (a_1 (-a_2^2 b_2^2 T_{11}+b_1^2 (T_{11}-a_2 b_2 T_{12})+a_2 b_1 b_2 (a_2 b_2 T_{12}+T_{11}))\nonumber\\ &\quad +a_2 (b_1-a_2 b_2)^2 (b_1 T_{12}-T_{11})))),\nonumber\\ m_7&=-(a_{1} b_{1}^2 \Gamma (\alpha +1))^{-1}((\bar{\rho }+\tilde{\rho })(\xi_{2} \bar{\rho }-2) (2 \xi_{2} \bar{\rho } (a_{1} (-a_{2}^2 b_{2}^2 T_{11}+b_{1}^2 (T_{11}-a_{2} b_{2} T_{12})\nonumber\\ &\quad +a_{2} b_{1} b_{2} (a_{2} b_{2} T_{12}+T_{11}))+a_{2} (b_{1}-a_{2} b_{2})^2 (b_{1} T_{12}-T_{11})),\nonumber\\ &\quad -4a_{1}(-a_{2}^2 b_{2}^2 T_{11}+b_{1}^2 (T_{11}-a_{2} b_{2} T_{12})+a_{2} b_{1} b_{2} (a_{2} b_{2} T_{12}+T_{11}))\nonumber\\ &\quad+(4 a_{2}-B_{21}) (-(b_{1}-a_{2} b_{2})^2) (b_{1} T_{12}-T_{11}))),\end{align*} $$

(4.6) $$ \begin{align} n_1&=(a_1^2 b_1^3 \Gamma (\alpha +1))^{-1}(8 a_2 (\bar{\rho }+\tilde{\rho }) (b_1-a_2 b_2)^2 (b_1 T_{22}-T_{21}) (a_2 (b_2 (6 a_1+B_{21})-6 b_1)\nonumber\\ &\quad +6 a_2^2 b_2-b_1 B_{21})), \nonumber\\ n_2&=(a_1^2 b_1^3 \Gamma (\alpha +1))^{-1}(2 a_2 (\bar{\rho }+\tilde{\rho }) (r_2 \bar{\rho }-2) (b_1-a_2 b_2)^2 (b_1 T_{22}-T_{21}), \nonumber\\ &\quad \times (r_2 \bar{\rho }(a_2 (18 a_1 b_2-18 b_1+b_2 B_{21})+18 a_2^2 b_2-b_1 B_{21})\nonumber\\ &\quad +2 a_2 (-18 a_1 b_2+18 b_1+b_2 B_{21})-36 a_2^2 b_2-2 b_1 B_{21})),\nonumber\\n_3&=(a_1^2 b_1^3\Gamma (\alpha +1))^{-1}(8 a_2 (\bar{\rho }\!+\!\tilde{\rho }) (b_1\!-\!a_2 b_2)^2 (b_1 T_{22}\!-\! T_{21}) (r_2 \bar{\rho } (a_2 (b_2 (9 a_1+B_{21})-9 b_1)\nonumber\\ &\quad +9 a_2^2 b_2-b_1 B_{21})+a_2 (18 b_1-b_2 (18 a_1+B_{21}))-18 a_2^2 b_2+b_1 B_{21})),\nonumber\\ n_4&=(a_1^2 b_1^3 \Gamma (\alpha +1))^{-1}(2 a_2 (\bar{\rho }+\tilde{\rho }) (r_2 \bar{\rho }-2)^2 (b_1-a_2 b_2)^2 (b_1 T_{22}-T_{21})\nonumber\\ &\quad \times (3 a_2 r_2 \bar{\rho } (b_2 (a_1+a_2)-b_1)+a_2 (-6 a_1 b_2+6 b_1+b_2 B_{21})-6 a_2^2 b_2-b_1 B_{21})), \nonumber\\ n_5&=-((a_{1} b_{1}^2 \Gamma (\alpha +1))^{-1}(2 (\bar{\rho }+\tilde{\rho }) (4a_{1} (-a_{2}^2 b_{2}^2 T_{21}+b_{1}^2 (T_{21}-a_{2} b_{2} T_{22})\nonumber\\ &\quad +a_{2} b_{1} b_{2} (a_{2} b_{2} T_{22}+T_{21}))+(4 a_{2}+B_{21}) (b_{1}-a_{2} b_{2})^2 (b_{1} T_{22}-T_{21}))),\nonumber\\ n_6&=-(a_1 b_1^2 \Gamma (\alpha +1))^{-1}((\bar{\rho }+\tilde{\rho }) (r_2 \bar{\rho } (8 a_1 (-a_2^2 b_2^2 T_{21}+b_1^2 (T_{21}-a_2 b_2 T_{22})\nonumber\\ &\quad +a_2 b_1 b_2 (a_2 b_2 T_{22}+T_{21}))+(8 a_2+B_{21}) (b_1-a_2 b_2)^2 (b_1 T_{22}-T_{21}))\nonumber\\ &\quad -16 (a_1 (-a_2^2 b_2^2 T_{21}+b_1^2 (T_{21}-a_2 b_2 T_{22})+a_2 b_1 b_2 (a_2 b_2 T_{22}+T_{21}))\nonumber\\ &\quad +a_2 (b_1-a_2 b_2)^2 (b_1 T_{22}-T_{21})))),\nonumber\\ n_7&=-(a_{1} b_{1}^2 \Gamma (\alpha +1))^{-1}((\bar{\rho }+\tilde{\rho }) (\xi_{2} \bar{\rho }-2) (2 \xi_{2} \bar{\rho } (a_{1} (-a_{2}^2 b_{2}^2 T_{21}+b_{1}^2 (T_{21}-a_{2} b_{2} T_{22})\nonumber\\ &\quad +a_{2} b_{1} b_{2} (a_{2} b_{2} T_{22}+T_{21}))+a_{2} (b_{1}-a_{2} b_{2})^2 (b_{1} T_{22}-T_{21}))\nonumber\\ &\quad -4a_{1}(-a_{2}^2 b_{2}^2 T_{21}+b_{1}^2 (T_{21}-a_{2} b_{2} T_{22})+a_{2} b_{1} b_{2} (a_{2} b_{2} T_{22}+T_{21})) \nonumber\\ &\quad +(4 a_{2}-B_{21}) (-(b_{1}-a_{2} b_{2})^2) (b_{1} T_{22}-T_{21}))), \end{align} $$

with

$$ \begin{align*} T_{11}&=\frac{B_{21}}{-4 B_{21}+B_{21}\bar{\rho}\xi_{2}},\\[4pt] T_{12}&=\frac{-2+\bar{\rho}\xi_{2}}{-4 B_{21}+B_{21}\bar{\rho}\xi_{2}}, \\[4pt] T_{21}&=\frac{-B_{21}}{-4 B_{21}+B_{21}\bar{\rho}\xi_{2}}, \\[4pt] T_{22}&=\frac{-2}{-4 B_{21}+B_{21}\bar{\rho}\xi_{2}}, \end{align*} $$

and $B_{21}$ is given in (4.3). Applying centre manifold theorem [9, 35, 53] to the system (4.5), one may obtain the one-dimensional system

(4.7) $$ \begin{align} U_{n+1}=-U_{n}+\tilde{\rho} m_3 U_{n}+(m_1+\tilde{\rho}h_1m_4+\tilde{\rho}m_6)U_{n}^{2} +(h_1m_2+m_5+\tilde{\rho}m_4)U_{n}^{3} \end{align} $$

with

(4.8) $$ \begin{align} h_1=\frac{m_1}{1-\lambda_2} \quad \text{and} \quad h_2=\frac{-m_3}{1+\lambda_2}. \end{align} $$

Now, let us apply Theorem 2.4 to (4.7). Considering the right-hand side of (4.7) as $f( U_n,\tilde {\rho })$ , we calculate that

$$ \begin{align*} &\frac{\partial{f}}{\partial{U_n}}( 0,0)=-1,\\ &\frac{1}{2}\bigg( \frac{\partial^2{f}}{\partial{U_n^2}}(0,0)\bigg)^{2} +\frac{1}{3} \frac{\partial^3{f}}{\partial{U_n^3}}(0,0)=2m_1^{2}+2(h_1 m_2+m_5)\neq 0, \\ &\frac{\partial^2{f}}{\partial{U_n} \partial{\tilde{\rho}}}=m_3 \neq 0. \end{align*} $$

So, from Theorem 2.4, we may conclude that there exists a flip bifurcation for the system (1.8) under the conditions (i), (ii), (iii) in Theorem 4.2. So the proof of Theorem 4.2 is complete.

4.2. Neimark–Sacker bifurcation

This section focuses on the construction of the Neimark–Sacker bifurcation at the point

$$\begin{align*}E_{*}=( M_{*}, N_{*} )=\bigg(\frac{a_1 b_2}{b_1-a_2 b_2}, \frac{a_1 b_1 (b_1- b_2(a_1+a_2))}{(b_1-a_2 b_2)^{2}}\bigg)\end{align*}$$

with the help of condition (v) in Lemma 2.2.

Theorem 4.3. If

(4.9) $$ \begin{align} \rho=\frac{1}{\xi_1} \quad \text{and} \quad \frac{\xi_2}{4}<\xi_1, \end{align} $$

then there exists a Neimark–Sacker bifurcation for the equilibrium point $E_*$ .

Proof. It is known that if the characteristic roots of the Jacobian matrix at $E_*$ are complex conjugate with modulus 1, then there may be a Neimark–Sacker bifurcation at this point. In our problem, the characteristic equation is given by

$$ \begin{align*} \lambda^2-(2- \rho \xi_2)\lambda+(\xi_1 \xi_2 \rho^2 - \xi_2 \rho +1)=0. \end{align*} $$

Under (4.9), conditions given in Lemma 2.2(v) are satisfied and the proof is complete.

To show the existence of the Neimark–Sacker bifurcation, we apply Theorem 2.6. For this aim, let us obtain the characteristic roots and the bifurcation parameter explicitly under the condition $ b_1(a_1 - a_2) + a_2 b_2(a_1 + a_2)> 0.$ The complex characteristic roots are in the form

$$ \begin{align*} \lambda_{1,2}(\rho)=\alpha(\rho)\pm i \beta(\rho)=r(\rho) e^{\pm i \theta(\rho)} \end{align*} $$

with

$$ \begin{align*} \alpha(\rho)&=\frac{2-\xi_2 \rho}{2},\qquad \beta(\rho)=\frac{\rho\sqrt{\xi_2(4 \xi_1-\xi_2)}}{2}, \\ r(\rho)&=\sqrt{ \alpha^{2}(\rho)+\beta^{2}(\rho)},\quad \theta(\rho)=\arctan\bigg(\frac{\beta(\rho)}{\alpha(\rho)}\bigg). \end{align*} $$

Note that when $\rho ={1}/{\xi _1}$ ,

$$ \begin{align*} \lvert\lambda_{1,2}\rvert= r(\rho)=1,\quad \theta(\rho)=\arctan \bigg( \frac{\sqrt{\xi_2(4 \xi_1-\xi_2)}}{2-\xi_2 \rho}\bigg)=\theta_{0}, \end{align*} $$

and our critical bifurcation parameter is $\rho _c={1}/{\xi _1}$ and

$$ \begin{align*} \frac{d(r(\rho))}{d \rho}\bigg \rvert_{\rho_c={1}/{\xi_1}}=\frac{\xi_2}{2}>0 , \end{align*} $$

because of the condition $a_1>a_2$ . Furthermore, since ${\xi _2}>0$ , it is seen that $( \lambda _{1,2}(\rho _c))^k=e^{\pm i k \theta =0} \neq 1$ for $k=1,2,3,4$ . So, we conclude that Neimark–Sacker bifurcation occurs as for the fixed point $E_*$ as $\rho $ passes $\rho _c={1}/{\xi _1}$ from left to right.

To construct the normal form of this bifurcation, let us consider the system (4.2) again as follows:

(4.10) $$ \begin{align} \begin{pmatrix} X_n\\ Y_n \end{pmatrix} = \begin{pmatrix} B_{11} & B_{12} \\ B_{21} & 1 \end{pmatrix} \begin{pmatrix} X_n\\ Y_n \end{pmatrix} + \begin{pmatrix} f_1(X_n,Y_n,\tilde{\rho},\bar{\rho})\\ g_1(X_n,Y_n,\tilde{\rho},\bar{\rho})\end{pmatrix} , \end{align} $$

where $B_{11}, B_{12}, B_{21}$ and $f_1(X_n,Y_n,\tilde {\rho },\bar {\rho })$ , $g_1(X_n,Y_n,\tilde {\rho },\bar {\rho })$ are given in (4.3) and (4.4), respectively. With this system, we transform the fixed point $E_*$ to the origin and the bifurcation parameter $\rho _c=\bar {\rho }$ to zero.

Now, let us consider the following Jacobian matrix:

$$ \begin{align*} J(\bar{\rho})=J\bigg( \frac{1}{\xi_1}\bigg)= \begin{pmatrix} j_{11} & j_{12} \\ j_{21} & j_{22} \end{pmatrix} \quad \text{where} \quad \end{align*} $$

$$ \begin{align*} j_{11}&=\frac{1}{b_{1} (b_{1}-a_{2} b_{2})^2 (b_{1}-b_{2} (a_{1}+a_{2}))}(b_{1}^2 b_{2} (-a_{1}^2+2 a_{1} a_{2} (b_{2}+1)+a_{2}^2 (3 b_{2}-1))\\ &\quad -a_{2}^2 b_{2}^3 (a_{1}+a_{2})^2-b_{1}^3 b_{2} (a_{1}+3 a_{2})-a_{2} b_{1} b_{2}^2 (a_{1}+a_{2}) (2 a_{1}+a_{2} (b_{2}-2))+b_{1}^4),\\ j_{12}&= -\frac{b_{2} (b_{1} (a_{1}-a_{2})+a_{2} b_{2} (a_{1}+a_{2}))}{b_{1} (b_{1}-a_{2} b_{2}) (b_{1}-b_{2} (a_{1}+a_{2}))},\\ j_{21} &=\frac{b_{1} (a_{1}-a_{2})+a_{2} b_{2} (a_{1}+a_{2})}{b_{1}-a_{2} b_{2}}, \\ j_{22}&=1. \end{align*} $$

We now create a matrix

$$ \begin{align*} V=\begin{pmatrix} -j_{12} && 0 \\ j_{11}- \mu && \omega \end{pmatrix} , \end{align*} $$

whose columns are made up of the real and imaginary parts of an eigenvector that correspond to $\lambda _1=\mu -i \omega $ , where

$$ \begin{align*} \mu=1-\frac{\xi_2}{2\xi_1} \quad \text{and} \quad \omega=\frac{\sqrt{\xi_2( 4 \xi_1 - \xi_2 )}}{2\xi_1}. \end{align*} $$

Next, we apply the transformation

$$ \begin{align*} \begin{pmatrix} X_n\\ Y_n \end{pmatrix}=V \begin{pmatrix} U_n\\ V_n \end{pmatrix} \end{align*} $$

to the system (4.10) and procure the following system:

$$ \begin{align*} \begin{pmatrix} U_{n+1}\\ V_{n+1} \end{pmatrix}= \begin{pmatrix} \mu && -\omega \\ \omega && \mu \end{pmatrix} \begin{pmatrix} U_{n}\\ V_{n} \end{pmatrix}+ \begin{pmatrix} H^{1} ( U_{n}, V_{n})\\ H^{2} ( U_{n}, V_{n}) \end{pmatrix} , \end{align*} $$

where

$$ \begin{align*} H^{1} ( U_{n}, V_{n})&=h_{11} U_{n}^{3}+h_{12} U_{n}^{2} V_{n}+h_{13} U_{n}^{2}+h_{14} U_{n} V_{n} ,\\ \nonumber H^{2} ( U_{n}, V_{n})&=h_{21} U_{n}^{3}+h_{22} U_{n}^{2} V_{n}+h_{23} U_{n}^{2}+h_{24} U_{n} V_{n} , \end{align*} $$

with

$$ \begin{align*} h_{11}&=-\frac{1}{a_{1}^2 b_{1}^6 \xi_{1}^4 \Gamma (\alpha +1)^4}(2 a_{2} b_{2}^2 (b_{1}-a_{2} b_{2})^2 (b_{2} (a_{1} (b_{1}-2 a_{2} b_{2})+2 a_{2} (b_{1}-a_{2} b_{2}))\\ &\quad +b_{1} (\mu -1) \xi_{1} \Gamma (\alpha +1) (b_{1}-a_{2} b_{2}))),\\ h_{12}&=\frac{2 a_{2} b_{2}^2 \omega (b_{1}-a_{2} b_{2})^3}{a_{1}^2 b_{1}^5 \xi_{1}^3 \Gamma (\alpha +1)^3},\\ h_{13}&=-\frac{b_{2}}{a_{1} b_{1}^4 \xi_{1}^3 \Gamma (\alpha +1)^3}((b_{2} (a_{1} (a_{2}^2 b_{2}^2-2 a_{2} b_{1} b_{2}+3 b_{1}^2)+a_{2} (b_{1}-a_{2} b_{2})^2)\\ &\quad +b_{1} (\mu -1) \xi_{1} \Gamma (\alpha +1) (b_{1}-a_{2} b_{2})^2)),\\ h_{14}&=-\frac{b_{2} \omega (b_{1}-a_{2} b_{2})^2}{a_{1} b_{1}^3 \xi_{1}^2 \Gamma (\alpha +1)^2} , \end{align*} $$

and

$$ \begin{align*} h_{21}&=\frac{1}{a_{1}^2 b_{1}^5 \xi_{1}^4 \Gamma (\alpha +1)^4}(2 a_{2} b_{2}^2 (b_{1}-a_{2} b_{2})^2 (b_{2} (a_{1} (b_{1}-2 a_{2} b_{2})+2 a_{2} (b_{1}-a_{2} b_{2}))\\ &\quad +b_{1} (\mu -1) \xi_{1} \Gamma (\alpha +1) (b_{1}-a_{2} b_{2}))),\\ h_{22}&=\frac{2 a_{2} b_{2}^2 \omega (a_{2} b_{2}-b_{1})^3}{a_{1}^2 b_{1}^4 \xi_{1}^3 \Gamma (\alpha +1)^3},\\ h_{23}&=\frac{b_{2}}{a_{1} b_{1}^3 \xi_{1}^3 \Gamma (\alpha +1)^3}((b_{2} (a_{1} (b_{1}-a_{2} b_{2})^2+a_{2} (a_{2}^2 b_{2}^2-2 a_{2} b_{1} b_{2}-3 b_{1}^2))\\ &\quad +b_{1} (\mu -1) r_1 \Gamma (\alpha +1) (b_{1}-a_{2} b_{2})^2)),\\ h_{24}&=\frac{b_{2} \omega (b_{1}-a_{2} b_{2})^2}{a_{1} b_{1}^2 \xi_{1}^2 \Gamma (\alpha +1)^2}. \end{align*} $$

However, the direction of this bifurcation is calculated with the help of the formula [21]

(4.11) $$ \begin{align} a\bigg(\frac{1}{\xi_1}\bigg)=-Re\bigg(\frac{(1-2 \lambda)\bar{\lambda}^{2}}{1-\lambda} \zeta_{20} \zeta_{11}\bigg)- \frac{1}{2}\lvert\zeta_{11}\rvert^{2}- \lvert\zeta_{02}\rvert^{2}+Re(\bar{\lambda} \zeta_{21}) \end{align} $$

with the following functions, which are all calculated at the point $(0,0,{1}/{\xi _1})$ :

$$ \begin{align*} \zeta_{20}&=\tfrac{1}{8}[(H^{1}_{UU}-H^{1}_{VV}+2H^{2}_{UV})+i(H^{2}_{UU}-H^{2}_{VV}-2H^{1}_{UV})],\\ \zeta_{11}&=\tfrac{1}{4}[(H^{1}_{UU}+H^{1}_{VV})+i(H^{2}_{UU}+H^{2}_{VV})],\\ \zeta_{02}&=\tfrac{1}{8}[(H^{1}_{UU}-H^{1}_{VV}-2H^{2}_{UV})+i(H^{2}_{UU}-H^{2}_{VV}+2H^{1}_{UV})],\\ \zeta_{21}&=\tfrac{1}{16}[(H^{1}_{UUU}+H^{1}_{UVV}+H^{2}_{UUV}+H^{2}_{VVV}) +i(H^{2}_{UUU}+H^{2}_{UVV}-H^{1}_{UUV}-H^{1}_{VVV})], \end{align*} $$

and it is said that the bifurcation is supercritical if $ a({1}/{\xi _1})<0$ and subcritical if $a({1}/{\xi _1})>0$ . As a result of the above discussion, we state the following theorem.

Theorem 4.4. If ${\xi _2}/{4}<\xi _1$ and $ a({1}/{\xi _1})\neq 0$ , then system (1.8) has a Neimark–Sacker bifurcation at the fixed point $E_{*}=(M_{*},N_{*})$ and this bifurcation is supercritical if $ a({1}/{\xi _1})<0$ and subcritical if $a({1}/{\xi _1})>0$ .

5. Examples

This section provides two numerical examples to strengthen and broaden the conclusions drawn from the theory in the other sections.

Example 5.1. This example is constructed by choosing the parameters in [7]. Here, we obtain the system

$$ \begin{align*} \begin{cases} M_{n+1}=M_{n}+ 2.11101\bigg(M_{n} (1-M_{n} )-\dfrac{M_{n} N_{n} }{26+39.2308 M_{n}}\bigg),\\ N_{n+1}=N_{n}+ 2.11101\bigg(\dfrac{8 M_{n} N_{n} }{26+39.2308 M_{n}}-0.12 N_{n}\bigg) \end{cases} \end{align*} $$

with the initial conditions

$$ \begin{align*} M_0=0.1 \quad \text{and} \quad N_0=0.1. \end{align*} $$

The equilibrium point for this system is determined to be $E_*=(0.947664,3.30649)$ . As the first condition in Theorem 3.3 holds, we conclude that this point is a sink as illustrated in Figure 1(a) and the phase portrait of the system is given in Figure 1(b).

Figure 1 (a) Stability of $E_*=(0.947664,3.30649)$ of system (1.8) and (b) phase portrait for $(M_{n},N_{n}).$

However, since the conditions of Theorem 4.2 are satisfied, we have established that there exists a flip bifurcation when $\rho $ exceeds the critical value $1.93915$ . This can also be observed in Figures 2(a) and 2(b), where the flip bifurcation takes place within the range $[0, 3],$ which includes the critical value $\rho _c=1.93915$ .

Figure 2 Bifurcation diagrams in $\rho \in [0,3]$ with initial condition $(M_{0},N_{0})=(0.1,0.1).$

Example 5.2. Let us consider the following system:

$$ \begin{align*} \begin{cases} M_{n+1}=M_{n}+ 1.5957 \bigg(M_{n} (1-M_{n} )-\dfrac{M_{n} N_{n} }{26+39.2308 M_{n}}\bigg),\\ N_{n+1}=N_{n}+ 1.5957\bigg(\dfrac{8 M_{n} N_{n} }{26+39.2308 M_{n}}-0.06 N_{n}\bigg) \end{cases} \end{align*} $$

with the initial conditions

$$ \begin{align*} M_0=0.1 \quad \text{and} \quad N_0=0.1. \end{align*} $$

The equilibrium $E_*=(0.276294,26.6608)$ of this system is a sink, as it satisfies the first condition of Theorem 3.3. This can be seen in Figure 3(a) and the phase portrait of the system is given in Figure 3(b). Meanwhile, since the conditions in Theorem 2.6 are true, there exists a Neimark–Sacker bifurcation at $\rho _{c}=1.83218$ . This bifurcation is supercritical in view of (4.11), which is calculated as $-0.000042416 < 0$ . This behaviour is presented in Figures 4(a) and 4(b), which occurs within the range $[0, 3]$ .

Figure 3 (a) Stability of $E_*=(0.276294,26.6608)$ of system (1.8) and (b) phase portrait for $(M_{n},N_{n}).$

Figure 4 Bifurcation diagrams in $\rho \in [0,3]$ with initial condition $(M_{0},N_{0})=(0.1,0.1).$

6. Conclusion

This paper explores the bifurcation properties of a fractional-order prey–predator model influenced by the wind, which is one of the most important abiotic factors often overlooked in ecological studies. The model was discretized using a piecewise discretization approach, allowing a detailed analysis of its discrete-time dynamics. We choose $\rho =h^\alpha $ as the bifurcation parameter that consists of both the discretization parameter h and the fractional order $\alpha $ of the system. The local stability of fixed points was investigated, and it was shown using the centre-manifold theorem and bifurcation theory that the system experiences both flip and Neimark–Sacker bifurcations at a positive fixed point. Numerical simulations were provided to illustrate and validate the theoretical results.

This work underscores the importance of including the abiotic factors such as wind in ecological models and demonstrates the advantages of fractional-order systems for capturing complex dynamics. Future research could focus on extending the model to include additional ecological interactions or exploring control strategies for practical applications. Further, validating the theoretical findings with empirical data could enhance the model’s relevance and applicability to the real-world ecological scenarios.