Dynamic Analysis of a Fractional-Modified Leslie-Gower Model with Predator Harvesting

1. Introduction

Predation is a crucial interaction that occurs between various species in nature, enabling the existence of most species in our ecosystem and supporting the rich biodiversity of complex ecosystems [1]. Mathematical models and their dynamics provide an apt explanation for the complexities arising from predator-prey relationships. These models also offer methods to optimally manage renewable resources and establish coexistence conditions for predators and prey [2]. The most fundamental model proposed for this ecosystem was the Lotka-Volterra model, introduced in the early 20th century [3,4]. While this model captured the oscillating behavior in populations of predators and their prey, its simplicity is unable to address most real-world scenarios, leading researchers to propose many modifications of the Lotka-Volterra type model [5].

An alternative to the Lotka-Volterra model is the Leslie model proposed by Leslie [6]. This model suggests that the number of prey and the carrying capacity of the predator’s environment are proportional, a factor not included in the Lotka-Volterra model. Leslie’s model addresses the issue of unbounded growth in both predator and prey populations, and it has been shown to possess globally asymptotically stable and unique positive equilibrium for any permissible parameters [7]. The model’s stability with Holling functional responses was later shown by May [8], and global stability of the unique interior equilibrium was proved for a Leslie-Gower predator-prey model with feedback controls by Chen et al. [9]. The uniqueness of limit cycles and Hopf bifurcation for this model were also discussed in [10,11].

While predator-prey models are widely used to understand complex interactions between different species in an ecosystem, the primary limitation of the Leslie-Gower predator-prey model is that the predator cannot switch to alternative food sources if the prey population is at low densities [12]. To address this limitation, Aziz-Alaoui and Daher [13] introduced a provisional alternative food source parameter, leading to the modified Leslie-Gower predator-prey model. This model has been utilized in various applications, such as modeling the biological control of the prickly-pear cactus by the moth Cactoblastis cactorum [14] and predator-prey mite outbreak interactions on fruit trees [15,16]. For further details on the modified Leslie-Gower model and its applications, see [17,18] and the references therein. Additionally, the Allee effect, introduced by W.C. Allee in [19], has been incorporated into the Leslie-Gower type models, and the stability dynamics and bifurcation analysis of predator-prey systems subject to Allee effects are discussed in [20-23].

Harvesting biological species is a necessary practice for economic development and resource utilization. However, human population growth and natural resource exploitation have resulted in ecological imbalances and species extinction. Although harvesting can minimize damage, it can also lead to predator species extinction [24]. Therefore, it is essential to reinforce scientific management of harvesting practices. Clark [2] discussed the problem of non-selective harvesting of two ecologically independent populations with the logistic growth law. Multi-species harvesting models have been studied in detail by Chaudhuri [25,26], Mesterton-Gibbons [27], and Kar and Chaudhuri [28,29]. Non-selective harvesting models of prey-predator fisheries have also been studied in detail by Chaudhuri and Ray [30] and Kar et al., [31]. Scientific management of harvesting practices can help minimize the impact on predator-prey dynamics and promote ecological stability.

The Michaelis-Menten-type predator-prey model proposes that the per capita growth rate of predators is directly affected by the ratio of prey to predator abundance. Experimental and observational data support this model’s effectiveness, particularly for predators that must compete for prey. For further details, see [32]-[33]. This model has been extensively studied ([34]-[35]), revealing rich dynamics such as stable limit cycles, multiple attractors, and deterministic extinction. Various parametric values result in the existence of hyperbolic, parabolic, and elliptic orbits near the origin and combinations of such orbits.

Continuous predator-prey models have been extensively used to study the interactions between two species. However, these models have some limitations. Specifically, they assume that the subject species have continuous and overlapping generations, which is not always the case. For instance, salmon have an annual spawning season and are born at the same time each year. In such cases, discrete-time models are more appropriate since they provide a better description of the population dynamics. In fact, discrete models often yield richer and more realistic results than continuous models. Moreover, since many continuous models cannot be solved analytically, using difference equations for approximation and finding solutions is a practical approach. Discrete-time models are widely used in population biology and complex ecosystems to examine taxonomic groups of organisms and species with the passage of time. These models are also well-suited for describing the chaotic behavior of nonlinear dynamics.

Fractional calculus is a branch of mathematics that studies the generalization of traditional calculus to include non-integer order derivatives and integrals. The origins of fractional calculus can be traced back to the seventeenth century, with the work of mathematicians such as Leibniz and Euler. However, it was not until the nineteenth century that the concept of fractional differentiation and integration was introduced by Liouville, who investigated the fractional calculus of analytic functions. In the twentieth century, fractional calculus gained further attention as it proved to be useful in a wide range of fields, including physics, engineering, and biology. One of the most significant developments in the field occurred in the 1960s and 1970s, when mathematicians and scientists realized that fractional calculus could be used to model many physical systems that exhibit anomalous diffusion, such as transport in porous media and electrical conduction in disordered systems. Today, fractional calculus continues to be an active area of research, with applications in various fields such as signal processing, control theory, and finance. The study of fractional calculus involves understanding the properties and behaviors of fractional derivatives and integrals, as well as developing techniques for their computation and application. The history of fractional calculus highlights the importance of exploring new and unconventional mathematical concepts, as they can often lead to significant breakthroughs in science and technology. Fractional derivatives and integrals are powerful tools for modeling non-locality and memory effects, which are not adequately captured by integer-order calculus. In contrast to integer-order calculus, fractional derivatives and integrals possess a non-local character and represent a suitable approach for modeling phenomena that depend on memory.

In this paper, we explore the fractional modified Leslie-Gower predator-prey model with Michaelis-Menten type prey harvesting. We utilize the Caputo fractional-order derivative operator due to its numerous advantages as previously discussed. The updated system, based on the model proposed in [36], is presented below: $$\begin{array}{rlr}\frac{dx}{dt}& =& x[t](1-x[t]–\frac{ay[t]}{m+x[t]}-\frac{k}{c+x[t]}),\\ \frac{dy}{dt}& =& ry[t](1-\frac{by[t]}{m+x[t]}).\end{array}$$

In this model, the temporal dynamics with memory can only be captured by the fractional derivative. The fractional integral of order $\beta \in R^{+}$ of the function $f(t),t>0,$ is defined as $$I^{\beta }f(t)={\int }_0{}^t\frac{(t-s{)}^{\beta -1}}{\mathrm{\Gamma }(\beta )}f(s)ds,$$ and the fractional derivative of order $\alpha \in (n-1,n)$ of $f(t),t>0,$ is defined as $${}^C{D}_t^{\alpha }f(t)=I^{n-\alpha }{D}_t{}^{n}f(t),$$ where $D_t=\frac{d}{dt}.$ Therefore,

$${}^C{D}_t^{\alpha }f(t)=\frac{1}{\mathrm{\Gamma }(n-\alpha )}{\int }_0{}^t(t-s{)}^{n-\alpha -1}{f}^{(n)}(s)ds.$$

Using the caputo derivative, instead of the classical one, we get, $$\begin{array}{}\text{(1)}& \{\begin{array}{l}{}^C{D}_t^{\theta }x(t)=x(1-x–\frac{ay}{m+x}-\frac{k}{c+x}),\\ {}^C{D}_t^{\theta }y(t)=ry(1-\frac{by}{m+x}),\end{array}\end{array}$$ where, ${}^C{D}_t^{\theta }$ is the fractional order derivative operator of Caputo type of order $\theta$.

The investigation presented in this paper focuses on a modified Leslie-Gower predator-prey model with fractional discrete-time Michaelis-Menten type prey harvesting. The analysis reveals the occurrence of nonnegative interior fixed points and their stability dynamics. Using the center manifold theorem and bifurcation theory, we derive conditions for flip and Neimark-Sacker bifurcations. The theoretical findings are consistently demonstrated through numerical simulations conducted with a computer package. Overall, our study highlights the importance of incorporating fractional calculus in predator-prey systems with harvesting and sheds light on the complex dynamics that arise in such models.

2. Existence and uniqueness of the solution for the fractional system (1)

We will now establish the existence and uniqueness of the solution for the fractional system discussed earlier. To do this, we consider the initial value problem (IVP): $${}^C{D}_t^{\theta }x(t)=f_x(t,x),,$$ with initial condition $x(0)=x_0$.

By utilizing the Volterra-type integral equation for the aforementioned initial value problem (IVP), we derive: $$x(t)=x_0+\frac{1}{\mathrm{\Gamma }(\theta )}{\int }_0{}^t(t-\mu {)}^{\theta -1}{f}_x(\mu ,x)d\mu ,.$$

To ensure the existence of a solution for this IVP, we refer to a lemma proposed by Katugampola [6]:

2.1 Existence

We will now examine the existence of a solution for the IVP based on the theorem presented by [36,37]:

Note that $P$ is a compact set since it is closed and bounded in ${\mathbb{R}}^2$. The function $f_x$ is continuous on $P$ and $M_x$ is well-defined since $P$ is compact. Hence, the conditions of the theorem are satisfied and a solution $x$ exists for the IVP.

2.2 Uniqueness

Katugampola [6] discusses the uniqueness of solutions for the fractional system, and we state the following theorem:

2.3 Process of Discretization

Let us discuss the process of discretization, which has been previously introduced in the literature [5,6, 37,38], for differential equations involving fractional orders. In this study, we will employ this technique for modeling the fractional glycolysis. We will discretize the system using a method that involves arguments with piecewise constants.

To begin, we consider $t$ satisfying the condition $0\le t<h$, which implies that $0\le \frac{t}{h}<1$. This allows us to express the discretization of the system with the following equations:

$$\begin{array}{rlr}{}^C{D}_t^{\theta }x(t)& =& x\left(\left[\frac{t}{h}\right]h\right)(1-x\left(\left[\frac{t}{h}\right]h\right)–\frac{ay\left(\left[\frac{t}{h}\right]h\right)}{mx\left(\left[\frac{t}{h}\right]h\right)}-\frac{k}{c+x\left(\left[\frac{t}{h}\right]h\right)}),\\ {}^C{D}_t^{\theta }y(t)& =& ry\left(\left[\frac{t}{h}\right]h\right)(1-\frac{by\left(\left[\frac{t}{h}\right]h\right)}{m+x\left(\left[\frac{t}{h}\right]h\right)}),\end{array}$$ where $[\cdot ]$ denotes the integer part function. Using the above equations and the condition on $t$, we can simplify the solution to obtain:

$$\begin{array}{rlr}{x}_1& =& x_0+\frac{h^{\theta }}{\mathrm{\Gamma }(1+\theta )}{x}_0(1-x_0–\frac{a{y}_0}{m+x_0}-\frac{k}{c+x_0}),\\ y_1& =& y_0+\frac{h^{\theta }}{\mathrm{\Gamma }(1+\theta )}r{y}_0(1-\frac{b{y}_0}{m+x_0}).\end{array}$$

Next, we consider $h\le t<2h$, which means $1\le \frac{t}{h}<2$. Using this, we get:

$$\begin{array}{rlr}{x}_2& =& x_1+\frac{(t-h{)}^{\theta }}{\mathrm{\Gamma }(1+\theta )}{x}_1(1-x_1–\frac{a{y}_1}{m+x_1}-\frac{k}{c+x_1}),\\ y_2& =& y_1+\frac{h^{\theta }}{\mathrm{\Gamma }(1+\theta )}r{y}_1(1-\frac{b{y}_1}{m+x_1}).\end{array}$$

Repeating this process $n$ times yields: $$\begin{array}{rlr}{x}_{n+1}& =& x_n+\frac{(t-nh{)}^{\theta }}{\mathrm{\Gamma }(1+\theta )}{x}_n(1-x_n–\frac{a{y}_n}{m+x_n}-\frac{k}{c+x_n}),\\ y_{n+1}& =& y_n+\frac{(t-nh{)}^{\theta }}{\mathrm{\Gamma }(1+\theta )}r{y}_n(1-\frac{b{y}_n}{m+x_n}),\end{array}$$ where, $nh\le t<(n+1)h$. For $t\to (n+1)h$, system is reduced to $$\begin{array}{}\text{(3)}& \{\begin{array}{l}{x}_{n+1}=x_n+\frac{h^{\theta }}{\mathrm{\Gamma }(1+\theta )}{x}_n(1-x_n–\frac{a{y}_n}{m+x_n}-\frac{k}{c+x_n}),\\ y_{n+1}=y_n+\frac{h^{\theta }}{\mathrm{\Gamma }(1+\theta )}r{y}_n(1-\frac{b{y}_n}{m+x_n}).\end{array}\end{array}$$

2.4 Stability of Fixed Points

In order to find the fixed points, we solve the following system. $$\begin{array}{rlr}x& =& x+\frac{h^{\theta }}{\mathrm{\Gamma }(1+\theta )}x(1-x–\frac{ay}{m+x}-\frac{k}{c+x}),\\ y& =& y+\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}ry(1-\frac{by}{m+x}).\end{array}$$ The fixed points are

$$\begin{array}{rl}& \{E_0,E_1,E_2,E_3,E_m,E_p\}\\ & =\{\{0,0\},\{\frac{1}{2}(1-c-\sqrt{(1+c{)}^2-4h}),0\},\{\frac{1}{2}(1-c+\sqrt{(1+c{)}^2-4h}),0\},\{0,\frac{m}{\beta }\},\{x_m,y_m\},\{x_p,y_p\}\},\end{array}$$

where $$\begin{array}{rlr}{x}_m& =& \frac{1-\frac{a}{b}-c-\sqrt{{(\frac{a}{b}-1-c)}^2-4k}}{2},\\ x_p& =& \frac{1-\frac{a}{b}-c+\sqrt{{(\frac{a}{b}-1-c)}^2-4k}}{2},\\ y_m& =& \frac{m+x_m}{b},\\ y_p& =& \frac{m+x_p}{b}.\end{array}$$ The Jacobian matrix of the system is: $$J=\left(\begin{array}{cc}1+\frac{h^{\theta }x(-1+\frac{k}{(c+x{)}^2}+\frac{ay}{(m+x{)}^2})}{\mathrm{\Gamma }[1+\theta ]}+\frac{h^{\theta }(1-x-\frac{k}{c+x}-\frac{ay}{m+x})}{\mathrm{\Gamma }[1+\theta ]}& -\frac{a{h}^{\theta }x}{(m+x)\mathrm{\Gamma }[1+\theta ]}\\ \frac{b{h}^{\theta }r{y}^2}{(m+x{)}^2\mathrm{\Gamma }[1+\theta ]}& 1-\frac{b{h}^{\theta }ry}{(m+x)\mathrm{\Gamma }[1+\theta ]}+\frac{h^{\theta }r(1-\frac{by}{m+x})}{\mathrm{\Gamma }[1+\theta ]}\end{array}\right).$$ To determine the properties and topology of the fixed points mentioned earlier, we will utilize the following lemma, which can be found in textbooks on discrete dynamical systems.

It should be noted that if ${\lambda }_{1,2}$ are the eigenvalues of the 2 by 2 Jacobian matrix, the following lemma holds:

To determine the stability of the fixed points, we will evaluate $J$ at each of the boundary fixed points $E_0$, $E_1$, $E_2$, and $E_3$, as well as the positive interior fixed points $E_m$ and $E_p$. $$J_0=J{|}_{(0,0)}=\left(\begin{array}{cc}1+\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}(1-\frac{k}{c})& 0\\ 0& 1+\frac{h^{\theta }r}{\mathrm{\Gamma }[1+\theta ]}\end{array}\right)$$ and the corresponding eigenvalues are $\{1+\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}r,1+\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\left(\frac{c-k}{c}\right)\}.$ Thus the fixed point is saddle if $k>c,$ unstable if $k<c$ and non- hyperbolic if $k=c.$ $$J_{1,2}=J{|}_{E_1,E_2}=\left(\begin{array}{cc}1+\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}{x}_{1,2}(-1+\frac{k}{(c+x_{1,2}){}^2})& -\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{x_{1,2}a}{m+x_{1,2}}\\ 0& 1+\frac{h^{\theta }r}{\mathrm{\Gamma }[1+\theta ]}\end{array}\right)$$ and the corresponding eigenvalues are $\{1+\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}r,1+\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}{x}_{1,2}(\frac{k}{(c+x_{1,2}){}^2}-1)\}.$ Thus, for both boundary fixed points $E_1$ and $E_2$, if $k<{(c+x_{1,2})}^2$, i-e, if $0<c<1\mathrm{\&}\mathrm{\&}c<k\le \frac{(1+c{)}^2}{4}$, then the fixed point is saddle. If $k>{(c+x_{1,2})}^2,$ i-e, $(0<c\le 1\mathrm{\&}\mathrm{\&}0<k<c)\parallel (c>1\mathrm{\&}\mathrm{\&}0<k\le \frac{(1+c{)}^2}{4})$, then the fixed point is unstable. If $k={(c+x_{1,2})}^2,$ the the fixed point may undergo transcritical or fold bifurcation and if $k=(1-\frac{2\mathrm{\Gamma }[1+\theta ]}{h^{\theta }{x}_{1,2}}){(c+x_{1,2})}^2,$ then the fixed point undergoes period- douling bifurcation. The possibility of Neimark-Sacker bifurcation is mute for these fixed points. $$J_3=J{|}_{(0,\frac{m}{\beta })}=\left(\begin{array}{cc}1+\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}(1-\frac{k}{c}-\frac{a}{b})& 0\\ \frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{r}{b}& 1-\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}r\end{array}\right)$$ and the corresponding eigenvalues are $\{1+\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}(1-\frac{k}{c}-\frac{a}{b}),1-\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}r\}$. Therefore if $k<c(1-\frac{a}{b})$ the the fixed point is saddle, if $k>c(1-\frac{a}{b})$ the fixed point is stable, if $k=c(1-\frac{a}{b})$ the fixed point may undergo transcritical or fold bifurcation and if $k=c(1-\frac{a}{b}+\frac{2\mathrm{\Gamma }[1+\theta ]}{h^{\theta }}),$ the fixed point undergoes period-douling bifurcation. The possibility of Neimark-Sacker bifurcation is mute for this fixed point also.

2.5 Interior Fixed Points

For the interior fixed points the analysis of stability is the same and we will perform it under the generic variable name $x_{m,p}=x$ and $y_{m,p}=y$. Let the Jacobian matrix of the system evaluated at $E_{m,p}$ be $$\overline{J}=\left(\begin{array}{cc}1+\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}x(-1+\frac{k}{(c+x{)}^2}+\frac{ay}{(m+x{)}^2})& -\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{ax}{(m+x)}\\ \frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{ry}{(m+x)}& 1-\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}r\end{array}\right).$$

3. Bifurcations

Bifurcation is a fundamental concept in dynamical systems theory that describes the qualitative changes that occur in a system as a parameter is varied. The basic idea is that the behavior of the system can change abruptly as the parameter crosses a certain threshold value. There are several types of bifurcations that can occur, including saddle-node bifurcation, pitchfork bifurcation, transcritical bifurcation, Hopf bifurcation, and period-doubling bifurcation. Each type of bifurcation has its own characteristic behavior and can lead to different types of dynamical behavior, such as limit cycles, chaos, and multistability. Understanding bifurcations is crucial for analyzing the behavior of complex systems and predicting how they will respond to changes in their environment or parameters [39,40].

3.1 Period – Doubling Bifurcation

Let $M_{PD}=\{(a,b,c,h,k,m,r,\theta )\in {\mathbb{R}}_{+}^7:h\in \{h_1,h_2\},h\notin \{{\left(\frac{2\mathrm{\Gamma }(1+\theta )}{r+x-\frac{axy}{(m+x{)}^2}}\right)}^{\frac{1}{\theta }},{\left(\frac{4\mathrm{\Gamma }(1+\theta )}{r+x-\frac{axy}{(m+x{)}^2}}\right)}^{\frac{1}{\theta }}\}\}$. Let $h$ be the bifurcation parameter. Then, variation of parameters in small neighborhood of $M_{PD}$ gives emergence of period doubling bifurcation. Let $k\in M_{PD}$ and let $H$ be the perturbation of $h$ where $|H|<<<1$. The system can be written as $$\left(\begin{array}{c}x\\ y\end{array}\right)\to \left(\begin{array}{c}x+\frac{(h+H{)}^{\theta }}{\mathrm{\Gamma }[1+\theta ]}x(1-x–\frac{ay}{m+x}-\frac{k}{c+x})\\ y+\frac{(h+H{)}^{\theta }}{\mathrm{\Gamma }[1+\theta ]}ry(1-\frac{by}{m+x})\end{array}\right).$$ Translate the fixed point to $(0,0).$ Define $X=x-\overline{x}$ and $Y=y-\overline{y}$. Then, the above system can be written as $$\left(\begin{array}{c}X\\ Y\end{array}\right)\to \left(\begin{array}{c}X+\frac{(h+H{)}^{\theta }(X+\overline{x})(1-X-\overline{x}-\frac{k}{c+X+\overline{x}}-\frac{a(Y+\overline{y})}{m+X+\overline{x}})}{\mathrm{\Gamma }[1+\theta ]}\\ Y+\frac{(h+H{)}^{\theta }r(Y+\overline{y})(1-\frac{b(Y+\overline{y})}{m+X+\overline{x}})}{\mathrm{\Gamma }[1+\theta ]}\end{array}\right).$$ If we define,
$a_{11}=1+\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\overline{x}(-1+\frac{k}{{(c+\overline{x})}^2}+\frac{a\overline{y}}{{(m+\overline{x})}^2});$
$a_{12}=-\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{a\overline{x}}{(m+\overline{x})};$
$b_1=\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}(-1+\frac{k}{{(c+\overline{x})}^2}+\frac{a\overline{y}}{{(m+\overline{x})}^2}+\overline{x}(-\frac{k}{{(c+\overline{x})}^3}-\frac{a\overline{y}}{{(m+\overline{x})}^3}));$
$b_2=-\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{am}{{(m+\overline{x})}^2};$
$b_3=\frac{h^{-2+\theta }(-1+\theta )\theta \left(\overline{x}(-1+\frac{k}{{(c+\overline{x})}^2}+\frac{a\overline{y}}{{(m+\overline{x})}^2})\right)}{2\mathrm{\Gamma }[1+\theta ]}$;
$b_4=-\frac{a{h}^{-1+\theta }\theta \overline{x}}{(m+\overline{x})\mathrm{\Gamma }[1+\theta ]}$;
$b_5=0;$
$b_6=\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}(-\frac{k}{{(c+\overline{x})}^3}-\frac{a\overline{y}}{{(m+\overline{x})}^3}+\overline{x}(\frac{k}{{(c+\overline{x})}^4}+\frac{a\overline{y}}{{(m+\overline{x})}^4}));$
$b_7=\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{am}{{(m+\overline{x})}^3};$
$b_8=\frac{h^{-1+\theta }\theta (-1+\frac{k}{{(c+\overline{x})}^2}+\frac{a\overline{y}}{{(m+\overline{x})}^2}+\overline{x}(-\frac{k}{{(c+\overline{x})}^3}-\frac{a\overline{y}}{{(m+\overline{x})}^3}))}{\mathrm{\Gamma }[1+\theta ]};$
$b_9=-\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{a{h}^{-1}m\theta }{{(m+\overline{x})}^2};$
$b_{10}=\frac{h^{-2+\theta }(-1+\theta )\theta \left(\overline{x}(-1+\frac{k}{{(c+\overline{x})}^2}+\frac{a\overline{y}}{{(m+\overline{x})}^2})\right)}{2\mathrm{\Gamma }[1+\theta ]};$
$b_{11}=-\frac{a{h}^{-2+\theta }(-1+\theta )\theta \overline{x}}{2(m+\overline{x})\mathrm{\Gamma }[1+\theta ]};$
$b_{12}=0;$
$b_{13}=0;$
$b_{14}=0;$
Then, $$\begin{array}{rl}{f}_{PD}^1(X,Y,H)=& b_1{X}^2+b_{2}XY+b_{3}HX+b_{4}HY+b_5{Y}^2+b_6{X}^3+b_7{X}^{2}Y+b_{8}H{X}^2+b_{9}HXY\\ & +b_{10}{H}^{2}X+b_{11}{H}^{2}Y+b_{12}X{Y}^2+b_{13}H{Y}^2+b_{14}{Y}^3+O\left(|H+X+Y{|}^4\right),\end{array}$$ and $$X\to a_{11}X+a_{12}Y+f_{PD}^1(X,Y,H).$$ Similarly, we can write $Y\to a_{21}X+a_{22}Y+f_{PD}^2(X,Y,H)$ by defining
$a_{21}=\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{br{\overline{y}}^2}{{(m+\overline{x})}^2}$;
$a_{22}=1-\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}r;$
$b_{15}=-\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{br{\overline{y}}^2}{{(m+\overline{x})}^3};$
$b_{16}=\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{2br\overline{y}}{{(m+\overline{x})}^2};$
$b_{17}=\frac{h^{-1+\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{br\theta {\overline{y}}^2}{{(m+\overline{x})}^2};$
$b_{18}=-\frac{h^{-1+\theta }}{\mathrm{\Gamma }[1+\theta ]}r\theta ;$
$b_{19}=-\frac{b{h}^{\theta }r}{(m+\overline{x})\mathrm{\Gamma }[1+\theta ]};$
$b_{20}=\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{br{\overline{y}}^2}{{(m+\overline{x})}^4};$
$b_{21}=-\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{2br\overline{y}}{{(m+\overline{x})}^3};$
$b_{22}=-\frac{h^{-1+\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{br\theta {\overline{y}}^2}{{(m+\overline{x})}^3}$;
$b_{23}=\frac{h^{-1+\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{2br\theta \overline{y}}{{(m+\overline{x})}^2};$
$b_{24}=\frac{h^{-2+\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{br(-1+\theta )\theta {\overline{y}}^2}{2{(m+\overline{x})}^2};$
$b_{25}=-\frac{h^{-2+\theta }}{2\mathrm{\Gamma }[1+\theta ]}r(-1+\theta )\theta ;$
$b_{26}=\frac{h^{\theta }}{\mathrm{\Gamma }[1+\theta ]}\frac{br}{{(m+\overline{x})}^2};$
$b_{27}=-\frac{b{h}^{-1+\theta }r\theta }{(m+\overline{x})\mathrm{\Gamma }[1+\theta ]}$;
$b_{28}=0;$
and
$$\begin{array}{rl}{f}_{PD}^2(X,Y,H)=& b_{15}{X}^2+b_{16}XY+b_{17}HX+b_{18}HY+b_{19}{Y}^2+b_{20}{X}^3+b_{21}{X}^{2}Y+b_{22}H{X}^2+b_{23}HXY\\ & +b_{24}{H}^{2}X+b_{25}{H}^{2}Y+b_{26}X{Y}^2+b_{27}H{Y}^2+b_{28}{Y}^3+O\left(|H+X+Y{|}^4\right).\end{array}$$ Thus, the system can be written as $$\left(\begin{array}{c}X\\ Y\end{array}\right)\to \left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right)\left(\begin{array}{c}X\\ Y\end{array}\right)+\left(\begin{array}{c}{f}_{PD}^1(X,Y,H)\\ f_{PD}^2(X,Y,H)\end{array}\right).$$ Let $T=\left(\begin{array}{cc}{a}_{12}& a_{12}\\ -1-a_{11}& {\lambda }_2–a_{22}\end{array}\right)$ be an invertable transformation matrix such that $\left(\begin{array}{c}X\\ Y\end{array}\right)\to T\left(\begin{array}{c}\mu \\ \eta \end{array}\right)$
The transformed system is given by $$\left(\begin{array}{c}\mu \\ \eta \end{array}\right)\to \left(\begin{array}{cc}-1& 0\\ 0& {\lambda }_2\end{array}\right)\left(\begin{array}{c}\mu \\ \eta \end{array}\right)+\left(\begin{array}{c}{g}_{PD}^1(X,Y,H)\\ g_{PD}^2(X,Y,H)\end{array}\right)$$ where $$\begin{array}{rl}{g}_{PD}^1(\mu ,\eta ,H)=& d_1{\mu }^2+d_2\mu \eta +d_{3}H\mu +d_{4}H\eta +d_5{H}^2+d_6{\eta }^2+d_7{\mu }^3+d_8{\mu }^2\eta +d_{9}H{\mu }^2+d_{10}H\mu \eta \\ & +d_{11}{H}^2\mu +d_{12}{H}^2\eta +d_{13}\mu {\eta }^2+d_{14}H{\eta }^2+d_{15}{\eta }^3+O\left(|H+\mu +\eta {|}^4\right)\\ g_{PD}^2(\mu ,\eta ,H)=& d_{16}{\mu }^2+d_{17}\mu \eta +d_{18}H\mu +d_{19}H\eta +d_{20}{H}^2+d_{21}{\eta }^2+d_{22}{\mu }^3+d_{23}{\mu }^2\eta +d_{24}H{\mu }^2\\ & +d_{25}H\mu \eta +d_{26}{H}^2\mu +d_{27}{H}^2\eta +d_{28}\mu {\eta }^2+d_{29}H{\eta }^2+d_{30}{\eta }^3+O\left(|H+\mu +\eta {|}^4\right)\end{array}$$ and $$\begin{array}{rl}{d}_1=& a_{12}^2{c}_1(-b_{15}+b_1{c}_2)+(-1-a_{11})a_{12}{c}_1(-b_{16}+b_2{c}_2)+(-1-a_{11}){}^2{c}_1(-b_{19}+b_5{c}_2);\\ d_2=& 2a_{12}^2{c}_1(-b_{15}+b_1{c}_2)+(-1-a_{11})a_{12}{c}_1(-b_{16}+b_2{c}_2)+a_{12}{c}_1(-b_{16}+b_2{c}_2)(-a_{22}+{\lambda }_2)\\ & +2(-1-a_{11})c_1(-b_{19}+b_5{c}_2)(-a_{22}+{\lambda }_2);\\ d_3=& a_{12}{c}_1(-b_{17}+b_3{c}_2)+(-1-a_{11})c_1(-b_{18}+b_4{c}_2);\\ d_4=& a_{12}{c}_1(-b_{17}+b_3{c}_2)+c_1(-b_{18}+b_4{c}_2)(-a_{22}+{\lambda }_2);\\ d_5=& 0;\\ d_6=& a_{12}^2{c}_1(-b_{15}+b_1{c}_2)+a_{12}{c}_1(-b_{16}+b_2{c}_2)(-a_{22}+{\lambda }_2)+c_1(-b_{19}+b_5{c}_2){(-a_{22}+{\lambda }_2)}^2;\\ d_7=& a_{12}^3{c}_1(-b_{20}+b_6{c}_2)+(-1-a_{11})a_{12}^2{c}_1(-b_{21}+b_7{c}_2)+(-1-a_{11}){}^2{a}_{12}{c}_1(-b_{26}+b_{12}{c}_2)\\ & +(-1-a_{11}){}^3{c}_1(-b_{28}+b_{14}{c}_2);\\ d_8=& 3a_{12}^3{c}_1(-b_{20}+b_6{c}_2)+2(-1-a_{11})a_{12}^2{c}_1(-b_{21}+b_7{c}_2)+{(-1-a_{11})}^2{a}_{12}{c}_1(-b_{26}+b_{12}{c}_2)\\ & +a_{12}^2{c}_1(-b_{21}+b_7{c}_2)(-a_{22}+{\lambda }_2)+2(-1-a_{11})a_{12}{c}_1(-b_{26}+b_{12}{c}_2)(-a_{22}+{\lambda }_2)\\ & +3{(-1-a_{11})}^2{c}_1(-b_{28}+b_{14}{c}_2)(-a_{22}+{\lambda }_2);\\ d_9=& a_{12}^2{c}_1(-b_{22}+b_8{c}_2)+(-1-a_{11})a_{12}{c}_1(-b_{23}+b_9{c}_2)+(-1-a_{11}){}^2{c}_1(-b_{27}+b_{13}{c}_2);\\ d_{10}=& 2a_{12}^2{c}_1(-b_{22}+b_8{c}_2)+(-1-a_{11})a_{12}{c}_1(-b_{23}+b_9{c}_2)+a_{12}{c}_1(-b_{23}+b_9{c}_2)(-a_{22}+{\lambda }_2)\\ & +2(-1-a_{11})c_1(-b_{27}+b_{13}{c}_2)(-a_{22}+{\lambda }_2);\\ d_{11}=& a_{12}{c}_1(-b_{24}+b_{10}{c}_2)+(-1-a_{11})c_1(-b_{25}+b_{11}{c}_2);\\ d_{12}=& a_{12}{c}_1(-b_{24}+b_{10}{c}_2)+c_1(-b_{25}+b_{11}{c}_2)(-a_{22}+{\lambda }_2);\\ d_{13}=& 3a_{12}^3{c}_1(-b_{20}+b_6{c}_2)+(-1-a_{11})a_{12}^2{c}_1(-b_{21}+b_7{c}_2)+2a_{12}^2{c}_1(-b_{21}+b_7{c}_2)(-a_{22}+{\lambda }_2)\\ & +2(-1-a_{11})a_{12}{c}_1(-b_{26}+b_{12}{c}_2)(-a_{22}+{\lambda }_2)+a_{12}{c}_1(-b_{26}+b_{12}{c}_2){(-a_{22}+{\lambda }_2)}^2\\ & +3(-1-a_{11})c_1(-b_{28}+b_{14}{c}_2){(-a_{22}+{\lambda }_2)}^2;\\ d_{14}=& a_{12}^2{c}_1(-b_{22}+b_8{c}_2)+a_{12}{c}_1(-b_{23}+b_9{c}_2)(-a_{22}+{\lambda }_2)+c_1(-b_{27}+b_{13}{c}_2)(-a_{22}+{\lambda }_2){}^2;\\ d_{15}=& a_{12}^3{c}_1(-b_{20}+b_6{c}_2)+a_{12}^2{c}_1(-b_{21}+b_7{c}_2)(-a_{22}+{\lambda }_2)+a_{12}{c}_1(-b_{26}+b_{12}{c}_2)(-a_{22}+{\lambda }_2){}^2\\ & +c_1(-b_{28}+b_{14}{c}_2)(-a_{22}+{\lambda }_2){}^3;\\ d_{16}=& a_{12}^2{c}_1(-b_{15}+b_1{c}_3)+(-1-a_{11})a_{12}{c}_1(-b_{16}+b_2{c}_3)+(-1-a_{11}){}^2{c}_1(-b_{19}+b_5{c}_3);\\ d_{17}=& 2a_{12}^2{c}_1(-b_{15}+b_1{c}_3)+(-1-a_{11})a_{12}{c}_1(-b_{16}+b_2{c}_3)+a_{12}{c}_1(-b_{16}+b_2{c}_3)(-a_{22}+{\lambda }_2)\\ & +2(-1-a_{11})c_1(-b_{19}+b_5{c}_3)(-a_{22}+{\lambda }_2);\\ d_{18}=& a_{12}{c}_1(-b_{17}+b_3{c}_3)+(-1-a_{11})c_1(-b_{18}+b_4{c}_3);\\ d_{19}=& a_{12}{c}_1(-b_{17}+b_3{c}_3)+c_1(-b_{18}+b_4{c}_3)(-a_{22}+{\lambda }_2);\\ d_{20}=& 0;\\ d_{21}=& a_{12}^2{c}_1(-b_{15}+b_1{c}_3)+a_{12}{c}_1(-b_{16}+b_2{c}_3)(-a_{22}+{\lambda }_2)+c_1(-b_{19}+b_5{c}_3)(-a_{22}+{\lambda }_2){}^2;\\ d_{22}=& a_{12}^3{c}_1(-b_{20}+b_6{c}_3)+(-1-a_{11})a_{12}^2{c}_1(-b_{21}+b_7{c}_3)\\ & +(-1-a_{11}){}^2{a}_{12}{c}_1(-b_{26}+b_{12}{c}_3)+(-1-a_{11}){}^3{c}_1(-b_{28}+b_{14}{c}_3);\\ d_{23}=& 3a_{12}^3{c}_1(-b_{20}+b_6{c}_3)+2(-1-a_{11})a_{12}^2{c}_1(-b_{21}+b_7{c}_3)+(-1-a_{11}){}^2{a}_{12}{c}_1(-b_{26}+b_{12}{c}_3)\\ & +a_{12}^2{c}_1(-b_{21}+b_7{c}_3)(-a_{22}+{\lambda }_2)+2(-1-a_{11})a_{12}{c}_1(-b_{26}+b_{12}{c}_3)(-a_{22}+{\lambda }_2)\\ & +3(-1-a_{11}){}^2{c}_1(-b_{28}+b_{14}{c}_3)(-a_{22}+{\lambda }_2);\\ d_{24}=& a_{12}^2{c}_1(-b_{22}+b_8{c}_3)+(-1-a_{11})a_{12}{c}_1(-b_{23}+b_9{c}_3)+(-1-a_{11}){}^2{c}_1(-b_{27}+b_{13}{c}_3);\\ d_{25}=& 2a_{12}^2{c}_1(-b_{22}+b_8{c}_3)+(-1-a_{11})a_{12}{c}_1(-b_{23}+b_9{c}_3)+a_{12}{c}_1(-b_{23}+b_9{c}_3)(-a_{22}+{\lambda }_2)\\ & +2(-1-a_{11})c_1(-b_{27}+b_{13}{c}_3)(-a_{22}+{\lambda }_2);\\ d_{26}=& a_{12}{c}_1(-b_{24}+b_{10}{c}_3)+(-1-a_{11})c_1(-b_{25}+b_{11}{c}_3);\\ d_{27}=& a_{12}{c}_1(-b_{24}+b_{10}{c}_3)+c_1(-b_{25}+b_{11}{c}_3)(-a_{22}+{\lambda }_2);\\ d_{28}=& 3a_{12}^3{c}_1(-b_{20}+b_6{c}_3)+(-1-a_{11})a_{12}^2{c}_1(-b_{21}+b_7{c}_3)+2a_{12}^2{c}_1(-b_{21}+b_7{c}_3)(-a_{22}+{\lambda }_2)\\ & +2(-1-a_{11})a_{12}{c}_1(-b_{26}+b_{12}{c}_3)(-a_{22}+{\lambda }_2)+a_{12}{c}_1(-b_{26}+b_{12}{c}_3){(-a_{22}+{\lambda }_2)}^2\\ & +3(-1-a_{11})c_1(-b_{28}+b_{14}{c}_3)(-a_{22}+{\lambda }_2){}^2;\\ d_{29}=& a_{12}^2{c}_1(-b_{22}+b_8{c}_3)+a_{12}{c}_1(-b_{23}+b_9{c}_3)(-a_{22}+{\lambda }_2)+c_1(-b_{27}+b_{13}{c}_3)(-a_{22}+{\lambda }_2){}^2;\\ d_{30}=& a_{12}^3{c}_1(-b_{20}+b_6{c}_3)+a_{12}^2{c}_1(-b_{21}+b_7{c}_3)(-a_{22}+{\lambda }_2)+a_{12}{c}_1(-b_{26}+b_{12}{c}_3)(-a_{22}+{\lambda }_2){}^2\\ & +c_1(-b_{28}+b_{14}{c}_3)(-a_{22}+{\lambda }_2){}^3,\end{array}$$ and $$X=a_{12}(\mu +\eta ),Y=-(1+a_{11})\mu +({\lambda }_2–a_{22})\eta .$$ Define $M_C=\{(\mu ,\eta ,H)|\eta =\mathcal{G}(\mu ,H)\},$ $|\mu |<{\delta }_1,$ $|H|<{\delta }_2,$ $f_c(0,0)=0,$ $D{f}_c(0,0)=0$ and $$\mathcal{F}(\mathcal{G}(\mu ,H),H)=\mathcal{G}(-\mu +g_{PD}^1(\mu ,\mathcal{G}(\mu ,H),H),H)–{\lambda }_2\mathcal{G}(\mu ,H)–g_{PD}^2(\mu ,\mathcal{G}(\mu ,H),H)$$ Assume $f_c(\mu ,H)=m_1{\mu }^2+m_{2}H\mu +m_3{H}^2+O\left(|\mu +H{|}^3\right).$ By solving and comparing coefficients, we get